Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab.pdf), страница 101

Thebound.H class is used to bound variables within certain limits.#include#include#include#include#include#include#include#include#include#include#include// * * ** * * //"fvCFD.H""psiThermo.H""RASModel.H""upwind.H""gaussConvectionScheme.H""bound.H""simpleControl.H""totalPressureCompFvPatchScalarField.H""totalPressureCorrectorCompFvPatchScalarField.H""totalVelocityFvPatchVectorField.H""orthogonalSnGrad.H"* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *int main(int argc, char *argv[]){###include "setRootCase.H"include "createTime.H"include "createMesh.H"simpleControl simple(mesh);#include "createFields.H"Listing 16.3 The include files used in simpleFoamCompressibleThe createFields.H now includes the definition of the density field and othervariables and constants related to compressible flow physics.

The psiTermo class,depicted in Listing 16.4, provides access to the thermophysical relations that arepart of the OpenFOAM® library [27], such as the perfect gas law described inEq. (16.4).67616Fluid Flow Computation: Compressible FlowsInfo<< "Reading thermophysical properties\n" << endl;autoPtr<psiThermo> thermo(psiThermo::New(mesh));Info<< "Calculating field rho\n" << endl;volScalarField rho(IOobject("rho",runTime.timeName(),mesh,IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),thermo->rho());volScalarField& p = thermo->p();volScalarField& h = thermo->he();thermo->correct();volScalarField gammaGas(IOobject("gammaGas",runTime.timeName(),mesh,IOobject::NO_READ,IOobject::AUTO_WRITE),thermo->Cp() / thermo->Cv());volScalarField RGas(IOobject("RGas",runTime.timeName(),mesh,IOobject::NO_READ,IOobject::AUTO_WRITE),thermo->Cp() - thermo->Cv());Listing 16.4 An excerpt of the psiTermo class16.9Computational Pointers677The pressure and enthalpy fields are defined in the thermo class, displayed inListing 16.5, and accessed as references in createFields.H.volScalarField& p = thermo->p();volScalarField& h = thermo->he();const volScalarField::GeometricBoundaryField& pbf=p.boundaryField();wordList pbt = pbf.types();volScalarField pp(IOobject("pp",runTime.timeName(),mesh,IOobject::NO_READ,IOobject::AUTO_WRITE),mesh,dimensionedScalar("zero", p.dimensions(), 0.0),pbt);Listing 16.5 Defining the pressure and enthalpy field in thermo classThe velocity is defined as in the incompressible version but the mass flux nowincludes density in its definition (Listing 16.6).Info << "Reading field U\n" << endl;volVectorField U(IOobject("U",runTime.timeName(),mesh,IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);surfaceScalarField mDot(IOobject("mDot",runTime.timeName(),mesh,IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),linearInterpolate(rho*U) & mesh.Sf());Listing 16.6 Defining the velocity and mass flux fields67816Fluid Flow Computation: Compressible FlowsIn compressible flows, setting physical limits on some variables such as densityand pressure can enhance robustness, especially during the first few iterations.

Thisprevents variables from assuming non-physical values (like negative densities orpressures). Therefore bounds can be set as part of the case definition, as shown inListing 16.7, and read in createFields.H.dimensionedScalar rhoMax(simple.dict().lookup("rhoMax"));dimensionedScalar rhoMin(simple.dict().lookup("rhoMin"));dimensionedScalar pMax(simple.dict().lookup("pMax"));dimensionedScalar pMin(simple.dict().lookup("pMin"));Listing 16.7 Setting lower and upper bounds for the density and pressure fieldsThe momentum equation is defined with a slightly modified syntax that accountsfor density and thermophysical property relations. The syntax of the linearizedformula is given in Listing 16.8.// Solve the Momentum equationfvVectorMatrix UEqn(fvm::ddt(rho,U)+ fvm::div(mDot, U)- fvm::laplacian(thermo->mu, U));UEqn.relax();solve(UEqn == -fvc::grad(p));Listing 16.8 Syntax used to solve the momentum equationThe first instruction defines the finite volume discretization of the momentumequation in a vector form (the three components of the velocity are solved in asegregated manner despite the vectorial implementation).

The system is implicitlyrelaxed and then solved with an iterative solver.Once the momentum equation is solved, a new guess of the velocity field isobtained. This velocity field does not necessarily satisfy the continuity equation.16.9Computational Pointers679To enforce mass conservation, assembly of the pressure correction equation is nowrequired to correct the velocity. Following Eq. (16.19) the syntax used for thatpurpose is shown in Listing 16.9.pp = scalar(0.0)*pp;pp.correctBoundaryConditions();surfaceScalarField phid(IOobject("phid",runTime.timeName(),mesh,IOobject::NO_READ,IOobject::NO_WRITE),mDot*drhodp/(rhofd));fvMatrix<scalar> ppCompEqn(fvm::ddt(thermo->psi(),pp) +fv::gaussConvectionScheme<scalar>(mesh,phid,tmp<surfaceInterpolationScheme<scalar> >(new upwind<scalar>(mesh,phid))).fvmDiv(phid, pp)- fvm::laplacian(pDiff, pp)+ fvc::div(mDot) + fvc::ddt(rho))Listing 16.9 Syntax used to assemble the pressure correction equationAs for the incompressible case, in order to avoid checker boarding, the mDotmass flux field is evaluated using the Rhie-Chow interpolation but now taking intoaccount also the density field, evaluated based on the thermo model as shown inListing 16.10.68016Fluid Flow Computation: Compressible Flowsrho = thermo->rho();Foam::fv::orthogonalSnGrad<scalar> faceGradient(mesh);surfaceVectorField gradp_avg_f = linearInterpolate(fvc::grad(p));surfaceVectorField gradp_f = gradp_avg_f - (gradp_avg_f & ed)*ed +(faceGradient.snGrad(p))*ed;surfaceVectorField U_avg_prevIter_f = linearInterpolate(U.prevIter());surfaceVectorField U_avg_f = linearInterpolate(U);surfaceScalarField rhofd = upwind<scalar>(mesh,mDot).interpolate(rho);surfaceScalarField rhof("srho",fvc::interpolate(rho));surfaceScalarField DUf("srUA",fvc::interpolate(DU,"interpolate((1|A(U)))"));volScalarField dt(IOobject("dt",runTime.timeName(),mesh,IOobject::NO_READ,IOobject::NO_WRITE),mesh ,dimensionedScalar("dt" ,runTime.deltaT().dimensions(),runTime.deltaT().value()),zeroGradientFvPatchScalarField::typeName);surfaceScalarField dt_f = linearInterpolate(dt);surfaceScalarField drhodp = linearInterpolate(thermo->psi());scalar UURF = mesh.equationRelaxationFactor("U");// Rhie-Chow interplationmDot = rhof*((U_avg_f & mesh.Sf()) - ( (DUf*( gradp_f - gradp_avg_f))& mesh.Sf() ));+ (scalar(1) - UURF)*(mDot.prevIter() ( (rhof*U_avg_prevIter_f) & mesh.Sf() ) )+rhof*(DUf/dt_f)*(mDot.prevIter() ( (rhof*U_avg_prevIter_f) & mesh.Sf() ) );Listing 16.10 Calculation of mass fluxes at cell faces using the Rhie-Chow interpolationIt is worth mentioning that density is interpolated to faces using an upwindscheme in order to mimic the hyperbolic behavior of compressible flows.The pressure correction equation is fully set and is solved using the syntaxdisplayed in Listing 16.11.16.9Computational Pointers681ppEqn.solve();Listing 16.11 Syntax for solving the pressure correction equationAfter solving the pressure correction equation, variables that depend on pressurecorrection are updated.

For the mass flux field this is performed using the syntax inListing 16.12, this is similar to the incompressible flux correction.mDot += ppEqn.flux();Listing 16.12 Syntax to update the mass flux fieldWhere again the flux() function in Listing 16.12 updates the fluxes using directlythe matrix coefficients and cell values. A simplified version of the flux() function isshown in Listing 16.13.for (label face=0; face<lowerAddr.size(); face++){mDotPrime[face] =upperCoeffs[face]*pp[upperAddr[face]]- lowerCoeffs[face]*pp[lowerAddr[face]];}return mDotPrime;Listing 16.13 A simplified version of the flux( ) function where the flux correction mDotPrime iscomputedIn Listing 16.13 the correction flux mDotPrime is basically evaluated by performing a loop over the faces using the upper and lower coefficients of the matrixand multiplying these coefficients with the corresponding cell values.Finally the velocity, density and pressure at cell centroids are updated usingEqs.

(16.20), (16.21), and (16.22), as shown in Listing 16.14, where the variablealphaP is the explicit relaxation factor for pressure and density updates kp , necessary for a stable SIMPLE solver.68216Fluid Flow Computation: Compressible Flowsscalar alphaP = mesh.equationRelaxationFactor("pp");mDot += ppCompEqn.flux();p += alphaP*pp;p.correctBoundaryConditions();rho += alphaP*pp*thermo->psi();boundMinMax(p, pMin, pMax);boundMinMax(rho, rhoMin, rhoMax);U -= fvc::grad(pp)*DU;U.correctBoundaryConditions();Listing 16.14 Update of the velocity and pressure fields at cell centroidsIn order to account for compressibility effects, the energy equation is introduced®and the related temperature is calculated.

In OpenFOAM , the energy equationexpressed in terms of specific static enthalpy h ¼ Cp T , given by Eq. (3.61), issolved as depicted in Listing 16.15.fvScalarMatrix hEqn(fvm::ddt(rho,h)+ fvm::div(mDot, h)- fvm::laplacian(turbulence->alphaEff(), h)+ fvm::SuSp(-fvc::div(mDot),h)==fvc::div(mDot/fvc::interpolate(rho)*fvc::interpolate(p))- p*fvc::div(mDot/fvc::interpolate(rho)));hEqn.relax();hEqn.solve();h.correctBoundaryConditions();thermo->correct();gammaGas = thermo->Cp()/ thermo->Cv();gammaGas.correctBoundaryConditions();RGas = thermo->Cp() - thermo->Cv();RGas.correctBoundaryConditions();Listing 16.15 Solving the energy equation16.9Computational Pointers683Once the energy equation is solved, the new enthalpy is used to update thetemperature and gas properties (e.g., specific heats).In addition to the main solver, new total pressure and total temperatureboundary conditions are implemented for subsonic inlet patches, these are oftenused boundary conditions for the simulation of compressible flows.