Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 61
Текст из файла (страница 61)
In addition, DD and UUdenote the nodes downstream and upstream of D and U, respectively. The corresponding values at these locations are denoted by /DD ; /D ; /C ; /U and /UU ,respectively. This notation is displayed graphically in Fig. 11.10 for the cases whenthe velocity at the face is either positive (→) or negative (←).11.5.1 Second Order Upwind SchemeA second order scheme, requires the use of a linear profile, as was the case with thecentral difference scheme. However, instead of a symmetric profile, an upwindbiased stencil is used [10]. As depicted in Fig.
11.11, the linear profile is constructed by employing the / values at nodes C and U. Therefore the value at theface is actually calculated by extrapolation rather than interpolation.39011 Discretization of the Convection TermfDDDCUUUUUUDCDDfvfFig. 11.11 The Second Order Upwind (SOU) scheme profile11.5.2 The Interpolation ProfileStarting with the linear profile/ð xÞ ¼ k0 þ k1 ðx xC Þð11:83Þand fitting it to the nodal values at xC and xU where the / values are /C and /U ,respectively, the profile becomes/ ð xÞ ¼ / C þ/C /Uð x xC ÞxC xUð11:84ÞUsing the above equation, the / value at face f, shown in Fig.
11.11, is given by / /U /f ¼ / xf ¼ /C þ Cxf xCxC xUð11:85Þwhich, for a uniform grid reduces to31/f ¼ /C /U22ð11:86Þ11.5.3 The Discretized EquationUsing this profile to approximate the interface values in the discretized onedimensional convection diffusion equation (Eq. 11.16), the fluxes at the faces areobtained as11.5Higher Order Upwind Schemes3131/ / km_ e ; 0k / /m_ e /e ¼km_ e ; 0k2 C 2 W2 E 2 EE3131/ / km_ w ; 0k / /m_ w /w ¼km_ w ; 0k2 C 2 E2 W 2 WW391ð11:87ÞSubstitution of these fluxes in Eq. (11.16), the discretized form is transformed to313131/C /W km_ e ; 0k /E /EE km_ e ; 0k þ/C /E km_ w ; 0k22222231/ S/ S/W /WW km_ w ; 0k Cð/E /C Þ Cð/C /W Þ ¼ 022dx edx wð11:88Þwhich can be modified into the formaC /C þ aE /E þ aW /W þ aEE /EE þ aWW /WW ¼ 0ð11:89ÞwhereSe 311 km_ e ; 0k km_ w ; 0k aEE ¼ FluxFee ¼ km_ e ; 0k22dxe 2S311waW ¼ FluxFw ¼ C/w km_ w ; 0k km_ e ; 0k aWW ¼ FluxFww ¼ km_ w ; 0k22dxw 2XFluxCfaC ¼aE ¼ FluxFe ¼ C/ef nbðC ÞSeSw 33þ C/wþ km_ e ; 0k þ km_ w ; 0k2dxedxw 2¼ ðaE þ aW þ aEE þ aWW Þ þ ðm_ e þ m_ w Þ¼ C/eð11:90Þ11.5.4 Truncation ErrorFollowing a procedure similar to the one used with the central difference scheme,the truncation error is found to be31 3 ivTE ¼ Dx2 /000C Dx /C þ 84indicating second order accuracy.ð11:91Þ39211 Discretization of the Convection Term11.5.5 Stability AnalysisTo check the stability of the SOU scheme, the discretized convective fluxes aresubstituted in Eq.
(11.75) resulting in the following RHSconvection term:RHSconvection3131/ /¼ /C /W km_ e ; 0k þkm_ e ; 0k222 E 2 EE3131/ / km_ w ; 0k þ/ /km_ w ; 0k2 C 2 E2 W 2 WWð11:92ÞThe rate of change of this term with respect to /C is found as@ ðRHSconvection Þ33¼ km_ e ; 0k km_ w ; 0k@/C22ð11:93Þwhich is always negative indicating a stable scheme, i.e., the solution error willalways remain under control. It should be kept in mind that this is true for theconditions stated in the derivations and not for the general case when the velocity isnot constant.11.5.6 The QUICK SchemeThe Quadratic Upstream Interpolation for Convective Kinematics (QUICK) schemewas developed by Leonard [2].
The method is based on interpolating the value ofthe dependent variable at each face of the element by using a quadratic polynomialbiased toward the upstream direction, as shown in Fig. 11.12. The interpolatedvalue is used to calculate the convective term in the governing equations for thedependent variable.
The calculation of the values of the dependent variable at a cellface is detailed next.DfUUUDDCUUUDCfvfFig. 11.12 The QUICK scheme profileDD11.5Higher Order Upwind Schemes39311.5.7 The Interpolation ProfileIn its original form, the QUICK scheme required the use of a general multi dimensionalsecond degree polynomial for the calculation of the unknown variable /. However it isgenerally sufficient to treat the flow as locally one dimensional [11] and to use a secondorder one dimensional profile in each coordinate direction. For the one dimensionalconvection-diffusion problem under discussion, the value of / is computed using/ ¼ k0 þ k1 x þ k2 x2 ;ð11:94Þsubject to8< /U/ ¼ /C:/Dat x ¼ xUat x ¼ xCat x ¼ xDð11:95ÞApplying the conditions given by Eq.
(11.95), the profile for / is found to be/ ¼ /U þðx xU Þðx xC Þðx xU Þðx xD Þð/ /U Þ þð/ /U ÞðxD xU ÞðxD xC Þ DðxC xU ÞðxC xD Þ Cð11:96ÞFor the case of a uniform grid, the value of / at the cell face f reduces to/f ¼/C þ /D /D 2/C þ /U28ð11:97ÞThus, the convective fluxes at the element faces can be computed as313313/C /W þ /E km_ e ; 0k /E /EE þ /C km_ e ; 0k488488313313/C /E þ /W km_ w ; 0k /W /WW þ /C km_ w ; 0km_ w /w ¼488488m_ e /e ¼ð11:98ÞSubstituting back in the one dimensional convection-diffusion equation[Eq. (11.16)], the discretized form becomes313313/C /W þ /E km_ e ; 0k /E /EE þ /C km_ e ; 0k488488313313þ/C /E þ /W km_ w ; 0k /W /WW þ /C km_ w ; 0k ð11:99Þ488488SSð/ /C Þ C/ð/ / W Þ ¼ 0 C/dx e Edx w C39411 Discretization of the Convection Termwhich can be modified into the formaC /C þ aE /E þ aW /W þ aEE /EE þ aWW /WW ¼ 0ð11:100Þwith the coefficients given bySe 331 km_ e ; 0k þ km_ e ; 0k km_ w ; 0kdxe 488331/ Sw km_ w ; 0k þ km_ w ; 0k km_ e ; 0kaW ¼FluxFw ¼ Cwdxw 48811aEE ¼FluxFee ¼ km_ e ; 0k aWW ¼ FluxFww ¼ km_ w ; 0k88XFluxCfaC ¼aE ¼FluxFe ¼ C/eð11:101Þf nbðC ÞSeSw 3333þ C/wþ km_ e ; 0k km_ e ; 0k þ km_ w ; 0k km_ w ; 0kdxedxw 4848¼ ðaE þ aW þ aWW þ aEE Þ þ ðm_ e þ m_ w Þ¼ C/e11.5.8 Truncation ErrorAgain following the procedure described above, the truncation error is found to beTE ¼13Dx3 /ivC Dx4 /vC þ 16128ð11:102Þwhich is clearly third order accurate.11.5.9 Stability AnalysisTo check the stability of the QUICK scheme, the discretized convective fluxes aresubstituted in Eq.
(11.75) resulting in the following RHSconvection term:313313/C /W þ /E km_ e ; 0k þ/E /EE þ /C km_ e ; 0k488488313313/C /E þ /W km_ w ; 0k þ/W /WW þ /C km_ w ; 0k488488RHSconvection ¼ ð11:103Þ11.5Higher Order Upwind Schemes395The rate of change of this term with respect to /C is found to be given by@ ðRHSconvection Þ333¼ km_ e ; 0k km_ w ; 0k ðm_ e þ m_ w Þ@/C888ð11:104Þwhich for a uniform velocity is always negative indicating a stable scheme.However this does not guarantee solution boundedness especially in the generalcase of nonuniform velocity.11.5.10The FROMM SchemeThe FROMM scheme [12] fits a linear profile between the far Upwind (U) andDownwind (D) nodes straddling the interface. As depicted in Fig. 11.13, instead ofa symmetric profile, an upwind biased stencil is used to calculate the value of thedependent variable / at face f. Based on the adopted profile, /U , /C , and /D areassumed to be collinear.11.5.11The Interpolation ProfileStarting with the linear profile/ð xÞ ¼ k0 þ k1 ðx xC Þð11:105Þand fitting it to the nodal values at xD and xU where the / values are /D and /U ,respectively, the profile becomes/ð xÞ ¼ /U þ/D /Uð x xU ÞxD xUCð11:106ÞfDDDUUUUUUDCfvfFig.
11.13 The FROMM scheme profileDD39611 Discretization of the Convection TermUsing the above equation, the / value at upwind node C is obtained as/C ¼ /U þ/D /UðxC xU ÞxD xUð11:107ÞFor the case of a uniform grid the above equation becomes/C ¼/D þ /U2ð11:108ÞThe required value at the face f, shown in Fig. 11.13, is given by / /U xf xC/f ¼ / xf ¼ /U þ Dxf xU ¼ /C þð/ /U ÞxD xUxD xU Dð11:109Þwhich, for a uniform grid reduces to/f ¼ /C þ/D /U4ð11:110ÞThe final expression for /f given in Eq.
(11.109) was obtained by invokingEq. (11.107).11.5.12The Discretized EquationUsing this profile to approximate the interface values in the discretized onedimensional convection diffusion equation (Eq. 11.16), the fluxes at the faces areobtained as1111m_ e /e ¼ /C /w þ /E km_ e ; 0k /E /EE þ /C km_ e ; 0k44441111m_ w /w ¼ /C /E þ /W km_ w ; 0k /W /WW þ /C km_ w ; 0k4444ð11:111ÞSubstitution of these fluxes in Eq. (11.16), the discretized form is transformed to1111/C /W þ /E km_ e ; 0k /E /EE þ /C km_ e ; 0k44441111þ /C /E þ /W km_ w ; 0k /W /WW þ /C km_ w ; 0k4444SSð/ /C Þ C/ð/ /W Þ ¼ 0 C/dx e Edx w Cð11:112Þ11.5Higher Order Upwind Schemes397which can be modified into the formaC /C þ aE /E þ aW /W þ aEE /EE þ aWW /WW ¼ 0ð11:113ÞwhereSe 11þ km_ e ; 0k km_ e ; 0k km_ w ; 0k4dxe 4S11waW ¼ FluxFw ¼ C/wþ km_ w ; 0k km_ w ; 0k km_ e ; 0k4dxw 411aEE ¼ FluxFee ¼ km_ e ; 0k aWW ¼ FluxFww ¼ km_ w ; 0k44XaC ¼FluxCfaE ¼ FluxFe ¼ C/ef nbðCÞSeSw11þ C/wþ km_ e ; 0k km_ e ; 0k þ km_ w ; 0k km_ w ; 0k44dxedxw¼ ðaE þ aW þ aEE þ aWW Þ þ ðm_ e þ m_ w Þ¼ C/eð11:114Þ11.5.13Truncation ErrorFollowing a procedure similar to the one used with the central difference scheme,the truncation error can be derived to be of OðDx2 Þ, i.e.,TE ¼ O Dx2ð11:115Þindicating second order accuracy.11.5.14Stability AnalysisTo check the stability of the FROMM scheme, the discretized convective fluxes aresubstituted in Eq.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
