Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 98
Текст из файла (страница 98)
It clearly shows that pressure acts on both the density fieldthrough the equation of state and the velocity field via the gradient in themomentum equation to enforce mass conservation. This dual role explains thesuccess of the pressure-based approach to predict fluid flow at all speeds. This facthowever did not deter workers in the density based track from using the artificialcompressibility technique [19] to develop methods capable of solving fluid flow atall speeds. To overcome degradation in performance due to the stiff matricesencountered in these methods, preconditioning of the resulting stiff matrices wasintroduced and several methods [20, 21] using this technique have recentlyappeared in the literature.16.2Introduction657Similarly several pressure-based methods [22] for predicting fluid flow at allspeeds following either a staggered grid approach [23] or a collocated variableformulation have been developed.
While in some methods primitive variables wereused, others employed the momentum components as dependent variables. Someworkers adopted the stream-wise directed density-retardation concept, which iscontrolled by Mach-number-dependent monitor functions [24, 25], to account forthe hyperbolic character of the conservation laws in the transonic and supersonicregimes. Other techniques used the first order upwind scheme for evaluating thedensity at the control volume faces at high Mach number values and the centraldifference scheme at low values [26].This chapter extends the collocated pressure-based technique developed in theprevious chapter allowing the simulation of fluid flows at all Mach number values.The adopted method is easy to implement, highly accurate, and does not require anyexplicit addition of damping terms to improve robustness or to properly resolveshock waves.16.3 The Conservation EquationsThe conservation equations for solving compressible flow problems include thecontinuity, momentum, and energy equations.
For a Newtonian fluid behaving as anideal gas, these equations can be expressed as@qþ r ðqvÞ ¼ 0@tð16:1Þno 2@½qv þ r fqvvg ¼ r flrvg rp þ r lðrvÞT rðlr vÞ þ f b@t3ð16:2ÞDcp Dp 2@qcp T þ r qcp vT ¼ r ½krT þ qT lW þ lU þ q_ Vþ@tDt 3Dtð16:3Þwhere the form of the energy equation adopted here is the one expressed in terms oftemperature. The above set of equations should be appended by an equation of staterelating density to pressure and temperature, i.e., q ¼ qðp; T Þ, which for an idealgas is given byq¼where R is the gas constant.pRTð16:4Þ65816Fluid Flow Computation: Compressible FlowsIn the derivations to follow, superscript n refers to values used at the beginningof an iteration, superscript indicates values updated once during an iteration, andsuperscript refers to values updated twice during the same iteration.16.4 Discretization of the Momentum EquationThe discretized momentum equation, Eq.
(16.2), over the control volume C shownin Fig. 16.1 is similar to its incompressible form given in Chap. 15. The only twodifferences are related to the interpolation of density to the interface and theadditional term ð2=3Þrðlr vÞ involving the bulk viscosity. Starting with thefirst difference, the density in compressible flows is no longer constant and since itis stored at the control volume centroids it has to be interpolated to find its value atthe control volume faces where it is needed for computing the mass flow rate.
Theuse of a linear interpolation profile (central difference) causes oscillation at highspeeds. Thus a bounded upwind biased scheme should be used. Any of the boundedconvective schemes presented in Chaps. 11 and 12 can be adopted for that purpose.The second difference is the additional term involving rðlr vÞ. This term hasnot been discretized so far and its discretized form is obtained by making use ofEq.
(2.85) based on which the volume integral of the gradient of a scalar quantity isF1F2Convectionf2Diffusionf1TransientCf6f3F3Source/Sinkf4f5F4F5Fig. 16.1 A schematic of a control volume C with its neighborsF616.4Discretization of the Momentum Equation659transformed into a surface integral and then into a summation of fluxes over thefaces of the control volume according toZZX½rðlr vÞdV ¼ðlr vÞdS ¼ðlr vÞf Sfð16:5ÞVCf nbðC Þ@VCThe divergence of the velocity vector at the faces is computed as" #@u ðnÞ@v ðnÞ@w ðnÞðnÞþþðlr vÞf ¼ lf@x f@y f@z fð16:6Þwhere the gradient of / ¼ u; v or w is interpolated linearly to the face ðnÞ ðnÞ ðnÞ@/@/@/¼ gCþ gF@x f@x C@x Fð16:7ÞThe final discretized form of the momentum equation is givenPby Eq.
(15.70)ðlr vÞf Sfwith its coefficients given by Eq. (15.71) with the term ð2=3Þf nbðC Þadded to the source term in that equation.As for incompressible flow problems, the algebraic equations are under relaxedand written in the form of Eq. (15.78), which is suitable for the derivation of thepressure correction equation.16.5 The Pressure Correction EquationThe pressure correction equation for compressible flows is obtained by a simpleextension of that for incompressible flows. The difference is related to variations indensity which are accounted for by defining a density correction field q0 andrelating it to the pressure-correction field through a pressure-density relation.
This,however, yields substantial differences in the treatment of boundary conditions, aswill be explained later in the chapter.For an ideal gas the relation between pressure and density is written asqRT ¼ pð16:8ÞUsing this relation, an equation relating density correction to pressure correctioncan be derived by expanding Eq. (16.8) via a Taylor series to yieldqjðpðnÞ þp0 Þ ¼ qjðpðnÞ Þ þ@q 0@q 01 0p ¼ q þ q0 ) q0 ¼p ¼p ¼ Cq p0 ð16:9Þ@p@pRT66016Fluid Flow Computation: Compressible FlowsThe corrected pressure, density, velocity, and mass flow rate fields are defined asp ¼ pðnÞ þ p0q ¼ q þ q0ð16:10Þv ¼ v þ v0m_ ¼ m_ þ m_ 0and the semi-discretized continuity equation can be written in terms of the correction fields as o X qC þ q0C qCVC þð16:11Þm_ f þ m_ 0f ¼ 0Dtf nbðCÞwherem_ f ¼ qf þ q0f vf þ v0f Sf¼ qf vf Sf þ qf v0f Sf þ q0f vf Sf þ q0f v0f Sf|fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}m_ fð16:12Þm_ 0fThe second order correction term q0f v0f Sf is usually neglected because it isconsiderably smaller than other terms.
This approximation does not influence theconvergence rate except during the first few iterations of the solution process. Inaddition, the final solution is not affected, since at the state of convergence, thecorrection fields vanish.Using the Rhie-Chow interpolation for vf and v0f ; m_ f and m_ 0f are respectivelyexpressed asðnÞðn Þð16:13Þm_ f ¼ qf vf Sf qf Dvf rpf rpf Sf|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}qf vf Sfandm_ 0f!m_ f¼ Sf Sf Cq; f p0f Sf þqf|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}q v0 Sfqf v0fqf Dvff¼rp0ff¼ qf Dvf rp0f Sf þqf Dvf rp0frp0f Sf þm_ fqf!q0f vf Sf Sf Cq; f p0f þ qf v0f Sf þ qf Dvf rp0f Sfð16:14Þ!m_ f Sf Cq; f p0f qf Hf ½v0 Sfqfwhere the second order term is neglected. Note the substitution of q0f with Cq; f p0f .16.5The Pressure Correction Equation661The underlined term in Eq.
(16.14) presents the same difficulties as for theincompressible algorithm and is usually dropped from the equation. Neglecting thisterm, the correction to the mass flow rate becomesm_ 0f¼qf Dvf rp0f Sf þ!m_ f Sf Cp; f p0fqfð16:15Þwhere the first term on the right hand side of Eq.
(16.15) is similar to that arising inincompressible flow while the second term is the new density correction contribution. This second term is important as it transforms the pressure correctionequation from an elliptic equation to a hyperbolic one capable of resolving shockwaves that may arise at supersonic and hypersonic speeds.
This allows the compressible SIMPLE algorithm to be used for predicting fluid flow at all speedswithout the need for any special preconditioning.More insight can be gainedthrougha simple normalization procedure whereby.Eq. (16.15) is divided by m_ f Sf Cp; f qf yielding a weighting factor of 1 for thep0f term, and a weighting factor proportional to 1 ðM 2 Þ (where M is the Machnumber of the flow) for the rp0f term, i.e., 2RT qf Dvf rp0f Sf þ p0fm_ 0f ¼ m_ f Sfð16:16ÞFor flows at low Mach number values, the rp0 correction term dominatesreturning the equation to an elliptic form as in the incompressible case. On the otherhand, for flows at very high Mach number values the p0f correction term can nolonger be neglected giving a hyperbolic character to the correction equation.