Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 58
Текст из файла (страница 58)
The computation of the flow field will be the subject of later chapters.11.1IntroductionAlthough the convection term looks simple, its discretization presents a number ofdifficulties that have occupied researchers for more than three decades now. Theirwork has shed light on the problems that hindered its discretization and resulted inthe development of a large number of convection schemes. The body of literature isso large that two chapters are devoted for the analysis of this term.
In this chapterthe basic concepts are introduced and High Order (HO) [1–3] upwind biasedschemes are discussed. The bounding of the convective flux, which made possiblethe development of non-oscillatory convection schemes of high order of accuracy,denoted by High Resolution (HR) schemes, will be discussed in the next chapter.For clarity, the presentation of new concepts will be initially introduced using aone dimensional grid, and then extended to multi dimensional non-orthogonal grids.The chapter starts with the discretization of the one-dimensional convection-diffusion© Springer International Publishing Switzerland 2016F. Moukalled et al., The Finite Volume Method in Computational Fluid Dynamics,Fluid Mechanics and Its Applications 113, DOI 10.1007/978-3-319-16874-6_1136536611 Discretization of the Convection Termproblem to highlight, through a stability criterion, the shortcomings of the centraldifference scheme [4, 5].
The upwind scheme [6] is shown to pass the stability test.The numerical diffusion error associated with low order schemes and the numericaldispersion error accompanying high order schemes are also explained. The treatmentof the class of HR schemes will be discussed in the next chapter.11.2Steady One Dimensional Convection and DiffusionInitially the derivations are performed for a very simple, one dimensional, steadyconvection-diffusion problem to learn as much as possible without being distractedby the complexity of the calculations. The governing equation for the case studiedcan be written asd ðqu/Þ dd/C/¼0dxdxdxð11:1ÞLuckily an analytical solution for the problem is available, and is used as a reference with which various numerical solutions are compared.11.2.1 Analytical SolutionThe continuity equation for this steady-state one dimensional problem of constantcross sectional area is given byd ðquÞ¼0dxð11:2Þimplying that ρu is constant.
Having this in mind and integrating with respect to x,Eq. (11.1) becomesqu/ C/d/¼ c1dxð11:3Þwhere c1 is a constant of integration the value of which depends on the boundaryconditions used. Rearranging, Eq. (11.3) is rewritten asd/ quc1¼ // /dx CCð11:4Þ11.2Steady One Dimensional Convection and Diffusion367Through a change of variable, Eq. (11.4) is transformed todU qu¼Udx C/ð11:5ÞwhereU¼quc1/ //CCð11:6ÞSeparating variables and integrating, the solution to Eq.
(11.5) is found to bequdU ququx¼ / dx ) LnðUÞ ¼ / x þ c3 ) U ¼ c2 eC/UCCð11:7Þwhere c2 is another constant of integration. Reverting back to the original variable,the general solution for / is given byquc2 C/ eC/ x þ c1/¼quð11:8ÞThus the analytical solution between the two points W and E shown in Fig. 11.1 andsubject to/ ¼ /W at x ¼ xW/ ¼ /E at x ¼ xEð11:9Þis obtained asxxW/ /WePeL L 1¼ Pee L 1/E /Wð11:10Þwhere PeL is the Péclet number (based on the length L), which represents the ratioof the advective transport rate of / to its diffusive transport rate, and is given byPeL ¼quLC/and L ¼ xE xW :ð11:11ÞEquation (11.10) is evaluated for different values of PeL and results are displayed inFig. 11.2. As shown, variations in / between W and E change from a linear profilefor pure diffusion problems to almost a step profile at high values of PeL .36811 Discretization of the Convection TermCWWWEEEuwwuwWWueWxCECxwSw SwueewEEyCxeS e SeeCFig.
11.1 Notation for a one dimensional grid system1.0EW0.80.6PeL=0PeL=10.4PeL=2PeL=40.2PeL=100.20.40.60.8PeL=1001.0x xWLFig. 11.2 Analytical solutions of the one dimensional convection and diffusion problem for various PeL11.2.2 Numerical SolutionThe discretization of Eq.
(11.1) starts by integrating the conservation equation overthe one dimensional element shown in Fig. 11.1 to yieldZr ðqv/Þ r C/ r/ dV ¼ 0ð11:12ÞVc11.2Steady One Dimensional Convection and Diffusion369where v = ui is the velocity vector. The conservation equation can be written interms of the convection and diffusion fluxes asZr J/;C þ J/;D dV ¼ 0 where J/;C ¼ qv/ and J/;D ¼ C/ r/: ð11:13ÞVcThen, using the divergence theorem, the volume integral is transformed into asurface integral givingZr J/;C þ J/;D dV ¼Z@VcVcJ/;C þ J/;D dS ¼Z qu/i C/@Vcd/i dS ¼ 0:dxð11:14ÞReplacing the surface integral by a summation of fluxes over the element faces,Eq.
(11.14) becomesX d/i Sf ¼ 0:qu/i C/dx ff nbðC Þð11:15ÞNoting that the surface vectors on opposite sides of the element have opposite signs,the expanded form of Eq. (11.15) for a constant cross section is obtained as d/d/DyDyðquDy/Þe C/ ðquDy/Þw C/¼0dxdxewð11:16ÞThe values of u at the cell faces are known, and those of the gradient can bediscretized in the way previously described. The question is how to proceed in thediscretization of the face values /e and /w in terms of the values at adjacent nodes.The method of specifying these face values is denoted in the literature by the“advection scheme”.11.2.3 A Preliminary Derivation: The Central Difference(CD) SchemeAt first sight the “obvious” answer would be a linear interpolation profile similar tothe one used for the diffusion term.
Hence, the value of / at a given face, say theright face e, will be computed as/ð xÞ ¼ k0 þ k1 ðx xC Þð11:17Þ37011 Discretization of the Convection TermCwWeEuwwWWueuwWwCeueeEEEFig. 11.3 The Central Difference scheme profilewhere k0 and k1 are constants evaluated using the two nodes straddling face e. Thususing the fact that / ¼ /E at x ¼ xE and / ¼ /C at x ¼ xC , Eq. (11.17) evaluated atx ¼ xe gives/e ¼ /C þð/E /C Þðxe xC Þ:ð xE xC Þð11:18ÞThis is basically the central difference scheme that can be obtained by a Taylorseries expansion where terms involving derivatives of the second order and higherare neglected, which means it is second order accurate.For the uniform grid shown in Fig.
11.3, Eq. (11.18) becomes/e ¼/C þ /E2ð11:19ÞExample 1Derive the value of /e for a central difference scheme using a Taylorexpansion approachSolutionThe Taylor expansion about e can be written as@//C ¼ /ðDxC =2Þ þ ðDxC =2Þ þ O Dx2C|fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} @x e¼/e¼/e /E /C DxCdxe211.2Steady One Dimensional Convection and Diffusion371thus/e ¼ /C þDxC /E /Cdxe2for a uniform grid dxe ¼ DxC yielding/e ¼/C þ /E2Thus after the discretization of the diffusion term using a linear profile(Fig. 11.3), the term in the first square bracket of Eq. (11.16) becomesd/ðquDy/Þe CDydx/ð/E þ /C Þ/ Dy C¼ ðquDyÞeð/ /C Þ2dx e Eð11:20Þe¼ FluxCe /C þ FluxFe /E þ FluxVewhereDye ðquDyÞeþdxe2ðquDyÞe/ DyeþFluxFe ¼ Cedxe2FluxVe ¼ 0/FluxCe ¼ Ceð11:21ÞA similar expression for the term in the second square bracket of Eq. (11.16) canalso be derived and is given by d/ð/ þ /C ÞDy ðquDy/Þw C/ð/ C / W Þ¼ ðquDyÞw WDy C/dx2dx ww¼ FluxCw /C þ FluxFw /W þ FluxVwð11:22Þwhere nowDyw ðquDyÞwdxw2ðquDyÞw/ DywFluxFw ¼ Cwdxw2FluxVw ¼ 0/FluxCw ¼ Cwð11:23Þ37211 Discretization of the Convection TermSubstitution of these values into the convection-diffusion equation yieldsaC /C þ aE /E þ aW /W ¼ 0ð11:24ÞwhereDye ðquDyÞeþdxe2ðquDyÞw/ DywaW ¼ FluxFw ¼ Cwdxw2 ðquDyÞeðquDyÞw/ Dye/ DywaC ¼ FluxCe þ FluxCw ¼þ Ceþ Cwþ 2dxe2dxw/aE ¼ FluxFe ¼ Ceð11:25ÞAs the problem is one dimensional Dye ¼ Dyw and, without loss of generality, canbe set equal to 1.
Moreover, from continuity u is constant and thusðquDyÞe ðquDyÞw ¼ 0. Assuming a uniform diffusion coefficient, the coefficients ofthe discretized equation can be simplified toC/ðquÞeþx E xC2/CðquÞwaW ¼ xC xW2aC ¼ ð aE þ aW ÞaE ¼ ð11:26ÞSubstituting back in Eq. (11.24), the value for /C is found as/C /WaE¼/E /W aE þ aWð11:27ÞIf the grid is assumed to be uniform, then the above equation can be written in termsof PeL as/C /W 1PeL1¼/E /W 22ð11:28Þwhere L is ðxE xW Þ, which is the size of two elements.
The analytical solution forthe problem can be obtained from Eq. (11.10) by setting ðx xW Þ=L to 0.5 and isgiven byPeL/C /We 2 1¼ Pe/E /W e L 1ð11:29ÞThe two solutions are compared in Fig. 11.4 as PeL varies from −10 to +10.Figure 11.4 demonstrates that at low values of PeL the numerical and analytical11.2Steady One Dimensional Convection and Diffusion373Fig.
11.4 Comparison of numerical (obtained with the CD scheme) and analytical solutions forthe one dimensional convection and diffusion problemsolutions are very close to each other. However as PeL is increased beyond a certainvalue the central difference (CD) numerical solution greatly departs from the analytical solution and becomes unbounded experiencing unphysical behavior.Whereas the analytical solution approaches asymptotically the values 0 and 1 forpositive and negative values of PeL , respectively, the CD solution decreases linearlyfrom +∞ to −∞ as PeL increases from −∞ to +∞.
This solution indicates thatsome of the assumptions used in the discretization of the equation are unrealistic orunphysical. What is the reason for this behavior?As depicted in Fig. 11.5, whereas diffusion at point C is equally affected byconditions upstream and downstream of C (Fig. 11.5a), the advection process is ahighly directional process transporting properties only in the direction of the flow(Fig. 11.5b). Therefore the linear profile approximation, which assigns equal weightto both the upwind and downwind nodes, is a good approximation for the diffusionterm (Fig.
11.5a). However it cannot describe the directional preference of convection, for which a step profile is more appropriate (Fig. 11.5b), and is the cause ofthe problem.The combined convection-diffusion zone of influence and the more relevantprofile in this case is schematically depicted in Fig. 11.5c. This zone of influenceapproaches the diffusion region displayed in Fig. 11.5a and the advection regiondepicted in Fig.
11.5b at low and high values of the Péclet number, respectively.37411 Discretization of the Convection Term(a)diffusion=1diffusionWCE=0x=0(b)x =1convection=1convectionuW=0CEx=0(c)x =1convection + diffusion=1uWconvection + diffusionC=0Ex=0x =1Fig. 11.5 Zone of influence of a diffusion, b convection, and c combined convection anddiffusion terms and expected problem solution11.2Steady One Dimensional Convection and Diffusion375Therefore, as long as diffusion is the dominant transfer mechanism, the use of alinear profile yields physical results.
However, once convection overwhelms diffusion, unphysical results are obtained. The value of Péclet number at which thisoccurs can be easily calculated. Assuming the flow to be in the positive x direction,it is noted that there is a possibility for the aE coefficient to become positive, thusleading to unphysical results (if the flow is in the negative x direction, then aW maybecome positive) when,/CeðquÞe dxeDye ðquDyÞe[0 )þ[ 2:/dxe2Ceð11:30ÞDefining the cell Péclet number asPe ¼qudx:C/ð11:31Þwhich, for a uniform grid, is half PeL , then Eq.
(11.30) can be rewritten asPe [ 2:ð11:32ÞThus for a cell Péclet number (Pe) larger than 2 the discretization process becomesinconsistent as now an increase in the neighboring value will lead to a decrease inthe value at C. This in turn will lead to a further increase in the neighboring value,and the error is magnified.This situation can be avoided by decreasing the grid size so that the cell Pécletnumber is smaller than 2. For many practical situations, however, the increase instorage and computational requirements may be too large to afford. Moreover forpurely convected flows (e.g., Euler flows) this is simply not feasible.
Therefore aremedy is needed.11.2.4 The Upwind SchemeWhen examining the discretization procedure described above, it is noticed that thereason for obtaining these positive coefficients is because of the adopted linearsymmetric profile. A linear symmetric profile gives equal weights to the two nodessharing the face with no directional preference, which is appropriate fornon-directional phenomena with an elliptic type term, such as the diffusion term.It is simply not adequate for the convection term [4].A scheme that is more compatible with the advection process is the upwindscheme [4, 6] schematically displayed in Fig.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
