Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 39
Текст из файла (страница 39)
Themain coefficient and source term are calculated asaC ¼ FLuxCe þ FLuxCw þ FLuxCn þ FLuxCs¼ 40 þ 0:002 þ 200 þ 200 ¼ 440:002bC ¼ 0Substituting, the discretized algebraic equation is obtained as440:002T5 40T6 0:002T4 200T8 200T2 ¼ 0Element #6The needed diffusion terms are calculated asgDiffw ¼Dy60:1Dx60:3Dx60:3¼¼¼¼ 0:4 gDiffn ¼¼ 3 gDiffs ¼¼3dx56 0:25dy69 0:1dy36 0:1The interface conductivities arekw ¼ k56 ¼ 102kn ¼ k69 ¼ 102ks ¼ k36 ¼ 102The general form of the equation is written asaC T6 þ aW T5 þ aN T9 þ aS T3 ¼ bCwhereaW ¼ FluxFw ¼ kw gDiffw ¼ 102 0:4 ¼ 40aN ¼ FluxFn ¼ kn gDiffn ¼ 102 3 ¼ 300aS ¼ FluxFs ¼ ks gDiffs ¼ 102 3 ¼ 300The east coefficient does not appear in the equation as its represents a zeroflux boundary condition.
Moreover the main coefficient and source term aregiven byaC ¼ FLuxCw þ FLuxCn þ FLuxCs¼ 40 þ 300 þ 300 ¼ 640bC ¼ FLuxVe ¼ 0Substituting, the discretized algebraic equation is obtained as640T6 40T5 300T9 300T3 ¼ 02348Spatial Discretization: The Diffusion TermElement #7The needed diffusion terms are calculated asgDiffe ¼Dy70:1¼¼ 0:667dx78 0:15gDiffw ¼Dy70:1¼¼2dx197 0:05gDiffs ¼Dx70:1¼¼1dy47 0:1The interface conductivities areke ¼ k78 ¼ 3 103kw ¼ k19 ¼ 103kn ¼ k18 ¼ 103ks ¼ k47 ¼ 103The general form of the equation is written asaC T7 þ aE T8 þ aS T4 ¼ bCwhereaE ¼ FluxFe ¼ ke gDiffe ¼ 3 103 0:667 ¼ 0:002aS ¼ FluxFs ¼ ks gDiffs ¼ 103 1 ¼ 0:001The west and north coefficients do not appear in the equation as theirinfluence is integrated through the boundary conditions asReq ¼h1 ðk18 =dy718 Þ20ð103 =0:05ÞDx7 ¼0:1 ¼ 0:001998h1 þ ðk18 =dy718 Þ20 þ 103 =0:05FluxCw ¼ kw gDiffw ¼ 103 2 ¼ 0:002 FluxVw ¼ kw gDiffw T19 ¼ 0:64FluxVn ¼ Req T1 ¼ 0:5994FluxCn ¼ Req ¼ 0:001998The main coefficient and source term can now be calculated and are given byaC ¼ FLuxCe þ FLuxCw þ FLuxCn þ FLuxCs¼ 0:002 þ 0:002 þ 0:001998 þ 0:001 ¼ 0:006998bC ¼ FLuxVw FLuxVn ¼ 0:64 þ 0:5994 ¼ 1:2394Substituting, the discretized algebraic equation is obtained as0:006998T7 0:002T8 0:001T4 ¼ 1:2394Element #8The needed diffusion terms are calculated asgDiffe ¼Dy80:1Dy80:1¼¼¼ 0:4 gDiffw ¼¼ 0:667dx89 0:25dx78 0:15gDiffs ¼Dx80:2¼¼2dy58 0:18.4 The Interface Diffusivity235The interface conductivities arekw ¼ k78 ¼ 3 103ke ¼ k89 ¼ 102kn ¼ k17 ¼ 102ks ¼ k58 ¼ 102The general form of the equation is written asaC T8 þ aE T9 þ aW T7 þ aS T5 ¼ bCwhereaE ¼ FluxFe ¼ ke gDiffe ¼ 102 0:4 ¼ 40aW ¼ FluxFw ¼ kw gDiffw ¼ 3 103 0:667 ¼ 0:002aS ¼ FluxFs ¼ ks gDiffs ¼ 102 2 ¼ 200The north coefficient does not appear in the equation as its influence isintegrated through the boundary condition asReq ¼h1 ðk17 =dy817 Þ20ð102 =0:05ÞDx8 ¼0:2 ¼ 3:9604h1 þ ðk17 =dy817 Þ20 þ 102 =0:05FluxCn ¼ Req ¼ 3:9604FluxVn ¼ Req T1 ¼ 1188:12The main coefficient and source term can now be calculated and are givenbyaC ¼ FLuxCe þ FLuxCw þ FLuxCn þ FLuxCs¼ 40 þ 0:002 þ 3:9604 þ 200 ¼ 243:9624bC ¼ FluxVn ¼ 1188:12Substituting, the discretized algebraic equation is obtained as243:9624T8 40T9 0:002T7 200T5 ¼ 1188:12Element #9The needed diffusion terms are calculated asgDiffw ¼Dy90:1¼ 0:4¼dx89 0:25gDiffs ¼Dx90:3¼3¼dy69 0:12368Spatial Discretization: The Diffusion TermThe interface conductivities arekw ¼ k89 ¼ 102kn ¼ k16 ¼ 102ks ¼ k69 ¼ 102The general form of the equation is written asaC T9 þ aW T8 þ aS T6 ¼ bCwhereaW ¼ FluxFw ¼ kw gDiffw ¼ 102 0:4 ¼ 40aS ¼ FluxFs ¼ ks gDiffs ¼ 102 3 ¼ 300The east and north coefficients do not appear in the equation as theirinfluence is integrated through the boundary conditions asReq ¼h1 ðk16 =dy916 Þ20ð102 =0:05ÞDx9 ¼0:3 ¼ 5:9406h1 þ ðk16 =dy916 Þ20 þ 102 =0:05FluxCn ¼ Req ¼ 5:9406FluxVn ¼ Req T1 ¼ 1782:18FluxCe ¼ FluxVe ¼ 0The main coefficient and source term can now be calculated and are givenbyaC ¼ FLuxCw þ FLuxCn þ FLuxCs ¼ 40 þ 5:9406 þ 300 ¼ 345:9406bC ¼ FLuxVn ¼ 1782:18Substituting, the discretized algebraic equation is obtained as345:9406T9 40T8 300T6 ¼ 1782:18Summary of equationsThe discretized algebraic equations are given by0:005T1 0:002T2 0:001T4 ¼ 0:64240:002T2 40T3 0:002T1 200T5 ¼ 208.4 The Interface Diffusivity237340T3 40T2 300T6 ¼ 300:006T4 0:002T5 0:001T7 0:001T1 ¼ 0:64440:002T5 40T6 0:002T4 200T8 200T2 ¼ 0640T6 40T5 300T9 300T3 ¼ 00:006998T7 0:002T8 0:001T4 ¼ 1:2394243:9624T8 40T9 0:002T7 200T5 ¼ 1188:12345:9406T9 40T8 300T6 ¼ 1782:18Solution of EquationsTo solve the above system of equations via the Gauss-Seidel method, theequations are rearranged intoT1 ¼ ð0:002T2 þ 0:001T4 þ 0:64Þ=0:005T2 ¼ ð40T3 þ 0:002T1 þ 200T5 þ 20Þ=240:002T3 ¼ ð40T2 þ 300T6 þ 30Þ=340T4 ¼ ð0:002T5 þ 0:001T7 þ 0:001T1 þ 0:64Þ=0:006T5 ¼ ð40T6 þ 0:002T4 þ 200T8 þ 200T2 Þ=440:002T6 ¼ ð40T5 þ 300T9 þ 300T3 Þ=640T7 ¼ ð0:002T8 þ 0:001T4 þ 1:2394Þ=0:006998T8 ¼ ð40T9 þ 0:002T7 þ 200T5 þ 1188:12Þ=243:9624T9 ¼ ð40T8 þ 300T6 þ 1782:18Þ=345:9406The solution is found iteratively with the latest available values usedduring the solution process.
Starting with a uniform initial guess of 300 K,Table 8.1 shows the results during the first two iterations along with the finalconverged solution.2388Spatial Discretization: The Diffusion TermTable 8.1 Summary of results obtained using the Gauss-Seidel iterative methodIterTlT2T3T4T5T6T7T8T9012⋮Solution300308.000309.633…312.490300300.083300.131…305.254300300.098300.146…305.254300308.000309.428…311.944300300.037300.078…305.154300300.048300.094…305.154300306.859307.072…308.867300300.031300.071…305.054300300.045300.090…305.054Values of T along the bottom boundaryThese can be calculated from the specified boundary condition, in this casea specified flux, askb@TTC Tbqb¼ qb ) Tb ¼ TC þ ðyC yb Þ¼ qb ) kbyC ybkb@yUsing the above equation, the temperatures are calculated asT10 ¼ T1 ¼ 312:490qb100ðy2 y11 Þ ¼ 305:254 þð0:05 0Þ ¼ 305:304T11 ¼ T2 þk11100qb100ðy3 y12 Þ ¼ 305:254 þð0:05 0Þ ¼ 305:304T12 ¼ T3 þk12100Values of T along the right boundaryThe above equation with the heat flux set to zero can be used to calculatethe temperature along the zero flux boundary, leading tokb@T¼ 0 ) Tb ¼ TC@yThe boundary temperatures are therefore given byT13 ¼ T3 ¼ 305:254 T14 ¼ T6 ¼ 305:154 T15 ¼ T9 ¼ 305:054Values of T along the top boundaryUsing Eq.
(8.45) the temperature values along the top boundary arecomputed ashb T1 þ ðk18 =dy718 ÞT7 20 300 þ ð103 =0:05Þ308:867¼¼ 300:009hb þ ðk18 =dy718 Þ20 þ ð103 =0:05Þhb T1 þ ðk17 =dy817 ÞT8 20 300 þ ð102 =0:05Þ305:054¼¼ 305:004¼hb þ ðk17 =dy817 Þ20 þ ð102 =0:05Þhb T1 þ ðk16 =dy916 ÞT9 20 300 þ ð102 =0:05Þ305:054¼¼ 305:004¼hb þ ðk16 =dy916 Þ20 þ ð102 =0:05ÞT18 ¼T17T168.4 The Interface Diffusivity239Total heat transfer along the left boundaryThe total heat transfer along the west boundary is given byQLeft ¼ q21 Dy21 þ q20 Dy20 þ q19 Dy19T1 T21T4 T20T7 T19¼ k21Dy21 þ k20Dy20 þ k19Dy19dx211dx204dx197103 0:1½312:49 þ 311:944 þ 308:867 3 320¼0:05¼ 0:053398 WTotal heat transfer along the top boundaryAlong the top boundary, the total heat transfer is computed asQtop ¼ q18 Dx18 þ q17 Dx17 þ q16 Dx16¼ h1 ½Dx18 ðT18 T1 Þ þ Dx17 ðT17 T1 Þ þ Dx16 ðT16 T1 Þ¼ 20½0:1ð300:009 300Þ þ 0:2ð305:004 300Þ þ 0:3ð305:004 300Þ¼ 50:058 WTotal heat transfer along the bottom boundaryKnowing the heat flux, the total amount of heat transfer is found asQBottom ¼ qb Dy ¼ 100 0:5 ¼ 50 WCheck for energy conservationFor conservation, the sum of heat entering and leaving the domain shouldbe equal to zero.
The sum is computed asQTop þ QLeft þ QBottom ¼ 50:058 0:053398 50 ¼ 0:004602 0The slight deviation from zero is due to the use of a limited number ofdigits during computations.8.5Non-Cartesian Orthogonal GridsLet us now consider the case of an orthogonal grid that is not oriented along thex and y axes. As shown in Fig. 8.8, such a grid can be obtained by rotating theCartesian grid of Fig. 8.1 by some angle.The discretized equation for this grid should be exactly the same as the oneobtained for the Cartesian grid.
Moreover for similar boundary conditions, the samesolution should be obtained.2408Spatial Discretization: The Diffusion TermNEyxNSnNW( x )CESeCtSwSsWESeSwSEEnCC( x )eWSy(yxWy )C( x )wd CE = ( x )e nd CE = ( x )eSWSe = Se nSw = SexSe = ( y)CSw = ( y)CFig. 8.8 Example of non-Cartesian orthogonal gridsConsider again the steady state conduction Eq.
(8.1)r J/;D ¼ Q/ð8:58ÞAs before its discretized form is given byXf nbðCÞ C/ r/ f Sf ¼ Q/C VCð8:59ÞConsidering the discretization along face e, the following is obtained:Je/;D Se ¼C/e ðr/ nÞe Se ¼C/e @/Se@n eð8:60Þwhere nowðr/ nÞe ¼@/@nð8:61Þeis the gradient of ϕ at face e along the n direction. Assuming again a linear profilefor ϕ along the n coordinate axis, the gradient can be written as @// /C¼ EdCE@n eð8:62ÞThe discretization of other terms proceeds as for a Cartesian grid, leading to thesame final discretization equation.8.6 Non-orthogonal Unstructured Grid8.68.6.1241Non-orthogonal Unstructured GridNon-orthogonalityIn the above configurations, the fluxes were normal to the face.
In general, structured curvilinear grids and unstructured grids are non-orthogonal [2–5].Thereforethe surface vector Sf and the vector CF joining the centroids of the elementsstraddling the interface are not collinear (see Fig. 8.9). In this case the gradientnormal to the surface cannot be written as a function of ϕF and ϕC, as it has acomponent in the direction perpendicular to CF.Thus while on orthogonal grids the gradient in the direction normal to theinterface yields @// /C/ /C¼ F¼ Fðr/ nÞf ¼dCF@n f krF rC kð8:63Þbecause CF and n (the unit vector normal to the surface) are aligned, onnon-orthogonal grids [6, 7], the gradient direction that yields an expressioninvolving ϕF and ϕC will have to be along the line joining the two points C and F.If e represents the unit vector along the direction defined by the line connectingnodes C and F thene¼rF rCdCF¼krF rC k dCFð8:64ÞTherefore, the gradient in the e direction can be written as @// /C/ /C¼ Fðr/ eÞf ¼¼ FdCF@e f krF rC kð8:65ÞThus to achieve the linearization of the flux in non-orthogonal grids, the surfacevector Sf should be written as the sum of two vectors Ef and Tf, i.e.,Sf ¼ Ef þ Tfð8:66ÞFig.