Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 95
Текст из файла (страница 95)
The basesolver, simpleFoam, will be presented first. This is followed by a number of versions, with each one adding more capabilities to the base code. These solvers can besummarized as follows:1. simpleFoam (not the OpenFOAM® built-in solver) is the base code thatincorporates the SIMPLE Algorithm in its most basic form.2. simpleFoamImproved extends the base code to allow for improved treatment ofrelaxation.3. simpleFoamTransient adds transient capabilities to the steady-state simpleFoam.4.
simpleFoamBuoyancy adds to the code the body force treatment.More versions will be covered in the chapters to follow, each one with extendedcapabilities, added by modifying the base code described in this chapter. A list ofthe versions that will be covered in the next chapters is given below.5. simpleFoamCompressible is the compressible version of simpleFoam (Chap. 16)6. simpleFoamTurbulent includes capabilities for treating turbulent flows(Chap. 17).simpleFoamBefore reviewing the simpleFoam code, some basic notational issues are addressed.The first step is to define, as shown in Listing 15.2, the geometric fields andparameters that will be initialized and used in the code.15.10Computational Pointers639#include "fvCFD.H"#include "orthogonalSnGrad.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** * * //int main(int argc, char *argv[]){###include "setRootCase.H"include "createTime.H"include "createMesh.H"Listing 15.2 The #include macro derivatives used to define the types of objects neededIn Listing 15.2, the #include macro directives outside the main function areneeded to define the types of objects that are then declared and used in the application.
The #include “fvCFD.H” contains a list of definitions for classes that are ingeneral necessary to build any application in OpenFOAM®. In the developedapplication an additional header, not present in the fvCFD.H header, necessary forthe SIMPLE solver implementation will be added.The use of the #include statements inside the main function is a compactingprocedure, with each declared statement representing a piece of the code moved tothe corresponding file name. For example, the statement #include “createMesh.H”just represents the code shown in Listing 15.3, which is necessary to instantiate themesh class.createMesh.H~~~~~~~~~~~Foam::Info<< "Create mesh for time = "<< runTime.timeName() << Foam::nl << Foam::endl;Foam::fvMesh mesh(Foam::IOobject(Foam::fvMesh::defaultRegion,runTime.timeName(),runTime,Foam::IOobject::MUST_READ));Listing 15.3 The code representing the #include createMesh.H file necessary to instantiate themesh class64015 Fluid Flow Computation: Incompressible FlowsOnce the necessary initialization has been performed, the next step is thedefinition of the proper fields or variables needed by the solver.
These are definedin file “createFields.H”. The first defined field, shown in Listing 15.4, is thepressure field (p).Info << "Reading field p\n" << endl;volScalarField p(IOobject("p",runTime.timeName(),mesh,IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);Listing 15.4 Script used to define the pressure fieldSince the solution is obtained by solving a pressure correction equation insteadof a pressure equation, a pressure correction field (pp) is also defined (Listing 15.5).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 15.5 Script used to define the pressure correction fieldIt is worth noting that in Listing 15.5 a different constructor is used to define thepressure correction field.
Since the pressure corrector represents the pressure itself,the same boundary conditions defined for the real pressure can be used for thepressure correction field without the need to define the same quantity twice.The list of pressure boundary types displayed in Listing 15.6 are now copiedunder the pbt variable in Listing 15.5 and directly used in the constructor of pp.15.10Computational Pointers641// Set pp boundary valuesforAll(pp.boundaryField(), patchi){if (isType<fixedValueFvPatchScalarField>(pp.boundaryField()[patchi])){fixedValueFvPatchScalarField& ppbound =refCast<fixedValueFvPatchScalarField>(pp.boundaryField()[patchi]);ppbound == scalarField(ppbound.size(),0.0);}}Listing 15.6 Script showing the declaration of the different pressure boundary typesThe corrector should reset to zero the correction field at every iteration andshould also apply a zero value at all boundaries for which a Dirichlet boundarycondition is used for the pressure.The velocity and corresponding mass flux fields must also to be defined.
Asdepicted in Listing 15.7, the velocity field is defined through an input file while themass flux field can be defined as a derived quantity.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(U) & mesh.Sf());Listing 15.7 Script used to define the velocity and mass flux fields64215 Fluid Flow Computation: Incompressible FlowsFinally the fluid thermo-physical properties should also be defined. For incompressible laminar flows this involves simply assigning a value to the kinematicviscosity nu, as shown in Listing 15.8.Info<< "Reading transportProperties\n" << endl;IOdictionary transportProperties(IOobject("transportProperties",runTime.constant(),mesh,IOobject::MUST_READ_IF_MODIFIED,IOobject::NO_WRITE));dimensionedScalar nu(transportProperties.lookup("nu"));Listing 15.8 Code used to define the fluid thermo-physical propertiesAfter defining all variables the implementation of the SIMPLE algorithm canproceed.
A while loop can be used for the cases when the stopping criterion is thenumber of SIMPLE iterations. For each single loop, the momentum and pressurecorrection equations are solved and updates of the variables are performed. Startingwith the momentum equation written as (the form solved in OpenFOAM®),r fvvg ¼ vr2 v rpð15:217Þwhere the kinematic viscosity v is defined asv¼lqð15:218Þand p represents the pressure divided by the density, i.e.,p¼static pressureqits solution is translated into the script shown in Listing 15.9.ð15:219Þ15.10Computational Pointers643// Solve the Momentum equationfvVectorMatrix UEqn(fvm::div(mDot, U)- fvm::laplacian(nu, U));UEqn.relax();solve(UEqn == -fvc::grad(p));Listing 15.9 Script to solve the momentum equationThe first instruction in Listing 15.9 defines the finite volume discretization of themomentum equation in vector form with its storage matrix (the three components ofthe velocity vector are solved in a segregated manner despite its vectorial implementation).
The system is then implicitly relaxed and solved via an iterative solver.Once the momentum equation is solved, a new guess for the velocity field isobtained. This velocity does not satisfy the continuity equation in general and theassembly of the continuity equation in the form of a pressure correction equation isnow required to correct the flow field. The pressure correction equation is written as f rp0 Þ ¼ r ðvÞr ðDð15:220Þ at an element centroid is computed aswhere a component of DD¼Vavcð15:221Þwith values at the faces obtained by interpolation.
Its implementation is translatedinto the following syntax (Listing 15.10):pp = scalar(0.0)*pp;pp.correctBoundaryConditions();fvScalarMatrix ppEqn(- fvm::laplacian(DUf, pp, "laplacian(pDiff,pp)")+ fvc::div(mDot));Listing 15.10 Script to implement the pressure correction equation64415 Fluid Flow Computation: Incompressible FlowsPwhere div(mDot) is basicallym_ f , while pp represents p′ and is reset to zero ateach iteration. Moreover, the DUf variable is the value of D at the cell faceobtained, as shown in Listing 15.11, by linear interpolation using the values at thenodes straddling the face. Further, .A() represents the diagonal terms in themomentum matrix divided by the volume.volScalarField DU = 1.0/UEqn.A();surfaceScalarField DUf("DUf",linearInterpolate(DU));Listing 15.11 Script to calculate the values of DUfListing 15.12 computes the mass flow rate at cell faces (mDot) using theRhie-Chow interpolation where the calculation of ∇pf and rpf is clearly shown.const surfaceVectorField ed = mesh.delta()()/mag(mesh.delta()());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_f = linearInterpolate(U);// Rhie-Chow interplationmDot = (U_avg_f & mesh.Sf()) - ( (DUf*( gradp_f - gradp_avg_f)) &mesh.Sf() );Listing 15.12 Script to compute the mass flow rate at cell faces using the Rhie-ChowinterpolationThe pressure correction equation is now fully set and is solved by executing thestatement in Listing 15.13.ppEqn.solve();Listing 15.13 Statement used to solve the pressure correction equationOnce the pressure correction equation is solved, the velocity, pressure, and massflow rate fields are updated using the obtained pressure correction field.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
