[Ford Prefect]

Hi,

Problem summary:
A scalar field that should be bounded between 0-1 is not bounded.


Description:
I have implemented a scalar transport in OpenFOAM and get unbounded results. All boundary conditions are zero-flux except for one surface with specified value. The underlying transport equations are interFoam (VOF solver), and I solve this in a sequential manner. The diffusivity is dependent on the position (the volume-fraction). So far I have tested the following in order to avoid this.

1. Lower the convergence criterion for the residual transport.
2. Decrease the time-step size
3. Use higher order discretization in time (Crank-Nicolson)
4. Use upwind for the discretization in space.
5. Set the diffusivity to constant throughout the domain.

The only thing that really seems to work is to decrease the time-step size. In the base case the Courant number is around 0.2, which causes unbounded results. However, for a Courant number of 0.02 the result is ok.

The other (1, 3-5) actions only seem to help somewhat. For instance, if I use a limitedLinear scheme then obviously I get a more unbounded result, but upwind still causes problems.

The questions:

- Do you have any suggestion how to enforce bounded solutions?
- How does time-step size come into play for boundedness? Is it the linearization method?

Extra:

I have also tested to vary the scalar field between some other range, e.g. 100-200 and the resulting unbounded results still occurs at the same position in terms of significant digits.

[LuckyTran]

Have you tried using the 01 schemes, e.g. limitedLinear01, for hard bounding between 0 and 1.

[FMDenaro]

Originally Posted by Ford Prefect

Hi,

Problem summary:
A scalar field that should be bounded between 0-1 is not bounded.


Description:
I have implemented a scalar transport in OpenFOAM and get unbounded results. All boundary conditions are zero-flux except for one surface with specified value. The underlying transport equations are interFoam (VOF solver), and I solve this in a sequential manner. The diffusivity is dependent on the position (the volume-fraction). So far I have tested the following in order to avoid this.

1. Lower the convergence criterion for the residual transport.
2. Decrease the time-step size
3. Use higher order discretization in time (Crank-Nicolson)
4. Use upwind for the discretization in space.
5. Set the diffusivity to constant throughout the domain.

The only thing that really seems to work is to decrease the time-step size. In the base case the Courant number is around 0.2, which causes unbounded results. However, for a Courant number of 0.02 the result is ok.

The other (1, 3-5) actions only seem to help somewhat. For instance, if I use a limitedLinear scheme then obviously I get a more unbounded result, but upwind still causes problems.

The questions:

- Do you have any suggestion how to enforce bounded solutions?
- How does time-step size come into play for boundedness? Is it the linearization method?

Extra:

I have also tested to vary the scalar field between some other range, e.g. 100-200 and the resulting unbounded results still occurs at the same position in terms of significant digits.

The Godunov theorem states that only first order linear scheme is monotone.
That is not an issue caused directly by the time step but is due to the requirement that all coefficient must be positive.
For a convection-diffusion scalar equation you can have an oscillation-free solution using a cell Peclet <2.
A further task requires introducing a non-linear scheme also for a linear problem.
For specific question on how using OF you should post in the proper subsection.

[Ford Prefect]

Originally Posted by LuckyTran

Have you tried using the 01 schemes, e.g. limitedLinear01, for hard bounding between 0 and 1.

No, this looks promising. I did not see that the alpha (VOF) function is using this scheme though. Is it fully conservative? I will have a look at the source code when I get back home, thank you!

Originally Posted by FMDenaro

The Godunov theorem states that only first order linear scheme is monotone.
That is not an issue caused directly by the time step but is due to the requirement that all coefficient must be positive.
For a convection-diffusion scalar equation you can have an oscillation-free solution using a cell Peclet <2.
A further task requires introducing a non-linear scheme also for a linear problem.
For specific question on how using OF you should post in the proper subsection.

Yes, and upwind was my first go-to test in order to enforce boundedness. As per my initial post .That did not work. However, changing the time-step did work. The scalar transport equation is in fact fully linear when I think about it, since I have also tested the problem for constant diffusivity.


Can it be related to the gradient calculation for the unstructured mesh?



The problem is not OpenFOAM specific. I can create a similar case in Fluent.

[sbaffini]

The ideal bounded setting is, or should be, first order Euler in time, first order upwind for convection and, if the grid is non orthogonal, zero non-orthogonal corrections for the diffusion part. But I am no expert in openfoam, and there certainly are much more points, in a code, where unboundedness can get in. In general, a code is made in such a way to not do that if all the main discretization strategies are first order, but I can't speak for OF.

In general, even with a first order discretization, depending from the specific equation form in use, a scalar can become unbounded if the underlying continuity equation is not exactly satisfied. That error acts as source term, for example, in the conservative formulation.

Have you tried an equivalent steady problem on the same grid, just with the scalar?

[Ford Prefect]

Originally Posted by sbaffini

In general, even with a first order discretization, depending from the specific equation form in use, a scalar can become unbounded if the underlying continuity equation is not exactly satisfied. That error acts as source term, for example, in the conservative formulation.

This is probably worth investigating. You mean that the incompressibility condition is not sufficiently satisfied and thus there is a possibility of unbounded solutions since the velocity field is not consistent? This could also explain why a smaller time-step is better in my case. I have not explicitly looked at the residuals but since OF has a minimum number of iterations for each variable it is possible that the residuals are much lower for small time-steps.



I will test to increase the demands on global continuity, thank you!


EDIT:







So if the equations are written in conservation form I believe your suggestion may be very important, depending on the sensitivity of the equations. The residual in the continuity equation is directly affecting the boundedness of the scalar equation. In all "text-book examples" as far as I can tell, which does not account for the solution of the velocity field.

[FMDenaro]

Originally Posted by Ford Prefect

No, this looks promising. I did not see that the alpha (VOF) function is using this scheme though. Is it fully conservative? I will have a look at the source code when I get back home, thank you!






Yes, and upwind was my first go-to test in order to enforce boundedness. As per my initial post .That did not work. However, changing the time-step did work. The scalar transport equation is in fact fully linear when I think about it, since I have also tested the problem for constant diffusivity.


Can it be related to the gradient calculation for the unstructured mesh?



The problem is not OpenFOAM specific. I can create a similar case in Fluent.

From a theoretical point of view, all coefficients should be positive in a first order scheme. The diffusive term is second order accurate, I suppose but if the cell Peclet number is O(1) it should work.



The first test I would check is resolving the transport equation without diffusivity using the first order upwind. In such a way no gradients are involved and the solution should be bounded. Try that and let us know.

[Ford Prefect]

Thanks to @sbaffini, I have now tested to lower the continuity (p_rgh) tolerance from 1e-7 to 1e-12 and for the current case this eliminated the problem with unbounded scalars.


I can run it with any time-step that the momentum solution admits and it is also possible to use higher order schemes for the scalar transport (e.g. limitedLinear) without getting unbounded results.


So in conclusion: it is extremely important to have a tight convergence for the continuity equation when dealing with scalar transport, which in hindsight seems perfectly reasonable.



I find it curious that this is not mentioned more often for multi-phase models since the volume-fraction should be bounded between 0 and 1.


Anyways, case closed (for now). Thank you for all valuable input!

[sbaffini]

Originally Posted by Ford Prefect

I find it curious that this is not mentioned more often for multi-phase models since the volume-fraction should be bounded between 0 and 1.

I learned this myself no more than a couple of months ago while testing my code for passive scalars and (as I use an immersed boundary code, which generally leaks at boundaries) the new algorithm to fix leaks. The most spectacular way to see this is to perform a computation with the scalar at convergence (say, initialized with the same inlet value and with zero flux bcs at walls) but not the velocity field, and then monitor the max and min scalar values in the domain. Enlightening. I am not into VOF, like at all, but that's the first thing that came into my mind and, like you, I was: "Wow, you guys could have written that somewhere, like in any textbook, when you present the non conservative form of the VOF equation".

If you look at the equations, it is trivial:



and indeed mentioned often out of context. But I have seen things much more trivial repeated at nauseam in many textbook, so, in the end, I'm pretty sure that's not why it is not mentioned, considering how peculiar the VOF equation is.

EDIT: Another lesson I learned was that I can't reliably test scalars using only 0-1 values, as the equation above wouldn't go below a minimum 0 but would go below any other non 0 minimum if the continuity is not perfectly satisfied.