Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 38
Текст из файла (страница 38)
8.5 Mixed typeboundary conditionbNNWSnSwWbCSbbhSsSWS8.3 Boundary Conditions223from which an equation for /b is obtained ash1 /1 þ C/b =dxb /C/b ¼h1 þ C/b =dxbð8:45ÞSubstituting /b back in Eq. (8.43), the flux equation is transformed to3h1 C/b =dxb/;D Sb 5 ð / 1 / C ÞJb S b ¼ 4h1 þ C/b =dxb|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}2ð8:46ÞReq¼ FluxCb /C þ FluxVbwhere nowFluxCb ¼ ReqFluxVb ¼ Req /1ð8:47ÞUsing the flux term given by Eq. (8.47) in the discretized equation of theboundary element C, the modified equation is found to beaC /C þ aW /W þ aN /N þ aS /S ¼ bCð8:48ÞwhereaE ¼ 0aW ¼ FluxFw ¼ C/w gDiffwaN ¼ FluxFn ¼ C/n gDiffnaS ¼ FluxFs ¼ C/s gDiffsXð8:49ÞaC ¼ FluxCb þFluxCf ¼ ðFluxCb þ FluxCw þ FluxCn þ FluxCs Þf nbðCÞ0bC ¼ Q/C VC @FluxVb þX1FluxVf Af ¼nbðCÞ8.3.4Symmetry Boundary ConditionAlong a symmetry boundary the normal flux to the boundary of a scalar variable /is zero.
Therefore a symmetry boundary condition is equivalent to a Neumannboundary condition with the value of the flux set to zero (i.e., FluxCb = FluxVb = 0).2248Spatial Discretization: The Diffusion TermThus the modified equation along a symmetry boundary condition can be deducedfrom Eqs. (8.40) and (8.41) by setting qb to zero and is given byaC /C þ aW /W þ aN /N þ aS /S ¼ bCð8:50ÞwhereaE ¼ 0aW ¼ FluxFw ¼ C/w gDiffwaN ¼ FluxFn ¼ C/n gDiffnaS ¼ FluxFs ¼ C/s gDiffsXaC ¼FluxCf ¼ ðaW þ aN þ aS Þf nbðCÞbC ¼ Q/C VC Xð8:51ÞFluxVff nbðCÞ8.4The Interface DiffusivityIn the above discretized equation, the diffusion coefficients C/e ; C/w ; C/n , and C/s havebeen used to represent the value of C/ at the e, w, n and s faces of the element,respectively.
When the diffusion coefficient C/ varies with position, its value will beknown at the cell centroids E, W, … and so on. Then a prescription for evaluating theinterface value is needed in terms of values at these grid points. The followingdiscussion is, of course, not relevant to situations of uniform diffusion coefficient.The discussion may become clearer if the energy equation is considered, inwhich case the diffusion coefficient represents the conductivity of the material usedand / denotes the temperature. Non-uniform conductivity occurs innon-homogeneous materials and/or when the thermal conductivity is temperaturedependent.
In the treatment of the general differential equation for /, the diffusioncoefficient C/ will be handled in the same way. Significant variations of C/ arefrequently encountered for example, in turbulent flow, where C/ may stand for theturbulent viscosity or the turbulent conductivity. Thus, a proper formulation fornon-uniform C/ is highly desirable.A simple approach for calculating the interface conductivity is the linear profileassumption for the variation of C/ , between point C and any of its neighbors. Thus,along the east face the value of C/ , is obtained asC/e ¼ ð1 ge ÞC/C þ ge C/Eð8:52Þ8.4 The Interface Diffusivity225where the interpolation factor ge is a ratio in terms of distances between centroidsgiven byge ¼dCedCe þ deEð8:53ÞHence for a Cartesian grid if the interface is midway between the grid points, gewould be 0.5, and C/e would be the arithmetic mean of C/C and C/E .
Similar coefficients can be defined for other faces asdCwdCw þ dwWdCngn ¼dCn þ dnNdCsgs ¼dCs þ dsSgw ¼ð8:54ÞTherefore it is sufficient to calculate the coefficients of each surface of an element only once. Moreover, the use of distances instead of volumes will lead to thesame interpolation factors as the grid here is Cartesian.This basic approach leads to rather incorrect implications in some cases and cannotaccurately handles, for example, abrupt changes of conductivity that may occur incomposite materials. Fortunately, a much better alternative of comparable simplicityis available. In developing this alternative, it is recognized that the local value ofconductivity at an interface is not of primary concern. Rather, the main objective is toobtain a good representation of the diffusion flux Jϕ,D at the interface [1].For the one dimensional problem shown in Fig. 8.6, it is assumed that theelement C is composed of a material having a thermal conductivity C/C , whileelement E is made of a material of thermal conductivity C/E .
For the non-homogeneous slab between points C and E, a steady one-dimensional analysis (withoutsources) results in (the flux on either side of the interface e is supposed to be thesame)Je/;D Se ¼/C /eðdxÞCeC/C¼/e /EðdxÞeEC/E/ // /¼ ðdxÞ C ðdxEÞ ¼ CðdxÞ ECeeECEþ ///CCCe( x )e( x )wwWCEeC( x)CFig. 8.6 Interpolation of properties at element facesEð8:55Þ2268Spatial Discretization: The Diffusion TermHence the effective conductivity for the slab is found to beðdxÞCE ðdxÞCe ðdxÞeE1¼þ) /¼//C/eCCCCEe1 geC/EþgeC/C!ð8:56ÞWhen the interface is halfway between C and E (ge = 0.5), Eq.
(8.56) reduces toC/e ¼2C/C C/EC/C þ C/Eð8:57Þwhich is the harmonic mean of C/C and C/E , rather than the arithmetic mean.It is important to note that the harmonic mean interpolation for discontinuousdiffusion coefficients is exact only for one-dimensional diffusion. Nevertheless, itsapplication for multi-dimensional situation has an important advantage. With thistype of interpolation, nothing special need to be done when treating conjugateinterfaces. Solid and fluid cells are simply treated as part of the same domain withdifferent diffusion coefficients stored at the cell centroids.
By calculating the facediffusivity as the harmonic-mean of the values at the centroids sharing the face, thediffusion flux at the conjugate interface is correctly computed.Example 1The heat conduction in the two-dimensional rectangular domain composed oftwo materials shown in Fig. 8.7 is governed by the following differentialequation:r ðkrT Þ ¼ 0where T represents temperature. For the thermal conductivities (k) andboundary conditions displayed in the figure:(a) Derive the algebraic equations for all the elements shown in the figure.(b) Using the Gauss-Seidel iterative method discussed in Chap. 4, solve thesystem of equations obtained and compute the cell values of T.(c) Compute the values of T at the bottom, right, and top boundaries.(d) Compute the net heat transfer through the top, bottom, and leftboundaries and check that energy conservation is satisfied.8.4 The Interface Diffusivity227Fig.
8.7 Conduction heat transfer in a two dimensional rectangular domainSolutionThe first step is to determine the geometric quantities that are needed in thesolution. The coordinates of the various nodes are found to bex21x10x11x12x13¼ x20 ¼ x19 ¼ 0¼ x1 ¼ x4 ¼ x7 ¼ x18 ¼ 0:05¼ x2 ¼ x5 ¼ x8 ¼ x17 ¼ 0:2¼ x3 ¼ x6 ¼ x9 ¼ x16 ¼ 0:45¼ x14 ¼ x15 ¼ 0:6y10y21y20y19y18¼ y11 ¼ y12 ¼ 0¼ y1 ¼ y2 ¼ y3 ¼ y13 ¼ 0:05¼ y4 ¼ y5 ¼ y6 ¼ y14 ¼ 0:15¼ y7 ¼ y8 ¼ y9 ¼ y15 ¼ 0:25¼ y17 ¼ y16 ¼ 0:3The distance between nodes are computed asdx211 ¼ dx204 ¼ dx197 ¼ 0:05dx12 ¼ dx45 ¼ dx78 ¼ 0:15dx23 ¼ dx56 ¼ dx89 ¼ 0:25dx313 ¼ dx614 ¼ dx915 ¼ 0:15dy101 ¼ dy112 ¼ dy123 ¼ 0:05dy14 ¼ dy25 ¼ dy36 ¼ 0:1dy47 ¼ dy58 ¼ dy69 ¼ 0:1dy718 ¼ dy817 ¼ dy916 ¼ 0:052288Spatial Discretization: The Diffusion TermThe size of the elements are given byDx1 ¼ Dx4 ¼ Dx7 ¼ 0:1Dx2 ¼ Dx5 ¼ Dx8 ¼ 0:2Dx3 ¼ Dx6 ¼ Dx9 ¼ 0:3Dy1 ¼ Dy2 ¼ Dy3 ¼ 0:1Dy4 ¼ Dy5 ¼ Dy6 ¼ 0:1Dy7 ¼ Dy8 ¼ Dy9 ¼ 0:1The computed volumes of the various cells are obtained asV1 ¼ V4 ¼ V7 ¼ 0:01V2 ¼ V5 ¼ V8 ¼ 0:02V3 ¼ V6 ¼ V9 ¼ 0:03The interpolation factors needed to find values at element faces are calculated asV10:01¼¼ 0:333V1 þ V2 0:01 þ 0:02V20:02ðge Þ2 ¼ ðge Þ5 ¼ ðge Þ8 ¼¼¼ 0:4V2 þ V3 0:02 þ 0:03V10:01ðgn Þ1 ¼ ðgn Þ2 ¼ ðgn Þ3 ¼¼¼ 0:5V1 þ V4 0:01 þ 0:01V100:01¼ 0:5ðgn Þ4 ¼ ðgn Þ5 ¼ ðgn Þ6 ¼¼V10 þ V15 0:01 þ 0:01ðge Þ1 ¼ ðge Þ4 ¼ ðge Þ7 ¼The thermal conductivity values over the domain are given byk1 ¼ k4 ¼ k7 ¼ k10 ¼ k21 ¼ k20 ¼ k19 ¼ k18 ¼ 103k2 ¼ k3 ¼ k5 ¼ k6 ¼ k8 ¼ k9 ¼ k11 ¼ k12 ¼ k13 ¼ k14 ¼ k15 ¼ k16 ¼ k17 ¼ 102while values at the element faces are found to bek1 k2103 102k12 ¼ ¼ 3 1031 ðge Þ1 k1 þ ðge Þ1 k2 ð1 0:333Þ 103 þ 0:333 102k12 ¼ k45 ¼ k78 ¼ 3 103k14 ¼ k47 ¼ 103k25 ¼ k36 ¼ k58 ¼ k69 ¼ 102Using the above values, the discretized algebraic equations for all elementsare derived next.Element #1The needed diffusion terms are calculated asDy10:1¼¼ 0:667dx12 0:15Dx10:1¼¼1gDiffn ¼dy14 0:1gDiffe ¼gDiffw ¼Dy10:1¼¼2dx211 0:058.4 The Interface Diffusivity229The interface conductivities areke ¼ k12 ¼ 3 103kn ¼ k14 ¼ 103kw ¼ k21 ¼ 103The general form of the equation is written asaC T1 þ aE T2 þ aN T4 ¼ bCwhereaE ¼ FluxFe ¼ ke gDiffe ¼ 3 103 0:667 ¼ 0:002aN ¼ FluxFn ¼ kn gDiffn ¼ 103 1 ¼ 0:001The west and south coefficients do not appear in the equation as theirinfluence is integrated through the boundary conditions asFluxCw ¼ kw gDiffw ¼ 103 2 ¼ 0:002FluxVs ¼ 0FluxVw ¼ kw gDiffw T21 ¼ 103 2 320 ¼ 0:64The main coefficient and source term can now be calculated and are given byaC ¼ FLuxCe þ FLuxCn þ FLuxCw ¼ 0:002 þ 0:001 þ 0:002 ¼ 0:005bC ¼ FluxVS FluxVw ¼ 0 þ 0:64 ¼ 0:64Substituting, the discretized algebraic equation is obtained as0:005T1 0:002T2 0:001T4 ¼ 0:64Element #2The needed diffusion terms are calculated asDy20:1Dy20:1¼¼¼ 0:4 gDiffw ¼¼ 0:667dx23 0:25dx12 0:15Dx20:2¼¼2gDiffn ¼dy25 0:1gDiffe ¼The interface conductivities areke ¼ k23 ¼ 102kn ¼ k25 ¼ 102kw ¼ k12 ¼ 3 1032308Spatial Discretization: The Diffusion TermThe general form of the equation is written asaC T2 þ aE T3 þ aW T1 þ aN T5 ¼ bCwhereaE ¼ FluxFe ¼ ke gDiffe ¼ 102 0:4 ¼ 40aW ¼ FluxFw ¼ kw gDiffw ¼ 3 103 0:667 ¼ 0:002aN ¼ FluxFn ¼ kn gDiffn ¼ 102 2 ¼ 200The south coefficient does not appear in the equation as its influence isintegrated through the boundary conditions asFluxVs ¼ qb j S1 ¼ qb j ðDx2 jÞ ¼ 100 0:2 ¼ 20The main coefficient and source term can now be calculated and are givenbyaC ¼ FLuxCe þ FLuxCn þ FLuxCw ¼ 40 þ 0:002 þ 200 ¼ 240:002bC ¼ FluxVS ¼ 20Substituting, the discretized algebraic equation is obtained as240:002T2 40T3 0:002T1 200T5 ¼ 20Element #3The needed diffusion terms are calculated asgDiffw ¼Dy30:1¼ 0:4¼dx23 0:25gDiffn ¼0:3Dx3¼3¼0:1dy36The interface conductivities arekw ¼ k23 ¼ 102kn ¼ k36 ¼ 102The general form of the equation is written asaC T3 þ aW T2 þ aN T6 ¼ bCwhereaW ¼ FluxFw ¼ kw gDiffw ¼ 102 0:4 ¼ 40aN ¼ FluxFn ¼ kn gDiffn ¼ 102 3 ¼ 3008.4 The Interface Diffusivity231The east and south coefficients do not appear in the equation as theirinfluence is integrated through the boundary conditions asFluxVe ¼ 0FluxVs ¼ qb j S1 ¼ qb j ðDx3 jÞ ¼ 100 0:3 ¼ 30The main coefficient and source term can now be calculated and are givenbyaC ¼ FLuxCw þ FLuxCn ¼ 40 þ 300 ¼ 340bC ¼ FluxVS FluxVe ¼ 30 þ 0 ¼ 30Substituting, the discretized algebraic equation is obtained as340T3 40T2 300T6 ¼ 30Element #4The needed diffusion terms are calculated asDy40:1Dy40:1¼ 0:667 gDiffw ¼¼2¼¼dx45 0:15dx204 0:05Dx40:1Dx40:1gDiffn ¼¼¼¼1gDiffs ¼¼1dy47 0:1dy14 0:1gDiffe ¼The interface conductivities areke ¼ k45 ¼ 3 103kw ¼ k20 ¼ 103kn ¼ k47 ¼ 103ks ¼ k14 ¼ 103The general form of the equation is written asaC T4 þ aE T5 þ aN T7 þ aS T1 ¼ bCwhereaE ¼ FluxFe ¼ ke gDiffe ¼ 3 103 0:667 ¼ 0:002aN ¼ FluxFn ¼ kn gDiffn ¼ 103 1 ¼ 0:001aS ¼ FluxFs ¼ ks gDiffs ¼ 103 1 ¼ 0:0012328Spatial Discretization: The Diffusion TermThe west coefficient does not appear in the equation as its influence isintegrated through the boundary condition asFluxCw ¼ kw gDiffw ¼ 103 2 ¼ 0:002FluxVw ¼ kw gDiffw T20 ¼ 103 2 320 ¼ 0:64The main coefficient and source term can now be calculated and are givenbyaC ¼ FLuxCe þ FLuxCw þ FLuxCn þ FLuxCs¼ 0:002 þ 0:002 þ 0:001 þ 0:001 ¼ 0:006bC ¼ FluxVw ¼ 0:64Substituting, the discretized algebraic equation is obtained as0:006T4 0:002T5 0:001T7 0:001T1 ¼ 0:64Element #5The needed diffusion terms are calculated asDy50:1¼¼ 0:4dx56 0:25Dx50:2gDiffn ¼¼¼2dy58 0:1gDiffe ¼Dy50:1¼¼ 0:667dx45 0:15Dx50:2¼gDiffs ¼¼2dy25 0:1gDiffw ¼The interface conductivities areke ¼ k56 ¼ 102kw ¼ k45 ¼ 3 103kn ¼ k58 ¼ 102ks ¼ k25 ¼ 102The general form of the equation is written asaC T5 þ aE T6 þ aW T4 þ aN T8 þ aS T2 ¼ bCwhereaE ¼ FluxFe ¼ ke gDiffe ¼ 102 0:4 ¼ 40aW ¼ FluxFw ¼ kw gDiffw ¼ 3 103 0:667 ¼ 0:002aN ¼ FluxFn ¼ kn gDiffn ¼ 102 2 ¼ 200aS ¼ FluxFs ¼ ks gDiffs ¼ 102 2 ¼ 2008.4 The Interface Diffusivity233All coefficients appear in the equation as this is an internal element.