Versteeg H., Malalasekra W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method
Подождите немного. Документ загружается.
6.4 THE SIMPLE ALGORITHM 187
a
I, j
v*
I, j
=∑a
nb
v*
nb
+ (p *
I, J−1
− p*
I, J
)A
I, j
+ b
I, j
(6.13)
Now we define the correction p′ as the difference between correct pressure
field p and the guessed pressure field p*, so that
p = p* + p′ (6.14)
Similarly we define velocity corrections u′ and v′ to relate the correct velo-
cities u and v to the guessed velocities u* and v*:
u = u* + u′ (6.15)
v = v* + v′ (6.16)
Substitution of the correct pressure field p into the momentum equations
yields the correct velocity field (u, v). Discretised equations (6.8) and (6.10)
link the correct velocity fields with the correct pressure field.
Subtraction of equations (6.12) and (6.13) from (6.8) and (6.10), respec-
tively, gives
a
i, J
(u
i, J
− u*
i, J
) =∑a
nb
(u
nb
− u*
nb
) + [(p
I−1, J
− p*
I−1, J
) − (p
I, J
− p*
I, J
)]A
i, J
(6.17)
a
I, j
(v
I, j
− v*
I, j
) =∑a
nb
(v
nb
− v*
nb
) + [(p
I, J−1
− p*
I, J−1
) − (p
I, J
− p*
I, J
)]A
I, j
(6.18)
Using correction formulae (6.14)–(6.16) the equations (6.17)–(6.18) may be
rewritten as follows:
a
i, J
u′
i, J
=∑a
nb
u ′
nb
+ (p ′
I−1, J
− p ′
I, J
)A
i, J
(6.19)
a
I, j
v ′
I, j
=∑a
nb
v′
nb
+ (p ′
I, J−1
− p ′
I, J
)A
I, j
(6.20)
At this point an approximation is introduced: ∑a
nb
u′
nb
and ∑a
nb
v′
nb
are
dropped to simplify equations (6.19) and (6.20) for the velocity corrections.
Omission of these terms is the main approximation of the SIMPLE
algorithm. We obtain
u ′
i, J
= d
i, J
(p ′
I−1, J
− p′
I, J
) (6.21)
v′
I, j
= d
I, j
(p′
I, J −1
− p ′
I, J
) (6.22)
where d
i , J
= and d
I, j
= (6.23)
Equations (6.21) and (6.22) describe the corrections to be applied to velo-
cities through formulae (6.15) and (6.16), which gives
u
i, J
= u*
i, J
+ d
i, J
(p ′
I−1, J
− p ′
I, J
) (6.24)
v
I, j
= v*
I, j
+ d
I, j
(p ′
I, J −1
− p ′
I, J
) (6.25)
Similar expressions exist for u
i+1, J
and v
I, j+1
:
u
i+1, J
= u*
i +1, J
+ d
i+1, J
(p ′
I, J
− p ′
I+1, J
) (6.26)
v
I, j+1
= v*
I, j+1
+ d
I, j+1
(p ′
I, J
− p ′
I, J+1
) (6.27)
A
I, j
a
I, j
A
i , J
a
i, J
ANIN_C06.qxd 29/12/2006 09:59 AM Page 187
188 CHAPTER 6 ALGORITHMS FOR PRESSURE---VELOCITY COUPLING
Figure 6.5 The scalar
control volume used for the
discretisation of the continuity
equation
where d
i+1, J
= and d
I, j+1
= (6.28)
Thus far we have only considered the momentum equations but, as men-
tioned earlier, the velocity field is also subject to the constraint that it should
satisfy continuity equation (6.3). Continuity is satisfied in discretised form
for the scalar control volume shown in Figure 6.5:
[(
ρ
uA)
i+1, J
− (
ρ
uA)
i, J
] + [(
ρ
vA)
I, j+1
− (
ρ
vA)
I, j
] = 0 (6.29)
A
I, j+1
a
I, j+1
A
i+1, J
a
i+1, J
Substitution of the corrected velocities of equations (6.24)–(6.27) into
discretised continuity equation (6.29) gives
[
ρ
i+1, J
A
i+1, J
(u*
i +1, J
+ d
i+1, J
(p ′
I, J
− p ′
I+1, J
)) −
ρ
i, J
A
i, J
(u*
i, J
+ d
i, J
(p ′
I−1, J
− p ′
I, J
))] + [
ρ
I, j+1
A
I, j+1
(v*
I, j+1
+ d
I, j+1
(p ′
I, J
− p ′
I, J+1
))
−
ρ
I, j
A
I, j
(v*
I, j
+ d
I, j
(p ′
I, J −1
− p ′
I, J
))] = 0 (6.30)
This may be rearranged to give
[(
ρ
dA)
i+1, J
+ (
ρ
dA)
i, J
+ (
ρ
dA)
I, j+1
+ (
ρ
dA)
I, j
]p ′
I, J
= (
ρ
dA)
i+1, J
p ′
I+1, J
+ (
ρ
dA)
i, J
p ′
I−1, J
+ (
ρ
dA)
I, j+1
p ′
I, J+1
+ (
ρ
dA)
I, j
p ′
I, J−1
+ [(
ρ
u*A)
i, J
− (
ρ
u*A)
i+1, J
+ (
ρ
v*A)
I, j
− (
ρ
v*A)
I, j+1
] (6.31)
Identifying the coefficients of p′, this may be written as
a
I, J
p′
I, J
= a
I+1, J
p ′
I+1, J
+ a
I−1, J
p ′
I−1, J
+ a
I, J+1
p ′
I, J+1
+ a
I, J−1
p ′
I, J−1
+ b ′
I, J
(6.32)
where a
I,J
= a
I+1, J
+ a
I−1, J
+ a
I, J+1
+ a
I, J−1
and the coefficients are given below:
a
I+1, J
a
I−1, J
a
I, J+1
a
I, J−1
b ′
I, J
(
ρ
dA)
i+1, J
(
ρ
dA)
i, J
(
ρ
dA)
I, j+1
(
ρ
dA)
I, j
(
ρ
u*A)
i, J
− (
ρ
u*A)
i+1, J
+ (
ρ
v*A)
I, j
− (
ρ
v*A)
I, j+1
ANIN_C06.qxd 29/12/2006 09:59 AM Page 188
6.4 THE SIMPLE ALGORITHM 189
Equation (6.32) represents the discretised continuity equation as an equa-
tion for pressure correction p′. The source term b′ in the equation is
the continuity imbalance arising from the incorrect velocity field u*, v*. By
solving equation (6.32), the pressure correction field p′ can be obtained at
all points. Once the pressure correction field is known, the correct pressure
field may be obtained using formula (6.14) and velocity components through
correction formulae (6.24)–(6.27). The omission of terms such as ∑a
nb
u ′
nb
in
the derivation does not affect the final solution because the pressure correc-
tion and velocity corrections will all be zero in a converged solution, giving
p* = p, u* = u and v* = v.
The pressure correction equation is susceptible to divergence unless some
under-relaxation is used during the iterative process, and new, improved,
pressures p
new
are obtained with
p
new
= p* +
α
p
p′ (6.33)
where
α
p
is the pressure under-relaxation factor. If we select
α
p
equal to 1
the guessed pressure field p* is corrected by p′. However, the corrections p′,
in particular when the guessed field p* is far away from the final solution, is
often too large for stable computations. A value of
α
p
equal to zero would
apply no correction at all, which is also undesirable. Taking
α
p
between 0 and
1 allows us to add to guessed field p* a fraction of the correction field p′ that
is large enough to move the iterative improvement process forward, but
small enough to ensure stable computations.
The velocities are also under-relaxed. The iteratively improved velocity
components u
new
and v
new
are obtained from
u
new
=
α
u
u + (1 −
α
u
)u
(n−1)
(6.34)
v
new
=
α
v
v + (1 −
α
v
)v
(n−1)
(6.35)
where
α
u
and
α
v
are the u- and v-velocity under-relaxation factors, u and v
are the corrected velocity components without relaxation, and u
(n−1)
and v
(n−1)
represent their values obtained in the previous iteration. After some algebra
it can be shown that with under-relaxation the discretised u-momentum
equation takes the form
u
i, J
=∑a
nb
u
nb
+ (p
I−1, J
− p
I, J
)A
i, J
+ b
i, J
+ (1 −
α
u
) u
i, J
(n −1)
(6.36)
and the discretised v-momentum equation
v
I, j
=∑a
nb
v
nb
+ (p
I, J−1
− p
I, J
)A
I, j
+ b
I, j
+ (1 −
α
v
) v
I, j
(n −1)
(6.37)
The pressure correction equation is also affected by velocity under-relaxation,
and it can be shown that d-terms of the pressure correction equation become
d
i, J
= , d
i+1, J
= , d
I, j
= , d
I, j+1
=
Note that in these formulae a
i, J
, a
i+1, J
, a
I, j
and a
I, j+1
are the central coeffici-
ents of discretised velocity equations at positions (i, J ), (i + 1, J ), (I, j) and
(I, j + 1) of a scalar cell centred around P.
A
I, j+1
α
v
a
I, j+1
A
I, j
α
v
a
I, j
A
i+1, J
α
u
a
i+1, J
A
i, J
α
u
a
i, J
J
K
L
a
I, j
α
v
G
H
I
a
I, j
α
v
J
K
L
a
i, J
α
u
G
H
I
a
i, J
α
u
ANIN_C06.qxd 29/12/2006 09:59 AM Page 189
A correct choice of under-relaxation factors
α
is essential for cost-
effective simulations. Too large a value of
α
may lead to oscillatory or even
divergent iterative solutions, and a value which is too small will cause
extremely slow convergence. Unfortunately, the optimum values of under-
relaxation factors are flow dependent and must be sought on a case-by-case
basis. The use of under-relaxation will be discussed further in Chapters 7
and 8.
The SIMPLE algorithm gives a method of calculating pressure and
velocities. The method is iterative, and when other scalars are coupled to the
momentum equations the calculation needs to be done sequentially. The
sequence of operations in a CFD procedure which employs the SIMPLE
algorithm is given in Figure 6.6.
190 CHAPTER 6 ALGORITHMS FOR PRESSURE---VELOCITY COUPLING
Figure 6.6 The SIMPLE
algorithm
Assembly of a
complete method
6.5
ANIN_C06.qxd 29/12/2006 09:59 AM Page 190
6.6 THE SIMPLER ALGORITHM 191
The SIMPLER (SIMPLE Revised) algorithm of Patankar (1980) is an
improved version of SIMPLE. In this algorithm the discretised continuity
equation (6.29) is used to derive a discretised equation for pressure,
instead of a pressure correction equation as in SIMPLE. Thus the intermedi-
ate pressure field is obtained directly without the use of a correction. Velocities
are, however, still obtained through the velocity corrections (6.24)–(6.27) of
SIMPLE.
The discretised momentum equations (6.12)–(6.13) are rearranged as
u
i, J
=+(p
I−1, J
− p
I, J
) (6.38)
v
I, j
=+(p
I, J −1
− p
I, J
) (6.39)
In the SIMPLER algorithm pseudo-velocities û and W are now defined as
follows:
û
i, J
= (6.40)
W
I, j
= (6.41)
Equations (6.38) and (6.39) can now be written as
u
i, J
= û
i, J
+ d
i, J
(p
I−1, J
− p
I, J
) (6.42)
v
I, j
= W
I, j
+ d
I, j
(p
I, J−1
− p
I, J
) (6.43)
The definition for d, introduced in the developments of section 6.4, is
applied in (6.42)–(6.43). Substituting for u
i, J
and v
I, j
from these equations
into the discretised continuity equation (6.29), using similar forms for u
i+1, J
and v
I, j+1
, results in
[
ρ
i+1, J
A
i+1, J
(û
i+1, J
+ d
i+1, J
(p
I, J
− p
I+1, J
)) −
ρ
i, J
A
i, J
(û
i, J
+ d
i, J
(p
I−1, J
− p
I, J
))] + [
ρ
I, j+1
A
I, j+1
(W
I, j+1
+ d
I, j+1
(p
I, J
− p
I, J+1
))
−
ρ
I, j
A
I, j
(W
I, j
+ d
I, j
(p
I, J−1
− p
I, J
))] = 0 (6.44)
Equation (6.44) may be rearranged to give a discretised pressure equation
a
I, J
p
I, J
= a
I+1, J
p
I+1, J
+ a
I−1, J
p
I−1, J
+ a
I, J+1
p
I, J+1
+ a
I, J−1
p
I, J−1
+ b
I, J
(6.45)
where a
I, J
= a
I+1, J
+ a
I−1, J
+ a
I, J+1
+ a
I, J−1
and the coefficients are given below:
a
I+1, J
a
I−1, J
a
I, J+1
a
I, J−1
b
I, J
(
ρ
dA)
i+1, J
(
ρ
dA)
i, J
(
ρ
dA)
I, j+1
(
ρ
dA)
I, j
(
ρ
ûA)
i, J
− (
ρ
ûA)
i+1, J
+ (
ρ
WA)
I, j
− (
ρ
WA)
I, j+1
Note that the coefficients of equation (6.45) are the same as those in the
discretised pressure correction equation (6.32), with the difference that
the source term b is evaluated using the pseudo-velocities. Subsequently, the
discretised momentum equations (6.12)–(6.13) are solved using the pressure
∑a
nb
v
nb
+ b
I, j
a
I, j
∑a
nb
u
nb
+ b
i, J
a
i, J
A
I, j
a
I, j
∑a
nb
v
nb
+ b
I, j
a
I, j
A
i, J
a
i, J
∑a
nb
u
nb
+ b
i, J
a
i, J
The SIMPLER
algorithm
6.6
ANIN_C06.qxd 29/12/2006 09:59 AM Page 191
192 CHAPTER 6 ALGORITHMS FOR PRESSURE---VELOCITY COUPLING
Figure 6.7 The SIMPLER
algorithm
field obtained above. This yields the velocity components u* and v*. The
velocity correction equations (6.24)–(6.27) are used in the SIMPLER algo-
rithm to obtain corrected velocities. Therefore, the p′-equation (6.32) must
also be solved to obtain the pressure corrections needed for the velocity
corrections. The full sequence of operations is described in Figure 6.7.
ANIN_C06.qxd 29/12/2006 09:59 AM Page 192
6.8 THE PISO ALGORITHM 193
The SIMPLEC (SIMPLE-Consistent) algorithm of Van Doormal and Raithby
(1984) follows the same steps as the SIMPLE algorithm, with the difference
that the momentum equations are manipulated so that the SIMPLEC velocity
correction equations omit terms that are less significant than those in SIMPLE.
The u-velocity correction equation of SIMPLEC is given by
u ′
i, J
= d
i, J
(p ′
I−1, J
− p ′
I, J
) (6.46)
where d
i, J
= (6.47)
Similarly the modified v-velocity correction equation is
v ′
I, j
= d
I, j
(p ′
I, J−1
− p ′
I, J
) (6.48)
where d
I, j
= (6.49)
The discretised pressure correction equation is the same as in SIMPLE, except
that the d-terms are calculated from equations (6.47) and (6.49). The sequence
of operations of SIMPLEC is identical to that of SIMPLE (see section 6.5).
The PISO algorithm, which stands for Pressure Implicit with Splitting of
Operators, of Issa (1986) is a pressure–velocity calculation procedure developed
originally for non-iterative computation of unsteady compressible flows. It has
been adapted successfully for the iterative solution of steady state problems.
PISO involves one predictor step and two corrector steps and may be seen as
an extension of SIMPLE, with a further corrector step to enhance it.
Predictor step
Discretised momentum equations (6.12)–(6.13) are solved with a guessed or
intermediate pressure field p* to give velocity components u* and v* using
the same method as the SIMPLE algorithm.
Corrector step 1
The u* and v* fields will not satisfy continuity unless the pressure field p* is
correct. The first corrector step of SIMPLE is introduced to give a velocity
field (u**, v**) which satisfies the discretised continuity equation. The result-
ing equations are the same as the velocity correction equations (6.21)–(6.22)
of SIMPLE but, since there is a further correction step in the PISO algo-
rithm, we use a slightly different notation:
p** = p* + p′
u** = u* + u′
v** = v* + v′
These formulae are used to define corrected velocities u** and v**:
u**
i, J
= u*
i, J
+ d
i, J
(p′
I−1, J
− p′
I, J
) (6.50)
v**
I, j
= v*
I, j
+ d
I, j
(p′
I, J−1
− p′
I, J
) (6.51)
As in the SIMPLE algorithm equations (6.50)–(6.51) are substituted into
the discretised continuity equation (6.29) to yield pressure correction
A
I, j
a
I, j
−∑a
nb
A
i, J
a
i, J
−∑a
nb
The SIMPLEC
algorithm
6.7
The PISO
algorithm
6.8
ANIN_C06.qxd 29/12/2006 09:59 AM Page 193
194 CHAPTER 6 ALGORITHMS FOR PRESSURE---VELOCITY COUPLING
equation (6.32) with its coefficients and source term. In the context of the
PISO method equation (6.32) is called the first pressure correction equation.
It is solved to yield the first pressure correction field p′. Once the pressure
corrections are known, the velocity components u** and v** can be obtained
through equations (6.50)–(6.51).
Corrector step 2
To enhance the SIMPLE procedure PISO performs a second corrector step.
The discretised momentum equations for u** and v** are
a
i, J
u**
i, J
=∑a
nb
u*
nb
+ (p**
I−1, J
− p**
I, J
)A
i, J
+ b
i, J
(6.12)
a
I, j
v **
I, j
=∑a
nb
v*
nb
+ (p**
I, J−1
− p**
I, J
)A
I, j
+ b
I, j
(6.13)
A twice-corrected velocity field (u***, v***) may be obtained by solving the
momentum equations once more:
a
i, J
u***
i, J
=∑a
nb
u**
nb
+ (p***
I−1, J
− p***
I, J
)A
i, J
+ b
i, J
(6.52)
a
I, j
v***
I, j
=∑a
nb
v**
nb
+ (p***
I, J−1
− p***
I, J
)A
I, j
+ b
I, j
(6.53)
Note that the summation terms are evaluated using the velocities u** and
v** calculated in the previous step.
Subtraction of equation (6.12) from (6.52) and (6.13) from (6.53) gives
u
i, J
*** = u**
i , J
++d
i, J
(p ″
I−1, J
− p ″
I, J
) (6.54)
v
I, J
*** = v**
I, j
++d
I, j
(p ″
I, J−1
− p″
I, J
) (6.55)
where p″ is the second pressure correction so that p*** may be obtained by
p*** = p** + p″ (6.56)
Substitution of u*** and v*** in the discretised continuity equation (6.29)
yields a second pressure correction equation
a
I, J
p ″
I, J
= a
I+1, J
p ″
I+1, J
+ a
I−1, J
p ″
I−1, J
+ a
I, J+1
p ″
I, J+1
+ a
I, J−1
p ″
I, J−1
+ b ″
I, J
(6.57)
with a
I, J
= a
I+1, J
+ a
I−1, J
+ a
I, J+1
+ a
I, J−1
, and the neighbour coefficients are
as follows:
a
I+1, J
a
I−1, J
a
I, J+1
a
I, J−1
b ″
I, J
(
ρ
dA)
i+1, J
(
ρ
dA)
i, J
(
ρ
dA)
I, j+1
(
ρ
dA)
I, j
i, J
∑a
nb
(u
nb
** − u*
nb
) −
i+1, J
∑a
nb
(u
nb
** − u*
nb
)
+
I, j
∑a
nb
(v
nb
** − v *
nb
) −
I, j+1
∑a
nb
(v
nb
** − v*
nb
)
In the derivation of (6.57) the source term
[(
ρ
Au**)
i, J
− (
ρ
Au**)
i+1, J
+ (
ρ
Av**)
I, j
− (
ρ
Av**)
I, j+1
]
is zero since the velocity components u** and v** satisfy continuity.
J
K
L
D
E
F
ρ
A
a
A
B
C
D
E
F
ρ
A
a
A
B
C
D
E
F
ρ
A
a
A
B
C
D
E
F
ρ
A
a
A
B
C
G
H
I
∑a
nb
(v
nb
** − v*
nb
)
a
I, j
∑a
nb
(u
nb
** − u*
nb
)
a
i, J
ANIN_C06.qxd 29/12/2006 09:59 AM Page 194
6.8 THE PISO ALGORITHM 195
Figure 6.8 The PISO algorithm
Equation (6.57) is solved to obtain the second pressure correction field p″,
and the twice-corrected pressure field is obtained from
p*** = p** + p″=p* + p′+p″ (6.58)
Finally the twice-corrected velocity field is obtained from equations
(6.54)–(6.55).
In the non-iterative calculation of unsteady flows the pressure field p***
and the velocity field u*** and v*** are considered to be the correct u, v and
p. The sequence of operations for an iterative steady state PISO calculation
is given in Figure 6.8.
ANIN_C06.qxd 29/12/2006 09:59 AM Page 195
The PISO algorithm solves the pressure correction equation twice so the
method requires additional storage for calculating the source term of the sec-
ond pressure correction equation. As before, under-relaxation is required
with the above procedure to stabilise the calculation process. Although
this method implies a considerable increase in computational effort it has
been found to be efficient and fast. For example, for a benchmark laminar
backward-facing step problem Issa et al. (1986) reported a reduction of CPU
time by a factor of 2 compared with standard SIMPLE.
The PISO algorithm presented above is the adapted, steady state version
of an algorithm that was originally developed for non-iterative time-
dependent calculations. The transient algorithm can also be applied to steady
state calculations by starting with guessed initial conditions and solving as a
transient problem for a long period of time until the steady state is achieved.
This will be discussed in Chapter 8.
The SIMPLE algorithm is relatively straightforward and has been success-
fully implemented in numerous CFD procedures. The other variations of
SIMPLE can produce savings in computational effort due to improved con-
vergence. In SIMPLE, the pressure correction p′ is satisfactory for correct-
ing velocities but not so good for correcting pressure. Hence the improved
procedure SIMPLER uses the pressure corrections to obtain velocity cor-
rections only. A separate, more effective, pressure equation is solved to yield
the correct pressure field. Since no terms are omitted to derive the discre-
tised pressure equation in SIMPLER, the resulting pressure field corre-
sponds to the velocity field. Therefore, in SIMPLER the application of the
correct velocity field results in the correct pressure field, whereas it does not
in SIMPLE. Consequently, the method is highly effective in calculating the
pressure field correctly. This has significant advantages when solving the
momentum equations. Although the number of calculations involved in
SIMPLER is about 30% larger than that for SIMPLE, the fast convergence
rate reportedly reduces the computer time by 30–50% (Anderson et al., 1984).
Further details of SIMPLE and its variants may be found in Patankar (1980).
SIMPLEC and PISO have proved to be as efficient as SIMPLER in
certain types of flows but it is not clear whether it can be categorically stated
that they are better than SIMPLER. Comparisons have shown that the
performance of each algorithm depends on the flow conditions, the degree
of coupling between the momentum equation and scalar equations (in
combusting flows, for example, due to the dependence of the local density on
concentration and temperature), the amount of under-relaxation used, and
sometimes even on the details of the numerical technique used for solving
the algebraic equations. A comprehensive comparison of PISO, SIMPLER
and SIMPLEC methods for a variety of steady flow problems by Jang et al.
(1986) showed that, for problems in which momentum equations are not
coupled to a scalar variable, PISO showed robust convergence behaviour and
required less computational effort than SIMPLER and SIMPLEC. It was
also observed that when the scalar variables were closely linked to velocities,
PISO had no significant advantage over the other methods. Iterative methods
using SIMPLER and SIMPLEC have robust convergence characteristics
in strongly coupled problems, and it could not be ascertained which of
SIMPLER or SIMPLEC was superior.
196 CHAPTER 6 ALGORITHMS FOR PRESSURE---VELOCITY COUPLING
General
comments on
SIMPLE, SIMPLER,
SIMPLEC and PISO
6.9
ANIN_C06.qxd 29/12/2006 09:59 AM Page 196