Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 40
Текст из файла (страница 40)
8.9 An element in anon-orthogonal mesh systemSft neCfFd CF2428Spatial Discretization: The Diffusion Termwith Ef being in the CF direction to enable writing part of the diffusion flux as afunction of the nodal values ϕF and ϕC such thatEf ez}|{ðr/Þf Ef|fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}ðr/Þf Sf ¼orthogonallike contributionþ @/¼ Efþðr/Þf Tf@e f¼ Efðr/Þf T|fflfflfflfflfflffl{zfflfflfflfflfflffl}nonorthogonal like contributionð8:67Þ/F /Cþ ðr/Þf TfdCFThe first term on the right hand side of Eq.
(8.67) represents a contributionsimilar to the contribution on orthogonal grids, i.e., involving ϕF and ϕC, while thesecond term on the right hand side is called cross-diffusion or non-orthogonaldiffusion [8] and is due to the non-orthogonality of the grid. Different options forthe decomposition of Sf are available and are discussed next.8.6.2Minimum Correction ApproachAs shown in Fig. 8.10, the decomposition of Sf is done in such a way as to keep thenon-orthogonal correction in Eq.
(8.67) as small as possible, by making Ef and Tforthogonal. As the non-orthogonality increases, the contribution to the diffusionflux from ϕF and ϕC decreases. In this case the vector Ef is computed asEf ¼ e Sf e ¼ Sf cos h eð8:68ÞFig. 8.10 Decomposing Sfvia the minimum correctionapproachSfneCfEfTfF8.6 Non-orthogonal Unstructured Grid8.6.3243Orthogonal Correction ApproachThis approach, schematically depicted in Fig. 8.11, keeps the contribution of theterm involving ϕF and ϕC the same as on an orthogonal mesh irrespective of thedegree of grid non-orthogonality.
To achieve this, Ef is defined asEf ¼ Sf eð8:69ÞFig. 8.11 Decomposing Sfvia the orthogonal correctionapproachSfTfFneEffC8.6.4Over-Relaxed ApproachIn this approach, the importance of the term involving ϕF and ϕC is forced toincrease as grid non-orthogonality increases. As shown in Fig. 8.12, this is achievedby selecting Tf to be normal to Sf. Mathematically Ef is computed asEf ¼Sfe¼cos h!S2fS f Sfe¼eSf cos he Sfð8:70ÞTo summarize, the diffusion flux at an element face of a non-orthogonal gridcannot be written solely in terms of the values at the nodes straddling the face.Fig. 8.12 Decomposing Sfvia the over-relaxed approachSfneCfTfEfF2448Spatial Discretization: The Diffusion TermA term that accounts for non-orthogonality has to be added. This term is denoted inthe literature by “cross diffusion” and is computed asðr/Þf Tf ¼ ðr/Þf Sf Ef8ðr/Þf ðn cos heÞSf>>>< ðr/Þ ðn eÞSff¼>>> ðr/Þ n 1 e Sf:fcos hminimum correctionnormal correctionð8:71ÞoverrelaxedFor orthogonal meshes n and e are collinear, the angle θ shown in Fig.
8.9 iszero, and the cross-diffusion term is zero. When cross diffusion is not zero, andsince it cannot be written as a function of ϕF and ϕC, it is added as a source term inthe element algebraic equation.All approaches described above are correct and satisfy Eq. (8.59). The differencebetween these methods is in their accuracy and stability on non-orthogonal meshes.The over-relaxed approach has been found to be the most stable even when the gridis highly non-orthogonal. Nevertheless, the final general form of the discretizeddiffusion term is the same for all three approximations.8.6.5Treatment of the Cross-Diffusion TermThe cross diffusion term cannot be expressed in terms of nodal values. Due to thisfact, it is treated in a deferred correction manner by computing its value using thecurrent gradient field and adding it as a source term on the right hand side of thealgebraic equation.
Gradients are computed at the main grid points, as describednext, and values at the interfaces are obtained by interpolation.8.6.6Gradient ComputationIn discretizing the diffusion term over a one-dimensional or orthogonalmulti-dimensional computational domain, it was shown that the gradient of ϕ can beexplicitly written as a function of the cell nodal ϕ values. In non-orthogonaldomains however, the computation of the diffusion flux was found to be morecomplex. The non-orthogonal component of the gradient could not be linearizedand written as a function of nodal values, but rather had to be moved to the righthand side and evaluated explicitly. This means that the gradient has to be evaluatedin order to incorporate its non-orthogonal contribution in the discretization equation. One widely used approach for computing the gradient at a cell is the8.6 Non-orthogonal Unstructured Grid245Green-Gauss or gradient theorem, which states that for any closed volume V,surrounded by a surface ∂V the following holds:IZr/ dV ¼/ dSð8:72Þ@VVwhere dS is the outward pointing incremental surface vector.
In order to obtain adiscrete version of this equation, the mean value theorem is applied according towhich the integral on the left hand side of Eq. (8.72) and the average gradient overthe volume V are related byZr/V ¼ r/ dVð8:73ÞVCombining Eqs. (8.72) and (8.73), the average gradient over element C shown inFig. 8.9 is found to be1r/C ¼VCI/f Sfð8:74Þ@VCNext the integral over a cell face is approximated by the face centroid valuetimes the face area. Thus r/C , or simply ∇ϕC, is computed asr/C ¼1 X/ SfVC f nbðCÞ fð8:75ÞThe gradient at the face of an element can then be obtained as the weightedaverage of the gradients at the centroids straddling the interface and is given byr/f ¼ gC r/C þ gF r/Fð8:76Þwhere gC and gF are geometric interpolation factors related to the position of theelement face f with respect to the nodes C and F.8.6.7Algebraic Equation for Non-orthogonal MeshesSplitting the surface vector Sf into the two Ef and Tf vectors and substituting itsequivalent expression into the semi-discretized equation of the diffusion fluxes, yield2468X Jf/;D Sf¼f nbðCÞSpatial Discretization: The Diffusion TermX C/ r/ f Ef þ Tff nbðCÞ¼X C/ r/ f Ef þX C/ r/ f Tff nbðCÞf nbðCÞ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}Orthogonal Linearizeable PartNonOrthogonal NonLinearizeable PartX / ð/ F / C ÞX Cf Ef C/ r/ f TfþdCFf nbðCÞf nbðCÞX X /Cf gDifff ð/C /F Þ þ C/ r/ f Tf¼¼f nbðCÞ0X¼@f nbðCÞ1FluxCf A/C þf nbðCÞX X FluxFf /F þFluxVff nbðCÞf nbðCÞð8:77Þwhere again gDifff is a geometric diffusion coefficient defined asgDifff ¼EfdCFð8:78ÞUsing the above form of the diffusion fluxes and expanding, the final form of thediscretized diffusion equation over unstructured/structured non-orthogonal grid isobtained asXaF /F ¼ bCð8:79ÞaC /C þFNBðCÞwhereaF ¼ FluxFf ¼ C/f gDifffXXX /aC ¼FluxCf ¼ FluxFf ¼Cf gDiffff nbðCÞbC ¼Q/C VCX f nbðCÞf nbðCÞFluxVf ¼f nbðCÞQ/C VCþX f nbðCÞC r//f Tfð8:80ÞNote the change in sign of the non-orthogonal term on the right hand side of theequation.Example 2For the polygonal element C and its neighbor F shown in Fig.
8.13, thesolution at any point satisfies ϕ = x2 + y2 + x2y2. If the volume of cell C isVC = 8.625 calculate the following:8.6 Non-orthogonal Unstructured Grid2471. The gradient of ϕ at (i.e., ∇ϕC) both numerically and analytically.2. The analytical value of ∇ϕF.3. Interpolate between the numerical value of ∇ϕC and the analytical valueof ∇ϕF to find an approximate value for r/f1 . Compare with the analytical value of r/f1 .4. Express r/f1 · Sf1 in terms of ϕC and ϕF using(a) The minimum correction approach(b) The orthogonal correction approach(c) The over-relaxed approachySf(0.5,3)Sf2F ( 4.25,3.5)( 2.5,4)f2e13f1n1f3f5f4Sff5 ( 2.75,0.75)Sf4( 2,0)1(3.5,1.5)C (1.75,2 )(0,1)Sf5xFig.
8.13 A polygonal element with its geometric entitiesSolution1. Knowing that ϕ = x2 + y2 + x2y2 at any point in the domain, the value of thegradient at any point can be calculated asr/ ¼ @/@/iþj ¼ 2x þ 2xy2 i þ 2y þ 2yx2 j@x@y2488Spatial Discretization: The Diffusion TermTherefore the analytical value of ∇ϕC is given byr/C ¼ 17:5i þ 16:25jThe numerical value of ∇ϕC can be computed using1 Xr/C ¼/ SfVC f nbðCÞ fTherefore the values of ϕf and Sf are required.
Using the given analyticalexpression, the values of ϕ at the various locations are found as/C ¼ 1:752 þ 22 þ 1:752 22 ¼ 19:3125/F ¼ 4:252 þ 3:52 þ 4:252 3:52 ¼ 251:578125/f1 ¼ 32 þ 2:752 þ 32 2:752 ¼ 84:625/f2 ¼ 1:52 þ 3:52 þ 1:52 3:52 ¼ 42:0625/f3 ¼ 0:252 þ 22 þ 0:252 22 ¼ 4:3125/f4 ¼ 12 þ 0:52 þ 12 0:52 ¼ 1:5/f5 ¼ 2:752 þ 0:752 þ 2:752 0:752 ¼ 12:37890625The surface vector is given bySf ¼ ðDyi DxjÞwhere the correct sign is chosen such that the vector is pointing outward.
Thisis done by computing the surface vector as Sf = Δyi − Δxj and then performingthe dot product of Sf with the distance vector dCF (pointing from C to F). If theproduct is positive then the direction of Sf is correct. If the product is negativethen the sign of Sf should be reversed, as shown below.The distance vectors are computed asdCf1 ¼ ð3 1:75Þi þ ð2:75 2Þj ¼ 1:25i þ 0:75jdCf2 ¼ ð1:5 1:75Þi þ ð3:5 2Þj ¼ 0:25i þ 1:5jdCf3 ¼ ð0:25 1:75Þi þ ð2 2Þj ¼ 1:5idCf4 ¼ ð1 1:75Þi þ ð0:5 2Þj ¼ 0:75i 1:5j8.6 Non-orthogonal Unstructured Grid249dCf5 ¼ ð2:75 1:75Þi þ ð0:75 2Þj ¼ i 1:25jwhile the tentative surface vectors are given bySf1 ¼ ð4 1:5Þi ð2:5 3:5Þj ¼ 2:5i þ jSf2 ¼ ð4 3Þi ð2:5 0:5Þj ¼ i 2jSf3 ¼ ð3 1Þi ð0:5 0Þj ¼ 2i 0:5jSf4 ¼ ð1 0Þi ð0 2Þj ¼ i þ 2jSf5 ¼ ð1:5 0Þi ð3:5 2Þj ¼ 1:5i 1:5jPerforming the dot products, the obtained values areSf1 dCf1 ¼ ð2:5i þ jÞ ð1:25i þ 0:75jÞ ¼ 3:875 [ 0Sf2 dCf2 ¼ ði 2jÞ ð0:25i þ 1:5jÞ ¼ 3:25\0Sf3 dCf3 ¼ ð2i 0:5jÞ ð1:5iÞ ¼ 3\0Sf4 dCf4 ¼ ði þ 2jÞ ð0:75i 1:5jÞ ¼ 3:75\0Sf5 dCf5 ¼ ð1:5i 1:5jÞ ði 1:25jÞ ¼ 3:375 [ 0Therefore the correct values of the surface vectors should beSf1 ¼ 2:5i þ j Sf2 ¼ i þ 2j Sf3 ¼ 2i þ 0:5jSf4 ¼ i 2jSf5 ¼ 1:5i 1:5jThus, the numerical value of the gradient at C is obtained asr/C ¼84:625ð2:5i þ jÞ þ 42:0625ði þ 2jÞ þ 4:3125ð2i þ 0:5jÞ18:625þ1:5ði 2jÞ þ 12:37890625ð1:5i 1:5jÞ¼ 20:63111i þ 17:31454jwhich is close to the analytical value.2.