Numerical instability--- fvm::div environment
 

Hello Foamers,

hope it's the right section to post.
Recently I faced a problem while implementing an advection-diffusion equation. Basically the equation is:



when I implement it in OpenFOAM:

-This formulation is diverging:

-This formulation is converging:

Basically just the minus sign is inside the brackets. I guess the problem is related on how OpenFOAM is discretizing the matrices, but I was wandering if anyone could give me a more detailed explanation on what happens inside the openFOAM codes.

Thanks, and have a nice day,

Andrea.

Are you sure the divergence and Laplacian should have the seme sign? I'm not sure if there exists a stable solution for your equation. The discretion depends on the scheme you chose

Hi, thanks for replying.
I am sure it's correct, it's not the classical advection-difffusion you would implement for a real flow. It comes from the derivation of an adjoint method formulation, for that reason it appears a minus in both. It has no physical sense, but it's used to compute the sensitivity for the adjoint equations.

Thanks, Andrea.

A ok. For a central scheme div(-phi,T) and -div(phi,T) should give the same result. For an upwind biased scheme the two formulation are different. You can have a look here https://www.openfoam.com/documentati...ivergence.html.

A good book describing how the different operators in OF a discretized is https://www.springer.com/de/book/9783319168739. Different surface interpolation schemes used to compute the value at the face center in order to dicretize the divergence in a finite volume framework can be found in $FOAM_SRC/finiteVolume/interpolation/surfaceInterpolation


Hope this helps

Hello, thanks for all the suggestions. I'll have a look on the documentation and the book you suggested.
Anyway, my scheme is:

so as you said, is not centered, thus the solution is different in the 2 cases. If you have any explanation on this particular behavior I would be glad to know.

So in Chapter 11.2 in the book I suggest to read, you find that for a 1D convective diffusion equation a upwind scheme is always bounded for any Pe number, a central scheme becomes unstable at Pe > and a downwind solution is completely unbounded at Pe = 4.



In your case however you have a convective diffusion equation with a negative diffusivity. By doing div(-phi,T) you actually transform the upwind scheme to a downwind scheme. So maybe for your equation a downwind scheme is always bounded. You can try to repeat the analysis of the book and see the outcome. I haven't done it by myself.


Best


Michael

Hi, sorry for the silence.

I had a look to the book you suggested me, so to sum up we can say that:

- This is the equation that I am trying to solve



that traslated in OpenFOAM would become:

I am trying to solve it with an upwind scheme, the solution diverges.

Changing it to:

what we do is actually removing the minus sign and incorporate it into the values of phi. With the same upwind scheme is converging.

From the math point of view, now the equation is again a "common advection-diffusion equation" with coherent signs.

You were suggesting me that a downwind solution is completely unbounded at Pe = 4. while an upwind is bounded for any Pe. In the case where:
my upwind scheme was changing to a downwind scheme due to the minus sign in front of the advection term, thus the solution becomes unbounded. Is that right what you were suggesting me?

In the case the answer is yes, I would like to ask why the minus changes the upwind to a downwind, I mean, isn't it depending on the spatial point I choose for discretizing the derivative of the advection part?

thanks for your time,

Andrea

Hello Andrea,


so in an upwind scheme it is decided based on the sign of phi which cell center value is used to compute the face quantity. The cell center value is used from the cell very the flux come from. So if you write the eqution:


the cell enters from which the flux comes from are used to calculate phi * T at the faces. After that the whole thing is multiplied by -1. You get then -phi * T



the cell enters where the flux is directed are used to calculate -phi * T at the faces. The reason is that you invert the sign of the flux before you evaluate the upwind scheme.



Hope the explanation classifies some things


Michael

Hi Michael,

thanks again for answering. I think I understood your point and I will resume it to have your assurance on my understanding, and also to help some interested reader.
Referring to the book you suggested me, from page 366 to 380, where there is the analysis on the upwind and downwind scheme.

Using the same nomenclature: phi_W --- phi_w --- phi_C --- phi_e --- phi_E
where capital letters stands for the cell centers and phi_w and phi_e are the values at the faces where phi is evaluated.


CASE 1
In this case what happens is that I am evaluating the quantity :
where, using an upwind scheme I obtain :

phi_e = phi_C and phi_w = phi_W

then of course I use these values multiplied by T and I obtain my advection term.

Later everything is multiplied by -1, thus what was an upwind scheme, in some manner is seen by the physic of the problem as a downwind scheme (an upwind scheme used with "negative velocity"). Thus the solution diverges since it seems to be unbounded.

CASE 2
in this case the minus sign is "absorbed" by the value of the phi, thus with the same procedure described above I obtain an upwind scheme that is coherent with the advection velocity and thus always bounded.



I know that my explanation is far from being clear, but I hope that I made the point of the discussion.

Thanks for you time,

Andrea

Hm i think for the case fvm::(-phi,T) you have a downwind scheme which is in your case stable since your advection diffusion equation had a minus in front of the advection term

Maybe I missed that in openFOAM in both cases I use an upwind discretization scheme. May I ask why it would be a downwind scheme? I mean, I am using an upwind as discretization, but with a minus inside:

So I can agree that results in a downwind scheme at the end. But can you explain me once again why it results as a stable solution? evenif I studied on the book that a downwind Is unbounded as Pe approaches to 4?

A downwind scheme is unbounded for your equation is so the stability conditions change.

Andrea, there are several reasons that your results blow up.

1. For the advection, only a first order scheme (upwind/downwind) produces a diagonally equal matrix system. Meanwhile, fvc::div(phi,T,upwind) = -fvc::div(-phi,T,downwind), fvc::div(-phi,T,upwind) = -fvc::div(phi,T,downwind) so please try -fvc::div(phi,T,downwind). However, I prefer fvc::div(-phi,T,upwind) since it looks more physical.

2. You are solving a steady-state equation. The laplacian term produces a diagonally equal matrix for orthogonal meshes. If your advection uses any kind of high order scheme. It violates the diagonal equality and create unbounded values of T. Try relax this system to boost diagonal dominance.

Hi,

Thanks for replying. I don't think OpenFOAM has a downwind scheme already implemented, so I'll try to implement it and I'll try to see if -fvc(phi,T, downwind) works as wells as fvc(-phi,T,upwind).

In any case I have only first order schemes for all the terms in the equations, so it should work.

In any case the problem was that with -fvc(phi,T,upwind) I was trying to solve an unphysical problem?

Thanks,

Andrea

If you solve UT, you should use div(phi,T), or -div(phi_neg,T,DOWNWIND).
If you solve -UT, you should implement div(phi_neg,T), or -div(phi_pos,T,DOWNWIND).

In any way the first one makes more sense, since the latter one can only use downwind scheme. Downwind is implemented by OpenFOAM by keyword "downwind".

Hi,

thanks for the clarification.
In any case my question is still why this phenomena happens? why -fvc(phi,T,upwind) is unbounded? Is that related to a physical meaning or is as MAletto was suggesting a problem about the stability of the advection-diffusion equation?

Sorry if I ask this again, but I need to be sure on the inner cause of the problem.

Please see this image for why it is unbounded

http://www.cfd-china.com/assets/uplo...0320103052.jpg