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

Why we subtract continuity Eqn into momentum Eqn instead of adding it?

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

Like Tree1Likes
  • 1 Post By Cyp

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 27, 2015, 04:14
Default Why we subtract continuity Eqn into momentum Eqn instead of adding it?
  #1
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17
sharonyue is on a distinguished road
Hey guys,

This piece of code is quite common in OpenFOAM;
Code:
 fvm::ddt(alpha1, U1)  //mom eqn
         -fvm::Sp(fvc::ddt(alpha1),U1) //continuity eqn
          + fvm::div(alpha1, U1) //mom eqn
          - fvm::Sp(fvc::div(alpha1), U1) ////continuity eqn
Why dont we add continuity eqn?...like this:..

Code:
 fvm::ddt(alpha1, U1)  //mom eqn
         +fvm::Sp(fvc::ddt(alpha1),U1) //continuity eqn
          + fvm::div(alpha1, U1) //mom eqn
          +fvm::Sp(fvc::div(alpha1), U1) ////continuity eqn
continuity is zero right?
sharonyue is offline   Reply With Quote

Old   April 27, 2015, 06:38
Default
  #2
Member
 
ali alkebsi
Join Date: Jan 2012
Location: Strasbourg, France
Posts: 82
Rep Power: 14
kebsiali is on a distinguished road
I dont really know
normally we seek diagnal dominance, so i would say so
but thats the contrary of what is done here ??????
kebsiali is offline   Reply With Quote

Old   April 27, 2015, 12:31
Default
  #3
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 838
Rep Power: 17
sharonyue is on a distinguished road
Yep...Wish someone who can explain this. lol
sharonyue is offline   Reply With Quote

Old   April 27, 2015, 13:32
Default
  #4
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 18
Cyp is on a distinguished road
Quote:
Originally Posted by sharonyue View Post
Yep...Wish someone who can explain this. lol
Hi,


Actually, we subtract a momentum transfer term because the real equation (left hand side of NS) we want to solve is

\rho \frac{\partial \mathsf{U}}{\partial t} + \rho  \mathsf{U} . \nabla  \mathsf{U} = (...)

which is written in a non-conservative form. In FVM, however, we prefer conservative form,
\frac{\partial \rho \mathsf{U}}{\partial t} + \nabla.  \left( \rho  \mathsf{U}   \mathsf{U} \right)= (...)

These two equations are not the same if the fluid is not incompressible or if there is phase change. You can notice that you can switch from non-conservative to conservative form with,
\rho \frac{\partial \mathsf{U}}{\partial t} + \rho  \mathsf{U} . \nabla  \mathsf{U}  = \frac{\partial \rho \mathsf{U}}{\partial t} + \nabla.  \left( \rho  \mathsf{U}   \mathsf{U}  \right) -  \left(  \frac{\partial \rho }{\partial t} + \nabla.  \left( \rho  \mathsf{U}  \right) \right)\mathsf{U}

So adding the continuity equation do not have any sense.

Cheers,
wyldckat likes this.
Cyp is offline   Reply With Quote

Old   April 28, 2015, 03:21
Default
  #5
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 Cyp View Post
Hi,


Actually, we subtract a momentum transfer term because the real equation (left hand side of NS) we want to solve is

\rho \frac{\partial \mathsf{U}}{\partial t} + \rho  \mathsf{U} . \nabla  \mathsf{U} = (...)

which is written in a non-conservative form. In FVM, however, we prefer conservative form,
\frac{\partial \rho \mathsf{U}}{\partial t} + \nabla.  \left( \rho  \mathsf{U}   \mathsf{U} \right)= (...)

These two equations are not the same if the fluid is not incompressible or if there is phase change. You can notice that you can switch from non-conservative to conservative form with,
\rho \frac{\partial \mathsf{U}}{\partial t} + \rho  \mathsf{U} . \nabla  \mathsf{U}  = \frac{\partial \rho \mathsf{U}}{\partial t} + \nabla.  \left( \rho  \mathsf{U}   \mathsf{U}  \right) -  \left(  \frac{\partial \rho }{\partial t} + \nabla.  \left( \rho  \mathsf{U}  \right) \right)\mathsf{U}

So adding the continuity equation do not have any sense.

Cheers,
hey Cyprien,

Yeah, The problem is why we want to solve the first Eqn in your posts. The reason why I ask this kind of problem is that: "openfoam22x is solving non-conservative equations but 23x is solving conservative eqns."for two-fluid model.

So I think the thing deeper is :

1. Why we solve Eqn 2 instead of Eqn 1. I saw some posts that Hrv and Henry explain this, They just said the conservative eqn is stable and good. I dont know the reason.

2. Eqn 2 and Eqn 1 should predicts the same results rite? But sometimes it does not.


Best,
sharonyue is offline   Reply With Quote

Old   April 28, 2015, 12:39
Default
  #6
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 18
Cyp is on a distinguished road
Hi Dongyue,

Eq 1 and Eq 2 are equivalent only if the continuity equation (equal to zero) is satisfy. I guess that since we deal with an operator-splitting scheme (first we solve the momentum, then the pressure, then the correction, then we iterate), the continuity equation is not fully respected, and that why we add to add this term, to avoid spurious momentum transfer due to an ill-evaluation of the continuity.


Cheers,
Cyp is offline   Reply With Quote

Old   June 3, 2015, 13:00
Default
  #7
Member
 
Bruno Blais
Join Date: Sep 2013
Location: Canada
Posts: 64
Rep Power: 12
blais.bruno is on a distinguished road
Quote:
Originally Posted by sharonyue View Post
hey Cyprien,

Yeah, The problem is why we want to solve the first Eqn in your posts. The reason why I ask this kind of problem is that: "openfoam22x is solving non-conservative equations but 23x is solving conservative eqns."for two-fluid model.

So I think the thing deeper is :

1. Why we solve Eqn 2 instead of Eqn 1. I saw some posts that Hrv and Henry explain this, They just said the conservative eqn is stable and good. I dont know the reason.

2. Eqn 2 and Eqn 1 should predicts the same results rite? But sometimes it does not.


Best,

Equation 1 and 2 are not equivalent numerically because they do not allow for the same Riemann invariant. If you take the example of a shockwave or any other type of discontinuity, as a Cauchy-Riemann classical problem, you find that the wave propagation / shock propagation/ discontinuity propagation will be at the wrong velocity if you do not formulate your problem using the conserved variables.

I am not familiar with these issues in the context of two phases flows, but for compressible flows this is an extremely critical issue and you always need to formulate the problem in conservative variables (rho, rho * u, rho * e). I believe that since these multiphase flows contain large amount of discontinuity in the phases (like air on top of water or whatever), a conservative formulation is a lot more appropriate if you want to obtain the right propagation of the interface.

There is a very very nice book by Euleterio Toro on hyperbolic systems that discusses these issues if I remember.
blais.bruno is offline   Reply With Quote

Old   June 4, 2015, 06:20
Default
  #8
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 14
danny123 is on a distinguished road
Hi,

I am confused with these statements. Sharonyue refers to the alpha1Eq, not the momentum equation. In interDymFoam, I have the following code:

Code:
        fvScalarMatrix alpha1Eqn
        (
            #ifdef LTSSOLVE
            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
            #else
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            #endif
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
        );
which is the same, no?

I changed to:

Code:
        fvScalarMatrix alpha1Eqn
        (
            #ifdef LTSSOLVE
            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
            #else
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            #endif
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
          - fvc::div(phi)*fvm::Sp(1,alpha1)
        );
It does not make a whole lot of difference in my case, but the idea is that there are large areas, where alpha1 does not change (only around the interface). In those domains I want interFoam to behave like a single phase solver, so any incosistency in mass conservation should not affect alpha1Eq (d alpha1/dt is 0 if there is no alpha change in space). If alpha1 is constant, the last 2 terms become equivalent. Because div(phi) has to be 0, the ultimate result should be identical. You basically eliminate solving alpha1Eq, where it is not needed (the equation system is overdetermined).

By the way, for compressible solvers, this code is wrong (div phi is not zero, but div rho phi.

Regards,

Daniel
danny123 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
Micro Scale Pore, icoFoam gooya_kabir OpenFOAM Running, Solving & CFD 2 November 2, 2013 14:58
How to write k and epsilon before the abnormal end xiuying OpenFOAM Running, Solving & CFD 8 August 27, 2013 16:33
Upgraded from Karmic Koala 9.10 to Lucid Lynx10.04.3 bookie56 OpenFOAM Installation 8 August 13, 2011 05:03
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 03:58
Could anybody help me see this error and give help liugx212 OpenFOAM Running, Solving & CFD 3 January 4, 2006 19:07


All times are GMT -4. The time now is 10:00.