Understanding Terms in Compressible pEqn

AMK53 · November 29, 2017, 09:47

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?

AMK53 · November 29, 2017, 12:54

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.

AMK53 · December 2, 2017, 13:28

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.

Ahmedshoeb · February 19, 2018, 12:17

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 ?

janGaertner · October 16, 2018, 10:59

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)

kaifu · September 18, 2019, 02:50

Quote:

Originally Posted by AMK53

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

November 29, 2017, 09:47	Understanding Terms in Compressible pEqn	#1
AMK53 New Member Join Date: Oct 2017 Posts: 12 Rep Power: 9	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) + rhorAUffvc::ddtCorr(rho, U, phi) + rhorAUffvc::interpolate(rho)(g & mesh.Sf()) ); Code: // Setup pressure correction equation fvScalarMatrix pDDtEqn ( fvc::ddt(rho) + psicorrection(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 psicorrection(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? kaifu, massive_turbulence, lukasf and 2 others like this.

November 29, 2017, 12:54		#2
AMK53 New Member Join Date: Oct 2017 Posts: 12 Rep Power: 9	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 "psifvm::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. Kummi likes this.

December 2, 2017, 13:28		#3
AMK53 New Member Join Date: Oct 2017 Posts: 12 Rep Power: 9	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] = psicorrection(fvm::ddt(P))$ Adding all the terms, we get: Code: fvc::ddt(rho) + psicorrection(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 psicorrection(fvm::ddt(p)) and I just don't know? Any help is appreciated. mickbatti, bharat_aero, peyman.havaej and 4 others like this.

February 19, 2018, 12:17	Have you been able to solve for the missing term ?	#4
Ahmedshoeb New Member Ahmed Shoeb Join Date: Sep 2017 Location: Germany Posts: 3 Rep Power: 9	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 ?

October 16, 2018, 10:59		#5
janGaertner New Member Jan Gaertner Join Date: Nov 2017 Posts: 20 Rep Power: 9	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) Zhiheng Wang and lukasf like this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
How rho_ and psi_ are calculated in compressible solvers of OpenFOAM?	JBUNSW	OpenFOAM Programming & Development	17	February 12, 2021 09:59
Flux discretization of 2D Compressible Navier-Stokes	selig5576	Main CFD Forum	7	March 29, 2017 23:25
Simulation of high pressure diesel injector - all phases compressible with cavitation	fivos	CFX	4	July 30, 2015 07:48
One phase is compressible, can I use compressibleinterfoam?	sharonyue	OpenFOAM Running, Solving & CFD	2	March 31, 2014 10:40
momentum under-relaxation for compressible flow with SIMPLE	Mihai ARGHIR	Main CFD Forum	0	April 7, 2000 05:58