Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 70
Текст из файла (страница 70)
For that purpose the multi dimensional advection equation with a sourceis considered. Again the velocity field is assumed to be known and the conservationequation in vector form is given byr ðqv/Þ ¼ Q/ð12:67Þ12.7UUDeferred Correction for HR SchemesU457FCfvfSfDDFig. 12.16 Two three-dimensional elements, which represent part of an unstructured grid systemIntegrating over the element of volume VC shown in Fig.
12.16, applying thedivergence theorem, and replacing the surface integral by a summation over theelement faces, Eq. (12.67) becomesXðqv/Þf Sf ¼ Q/C VCð12:68Þf nbðCÞNoticing that the mass flow rate at a cell face is given bym_ f ¼ ðqvÞf Sfð12:69Þthen Eq. (12.68) can be rewritten asXm_ f /f ¼ Q/C VCð12:70Þf nbðC ÞThe value of /f is obtained using any of the advection schemes presented earlierkeeping in mind that, in order to be able to solve for the unknown values at the mainnodes, the algebraic form of the discretized equation should look likeaC /C þXaF / F ¼ b Cð12:71ÞF NBðC ÞThe difficulty lies in the instability that arises when expressing /f in terms ofnodal values as described next.45812 High Resolution Schemes12.7.1 The Difficulty with the Direct Use of Nodal ValuesThe difficulty that arises when explicitly expressing /f in terms of neighboringvalues will be explained by discretizing the convection flux using a TVD scheme.Referring to Fig.
12.16, the convective flux at a face f is written as23rfþzfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflffl{67671/C /Um_ f ; 06m_ f /f ¼ 6/C þ wð/F /C Þ772/F /C4523rfzfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{671/F /DD67 6/F þ wð/C /F Þ7m_ f ; 0452/C /Fð12:72Þ1 ¼ /C þ w rfþ ð/F /C Þ m_ f ; 021 /F þ w rf ð/C /F Þ m_ f ; 02Substituting Eq. (12.72) into Eq. (12.70) yieldsXaC /C þaF / F ¼ b Cð12:73ÞF NBðC Þwhere 1 1 aF ¼ FluxFf ¼ m_ f ; 0 þ w rfþ m_ f ; 0 þ w rf m_ f ; 02XX2XaC ¼FluxCf ¼ aF þm_ ff nbðC ÞF NBðC Þð12:74ÞF NBðC ÞbC ¼ Q/C VCTo see the weakness in this formulation, Eq.
(12.67) is simplified to the onedimensional problem depicted in Fig. 12.8. Assuming the flow to be in the positivedirection and using the terminology displayed in the figure, Eq. (12.73) becomesaC /C þ aE /E þ aW /W ¼ bCð12:75Þ12.7Deferred Correction for HR Schemes459where1 aE ¼ FluxFe ¼ w reþ m_ e21 aW ¼ FluxFw ¼ 1 þ w rw m_ e2PaC ¼FluxCf ¼ ðaE þ aW Þð12:76Þf nbðC ÞbC ¼ Q/C VCSince 0 wðr Þ 2; the aE and aW coefficientswill be of opposite signs (except for the UPWIND scheme, i.e., when w rf ¼ 0), thereby violating one of the basicrules for stability and causing convergence difficulties of the iterative procedure.Following a similar approach with the NVF leads to the same shortcomings.A remedy, which was presented in the previous chapter, is the deferred correction(DC) procedure, in which the coefficients are based on the upwind scheme, while thedifference between the HR and upwind schemes is added as a source term in thealgebraic equation.
The DC procedure is simple to implement and can be used instructured and unstructured grid systems, however as the difference between the cellface values calculated with the upwind scheme and that calculated with the HR schemebecomes larger, the convergence rate diminishes. This effect can be easily estimated onan NVD; the difference between the UPWIND line and that of the chosen HR schemeis the normalized difference between the cell face values. The larger this difference is,the lower the convergence rate will be. This has enticed researchers to look for othertechniques for implementing HR schemes that are more implicit not affecting theconvergence rate.
Two of these techniques are described in the next section.12.8The DWF and NWF MethodsSeveral techniques have been developed to overcome the reduction in the convergence rate associated with the use of the Deferred Correction (DC) procedure forthe implementation of HR schemes. Two of these methods are presented below,namely the Downwind Weighing Factor (DWF) method of Leonard and Mokhtari[25] (implemented in OpenFOAM®) and the Normalized Weighing Factor(NWF) method of Darwish and Moukalled [29].The implementation details for both methods are presented in the context ofsolving the convection equation (Eq.
12.67) over the three-dimensional unstructured grid system shown in Fig. 12.16. The discretized equation for the element ofvolume VC shown in Fig. 12.16 can be written asXð12:77Þm_ f /f ¼ Q/C VCf nbðC ÞThe / values at cell faces are computed using a HR scheme and the objective ofthe various methods is to incorporate these values in the discretized equation in themost effective manner.46012 High Resolution Schemes12.8.1 The Downwind Weighing Factor (DWF) MethodThe Downwind Weighing Factor (DWF) [25] defined asDWFf ¼~ /~//f /CfC¼~/D /C1 /Cð12:78Þis used to rewrite the face value such that/f ¼ DWFf /D þ 1 DWFf /C ¼ /C þ DWFf ð/D /C Þð12:79Þ~ betweenthereby redistributing the HR scheme estimate /f or the normalized value /fthe Upwind and Downwind nodes.
The effect is a reduced stencil for the discretizedcoefficients. Since the value of /f computed using a HR scheme lies between /Cand /D ; the value of DWFf is always between 0 and 1, i.e., 0 DWFf 1:Now rather than computing the DWFf explicitly from the computed /f value,the DWFf can be expressed directly from the functional relationships of the HRscheme. Table 12.1 presents such relationships for several HO and HR schemes onuniform grid.Comparing the TVD formulation given by Eq. (12.30) with Eq. (12.79), it isclear that1 DWFf ¼ w rfð12:80Þ2As shown in the previous section on TVD schemes, the coefficients obtained fromsuch implementation will not be diagonally dominant and the formulation is thus notstable for many flow configurations.
As the method is used in OpenFOAM®, thereasons behind the numerical difficulties that are generally experienced by the codewhen solving convection dominated flow problems are now clear.For completeness, the analysis of the implementation via the DWF is presented.Starting with Eq. (12.79), the convection flux in the general case is written in the form ihm_ f /f ¼ m_ f ; 0 DWFfþ /F þ 1 DWFfþ /C ih m_ f ; 0 DWFf /C þ 1 DWFf /Fð12:81Þwith/f /C/F /C/f /FDWFf ¼/C /FDWFfþ ¼ð12:82Þ12.8The DWF and NWF Methods461Table 12.1 Functional relationships of the DWFf for some HO and HR schemes on uniform gridSchemeDownwind weighing factor-NVFUpwindSOUDWFf ¼ 0~/DWFf ¼ C ~2 1/C12CDDWFf ¼FROMM1DWFf ¼ ~4 1/CQUICKDWFf ¼DownwindMINMODDWFf ¼ 18~/1>C>~ 1 0/> C>>22~< 1 /CDWFf ¼11 ~>> /C 1>>22>:0elsewhere8<1~ 10/CDWFf ¼ 2:0 elsewhere8~>1/>C~ 2> 0/>C>3<2 1 /~CDWFf ¼2 ~>>> /C 11>>3:0elsewhere8~2/>C>~ 1>0/C>>~6>1/C>>>11 ~5<1þ /C DWFf ¼ 4 8 1 /66~>C>>>>5>~ 1>1/>C>6:0elsewhere8~2/>C>~ 1>0/C>~>51/>C>>>111>~> /C >>>252<DWFf ¼ 111 ~5> /C >>4 þ >26~>81/>C>>>>5 ~>> /C 11>>6:0elsewhere8~/>C>~ 1>0/>C> 1/~4>C>>>3< 1 1 ~ /C DWFf ¼ 4 1 /44~>C>>>>3 ~>> /C 11>>4:0elsewhere11þ 4 8 1/~CBounded CDOSHER [30]SMARTSTOICMUSCL46212 High Resolution SchemesThe fluxes in Eq.
(12.77) can now be expressed asFluxFf ¼ m_ f ; 0DWFfþ m_ f ; 0 1 DWFf FluxCf ¼ m_ f ; 0 1 DWFfþ m_ f ; 0DWFfð12:83Þand the discretized equation becomesaC /C þXaF / F ¼ b Cð12:84ÞF NBðC ÞwithaF ¼ FluxFfXaC ¼FluxCf ¼ f nbðC ÞXF NBðCÞXaF þm_ ff nbðC Þð12:85ÞbC ¼ Q/C VCThe coefficients in Eq. (12.85) result in a highly unstable system of equations,thus requiring substantial relaxation. This can be demonstrated on a simpleone-dimensional mesh (Fig. 12.8). Without loss of generality, a positive flow fieldis assumed, which reduces Eq. (12.84) toaC /C þ aE /E þ aW /W ¼ bCð12:86ÞwhereaE ¼ FluxFe ¼ m_ e DWFeþaW ¼ FluxFw ¼ m_ w 1 DWFwXaC ¼FluxCf ¼ m_ e DWFeþ þ m_ w 1 DWFw þ m_ e þ m_ wð12:87Þf nbðC ÞThe continuity constraint implies thatm_ e þ m_ w ¼ 0 ) m_ w ¼ m_ eð12:88Þthus the coefficients becomeaE ¼ FluxFe ¼ m_ e DWFeþaW ¼ FluxFw ¼ m_ e 1 DWFwXaC ¼FluxCf ¼ m_ e DWFeþ þ DWFw 1f nbðCÞð12:89Þ12.8The DWF and NWF Methods463In the above equation the aE and aW coefficients have opposite signs, a violation toone of the basic coefficient rules.
Furthermore for values of the DWF factors that arelarger than 0.5 (i.e., DWF [ 0:5), the diagonal coefficient aC becomes negativeresulting in a system not solvable by iterative means. This would occur whenever/f [ 0:5ð/C þ /D Þ a situation common to all HR schemes for /C [ 0:5: In fact theDWF moves much of the HR flux influence onto the downwind value causing theabove mentioned issues, a situation that resembles in effect the central differencescheme.12.8.2 The Normalized Weighing Factor (NWF) MethodThe Normalized Weighing Factor (NWF) method was developed [29] to addressthe shortcomings of the DWF method.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
