CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Gravity source term causes instability?

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

Like Tree1Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 20, 2023, 22:53
Default Gravity source term causes instability?
  #1
New Member
 
Join Date: May 2023
Posts: 15
Rep Power: 2
ALVRM is on a distinguished road
I have built a pressure projection FVM solver for incompressible flow, and it works as intended if I use a small gravitational force, like 0.01 * 9.81. But when I use a larger value, like the actual acceleration of gravity of 9.81, the solver becomes very unstable and blows up.

I built the code without any third part libraries, so my question is conceptual... I would like to know what stability condition I am missing... My code is fully implicit for both the advection and diffusion term. I was using Crank Nicolson semi-implicit but realized that fully implicit is more stable.

I tried a pressure correction approach and it still failed at 9.81 source value. I might try to damp the pressure corrector, as I think that the large source term is causing the pressure to over correct and leading to large errors.

Any suggestions would be helpful. My ultimate goal is to create an unconditionally stable incompressible solver for natural convection.
ALVRM is offline   Reply With Quote

Old   May 21, 2023, 10:31
Default
  #2
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,120
Blog Entries: 29
Rep Power: 38
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Do you have variable density and/or temperature?

Are you using a boussinesq like term or full gravity?

Usually, for full gravity, one splits it in a constant density term absorbed by the pressure and a variable density one, which should be more benign. Fully coupled scheme then make it implicit as well
sbaffini is offline   Reply With Quote

Old   May 21, 2023, 13:16
Default
  #3
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,598
Rep Power: 70
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by ALVRM View Post
I have built a pressure projection FVM solver for incompressible flow, and it works as intended if I use a small gravitational force, like 0.01 * 9.81. But when I use a larger value, like the actual acceleration of gravity of 9.81, the solver becomes very unstable and blows up.

I built the code without any third part libraries, so my question is conceptual... I would like to know what stability condition I am missing... My code is fully implicit for both the advection and diffusion term. I was using Crank Nicolson semi-implicit but realized that fully implicit is more stable.

I tried a pressure correction approach and it still failed at 9.81 source value. I might try to damp the pressure corrector, as I think that the large source term is causing the pressure to over correct and leading to large errors.

Any suggestions would be helpful. My ultimate goal is to create an unconditionally stable incompressible solver for natural convection.



There are many possible issues in the projection method for bouyancy-driven flows. I assume that you validated carefully your code in standard solution without bouyancy.

You should detail how you solve the system, are you solving the internal energy equation with the Bousinnesq coupling?
How do you discretize the source term in time?
Is the non linear term linearized in the implicit scheme?
FMDenaro is offline   Reply With Quote

Old   May 21, 2023, 15:38
Default
  #4
New Member
 
Join Date: May 2023
Posts: 15
Rep Power: 2
ALVRM is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
Do you have variable density and/or temperature?

Are you using a boussinesq like term or full gravity?

Usually, for full gravity, one splits it in a constant density term absorbed by the pressure and a variable density one, which should be more benign. Fully coupled scheme then make it implicit as well

I am using a boussinesq like term:

\frac{(\rho_{cell} - \rho_{rel})}{\rho_{cell}} 9.81

The density is calculated at each time step based on temperature, using the equation of state for air. I solve an energy equation for the temperature as part of the procedure to get the density.

Can you describe what you mean by "make it implicit as well"?
ALVRM is offline   Reply With Quote

Old   May 21, 2023, 16:02
Default
  #5
Senior Member
 
Lucky
Join Date: Apr 2011
Location: Orlando, FL USA
Posts: 5,493
Rep Power: 64
LuckyTran has a spectacular aura aboutLuckyTran has a spectacular aura aboutLuckyTran has a spectacular aura about
You have now a source term which may be discretized differently than advection terms or diffusion terms. Is \rho_{cell} also implicit?

Gravity waves do go into the Courant number and can influence stability but somehow I doubt that is the issue you have having.
LuckyTran is offline   Reply With Quote

Old   May 21, 2023, 16:04
Default
  #6
New Member
 
Join Date: May 2023
Posts: 15
Rep Power: 2
ALVRM is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
There are many possible issues in the projection method for bouyancy-driven flows. I assume that you validated carefully your code in standard solution without bouyancy.

You should detail how you solve the system, are you solving the internal energy equation with the Bousinnesq coupling?
How do you discretize the source term in time?
Is the non linear term linearized in the implicit scheme?
I solve a temperature equation explicitly:

\frac{\partial T}{\partial t} = -u \nabla T + \alpha \nabla^{2}T + Q

Then I use the equation of state for air to get a density field.

Then the fully implicit momentum equation:

\frac{\partial u^{*}}{\partial t} + u^{n}\nabla u^{*} - \nu \nabla^{2}u^{*} = \frac{\Delta\rho}{\rho}  g

I make the convection term linear by using the previous time step velocity.

In matrix form it looks like this:

\left[ \frac{1}{\Delta t} \right] \left[ u^{*} \right]
+
\left[ u^{n}\nabla \right] \left[ u^{*} \right]
+
\left[ -\nu \nabla^{2} \right] \left[ u^{*} \right]
=
\left[ \frac{\Delta\rho}{\rho}  g + \frac{ u^{n}}{\Delta t} \right]

I then solve the pressure poisson equation with zero divergence for the u^{n+1} term, and solve for u^{n+1} using the gradient pressure.
ALVRM is offline   Reply With Quote

Old   May 21, 2023, 16:24
Default
  #7
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,598
Rep Power: 70
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Why do you use the state equation?
When you solve the temperature field, the buoyancy term in the momentum is expressed directly by the temperature.

Be also careful, the time step you must use is dictated by the explicit scheme in the temperature equation, your implicit scheme cannot use an arbitrary high time step. The stability in the temperature equation depends on the Peclet and cfl but for the momentum equation the source term has influence on the stability.

Before introducing the bouyancy term what kind of test case have you resolved to assess the omothermal case?
FMDenaro is offline   Reply With Quote

Old   May 21, 2023, 16:26
Default
  #8
New Member
 
Join Date: May 2023
Posts: 15
Rep Power: 2
ALVRM is on a distinguished road
Quote:
Originally Posted by LuckyTran View Post
You have now a source term which may be discretized differently than advection terms or diffusion terms. Is \rho_{cell} also implicit?

Gravity waves do go into the Courant number and can influence stability but somehow I doubt that is the issue you have having.
Can you explain how to make \rho implicit? Or point me to a reference that explains this?
ALVRM is offline   Reply With Quote

Old   May 21, 2023, 16:31
Default
  #9
New Member
 
Join Date: May 2023
Posts: 15
Rep Power: 2
ALVRM is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Why do you use the state equation?
When you solve the temperature field, the buoyancy term in the momentum is expressed directly by the temperature.

Be also careful, the time step you must use is dictated by the explicit scheme in the temperature equation, your implicit scheme cannot use an arbitrary high time step. The stability in the temperature equation depends on the Peclet and cfl but for the momentum equation the source term has influence on the stability.

Before introducing the bouyancy term what kind of test case have you resolved to assess the omothermal case?
I have two stability constraints, they come from the advection and diffusion terms in the temperature equation. The issue is the gravity source term. It works for small source values...

From what I understand, the bousenessq term is derived from the equation of state anyway. I just use density variation instead of temperature variation. I can substitute the actual bousenessq term in my equation and the same results occur.
ALVRM is offline   Reply With Quote

Old   May 22, 2023, 06:12
Default
  #10
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,120
Blog Entries: 29
Rep Power: 38
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by ALVRM View Post
I am using a boussinesq like term:

\frac{(\rho_{cell} - \rho_{rel})}{\rho_{cell}} 9.81

The density is calculated at each time step based on temperature, using the equation of state for air. I solve an energy equation for the temperature as part of the procedure to get the density.

Can you describe what you mean by "make it implicit as well"?
First, you can split the gravity term in:

\rho \textbf{g} = \nabla \left(\rho_0 \textbf{g} \cdot \textbf{r}\right) +\left(\rho - \rho_0\right) \textbf{g}

where the first term, with a suitably chosen value for the reference density \rho_0, can be included in the pressure gradient. The second term then will be less troublesome.


Second, a coupled solver will typically handle all the equations together. For example, a preconditioned density based code has p, u, v, w, T as vector of independent variables and the matrix coefficients are 5x5 blocks, so a momentum source term depending on density has an explicit dependence from the independent variables and I can treat it implicitly with no particular reasoning.

In your case, density is function of T, whose equation is not solved with momentum but segregated. The stability analysis here is more complex, and I have no direct experience, but solving for the temperature first should already give you a better stability.
FMDenaro likes this.
sbaffini is offline   Reply With Quote

Old   May 22, 2023, 07:28
Default
  #11
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,120
Blog Entries: 29
Rep Power: 38
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by ALVRM View Post
From what I understand, the bousenessq term is derived from the equation of state anyway. I just use density variation instead of temperature variation. I can substitute the actual bousenessq term in my equation and the same results occur.
They are different. With Boussinesq, your density is constant everywhere except the gravity term, especially relevant being the pressure equation. Also, your EOS in the gravity term is fixed to a specific model of linear dependence of density from temperature. They are equal only if you make all the Boussinesq assumptions but just retain the density form (which, at that point, is useless).

With Boussinesq you can use a constant density code. Without it, you need a variable density one.
sbaffini is offline   Reply With Quote

Old   May 22, 2023, 08:04
Default
  #12
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,598
Rep Power: 70
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Solve first for the temperature, then you have the temperature solutions both at tn and tn+1 to insert into the momentum by means of Bousinnesq.

However, check that the stabilty constraint is satisfied for the temperature equation.
FMDenaro is offline   Reply With Quote

Old   May 23, 2023, 01:35
Default
  #13
Member
 
EM
Join Date: Sep 2019
Posts: 53
Rep Power: 6
gnwt4a is on a distinguished road
Quote:
Originally Posted by ALVRM View Post
I would like to know what stability condition I am missing...

*important*: since u have incompressible flow, your buoyancy term *must* obey the pressure compatibility condition otherwise your discrete system has no solution. once fulfilled, the B-V frequency below is the next consideration.

Brunt-Vaisala instability. Even if u do not explicitly impose stable stratification, thermal flows will develop regions (at least transiently) where these occur. There is no steady state solution and codes will never converge to a steady-state. If timestepping is used, u must time-resolve those gravity waves or up it goes. these numerical effects have been known for quite a while.

Last edited by gnwt4a; May 23, 2023 at 03:43. Reason: Added compatibility precondition
gnwt4a is offline   Reply With Quote

Old   May 23, 2023, 03:22
Default
  #14
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,598
Rep Power: 70
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by gnwt4a View Post
Brunt-Vaisala instability. Even if u do not explicitly impose stable stratification, thermal flows will develop regions (at least transiently) where these occur. There is no steady state solution and codes will never converge to a steady-state. If timestepping is used, u must time-resolve those gravity waves or up it goes. these numerical effects have been known for quite a while.
Actually, you have to consider in the code that a physical rest condition must be saved for a stable temperature distribution T(z). Thus, you must introduce a gravity forcing only for the difference to such condition.
FMDenaro is offline   Reply With Quote

Old   May 23, 2023, 07:24
Default
  #15
New Member
 
Join Date: Nov 2013
Posts: 5
Rep Power: 11
jpanchog is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
They are different. With Boussinesq, your density is constant everywhere except the gravity term, especially relevant being the pressure equation. Also, your EOS in the gravity term is fixed to a specific model of linear dependence of density from temperature. They are equal only if you make all the Boussinesq assumptions but just retain the density form (which, at that point, is useless).

With Boussinesq you can use a constant density code. Without it, you need a variable density one.

Regarding your last statement: is that really so? shouldn't be also possible to use an EoS for the buoyancy term rho*g in an incompressible (constant density) setting?
I guess the results should be pretty similar (if not the same) if a incompressible solver is used together with an EoS, or a linearized form of it... at least for very small temperature differenecs
jpanchog is offline   Reply With Quote

Old   May 23, 2023, 08:47
Default
  #16
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,120
Blog Entries: 29
Rep Power: 38
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by jpanchog View Post
Regarding your last statement: is that really so? shouldn't be also possible to use an EoS for the buoyancy term rho*g in an incompressible (constant density) setting?
I guess the results should be pretty similar (if not the same) if a incompressible solver is used together with an EoS, or a linearized form of it... at least for very small temperature differenecs

Boussinesq is not an arbitrary EOS, it is, specifically, a linearization of a density depending only from temperature around a given temperature. It is arbitrary in the thermal expasion coefficient, indeed, not the rest. That's where your freedom is in using Boussinesq. So, if you are using a general density form but only in the gravity term, it certainly is a mistake. Is it the cause of your instability? I don't know. But it might be.
sbaffini is offline   Reply With Quote

Old   May 23, 2023, 10:48
Default
  #17
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,598
Rep Power: 70
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by jpanchog View Post
Regarding your last statement: is that really so? shouldn't be also possible to use an EoS for the buoyancy term rho*g in an incompressible (constant density) setting?
I guess the results should be pretty similar (if not the same) if a incompressible solver is used together with an EoS, or a linearized form of it... at least for very small temperature differenecs



And how do you think you can compute rho=p/RT ? How do you evaluate "p" in the incompressible flow model?
FMDenaro is offline   Reply With Quote

Old   May 23, 2023, 12:32
Default
  #18
New Member
 
Join Date: Nov 2013
Posts: 5
Rep Power: 11
jpanchog is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
And how do you think you can compute rho=p/RT ? How do you evaluate "p" in the incompressible flow model?


I'm really not sure if it makes sense, but wouldn't it be possible to just use for the pressure appearing in the EoS a constant pressure (the thermodynamic pressure) , as it is done in the low-Mach formulation of the governing equations?
jpanchog is offline   Reply With Quote

Old   May 23, 2023, 12:47
Default
  #19
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,598
Rep Power: 70
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by jpanchog View Post
I'm really not sure if it makes sense, but wouldn't it be possible to just use for the pressure appearing in the EoS a constant pressure (the thermodynamic pressure) , as it is done in the low-Mach formulation of the governing equations?
You donít have any thermodynamics law in the incompresdible flow model.
FMDenaro is offline   Reply With Quote

Old   May 23, 2023, 14:32
Default
  #20
New Member
 
Join Date: May 2023
Posts: 15
Rep Power: 2
ALVRM is on a distinguished road
Quote:
Originally Posted by gnwt4a View Post
*important*: since u have incompressible flow, your buoyancy term *must* obey the pressure compatibility condition otherwise your discrete system has no solution. once fulfilled, the B-V frequency below is the next consideration.

Brunt-Vaisala instability. Even if u do not explicitly impose stable stratification, thermal flows will develop regions (at least transiently) where these occur. There is no steady state solution and codes will never converge to a steady-state. If timestepping is used, u must time-resolve those gravity waves or up it goes. these numerical effects have been known for quite a while.
For the pressure poisson solver, I use the boundary condition of u^{*} = 0 and \frac{\partial P}{\partial \hat{n}} = 0. I set an arbitrary pressure at one point, P = 0, and I always get a unique solution for the poisson system of equations. I use P = 0 because I'm only interested in the pressure gradient, not absolute pressure. This is typical for walls, correct? I tried to find literature on "pressure compatibility condition" and I cannot find anything useful. Do you have a source for this?

Brunt-Vaisala instability might be what I am experiencing... I looked at the Wikepedia article for it and I understand the general concept. Is there a way to eliminate this stability condition?

Last night, I ran the simulation for 12 hours (around 100,000 iterations & 1,072 simulated flow seconds). It's perfectly stable with a body force term of (\rho - \rho_{ref}) 0.0981. By keeping the gravity term low, I notice the lack of "swirling" or curling in the simulation. When I bump up the gravity term to the correct 9.81 value, I get swirling effects and eventual very high velocities that blow up. The simulation usually fails within 30 seconds.
ALVRM is offline   Reply With Quote

Reply

Tags
finite volume method, implicit, incompressible flow, pressure projection

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
[foam-extend.org] Problems installing foam-extend-4.0 on openSUSE 42.2 and Ubuntu 16.04 ordinary OpenFOAM Installation 19 September 3, 2019 18:13
[Other] How to use finite area method in official OpenFOAM 2.2.0? Detian Liu OpenFOAM Meshing & Mesh Conversion 4 November 3, 2015 03:04
[swak4Foam] Swak4FOAM 0.2.3 / OF2.2.x installation error FerdiFuchs OpenFOAM Community Contributions 27 April 16, 2014 15:14
[swak4Foam] Error bulding swak4Foam sfigato OpenFOAM Community Contributions 18 August 22, 2013 12:41
[swak4Foam] build problem swak4Foam OF 2.2.0 mcathela OpenFOAM Community Contributions 14 April 23, 2013 13:59


All times are GMT -4. The time now is 07:58.