Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 47
Текст из файла (страница 47)
Theinput to the function is theInterpolationScheme, the element gradient grad, theelement array phi and mdot for cases where an upwind or downwind scheme isused. The function is shown in Listing 9.3.2949Gradient Computationfunction grad_f=cfdInterpolateGradientsFromElementsToInteriorFaces(theInterpolationScheme,grad,phi,mdot_f)%===================================================% written by the CFD Group @ AUB, Fall 2006%===================================================theMesh= cfdGetMesh;numberOfInteriorFaces = theMesh.numberOfInteriorFaces;iOwners = [theMesh.faces(1:numberOfInteriorFaces).iOwner]';iNeighbours = [theMesh.faces(1:numberOfInteriorFaces).iNeighbour]';gf = [theMesh.faces(1:numberOfInteriorFaces).gf]';if(strcmp(theInterpolationScheme,'Average')==1)grad_f(:,1) = (1-gf).*grad(iNeighbours,1) + gf.*grad(iOwners,1);grad_f(:,2) = (1-gf).*grad(iNeighbours,2) + gf.*grad(iOwners,2);grad_f(:,3) = (1-gf).*grad(iNeighbours,3) + gf.*grad(iOwners,3);elseif(strcmp(theInterpolationScheme,'Upwind')==1)pos = zeros(size(mdot_f));pos((mdot_f>0))=1;%grad_f(:,1) = pos.*grad(iNeighbours,1) + (1-pos).*grad(iOwners,1);grad_f(:,2) = pos.*grad(iNeighbours,2) + (1-pos).*grad(iOwners,2);grad_f(:,3) = pos.*grad(iNeighbours,3) + (1-pos).*grad(iOwners,3);elseif(strcmp(theInterpolationScheme,'Downwind')==1)pos = zeros(size(mdot_f));pos((mdot_f>0))=1;%grad_f(:,1) = (1-pos).*grad(iNeighbours,1) + pos.*grad(iOwners,1);grad_f(:,2) = (1-pos).*grad(iNeighbours,2) + pos.*grad(iOwners,2);grad_f(:,3) = (1-pos).*grad(iNeighbours,3) + pos.*grad(iOwners,3);elseif(strcmp(theInterpolationScheme,'Average:Corrected')==1)grad_f(:,1) = (1-gf).*grad(iNeighbours,1) + gf.*grad(iOwners,1);grad_f(:,2) = (1-gf).*grad(iNeighbours,2) + gf.*grad(iOwners,2);grad_f(:,3) = (1-gf).*grad(iNeighbours,3) + gf.*grad(iOwners,3);d_CF=[theMesh.elements(iNeighbours).centroid]'[theMesh.elements(iOwners).centroid]';dmag=cfdMagnitude(d_CF);-e_CF(:,1)=d_CF(:,1)./dmag;e_CF(:,2)=d_CF(:,2)./dmag;e_CF(:,3)=d_CF(:,3)./dmag;local_grad_mag_f = (phi(iNeighbours)-phi(iOwners))./dmag;local_grad(:,1) = local_grad_mag_f.*e_CF(:,1);local_grad(:,2) = local_grad_mag_f.*e_CF(:,2);local_grad(:,3) = local_grad_mag_f.*e_CF(:,3);local_avg_grad_mag = dot(grad_f',e_CF')';local_avg_grad(:,1)=local_avg_grad_mag.*e_CF(:,1);local_avg_grad(:,2)=local_avg_grad_mag.*e_CF(:,2);local_avg_grad(:,3)=local_avg_grad_mag.*e_CF(:,3);grad_f = grad_f - local_avg_grad + local_grad;elsetheInterpolationSchemegrad,phi,mdotexit;endListing 9.3 Function used to interpolation the element gradient field to the faces9.5 Computational Pointers295It should be noted that theInterpolationScheme “Average:Corrected” implements the face gradient correction technique used to get a more accurate gradientrepresentation along the CF direction, i.e., Eq.
(9.33).9.5.2OpenFOAM®In OpenFOAM® [8] several types of gradient evaluation techniques are defined: thestandard Green-Gauss method, the second order least square method, and the fourthorder least square method. Attention will be focussed here on the Green-Gaussmethod. Nonetheless the discretization in all cases is performed explicitly and isthus part of the fvc operator namespace. The Green-Gauss gradient is defined at thecell centre, as in Eq.
(9.4), and its source code in OpenFOAM is located in thedirectory “src/finiteVolume/finiteVolume/gradSchemes/gaussGrad”.Implementation-wise, the gradient evaluation is performed in the following twosteps:• Face interpolation of the variable• Green-Gauss formula evaluationThe interpolation to the face is in the calcGrad routine listed below (Listing 9.4).Foam::fv::gaussGrad<Type>::calcGrad(const GeometricField<Type, fvPatchField, volMesh>& vsf,const word& name) const{typedef typename outerProduct<vector, Type>::type GradType;tmp<GeometricField<GradType, fvPatchField, volMesh> > tgGrad(gradf(tinterpScheme_().interpolate(vsf), name));GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();correctBoundaryConditions(vsf, gGrad);return tgGrad;}Listing 9.4 Script to interpolate variable values to cell facesAs a first step, values are interpolated to the faces and stored in the generic field“vsf” (the interpolation class will be described with more details in Chap.
11). Thenthe gradient is computed based on the Green-Gauss formula in the gradf routineshown in Listing 9.5.2969Gradient ComputationFoam::fv::gaussGrad<Type>::gradf(const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,const word& name){Info << "Calculating gradient for " << name << endl;//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 9.5 Routine used to compute the gradient using the Green-Gauss method9.5 Computational Pointers297The sum over cell faces is performed using the LDU addressing.
As such the“for” loop that evaluates the sum over the faces of the cell is based only on theglobal face numbering and uses the upper and lower addressing vectors to add orsubtract (Listing 9.6) the flux to cell values. After all fluxes have been processed,the net value is divided by the volume to yield the gradient, as in Eq. 9.4.igGrad[owner[facei]] += Sfssf;igGrad[neighbour[facei]] -= Sfssf;Listing 9.6 Adding or subtracting fluxes to cell valuesThe type of gradient is defined in fvSchemes shown in Listing 9.7.gradSchemes{defaultgrad(phi)}none;Gauss;Listing 9.7 Defining the gradient calculation methodThe face interpolation scheme is defined as displayed in Listing 9.8.interpolationSchemes{interpolate(phi) linear;}Listing 9.8 Defining the interpolation method used in calculating face valuesThe idea of evaluating the gradient based on a generic interpolation scheme thathas to be defined by dictionary, allows computing the gradient with theGreen-Gauss formula in several ways by just changing the interpolation scheme.If skew correction, as defined in Eq.
(9.8), is required then the above interpolation scheme definition should be replaced by (Listing 9.9)interpolationSchemes{interpolate(phi) skewCorrected linear;}Listing 9.9 Defining the interpolation scheme to calculate face values with skew correction2989Gradient ComputationA more compact syntax can be used for defining the gradient in which theinterpolation type is specified directly under gradSchemes shown in Listing 9.10.gradSchemes{defaultgrad(phi)}none;Gauss linear;Listing 9.10 Defining the gradient calculation method: compact syntaxIn this case the interpolation method is defined directly in the gradient dictionary. The choice of the syntax type is up to the user keeping in mind that a separatedefinition of the interpolation scheme helps clarifying the various steps of thegradient calculation.9.6ClosureThis chapter presented the discretization details of two methods for computing thegradient at the centroids of control volume meshes in general non-orthogonal gridsystems.
One method is based on the Green-Gauss theorem and the second on theLeast-Square reconstruction approach. Chapter 10 will be devoted to methods usedfor solving systems of algebraic equations.9.7ExercisesExercise 1In the configuration depicted in Fig. 9.22, the values of the variable / and its gradientr/ at C are known. Using these known values estimate the value of / at F.F( 4,5.5 )C( 2,1)C(= 10.3)C = 2i + 3.5 jFig. 9.22 A two dimensional element of centroid C and neighbor F9.7 Exercises299Exercise 2For the two cells shown in Fig. 9.23 the value of some scalar ϕ at their centroidsC and F can be calculated from/ðx; yÞ ¼ 100 x2 y þ y2 xThe gradient at C and F are also computed to be the exact gradient of thefunction at these points such thatr/ðx; yÞ ¼ 100 2xy þ y2 i þ 100 x2 þ 2yx jWith C and F located at (0.32, −0.1) and (0.52, 0.31), respectively, find the value ofthe gradient at face f (0.43, 0.1) numerically (not from the expression for thegradient) usinga. A simple averaging between C and F.b.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
