Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 88
Текст из файла (страница 88)
(2.85) and computed asZZXrp dV ¼p dS ¼pf S fð15:68Þ@VCVCf nbðCÞThe body force term is integrated directly over the control volume to yieldZf b dV ¼ ðf b ÞC VCð15:69ÞVCUsing a first order Euler scheme for the discretization of the unsteady term, a HRscheme for the convection term implemented via the deferred correction approach,and decomposing the diffusion flux into an implicit part aligned with the grid and anexplicit cross diffusion part, the discretized momentum equation is written in vectorform asXavC vC þavF vF ¼ bvCð15:70ÞFNBðCÞ59015 Fluid Flow Computation: Incompressible Flowswhere the coefficients are given byXavC ¼ FluxCC þðFluxCf Þf nbðCÞavFbvC¼ FluxFf¼ FluxVC XXFluxVf þf nbðCÞð15:71Þlf ðrvÞTf Sf ðrpÞC VCf nbðCÞwith the face fluxes calculated usingEfFluxCf ¼ m_ f ; 0 þ lfdCF|fflfflfflffl{zfflfflfflffl}|fflffl{zfflffl}convectioncontributiondiffusioncontributionEfFluxFf ¼ m_ f ; 0 lfdCF|fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}|fflffl{zfflffl}convectioncontributionFluxVf ¼ lf ðrvÞf diffusioncontributionTf þ m_ f ðvHRfð15:72Þ vUf Þand the element fluxes computed fromqC VCDtq VCFluxVC ¼ C vC ðf b ÞC VC|fflfflfflfflfflfflfflDt{zfflfflfflfflfflfflffl} |fflfflfflffl{zfflfflfflffl}FluxCC ¼transientcontributionð15:73Þsource termcontributionEven though the algebraic form of the momentum equation, Eq.
(15.70), islinear, its coefficients depend on the velocity and pressure fields. This nonlinearityis handled by an iterative process during which the coefficients are calculated at thestart of every iteration based on values of the dependent variables obtained in theprevious iteration. This change in the values of the coefficients results in largechanges in v and affects the rate of convergence to the degree of even causingdivergence. To slow down the changes, under-relaxation can be applied when thetransient time steps used are large.
Denoting the under relaxation factor by λv andadopting Patankar’s implicit relaxation approach, the under relaxed momentumequation can be written asXavC1 kv v ðnÞvvvþav¼bþaC vCCFFCkvkvFNBðCÞð15:74Þ15.5General Derivation591By redefining avC and bvC such thatavCavCkvbvCbvC1 kv v ðnÞþaC vCkvð15:75Þthe under-relaxed momentum equation can be rewritten asavC vC þXavF vF ¼ bvCð15:76ÞFNBðCÞFor the derivation of the collocated pressure correction equation, the pressuregradient is taken out of the bvC source term and displayed explicitly to yield^bvC ¼ VC ðrpÞC þ bCvð15:77ÞSubstituting back in Eq.
(15.76), the momentum equation becomesvC þX av^vbVCFCv¼ðrpÞþFCvvv :aaaCCCFNBðCÞð15:78ÞDefining the following vector operators:HC ½v ¼X avFv vFaCf NBðCÞ^vbCavCVCDvC ¼ vaCBvC ¼ð15:79ÞEquation (15.78) is reformulated asvC þ HC ½v ¼ DvC ðrpÞC þ BvC ;a form that will be useful in later derivations.ð15:80Þ59215 Fluid Flow Computation: Incompressible Flows15.5.2 The Collocated Pressure Correction EquationAs in the case of a staggered grid, starting with guessed values or values obtainedfrom the previous iteration vðnÞ ; m_ ðnÞ ; pðnÞ , the momentum equation, Eq.
(15.80),is first solved to obtain a momentum conserving velocity field v*. Thus the obtainedsolution satisfiesvC þ HC ½v ¼ DvC rpðnÞ þ BvCCð15:81Þwhile the final solution should satisfy Eq. (15.80). The difference between these twoequations is that the velocity field in Eq. (15.80) satisfies both the momentum andcontinuity equations while the one in Eq. (15.79) does not necessarily satisfy thecontinuity equation because of the linearization in which pressure and velocity arebased on the previous iteration values. Therefore corrections to the velocity, massflow rate, and pressure fields are needed to enforce mass conservation. Denotingthese corrections by ðv0 ; p0 ; m_ 0 Þ the relations between the exact and computed fieldscan be written asv ¼ v þ v 0p ¼ pðnÞ þ p0m¼m þmð15:82Þ0Substituting the mass flow rate given by Eq.
(15.82) into Eq. (15.25), the continuity equation becomesXXð15:83Þm_ 0f ¼ m_ f where m_ f ¼ qf vf Sff nbðCÞf nbðCÞwith the face velocity computed using the Rhie-Chow interpolation asðnÞðnÞvf ¼ vf Dvf rpf rpf :ð15:84ÞWhen the computed mass flow rate field is conservative, the RHS of Eq. (15.83)is zero yielding a zero correction field. On the other hand, an incorrect velocity fieldleads to an imbalance in mass and a nonzero value of the RHS of Eq. (15.83)implying the need for a correction field for conservation to be enforced.Mass flow rate corrections can be written in terms of velocity corrections, whichcan be derived by subtracting Eq.
(15.81) from Eq. (15.80) to yieldv0 C þ HC ½v0 ¼ DvC ðrp0 ÞCð15:85Þ15.5General Derivation593A similar equation holds for element F and is given byv0 F þ HF ½v0 ¼ DvF ðrp0 ÞFð15:86ÞThe mass flow rate correction at a cell face can be expressed asm_ 0f ¼ qf v0 f Sfð15:87Þwhere the face velocity correction is obtained by subtracting Eq. (15.84) fromEq. (15.60) to givev0 f ¼ v0 f Dvf ðrp0f rp0f Þð15:88ÞSubstitution of Eqs.
(15.87) and (15.88) in Eq. (15.83), leads to the followingform of the pressure correction equation:X X X qf v0 f Sf þqf Dvf rp0f Sf qf Dvf ðrp0 Þf Sff nbðCÞf nbðCÞf nbðCÞ¼Xm_ fð15:89Þf nbðCÞIn this equation the underlined part represents the effects of the neighboringvelocity corrections on the velocity correction of the element under consideration.This influence becomes clearer by interpolating Eqs. (15.85) and (15.86) to the faceyielding the following equivalent expression for the underline terms:v0 f þ Hf ½v0 ¼ Dvf rp0 f ) v0 f þ Dvf rp0 f ¼ Hf ½v0 :ð15:90ÞSubstituting Eq.
(15.90) in Eq. (15.89), the pressure correction equation isrewritten asX qf Dvf ðrp0 Þf Sff nbðCÞX¼f nbðCÞX m_ f þqf Hf ½v0 Sfð15:91Þf nbðCÞor more explicitly in the formX f nbðCÞqf Dvf ðrp0 Þf Sf11X avF 0 A@qf @ Sf A :m_ f þ¼v vFaCf nbðCÞf nbðCÞFNBðCÞXX0 0ð15:92ÞIn Eq.
(15.91) or (15.92) the treatment of the underlined term is critical torendering the equation solvable. In the original SIMPLE algorithm it is neglected,59415 Fluid Flow Computation: Incompressible Flowsthus linking the velocity correction at a point directly to pressure corrections.Because this is a correction equation the modification or dropping of the term willnot affect the final solution, since at convergence the corrections become zero.However it will affect the convergence rate in that the larger is the neglected termthe higher will be the error present in the approximation at each iteration.The remaining terms in Eq.
(15.91) or (15.92) can be easily treated. The coefficients of the pressure correction equation are obtained as per the discretization ofthe diffusion term in Chap. 8, specifically the treatment of anisotropic diffusion.Thus the term on the LHS is modified into a gradient dot product of the formTDvf ðrp0 Þf Sf ¼ ðrp0 Þf Dvf Sf T¼ ðrp0 Þf Dvf Sfð15:93Þ¼ ðrp0 Þf S0 fThe expanded expression of S0 f is given by2S0 f ¼ DvfTDuf6 0 Sf ¼ 400Dvf032 3 230Duf SxfSxf6yy 70 754 Sf 5 ¼ 4 Dvf Sf 5zSfDwfDwf Szfð15:94ÞWorking with S0 f , the discretization of the pressure correction gradient termproceeds as usual resulting inðrp0 Þf S0 f ¼ ðrp0 Þf Ef þ ðrp0 Þf TfEf 0¼pF p0C þ ðrp0 Þf TfdCFð15:95Þwhere the following decomposition of S0 f was used:S0 f ¼ Ef þ Tf :ð15:96ÞThe type of decomposition could be any of those reviewed in Chap.
8, as will bedetailed later. The underlined term, arising due to grid non-orthogonality, can eitherbe neglected or retained. If neglected, it will not affect the final solution as it is acorrection term. If retained, then it will be treated explicitly with an internal loop(non-orthogonal loop in OpenFOAM®). As the solution starts with a zero pressurecorrection field at every iteration, the term has to be updated iteratively whilesolving the equation.Dropping the non-orthogonal contribution, the linearized term of the pressurecorrection equation becomes15.5General Derivation595Ef 0pF p0CdCF¼ Df p0F p0Cðrp0 Þf S0 f ¼ð15:97ÞSubstituting back in Eq.
(15.91) the algebraic form of the pressure correctionequation is obtained as0apC p0C þX00apF p0F ¼ bCpð15:98ÞFNBðCÞwith the coefficients given by0apF ¼ FluxFf ¼ qf DfXX p00apC ¼ FluxFf ¼ aFf nbðCÞ0bpC¼XFNBðCÞFluxVf þf nbðCÞ¼Xf nbðCÞX qf Hf ½v0 Sff nbðCÞm_ f þX qf Hf ½v0 Sfð15:99Þf nbðCÞNote that different approximations to the underlined terms in Eq. (15.99) resultin different algorithms.
In the original SIMPLE algorithm these terms are simplyneglected.Finally the mass flow rate m_ f in Eq. (15.99), is the one computed after solvingthe momentum equation using as usual the Rhie-Chow interpolation technique withthe latest velocity field, i.e.,ðnÞðnÞm_ f ¼ qvf Sf ¼ qvf Sf Dvf rpf rpf Sf :ð15:100ÞFollowing the calculation of the pressure correction field, the pressure andvelocity at the element centroids and the mass flow rate at the element faces are allcorrected.
As mentioned above, the underlined term in Eq. (15.99) is neglected inthe SIMPLE algorithm resulting in large pressure correction values that may slowthe rate of convergence or cause divergence. To increase robustness and improvethe convergence behavior, pressure correction values obtained from Eq. (15.98) areexplicitly under relaxed. No under relaxation is used when updating the velocityand mass flow rate fields since the pressure correction will ensure mass conservation for these fields. Denoting the under relaxation factor by λp, the followingcorrection equations are used:59615 Fluid Flow Computation: Incompressible Flows0vC ¼ vC þ v C_ f þ m_ 0fm_ f ¼mpC¼ðnÞpCþv0 C ¼ DvC ðrp0 ÞCm_ 0f ¼ qf Dvf rp0f Sfð15:101Þkp p0C15.5.3 Calculation of the Df TermThe type of decomposition suggested in Eq. (15.96) could be any of those reviewedin Chap.
8 with different approaches leading to different expressions for Df asderived below.15.5.3.1 Minimum Correction ApproachFor this approach Ef is obtained by substituting S0 f for Sf in Eq. (8.68) leading toEf ¼ eCF S0 f eCFð15:102Þwhere eCF is a unit vector in the CF direction. Combining Eqs. (15.94), (15.102),and (8.64), Ef becomesEf ¼yzxdCFDuf Sxf þ dCFDvf Syf þ dCFDwf Szf2dCFdCFð15:103ÞUsing Eq. (15.103), the following expression for Df is derived:Df ¼yzxdCFDuf Sxf þ dCFDvf Syf þ dCFDwf SzfEf¼ 2 y 2 2dCFdxþ dþ dzCFCFð15:104ÞCF15.5.3.2 Orthogonal Correction ApproachThe definition of Ef in this case is obtained from Eq.
(8.69) and written asEf ¼ S0f eCFð15:105ÞCombining Eqs. (15.105), (15.94), and (8.64), Df is found to beDf ¼EfdCFvffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2uu u x 2yu Df Sf þ Dvf Sf þ Dwf Szf¼ t 2 y 2 2zxdCFþ dCF þ dCFð15:106Þ15.5General Derivation59715.5.3.3 Over-Relaxed ApproachIn this method, Ef is computed from Eqs. (8.64) and (8.70) asEf ¼S0 f S0 fdCFdCF S0 fð15:107ÞCombining Eq.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
