[wildfire230]

Hi all,

I am using a modified version of scalarTransportFoam which allows for variable diffusivity (although the problem I am encountering is also present for constant diffusivity). First, I solve the steady state flow profile using simpleFoam, which shows excellent convergence, approaches velocity profile I expect, and shows very small divergence errors, i.e. max(fvc::div(phi)) < 10^(-10).

This velocity field includes some recirculation near the boundary, where I impose a fixed uniform value of the scalar. Due to this convective term, my scalar concentration goes negative near the boundary, which is unphysical. I am trying to figure out exactly what numerical schemes and solvers I should be using for this problem. Currently I am using:
C
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-10;
relTol 0;
}

where C is my scalar concentration, and
div(phi,C) Gauss linearUpwind Gauss linear;
but I have also tried many alternatives including bounded Gauss upwind, Gauss linear, Gauss upwind, etc...

I would really appreciate any help!

Thanks

[invadoria]

I recommend you to use limited scheme for the divergence term to limit the scalar variable C in the range of 0 to 1. Otherwise negative values may be encountered due to the truncation errors arising from the discretization scheme.

[wildfire230]

Awesome, thanks for this suggestion. I didn't know there were specific schemes for this case, but from what I can gather from the documentation it seems like:

div(phi,C) Gauss limitedLinear01 1;

is the proper choice. Does that make sense? It does seem to result in smaller negative values for my scalar concentration, but it still doesn't strictly enforce positivity. Is there something I'm missing?

[invadoria]

I had used Gamma01 method to strictly bound the scalar variable.

[wildfire230]

Thanks, I didn't find the Gamma01 method in the fvSchemes documentation, but when I try:

div(phi,C) Gauss Gamma01 1;

in my divSchemes, I still find values of C that are negative. Am I still missing something?

Edit: following hjasak's instructions here: Schemes for scalar in LES I also switched to

laplacian(D,C) Gauss linear limited 1;

and switched to Euler for ddt scheme. I still end up with negative values for my C. Also they are not trivially small values, but almost the same order of magnitude as my positive values.

Thanks again

[akidess]

Post your code.

[wildfire230]

Here is an example code that illustrates the problem. I am using OpenFOAM 5.0. Simply unzip the directory, run blockMesh, and then simpleFoam, which will converge in about 0.525 seconds. Then copy the T initial condition into the converged directory: "cp 0/T 0.525". Then run scalarTransportFoam. If you use paraView to look at the final result, you will see that there are negative temparature values, which are non-physical. They're not that dramatic, but they also are not negligible, and I cannot eliminate them regardless of which fvSchemes I choose.

I would greatly appreciate any help. This same problem has been a hindrance on multiple different simulations at this point.

Thanks!

Edit: Also I should note that these negative values do not appear to go away with a more refined grid.

[wildfire230]

I haven't done a comprehensive set of tests, but these negative values from scalarTransportFoam do not appear to go away with a finer mesh or smaller time step. I think it must be an issue with the selection of schemes, but I do not understand the details well enough.

[wildfire230]

Hi all, I am still trying to solve this problem. Does anyone know the requirements on grid size for stability of the scalarTransportFoam solver for a given diffusivity coefficient and velocity field? For example, here: https://en.wikipedia.org/wiki/Numeri...usion_equation there is a condition listed as h < 2a/u, where a is the diffusivity coefficient, h is the grid size, and u is the fluid velocity. However, I believe this is specifically for the case of forward Euler centered difference.


Does anyone know the restrictions on grid size for scalarTransportFoam for a given diffusivity coefficient. Specifically, in my case the diffusivity is very small, so that the Peclet number is large.

[randolph]

Originally Posted by wildfire230

Hi all, I am still trying to solve this problem. Does anyone know the requirements on grid size for stability of the scalarTransportFoam solver for a given diffusivity coefficient and velocity field? For example, here: https://en.wikipedia.org/wiki/Numeri...usion_equation there is a condition listed as h < 2a/u, where a is the diffusivity coefficient, h is the grid size, and u is the fluid velocity. However, I believe this is specifically for the case of forward Euler centered difference.


Does anyone know the restrictions on grid size for scalarTransportFoam for a given diffusivity coefficient. Specifically, in my case the diffusivity is very small, so that the Peclet number is large.

Hi,

I have encountered the same issue. I am having the scalar negative. I think its the solver issue more than the scheme. I am planning dig into the multiphase solver to check the implementation. Have you gotta any updates?

Thanks,
Rdf

[piu58]

> This velocity field includes some recirculation near the boundary, where I impose a fixed uniform value of the scalar. Due to this convective term, my scalar concentration goes negative near the boundary, which is unphysical.

I assume this is an outlet with a zero gradient b.c. for velocity.

If you give zero gradient than it has to be there. The simulation is overstrained to enforce this. If possible, extend the outflow area to give the flow room for homogenizing.

[randolph]

After some playing around with the code, it looks like this problem is associated with the boundary condition when the flow is pointing outwards. If you apply a fixedValues or a fixedGradient BC (centrally a fixedValues types too), the outgoing flux on that face will be calculated with the values that you specified which potentially will be larger than the quantity in that cell and produce a negative value. This is also should be the reason why there is the inletOutlet BC in the OpenFOAM (zeroGradient for outflow and fixedValue for inflow). Someone correct me if I am wrong.

Now for the most of the problem that you can avoid the inflow (reverse flow), you can just take care of this problem when you compute the flow field.

[Hughtong]

Originally Posted by randolph

Hi,

I have encountered the same issue. I am having the scalar negative. I think its the solver issue more than the scheme. I am planning dig into the multiphase solver to check the implementation. Have you gotta any updates?

Thanks,
Rdf

Hi, I'm dealing with a periodic channel in the streamwise direction and periodic temperature BCs are applied, and the results got the negative, which is totally non-physical, but with the process of nondimensionless, the profile of temperature with the function of yplus can match with DNS database