Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 99
Текст из файла (страница 99)
Thiscombined behavior allows the prediction of fluid flow at all speeds.Substitution of Eq. (16.14) in the continuity equation, Eq. (16.11), yields thecompressible form of the pressure correction equation and is written as(!)Xm_ fVC0 v00Cq pC þqf Df rpf Sf þC q pfDtqff nbðC Þ01XXqqCVC þ¼ @ CHf ½v0 Sfm_ f A þDtf nbðC Þf nbðCÞð16:17ÞAgain the treatment of the underlined term yields the different variants of theSIMPLE family of algorithms.
Dropping the underlined term, the pressure correction equation for the SIMPLE algorithm is obtained as66216Fluid Flow Computation: Compressible Flows!Xm_ f 0þCq pf qf Dvf ðrp0 Þf Sfqff nbðC Þf nbðC Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}transientlike termXVC Cq 0pDt C|fflfflfflffl{zfflfflfflffl}convectionlike termqCdiffusionlike termX qCVC ¼m_ fDtf nbðC Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}ð16:18Þsourcelike termEquation (16.18) can be obtained directly by substituting Eq.
(16.15) in Eq. (16.11).It is worth stressing that the convection-like term came about naturally duringthe derivation of the pressure correction equation and its presence is critical for theability of the algorithm to resolve flows at all speeds. Moreover, for flows at highMach number values, density correction is convected (i.e., exhibiting a hyperbolicbehavior) and the mathematical operator describing this phenomenon is the firstorder divergence operator. Thus, contrary to incompressible flows where only thediffusion-like term is present implying that the pressure correction equation exhibitsan elliptic behavior, pressure correction solutions of the form p0 þ C can no longersatisfy the equation.
This indicates that while for incompressible flows any pressurevalue can be set as a boundary condition without affecting the solution, for compressible flows it is important to define the exact value of pressure at the boundariesbecause the chosen value will affect the final solution.It should also be noted that because at convergence the correction field is zero,the order of the scheme used to discretize the convection-like term is of no consequence on the accuracy of the final results. However, this is not the case for m_ fwhere the use of high order schemes in its evaluation does improve the shockcapturing properties of the algorithm.
Thus, to enhance robustness, it is helpful touse an upwind scheme for the discretization of the convection like term. Further,neglecting the non-orthogonal contribution of the diffusion-like term as explained inChap. 15, the pressure correction equation and its coefficients become0aCp p0C þX0aFp p0F ¼ bCp0FNBðC ÞC q; faFp ¼ qF Df m_ f ; 0 qf!XXVC CqCq; f p0þaC ¼qf Dfm_ f ; 0 þDtqff nbðC Þf nbðC Þ01 XX0q qCVC þbpC ¼ @ Cqf Dvf rp0f Tfm_ f A þDtf nbðC Þf nbðC Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}0nonorthogonal term usually neglectedð16:19Þ16.5The Pressure Correction Equation663Following the calculation of the pressure-correction field p′ by solvingEq.
(16.19) the velocity, pressure, density, and mass flow rate fields are correctedusing the following equations:0vC ¼ vC þ vCv0C ¼ DvC ðrp0 ÞCð16:20ÞðnÞpC ¼ pC þ kp p0Cð16:21Þq0qC ¼ qC þ k Cq pCm_ f¼m_ fþm_ 0fm_ 0fð16:22Þ¼0vqf Df rpf Sf þ!m_ f Sf Cq; f p0fqfð16:23Þwhere kq is the under relaxation factor for density.16.6 Discretization of The Energy EquationThe discretization of the unsteady, convection, and diffusion terms of the energyequation, Eq. (16.3), follows the general procedures described in previous chaptersand will not be repeated.16.6.1 Discretization of the Extra TermsThe focus here will be on the discretization over the control volume C shown inFig.
(16.1), of the new terms appearing on the right hand side of Eq. (16.3). Theseare specific to the energy equation and have not been handled during the discretization of the general scalar equation. Many of the terms are treated as source termsand evaluated at the centroid of the element during their integration to ensure asecond order accurate discretization.16.6.1.1 The Specific Heat TermThe discretization of the term involving the specific heat proceeds as follows:664ZVC16Fluid Flow Computation: Compressible FlowsDcp ðnÞ DcpdV ¼ qC TCqTVCDtDt C" ðnÞ ðnÞ ðnÞ #cðpnÞ cp ðnÞ @cp @cp @cpþ uC¼ qC TCþ vCþ wCVCDt@x C@y C@z Cð16:24Þ16.6.1.2 The Substantial Derivative TermThe discretization of the substantial derivative of the pressure term is performed asZDpdV ¼DtVc Dpp pC@p @p @pþ uVC ¼ CþvþwVCCCCDt C@x C@y C@z CDtð16:25Þ16.6.1.3 The Dissipation TermThe discretized form of the dissipation term involving the bulk viscosity is obtained asZlWdV ¼ðnÞlC WC VCVC¼ðnÞlC@u@x 2@v@wþþVC@y C@z CCð16:26Þ16.6.1.4 The Viscous Dissipation TermThe discretization of the viscous dissipation term is performed in a way similar tothe term involving the bulk viscosity and is given byZðnÞlUdV ¼ lC UC VCVCðnÞ¼ lC8 " # 2 9>>@u 2@v 2@w 2@u @v>>>þþþ>>>=< 2 @x þ @y þ @z@y@xCC 2 2>>@u@w@v@w>>:þþþ@z C@x C@z C@y CC>>>>;VCð16:27Þ16.6Discretization of The Energy Equation66516.6.1.5 The Source/Sink TermThe term involving the heat source/sink per unit volume is discretized asZð16:28Þq_ V dV ¼ ðq_ V ÞC VCVCThe discrete forms of the above terms are substituted into the energy equation toyield the algebraic form of the energy equation as described next.16.6.2 The Algebraic Form of the Energy EquationAssuming a first order Euler scheme for the discretization of the unsteady term anda high resolution scheme for the discretization of the convection term applied in thecontext of a deferred correction approach, the final algebraic form of the energyequation can be written asaTC TC þXaTF TF ¼ bTCð16:29ÞFNBðCÞwhere the coefficients are given by Ef m_ f ; 0 cp fdCFXX c aTC ¼ aC aTF þm_ f cp f þ aCp ; 0aTF ¼ kfFNBðCÞf nbðCÞ qVCC cpqC cp C VCDcpcCaC ¼aCp ¼ qC VCaC ¼DtDtDt CXX HRTbC ¼kf ðrT Þf Tf m_ f cp f Tf TfU þ aC TCf nbðC Þf nbðCÞ c þ TC aCp ; 0 þDpDt2þ lC WC þ UC þ ðq_ V ÞC VC3Cð16:30ÞSimilar to other variables, under relaxation of the energy equation is usuallyrequired.66616Fluid Flow Computation: Compressible Flows16.7 The Compressible SIMPLE AlgorithmThe various elements of the collocated compressible SIMPLE algorithm are displayed in Fig.
16.2 and can be summarized as follows:1. To compute the solution at time t þ Dt, start with the solution at time t forpressure, velocity, density, temperature, and mass flow rate fields pðnÞ ; vðnÞ ; qðnÞ ;T ðnÞ , and m_ ðnÞ , respectively, as the initial guess.set initial guess m( n ) ,v ( n ) ,(n),and p( n )at time t + t toconverged values at timetassemble and solve momentumequation for v *Compute*usingtheequation of stateand*fm usingthe Rhie Chowinterpolationassemble and solve pressurecorrectionequation for prepeatrepeatcorrect m *f ,v * ,*,and p( n )toobtain m **f ,v ** ,**,and p*Assembleand solveenergyequation for Tset m (fn ) = m **f ,v ( n ) = v ** ,(n)=**,and p( n ) = p*noconverged ?yesset solution at time t + t tobeequal tothe converged solutionadvanceintimeset t = t + tnotimeexceeded ?yesstopFig.
16.2 A flow chart of the SIMPLE algorithm for compressible fluid flow16.7The Compressible SIMPLE Algorithm6672. Solve the momentum equation given by Eq. (16.2) to obtain a new velocityfield v .3. Use the equation of state to calculate a new density field q .4. Update the mass flow rate at the control volume faces using the Rhie-Chowinterpolation technique (Eq. 16.13) to obtain a momentum satisfying mass flowrate field m_ .5.
Using the new mass flow rates calculate the coefficients of the pressure correction equation and solve it (Eq. 16.19) to obtain a pressure correction field p0 .6. Update the pressure, density, and velocity fields at the control volume centroids and the mass flow rate at the control volume faces to obtaincontinuity-satisfying fields using Eqs. (16.20)–(16.23).7. Solve the energy equation to obtain a new temperature field T .8. Set v ; m_ ; q ; T , and p as the initial guess for velocity, mass flow rate,density, temperature, and pressure.9. Go back to step 2 and repeat until convergence.10. Set the solution at time t þ Dt to be equal to the converged solution.11. Advance to the next time step by setting the current time t to t þ Dt.12.
Go to step 1 and repeat until the last time step is reached.16.8 Boundary ConditionsGenerally, there is no difference in the implementation of boundary conditions forthe momentum equation between incompressible and compressible flow. Thereforethe required modifications in the momentum equation are those discussed inChap. 15 and will not be repeated here. However substantial differences do arisewith the pressure correction equation and will form the main subject of this section.The boundary conditions for the energy equation follow the ones described for ageneral scalar variable / and also will not be repeated here (inlet, outlet, Dirichlet,Van-Neumann, and symmetry conditions).For a boundary cell, such as the one shown in Fig.
16.3, the continuity equationis written as X qC þ q0C qCVC þm_ f þ m_ 0f þ ðm_ b þ m_ 0b Þ ¼ 0|fflfflfflfflfflffl{zfflfflfflfflfflffl}Dtf nbðCÞð16:31Þboundary facewhere the contribution of the boundary face is separately displayed with m_ b representing the boundary mass flux and m_ 0b its correction. Moreover, expressing theRhie-Chow interpolation at the boundary faces by adopting the same approach thatwas used for incompressible flow, the velocity, mass flow rate, and mass flow ratecorrection at a boundary face are expressed as66816Fluid Flow Computation: Compressible FlowsFig. 16.3 A boundarycontrol volumeboundaryelementebCbSbnboundaryfacevb|{z}boundary faceðnÞðnÞ¼ vC DC rpb rpC|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}ð16:32Þboundary RhieChowðnÞðnÞðnÞm_ b ¼ qb vC Sb qb DvC rpb rpC Sbð16:33Þ m_ bm_ 0b ¼ qb DC p0b p0C þCq; b p0bqbð16:34ÞThe only difference between these expressions and those presented in Chap.
15is in the correction equation of the mass flow rate, which has an additional termrelated to density correction. Moreover, there is no difference in the implementationof boundary conditions at a wall and a symmetry plane between incompressible andcompressible flows. Therefore the modifications to the boundary elements presented in the previous chapter for wall and symmetry boundary conditions in thepressure correction equation are applicable here with no need to be repeated.The boundary conditions left to be discussed here are those applicable at the inletand outlet.