Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 60
Текст из файла (страница 60)
To calculate this truncation error, the followingTaylor series expansions with respect to the exact /e and /w values are needed:Dx 0Dx2 00 Dx3 000 Dx4 ivDx5 v/w þ/w /w þ/w / þ 28483843840 wDx 0Dx2 00 Dx3 000 Dx4 ivDx5 v/w þ/C ¼ /w þ/w þ/w þ/w þ/ þ 28483843840 wDx 0 Dx2 00 Dx3 000 Dx4 ivDx5 v/e þ/e /e þ/e / þ /C ¼ /e 28483843840 eDx 0 Dx2 00 Dx3 000 Dx4 ivDx5 v/e þ/E ¼ /e þ/e þ/e þ/e þ/ þ 28483843840 e/W ¼ /w ð11:64ÞUsing these expansions, the average values at the faces are obtained as1Dx2 00 Dx4 iv/ þ/ þ ð/E þ /C Þ ¼ /e þ28 e 384 e2Dx 00Dx4 iv/e ¼ ð/C 2/e þ /E Þ / þ )4192 e1Dx2 00 Dx4 iv/ þ/ þ ð/C þ /W Þ ¼ /w þ28 w 384 wDx2 00Dx4 iv/w ¼ ð/W 2/w þ /C Þ / þ )4192 wð11:65ÞSubtracting the above two terms given in Eq. (11.65), their difference is found to be11Dx2 00/e /00wð/E þ /C Þ ð/C þ /W Þ ¼ ð/e /w Þ þ228Dx4 ivþ/ /ivw þ 384 eð11:66ÞFurther, the expanded forms of the exact /e and /w with respect to /C are given byDx 0Dx2 00 Dx3 000 Dx4 ivDx5 v/C þ/C þ/C þ/C þ/ þ 28483843840 CDx 0Dx2 00 Dx3 000 Dx4 ivDx5 v/C þ/w ¼ /C /C /C þ/C / þ 28483843840 C/e ¼ /C þð11:67ÞIn addition, /E and /W can be expanded in terms of /C asDx2 00 Dx3 000 Dx4 iv Dx5 v/ þ/ þ/ þ/ þ 2 C6 C24 C 120 CDx2 00 Dx3 000 Dx4 iv Dx5 v/W ¼ /C Dx/0C þ/ / þ/ / þ 2 C6 C24 C 120 C/E ¼ /C þ Dx/0C þð11:68Þ11.3Truncation Error: Numerical Diffusion and Anti-Diffusion385Using the above two sets of equations, the following is obtained:Dx2 00 Dx3 000 Dx4 ivDx5 v/ / / / þ 4 C24 C 192 C 1920 C2345Dx 00 Dx 000 Dx iv Dx vþ /C þ Dx/0C þ/ þ/ þ/ þ/ þ 2 C6 C24 C 120 CDx2 00 Dx3 000 7Dx4 iv 15Dx5 v/ þ/ þ/ þ/ þ ¼4 C8 C192 C1920 C2345Dx 00 Dx 000 7Dx iv 15Dx v/W 2/w þ /C ¼/ / þ/ / þ 4 C8 C192 C1920 C/C 2/e þ /E ¼ /C 2/C Dx/0C ð11:69ÞThen /00e /00w is computed as Dx3 000 Dx5 vDx2 00/e /00w ¼/ þ/ þ 88 C 128 Cð11:70ÞSubstituting back, Eq.
(11.66) is transformed to11Dx3 000 Dx5 v/ þ/ þ ð/E þ /C Þ ð/C þ /W Þ ¼ /e /w þ228 C 128 Cð11:71ÞDividing throughout by Dx, the truncation error associated with the calculation ofthe gradient is obtained asTE ¼Dx2 000 Dx4 v/ þ/ þ 8 C 128 Cð11:72Þwhich indicates that the method is second order accurate.11.4Numerical StabilityThe confusion related to truncation error (first order for the physical solution andsecond order for the unphysical one) has led many workers to extrapolate that sincecentral differencing of the diffusion term is so accurate, then central differencing ofthe convection term should also be similarly accurate.
Of course, as seen in theprevious sections, central differencing of the convection term can lead to unphysicalsolutions. The cause for this behavior is that the central difference scheme has noinherent convective stability when applied to derivatives of odd order like theconvection term.Leonard [9] advanced the convective stability concept as an explanation tooscillation in numerical solutions generated by applying the central differencescheme to convection dominated flows.
Leonard used the general one dimensionalunsteady convection-diffusion equation with constant velocity given by38611 Discretization of the Convection TermCeEwWwWeuwCueEFig. 11.9 Insensitivity of CD convection term to the values of /C@ ðq/Þ@ ðqu/Þ @@/¼þC/þ Q/@t@x@x@xð11:73Þto describe this concept. When applying this equation over the element of centroidC shown in Fig. 11.9, its left hand side (LHS) represents the rate of change of /Cwithin the control cell per unit time, and its right hand side (RHS) represents the netinflux across the element surface and source terms within the element that areaffecting the value of /C .
If there were numerical errors in the RHS then the valueof /C calculated from Eq. (11.73) would either increase or decrease depending onthe scheme used in the discretization process. In an unstable scheme, a smalldeviation from the correct value of /C gives a corresponding increase/decrease inthe net influx represented by the RHS. When an iterative procedure is used as part ofthe solution mechanism, an increase/decrease in the net influx will furtherincrease/decrease the value of /C at each subsequent step of the iterative process. Ina stable scheme this change in /C due to errors in the RHS should feed backnegatively into the RHS as a self correction device.
Clearly, for this kind ofnumerical stability, the RHS should satisfy@ ðRHSÞ\0@/Cð11:74Þindicating that the sensitivity to /C of the combination of the modeled terms on theright hand side of Eq. (11.73) should be negative. In this case, an increase/decreasein /C will correspond to a decrease/increase in the influx, which will in turn pushes/C downward/upward towards its correct value. However, stability should not beconfused with boundedness or accuracy.
A stable numerical scheme could actuallybe unbounded giving rise to over/under shoots and oscillations/wiggles or be verydiffusive and give results that are of low accuracy. The stability here refers tocontrolling the numerical error to remain bounded in order for it not to increaseindefinitely as was the case with the central difference scheme where the normalized11.4Numerical Stability387value of /C (Fig. 11.4) was found to vary from +∞ to −∞ while varying Pe from−∞ to +∞, even though the correct value of /C should vary between 1 and 0.
It sohappens that the upwind scheme possesses both the boundedness and stabilitycharacteristics.With this in mind, the general discretized form of the RHS of Eq. (11.73) isgiven byRHS ¼ ðquDyÞe /e þ ðquDyÞw /wDyDyþ C/ð/ E / C Þ C /ð/C /W Þ þ Q/C VCdx edx wð11:75ÞThus using the central difference scheme, Eq. (11.75) becomesð/ þ / C Þð/ þ /W ÞRHSCD ¼ ðquDyÞe Eþ ðquDyÞw C22/ Dy/ Dyþ Cð/E /C Þ Cð/C /W Þ þ Q/C VCdx edx wð11:76ÞAnalyzing the central difference scheme using the above criteria, it is found that forthe diffusion term the sensitivity is given by@ RHSDiffCD@/C¼ 2C/Dydxð11:77Þwhich is negative since C/ is positive, indicating a stable scheme. However, for theconvective term the sensitivity equation gives@ RHSConv1CD¼ ðm_ e þ m_ w Þ2@/Cð11:78Þwhich is equal to zero for steady flows but not necessarily for unsteady flows.
Forunsteady situations, its value will be positive for decelerating flows. In a generalflow, such regions act as wiggle sources and can easily lead to a total numericalcatastrophe when the Péclet number is large enough. Even for steady flows, as itsvalue is zero, it cannot feed back into the equation to act as a self correction device.Equation (11.78) also indicates that for steady flow the net convective flux computed with the CD scheme is independent of the value of /C . Therefore, as shownin Fig. 11.9, the different possible values of /C will result in the same net convective flux over the element of centroid C.38811 Discretization of the Convection TermWith the upwind scheme, the RHS can be obtained from Eq. (11.36) asRHSUpwind ¼ km_ e ; 0k/C þ km_ e ; 0k/E km_ w ; 0k/C þ km_ w ; 0k/Wð11:79ÞSSþ C/ð/ E / C Þ C /ð/ /W Þ þ Q/Cdx edx w CThe sensitivity is expected to be negative since the scheme was found to be stablefor all Péclet number, indeed@ RHSConvUpwind@/C¼ km_ e ; 0k km_ w ; 0kð11:80Þwhich is negative or equal to zero for all flows.
It will be equal to zero when bothmass flow rates are negative, a situation that does not arise in a one dimensionalsituation of constant cross-sectional area. When added to the false diffusion introduced by the first order approximation, Eq. (11.80) indicates that the scheme is verystable. However, this stability is achieved at the expense of accuracy, as wasdemonstrated in the previous section.For the downwind scheme the RHS is given byRHSDownwind ¼ km_ e ; 0k/E þ km_ e ; 0k/C km_ w ; 0k/W þ km_ w ; 0k/C/ S/ Sþ Cð/ /C Þ þ Cð/ /C Þdx e Edx w Wð11:81ÞThe sensitivity is given by@ RHSConvDownwind¼ km_ e ; 0k þ km_ w ; 0k@/Cð11:82Þwhich is always positive or equal to zero for all flows. When this is added to theanti-diffusion effects it gives a highly unstable scheme.11.5Higher Order Upwind SchemesThe previous sections demonstrated that both the upwind and central differenceschemes have severe limitations, the former because of its poor accuracy due tonumerical diffusion, and the latter because of its instability also known as numericaldispersion error.
These shortcomings have provoked a great deal of research toimprove the accuracy and stability of advection schemes by using higher order11.5Higher Order Upwind Schemes389DCDDvfUUUUUDCUDDfCDUvfDDUUDDDCUUUfFig. 11.10 A schematic showing the UU, U, C, D and DD node locations used in describingconvection schemesupwind biased interpolation profiles. These higher-order schemes aim at producingat least a second order accurate solutions, while being unconditionally stable.To underline the conservation property which associates fluxes with cell faces,not nodes, in what follows D, C, and U will be used to denote the Downwind,Upwind, and far Upwind nodes at any particular face.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
