Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 94
Текст из файла (страница 94)
(15.190) in Eq. (15.186), an expression relating λv and λp isobtained askp 1 kvð15:191ÞExperience has shown that the performance of the SIMPLE algorithm with itsunder relaxation factors satisfying Eq. (15.191) is similar to that of the SIMPLECalgorithm.63015 Fluid Flow Computation: Incompressible Flows15.9 Treatment of Various Terms with the Rhie-ChowInterpolation15.9.1 Treatment of the Under-Relaxation TermThe use of the collocated grid method with the Rhie-Chow interpolation resulted insolutions that are dependent on the value of under relaxation factor in the momentumequation. To eliminate this dependence, a modification to the Rhie-Chow interpolation is required.
The under relaxed momentum equation is written asX1 v1 kv v ðnÞvvav¼avþbVrpþaC vCCCCF FCkv CkvFNBðCÞð15:192Þwhere bvC is the source term of the momentum equation from which the pressureðnÞand under relaxation source terms are extracted and vC is the previous iterationvalue of velocity at cell centroid C. The corresponding under relaxed momentumequation using a staggered grid formulation can be expressed asX1 v1 kv v ðnÞvvav¼avþbVrpþaf vffnbffnbfkv fkvnbNBðf Þð15:193ÞThe Rhie-Chow interpolation method mimics the staggered grid formulation byforming a pseudo-momentum equation at the cell face. It is because of thisbehavioral imitation that the Rhie-Chow interpolation is successful. Therefore as aguiding principle, the yardstick to any modification to the Rhie-Chow interpolationshould be whether the modified formulation is similar to the staggered grid formulation.
Therefore, the form of the under relaxed equation using the Rhie-Chowinterpolation should be given byX1 v1 kv v ðnÞvva vf ¼ anb vnb þ bf Vf rpf þaf vfkv fkvnbNBðf Þð15:194ÞThe average of the first term on the right hand side is obtained asXnb NBðf Þ0avnb vnb þ bvf ¼ gC @XF NBðCÞ10avF vF þ bvC A gF @X1avN vN þ bvF AN NBðFÞ1 v1 kv v ðnÞ¼ gC v aC vC þ VC rpC aC vCkvk11 kv v ðnÞavþ gF v avF vF þ VF rpF F Fkvkv11kv ðnÞaf vf¼ v avf vf þ Vf rpf kvkð15:195Þ15.9Treatment of Various Terms with the Rhie-Chow Interpolation631Substituting Eq. (15.195) into Eq. (15.194), the extended Rhie-Chow interpolated cell face velocity vf is obtained asðnÞðnÞvf ¼ vf Dvf rpf rpf þ ð1 kv Þ vf vfð15:196ÞNot accounting for the effect of under-relaxation on the face velocity results insolutions that depend on the under relaxation factor.15.9.2 Treatment of the Transient TermWhen solving a transient problem with a backward Euler transient scheme thediscretized momentum equation can be written asX avC vC ¼ avF vF þ bvC VC rpC þ aC vCð15:197ÞFNBðCÞwhere bvC is the source term of the momentum equation from which the pressureand transient source terms are extracted.
The equivalent equation for the staggeredgrid variable arrangement has a similar form given byavf vf ¼ X avnb vnb þ bvf Vf rpf þ af vfð15:198ÞnbNBðf ÞUsing the Rhie-Chow interpolation method, a pseudo cell-face equation will beconstructed asXavf vf ¼ avnb vnb þ bvf Vf rpf þ af vfð15:199ÞnbNBðf ÞThe average of the first term on the right hand side is obtained as01XX avnb vnb þ bvf ¼ gC @avF vF þ bvC AnbNBðf ÞFNBðCÞ0 gF @X 1avN vN þ bvF ANNBðFÞ¼ gC avC vC þ VC rpC aC vCþ gF avF vF þ VF rpF aF vF¼ avf vf þ Vf rpf af vfð15:200Þ63215 Fluid Flow Computation: Incompressible FlowsSubstituting into Eq. (15.84), the extended Rhie-Chow interpolated cell facevelocity vf is obtained af Dvf vf vfvf ¼ vf Dvf rpf rpf þVfð15:201ÞNot accounting for the effect of the unsteady term on the face velocity results insolutions that are time step dependent and have an oscillatory behavior for smalltime step.
This correction is valid only for the first order Euler discretization. In caseof more accurate time discretization schemes similar corrections can be performedfollowing the same principles.15.9.3 Treatment of the Body Force TermWhen treating body forces in the staggered grid arrangement, the stencil of the bodyforce term is exactly that of the pressure gradient term.
In the case of a collocatedgrid arrangement, the body force, velocity, and momentum variables are calculatedat the same location. Thus, in order to have a discretization of the body force thatretains a similar stencil as the pressure, a redistribution of the body force term isneeded. The discretized momentum equation is written asavC vC ¼ XavF vF þ bvC VC ðrpÞC þ VC BvCð15:202ÞFNBðCÞwhere the double bar indicates two averaging steps. The first step is to compute BvC(Fig. 15.28) at the cell face asBvf ¼ gC BvC þ ð1 gC ÞBvFpxstaggered gridCfFBFig. 15.28 Treatment of body force and pressure gradient on a staggered gridð15:203Þ15.9Treatment of Various Terms with the Rhie-Chow InterpolationFig. 15.29 Standard andimproved Rhie-Chowtreatment of body forces633pxstandard treatmentCBBCimproved treatmentpxwhile the second (Fig. 15.29) is to get an average of these face values at the cellcentre.The average values at cell center can best be derived [26] by considering the onedimensional situation depicted in Fig.
15.30.For the case of a stationary fluid, the pressure gradient should be in equilibriumwith the body forces leading to0 ¼ rpf þ Bvfð15:204ÞExpanding the above equation, a relation between the pressures at C and F canbe written aspC ¼ pF þ Bf dyFyfd CFCFig. 15.30 One dimensional stationary fluidBfð15:205Þ63415 Fluid Flow Computation: Incompressible Flowsor more generally aspC ¼ pF þ Bf dCFð15:206Þwhere Bf is the magnitude of Bf given byBf ¼ qf gð15:207ÞFor incompressible flows, the variation with temperature of the densityappearing in the body force term is modeled using the Boussinesq approximation asgiven by Eq. (3.101).Again for cell C the pressure gradient should be in equilibrium with the bodyforces, resulting in0 ¼ rpC þ BvC ) rpC ¼ BvCð15:208ÞHowever the pressure gradient for cell C is computed asPpf SfrpC ¼fVP CðgC pC þ ð1 gC ÞpF ÞSfð15:209Þf¼VCsubstituting from Eq.
(15.206) givesPrpC ¼fP¼gC pF þ Bvf dCF þ ð1 gC ÞpF SfpF S ffVCP¼ pFPfþfþfVCVCgC Bvf dCF SfPSfgC Bvf df Sfð15:210ÞVCVC|{z}0P vgC Bf df Sf¼f¼BvCVCwhich implies thatPBvC ¼fgC Bvf df SfVCð15:211Þ15.9Treatment of Various Terms with the Rhie-Chow Interpolation635The second requirement is that the cell-face velocity be similar to that of thestaggered arrangement equation:avf vf ¼ Xavnb vnb þ bvf Vf rpf þ Vf Bvfð15:212ÞnbNBðf Þwhere bvC is the source term given in Eq. (15.71) from which the pressure and bodyforce terms are extracted.
The averaging of the coefficients yieldsX0X avnb vnb þ bvf ¼ gC @nbNBðf ÞavF vF þ bvC AFNBðCÞ0 gF @1X 1avN vN þ bvF ANNBðFÞhi¼ gC avC vC þ VC rpC VC BvChiþ gF avF vF þ VF rpF VF BvFð15:213Þ¼ avf vf þ Vf rpf Vf Bvfand substituting into Eq. (15.179), the extended Rhie-Chow interpolated cell facevelocity vf is obtained asvf ¼ vf Dvf rpf rpf þ Dvf Bvf Bvfð15:214Þwhere Bvf is calculated asBvf ¼ gC BvC þ ð1 gC ÞBvFð15:215ÞThe above additional treatment of the cell face velocity increases the overallrobustness of the solution procedure for situations where variations in body forcesare important (e.g., free-surface flows).63615 Fluid Flow Computation: Incompressible Flows15.9.4 Combined Treatment of Under-Relaxation, Transient,and Body Force TermsIn general all three terms described above should be dealt with together.
Thisnecessitates modifying the Rhie-Chow interpolation to account for all three effects.Fortunately this can easily be derived by using the principle of superpositionleading to the following interface velocity:vf ¼ vf Dvf rpf rpf þ Dvf Bvf Bvfþaf Dvf Vfvfvfð15:216ÞðnÞðnÞþ ð 1 k Þ vf vfvwhere in calculating Dvf the under relaxed value of the avC coefficient is used.15.10 Computational Pointers15.10.1 uFVMIn uFVM, the pressure correction equation is implemented in one script file denotedby cfdAssembleMdotTerm. Listing 15.1 shows the core of the algorithm wherebythe coefficients of the pressure correction equation are assembled by linearizing thefluxes at each of the interior faces.
In addition, the mdot field (i.e., the mass flowrate at the faces) is calculated based on Eq. (15.216), which is subdivided into 9terms (i.e., terms I trough IX) and assembled step by step to produce the cell facevelocity. The terms into which the velocity at the face is decomposed are asfollows:term I: the interpolated velocity field vf ,terms II and III: the face and average pressure gradients Dvf rpf rpf ,terms IV and V: the average and redistributed body forces Dvf Bvf Bvf ,terms VI and VII: the transient fluxesaf Dvf Vfvf vf , andðnÞðnÞterms VIII and IX: the relaxation correction term ð1 kV Þðvf vf Þ.15.10Computational Pointers637%% assemble term I%density_f [v]_f.Sf%U_bar_f = (dot(vel_bar_f(:,:)',Sf(:,:)'))';local_FLUXVf =local_FLUXVf+ density_f.*U_bar_f;%% Assemble term II and linearize it%- density_f ([DPVOL]_f.P_grad_f).Sf%DUSf = [DU1_f.*Sf(:,1),DU2_f.*Sf(:,2),DU3_f.*Sf(:,3)];geoDiff =( dot(Sf(:,:)',DUSf') ./ dot(CN(:,:)',Sf(:,:)'))';local_FLUXCf1 = local_FLUXCf1 + density_f.*geoDiff;local_FLUXCf2 = local_FLUXCf2 - density_f.*geoDiff;local_FLUXVf = local_FLUXVf - density_f.*dot(p_grad_f',T')';%% assemble term III%density_f ([P_grad]_f.([DPVOL]_f.Sf))%local_FLUXVf = local_FLUXVf +density_f.*dot(p_grad_bar_f(iFaces,:)',DUSf(iFaces,:)')';%% assemble terms IV and V%density_f [DBVOL]_f.([B]_f -[[B]]_f).S_f%local_FLUXVf = local_FLUXVf +density_f[iFace]*FVVectorDotProduct(FVTensorVectorDotProduct(DB_f,FVVectorSubstract(FVMakeVector(bf1_bar_f[iFace],bf2_bar_f[iFace]),FVMakeVector(bf1_redistributed_f[iFace],bf2_redistributed_f[iFace]))),S_f) ;%% assemble terms VI and VII%[Dt]_f (U_Old_f -[v_old]_f.S_f)%U_bar_old_f = [velx_old_bar_f[iFace] vely_old_bar_f[iFace])] *S_f'local_FLUXVf = local_FLUXVf + DT_f*(mdot_old_f[iFace] density_old_f[iFace]*U_bar_old_f);%% assemble terms VIII and IX%(1-URF)(U_f -[v]_f.S_f)%local_FLUXVf = local_FLUXVf + (1.0-Mdot_URF)*(mdot_previous_f density_f.*U_bar_f);%% assemble the flow term dot for the face%local_mdot_f = local_FLUXCf_1*(pressure[iElement1]+ Pref) +local_FLUXCf_2*(pressure[iElement2]+ Pref) + local_FLUXVf;%%% Assemble in Global Fluxes%Listing 15.1 Script used for the calculation of the mass flow rates and coefficients of the pressurecorrection equation in uFVM63815 Fluid Flow Computation: Incompressible FlowstheFluxes.FLUXC1f(iFaces,1) = local_FLUXCf1;theFluxes.FLUXC2f(iFaces,1) = local_FLUXCf2;theFluxes.FLUXVf(iFaces,1) = local_FLUXVf;%theFluxes.FLUXTf(iFaces,1) = theFluxes.FLUXC1f.*pressureC(iFaces)+ theFluxes.FLUXC2f.*pressureN(iFaces) + theFluxes.FLUXVf(iFaces);%mdot_f = theFluxes.FLUXTf(iFaces);Listing 15.1 (continued)15.10.2 OpenFOAM®The numerical techniques introduced so far are used in what follows to developOpenFOAM® [27] applications for solving the incompressible Navier-Stokesequations.15.10.2.1 Pressure Correction SIMPLE SolversBased on the SIMPLE algorithm, a number of solvers will be constructed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
