# RhoSimpleFoam

rhoSimpleFoam

rhoSimpleFoam is a steady-state solver for compressible, turbulent flow, using the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm. In the newer releases it also includes an option to use the SIMPLEC (Semi-Implicit Method for Pressure Linked Equations Consistent) algorithm. It is a pressure based compressible solver. That means that a pressure equation is solved and the density is related to the pressure via an equation of state.

## 1 Solution Strategy

The solver follows a segregated solution strategy. This means that the equations for each variable characterizing the system (the velocity $\bold u$, the pressure $p$, the energy (either internal energy or enthalpy depending on the choice of the user) and the variables characterizing turbulence) is solved sequentially and the solution of the preceding equations is inserted in the subsequent equation. The non-linearity appearing in the momentum equation (the face flux $\phi_f$ which is a function of the velocity and density $\rho_f$ ) is resolved by computing it from the velocity and pressure values of the preceding iteration. Note that the density is related to the pressure via an equation of state which can be selected by the user. The dependence from the pressure is introduced to avoid a decoupling between the momentum and pressure equations and hence the appearance of high frequency oscillation in the solution (check board effect). The first equation to be solve is the momentum equation. It delivers a velocity field $\bold u^*$ which is in general not divergence free, i.e. it does not satisfy the continuity equation. After that the energy equation is solved. Then that the momentum and the continuity equations are used to construct an equation for the pressure. The aim is to obtain a pressure field $p^{n}$, which, if inserted in the momentum equation, delivers a divergence free velocity field $\bold u$. After correcting the velocity field, the equations for turbulence are solved. The above iterative solution procedure is repeated until convergence.

A good book which describes the derivation of the pressure equation for compressible cases is the one of [1].

The source code can be found in rhoSimpleFoam.C



int main(int argc, char *argv[])
{
(
"Steady-state solver for compressible turbulent flow."
);

#include "postProcess.H"

#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "initContinuityErrs.H"

turbulence->validate();

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

// Pressure-velocity SIMPLE corrector
#include "UEqn.H"
#include "EEqn.H"

if (simple.consistent())
{
#include "pcEqn.H"
}
else
{
#include "pEqn.H"
}

turbulence->correct();

runTime.write();

runTime.printExecutionTime(Info);
}

Info<< "End\n" << endl;

return 0;
}


## 2 Equations

In the next section the equation to be solved are described. In contrast to incompressible fluids, where the density is constant, in compressible flows the density varies in space and time, i.e., it the density became a thermodynamic property. For a pure mixture the thermodynamic state is fully described by to variables, i.e., $\rho = \rho(p,T)$ the density is a function of the pressure and the temperature. For this reason in compressible flows we need to solve besides the equations for the momentum and the pressure (which is derived from the continuity equations) also an equation which gives us the evolution of the temperature. We will see later that the temperature is obtained via a thermodynamic relation after we solved the energy equation. After that, the density $\rho$ can be obtained from the other two variables via an equation of state.

### 2.1 Momentum Equation

The equation of motion in rhosimpleFoam are written for a moving frame of reference. They are however formulated for the absolute velocity (the derivation of the equations of motion can be found in https://openfoamwiki.net/index.php/See_the_MRF_development and also in https://diglib.tugraz.at/download.php?id=581303c7c91f9&location=browse. Some additional information can be found in https://pingpong.chalmers.se/public/pp/public_courses/course07056/published/1497955220499/resourceId/3711490/content/UploadedResources/HakanNilssonRotatingMachineryTrainingOFW11-1.pdf):

 - $\nabla \cdot ( \rho \bold{u_R} \otimes \bold{u_I}) + \rho \bold{\Omega} \times \bold{u_I} = \nabla p + \nabla \cdot (\bold{R_{vis}} + \bold{R_{tur}} )$ (1)

$\bold{u_R}$ represents the velocity in the relative reference frame, $\bold{u_I}$ is the velocity described in the inertial reference frame, $\bold{\Omega}$ is the angular velocity of the moving reference frame, $\bold{R_{vis}} + \bold{R_{tur}}$ is the sum of the viscous and turbulent stress tensor. It is evident that the above equation can be interpreted as the Favre averaged Navier Stokes equation rather than the Reynolds averaged Navier stokes equation. Favre averaging is a transformation which avoids the appearance of correlation involving the density fluctuations in the continuity and momentum equation. For a more detailed discussion see the thread https://www.cfd-online.com/Forums/openfoam-solving/193246-favre-reynolds-average-openfoam.html or the book of [2].

The source code can be found in the file UEqn.H:


// Solve the Momentum equation

MRF.correctBoundaryVelocity(U);

tmp<fvVectorMatrix> tUEqn
(
fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

fvOptions.constrain(UEqn);

fvOptions.correct(U);

Note that the above equation differs from the equations without considering a moving frame of reference in two terms: $\rho \bold{u_R}$ which if interpolated on the face of the control volume gives the relative face flux and the term $\rho \bold{\Omega} \times \bold{u_I}$ which is simply a explicit source term in the momentum equation. We will see later in the pressure equation how the relative face flux is calculated.

The source code of the acceleration resulting from the description in a moving frame of reference can be found in the following src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C


Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::DDt
(
const volVectorField& U
) const
{
tmp<volVectorField> tacceleration
(
new volVectorField
(
IOobject
(
"MRFZoneList:acceleration",
U.mesh().time().timeName(),
U.mesh()
),
U.mesh(),
dimensionedVector(U.dimensions()/dimTime, Zero)
)
);
volVectorField& acceleration = tacceleration.ref();

forAll(*this, i)
{
}

return tacceleration;
}

Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::DDt
(
const volScalarField& rho,
const volVectorField& U
) const
{
return rho*DDt(U);
}

The calculation of the Coriolis force is done in the file src/finiteVolume/cfdTools/general/MRF/MRFZone.C



(
const volVectorField& U,
volVectorField& ddtU
) const
{
if (cellZoneID_ == -1)
{
return;
}

const labelList& cells = mesh_.cellZones()[cellZoneID_];
vectorField& ddtUc = ddtU.primitiveFieldRef();
const vectorField& Uc = U;

const vector Omega = this->Omega();

forAll(cells, i)
{
label celli = cells[i];
ddtUc[celli] += (Omega ^ Uc[celli]);
}
}

The way the turbulent stresses enter into the momentum equation is well described here: https://www.openfoam.com/documentation/guides/latest/doc/guide-turbulence-ras.html.

### 2.2 Energy Equation

The derivation of the Favre agraged energy equation for compressible flows can be found in the book of [3]. A description of the terms resulting in the energy equation from the favre averaging can also be found in https://www.cfd-online.com/Wiki/Favre_averaged_Navier-Stokes_equations.

$\frac{\partial}{\partial t} (\overline{\varrho} \tilde{e_0})+\nabla\cdot\left(\overline{\varrho}\tilde{e_0}\widetilde{\mathbf{u}}+\tilde{\mathbf{u}}\overline{p} + \overline{\varrho \mathbf{u''}e_0''} + \mathbf{q_{vis}} + \overline{\mathbf{\tau \cdot u}} \right) = 0$

The total energy $\tilde{e_0}$ is defined as follows:

$\tilde{e_0} = \tilde{e} + 0.5 \tilde{\mathbf{u}}\cdot \tilde{\mathbf{u}} + \tilde{k}$. The overbar denotes the Reynolds averaging, the tilde the Favre averaging and the two $''$ denote the fluctuation due to the Favre averaging:

$\phi = \tilde{\phi} + \phi''$.

$\tau$ is the viscous stress tensor.

The equation solved can be found in the file EEqn.H


{
volScalarField& he = thermo.he();

fvScalarMatrix EEqn
(
fvm::div(phi, he)
+ (
he.name() == "e"
? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
: fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
)
- fvm::laplacian(turbulence->alphaEff(), he)
==
fvOptions(rho, he)
);

EEqn.relax();

fvOptions.constrain(EEqn);

EEqn.solve();

fvOptions.correct(he);

thermo.correct();
}

# 3 References

1. Moukalled, F., L. Mangani, and M. Darwish. "The finite volume method in computational fluid dynamics." An Advanced Introduction with OpenFOAM and Matlab (2016):
2. Wilcox, David C. Turbulence modeling for CFD. Vol. 3. La Canada, CA: DCW industries, 2006.
3. Wilcox, David C. Turbulence modeling for CFD. Vol. 3. La Canada, CA: DCW industries, 2006.