Versteeg H., Malalasekra W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method
Подождите немного. Документ загружается.
5.9 HIGHER-ORDER DIFFERENCING SCHEMES 157
For example, when u
w
> 0 and u
e
> 0 a quadratic fit through WW, W and
P is used to evaluate
φ
w
, and a further quadratic fit through W, P and E to
calculate
φ
e
. For u
w
< 0 and u
e
< 0 values of
φ
at W, P and E are used for
φ
w
,
and values at P, E and EE for
φ
e
. It can be shown that for a uniform grid the
value of
φ
at the cell face between two bracketing nodes i and i − 1 and
upstream node i − 2 is given by the following formula:
φ
face
=
φ
i −1
+
φ
i
−
φ
i −2
(5.45)
When u
w
> 0, the bracketing nodes for the west face w are W and P, the
upstream node is WW (Figure 5.17) and
φ
w
=
φ
W
+
φ
P
−
φ
WW
(5.46)
When u
e
> 0, the bracketing nodes for the east face e are P and E, the
upstream node is W, so
φ
e
=
φ
P
+
φ
E
−
φ
W
(5.47)
The diffusion terms may be evaluated using the gradient of the approximat-
ing parabola. It is interesting to note that on a uniform grid this practice gives
the same expressions as central differencing for diffusion, since the slope
of the chord between two points on a parabola is equal to the slope of the
tangent to the parabola at its midpoint. If F
w
> 0 and F
e
> 0, and if we use
equations (5.46)–(5.47) for the convective terms and central differencing for
the diffusion terms, the discretised form of the one-dimensional convection–
diffusion transport equation (5.9) may be written as
F
e
φ
P
+
φ
E
−
φ
W
− F
w
φ
W
+
φ
P
−
φ
WW
= D
e
(
φ
E
−
φ
P
) − D
w
(
φ
P
−
φ
W
)
which can be rearranged to give
D
w
− F
w
+ D
e
+ F
e
φ
P
= D
w
+ F
w
+ F
e
φ
W
+ D
e
− F
e
φ
E
− F
w
φ
WW
(5.48)
1
8
J
K
L
3
8
G
H
I
J
K
L
1
8
6
8
G
H
I
J
K
L
6
8
3
8
G
H
I
J
K
L
D
E
F
1
8
3
8
6
8
A
B
C
D
E
F
1
8
3
8
6
8
A
B
C
G
H
I
1
8
3
8
6
8
1
8
3
8
6
8
1
8
3
8
6
8
Figure 5.17 Quadratic profiles
used in the QUICK scheme
ANIN_C05.qxd 29/12/2006 04:36PM Page 157
158 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS
This is now written in the standard form for discretised equations:
a
P
φ
P
= a
W
φ
W
+ a
E
φ
E
+ a
WW
φ
WW
(5.49)
where
a
W
a
E
a
WW
a
P
D
w
+ F
w
+ F
e
D
e
− F
e
− F
w
a
W
+ a
E
+ a
WW
+ (F
e
− F
w
)
For F
w
< 0 and F
e
< 0 the flux across the west and east boundaries is given
by the expressions
φ
w
=
φ
P
+
φ
W
−
φ
E
(5.50)
φ
e
=
φ
E
+
φ
P
−
φ
EE
Substitution of these two formulae for the convective terms in the discretised
convection–diffusion equation (5.9) together with central differencing for
the diffusion terms leads, after rearrangement as above, to the following
coefficients:
a
W
a
E
a
EE
a
P
D
w
+ F
w
D
e
− F
e
− F
w
F
e
a
W
+ a
E
+ a
EE
+ (F
e
− F
w
)
General expressions, valid for positive and negative flow directions, can be
obtained by combining the two sets of coefficients above.
The QUICK scheme for one-dimensional convection–diffusion
problems can be summarised as follows:
a
P
φ
P
= a
W
φ
W
+ a
E
φ
E
+ a
WW
φ
WW
+ a
EE
φ
EE
(5.51)
with central coefficient
a
P
= a
W
+ a
E
+ a
WW
+ a
EE
+ (F
e
− F
w
)
and neighbour coefficients
a
W
a
WW
a
E
a
EE
D
w
+
α
w
F
w
+
α
e
F
e
−
α
w
F
w
D
e
−
α
e
F
e
− (1 −
α
e
)F
e
(1 −
α
e
)F
e
+ (1 −
α
w
)F
w
− (1 −
α
w
)F
w
1
8
3
8
1
8
6
8
3
8
1
8
1
8
6
8
1
8
1
8
6
8
3
8
1
8
3
8
6
8
1
8
3
8
6
8
1
8
3
8
1
8
6
8
ANIN_C05.qxd 29/12/2006 04:36PM Page 158
5.9 HIGHER-ORDER DIFFERENCING SCHEMES 159
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
Using the QUICK scheme solve the problem considered in Example 5.1 for
u = 0.2 m/s on a five-point grid. Compare the QUICK solution with the
exact and the central differencing solution.
As before, the five-node grid introduced in Example 5.1 is used for the
discretisation. With the data of this example and u = 0.2 m/s we have F = F
e
= F
w
= 0.2 and D = D
e
= D
w
= 0.5 everywhere so that the cell Peclet number
becomes Pe
w
= Pe
e
=
ρ
u
δ
x/Γ=0.4. The discretisation equation with the
QUICK scheme at internal nodes 3 and 4 is given by (5.51) together with its
coefficients.
In the QUICK scheme the
φ
-value at cell boundaries is calculated with
formulae (5.46)–(5.47) that use three nodal values. Nodes 1, 2 and 5 are
all affected by the proximity of domain boundaries and need to be treated
separately. At the boundary node 1,
φ
is given at the west (w) face (
φ
w
=
φ
A
),
but there is no west (W) node to evaluate
φ
e
at the east face by (5.47). To
overcome this problem Leonard (1979) suggested a linear extrapolation to
create a ‘mirror’ node at a distance
δ
x/2 to the west of the physical bound-
ary. This is illustrated in Figure 5.18.
Figure 5.18 Mirror node
treatment at the boundary
Example 5.4
Solution
It can be easily shown that the linearly extrapolated value at the mirror
node is given by
φ
0
= 2
φ
A
−
φ
P
(5.52)
The extrapolation to the ‘mirror’ node has given us the required W node for
the formula (5.47) that calculates
φ
e
at the east face of control volume 1:
φ
e
=
φ
P
+
φ
E
− (2
φ
A
−
φ
P
)
=
φ
P
+
φ
E
−
φ
A
(5.53)
At the boundary nodes the gradients must be evaluated using an expression
consistent with formula (5.53). It can be shown that the diffusive flux through
the west boundary is given by
2
8
3
8
7
8
1
8
3
8
6
8
ANIN_C05.qxd 29/12/2006 04:36PM Page 159
160 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS
Γ
A
= (9
φ
P
− 8
φ
A
−
φ
E
) (5.54)
where D
A
* =
The superscript * is used to indicate that, in the QUICK scheme, the diffu-
sion conductances at boundary nodes and interior nodes have the same value,
i.e. D
A
* = D =Γ/
δ
x. This aspect is different from the discretisation schemes
we have discussed thus far. These used the half-cell approximation, so the
diffusive conductance at the boundary cell was always D
A
= 2D = 2Γ/
δ
x.
The discretised equation at node 1 is
F
e
φ
P
+
φ
E
−
φ
A
− F
A
φ
A
= D
e
(
φ
E
−
φ
P
) − (9
φ
P
− 8
φ
A
−
φ
E
) (5.55)
At control volume 5, the
φ
-value at the east face is known (
φ
e
=
φ
B
) and the
diffusive flux of
φ
through the east boundary is given by
Γ
B
= (8
φ
B
− 9
φ
P
+
φ
W
) (5.56)
where D
B
* =
At node 5 the discretised equation becomes
F
B
φ
B
− F
w
φ
W
+
φ
P
−
φ
WW
= (8
φ
B
− 9
φ
P
+
φ
W
) − D
w
(
φ
P
−
φ
W
) (5.57)
Since a special expression is used to evaluate
φ
at the east face of control
volume 1 we must use the same expression for
φ
to calculate the convective
flux through the west face of control volume 2 to ensure flux consistency. So
at node 2 we have
F
e
φ
P
+
φ
E
−
φ
W
− F
w
φ
W
+
φ
P
−
φ
A
= D
e
(
φ
E
−
φ
P
) − D
w
(
φ
P
−
φ
W
) (5.58)
The discretised equations for nodes 1, 2 and 5 are now written to fit into the
standard form to give
a
P
φ
P
= a
WW
φ
WW
+ a
W
φ
W
+ a
E
φ
E
+ S
u
(5.59)
with
a
P
= a
WW
+ a
W
+ a
E
+ (F
e
− F
w
) − S
P
and
J
K
L
2
8
3
8
7
8
G
H
I
J
K
L
1
8
3
8
6
8
G
H
I
D*
B
3
J
K
L
1
8
3
8
6
8
G
H
I
Γ
δ
x
D*
B
3
∂φ
∂
x
D
A
*
3
J
K
L
2
8
3
8
7
8
G
H
I
Γ
δ
x
D
A
*
3
∂φ
∂
x
ANIN_C05.qxd 29/12/2006 04:36PM Page 160
5.9 HIGHER-ORDER DIFFERENCING SCHEMES 161
Table 5.10
Node a
W
a
E
a
WW
S
u
S
P
a
P
1 0 0.592 0 1.583
φ
A
−1.583 2.175
2 0.7 0.425 0 −0.05
φ
A
0.05 1.075
3 0.675 0.425 −0.025 0 0 1.075
4 0.675 0.425 −0.025 0 0 1.075
5 0.817 0 −0.025 1.133
φ
B
−1.133 1.925
The matrix form of the equation set is
G 2.175 −0.592 0 0 0 JG
φ
1
JG1.583J
H−0.7 1.075 −0.425 0 0 KH
φ
2
KH−0.05 K
H 0.025 −0.675 1.075 −0.425 0 KH
φ
3
K = H 0 K (5.60)
H 0 0.025 −0.675 1.075 −0.425KH
φ
4
KH0 K
I 0 0 0.025 −0.817 1.925LI
φ
5
LI0 L
The solution to the above system is
G
φ
1
JG0.9648J
H
φ
2
KH0.8707K
H
φ
3
K = H0.7309K (5.61)
H
φ
4
KH0.5226K
I
φ
5
LI0.2123L
Comparison with the analytical solution
Figure 5.19 shows that the QUICK solution is almost indistinguishable from
the exact solution. Table 5.11 confirms that the errors are very small even
with this coarse mesh. Following the steps outlined in Example 5.1 the cen-
tral differencing solution is computed with the data given above. The sum of
absolute errors in Table 5.11 indicates that the QUICK scheme gives a more
accurate solution than the central differencing scheme.
Node a
WW
a
W
a
E
S
P
S
u
10 0 D
e
+ D
A
* − F
e
− D
A
* + F
e
+ F
A
D
A
* + F
e
+ F
A
φ
A
20 D
w
+ F
w
+ F
e
D
e
− F
e
F
w
− F
w
φ
A
5 − F
w
D
w
+ D
B
* + F
w
0 − D
B
* − F
B
D
B
* − F
B
φ
B
Substitution of numerical values gives the coefficients summarised in
Table 5.10.
D
E
F
8
3
A
B
C
D
E
F
8
3
A
B
C
6
8
1
3
1
8
1
4
1
4
3
8
1
8
7
8
D
E
F
2
8
8
3
A
B
C
D
E
F
2
8
8
3
A
B
C
3
8
1
3
ANIN_C05.qxd 29/12/2006 04:36PM Page 161
5.9.2 Assessment 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 Taylor series
truncation error is third-order on a uniform mesh. The transportiveness
property is built into the scheme as the quadratic function is based on two
upstream and one downstream nodal values. If the flow field satisfies con-
tinuity the coefficient a
P
equals the sum of all neighbour coefficients, which
is desirable for boundedness.
On the downside, the main coefficients (E and W) are not guaranteed
to be positive and the coefficients a
WW
and a
EE
are negative. For example, if
u
w
> 0 and u
e
> 0 the east coefficient becomes negative at relatively modest
cell Peclet numbers (Pe
e
= F
e
/D
e
> 8/3). This 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.
Another notable feature is the fact that the discretised equations involve
not only immediate-neighbour nodes but also nodes further away. Tri-diagonal
matrix solution methods (see Chapter 7) are not directly applicable.
162 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS
Table 5.11
Node Distance Analytical QUICK Difference CD Difference
solution solution solution
1 0.1 0.9653 0.9648 0.0005 0.9696 0.0043
2 0.3 0.8713 0.8707 0.0006 0.8786 0.0073
3 0.5 0.7310 0.7309 0.0001 0.7421 0.0111
4 0.7 0.5218 0.5226 −0.0008 0.5374 0.0156
5 0.9 0.2096 0.2123 −0.0027 0.2303 0.0207
∑ Absolute error 0.0047 0.059
Figure 5.19 Comparison of
QUICK solution with the
analytical solution
ANIN_C05.qxd 29/12/2006 04:36PM Page 162
5.9 HIGHER-ORDER DIFFERENCING SCHEMES 163
5.9.3 Stability problems of the QUICK scheme and remedies
Since the QUICK scheme in the form presented above can be unstable due
to the appearance of negative main coefficients, it has been reformulated
in different ways that alleviate stability problems. These formulations all
involve placing troublesome 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.
Some of the better known practical approaches are described in Han et al.
(1981), Pollard and Siu (1982) and Hayase et al. (1992). The last authors
generalised the approach for rearranging QUICK schemes and derived a
stable and fast converging variant.
The Hayase et al. (1992) QUICK scheme can be summarised as follows:
φ
w
=
φ
W
+ [3
φ
P
− 2
φ
W
−
φ
WW
] for F
w
> 0
φ
e
=
φ
P
+ [3
φ
E
− 2
φ
P
−
φ
W
] for F
e
> 0 (5.62)
φ
w
=
φ
P
+ [3
φ
W
− 2
φ
P
−
φ
E
] for F
w
< 0
φ
e
=
φ
E
+ [3
φ
P
− 2
φ
E
−
φ
EE
] for F
e
< 0
The discretisation equation takes the form
a
P
φ
P
= a
W
φ
W
+ a
E
φ
E
+ D (5.63)
The central coefficient is
a
P
= a
W
+ a
E
+ (F
e
− F
w
)
and
a
W
a
E
D
D
w
+
α
w
F
w
D
e
− (1 −
α
e
)F
e
(3
φ
P
− 2
φ
W
−
φ
WW
)
α
w
F
w
+ (
φ
W
+ 2
φ
P
− 3
φ
E
)
α
e
F
e
+ (3
φ
W
− 2
φ
P
−
φ
E
)(1 −
α
w
)F
w
+ (2
φ
E
+
φ
EE
− 3
φ
P
)(1 −
α
e
)F
e
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
1
8
1
8
1
8
1
8
1
8
1
8
1
8
1
8
ANIN_C05.qxd 29/12/2006 04:36PM Page 163
The advantage of this approach is that the main coefficients are positive
and satisfy the requirements for conservativeness, boundedness and trans-
portiveness. The allocation to the source term of the part of the discretisa-
tion that contains negative coefficients is called deferred correction and
relies on the scheme being applied as part of an iterative loop structure. At
the nth iteration the source term is evaluated using values known at the end
of the previous (n − 1)th iteration, i.e. ‘correction’ of the main coefficients is
‘deferred’ by one iteration. After a sufficiently large number of iterations
the correction ‘catches up’ with the rest of the solution, so all variations
of QUICK, including the one developed by Hayase et al., will give the same
converged solution.
5.9.4 General comments on the QUICK differencing scheme
The QUICK differencing scheme has greater formal accuracy than the
central differencing or hybrid schemes, and it retains the upwind-weighted
characteristics. The resultant false diffusion is small, and solutions achieved
with coarse grids are often considerably more accurate than those of the
upwind or hybrid schemes. Figure 5.20 shows a comparison between
upwind and QUICK for the two-dimensional test case considered in section
5.6.1. It can be seen that the QUICK scheme matches the exact solution
much more accurately than the upwind scheme on a 50 × 50 grid.
164 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS
Figure 5.20 Comparison of
QUICK and upwind solutions
for the 2D test case considered in
section 5.6.1
TVD schemes
The QUICK scheme can, however, give (minor) undershoots and over-
shoots, as is evident in Figure 5.20. In complex flow calculations, the use of
QUICK can lead to subtle problems caused by such unbounded results: for
example, they could give rise to negative turbulence kinetic energy (k) in k–
ε
model (see Chapter 3) computations. The possibility of undershoots and
overshoots needs to be considered when interpreting solutions.
Schemes of third-order and above have been developed for the discretisation
of convective terms with varying degrees of success. Implementation of
boundary conditions can be problematic with such higher-order schemes.
The fact that the QUICK scheme and other higher-order schemes can give
undershoots and overshoots has led to the development of second-order
5.10
ANIN_C05.qxd 29/12/2006 04:36PM Page 164
5.10 TVD SCHEMES 165
schemes that avoid these problems. The class of TVD (total variation dimin-
ishing) schemes has been specially formulated to achieve oscillation-free
solutions and has proved to be useful in CFD calculations. TVD is a
property used in the discretisation of equations governing time-dependent
gas dynamics problems. More recently, schemes with this property have
also become popular in general-purpose CFD solvers. Fundamentals of the
development of TVD methodology involves a fair amount of mathematical
background. However, the ideas behind TVD schemes can be easily illus-
trated in the context of the discretisation practices presented in the previous
sections by considering the basic properties of standard schemes and their
deficiencies.
As discussed earlier, the basic upwind differencing scheme is the most
stable and unconditionally bounded scheme, but it introduces a high level of
false diffusion due to its low order of accuracy (first-order). Higher-order
schemes such as central differencing and QUICK can give spurious oscilla-
tions or ‘wiggles’ when the Peclet number is high. When such higher-order
schemes are used to solve for turbulent quantities, e.g. turbulence energy
and rate of dissipation, wiggles can give physically unrealistic negative values
and instability. TVD schemes are designed to address this undesirable oscil-
latory behaviour of higher-order schemes. In TVD schemes the tendency
towards oscillation is counteracted by adding an artificial diffusion fragment
or by adding a weighting towards upstream contribution. In the literature
early schemes based on these ideas were called flux corrected transport
(FCT) schemes: see Boris and Book (1973, 1976). Further work by Van Leer
(1974, 1977a,b, 1979), Harten (1983, 1984), Sweby (1984), Roe (1985), Osher
and Chakravarthy (1984) and many others has contributed to the develop-
ment of present-day TVD schemes. In the next section we explain the funda-
mentals of the TVD methodology.
5.10.1 Generalisation of upwind-biased discretisation schemes
Consider the standard control volume discretisation of the one-dimensional
convection–diffusion equation (5.3). Discretisation of the diffusion terms
using the central differencing practice is standard and does not require
any further consideration. It is the discretisation of the convective flux term
that requires special attention. We assume that the flow is in the positive
x-direction, so u > 0, and develop the TVD concept as a generalisation of
upwind-biased expressions for the value of transported quantity
φ
at the east
face of a one-dimensional control volume.
The standard upwind differencing (UD) scheme for the east face value of
φ
e
gives
φ
e
=
φ
P
(5.64)
A linear upwind differencing (LUD) scheme, which involves two upstream
values, yields the following expression for
φ
e
:
φ
e
=
φ
P
+
=
φ
P
+ (
φ
P
−
φ
W
) (5.65)
1
2
δ
x
2
(
φ
P
−
φ
W
)
δ
x
ANIN_C05.qxd 29/12/2006 04:36PM Page 165
This can be thought of as a second-order extension of the original UD
estimate (5.64) of
φ
e
with a correction based on an upwind-biased estimate
(
φ
P
−
φ
W
)/
δ
x of the gradient of
φ
multiplied by the distance
δ
x/2 between
node P and the east face. Another way of looking at this is to recall that our
aim is to construct expressions for convective flux F
e
φ
e
. Hence, for positive
flow direction, the convective flux discretisation by means of the LUD scheme
can be thought of as the sum of the basic UD convective flux F
e
φ
P
plus an
additional flux contribution F
e
(
φ
P
−
φ
W
)/2 to improve the order of accuracy.
The QUICK scheme (5.47) can be similarly rearranged in the form of the
UD estimate plus a correction:
φ
e
=
φ
P
+ [3
φ
E
− 2
φ
P
−
φ
W
] (5.66)
The central differencing (CD) scheme can be written as follows:
φ
e
=
=
φ
P
+ (
φ
E
−
φ
P
) (5.67)
We consider a generalization of the higher-order schemes in the following
form:
φ
e
=
φ
P
+
ψ
(
φ
E
−
φ
P
) (5.68)
where
ψ
is an appropriate function.
In choosing this form we express the convective flux at the east face as the
sum of the flux F
e
φ
P
that is obtained when we use UD and an additional con-
vective flux F
e
ψ
(
φ
E
−
φ
P
)/2. The extra contribution is connected in some
way to the gradient of the transported quantity
φ
at the east face, as indicated
by its central difference approximation (
φ
E
−
φ
P
). It is easy to see that the
central difference scheme (5.68) leads to function
ψ
= 1, but in sections
5.3–5.5 we have established that an additional convective flux based on this
choice of
ψ
leads to wiggles in the solution if the grid is too coarse due to lack
of transportiveness. The upwind scheme (5.64) corresponds to function
ψ
= 0, but this choice of
ψ
gave rise to false diffusion. Looking at the higher-
order schemes we find that the LUD scheme (5.65) may be rewritten as
φ
e
=
φ
P
+ (
φ
E
−
φ
P
) (5.69)
Hence, for LUD the function is
ψ
= (
φ
P
−
φ
W
)/(
φ
E
−
φ
P
).
After some algebra the QUICK expression (5.66) can be rewritten as
φ
e
=
φ
P
+ 3 + (
φ
E
−
φ
P
) (5.70)
By comparing equation (5.70) with equation (5.68) it can be seen that the
appropriate function for the QUICK scheme is
ψ
= 3 +
1
4
D
E
F
φ
P
−
φ
W
φ
E
−
φ
P
A
B
C
J
K
L
1
4
D
E
F
φ
P
−
φ
W
φ
E
−
φ
P
A
B
C
G
H
I
1
2
D
E
F
φ
P
−
φ
W
φ
E
−
φ
P
A
B
C
1
2
1
2
1
2
(
φ
P
+
φ
E
)
2
1
8
166 CHAPTER 5 FINITE VOLUME METHOD FOR C---D PROBLEMS
ANIN_C05.qxd 29/12/2006 04:36PM Page 166