Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 64
Текст из файла (страница 64)
Theowner and neighbor indices for the interior faces are first retrieved and an upwindindex is defined (pos). Then the coefficients for the upwind scheme are computed.theMdotName = ['Mdot' theFluidTag];mdotField = cfdGetMeshField(theMdotName,'Faces');mdot_f = mdotField.phi(iFaces);iOwners = [theMesh.faces(iFaces).iOwner];iNeighbours = [theMesh.faces(iFaces).iNeighbour];pos = zeros(size(mdot_f));pos((mdot_f>0))=1;theFluxes.FLUXC1f(iFaces,1) = mdot_f.*pos;theFluxes.FLUXC2f(iFaces,1) = mdot_f.*(1-pos);theFluxes.FLUXVf(iFaces,1) = 0;Listing 11.1 Convection scheme class declarationThe implementation of High Order schemes is performed after the upwindassembly using the deferred correction method. The assembly of the QUICKscheme (cfdAssembleConvectionTermDCQUICK) is shown in Listing 11.2.iUpwind = pos.*iOwners + (1-pos).*iNeighbours;%get the upwind gradient at the interior facesphiGradCf = phiGrad(iUpwind,:,iComponent);%interpolated gradient to interior facesiOwners = [theMesh.faces(iFaces).iOwner]';iNeighbours = [theMesh.faces(iFaces).iNeighbour]';pos = zeros(size(mdot_f));pos(mdot_f>0)=1;%phiGradf=cfdInterpolateGradientsFromElementsToInteriorFaces('Average:Corrected',phiGrad,phi);rc = [theMesh.elements(iUpwind).centroid]';rf = [theMesh.faces(iFaces).centroid]';rCf = rf-rc;%corr = mdot_f .* dot(phiGradCf'+phiGradf',rCf')'*0.5;%theFluxes.FLUXTf(iFaces) = theFluxes.FLUXTf(iFaces) + corr;Listing 11.2 Assembly of the QUICK scheme11.9Computational Pointers413The deferred correction value is computed for all interior faces as10CB1UPWIND/C Cm_ f /HR¼ m_ f Bf /fA@/C þ 2 r/C þ r/f dCf |{z}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} UPWINDQUICK¼ m_ fð11:162Þ1r/C þ r/f dCf211.9.2 OpenFOAM®In OpenFOAM® [18] the convection term can be evaluated either explicitly usingfvc::div(mDot, phi) or implicitly via the fvm::div(mDot,phi) function.
The fvc::div(mDot,phi) returns a field in which the divergence of phi is evaluated in eachcell. The field is added to the right hand side of the system of equations. The fvm::div(mDot,phi) returns the fvMatrix, a matrix of coefficients that are evaluatedbased on the linearization of the face fluxes. The matrix of coefficients is then addedto the left hand side of the system of equations.The scripts of fvm::div and fvc::div functions can be found in the file FOAM_SRC/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionSheme.C.
As already stated, the convection class is based on thetype of interpolation scheme at the face and accordingly the declaration of the classis performed as displayed in Listing 11.3.template<class Type>class gaussConvectionScheme:public fv::convectionScheme<Type>{// Private datatmp<surfaceInterpolationScheme<Type> > tinterpScheme_;Listing 11.3 Convection scheme class declarationIn the above declaration the private member tinterpScheme_ describes the faceinterpolation type on which the divergence operator is based. Additionally to theabove mentioned function, the gaussConvectionScheme class defines two auxiliary functions named interpolate and flux, as shown in Listing 11.4, which wrapthe surfaceInterpolation class used to perform interpolation from volume fields tosurface fields.41411 Discretization of the Convection Termtemplate<class Type>tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >gaussConvectionScheme<Type>::interpolate(const surfaceScalarField&,const GeometricField<Type, fvPatchField, volMesh>& vf) const{return tinterpScheme_().interpolate(vf);}template<class Type>tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >gaussConvectionScheme<Type>::flux(const surfaceScalarField& faceFlux,const GeometricField<Type, fvPatchField, volMesh>& vf) const{return faceFlux*interpolate(faceFlux, vf);}Listing 11.4 The gaussConvectionScheme class that defines the interpolate and flux functionsFor the explicit discretization of convection, the operator defines a field wherethe divergence of the cell face fluxes are stored.
The code used by the fvc::divoperator for the explicit evaluation of the divergence of a field over a specificvolume is as follows (Listing 11.5):template<class Type>tmp<GeometricField<Type, fvPatchField, volMesh> >gaussConvectionScheme<Type>::fvcDiv(const surfaceScalarField& faceFlux,const GeometricField<Type, fvPatchField, volMesh>& vf) const{tmp<GeometricField<Type, fvPatchField, volMesh> > tConvection(fvc::surfaceIntegrate(flux(faceFlux, vf)));tConvection().rename("convection(" + faceFlux.name() + ',' + vf.name() + ')');return tConvection;}Listing 11.5 Script used for the explicit evaluation of the divergence operator11.9Computational Pointers415The calculation of the divergence of a field involves the following three steps:1.
Evaluating values at the faces of the element.2. Multiplying the value at the face with the mass flux at the face (i.e., faceFlux).3. Summing the contributions of all cell faces and dividing the sum by the cellvolume.First the face value of the generic variable vf is evaluated. Then the estimateobtained is multiplied by the corresponding face mass flux value using auxiliaryfunctions. Finally the sum over the faces of the cell is performed using thesurfaceIntegrate function.
The implementation of this function can be foundunder FOAM_SRC/finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.C andits script is given in Listing 11.6.template<class Type>void surfaceIntegrate(Field<Type>& ivf,const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf){const fvMesh& mesh = ssf.mesh();const labelUList& owner = mesh.owner();const labelUList& neighbour = mesh.neighbour();const Field<Type>& issf = ssf;forAll(owner, facei){ivf[owner[facei]] += issf[facei];ivf[neighbour[facei]] -= issf[facei];}forAll(mesh.boundary(), patchi){const labelUList& pFaceCells =mesh.boundary()[patchi].faceCells();const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];forAll(mesh.boundary()[patchi], facei){ivf[pFaceCells[facei]] += pssf[facei];}}ivf /= mesh.V();}Listing 11.6 Script of the surfaceIntegrate function41611 Discretization of the Convection TermThe implicit discretization of the convection term is performed using the fvm::div operator.
It does so by implicitly expressing the dependent variable value at theface as function of the owner and neighbor values according to- /O þ ð1 -Þ /N/f ¼ |{z}|fflfflfflffl{zfflfflfflffl}ownercoefficientð11:163Þneighbourcoefficientwhere subscripts O and N refer to owner and neighbor, respectively, and - representsthe weight assigned to the contribution of the owner to the value at the face.
The waythe weight is calculated will be described in the next chapter. The coefficients of /Oand /N are then used to assemble the matrix of coefficients. The script used toperform implicit discretization of the divergence operator reads (Listing 11.7)template<class Type>tmp<fvMatrix<Type> >gaussConvectionScheme<Type>::fvmDiv(const surfaceScalarField& faceFlux,const GeometricField<Type, fvPatchField, volMesh>& vf) const{tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);const surfaceScalarField& weights = tweights();tmp<fvMatrix<Type> > tfvm(new fvMatrix<Type>(vf,faceFlux.dimensions()*vf.dimensions()));fvMatrix<Type>& fvm = tfvm();Listing 11.7 Implicit calculation of the divergence of a fieldwhere in a first place an fvMatrix is defined and then, as shown in Listing 11.8, it isproperly filled with the corresponding coefficients asfvm.lower() = -weights.internalField()*faceFlux.internalField();fvm.upper() = fvm.lower() + faceFlux.internalField();fvm.negSumDiag();Listing 11.8 Properly assembling the weights contributions to coefficients11.9Computational Pointers417clippedLinear<type>CoBlended<type>downwind<type>harmonic<scalar>fixedBlended<type>limitedSurfaceInterpolationScheme<type>limiterBlended<type>limitWith<type>linear<type>localBlended<type>refCountsurfaceInterpolationScheme<type>localMax<type>localMin<type>midPoint<type><GType>fieldScheme<type><scalar>surfaceInterpolationScheme<scalar>surfaceInterpolationScheme<GType>weighted<type>skewCorrected<type>outletStabilised<type>reverseLinearFig.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
