CFD Online Discussion Forums - Understanding Terms in Compressible pEqn

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - Understanding Terms in Compressible pEqn (https://www.cfd-online.com/Forums/openfoam-solving/196243-understanding-terms-compressible-peqn.html)

Understanding Terms in Compressible pEqn

I am trying to understand some of the terms in the compressible version of the pEqn (rhoPimpleFoam) and had some questions that I have not been able to figure out on my own.

My first question has to do with the pressure equation itself. In "Computational Methods for Fluid Dynamics" chapter 10, the pressure-velocity correction is explained as inputting the predicted velocity and density into the continuity equation with a small correction parameter ( $\dot{m}'$ , $\rho'$ ):

$\frac{\rho' \Delta \Omega}{\Delta t} + \sum (\dot{m}* + \dot{m}') = 0$

Where (in OpenFOAM notation)

$\dot{m}* = \rho (A^{-1} H - A^{-1}\nabla P + A^{-1}\rho g) S$

and according to Ferziger

$\dot{m}' = -(\rho^{m-1} S \Delta\Omega) A^{-1} \nabla{p'}_{n} + \frac{Cp\dot{m}*}{\rho^{m-1}}P'$

In OpenFOAM programming, this is implemented as:

Code:

// Calculate flux of HbyA and add contribution from gravity term

surfaceScalarField phiHbyA

(

    "phiHbyA",

    fvc::interpolate(rho)*fvc::flux(HbyA)

  + rhorAUf*fvc::ddtCorr(rho, U, phi)

  + rhorAUf*fvc::interpolate(rho)*(g & mesh.Sf())

);

Code:

// Setup pressure correction equation

fvScalarMatrix pDDtEqn

(

    fvc::ddt(rho) + psi*correction(fvm::ddt(p))

  + fvc::div(phiHbyA)

  ==

    fvOptions(psi, p, rho.name())

);



while (pimple.correctNonOrthogonal())

{

  fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));



  pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

  

  // Correct velocity flux by (subtracting?) flux from P Equation

  if (pimple.finalNonOrthogonalIter())

  {

     phi = phiHbyA + pEqn.flux();

  }

}

What I don't understand is the following:

1. Why is "fvc::div(phiHbyA) - fvm::laplacian(rhorAUf, P)" is included instead of "phiHbyA - fvm::grad(rhorAUf,P)"? Is there a particular reason the divergence of those terms taken? "(phiHbyA - fvm::grad(rhoAUf,P))*S" should be the correct velocity flux, so taking the divergence seems unnecessary to me. And even if we take the divergence, shouldn't we have to take the divergence of every term in that equation including the 'ddt' terms?

2. Is psi*correction(fvm::ddt(p)) the equivalent $\dot{m}'$ term? Or is this simply the density ( $\rho'$ ) correction?

My second question has to do with how the density is calculated. In rhoPimpleFoam, the density is initially predicted by "rhoEqn.H" prior to calculating the velocity terms:

Code:

if (pimple.nCorrPIMPLE() <= 1)

        {

            #include "rhoEqn.H"

        }



        // --- Pressure-velocity PIMPLE corrector loop

        while (pimple.loop())

        {

            #include "UEqn.H"

            #include "EEqn.H"



            // --- Pressure corrector loop

            while (pimple.correct())

            {

                if (pimple.consistent())

                {

                    #include "pcEqn.H"

                }

                else

                {

                    #include "pEqn.H"

                }

            }



            if (pimple.turbCorr())

            {

                turbulence->correct();

            }

        }

But the code then solves again the "rhoEqn.H" after the PHI term is updated:

Code:

while (pimple.correctNonOrthogonal())

{

  fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));



  pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

  

  // Correct velocity flux by (subtracting?) flux from P Equation

  if (pimple.finalNonOrthogonalIter())

  {

     phi = phiHbyA + pEqn.flux();

  }

}



#include "rhoEqn.H"

#include "compressibleContinuityErrs.H"

And then it is corrected again afterwards by:

Code:

thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;

rho = thermo.rho();

Is there a reason why the density is corrected so many times? And shouldn't the equation state act as the constraint on the density instead of solving the continuity equation again?

In regard to my first question, you can always write the continuity equation in differential form instead of integral form:

$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho U) = 0$

If we break it into the (.)* and (.)' terms, we get for the time derivative:

$\frac{\partial \rho}{\partial t} = \frac{\partial \rho^*}{\partial t} + \frac{\partial \rho'}{\partial t}$

and for the divergence term

$\nabla \cdot (\rho U) = \nabla \cdot (\rho^*U*) + \nabla \cdot (\rho^*U' + U^*\rho' + U'\rho') = 0$

The first term in the time derivative should be "fvc::ddt(rho)" while the second term would then be approximated by "psi*fvm::ddt(p)", where p here is the correction pressure.

For the divergence terms, the first term should be "fvc::div(phiHbyA)-fvc::div(laplacian(rhorAUf, p))", but in rhoPimpleFoam, it is written as "fvc::div(phiHbyA)-fvm::div(laplacian(rhorAUf, p))".

Why is that? And what about the remaining terms in the divergence term (i.e. all the u' and rho' terms)? I think I am missing something, as I am not sure how the pressure equation is derived in rhoPimpleFoam.

I think I have (semi) figured out my first question. If we re-write the continuity equation as:

$\frac{\partial \rho^{*}}{\partial t} + \nabla \cdot (\rho^{*} u^{*}) + \frac{\partial \rho^{'}}{\partial t} + \nabla \cdot (\rho^{*} u^{'}) + + \nabla \cdot (\rho^{'} u^{*}) = 0$

Where (*) is the predicted values and (') are the corrected values. We can then break each term down to:

$\frac{\partial \rho^{*}}{\partial t} = fvc::ddt(rho)$

$\nabla \cdot (\rho^{*} u^{*}) + \nabla \cdot (\rho^{*} u^{'}) = \nabla \cdot \big[\rho^{*}A^{-1}(H-\nabla P^{*} + \rho^{*}g)\big] - \nabla \cdot \big[\rho^{*} A^{-1} \nabla P^{'}\big]$

But from the definition of $P_{new} = P^{*} + P^{'}$ , we can rewrite the above equation as:

$\nabla \cdot \big[\rho^{*}A^{-1}(H-\nabla P^{*} + \rho^{*}g)\big] - \nabla \cdot \big[\rho^{*} A^{-1} \nabla P^{'}\big] = \nabla \cdot\big[\rho^{*} A^{-1}(H+\rho^{*}g)\big] - \nabla \cdot \big[\rho^{*} A^{-1} \nabla P\big]$

This simply equals (in OpenFOAM notation):

Code:

fvc::div(phiHbyA) - fvm::laplacian(rhorAUf,P)

And for the density correction:

$\frac{\partial \rho^{'}}{\partial t} = psi \frac{\partial P^{'}}{\partial t} = psi\big[\frac{\partial P}{\partial t}-\frac{\partial P^{*}}{\partial t}\big] = psi*correction(fvm::ddt(P))$

Adding all the terms, we get:

Code:

fvc::ddt(rho) + psi*correction(fvm::ddt(p))

  + fvc::div(phiHbyA) - fvm::laplacian(rhorAUf, p)

Which is exactly what's given in rhoPimpleFoam solver. The only remaining term that is left that I can't figure out is the following:

$\nabla \cdot (\rho^{'} u^{*}) = \nabla \cdot (psi P^{'} u^{*})$

There doesn't seem to be any additional terms in the OpenFOAM pressure equation to account for this. Is this somehow hidden inside psi*correction(fvm::ddt(p)) and I just don't know?

Any help is appreciated.

Have you been able to solve for the missing term ?

Hello

thank you for your detailed explanation, have you been able to find a solution to account for the missing term ?

Also could you please explain the role of the the fvc::ddtCorr(rho, U, phi) term ?

Hello,

the term $\nabla \cdot ( \psi P' u^*)$ is found when you use the transonic option.

As I understand it they do not use the correction term here and instead the fvm description as there is no

Code:

fvc::div(phi,rho)

Quote:

Originally Posted by AMK53 (Post 673314)

#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
[/CODE]And then it is corrected again afterwards by:

Code:

thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;

rho = thermo.rho();

Is there a reason why the density is corrected so many times? And shouldn't the equation state act as the constraint on the density instead of solving the continuity equation again?

When it goes to convergece, both of density update by thermo model and rhoEqn.H are all the same. I guess during calculaltion, this method could bring robustness to the solver.