Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 79
Текст из файла (страница 79)
(13.1) over the element C shown in Fig. 13.2 is obtained by substituting Eq. (13.104) in Eq. (13.57)and is given byDtVCDtDtVCðqC /C Þ 1 þþ ðq / Þ1þDt þ DtDt C CDt þ Dt DtDt þ DtDtVCþ ðq / Þ þLð/C Þ ¼ 0Dt þ Dt Dt C Cð13:105ÞThe linearization coefficients for the SOUE scheme with non uniform time steps areinferred to be11q VCþFluxC ¼Dt Dt þ Dt C11Dt =Dtq VCþFluxC ¼ þDt Dt þ Dt Dt þ Dt CDt =Dtq VC /FluxV ¼CDt þ Dt Cð13:106ÞSimilar to the constant time step case, the method is implicit as it requires solving asystem of equations to obtain the ϕ field at every time step. The uniform time stepform of the equation given by Eq.
(13.80) is obtained by setting Dt ¼ Dt ¼ Dt inEq. (13.105).13.513.5.1Computational PointersuFVMThe discretization of the transient term in uFVM follows the finite volume methodimplemented within an implicit framework. The assembly of the transient fluxesresulting from the first order backward Euler scheme is shown in Listing 13.1.52613 Temporal Discretization: The Transient Term%theDensityField = cfdGetMeshField(['Density' theFluidTag]);density = theDensityField.phi(iElements);density_old = theDensityField.phi_old(iElements);volumes = [theMesh.elements(iElements).volume]’;theFluxes.FLUXCE(iElements)theFluxes.FLUXCEOLD(iElements)theFluxes.FLUXTE(iElements)theFluxes.FLUXTEOLD(iElements)=volumes .* density / dt;= - volumes .* density_old / dt;= theFluxes.FLUXCE .* phi;= theFluxes.FLUXCEOLD .* phi_old ;Listing 13.1 Assembly of the transient fluxes resulting from the implicit Euler scheme13.5.2OpenFOAM®In OpenFOAM®, explicit and implicit time derivatives [21] are defined via thenamespaces fvm, fvc and the corresponding functions fvc::ddt(rho, phi) and fvm::ddt(rho, phi), respectively.
Moreover the first and second order upwind Eulerschemes in addition to the second order Crank-Nicholson scheme are available,with the latter implemented following the two-step approach.The files of the transient schemes are located in the directory “$FOAM_SRC/finiteVolume/finiteVolume/ddtSchemes”. A base class denoted by ddtScheme<Type> is defined from which all time discretization schemes have to be derived.The first order Euler scheme is implemented in the class EulerDdtScheme.
Theclass is declared on top of the base class ddtScheme <Type>, as shown in Listing13.2.template<class Type>class EulerDdtScheme:public ddtScheme<Type>Listing 13.2 Decleration of the EulerDdtScheme classThe implementation of the associated fvc and fvm namespaces are defined in thefile EulerDdtScheme.C. The implicit evaluation of the Euler scheme is defined viathe following function in Listing 13.3:template<class Type>tmp<fvMatrix<Type> >EulerDdtScheme<Type>::fvmDdt(const volScalarField& rho,const GeometricField<Type, fvPatchField, volMesh>& vf)Listing 13.3 Definition of the function needed for implicit solution using the Euler scheme13.5Computational Pointers527As depicted in Listing 13.4, the first step in this function is the definition of thefvMatrix in which only the diagonal coefficient vector is filled.{tmp<fvMatrix<Type> > tfvm(new fvMatrix<Type>(vf,vf.dimensions()*dimVol/dimTime));fvMatrix<Type>& fvm = tfvm();Listing 13.4 Definition of the fvMatrixAfter allocating the needed space for storing the diagonal coefficients andsources, the values to be stored are defined and computed.
As shown in Listing 13.5,this is accomplished by first defining the reciprocal of the time step as rDeltaT, thencalculating at ¼ qC VC =Dt and storing its value in the fvm.diag() vector. The sourcecontribution is computed as the product of aot ¼ qoC VC =Dt and the old value of vfand stored in the fvm.source() vector, where vf is the generic variable used whileapplying the time scheme.scalar rDeltaT = 1.0/mesh().time().deltaTValue();fvm.diag() = rDeltaT*rho*mesh().V();fvm.source()=rDeltaT*rho.oldTime()*vf.oldTime().internalField()*mesh().V();return tfvm;}Listing 13.5 Calculation of the terms added to the diagonal and source vectorsThe script in Listing 13.6 shows that OpenFOAM®allows explicit evaluation ofthe unsteady term using the Euler scheme.
In this case, a GeometricField objectcontaining the value of ðq/ qo /o Þ=Dt is returned (here the value is per unitvolume).tmp<GeometricField<Type, fvPatchField, volMesh> >EulerDdtScheme<Type>::fvcDdt(const GeometricField<Type, fvPatchField, volMesh>& vf)return tmp<GeometricField<Type, fvPatchField, volMesh> >(new GeometricField<Type, fvPatchField, volMesh>(ddtIOobject,rDeltaT*(vf*rho - vf.oldTime()*rho.oldTime())));Listing 13.6 Explicit calculation of the unsteady term using the Euler scheme52813 Temporal Discretization: The Transient TermThe SOUE scheme is implemented in OpenFOAM® under the classbackwardDdtScheme. The definition of the class is on top of the base class, asshown in Listing 13.7.template<class Type>class backwardDdtScheme:public fv::ddtScheme<Type>{// Private Member Functions//- Return the current time-stepscalar deltaT_() const;//- Return the previous time-stepscalar deltaT0_() const;Listing 13.7 Script used to define the backwardDdtScheme class for the implementation of theSOUE schemeIn this case OpenFOAM® uses information from the current and the previoustime steps.
In the general case, the time steps are different necessitating the use oftwo variables to store their values.The implicit time discretization is defined in a way similar to the first ordertransient scheme using the function shown in Listing 13.8 through Listing 13.10.template<class Type>tmp<fvMatrix<Type> >backwardDdtScheme<Type>::fvmDdt(const volScalarField& rho,const GeometricField<Type, fvPatchField, volMesh>& vf)Listing 13.8 Script used to define the implicit discretization using the SOUE schemeThe part of the function displayed in Listing 13.9 calculates the transient coefficients that multiply the current, old, and old-old values of the dependent variableand stores the contribution to the diagonal coefficients in fvm.diag() vector.scalar rDeltaT = 1.0/deltaT_();scalar deltaT = deltaT_();scalar deltaT0 = deltaT0_(vf);scalar coefft= 1 + deltaT/(deltaT + deltaT0);scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));scalar coefft0 = coefft + coefft00;fvm.diag() = (coefft*rDeltaT)*rho.internalField()*mesh().V();Listing 13.9 Script used to calculate the unsteady coefficients and to store the contribution to thediagonal coefficients in the fvm.diag() vector13.5Computational Pointers529In the last part of the function shown in Listing 13.10, the contribution of theunsteady term is computed and stored in the fvm.source() vector.fvm.source() = rDeltaT*mesh().V()*(coefft0*rho.oldTime().internalField()*vf.oldTime().internalField()- coefft00*rho.oldTime().oldTime().internalField()*vf.oldTime().oldTime().internalField());Listing 13.10 Script used to calculate and store contribution to the source in the fvm.source()vectorBy comparing the coefficients of the SOUE scheme with the ones given inEq.
(13.100) it is easily seen that OpenFOAM® adopts a finite difference approachfor the discretization of the unsteady term whereby the time derivative is approximated via Taylor series expansions.13.6ClosureThe chapter covered the discretization of the transient term in the unsteady conservation equation.
For that purpose, two general methodologies were discussed.One method is based on a finite difference discretization while the other follows afinite volume approach in which the conservation equation is integrated over atemporal element. The first order fully implicit and fully explicit transient schemeswere presented. The formulation of higher order approximations was also investigated. This included the CN and the SOUE schemes for uniform and non uniformtime steps.
The next chapter is devoted to the discretization of the source term,relaxation of the algebraic system of equations, and other related details.13.7ExercisesExercise 1Transient heat transfer for the one dimensional body shown in Fig. 13.21 is governed by the following energy equation:@ qcp T@@T¼k@t@x@x53013 Temporal Discretization: The Transient TermThe body is insulated at one end while subjected to convective heat transfer at thesecond end.
Other parameters include TR = 330 K, hR = 400 W/m2K, k = 55 W/mK,q ¼ 7000 kg=m3 , and cp ¼ 400 J=Kg K.xTb1T1xTb2T2hRTRInsulatedFig. 13.21 One-dimensional domain used for Exercise 1(a) Compute the temperature field using the Euler Explicit method for three timesteps. Note that the initial temperature is Ti ¼ 273 K with Dt ¼ 20 s andDx ¼ 0:015 m.(b) Repeat part a using an implicit Euler scheme.(c) Explain the difference in temperatures between the two methods at timet = 60 s.Exercise 2The body described in Exercise 1 is now insulated at one end while subjected to aDirichlet boundary condition at the second end. The initial and boundary conditionsare Ti ¼ 273 K, Tb2 ¼ 330 K while values of other parameters are given byDx ¼ 0:015 m; k ¼ 55 W=mK; q ¼ 7000 kg=m3 and cp ¼ 400 J=Kg K:xTb1T1xT2Tb2InsulatedFig.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
