Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 37
Текст из файла (страница 37)
(8.5), the algebraic form of the diffusion equation isobtained asaC /C þ aE /E þ aW /W þ aN /N þ aS /S ¼ bCð8:18ÞwhereaE ¼ FluxFe ¼ C/e gDiffeaW ¼ FluxFw ¼ C/w gDiffwaN ¼ FluxFn ¼ C/n gDiffnaS ¼ FluxFs ¼ C/s gDiffsaC ¼ FluxCe þ FluxCw þ FluxCn þ FluxCsð8:19Þ¼ ðaE þ aW þ aN þ aS ÞbC ¼ Q/C VC ðFluxVe þ FluxVw þ FluxVn þ FluxVs Þor, more compactly, asaC /C þXaF /F ¼ bCð8:20ÞFNBðCÞwithaF ¼ FluxFf ¼ C/f gDifffXaC ¼FluxCff nbðCÞbC ¼Q/C VCXð8:21ÞFluxVff nbðCÞwhere the subscript F denotes the neighbors of element C (E, W, N, S), and thesubscript f denotes the neighboring faces of element C (e, w, n, s).2168.28Spatial Discretization: The Diffusion TermComments on the Discretized EquationA proper discretization method should result in a discretized algebraic equation thatreflects the characteristics of the original conservation equation.
The properties ofthe discretization techniques were introduced in the previous chapter and twoadditional rules that the coefficients of the discretized equation have to satisfy arepresented next.8.2.1The Zero Sum RuleLooking back at the discretization process, the first major approximation made wasthe linear profile assumption for the variation of ϕ between the centroids of theelements straddling the element face. The reader might ask why to use a first orderrather than a higher order profile. To answer this question, a one dimensionalconfiguration with no source term is considered.
Under these conditions the discretized equation reduces toaC /C þ aE /E þ aW /W ¼ 0ð8:22ÞwhereaE ¼ C/e gDiffeaW ¼ C/w gDiffwaC ¼ ð aE þ aW Þð8:23ÞIn the absence of any source or sink within this one dimensional domain, thetransfer of ϕ occurs by diffusion only and is governed by Fourier’s law (ellipticequation), i.e., in the direction of decreasing ϕ. As such, the value of ϕe or ϕw shouldlie between the values ϕC and ϕE or ϕC and ϕW, respectively, which is guaranteed bythe linear profile. A second order profile (e.g., a parabolic profile) may result in avalue at the face that is higher or lower than the values at the centroids of the cellsstraddling the face, which is unphysical. The same is true for other higher orderprofiles.
If the discretization scheme is to guarantee physical results, then a linearprofile should be used. Moreover, the adopted profile becomes less important as thesize of the element decreases, since all approximations are expected to yield thesame analytical solution in the limit when the size of the element approaches zero.Furthermore, in the absence of any source term, the multi-dimensional heat conduction equation reduces to,r C/ r/ ¼ 0ð8:24ÞThis implies that ϕ and ϕ + constant are solutions to the conservation equation.A consistent discretization method should reflect this property through its discretized equation and the discretized equation should fulfill8.2 Comments on the Discretized EquationaC / C þX217aF /F ¼ 0FNBðCÞaC ð/C þ constantÞ þXFNBðCÞ9>>=aF ð/F þ constantÞ ¼ 0 >>;) aC þXaF ¼ 0FNBðCÞð8:25Þwhich is actually satisfied by the discretized equation.
The above equation, which isvalid in the presence or absence of a source/sink term as revealed by Eq. (8.19),may be written asXaC ¼ aFð8:26ÞFNBðCÞor asX aF¼ 1aFNBðCÞ Cð8:27ÞThus ϕC can be viewed as the weighted sum of its neighbors, and in the absenceof any source term it should always be bounded by these neighboring ϕF values.When a source term is present, i.e., when S/C 6¼ 0; /C does not need to be boundedin this manner, and can over/undershoot the neighboring values, but this is perfectlyphysical.
The extent of over/under shoot is determined by the magnitude of S/C withrespect to the size of the coefficients at the neighboring nodes (aF).8.2.2The Opposite Signs RuleThe above derivations demonstrated that the coefficients aC and aF are of oppositesigns. This is of physical significance implying that as the value of ϕF isincreased/decreased, the value of ϕC is expected to increase/decrease. This isbasically related to boundedness suggesting that a sufficient condition for thisproperty to be satisfied is for the neighboring and main coefficients to be of oppositesigns.
If not, then the boundedness property may not be enforced.8.3Boundary ConditionsIt is well known that the analytical solution to any ordinary or partial differentialequation is obtained up to some constants that are fixed by the applicable boundaryconditions to the situation being studied.
Therefore, using different boundary2188Spatial Discretization: The Diffusion Termconditions will result in different solutions even though the general equationremains the same. Numerical solutions follow the same constrain, necessitatingcorrect and accurate implementation of boundary conditions as any slight change inthese conditions introduced by the numerical approximation leads to a wrongsolution of the problem under consideration.Boundary conditions will be discussed as deemed relevant to the terms beingdiscretized. For conduction/diffusion problems Dirichlet, Neumann, mixed, andsymmetry boundary condition types are encountered, which are detailed next.Boundary conditions are applied on boundary elements, which have one or morefaces on the boundary.
Discrete values of ϕ are stored both at centroids of boundarycells and at centroids of boundary faces.Let C denotes the centroid of the boundary element shown in Fig. 8.3 with oneboundary face of centroid b and of surface vector Sb pointing outward. As before,the discretization process over cell C yieldsX J/;D S f ¼ Q/C VCð8:28Þf nbðCÞThe fluxes on the interior faces are discretized as before, while the boundary fluxis discretized with the aim of constructing a linearization with respect to ϕC, thus Sb ¼ FluxTbJ/;Db¼ C/b ðr/Þb Sbð8:29Þ¼ FluxCb /C þ FluxVbThe specification of boundary conditions involves either specifying the unknownboundary value ϕb, or alternatively, the boundary flux Jb/;D . Using Eq.
(8.18), thediscretized equation at a boundary element for the different boundary conditiontypes of diffusion problems are derived next.8.3.1Dirichlet Boundary ConditionA Dirichlet boundary condition is a type of boundary condition that specifies thevalue of ϕ at the boundary, i.e.,/b ¼ /specifiedð8:30Þ8.3 Boundary Conditions219( x)Fig. 8.3 Value specifiedboundary conditionbNNWSnSwWbCSbb =specifiedSsSWSFor this caseFluxTb ¼ C/b ðr/Þb SbkSb k¼ C/bð/ /C ÞkdCb k b¼ FluxCb /C þ FluxVbð8:31ÞyieldingFluxCb ¼ C/b gDiffb ¼ abFluxVb ¼ C/b gDiffb /b ¼ ab /bð8:32ÞwithgDiffb ¼SbdCbð8:33ÞThus for element C shown in Fig. 8.3 the aE coefficient is zero reducing thediscretized equation toaC /C þ aW /W þ aN /N þ aS /S ¼ bCð8:34ÞwhereaE ¼ 0aW ¼ FluxFw ¼ C/w gDiffwaN ¼ FluxFn ¼ C/n gDiffnaS ¼ FluxFs ¼ C/s gDiffsXaC ¼ FluxCb þFluxCf ¼ FluxCb þ ðFluxCw þ FluxCn þ FluxCs Þf nbðCÞð8:35Þ2208andSpatial Discretization: The Diffusion Term0bC ¼Q/C VC @FluxVb þX1FluxVf Að8:36Þf nbðCÞThe following important observations can be made about the discretizedboundary equation:1.
The coefficient ab is larger than other neighbor coefficients because b is closer toC and consequently has a more important effect on ϕC.2. The coefficient aC is still the sum of all neighboringP coefficients including ab.This means that for the boundary elementjaF j=jaC j\1 giving theF NBðC Þsecond necessary condition to satisfy the Scarborough criterion, thus guaranteeing, at any one iteration, the convergence of the linear system of equationsvia an iterative solution method.3. The abϕb product (=FluxVb) is now on the right hand side of the equation, i.e.,part of bC, because it contains no unknowns.8.3.2Von Neumann Boundary ConditionIf the flux (or normal gradient to the face) of ϕ is specified at the boundary(Fig. 8.4), then the boundary condition is denoted by a Neumann boundary condition.
In this case the specified flux is given by C/ r/ b i ¼ qbð8:37Þwhich is, in effect, the flux J/;Db , since Sb ¼ C/ r/ b kSb ki ¼ qb kSb kJ/;Db¼ FluxCb /C þ FluxVbð8:38Þwhere nowFluxCb ¼ 0FluxVb ¼ qb Sb¼ qb ðDyÞCð8:39ÞHere the flux components are assumed to be positive when they act in the samedirection as the coordinate system used.8.3 Boundary Conditions221( x)Fig. 8.4 Flux specifiedboundary conditionbNNWSnbSwWbCnb = qspecifiedSbSsSWSThus qb may directly be included in Eq. (8.18) to yield the following discreteequation for the boundary cell C:aC /C þ aW /W þ aN /N þ aS /S ¼ bCð8:40ÞwhereaE ¼ 0aW ¼ FluxFw ¼ C/w gDiffwaN ¼ FluxFn ¼ C/n gDiffnaS ¼ FluxFs ¼ C/s gDiffsXaC ¼FluxCf ¼ ðaW þ aN þ aS Þf nbðCÞ0bC ¼ Q/C VC @FluxVb þXð8:41Þ1FluxVf Af nbðCÞThe following important points can be made about the above discretizedequation:1. A Von Neumann boundary condition does not result in a dominant aCcoefficient.2.
If both qb and S/C are zero, then /C will be bounded by its neighbors. Otherwise,/C can exceed (or fall below) the neighbor values of /, which is admissible. If/ is temperature, for example, then qb represents the heat flux applied at theboundary. Therefore if heat is added at the boundary, then the temperature in theregion close to the boundary is expected to be higher than that in the interior.2228Spatial Discretization: The Diffusion Term3.
Once /C is computed, the boundary value /b may be computed using/b ¼C/b gDiffb /C qbð8:42ÞC/b gDiffb4. Finally the Von Neumann condition can be considered as a natural boundarycondition for the Finite Volume method since, for the case where the specifiedflux is zero, nothing needs to be done in terms of discretization for the face,while a specified value of zero (Dirichlet condition) will still require the discretization to be carried out.8.3.3Mixed Boundary ConditionThe mixed boundary condition, schematically depicted in Fig. 8.5, refers to thesituation where information at the boundary is given via a convection transfercoefficient (h1 ) and a surrounding value for /ð/1 Þ asJ/;D Sb ¼ C/ r/ b iSb ¼ h1 ð/1 /b ÞðDyÞCbð8:43Þwhich can be rewritten asC/b Sb/b /C¼ h1 ð/1 /b ÞSbdxbð8:44Þ( x)Fig.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
