Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 71
Текст из файла (страница 71)
It operates by linearizing the normalizedinterpolation profile such that~ ¼ ‘/~ þk/fCð12:90Þwhere ‘ and k are constants that represent the slope and intercept of the linearfunction within any interval of /f ; with the number of intervals depending on theHR scheme used. This is an exact representation for nearly all HR schemes.For example, by equating Eq. (12.90) to the NVF form of the MINMOD scheme(Eq. 12.14), the values of ‘ and k are deduced to be83>>;0>>< 2 ½‘; k ¼1 1>;>>>: 2 2½1; 0~ \0\/C121 ~ /C \12elsewhereð12:91ÞIn a second step Eq.
(12.90) is rewritten as/f /U/ /U¼‘ Cþk/D /U/D /Uð12:92Þand transformed to/f ¼ ‘ð/C /U Þ þ kð/D /U Þ þ /U ¼ ‘/C þ k/D þ ð1 ‘ kÞ/U ð12:93Þwhere /U ; /D ; and /C are the values at the U, D, and C nodes whose locationsdepend on the flow direction. The values of ‘ and k for a number of HO and HRschemes are listed in Table 12.2 (for unstructured and/or structured uniform grid).Since for an unstructured grid the U location is virtual, the term involving /U istreated in a deferred correction fashion. However the value of the resulting deferred46412 High Resolution SchemesTable 12.2 NVF values of NWF ½‘; k factors for some HO and HR schemesSchemeUniform grid (NVF)UpwindSOU½‘; k ¼ ½1; 03½‘; k ¼ ; 021 1½‘; k ¼ ;2 2 1½‘; k ¼ 1;43 3½‘; k ¼ ;4 883>>;0>>< 2 ½‘; k ¼1 1>;>> 2 2>:½1; 083>>;0>< 2½‘; k ¼>½0; 1>>:½1; 08>>½2; 0>>> >>1<1;½‘; k ¼4>>>>>½0;1>>:½1; 08>>½4; 0>>>>>< 3 3;½‘; k ¼4 8>>>>>½0; 1>>:½1; 0CDFROMMQUICKMINMODOSHERMUSCLSMART~ \0\/C121 ~ /C \12elsewhere~ \20\/C32 ~ /C \13elsewhere~ \10\/C41 ~3 /C \443 ~ /C \14elsewhere~ \10\/C61 ~5 /C \665 ~ /C \16elsewherecorrection source term ði.e., ð1 ‘ kÞ/U l 0; k 0Þ is smaller than the one thatwould be obtained with the standard deferred correction treatment (i.e., /U ).
Thusthe NWF requires less under-relaxation than the standard DC method and thusallows for faster convergence.Starting with Eq. (12.93), the convection flux in the general case (i.e., multidimensional unstructured grid) can be written in the formh iþþþþm_ f /f ¼k m_ f ; 0 k ‘þ/þk/þ1‘kCFffff /Uh i k m_ f ; 0 k ‘f /F þ kf /C þ 1 ‘f kf /Uð12:94Þ12.8The DWF and NWF Methods465and linearized to yieldFluxFf ¼k m_ f ; 0 k kfþ km_ f ; 0 k ‘f_ f ; 0 k kfFluxCf ¼k m_ f ; 0 k ‘þf kmþþ_FluxVf ¼k m_ f ; 0 k 1 ‘þkkm;0k1‘k/fUffff /Uð12:95ÞSubstituting Eq. (12.95) into Eq.
(12.77) the algebraic equation is obtained asXaC /C þaF / F ¼ b Cð12:96ÞF NBðC Þwhere now _ f ; 0aF ¼ FluxFf ¼ kfþ m_ f ; 0 ‘f mXX _ f ; 0 k m_ f ; 0aC ¼FluxCf ¼‘þf mff nbðC Þf nbðCÞXbC ¼ Q/C VC FluxVfð12:97Þf nbðCÞ|fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}bDCCX h¼ Q/C VC iþþ_ f ; 0_ f ; 0 1 ‘1 ‘þf kf /U mf kf /U mf nbðCÞ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}bDCCThe NWF was originally developed for use on structured grids with its formulation in that context allowing for a full implicit treatment of HR schemes. Thefull implicitness of the method on structured grid is the result of /U being an actualnode in the computational domain that can be resolved in the algebraic equation,which now has a larger stencil that includes the far nodes EE and WW.
For the onedimensional structured grid displayed in Fig. 12.8, the NWF form of the algebraicequation becomesXaC / C þaF / F ¼ bCð12:98ÞF E;W;EE;WWwhereþ_ w ; 0k 1 ‘þaE ¼ FluxFe ¼ km_ e ; 0kkeþ km_ e ; 0k‘e þ kmw kwþ_ e ; 0k 1 ‘ þaW ¼ FluxFw ¼ km_ w ; 0kkwþ km_ w ; 0k‘w þ kme keaEE ¼ FluxFee ¼ km_ e ; 0k 1 ‘e keð12:99ÞaWW ¼ FluxFww ¼ km_ w ; 0k 1 ‘w kwX_ w ; 0k‘þ_ e ; 0kke km_ w ; 0kkwFluxCf ¼ km_ e ; 0k‘þaC ¼e þ kmw kmf nbðC Þ¼ ðaE þ aW þ aEE þ aWW Þ þ ðm_ e þ m_ w Þ46612 High Resolution SchemesIn the NWF reformulation of the HR schemes since the value of ‘ is greater thanthat of k (see Fig. 12.6), except in a narrow region of the NVD close to the Downwindline as explained next, the value of aC is always positive and instability does notarise. Along the Downwind line of the NVD, where ð‘; k Þ ¼of zero ð0; 1Þ; a valuefor the aC coefficient is obtained.
In this case ð‘; kÞ is set to L; 1 L/f where L isusually set to the value of ‘ in the previous interval of the composite scheme.This basically allows the NWF to be much more robust than the DWF as itguarantees positive aC coefficients.12.8.2.1The NWF Method in the Context of the TVDWith the exception of the MUSCL Van Leer limiter, the limiters of the TVDformulation of all HR schemes presented earlier appear as straight lines on Sweby’sdiagram and as such can be written as a set of linear equations of the form w rf ¼ mrf þ nð12:100Þwhere m and n are constants (slope and intercept of the linear function and dependon geometric quantities only) within any interval of w rf ; with the number ofintervals depending on the HR TVD scheme used.
For example, by equatingEq. (12.100) to the TVD form of the MINMOD scheme (Eq. 12.44), the values ofm and n are deduced to be8< ½1; 0 0\rf \1MINMOD ½m; n ¼ ½0; 1 rf 1ð12:101Þ:½0; 0 rf 0Substituting Eq. (12.100) in Eq. (12.30), the interface value is found to be1mrf þ n ð/D /C Þ21/C /Um¼ /C þþ n ð/D /C Þ2/D /C1111¼ 1 þ m n /C þ n/D m/U2222/f ¼ /C þð12:102Þwhere /U ; /D ; and /C are again the values at the U, D, and C nodes whoselocations depend on the flow direction. The values of m and n for a number of HOand HR schemes are listed in Table 12.3 for uniform grids.
Moreover Eq. (12.102)has the same form as Eq. (12.93) with11‘ ¼ 1 þ m n and221k¼ n2ð12:103ÞThus an approach similar to that of the NVF-NWF can be used for the implementation of the TVD-NWF.12.9Boundary Conditions467Table 12.3 TVD values of NWF ½m; n factors for some HO and HR schemesSchemeUniform grid (NVF)UpwindSOUCDFROMM½m; n ¼ ½0; 0½m; n ¼ ½1; 0½m; n ¼ ½0; 11 1½m; n ¼ ;2 21 3½m; n ¼ ;4 4½m; n ¼ ½0; 28< ½1; 0½m; n ¼ ½0; 2:½0; 08>½2; 0>>>>< 1 1;½m; n ¼2 2>>>>> ½0; 2:½0; 08>>> ½2; 0>>>><½0; 1½m; n ¼>½1; 0>>>>>½0; 2>:½0; 0QUICKDOWNWINDOSHERMUSCLSUPERBEE12.90\rf \2rf 2r00\rf \131 rf \33rf 3rf 010\rf \21 rf \121 rf \2rf 2rf 0Boundary ConditionsThe Boundary conditions for the convection term are generally much simpler thanfor the diffusion term.
The particulars of the implementation of the followingboundary conditions: “Inlet”, “Outlet”, “Wall”, and “Symmetry” are now detailed.Typical boundary elements are shown in Fig. 12.17. A boundary cell, as mentioned earlier in this book, is one that has one or more faces on the boundary. Discretevalues of / are stored both at centroids of boundary cells and of boundary faces.Let C denotes the centroid of the boundary element with one boundary face ofcentroid b and of surface vector Sb pointing outward (Fig. 12.17).
As before, thediscretization process over cell C of a pure convection problem in a multidimensional domain yieldsX J/;C S f ¼ 0ð12:104Þf nbðC ÞThe fluxes on the interior faces are discretized as before. Independent of theboundary condition type, the boundary flux Jb/;C may be written using the boundaryface centroid value as46812 High Resolution SchemesFig. 12.17 Boundaryelements with one or twoboundary facesSbCbJb/;C ¼ ðqv/Þbð12:105ÞJb/;C Sb ¼ ðqv/Þb Sb ¼ m_ b /bð12:106Þsuch thatThus the discretized equation of the boundary cell is expressed asXðqv/ SÞf þðqv/ SÞb ¼ 0ð12:107Þf nbðC Þwhere subscript f refers to interior faces and subscript b to the boundary face.
Thespecification of boundary conditions involves either specifying the unknownboundary value /b ; or alternatively, the boundary flux Jb/;C : Using Eq. (12.107), thediscretized equations at a boundary element for the different boundary conditiontypes of convection problems are derived next.12.9.1 Inlet Boundary ConditionAt inlet to a domain (Fig. 12.18), the value of / is usually specified. Since thevelocity field is assumed to be known, then the convective flux at inlet is also known.Therefore the boundary flux is moved to the right hand side and treated as a sourceterm. With this modification Eq. (12.107) becomes12.9Boundary Conditions469Fig. 12.18 Inlet boundarycondition for the convectionfluxCebbnSbXðqv SÞf /f ¼ ðqv SÞb /b ¼ m_ b /bð12:108Þf nbðCÞIf a HR scheme is used to discretize the convection flux at interior faces and isimplemented via a deferred correction approach, then the modified algebraicequation for the boundary element can be written asXaC /C þaF / F ¼ b Cð12:109ÞF NBðC ÞwhereaF ¼ FluxFf ¼ km_ f ; 0 kXXFluxCf ¼km_ f ; 0 kaC ¼f nbðC Þ¼f nbðC ÞXaF þF NBðC ÞbC ¼ Xf nbðC ÞXm_ ff nbðCÞFluxVf ¼ m_ b /b ð12:110ÞXUm_ f /HRf /ff nbðCÞ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}bDCcwhere F refers to interior neighboring nodes of the C grid point and f refers tointerior faces of the boundary element.47012 High Resolution Schemes12.9.2 Outlet Boundary ConditionAt outlet from the domain (Fig.
12.19) no information downstream of the boundarygrid point is available. However, being a directional phenomenon, the value of / atthe boundary is highly dependent on upstream locations. In fact, the upwind andSOU schemes, for example, do not require any information at outlet since its valuecan be expressed as a function of values at upstream nodes. The treatment that hasproven to be very effective at an outlet boundary condition is to assume the /profile to be fully developed, which is equivalent to assuming that the normalgradient to the face is zero [i.e., ðr/ nÞb ¼ ð@/=@nÞb ¼ 0].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
