Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 72
Текст из файла (страница 72)
The usual practice atan outlet (Fig. 12.19) is to apply the upwind scheme ð/b ¼ /C Þ; which automatically results in a zero normal gradient.Discretizing the convection flux at interior faces using a HR scheme implemented through a deferred correction approach, the modified algebraic equation forthe boundary element can be written asXaC /C þaF / F ¼ b Cð12:111ÞF NBðC ÞSbnbebCFig. 12.19 Outlet boundary condition for the convection flux12.9Boundary Conditions471whereaF ¼ FluxFf ¼ km_ f ; 0 kXXFluxCf ¼km_ f ; 0 kaC ¼f nbðC Þ¼f nbðC ÞXaF þF NBðC ÞbC ¼ XX f nbðC ÞFluxVf ¼ f nbðC Þm_ f þ m_ bð12:112ÞXUm_ f /HRf /ff nbðC Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}bDCcwhere f refers to the interior faces of the boundary element, and C and F refer toowner and neighbor, respectively.12.9.3 Wall Boundary ConditionAs shown in Fig. 12.20, the normal velocity at a wall is zero.
As such the convection flux is zero and does not appear in the algebraic equation.Again adopting a HR scheme with a deferred correction approach, the modifiedalgebraic equation for the boundary element can be written asaC /C þXaF / F ¼ b Cð12:113ÞF NBðC ÞFig. 12.20 Wall boundarycondition for the convectionfluxSbnbebC47212 High Resolution SchemeswhereaF ¼ FluxFf ¼ km_ f ; 0 kXXFluxCf ¼km_ f ; 0 kaC ¼f nbðC Þf nbðC ÞX¼aF þFNBðC ÞXbC ¼ f nbðC ÞXm_ ff nbðC ÞFluxVf ¼ Xð12:114ÞUm_ f /HRf /ff nbðC Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}bDCcAgain f refers to the interior faces of the boundary element and C and F refer tothe owner and neighbor nodes, respectively.12.9.4 Symmetry Boundary ConditionNo flow crosses a symmetry boundary.
Therefore it is treated in a similar fashion toa wall boundary condition with the convection flux normal to a symmetry boundaryset to zero.12.1012.10.1Computational PointersuFVMSimilar to HO schemes discussed in Chap. 11, HR schemes are implemented inuFVM using the deferred correction method. The main difference in the implementation stems from the use of the NVF or TVD relations rather than the calculation of the face value using directly the gradient.
The implementation of theSTOIC HR [31] scheme using the NVF formulation is shown in Listing 12.1(cfdAssembleConvectionTermDCSTOIC).It starts with the retrieval of the needed fields, followed by setting the upwindand downing indices for all interior faces.
Then for each face, /C and /D are~ is computed and used to calculate the faceidentified, /U is constructed, and /Cvalue from the NVF relationship of the adopted HR scheme. Finally /f is recon~ and used in the deferred correction method.structed from /f12.10Computational Pointers473theFluidTag = cfdGetFluidTag(theEquationName);theMdotName = ['Mdot' theFluidTag];theMdotField = cfdGetMeshField(theMdotName,'Faces');mdot_f = theMdotField.phi(iFaces);iOwners = [theMesh.faces(iFaces).iOwner]';iNeighbours = [theMesh.faces(iFaces).iNeighbour]';pos = zeros(size(mdot_f));pos(mdot_f>0)=1;% find indices of U and D cellsiUpwind=pos.*iOwners + (1-pos).*iNeighbours;iDownwind = (1-pos).*iOwners +pos.*iNeighbours;% find phi_C, phi_D and calculate phi_Uphi_C = phi(iUpwind,iComponent);phi_D = phi(iDownwind,iComponent);rCD=[theMesh.elements(iDownwind).centroid]'[theMesh.elements(iUpwind).centroid]';phi_U = phi_D - 2*dot(phiGrad(iUpwind,:,iComponent)',rCD')';SMALL= 1e-6;% calculate phi_tildaCnominator = phi_C-phi_U;denominator = phi_D-phi_U;divideLoc = find(~((denominator<SMALL) & (denominator>-SMALL)));phi_tildaC = ones(size(phi_C));phi_tildaC(divideLoc) = nominator(divideLoc)./denominator(divideLoc);Listing 12.1 Implementation of the STOIC HR scheme47412 High Resolution Schemes% get phi_tildaf from STOIC functionphi_tildaf = zeros(size(phi_tildaC));% lower UPWIND sectionphi_tildaf0)=phi_tildaf.*(phi_tildaC);+(phi_tildaC<=% intermediate sectionphi_tildaf = phi_tildaf0.2).*(3*phi_tildaC);+(phi_tildaC>0).*(phi_tildaC<% CDS sectionphi_tildaf=phi_tildaf+0.5).*(0.5*phi_tildaC + 0.5);(phi_tildaC>=0.2).*(phi_tildaC<(phi_tildaC>=0.5).*(phi_tildaC<phi_tildaf=phi_tildaf+(phi_tildaC1) .*(ones(size(phi_tildaC)));>=5/6).*(phi_tildaC<% SMART sectionphi_tildaf=phi_tildaf+5/6).*(0.75*phi_tildaC + 3/8);% DOWNWIND section% upper UPWIND sectionphi_tildaf1)=phi_tildaf.*(phi_tildaC);+(phi_tildaC% calculate phi_fphi_f = phi_tildaf.*(phi_D-phi_U) + phi_U;% calculate correctioncorr = mdot_f .* (phi_f - phi_C);% apply deferred correctiontheFluxes.FLUXTf(iFaces) = theFluxes.FLUXTf(iFaces) + corr;Listing 12.1 (continued)>=12.10Computational Pointers12.10.2475OpenFOAM®As discussed in Chap.
11 the discretization of the convection term in OpenFOAM®[32] is accomplished through the base class surfaceInterpolationScheme. Thisclass is inherited by any face interpolation algorithm. High resolution schemes, asdiscussed in this chapter, are special face interpolation algorithms. Thus, as displayed in Fig. 12.21, OpenFOAM® groups all TVD schemes in a base class,derived from the surfaceInterpolationScheme class, denoted by thelimitedSurfaceInterpolationScheme class.As shown in Listing 12.2, this class consists of the following three main functions: the virtual weights function representing the main virtual function of thesurfaceInterpolationScheme class, a local weights function defined with threearguments, and a new virtual base function called limiter.SurfaceInterpolationScheme<type>LimitedSurfaceInterpolationScheme<type>blended<type>LimitedScheme<type>PhiScheme<type>upwind<type>Fig.
12.21 UML showing the base class that groups all TVD schemes47612 High Resolution Schemestemplate<class Type>class limitedSurfaceInterpolationScheme:public surfaceInterpolationScheme<Type>{// Member Functions//- Returns the interpolation weighting factorsvirtual tmp<surfaceScalarField> limiter(const GeometricField<Type, fvPatchField, volMesh>&) const = 0;//- Returns the interpolation weighting factors for the given//field, by limiting the given weights with the given limitertmp<surfaceScalarField> weights(const GeometricField<Type, fvPatchField, volMesh>&,const surfaceScalarField& CDweights,tmp<surfaceScalarField> tLimiter) const;//- Return the interpolation weighting factors for the given fieldvirtual tmp<surfaceScalarField> weights(const GeometricField<Type, fvPatchField, volMesh>&) const;Listing 12.2 The limitedSurfaceInterpolationScheme class showing its three main functionsThe script of the virtual weights function, shown in Listing 12.3, is given bytemplate<class Type>tmp<surfaceScalarField>limitedSurfaceInterpolationScheme<Type>::weights(const GeometricField<Type, fvPatchField, volMesh>& phiListing 12.3 Script showing the implementation of the virtual weights function12.10Computational Pointers477) const{return this->weights(phi,this->mesh().surfaceInterpolation::weights(),this->limiter(phi));}Listing 12.3 (continued)It is clear that executing the script instantiates the additional local weightsfunction.
The three main arguments of the local weights function are the flux phi,the linear interpolation weights (representing the central differencing weights), andthe returning object of the virtual base class limiter. These are clearly shown inListing 12.4 where the local weights function is defined.template<class Type>tmp<surfaceScalarField>limitedSurfaceInterpolationScheme<Type>::weights(const GeometricField<Type, fvPatchField, volMesh>& phi,const surfaceScalarField& CDweights,tmp<surfaceScalarField> tLimiter) const{surfaceScalarField& Weights = tLimiter();scalarField& pWeights = Weights.internalField();forAll(pWeights, face){pWeights[face] =pWeights[face]*CDweights[face]+ (1.0 - pWeights[face])*pos(faceFlux_[face]);}Listing 12.4 Script showing the implementation of the local weights function47812 High Resolution SchemesThe calculation of the interpolation weights according to the TVD formulationimplemented in Listing 12.4 may seem, from a first look, a bit unclear.
Confusionmay be eliminated by considering the case of a positive flux, for example, for whichthe weight - is calculated as- ¼ -CD w þ ð1 wÞð12:115ÞFor a uniform grid -CD ¼ 1=2 yielding-¼wwþ ð1 w Þ ¼ 1 22ð12:116ÞRecalling Eq. (11.164) and substituting the weight yieldsww/f ¼ /N þ -ð/O /N Þ ¼ /N þ 1 ð/O /N Þ ¼ /O þ ð/N /O Þ22ð12:117Þwhich is Eq. (12.30) in the TVD formulation.So far the implementation of HR schemes in OpenFOAM® following the TVDformulation has been discussed along with its integration in the standard convectiondiscretization procedure by properly defining the interpolation weights.
Due to themany TVD schemes that can be used, OpenFOAM® has introduced a base virtualfunction denoted by limiter for the implementation of these schemes. The last stepof the computational pointer is to describe how this class is defined and whichclasses, this limiter base class, engage.All TVD limiters are organized through a class named LimitedScheme thatinherits and defines the limiter function of the limitedSurfaceInterpolationScheme class. The definition of this class is based on nested template classesdefinition in which the derived class is described with a template argument, asshown in Listing 12.5.template<class Type, class Limiter, template<class> class LimitFunc>class LimitedScheme:public limitedSurfaceInterpolationScheme<Type>,public Limiter{Listing 12.5 The LimitedScheme class with the Limiter template12.10Computational Pointers479Here Limiter is not a function but just a template definition.
Additionally thevirtual limiter class is now specialized (Listing 12.6).//- Return the interpolation weighting factorsvirtual tmp<surfaceScalarField> limiter(const GeometricField<Type, fvPatchField, volMesh>&) const;Listing 12.6 The virtual limiter classAs depicted in Listing 12.7 the definition of the limiter function is practicallylinked to an auxiliary function named calcLimiter.template<class Type, class Limiter, template<class> class LimitFunc>Foam::tmp<Foam::surfaceScalarField>Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter(const GeometricField<Type, fvPatchField, volMesh>& phi) const{tmp<surfaceScalarField> tlimiterField(new surfaceScalarField(IOobject(limiterFieldName,mesh.time().timeName(),mesh),mesh,dimless));Listing 12.7 The script used to call the calcLimiter function48012 High Resolution SchemescalcLimiter(phi, tlimiterField());return tlimiterField;Listing 12.7 (continued)The core of the calcLimiter function is to evaluate the TVD limiter and to storeit in the tlimiterField, as shown in Listings 12.7 and 12.8.template<class Type, class Limiter, template<class> class LimitFunc>void Foam::LimitedScheme<Type, Limiter, LimitFunc>::calcLimiter(const GeometricField<Type, fvPatchField, volMesh>& phi,surfaceScalarField& limiterField) const{const fvMesh& mesh = this->mesh();tmp<GeometricField<typenamevolMesh> >Limiter::phiType,fvPatchField,tlPhi = LimitFunc<Type>()(phi);constvolMesh>&GeometricField<typenameLimiter::phiType,fvPatchField,Limiter::gradPhiType,fvPatchField,lPhi = tlPhi();tmp<GeometricField<typenamevolMesh> >tgradc(fvc::grad(lPhi));const GeometricField<typename Limiter::gradPhiType, fvPatchField,volMesh>&gradc = tgradc();constsurfaceScalarField&mesh.surfaceInterpolation::weights();Listing 12.8 Script used with the calcLimiter functionCDweights=12.10Computational Pointers481const labelUList& owner = mesh.owner();const labelUList& neighbour = mesh.neighbour();const vectorField& C = mesh.C();scalarField& pLim = limiterField.internalField();forAll(pLim, face){label own = owner[face];label nei = neighbour[face];pLim[face] = Limiter::limiter(CDweights[face],this->faceFlux_[face],lPhi[own],lPhi[nei],gradc[own],gradc[nei],C[nei] - C[own]);}Listing 12.8 (continued)In the calcLimiter function the following steps are performed:•••••Storing a copy of the field to be interpolated in lPhiEvaluating the gradient of the lPhi field using fvc::grad(lPhi)Collecting the central differencing weightsCollecting the cell centersEvaluating the limiter by calling the nested template class: pLim[face] = Limiter::limiter.As an example of a Limiter::limiter function, the TVD formulation of the SUPERBEEgiven in Eq.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
