Department of Mechanical & Aerospace Engineering CARLETON UNIVERSITY
AERO 4304: Computational Fluid Dynamics Winter 2013
Lecture 14: Finite Volume Method II: Steady Convection-Diffusion 1
Lect Lectur ure e summ ummary
In this lecture, we will discuss • Derivation of the integral form of the convection-diffusion equation applied to one-dimen one-dimensional sional,, steady-st steady-state ate conve convection ction-diffu -diffusion sion problems: problems: • FVM applied introduction and worked example • Central, upwind, hybrid, power-law, and higher-order differencing schemes • Assessment of the various differencing schemes
2
Integ Integra rall con convection ection-di -diffu ffusio sion n equa equatio tion n
In proble problems ms where where fluid fluid flow flow plays plays an importa important nt role, role, conve convect ction ion will will accomp accompan any y diffusi diffusion. on. Becaus Becausee diffusi diffusion on will will alway alwayss occur occur alongs alongside ide conv convectio ection n in nature nature,, we consider consider in this lecture lecture the combined combined effects effects of convecti convection on and diffusion. diffusion. Recall Recall the integral form of the general conservation equation considered in the last lecture:
∂ ∂t
∆t
CV
ρφdV dt+ dt + ∆t
n · (ρφu)dAdt = dAdt =
A
∆t
n · (Γ grad (Γ grad φ) φ)dAdt+ dAdt +
A
∆t
S φ dV dt
CV
To obtain the steady convection-diffusion equation, drop the unsteady term on the left-hand side and eliminate the time integration. This yields
A
n · (ρφu)dA = dA =
A
n · (Γ grad (Γ grad φ) φ)dA +
S φ dV.
(1)
CV
The term on the left-hand side is the convective term, the first term on the right-hand side is the diffusive term, and the second term on the right-hand side is the source term. Together, they represent flux balance in the control volume. The principal difficulty in discretizing the convective term is calculating the value of the transported property φ at control volume faces and calculating the fluxes of
AERO 4304
2
Lecture 10
across these these boundar b oundaries. ies. In the previous previous lecture, lecture, we used the central central differencin differencingg φ across method to obtain discretized equations for the diffusive and source terms on the righthand side of Eqn. 1, and it worked worked quite well. Howev However, er, the conve convectiv ctivee term is not so simple simple.. In diffusi diffusion, on, the diffusion diffusion process process affects affects the quant quantit ity y of a transpo transporte rted d quantity along its gradient in all directions, whereas convection affects the quantity only in the flow direction. This difference presents itself as a strict upper-limit on the grid size when central differencing is used. The restriction is a function of the relative strengths of convection and diffusion. A final note before we begin the discussion of the different discretization methods: the velocity u in the convective term is the velocity vector at the boundary faces of the control control volume. volume. These These velocitie velocitiess are generally generally unknown; unknown; they a part of the CFD solution solution and what we are intereste interested d in finding. Howev However, er, in this lecture, lecture, we will not discuss their evaluation; we will assume that they are somehow known. The method for computing these velocities will be discussed in Lecture 11.
3
FVM FVM disc discre reti tiza zati tion on of the the con convecti ection on-d -diff iffus usio ion n equation
In the absence of sources, the steady convection-diffusion equation of a property φ in φ in a 1D flow field u field u is governed by
ρuφdA = ρuφdA =
A
Γ
dφ dA dx
(2)
A
and must also satisfy continuity so
ρudA = ρudA = 0.
(3)
A
This equation equation is integrat integrated ed in the 1D control control volume volume shown in Fig. 1 surrounding surrounding node P . neighbouri uring ng nodes nodes are identifi identified ed as W and E and P . The neighbo E and the CV faces are denoted w and and e e.. Integration of Eqn. 2 within this control volume yields
dφ (ρuA) ρuA)e − (ρuA) ρuA)w = ΓA dx
dφ − ΓA dx e
(4)
w
and integration of Eqn. 3 yields (ρuA) ρuA)e − (ρuA) ρuA)w = 0
(5)
To obtain a discretized form of the convection-diffusion equation, we must approximate imate the terms terms in Eqn. 4. It is conveni convenien entt to define define two two variabl ariables es F and D to
AERO 4304
3
Lecture 10
δ xwP
δ x Pe ue
uw w
e P
W
E
δ xwe
Figure 1: Control volume around node P for discretizing the convection-diffusion equation. represent the convective mass flux per unit area and the diffusion conductance at cell faces: Γ (6) F = ρu and D = δx The cell-face values of the variables F and D can be written as F w = (ρu)w , Γw Dw = , δxW P
F e = (ρu)w Γe De = δx P E
For now, let’s assume that A w = A e = A and use central differencing to represent the contribution of the diffusion terms on the right hand side. The integrated convectiondiffusion equation (Eqn. 4) can now be written as F e φe − F w φw = D e (φE − φP ) − Dw (φP − φW )
(7)
and the integrated continuity equation (Eqn. 5) as F e − F w = 0
(8)
Assuming that the velocity field is “somehow” known (see above note) takes care of the F e and F w . In order to solve Eqn. 7 we need to calculate the transported property φ at the e and w faces. Schemes for this purpose are assessed in the following sections.
3.1
Centered differencing
The central differencing scheme has been used already to represent the diffusion terms that appear on the right-hand side of Eqn. 7. It seems logical that we try the same
AERO 4304
4
Lecture 10
approach to discretize the convective term on the left-hand side. For a uniform grid, we can write the cell-face property values as φe = (φP + φE )/2 φw = (φW + φP )/2
(9) (10)
Substituting these into Eqn. 7 yields F e F w (φP + φE ) − (φW + φP ) = De (φE − φP ) − Dw (φP − φW ) 2 2
(11)
which can be rearranged to give
F w Dw + 2
F e + De − 2
F w + (F e − F w ) φP = Dw + 2
F e φW + De − 2
φE
(12) If we identify the coefficients of the φP , φW , and φE terms as aP , aW , and aE , respectively, we write the central differencing expressions for the discretized convectiondiffusion as (13) aP φP = a W φW + aE φE where
F w F e , ae = D e + , and aP = aw + aE + (F e − F w ) 2 2 It is clear that Eqn. 13 for steady convection-diffusion problems takes the same form as the discretized equation for pure diffusion problems. The difference is that now the coefficients contain extra terms to account for convection. To solve a 1D convection-diffusion problem, we write discretized equations in the form of Eqn. 13 for all cells, which results in an algebraic system of equations for the transported property φ. To illustrate the solution of a 1D steady convection-diffusion problem, let’s consider a worked example. aw = D w +
Example 1: A property φ is transported by means of convection and diffusion
through the one-dimensional domain sketched in Fig. 2. Boundary conditions are φ = 1 at x = 0 and φ = 0 and x = L. Using five equally-spaced points and the central differencing scheme, plot the distribution of φ as a function of x for (i) Case 1: u = 0.1 m/s and (ii) Case 2: u = 2.5 m/s. Compare with the analytical solution exp(ρux/Γ) − 1 φ − φ0 = exp(ρuL/Γ) − 1 φL − φ0
(14)
The following data apply: L = 1.0 m, ρ = 1.0 kg/m3 , Γ = 0.1 kg/m/s. Solution: The method of solution is illustrated for the simple grid shown in Fig.
3. The domain is divided into five equally-sized control volumes with δx = 0.2 m.
AERO 4304
5
Lecture 10
u ϕ = 0
ϕ = 1 x = 0
x = L
Figure 2: One-dimensional domain for Example 1. Note that F = ρu, D = Γ/δx, F e = F w = F , and De = Dw = D everywhere. The boundaries are denoted A and B. The discretization of the convection-diffusion given in Eqn. 13 is applicable in interior nodes 2, 3, 4, but special treatment of the boundary volumes is needed because of they are adjacent to the domain boundaries. We use central differencing for both the diffusion terms and the convective flux through the east face of volume 1. The value of φ on the west face of control volume 1 is given ( φw = φ A = 1) so we do not need to make any approximations in the convective flux term at this boundary. This yields the following equation for control volume 1: F e (φP + φE ) − F AφA = De (φE − φP ) − DA(φP − φA). 2
(15)
For control volume 5, the value of φ on the east face is known (φw = φB = 0). We obtain F w (φP + φW ) = DB (φB − φP ) − Dw (φP − φW ). (16) F B φB − 2 Rearrangement of the above two equations for CVs 1 and 5, noting that D A = D B = 2Γ/δx = 2D and F A = F B = F , gives the discretized equations at the boundary nodes as (17) aP φP = aW φW + aE φE + S u where aP = aW + aE + (F e − F w ) − S P
(18)
and the coefficients a P , aW , a E are as given in Table 3.1. Here, to introduce the boundary conditions, we have suppressed the link to the boundary side (hence the zero aw term for CV 1 and the zero ae term in CV 5). The boundary flux is entered through the source terms. We can now construct the algebraic system of equations for Case 1 and Case 2. δx u
1
ϕ = 1
u
A
5
ϕ = 0 B
x
x
Figure 3: Grid used for discretiztion in Example 1.
AERO 4304
6
Lecture 10
Node aW aE S p 1 0 D − F/2 −(2D + F ) 2,3,4 D + F/2 0 D − F/2 5 D + F/2 0 −(2D − F ) (2D − F )φB
S u (2D + F )φA 0
Table 1: Coefficients for Case 1 of Example 1. Node aW aE S u S p aP = aW + aE − S p 1 0 0.45 1.1φA -1.1 1.55 2 0.55 0.45 0 0 1.0 3 0.55 0.45 0 0 1.0 4 0.55 0.45 0 0 1.0 5 0.55 0 0.5φB -0.9 1.45 Case 1: With u = 0.1 m/s, F = ρu = 0.1, D = Γ/δx = 0.5 and the coefficients
are as given in Table 3.1. The matrix form of the equation set, with φ A = 1 and φB = 0 is
1.55 −0.55 00 0
0 0 0 −0.45 1.0 −0.45 0 0 0 −0.55 1.0 −0.45 0 −0.55 1.0 −0.45 0 0 −0.55 1.45
φ φ φφ
1 2 3 4
φ5
1.1 0 = 00 .
(19)
0
The solution to this equation set is
φ φ φφ
1 2 3 4
φ5
0.9421 0.8006 = 0.6276 0.4163
(20)
0.1579
Comparison to the analytical solution is accomplished by plotting the simulated and analytical distributions of φ in Fig. 4. Considering the coarseness of the grid in this example, reasonable agreement to the analytical solution is achieved using central differencing. Case 2: With u = 2.5 m/s, F = ρu = 2.5, D = Γ/δx = 0.5 gives the coefficients shown in Table 3.1. When these are formed into matrix system of equations and solved as in Case 1, the resulting solution is shown in Fig. 5. In this case, the central differencing scheme produces a solution that appears to oscillate about the exact solution. These oscillations indicate that central differencing of the convective term does not yield an accurate solution under the conditions of Case 2. If Case 2 is repeated with 20 control volumes instead of 5, the results shown in Fig. 6 is obtained. This time, the simulated results do not oscillate and are in much better agreement with the analytical solution.
AERO 4304
Lecture 10
7
Figure 4: Comparison of analytical and simulated solution to Case 1 in Example 1.
Table 2: Coefficients for Case 2 of Example 1. Node aW aE S u S p aP = a W + aE − S p 1 0 -0.75 3.5φA -3.5 2.75 2 1.75 -0.75 0 0 1.0 3 1.75 -0.75 0 0 1.0 4 1.75 -0.75 0 0 1.0 5 1.75 0 -1.5φB 1.5 0.25
Figure 5: Comparison of analytical and simulated solution to Case 2 in Example 1 with the domain discretized with 5 control volumes.
AERO 4304
Lecture 10
8
Figure 6: Comparison of analytical and simulated solution to Case 2 in Example 1 with the domain discretized with 20 control volumes.
3.2
Properties of discretization schemes
The above example shows that central differencing failed under the conditions of Case 2 when the grid was too coarse (5 control volumes) but produced good results when the grid was refined to 20 control volumes. It suggests that the convective term is sensitive to the properties of the discretization scheme. In this section, let’s take a more in-depth look at the properties of discretization schemes. In theory, if we have an infinite number of control volumes, all schemes will achieve results that are indistinguishable from the exact results. But in practice, only a finite number of cells can be used and the results obtained from a particular discretization scheme will only be physically realistic if the scheme possesses certain properties. The most important ones are: • Conservativeness • Boundedness • Transportiveness 3.2.1
Conservativeness
The conservation equations that we have discretized yield an algebraic system of equations that involve the fluxes of the transported property φ through the control volume faces. To ensure conservation of φ for the whole solution domain, the flux of φ leaving the control volume across a certain face must equal the flux of φ entering
AERO 4304
9
Lecture 10
ϕ1
Gradient = (ϕ2 - ϕ1)/δx
ϕ3
ϕ2 q A
1
2
δ x
δ x/ 2
ϕ4
3
δ x
4
δ x
q B
δx/ 2
Figure 7: Example of conservative specification of diffusion fluxes the adjacent control volume through the same face. To achieve this, the flux must be represented by one and the same expression in each control volume. The best way to illustrate this property is with an example. Consider the steady-state 1D diffusion problem without sources shown in Fig. 7. The fluxes across the domain boundaries are q A and q B . Assume that Γ is constant. Let us consider four control volumes and apply central differencing to calculate the diffusive flux across the cell faces. An overall flux balance for the whole domain is obtained by summing the net flux through the each control volume, taking into account the boundary fluxes for the control volumes around nodes 1 and 4:
(φ − φ ) (φ − φ ) (φ − φ ) Γ − q + Γ −Γ δx δx (φ − φ ) (φ − φ ) (φ δx− φ ) 2
1
3
2
2
1
+
4
3
= 0
A
Γ
4
3
δx
−Γ
3
2
δx
+ q B − Γ
δx
(21)
Since Γ is constant, the fluxes across control volume faces cancel out, and the expression simplifies to q B − q A = 0. Therefore, Eqn. 21 express overall conservation of φ throughout the whole domain, as well as conservation across each control boundary. The scheme is there conservative . Non-conservative schemes can occur when the flux across control volume faces are interpolated inconsistently on either side of the face. For example, consider Fig. 8. The variation of φ in control volume 2 is represented by a quadratic function is fitted based on the values of φ at 1, 2, and 3, and similarly, a quadratic function based on the values at 2, 3, and 4 is used to represent the variation of φ in control volume 3. Fig. 8 shows that the resulting quadratic functions can be quite different, and that the gradients at the face between volumes 2 and 3 are unequal. If this is the case, the two fluxes are not equal and so they do not cancel out when summed. Therefore, overall conservation is not satisfied, and the scheme is non-conservative . Figure 8 does not imply that all quadratic interpolation schemes are always non-conservative; later in this lecture, the QUICK scheme will be described; it is a popular quadratic discretization scheme that is consistent.
AERO 4304
10
Lecture 10
Gradient of function 3
Gradient of function 2
Quadratic function 3
Quadratic function 2
ϕ1
ϕ3
ϕ4
3
4
ϕ2 q A
1
2
δx
δx/ 2
δx
δx
q B
δx/ 2
Figure 8: Example of non-conservative specification of diffusion fluxes 3.2.2
Boundedness
For reasons that were described earlier in the course, the algebraic system of equations that result from the discretization at each nodal point are normally solved with iterative methods. These methods start the solution process from a guessed distribution of φ and performs successive updates until a converged solution is achieved. Boundedness is a property of the discretization method that determines whether the method will achieve a converged solution. The Scarborough criterion already brought up provides a sufficient condition for convergence:
|a
nb
|aP | ′
|
≤ 1 for all grid cells < 1 for at least one grid cell
where a P is the net coefficient of the central node P (i.e. aP − S p ) and the summation operator is taken over the neighbouring nodes ( nb). If the differencing scheme produces coefficients that satisfy the above criterion, the resulting matrix is diagonallydominant. To achieve this, the value (aP − S p ) should be large, so it is common practice to linearize the source terms such that S p is always negative and −S p adds to a P . Another essential requirement for boundedness is that all the coefficients must have the same sign (usually all positive). Physically, this means that an increase of φ at one node should produce an increase in φ at the neighbouring nodes. If a scheme does not produce coefficients with all the same sign, it is possible to contain “wiggles” in the solution, or not converge at all. Such was the case in Case 2 of Example 1; most of the east coefficients were negative (see Table 3.1) and the solution contained large undershoots and overshoots. ′
AERO 4304
11
Lecture 10
Pe = 0 Pe > 0 Direction of flow
Pe
W
P
8
E
Figure 9: Distribution of φ in the vicinity of a source located at P for different Peclet numbers 3.2.3
Transportiveness
The transportiveness of a discretization scheme relates to how well the directionality of influence in the physical problem is represented in the discretization scheme. Transportiveness can be illustrated by considering a constant source of φ at a point P as shown in Fig. 9. We define the non-dimensional cell Peclet number as a measure of the relative strengths of convection and diffusion : P e =
F ρu = Γ/δx D
(22)
where δx is a characteristic length (e.g. cell width). The contour lines shown in Fig. 9 indicate the general shape of contours of constant φ for different values of the Peclet number. To illustrate the effect of the Peclet number on the influence of the upstream node P on the downstream node E , let’s consider two extremes: • pure diffusion (P e = 0) • pure convection (P e → ∞) In the case of pure diffusion (P e = 0), the fluid is stagnant and the contours of φ will be concentric circles about P because diffusion will spread φ evenly in all directions. Conditions at E will be influence by P and equally by conditions further downstream. As the Peclet number increases, the contours change shape from circular to elliptical and are shifted in the direction of flow as indicated in Fig. 9. Influence is strongly biased to the upstream direction, so that E is strongly affected by conditions at P , but there is a weak or no influence of E on the conditions at P . In the case of pure convection (P e → ∞), the elliptical contours are completely stretched out in the flow direction. All of property φ emanating from P is convected downstream towards E . Thus the conditions at E are influenced only by the upstream conditions at P , and because there is no diffusion, φE is equal to φP . For a discretization scheme to have the transportiveness property, the relationship between the Peclet number and the directionality of influence must be present within the scheme.
AERO 4304
4 4.1
Lecture 10
12
Assessment of convection-diffusion discretization schemes Centered differencing
Conservativeness. Centered differencing was used in Section 3.1 to discretize the
convection-diffusion equation. The discussion in Section 3.2.1 clearly indicates that the expressions for evaluating the convective and diffusive fluxes at the control faces are consistent, and hence the scheme is conservative. Boundedness. The internal coefficients of the discretized transport equations are
aW = Dw +
F w F e , aE = De − , aP = aW + aE + (F e − F w ). 2 2
A steady one-dimensional flow field is also goverend by the discretized continuity equation, which states that (F e − F w ) is zero when the flow field satisfies continuity. Thus the expression for aP becomes equal to aP = aW +aE . Therefore, the coefficients of the central differencing scheme satisfy the Scarborough cirterion (Eqn. 22). With aE = D e −F e /2 the convective contribution to the east coefficient is negative; if the convection domainates it is possible for aE to be negative, which violates the conditions required for boundedness. Given that F w > 0 and F e > 0 (i.e. the flow is unidirectional), for a E to be positive and the scheme to bounded, D e and F e must satisfy the following condition: F e /De = P ee < 2
(23)
If P ee is greater than 2, the east coefficient will be negative and the solution will be unbounded, leading to oscillations and a physically impossible solution. In the example in Section 3.1, Case 2 takes a Peclet number of 5, so the above condition is violated. The consequences were evident in the results in Fig. 5, which shows “undershoots” and “overshoots” in the solution. When the number of control volumes is increased to 20, the Peclet number reduces to below 2 and the discretization gives answers close to the analytical solution. Transportiveness. The central differencing scheme introduces influencing at node
P from the directions of all its neighbours to calculate the convective and diffusive flux. Thus the scheme does not recognize the direction of flow or the strength of convection relative to diffusion. Therefore, it does not possess the transportiveness property at high Peclet numbers. Accuracy. The Taylor series truncation error of the central differencing scheme is
second order. The requirement for positive coefficients in the central differencing scheme implies that the scheme will be stable and accurate only if P e = F/D < 2.
AERO 4304
Lecture 10
13
It is important to note that the cell Peclet number is a combination of fluid properties (ρ and Γ), and flow property (u), and a property of the computational grid (δx). So, for given values of ρ and Γ, it is only possible to satisfy the restriction given by Eqn. 23 if the velocity is small, as in diffusion-dominated low-Reynolds number flows, or if the grid spacing is small. Owing to this limitation, central differencing is not normally a suitable discretization practice for general-purpose flow calculations, and other discretization schemes have been proposed that possess more favourable properties. These are discussed presently.
4.2
Upwind differencing
One of the main deficiencies of the centered differencing scheme is its lack of transportiveness; that is, its inability to identify flow direction. The value of property φ at a west cell face is always influenced by both φP and φW in centered differencing. However, in a strongly convective flow from west to east, such a treatment is unsuitable because the west face should receive much strongly influencing from node W than from node P , as was illustrated in Fig. 9. The upwind differencing scheme takes the flow direction into account when determining the value of φ at a cell face by setting it equal to the value of the upstream node. The definition of “upstream” varies depending on the value of ue and uw , which results in coefficients that are functions of the values of F e and F w . Compactly, the coefficients for the upwind differencing method applied to the one-dimensional convection-diffusion equation are aP = aW + aE + (F e − F w ) aW = D w + max(F w , 0) aE = D e + max(0, −F e )
(24)
Conservativeness. The upwind differencing scheme uses consistent expressions to
calculate fluxes at the cell faces, and can be easily shown to be conservative. Boundedness. The coefficients of the discretized equation are always positive and
satisfy the requirements for boundedness. When the flow satisfies the continuity equation, the (F e −F w ) term in the central coefficient term aP is zero and so aP = aE +aW , which is desirable for stable iterations. All the coefficients are positive and the coefficient matrix is diagonally-dominant (i.e. it satisfies the Scarborough criterion), and hence no oscillations or “wiggles” occur in the solution. Transportiveness. The scheme considers the directionality of flow, so transportive-
ness is built into the formulation. Accuracy. Because upwind differencing is based on the backwards differencing for-
mula, it is only first-order accurate on the basis of the Taylor series truncation error.
AERO 4304
14
Lecture 10
The resulting error has a diffusion-like appearance and is referred to as false diffusion or numerical diffusion . In practice, grid refinement can reduce the amount of numerical diffusion in an upwind differencing scheme, but the level of refinement that is required to reduce the error to an acceptable level can be prohibitively expensive. Therefore, upwind differencing is not entirely suitable for accurate flow calculations. In certain cases, the numerical diffusion inherent in upwind differencing schemes can be exploited to “damp” the solution and ease convergence, particularly in flows that are sensitive to the initial conditions. If oscillatory-type errors are preventing convergence of a particular flow with a higher-order accuracy discretization scheme, switching to a first-order-accurate scheme such as upwind differencing can sometimes damp the oscillations and allow the solution to converge to a “smoother” flow field, after which a higher-order accuracy scheme can be resumed to obtain a final solution with higher overall accuracy. This approach is not always applicable, however, especially if the solution field is transient.
4.3
Hybrid differencing
The hybrid differencing scheme was developed to combine the higher accuracy of the centered differencing scheme and the improved boundedness and transportiveness properties of the upwind differencing scheme. The central differencing scheme is employed in control volumes with small local Peclet numbers ( P e < 2), while the upwind scheme is used for control volumes with larger Peclet numbers ( P e ≥ 2). Compactly, the coefficients for the hybrid differencing method applied to the onedimensional convection-diffusion equation are aP = aW + aE + (F e − F w ) aW
F = D + max F , D + ,0 2 F . w
w
w
w
aE = D e + max −F e , De −
e
2
(25)
,0
Conservativeness. The hybrid differencing scheme exploits the favourable proper-
ties of the upwind and centered differencing schemes. Because both of those schemes are conservative, the hybrid scheme is also conservative. It switches to upwind differencing when centered differencing produces inaccurate results at high Peclet numbers. Boundedness. Since the coefficients of the hybrid scheme are always positive, it is
unconditionally bounded. Transportiveness. The hybrid scheme possess the transportiveness property by us-
ing an upwind formulation for large values of cell Peclet number.
AERO 4304
15
Lecture 10
Accuracy. Because hybrid differencing is conservative, bounded, and possesses the
transportiveness property, it has been widely used in various CFD procedures and has been very useful in predicting practical engineering flows. However, its disadvantage is that it suffers from lower accuracy; in terms of the Taylor series truncation error, it is only first-order accurate.
4.4
Power-law differencing
Because the hybrid scheme is only first-order accurate, alternate schemes were developed that retain its desirable features while improving the accuracy. The power-law scheme is a more accurate approximation to the one-dimensional exact solution of the convection-diffusion equation and hence produces better results than the hybrid scheme. In the power-law scheme, the diffusion is set to zero when the cell Peclet number exceeds 10. If 0 < P e < 10, the flux is evaluated by using a polynomial expression. For example, the net flux per unit area at the west control volume face is evaluated using q w = F w [φW − β w (φP − φW )] for 0 < Pe < 10
(26)
where β w = (1 − 0.1P ew )5 /P ew and q w = F w φW for P e > 10.
(27)
The coefficients of the discretized equation using the power-law scheme for the onedimensional convection-diffusion equation are given by aP = aW + aE + (F e − F w ) aW = D w · max 0, (1 − 0.1|P ew |)5 + max[F w , 0] . aE
= D · max 0, (1 − 0.1|P e |) + max[−F , 0]
(28)
5
e
e
e
Properties of the power-law scheme. The power-law differencing scheme has
similar properties as the hybrid scheme. The power-law scheme is more accurate for one-dimensional problems because it attempts to represent the exact solution more closely. The scheme has been used in practical flow calculations as an alternative for the hybrid scheme. It is available for use in some commercial codes, such as FLUENT version 4.22, which uses the power-law scheme as the default for flow calculations.
4.5
Higher-order differencing schemes: the QUICK method
The hybrid, power-law, and upwind schemes are only first-order accurate in terms of the Taylor series truncation error. The use of upwind quantities ensures that the schemes are very stable, bounded, and obey the transportiveness requirements, but
AERO 4304
16
Lecture 10
also makes them suffer from numerical diffusion-type errors. These errors can be minimised by employing higher-order methods of discretization. To achieve higher-order accuracy, the computational stencil needs to be enlarged by incorporating neighbouring points and bringing in a wider influence onto the point of interest. The centered differencing scheme does this but it is not bounded for large Peclet numbers and does not possess the transportiveness property. Therefore, alternate formulations are sought that are more accurate than upwind schemes but still preserve the stability and sensitivity to flow direction. A widely-used approach is the quadratic upwind differencing scheme developed by Leonard (1979), which is termed the quadratic upstream interpolation for convective kinetics (QUICK) scheme. This scheme uses a three-point upstream-weighted quadratic interpolation for cell face values. The face value of φ is obtained from a quadratic function passing through two bracketing nodes on each side of the face plus a node on the upstream side, as shown in Fig. 10. The resulting expression for the face values is 6 3 1 (29) φface = φi 1 + φi − φi 2 8 8 8 For example, when uw > 0 and ue > 0, a quadratic fit through W W , W , and P is used to evaluate φ w , and a quadratic fit through W , P , and E is used to evaluate φ e . For uw < 0 and u e < 0, a quadratic fit through W , P , and E is used to evaluate φ w , and a quadratic fit through P , E , and EE is used to evaluate φe . The diffusion terms may be evaluated using the gradient of the appropriate parabola. It is interesting to note that on a uniform grid this practice gives the same expressions as central differencing for diffusion. General expressions of the QUICK scheme for one-dimensional convectiondiffusion problems valid for positive and negative flow directions is summarized as follows: (30) a pφP = aW φW + aE φE + aW W φW W + aEE φEE −
−
with a central coefficient aP = aW + aE + aW W + aEE + (F e − F w )
(31)
and neighbour coefficients 6 1 3 aW = D w + αw F w + αe F e + (1 − αw )F w 8 8 8 1 aW W = − αw F w 8 . 3 6 1 aE = D e − αe F e − (1 − αe )F e − (1 − αw )F w 8 8 8 1 aEE = (1 − αe )F e 8
(32)
where αw = 1 for F w > 0 and αe = 1 for F e > 0 αw = 0 for F w < 0 and αe = 0 for F e < 0
(33)
AERO 4304
Lecture 10
17
Figure 10: Schematic illustrating the quadratic interpolation between nodes in the QUICK scheme Properties of the QUICK scheme. The scheme uses consistent quadratic profiles—
the cell face values of fluxes are always calculated by quadratic interpolation between two bracketing nodes and an upstream node—and is therefore conservative. Since the scheme is based on a quadratic function, its accuracy in terms of the Taylor series truncation error is third-order on a uniform mesh. The transportiveness property is built into the scheme as well, since the quadratic function is based on two upsteam and one downstream nodal values. Additionally, if the flow field satisfies continuity, the coefficient aP equals the sum of all neighbour coefficients, which is desirable for boundedness. There are certain undesirable features of the QUICK scheme. Firstly, the main coefficients (E and W ) are not guaranteed to be positive and the coefficients aW W
AERO 4304
18
Lecture 10
and aEE are negative. For example, if u w > 0 and u e > 0, the east coefficient becomes negative at relatively modest cell Peclet numbers (P ee = F e /De > 8/3), which gives rise to stability problems and unbounded solutions under certain flow conditions. Similarly the west coefficient can become negative when the flow is in the negative direction. The QUICK scheme is therefore conditionally stable. Secondly, because the discretized equations involve not only immediate-neighbour nodes but also nodes further away, the algebraic system of equations is not triadiagonal and the efficient tridiagonal-matrix methods are not directly applicable. 4.5.1
Reformulated QUICK schemes
Because the QUICK scheme as presented above can be unstable due to the appearance of negative main coefficients, it has been reformulated in various different ways to alleviate the stability constraints. For the most part, the reformulations place negative coefficients in the source term so as to retain positive main coefficients. The contributing part is appropriately weighted to give better stability and positive coefficients as far as possible. Hayase et al. (1992) generalized the approach for re-arranging QUICK schemes and derived a stable and fast-converging variant. Their scheme can be summarized as follows: 1 φw = φ W + [3φP − 2φW − φW W ] for F w > 0 8 1 φe = φ E + [3φE − 2φP − φW ] for F e > 0 8 (34) 1 φw = φ W + [3φW − 2φP − φE ] for F w < 0 8 1 φe = φ E + [3φP − 2φE − φEE ] for F e < 0 8 The discretized equation takes the form aP φP = aW φW + aE φE + S
(35)
where the central coefficient is aP = aW + aE + (F e − F w )
(36)
and the neighbouring coefficients and source term are aW = D w + αw F w aE = D e − (1 − αe )F e 1 1 . (37) S = (3φP − 2φW − φW W )αw F w + (φW + 2φP − 3φE )αe F e 8 8 1 1 + (3φW − 2φP − φE )(1 − αw )F w + (2φE + φEE − 3φP )(1 − αe )F e 8 8 The advantage of this approach is that the coefficients are always positive and now satisfy the requirements for conservativeness, boundedness, and transportiveness.
AERO 4304
Lecture 10
19
Figure 11: Comparison of QUICK and upwind differencing schemes against the exact solution of the convection-diffusion equation 4.5.2
General comments on the QUICK differencing scheme
The QUICK differencing scheme has greater formal accuracy than the central differencing or hybrid schemes while also retaining the upwind-weighted characteristics that make it conservative and possess transportiveness property. The resultant numerical diffusion is small and solutions achieved with coarse grids are often considerably more accurate than those achieved with upwind or hybrid schemes. Figure 11 shows a comparison between upwind and QUICK against the exact solution of the convection-diffusion equation. It can be seen that the QUICK scheem matches the exact solution much more accurately than the upwind scheme on a 50 × 50 grid. Figure 11 also illustrates that the QUICK scheme can give minor undershoots and overshoots in the solution. In complex flow calculations, such oscillations can
AERO 4304
Lecture 10
20
cause QUICK schemes to give unbounded results, and such a possibility needs to be considered when interpreting solutions. Classes of QUICK schemes with limiter functions that are specially formulated to achieve oscillation-free solutions have been proposed but are beyond the scope of this course. Interested readers can refer to Hirsch (1990).
References Hayase, T., Humphrey, J.A.C & Greif, R. 1992 A consistenly formulated
QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures. Journal of Computational Physics 98, 108–118. Hirsch, C. 1990 Numerical Computation of Internal and External Flows , , vol. 2.
John Wiley & Sons, Chichester, England. Leonard, B.P. 1979 A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Computer Methods in Applied Mechanics and Engineering 19, 59–98.
