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

Understanding Terms in Compressible pEqn

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

Like Tree15Likes
  • 5 Post By AMK53
  • 1 Post By AMK53
  • 7 Post By AMK53
  • 2 Post By janGaertner

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 29, 2017, 09:47
Default Understanding Terms in Compressible pEqn
  #1
New Member
 
Join Date: Oct 2017
Posts: 12
Rep Power: 8
AMK53 is on a distinguished road
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 is offline   Reply With Quote

Old   November 29, 2017, 12:54
Default
  #2
New Member
 
Join Date: Oct 2017
Posts: 12
Rep Power: 8
AMK53 is on a distinguished road
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.
Kummi likes this.
AMK53 is offline   Reply With Quote

Old   December 2, 2017, 13:28
Default
  #3
New Member
 
Join Date: Oct 2017
Posts: 12
Rep Power: 8
AMK53 is on a distinguished road
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.
AMK53 is offline   Reply With Quote

Old   February 19, 2018, 12:17
Default Have you been able to solve for the missing term ?
  #4
New Member
 
Ahmed Shoeb
Join Date: Sep 2017
Location: Germany
Posts: 3
Rep Power: 8
Ahmedshoeb is on a distinguished road
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 ?
Ahmedshoeb is offline   Reply With Quote

Old   October 16, 2018, 10:59
Default
  #5
New Member
 
Jan Gaertner
Join Date: Nov 2017
Posts: 17
Rep Power: 8
janGaertner is on a distinguished road
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.
janGaertner is offline   Reply With Quote

Old   September 18, 2019, 02:50
Default
  #6
Member
 
Kai
Join Date: May 2010
Location: Shanghai
Posts: 61
Blog Entries: 1
Rep Power: 15
kaifu is on a distinguished road
Quote:
Originally Posted by AMK53 View Post
#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.
__________________
Kai
kaifu is offline   Reply With Quote

Reply

Tags
compressible, openfoam, peqn, pressure correction

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


All times are GMT -4. The time now is 04:11.