Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 35
Текст из файла (страница 35)
The scriptused is shown in Listing 7.47 in which one of the several constructors of the classGeometricField is used (other constructors can be found in the header file of thespecific class GeometricField.H).volScalarField T(IOobject("T",runTime.timeName(),mesh,IOobject::NO_READ,IOobject::AUTO_WRITE),mesh,dimensionedScalar("DTVol", dimensionSet(0,0,0,1,0,0,0), 300.0),"zeroGradient");volVectorField U(IOobject("U",runTime.timeName(),mesh,IOobject::NO_READ,IOobject::AUTO_WRITE),mesh,dimensionedVector("U",dimensionSet(0,1,-1,0,0,0,0),vector::zero),"zeroGradient");Listing 7.47 Script used to define the two variables U and T in OpenFOAM®As can be seen in the code, both fields are linked to the mesh.
Both requirespecifying the following four arguments: (i) field name, (ii) field dimension,(iii) initialization, and (iv) boundary conditions. Specific boundary conditions haveto be defined for the field and in this case a zero order extrapolation (i.e.,“zeroGradient”) is specified for the full set of boundaries and for both variables.7.2 OpenFOAM®199Similar constructors can be used for variables defined at faces of the mesh. Forexample, the script needed to define the mass flux (volume flow rate) field at cellfaces is shown in Listing 7.48.surfaceScalarField mdot(IOobject("mdot",runTime.timeName(),mesh,IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),linearInterpolate(U) & blockMesh.Sf());Listing 7.48 Script used to define the mass flux field at mesh facesIn this case the field is constructed with a scalar product between velocity(interpolated at the face) and face area vectors.
This field represents the volumeflow rate at control volume faces. It also represents the mass flow rate field forincompressible flows with a density value of one.Once fields are defined, access to specific data in different parts of the mesh ispossible as shown in the scripts below.7.2.2 InternalField DataInternal field data can be accessed (Listing 7.49) using// internal accessforAll(T.internalField(), cellI){scalar cellT = T.internalField()[cellI];}// internal accessforAll(U.internalField(), cellI){vector cellU = U.internalField()[cellI];}Listing 7.49 Accessing internal field data2007 The Finite Volume Mesh in OpenFOAM® and uFVM7.2.3 BoundaryField DataBoundary field data can be accessed (Listing 7.50) usingconst volVectorField::GeometricBoundaryField& UBoundaryList =U.boundaryField();// boundary accessforAll(UBoundaryList, patchI){const fvPatchField<vector>& fieldBoundary = UBoundaryList[patchI];forAll(fieldBoundary, faceI){vector faceU = fieldBoundary[faceI];}}Listing 7.50 Accessing boundary field dataor in a more compact form via (Listing 7.51)// boundary accessforAll(T.boundaryField(), patchI){forAll(T.boundaryField()[patchI], faceI){scalar faceT = T.boundaryField()[patchI][faceI];}}Listing 7.51 A more compact script to access boundary field data7.2.4 lduAddressingOpenFOAM® uses exclusively face addressing in its discretization loops andcoefficients storage.
It also uses arbitrary polyhedral elements in its meshes.A polyhedron element can have any number of faces with a neighbor elementassociated with each interior face. The storage of coefficients in OpenFOAM® isbased on the face addressing scheme. In this approach, coefficients are storedfollowing the interior face ordering, with access to elements and their coefficients,based on the owner/neighbor indices associated with interior faces.
As mentioned inthe previous section, the element with the lower index is the owner while theneighbor is the element with the higher index. For boundary faces the owner isalways the cell to which the face is attached while no neighbor is defined by settingthe neighbour index to −1. The list of owner or neighbor indices thus define theorder in which the element-to-element coefficients are assembled for the variousintegral operators.7.2 OpenFOAM®201The above described scheme is denoted by lduAddressing and is implementedin the lduMatrix class displayed in Fig.
7.10 that includes 5 arrays representing thediagonal, upper, and lower coefficients and the lower and upper indices of the faceowner and neighbor, respectively. In this scheme, owners represent the lower triangular part of the matrix (lower addressing) while neighbors refer to the uppertriangular part (upper addressing). It is worth mentioning that given a face for anowner element the lower and upper addressing provide respectively the column androw where the coefficient of the face flux is stored in the matrix while it is theopposite for a neighbor cell. For the domain shown on the upper left side ofFig.
7.10, the owner and neighbor of each face are displayed in the lowerAddr()and upperAddr() array, respectively. The face number is the position of the owneror neighbor in the array plus one (since numbering starts with 0). Consideringinternal face number 4, for example, its related information is stored in the fifth row(in C++ indexing starts from 0) of the lower(), upper(), lowerAddr(), andupperAddr() array, respectively.
The stored data indicate that its owner is elementnumber 2 (lowerAddr() array), its neighbor is element number 4 (upperAddr()array), the coefficient multiplying /4 in the algebraic equation for element number 2Fig. 7.10 The lduAddressing and lduMatrix2027 The Finite Volume Mesh in OpenFOAM® and uFVMis stored in the fifth row of the array upper(), and the coefficient multiplying /2 inthe algebraic equation for element number 4 is stored in the fifth row of the arraylower().Thus the lduAdressing provides information about the addresses of the offdiagonal coefficients in relation to the faces to which they are related. This meansthat while the computational efficiency for various operations on the matrix is highwhen they are mainly based on loops over all faces of the mesh, direct access to aspecific row-column matrix element is difficult and inefficient.
One example is thesummation of the off-diagonal coefficients for each row given byXac ¼anð7:1ÞnnbðCÞIn this case the use of face addressing does not allow direct looping over theoff-diagonal elements of each row and performing such operation requires loopingover all elements because it can only be done, as shown in Listing 7.52, following aface based approach.for (label faceI=0; faceI<l.size(); faceI++){ac[l[faceI]] -= Lower[faceI];ac[u[faceI]] -= Upper[faceI];}Listing 7.52 Summation of the off-diagonal coefficientsIn Listing 7.52 ac is the sum of the off-diagonals coefficients, l and u are theupper and lower addressing while Lower and Upper stores the correspondingcoefficients.
As can be noticed the summation is performed looping through all thefaces and the summation of each row is not sequential, depending only on theowner-neighbor numeration of the mesh.So in general the lduAddressing introduces a more complex handling of thematrix operations and it will be reflected also in the implementation of linear solversbut with the advantage of higher computational speed.7.2.5 Computing the GradientThe script used to compute the Green-Gauss gradient in OpenFOAM® is shown inListing 7.53.7.2 OpenFOAM®Foam::fv::gaussGrad<Type>::gradf(const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,const word& name){typedef typename outerProduct<vector, Type>::type GradType;const fvMesh& mesh = ssf.mesh();tmp<GeometricField<GradType, fvPatchField, volMesh> > tgGrad(new GeometricField<GradType, fvPatchField, volMesh>(IOobject(name,ssf.instance(),mesh,IOobject::NO_READ,IOobject::NO_WRITE),mesh,dimensioned<GradType>("0",ssf.dimensions()/dimLength,pTraits<GradType>::zero),zeroGradientFvPatchField<GradType>::typeName));GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();const labelUList& owner = mesh.owner();const labelUList& neighbour = mesh.neighbour();const vectorField& Sf = mesh.Sf();//Field<GradType>& igGrad = gGrad;const Field<Type>& issf = ssf;forAll(owner, facei){GradType Sfssf = Sf[facei]*issf[facei];igGrad[owner[facei]] += Sfssf;igGrad[neighbour[facei]] -= Sfssf;}forAll(mesh.boundary(), patchi){const labelUList& pFaceCells =mesh.boundary()[patchi].faceCells();const vectorField& pSf = mesh.Sf().boundaryField()[patchi];const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];forAll(mesh.boundary()[patchi], facei){igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];}}igGrad /= mesh.V();gGrad.correctBoundaryConditions();return tgGrad;}Listing 7.53 Script used to compute a gradient field in OpenFOAM®2032047 The Finite Volume Mesh in OpenFOAM® and uFVMThe small introduction reported above is neither intended as a replacement ofOpenFOAM® user’s manual [1] nor a replacement of a C++ manual.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
