Issa 1986
Issa 1986
Issa 1986
AND
INTRODUCTION
In a companion paper [l], a method (PISO) is presented for the solution of the
implicitly discretised fluid flow equations by splitting of operations. The method
dispenses with outer’ iteration and is equally applicable to both compressible and
incompressible flows. In the aforementioned paper, the splitting technique is presen-
ted and is rhen assessed for accuracy and stability in relation to a linearised form of
’ The term “outer iteration” is used to distinguish the iteration process on the coupled sets of
equations for different variables from the iteration process which may be used to solve the simultaneous
nodal tinite-difference equations for a single variable.
66
no2 l-999 I!86 $3.00
APPLICATIONS OF PISO &54
the equations. It was found that for such a system, the temporal error introduced
by the splitting scheme (herein called splitting error) is of higher order in 6r than
the errors embodied in most of the temporal difference schemes currently used in
time-discretising the original differential equations (herein called temporal dis-
cretisation error). This not only would enable the control of the splitting error by
the time-step size 6t (rather than by recourse to iteration), but should also lead, a-t
least in the case of low-order temporal difference schemes, to the attainment of
time-accurate solutions at 6t values comparable to those dictated by the accuracy of
the difference scheme. Also, depending on the spatial difference scheme used, the
method can be shown to retain some of the stability endowed by implicit differen-
cing.
In principle, therefore, PISO should be efficient for time-dependent calculations
as iteration is disposed of without paying the penalty of having to reduce Or in
order to reduce splitting errors. At the same time, because of the ability to cope
with large 6t, the method should also be useful for applications to steady-stare
problems. Whether these merits are preserved when it is applied to the full non-
linear equations which govern fluid flow is the subject of the present paper.
First, it is useful to state the objectives which need to be fulfilled in order to
validate the procedure; these are as follows.
(i) To demonstrate that the method is capable of handling both com-
pressible and incompressible flows.
(iii To show, at least for one particular combination of spatial and temporal
discretisation schemes and a given spatial mesh size, that the temporal errors
introduced by the splitting procedure vanish with 6t and do so at least as fast as the
discretisation errors due to the temporal differencing scheme (first-order-accurate in
the present case). The implication here is that PISO should not hinder the
achievement of a time-accurate solution at the 6r values dictated by the accuracy of
the difference scheme.
(iii) To verify the saving in computing effort achieved by dispensing with
iteration when computing time-accurate solutions.
(iv) To demonstrate the ability to handle large time-steps which renders the
method also useful for steady-state calculations.
To accomplish these objectives, comparisons are made against computations with
an existing iterative method employing the same spatial and temporal difference
practices. This method yields the exact solution (to within the iteration tolerance
imposed) to the discretised equations over one time-step, i.e., the errors in the
calculated fields are solely due to the spatial and temporal discretisation practices
employed. By refinement of bt, the value at which the particular temporal difference
scheme achieves a time-accurate solution can thus be determined. Since PISO is
intended as an alternative to iteration, a comparison between the computing efforts
required by the two will serve as a measure of relative performance.
Other non-iterative, time-marching schemesexist which are either semi-implicit,
68 ISSA, GOSMAN, AND WATKINS
such as those in [2, 31, or fully implicit and based on time-splitting (e.g., ADI) of
the spatial fluxes, such as those of [7-91. The first group of methods is subject to
well-known restrictions on 6t imposed by stability considerations and is hence at a
disadvantage. The second group has to be based on the compressible form of the
continuity equation, which restricts their application to that particular class of
flows. Furthermore, the latter methods in general require the solution of block
simultaneous sets of equations, which can be a far more complex affair than the
solution of scalar ones as is the case with PISO.
The best established methods using pressure and velocity as working variables
are those found in [5, 61, namely, the SIMPL,E and SIMPLER schemes, both of
which are iterative. Both were developed originally for steady-state incompressible
flow although extensions to time-dependent, compressible flow have been made, for
example, in [4]. Of the two, SIMPLER is reported to be the most efficient and
robust. Other variants on these schemes exist, in particular the PUP method in
[12], which shares some features with the incompressible version of PISO (see [l]
for a discussion of the similarities and differences). A more recent method which
also uses the same working variables, but employs block iteration, is reported in
[lo] to be efficient and robust for steady-state incompressible flow.
In what follows, the incorporation of the PISO methodology into a finite-volume
procedure is presented. The scheme employs a staggered grid for the storage of
velocity and pressure, and uses backwards differencing for the representation of the
temporal variations (i.e., the Euler implicit scheme). The method is applied to both
compressible and incompressible flows to demonstrate its versatility. A comparison
is made against computations performed with fully iterative methods based on the
SIMPLE algorithm [h-6] and employing an identical difference scheme in order to
assessthe relative efficiency and stability of PISO in. transient flows. A study is also
made of the performance of PISO for a steady-state incompressible flow
calculation.
The incompressible case considered is that of a laminar flow (Re = 100) through
a suddenly expanding circular pipe with an open outlet. The compressible case is
that of laminar flow in a similar pipe but with a closed end, where a peak Re of
1000 and a peak Mach number of 0.2 are attained.
(1)
APPLICATIONS OF PISO
where p is the laminar viscosity and Pr is the Prandtl number. The total energy e is
related to the temperature by
e = c,. T+ fz4$4, (4)
where C,. is the constant volume specific heat. The equation of state taken here is
that of a perfect gas
P= RpT (jj
where 4 now stands for any of the dependent variables U: t’, and r, r is the “dif-
fusion” coefficient of property 4, and S, contains ail the remaining terms present in
the parent equations, as well as terms arising from the coordinate transformation.
The continuity equation becomes
Discretisntion
A staggered grid arrangement in which velocity nodes are located in between
pressure ones, as depicted in Figs. 1 and 2, is used. The discretisation of the trans-
port equations is effected by a finite-volume technique for which purpose control
volumes (indicated in Figs. 1 and 2) surrounding each variable location are defined.
The formulation of the difference equations follows standard practices employed in
earlier work [4-61, hence only a brief description of it is provided below. The
spatial variations are approximated by hybrid upwind/central difference formulae
giving first- to second-order spatial accuracy. It should be stressed, however, that
the applicability of PISO is in no way restricted to this choice, as its application in
[14] to a nine-point skew upwind scheme verifies. The temporal difference scheme
used is the Euler implicit formula which, although of only first-order accuracy, is
chosen mainly because of the simplicity of its implementation. Here again, the
validity of PISO is in principle not restricted to this choice, as inspection of the
methodology in Cl] will reveal.
The derivation of the difference equations is now illustrated by considering
Eq. (6) as an example. Integration of this equation over a control volume such as
those shown in Fig. 2 gives
where the subscripts refer to the locations indicated in Fig. 2, and superscripts n
and y1+ 1 denote successivetime levels. As implied by the Euler implicit temporal
difference scheme used, all the spatial fluxes (which appear non-superscripted) in
Eq. (8) are evaluated at time n + 1. These spatial fluxes are now approximated by
APPLICATIONS OF PISO
the hybrid upwind/centred difference scheme mentioned earlier. Taking the flux F:,.
at the n:-face of the control volume as an example: it can be represented as
where M,, is the mass flux through the cell face w and is defined as
and a is the cell-face area. The quantity CIis a weighting factor which depends on
the cell Peclet number Pe as
1 2’
Lx,=- l+- for lPej d 2
2( Pe )
=1 for Pe>2
=o for Pe < --2 ill)
and Pe is given by
p116.X
Pe= -
( r i I,..
and !?, is the volume integral of the source term S,. The operator H’ relates to the
4 values prevailing at the surrounding nodes and is given by
H’(~)=A,~,+A,,~~,+A,~,+A,~, (17)
where p is given by
and the mass fluxes are evaluated at the n + 1 time level. The pressure equation can
now be derived from a combination of the discretised continuity relation (19) and
the momentum equations which are given by Eq. (15) when d stands for u and a;
this procedure is outlined below.
The velocities u,, u,, etc., in Eq. (19) are to be replaced by expressions obtained
from the momentum equations. Taking the velocity U, as an example, Eq. (15) for
this quantity can be cast into the form
u, = [HL - (PO- PE) a, + C,]/(B - A,) (21)
where the pressure gradient term is now written out explicitly, and the quantity C
contains all other terms in the parent equation. The mass flux at the e-face of the
cell can now be evaluated from Eq. (21) as
(pus), = A, - D,(Po - PE) + c, (22 1
where the following have been introduced:
(23)
(24)
(25)
APPLICATIONS OF PISO 73
Similar expressions for the mass fluxes at faces W,n, and s of the cell can be derived,
ail of which may be substituted into Eq. (19) to yield
D,P, = J(P) - div &- div c - fl(p’+’ -p”) i26)
and
D,,=D,fD,,+D,,+D,. (28)
Equation (26) corresponds to the pressure equation (12) in [I], and its form is
similar to the general 4 equation (15).
METHOD OF SOLUTION
The system of Eqs. (15) for u. c’, and e (the latter of which becomes redundant in
the case of incompressible flow) together with Eq. (26) for the pressure are solved
by the PISO algorithm, a full description of which appears in [l]. The non-
linearity in Eq. (15) arising from the dependency of the coefficients .4 and B on the
field variables themselves is handled by evaluating these coefficients from the old
time level values. Although this practice is only first-order accurate in time it is of
the same order of accuracy as the temporal difference scheme, and is therefore con-
sistent with it. Furthermore, for the same reason, all the computations for the com-
pressible flow case presented herein were made with the two-stage version of the
PISO algorithm presented in [l]. Indeed, computations with the three-stage
variant that were performed did not bring any improvement in accuracy while
requiring substantially higher computing effort.
The set of nodal simultaneous algebraic equations for each variable are solved by
a line successiveover-relaxation procedure to a specified convergence level. Recent
work [15] shows that considerable saving can be made if Stone’s strongly impkit
procedure [I131 is used instead, especially for the pressure equation. Further saving
can be made in the case of steady-state calculations, when temporal accuracy is of
no consequence, by relaxing the convergence tolerance on these sets of equations.
However, the pressure equations must still be converged to at least 10% of the
initial residuals, otherwise the method is forced to take more time-steps to arrive to
the final solution, to the detriment of effkiency.
For the purpose of evaluating the performance of PISO, a comparison is made
against calculations carried out with an iterative method. The method chosen is
that based on the well-established SIMPLE algorithm in [S, 61: which is com-
prehensively documented in the literature, The method was initially developed for
steady-state problems, but has been extended to compressible transient fiows (as
74 ISSA, GOSMAN, AND WATKINS
in [S]). For the steady-state case, the temporal derivative terms in the transport
equations are suppressed and iteration replaces time-marching with under-
relaxation used to procure stability and convergence. This under-relaxation is
related to the time step St of an equivalent transient calculation as shown in [12]
and in the Appendix. For unsteady-state problems, iteration is employed at each
time step to satisfy all the equations simultaneously. Under-relaxation must be used
here also in order to ensure the convergence of the iteration process. As is the case
with the PISO computations, a line successive over-relaxation scheme was
employed for each set of algebraic simultaneous equations.
RESULTS
Incompressible Flow
The incompressible flow case chosen is that of an axisymmetric laminar flow in a
duct with a sudden enlargement at entry (Fig. 3). The ratio of the diameter at inlet
to that of the duct is I:2 and the length of the duct is 4 times its diameter. The
time-varying velocity profile at inlet is assumed to be spatially uniform across inlet
and the flow is started from rest with an inlet velocity rising linearly from zero at
t = 0 to its final steady value of I/ at t = L/V: where L is the length of the duct. The
Reynolds number (based on the peak velocity V and the duct diameter) is 100. At
the outlet, zero velocity gradients are assumed, while at the walls? the usual no-slip
conditions are imposed. The mesh size is 20 x 20 (uniform in both directions)
throughout and is kept unaltered. The tests carried out are designed to serve three
purposes. The first is to demonstrate that the accuracy of the PISO method is at
least as good as the accuracy of the temporal difference scheme used. The second is
to illustrate the stability of the schemefor large at, which makes it equally suited to
the calculation of steady-state problems. And third, to show that the first objective
is achieved at substantially lower computing effort than with existing methods
utilising sequential iteration, as the one compared here.
10
08
"
ilstC&
0.6
0.4
02
1 I I I I i
0
-0L
-0 8
-1.2
I
-1. 6 I I I I I
values of &I&,, as well as that for the time-accurate solution obtained with the
iterative method. The latter curve is indistinguishable from that pertaining to
6t/6to = 1, verifying that the PISO solution is at least as accurate as the temporal
difference scheme used, and that this algorithm allows the use of the maximum
valued of bt permitted by the accuracy of that difference scheme.
The benefits of avoidance of iteration were verified by the finding that the ratio of
the computing time required by the two methods (with 6t/6to = 1) was 9.3 in favour
of PISO.
Figure 4 also shows that PISO remains stable for btjbt, as large as 20. This is
very useful when only the steady-state solution is sought (in which case temporal
accuracy is of no consequence). Taking large time-steps permits rapid arrival at the
steady state with minimum computing effort, as exemplified in Fig. 6. This figure
shows the computing times (on a CDC 6400 machine) required respectively by
PISO and the steady-state version of SIMPLE (in which time derivatives are omit-
ted) to reach that state. In the iterative calculations the under-relaxation factor i
(imposed on the velocities only) was varied to determine the optimum value, i.e.,
that yielding the minimum computing time. A similar exercise was then performed
with PISO, using now dt as the controlling parameter; thus the abscissa in Fig. 6 is
either 2 or 6t (made dimensionless by dividing by the convection time scale 6x/V)
according to the method.
Two facts emerge from an examination of this plot. The first is that PISO arrives
at the solution at a much lower level of computing effort than the other method,
whatever the value of ,? or ht. The second is that the PISO curve is remarkably flat
over a wide range of values, thus indicating stability and robustness. In contrast,
00 0.5 1 1.0
CPU
I seconds1
100
the performance of SIMPLE is fairly sensitive to the relaxation parameter and the
scheme diverges beyond a value of 0.55 (which is equivalent to GrV/Ex of 1.2 for a
time-marching method).
It should. however, be pointed out that the SIMPLE scheme used for these com-
parisons is not the most efficient of the existing iterative methods. For example,
several versions of the same scheme have been proposed and tested in [12], some
of which show marked improvement in efftciency over the original one. Other
techniques, notably SIMPLER and PUP in [6] and [12], respectively. are also
reported to be much faster as well as more stable.
A qualitative assessment of the results gleaned from [12] shows that rhc
improvement in speed over SIMPLE by the best of these methods for steady-state
calculations is about the same as that achieved by PISO for steady-state problems
(including ones not presented here). This is not surprising as the latter shares some
common features with the others as explained in [ 1]. However, the improvement in
speed attained by these other methods over SIMPLE still does not match the
margin of saving achieved by PESO in time-evolving flows, thanks to its non-
iterative strategy.
ln this flow case, the geometry of the duct is similar to that of the previous case.
with the exception that the downstream end is now closed, so that the fluid mass
inside the duct varies with time. The velocity profile at inlet, which is again assumed
to be uniform in the radial direction, is taken to vary sinusoidally with time with a
peak velocity V corresponding to a Reynolds number of 1000 and a Mach number
of 0.2. The period of the cycle is taken as twice the time taken for a fluid particle
traveling at the peak velocity to traverse the length of the duct; the computations
are carried out over a full cycle of velocity variations starting from rest. The density
and temperature at inlet are assumed to be uniform and constant with time and the
duct wall temperatures are taken to be constant at the same value as that at inlet. .4
uniform mesh of 20 x 20 is used throughout.
As in the incompressible case, calculations were first performed with the time-
marching, compressible version of SIMPLE using various values of c‘it to find the
value Sr, for which temporal truncation errors are negligible. Here again under-
relaxation was necessary to ensure convergence and values of A = 0.5 were used for
the velocities and temperature (but not pressure). The exercise was then repeated
with PISO using step sizes that are multiples of at,. The results of those com-
putations are displayed in Figs. 7 to 10, which show the predicted transience of the
axial velocity (normalised by the peak inlet value) and of the pressure (normalised
by the pressure at inlet) at two locations in the duct. The first point (Figs. 7 and 8)
is located at one-third of the duct length downstream of inlet and at 80% of the
duct radius, this being the location where some of the highest reverse velocities
occur. The second point (Figs. 9 and 10) is located on the centreline halfway down
the length of the duct. In all these figures, the time coordinate appears normaiised
78 ISSA, GOSMAN, AND WATKINS
by the period of the cyclic inlet-velocity variations. These figures indicate that the
solutions given by PISO for different values of the time-step size (st/st,) converge
to an acceptable level of temporal accuracy at the same value of 6t as that required
by the iterative method, i.e., dt,. This verifies that the accuracy of PISO for com-
pressible flow is as good as the accuracy of the temporal difference scheme
employed. The most significant finding emerging from these calculations was.
however, that the PISO computations required only 0.19 of the computing effort
demanded by the iterative scheme.
OL 06 08 IO
Normahsed Time
Normollsed Time
The figures also demonstrate that time-step sizes of up to 20 St, can be taken
with PISO without any signs of instability, whereas solutions could not be obtained
with the iterative method for values greater than 3 tit,: as the iteration process then
failed to converge with the chosen set of values of the under-relaxation parameter A.
Here again, the stability of PISO for large 6t renders it useful for the computation
of steady-state compressible flows.
CONCLUSIONS
---- 2
- 20
- lteratlve lbtzbt,l
0.2 OL 06 08 10
Normallsed TI me
551’6? l-6 FIG. 10. Predicted pressure transience on centreline jcompressible case j.
80 ISSA,GOSMAN,AND WATKINS
regarding the accuracy and stability of the algorithm. This is achieved through a
comparison against calculations performed with an existing iterative method. The
most important fact emerging from these comparisons is that PISO is many times
faster than its iterative counterpart for transient flows, whether compressible or
incompressible. Furthermore, it exhibits stable behaviour for large time-step sizes
which makes it a reliable technique for steady-state calculations.
Recent work Cl.51 where the PJSO method was applied to turbulent flows (for
which the--k-s model of turbulence was used) shows that the same significant
saving over iterative schemescan be achieved for such flows also. This is the case
provided that the splitting procedure for the source terms in the k and F equations
outlined in [l] is implemented in conjunction with PISO.
Consider a model implicitly discretised scalar equation (such as Eq. (15) in the
text) for time-dependent flow in the absence of sources and for a constant density
fluid. It can be written as
(B--A,)$h”,+‘=H’((LY+‘)+Bqg. (A.1)
The corresponding equation in steady-state form in which the time derivative terms
are suppressed would be obtained if the term B in Eq. (A.1) is set to zero and the
superscripts now stand as iteration counters. An iterative procedure based on this
form of the equation invariably requires under-relaxation in order to stabilise the
computations. A standard practice adopted in introducing under-relaxation is to
base the new iteration level value (i.e., at rz+ 1) on
4 g+‘=nq5,+(1-l)qg lA.2)
where d is the relaxation factor and the superscript denotes the value that would be
obtained if no relaxation was used. Thus the equation actually solved in an
iterative, under-relaxed procedure becomes
--:4 ;;+‘=H’(&“+‘)-~i1&
1-I (A.3)
Equation (A.3) is obtained from the substitution of relation (A.2) into the form of
E,q.(A.1) in which the term B is set to zero. It is not difficult to discern the
similarities between the time derivative terms in Eq. (A.l) and the damping terms
due to under-relaxation in Eq. (A.3). Equating these terms gives
APPLICATIONS OF PISO 81
Consider, for the moment, the case of convection-dominated flow (which is often
the case); the coefficient ‘4, is given approximately by
- .4, = purEr + pw6.u. +..6)
($+$)=A. (A.7)
It is evident from (A.7) that the relaxation parameter h is related to the local
Courant numbers, uSt/Gx, dt/6r, of the equivalent time-dependent procedure (a
similar result was also obtained in [ 121). For example, if o is zero, the case of E.= 0
corresponds to 6t = 0. while the case of A = 1 corresp0nd.sto the case of 6; = (=:‘.The
case of A = 0.5 corresponds to a local Courant number ufit,iS.u of 1.
Similar findings can be reached for diffusion-dominated flows, where now the
relaxation factor i becomes related to the parameters 12t/(pkx2) and F6rjipbr’). In
all cases,therefore, there is a definite relationship between the value of the time-step
size in a time-marching procedure and the under-relaxation factor employed rn
steady-state algorithms, a relationship which relates local flow properties with the
mesh size and 6r.
ACKNOWLEDGMENT
This work was carried out under the financial support of the Science and Engineering Research
Council.
REFERENCES
10. J. P. VAN DOORMAALAND G. D. RAITHB~,in “Proceedings 10th 04.4C.S World Congress on System
Simulation and Sci. Comp., 1981,” p. 128.
11. B. E. LAL~NDER AND D. B. SPALDING “Mathematical Models of Turbulence,” Academic Press, Lon-
don/New York, 1972.
12. G. D. RAITHBYAND G. E. SCHNEIDER, Numer. Heat Transfer 2 (1979), 417.
13. H. L. STONE,SIAM J. Numer. Anal. 5 (1968). 530.
14. R. W. BENODEKAR, A. J. H. GODDARD. A. D. GOSMAN,AND R. I. Issa, in “Proceedings, A.SIZfE
Fluids Engineering Congress, 1983.”
15. A. D. GOSMAN,Y. Y. Tsur, AND A. P. WATKINS, SAE paper 840229, 1984.