Gravity source term causes instability?

ALVRM · May 20, 2023, 23:53

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.

sbaffini · May 21, 2023, 11:31

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

FMDenaro · May 21, 2023, 14:16

Quote:

Originally Posted by ALVRM

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?

ALVRM · May 21, 2023, 16:38

Quote:

Originally Posted by 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

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"?

LuckyTran · May 21, 2023, 17:02

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.

ALVRM · May 21, 2023, 17:04

Quote:

Originally Posted by FMDenaro

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.

FMDenaro · May 21, 2023, 17:24

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?

ALVRM · May 21, 2023, 17:26

Quote:

Originally Posted by LuckyTran

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 · May 21, 2023, 17:31

Quote:

Originally Posted by FMDenaro

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.

sbaffini · May 22, 2023, 07:12

Quote:

Originally Posted by ALVRM

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

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.

sbaffini · May 22, 2023, 08:28

Quote:

Originally Posted by ALVRM

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.

FMDenaro · May 22, 2023, 09:04

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.

gnwt4a · May 23, 2023, 02:35

Quote:

Originally Posted by ALVRM

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.

FMDenaro · May 23, 2023, 04:22

Quote:

Originally Posted by gnwt4a

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.

jpanchog · May 23, 2023, 08:24

Quote:

Originally Posted by sbaffini

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

sbaffini · May 23, 2023, 09:47

Quote:

Originally Posted by jpanchog

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.

FMDenaro · May 23, 2023, 11:48

Quote:

Originally Posted by jpanchog

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?

jpanchog · May 23, 2023, 13:32

Quote:

Originally Posted by FMDenaro

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?

FMDenaro · May 23, 2023, 13:47

Quote:

Originally Posted by jpanchog

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.

ALVRM · May 23, 2023, 15:32

Quote:

Originally Posted by gnwt4a

*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.

May 20, 2023, 23:53	Gravity source term causes instability?	#1
ALVRM New Member Join Date: May 2023 Posts: 15 Rep Power: 3	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.

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 19:13
[Other] How to use finite area method in official OpenFOAM 2.2.0?	Detian Liu	OpenFOAM Meshing & Mesh Conversion	4	November 3, 2015 04:04
[swak4Foam] Swak4FOAM 0.2.3 / OF2.2.x installation error	FerdiFuchs	OpenFOAM Community Contributions	27	April 16, 2014 16:14
[swak4Foam] Error bulding swak4Foam	sfigato	OpenFOAM Community Contributions	18	August 22, 2013 13:41
[swak4Foam] build problem swak4Foam OF 2.2.0	mcathela	OpenFOAM Community Contributions	14	April 23, 2013 14:59

May 21, 2023, 11:31		#2
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,195 Blog Entries: 29 Rep Power: 39	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

May 21, 2023, 17:02		#5
LuckyTran Senior Member Lucky Join Date: Apr 2011 Location: Orlando, FL USA Posts: 5,761 Rep Power: 66	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.

May 21, 2023, 17:24		#7
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,896 Rep Power: 73	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?

May 22, 2023, 09:04		#12
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,896 Rep Power: 73	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.