Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 41
Текст из файла (страница 41)
The analytical value of ∇ϕF is easily calculated asr/F ¼ 112:625i þ 133:4375j2508Spatial Discretization: The Diffusion Term3. The interpolated value of r/f1 is obtained usingr/f1 ¼ gf1 r/F þ 1 gf1 r/Cwheregf1 ¼dCf1dCf1 þ df1 FPerforming the calculations, the interpolation factor is found to be given byqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffixf1 xC þ yf1 yC ¼ ð3 1:75Þ2 þð2:75 2Þ2 ¼ 1:4577qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 2ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi22dfF ¼xF xf1 þ yF yf1 ¼ ð4:25 3Þ þð3:5 2:75Þ ¼ 1:4577dCf1 ¼gf 1 ¼1:4577¼ 0:51:4577 þ 1:4577leading to the following value for r/f1 :r/f1 ¼ 66:628055i þ 75:37602jThe analytical value of r/f1 is easily obtained asr/f1 ¼ 51:375i þ 55jAgain values are close.4. The general form of r/f1 Sf1 is given byð/ /F Þr/f1 Sf1 ¼ Ef1 Cþ r/f1 Tf1dCF|fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} nonorthogonallikeorthogonallikecontributioncontributionwithdCF ¼ ð4:25 1:75Þi þ ð3:25 2Þj ¼ 2:5i þ 1:5j ) dCF ¼¼ 2:9155pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2:52 þ 1:52The unit vector e1 in the direction of dCF is calculated usinge1 ¼dCF¼ 0:8575i þ 0:5145jdCF8.6 Non-orthogonal Unstructured Grid(a) The minimum correction approachUsing this approach, the expression for Ef1 is computed asEf1 ¼ e1 Sf1 e1 ¼ 2:279i þ 1:368j ) Ef1 ¼ 2:658The normal component is computed asTf1 ¼ Sf1 Ef1 ¼ 0:221i 0:368j ) r/f1 Tf1 ¼ 13:014The diffusion flux at the face becomesr/f1 Sf1 ¼ 0:9122ð/C /F Þ þ 13:014(b) The orthogonal correction approachIn this approach, the expression for Ef1 is computed asEf1 ¼ Sf1 e1 ¼ 2:309i þ 1:385j ) Ef1 ¼ 2:693The normal component is found to beTf1 ¼ 0:191i 0:385j ) r/f1 Tf1 ¼ 16:294The diffusion flux at the face becomesr/f1 Sf1 ¼ 0:924ð/C /F Þ þ 16:294(c) The over-relaxed approachThe expression for Ef1 is computed asEf1 ¼Sf1 Sf1e1 ¼ 2:339i þ 1:403j ) Ef1 ¼ 2:728e 1 Sf 1The normal component is found to beTf1 ¼ 0:161i 0:403j ) r/f1 Tf1 ¼ 19:649The diffusion flux at the face becomesr/f1 Sf1 ¼ FluxCf1 /C þ FluxFf1 /F þ FluxVf1¼ 0:936ð/C /F Þ þ 19:6492512528.6.88Spatial Discretization: The Diffusion TermBoundary Conditions for Non-orthogonal GridsThe treatment of boundary conditions for non-orthogonal grids is similar to that fororthogonal grids with some minor differences related to the non-orthogonal diffusion contribution.
This is outlined next.8.6.8.1Dirichlet Boundary ConditionFor the case of a Dirichlet boundary condition, i.e., when ϕ is specified by the userat the boundary as shown in Fig. 8.14, the boundary discretization proceeds as inorthogonal-grids. However, there is a need now to account for the cross-diffusion,which arises on boundary faces as on interior faces.
This happens whenever thesurface vector is not collinear with the vector joining the centroids of the elementand boundary face. The diffusion flux along the boundary face is discretized as Sb ¼ C/b ðr/Þb Sb ¼ C/b ðr/Þb ðEb þ Tb ÞJ/;Db/ /b /CEb C/b ðr/Þb Tb¼ CbdCbð8:81Þ¼ FluxCb /C þ FluxVbwhereFluxCb ¼ C/b gDiffbFluxVb ¼ C/b gDiffb /b C/b ðr/Þb TbEbgDiffb ¼dCbð8:82ÞSubstituting into Eq. (8.79), the modified coefficients are obtained asXaF ¼ FluxFf aC ¼FluxCf þ FluxCbf nbðCÞbC ¼Q/C VC FluxVb Xð8:83ÞFluxVff nbðCÞFig.
8.14 Dirichlet boundaryfor a non-orthogonal meshSfneCd Cbb8.6 Non-orthogonal Unstructured Grid8.6.8.2253Neumann Boundary ConditionThe Neumann type condition for non-orthogonal grids follows that for orthogonalgrids. In this case the user-specified flux at the boundary is just added as a sourceterm as with orthogonal grid. The algebraic equation at the boundary is given byEq.
(8.40) with the modified coefficients specified by Eq. (8.41).8.6.8.3Mixed Boundary ConditionFor the mixed boundary condition case (Fig. 8.5), denoting the convection transfercoefficient by h∞ and the surrounding value of ϕ by /1 , the diffusion flux at theboundary can be written as/ /b /CEb C/b ðr/Þb TbS¼CJ/;DbbbdCbð8:84Þ¼ h1 ð/1 /b ÞSbfrom which an equation for ϕb is obtained as/b ¼h1 Sb / 1 þC/b Eb/ C/b ðr/Þb TbdCb CC/ Ebh1 Sb þ bdCbSubstituting ϕb back in Eq. (8.84), the flux equation is transformed to32C/b Eb76 h1 Sbh1 Sb C/b ðr/Þb TbdCb 76ð//ÞJb/;D Sb ¼ 671C4C/ Eb 5C/ Ebh1 Sb þ bh1 Sb þ bdCbdCb|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}abð8:85Þð8:86ÞS/;CDb¼ FluxCb /C þ FluxVbwhere nowC/b EbdCbFluxCb ¼C/ Ebh1 Sb þ bdCbh1 SbFluxVb ¼ FluxCb /1 ð8:87Þh1 Sb C/b ðr/Þb TbC/ Ebh1 Sb þ bdCb2548Spatial Discretization: The Diffusion Termand the modified coefficients for the boundary element are obtained asEfaF ¼ FluxFf C/fdCfXaC ¼ FluxCb þFluxCff nbðCÞbC ¼ Q/C VC FluxVb Xð8:88ÞFluxVff nbðCÞ8.7SkewnessTo evaluate many of the terms constituting the general discretized equation of aquantity ϕ, it is necessary to estimate its value at element faces.
An estimated valueat the face should be the average value of the entire face. At different steps of thediscretization process, linear variation of variables between nodes is assumed. Ifthis is extended to variation of variables along the face, then the average value ofany variable ϕ should be found at the face centroid. The common practice is to use alinear interpolation profile and to estimate the value at the face at the intersectionbetween the face and the line connecting the two nodes straddling the face.
Whenthe grid is skewed the line does not necessarily pass through the centroid of the face[9, 10]. An example is depicted in Fig. 8.15 where the intersection point of segment[CF] with the face is at point, f ′, which does not coincide with the face centroid, f.To keep the overall accuracy of the discretization method second order, all faceintegrations need to take place at point f. Thus a correction for the interpolated valueat f′ is needed in order to get the value at f.The skewness correction is derived by expressing the value of ϕ at f in terms ofits value and the value of its derivative at f ′ via a Taylor expansion such that/f ¼ /f 0 þ ðr/Þf 0 df 0 fð8:89Þwhere df ′f is a vector from the intersection point f ′ to the face centre f.Fig.
8.15 Non-conjunctionalelementsSfFffeC8.8 Anisotropic Diffusion8.8255Anisotropic DiffusionThe diffusion equation presented so far assumed the material has no preferreddirection for transfer of ϕ with the same diffusion coefficient in all directions, i.e.,the medium was assumed to be isotropic.
For the case when the diffusion coefficientof the medium is direction dependent, diffusion is said to be anisotropic [11–16]. Asmentioned in Chap. 3, some solids are anisotropic for which the semi-discretizeddiffusion equation becomesX j/ r/ f Sf ¼ S/C VCð8:90Þf nbðCÞwhere κϕ is a second order symmetric tensor. Assuming a general three dimensionalsituation, the term on the left hand side, through some mathematical manipulation[17], can be rewritten as2j/ r/fj/11j/126 / S f ¼ 64 j21j/22j/31j/3223@/7j/13 66 @x 7767@/67j/23 75 6 @y 7 Sf67j/33 f 4 @/ 53@zð8:91ÞfPerforming matrix multiplication, Eq. (8.91) becomes3@// @// @/þ j12þ j132 3@x@y@z 77 Sx6766 / @/6 7 // @// @/ 77 6 y7j r/ f Sf ¼ 66 j21 @x þ j22 @y þ j23 @z 7 4 S 5767 Sz654 / @/f/ @// @/þ j32þ j33j31@x@y@z f2/6 j11ð8:92ÞFurther manipulations yieldj/ r/f2j/11j/21j/3132Sx36 /76 7// 7 6 y76j4 12 j22 j32 5 4 S 5fj/13 j/23 j/33 f Sz fh T i¼ ðr/Þf j/ S ¼ ðr/Þf S0f Sf ¼ @/@x@/@y@/@zð8:93ÞfSubstituting Eq.
(8.93) into Eq. (8.90), the new form of the diffusion equationbecomes2568XSpatial Discretization: The Diffusion Termðr/Þf S0f ¼ Q/C VCð8:94Þf nbðCÞIt is obvious that in this form the discretization procedure described above can beapplied by simply setting C/ to 1 and replacing Sf by S0f . Therefore the same codecan be used to solve isotropic and anisotropic diffusion problems.8.9Under-Relaxation of the Iterative Solution ProcessFor a general diffusion problem, C/ may be a function of the unknown dependentvariable ϕ and the grid may be highly non-orthogonal with a large cross diffusionterm, which is treated using a deferred correction approach.
Therefore large variations in ϕ between iterations result in large source terms and large changes in thecoefficients, which may cause divergence of the iterative solution procedure. Thisdivergence is usually due to the non-linearity introduced by the coefficients and thecross-diffusion term, which makes the source term highly affected by the currentsolution field which is not yet converged.