Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 63
Текст из файла (страница 63)
The QUICK scheme profile is shown to be much sharper and moreaccurate than the upwind profile however it is infected with over/undershoots nearthe sharp gradient. As mentioned earlier this error is called the dispersion error,which causes the generation of maxima/minima in the solution domain and is acharacteristic of all High Order (HO) schemes.40411 Discretization of the Convection Term11.6.1 Error SourcesBased on the discussions in the previous sections the sources of numerical error inthe discretization of the convective flux can be divided into numerical diffusion andnumerical dispersion.Numerical diffusion, which causes smearing of sharp gradients (Fig.
11.18a) canalso be divided into stream wise and cross stream numerical diffusion. Stream wisenumerical diffusion can be reduced by increasing the order of the interpolationprofile, as shown in Fig. 11.18b, resulting in sharper profiles but inducingunder/overshoots in the presence of large gradients. Cross-stream numerical diffusion, shown in Fig.
11.17b, is caused by the one dimensional nature of theassumed profile and can be reduced either by interpolating in the direction of theflow (i.e., multi-dimensional profiles) [15, 16] or by using one-dimensional higherorder interpolation profiles (Fig.
11.18b).Numerical dispersion error manifests itself through oscillation in the generatedprofile in the presence of large gradients rendering the solution unbounded. Asshown in Figs. 11.18b and 11.19 and the results reported in Fig. 11.14b, it existswith all assumed interpolation profiles except the upwind scheme. It is the result ofthe unphysical behavior of the assumed interpolation profile.(a)CfDUDD(b)CUfDDDFig.
11.18 a Smearing of sharp profiles by numerical diffusion and b wiggles in the computedprofile due to numerical dispersion11.6Steady Two Dimensional Advection405fCDUDDFig. 11.19 Numerical dispersion error causing oscillations in the presence of a large gradientAn evaluation of this error can be obtained by using a simplified version ofEq. (11.73) in which diffusion and sources are neglecting. If further, velocity anddensity fields are considered constants, Eq. (11.73) simplifies to@/@/þu¼0@t@xð11:130ÞAssuming an exact solution of the form/ðx; tÞ ¼ /ðtÞe jkxð11:131Þwhere j is the imaginary number defined by j2 = −1.
Then the exact value of thegradient becomes@/ðx; tÞ ¼ jk/ðtÞe jkx ¼ jk/ðx; tÞ:@xð11:132ÞWith an interpolation profile, the numerical approximation of the gradient is writtenin terms of / values at locations MDx; ðM þ 1ÞDx; . . .; Dx; 2Dx; . . .; NDx asNN@1 X1 X/ðx; tÞ an /ðx þ nDx; tÞ ¼an /ðtÞe jkðxþnDxÞ@xDx n¼MDx n¼Mð11:133Þwhich, upon substitution of the assumed solution, becomesNN@1 X1 X/ðx; tÞ an /ðtÞe jkðxþnDxÞ ¼an e jknDx /ðx; tÞ:@xDx n¼MDx n¼Mð11:134ÞBy comparing the exact and numerical solutions, it is found thatk¼Nj Xan e jknDx :Dx n¼Mð11:135Þ40611 Discretization of the Convection TermIn general k is an imaginary number that can be written ask ¼ Reðk Þ þ j ImðkÞð11:136Þwhere Re and Im refer to the real and imaginary part, respectively.
Inserting k in theexact solution, the approximate solution is obtained as/ðx; tÞ ¼ /ðtÞe jkxðk ÞxImðk Þx /ðtÞe j½ReðkÞþjImðkÞx ¼ /ðtÞ |fflffle jReffl{zfflfflffl} e|fflfflffl{zfflfflffl}ð11:137ÞPhaseAmplitudeDispersion DissipationTherefore, the numerical solution may include both diffusion (or dissipative) anddispersion errors. If k is real, only dispersion error occurs. However if k is complex,then both types of errors will arise.
Based on this analysis, the value of k for the upwindand CD schemes can be checked. For the upwind scheme, the gradient is computed as@/ /e /w /C /W e jkx ejkðxdxÞ’¼¼/ðx; tÞ@xDxDxDx1 cosðkdxÞ þ j sinðkdxÞ/ðx; tÞ¼DxsinðkdxÞ1 cosðkdxÞ¼ jk/ðx; tÞ ) k ¼jDxDxð11:138ÞIt is clear that the upwind scheme gives rise to both types of errors. For the centraldifference scheme, the gradient is given by@/ /e /w /E /W e jkðxþdxÞ e jkðxdxÞ’¼¼/ðx; tÞ@xDx2Dx2Dx1sinðkdxÞj sinðkdxÞ/ðx; tÞ ¼ jk/ðx; tÞ ) k ¼¼DxDxð11:139ÞSince k is imaginary, then only numerical dispersion error arises.
This dispersionerror causes oscillations and under/overshoots in the solution.Having developed a better understanding of numerical dispersion, it is desirableto develop convective non-oscillatory schemes of high order of accuracy. This haskept researchers busy for an extended period of time until the bounding of theconvective flux became understood. The development of such schemes will bedetailed in the next chapter.11.7High Order Schemes on Unstructured GridsAs for structured grids, the functional relationships of HO schemes are defined asfunctions of the values at the U, C, and D nodes.
While the C and D nodes arereadily available for any interior face (Fig. 11.20a), defining the U node is not11.7High Order Schemes on Unstructured Grids407(a)CFdCF(c)(b)?UCDfDCf?UdCDdCDFig. 11.20 a An element in an unstructured grid; the upwind node for HO and HR convectionschemes in unstructured grids when the velocity at the element face is b negative or c positivestraightforward (Fig.
11.20b, c) in an unstructured grid. A simple way around thisdifficulty is to simply redefine HO schemes in terms of the gradients at the C andD nodes, or a combination of thereof. Another approach is to reconstruct a pseudoU node, which will be detailed in the next chapter and used in developing HRschemes.11.7.1 Reformulating HO Schemes in Terms of GradientsThis approach relies on profiles developed over structured grids and is bestexplained by rewriting the QUICK [2] scheme using the new terminology and thengeneralizing results for any second order profile developed using three points.
Thefunctional relationship of the QUICK [2] scheme can be written as1 /D /U1/f ¼ /C þþ ð/D /C Þ442ð11:140ÞBy approximating the gradients at the locations C and f in the dCF or f direction asshown in Fig. 11.19, then the gradients at the centroids C and f are computed as@/C /D /U¼@f2Df@/f /D /C¼@fDfð11:141Þ40811 Discretization of the Convection TermUsing Eq. (11.141), Eq. (11.140) can be recast into/f ¼ /C þ1 /D /C Df 1 /D /C Dfþ22222DfDfð11:142Þor/f ¼ /C þ1 @/C Df 1 @/f Dfþ2 @f 22 @f 2ð11:143ÞDenoting the vector between C and f by dCf , in vector form the above equationbecomes11/f ¼ /C þ r/C dCf þ r/f dCf22ð11:144Þwhich is quite suitable for use in the context of unstructured grids since it onlyrequires information related to gradient at the C and f locations.
As long as thecomputation of these gradients is second order accurate, the way they are calculatedbecomes immaterial. This gives higher flexibility, as compared to the originalformulation, over unstructured grids. Furthermore, Eq. (11.144) suggests that aprofile based on three points can be written as/f ¼ a/C þ br/C dCf þ cr/f dCfð11:145Þwhere the constants a, b, and c are determined by equating /f to the profileobtained over structured grids. The general discretized form of Eq.
(11.145) can befound once and then used in all subsequent derivations. This is derived by firstsubstituting the approximate values of the gradients using Eq. (11.141) to yield/f ¼ a/C þ br/C dCf þ cr/f dCf/ /U Df/ /C Dfþc D¼ a/C þ b D22Df2Df 2ð11:146Þand then after some algebraic manipulations the final form is obtained ascb cb/f ¼ a /C þþ /D /U24 44ð11:147ÞUsing the above equation, the calculation of the a, b, and c coefficients is easilyaccomplished. For example the SOU profile in terms of the gradients can be foundas follows:11.7High Order Schemes on Unstructured Grids4098b 1>9>¼ )b¼2>>cb cb>>4 2>/f ¼ a /C þþ /D /U >=<24 24) b þ c ¼ 0 ) c ¼ 1>>314 2>>;>¼ /C /U>>>22:a c ¼ 3 ) a ¼ 12 2Thus the equivalent form of the SOU scheme is given by/f ¼ /C þ 2r/C r/f dCfð11:148Þð11:149ÞFollowing the same procedure, the gradient forms of the schemes presented earlierare found to beUpwind scheme :/f ¼ /Cð11:150ÞCentral difference :/f ¼ /C þ r/f dCfð11:151ÞSOU scheme :/f ¼ /C þ 2r/C r/f dCfð11:152ÞFROMM scheme :/f ¼ /C þ r/C dCfð11:153ÞQUICK scheme :/f ¼ /C þDownwind scheme :/f ¼ /C þ 2r/f dCf1r/C þ r/f dCf2ð11:154Þð11:155ÞA full chapter was devoted to the calculation of the gradient at the element centroidand faces and the reader is referred to Chap.
9 for the calculation of r/C and r/f .11.8The Deferred Correction ApproachThe Deferred Correction (DC) procedure of Khosla and Rubin [17] is a compactingtechnique that enables the use of HO schemes in codes initially written for loworder schemes without violating any of the stability rules. The approach is applicable on any type of structured or unstructured grid systems.
The method is basedon writing the convection flux at a cell face f calculated using a HO scheme asUHOU__m_ f /HO¼m/þm//ð11:156Þffffff|fflffl{zfflffl} |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}implicitexplicitwhere the superscripts U and HO refer to the Upwind and High Order scheme,respectively. By expressing the convection flux in this way, the first term on theright hand side (RHS) is implicitly evaluated by expressing it in terms of nodal41011 Discretization of the Convection Termvalues, while the second term on the RHS is evaluated explicitly using the latestavailable / values, i.e., values from the previous iteration in an iterative solutionprocedure. In terms of nodal values, Eq.
(11.156) is expressed as_ f ; 0 k /C þ k m_ f ; 0 k /F¼ k m_ f ; 0 k /C k m_ f ; 0 k /F þ m_ f /HRm_ f /HOff k m¼ FluxCf /C þ FluxFf /F þ FluxVf|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflffl{zfflfflffl}implicitexplicitð11:157ÞwhereFluxCf ¼k m_ f ; 0 kFluxFf ¼ k m_ f ; 0 kFluxVf ¼m_ f /HRfð11:158Þ FluxCf /C FluxFf /FSubstituting the convection flux given by Eqs. (11.157) and (11.158) inEq. (11.156) and rearranging, the final form of the algebraic equation is obtained asXaC /C þaF /F ¼ bCð11:159ÞF NBðC ÞwhereaF ¼ FluxFf ¼ km_ f ; 0 kXXX FluxCf ¼km_ f ; 0 k¼m_ f þ km_ f ; 0kaC ¼f nbðC Þ¼XF NBðC Þf nbðC ÞaF þXf nbðC Þð11:160Þm_ ff nbðC ÞandbC ¼ Q/C VC XFluxVff nbðC Þ|fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}¼ Q/C VC XbDCC /Um_ f /HOffð11:161Þf nbðC Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}bDCCthe bDCC represents the extra source term arising due to the DC procedure. Moreover,the DC technique results in an equation for which the coefficient matrix is alwaysdiagonally dominant since it is formed using the upwind scheme.11.8The Deferred Correction Approach411Example 2Derive the coefficients for a deferred correction implementation of theQUICK schemeSolutionThe conservation equation is discretized into the following algebraic equationin every cell:XaF / F ¼ b CaC /C þF NBðC ÞFor the quick scheme the value at a cell face is computed as313/f ¼ /C /U þ /D488The coefficients of the algebraic system of equations in a deferred correctionapproach are based on the UPWIND convection scheme.
As such thesecoefficients are given byaF ¼ km_ f ; 0 kXXaF þaC ¼ m_ fF NBðC Þf nbðC ÞThe difference between the UPWIND and QUICK schemes is explicitlyaccounted for in the source term which is modified tobC ¼S/C VC393/ / þ /m_ f4 C 8 U 8 Df nbðC ÞXThis compacting procedure is simple to implement, however as the differencebetween the cell face values calculated with the upwind scheme and that calculatedwith the HO scheme becomes larger, the convergence rate diminishes.11.9Computational Pointers11.9.1 uFVMThe assembly of the convection term in uFVM is performed in cfdAssembleConvectionTerm, where the assembly for the upwind scheme is performedfirst for interior faces (cfdAssembleConvectionTermInterior), then for boundary facesby looping over the various boundary patches (cfdAssembleConvectionTermInletBC,cfdAssembleConvectionTermOutletBC, etc.).41211 Discretization of the Convection TermThe core of the interior faces assembly routine is shown in Listing 11.1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
