CHAPTER 7

PISO ALGORITHM FOR UNSTEADY FLOWS

7.1 Introduction
Several authors have presented works based on predictor-corrector operator-splitting
for unsteady flows. Issa first proposed the PISO (pressure implicit splitting of operators)
algorithm (Issa 1986) followed by its variants and its application to various flows (Issa et al.
1986, Issa and Javareshkian 1996). Though the PISO algorithm was formulated for both
incompressible and compressible flows, there is little evidence of its application to high-
speed compressible flows, especially supersonic flows. It is the objective of the present paper
to develop a pressure-based algorithm for the efficient computation of unsteady flows at all
Mach numbers (ranging from incompressible to hypersonic) in complex three-dimensional
geometries.
Here, two algorithms based on operator splitting for unsteady all-speed flows are
developed and assessed. First, the PISO algorithm of Issa (1986) is re-formulated in the
context of a three-dimensional flow solver employing multi-block, body-fitted (curvilinear),
structured grids. It is found that although the algorithm works well for unsteady simulation of
incompressible and subsonic flows, it is not stable for unsteady supersonic flows. To rectify
this, a variant of the original PISO algorithm is developed in the present work. The predictor-
corrector sequence of the proposed algorithm is similar to that of the original PISO scheme.
However, the density-velocity-pressure coupling follows the compressible extension of the
original SIMPLE algorithm for steady flows (e.g., Karki and Patankar 1989, Thakur and
Shyy).

w e
W P E
i-1/2 i+1/2
i-1 i i+1

Figure 7.1. Schematic of the control volume with contravariant velocities at interfaces.

STREAM Version 4.5.2 129

 7.2 PISO Algorithm
In this section, the discretized equations for the algorithm will be derived for a
compressible flow using curvilinear coordinates. For clarity and brevity, the derivation is
presented for a two-dimensional situation. However, there is no loss of generality; the
extension to three-dimensional flows is straightforward.
First consider the discretized equation governing u-momentum using the fully
implicit backward-Euler time-stepping scheme. In curvilinear coordinates, this equation can
be written as
 ρ n +1 J  n +1 ρ nu n J
A
 P +
∆t 
 u = H n
( u n +1
) +
∆t
− f11 ( pe − pw ) − f 21 ( pn − ps )
n n
(7.1)


where AP is the nodal coefficient involving only the spatial discretization. The term H n
contains all terms resulting from the discretization of the convective and diffusive fluxes as
well as any source terms with the exception of the pressure gradient terms. Note that the
superscript n denotes the old (or the present) time level and n+1 denotes the new time level
being computed. The implicit solution of the above equation requires several iterations at
every time step, which is very inefficient. The unsteady algorithm developed in the present
work seeks to eliminate this iterative solution procedure by employing a series of predictor-
corrector steps. The solution field at the predictor stage is denoted by the superscript *, that
in the successive corrector steps by **, ***, etc. The solution obtained at the last corrector
step is considered as the solution at the new time level (n+1).
7.2.1 Momentum Predictor Step

In the predictor step, the u- and v-momentum equations are solved implicitly and can
be written as follows:
 ρnJ  * ρ nu n J
A
 P +
∆t 
 u = H n
( u *
) +
∆t
− f11 ( pe − pw ) − f 21 ( pn − ps )
n n
(7.2)


 ρnJ  * ρ nvn J
 v = H (v ) + − f12 ( pe − pw ) − f 22 ( pn − ps )
n n
 AP +
n *
(7.3)
 ∆t  ∆t

Note that the above equations use the pressure field at the time level t n . Similarly the scalar
equations can be written.
7.2.2 First Momentum Corrector Step

In the first corrector stage of the overall splitting process, a new velocity field, ( u **
and v** ), is sought, along with a corresponding pressure field, p* . Towards this, the velocity
field obtained from the predictor step ( u * and v* ) is used to write the momentum equations
for the first corrector step as follows:
 ρ n J  ** ρ nu n J
A
 P +
∆t 
 u = H n
( u *
) +
∆t
− f11 ( pe − pw ) − f 21 ( pn − ps )
* *
(7.4)


STREAM Version 4.5.2 130

  ρnJ  ** ρ nvn J
A
 P +
∆t
 v = H n
( v *
) +
∆t
− f12 ( pe − pw ) − f 22 ( pn − ps )
* *
(7.5)
 

Note that we have used ρ n on the left hand side of the above equations
In Eqs. (7.4)-(7.5), p* is yet to be determined. To determine it, the above momentum
equations and the continuity equation are used to formulate an equation governing pressure
correction, p′ , which is defined as:
p′ = p * − p n (7.6)
First, we rewrite the momentum equations in terms of p′ . Then, subtracting Eqs. (7.2) and
(7.3) from Eqs. (7.4) and (7.5), respectively, we obtain:
1
u ** = u * − n  f11 ( pe′ − pw′ ) − f 21 ( pn′ − ps′ )  (7.7)
AP
1
v** = v* −  f12 ( pe′ − pw′ ) − f 22 ( pn′ − ps′ )  (7.8)
APn 
where we have defined
ρnJ
APn = AP + (7.9)
∆t
Eqs. (7.7) and (7.8) will be used to correct u * and v* to yield u ** and v** once the p′ field
is solved for. To obtain an equation for p′ , we first need to formulate the interface
contravariant velocities, which in turn require interface Cartesian velocity components. This
will be described next.
7.2.2.1 Formulation of Cell-Face Contravariant Velocities

We will illustrate the formulation of the contravariant velocities along the i-direction,
i.e., U * . Due to the finite volume approach being employed here, the contravariant
velocities, multiplied by the density, represent the mass flux at the faces of control volumes.
Thus, for example, for a control volume located arou9nd the point P (index i) shown in
Figure 1, the contravariant velocities are required at the east (i+1/2) and west (i-1/2) faces.
Thus, in the algorithm described here, the contravariant velocities are only stored at the
control volume interfaces, as shown in Figure 7.1.
In this section, we will describe the estimation of the contravariant velocity on the
east face, i.e., U i +1/ 2 . First, consider the Cartesian velocity components; rewrite Eqs. (7.2)
and (7.3) as:
G u* f
ui* = i n − 11n ( pe − p w )
n
(7.10)
AP AP

Giv* f12
( )
n
vi* = − pe − p w (7.11)
APn APn

STREAM Version 4.5.2 131

where
ρ nu n J
Giu = H ( u * ) + − f 21 ( pn − ps )
n
(7.12)
∆t
ρ nvn J
Giv = H ( v* ) + − f 22 ( pn − ps )
n
(7.13)
∆t
Similarly, for V * (i.e., the contravariant velocity along the j-direction), write Eqs. (7.2) and
(7.3) as:
G uj * f 21
n ( n
p − ps )
n
u *j = n
− (7.14)
A P AP

G vj * f 22
n ( n
p − ps )
n
v*j = n
− (7.15)
A P AP

Now, let us consider the east face (i+1/2) and write ui +1/ 2 and vi +1/ 2 similar to ui and vi in
Eqs. (7.10) and (7.11):
Giu+*1/ 2 f
( )
n
ui*+1/ 2 = n
− n11 pi +1 − p i (7.16)
APi+1/ 2 APi+1/ 2

Giv+*1/ 2 f
( )
n
v *
i +1/ 2 = n − n12 pi +1 − p i (7.17)
APi+1/ 2 APi+1/ 2

Note that Gi +1/ 2 in Eqs. (7.16) and (7.17) are evaluated using Eqs. (7.10)-(7.11) instead of Eqs.
(7.12)-(7.13), as follows:

Giu+*1/ 2 1  Giu* Giu+*1 

=  + 
APni+1/ 2 2  APni APni +1 

1 * * 1  f   f  
=
2
( ui + ui +1 ) +  11  ( pi +1/ 2 − pi −1/ 2 ) +  11  ( pi +3/ 2 − pi +1/ 2 ) 
2  AP i
(7.18)
 AP i +1 

Giv+*1/ 2 1  Giv* Giv+*1 

=  + 
APni+1/ 2 2  APni APni +1 

1 * * 1  f   f  
=
2
( vi + vi +1 ) +  12  ( pi +1/ 2 − pi −1/ 2 ) +  12  ( pi +3/ 2 − pi +1/ 2 ) 
2  AP i
(7.19)
 AP i +1 
Now we can formulate the east-interface contravariant velocity as:

STREAM Version 4.5.2 132

 U i*+1/ 2 = f11ui*+1/ 2 + f12 vi*+1/ 2 (7.20)

f11 * * f12 * * f i +21/ 2 fi +21/ 2

=
2
( i i+1 ) 2 ( i i+1 ) A
u + u + v + v − ( pi +1 − pi )
n
+
2 A
( pin+1/ 2 − pin−1/ 2 ) + ( pin+3/ 2 − pin+1/ 2 ) 
 
Pi +1 / 2 Pi+1 / 2

7.2.2.2 Formulation of Equation Governing p′

The next step is to formulate the equation for p′ . Towards this end, let us first
write the momentum equations for the i+1/2 face:
G u* f
( )
n
predictor: ui*+1/ 2 = in+1/ 2 − n11 pi +1 − p i (7.21)
APi+1/ 2 APi+1/ 2

Giu+*1/ 2 f
( )
*
corrector: ui**+1/ 2 = n
− n11 pi +1 − p i (7.22)
APi+1/ 2 APi+1/ 2

Now, subtracting Eq. (7.21) from Eq. (7.22), and doing the same for vi +1/ 2 , we get

f11
ui**+1/ 2 = ui*+1/ 2 − (
APni+1/ 2
pi′+1 − p′i ) (7.23)

f12
vi**+1/ 2 = vi*+1/ 2 − (
APni+1/ 2
pi′+1 − p′i ) (7.24)

Multiplying Eqs.(7.23) and (7.24) with f11 and f12 , respectively, and then subtracting the
former from the latter, we get the following:
(
U i**+1/ 2 = U i*+1/ 2 − Fi +1/ 2 pi′+1 − p′i )
(7.25)

where the coefficient F is defined as follows:

 f2+ f2
F =  11 n 12  (7.26)
 AP 
Similarly, at the i-1/2 interface, we can write:

(
U i**−1/ 2 = U i*−1/ 2 − Fi −1/ 2 pi′ − p′i−1 ) (7.27)

Now, we invoke the continuity equation for this stage, which is

( ρ *U ** ) J
− ( ρ *U ** )  = − ( ρ * − ρ n ) (7.28)
 i +1/ 2 i −1/ 2  ∆t
The density ρ * (which is an unknown so far) is expressed as

ρ* = ρ n + ρ′ (7.29)

STREAM Version 4.5.2 133

Using the above equations, the continuity equation, Eq. (7.28), can be written as follows:

(ρ n
+ ρ ′)
i +1/ 2
U i*+1/ 2 − ( ρ n + ρ ′ )
i −1/ 2
U i*−1/ 2 −
J (7.30)
( ρ + ρ ′)
n
i +1/ 2
Fi +1/ 2 ( pi′+1 − pi′ ) − ( ρ n + ρ ′ )
i −1/ 2
Fi −1/ 2 ( pi′ − pi′−1 ) = −
∆t
( ρ* − ρ n )

In Eq. (7.30), the terms involving the product of ρ ′ and p′ are neglected. Thus, the equation
governing p′ can be written as:

( ρ n F )i+1/ 2 ( pi′+1 − pi′ ) − ( ρ n F )i −1/ 2 ( pi′ − pi′−1 ) + ρi′+1/ 2U i*+1/ 2 − ρi′−1/ 2U i*−1/ 2 + ∆ρti J =
′
(7.31)
(ρ U ) −(ρ U )
n *
i +1/ 2
n *
i −1/ 2

where the term involving density difference is treated implicitly by invoking the equation of
state:
p′
ρ′ = (7.32)
RT
Moreover, when evaluating ρ ′U * at the control volume interface, an upwind scheme is
employed, so that at the east (i+1/2) interface, for example, we have
 pi′ *
 U i +1/ 2 , if U i*+1/ 2 > 0
p ′ 
ρi′+1/ 2U i*+1/ 2 = i +1/ 2 U i*+1/ 2 =  RT (7.33)
RT  p ′
i +1
U *
, if U i +1/ 2 < 0
*
 RT i +1/ 2

The final discretized form of the equation governing the p′ field can be written as:
AP p′P = AE p′E + AW pW′ + AN p′N + AS pS′ + b (7.34)
where
1
AE = ( ρ n F ) + max ( −U i +1/ 2 , 0 ) (7.35)
i +1/ 2 RT
1
AW = ( ρ n F ) + max (U i −1/ 2 , 0 ) (7.36)
i −1/ 2 RT
1 1 J
AP = ( ρ n F ) + (ρnF ) + max (U i +1/ 2 , 0 ) + max ( −U i −1/ 2 , 0 ) + (7.37)
i +1/ 2 i −1/ 2 RT RT RT ∆t
b = ( ρ nU * ) − ( ρ nU * ) (7.38)
i +1/ 2 i −1/ 2

and the coefficients AN and AS have a form similar to AE and AW along with the
corresponding contributions to the coefficient AP .

STREAM Version 4.5.2 134

 7.3 Second Momentum Corrector Step
The u-momentum equation for the second corrector stage is written as:
 ρnJ  *** ρ nu n J
A
 P +
∆t
 u = H *
( ) ∆t − f11 ( pe − pw )* − f 21 ( pn − ps )**
u **
+ (7.39)
 

Note that, similar to Eq. (7.4), we have used ρ n in the term on the left hand side of Eq. (7.39)
Let us define the pressure correction at the second stage as
p′′ = p** − p* (7.40)
Subtracting the first corrector equation, Eq. (7.4), from Eq. (7.39), we get

APn ( u *** − u ** ) =  H * ( u ** ) − H n ( u * )  − f11 ( pe′′ − pw′′ ) − f 21 ( pn′′ − ps′′) (7.41)

from which we obtain

u *** =
1 **
n
AP
{
u +  H * ( u ** ) − H n ( u * )  − f11 ( pe′′ − pw′′ ) − f 21 ( pn′′ − ps′′) } (7.42)

Next, we will derive the contravariant velocities (or fluxes) at the control volume
interfaces for the second corrector stage. First, let us rewrite Eqs. (7.39) and (7.4) as follows:
1
u *** = n uˆ * − f11 ( pe − pw ) 
**
2nd Corrector: (7.43)
AP  

1  n
uˆ − f11 ( pe − pw ) 
*
1st Corrector: u ** = n 
(7.44)
AP 

where
ρ nu n J
uˆ * = H * ( u** ) + − f 21 ( pn − ps )
*
(7.45)
∆t
ρ nu n J
uˆ n = H n ( u* ) + − f 21 ( pn − ps )
n
(7.46)
∆t
Now, let us rewrite Eqs. (7.43) and (7.44) for the interface i+1/2 as follows:
1  *
− f11 ( pi +1 − pi ) 
**
2nd Corrector: +1/ 2 =
ui*** uˆ
n  i +1/ 2
(7.47)
AP 

1  n
− f11 ( pi +1 − pi ) 
*
1st Corrector: ui**+1/ 2 = uˆ
n  i +1/ 2
(7.48)
AP 

Subtracting Eq. (7.48) from Eq. (7.47) and rearranging, we get

+1/ 2 = ui +1/ 2 + ( ui +1/ 2 − ui + /12 ) − f11 ( pi +1 − pi )

ui*** **
ˆ* ˆn ′′ ′′ (7.49)

STREAM Version 4.5.2 135

Similarly, we can write the expression for vi***
+1/ 2 as

+1/ 2 = vi +1/ 2 + ( vi +1/ 2 − vi + /12 ) − f12 ( pi +1 − pi )

vi*** **
ˆ* ˆn ′′ ′′ (7.50)

Then, multiplying Eqs. (7.49) and (7.50) with f11 and f12 , respectively, the contravariant
velocity for the second corrector stage at the i+1/2 interface can be written as:
+1/ 2 = U i +1/ 2 + U i +1/ 2 − Fi +1/ 2 ( pi +1 − pi )
U i*** ** ˆ ′′ ′′ (7.51)

where
 f   f 
Uˆ i +1/ 2 =  11n 
A
( uˆi*+1/ 2 − uˆin+1/ 2 ) +  12n 
A
( vˆi*+1/ 2 − vˆin+1/ 2 ) (7.52)
 P i +1/ 2  P i +1/ 2
and F is defined in Eq. (7.26). The contravariant velocity at the i-1/2 interface can be
obtained in a similar manner.
Next, we invoke the continuity equation for the second corrector stage:
J
= − ( ρ ** − ρ n )
i +12
 ρ **U ***  (7.53)
i −1/ 2 ∆t
The density for the second stage, ρ ** , is expressed as
ρ ** = ρ * + ρ ′′ (7.54)
Using the above equations, the continuity equation, Eq. (7.53), can be written as:

( ρ * + ρ ′′ ) U ** 
 i −1/ 2 + ( ρ + ρ ) U  i −1/ 2 −
i +1/ 2 i +1/ 2
 * ′′ ˆ 

(7.55)
( ρ + ρ ′′)i +1/ 2 Fi +1/ 2 ( pi′′+1 − pi′′) + ( ρ + ρ ′′)i−1/ 2 Fi −1/ 2 ( pi′′− pi′′−1 ) = − ∆Jt ( ρ ** − ρ n )
* *

In Eq. (7.55), the terms involving the product of ρ ′′ and p′′ are neglected. Thus, the equation
governing p′′ can be written as:

( ρ * F )i+1/ 2 ( pi′′+1 − pi′′) − ( ρ * F )i−1/ 2 ( pi′′− pi′′−1 ) + ρi′′+1/ 2U i**+1/ 2 − ρi′′−1/ 2U i**−1/ 2 + ( i ∆t i ) J =
ρ ′′+ ρ ′
(7.56)
( )
i +1/ 2
 ρ U + Uˆ 
* **
  i −1/ 2
where the term involving density difference is treated implicitly by invoking the equation of
state:
p′′
ρ ′′ = (7.57)
RT
Moreover, the ρ ′′U ** term at control volume interfaces is evaluated using an upwind
scheme, similar to Eq. (7.33).
The final discretized form of the equation governing the p′′ field can be written as:
AP p′′P = AE p′′E + AW pW′′ + AN p′′N + AS pS′′ + b (7.58)

STREAM Version 4.5.2 136

where
1
AE = ( ρ * F ) + max ( −U i*+1/ 2 , 0 ) (7.59)
i +1/ 2 RT
1
AW = ( ρ * F ) + max (U i*−1/ 2 , 0 ) (7.60)
i −1/ 2 RT
1 1 J
AP = ( ρ * F ) + ( ρ *F ) + max (U i*+1/ 2 , 0 ) + max ( −U i*−1/ 2 , 0 ) + (7.61)
i +1/ 2 i −1/ 2 RT RT RT ∆t
ρ′
b = ( ρ *U ** )
i +1/ 2
− ( ρ *U ** )
i −1/ 2
(
+ ρ *Uˆ ) i +1/ 2
(
− ρ *Uˆ ) i −1/ 2
+
∆t
J (7.62)

and the coefficients AN and AS have a form similar to AE and AW along with the
corresponding contributions to the coefficient AP .

7.4 Summary of the Algorithm

• The algorithm described in the preceding section can be summarized as
follows:
• Predictor Eq. (7.2) and (7.3)
• Compute contravariant velocities: Eq. (7.20)
• Compute p′ : Eq. (7.34)
• Update velocities using p′ : u ** and v** : Eq. (7.23)-(7.24); U ** and V ** : Eq.
(7.27)
• Update pressure and density: p* from Eq. (7.6) and ρ * from equation of state
• Compute Û and Vˆ : Eq. (7.52)
• Compute p′′ : Eq. (7.58)
• Update velocities using p′′ : u *** and v*** : Eq. (7.42); U *** and V *** : Eq.
(7.51)
• Update pressure and density: p** from Eq. (7.40) and ρ ** from equation of
state
u , v , U , V *** , p** , ρ ** are considered the values at the new time level (n+1) and we
*** *** ***

proceed to step 1 for the next time level.

7.5 Test Cases
7.5.1 Evolution of Steady Incompressible Flow Past a Cylinder

The laminar, incompressible flow past an impulsively started circular cylinder in a

quiescent fluid is used as the first test case to assess the present algorithm. The Reynolds
number of the flow, based on the cylinder diameter is Re=40. A single block C-type grid
(shown in Figure 7.2) consisting of 101¥201 nodes, is used to discretize the domain. A

STREAM Version 4.5.2 137

second-order upwind scheme is used for the convection terms in the momentum equation.
The time step size used is 1.0 × 10−5 . The main characteristic of the flow field is a pair of
counter-rotating recirculation regions behind the cylinder. The length of these recirculation
regions grows from zero to a constant value when the steady state is reached. Experimental
data is used to compare the evolution of the recirculation region as we march in time, as
shown in Figure 7.3.

Figure 7.2. Computation of steady flow past a circular cylinder; velocity vectors showing the
recirculation region

STREAM Version 4.5.2 138

 2.5

Recirculation Length 2

1.5

PISO Computation
Experiment

0.5

0
0 5 10 15 20 25
Time
Figure 7.3. Development of the recirculation region.

7.5.2 Unsteady Incompressible Flow Past a Circular Cylinder

Although the history of development of the recirculation regions behind the cylinder
in the previous test case provides a validation of the algorithm with experimental data, the
next case allows the validation for a truly unsteady flow. This case involves the flow past a
circular cylinder at Re=100. To demonstrate the performance of the algorithm for a
multiblock grid, the domain is discretized using eight blocks, as shown in Figure 7.4. The
grid consists of 100,000 nodes. The time step size used is 1.0 × 10 −5 .
To validate the algorithm for time-dependent computations, the unsteady problem of
flow over a cylinder is studied. The Reynolds number investigated is 100, at which the flow
is laminar and exhibits the well-known vortex shedding phenomenon. A perturbation is given
to the initial flow field to quickly reach the periodic oscillatory flow pattern. The influence of
the imposed boundary conditions on the time-dependent behavior is studied. The
computational domain and the imposed boundary conditions are shown in Figure 7.4.
Williamson (1989) has performed experiments with an open domain. On the left boundary,
the steady-state inlet condition with constant horizontal velocity is adopted. On the right, the
outlet condition employing zero gradient extrapolation of the velocity and specified pressure
is invoked. Several conditions on the top and bottom boundaries are tested. The inlet
boundary condition imposes the velocity, no-slip boundary condition sets all the velocity
components to zero, the slip boundary condition imposes a zero gradient for the velocity, and

STREAM Version 4.5.2 139

the outlet boundary condition specifies the pressure. Table 1 presents the Strouhal number
corresponding to each case. With the same boundary conditions, the computed and measured
values of the Strouhal number are in excellent agreement.

Figure 7.4. Computational domain and imposed boundary conditions

Table 7-1. Influence of top and bottom boundary conditions on the Strouhal number.
Strouhal number (fD/U∞)

No-slip 0.1267

Slip 0.1658

Inlet 0.1750

Outlet 0.1640

Exp. (Williamson
0.1640
1989)

Figure 7.5 shows the periodic behavior of the cross-stream (v) component of the
velocity obtained by probing at a single point downstream of the cylinder. The v component
of the velocity alternates due to the vortex shedding behind the cylinder. Vortices form on the
surface of the cylinder and shed downstream forming the well known von Karman vortex
street as shown in Figure 7.6.

STREAM Version 4.5.2 140

 Flow over a cylinder, Re=100
0.8
v-velocity
0.6

0.4

0.2

v/U∞ 0

-0.2

-0.4

-0.6

-0.8
0 10 20 30 40 50 60
t/T
Figure 7.5Periodic oscillation of the cross-stream (v) component of velocity.

Figure 7.6. Spanwise vorticity distribution. Snapshot of the vortex shedding and the von
Karman vortex street is shown.

STREAM Version 4.5.2 141

7.5.3 Subsonic Flow Over a Bump

An inviscid flow in a channel with a bump is used to assess the algorithm for
subsonic compressible flows. The grid is composed of three blocks: one upstream of the
bump, one over the bump, and one downstream of the bump. The grid size is 101¥41. The
time step size used is 1.0 × 10−5 . The boundary conditions involve fixing the total pressure,
total temperature and the flow angle at the inlet, and the static pressure at the outlet. The
assigned boundary conditions yield a Mach number M=0.5 at the inlet, for which the flow
remains subsonic throughout the channel. The solution for this case is symmetrical with no
shock waves in the domain, as seen by the Mach number contours plotted in Figure 7.7.

Figure 7.7 Pressure contours for subsonic flow over a bump.

7.5.4 Supersonic Flow Past a Wedge

This case involves a supersonic flow past a wedge placed in a channel. The lower half
of the geometry is shown in Figure 7.8. Two grids are used consisting of 101¥41 and 201¥41
nodes, respectively.. The time step size used is 1.0 × 10−5 . The boundary conditions involve
fixing the magnitude of all the variables at the inlet and the static pressure at the outlet. The
assigned boundary conditions yield a Mach number M=2.3 at the inlet, for which the flow
remains supersonic throughout the channel. The solution contains a shock wave generated at
the leading edge of the wedge and subsequently reflected at the lower channel wall and then
the upper (wedge) boundary (Figure 7.9).
leading edge of the wedge

10.94 0

M 1 + 2.9 M 2 + 2.378
p 1 + 0.714 p 2 + 1.529
M 3 + 1.942 M 4 + 1.55
p 3 + 2.934 p 4 + 5.2

Figure 7.8 Schematic of the supersonic flow past a wedge.

STREAM Version 4.5.2 142

 (a) 101¥21 grid.

(b) 201¥41 grid.

Figure 7.9. Pressure contours for the supersonic wedge flow.

7.6 References
Jones, W. P. and Launder, B. E., 1972, “The Prediction of Laminarization with a Two-Equation
Model of Turbulence,” Int. J. Of Heat Transfer, Vol. 15, pp 301-314.
Patankar, S. V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.
Rhie, C. L. and Chow, W. L., 1983, “A Numerical Study of the Turbulent Flow Past an Isolated
Airfoil with Trailing Edge Separation,” AIAA J., Vol. 21, pp1525-1532.
Shyy, W., Thakur, S. S., Ouyang, H., Liu, J. And Blosch, E. 1997, Computational Techniques for
Complex Transport Phenomena, Cambridge University Press, New York.
Thakur, S. And Shyy, W., 1993, “Some Implementational Issues of Convection Schemes for Finite
Volume Formulations,” Num. Heat Transfer, Part B, Vol. 24, pp 31-35
Thakur, S., Shyy, W., Udaykumar, H. S. and Hill, L., 1998, “Multiblock Interface Treatments in a
Pressure-Based Solver,” Num . Heat Transfer, Part B, Vol. 33, pp 367-396.
Thakur, S and Wright, J. 1999, “STREAM: A Computational Fluid Dynamics and Heat Transfer
Code for Complex Geometries. Part 1: Theory. Part 2: User’s Guide,” Streamline Numerics, Inc.,
Gainesville, Florida

STREAM Version 4.5.2 143