计算流体力学(CFD)文档——4. The scalar-transport equation

4. THE SCALAR-TRANSPORT EQUATION SPRING 20114.1 The scalar-transport equation 4.2 Control-volume notation4.3 The steady-state 1-d advection-diffusion equation 4.4 Discretising diffusion4.5 Discretising the source term4.6 Assembling the algebraic equations 4.7 Extension to 2 and 3 dimensions 4.8 Discretising advection (part 1) 4.9 Discretisation properties4.10 Discretising advection (part 2)4.11 Implementation of advanced advection schemes 4.12 Implementation of boundary conditions 4.13 Solution of the algebraic equations Summary Examples4.1 The Scalar-Transport EquationThis is a generic equation for transported physical quantities (momentum, energy,...) For an arbitrary control volume V := + inside V y of boundar out inside V SOURCE FLUX CHANGE OF RATEThe total flux comprises advection (transport with the flow) and diffusion (molecular, or due to turbulent fluctuations). The resulting advection-diffusion equation for concentration φ may be written, for each computational control volume (or “cell”) as:sourcediffusion advection change of rate SV A C V t faces=−φ+φ∑)()(d d(Any non-advective fluxes not described by gradient diffusion are assumed to be incorporated in the source term – see Section 2.2.5.)4.2 Control-Volume NotationThe commonest configurations are cell-centred storage or cell-vertex storage.nThis course focuses on structured meshes using cell-centred storage. (Unstructured meshes will be discussed briefly at the end of the course, but the 2nd edition of Versteeg and Malalasekera’s book gives a much better description.)A typical 3-d control volume is shown right. Relativeto the cell centre (point P) the coordinate directions Array are commonly denoted west, east, south, north,bottom, top with:•lower case w, e, s, n, b, t used for cell faces;•upper case W, E, S, N, B, T for adjacent nodes.For a Cartesian mesh these would usually correspondto ±x, ±y, ±z directions respectively.Cell-face areas will be denoted A w, A e, A s, A n, A b, A t.Cell volumes will be denoted V. In 2 dimensions onecan think of a single layer of cells with unit depth.When referring to the entire set of control volumes (as opposed to looking at a representative one) it is common to switch between geographic and ijk notation, so thatφP≡φijk, φE≡φi+1 jk , etc.4.3 The Steady-State 1-D Advection-Diffusion EquationConsider first the steady-state, 1-d advection-diffusion equation. This is worthwhile because: • it greatly simplifies the analysis; • it permits a hand solution of the discretised equations; • subsequent generalisation to 2 and 3 dimensions is straightforward; • in practice, discretisation of fluxes is generally carried out coordinate-wise; • many important theoretical problems are 1-d.Integral (Control-Volume) FormConservation for one control volume givessource flux flux w e =− (1) where “flux ” means the rate at which something passesthrough a given area (here, a cell face).If φ is the amount per unit mass, then the total flux hasadvective and diffusive parts:advection:φ)(uA diffusion: xA d d φ−If S is the source per unit length then the advection-diffusion equation for φ isx S x A uA x A uA w e d d d d = φ−φ−φ−φ(2)Conservative Differential FormDividing by x and taking the limit as x → 0 gives a corresponding differential equation:S xA uA x =φ−φ)d d (d d(3)Non-Conservative Differential FormMass conservation implies that uA = constant and hence (3) can also be writtenS xA x x uA =φ−φd d (d d d d (4)Note : • This system is quasi-1-d in the sense that the cross-sectional area A may vary. Tosolve a truly 1-d problem, set A = 1. The differential equation is thenS u x =−φ(d d• In this instance, , u , and S are assumed to be known. In the general CFD problem,u is itself the subject of a separate transport equation.eflux4.4 Discretising DiffusionGradient diffusion is usually discretised by central differencing :−→φ−A x A e e )(d dor)(d d P E e eD x A φ−φ−→φ− where D ≡(6)D is a diffusive transfer coefficientA similar expression is used for the west face:)(d d W P w wD x A φ−φ−→φ−Note : • In the finite-volume method, fluxes are required at cell faces , not nodes. • This approximation for (d φ/d x )e is second-order accurate in x (see later). • If the diffusivity varies then its cell-face value must be obtained by interpolation. • Equal weighting is applied to the two nodes on either side of the cell face, consistentwith diffusion acting equally in all directions. Later on, we shall see that this contrasts with advection, which has a directional bias.4.5 Discretising the Source TermWhen the source is proportional to the amount of fluid, the total source strength for the cell is (source per unit volume ) × (volume ) V S ×=(In 1-d problems V is the cell length, x ).S may depend on the solution φ as in the example above. The source term is conveniently broken down into solution-independent and solution-dependent parts in the linear form,≤φ+=P P P P s s b source (7) The reason for this particular form will become apparent later.4.6 Assembling the Algebraic EquationsAs seen in the above classroom examples, when there is no flow (u = 0) the steady-state diffusion problem discretises as follows.p P P W P w P E e w e s b D D sourceflux flux φ+=φ−φ+φ−φ−=−)()( or, collecting multiples of φP , φE and φW together: P E e P p e w W w b D s D D D =φ−φ−++φ−)(This is commonly written P E E P P W w b a a a =φ−φ+φ−or, in a notation that generalises to 2 and 3 dimensions (and to the case with flow): ∑=φ−φnodesadjacentP F F P P b a a(8)(8) is a canonical form for the discretised scalar-transport equation.There will be a discretised equation of this form for each variable and for each control volume. For one variable φ, if the nodal values are assembled into a vector then the set of algebraic equations takes the form= φ−−M M M M M M M M O O O O O O O O OO P P EP W b a a a 00 orFor a 1-d system this is tri-diagonal . If the coefficients are constants then it can be solved directly by Gaussian elimination or very efficiently on a computer by the tri-diagonal matrix algorithm (see the Appendix). If the elements of the matrix are dependent on the solution then it must be solved iteratively.4.7 Extension to 2 and 3 Dimensions For a multi-dimensional flow the net flux out of a cell can be obtained by summing the outward fluxes through all faces. For quadrilateral elements (in 2d) or hexahedral elements (in 3d) then the net flux out of a cell is simply the sum of the net fluxes through opposing sides and the general conservation equation may be written: source flux flux flux flux flux flux b t s n w e =−+−+−)()()( (9)The discretised equations are still assembled in the same matrix formP FF F P P b a a =φ−φ∑(10) with the summation extended to include nodes in the other directions.Combining the individual equations (10) from all control volumes gives a matrix equation for the set of nodal values. In two dimensions this gives the banded matrix system shown. Further bands appear in three dimensions.Thus, (8) or (10) has the same form in 1-, 2- or 3-d problems, but in multiple dimensions the summation is extended to the other (S , N , B , T ) nodes and there is a corresponding increase in the number of non-zero diagonals in the assembled matrix equation. This makes the matrix equation much harder to solve (see Section 4.13).Equation (10) formally describes the discretisation for a single control volume in an unstructured mesh. However, since the nodes do not have a simple ijk indexing, the resulting matrix and solution method for unstructured meshes are much more complex.ΦbΦbe flux s4.8 Discretising Advection (Part 1)Purely diffusive problems (u = 0) are of passing interest only. In typical engineering flow problems, advection fluxes far exceed diffusive fluxes because the Reynolds number (Re UL / = ratio of inertial forces [mass × acceleration] to viscous forces) is very large.The 1-d steady-state advection-diffusion equation issource flux flux w e =− where, with mass flux C (= uA ):xA C flux d d φ−φ= .Discretising the diffusion and source terms as before, the equation becomessourcediffusion advection s b D D C C PP P W P w P E e w w e e φ+=φ−φ−φ−φ−φ−φ)]()([][ (11)φe and φw have yet to be approximated. The problem is … how to approximate these face values in terms of the values at adjacent nodes . A method of specifying these face values in order to calculate advective fluxes is called an advection scheme or advection-differencing scheme .4.8.1 Central DifferencingIn central differencing the cell-face value is approximated by the average of values at the nodes on either side:)(21E P e φ+φ=φThis is second-order accurate for φe in terms of ∆x (see later). Substituting into (11) (with a similar expression for φw ) givesP P P W P w P E e P W w E P e s b D D C C φ+=φ−φ+φ−φ−φ+(φ−φ+(φ)()())2121 or, collecting terms,P E e e P P w e w e W w w b D C s D D C C D C =φ+−−)φ−++−+)φ+−)(((21212121In out canonical notation: P FF F P P b a a =φ−φ∑where, in one dimension:)(2121w e P W E P ee E w w W C C s a a a D C a D C a −+−+=+−=+= (12)(By mass conservation, C e – C w = 0, so that the expression for a P can be simplified.)The graphs below show the solution of an advection-diffusion problem with no sources and constant diffusivity for the two combinationsPe = 1/2 (advection « diffusion) (equation: 024345=φ−φ+φ−E P W )Pe = 4 (advection » diffusion)(equation: 023=φ+φ+φ−E P W )where the Peclet number Pe is defined byi.e.Pe diffusion advection D C ==(13)Pe = 1/2Pe = 4In the first case the solution is good (consistent with second-order accuracy).In the second case there are pronounced “wiggles” in what should be a perfectly smoothsolution. What is wrong?eMathematically, when the cell Peclet number Pe is bigger than 2, the a E coefficient becomes negative ,meaning that, for example, an increase in φE wouldlead to a decrease in φP . This is impossible for a quantity that is simply advected and diffused.Physically, the advection process is directional ; ittransports properties only in the direction of the flow.However, the central-differencing formula assigns equal weight to both upwind and downwind nodes.4.8.2 Upwind DifferencingIn upwind differencing φface is taken as the value at whichever is the upwind node; i.e. in one dimension:<φ>φ=φ)0if ()0if (u u EP eThis is only first-order accurate in x (see later) but acknowledges the directional nature of advection. The alternatives can be summarised as U face φ=φwhere subscript U denotes whichever is the upwind node for that face.With this scheme, the algebraic equation for each control volume takes the canonical form∑=φ−φFP F F P P b a a where:PW E P w w W ee E s a a a D C a D C a −+=+=+−=)0,max()0,max((14)If the “max” bit confuses you, consider separately the two cases where the mass flux is positive (flow from left to right) or negative (flow from right to left).When applied to the same advection-diffusion problem as the central-differencing scheme it is found that: • when Pe = ½ the upwind-differencing scheme is not as accurate as centraldifferencing; this is to be expected from its order of accuracy; • when Pe = 4 the upwind-differencing solution is not particularly accurate, but the“wiggles” have disappeared.In all cases, however, both a W and a E are unconditionally positive.So there is a pay-off – accuracy versus boundedness (absence of wiggles). The next sections examine the desirable properties of discretisation schemes, the constraints that they impose upon the matrix coefficients and more advanced advection schemes that are both accurate andbounded.diffusion only advection+diffusionφφe4.9 Discretisation Properties(i) ConsistencyAn approximation is consistent if the discretised equations are equivalent to the continuum equations in the limit as the grid size tends to zero. e.g. by the definition of a derivative,is a consistent approximation forx∂φ∂.(ii) ConservativenessA scheme is conservative if fluxes are associated with faces, not particularcells, so that what goes out of one cell goes into the adjacent cell;This is automatically built into the finite-volume method.(iii) TransportivenessAn advection scheme is transportive if it is upstream-biased. In practice this means higher weighting to nodes on the upstream side of a face.(iv) BoundednessA flux-differencing scheme is bounded if, in an advection-diffusion problem without sources: •the value of φ at a node always lies between the maximum and minimum values at surrounding nodes;•φ = constant is a possible solution.This imposes conditions on the matrix coefficients a P, a W, a E, ... If there are no sources then=φ−φ∑F FPPaa(15) Suppose that φ is only non-zero at one adjacent node, F. Then=φ−φFFPPaa orFPFP aaφ=φSince φP must lie between 0 and φF, this requires that10≤≤PFaaHence, a F and a P must have the same sign (invariably positive in practice). Thus, we require: FaFallfor≥ (“positive coefficients”) (16) (It is the contravening of the positivity condition that leads the central-differencing scheme to produce “wiggles”.)If (15) is also to admit the solution φ = constant then this can be divided out to yield ∑=FPaa (“sum of neighbouring coefficients”) (17)(v) StabilityA solution method (not advection scheme) is stable if small errors do not grow in the course of the calculation. This determines whether it is possible to obtain a converged solution – it says nothing about its accuracy or whether it contains wiggles. Stability is heavily influenced by how the source term is discretised.Any source term should be linearised as P P P s b φ+; the complete equation for one cell is then P P P F F P P s b a a φ+=φ−φ∑ If the solution-dependent part of the source (s P φP ) is transferred to the LHS then the diagonal coefficient is modified to read P F P s a a −=∑ Numerical stability requires negative feedback ; otherwise, an increase in φ would lead to an increase in the source, which would lead to a further increase in φ and so on. Thus: s P ≤ 0 (“negative-slope linearisation of the source term ”) If this condition and the positivity of the a F is maintained then ∑≥F P a a (“diagonal dominance ”)The last condition is, in fact, a necessary requirement of many matrix solution algorithms.In summary, boundedness and stability place the following constraints on the discretisation of flux and source terms:positive coefficients : F a F all for 0≥ negative-slope linearisation of the source term : 0,≤φ+=P P P P s s b source sum of neighbouring coefficients : ∑−=FP F P s a a(vi) OrderOrder is a measure of accuracy. Formally, it defines how fast the error in a numerical approximation diminishes as the grid spacing gets smaller.If, on a uniform grid of spacing x , the truncation error is proportional to x n as x 0 then that scheme is said to be of order n .Order can be established formally by a Taylor-series expansion about a cell face . e.g. for the nodes either side of the east face:L +∂φ∂1+ ∂φ∂+∂φ∂+φ=φ333222!3!21x x x e e e e E (18)(a)L +∂φ∂1− ∂φ∂+∂φ∂−φ=φ333222!3!21x x x e e e eP (18)(b)Subtracting (18)(a) – (b) gives:L +∂φ∂1++∂φ∂+=φ−φ333300x x x e e P E whence, dividing by x :)(2x O x e+∂φ∂=As the leading error term is O (∆x 2is a second-order approximation for ex∂φ∂.Alternatively, adding (18)(a) – (b) gives:L +∂φ∂++φ=φ+φ22202x e e E P whence:)()(221x O e E P +φ=φ+φAs the error term is O (x 2), the central differencing formula )(21E P φ+φ is a second-order approximation for φe . On the other hand, the Upwind-differencing approximations φP or φE (depending on the direction of the flow) are formally first-order accurate.Higher accuracy can be obtained by using more nodes in an approximation – one node permits schemes of at most first-order accuracy, two permit second-order accuracy and so on.Schemes of low-order accuracy – e.g. Upwind – lead to substantial numerical diffusion in 2- and 3-d calculations when the velocity vector is not aligned with the grid lines.Notes . (1) Order is an asymptotic concept; i.e. it refers to behaviour as x 0. In this limit onlythe first non-zero truncation term in the Taylor series is important. However, the full expansion includes terms of higher order which may be non-negligible for finite x . (2) Order refers to the theoretical truncation error in the approximation, not thecomputer’s round-off error (the accuracy with which it can store floating-point numbers). (3) The more accurate a scheme then, in principle, the greater the reduction in numericalerror as the grid is made finer or, conversely, the less nodes required for a given accuracy. However, high-order schemes tend to require more computational calculations and often have boundedness or stability problems. Also, the law of diminishing returns applies when the truncation error becomes of similar size to the floating-point round-off error.4.10 Discretising Advection (Part 2)With an understanding of the desirable properties of a flux-differencing scheme it is now possible to examine more advanced schemes.4.10.1 Exponential SchemeThe 1-d advection diffusion equation issource flux flux w e ≡− or, equivalently,S xA uA x =φ−φ)d d (d d (19) If there are no sources (S = 0) then the total flux must be constant:constant xA uA flux =φ−φ=d dIf , u , A and are constant this equation can be solved exactly, with boundary conditions φ = φP and φ = φE at adjacent nodes, to give (see the Examples):−φ−φ=1P P e E P e e e e Cflux wheret coefficien transfer diffusive fluxmass ====D uA C Pe DC==.With a similar expression for the west face, one obtains ∑φ−φ=−F F P P w e a a flux flux where:W E P E W a a a e Ca e Ce a +=−=−=,1,1Pe Pe Pe (20)Assessment • Conservative by construction. • Transportive, because there is a larger weighting on the upwind node. • Bounded: all a F are positive and a P is the sum of the neighbouring coefficients.To see the last two of these you will have to consider separately the cases C > 0 (for which e Pe > 1) and C < 0 (for which e Pe < 1).This scheme – by construction – gives the exact solution for zero sources and constant velocity and diffusivity , but this is something we could have found analytically anyway. The scheme has never really found favour because: • the scheme isn’t exact when u or vary, or if there are sources, or in 2-d or 3-d flow; • exponentials are extremely expensive to compute.eflux4.10.2 Hybrid Scheme (Spalding (1972)This is based on a piecewise-linear approximation to the exponential scheme. It amounts to: • central differencing if 2Pe ≤;• upwind differencing (with zero diffusion!) if 2Pe >. ∑φ−φ=−F F P P w e a a flux fluxwhere (if u > 0 and the mass flux C and diffusive transport coefficient D are constant):WE P E W E W a a a a C a D C DC aD C a +=>==≤≡+−=+=2Pe if 0,2/Pe if ,2121 (21)AssessmentThe scheme is conservative, transportive and bounded.The hybrid scheme remained extremely popular in commercial codes for a long time because it was stable and robust. However, most flows of interest operate in the high-advection/low-diffusion regime, where this scheme amounts to first-order upwinding (with no diffusion). Modern CFD practitioners seek much higher accuracy.Patankar also developed a power-law approximation to the exponential scheme, to overcome the heavy-handed switch-off of diffusion at Pe = 2. However, “powers” are as computationally expensive as exponentials.4.10.3 QUICK (QUadratic Interpolation for Convective Kinematics – Leonard, 1979)Fits a quadratic polynomial through 3 nodes to get 3rd -order accuracy.For each cell face, QUICK uses the nodes either side of thecell face, plus a further upwind node depending on thedirection of the flow as shown right.To emphasise the conservation property which associates fluxes with cell faces , notnodes, we shall, in future, for all such three-point schemes use the notation φD , φU and φUUfor the Downwind, Upwind and Upwind-Upwind nodes at any particular face.By fitting a quadratic polynomial to these nodes (see the Examples) the QUICK scheme gives:D UUU face φ+φ+φ−=φ834381 (22)u>0u<0For example, if u > 0 on the east face then:E P W e φ+φ+φ−=φ834381 whereas, if u < 0:P E E e φ+φ+φ−=φ834381Assessment • 3rd -order accurate. • Conservative by construction. • Transportive (upwind bias in the selection of the third node and relative weightings). • Not bounded; (for example, if u > 0 then a E is negative – see the Examples).Despite boundedness not being guaranteed (which can be a major problem in turbulent flows, where k and ε are required to be positive – see later) the high-order accuracy of the QUICK scheme make it popular and widely-used.4.10.4 Flux-Limited SchemesHitherto we have only seen schemes where the matrix coefficients a F are constants (i.e. independent of the solution φ). The only unconditionally-bounded scheme of this type is first-order upwind differencing. Schemes such as QUICK, which fit a polynomial through several points, are prone to generate cell-face values which lie outside the interpolating values φD , φU , φUU . To prevent this, modern schemes employ solution-dependent limiters , which enforce boundedness whilst trying to retain high-order accuracy wherever possible.For three-point schemes, φ is said to be: monotonic increasing if D U UU φ<φ<φ, monotonic decreasing if D U UU φ>φ>φ. A necessary condition for boundedness is that the schemes must default to first-order upwinding (i.e. φface = φU ) if φ is not locallymonotonic (either increasing or decreasing).Monotonic variation in φ may be gauged by whether the changes in φ between successive pairs of nodes have the same sign; i.e. 0))((>φ−φφ−φ⇔UU U U D monotonicSuch schemes can then be written (in the notation of Versteeg and Malalasakera) as the sum of the upstream value (φU ) and a solution-dependent fraction of the difference between downstream and upstream nodal values:φφ−φ+φ=φotherwise monotone if ))(UU D U face rwhere r is the ratio of successive differences (>0 where monotonic):U D UU U r φ−φφ−φ=monotonicThe limiter (r ) is given below for various schemes used at the University of Manchester.(r ) is taken as 0 if non-monotonic (i.e., r < 0).The choice of these examples is (obviously!) parochial. Many other equally-good schemes exist (see Versteeg and Malalasekera, 2nd edition, for a list). The important points about these schemes are that they are (a) bounded, and (b) non-linear (i.e. the resulting matrix elements depend on, and hence change with, the solution φ). The last property means that an iterative numerical solution is inevitable.4.11 Implementation of Advanced Advection SchemesThe general steady-state scalar transport equation issourcediffusion advection SV A C faces=−φ∑)( (23)If C and D are the outward mass flux and diffusive transport coefficient on each cell face, then, with the standard discretisation for diffusion and sources, (23) becomes P P P facesF P face s b D C φ+=φ−φ+φ∑)]([where F denotes an adjacent node and P is the index of the cell-centre node. Since 0=∑C by mass conservation, it is convenient to subtract P C φ∑ (which is 0) from both sides:P P P facesF P P faces b D C φ+=φ−φ+φ−φ∑)]()([(24)An advection scheme (Upwind, Central, QUICK, …) is needed to specify the cell-face value φface . Because many matrix-solution algorithms require positive coefficients and diagonal dominance, it is common practice to separate φface into the Upwind part plus a correction; i.e. )(U face U face φ−φ+φ=φThen, for the advective flux:)())(0,max()()()(U face F P correctionU face upwindP U P face C C C C C φ−φ+φ−φ−=φ−φ+φ−φ=φ−φ443442143421 (25)The first part of (25) gives rise to a positive matrix coefficient D C a F +−=)0,max(whilst the latter part is transferred to the RHS of the equation as what is called a deferred correction ; (“deferred” because it won’t be updated until the next iteration): 444344421correctiondeferred faces U face P P P FF P F C s b a ∑∑φ−φ−φ+=φ−φ)()(4.12 Implementation of Boundary ConditionsThe most common types of boundary condition are: • value φ specified (Dirichlet boundary conditions); e.g. φ = 0: velocity at a wall, or temperature fixed at some surface; • gradient ∂φ/∂n specified (Neumann boundary conditions). e.g. ∂φ/∂n = 0 on a symmetry plane, or at an outflow boundary.In the finite-volume method, both types of boundary condition can be implemented by transferring the boundary flux to the source term .For a cell abutting a boundary: source flux flux boundary boundarynot=+∑⇓boundary P F F boundarynot P P flux b a a −=φ−φ∑Thus, there are two modifications:• the a F coefficient in the direction of the boundary is set to zero; • the outward boundary flux is subtracted from the source terms.If flux (φ) is specified on the boundary, then this is immediate. If φ itself is fixed on the boundary node B then, with x the width of the cell:)(2)(P B D A flux φ−φ−=−=φTo subtract this flux from the source term requires a simple change of coefficients: D s s D b b P P B P P 2,2−→φ+→ (26)(You should revisit the classroom examples of Sections 4.4 and 4.5 to see this in action.)b o u n d a r yboundary4.13 Solution of the Algebraic EquationsThe discretisation of a single scalar-transport equation over a single control volume produces an algebraic equation of the form P F F P P b a a =φ−φ∑where the summation is over adjacent nodes. Combining the equations for all control volumes produces a set of simultaneous equations, i.e. a matrix equation b A =Φwhere Φ is the vector of nodal values. Matrix A is sparse (i.e. has only a few non-zero elements). Many algebraic methods have been used to tackle this problem; a few that are suitable for structured grids are mentioned below.4.13.1 Matrix Solution AlgorithmsGaussian EliminationThis is a direct (i.e. non-iterative) method. It consists of a sequence of row operations to obtain zeros below the main diagonal (upper-triangular matrix), followed by back-substitution.In general, it is inefficient because it tends to fill in sparse matrices (tridiagonal systems are an important exception), whilst for fluid-flow problems the matrix elements vary with the solution, so that an iterative solution is necessary anyway.Gaussian elimination is OK for small hand calculations, but not recommended for large systems of equations.Gauss-SeidelRearrange the equation foreach control volume as aniterative update for thecentral node:)*(1∑φ+=φF F P PP a b a(the asterisk * denotes the current value). Thenrepeatedly cycle through the entire set of equationsuntil convergence is achieved.Gauss-Seidel is very simple to code and is often used for unstructured grids. However, it tends to converge slowly for large matrices and may require substantial under-relaxation (see below).。