|
[Sponsors] |
Numerical instability--- fvm::div environment |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
February 21, 2020, 14:45 |
Numerical instability--- fvm::div environment
|
#1 |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
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: 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. |
|
February 22, 2020, 15:28 |
|
#2 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 |
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
|
|
February 22, 2020, 15:32 |
|
#3 | |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
Quote:
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. |
||
February 22, 2020, 15:52 |
|
#4 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 |
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.
|
|
February 23, 2020, 11:30 |
|
#5 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 |
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 |
|
February 23, 2020, 14:46 |
|
#6 |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
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; |
|
February 24, 2020, 04:58 |
|
#7 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 |
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 |
|
March 10, 2020, 07:24 |
|
#8 |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
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: Code:
-fvm::div(phi, Ta) - fvm::laplacian(alpha_eff, Ta) Changing it to: Code:
fvm::div(-phi, Ta) - fvm::laplacian(alpha_eff, Ta) 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) 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 |
|
March 10, 2020, 09:22 |
|
#9 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 |
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) Code:
fvm::(-phi,T) Hope the explanation classifies some things Michael |
|
March 11, 2020, 11:27 |
|
#10 |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
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) Code:
fvm::(phi,T) 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) 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 |
|
March 11, 2020, 11:55 |
|
#11 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 |
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
|
|
March 11, 2020, 12:02 |
|
#12 |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
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) |
|
March 11, 2020, 12:30 |
|
#13 |
Senior Member
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 615
Rep Power: 15 |
A downwind scheme is unbounded for
Code:
fvm::div(phi,T) - fvm:: Laplacian (alpha,T) Code:
-fvm::div(phi,T) - fvm:: Laplacian (alpha,T) |
|
March 16, 2020, 19:46 |
|
#14 |
Senior Member
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17 |
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 |
|
March 17, 2020, 05:57 |
|
#15 |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
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 |
|
March 17, 2020, 20:41 |
|
#16 | |
Senior Member
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17 |
Quote:
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 |
||
March 18, 2020, 05:27 |
|
#17 |
New Member
Andrea Trotta
Join Date: Aug 2019
Posts: 14
Rep Power: 6 |
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. |
|
March 19, 2020, 22:33 |
|
#18 |
Senior Member
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17 |
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 |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
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 |