CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Numerical instability--- fvm::div environment

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By mAlletto

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 21, 2020, 14:45
Default Numerical instability--- fvm::div environment
  #1
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
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:

-\boldsymbol{v} \cdot \nabla (T) - \alpha \nabla^2 (T) = 0

when I implement it in OpenFOAM:

-This formulation is diverging:
Code:
 -fvm::div(phi, Ta) - fvm::laplacian(alpha_eff, Ta)

-This formulation is converging:
Code:
fvm::div(-phi, Ta) - fvm::laplacian(alpha_eff, Ta)

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.
andrea5 is offline   Reply With Quote

Old   February 22, 2020, 15:28
Default
  #2
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
mAlletto is offline   Reply With Quote

Old   February 22, 2020, 15:32
Default
  #3
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
Quote:
Originally Posted by mAlletto View Post
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.
andrea5 is offline   Reply With Quote

Old   February 22, 2020, 15:52
Default
  #4
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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.
mAlletto is offline   Reply With Quote

Old   February 23, 2020, 11:30
Default
  #5
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
Tom Lauriks likes this.
mAlletto is offline   Reply With Quote

Old   February 23, 2020, 14:46
Default
  #6
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
Hello, thanks for all the suggestions. I'll have a look on the documentation and the book you suggested.
Anyway, my scheme is:

Code:
  div(-phi,Ta)	    bounded Gauss upwind;
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.
andrea5 is offline   Reply With Quote

Old   February 24, 2020, 04:58
Default
  #7
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
mAlletto is offline   Reply With Quote

Old   March 10, 2020, 07:24
Default
  #8
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
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

-\boldsymbol{v} \cdot \nabla (T) - \alpha \nabla^2 (T) = 0

that traslated in OpenFOAM would become:

Code:
 -fvm::div(phi, Ta) - fvm::laplacian(alpha_eff, Ta)
I am trying to solve it with an upwind scheme, the solution diverges.

Changing it to:

Code:
 fvm::div(-phi, Ta) - fvm::laplacian(alpha_eff, Ta)
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:
Code:
 -fvm::div(phi, Ta) - fvm::laplacian(alpha_eff, Ta)
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
andrea5 is offline   Reply With Quote

Old   March 10, 2020, 09:22
Default
  #9
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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:


Code:
 -fvm::(phi,T)
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



Code:
 fvm::(-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
mAlletto is offline   Reply With Quote

Old   March 11, 2020, 11:27
Default
  #10
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
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
Code:
-fvm::(phi,T)
In this case what happens is that I am evaluating the quantity :
Code:
 fvm::(phi,T)
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
Code:
fvm::(-phi,T)
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
andrea5 is offline   Reply With Quote

Old   March 11, 2020, 11:55
Default
  #11
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
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
mAlletto is offline   Reply With Quote

Old   March 11, 2020, 12:02
Default
  #12
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
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:

Code:
fvm::(-phi,T)
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?
andrea5 is offline   Reply With Quote

Old   March 11, 2020, 12:30
Default
  #13
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15
mAlletto will become famous soon enough
A downwind scheme is unbounded for
Code:
fvm::div(phi,T) - fvm:: Laplacian (alpha,T)
your equation is
Code:
-fvm::div(phi,T) - fvm:: Laplacian (alpha,T)
so the stability conditions change.
mAlletto is offline   Reply With Quote

Old   March 16, 2020, 19:46
Default
  #14
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17
sharonyue is on a distinguished road
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.
__________________
My OpenFOAM algorithm website: http://dyfluid.com
By far the largest Chinese CFD-based forum: http://www.cfd-china.com/category/6/openfoam
We provide lots of clusters to Chinese customers, and we are considering to do business overseas: http://dyfluid.com/DMCmodel.html
sharonyue is offline   Reply With Quote

Old   March 17, 2020, 05:57
Default
  #15
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
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
andrea5 is offline   Reply With Quote

Old   March 17, 2020, 20:41
Default
  #16
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17
sharonyue is on a distinguished road
Quote:
Originally Posted by andrea5 View Post
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".
__________________
My OpenFOAM algorithm website: http://dyfluid.com
By far the largest Chinese CFD-based forum: http://www.cfd-china.com/category/6/openfoam
We provide lots of clusters to Chinese customers, and we are considering to do business overseas: http://dyfluid.com/DMCmodel.html
sharonyue is offline   Reply With Quote

Old   March 18, 2020, 05:27
Default
  #17
New Member
 
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6
andrea5 is on a distinguished road
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.
andrea5 is offline   Reply With Quote

Old   March 19, 2020, 22:33
Default
  #18
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17
sharonyue is on a distinguished road
Quote:
Originally Posted by andrea5 View Post
Hi,

why -fvc(phi,T,upwind) is unbounded?
Please see this image for why it is unbounded

__________________
My OpenFOAM algorithm website: http://dyfluid.com
By far the largest Chinese CFD-based forum: http://www.cfd-china.com/category/6/openfoam
We provide lots of clusters to Chinese customers, and we are considering to do business overseas: http://dyfluid.com/DMCmodel.html
sharonyue is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Numerical instability using "small" time step Pier84 OpenFOAM Running, Solving & CFD 7 May 6, 2020 12:08
TVD scheme Numerical Instability Tempa Main CFD Forum 0 June 29, 2011 02:54
Numerical Instability PattiMichelle Phoenics 2 December 26, 2005 20:49
NUMERICAL INSTABILITY OVER THE SWIRLING VANES Dan Main CFD Forum 2 February 9, 2005 09:12
Problems of numerical instability because of high gradients Xiangyang Ye Main CFD Forum 4 September 28, 1998 04:48


All times are GMT -4. The time now is 19:32.