Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 42
Текст из файла (страница 42)
To promote convergence and stabilize theiterative solution process, slowing down the changes in ϕ between iterations ishighly desirable and is enforced by a technique called under-relaxation. There aremany ways of introducing under-relaxation. One of the practices will be describedhere, others in Chap. 14. Derivations will be performed on the general discretizationequation of the formXaC /C þaF /F ¼ bCð8:95ÞFNBðCÞEquation (8.95) can be written as/C ¼PaF / F þ bCFNBðCÞaCð8:96ÞLet /C represents the value of ϕC from the previous iteration. If /C is added toand subtracted from the right hand side, then Eq. (8.96) becomesP10aF / F þ bCCB FNBðCÞð8:97Þ/C ¼ /C þ @ /C AaCwhere the expression between the parentheses represents the change in ϕC producedby the current iteration.
This change can be modified by the introduction of arelaxation factor λϕ, such that8.9 Under-Relaxation of the Iterative Solution Process0B/C ¼ /C þ k/ @P257aF /F þ bCFNBðCÞaC1C /C Að8:98ÞoraCk/ þ/ CXaF / F ¼ b C þ1 k/ a CFNBðCÞk//Cð8:99ÞAt first, it should be noted that at convergence ϕC and /C become equal independent of the value of relaxation factor used. This is indeed reflected by Eq. (8.98),which shows that the converged value ϕC does satisfy the original equation(Eq. 8.95). This property should be satisfied by any relaxation scheme.Depending on the value of the relaxation factor λϕ, the equation may be eitherunder relaxed (i.e., 0 < λϕ < 1) or over-relaxed (i.e., λϕ > 1).
In CFD applications,under relaxation is usually used. A value of λϕ close to 1 implies little underrelaxation, while a value close to 0 produces heavy under relaxation effects withvery small changes in ϕC from iteration to iteration.The optimum under relaxation factor is problem dependent and is not governedby any general rule. The factors affecting λϕ values include the type of problemsolved, the size of the system of equation (i.e., number of grid points in the domain),the grid spacing and its expansion rate, and the adopted iterative method, amongothers. Usually, values of λϕ are assigned based on experience or from preliminarycalculations. Moreover, it is not necessary to use the same under-relaxation valuethroughout the computational domain and values may vary from iteration toiteration.Equation (8.99) can be recast into the form of Eq.
(8.95), where now the centralcoefficient ac becomesaCaCð8:100Þk/while the source term is added to produce a new term asbCbC þ1 k/ a Ck//Cð8:101ÞThis method of relaxation plays an important role in stabilizing the solution ofnon-linear problems.2588.108.10.18Spatial Discretization: The Diffusion TermComputational PointersuFVMIn uFVM the implementation of the diffusion term discretization for interior faces isperformed in function cfdAssembleDiffusionTermInterior, the core of which isshown in Listing 8.1.gamma_f = cfdInterpolateFromElementsToFaces('Average',gamma);gamma_f = [gamma_f(iFaces)];geoDiff_f = [theMesh.faces(iFaces).geoDiff]';Sf = [theMesh.faces(iFaces).Sf]';Tf = [theMesh.faces(iFaces).T]';iOwners = [theMesh.faces(iFaces).iOwner]';iNeighbours = [theMesh.faces(iFaces).iNeighbour]';theFluxes.FLUXC1f(iFaces,1) = gamma_f .* gDiff_f;theFluxes.FLUXC2f(iFaces,1) = -gamma_f .* gDiff_f;theFluxes.FLUXVf(iFaces,1) = gamma_f .* dot(grad_f(:,:)',Tf(:,:)')';theFluxes.FLUXTf(iFaces,1)=theFluxes.FLUXC1f(iFaces).*phi(iOwners)+theFluxes.FLUXC2f(iFaces).*phi(iNeighbours)+theFluxes.FLUXVf(iFaces);Listing 8.1 Assembly of the diffusion term at interior facesIn the above, FLUXC1f and FLUXC2f are equivalent to FluxCf and FluxFf inthe text and are the coefficients of the owner and neighbor element, respectively.Moreover, gamma_f and gDiff_f are arrays defined over all interior faces, andusing the dot notation, the expression gamma_f.* gDiff_f returns an arraycontaining the product of the respective elements of the gamma_f and gDiff_farrays.
The dot notation in Matlab® allows for efficient operations over arrays andfor the writing of a clearer code and is used in uFVM whenever practical.The non-orthogonal term is stored in the FLUXVf coefficient and is equal togamma_f .* dot(grad_f(: , :)’,Tf( :, :)’)’.The total flux passingthrough face f, FLUXTf, is assembled as in Eq. (8.15).The diffusion term is setup in cfdProcessOpenFoamMesh.m and Listing 8.2shows the computation of the gDiff coefficient.8.10Computational Pointers259theMesh.faces(iFace).dCF = element2.centroid - element1.centroid;eCF = dCF/cfdMagnitude(dCF);E = theFace.area*eCF;theMesh.faces(iFace).gDiff = cfdMagnitude(E)/cfdMagnitude(dCF);theMesh.faces(iFace).T = theFace.Sf - E;Listing 8.2 Computing the geometric diffusion coefficientThe effects of boundary conditions on the equations of boundary elementsshould be accounted for and Listing 8.3 shows the assembly of the diffusion term atboundary faces for the case of a Dirichlet boundary condition type.theMesh = cfdGetMesh;numberOfElements = theMesh.numberOfElements;numberOfInteriorFaces = theMesh.numberOfInteriorFaces;%theBoundary = theMesh.boundaries(iPatch);numberOfBFaces = theBoundary.numberOfBFaces;%iFaceStart = theBoundary.startFace;iFaceEnd = iFaceStart+numberOfBFaces-1;iBFaces = iFaceStart:iFaceEnd;%iElementStart = numberOfElements+iFaceStart-numberOfInteriorFaces;iElementEnd = iElementStart+numberOfBFaces-1;iBElements = iElementStart:iElementEnd;geodiff = [theMesh.faces(iBFaces).geoDiff]';Tf = [theMesh.faces(iBFaces).T]';iOwners = [theMesh.faces(iBFaces).iOwner]';theFluxes.FLUXC1f(iBFaces) =gamma(iBElements).*geodiff;theFluxes.FLUXC2f(iBFaces) = - gamma(iBElements).*geodiff;theFluxes.FLUXVf(iBFaces)=gamma(iBElements).*dot(grad(iBElements,:)’,Tf(:,:)’)’;-Listing 8.3 Assembly of the diffusion term at boundary faces for a Dirichlet ConditionThe implementation of other boundary condition types can be reviewed in filecfdAssembleDiffusionTerm.m.Once the linearized coefficients are computed for each of the interior andboundary faces, then assembly into the global (sparse) matrix can proceed.
Thisapproach is used for the discretization of all flux terms in uFVM and the assembly isperformed in the cfdAssembleIntoGlobalMatrixFaceFluxes function.Listing 8.4 shows the assembly into the sparse LHS matrix and the RHS vector of thecoefficients at interior faces.2608Spatial Discretization: The Diffusion Term%% Assemble fluxes of interior faces%for iFace = 1:numberOfInteriorFacestheFace = theMesh.faces(iFace);iOwner= theFace.iOwner;iOwnerNeighbourCoef= theFace.iOwnerNeighbourCoef;iNeighbour= theFace.iNeighbour;iNeighbourOwnerCoef = theFace.iNeighbourOwnerCoef;%% assemble fluxes for owner cellac(iOwner) = ac(iOwner)+ vf_f(iFace)*theFluxes.FLUXC1f(iFace);anb{iOwner}(iOwnerNeighbourCoef) = anb{iOwner}(iOwnerNeighbourCoef)+ vf_f(iFace)*theFluxes.FLUXC2f(iFace);bc(iOwner) = bc(iOwner)- vf_f(iFace)*theFluxes.FLUXTf(iFace);%% assemble fluxes for neighbour cellac(iNeighbour) = ac(iNeighbour)- vf_f(iFace)*theFluxes.FLUXC2f(iFace);anb{iNeighbour}(iNeighbourOwnerCoef)=anb{iNeighbour}(iNeighbourOwnerCoef)- vf_f(iFace)*theFluxes.FLUXC1f(iFace);bc(iNeighbour)= bc(iNeighbour)+ vf_f(iFace)*theFluxes.FLUXTf(iFace);endListing 8.4 Assembly into LHS spare matrix and RHS arrayFor each interior face the coefficients are assembled into both the owner andneighbor element equations.
For the owner the coefficients are (ac(iOwner),anb(iOwner)(iOwnerNeighbourCoef), and bc(iOwner)), while for theneighbor the coefficients become: ac(iNeighbour), anb(iNeighbour)(iNeighbourOwnerCoef), and bc(iNeighbour)). The different signs usedin the assembly of the two equations is due to the fact that the face surface vector ispointing into the neighbor cell and out of the owner cell.8.10.2OpenFOAM®In OpenFOAM® [18] the diffusion term can be evaluated explicitly using “fvc::laplacian(gamma, phi)”or implicitly via “fvm::laplacian(gamma,phi)” namespacesand functions.
The “fvc::laplacian(gamma,phi)” returns a field in which the laplacianof a generic field ϕ (phi) is evaluated at each cell. The field is added to the right handside of the system of equations. The “fvm::laplacian(gamma,phi)” returns instead anfvMatrix of coefficients evaluated as per Eq. (8.19), which is added to the left hand sideof the system of equations, in addition to a field containing the non-orthogonal termsthat is added to the right hand side of the system.
The definition of the laplacianoperator is located in the directory “$FOAM_SRC/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme” in the files “gaussLaplacianScheme.C”,“gaussLaplacianScheme.H”, and “gaussLaplacianSchemes.C”.8.10Computational PointersThe “fvm::laplacian” definition displayed in Listing 8.5 readstemplate<>Foam::tmp<Foam::fvMatrix<Foam::Type> >Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian(const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,const GeometricField<Type, fvPatchField, volMesh>& vf){const fvMesh& mesh = this->mesh();GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf(gamma*mesh.magSf());tmp<fvMatrix<Type> > tfvm = fvmLaplacianUncorrected(gammaMagSf,this->tsnGradScheme_().deltaCoeffs(vf),vf);fvMatrix<Type>& fvm = tfvm();if (this->tsnGradScheme_().corrected()){if (mesh.fluxRequired(vf.name())){fvm.faceFluxCorrectionPtr() = newGeometricField<Type, fvsPatchField, surfaceMesh>(gammaMagSf*this->tsnGradScheme_().correction(vf));fvm.source() -=mesh.V()*fvc::div(*fvm.faceFluxCorrectionPtr())().internalField();}else{fvm.source() -=mesh.V()*fvc::div(gammaMagSf*this->tsnGradScheme_().correction(vf))().internalField();}}return tfvm;}Listing 8.5 Calculation of the Laplacian operator2612628Spatial Discretization: The Diffusion Termwith the following additional functions (Listing 8.6):template<class Type, class GType>tmp<fvMatrix<Type> >gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected(const surfaceScalarField& gammaMagSf,const surfaceScalarField& deltaCoeffs,const GeometricField<Type, fvPatchField, volMesh>& vf){tmp<fvMatrix<Type> > tfvm(new fvMatrix<Type>(vf,deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()));fvMatrix<Type>& fvm = tfvm();fvm.upper()= deltaCoeffs.internalField()*gammaMagSf.internalField();fvm.negSumDiag();forAll(vf.boundaryField(), patchi){const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];const fvsPatchScalarField& pDeltaCoeffs =deltaCoeffs.boundaryField()[patchi];if (pvf.coupled()){fvm.internalCoeffs()[patchi] =pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);fvm.boundaryCoeffs()[patchi] =-pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);}else{fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();}}return tfvm;}Listing 8.6 Additional functions needed for the calculation of the Laplacian operatorIt is worth noting that in OpenFOAM® the assembly takes place directly into theglobal coefficients, which as described earlier in Chaps.
5, 6, and 7 are stored inthree arrays, namely fvm.upper(), fvm.lower(), and fvm.diag(). Themain part of the discretization is defined in the fvmLaplacianUncorrectedfunction, where Eq. (8.19) is evaluated by first defining an fvMatrix objectand then filling its upper triangle by the extra diagonal coefficients (this approachrelies on the fact that the Laplacian operator returns a symmetric matrix).8.10Computational Pointers263The “deltaCoeffs.internalField()” represents the gDiff field whilegamma indicates the diffusion field. The diagonal coefficients are evaluated infvm.negSumDiag(),where the negative sum of the upper and lower coefficientsare assembled into the diagonal coefficients as per Eq. (8.18).The “fvm.negSumDiag()” method is defined in lduMatrixOperations.C as (Listing 8.7)void Foam::lduMatrix::negSumDiag(){const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();scalarField& Diag = diag();const labelUList& l = lduAddr().lowerAddr();const labelUList& u = lduAddr().upperAddr();for (register label face=0; face<l.size(); face++){Diag[l[face]] -= Lower[face];Diag[u[face]] -= Upper[face];}}Listing 8.7 Calculation of the negative sum of the upper and lower coefficientsThe boundary conditions are implemented exactly as in Eqs.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
