Fundamentals of compressible

and viscous flow analysis - Part III

Lectures 6, 7, 8 and Tables

Instantaneous and averaged temperature contours in a shock-boundary layer interaction.

Taken from (Pasquariello et al., Journal of Fluid Mechanics 2017).

Christophe Corre (LMFA - ECL)

January 17, 2019

 Lecture 6
2D compressible flows

Contents
6.1 Oblique shockwaves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 274
6.1.1 Description of the problem / Flow configuration . . . . . . . . . . . . . 274
6.1.2 Solution method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 275
6.1.3 Oblique shock relationships . . . . . . . . . . . . . . . . . . . . . . . . . 278
6.1.4 Physical interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 279
6.2 Weak shock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282
6.2.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282
6.2.2 θ − β − M relationship for a weak shock . . . . . . . . . . . . . . . . . . 284
6.2.3 Weak shock relationships . . . . . . . . . . . . . . . . . . . . . . . . . . 284
6.3 Supersonic (isentropic) compression . . . . . . . . . . . . . . . . . . . 285
6.4 Supersonic expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 287
6.4.1 Flow description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 287
6.4.2 Prandtl-Meyer function . . . . . . . . . . . . . . . . . . . . . . . . . . . 289
6.5 Computation of aerodynamic forces . . . . . . . . . . . . . . . . . . . 292
6.5.1 Exact computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292
6.5.2 Approximate calculation for a thin airfoil . . . . . . . . . . . . . . . . . 295
6.5.3 Numerical simulation of the supersonic flow over a diamond-shaped airfoil 299
6.6 Exercises and problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 301
6.6.1 Exercise #1 : θ − β − M relationship for a weak shock . . . . . . . . . . 301
6.6.2 Exercise #2 : normal velocity jump through a weak shock . . . . . . . . 305
6.6.3 Exercise #3 : velocity (magnitude) jump through a weak shock . . . . . 305
6.6.4 Exercise #4 : regular shock reflection on a solid wall . . . . . . . . . . . 306
6.6.5 Problem #1 : analysis of the supersonic flow around a projectile . . . . 312
 274 2D compressible flows

6.1 Oblique shockwaves

6.1.1 Description of the problem / Flow configuration
Let us consider now steady 2-dimensional (2D) supersonic flows. Let us also consider two
~ 1 ) and (ρ2 , p2 , U
physical states, (ρ1 , p1 , U ~ 2 ). Denoting Ux and Uy the components of the velocity
vector in the Cartesian frame of reference (~i,~j), we have :

~ 1 = (Ux )1~i + (Uy )1~j


 U
 ~
U2 = (Ux )2~i + (Uy )2~j

We wish to determine the relationships between states 1 and 2 when they are assumed to be
respectively the upstream state and the downstream state of an oblique shockwave such as the
one depicted in Fig.6.1.

n
u2
v1 V2
Etat 1 Etat 2
V1 θ
ρ1 p1 M 1 =|| V1 || v2 ρ2 p2 M 2 =|| V2 ||
a1 u1 β a2

onde de choc
oblique
Figure 6.1: Oblique shockwave (=”onde de choc oblique”) separating two constant states 1
(upstream) and 2 (downstream).

Two essential angles appear in Fig.6.1 :

-• angle β is the shock angle. It is the angle between the direction of the incident or upstream
flow and the oblique shock.

-• angle θ is the flow deviation. It is the angle between the direction of the upstream flow
and the direction of the downstream flow.

As will be seen next, the oblique shockwave appears in the flow as a consequence of the supersonic
flow turning.
Let us introduce for the following shock analysis a reference frame attached to the shock defined
as (~n, ~t), with ~n the unit normal vector to the oblique shock oriented from the upstream to the
downstream state and ~t the unit normal vector tangential to the shock. In this reference frame,
the velocity vector is decomposed using the normal component u and the tangential component
6.1. Oblique shockwaves 275

v so that the velocity vectors can be also expressed as :

(
~ 1 = u1~n + v1~t
U
~ 2 = u2~n + v2~t
U

As explained below, this decomposition in the frame of reference linked to the shockwave is
particularly useful to derive the shock relationships between state 1 and state 2.

6.1.2 Solution method

The flow under study satisfies the 2D steady Euler equations. These governing equations are
expressed in their integral form :
 ˆ

 ρU~ ·N ~ dS = 0



 S

 ˆ



[ρU~ (U
~ ·N~ ) + pN
~ ]dS = 0

 S
ˆ







 ρH U ~ ·N~ dS = 0
S

where the notation N ~ has been retained to denote the unit normal vector to the surface S in
order to avoid confusion with ~n the unit normal vector to the shockwave oriented towards the
downstream state.
The decomposition of the velocity vector in the frame of reference (~n, ~t) attached to the shock-
wave is used to express U~ = u~n + v~t.
The surface S (or more precisely in 2D the contour S) is chosen as described in Fig.6.2 and is
therefore such that :

-• either a face (top and bottom faces on both sides of the shockwave included in the control
~ ·N
surface) is aligned with the flow direction so that U ~ =0

~ = ±~n (with the minus sign

-• or a face is parallel to the oblique shockwave so that N
along the upstream left face and the plus sign along the downstream right face) and
~ ·N
U ~ = (u~n + v~t) · N
~ = ±u

Taking into account these properties yields the following form of the governing equations applied
to the surface S displayed in Fig.6.2 :
 ˆ ˆ



 ρ (−u) dS + ρ (u) dS = 0

 Σ1 Σ2

 ˆ ˆ



[ρU~ (−u) + p (−~n)]dS + ~ u + p ~n]dS = 0
[ρU

 Σ1 Σ2

ˆ ˆ







 ρH (−u)dS + ρH udS = 0
Σ1 Σ2
 276 2D compressible flows

Figure 6.2: Choice of a relevant control surface for the study of an oblique shockwave. The
control volume includes the shockwave and its faces are chosen so that they are either aligned
with the flow direction (top and bottom faces on both sides of the shockwave) or parallel to
the shockwave. Σ1 denotes the upstream (left) face while Σ2 denotes the downstream (right)
face, both parallel to the shockwave.
6.1. Oblique shockwaves 277

Since the flow state along Σ1 is state 1 while the flow state along Σ2 is state 2, the mass
conservation yields, denoting L the length of the edges Σ1 and Σ2 of the contour S :

−ρ1 u1 L = ρ2 u2 L

which simplifies into ρ1 u1 = ρ2 u2 .

The momentum conservation equation yields :

− ρ1 u1 (u1~n + v1~t) + p1~n L + ρ2 u2 (u2~n + v2~t) + p2~n L = 0

Projecting this vector equation on ~n and ~t yields :

 ρ1 u21 + p1 = ρ2 u22 + p2


ρ1 u1 v1 = ρ2 u2 v2


Finally, the energy conservation yields :

ρ1 u1 H1 = ρ2 u2 H2

Taking into account the mass conservation identity when expressing the tangential momentum
conservation equation and the energy equation yields the following system of equations to link
states 1 and 2 on both sides of the oblique shockwave :

 ρ1 u1 = ρ2 u2
p1 + ρ1 u21 = p2 + ρ2 u22


(6.1)
v = v2
 1


H1 = H2

We can thus observe that through an oblique shockwave :

-• the tangential (to the shock) component of velocity v is conserved

-• the normal (to the shock) component of velocity u satisfies exactly the set of conservations
laws established for a 1D shockwave in Lecture 4. This means the oblique shock can be
analysed as a 1D shock in the direction ~n. The shock relationships established in 1D can
be applied but using the normal Mach number associated to the upstream state, denoted
(Mn )1 , instead of the Mach number M1 .

The upstream Mach number M1 is defined as :

p
u21 + v12 w1
M1 = =
a1 a1

~ ||.
where w denotes the magnitude or norm of the velocity vector : w = ||U
The normal Mach number (Mn )1 is defined as :

u1
(Mn )1 =
a1
 278 2D compressible flows

6.1.3 Oblique shock relationships

From the geometry of the flow configuration (see Fig.6.1), the shock angle β can be linked to
the velocity components as follows :

u1 u1 a1 (Mn )1
sin (β) = = =
~ 1 ||
||U ~ 1 ||
a1 ||U M1

that is
(Mn )1 = M1 sin (β)

Replacing the Mach number M1 with the normal Mach number (Mn )1 in the 1D shock relation-
ships yields immediately the following oblique shock relationships for the velocity, density and
pressure jumps :

u1 ρ2 (γ + 1)M12 sin2 (β)

= = (6.2)
u2 ρ1 (γ − 1)M12 sin2 (β) + 2

p2 2γ
=1+ (M 2 sin2 (β) − 1) (6.3)
p1 (γ + 1) 1

Note the temperature jump is such that (T2 /T1 ) = (p2 /p1 ) · (ρ1 /ρ2 ), where the perfect gas equa-
tion of state p = ρrT has been used to express T using p and ρ.

The following remarks can be made regarding these oblique shock relationships :

-• in the purely 1D case, it was established a necessary condition for the shockwave to appear
is that the upstream Mach number is larger than 1, that is the upstream flow is supersonic.
The corresponding condition for an oblique shockwave is that the upstream normal Mach
number must be larger than 1 for the oblique shockwave to appear (this is a necessary
condition but not a necessary and sufficient condition : a supersonic flow can develop
without shockwaves occurring in the flow). This condition reads M1 sin(β) ≥ 1 and yields
a minimal shock angle value given by :
 
−1 1
β ≥ sin
M1

The maximum shock angle value corresponds to the 1D shock or normal or straight shock,
π
such that β = .
2
-• Since (Mn )1 ≥ 1, it follows that u1 /u2 > 1, which means the flow deviation takes place to-
wards the oblique shock or, in other words, the flow turns in the direction of the shockwave
(see Fig.6.1) and not away from the shockwave.

-• Relationships (6.2) and (6.3) also indicate that, across the oblique shockwave, the velocity
decreases while the pressure and the density increase - the temperature also increases
through the shockwave.
6.1. Oblique shockwaves 279

-• Regarding the Mach number, the 1D formula between M2 and M1 can be applied replacing
respectively M1 with (Mn )1 = u1 /a1 and M2 with (Mn )2 = u2 /a2 :

γ−1
1+( )(Mn )21
(Mn )22 = 2 (6.4)
γ−1
γ(Mn )21 − ( )
2

Given the definition of the shock angle β and the flow deviation θ, we can write :

(Mn )2 = M2 sin(β − θ)

so that M2 (which will be such that M2 < M1 ) can be computed from (6.4) provided the de-
flection angle θ is known (and the shock angle is also known of course in order to compute (Mn )1 .

The shock angle β, the flow deviation θ and the upstream Mach number M1 are linked by the
following relationship. Since :
u1 u2
tan(β) = and tan(β − θ) =
v1 v2
we can write :
u2 u2 u1 v1
tan(β − θ) = = × ×
v2 u1 v1 v2
|{z} |{z} |{z}
known from (6.2) =tan(β) =1

so that we eventually obtain the so-called ”θ-β-M ” relationship or identity, which links the
upstream Mach number M1 , the shock angle β and the flow deviation θ :

(γ − 1)M12 sin2 (β) + 2

tan(β − θ) = tan(β)[ ] (6.5)
(γ + 1)M12 sin2 (β)

The θ − β − M relationship (6.5) plays an essential role in the analysis of oblique shockwaves. It
is graphically displayed in Appendix A (see plot III) to allow for an easier use. This graphical
representation is analyzed in the next section.

6.1.4 Physical interpretation

The main characteristics of the oblique shock configurations are reported in Fig. 6.3(a). For
the inviscid flow under study, a streamline can be replaced with a wall boundary since the wall
boundary condition for an inviscid flow is a slip condition (the fluid velocity is purely tangential
to the wall). Therefore, the oblique shockwave displayed on the left side of Fig.6.3 can also be
put in the form of the flow configuration displayed on the right side of this same figure, that is
Fig.6.3(b).

This flow configuration corresponds to a supersonic flow (M1 > 1) over a wedge forming an angle
θ with the incoming flow direction. When turning towards its incident direction, the supersonic
flow generates an oblique shockwave. In this ”supersonic flow over a wedge” configuration, the
flow deviation angle θ is known since imposed by the wedge geometry and the incoming flow
Mach number M1 is also known. The θ − β − M relationship allows to compute the shock angle
 280 2D compressible flows

M
2
β M
M θ 2
1
M
1 β
θ

(a) (b)

Figure 6.3: General oblique shockwave configuration in a free flow and equivalent flow
configuration in the presence of a wall boundary (equivalent to a streamline).

β from θ and M1 .

The standard approach for computing the physical state downstream of an oblique shock with
a known upstream physical state and prescribed flow deviation θ is the following one :

-• from M1 and θ, the shock angle β is computed (6.5), more precisely using the graphical
representation of this relationship (see Appendix A plot III), which will be detailed next.

-• with β known, the normal Mach number can be computed as (Mn )1 = M1 sin(β)

-• the shock relationships, of the form (.)2 /(.)1 = f (M1 sin(β)), can then be applied to obtain
state 2 from state 1.

-• to obtain M2 , a special attention must be paid to the fact that formula (6.4) provides
(Mn )2 = M2 sin(β − θ) so that M2 is obtained as (Mn )2 / sin(β − θ).

The relationship (6.5) is plotted in Fig.6.4 for a given upstream Mach number M1 > 1.

Several important comments must be made on this plot :

-• for every flow deviation angle θ below the maximum deviation angle θmax , there are two
possible values for the shock angle β : β1 < β2 .
On the interval β ∈ [sin−1 (1/M1 ), π/2], the quantity sin(β) is an increasing function of β
and the strenght of the shock, defined as the static pressure ratio p2 /p1 , increases with β.
The solution corresponding to β1 is called the weak shock solution, while the solution
corresponding to β2 is called the strong shock solution.
Fig. 6.5) displays both solutions for a fixed deviation angle (fixed wedge angle) θ below
the maximum angle θmax (M1 ).
Both flow configurations are physically admissible.
In the strong shock solution, the downstream Mach number M2 is always such that the
flow becomes subsonic after the shock that is M2 < 1.
In the weak shock solution, the flow remains generally supersonic after the shock, M2 >
1, except when the deviation angle θ becomes close the maximum angle θmax (see also
Fig.6.4).
The occurrence of one shock solution (weak shock instead of strong shock) rather than the
other (strong shock instead of weak shock) may depend on the downstream flow conditions
(pressure level). For instance a strong shock can develop in the inlet of a ramjet or a
6.1. Oblique shockwaves 281

Déflection du choc θ

M =1 Μ donné
θmax 2 1

Μ >1 Μ <1
2 2

β β
1 2
0
-1
sin ( 1/M ) π/2 Inclinaison β
1

Figure 6.4: Graphical representation of the θ − β − M relationship (6.5) for a given upstream
Mach number M1 .

scramjet depending on the downstream pressure level. However, in general, the weak shock
solution is the one which is physically selected - because it is the most stable solution.
Consequently, except in a specific context motivating the selection of the strong shock
solution, it will always be assumed the weak shock solution is the one occurring in the flow
so that the value of the shock angle corresponds to the lower value β1 .

Choc fort M <1

2 Choc faible M >1
2

M >1 β2 M >1 β1
1 θ < θmax 1 θ < θmax

Figure 6.5: Flow configuration for a strong shock (left) and a weak shock (right) solution in
the case of the supersonic flow over a wedge.

-• for a deviation angle θ such that θ > θmax , there is no physical solution corresponding to an
oblique shock wave, as observed in Fig.6.4. In practice, the physical solution developing in
that case is a curved shockwave which is detached from the wedge corner (see Fig.6.6). The
flow downstream of the shockwave is no longer uniform in that case : a subsonic pocket
forms close to the corner where the shockwave is close to a normal or straight (1D) shock
while the flow remains supersonic further away from the wall where an oblique shockwave
is recovered.

Computing a detached shock solution requires the use of approximate solution technique
 282 2D compressible flows

Choc détaché
M>1

θ > θmax
M >1 M<1
1

Figure 6.6: Detached shock solution for the supersonic flow (M1 > 1) over a wedge with an
angle θ greater than the limit angle θmax (M1 ).

for the 2D Euler equations. There is no exact explicit relationship giving θmax but the
following approximate formula can be used :

4 (M12 − 1)3/2
θmax ≈ √
3 3(γ + 1) M12

The oblique shock angle βmax corresponding to the maximum flow deviation angle θmax
is given by :

1
q
sin2 (βmax ) = [(γ + 1)M12 − 4 + (γ + 1)((γ + 1)M14 + 8(γ − 1)M12 + 16)]
4γM12

Figure 6.7, taken from [3], displays shadowgraphs taken in 1964 by Freeman, Cash and
Bedder at the National Physical Laboratory (USA) for the flow at M∞ = 8.8 over bodies
of revolution with various radius of curvature. For the sharper bodies, an attached shock
solution exists because the flow deviation θ at the body leading edge is below θmax (M∞ )
while for the blunt bodies such an attached shock solution no longer exists and is replaced
with a detached bow shock solution.

-• if θ = 0 then β = sin−1 (1/M1 ) and the ”shock” corresponds to the weakest possible shock
which is actually no longer a shock (finite variation of the flow properties) but a Mach line
through which the flow undergoes an infinitesimal variation, the pressure p for instance
increasing to p + dp through the Mach line.
The angle µ = sin−1 (1/M1 ), known as the Mach angle, is the angle between the incident
flow direction and the Mach line (see Fig.6.8).

6.2 Weak shock

6.2.1 Description
In this section, the flow deviation angle θ is assumed to be small, more precisely small enough
to justify the validity of a simplified form for the θ − β − M relationship.
As displayed in Fig.6.9 for the case of a supersonic flow over a wedge forming an angle θ with
the incident (horizontal) flow direction, if the flow deviation angle θ is sufficiently small, the
shock angle β will be close to the Mach angle µ (which corresponds to the limiting case θ = 0).
6.2. Weak shock 283

Figure 6.7: Attached oblique shock and detached bow shock solutions for the hypersonic
(M∞ = 8.8) flow over sharp and blunt bodies of revolution. Picture taken from [3].

Ligne de Mach

-1
µ = sin (1 / M 1)
M >1
1

Figure 6.8: Mach line configuration

 284 2D compressible flows

The shock is then said to be weak, where the adjective ”weak” must be understood here in the
absolute sense of the work and not ”weak” as opposed to a ”strong” shock when the β1 solution
is selected for the oblique shock angle rather than the β2 solution.
Put in different terms, the weak shock can be considered as weak when θ is small enough; this
is because a small value of β, close to the minimum value µ, yields a low value for the shock
strength, that is the static pressure jump p2 /p1 given by (6.3).

M >1
1

β=µ+ε
µ
θ ε <<1
( θ petit )
Figure 6.9: Weak shock configuration (weak in the sense : small (finite) variations of the flow
properties through the shock.

Using the fact β is close to the Mach angle µ, that is β = µ +  with  << 1, and the fact that
θ is small, a simplified form of the θ − β − M relationship can be derived.

6.2.2 θ − β − M relationship for a weak shock

Let us start from the general θ − β − M relationship (6.5) :

(γ − 1)M12 sin2 (β) + 2

tan(β − θ) = tan(β)[ ]
(γ + 1)M12 sin2 (β)

1
Taking advantage of θ << 1 and tan(β) = tan(µ)+O(), with tan(µ) = p 2 , the following
M1 − 1
simplified relationship can be derived (see Section 6.6 for a detailed demonstration proposed as
an exercise) :
M2
 
2 2 γ+1
M1 sin (β) − 1 ≈ p 1 θ (6.6)
2 M12 − 1

It is particularly interesting to notice it is possible, for a weak shock, to express the shock Mach
angle β as an explicit function of the flow deviation angle θ and the upstream Mach number
M1 .

6.2.3 Weak shock relationships

Formula (6.6) is now exploited in order to derive simplified shock relationships valid for weak
shocks. Since the upstream normal Mach number (Mn )1 = M1 sin (β) is explicitely known
from M1 and θ, the pressure jump for instance (the jump of any flow quantity in fact) can be
expressed as a function of M1 and θ in the case of a weak shock. Injecting the relationship
(Mn )21 = f (M12 , θ) given by (6.6) into the general shock relationships of the form (·)2 /(·)1 =
6.3. Supersonic (isentropic) compression 285

g((Mn )21 ) yields immediately :

p2 − p1 ∆p γM 2
= = p 21 θ (6.7)
p1 p1 M1 − 1

Let us emphasize again the relationship (6.7) can be used, instead of the general relationship
(6.3) only for weak (weak) shock, associated with small values of the flow deviation angle θ.
How small the value of θ must be for (6.7) to provide a good approximation of (6.3) is explained
in one of the exercise proposed in Section 6.6.
To compute the normal velocity (u) jump, we proceed in the similar way, injecting (6.6) in (6.2)
to obtain eventually (the detail of the calculation is proposed as another exercise in Section
6.6) :
u1 M2
= 1 + p 21 θ (6.8)
u2 M1 − 1

The variation of the velocity norm or amplitude w = ||U ~ || through a weak oblique shock can
also be explicitely computed as follows (see again Section 6.6 for the detailed calculation) :

w2 − w1 ∆w −θ
= =p 2 (6.9)
w1 w1 M1 − 1

Note the flow deviation angle θ is counted positively so that the static pressure variation through
the weak shock is a static pressure increase while the velocity (normal velocity or velocity
magnitude) decreases through the weak shock. Note that density and temperature, like the
pressure, increase through a weak shock : though weak, the weak shock remains a shock which
compresses the flow.

6.3 Supersonic (isentropic) compression

Up to now, we have considered supersonic flows over a wedge forming a finite angle θ with the
incident or incoming flow direction. An oblique shock develops in the flow (if θ < θmax (M1 ))
from the wedge corner, which compresses the flow in an adiabatic way but definitely not in an
isentropic way (see Fig.6.10 (a)).
If the wedge is replaced with a series of small wall segments where each segment displays a small
deviation angle ∆θ with respect to the previous one, a weak shock forms at each change of slope
for the wall, with small (but finite) pressure and velocity variations through each weak shock
(see Fig. 6.10 (b)). Each weak shock separates two flow regions in which the flow properties are
uniform. The flow state in a region located upstream of one of the weak shocks determines the
downstream flow state which will become the upstream flow state of the next successive weak
shock.

Pushing the analysis to the limit, we can imagine the small finite slope variation ∆θ becomes
infinitesimal, equal to dθ (see Fig.6.10 (c)). In that case, there is no longer a finite variation of
the flow properties at each change of slope but an infinitesimal variation only, corresponding to
the variation through a Mach line. Along each Mach line the flow properties are constant but
across the Mach line an infinitesimal variation exists which eventually yields a finite variation
 286 2D compressible flows

lignes de Mach
chocs faibles

M M M
1 1 1

θ ∆θ

(a) (b) (c)

Figure 6.10: From an oblique shock configuration (left) throug a series of successive weak
shocks (middle) to a series of Mach lines (right).

when integrated over all the Mach lines created between the upstream state 1 and the down-
stream state 2. The flow across each Mach line can be considered isentropic so that the overall
flow, which is a compression, corresponds to a supersonic isentropic compression.

Going to the limit ∆θ → 0 for the flow deviation translates for instance into the fact the small
finite variation of the velocity magnitude expressed by (6.9) becomes a differential relationship :

dw −dθ
=√ (6.10)
w M2 − 1

Similarly, the infinitesimal pressure increase (dθ > 0) is given by :

dp γM 2
=√ dθ (6.11)
p M2 − 1

Note that in formulae (6.10) and (6.11) the quantity M denotes the local Mach number just
upstream of the current Mach line through which the infinitesimal flow variation is computed,
with dθ corresponding to the infinitesimal flow deviation through the Mach line.

Mach line concept

Let us go back for a moment to the concept of Mach line and let us imagine a source of
perturbation moving in space with a supersonic velocity, say equal to twice the local speed of
sound V = 2a. In other words, this source of perturbation propagates in the flow with a velocity
which is twice the velocity of the perturbations generated by this source since a perturbation
propagates with the speed of sound a.
During a given time t, the source moves from a point A to a point B, covering say a unit distance
(see Fig.6.11). When the source arrives at point B, the perturbations initially created at point
A have been propagated isotropically on a distance 1/2 since a = V /2; a similar analysis can be
performed for perturbations created and emitted at intermediate positions between A and B.
The resulting flow picture is displayed in Fig. 6.11.
The envelope of the perturbations correspond to the wave designed as the Mach line. The
6.4. Supersonic expansion 287

B A
Vt
at

Figure 6.11: Development of a Mach cone.

angle between this ligne and the direction of propagation for the source is given by :

at 1
sin(µ) = =
Vt M

hence µ = sin−1 (1/M ) and µ corresponds to the Mach angle.

The practical generation of a Mach cone can be observed in Fig.6.12. The shadowgraph dis-
played on this picture shows a sphere shot past a plate perforated by a line of small holes.
The pressure of the bow wave in front on the sphere produces below the plate a Mach cone,
corresponding to the envelope of a series of expanding acoustic waves . The angle of the Mach
cone can be measured to find the flow Mach number, which is seen to be equal to 3.

If the ”perturbations” propagating in the flow are no longer weak perturbations of acoustic type
but much stronger perturbations such as, for instance, those produced by a cone moving into
still air with a velocity larger than the local speed of sound hence a Mach number M > 1,
an oblique shockwave develops in the flow, which forms an angle β > µ with the incoming
flow direction (see Fig. 6.13). Note the symmetric flow over a cone is readily deduced from
the previous flow over a wedge : the wall boundary is replaced with the streamline going from
infinity to the stagnation point. Note also the cone could be non-symmetric (see Fig. 6.13) : in
that case the upper and lower attached shocks (assuming θ and θ0 less than θmax (M1 )) can be
analyzed independently from one another.

6.4 Supersonic expansion

6.4.1 Flow description
It has been observed in the previous sections that when a supersonic flow takes place over a
concave wall surface, a compression develops. If the flow must turn rapidly, for instance on a
 288 2D compressible flows

Figure 6.12: Sphere flying at a supersonic speed over a perforated plate. Shadowgraph from
U.S. Army Ballistic Research Laboratory, taken from [3]. The Mach cone can be observed
when looking closely at the lower part of the picture, below the flat plate.

µ
Μ >1
1 β

β Μ >1 θ
1 β
µ
β’
θ’

Figure 6.13: Left : symmetric supersonic flow over a sharp body. Right : non-symmetric
supersonic flow.
6.4. Supersonic expansion 289

wedge imposing a brutal change of direction to the supersonic flow, this compression is no longer
an isentropic process but becomes an oblique shock.
If the supersonic flow takes place over a convex wall surface, Mach lines form as in the case of
the supersonic flow over a concave wall but since these Mach lines are now diverging from one
another, instead of converging to a focal point the flow structure as observed for a compression
over a concave surface, the flow corresponds to a supersonic expansion, which is an isentropic
process (see Fig.6.14). Through each Mach line, an infinitesimal variation of the flow properties
takes place : elementary increase of the velocity, elementary decrease of pressure (as well as
density and temperature). When integrating over the Mach lines, the net global result of the
flow process is : M2 > M1 and p2 < p1 , with upstream state 1 and downstream state 2 connected
through an isentropic process.
If the flow takes place on a convex corner, with a brutal change of flow direction imposing to
the incoming flow a deviation away from its upstream direction, the isentropic expansion is
centered on the corner, with an expansion fan (set of Mach lines) issued from the corner (see
Fig.6.14). This qualitative description being established, it remains to determine how state 2
can be computed from state 1.

µ
1

M1 Μ1

µ1 µ
2

M2
µ2
M2

Figure 6.14: Left : supersonic expansion over a convex wall with a regular change of slope.
Right : supersonic expansion over a convex corner (expansion centered on the corner).

6.4.2 Prandtl-Meyer function

The state downstream of the expansion results from the superposition of the elementary or in-
finitesimal variation of the flow properties through each Mach line of the expansion fan. The
identity (6.10) previously established in 6.3 for the analysis of a supersonic compression also ap-
plies to the isentropic supersonic expansion process, provided an adequate convention is adopted
regarding the sign of the flow deviation angle θ. If θ is still counted positively in the case of
a convex corner (see Fig. 6.15), the expected elementary increase of the velocity magnitude w
produced by an elementary flow deviation dθ is given by :

dw dθ
=√ (6.12)
w M2 − 1
 290 2D compressible flows

When compared with the formula (6.10) previously established for a supersonic compression :

dw −dθ
=√
w M2 − 1

we note the only difference lies in the minus sign in front of the elementary flow deviation dθ
(always counted positively). For a concave surface generating a compression, the elementary ve-
locity change is a decrease of the magnitude while for a convex surface generating an expansion,
the elementary velocity change is an increase of the magnitude.

Identity (6.12) can also be expressed as :

p dw
dθ = M2 − 1 (6.13)
w

In order to account for a finite change θ of flow direction (away from the incoming flow direc-
tion), the above relationship must be integrated, which implies the relative variation of velocity
magnitude should be expressed as a function of dM and M . If this is achieved, it will be pos-
sible to integrate from the zero deviation angle corresponding to the upstream flow with Mach
number M1 to the θ variation angle corresponding to the downstream flow with Mach number
M2 (see Fig. 6.15). From the known deviation angle θ and upstream Mach number M1 , it will
then be possible to compute the Mach number M2 upstream of the expansion and with M1 and
M2 known, it will be possible to compute any flow quantity ratio ϕ2 /ϕ1 taking advantage of the
fact the flow is isentropic.

µ
1

Μ1

θ (>0)
M2

Figure 6.15: Flow configuration for a supersonic centered expansion over a convex corner.

In order to express dw/w as a function of dM and M , we use the definition w = aM as a starting

point. Differentiating this identity yields :

dw da dM
= + (6.14)
w a M
6.4. Supersonic expansion 291

It remains to express da/a using dM and M . Since the speed of sound of an ideal gas depends
on the temperature only and since there is a connection between the static temperature and the
local Mach number in an isentropic flow through the total temperature, we write :

a2 T γ − 1 2 −1
2 = = (1 + M )
a0 T0 2

γ−1 2
Let us denote f (M ) = 1 + M . Since a2 /a20 = 1/f (M ), it is easy to obtain :
2
da 1 f 0 (M )
=− dM
a 2 f (M )

so that, replacing in the identity (6.14) yields :

dw 2f (M ) − M f 0 (M ) dM
=
w 2f (M ) M

Using eventually the expression of f (M ) and f 0 (M ), the following explicit relationship between
dw/w and dM/M is obtained :

dw 1 dM
= (6.15)
w γ−1 2 M
1+ M
2

Inserting (6.15) into (6.13) yields :

√
M 2 − 1 dM
dθ = (6.16)
γ−1 2 M
1+ M
2

The Prandtl-Meyer function is defined as :

ˆ √
M 2 − 1 dM
ν(M ) = (6.17)
γ−1 2 M
1+ M
2
Integrating (6.16) from state 1 to state 2, that is from M1 to M2 for the Mach number and from
0 to θ for the deviation angle, yields :

ν(M2 ) = θ + ν(M1 ) (6.18)

Note the Prandtl-Meyer function can be analytically defined (from atan functions) but it is
tabulated in Table IV of Appendix A for a (much) more practical use.
Knowing M1 , the corresponding value ν(M1 ) can be obtained from Table IV. Adding θ to ν(M1 )
gives the value ν(M2 ) (greater than ν(M1 )) corresponding to a Mach number M2 read again
from Table IV (and which is also greater than M1 ).

Once M2 is known, state 2 can be entirely defined as follows :

 292 2D compressible flows

-• since the total temperature T0 is constant through the expansion (adiabatic and reversible
process) it is possible to write :
T0 γ−1 2 T0 γ−1 2
=1+ M1 and =1+ M2
T1 2 T2 2
so that the static temperature ratio through the expansion is computed as :

γ−1 2
T2 1+ M1
= 2 (6.19)
T1 γ−1 2
1+ M2
2
Note that since M2 > M1 , it follows that T2 < T1 : the static temperature decreases
through an expansion and so do the pressure and the density.

-• since the expansion is isentropic, the total pressure is also constant through the expansion
so that :
 γ
γ − 1 2 (γ−1)

  γ 1 + M
p2 T2 (γ−1)  2 1
= = (6.20)
p1 T1 γ − 1 2
1+ M2
2
Note again that p2 < p1 : the flow is expanding.

6.5 Computation of aerodynamic forces

6.5.1 Exact computation
Using the oblique shock theory and the description of isentropic expansions (and compressions)
provided in the previous sections, it is possible to compute exactly the aerodynamic force applied
by a supersonic flow over an airfoil positioned in this flow. Let us consider for instance the
diamond-shaped airfoil displayed in Fig. 6.16.

From what we learned about supersonic flows in this Lecture, we know an attached oblique
shock develops from the leading edge of the airfoil (provided it is sharp enough for the flow
deviation θ, corresponding here to the angle , to be less tha θmax (M1 )). More precisely, an
oblique shock develops from the upper leading edge (concave) corner made by the streamline
connecting the upstream flow to the airfoil leading edge and the AB segment or edge of the
airfoil. Symmetrically, the same oblique shock develops from the lower leading edge because
of the flow deviation imposed by the (concave) corner made by the same streamline going to
the leading edge and the edge AD of the airfoil. On the convex corner centered on the vertex
B of the airfoil, made by edges AB and BC, a supersonic (centered) expansion fan is created
which makes the flow turn with an angle θ = 2. The airfoil being symmetric and the flow at
zero incidence, the same expansion is also created at point D on the lower part (intrados) of the
airfoil. Finally, oblique shocks are also created at the trailing edge where the flow turns again
to align with the horizontal direction.
Let us denote respectively 1 the uniform supersonic flow upstream of the airfoil, 2 the (still
supersonic) state downstream of the leading edge oblique shocks and 3 the (supersonic) state
downstream of the centered expansions. State 4 will be the uniform flow downstream of the
 6.5. Computation of aerodynamic forces 293

Etat 2
B Etat 3
p
p2 3
Etat 1 t C
A ε
c Etat 4
M1 > 1 ε
p2 p3
Etat 2 D Etat 3
j

i
Figure 6.16: Symmetric diamond-shaped airfoil in a supersonic flow. ”Etat” stands for ”State”
(physical state).

trailing edge oblique shocks.

Since the flow is assumed inviscid, the aerodynamic force applied by the flow onto the airfoil is
exclusively a pressure force defined as :
ˆ
~
F =− p~ndl (6.21)
airfoil

where ~n denotes the unit normal vector to the airfoil surface pointing toward the flow. Note
that formula (6.21) gives more precisely the force applied on the airfoil per unit length of airfoil
span in the direction normal to the (Oxy) plane.
The aerodynamic force F~ is classically decomposed into a drag component D (projection on the
incoming flow direction) and a lift component L (projection on the direction orthogonal to the
incoming flow direction) :
F~ = D~i + L~j (6.22)
where ~i et ~j are the basis vectors of the Cartesian frame of reference, with ~i aligned with the
incoming flow direction and ~j normal to this direction.
In order to compute F~ , the unit normal vector ~n to the 4 edges or faces of the diamond-shaped
airfoil can be defined as a functionof the angle
. Elementary geometry indicates the components
−sin()
of the normal to the edge AB is and similarly for other faces. It follows that :
cos()
ˆ        
−sin() sin() sin() −sin()
− p~ndl = −[p2 d + p3 d + p3 d + p2 d]
airfoil cos() cos() −cos() −cos()
where d denotes the length of edges or faces AB, BC, CD and DA. Note that p2 is the uni-
form pressure downstream of the leading edge oblique shocks while p3 is the uniform pressure
 294 2D compressible flows

downstream of the centered expansions attached to B and D. After simplification, it comes :

 
~ −2sin()(p3 − p2 )d
F =
0

and since sin() = (t/2)/d (with t the airfoil thickness), the drag and lift are eventually given
by :

D = (p2 − p3 )t
(6.23)
L=0

The lift of a symetric airfoil positioned in a flow with zero incidence (that is aligned with the
airfoil chord) is naturally zero since the stresses applied on the upper part of the airfoil cancel
out with the stresses applied on the lower part of the airfoil.
The drag contribution in the present case is said to be a wave drag as it is created by the occur-
rence of oblique shockwaves and expansion waves along the airfoil surface which yield pressure
levels p2 and p3 distinct from each other and also distinct from the upstream pressure level p1 .
The pressure p2 applied on the first half of the airfoil is obtained by applying the oblique shock
relationships between state 1 and state 2. The pressure p3 applied on the second half of the
airfoil is obtained as the downstream state of a supersonic expansion turning the flow of an angle
θ = 2 with an upstream state corresponding to state 2.

It is customary in fluid mechanics and aerodynamics to work with non-dimensional forces. The
non-dimensional drag coefficient and the non-dimensional lift coefficient are thus defined by
dividing respectively the drag force and the lift force (per unit length in the z direction) by a
reference force (per unit length) computed from the upstream dynamic pressure q1 = 12 ρ1 w12
and the airfoil chord c (let us recall w1 denotes the magnitude of the velocity). Since it is
1
also possible to express the dynamics pressureas q1 = γp1 M12 , the following expressions can
2
eventually be obtained for the drag and lift coefficients :

2L 2D
CL (M1 ) = , CD (M1 ) = (6.24)
γp1 M12 c γp1 M12 c

Application
Let us consider the flow at M1 = 3 over a symmetric diamond-shaped airfoil such that  = 3◦ .
The upstream static pressure is assumed equal to p1 = 1 atm. Let us compute the drag coefficient
of this airfoil.
State 2 is the downstream state of an oblique shock characterized by an upstream Mach number
M1 = 3 and a flow deviation angle θ =  = 3◦ . Using plot III of Appendix A, which graphically
translates the θ − β − M identity, we obtain the value β = 22◦ for the shock angle. Let us
recall the selected value corresponds to the lower value of β. It follows the normal Mach number
(Mn )1 is such that (Mn )21 = M12 sin2 (β) = 1.263. Applying the pressure jump relationship (6.3)
yields p2 = 1.307 atm . The downstream normal Mach number (Mn )2 is obtained from (6.4)
and the downstream Mach number M2 is computed taking into account M2 = (Mn )2 /sin(β − θ)
so that M2 = 2.745. Note that the flow downstream of the oblique shock is still supersonic.
State 3 is computed as the downstream state of a supersonic expansion characterized by the
upstream state 2 and a flow deviation angle equal to θ = 2 = 6◦ . Reading the Prandtl-Meyer
6.5. Computation of aerodynamic forces 295

table IV in Appendix A yields ν(M2 ) = 44.5◦ . Consequently, ν(M3 ) = 50.5◦ and reading again
the Prandtl-Meyer table IV yields M3 = 3.039. It can be checked the flow is indeed accelerated
through the expansion wave. Applying formula (6.20) yields :
 γ
γ − 1 2 (γ−1)

p 3  1 + 2 M2 
=
p2 γ − 1 2
1+ M3
2

so that p3 = 0.837 atm .

The drag coefficient is computed as :

D (p2 − p3 ) t (p2 − p3 )
CD = = = tan()
q1 c q1 c q1
1
or else, since the dynamic pressure q1 is such that q1 = γp1 M12 :
2
 
2 p2 p3
CD = − tan ()
γM12 p1 p1

With the known values of p1 , M1 ,  and γ (p1 = 1 atm, M1 = 3,  = 3◦ and γ = 1.4) and the
computed values of p2 and p3 the exact value of the drag coefficient is found equal to :

CD = 3.91 × 10−3

Assuming an airfoil chord c = 0.2 m, the corresponding value of the drag force is :
1
D= × 1.4 × 101325 × 32 × 0.2 × 0.00391 ≈ 499 N
2
for a unit length in the z-direction (airfoil span).

6.5.2 Approximate calculation for a thin airfoil

Pressure coefficient for a thin airfoil
Let us assume the diamond-shaped airfoil is such that its thickness t is small with respect ot
its chord  and that is is also positioned with an incidence angle which remains small. Under
these assumptions, the flow deviation angle (counted with respect to the incident flow direction)
remains small at any location in the supersonic flow around the airfoil. Consequently, the
relationship (6.7) established for computing the pressure increase through a weak shock can be
applied as well as the relationship (6.11) established for computing the pressure decrease through
a supersonic expansion. Both formulae can be actually gathered into the single expression :

p − p1 γM 2
= p 21 θ (6.25)
p1 M1 − 1

where the reference value for the static pressure is the upstream farfield pressure p1 and the
reference value for the flow direction is the upstream farfield flow direction assumed horizontal.
 296 2D compressible flows

Consequently, θ is the local inclination of the flow with respect to the horizontal.

The pressure coefficient is defined as :

p − p1
Cp = (6.26)
q1

1
with the reference dynamic pressure q1 such that q1 = γp1 M12 .
2
Inserting (6.25) into (6.26), a key (simple) relationship is obtained between the local pressure
coefficient at the airfoil surface and the local flow deviation angle in the case of a supersonic
flow over a thin airfoil :
2θ
Cp = p 2 (6.27)
M1 − 1

The so-called thin-airfoil theory makes use of the approximation formula (6.27) yielding the local
pressure coefficient Cp to obtain simple formulae for the lift and drag coefficient of the airfoil.

Drag and lift for a thin airfoil

Let us assume the geometry of the airfoil under study is defined by a general shape function
y = yU (x) for the upper surface of the airfoil and a general shape function y = yL (x) for the
lower surface of the airfoil (see Fig. 6.17). It is assumed the upper and lower line (the airfoil is
defined in the Oxy plane and extends on a unit length in the z direction) do not cross and are
such that yU (0) = yU (0) and yL (c) = yU (c).

yU (x)

x x=0 x=c
yL (x)

Figure 6.17: Schematic view of the thin airfoil geometry. Note the orientation of the closed
integration contour when computing the aerodynamic force F~ .

The ~
ˆ aerodynamic force applied on the airfoil is still given by the general expression : F =
− p~ndl but we wish to derive an approximate formula for F~ which makes use only of the
airfoil
upstream flow conditions (state 1) and the airfoil geometry defined by yU (x) and yL (x).
p
The elementary length dl along the airfoil contour is defined
r as dl = dx2 + dy 2 . Along the
dy dyU dyU 2
upper airfoil surface, we have : = so that dl = (1 + ( ) )dx. Similarly, for the
dx dx dx
6.5. Computation of aerodynamic forces 297

r
dyL 2
lower airfoil surface, we obtain dl = (1 + ( ) )dx.
dx
The unit normal vector to the curve of equation y = yU (x) is given by :

dyU
!
1 −( )
~n = r dx
dyU 2 1
1+( )
dx
and a similar formula holds for the lower airfoil surface.
It is therefore possible to express F~ as :

ˆ c dyU
! ˆ 0 dyL
!
~ − −
F = −[ pU (x) dx dx + pL (x) dx dx]
0 1 c 1

where pU (x) and pL (x) denote respectively the pressure distribution along the upper and the
lower airfoil surface. From the decomposition of F~ into drag and lift it follows that :
ˆ c
 dyU dyL

 D = (pU − pL )dx

 0 dx dx
 ˆ c

 L= (pL − pU )dx

0

The pressure distributions pU and pL can be expressed respectively as a function of yU (x) and
yL (x) introducing the pressure coefficient and making use of the thin-airfoil formula (6.27). Let
us denote the upper and lower pressure coefficient :
pU − p1 pL − p1
(Cp )U = , (Cp )L =
q1 q1
We can then write :

pL − pU = q1 (Cp )L + p1 − q1 (Cp )U − p1 = q1 ((Cp )L − (Cp )U )

so that the lift force is given by :

ˆ c
L = q1 [(Cp )L − (Cp )U ]dx
0

In order to compute the drag force, we can also write :

dyU dyL dyU dyL dyU dyL
pU − pL = q1 [(Cp )U − (Cp )L ] + p1 ( − )
dx dx dx dx dx dx
The last term in this expression does not contribute to drag since its integral over x from 0 to
c is zero in the case of the closed airfoil under study. It comes eventually :
ˆ c
dyU dyL
D = q1 [(Cp )U − (Cp )L ]dx
0 dx dx
 298 2D compressible flows

Applying now the thin-airfoil approximate formula (6.27) which links the local pressure coeffi-
cient and the local flow deviation with respect to the upstream farfield incoming flow direction,
we obtain :
2 dyU 2 dyL
(Cp )U = p 2 ( ) , (Cp )L = p 2 (− )
M1 − 1 dx M1 − 1 dx
dyU
where we have used the fact the local flow inclination is given by for the upper airfoil
dx
dyL
surface and by − for the lower airfoil surface.
dx
The following thin-airfoil (approximate) expressions are finally obtained for drag and lift :
ˆ c
2q1 dyL 2 dyU 2
D=p 2 [( ) +( ) ]dx (6.28)
M1 − 1 0 dx dx
ˆ c
−2q1 dyL dyU
L= p [ + ]dx (6.29)
M12 −1 0 dx dx
Non-dimensional expressions are readily obtained for CD and CL , dividing D and L by q1 c :
ˆ
2 1 c dyL 2 dyU 2


 CD = p 2 [( ) +( ) ]dx
M1 − 1 c 0 dx dx



ˆ c (6.30)

 −2 1 dyL dyU 4 1
 CL = p 2 [ + ]dx = p 2 (y(0) − y(c))


M1 − 1 0c dx dx M1 − 1 c
with y(0) = yU (0) = yL (0) and y(c) = yU (c) = yL (c) for a closed airfoil surface.

Application
The thin-airfoil theory is applied to the supersonic flow at M1 = 3 over the diamond-shape airfoil
studied in the previous section (such that  = 3◦ ). The only information needed to compute CD
and CL is the definition of the airfoil upper surface yU (x) and the airfoil lower surface yL (x).
Since the airfoil contour is closed and such that y(0) = y(c) = 0 it can be concluded from the
thin airfoil theory that CL = 0, which also corresponds to the exact result.
To compute the airfoil drag coefficient, the slope of the airfoil contour is needed. Whatever the
face of the symmetric diamond-shaped airfoil under study, we have :
dyL 2 dyU 2
( ) =( ) = tan2 ()
dx dx
Inserting into (6.30) yields :
ˆ c
D 2 1
CD = =p 2 2tan2 ()dx
q1 c M1 − 1 c 0

hence
4tan2 ()
CD = p 2
M1 − 1

An immediate calculation yields CD = 3.88 × 10−3 , a thin-airfoil approximate prediction which

differs from the exact result CD = 3.91 × 10−3 , with a relative error below 1%. This excellent
agreement is ensured by the validity of the thin-airfoil theory in the present case : with  = 3◦
the thickness-to-chord ratio is low enough to ensure the airfoil is indeed thin.
6.5. Computation of aerodynamic forces 299

6.5.3 Numerical simulation of the supersonic flow over a diamond-shaped

airfoil
The exact shock-expansion theory or the thin-airfoil theory makes it possible to compute the
aerodynamic force applied on an airfoil or a wing in a supersonic flow, at least as long as the
geometries under study remain simple enough. For more complex geometries, the numerical
simulation of the supersonic flow becomes an essential tool of analysis. Some key ingredients
of such a discrete solution of the Euler equations are briefly reviewed in the present section.
Figure 6.18 displays an example of discretization grid used for a finite-volume calculation of the
supersonic flow over a diamond-shaped airfoil.

0.5

0.4

0.3

0.2

0.1
Y

-0.1

-0.2

-0.3

-0.4

-0.5
0 0.5 1
X

Figure 6.18: Partial view of the computational grid used to compute a discrete (finite volume)
solution of the supersonic flow over a diamond-shaped airfoil.

The integral formulation of the Euler equations is applied with each cell of the grid considered
as a control volume. The 2D Euler equations take the integral form :
ˆ ˛
d
w
~ dΩ + (f~nx + ~g ny ) dΓ = 0 (6.31)
dt Ω ∂Ω

where Ω is a closed domain with closed surface ∂Ω and ~n denotes the unit normal vector to
∂Ω with Carteisian components (nx , ny ). The vector w ~ is the vector of conservative varialbles
T ~
(w = (ρ, ρUx , ρUy ρE) while f and ~g denote the physical flux vectors in the x and y directions,
defined as f~ = (ρUx , ρUx2 + p, ρUx Uy , ρUx H)T and ~g = (ρUy , ρUx Uy , ρUy2 + p, ρUy H)T .
In a finite-volume approach, the volume (in 3D) or the surface (in 2D) Ω is a cell Ωi,j of the
grid, with surface |Ωi,j,k |. The integral form reads in that case :

|Ωi,j,k |(w
~ t )i,j +|Γi+ 1 ,j |F~ i+ 1 ,j − |Γi− 1 ,j |F~ i− 1 ,j
2 2 2 2 (6.32)
~
+|Γi,j+ 1 |G ~
2
i,j+ 1 − |Γi,j− 1 |Gi,j− 1 = 0
2 2 2
 300 2D compressible flows

where the overlined quantities denote the following exact averages :

ˆ
1 d
(w
~ t )i,j = w
~ dxdy
|Ωi,j | dt Ωi,j
ˆ
1
F~ i+ 1 ,j = (f nx + gny ) dΓ
2 |Γi+ 1 ,j | Γ 1
2 i+ ,j 2

and similarly for F~ i− 1 ,j and G

~
i,j± 1 .
2 2

Building a finite-volume (FV) approximation of the exact conservation equation (6.32) means
the averaged physical solution in each cell is now evolved using the following relationship :

|Ωi,j |(w
~ t )i,j ~ 1 1 − |Γ 1 |H
+|Γi+ 1 ,j |H ~1 1
2 i+ 2 ,j i− 2 ,j i− 2 ,j
~ 2 1 − |Γ ~ 2 (6.33)
+|Γi,j+ 1 |H i,j+ i,j− 1 |H i,j− 1
=0
2 2 2 2

~ 1 1 and H
where the numerical flux H ~ 2 1 approximate the physical fluxes on the faces Γ 1 ,
i± ,j i,j± i+ ,j 2
2 2
Γi,j+ 1 of the control cell, using the available averaged state w ~ in cell (i, j) and in the neighboring
2
cells (i ± 1, j), (i, j ± 1) but also (i ± 1, j ± 1) (see Fig. 6.19).

Γj,k+1/2

Γj+1/2,k
j,k

j-1,k

j,k-1
Figure 6.19: Control cell for the FV discretization of the Euler equations. Cell (j, k) also
corresponds to cell (i, j) in the text (so that j → i and k → j).

0
If an initial state w
~ is defined in each cell, formula (6.33) makes it possible to march in time
until a steady state is reached for the approximate discrete solution of the Euler equations. It
remains of course to choose a numerical formula for H ~ 1 and H
~ 2 (Roe numerical flux, Jameson
numerical flux, AUSM numerical flux or others), selecting a flux which offers good properties
both in terms of accuracy and robustness / stability. We will not further detail here a choice of
numerical flux. The scheme used in the present case is a second-order (in space) Roe scheme.
It is applied to the solution of the supersonic flow at M∞ = 1.5, M∞ = 2 and M∞ = 3 over
a diamond-shaped airfoil with unit chord length and a thickness such that t = 0.07 c. The
calculations are performed in a grid containing 120 points in the i direction of the grid and 25
points (only) in the perpendicular direction j. The pressure contours obtained for each of the 3
flow configurations are displayed in Fig. 6.20. The fact the computational or discretization grid
6.6. Exercises and problem 301

M1 (CD )thin airfoil (CD )simulation

1.5 0.01128 0.01118
2.0 0.01750 0.01742
3.0 0.00693 0.00687

Table 6.1: Comparison between the thin-airfoil solution and the discrete solution for the 3
supersonic flow configurations under study.

is rather coarse in the j direction is visible on the capture of the oblique shocks which develop
from the leading edge of the airfoil: when moving away from the airfoil surface (where the
grid is rather well refined in the j direction) the grid rapidly becomes less and less refined so
that the shocks are less and less well resolved with the numerical dissipation that significantly
diffuses or spreads the flow discontinuities.
The same calculations are also performed in a grid which is refined along the j direction with 3
times more points in this direction. It is clearly observed in Fig.6.21 this grid refinement reduces
the effects of the numerical dissipation on the shock capturing properties of the numerical solver.
Note the pressure levels appearing in the legend of Fig. 6.20 correspond to a non-dimensional
1
static pressure, such that the non-dimensional upstream pressure is equal to p1 = .
γM12
The numerical simulation gives access to the values of the conservative and primitive variables
in every cell of the computational grid. It is therefore possible in particular to estimate the
pressure at the airfoil surface and perform a numerical integration to estimate the aerodynamic
force applied on the airfoil and deduce a discrete estimate of the drag. The computed values
for the 3 tested values of M∞ using the refined 120 × 150 grid are provided in Table 6.1 along
with the estimated values provided by the thin-airfoil theory. A very good agreement between
the two series of results can be observed.
It is also possible to study graphically the computed shock angle for the leading edge oblique
shockwaves in the FV simulations. From the plot III provided in Appendix A, it can be checked
that for a flow deviation of θ = 4◦ , the shock angle decreases with the Mach number to go from
β ≈ 46.5◦ at M1 = 1.5 to β ≈ 33.5◦ at M1 = 2 and β ≈ 22.25◦ at M1 = 3. Looking at the
leading edge oblique shocks on the plots provided in Fig.6.21, it is possible to estimate the shock
angle β for the FV numerical simulation of the flow over the airfoil. We find for M1 varying
from 1.5 to 3 the successive values β ≈ 46.5◦ , β ≈ 33.6◦ and β ≈ 21.8◦ which are in very good
agreement with the exact shock theory.

6.6 Exercises and problem

6.6.1 Exercise #1 : θ − β − M relationship for a weak shock
• Show that for a weak shock (weak being used here in the ”absolute” meaning of the word and
not as opposed to a strong shock solution) the general θ − β − M relationship simplifies into :
γ+1 M2
M12 sin2 (β) − 1 = ( ) p 21 θ
2 M1 − 1
Note : this exercise is an illustration of a technique currently applied by the engineer : identify
small parameters and simplify complex relationships through relevant limited developments in-
 302 2D compressible flows

M=1.5

0.5
Pression
0.39269
0.3829
0.37311
0.363321
0.353531
Y

0.343741
0 0.333952
0.324162
0.314372
0.304583
0.294793
0.285003
0.275214
0.265424
0.255634
-0.5

0 0.5 1 1.5
X

M=2

0.5
Pression
0.227189
0.220782
0.214375
0.207968
0.201561
Y

0.195154
0 0.188747
0.18234
0.175933
0.169525
0.163118
0.156711
0.150304
0.143897
0.13749
-0.5

0 0.5 1 1.5
X

M=3

0.5
Pression
0.110409
0.106364
0.102319
0.0982735
0.0942285
Y

0.0901834
0 0.0861383
0.0820933
0.0780482
0.0740032
0.0699581
0.065913
0.061868
0.0578229
0.0537778
-0.5

0 0.5 1 1.5
X

Figure 6.20: Pressure contours of the discrete solution (120 × 50 grid) for the supersonic flow
over a diamond-shaped airfoil with t/c = 0.07 for a far-field Mach number equal to 1.5 (top), 2
(middle) and 3 (bottom).
6.6. Exercises and problem 303

0.5
Y

-0.5

0 0.5 1 1.5
X

0.5
Y

-0.5

0 0.5 1 1.5
X

0.5
Y

-0.5

0 0.5 1 1.5
X

Figure 6.21: Pressure contours of the discrete solution (refined 120 × 150 grid) for the
supersonic flow over a diamond-shaped airfoil with t/c = 0.07 for a far-field Mach number
equal to 1.5 (top), 2 (middle) and 3 (bottom).
 304 2D compressible flows

volving these small parameters.

 Let us start from the general θ − β − m identity (6.5) :

(γ − 1)M12 sin2 (β) + 2

tan(β − θ) = tan(β)[ ]
(γ + 1)M12 sin2 (β)

If the shock is assumed to be weak, that is such that weak pressure variations take place through
the shock, this corresponds necessarily to a flow where the flow deviation angle θ is small so
that a limited (first-order) development can be performed for tan(β − θ) using θ as a small
parameter :
θ
tan(β − θ) ≈ tan(β) − θ[tan(β)]0 = tan(β) −
cos2 (β)
Inserting into the general θ − β − M identity yields :

θ 2 − 2M12 sin2 (β)

tan(β) − = tan(β)[1 + ]
cos2 (β) (γ + 1)M12 sin2 (β)

Simplifying tan(β) in the left-hand-side and the right-hand-side and isolating the quantity
M12 sin2 (β) − 1 in the left-hand-side yields :

(γ + 1) 2
M12 sin2 (β) − 1 = M1 tan(β)θ
2

Now, if the shock is weak (in the absolute sense), this means the shock angle β is close to the
Mach angle µ that is β = µ +  with  << 1.
It is therefore possible to write :

tan(β)θ = tan(µ)θ + H.O.T

where the Higher Order Terms are at least of order 2 since resulting from the product of the
small parameter θ with the successive powers of the small parameter .
Since, by definition of the Mach angle, sin(µ) = 1/M1 we can also compute :
q
tan(µ) = 1/ M12 − 1

so that
θ
tan(β)θ ≈ p 2
M1 − 1

and replacing in the expression of M12 sin2 (β) yields :

(γ + 1) M2
M12 sin2 (β) − 1 = p 1 θ
2 M12 − 1

which corresponds to the expected weak shock θ − β − M relationship. 

6.6. Exercises and problem 305

6.6.2 Exercise #2 : normal velocity jump through a weak shock

• Show that for a weak shock the normal velocity ratio u1 /u2 through the shock is such that :

u1 M2
= 1 + p 21 θ
u2 M1 − 1

Both the general relationship giving u1 /u2 for an oblique shock which is not necessarily weak
and the θ − β − M identity for a weak shock will be assumed known.

 Since u1 /u2 is given, in the general case, by :

u1 (γ + 1)(Mn )21
=
u2 (γ − 1)(Mn )21 + 2
and since the upstream normal Mach number is given by :

(Mn )21 = 1 + φ(M12 )θ

γ+1 M2
for a weak shock (θ assumed small), with φ(M12 ) = p 1 , it follows that the normal
2 M12 − 1
velocity jump for a weak shock is such that :

u1 (γ + 1)(1 + φ(M12 )θ) 1 + φ(M12 )θ

= =
u2 2 + (γ − 1)(1 + φ(M12 )θ) γ−1
1+ φ(M12 )θ
γ+1
where the weak shock normal Mach number expression has been inserted into the general normal
velocity jump.
Since θ is a small parameter and since for x small, a first-order limited development yields
1/(1 + x) ≈ 1 − x + H.O.T., it is possible to write :
u1 γ−1 γ−1
= (1 + φ(M12 )θ)(1 − φ(M12 )θ) = 1 + φ(M12 )(1 − )θ + H.O.T.
u2 γ+1 γ+1
where the limited development is stopped at first-order.
Replacing φ with its expression yields :

u2 M2
= 1 + p 21 θ 
u1 M1 − 1

6.6.3 Exercise #3 : velocity (magnitude) jump through a weak shock

• Show that for a weak shock the jump of the velocity magnitude w through the shock is given
by the simplified expression :
w2 −θ
−1= p 2
w1 M1 − 1
The identity established in the previous exercise for the normal velocity jump through a weak
shock will be assumed to be known.
 306 2D compressible flows

 Let us start from the definition of the velocity magnitude :

w2 ~ 2 ||
||U a 2 M2 a2 u1 1 1 M2 (Mn )1 u2
= = = M2 u2 =
w1 ~
||U1 || a M
1 1 u2 a1 M u
1 1 (M n )2 M1 u1

Since
(Mn )2 = M2 sin(β − θ) and (Mn )1 = M1 sin(β)
it can be deduced the velocity magnitude ratio through a shock satisfies the following identity :
w2 sin(β) u2
=
w1 sin(β − θ) u1
Assuming now that θ is a small parameter, it is possible to write :

sin(β − θ) ≈ sin(β) − θcos(β)

hence
sin(β)/sin(β − θ) = 1/(1 − θcotan(β)) ≈ 1 + θcotan(β) + H.O.T.
or else, since β is close to µ for a weak shock :
q
sin(β)/sin(β − θ) ≈ 1 + θcot(µ) = 1 + θ M12 − 1

In the previous exercise, the ratio u1 /u2 has been computed as a function of the small parameter
θ. This identity can also be written as :
q
u2 /u1 = 1 − (M12 / M12 − 1)θ

where we have used 1/(1 + x) ≈ 1 − x for x small, x being here θ.

Inserting these various identities in the formula giving w2 /w1 eventually yields :

w2 M2
q
= (1 + θ M12 − 1)(1 − p 21 θ)
w1 M1 − 1

which can be expanded and simplified, retaining only the first-order term in θ :
w2 θ
=1− p 2 
w1 M1 − 1

6.6.4 Exercise #4 : regular shock reflection on a solid wall

Let us consider the flow configuration displayed in Fig. 6.22. A horizontal supersonic flow at
M1 > 1 takes place in a channel. An oblique shock with a shock angle β1 develops when the
flow is deviated by a wedge forming an angle θ with the lower wall of the channel. This oblique
shockwave impacts the upper wall of the channel where it is reflected. The incident shockwave
makes the upstream flow (state 1) turn of an angle θ in order to be aligned with (and slip along)
the wedge (state 2). The reflected shock realigns the flow in state 2 with the upper horizontal
wall of the channel (state 3). The angle between the upper wall of the channel and the incident
shockwave corresponds to β1 . The shock reflection is not specular, that is the angle Φ between
6.6. Exercises and problem 307

the reflected shockwave and the upper wall of the channel is not equal to β1 . The present exer-
cise analyzes such an oblique shock reflection in the case of an incoming supersonic flow defined
by M1 = 2.8, p1 = 1 atm, T1 = 300 K and with θ = 16◦ .
• Compute the angle Φ between the reflected shock and the upper wall of the channel and show
the shock reflection is not a specular process.

Figure 6.22: Regular reflection of an oblique shockwave on a solid wall.

 Using M1 = 2.8 and θ = 16◦ with the plot IV of Appendix A, the incident oblique shock angle
can be estimated as β1 ≈ 35◦ (the weak shock solution is selected rather than the strong shock
solution). The normal Mach number (Mn )1 can thus be computed : Mn1 = M1 sin(β) = 1.606.
Applying the oblique shock relationships yields p2 /p1 = 2.85, T2 /T1 = 1.388 and Mn2 = 0.6684.
Since Mn2 = M2 sin(β1 − θ) it follows that M2 = 0.6684/ sin(19◦ ) = 2.053.
At this stage the physical state 2 downstream of the incident shock has been entirely defined.
Now, this same physical state 2 is also the upstream state of the reflected shock. The reflected
oblique shock is such that the upstream Mach number is M2 = 2.053 and the flow deviation
angle is again θ = 16◦ since the flow is realigned with the upper channel boundary when crossing
the reflected shockwave. Using plot III of Appendix A provides β2 = 45.5◦ for the shock angle
of the reflected shockwave. Since Φ = β2 − θ it follows that Φ = 45.5 − 16 = 29.5◦ . Consequently
Φ 6= β1 since β1 = 35◦ . 

• Compute state 3 downstream of the reflected shock.

 When computing state 3, we make use of the upstream normal Mach number, which is
defined with respect to the normal direction to the oblique shock under study. Consequently,
the normal Mach number to consider when applying the shock relationships through the re-
flected shock is the normal Mach number defined as Mn2 = M2 sin(β2 ) = 1.46. It is of course
distinct from the normal Mach number (Mn )2 computed in the previous question which was the
normal Mach number downstream of the incident oblique shock hence using the normal velocity
relative to the incident shock.
The shock relationships applied with (Mn )2 = 1.46 yield p3 /p2 = 2.32, T3 /T2 = 1.294 and
 308 2D compressible flows

Mn3 = 0.7157 for the downstream normal Mach number relative to the reflected oblique shock-
wave. By definition, Mn3 = M3 sin(β2 − θ) hence M3 = 1.45. The pressure p3 is such that :
p3 p2
p3 = × × p1 = 2.32 × 2.85 × 1 = 6.54 atm
p2 p1
Similarly, the temperature T3 is obtained from :
T3 T2
T3 = × × T1 = 1.294 × 1.388 × 300 = 538.8 K
T2 T1
Note again the oblique shock reflection is not specular because even though the incident and
the reflected shock produce the same deviation of the flow (with an angle θ) the strenght of
these shocks is not the same since the upstream pressure is different for the incident and for the
reflected shock.

This shock reflection has been computed using the Ansys Fluent software on a mixed unstruc-
tured grid made of quadrangles and triangles. The computed Mach contours are displayed in
Fig. 6.23. One can clearly observe the three uniform flow regions : 1 with M1 = 2.8 upstream of
the incident shockwave, 2 with M2 = 2.053 between the incident shock and the reflected shock,
3 with M3 = 1.45 downstream of the reflected shock.

1.5

mach-number: 1.5 1.63 1.76 1.89 2.02 2.15 2.28 2.41 2.54 2.67 2.8

1
Y

0.5

0
0 0.5 1 1.5 2
X

Figure 6.23: Mach contours for a regular shock reflection. Incident flow such that M1 = 2.8,
p1 = 1 atm, T1 = 300 K and wedge on the lower channel wall with θ = 16◦ .

Plotting also the streamlines (see Fig.6.24) make very clear the effect of the oblique shocks on
the flow direction : the flow turns upward with θ = 16◦ when crossing the incident oblique shock
while it turns (back) downward with θ = 16◦ when crossing the reflected oblique shock.
6.6. Exercises and problem 309

1.05

0.95

0.9
Y

0.85

0.8

0.75

1.4 1.6 1.8 2

X

Figure 6.24: Mach contours and streamlines for the regular oblique shock reflection problem.

Finally the computed pressure and Mach number distributions displayed in Fig. 6.25 and Fig.
6.26 confirm the good agreement between the numerical discrete prediction of the flow and the
analytical results obtained in the above analysis.

Singular reflection on a solid wall

The adjective ”regular” has been used several times when describing the previous oblique shock
reflection. Indeed, for the reflected oblique shock to exist, the flow deviation angle θ must re-
main below the limit angle θmax associated with the Mach number M2 upstream of the reflected
oblique shock (see Fig. 6.27). Similarly of course the incident oblique shock exists only if θ is
below θmax (M1 ).
It is possible to find flow configurations such that θ < θmax (M1 ) but with θ > θmax (M2 ). In
that case, the incident oblique shock cannot be reflected as an oblique shock starting from the
point of impact of the incident shock on the channel upper surface. In order to make the flow
realign with this upper surface, a normal shockwave is formed close to the upper wall. When
moving away from this wall, the normal shock is going to curve slightly and to interact with the
incident oblique shock, producing a curved reflected shock from a so-called triple point (see Fig.
6.28).

This singular shock reflection phenomenon is illustrated, keeping the geometry displayed in Fig.
6.22 but with an upstream Mach number set now equal to M1 = 1.8 instead of M1 = 2.8 for the
regular case.
For this flow configuration, M2 is found slightly below 1.2. It can be checked on plot III of
Appendix A that θmax (M1 = 1.8) is greater than θ = 16◦ while θmax (M2 = 1.2) is lower than
this same deviation angle value θ = 16◦ . We are therefore meeting the conditions of occurrence
 310 2D compressible flows

6
Pression normalisée par pa

1.3 1.4 1.5 1.6 1.7 1.8 1.9

X

Figure 6.25: Regular oblique shock reflecion. Pressure distribution (in atm) along the
horizontal line y = 0.9 (for a unit channel height). Note the computed shocks are not pure
discontinuities because of the numerical dissipation which tends to spread the shockwaves.
Refining (locally) the grid would improve the discrete solution accuracy with steeper
shockwaves.

2.5
Nombre de Mach

1.5

1.3 1.4 1.5 1.6 1.7 1.8 1.9

X

Figure 6.26: Regular oblique shock reflecion. Mach number distribution along the horizontal
line y = 0.9.
6.6. Exercises and problem 311

Figure 6.27: Conditions of occurrence for a singular oblique shock reflection.

Figure 6.28: Schematic view of a singular shock reflection.

 312 2D compressible flows

of a singular shock reflection. The numerical simulation of the flow in the channel with M1 = 1.8
relies on the discretization of the 2D Euler equations. It allows to observe, as expected from
the theoretical analysis, a normal shock develops on the upper channel wall (see Fig. 6.29 or
6.30). Figure 6.31 also illustrates the horizontal flow direction is preserved on the upper part of
the flow while further away from the upper wall the incident oblique shock and the ”reflected”
(slightly) curved shock turn the flow successively upward with an angle θ then downward with
an angle θ to realign it with the horizontal direction.

mach-number: 0.33 0.47 0.61 0.75 0.89 1.03 1.17 1.31 1.45 1.59 1.73

0.8

0.6
Y

0.4

0.2

0 0.2 0.4 0.6 0.8 1 1.2

X

Figure 6.29: Numerical simulation of a singular shock reflection on a wall. Mach contours.

6.6.5 Problem #1 : analysis of the supersonic flow around a projectile

Questions
Let us consider a projectile moving at Mach number M∞ = 2 in air at rest with a pressure level
p∞ = 1 atm and let us focus on the analysis of the steady flow in the symmetry plane of the
projectile, the geometry of which is schematically displayed in Fig.6.32. The analysis will be
performed in the frame of reference attached to the projectile.

1) • Give a qualitative description of the flow along the projectile surface. This description shall
include a schematic drawing of the Mach lines or centered expansions occurring in the flow at
points A, B, C, D. The variation (increasing / decreasing / constant) of the Mach number and
of the static pressure in the various flow regions must also be provided.

2) • Compute, explaining your thought process, the pressure and the Mach number at points
A, B, C, D.
6.6. Exercises and problem 313

pression_atm: 1.2 1.545 1.89 2.35 2.695 3.04 3.385

0.8

0.6
Y

0.4

0.2

0 0.2 0.4 0.6 0.8 1 1.2

X

Figure 6.30: Numerical simulation of a singular shock reflection on a wall. Static pressure
contours.

0.8

0.6
Y

0.4

0.2

0 0.2 0.4 0.6 0.8 1 1.2

X

Figure 6.31: Numerical simulation of a singular shock reflection on a wall. Streamlines

superimposed with the Mach contours.
 314 2D compressible flows

Figure 6.32: Geometry of the projectile studied in questions 1) and 2). Different scales are
used along the axes x and y. Point B is located at mid-distance of points A and C (distance
measured along the circular arc AC).

3) • Apply again the qualitative description of question 1) when the geometry of the projectile
is modified as indicated in Fig.6.33. The values of the static pressure and the Mach number at
points A, B, C and D are not requested in that case.

Figure 6.33: Geometry of the projectile studied in question 3).

Answers
1)  Since the flow is not deviated from its incoming horizontal direction when passing at
point A, the Mach line issued from A is inclined by the Mach angle µA = sin−1 (1/M∞ ) =
sin−1 (1/2) = 30. The flow deviation produced by the concave wall boundary AC induces an
isentropic compression of the flow, with a series of converging Mach lines issues from the points
forming the AC arc. Across this isentropic compression, the pressure gradually increases while
the Mach number gradually decreases. When the flow reaches point C, it must turn to align with
6.6. Exercises and problem 315

the horizontal wall boundary CD. This flow deviation, away from the incoming flow direction,
corresponds to an isentropic expansion, through which the Mach number increases while the
pressure decreases. From the last Mach line of the expansion fan centered on C, the flow is
horizontal and uniform with a constant Mach number and constant pressure. This qualitative
description of the flow over the projectile is summarized in Fig.6.34. 

Figure 6.34: Schematic view of the key flow features for the supersonic flow over the projectile
geometrically described in Fig.6.32. ”compression isentropique” = ”isentropic compression”;
”détente centrée” = ”centered expansion”; ”lignes de Mach convergentes” = ”converging Mach
lines”.

2)  The Mach number at point A is MA = M∞ = 2 and the pressure at point A is pA = p∞ =

1 atm. To compute the Mach number at point B and at point C, immediately upstream of the
centered expansion,the Prandtl-Meyer relationship is used :

ν(M ) = ν(MA ) − θ

where θ denotes the flow deviation with respect to the horizontal direction (the minus sign
is associated with an isentropic compression : the Mach number decreases in an isentropic
compression). Let θC be the flow inclination at point C, just upstream of the centered expansion
fan. This inclination is also the angle between the segments OA and OC. Considering the right
 316 2D compressible flows

triangle with vertices O, C and the intersection between the horizontal ligne going through C
and D and the segment OA, we can write :

sin(θC ) = L/R = 1/3.236

hence θC = 18◦ . It follows that the flow deviation at point B is θB = 9◦ .

Therefore ν(MB ) = ν(MA ) − θB = 26.5 − 9 = 17.5◦ hence MB = 1.689 and similarly ν(MC ) =
ν(MA ) − θC = 26.5 − 18 = 8.5◦ so that MC = 1.383. The pressure at point B is obtained by
expressing the fact the flow is isentropic between A and B :
!
pB 1 + γ−1
2 MA
2
=
pA 1 + γ−1
2 MB
2

A numerical application yields pB ≈ 1.6 atm. The pressure at point C is similarly obtained :
pC ≈ 2.52 atm.

Through the centered expansion fan, the flow turns (still isentropically) with an angle of 18◦ in
order to be aligned again with the horizontal direction. The uniform state downstream of the
expansion, in particular at point D, is such that it satisfies :

ν(MD ) = ν(MC ) + 18◦ = 26.5◦

so that we recover MD = 2 and pD = pA which is logical since the flow is assumed to be entirely
isentropic (this is an idealized view, neglecting the viscous effects).

3)  For the new geometry described in Fig. 6.33, a brutal deviation of the flow occurs at point
A, with an angle θ = 18◦ . Using the θ − β − M plot provided in Appendix A III, we can check
that the maximum deviation angle θmax associated with an incident supersonic flow at M∞ = 2
is equal to θmax = 23◦ . Since the deviation imposed to the flow by the geometry is smaller than
this limit value, we can deduce that an attached oblique shock appears at point A. Across this
oblique shock, the flow velocity decreases while the pressure increases, according to the oblique
shock jump conditions. When the flow evolves along the convex arc AC, a supersonic expansion
takes place which produces a gradual increase of the Mach number and a gradual decrease of
the pressure. The maximum value of the Mach number in the expansion (and correspondingly
the minimum value of the pressure) is reached at the level of the Mach line issued from point C.
Downstream of this Mach line, the flow is uniform with M = MC and p = pC . This qualitative
description of the flow is summarized in Fig. 6.35.
Fig. 6.36 also provides a photograph of the flow over an ogive cylinder flying at Mach 1.7.
6.6. Exercises and problem 317

Figure 6.35: Schematic view of the key flow features for the projectile geometrically described
in Fig. 6.33.
 318 2D compressible flows

Figure 6.36: Ogive cylinder in free flight at M = 1.7. Photograph from Aerodynamics Range
Facility, U.S. Army Ballistic Research Laboratory. The details of the supersonic flow over the
body of 2-centimeter radius can be observed : i) initially straight bow wave evolving into a
finite coincal tip; ii) weak shockwaves generated by the roughness on the nose; iii) growing
turbulent boundary layer over the cylindrical midsection; iv) expansion fans at the start of the
2-degree coincal boattail and at the base; v) transient shockwaves from the unsteady turbulent
wake. Taken from [3].
 Lecture 7
Laminar boundary layer

Contents
7.1 High-Reynolds number flows along wall boundaries . . . . . . . . . . 319
7.1.1 General qualitative description . . . . . . . . . . . . . . . . . . . . . . . 320
7.1.2 Boundary layer equations . . . . . . . . . . . . . . . . . . . . . . . . . . 322
7.1.3 Key phenomena for a developing boundary layer . . . . . . . . . . . . . 328
7.1.4 Some useful quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . 338
7.1.5 Application to the detection of transition to turbulence . . . . . . . . . 341
7.2 Particular solutions of the boundary layer equations . . . . . . . . . 342
7.2.1 Blasius solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342
7.2.2 Falkner-Skan solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 350
7.3 Integral methods for boundary layer analysis . . . . . . . . . . . . . . 354
7.3.1 von Kármán integral equation . . . . . . . . . . . . . . . . . . . . . . . . 354
7.3.2 Polhausen’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 357
7.3.3 Walz-Thwaites’ method . . . . . . . . . . . . . . . . . . . . . . . . . . . 359
7.3.4 Simple integral method . . . . . . . . . . . . . . . . . . . . . . . . . . . 360
7.4 Exercises and problems . . . . . . . . . . . . . . . . . . . . . . . . . . . 362
7.4.1 Exercise #1 : boundary layer developing along the walls of a wind tunnel 362
7.4.2 Exercise # 2 : Polhausen polynomial and von Kármán equation . . . . 366
7.4.3 Exercise #3 : application of the Walz-Thwaites method . . . . . . . . . 367
7.4.4 Problem #1 : study of the boundary layer developing along an airfoil . 368

7.1 High-Reynolds number flows along wall boundaries

Before devoting Lecture 3 to energy conservation and Lectures 4 to 6 to the study of inviscid
compressible flows, we had first touched in Lecture 2, devoted to momentum conservation, the
concept of boundary layer (see in particular Section 2.3). We review some of these key ideas at
the start of this Lecture, in the context of viscous incompressible flows.
 320 Laminar boundary layer

7.1.1 General qualitative description

In the previous lectures, air was considered as an inviscid fluid, i.e. such that the viscous
constraints are assumed to be zero. The Cauchy stress tensor reduces in that case to :

~~σ = −p~~1

where ~~1 is the unit tensor. The governing equations of this inviscid flow reduce in that case
to the Euler equations. Along with this modeling assumption on the stress tensor, comes the
physical wall boundary for such an inviscid flow : when no viscosity effects exist in the flow, the
fluid slips along a solid wall with a velocity boundary condition at a solid wall S expressed as :
~ (~x ∈ S, t) · ~n = 0,
U

with ~n the unit normal vector to the wall.

Let us consider an inviscid or ideal fluid flow along a flat plate. The plate is assumed horizontal
and the incoming flow is also assumed horizontal with uniform velocity U∞~i or ue~i with ~i the
unit normal vector aligned with the horizontal direction. Since the inviscid fluid is slipping along
the flat plate, the flat plate is merely a streamline for the incoming flow which is not perturbed
by the flat plate. The net result of the inviscid flow assumption is that no drag is created when
the fluid is flowing along the flat plate which is obviously a non-physical result. The typical
behaviour of an inviscid flow along a flat plate is depicted in the left part of Fig.7.1.

In reality, air behaves as a Newtonian fluid with dynamic viscosity µ 6= 0, meaning the stress
tensor takes the following form for an incompressible flow :
~~
~~σ = −p~~1 + 2µD

~~
with D the deformation tensor such that :
→
−
→
−

1 ~ ~  ~ ~ T

D= ∇U + ∇U
2

or equivalently :  
1 ∂Ui ∂Uj
Dij = +
2 ∂xj ∂xi
Because of the viscosity, the real fluid sticks to the wall. A no-slip boundary condition is satisfied
at a wall, which states in the case of a fixed wall boundary S :
~ (~x ∈ S, t) = ~0
U

The fluid x-velocity component near the wall displays the following typical behaviour : it in-
creases from u = 0 at the wall to the value of the external velocity ue , which is reached at a
distance δ from the plate wall (see the right part of Fig.7.1).

The flow region between the plate wall and the inviscid flow region where u takes the velocity ue
it would take down to the wall in the case of an inviscid flow is called the boundary layer. In the
boundary layer, viscous effects are significant, even for a fluid such as air which displays a small
7.1. High-Reynolds number flows along wall boundaries 321

CAS IDÉAL CAS RÉEL

approximation
ue de fluide parfait valable

ue

couche limite u effets visqueux

d’épaisseur δ importants

condition de glissement u=0 (condition d’adhérence)

Figure 7.1: Incompressible flow over a flat plate. Left : case of an inviscid flow. Right : case of
a viscous Newtonian flow. ”Condition de glissement” =”slip condition”. ”Approximation de
fluide parfait valable” =”valid ideal fluid approximation”. ”Effets visqueux importants” =
”significant viscous effects”. ”Condition d’adhrence”=”no-slip condition”.

viscosity (µ ≈ 1.8 × 10−5 kg / m / s for standard pressure and temperature conditions) because
the velocity gradient is significant in the boundary layer where u goes from 0 to ue on a small
distance δ, corresponding to the boundary layer thickness. The boundary layer is also known as
being a shear layer, in which the shear stress (as opposed to the normal stress, see Lecture #
2) is large. The fact the boundary layer thickness δ is small at large Reynolds number will be
explained now from simple physical considerations.

The boundary layer development along a flat plate is governed by two key physical phenomena :
-• the advection of momentum by the flow, with a typical velocity ue . The time needed for
the flow to move on a distance x along the plate is thus such that t ∼ x/ue .

-• the diffusion of momentum by the viscous diffusivity or kinematic viscosity ν = µ/ρ. Since
∂u ∂2u
pure diffusion effects are typically described by an equation of the form = ν 2 , the
∂t ∂y
boundary layer is diffused during the time t on a distance δ in the normal direction to the
plate which is such that 1/t ∼ ν/δ 2 or t ∼ δ 2 /ν.
Eliminating t between the convective and the diffusive estimation yields :
δ2 x
∼
ν ue
Let us introduce the local Reynolds number along the plate defined as :
ue x
Rex = (7.1)
ν

where x is the distance from the plate leading edge (see Fig.7.2).
The evolution of δ can also be expressed using Rex :

δ(x) 1
∝√ (7.2)
x Rex
 322 Laminar boundary layer

which shows that for large Reynolds number flows the thickness of the boundary layer is small
with respect to the characteristic length of the problem, for instance the plate length.

frontière extérieure de la couche limite

ue

δ (x)

x
bord d’attaque de la plaque

Figure 7.2: Typical thickness evolution for a boundary layer developing along a flat plate.
”Bord d’attaque de la plaque” = ”plate leading edge”. ”Frontière extérieure de la couche
limite” = ”external or upper boundary of the boundary layer”.

To better understand how thin a boundary layer can be, let us consider an incompressible flow
of air, assuming a constant kinematic viscosity ν ≈ 10−5 m2 / s. Let us also assume we are
interested in a low-speed flow, which reaches ue ≈ 10 m/s. The quantity ue /ν, corresponding to
a Reynolds number per unit length, is such that ue /ν = 106 m−1 . The thickness of the boundary
layer developing in these conditions along a plane surface of length L = 1 m will be such that :

δ(L) 1
∼
L ReL

with ReL = ue L/ν = 106 . Consequently, δ(L)/L ≈ 10−3 which means the thickness of the
boundary layer at the trailing edge of the 1-meter long plate is of the order of 1 mm, hence
much smaller than the plate length. We keep in mind from now on that a large Reynolds
number flow along a solid boundary implies the development of a thin boundary layer along this
solid boundary.
Note that (7.2) can also be expressed as :

r
νx
δ(x) ∝ (7.3)
ue

which indicates the outer (upper) boundary of the boundary layer displays a parabolic profile,
as schematically displayed in Fig.7.2.
Figure 7.3 provides an example of a boundary layer visualization.

7.1.2 Boundary layer equations

Let us build a set of simplified governing equations to describe the steady incompressible high-
Reynolds number flow along a flat plate. The flow under study is sketched in Fig.7.4.
7.1. High-Reynolds number flows along wall boundaries 323

Figure 7.3: Laminar boundary layer (Re = 6000 based on the diameter of the body)
developing around an axisymmetric body known as a Rankine ogive. ONERA photograph
taken by Werl (1962). Streamlines are made visible by tiny air bubbles in water, illuminated
by a sheet of light in the ogive mid-plane. The boundary layer region corresponds to the white
area close to the ogive wall. Picture taken from [3].

frontière de la
couche limite
région d’écoulement
de fluide parfait

couche limite
x
O
L

Figure 7.4: Laminar boundary layer developing along a flat plate. ”Région d’écoulement de
fluide parfait” = ”Inviscid flow region”. ”Frontière de la couche limite” = ”boundary of the
boundary layer”.
 324 Laminar boundary layer

From the Navier-Stokes equations to the Prandt equations

We know the full model allowing to describe the flow displayed in Fig.7.4 corresponds to the
Navier-Stokes equations, which take the following form for the 2D steady incompressible flow
under study :
 ∂u ∂v
 + =0
 ∂x ∂y


∂2u ∂2u

 ∂u ∂u 1 ∂p
u +v =− + ν( 2 + 2 ) (7.4)
 ∂x ∂y ρ ∂x ∂x ∂y
2 2

 u ∂v + v ∂v = − 1 ∂p + ν( ∂ v + ∂ v )



∂x ∂y ρ ∂y ∂x2 ∂y 2
with u = Ux the x-velocity component, v = Uy the y-velocity component, p the static pressure
and ρ the constant density, ν the constant kinematic viscosity.

We have established in the previous section that two different length scales exist in the flow :
δ and L with the boundary layer thickness δ much smaller than the plate length L for a high-
Reynolds number flow. Still, we can imagine in a first step that we are careless enough to
disregard this observation and retain only L as characteristic length of the problem. Using
also the constant upstream velocity U∞ as reference velocity, we introduce the following non-
dimensional variables : 
 x = Lx̄, y = Lȳ;
u = U∞ ū, v = U∞ v̄;

p = ρU∞2 p̄.

The non-dimensional formulation of the Navier-Stokes equations (7.4) reads in that case :
 ∂ ū ∂v̄
 + =0
∂ x̄ ∂ ȳ



1 ∂ 2 ū ∂ 2 ū

 ∂ ū ∂ ū ∂ p̄
ū + v̄ =− + ( + 2) (7.5)
 ∂ x̄ ∂ ȳ ∂ x̄ Re ∂ x̄2 ∂ ȳ
2 2 v̄


 ∂v̄ ∂v̄ ∂ p̄ 1 ∂ v̄ ∂
 ū + v̄ =− + ( 2 + 2)

∂ x̄ ∂ ȳ ∂ ȳ Re ∂ x̄ ∂ ȳ

with the characteristic Reynolds number of the flow defined as Re = U∞ L/ν.

For the high-Reynolds number flow under study, Re >> 1, all the viscous or diffusive terms
in the right-hand-side of the momentum equations become negligible so that the set governing
equations simplifies into the Euler equations :

∂ ū ∂v̄


 + =0
 ∂∂ x̄ ∂ ȳ



ū ∂ ū ∂ p̄
ū + v̄ =− (7.6)
 ∂ x̄ ∂ ȳ ∂ x̄
∂v̄ ∂v̄ ∂ p̄


=−

 ū
 + v̄
∂ x̄ ∂ ȳ ∂ ȳ

However this set of equations is not consistent with the no-slip (zero velocity) condition along the
fixed wall plate boundary. There is necessarily a flow region of small thickness when Re >> 1,
close to the wall, in which the x-velocity component varies from 0 at the wall (no-slip condition)
7.1. High-Reynolds number flows along wall boundaries 325

to U∞ (or ue ) in the inviscid flow region. Having identified the need to introduce two distinct
characteristic lengths, L in the x direction (along the plate) and δ << L in the y direction
(normal to the plate), we also assume it is relevant to introduce two distinct characteristic
velocities, one for the x-velocity and the other one for the y-velocity. Thus, u is made non-
dimensional using U∞ (or ue ) while v is made non-dimensional using the (yet unknown) velocity
V0 . The Navier-Stokes equations are now made non-dimensional using the following variables :

 x = Lx̄, y = δ ỹ, δ << L
u = U∞ ū, v = V0 ṽ,

p = ρU∞2 p̄

Note also the boundary conditions expressed using the non-dimensional variables :

-• no-slip condition at the wall (y = 0 or ỹ = 0) : ū(x̄, ỹ) = 0 and ṽ(x̄, ỹ) = 0

-• matching between the solution at the upper boundary of the boundary layer and the outer
inviscid solution :
lim ū(x̄, ỹ) = ūe (x̄) , lim p̄(x̄, ỹ) = p̄e (x̄)
ỹ→∞ ỹ→∞

Inserting the above non-dimensional variables into (7.4) yields :


 ∂ ū + L V0 ∂ṽ = 0

∂ x̄ δ U∞ ∂ ỹ







1 ∂ 2 ū L ∂ 2 ū

 ∂ ū L V ∂ ū
0 ∂ p̄
ū + ṽ =− + ( 2 + ( )2 2 )

 ∂ x̄ δ U∞ ∂ ỹ ∂ x̄ Re ∂ x̄ δ ∂ ỹ



1 ∂ 2 ṽ L ∂ 2 ṽ

∂ṽ L V0 ∂ṽ L U∞ ∂ p̄

( 2 + ( )2 2 )

 ū + ṽ =− +


∂ x̄ δ U∞ ∂ ỹ δ V0 ∂ ỹ Re ∂ x̄ δ ∂ ỹ
Let us first analyze the mass conservation equation or continuity equation. Since this equation
expresses the fact the divergence of the velocity vector is zero in an incompressible flow, it is
not possible to simplify the equation by suppressing one of its two terms. Since each factor
(depending on L, δ, U∞ , V0 or Re) multiplying a term in one of the governing equations is an
order of magnitude for this term, only 3 options can be considered : i) this factor is >> 1 so
that the term it multiplies is dominant with respect to the other terms; ii) this factor is << 1
so that the term it multiplies is dominated by the other terms in the equation; iii) this factor
is equal to 1. For the continuity equation, the only possible choice is therefore LV0 /(δU∞ ) = 1
which imposes the following relationship between δ/L and V0 /U∞ :

δ V0
= (7.7)
L U∞

This is a first important information : since the order of magnitude V0 of the y-velocity with
respect to U∞ is the same as the order of magnitude of the boundary layer thickness δ with
respect to L and since we know δ << L for a high-Reynolds number flow (Re >> 1), it follows
that V0 << U∞ when Re >> 1.
 326 Laminar boundary layer

Looking now at the x-momentum equation, we know that the second derivative of u in the
x-direction is negligible with respect to the second derivative of u in the y-direction since the
ratio between the factors associated with these two terms is 1/Re << 1. Moreover, the factor
(1/Re) ∗ (L/δ)2 must be such that :
 2
1 L
=1 (7.8)
Re δ

Indeed, if this factor were much larger than 1, it would mean the only dominant (hence remain-
ing) term in the x-momentum equation would be the second derivative of u with respect to y
which would be physically meaningless (no convection, no pressure effect).
Reversely, if this factor were much small than 1, it would mean there is no diffusion term at all
in the x-momentum equation and we would go back to the issue faced with the Euler system of
equations not being consistent with the no-slip boundary condition at the plate wall.
Using (7.7) and (7.8) it is easy to establish the expressions of δ and V0 (respectively order of
magnitude for the boundary layer thickness and for the y-velocity in the boundary layer) :

L U∞
δ=√ , V0 = √
Re Re

Taking into account these expressions in (7.4) yields the following non-dimensional form :
 ∂ ū ∂ṽ
 + =0
∂ x̄ ∂ ỹ



1 ∂ 2 ū ∂ 2 ū

 ∂ ū ∂ ū ∂ p̄
ū + ṽ =− + + 2 (7.9)
 ∂ x̄ ∂ ỹ ∂ x̄ Re ∂ x̄2 ∂ ỹ
1 ∂ 2 ṽ 1 ∂ 2 ṽ

 1
 ∂ṽ ∂ṽ ∂ p̄
(ū + ṽ ) = − + ( + )


Re ∂ x̄ ∂ ỹ ∂ ỹ Re2 ∂ x̄2 Re ∂ ỹ 2
System (7.9) still corresponds to the full system of the steady incompressible Navier-Stokes
equations, which can be now simplified taking into account the fact Re >> 1 :

∂ ū ∂ṽ


 + =0


 ∂ x̄ ∂ ỹ
∂ ū ∂ ū ∂ p̄ ∂ 2 ū

ū + ṽ =− + (7.10)

 ∂ x̄ ∂ ỹ ∂ x̄ ∂ ỹ 2

 ∂ p̄

 =0
∂ ỹ

A first remarkable result of this dimensional analysis performed for the high-Reynolds number
flow is the following one :

-• the pressure p̄ does not vary along ỹ in the boundary layer, that is in the normal direction
to the plate. At a given x̄ section along the plate, the pressure is constant through the
boundary layer and is equal in particular to its value at the upper boundary of the boundary
layer, that is the inviscid solution p̄e (x̄) for the non-dimensional pressure or pe (x) in the
dimensional case. The following identity is therefore satisfied :
∂ p̄ dp̄e
=
∂ x̄ dx̄
7.1. High-Reynolds number flows along wall boundaries 327

-• the outer pressure distribution p̄e and outer velocity distribution ue (x) are linked by the
non-dimensional Bernoulli relationship (valid for the inviscid steady incompressible flow
with zero body force which is associated with the outer region) :
1
p̄e + ū2e = constant
2
The Bernoulli identity can be differentiated to yield :
dp̄e dūe
+ ūe =0
dx dx
The non-dimensionsal boundary layer equations can thus be expressed as :
 ∂ ū ∂ṽ
 + =0
∂ x̄ ∂ ỹ

2 (7.11)
 ū ∂ ū + ṽ ∂ ū = ūe dūe + ∂ ū

∂ x̄ ∂ ỹ dx̄ ∂ ỹ 2

where the red term is known from a previous inviscid analysis of the flow. For instance, in the
simple case of the flow over a flat plate with zero outer pressure gradient, ūe is constant so that
its derivative with respect to x̄ is equal to zero.

This set of equations is completed with the boundary conditions :


• no-slip at the wall : ū(x̄, 0) = 0 , ṽ(x̄, 0) = 0
• matching with the outer inviscid solution : ū(x̄, ∞) = ūe (x̄)

Going back to dimensional variables, the so-called Prandtl equations describing the boundary
layer are obtained :

∂u ∂v

 + =0
∂x ∂y



∂2u


 ∂u ∂u due
 u +v = ue +ν 2


∂x ∂y dx ∂y
(7.12)





 Boundary conditions :



 • no-slip at the wall : u(x, 0) = v(x, 0) = 0
• matching with the outer inviscid solution : u(x, ∞) = ue (x)


In the following sections, we will look for exact and approximate solutions of the Prandtl equa-
tions. We will see exact solutions can be obtained in particular cases only (corresponding to
specific expressions of the outer velocity distribution ue (x)) while it is possible to derive an ap-
proximate solution method which can be applied whatever the form of the outer inviscid velocity
distribution ue (x).
Prior to this exact and approximate analysis, some mathematical remarks can be made regard-
ing the nature of the boundary layer or Prandtl equations. The set of equations (7.12) is of
parabolic nature : these equations describe an evolution problem (for the velocity field u) along
the x space direction. This parabolic nature of the Prandtl equations has a practical consequence
on how the equations can be solved. It can be proved indeed that if ū(x̄, ỹ) remains positive,
then (7.11) admits a unique solution as soon as it is completed with :
 328 Laminar boundary layer

-• an initial profile of ū in a given section x̄0 : ū(x̄0 , ỹ)

-• a boundary condition on ū at the upper boundary ỹ → ∞

-• a boundary condition on ū at the lower boundary ỹ = 0 (plate wall)

The two requested boundary conditions are ū(x̄, ∞) = ūe (x̄) and ū(x̄, 0) = 0. If a velocity profile
can be determined to initialize the boundary layer calculation, it will be possible to march in
space along the increasing x direction. Marching in space means the velocity profile in a section
can be computed using the velocity profiles in upstream sections. This is possible because of
the parabolic nature of the Prandtl equation : a profile in a given x section is influenced only
by the upstream flow and, once computed, influences the downstream profiles (see Fig. 7.5).

~ condition aux limites connue

y

condition
initiale

sens du calcul = sens de l’évolution

x
condition aux limites connue

Figure 7.5: Boundary conditions for the boundary layer or Prandtl equations seen as a
space-evolution problem (along x). ”condition aux limites connue” = ”known boundary
condition”; ”condition initiale” = ”initial condition”; ”sens du calcul” = ”direction of the
calculation or space integration”.

It will be observed in section 7.1.3 devoted to the analysis of a so-called detached boundary
layer that the condition ū(x̄, ỹ) > 0 is not always fulfilled. In that case, the boundary layer
problem can become ill-posed and it is no longer possible to solve the Prandtl equations using
a space-marching technique.

7.1.3 Key phenomena for a developing boundary layer

Skin friction

In the general case of the incompressible flow of a Newtonian fluid, the viscous stress tensor is
given by :
∂u ∂u ∂v
 
2µ µ( + )
~
~τ = 
 ∂x ∂y ∂x 
∂u ∂v ∂v 
µ( + ) 2µ
∂y ∂x ∂y
7.1. High-Reynolds number flows along wall boundaries 329

The Cauchy stress tensor ~

~σ gathering the pressure stresses and viscous stresses reads :
∂u ∂u ∂v
 
−p + 2µ µ( + )
~σ = −p~
~ ~1 + ~
~τ =  ∂x ∂y ∂x 

∂u ∂v ∂v 
µ( + ) −p + 2µ
∂y ∂x ∂y
In the boundary layer region, the stress tensor can be expressed using the previously defined
non-dimensional variables p̄, ū, ṽ, x̄, ỹ :
 2 ∂ ū 1 ∂ ū 1 ∂ṽ 
−p̄ + √ + )
~
~σ = ρU∞2  Re ∂ x̄ Re ∂ ỹ Re3/2 ∂ x̄ 
 1 ∂ ū 1 ∂ṽ 2 ∂ṽ 
√ + ) −p̄ +
Re ∂ ỹ Re3/2 ∂ x̄ Re ∂ ỹ

It is possible to simplify the expression of ~~σ by :

-• neglecting the viscous terms with respect to the pressure term in the expression of the
normal stresses (diagonal of the stress tensor) since the order of magnitude of these terms
is 2/Re, a term which becomes << 1 when Re >> 1 while the order of magnitude of the
pressure term remains unity.
-• keeping only the dominant viscous terms in the expression of the shear stresses, that is
neglecting the non-dimensional terms multiplied by a factor 1/Re3/2 with respect to√the
term multiplied by a factor 1/Re1/2 only. Note the term multiplied by the factor 1/ Re
must be kept in the expression of the stress tensor, even though it is << 1 when Re >> 1
because if this term is also neglected, this means the viscous friction is not taken into
account in the analysis so that we are precisely missing the estimation of the viscous
effects at the wall !
The stress tensor ~~σ can thus be simplified as follows when Re >> 1 :
 1 ∂ ū 
−p̄ √
~
~σ = ρU∞2 
 1 ∂ ū Re ∂ ỹ 

√ −p̄
Re ∂ ỹ
or else, going back to the dimensional form of the stress tensor :
∂u
 
−p µ
~
~σ =  ∂u
 ∂y 

µ −p
∂y
The local stress vector at each point of the flat plate is given by :
T~ = ~~σ · ~n
with ~n the unit normal vector to the flat plate pointing towards the fluid, that is ~n = (0, 1)T for
an horizontal flat plate :
  
∂u

 µ ∂y
T~ = 

y=0 
 
−pe (x)
 330 Laminar boundary layer

The wall stress vector has two components :

-• a normal stress component equal to the (outer inviscid) pressure pe (x)

-• a shear stress or tangential (to the wall) stress component given by :

 
∂u
τw = µ (7.13)
∂y y=0

This wall shear stress corresponds to the viscous friction of the flow over the plate and
depends directly, through the dynamic viscosity µ, on the slope of the x-velocity profile at
the plate wall.

This wall friction stress or skin friction stress can be made non-dimensional using the outer
(inviscid) dynamic pressure 21 ρu2e as a reference stress value. The local skin friction coefficient
Cf (x) is defined as :
τw(x)
Cf (x) = 1 2 (7.14)
2 ρue

The friction force or drag force which applies on the plate wall of length L in the main flow
direction is obtained by integrating the local stress from the leading edge (x = 0) to the trailing
edge (x = L) of the plate :
ˆ L
D= τw dx (7.15)
0

This friction drag (or drag due to viscous friction) can be made non-dimensional by dividing
D by the product of the reference dynamic pressure 12 ρu2e and the contact surface between the
fluid and the plate. For the 2D flow under study, this contact surface is actually the length of
the plate multiplied by a unit length in the z (spanwise) direction. The friction drag coefficient
is thus computed as :
ˆ
D 1 L
CD = 1 2 = Cf dx (7.16)
2 ρue L
L 0

The friction drag is a key contribution to the drag of modern aicrafts. The reduction of the
friction drag remains a major challenge in aeronautics.

Boundary layer separation

The separation or detachment (from the wall) of a boundary layer may occur when this boundary
layer develops along a wall in the presence of an adverse outer pressure gradient. The outer
pressure gradient dpe /dx is said adverse when it is positive that is when the pressure pe is
increasing following the main flow direction. An increasing outer pressure distribution pe along
x corresponds, from the Bernoulli equation, to a decreasing outer velocity distribution ue along
x. The net effect of the adverse pressure gradient is to slow down the flow in the boundary layer,
since this gradient acts on the boundary layer flow down to the plate wall. In the region close to
the wall, the velocity is small since close to zero (because of the no-slip condition). The adverse
dpe
(positive) outer pressure gradient acting on the flow in this near-wall region can become
dx
strong enough to reduce the momentum of the fluid particles moving close the wall up to a point
7.1. High-Reynolds number flows along wall boundaries 331

where the velocity u of these fluid particles becomes zero and eventually becomes negative (see
Fig. 7.6 left). A negative x-velocity near the wall corresponds to a so-called back-flow region or
recirculating flow configuration. The point where the y derivative of u at the wall becomes zero
is defined, for a 2D flow, as the separation point. Note also the boundary layer becomes thicker
under the effect of an adverse pressure gradient.

GRADIENT DE PRESSION DÉFAVORABLE GRADIENT DE PRESSION FAVORABLE

renversement du sens
de l’écoulement (recirculation)

la couche limite s’amincit

point de décollement

Figure 7.6: Influence of the outer pressure gradient dpe /dx on the development of the
boundary layer. Left : adverse pressure gradien. ”point de décollement” = ”separation point”;
”renversement du sens de l’écoulement” = ”creation of a backflow”. Right : favorable pressure
gradient. ”la couche limite s’amincit” = ”the boundary layer becomes thinner”.

The separation location is defined for a 2D flow as the section where the skin friction or wall
shear stress τw becomes zero : τw = µ( ∂u
∂y )y=0 = 0.

In the recirculating zone which follows the separation, flow instabilities tend to develop which
lead the flow to become turbulent. Separation and transition to turbulence (see below) should
however not be confused. A boundary layer can separate from the wall along which it is devel-
oping without becoming turbulent further downstream.

An ”adverse pressure gradient” configuration is found for instance when considering the bound-
ary layer developing along the wall of a diverging nozzle in which the flow is subsonic. For
such a flow configuration, the inviscid solution (neglecting the development of the boundary
layer) is known to be such that ue decreases along the main flow direction while pe increases
hence dpe /dx > 0. Depending on the level of this adverse pressure gradient and on the level of
momentum in the developing boundary layer, a separation of the boundary layer might occur
at some section of the diverging nozzle.

If the outer inviscid flow is accelerated, ue (x) increases with x, then due /dx > 0 and dpe /dx < 0.
This negative outer pressure gradient corresponds to a favourable pressure gradient since it con-
tributes to accelerate the flow in the boundary layer near the wall. The velocity profile in the
boundary layer becomes fuller and the boundary layer becomes thinner under the influence of a
 332 Laminar boundary layer

favourable outer pressure gradient (see Fig. 7.6 right).

Such a flow configuration can be found for instance in the case of a subsonic flow in a converging
nozzle : the inviscid flow solution in the nozzle is indeed such that ue (x) increases with x.

An example of separated laminar boundary layer is displayed in Fig.7.7. The flow over a blunt
axisymmetric body at Re = 6000 based on the body diameter can be compared with the flow
over the Rankine ogive pictured in Fig. 7.3, where the laminar boundary layer remains attached.
For the blunt body of Fig.7.7 the flow deceleration in the leading edge region is stronger than
the flow deceleration over the Rankine ogive. This stronger outer velocity decrease for ue (x)
corresponds to a stronger (positive) outer pressure gradient dpe /dx which leads to the detach-
ment or separation of the boundary layer, with a well visible recirculation bubble on the upper
part of the ogive - ”recirculation bubble” means that in this flow region of limited extend a back
flow occurs.

Figure 7.7: Separating boundary layer over an axisymmetric body displaying a blunter
geometry than the Rankine ogive pictured in Fig. 7.3. ONERA photograph taken by Werl
(1962). Picture taken from [3].

Boundary layer separation plays a very important role in aeronautics application, for instance
when considering the flow over an airfoil with increasing incidence. Indeed, when the incidence
angle α of the airfoil increases (α can be defined as the angle between the upstream flow direction
and the chord of the airfoil), the positive outer pressure gradient dpe /dx along the airfoil upper
surface (extrados) tends to increase. An example is provided in Fig.7.8 where the distribution
of the pressure coefficient Cp = (p − p∞ )/q∞ along a NACA0012 airfoil is plotted for a flow
at zero degree of incidence and at α = 10◦ of incidence. When the incidence increases from
0◦ to 10◦ , the minimum value of pressure is reached at point A and is followed by a pressure
7.1. High-Reynolds number flows along wall boundaries 333

recovery (pressure increase) which becomes much stronger and corresponds to an higher adverse
(positive) pressure gradient, leading to boundary layer separation.
For a critical value αc of the (positive) incidence angle (in general αc ∈ [10◦ , 20◦ ]), the boundary
layer separates or detaches from the upper airfoil surface, leading to the creation of a large
recirculation zone followed by a thick wake (see Fig.7.9). The location of the separation point
can be assessed by plotting (see Fig.7.10) the skin friction  coefficient
 distribution along the
∂u
airfoil : the point where Cf = 0 corresponds to τw = 0 hence = 0 (with y the normal
∂y y=0
distance to the airfoil surface) and this defines the separation point of the boundary layer (see
also Fig.7.6 left).
The boundary layer separation induces a brutal decrease of the airfoil lift along with an increase
of the airfoil drag : this so-called stall phenomenon is therefore particularly impacting for the
airfoil performance and introduces in practice a constraint on the flight envelope- meaning that
such large values of the incidence angle should simply be avoided.
coefficient de pression Cp=(p-p∝)/(0.5 ρ ∝ V 2∝)

1.2

0.8 incidence 0°

0.6

0.4 incidence 10°

0.2

0
gradient de pression defavorable (dp/dx > 0)
-0.2

-0.4

-0.6
A
0 0.25 0.5 0.75
x/c

Figure 7.8: NACA0012 airfoil. Pressure coefficient distributions at 0◦ and 10◦ of incidence.

Figure 7.11 displays the key features of the flow topology in the neighbourhood of an airfoil
trailing edge, when the boundary layer is detached on the upper surface of the airfoil. In this
example, the separation of the boundary layer and its effects remain ”local” since a closed
separation bubble develops between the separation point D, where the skin friction coefficient
goes through zero for the first time before becoming negative and going through zero again
at reattachment point R. The separation and reattachment points are connected through a
separating streamline which defines the boundary between the recirculation bubble and the flow
outside this bubble which displays no backflow and goes only in the downstream direction.

Another classical example of separated flow corresponds to the flow over the rear part of a body
(see Fig. 7.12). When looking first at the inviscid flow solution for this geometry, it can be
 334 Laminar boundary layer

Figure 7.9: Streamlines over the NACA0012 airfoil at : 0◦ of incidence (left) and 10◦ of
incidence (right).

0.4

incidence 10°
Cf

0.2
incidence 0°

point de decollement : τ w=0

0 0.25 0.5 0.75

x/c

Figure 7.10: Boundary layer separation on the NACA0012 airfoil. Separation point and skin
friction distribution.
7.1. High-Reynolds number flows along wall boundaries 335

Figure 7.11: Flow topology for an airfoil displaying a separation bubble at its extrados, near
its trailing edge.

observed the flow around the edges D and D0 generates strong positive pressure gradient in
the main flow direction, corresponding to a strong decrease of the inviscid flow velocity. When
applied on the boundary layer which develops along the body surface, these adverse pressure
gradients lead to the separation of the boundary layer, precisely at the location of these geometric
singularities (sharp edges D and D0 ). Downstream of these two points, the detached flow sheets
converge towards the axis of symmetry of the body and the two separating streamlines meet
at a reattachment point R where the fluid velocity is zero (stagnation condition). The flow
recirculates at low speed in the region limited by the wall of the body rear part and the separating
streamlines. Downstream of the reattachment point R, the boundary layers merge into a wake.

Figure 7.12: Flow topology over the rear part of a body.

Since the occurrence of boundary layer separation is a key limiting factor for the performance
of many aerodynamic devices (because of loss of lift or increase or drag or else loss of efficiency
for a propulsion system), the accurate prediction of separation remains an important challenge
in aerodynamics. We will see later in this Lecture how separation can be estimated using
 336 Laminar boundary layer

approximate semi-analytical methods such as the method of Thwaites. However, in the general
case, numerical simulation of flows hence CFD plays today a key role into this accurate prediction
of separation, be it through the numerical solution of the boundary layer equations coupled with
an inviscid flow model or through the numerical solution of the full Navier-Stokes equations.

Transition to turbulence

Let us consider a laminar boundary layer on a flat plate. The word ”laminar” refers to the fact
the flow is very regular and can be seen as ”sheets of fluid” (”laminae” in Latin) flowing on top of
each other. Quantitatively, we have seen the thickness of the laminar boundary layer developing
on a flat plate with zero external pressure gradient (case of the exact Blasius solution) increases
as the square root of the distance x measured from the leading edge of the plate. We have also
established the ratio δ(x)/x varies as the inverse of the squareroot of the local Reynolds number.
The experimental observation of a laminar boundary layer development on a flat plate indicates
that for Reynolds numbers between, typically, 3 × 105 and 3 × 106 , the boundary layer loses its
regular aspect (it is no longer ”laminar”) to display a chaotic behaviour and become a turbulent
boundary layer. Such a transition from the laminar regime to the turbulent regime is illustrated
in Fig. 7.13.

Figure 7.13: Left : local Reynolds number (Rex ) less than the critical Reynolds number (Rec );
the boundary layer (BL) is laminar. Right : Rex becomes larger than Rec ; the BL becomes
turbulent. Photograph from the University of Iowa Institute of Hydraulic Research.

The transition of the boundary layer to the turbulent regime is driven by the triggering of flow
instabilies : when the local Reynolds number exceeds a critical value Rec , which depends on
the flow configuration, small perturbations existing in the flow become unstable and their rapid
growth leads to the transition of the flow to a turbulent regime. These small perturbations can
be due for instance to solid wall irregularities (rough surface or at least non-smooth surface) or to
small temperature fluctuations in the ambient medium, . . . The exact mechanisms of transition
are quite well known for some flow configurations (in particular geometrically simple ones such
as the flow over a flat plate) but remain an active subject of research for many other cases. In
the case of a flat plate, the location of the transition starting point along the flat plate depends
on several factors such as the wall roughness, the level of turbulence pre-existant in the incoming
7.1. High-Reynolds number flows along wall boundaries 337

flow, the outer pressure gradient (adverse or favourable), . . . Once the critical Reynolds reached,
the stable laminar flow developing from the plate leading edge becomes unstable with first
two-dimensional waves of Tollmien-Schlichting type propagating in the direction of the mean
flow. The development of these waves can be predicted by the linear stability theory but the
remainder of the transition process is fundamentally non-linear. The Tollmien-Schlichting waves
soon display variations in the direction transverse to the main flow direction and 3D waves as
well as hairpin vortices appear in the flow. These vortices are going to break down and generate
fully 3D fluctatuions in the flow. The most intense fluctuations generate turbulent spots which
will merge in the downstream flow to form a fully turbulent flow. Note that some steps or stages
of this transition process can be bypassed due for instance to an adverse pressure gradient which
generates velocity profile particularly prone to becoming unstable. Figures 7.14 and 7.15 display
visualizations of this laminar-to-turbulent transition over a flat plate.

Figure 7.14: Top : flow at ReL = 20000 over a flat plate. Fully laminar boundary layer.
Bottom : flow at ReL = 100000 over the same flat plate. 2D Tollmien-Schlichting waves
appear in the flow. Photograph taken at ONERA by Werlé (1980) [3].

As will be detailed next, both in the present Lecture 7 and in the next Lecture 8 entirely
devoted to turbulent flows, the properties of a laminar boundary layer and that of a turbulent
boundary layer are very different, both from the viewpoint of skin-friction and separation. It
is therefore important for the engineer to be able to estimate the location of the laminar-to-
turbulent transition in order to apply a laminar description where it is indeed valid and a
turbulent description where needed.
Estimating in an accurate and reliable way the transition remains a difficult task because of
the sensitivity of the phenomenon to a number of physical factors. However, some empirical
 338 Laminar boundary layer

Figure 7.15: Natural transition on a slightly inclined plate. ReL = 100000 and flow with a 1◦
angle of attack. Photograph taken at ONERA by Werlé (1980) [3].

tools are available for the engineer which make use of various characteristic Reynolds numbers
of the boundary layer. For instance, in the case of the boundary layer developing on a flat plate
with zero outer pressure gradient, the transition location xtr (distance measured from the plate
leading edge) is sometimes defined as :

U∞ xtr
Rextr = ≈ 5 × 105 (7.17)
ν

This formula corresponds to typically observed values for flat plates where factors such as wall
roughness or external flow perturbations trigger the transition to turbulence. In the case of
very well controlled experiments, with a smooth plate and a minimum level of external flow
perturbations, the natural transition to turbulence can actually occur at much higher Reynolds
number, up to 5 × 106 in some cases. Better criteria than (7.17) for detecting the onset of
transition to turbulence have been developed using Reynolds numbers based on characteristic
boundary layer thicknesses which we will define now before coming back to their application in
empirical transistion criteria.

7.1.4 Some useful quantities

Boundary layer thickness
When establishing the (laminar) boundary layer equations in 7.1.2, the order of magnitude of
the boundary layer thickness has been introduced : for a flat √plate of length L it was found
the order of magnitude of the boundary layer thickness is L/ ReL with ReL = U∞ L/ν the
characteristic Reynolds number of the flow, assumed to be large (hence δ << L). We have also
explained that when looking at this high-Reynolds number flow, the boundary layer thickness
tends actually to zero. To illustrate this idea, the development of a boundary layer is displayed
in Fig.7.16 for the following flow conditions : length of the plate L = 20 m/s; upstream flow
velocity U∞ = 0.5 m/s, kinematic viscosity ν = 10−5 m2 /s. The Reynolds number ReL is equal
to 106 in that case which means the boundary layer could actually become turbulent at some
location along the plate but this would not change the key message conveyed by Fig.7.16. This
key message is the tremendous scale difference between the length of the plate, which is the
relevant length scale for the advection along x (hence the use of the non-dimensional variable
7.1. High-Reynolds number flows along wall boundaries 339

x̄ = x/L), and the thickness of the boundary layer. With x and y plotted with (almost) the same
scale in Fig.7.16, the boundary layer is barely visible on the picture. The red line corresponding
to the upper or outer boundary of the boundary layer is almost superimposed with the x-axis
(y = 0). The same boundary layer is displayed in Fig.7.17 but using an enlargement in the
y-direction which corresponds to scaling the normal distance y to the plate no longer with the
plate’s length L but with the boundary layer thickness δ - this defines the non-dimensional
variable ỹ = y/δ instead of ȳ = y/δ.

4
y

0 5 10 15 20
x

Figure 7.16: Boundary-layer developing on a flat plate of length L = 20 m with ReL = 106 .
Real scale.

0.3

0.2
y

0.1

2 4 6 8 10 12 14 16 18 20
x

Figure 7.17: Boundary layer developing on a flat plate of length L = 20 m with ReL = 106 .
Enlarging the scale in the y-direction is equivalent to introducing the non-dimensional variable
ỹ = y/δ with δ the characteristic boundary layer thickness.

The definition of the boundary layer retained to draw the outer boundary of the boundary layer
in Fig.7.16 and 7.17 is now clarified. The boundary layer thickness δ(x) at a section x along the
flat plate is defined as the distance to the plate wall at which the x-velocity reaches 99% of the
 340 Laminar boundary layer

value of the outer velocity ue (x) at this section :

u(x, δ99% (x))

= 0.99 (7.18)
ue (x)

In the case of a boundary layer flow over a flat plate with zero outer pressure gradient, dpe /dx =
0, the outer velocity is constant such that ue (x) = U∞ . Note this definition is not really satisfying
because too arbitrary : a threshold value of 95% or 98% or else 99.9% could be retained instead of
99% and this would affect the position of the upper boundary of the boundary layer. Moreover,
using the boundary layer thickness thus defined from a specific point of the velocity profile
provides no relevant global information on the boundary layer velocity distribution. This leads
to retain alternative thickness definition based on integral definitions which take into account
the full velocity profile in the boundary layer.

Displacement thickness
The displacement thickness δ ∗ (also denoted δ1 in some textbooks) corresponds to the loss of
massflow through a section of the boundary layer with respect to the inviscid flow configuration
(where the external or outer velocity ue (x) applies down to the wall because of the slip boundary
condition in the inviscid case). According to this definition, δ ∗ is such that :
ˆ δ ˆ δ
ρe ue dy − ρudy = ρe ue δ ∗
0 0
| {z } | {z }
massflow for the inviscid case massflow for the viscous case
where the notation δ(x) corresponds to δ99% (x) and has been retained for the sake of simplicity.
Assuming the flow is incompressible and homogeneous (ρ = constant), δ ∗ can also be computed
from the normalized x-velocity profile (that is u(x, y) in the boundary layer divided by the outer
inviscid velocity ue (x)) :
ˆ δ(x)  
u(x, y)
δ ∗ (x) = 1− dy (7.19)
0 ue (x)
Note that for an inviscid flow such that ue (x) as assumed here, the x-velocity is equal to ue (x)
for y ≥ δ(x) so that an alternative expression for δ ∗ (x) is :
ˆ +∞  
∗ u(x, y)
δ (x) = 1− dy (7.20)
0 ue (x)

Momentum thickness
The momentum thickness θ (alos denoted δ2 in some textbooks) corresponds to the loss of
momentum in the boundary layer with respect to the case of an inviscid flow, that is (for a flow
with constant density ρ) :
ˆ δ ˆ δ
ρudy × ue − ρu2 dy = ρu2e θ
0 0
so that θ is expressed as :
ˆ δ  
u u
θ= 1− dy (7.21)
0 ue ue
7.1. High-Reynolds number flows along wall boundaries 341

or else
ˆ +∞  
u(x, y) u(x, y)
θ(x) = 1− dy (7.22)
0 ue (x) ue (x)

Shape factor
The shape factor H is defined as the ratio between the displacement thickness and the momentum
thickness :
δ∗
H= (7.23)
θ
This parameter plays an important role in the analysis of the boundary layer (in particular
regarding its separation characteristics); we will see next that H appears in the von Kármán
equation used to provide an estimate of skin friction associated with the boundary layer.

7.1.5 Application to the detection of transition to turbulence

We have previously mentioned that the validity of the criterion (7.17) for predicting the transition
to turbulence of a laminar boundary layer is rather limited. More general criteria have been
proposed, which make use not only of the local Reynolds number Rex (with x the distance
measured from the starting point of development of the boundary layer, for instance distance
along the flat plate measured from the plate’s leading edge) but also of the Reynolds number
Reθ based on the momentum thickness θ(x). Michel has thus proposed in 1952, based on an
attempt to correlate a set of experimental results, to characterize the transition location xtr
using the (approximate of course) formula :

(Reθ )tr ≈ 2.9(Rex )0.4

tr (7.24)

Since Reθ = ue (x)θ(x)/ν and Rex = ue (x)x/ν, the identity (7.24) becomes an equation on
x = xtr such that :
ue (xtr )xtr 0.4
 
ue (xtr )θ(xtr )
≈ 2.9
ν ν
which can be solved once a correlation is available to express θ(x). We will see next how the von
Kármán can be used to this end. In the case of a flat plate with zero outer pressure gradient,
ue (x) = U∞ so that the equation to solve is :

ν 0.6
 
θ(xtr )
≈ 2.9
x0.4
tr U∞
This criterion has been further generalized in 1974 by Cebeci and Smith who proposed to identify
the transition location using the following empirical correlation :
 
22400 0.46
(Reθ )tr ≈ 1.174 1 + (Rex )tr (7.25)
(Rex )tr

Here again, a correlation is needed for θ(x) in order to obtain an equation on xtr which can then
be solved to yield an estimate of the transition location. The following section will provide such
a correlation in the particular case of the boundary layer flow over a flat plate with zero outer
pressure gradient, for which the exact Blasius solution can be obtained.
 342 Laminar boundary layer

7.2 Particular solutions of the boundary layer equations

The Prandtl equations, in their non-dimensional (7.11) or dimensional (7.12) form, with appro-
priate boundary conditions, do not possess a general closed-form solution. The key boundary
condition, specific to the boundary layer under study, is of course the matching between the so-
lution inside the boundary layer and the outside inviscid velocity distribution ue (x). Quasi-exact
solutions can be computed for some specific outer inviscid velocity distributions ue (x). Approx-
imate solutions using the von Kármán Equation can be obtained for any ue (x), in particular for
inviscid velocity distribution around airfoil or blade geometries for instance. Before analyzing
such an approximate solution strategy in Section 7.3.1, we will first focus in the present Section
on exact or quasi-exact solutions which can be obtained for some simple velocity distributions
ue (x) associated with simple geometries (flat plate, wedge). Note that these exact solutions,
though limited in their range of validity, are of interest in particular to assess the more general
approximate or discrete solution strategies.

7.2.1 Blasius solution

Calculation of the Blasius solution
Let us consider a half flat plate (half meaning only the upper surface of the plate will be con-
sidered) positioned with zero incidence in a viscous incompressible flow, with uniform upstream
far-field conditions (u = U∞ , v = 0, p = p∞ ).
The plate is semi-inifinite since extending from x = 0 at its leading edge to a downstream lo-
cation which goes to +∞. Therefore, there is no a priori characteristic length which can be
readily derived from the problem’s geometry. The outer inviscid solution is very simple since it
corresponds to the inviscid flow over a flat plate with no pressure gradient introduced by the
environment of the flat plate so that, because of the slip boundary condition at the plate wall
for an inviscid flow, ue (x) = U∞ . The non-dimensional Prandtl equations or boundary layer
system (7.11) reduces in that case to :

∂ ū ∂ṽ
+ =0


∂ x̄ ∂ ỹ



∂ 2 ū

 ∂ ū ∂ ū
ū + ṽ = (7.26)

 ∂ x̄ ∂ ỹ ∂ ỹ 2



 ỹ = 0 : ū = ṽ = 0
 ỹ → ∞ : ũ = 1

The solution for the non-dimensional x-velocity ū is a priori of the form ū = h(x̄, ỹ), where
p h is an
unknown function. Going back to dimensional variables yields u/U∞ = h(x/L, (y/L) U∞ L/ν),
where the relationship between δ and L has been used (it involves the inverse of the squareroot
of the Reynolds number based on L). Note however that for the semi-infinite flat plate L is
an arbitrary distance so that there is no reason why the solution u(x, y) should depend upon
L. Therefore, one rather looks for a solution of the form u(x, y) = U∞ × h(η) where η is a
non-dimensional parameter built from x̄ and ỹ but which does not depend on L. The assumed
form for the solution u(x, y) is therefore :

u ỹ y yp
= h(η) = h( √ ) = h( q ) = h( Rex )
U∞ x̄ νx x
U∞
7.2. Particular solutions of the boundary layer equations 343

The x-velocity profile is assumed to display a so-called self-similar form : each profile for a given
section x is similar to another profile in an upstream or downstream section but with a different
scaling factor in the y-direction, assumed equal to νx/U∞ (see Fig. 7.18).

0.15

0.1
y

0.05

2 2.5 3 3.5 4 4.5 5 5.5 6 6.5

x

Figure 7.18: Boundary layer developing on a flat plate of length L = 20 m with ReL = 106 .
Detailed view of the x-velocity profiles at different successive x srections illustrating the
self-similar nature of the velocity field u(x, y).

For reasons which will soon become clear, it is actually better to assume :

u(x, y) = U∞ f 0 (η)

that is to introduce the derivative f 0 (η) of a function f (η).

For the incompressible flow under study, the mass balance reduces to the incompressibility
condition :
∂u ∂v
+ =0
∂x ∂y
This identity implies there exists a stream function ψ(x, y) such that :

∂ψ ∂ψ
u(x, y) = , v(x, y) = −
∂y ∂x

It is straightforward to check the incompressibility condition is indeed satisfied when u and v

are computed in that way. From the definition of the stream function ψ :
ˆ ˆ
ψ = udy = U∞ f 0 (η)dy

√ x
Since η = (y/x) × Rex , it follows that dy = √ dη. Applying this change of variable yields :
Rex
ˆ ˆ ˆ
x
f (η)dη = U∞ νx f 0 (η)dη
0
p
ψ= Ux dy = U∞ √
Rex
 344 Laminar boundary layer

hence
p
ψ(x, y) = U∞ νxf (η)

The choice of using f 0 (η) to define u(x, y) appears clearly : it allows to draw a direct link be-
tween ψ and f (η) (otherwise it would be necessary to introduce the primitive of the function
which would turn the problem into an integro-differential equation, more complex to solve).

The expression of v(x, y) as a function of f (η) and its derivative follows from the definition :

∂ψ ∂ hp i
v(x, y) = − =− U∞ νxf (η)
∂x ∂x
In order to compute v(x, y) one must take into account the fact that η depends on x (and also
on y of course but x and y are independent variables) :
 1
U∞ 2
η(x, y) = y
νx

It follows that :
  1
 ∂η U∞ 2

 = = η/y
 ∂y
 νx
 r
 ∂η = −y U∞ 1 = − η



∂x ν 2x 2x
Taking this into account for computing v(x, y) yields :
√
p ∂( x) p ∂f (η)
v(x, y) = − U∞ ν f (η) − U∞ ν x ×
∂x | ∂x
{z }
df (η) ∂η
= ×
dη ∂x

or else : r
1 U∞ ν p η
v(x, y) = − f (η) + U∞ ν x × × f 0 (η)
2 x 2x
q
1 U∞ ν
Factoring 2 x in this expression yields :

r
1 U∞ ν
v(x, y) = (ηf 0 (η) − f (η))
2 x

This expression for v(x, y), along with u(x, y) = U∞ f 0 (η), can be injected into the x-momentum
equation of the Prandtl equations :

∂u ∂u ∂2u
u +v =ν 2
∂x ∂y ∂y

in order to turn this partial differential equation (PDE) into an ordinary differential equation
(ODE) for the function f (η) and its derivatives.
7.2. Particular solutions of the boundary layer equations 345

The various partial derivatives of u which appear in the momentum equation can be expressed
using f (η) and its successive derivatives :
∂u ∂η


 = U∞ f 00 (η)



 ∂x ∂x


 ∂u

∂η
= U∞ f 00 (η)

 ∂y ∂y



2
 ∂ u = U f 000 (η)( ∂η )2



 ∞
∂y 2 ∂y
Injecting these expressions into the x-momentum equation yields :
r
∂η 1 U∞ ν ∂η ∂η
U∞ f 0 (η)U∞ f 00 (η) + (ηf 0 (η) − f (η))U∞ f 00 (η) = ν U∞ f 000 (η)( )2
∂x 2 x ∂y ∂y
Expliciting finally the partial derivatives of η with respect to x and y leads to :
r
2 0 η 1 U∞ ν η η2
−U∞ f (η) f 00 (η) + (ηf 0 (η) − f (η)) × U∞ f 00 (η) = ν U∞ f 000 (η) 2
2x 2 x y y
q
Since η/y = Uν∞ x it also comes :

η 0 2
1 U∞ U2
2
−U∞ f (η) f 00 (η) + (ηf 0 (η)f 00 (η) − f (η)f 00 (η)) = ∞ f 000 (η)
2x 2 x x
so that the successive simplification of the first term and the second term followed by the
2 /x yields an equation depending only on the similarity parameter η :
simplification of U∞

1
f 000 (η) + f (η)f 00 (η) = 0
2

This third-order ODE must be completed with 3 boundary conditions before being numerically
solved. These 3 conditions can be derived from the wall boundary conditions and the far-field
condition. Since u(x, y = 0) = 0 it follows that U∞ f 0 (η = 0) = 0 or else f 0 (0) = 0. Since
v(x, y = 0) = 0 it follows that [ηf 0 (η) − f (η)]η=0 = 0 or else f (0) = 0. Finally, far from the
plate, when y → ∞ or η → ∞, we have u → U∞ so that f 0 (η) → 1.

Finding the function f (η) requires to solve the following problem, known as Blasius equation :
 000

 f (η) + 21 f (η)f 00 (η) = 0
 f (0) = 0

f 0 (0) = 0 (7.27)

 lim f 0 (η) = 1


η→∞

In order to solve this third-order ODE, it is transformed into a system of first-order ODE which
will not be detailed here. A shooting method is eventually used to account for the fact one
of the boundary conditions applies at infinity. In the numerical simulation, η will naturally
keep a finite value but its upper value will be such that η >> 1. The numerical integration is
346 Laminar boundary layer

η f (η) f 0 (η) f 00 (η)

0. 0 0 0.33206
0.2 0.00664 0.06641 0.33199
0.4 0.02656 0.13277 0.33147
0.6 0.05974 0.19894 0.33008
0.8 0.10611 0.26471 0.32739
1.0 0.16557 0.32979 0.32301
1.2 0.23795 0.39378 0.31659
1.4 0.32298 0.45627 0.30787
1.6 0.42032 0.51676 0.29667
1.8 0.52952 0.57477 0.28293
2.0 0.65003 0.62977 0.26675
2.2 0.78120 0.68132 0.24835
2.4 0.92230 0.72899 0.22809
2.6 1.07252 0.77246 0.20646
2.8 1.23099 0.81152 0.18401
3.0 1.39682 0.84605 0.16136
3.2 1.56911 0.87609 0.13913
3.4 1.74696 0.90177 0.11788
3.6 1.92954 0.92333 0.09809
3.8 2.11605 0.94112 0.08013
4.0 2.30576 0.95552 0.06424
4.2 2.49806 0.96696 0.05052
4.4 2.69238 0.97587 0.03897
4.6 2.88826 0.98269 0.02948
4.8 3.08534 0.98779 0.02187
5.0 3.28329 0.99155 0.01591
6.0 4.27964 0.99898 0.00240
7.0 5.27926 0.99992 0.00022
8.0 6.27923 1.00000 0.00001

Table 7.1: Table of values for the Blasius solution.

7.2. Particular solutions of the boundary layer equations 347

typically performed using a Runge-Kutta method and yields a set of discrete values describing
the evolution of f (η), f 0 (η), f 00 (η) when η varies from 0 (falt plate wall) to ”infinity” (large value
with respect to unity) (see Table 7.1).

The values of Table 7.1 define the Blasius solution for the laminar boundary layer developing over
a semi-infinite flat plate with zero outside (inviscid) pressure gradient (dpe /dx = 0). Plotting
f 0 (η) eas a function of η yields the normalized or non-dimensional x-velocity profile, u/U∞ .
The agreement between this theoretical prediction and experimental measurements is excellent
as can be observed in Fig. 7.19. In the experiment, water flows at 9 cm/s so that the local
Reynolds number in the section is 500 (δ ∗ = 5 mm). A slender colloidal cloud is produced by
the fine tellurium wire perpendicular to the plate at the left of the photograph when this wire
is subjected to an electric impulse during a few milliseconds. The cloud is advected along the
streamlines and captured by the photograph (taken by F.X. Wortmann [3]).

Ecoulement de fluide parfait

8

4
η

Paroi (condition d’adherence)

0
0 0.25 0.5 0.75 1
u/ue = f’(η)

Figure 7.19: Left : theoretical Blasius solution for the laminar boundary layer over a
semi-infinite flat plate with zero external pressure gradient. Right : experimental visualization
of the Blasius boundary layer profile on a flat plate.
 348 Laminar boundary layer

Exploitation of the Blasius solution

Looking at Fig.7.19 and at the table of discrete values 7.1, it is observed the boundary layer
thickness defined as δ99% , that is such that f 0 (η) goes above 0.99, satisfies :

η(x, y = δ99% ) ≈ 5

so that, taking into account the definition of the similarity parameter η :

 1
U∞ 2
δ99% ≈5
νx
The evolution of the boundary layer thickness for the Blasius solution is therefore such that :
 1
νx 2
δ99% (x) = 5
U∞

This parabolic evolution of δ(x) with x had already been qualitatively predicted but the Blasius
solution provides a quantitative result. This law of evolution can also be rewritten using the
local Reynolds number Rex = U∞ x/ν :

δ99% (x) 5
=√
x Rex

If the plate is not semi-inifinite but is of total length L, the boundary layer thickness at the end
of the plate (x = L) is such that :
δ99% (L) 5
=√
L ReL
For the case ReL = 106 previously reported (even though such a large value of the Reynolds
number corresponds no longer to a laminar boundary layer but rather to a turbulent one),
δ99% (L)
= 0.005. Since L = 20 m, it follows that δ(L) = 0.1 m or 10 cm.
L
For the Blasius boundary layer, such that ue (x) = U∞ , the boundary layer thickness is such
that : ˆ ∞ 
∗ u(x, y)
δ (x) = 1− dy
0 U∞
Introducing the expression of Ux /U∞ as a function of the similarity parameter η yields :
ˆ ∞
∗
1 − f 0 (η) dy


δ (x) =
0

The integration is performed over η rather than over y using the identity :
 1
∂y νx 2
dy = dη = dη
∂η U∞
It follows that :
 1 ˆ ∞  1
∗ νx 2
0

νx 2
δ (x) = 1 − f (η) dη = [(η − f (η))η→∞ − (η − f (η))η=0 ]
U∞ 0 U∞
7.2. Particular solutions of the boundary layer equations 349

or else, using the boundary condition f (η = 0) = 0 and introducing the local Reynolds number :

δ ∗ (x) 1
=√ (η − f (η))η→∞
x Rex

From the values reported in Table 7.1, it follows that (η − f (η))η→∞ ≈ 1.72 so that :

δ ∗ (x) 1.72
=√
x Rex

The momentum thickness is such that :

ˆ ∞  1 ˆ ∞
0 0 νx 2
f 0 (η) 1 − f 0 (η) dη



θ(x) = f (η) 1 − f (η) dy =
0 U∞ 0

which yields after an immediate calculation :

θ(x) 0.664
=√
x Rex

This result can be used for instance to estimate the location of transition for the Blasius boundary
layer using the Michel criterion (7.24). Indeed :

U∞ θ U∞ x θ 0.664 p
Reθ = = × = Rex × √ = 0.664 Rex
ν ν x Rex

Inserting into (7.24) yields : p 0.4

0.664 (Rex )tr = 2.9(Rex )tr
or else (Rex )0.1
tr = 4.3675 so that (Rex )tr ≈ 2525000.
If the Cebeci-Smith criterion is used :
 
0.5 22400
0.664(Rex )tr = 1.174 1 + (Rex )0.46
tr
(Rex )tr

so that X = (Rex )tr is solution of the non-linear equation :

X 1.04 − 1.768X − 39603.2 = 0

hence (Rex )tr ≈ 2025000.

The shape factor of the Blasius boundary layer is such that :

δ ∗ (x) 1.72
H(x) = = ≈ 2.59
θ(x) 0.664

The skin friction coefficient along the flat plate is such that :
 1
τw 2ν ∂u 2ν 00 ∂η 2ν 00 U∞ 2
Cf = 1 2 = 2 ( )y=0 = (f (η))η=0 ( )y=0 = f (0)
2 ρU ∞ U∞ ∂y U ∞ ∂y U ∞ νx
 350 Laminar boundary layer

or else
2f 00 (0) 0.664
Cf (x) = √ =√ (7.28)
Rex Rex
The drag coefficient of the upper surface of the flat plate is given by (7.16) :
ˆ
D 1 L
CD = 1 2 = Cf (x) dx
2 ρue L
L 0

Since Cf (x) is more naturally expressed as a function of Rex , it is customary to compute the
above integral with the following change of variable :
U∞
dRex = dx
ν
so that :
ˆ ReL ˆ ReL
ν 1
CD = Cf (Rex ) dRex = Cf (Rex ) dRex
U∞ L 0 ReL 0

Inserting (7.28) into this formula yields :

ˆ ReL
1 0.664 0.664  ReL
CD = 0.5
dRex = 2Re0.5
x 0
ReL 0 Rex ReL
hence the final expression of the drag coefficient for the upper surface of the flat plate :

1.328
CD = √
ReL

For the boundary layer displayed in Fig. 7.17, ReL = 106 so that the friction drag is equal to
CD = 0.001328. The drag D (per unit width or unit span) is computed as :
1 2
D = ρU∞ L × CD (ReL ) = 0.5 × 1.225 × 0.52 × 20 × 0.001328 ≈ 4 × 10−3 N
2
This very low value is a direct consequence of the very low upstream dynamic pressure for the
flow under study : q∞ = 12 ρU∞2 = 0.153 P a.

The same Reynolds number ReL = 106 can be achieved with the same fluid (air with ν =
10−5 m2 /s) for a flow at U∞ = 200 m/s on a plate of length L = 0.05 m. In that case, the
upstream dynamic pressure is much larger : q∞ = 21 ρU∞ 2 = 24500 P a. The resulting friction

drag (by unit span or width of the plate) is given by :

D = q∞ L × CD (ReL ) = 1.63 N

7.2.2 Falkner-Skan solution

Another self-similar solution of the boundary layer equations can be obtained for an outer
inviscid velocity distribution ue (x) a bit more general than the constant distribution of the
Blasius solution. Falkner and Skan proved in 1930 (Blasius solution was obtained in 1905) a
self-similar solution exists for problem (7.12) in the case where ue (x) is of the form :
 x m
ue (x) = U0 or ue (x) = Cxm
L
7.2. Particular solutions of the boundary layer equations 351

Using the theory of potential flows (not covered in this course), it can be proved such a velocity
distribution corresponds to the incompressible inviscid flow over a dihedron with dihedral angle
βπ such that m = β/(2 − β) (see Fig.7.20).

β π/2

β π/2

2mπ
Figure 7.20: Potential flow over a dihedron with dihedral angle βπ = .
m+1

Since a streamline is equivalent to a wall boundary in the case of an inviscid flow (because of the
slip boundary condition applying along a wall boundary), the flow over the dihedron can also be
interpreted as described in Fig.7.21. The velocity distribution ue (x) = Cxm with m > 0 (hence
β > 0) corresponds to the accelerating flow over a concave corner of angle mπ/(m + 1) with
respect to the horizontal incoming flow direction. The velocity distribution ue (x) = Cxm with
m < 0 (but larger than −1, hence β < 0) corresponds to the decelerating flow over a convex
corner of angle mπ/(m + 1) with respect to the horizontal incoming flow direction. The concave
corner / accelerating flow configuration corresponds to a decreasing outer pressure distribution,
hence dpe /dx < 0 that is a favourable pressure gradient. The convex corner / decelerating flow
configuration corresponds to an increasing outer pressure distribution, hence dpe /dx > 0 that
is an adverse pressure gradient dpe /dx > 0. Note also the two particular cases : m = 0 that is
β = 0 corresponds to the flat plate with zero outer pressure gradient (Blasius solution); m = 1
corresponds to a dihedral angle π (β = 1) that is a flat plate positioned orthogonally with respect
to the incoming flow, a configuration corresponding to the flow in the immediate vicinity of a
stagnation point. The evolution of the flow velocity ue (x) is displayed in Fig.7.22; note that
x is counted / measured from the corner. The decreasing ue (x) distribution for m = −0.2
corresponds to a convex corner with angle mπ/(m + 1) = −π/4 = −0.7854 rad that is −45◦ or
β = 2m/(m + 1) = −0.5. The value m = −0.0905 corresponds to a decreasing inviscid velocity
ue (x) over a convex corner with angle −0.1π or −17.91◦ below the horizontal (incoming flow)
direction. The value m = 0.1 corresponds to an increasing inviscid velocity ue (x) over a concave
corner with angle 0.0909π that is 16.4◦ above the horizontal (incoming flow) direction.

As done for the Blasius solution the non-dimensional velocity profile u(x, y)/ue (x) is√searched
in the form of a function of the non-dimensional variables x̄ = x/L and ỹ = (y/L) Re with
Re = ue (L)L/ν, that is Re = CLm+1 /ν. Since the dihedron or the corners are semi-infinite
(they extend from x = 0 at the corner to a downstream section at an infinite distance from
the corner) there is no reason for the solution to depend on L. Consequently, u(x, y)/ue (x) is
assumed to be a function of a single non-dimensional (similarity) parameter η built from x̄ and
 352 Laminar boundary layer

Figure 7.21: Inviscid flow configurations for the Falkner-Skan solution. Top : flow over a
concave corner with m > 0. Bottom : flow over a convex corner with m < 0.

2.5
m=-0.2
m=-0.0905
m=0
m=+0.1
m=1
2

1.5
ue(x)

0.5

0
0 0.5 1 1.5 2
x

Figure 7.22: Velocity distribution ue (x) = xm for several values of m.

7.2. Particular solutions of the boundary layer equations 353

ỹ and not dependent on L :

√
r
u(x, y) Cxm−1
= h(η) = h(ỹ x̄m−1 ) = h(y ) = f 0 (η)
ue (x) ν

Following a line of thought perfectly similar to the one detailed for the Blasius solution, the
Prandtl equations can be turned into the following third-order ODE completed with wall bound-
ary conditions and matching condition with the inviscid outer velocity field :
 m+1
 f 000 (η) + f (η)f 00 (η) + m(1 − f 0 (η)2 ) = 0
2





f (0) = 0 , f 0 (0) = 0 (7.29)





 0
f (∞) = 1

Note that for m = 0 the Blasius equation (7.27) is recovered. The Falkner-Skan equation (7.29)
can be numerically solved using the same shooting method as the one used for the Blasius
equation. The various velocity profiles displayed in Fig.7.23 are obtained depending on the
value of the angle β or, equivalently, the value of the exponent m for the velocity distribution.

courant de retour
β < −0.199 decollement
β=−0.199

β<0
et decroissant

β=0

β>0
et croissant u/u e
0 1
Figure 7.23: Velocity profile in the boundary layer developing over a convex or concave corner
of angle βπ/2 = mπ/(m + 1) corresponding to an outer inviscid velocity ue (x) = Cxm .
Falkner-Skan solution. ”croissant” = ”increasing”, ”décroissant” = ”decreasing”,
”décollement” = ”separation” or ”detachment”, ”courant de retour” = ”backflow”.

The computed velocity profile for β = 0 (m = 0) corresponds to the Blasius solution ( dp e

dx = 0
hence ue = constant = U∞ ). For β > 0 (m > 0), the pressure gradient dpe /dx is favourable
 354 Laminar boundary layer

since it tends to accelerate the flow in the boundary layer in the immediate vicinity of the wall
( dp e
dx < 0) : the boundary layer gets thinner ahd the skin friction increases with a fuller velocity
profile near the wall. For β < 0, the outer inviscid pressure gradient is adverse ( dp e
dx > 0) : it
contributes to slow down the flow in the immediate vicinity of the wall, eventually leading to
the detachment or separation of the boundary layer from the wall. The limit attached state
corresponds
 to  a vertical tangent at the wall associated with a zero value of the wall shear stress
∂u
τw = µ or equivalenty a zero skin friction coefficient Cf . This flow configuration is
∂y y=0
reached when β = −0.199, wich corresponds to the flow over a convex corner inclined by an
angle −17.9◦ below the horizontal incoming flow direction. Note that since a wall is equivalent
to a streamline for an inviscid flow, the convex corner flow configuration displayed at the bottom
of Fig.7.21 corresponds also to the flow over a flat plate inclined by an angle −17.9◦ below the
horizontal incoming flow direction (see Fig. 7.24). Beyond this critical value, the boundary
layer is detached from the wall and backflow (u < 0) appears in the flow. The detached flow is
no longer governed by the boundary layer equations : the numerical simulation of such a flow
configuration requires to solve the full Navier-Stokes equations. Note the solution obtained for
β = 1 can be used as an initial velocity profile when integrating the boundary layer equations
along an airfoil from the leading edge (stagnation region where the β = 1 solution approximates
the initial velocity profile) down to the trailing edge.

7.3 Integral methods for boundary layer analysis

7.3.1 von Kármán integral equation
One of the key objectives of a quantitative analysis of the boundary layer is to compute the wall
shear stress τw so as to estimate the friction drag and also to predict a possible detachment of
the boundary layer. This analysis is performed for a known outer velocity distribution ue (x),
issued for instance from a preliminary potential analysis of the flow. When this external (invis-
cid) velocity distribution is arbitrary (and no longer corresponds to the specific form of ue seen
for the Blasius solution - ue (x) = constant - or the Falkner-Skan solution - ue (x) = Cxm -), the
search for a self-similar solution generally fails. An alternative solution strategy is to settle for
an approximate solution of the Prandtl equations which could be obtained whatever the outer
velocity distribution ue (x). The integral von Kármán equation is precisely such a tool.

Let us consider the Prandt equations written for a general outer velocity field ue (x) and with
the diffusive term in the momentum equation expressed using the shear stress τ given by the
∂u
simplified expression τ = µ in the boundary layer because of the high-Reynolds number
∂y
nature of the flow :
∂u ∂v


 + =0 (i)
 ∂x ∂y


 u ∂u + v ∂u = ue due + 1 ∂τ


 (ii)
∂x ∂y dx ρ ∂y
Equations (i) and (ii) are completed with the wall boundary conditions (u(x, y = 0) = v(x, y =
0) = 0), the matching condition with the outer inviscid velocity field and farfield boundary
conditions expressing the fact that far from the wall the shear stress is zero (the flow becomes
7.3. Integral methods for boundary layer analysis 355

Figure 7.24: Top : interpretation of the Falkner-Skan solution with m < 0 as the flow over a
flat plate at incidence. Theoretical prediction : separated flow for β < −0.199 corresponding to
an incidence angle greater than 17.9◦ (below the incoming flow direction). Bottom :
experimental observation of laminar boundary layer separation over a flat plate at 20◦ of
incidence (photograph by Werlé (ONERA, 1974) taken from [3]).
 356 Laminar boundary layer

uniform) and the y-velocity component tends to zero :

τ (y → ∞) = 0 , v(y → ∞) = 0.
Let us form the combination (ii)+(u−ue )(i) from the continuity equation (i) and the momentum
equation (ii) and let us integrate the resulting equation from the wall y = 0 to y → ∞. Taking
into account the boundary conditions at the wall and when y → ∞ yields :
ˆ ∞ ˆ ∞
1 due ∂
τw = · (ue − u)dy + ( u(ue − u)dy)
ρ dx 0 ∂x 0
where τw denotes the shear stress at the wall (y = 0). Introducing the displacement thickness
and the momentum thickness, the above equation can be recast as :
τw 1 due ∗ dθ
2
= (δ + 2θ) +
ρue ue dx dx
τw
Finally, using the shape factor H, ratio of δ ∗ to θ, and the skin friction coefficient Cf = 1 2
,
2 ρue
the following standard form of the von Kármán integral equation can be obtained :
Cf dθ θ due
= + (H + 2) (7.30)
2 dx ue dx
In this relationship only the outer inviscid velocity ue (x) is known; the momentum thickness θ,
the shape factor H and the skin friction distribution Cf are to be determined.
Equation (7.30) is the starting point of an approximate solution technique for the boundary
layer equations. The key idea of this solution strategy is that since the exact solution u(x, y)
of the Prandlt equations satisfies (7.30), this equation cna be used as a ”filter” to obtain an
approximate expression for Cf from an assumed velocity distribution u(x, y)/ue which can be a
coarse approximation of the real (unknown) velocity profile in the boundary layer. The typical
methodology can be decomposed into the following steps :
-• a form is postulated for the boundary layer x-velocity profile :
 
u(x, y) y
=f
ue (x) δ(x)

The outer velocity ue (x) is known. The function f is defined in such a way it satisfies
some of the boundary conditions associated with the boundary layer. The boundary layer
thickness δ(x) is unknown.
-• the displacement thickness and the momentum thickness are computed from their respec-
tive definitions. For the displacement thickness, we have :
ˆ δ(x)
∗ u(x, y)
δ (x) = (1 − )dy
0 ue (x)
y
so that, using the change of variable : ξ = we can also write :
δ(x)
ˆ 1
∗
δ (x) = δ(x) (1 − f (ξ))dξ = k1 δ(x)
|0 {z }
=k1
7.3. Integral methods for boundary layer analysis 357

Similarly, the momentum thickness is such that :

ˆ 1
θ(x) = δ(x) f (ξ)(1 − f (ξ))dξ = k2 δ(x)
|0 {z }
=k2

By definition, the shape factor is such that :

δ ∗ (x) k1
H= =
θ(x) k2
The wall shear stress is such that :
 
∂u µue (x) 0
τw (x) = µ = f (ξ)|w
∂y w δ(x)

so that :
µue (x) 0
τw (x) = f (0)
δ(x)
where f 0 (0) can be computed since f has an assumed form. Consequently the skin-friction
coefficient reads :
Cf τw ν 0
= 2 = f (0)
2 ρue ue δ

-• Inserting these expressions for Cf , θ and H into the von Kármán integral equation yields
the following differential equation for the unknown function δ(x) :

νf 0 (0) dδ δ 2 due
= k2 δ + (2k2 + k1 )
ue (x) dx ue (x) dx

This is more precisely a first-order ODE for the unknown quantity δ 2 (x) :

k2 d(δ 2 ) 1 due 2 νf 0 (0)

+ (2k2 + k1 ) δ (x) =
2 dx ue (x) dx ue (x)

-• solving this ODE for a known outer velocity distribution ue (x) and an assumed function
f (ξ) yields δ 2 (x) hence δ(x) from which the skin friction coefficient Cf (x) can be computed.

Examples of application will be proposed as Exercises in Section 7.4.

7.3.2 Polhausen’s method

It has been mentioned in the previous description of the von Kármán equation methodology that
the function f chosen to approximate the velocity profile in the boundary layer must satisfy a
number of physical boundary conditions.
Let us assume f is defined as a polynomial function of the non-dimensional variable ξ(x, y) =
y/δ(x). The number of boundary conditions which can be satisfied by f (ξ) depends on the order
of this polynomial function. For instance, let us assume f (ξ) is defined as a polynomial function
of order 1, that is a linear function of ξ : f (ξ) = aξ + b. The boundary condition at the wall
 358 Laminar boundary layer

states that u = 0 for y = 0, that is f (0) = 0. The boudary condition at the upper boundary
of the boundary layer states u = ue for y = δ that is f (1) = 1. Consequently, a and b are such
that a = 1 and b = 0 hence f (ξ) = ξ.

More generally, the x-velocity u(x, y) in the boundary layer and its successive derivatives must
satisfy the following boundary conditions :

-• u(x, y = 0) = 0 expresses the no-slip condition on a fixed wall. This translates into :
f (0) = 0 .

-• u(x, y = δ(x)) = ue (x) expresses the matching between the flow in the boundary layer and
the outer inviscid flow. This translates into : f (1) = 1 .

∂u(x, y = δ(x))
-• = 0 expresses the fact there is no friction outside the boundary layer. It
∂y
can also be written as : f 0 (1) = 0 .

-• for a steady flow, the momentum equation in the boundary layer reads :

∂u ∂u due ∂2u
u +v = ue +ν 2
∂x ∂y dx ∂y
When written at the wall and taking into account the wall boundary conditions, we obtain
a boundary condition on the second derivative of u with respect to y at the wall :

∂ 2 u(x, 0) ue due
2
=−
∂y ν dx

If the boundary layer develops with a zero outer pressure gradient (case of the Blasius
∂2u
problem), the boundary condition reduces to : (y = 0) = 0. If ue 6= constant, this
∂y 2
condition on the x-velocity second derivative in the normal direction to the wall expresses
the influence of the outer pressure gradient (be it adverse or favourable) on the boundary
∂2u ue
layer development. Since 2
= 2 f 00 (ξ), the condition can also be written as :
∂y δ

δ 2 due
f 00 (0) = − = −Λ
ν dx

where Λ is a parameter expressing the influence of the outer pressure gradient on the
velocity profile in the boundary layer.

-• Writing the momentum equation at y = δ and using the boundary conditions u(x, δ) = ue
∂u(x, δ)
and = 0, another last boundary condition can be written :
∂y

∂ 2 u(x, y = δ)
=0
∂y 2

which translates into f 00 (1) = 0 .

7.3. Integral methods for boundary layer analysis 359

These 5 boundary conditions can now be exploited to determine the coefficients of a polynomial
of order 4, known as the Polhausen polynomail (initial proposed by Polhausen in 1921) :

u(x, y) Λ
= f (ξ) = [2ξ − 2ξ 3 + ξ 4 ] + ξ(1 − ξ)3 (7.31)
ue (x) 6

The first component of the polynomial, between brackets, corresponds to the case of a flow
without external pressure gradient. The second component, proportional to the parameter Λ,
represents a correction of the first component, which accounts for the existence of an external
pressure gradient.
The general methodology described in the previous section can be applied with the choice (7.31)
for f (ξ) and an approximation of the skin friction coefficient Cf (x) can thus be derived whatever
the external velocity distribution ue (x). This approach can be rather accurate for some flow
configurations (for instance the simple case of a flat plate, see section 7.4). However, it does not
necessarily provide an accurate prediction of the skin friction in the general case (that is for any
outer velocity distribution ue (x)); in particular it usually fails to predict the correct location of
the boundary layer separation in the case of an adverse pressure gradient (Λ < 0).

7.3.3 Walz-Thwaites’ method

This solution method had been developed in several stages :

-• in 1940, Holstein and Bohlen propose to retain the following non-dimensional parameter
to characterize the development of boundary layers :

θ2 due
λ= . (7.32)
ν dx

This parameter is relevant because it makes it possible to correlate the key non-dimensional
quantities apparearing in the von Kármán equation, assuming :

τw θ
= S(λ) and H = H(λ)
µue

where S(λ) and H(λ) are functions to be determined. Under these assumptions, the von
Kármán equation can be recast in the form :

d λ
ue (x) ( 0 ) = 2[S(λ) − λ(2 + H(λ))] = F (λ) (7.33)
dx ue (x)

-• In 1941, Walz shows that if F (λ) is linear then (7.33) can be analytically integrated and
θ2
yields : = φ(ue ).
ν
-• In 1949, after a careful analysis of a large number of experimental and theoretical results,
Thwaites shows that the following simple (linear) correlation holds for a wide range of
boundary layers :
F (λ) ≈ 0.45 − 6.0λ
 360 Laminar boundary layer

λ H(λ) S(λ) λ H(λ) S(λ) λ H(λ) S(λ)

+0.25 2.0 0.500 -0.008 2.64 0.208 -0.072 3.23 0.083
0.20 2.07 0.463 -0.016 2.67 0.195 -0.074 3.30 0.076
0.14 2.18 0.404 -0.024 2.71 0.182 -0.076 3.38 0.067
0.12 2.23 0.382 -0.032 2.75 0.168 -0.078 3.47 0.055
0.10 2.28 0.359 -0.040 2.81 0.153 -0.079 3.52 0.049
0.080 2.34 0.333 -0.048 2.87 0.138 -0.080 3.58 0.039
0.075 2.35 0.326 -0.052 2.90 0.130 -0.0804 3.59 0.035
0.064 2.39 0.313 -0.056 2.94 0.122 -0.0808 3.61 0.030
0.048 2.44 0.291 -0.060 2.99 0.113 -0.0812 3.63 0.024
0.032 2.49 0.268 -0.064 3.05 0.104 -0.0816 3.66 0.016
0.016 2.55 0.244 -0.068 3.13 0.094 -0.0818 3.69 0.011
0.0 2.61 0.220 -0.070 3.17 0.089 -0.0820 3.70 0.

Table 7.2: Walz-Thwaites correlations for H(λ) and S(λ).

Using the result previously established by Walz, Thwaites derived the following formula
giving the momentum thickness θ for a laminar boundary layer with a known external
velocity distribution ue (x) :
ˆ x
2 0.45ν
θ (x) = 6 u5e (t)dt (7.34)
ue (x) 0

where t is a mute integration variable.

θ2 due
-• Once θ2 (x) is computed from ue , the parameter λ = can be also computed.
ν dx
Thwaites also provided correlation tables for H(λ) and S(λ) (see Table 7.2).
It is therefore possible to compute the displacement thickness from δ ∗ = H(λ)θ and the
µue
wall shear layer from τw = S(λ). Note also that the value λ = 0 corresponds to the
θ
position of the maximum external velocity (due /dx = 0). The value λ = λsep = −0.082
corresponds to the separation location since S(λsep ) = 0.
Note the Walz-Thwaites method does not allow to determine the velocity profile u(x, y)
in the boundary layer. It is nonetheless one of the best simple methods to approximately
compute the properties of a boundary layer with known external inviscid velocity distri-
bution ue (x). It has been of course progressively replaced with the numerical simulation
of the full Navier-Stokes equations but remains nonetheless an interesting tool.

An example of application of the method is proposed in Section 7.4.

7.3.4 Simple integral method

A simple way to use von Kármán integral equation consists in assuming the shape factor H
remains roughly constant in the flow, which is rather well verified in practice as long as the
external pressure gradient is not too strong. This constant value is taken equal to the value
computed for a flat plate boundary layer while the skin friction coefficient Cf is assumed to
7.3. Integral methods for boundary layer analysis 361

follow a flat plate law of the form :

Cf b
=
2 Remθ

Using the Blasius solution yields m = 1 and b = 0.2205 in the above expression as well as
H = 2.591. Injecting into the von Kármán equation and multiplying this equation by (m +
(m+1)(H+2)
1)θm ue , the following ODE is obtained :

(m+1)(H+2)
d ue
[θu(H+2)
e ]m+1 = (m + 1)b
dx (ue /ν)m

It can integrated in a straightforward way to yield the following equation of evolution for the
momentum thickness :

ˆ x1
(H+2) (m+1) (H+2) (m+1)[ue (x)](H+2)(m+1)
[ue (x1 ) θ(x1 )] = [ue (x0 ) θ(x0 )] + (m + 1)b dx
x0 (ue (x)/ν)m
(7.35)
The formula (7.35) provides an approximate expression for the momentum thickness evolution
as a function of the external inviscid velocity distribution ue (x) and the value of θ(x) in an
initial section x0 .

In the case of a laminar boundary layer developing over a flat plate, we can choose x0 = 0 the
leading edge of the plate where the momentum thickness is equal to zero. In the case where
ue = constant, the identity (7.35) yields :
ˆ x
m+1 1
θ(x) = b(m + 1) dt
0 (ue /ν)m

so that
θ(x) [b(m + 1)]1/(m+1)
= ue x m/(m+1)
x ( )
ν
The Blasius formula is recovered after a numerical application for b and m :

θ(x) 0.664
=√
x Rex

A few analytical or semi-analytical methods have been presented, which allow to predict in a
quantitative way the key properties (friction, separation) of a laminar boundary layer developing
in an incompressible flow over a wall boundary. Note however many external or internal flows in
aerodynamics are laminar only on a fraction of the wall boundary extent (near the leading edge
of the airfoil for instance) so that a correct estimate both of the skin friction and the separation
location must take into account the transition of the boundary layer to turbulence. This topic
will be covered in the next Lecture.
 362 Laminar boundary layer

couche limite se développant

sur la paroi de la soufflerie

ue(0) ue(L) h

x=0 x=L

Figure 7.25: Uniform flow between the upper and lower walls in the midsection of the wind
tunnel. The thickness of the boundary layers displayed in the picture is not to scale.

7.4 Exercises and problems

7.4.1 Exercise #1 : boundary layer developing along the walls of a wind

tunnel

Questions

We wish to improve the operation of a wind tunnel by taking into account the development of
a boundary layer on the walls of the wind tunnel. We work in the midplane of the wind tunnel,
with height h = 1 m and length L = 6 m. The inflow velocity is ue = 70 m / s. In order to
ensure the good quality of the experiments performed in the wind tunnel, we wish to establish
a flow as uniform as possible on the greatest streamwise extent of the wind tunnel. Velocity
measurements indicate the flow is slightly accelerated in the outlet section with respect to the
inlet section.

1) a) • Explain qualitatively why the development of a boundary layer on the upper and lower
walls of the wind tunnel midsection can cause the measured flow acceleration between the inlet
and the outlet of the wind tunnel.
The flow of air is assumed uniform in a section of the wind tunnel, between the lower and upper
boundary layers (see Figure 7.25).

• Assuming the flow takes place with a constant density ρ, express the volume flowrate conser-
vation between the inlet section and outlet section of the wind tunnel and obtain a relationship
between the inlet and outlet velocities ue (0), ue (L), the displacement thickness in the outlet
section δ ∗ (L) and the height of the wind tunnel.

b) • Explain why a small divergence of the wind tunnel walls (see Fig.7.26) can ensure an almost
constant velocity ue along the whole length of the wind tunnel.

It is now assumed the evolution of the displacement thickness and momentum thickness is well
7.4. Exercises and problems 363

x=0 x=L

Figure 7.26: Wind tunnel with inclined walls. The angle given to the walls of the wind tunnel
is not to scale.

approximated by the Blasius relationships valid for a flat plate with no outer pressure gradient :
1.7208x 0.664x
δ ∗ (x) = √ , θ(x) = √
Rex Rex
ue x
with Rex = . The kinematic viscosity of air in the wind tunnel is ν = 1.44 × 10−5 m2 /s).
ν
• Compute the angle α which should be applied to the wind tunnel walls to ensure a constant
velocity flow over the wind tunnel length.

2) In reality, the initially laminar boundary layer developing over the wind tunnel walls becomes
turbulent downstream of the transition location xt .

a) The location of transition is estimated using the empirical formula of Cebeci and Smith :
22400
Reθ(xt ) = 1.174(1 + )Re0.46
xt
Rext
• Check the boundary layer developing over the wind tunnel walls is essentially turbulent.
Note : the plot displayed in Fig.7.27 can be used to analyze the proposed transisition criterion.

b) It is assumed the evolution of the turbulent boundary layer thickness over the wind tunnel
wall is well approximated by the formula :
0.37x
δ(x) = 1/5
Rex

• Assuming the velocity profile in the turbulent boundary layer is of the form (u/ue ) = (y/δ)1/7
(where y denotes the distance to the wall along which the boundary layer develops) , compute
the displacement thickness δ ∗ (x).

• Compute in that case the angle α yielding a constant velocity over the wind tunnel length.

Answers
 1) a) • The flow velocity is zero on the upper and lower walls of the wind tunnel. In the
boundary layer, the velocity u goes from 0 to ue (x) over a normal distance to the wall equal to
 364 Laminar boundary layer

1.09

1.08

1.07

1.06

1.05

1.04
f(z)

1.03

1.02

1.01

0.99

0.98

1E+06 2E+06 3E+06

z

1.768 39603.2
Figure 7.27: Graphical representation of the function f (z) = + .
z 0.04 z 1.04

δ(x), the thickness of the boundary layer. For the constant density flow under study, mass con-
servation reduces to volume flowrate conservation. With a developing boundary layer in which
the velocity is between 0 and ue (x), the external velocity ue must increase with x to balance
this velocity decrease near the walls.

• More precisely, mass conservation or flowrate conservation between the inlet and outlet sections
of the wind tunnel reads :
ˆ h
ρue (0)h = ρu(L, y)dy
ˆ0 δ(L) ˆ h−δ(L) ˆ h
= ρu(L, y)dy + ρue (L)dy + ρu(L, y)dy
0
ˆ δ(L) δ(L) h−δ(L)

= 2ρ u(L, y)dy + ρue (L)(h − 2δ(L))

0 ˆ δ(L)
= ρue (L)h − 2ρ (ue (L) − u(L, y))dy
0
= ρue (L)h − 2ρue (L)δ ∗ (L)

so that :
δ ∗ (L)
 
ue (0) = ue (L) 1 − 2 .
h

b) • If a small inclination is given to the wind tunnel walls, it is possible to conserve the volume
flowrate at the outlet through an increase of the wind tunnel section, without increasing the
velocity ue (L). This slight increase of the section can be adjusted so as to balance the change
7.4. Exercises and problems 365

of velocity induced by the boundary layer development. The flowrate conservation implies :
ˆ h0
ue (0)h = u(L, y)dy = ue (L)h0 − 2ue (L)δ ∗ (L)
0

so that :
h0 δ ∗ (L)
 
ue (0) = ue (L) −2
h h
h0 δ ∗ (L)
Choosing h0 such that −2 = 1 ensures ue (0) = ue (L). Consequently, α such that
h  h∗
∗

δ (L) −1 δ (L)
tan (α) = or α = tan ensures ue (L) = ue (0).
L L
δ ∗ (L) 1.7208
The law of evolution proposed for δ ∗ yields = √ . Since ν = 1.44 × 10−5 m2 /s,
L ReL
ue L
ue (0) = 70 m/s and L = 6 m, the Reynolds number ReL is such that ReL = = 29166667,
∗
ν
δ (L)
hence = 3.1863 × 104 and α = 1.826 × 10−2 deg (or 3.2 × 10−4 rad). This value of the
L
wall deviation angle is very small. Note also the assumption of a laminar boundary layer is not
realistic considering the very large value of ReL .

2) a) • Using the proposed law of evolution of the momentum thickness of the laminar boundary
layer, the Reynolds number based on this thickness is computed as :
ue θ(xt ) ue 0.664xt p
Reθ(xt ) = = × p = 0.664 Rext .
ν ν Rext
Denoting z = Rext , the proposed transition criterion reads :
22400 0.46
0.664z 0.5 = 1.174(1 + )z
z
or else, after a straightforward calculation :
1.768 39603.2
f (z) = + = 1.
z 0.04 z 1.05
It can be checked in Fig.7.27 that z = 2 × 106 is the root of this equation. It follows that
ue L xt xt
× = 2 × 106 with ReL ≈ 29.167 × 106 . This yields ≈ 0.07 or xt = 0.4 m. The
ν L L
transition of the boundary layer to turbulence occurs near the inlet section of the wind tunnel.

b) • In order to simplify the analysis, the boundary layer is assumed turbulent all over the wind
tunnel walls, with a thickness given by :
0.37x
δ(x) = 1/5
Rex
u y
and a velocity in the boundary layer assumed of the form : = ( )1/7 .
ue δ
The displacement thickness δ ∗ (x) is computed as follows :
ˆ δ ˆ 1
∗ u δ
δ = (1 − )dy = δ (1 − Y 1/7 )dY = .
0 ue 0 8
 366 Laminar boundary layer

0.04625x
so that δ ∗ (x) = 1/5
.
Rex
• Working in the same way as in 1) b) to determine the angle α to be applied to the wind
tunnel walls in order to ensure ue = constant yields :

δ ∗ (L)
 
α = tan−1 = 0.085 deg = 1.5 × 10−3 rad.
L

This angle is larger than the one obtained in 1) b) assuming a laminar boundary layer, which
makes sense since the turbulent boundary layer is thicker than the laminar boundary layer. 

7.4.2 Exercise # 2 : Polhausen polynomial and von Kármán equation

Questions

Let us consider a laminar boundary layer developing on a flat plate with zero external pressure
gradient.

1) The velocity profile in the boundary layer is assumed of the form :

 
u(x, y) y
=f = f (ξ) = a0 + a1 ξ + a2 ξ 2 + a3 ξ 3
ue (x) δ(x)

• Determine the coefficients of the third-order polynomial from the boundary conditions at the
plate wall and at the upper boundary of the boundary layer.

2) • Compute the displacement thickness δ ∗ , momentum thickness θ and shape factor H.

3) • Solve the von Kármán equation and derive an expression for the skin friction coefficient
Cf (x).

Answers

1)  The boundary conditions fulfilled by the function f (ξ) are :

-• u(x, y = 0) = 0 hence f (0) = 0. Using the assumed form for f (ξ) this condition yields
a0 = 0 hence f (ξ) = a1 ξ + a2 ξ 2 + a3 ξ 3 .

∂2u
 
-• = 0 from the momentum equation written at the wall, hence f 00 (0) = 0. Since
∂y 2 y=0
f 00 (ξ) = 2a2 + 6a3 ξ, this boundary condition yields a2 = 0 hence f (ξ) = a1 ξ + a3 ξ 3 .

-• at y = δ(x) that is ξ = 1, u = ue that is f (ξ = 1) = 1. This condition yields a1 + a3 = 1 .

∂u
-• at y = δ(x) , = 0 (no friction outside the boundary layer) hence f 0 (ξ = 1) = 0 that is
∂y
a1 + 3a3 = 0 .
7.4. Exercises and problems 367

Since a1 = −3a3 and a1 + a3 = 1 we obtain a3 = −1/2 and a1 = 3/2 so that eventually :

3 1
f (ξ) = ξ − ξ 3 
2 2

2)  The displacement thickness is obtained from :

ˆ δ ˆ 1 1
ξ3

∗ u 3 3
δ = (1 − )dy = δ (1 − f (ξ))dξ = 1 − ξ + = δ ≈ 0.375 δ
0 ue 0 2 2 0 8

Similarly we obtain : θ = (39/280) δ ≈ 0.14 δ.

The shape factor is such that : H = δ ∗ /θ ≈ 2.69. 

∂u
3)  The wall shear stress is such that τw = µ( )y=0 so that it can also be computed as :
∂y
µue 0 3 µue
τw = f (0) = .
δ 2 δ
With zero external pressure gradient, the von Kármán equation reduces to :
dθ
τw = ρu2e
dx
Injecting the above expressions for τw and θ yields :
dδ 3 µue
0.14ρu2e =
dx 2 δ
or else
d(δ 2 (x)) µ
2δdδ = ≈ 21.5 dx
dx ρue
which is readily integrated into :
δ 4.64
=√
x Rex
where the local Reynolds number is defined as Rex = ue x/ν.
The skin friction coefficient is obtained from :

τw 0.646
Cf = 1 2
≈√
2 ρue Rex

0.664
which is a good approximation (within 2% error) of the exact Blasius solution : Cf = √ . 
Rex

7.4.3 Exercise #3 : application of the Walz-Thwaites method

Question
Let us study the viscous flow over a convex corner (see Fig. 7.28). The inviscid solution is
known to be such that ue (x) = Cxm , with x the downstream distance along the wall counted
from the corner and m ∈] − 1, 0[ with the angle α (< 0) such that α = [m/(m + 1)] × π.
 368 Laminar boundary layer

α(<0)

u e(x)

x
Figure 7.28: Boundary layer separation on a convex corner.

Since ue (x) decreases with x along the inclined wall, an adverse pressure gradient applies on the
boundary layer which eventually leads to the separation of the boundary layer when the angle
α becomes large enough (in absolute value).

• Apply the Walz-Thwaites method to estimate the value of the angle α (or equivalently the
value of the exponent m) yielding a detached boundary layer.

Answer
Formula (7.34) is applied with ue (x) = Cxm :
ˆ x
θ2 0.45 0.45 x5m+1 0.45 1
= 6 6m C 5 t5m dt = 6m
=
ν C x 0 Cx (5m + 1) (5m + 1) Cxm−1
The value of the parameter λ is computed from :
θ2 due θ2 0.45m
λ= = Cmxm−1 =
ν dx ν (5m + 1)
Reporting now the tabulated values of the function S(λ) in Table 7.2, we note that separation
(τw = 0 for a 2D flow) corresponds to S(λ) = 0 and occurs for λ = −0.082. The corresponding
value of m is therefore such that :
0.45m
= −0.082
(5m + 1)
hence m = −0.0953 corresponding to an angle α ≈ −19◦ . This value is close to the exact value
α = −17.9◦ provided by the numerical solution of the Falkner-Skan equation.

7.4.4 Problem #1 : study of the boundary layer developing along an airfoil

Questions
Let us consider an airfoil along which a laminar boundary layer is developing. The external
inviscid velocity ue (x) can be approximated as :

ue (x) x x


 = 19 if 0 ≤ ≤ 0.1
 U∞
 c c

 ue (x) = 2 − x x


if 0.1 ≤ ≤1

U∞ c c
7.4. Exercises and problems 369

where U∞ denotes the upstream farfield velocity and c is the airfoil chord.
The Walz-Thwaites or Thwaites method is used to analyze the boundary layer which develops
on the airfoil with this outer velocity distribution.
x
1) • Compute the normalized momentum thickness θ/c for 0 ≤ ≤ 0.1.
c
U∞ c
Note : θ/c will be given as a function of the Reynolds number Re = where ν is the air
ν
kinematic viscosity.
x
2) • Show the normalized momentum thickness θ for 0.1 ≤ ≤ 1 is given by : :
c
 
 2
θ 0.075  49.522
= − 1

c Re
 x 6
2−
c
3) • Compute the location of the boundary layer separation point along the airfoil.

Answers
 1)  In Thwaites’ method, the momentum thickness of the boundary layer is directly linked
to the outer inviscid velocity profile through the relationship :
ˆ
0.45ν x 5
θ2 (x) = u (t)dt
u6e 0 e
If we use the non-dimensional coordinate x/c and the non-dimensional x-velocity component
ue /U∞ , this formula becomes :
6 ˆ x/c 
θ(x/c) 2 0.45 ue (t/c) 5
   
U∞
= d(t/c)
c Re ue (x/c) 0 U∞
U∞ c
where the Reynolds number Re is defined as Re = .
ν
For 0 ≤ x/c ≤ 0.1, ue /U∞ is such that :
ue (x) x
= 19
U∞ c
The momentum thickness for 0 ≤ x/c ≤ 0.1 is therefore such that :
 2 ˆ x/c ˆ x/c
θ 0.45 1 5 0.45 1 1
= (19(t/c)) d(t/c) = (t/c)5 d(t/c)
c Re (19(x/c)6 ) 0 Re 19 (x/c)6 0
The calculation of the integral is immediate :
ˆ x/c
1 x/c (x/c)6
(t/c)5 d(t/c) = (t/c)6 0 =
0 6 6
so that
 2
θ 0.45 1 1 (x/c)6 0.00395
= 6
=
c Re 19 (x/c) 6 Re
 370 Laminar boundary layer

for x/c ∈ [0, 0.1] or else

θ 0.06285
= √ .
c Re
still for x/c ∈ [0, 0.1]. 

2)  For 0.1 ≤ x/c ≤ 1, because of the new definition for the outer inviscid velocity distribution,
Thwaites’s formula yields :
 2 ˆ x/c 
ue (t/c) 5

θ 0.45 1
= d(t/c)
c Re (2 − x/c)6 0 U∞
A key point to underline is that the integral must be computed from the start of the boundary
layer development and takes therefore into account both the velocity distribution valid for x/c ∈
[0, 0.1] and the velocity distribution valid for x/c ∈ [0.1, 1] :
ˆ x/c  ˆ 0.1 ˆ x/c
ue (t/c) 5

5
d(t/c) = (19(t/c)) d(t/c) + (2 − (t/c))5 d(t/c)
0 U ∞ 0 0.1

It is immediate to compute :
ˆ x/c 
ue (t/c) 5 195 

0.1 1  x/c
d(t/c) = (t/c)6 0 − (2 − (t/c))6 0.1
0 U∞ 6 6

1 1
= 0.41268 − (2 − (x/c))6 + (2 − 0.1)6
6 6
1 1
= 0.41268 + 7.841 − (2 − (x/c))6 = 8.25368 − (2 − (x/c))6
6 6
Reinjecting into the formula giving (θ/c)2 yields :
 2  
θ 0.45 1 1 6
= × 8.25368 − (2 − (x/c))
c Re (2 − x/c)6 6

3.71416 1 0.075
= 6
−
Re (2 − x/c) Re
 
0.075 49.522
= −1
Re (2 − x/c)6

valid for 0.1 ≤ x/c ≤ 1. 

θ2 due
3) • Since the momentum distribution is available all along the airfoil, the parameter λ =
ν dx
can be computed in order to determine (estimate) the separation location for the boundary layer.
Let us recall that in the framework of Thwaites’ theory, the separation location corresponds to
the value λ = λsep = −0.082.
Using non-dimensionval variables, λ is readily expressed as :
 2
θ d(ue /U∞ )
λ= Re
c d(x/c)
7.4. Exercises and problems 371

For x/c ∈ [0, 0.1], it comes :

0.00395
λ= × Re × 19 = 0.075
Re
which means the boundary layer remains attached in this zone. This behaviour is expected since
the increasing outer inviscid velocity distribution corresponds to a negative pressure gradient
dpe /dx that is a pressure gradient favoring the attachment of the boundary layer.
For x/c ∈ [0.1, 1], it comes :
 
0.075 49.522
λ= − 1 × Re × (−1)
Re (2 − x/c)6

or else
3.71415
λ = 0.075 −
(2 − (x/c))6
which is negative for x/c ∈ [0.1, 1]. If boundary layer separation occurs, it will necessarily occur
in that region of the airfoil. Solving for λ = λsep = −0.082 yields :

3.71415
0.082 = − 0.075
(2 − (xsep /c))6
or
(2 − (xsep /c))6 = 23.657
that is 2 − (xsep /c) = 1.694 and eventually xsep /c = 0.306 ≈ 0.31. 
372 Laminar boundary layer
 Lecture 8
Turbulent boundary layer and
turbulence modeling

Contents
8.1 Turbulent boundary layer . . . . . . . . . . . . . . . . . . . . . . . . . 373
8.1.1 Averaging operators and Reynolds decomposition . . . . . . . . . . . . . 373
8.1.2 Reynolds-averaged Navier-Stokes (RANS) equations . . . . . . . . . . . 376
8.1.3 von Kármán Equation and turbulent boundary layer . . . . . . . . . . . 377
8.1.4 Velocity profile in the turbulent boundary layer . . . . . . . . . . . . . . 383
8.2 Turbulence modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 388
8.3 Exercises and problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 388
8.3.1 Exercise #1 : application of the von Kármán equation . . . . . . . . . . 388
8.3.2 Exercise #2 : friction over a flat plate . . . . . . . . . . . . . . . . . . . 390
8.3.3 Exercise # 3 : simplified model of a turbulent jet . . . . . . . . . . . . . 393
8.3.4 Problem # 1 : total drag of an airfoil . . . . . . . . . . . . . . . . . . . 395

8.1 Turbulent boundary layer

A turbulent flow displays irregular features : velocity components, pressure and density fluctuate
in random fashion. Note these fluctutations of the flow properties favor exchanges within the
flow : increased momentum exchange leads for instance to an increase of the wall shear stress
or skin friction. The turbulent motion is unsteady and three-dimensional. However, it can be
considered as the superposition of a mean or averaged flow, often steady, and random fluctuations
(see Figure 8.1).

8.1.1 Averaging operators and Reynolds decomposition

Several averaging processes can be defined depending on the flow under study. A general average
is the so-called ”ensemble average” defined as an average over a set of flow realizations. Let us
suppose for instance the flow over a flat plate is repeated N times and let us denote fn (~x, t) the
 374 Turbulent boundary layer

FLUCTUATIONS DANS UNE COUCHE LIMITE TURBULENTE

NOTION D’ECOULEMENT MOYEN

Figure 8.1: Top : velocity profiles (in the normal direction to a plate) across a turbulent
boundary layer over a flat plate plotted at successive times ti for a fixed position along the
plate. Instantaneous flow evolution appears random. Down : superposition of the above
profiles (left) showing a mean velocity profile (right) around which the flow fluctuates over
time.
8.1. Turbulent boundary layer 375

measurement of the physical flow quantity f (~x, t) during the nth flow realization. The ensemble
average is then defined as :
N
1 X
f (~x, t) = lim fn (~x, t)
N →∞ N
n=1

It is also possible of course to define a time average :

ˆ t+T
1
f (~x) = lim f (~x, t)dt
T →∞ T t

Such an average is very commonly used since a large number of turbulent flows can be considered
as statistically steady. In practice, T cannot be chosen infinite but it is enough to select T large
with respect to the characteristic time τ of the turbulent fluctuations so that the time average
is rather defined as :
ˆ
1 t+T
f (~x) = f (~x, t)dt with T >> τ
T t
A space average can also be relevant for other flows :
ˆ
1
f (t) = lim f (~x, t)dv
V →∞ V V

Once the average process selected - it can be the above defined time-average but what follows is
valid for other types of average processes - any instantaneous flow variable φ can be decomposed
into a mean or averaged part φ and a fluctuating part φ0 :

φ = φ + φ0

Such a decomposition is called a Reynolds decomposition since Reynolds was the first to intro-
duce it. Let us briefly review some key properties of the averaging.
The fluctuation φ0 is such that its average or mean value is zero. Indeed :

φ0 = φ − φ = φ − φ = φ − φ = 0

It is particularly interesting to compute the average of the product of two instantaneous quan-
tities :
φψ = (φ + φ0 )(ψ + ψ 0 ) = φ ψ + φψ 0 + φ0 ψ + φ0 ψ 0
so that
φψ = φ ψ + φψ 0 + ψφ0 + φ0 ψ 0

= φ ψ + φ ψ 0 +ψ φ0 +φ0 ψ 0
|{z} |{z}
=0 =0

and eventually :
φψ = φ ψ + φ0 ψ 0 (8.1)
This important result will be used next when deriving the averaged Navier-Stokes equations.
The averaging operator is assumed to commute with the time and space differentiation operators.
 376 Turbulent boundary layer

8.1.2 Reynolds-averaged Navier-Stokes (RANS) equations

Let us consider the Navier-Stokes equations written for an incompressible and homogeneous flow
(ρ = constant, ν = constant) :

∂ui
=0


∂xi

∂ui ∂(ui uj ) 1 ∂p ∂ 2 ui

 + =− +ν
∂t ∂xj ρ ∂xi ∂xj ∂xj


Let us introduce in these instantaneous equations governing the flow (be it in laminar or turbu-
lent regime) the Reynolds decomposition for velocity and pressure :

ui = ui + u0i


p = p + p0

and let us take the average of these governing equations. Taking into account the properties
reviewed above (commutation between the averaging operator and the differentiation operators,
identity (8.1)) , it is easy to derive the following Reynolds equations or Reynolds-Averaged
Navier-Stoke equations (in short RANS equations) :

∂ui

 ∂x = 0


i
∂ui ∂(ui uj ) ∂(u0i u0j ) 1 ∂p ∂ 2 ui
+ + =− +ν



∂t ∂xj ∂xj ρ ∂xi ∂xj ∂xj

∂ 2 ui 1 ∂τij
Since ν = with ~
~τ the viscous stress tensor computed using the mean velocity
∂xj ∂xj ρ ∂xj
fields, one can also write :

∂ui
=0


 ∂xi


(8.2)
 ∂ui ∂(ui uj ) 1 ∂p 1 ∂
(τij − ρu0i u0j )


 + =− +
∂t ∂xj ρ ∂xi ρ ∂xj


With respect to the case of a laminar flow, the only (but crucial !) change in the conservation laws
introduced by the averaging process lies in the stress tensor. It is now made of two contributions :
the Reynolds stress tensor defined as

τijR = −ρu0i u0j

must be added to the viscous stress tensor ~~τ . The Reynolds stress tensor expresses the influence
of the turbulent fluctuations on the mean flow motion. Note that system (8.2) is not closed
since, in 3D, it is made of 4 scalar equations linking the 4 mean flow unknowns ui , p but
also the components of the Reynolds stress tensor. Since this tensor is obviously symmetric, it
introduces 6 new unknowns (in 3D). In order to close (8.2), the Reynolds stress tensor must be
expressed in terms of mean or averaged quantities. A classical assumption, known as Boussinesq
8.1. Turbulent boundary layer 377

hypothesis, is to consider the Reynolds stress tensor is proportional to the mean flow deformation
~~
tensor D :
∂ui ∂uj
−ρu0i u0j = µt ( + )
∂xj ∂xi
or
∂ui ∂uj
−u0i u0j = νt ( + )
∂xj ∂xi
where µt is the dynamic turbulent viscosity or eddy viscosity and νt = µt /ρ is the kinematic
turbulent viscosity. It must be emphasized these relationships are nothing more than modelling
assumptions, which in practice have a limited validity. Note also that the eddy viscosity is not
a fluid property, as is the physical viscosity µ or ν, but a flow property. The eddy viscosity is a
flow field ν(~x, t). Section 8.2 will briefly describe how the eddy viscosity can be expressed from
the mean flow quantities in order to fully close the RANS equations. First, we will see how the
friction can be estimated in the turbulent regime using the tools (in particular the von Kármán
equation) presented in the previous Lecture # 7.

Following a line of thought similar to the one presented in the laminar case, the fact the thickness
of the turbulent boundary layer is small with respect to the characteristic dimensions of the flow
problem (for instance the length L of the flat plate along which the turbulent boundary layer
develops) can be exploited to derive a simplified set of equations (Prandtl equations) describing
the flow in the turbulent boundary layer :
 ∂u ∂v
 + =0
∂x ∂y







∂u ∂u due 1 ∂τ
u +v = ue +
∂x ∂y dx ρ ∂y








p = pe (x)

with
∂u
τ = τ lam + τ turb = µ − ρu0 v 0
∂y
Notations have been simplified with φ denoting now the mean value φ. It can be observed the
Prandtl equations for the mean velocity components in a turbulent boundary layer are formally
identical to the Prandtl equations derived for a laminar boundary layer and applied to the
instantaneous velocity components of the laminar flow with a key difference in the definition
of the stress tensor which contains two contributions in the turbulent case : the viscous or
”laminar” shear stress τ lam and the turbulent or Reynolds shear stress τ turb = −ρu0 v 0 .

8.1.3 von Kármán Equation and turbulent boundary layer

The integral equation of von Kármán (denoted VKE from now on) can be derived from the
turbulent Prandtl equations exactly as it has been derived in the laminar case since the stress
tensor is zero when y → +∞ whatever the flow regime. It is therefore possible to write :

Cf dθ 1 due
= + (H + 2) θ
2 dx ue dx
 378 Turbulent boundary layer

with the displacement thickness δ ∗ , the momentum thickness θ and the shape factor H defined
as previously done for the laminar boundary layer, keeping in mind the mean velocity component
u, here denoted u, is actually used in these definitions. The velocity distribution ue remains the
outer invisicd distribution - that is the velocity field solution of the Euler equations, assuming
there is no boundary layer developing along the wall. The skin friction coefficient is the non-
dimensional form of the full stress tensor :
τw
Cf = 1 2
2 ρue

with τw = τ (x, y = 0) = τ lam (x, y = 0) + τ turb (x, y = 0). A key difference with the lami-
nar case is that τw could be computed from the velocity field for the laminar boundary layer
∂u
(τw = µ( )y=0 ) while, for the turbulent boundary layer, an explicit formulation of Cf as a func-
∂y
tion of the mean velocity u, also denoted u, requires a turbulence closure for the Reynolds stress.

It is however possible to apply VKE without choosing a turbulence model if the assumed velocity
profile in the boundary layer is completed with a friction law for Cf . This is an important dif-
ference with respect to the solution of VKE in the laminar case where only an assumed velocity
profile is required.

The eddy structures existing in the turbulent boundary layer enhance the momentum transfer
between the upper region of the boundary layer where the velocity is close to the outer inviscid
velocity ue and the inner layers, close to the fixed wall, where the velocity goes to zero and
is therefore small. The macroscopic turbulent mixing is much more efficient than the diffusive
momentum mixing - which is the only mixing process taking place in the laminar boundary layer.
The momentum brought near the wall by the turbulent eddies has important consequences on
the physical features of the turbulent boundary layer :
-• i) the turbulent boundary layer is thicker than the laminar boundary layer

-• ii) the velocity profile in the turbulent boundary layer increases much more rapidly with
the normal distance to the wall (see Figure 8.2) so that the friction applied by a turbulent
boundary layer is much larger than the friction applied by a laminar boundary layer. This
can be perceived as a drawback of the turbulent boundary layer with respect to the laminar
one but it must be balanced with the next feature.

-• iii) because of its velocity profile (large value of the velocity very near the wall), the
turbulent boundary layer is much less likely to separate than the laminar boundary layer.
More precisely, more intense adverse pressure gradients are needed to induce the separation
of a turbulent boundary layer.
This balance between friction and separation is well exemplified when considering the flow over
a sphere. Figure 8.3 displays the laminar flow over a sphere while Fig. 8.4 displays the turbulent
flow over the same sphere. In the laminar case, the laminar boundary layer developing along
the sphere tends to separate early, before reaching the top or the bottom of the sphere as can
be observed in Fig.8.3. Once the boundary layer separates, a large wake develops downstream
of the separation location. The total drag of the sphere in that case is therefore the sum of a
small friction drag over the front of the sphere where the boundary layer is both attached and
laminar and a very large pressure drag (also said ”form drag”) over the rear part of the sphere
8.1. Turbulent boundary layer 379

u/u e

couche limite laminaire

couche limite turbulente

y 1/2
(Rex )
x
Figure 8.2: Typical x-velocity profiles in a boundary layer developing over a flat plate plotted
versus a non-dimensional distance to the wall. In the turbulent case, the x-velocity increases
much faster when moving away from the plate, inducing both a larger friction but also a better
resistance to separation.

where a wake flow is taking place.

In the turbulent case, the turbulent boundary layer developing along the sphere displays a
separation which is significantly delayed with respect to the laminar case because the turbulent
boundary layer needs a much stronger adverse pressure gradient to eventually separate. Note
that in that case the trip wire observed on the front part of the sphere forces the transition of
the boundary layer to the turbulent regime. The turbulent boundary layer separates on the rear
part of the spheren which leads to the formation of a narrow wake region. The total drag of
the sphere in that case is therefore the sum of a large friction drag over the region of the sphere
where the boundary layer is both attached and turbulent and a small pressure drag or form drag
over the rear part of the sphere since the turbulent wake remains of limited extent. Eventually,
it is found out that the total drag for the turbulent flow over the sphere is actually smaller than
the total drag for the laminar flow. This explains for instance why golf balls display small holes
on their surface, which ensure the flow over the ball will be indeed turbulent.

These remarks on the turbulent boundary layer shall guide the choice of a relevant assumed
velocity profile for solving the VKE. Let us indeed recall that an accurate velocity profile leads
to a better estimate of the friction. A commonly used velocity profile for a laminar boundary
layer developing with small to moderate outer pressure gradients is a power law of the form :
u y
= ( )1/n (8.3)
ue δ

with a typical value of 7 for n.

With respect to the laminar case, a friction law (also derived from experimental observations)
must be provided. It is generally of the form : Cf as a function of the momentum thickness θ or
the Reynolds number Reθ based on θ and ue , that is Reθ = ue θ/ν; or else Cf as a function of the
 380 Turbulent boundary layer

Figure 8.3: Laminar flow over a sphere at Re = 15000 (based on the sphere diameter and the
far-field velocity). ONERA photographs taken by Werlé (1980). Left : instantaneous flow.
Right : mean flow.

Figure 8.4: Turbulent flow over a sphere (equipped with a trip wire) at Re = 30000 (based on
the sphere diameter and the far-field velocity). ONERA photographs taken by Werlé (1980).
Left : instantaneous flow. Right : mean flow.
8.1. Turbulent boundary layer 381

boundary layer thickness δ or Reδ . For instance, under assumptions which are not detailed here,
the skin friction of the turbulent boundary layer developing along the inner wall of a cylindrical
pipe can be estimated as :
0.0464
Cf = 1/4
(8.4)
Reδ
K
Let us assume now a velocity profile of the form (8.3) and a friction law of the form Cf = 1/p
Reθ
(generalization of (8.4)). Let us also assume there is no outer pressure gradient so that VKE
reduces to :
dθ
Cf (Reθ ) =
dx
or else, using the assumed friction law :

dθ 1/p
θ = K/(ue /ν)1/p
dx
This is an ODE for the function θ(x) which can be easily solved (see also the Exercises at the
end of the Lecture) to obtain θ(x) hence Cf (x) that is the skin friction as a direction function
of the distance x along the plate.

If the outer pressure gradient is non-zero and if the power-law velocity profile is used to estimate
δ ∗ , θ hence H and if the friction law is taken into account, the VKE reads :

dθ 1 due Cf (Reθ )
+ [(H(n) + 2) ]θ =
dx ue dx 2

which is again an ODE on θ(x), which can be solved to obtain θ(x) hence Cf (x).

The skin friction coefficient Cf estimated for a turbulent boundary layer (case of a zero outer
pressure gradient) can be plotted as a function of the local Reynolds number Rex (based on the
distance x from the plate’s leading edge and the outer velocity ue (x)) along with the laminar
skin friction coefficient derived for the Blasius solution in the laminar case :
0.664
(Cf )lam = 1/2
Rex
Figure 8.5 displays both the laminar skin friction distribution and the turbulent skin friction
distribution. The boundary layer tends to remain laminar as long as the local Reynolds number
does not exceed a critical value between 3 × 105 and 3 × 106 . For Rex = Recrit , a transition
from the laminar regime to the turbulent regime takes place along the plate. This transition
corresponds to a strong increase of the skin friction coefficient which can be quantified using
the above formulae. Obviously it can be interesting to preserve the laminarity of the flow so
as to keep a low friction drag. For Rex = 106 for instance, (Cf )turb = 3.73 × 10−3 while
(Cf )lam = 0.664 × 10−3 , corresponding to a value 5 times smaller than the value reached in
the turbulent flow. Naturally, this reduction of the friction drag must be balanced with the
occurrence of undesirable flow separation.

Note that formula (8.4) for the friction law is absolutely not universal. Many other empirical
laws exist and can be used depending on the flow configuration under study. Some of these
 382 Turbulent boundary layer

0.01

Turbulent
Coefficient de frottement Cf

Transition

0.001
Laminaire

0.0001
5 6 7
10 10 10
Nombre de Reynolds local

Figure 8.5: Evolution of the skin friction coefficient along a flat plate.

friction laws take into account a non-zero outer pressure gradient through the value of the shape
factor H (which is directly impacted by dpe /dx). For instance, the Ludwieg and Tillman formula
(1949) expresses Cf as :
0.246 −0.678H
Cf = 10
Reθ0.268
It is clear however that using this formula in the VKE can make the solution of the equation
difficult to compute. In practice, it is often preferred to use a power law velocity profile and a
friction law depending on Reθ or Reδ only, with the outer pressure gradient effect accounted for
through due /dx and a modified value of the shape factor H.

Another simple way to proceed when solving the VKE is to follow the method described in
7.3.4 for the case of a laminar boundary layer, adapting the value of the shape factor H and the
coefficients m, b pf the friction law to the case of a flat plate turbulent boundary layer. A typical
value for the shape factor of a turbulent boundary layer developing along a flat plate with a
moderate outer pressure gradient is H = 1.4. The coefficients of the friction law expressed as a
function of the Reynolds number based on the momentum thickness can be taken equal to the
typical values : m = 1/5 and b = 0.0086. Applying then formula (7.35) allows to compute the
evolution of the momentum thickness θ(x) for a given outer inviscid velocity distribution ue (x).
The friction law Cf (x) can then be obtained and the evolution of the displacement thickness is
given by δ ∗ = θH with the shape factor H equal to the previously specified value.
When computing for instance the characteristic features of the boundary layer developing along
an airfoil, the boundary layer can be assumed laminar near the leading edge and the stagnation
point. Falkner-Skan solution allows to estimate the momentum thickness of the boundary layer
at the stagnation point. The integral formula (7.35) can be applied in the laminar case (with
values of H, m, b corresponding to the laminar boundary layer over a flat plate) from the
stagnation point x0 = 0 to the transition point x = xt which can be estimated using an empirical
8.1. Turbulent boundary layer 383

transition criterion. This transition point becomes the initial point for the analysis of the
turbulent boundary layer. Assuming the continuity of the momentum thickness in the transition
zone allows to estimate the initial value of the turbulent momentum thickness :

θturb (x0 ) = θturb (xt ) = θlam (xt )

The evolution of the momentum thickness downstream of xt , for instance from x = x0 = xt to

x = c with c the airfoil chord, is obtained by applying formula (7.35) with the values of H, m,
b corresponding to the turbulent boundary layer over a flat plate. Such an approach is detailed
in the problem proposed in 8.3.4.
We show in the following section how the skin friction can also be estimated from a more detailed
knowledge of the velocity profile in the turbulent boundary layer - more detailed than a rough
approximation using a power-law description.

8.1.4 Velocity profile in the turbulent boundary layer

For the laminar boundary layer developing over a flat plate with zero outer pressure gradient,
it was possible to obtain an exact solution, the Blasius solution. It is no longer possible to
obtain such an equivalent exact solution for the turbulent boundary layer because of the closure
problem introduced by the Reynolds stress tensor. It is nonetheless possible to derive a refined
description of the velocity profile inside a turbulent boundary layer from a thorough analysis
of the boundary layer equations. Let us consider therefore the x-momentum equation written
inside a turbulent boundary layer with zero outer pressure gradient :

∂u ∂u ∂ ∂u ∂
u +v = (ν ) + (−u0 v 0 )
∂x ∂y ∂y ∂y ∂y
| {z } | {z } | {z }
inertia laminar/viscous friction turbulent friction

Universal velocity profile

The turbulent boundary layer can be decomposed into several zones or layers / sub-layers ac-
cording to the dominant term in the previous equation :

-• Very close to the wall (assumed at y = 0), the viscous friction is dominant because the
mean velocity components and the velocity fluctuations tend to 0 at the wall. Therefore,
when y → 0, the x-momentum equation reduces to :

∂2u
µ ≈0
∂y 2

hence
u = Cy (8.5)
with C a constant precisely equal to τw /µ. This zone very close to the wall is called the
viscous sub-layer. In the viscous sub-layer, the characteristic scale for the velocity is the
so-called friction velocity uτ , defined from the wall shear stress as :
r
τw
uτ =
ρ
 384 Turbulent boundary layer

Similarly, a characteristic length scale ν/uτ , relevant for the viscous sub-layer can be
defined from the kinematic viscosity ν and the friction velocity so that the following non-
dimensional velocity u+ and distance to the wall y + can be introduced :
u uτ y
u+ = , y+ = .
uτ ν

Note that y + is similar to a Reynolds number based on the distance y to the wall, the
friction velocity uτ and the fluid kinematic viscosity. It is thus possible to rewrite (8.5) in
the non-dimensional form :
u+ = y + (8.6)
Experiments show that the relationship (8.6) remains valid until y + reaches a value between
5 and 10. For larger values of y + , the turbulent friction is no longer negligible with respect
to the viscous friction and the shape of the velocity profile must be modified accordingly.

-• When leaving the viscous sub-layer, the turbulent friction term becomes progressively
dominant with respect to the viscous friction while the inertia term remains negligible
because of the low velocities existing in this flow region which is still close to the plate
wall.

-• When one keeps moving away from the wall to reach the upper boundary of the boundary
layer, the inertia term and the turbulent friction term becomes the dominant terms in the
x-momentum equation while the viscous friction becomes negligible. When approaching
the upper boundary of the boundary layer, the mean velocity u tends to the outer inviscid
velocity ue . It is observed : i) the difference between u and ue in this zone is of the order
of magnitude of the friction velocity uτ , ii) the non-dimensional velocity deficit (u − ue )/uτ
depends only on the distance y to the wall normalized by the characteristic length in this
flow region which is the thickness δ of the boundary layer. It is thus possible to write :
ue − u
= g(η) (8.7)
uτ

where η = y/δ.

-• Finally, in the outer region, only the inertia term is dominant, the viscous and turbulent
friction terms being both negligible.

In the zone where the viscous and turbulent friction have the same order of magnitude, one can
write :
∂ ∂u ∂
(µ ) + (−ρu0 v 0 ) ≈ 0
∂y ∂y ∂y
∂u
so that the total friction τ = µ − ρu0 v 0 is approximately constant and equal to the wall shear
∂y
stress τw . This remark is at the root of the development of so-called turbulent wall laws. The
mean velocity in this zone which is still close to the wall remains of the order of the friction
velocity uτ with the characteristic length still given by ν/uτ . Therefore, the velocity profile in
this zone can be assumed of the form :

u+ = f (y + ) (8.8)
8.1. Turbulent boundary layer 385

Two main layers can thus be distinguished in the turbulent boundary layer : an outer or external
layer in which the characteristic length is the boundary layer thickness δ eand the velocity deficit
u − ue is of the order of the friction velocity uτ ; an inner or internal layer in which the fluid
viscosity effects are significant so that the characteristic length scale is ν/uτ while the mean
velocity is of the order of the friction velocity uτ . Note that the ”outer” layer is inside the
boundary layer (in its upper part) and should not be confused with the outer inviscid flow
region outside the boundary layer. The non-dimensional variable η = y/δ is used in the outer
layer : η = 1 at the upper boundary of the boundary layer whereas η → 0 when approaching
the interface with the inner layer.
The non-dimensional variable y + is used to describe the inner layer : y + → 0 when approaching
the wall whereas y + → ∞ when approaching the interface with the outer layer.
In order to explicitly determine the functions f and g which have been previously introduced,
let us use the fact the inner velocity profile given by (8.8) should match, when y + → ∞, the
outer velocity profile given by (8.7), when η → 0. This is expressed by :
ue
f (y + ) = − g(η) (8.9)
uτ

Taking the ratio of the non-dimensional variables y + and η yields :

y+ uτ δ
= = Reτ
η ν
where Reτ denotes the Reynolds number based on the friction velocity and the boundary layer
thickness. The matching condition (8.9) can be recast in the form :
ue
f (Reτ η) = − g(η)
uτ
Differentiating this relationship with respect to η and multiplying the result by η in order to
reintroduce y + yields :
y + f 0 (y + ) = −ηg 0 (η)
This functional identity is satisfied if and only if :

y + f 0 (y + ) = constant , −ηg 0 (η) = constant

An immediate integration yields the form of the velocity profile in the inner layer when y + → ∞ :
1
u+ = ln(y + ) + B (8.10)
κ
where κ denotes the von Kármán constant, the value of which has been experimentally found
such that κ ≈ 0.40 or 0.41. The value of the integration constant B lies between 5 and 5.2. The
relationship (8.10) is a logarithmic law which does not depend on the parameters of the external
layer which tends to make this relationship universal - at least as long as the external pressure
gradient remain moderate. The relationship (8.10) is typically valid from y + ≈ 50 to y + ≈ 500.

Similarly, the following expression is obtained for the velocity profile in the outer layer when
η→0:
ue − u 1
= − ln(η) + A (8.11)
uτ κ
 386 Turbulent boundary layer

with A ≈ 2.35. This logarithmic law for the velocity deficit applies in the external / upper
region of the turbulent boundary layer, typically from η = y + /Reτ with y + ≈ 500 up to the
upper limit of the boundary layer. Since Reτ depends on the global Reynolds number (based
on the outer velocity ue and the plate length L), the location of the interface between the in-
ner layer and the outer layer of the turbulent boundary layer depends on this global Reynolds
number. Assuming Re = O(106 ) - order of magnitude for the characteristic Reynolds number in
the context of external aerodynamic flows -, Reτ = O(103 ) and the interface between the inner
and outer layers lies around y ≈ 0.1δ. Note that, when expressed as a function of the external
variable η, the thickness of the viscous sub-layer is of the order of 0.001δ, thus a very thin region
(but crucial regarding friction).

At this stage, a velocity profile has been obtained in the inner layer when y + → 0 and when
y + → ∞ as well as in the outer layer when η → 0. We are still missing a profile for u when y + is
neither very small nor very large as well as a profile for the velocity deficit ue −u when η does not
tend to zero but on the contrary increases to reach η = 1 at the upper limit of the boundary layer.

The matching between the linear law valid in the viscous sub-layer and the logarithmic profile
valid in the upper part of the inner layer is achieved in a progressive way through a so-called
buffer layer. In this buffer layer, the velocity profile can be described by an implicit relationship
due to Spalding (1961) :

+ 1 1
y + = u+ + e−κB (eκu − 1 − κu+ − (κu+ )2 − (κu+ )3 ) (8.12)
2 6
In the logarithmic or inertial zone, which ensures the transition from the near-wall region where
the viscous friction is dominant to the region where the turbulent friction becomes dominant, the
velocity profile is given by (8.10) using internal variables and by (8.11) using external variables.
When moving away from the wall in the external layer, the velocity deficit no longer depends
only on η = y/δ but also varies as a function of the outer pressure gradient - outer meaning here
outside the boundary layer, in the inviscid flow region. Coles proposed in 1952 to introduce a
correction to the logarithmic law in the form :
1 Π
u+ = ln(y + ) + B + × W (η) (8.13)
κ κ
where the parameter Π, known as the wake parameter, depends on the outer inviscid pressure
gradient through the β parameter introduced by Clauser (1954) :

δ ∗ dpe
Π = Π(β) = Π( )
τw dx

The correction function W (η), known as the wake correction, has been proposed by Coles and
reads :
π
W (η) = 2sin2 ( η).
2
As for the wake parameter, its value is deduced from experimental observations and is around
0.5 for a flat plate with zero outer pressure gradient. In the general case, White proposes the
following empirical expression : l’expression empirique suivante :

Π ≈ 0.8(β + 0.5)0.75
8.1. Turbulent boundary layer 387

Note that formula (8.13) is a composite formulation mixing the internal scale y + and the external
scale η. Following Simpson (1970) and Purtell et al. (1981) it is possible to express the correction
of the velocity deficit law as :

ue − u 1 Π Π
= − ln(η) − × W (η) + 2 (8.14)
uτ κ κ κ

Fig. 8.6 summarizes the different regions of a turbulent boundary layer. Note first that the

30 région INTERNE

25
sous-couche visqueuse zone tampon sous-couche
(loi linéaire) (loi de inertielle
20 Spalding) (loi log)
zone de
u =u/uτ

zone de
vitesse sillage
15 (loi de
déficitaire
+

(loi log sillage)

pour la
10 vitesse
déficitaire)

5
région EXTERNE
0 -1
10 100 101 102 y=0.1 δ 103 104
y=δ
y+=uτ y / ν

Figure 8.6: Evolution of the mean velocity profile u (or u) through a turbulent boundary layer.

velocity profile is plotted using the internal variables (u+ , y + ), which produces a magnifying
effect on the internal region of the boundary layer which typically extends to y ≈ 0.1δ and in
particular on the viscous sub-layer, the thickness of which is about 1/1000 of the boundary layer
thickness.
Knowing the velocity profile in the turbulent boundary layer is very useful when building a
turbulence model since it is obviously a validation test-case for a turbulence closure. The
logarithmic law (8.10), also known as the wall law, is also very useful since it can be used in
the estimation of the skin friction of wall shear layer without needing to discretize the boundary
layer down to the wall.
This velocity profile can also be used as an assumed profile when solving the VKE. However,
the analytical solution of the differential equation obtained in that case can become complex so
that the use of a power law such as (8.3) is often preferred. Note such a power law correctly
describes the velocity profile computed in the above analysis over the whole external layer, that
is typically over 90% of the boundary layer thickness. This is illustrated in Fig. 8.7 where the
power-law velocity profile (8.3) (with the recommended value n = 7 for a flate plate with zero
outer inviscid pressure gradient) is plotted and compared with Coles’ wake law (8.13). This law
 388 Turbulent boundary layer

can be expressed using the outer variable η :

u uτ 1 Π
= ( ln(Reτ δ) + A + W (η))
ue ue κ κ

with Π = 0.5 for a boundary layer with zero outer pressure gradient. The ratio ue /uτ is
determined by noting that u = ue when δ = 1 so that :

1
ue /uτ = ln(Reτ ) + A + 2Π/κ
κ
Eventually, the wake law is entirely determined by the specification of the Reynolds number
Reτ : in Fig. 8.7 the value Reτ = 1 × 104 has been retained. The agreement between the simple
power-law approximation and the wake law is quite satisfactory in the outer layer.

0.9 Profil de vitesse théorique

Loi en puissance (n=7)
0.8

0.7

0.6
y/δ

0.5

0.4

0.3

0.2

0.1

0
0 0.25 0.5 0.75 1
u/ue

Figure 8.7: Approximation of the velocity profile in the turbulent boundary layer using a
power-law velocity distribution.

8.2 Turbulence modelling

This section will be made available in a future version of the Lecture notes.

8.3 Exercises and problem

8.3.1 Exercise #1 : application of the von Kármán equation
Let us consider a turbulent boundary layer developing along a flat plate with zero outer pressure
gradient. Experimental observations show that :
8.3. Exercises and problem 389

-• the velocity profiles (along y, at constant x) in this boundary layer are rather well approx-
 1/7
u(x, y) y
imated by the power law : =
ue δ(x)
-• the wall friction distribution can be approximated by the law (8.4), which links the skin
friction coefficient Cf and the Reynolds number based on the boundary layer thickness
δ(x).

• Determine the wall friction distribution for Cf (x) along the flat plate from these assumed
velocity profiles and skin friction law using the von Kármán equation.

 Let us denote ξ = y/δ(x) and u/ue = f (ξ) = ξ 1/7 . The displacement thickness and momentum
thickness are such that : ˆ 1
∗ δ(x)
δ (x) = δ(x) (1 − ξ 1/7 )dξ =
0 8
and ˆ 1
7δ(x)
θ(x) = δ(x) ξ 1/7 (1 − ξ 1/7 )dξ =
0 72
The shape factor is therefore equal to H = 9/7 ≈ 1.285.
With a zero outer pressure gradient or constant outer inviscid velocity, the von Kármán equation
reduces to :
dθ 7 dδ
Cf = 2 =
dx 36 dx
Using now the friction law (8.4) allows to turn this relationship into an ODE for δ(x) :

1/4 7 dδ
0.0464 = Reδ ×
36 dx
which can also be expressed using systematically Reδ instead of δ. Let us recall Reδ = ue δ/ν
while Rex = ue x/ν. Therefore we can write :

7 d(Reδ )
0.0464 = (Reδ )1/4
36 d(Rex )

or else
1/4
Reδ d(Reδ ) = 0.2386 d(Rex )
Hence :
4 5/4
Reδ = 0.2386Rex
5
or
5/4
Reδ = 0.2983Rex
This eventually leads to the following friction law :

0.0592
(Cf )turb = 1/5
(8.15)
Rex

which can be directly exploited to estimate the friction along a plate of given length L since
expressed as a function of x, the distance from the plate’s leading edge. 
 390 Turbulent boundary layer

0.9
Regime turbulent
0.8

0.7
y Regime laminaire
0.6

y/δ
0.5
Ue(x)
0.4

0.3

0.2

δ (x) 0.1
x 0
0 0.25 0.5 0.75 1
U/U e
x
Figure 8.8: Main features and velocity profiles for the boundary layer developig along a flat
plate with zero outer pressure gradient.

8.3.2 Exercise #2 : friction over a flat plate

We wish to asses the influence of flow regime (laminar or turbulent) over the friction applied
by a fluid flowing along a flat plate. Experiments show the velocity profiles along y (normal
direction to the plate) can be approximated as :
 1/n
u(x, y) y
= (8.16)
ue (x) δ(x)
with n = 1 (this is a rough approximation) in the laminar case and n = 7 in the turbulent case
(see Fig. 8.8) in the case where the flow takes place with zero outer pressure gradient so that
ue (x) = cste = ue .

1) • Express the von Kármán equation as :

dδ
Cf = C(n) (8.17)
dx
where C(n) is a coefficient depending on the power law exponent in (8.16) which will be detailed.

• Deduce from (8.17) a simple expression for the friction drag coefficient in the case of the flow
over a flat plate of length L.

due (x)
1)  If the outer pressure gradient is zero, that is also if = 0 (from the Bernoulli equation
dx
valid in the inviscid flow outside the boundary layer region), the von Kármán equation reduces
to :
Cf dθ
=
2 dx
The momentum thickness θ is defined as :
ˆ δ
u u
θ= (1 − )dy.
0 ue ue
8.3. Exercises and problem 391

From the assumed form of the velocity profile in the boundary layer, we can introduce the
y
non-dimensional quantity Y = and compute θ as :
δ
ˆ 1
θ=δ Y 1/n (1 − Y 1/n )dY
0

which yields :
n
θ= δ
(n + 1)(n + 2)
The skin friction coefficient can therefore be expressed as :

2n dδ
Cf = 
(n + 1)(n + 2) dx

 By definition of the drag coefficient :

ˆ L
1
CD = Cf (x)dx
L 0

Assuming δ(0) = 0, that is the thickness of the boundary layer is zero at the leading edge of the
plate (x = 0), the drag coefficient reads :

δ(L)
CD = C(n) 
L

2) Case of a laminar flow (n = 1)

By definition, the local skin friction coefficient in a laminar flow reads :
τw
Cf = 1 2
2 ρue

where τw is the local wall shear stress.

• Solve (8.17) to find δ(x) and compute the drag coefficient as a function of the Reynolds number
ue L
ReL = .
ν
 The wall shear stress reads : 
∂u 
τw = µ ( )
∂y y=0
so that
u 
∂(
) 
τw νue ue 
= ( y )
ρ δ ∂( ) 
δ y=0
| {z }
=1 from (8.16)

2ν dδ
Consequently, the skin friction coefficient reads Cf = and since Cf = C(n) , an ODE for
ue δ dx
δ(x) is obtained :
d(δ) 2ν
C(1)δ =
dx ue
 392 Turbulent boundary layer

Since C(1) = 1/3, it comes :

d(δ 2 ) 12ν
=
dx ue
which is readily integrated into : √
2 3
δ(x) = r
ue
νx
since δ(0) = 0. For x = L (trailing edge of the plate) we obtain :
√
δ(L) 2 3
=√
L ReL
so that after replacement in the CD formula obtained in 1) the following estimate is obtained
for the drag coefficient in the laminar case :
√
(2/ 3)
(CD )lam = √ = 1.155(ReL )−1/2 (8.18)
ReL

This estimate can be compared with the exact√ value derived from the Blasius solution which
was computed in Lecture # 7 : CD = 1.328/ ReL . Formula (8.18) underestimates the drag
coefficient by about 13%, which is quite reasonable considering the initial rough estimate (linear
profile (8.16) with n = 1) used to obtain (8.18). 

3) Case of a turbulent flow (n = 7)

In the turbulent case, there is no longer a direct relationship between Cf and τw since the friction
also includes a turbulent contribution which is not a priori known - a turbulence model would be
needed to express the Reynold stress. Therefore, one relies here on experimental data indicating
the thickness of the boundary layer developing over a flat plate with zero outer pressure gradient
can be estimated using a law of the form :
δ(x)  u x −1/10
e
=κ (8.19)
x ν
with typically κ = 0.092. It is assumed here the boundary layer is forced to become turbulent
at the plate’s leading edge so that δ(0) ≈ 0 when x = 0 at the plate’s leading edge.
• Compute the drag coefficient defined in 1) for a turbulent flow.

 In the turbulent case, the drag coefficient is given by :

δ(L)
CD = C(7)
L
so that using the proposed law of variation (8.19) for δ(x) it is immediate to express the boundary
layer thickness at x = L as :
δ(L) ue L −1/10
= κ( )
L ν
which yields after replacement in the drag expression :

(CD )tur = C(7)κ(ReL )−1/10

8.3. Exercises and problem 393

or else
(CD )tur = 0.0179(ReL )−1/10
Plotting both the laminar and turbulent drag coefficient would show the turbulent boundary
layer generates a friction which is much larger than the laminar boundary layer.

8.3.3 Exercise # 3 : simplified model of a turbulent jet

Question
The Reynolds averaged equations take the form :

∂U i
=0



 ∂xi

(8.20)
∂ ∂  0 0  1 ∂p ∂2U i



U iU j + ui uj + =ν



∂xj ∂xj ρ ∂xi ∂xi ∂xi

The viscous stress cannot be neglected with respect to the turbulent stress in the so-called viscous
sublayer region, very close to a solid wall. In the case of a free flow, such as a turbulent jet
created from a slot, the viscous stress can therefore be neglected with respect to the turbulent
stress everywhere in the jet. Moreover, the statistical flow invariance along the x3 direction
allows to further simplify (8.20) into :

∂U 1 ∂U 2
+ =0
∂x1 ∂x2
2
∂(U 1 ) ∂(U 1 U 2 ) ∂  0 0 ∂  0 0  1 ∂p
+ + u1 u1 + u u + =0 (8.21)
∂x1 ∂x2 ∂x1 ∂x2 1 2 ρ ∂x1
2
∂(U 2 U 1 ) ∂(U 2 ) ∂  0 0 ∂  0 0  1 ∂p
+ + u2 u1 + u u + =0
∂x1 ∂x2 ∂x1 ∂x2 2 2 ρ ∂x2

It is observed the variation of the statistical flow properties (mean values and Reynolds stress)
in the turbulent jet is much faster in the transverse direction x2 than in the axial direction x1 .
This leads to introduce d the slot height (in the x2 direction) as characteristic length for the
flow variation along x2 and d/ with  << 1 as the much larger characteristic length for the flow
variation along x1 . Let us also denote V1 the order of magnitude of the mean axial velocity U 1 ,
V2 the order of magnitude of the mean transverse velocity U 2 and v 0 the order of magnitude of
the turbulent velocity fluctuation (be it u01 or u02 ). Let us finally denote P the order of magnitude
of the mean pressure p. It will be assumed no mean axial pressure gradient exists outside the
jet region.
√
• Show that V2 = V1 , v 0 = V1 , P = ρV12 and establish that system (8.21) can be simplified
into :
∂U 1 ∂U 2
+ =0
∂x1 ∂x2
(8.22)
2
∂(U 1 ) ∂(U 1 U 2 ) ∂(u01 u02 )
+ =−
∂x1 ∂x2 ∂x2
 394 Turbulent boundary layer

Answer
√
 Let us show that V2 = V1 , v 0 = V1 , P = ρV12 and let us then establish that system (8.21)
can be simplified into (8.22).
Inserting the orders of magnitude for U 1 , U 2 , x1 and x2 into the continuity equation yields :

∂U 1 ∂U 2
+ =0
∂x ∂x
| {z1} | {z2}
=O(V1 /(d/)) =O(V2 /d)

Both terms must have the same order of magnitude in order to properly express the continuity
equation, that is the fact the divergence of the velocity field is zero, so that :
V1 V2
=
d d
which simplifies into V2 = V1 indicating the order of magnitude of the transverse velocity in
the jet is much smaller than the order of magnitude of the axial velocity in the jet.

We introduce in the following x1 -momentum equation the non-dimensional flow quantities with-
out changing their notations and indicate in blue the order of magnitude for each term coming
from the non-dimensionalization process :
2
V12 ∂(U 1 ) V12 ∂(U 1 U 2 ) (v 0 )2 ∂(u01 u01 ) (v 0 )2 ∂(u01 u02 ) P ∂p
+ + + + =0
d ∂x1 d ∂x2 d ∂x1 d ∂x2 ρd ∂x1

∂(u01 u01 ) ∂(u01 u02 )

Since  << 1, the term can be neglected with respect to the term .
∂x1 ∂x2
Moreover, since the turbulent stress must be both taken into account into the x1 -momentum
conservation equation but not be the only remaining term in the equation, it follows that
V12 (v 0 )2
=
d d
√
which simplifies into v 0 = V1 , indicating the turbulent velocity fluctuation is small with
respect to the mean (axial) velocity.
Finally, the axial pressure gradient is retained in the equation if P = ρV12 .
This leads to the following simplified x1 -momentum equation :
2
∂(U 1 ) ∂(U 1 U 2 ) 1 ∂p ∂(u01 u02 )
+ =− −
∂x1 ∂x2 ρ ∂x1 ∂x2

The same analysis is applied to the x2 -momentum equation :

2
2 V12 ∂(U 2 U 1 ) 2 V12 ∂(U 2 ) (v 0 )2 ∂(u02 u01 ) (v 0 )2 ∂(u02 u02 ) V12 ∂p
+ + + + =0
d ∂x1 d ∂x2 d ∂x1 d ∂x2 d ∂x2
or else, using (v 0 )2 = V12 :
2
2 V12 ∂(U 2 U 1 ) 2 V12 ∂(U 2 ) 2 V12 ∂(u02 u01 ) V12 ∂(u02 u02 ) V12 ∂p
+ + + + =0
d ∂x1 d ∂x2 d ∂x1 d ∂x2 d ∂x2
8.3. Exercises and problem 395

which shows that all the terms in the x2 -momentum equation can be neglected with respect
to the last term corresponding to the transverse pressure gradient. The simplified form of the
x2 -momentum equation is therefore :
∂p
=0
∂x2
Since it is also assumed there is no mean axial pressure gradient outside the jet region, it can
be concluded that pressure is uniform in the flow so that the pressure gradient term in the
x1 -momentum equation vanishes and the simplified governing equations of the turbulent plane
jet read :
∂U 1 ∂U 2
+ =0
∂x1 ∂x2
2
∂(U 1 ) ∂(U 1 U 2 ) ∂(u01 u02 )
+ =−
∂x1 ∂x2 ∂x2
corresponding indeed to (8.22). 

8.3.4 Problem # 1 : total drag of an airfoil

Questions
The pressure distribution is measured along a symmetric airfoil of chord c at zero degree incidence
in a high-Reynolds number flow. The flow is supposed incompressible and homogeneous (ρ =
constant). The far-field upstream pressure and velocity are denoted respectively p∞ , U∞ . The
objective of the problem is to show how the total drag of the airfoil can be estimated from this
measured pressure distribution. It can be proved - and will be admitted here - that the drag
coefficient of the airfoil can be computed using the following formula :
ˆ
D 2 u u
CD = = (1 − ) dy (8.23)
q∞ c c U
(Σ) ∞ U∞

where D denotes the total drag force, q∞ = ρU∞ 2 /2 is the dynamic pressure of the incoming

flow and the integration line (Σ) is a vertical straightline located downstream of the airfoil,
sufficiently far away from the airfoil for the velocity outside the wake to recover the upstream
value U∞ (see Fig.8.9).

Let us introduce the momentum thickness θs of the wake developing downstream of the airfoil.
The thickness θs is defined as :
ˆ
u u
θs = (1 − ) dy
(sillage) U∞ U∞

It can be immediately deduced from (8.23) that :

θ∞
CD = 2 (8.24)
c
 396 Turbulent boundary layer

ue
ue

ue = U

8
θ1 θ2 θ(x )
CL θbf θ

8
U sillage
8

0 Xc c x

bord d’attaque bord de fuite ( Σ)

Figure 8.9: Overview of the airfoil configuration. ”bord d’attaque” = leading edge. ”bord de
fuite” = trailing edge. ”sillage” = wake. The outer inviscid velocity profile is reported along
the extrados - the same profile is observed along the intrados since the airfoil is symmetric and
at zero incidence.

where θ∞ denotes the value of the wake momentum thickness in the downstream section (Σ).

1) The equation governing the evolution of θs in the wake of the airfoil is similar to the von
Kármán equation for the boundary layer, without naturally the wall friction term since the wake
flow is a free flow :
dθs θs due
+ (H + 2) =0 (8.25)
dx ue dx
where ue denotes in that case the inviscid velocity distribution outside the wake downstream of
the airfoil (whereas ue will denote the inviscid velocity distribution outside the boundary layer
along the airfoil).
• Integrate (8.25) from the trailing ledge of the airfoil (where the initial momentum thickness
in the wake will be denoted θbf and where the velocity ue is equal to Vbf × U∞ ) and the section
(Σ) in the farfield wake (where ue = U∞ ) assuming the shape factor constant, equal to a value
Hs which will be later specified.
• Deduce the expression of the total drag coefficient as a function of θbf , Vbf , Hs and c.

2) • Explain why the pressure measured along the airfoil surface can be also considered as
the pressure pe outside the boundary layer developing along the extrados and intrados of the
airfoil from the leading edge x = 0.
• Show that the pressure coefficient :
pe − p∞
Cp =
q∞
is such that :
qe
Cp = 1 − (8.26)
q∞
8.3. Exercises and problem 397

where the dynamic pressure qe is defined as qe = 12 ρu2e .

Expression (8.26) can be used to plot, from experimental values of the pressure, the distri-
bution qe /q∞ or ue /U∞ associated to the airfoil. Let us introduce the non-dimensional outer
ue x
velocity V (X) = and the non-dimensional position along the airfoil X = .
U∞ c
Two zones are considered along the upper part of the airfoil : first the interval [0, X̄] along which
the outer inviscid velocity is constant, V (X) = V (X̄) = V̄ and next the interval [X̄, 1] along
which the inviscid flow slows down from V̄ to Vbf . The non-dimensional outer inviscid velocity
is therefore such that : 
 V̄ , for X ∈ [0, X̄]
V (X) =
AX + B, for X ∈ [X̄, 1]


Vbf − V̄ V̄ − Vbf X̄
with A = < 0 and B = .
1 − X̄ 1 − X̄
Let us denote θ1 = θ(X̄) lthe momentum thickness at X = X̄ and θ2 the momentum thickness at
the airfoil trailing edge. The flow being symmetric, the momentum thickness distribution along
the intrados (lower surface) of the airfoil is the same as the momentum thickness distribution
along the extrados (upper surface) so that θbf = 2θ2 . The estimation of the airfoil total drag
requires therefore the calculation of the momentum thickness θ2 .

3) • Use the simple integral method described in 7.3.4 to obtain :

 m+1
θ1 b(m + 1) X̄
= (8.27)
c V̄ m (Re∞ )m

U∞ c
where the notations b, m are those of 7.3.4 and Re∞ = .
ν
It will be admitted that the application of the same simple integral method over the interval
[X̄, X] with X ≤ 1 yields :
 m+1  m+1  α  β
θ(X) θ1 V̄ b(m + 1) 1 V̄
= + (X − X̄) (V (X))1−m (1 − )
c c V (X) β(V (X) − V̄ ) (Re∞ )m V (X)
(8.28)
where α = (H + 2)(m + 1) and β = α − m + 1. Applying formula (8.28) with X = 1 provides θ2 .

4) Let us first assume that the boundary layer developing along the airfoil from the leading edge
X = 0 is laminar all along the airfoil with the outer inviscid flow such that V (X) = V̄ = 1.3
from X = 0 to X̄ = 0.5 and V (X) linearly decreasing from V̄ to Vbf = 0.5 between X̄ and X = 1.

• Compute the total drag of the airfoil in that case.

Note : formulae (8.27) and (8.28) will be applied with the standard values of the parameters
b, m and H given in 7.3.4 for the laminar case, U∞ = 60 m/s−1 , c = 1 m, ν = 1.5 × 10−5 m2 /s
and ρ = 1.2 kg/m3 . The averaged shape factor Hs in the wake will be computed as the average
between its value at the trailing edge, Hs = 2.59, and its value in the far-field wake Hs = 1.
 398 Turbulent boundary layer

θ2 due
5) The laminar boundary layer separates when the non-dimensional parameter − be-
ν dx
comes larger than a given value C (C = 0.068 if the self-similar Falkner-Skan solution is used).

• Apply the above separation criterion to show that the laminar boundary layer along the
airfoil separates at a non-dimensional position Xd which does not depend on Re∞ - it is not
asked to explicitely compute Xd -.

6) In practice, for the flow under study, the boundary layer becomes turbulente before the
separation can take place and this turbulent boundary layer remains attached to the airfoil until
the trailing edge.

• Show, using the Michel criterion given in class, that the transition occurs for a non-dimensional
position Xt less than X̄.

Note : the formula (8.27) will be adapted to estimate θ(X) with X between 0 and X̄.

7) The end of the analysis is simplified by supposing Xt = X̄ so that the boundary layer is
laminar from X = 0 to X = X̄ and turbulent from X = X̄ to the trailing edge X = 1.

• Compute in that case the total drag of the airfoil and comment the result with respect to
the one obtained in 4).

Note : the standard values of m, b and H given in 8.1.3 will be used when applying the simple
integral method in the turbulent case. The value Hs will be computed as the average between
its value at the trailing edge, Hs = 1.4, and its value in the far-field wake, Hs = 1. The fluid
properties and the upstream conditions are unchanged with respect to 4).

Answers
1)  The following equation
dθs θs due
+ (H + 2) =0
dx ue dx
is integrated between the airfoil trailing edge and the section (Σ) in the far-field wake. H is
assumed constant, equal to Hs . The above equation can be immediately rewritten as :
θs0 (x) u0 (x)
= −(Hs + 2) e
θs (x) ue (x)
It is then easy to integrate it between the trailing edge where θs = θbf , ue = Vbf × U∞ and the
section (Σ) where θs = θ∞ and ue = U∞ :
ˆ (Σ) 0 ˆ (Σ) 0
θs (x) ue (x)
dx = −(Hs + 2) dx
(bf ) θs (x) (bf ) ue (x)

(Σ) (Σ)
⇔ [ln(θs (x))](bf ) = −(Hs + 2) × [ln(ue )](bf )

θs (Σ) ue (Σ)
⇔ ln( ) = −(Hs + 2) × ln( )
θs (bf ) ue (bf )
8.3. Exercises and problem 399

Since θs (Σ) = θ∞ , θs (bf ) = θbf and ue (Σ) = U∞ , ue (bf ) = Vbf × U∞ , it comes :

θ∞ U∞
ln( ) = −(Hs + 2) × ln( ) = (Hs + 2)ln(Vbf )
θbf Vbf × U∞

Taking the exponential of the previous expression yields :

(Hs +2)
θ∞ = θbf Vbf

Injecting this expression of θ∞ into (8.24) yields :

θbf (Hs +2)

CD = 2 V (8.29)
c bf

which expresses the total airfoil drag as a function of the trailing edge wake momentum thickness
θbf , the airfoil chord c, the shape factor Hs and the trailing edge velocity Vbf . Formula (8.29)
is due to Squire and Young. For a symmetric airfoil at zero incidence, θbf = 2θ2 where θ2 is the
boundary layer momentum thickness at the trailing edge at the airfoil extrados (also equal to
its value at the airfoil intrados). Formula (8.29) can therefore also be expressed as :

θ2 (Hs +2)
CD = 4 V (8.30)
c bf

It remains to estimate θ2 using the VKE. 

2)  In the boundary layer the pressure does not vary in the direction normal to the wall.
Therefore, the pressure measured at the wall at a location x = x0 along the airfoil is also equal
to the outer inviscid pressure at the same location outisde the boundary layer, that is pe (x0 ).
Applying the Bernoulli theorem in the inviscid flow region outside the boundary layer yields :

1 1
pe + ρu2e = p∞ + ρu2∞
2 2
or else pe + qe = p∞ + q∞ so that (pe − p∞ )/q∞ = 1 − qe /q∞ and finally Cp = 1 − qq∞e .
The outer inviscid velocity distribution ue can be therefore known from wall pressure measure-
ments. Such a distribution will be next given and the simple integral method will be applied in
order to compute θ2 and deduce the airfoil drag. 

3) • Since ue is constant for x between 0 and c × X̄, the simplified form of the simple integral
method can be used, which states that for a boundary layer with zero outer pressure gradient
(that is also a constant outer inviscid velocity) the momentum thickness evolves from x = 0 as :
ˆ x
m+1 1
θ(x) = b(m + 1) m
dt
0 (ue /ν)

where b and m are constant coefficients. In the present case one takes x = X̄ ×c where θ(x) = θ1
so that :
ˆ X̄×c
m+1 1 X̄c
θ1 = b(m + 1) m
dt = b(m + 1) .
0 (U∞ V̄ /ν) (U∞ V̄ /ν)m
 400 Turbulent boundary layer

Dividing by cm+1 and introducing the Reynolds number Re∞ = U∞ c/ν yields (8.27) that is :

θ1 m+1 b(m + 1) X̄
( ) =
c V̄ m (Re∞ )m

Formula (8.28) is admitted - it results from the same simple integral method applied between
X̄ and X. Using this formula with X = 1 (θ(1) = θ2 ,V (1) = Vbf ) yields :
 m+1  m+1  α  β
θ2 θ1 V̄ b(m + 1) 1 V̄
= + (1 − X̄) (Vbf )1−m (1 − ) (8.31)
c c Vbf β(Vbf − V̄ ) (Re∞ )m Vbf

Applying formula (8.27) yields θ1 ; once θ1 is known, formula (8.31) can be applied to obtain θ2 . 

4)  The boundary layer developing along the extrados is assumed laminar for the leading edge
to the trailing edge. The standard values m = 1, b = 0.2205, H = 2.591 are used to compute
the momentum thickness. The outer inviscid distribution is such that V̄ = 1.3, X̄ = 0.5 and
Vbf = 0.5.
A first immediate calculation yields Re∞ = U∞ ν
c
= 4 × 106 . It is unlikely the laminar regime
will be preserved all along the airfoil for such a large value of the Reynolds number.
A second immediate calculation yields (θ1 /c)2 = 0.1696/Re∞ = 4.24 × 10−8 so that θ1 =
2.06 × 10−4 m.
Applying (8.31) to compute θ2 is also straightforward. Since α = 9.182, β = α = 9.182, it
comes :
θ2 2 θ1 1.3 0.441 1 1.3
( ) = ( )2 × ( )9.182 + × 0.5 × × (1 × (1 − ( )9.182 ))
c c 0.5 9.182 × (0.5 − 1.3) Re∞ 0.5
1095.749 193.91
so that (θ2 /c)2 = + = 3.224 × 10−4 or else θ2 = 1.7955 × 10−2 m.
Re∞ Re∞
The total drag is obtained from D = q∞ cCD with q∞ = 21 ρ∞ U∞ 2 and C given by (8.30). In
D
the CD formula, θ2 /c has just been computed and Vbf = 0.5. Also Hs = 12 (2.59 + 1) = 1.795. It
follows that CD = 5.17 × 10−3 and D = 12 × 1.2 × 602 × 1 × 0.00517 = 11.17 N .

5) The separation of the laminar boundary layer is assumed to take place when :

θ2 due
− ≥C
ν dx
with C = 0.068. The boundary layer is likely to separate in a flow region where the pressure
gradient is adverse, that is where it is positive along the flow direction x, also corresponding to
a decrease of the flow velocity. For the flow under study, this means separation is likely to take
place for X ∈ [X̄, 1]. Combining formula (8.31) with (8.27) yields (m = 1 in the laminar case) :

θ(X) 2 C(X)
( ) = . (8.32)
c Re∞
Using non-dimensional variables, the separation criterion reads :

θ2 due θ c2 U∞ dV θ
− = −( )2 × × = −( )2 V 0 (X)Re∞ ≥ C
ν dx c ν c dX c
8.3. Exercises and problem 401

or else, using (8.32), −C(X)V 0 (X) ≥ C. The separation location Xd given by this criterion is
therefore such that C(Xd )V 0 (Xd ) = −C and does not depend on the Reynolds number Re∞ .

6) Michel’s criterion estimates the transition location xt as :

Reθt = 2.9(Rext )0.4 (8.33)

where Reθt = ue θ(xt )/ν and Rext = ue xt /ν. It is suggested to look for Xt = xt /c in the airfoil
region between 0 and X̄. Therefore, θ(X) can be derived from :
ˆ x
1
θ(x)m+1 = b(m + 1) m
dt
0 (ue /ν)
or else ˆ X
θ b(m + 1) dT b(m + 1) X
( )m+1 = =
c Rem
∞ 0 V m (T ) Rem
∞ V̄ m
since V (X) = V̄ for X ∈ [0, X̄]. With m = 1 and b = 0.2205, it comes (θ(X)/c)2 = 0.441 ×
X/(Re∞ × V̄ ). The criterion (8.33) can be rewritten by introducing Re∞ and V̄ to obtain :
θt
(Re∞ V̄ )0.6 × () = 2.9Xt0.4
c
Replacing θt /c with its expression as a function of Xt yields :
√
(Re∞ V̄ )0.1 0.441Xt0.5 = 2.9Xt0.4
or else √
(2.9/ 0.441)10
Xt =
(Re∞ V̄ )
and eventually Xt = 0.485. Since the Reynolds number based on the chord c is approximately
4 × 106 , the transition Reynolds number is around 2 × 106 which corresponds to a realistic value
for the critical Reynolds number in an external flow. 

7)  We still have of course θ1 /c = 2.06 × 10−4 but the momentum thickness θ2 must be
estimated from :
 m+1  m+1  α
θ2 θ1 V̄
=
c c Vbf
 β
b(m + 1) 1 1−m V̄
+ (1 − X̄) (Vbf ) (1 − )
β(Vbf − V̄ ) (Re∞ )m Vbf
with m = 0.2, b = 0.0086, H = 1.4 so that α = 4.08 and β = 4.88.
A careful calculation yields :
 1.2  1.2
θ2 θ1 0.079665
= 49.328 + 0.2
c c Re∞
so that θ2 /c = 0.0159. The total drag coefficient follows from :
θ2 (Hs +2)
CD = 4 V
c bf
with Hs = 1.2 hence CD = 0.00692. The total drag is therefore equal to D = 21 × 1.2 × 602 ×
1 × 0.00692 = 14.95 N . This value is 30% larger than the drag estimated under the (unrealistic)
assumption of a fully laminar flow. 
402 Turbulent boundary layer
Bibliography

[1] P.M. Gerhart, A.L. Gerhart, and J.I. Hochstein. Munson, Young and Okiishi’s Fundamentals
of Fluid Mechanics, 8th Edition. Wiley Editors, 2015.

[2] P.K. Kundu, I.M. Cohen, and D.R. Dowling. Fluid Mechanics, 6th Edition. Academic Press,
2016.

[3] Milton Van Dyke. An Album of Fluid Motion, 4th Edition. The Parabolic Press, 1988.
404 BIBLIOGRAPHY
Useful tables

Several useful tables and abaqus for compressible flow analysis are provided hereafter :

I) a table of shock relationships.

For a given value of the upstream normal Mach number (this normal Mach number (Mn )1 is
equal to the upstream Mach number M1 in the case of a 1D or straight or normal shock), the
values of the static pressure jump through the shock, of the total pressure loss and of the down-
stream Mach number M2 are provided. The value of the downstream Mach number is valid only
for a 1D shock. All the values are computed for γ = 1.4.
The static pressure jump is computed from (4.19).
The downstream Mach number for a normal or 1D shock is computed from (4.17).
The loss of total pressure is computed from (4.21).

II) a table for an isentropic flow in a varying-section nozzle providing the static to
total pressure ratio (p/p0 ) and the local aera to sonic or critical area ratio (A/A∗ ) for a given
Mach number. These ratios are computed from the fundamental relationship (5.11), valid for
isentropic flows. Keep in mind that two possible entries exist for the Mach number when con-
sidering a given value for the ratio A/A∗ : a subsonic value and a supersonic value. Selecting
one value of the other depends on the detailed analysis of the nozzle flow, allowing to identity
which flow regime (subsonic / supersonic) can take place in the flow region under analysis.

III) an abacus for the relationship θ − β − M (6.5) :

(γ − 1)M12 sin2 (β) + 2
tan(β − θ) = tan(β)[ ]
(γ + 1)M12 sin2 (β)
The curves β = β(θ, M1 ) are plotted for fixed values of the upstream Mach number M1 .

IV) a table of values for the Prandtl-Meyer function.

For a given value of the Prandtl-Meyer function ν(M ) :
ˆ √ 2 r r
M − 1 dM γ+1 −1 γ−1 p
ν(M ) = = tan ( (M 2 − 1)) − tan−1 ( M 2 − 1)
γ−1 2 M γ−1 γ+1
1+ M
2
Table IV) provides the value of the corresponding Mach number M and of the corresponding
Mach angle sin−1 ( M1
).
 406 Useful tables

p2 (p0 )2 p2 (p0 )2
(M1 )n M2 (M1 )n M2
p1 (p0 )1 p1 (p0 )1
(1D shock) (1D shock)
1.00 1.0000 1.0000 1.0000 1.36 1.9912 0.9676 0.7572
1.01 1.0235 1.0000 0.9901 1.37 2.0230 0.9653 0.7527
1.02 1.0471 1.0000 0.9805 1.38 2.0551 0.9630 0.7483
1.03 1.0710 1.0000 0.9712 1.39 2.0874 0.9607 0.7440
1.04 1.0952 0.9999 0.9620 1.40 2.1200 0.9582 0.7397
1.05 1.1196 0.9999 0.9531 1.41 2.1528 0.9557 0.7355
1.06 1.1442 0.9998 0.9444 1.42 2.1858 0.9531 0.7314
1.07 1.1691 0.9996 0.9360 1.43 2.2190 0.9504 0.7274
1.08 1.1941 0.9994 0.9277 1.44 2.2525 0.9476 0.7235
1.09 1.2195 0.9992 0.9196 1.45 2.2863 0.9448 0.7196
1.10 1.2450 0.9989 0.9118 1.46 2.3202 0.9420 0.7157
1.11 1.2708 0.9986 0.9041 1.47 2.3544 0.9390 0.7120
1.12 1.2968 0.9982 0.8966 1.48 2.3888 0.9360 0.7083
1.13 1.3230 0.9978 0.8892 1.49 2.4234 0.9329 0.7047
1.14 1.3495 0.9973 0.8820 1.50 2.4583 0.9298 0.7011
1.15 1.3763 0.9967 0.8750 1.51 2.4934 0.9266 0.6976
1.16 1.4032 0.9961 0.8682 1.52 2.5288 0.9233 0.6941
1.17 1.4304 0.9953 0.8615 1.53 2.5644 0.9200 0.6907
1.18 1.4578 0.9946 0.8549 1.54 2.6002 0.9166 0.6874
1.19 1.4855 0.9937 0.8485 1.55 2.6363 0.9132 0.6841
1.20 1.5133 0.9928 0.8422 1.56 2.6725 0.9097 0.6809
1.21 1.5415 0.9918 0.8360 1.57 2.7090 0.9062 0.6777
1.22 1.5698 0.9907 0.8300 1.58 2.7458 0.9026 0.6746
1.23 1.5984 0.9896 0.8241 1.59 2.7828 0.8989 0.6715
1.24 1.6272 0.9884 0.8183 1.60 2.8200 0.8952 0.6684
1.25 1.6562 0.9871 0.8126 1.61 2.8575 0.8915 0.6655
1.26 1.6855 0.9857 0.8071 1.62 2.8951 0.8877 0.6625
1.27 1.7150 0.9842 0.8016 1.63 2.9330 0.8838 0.6596
1.28 1.7448 0.9827 0.7963 1.64 2.9712 0.8799 0.6568
1.29 1.7748 0.9811 0.7911 1.65 3.0096 0.8760 0.6540
1.30 1.8050 0.9794 0.7860 1.66 3.0482 0.8720 0.6512
1.31 1.8354 0.9776 0.7809 1.67 3.0870 0.8680 0.6485
1.32 1.8661 0.9758 0.7760 1.68 3.1261 0.8639 0.6458
1.33 1.8970 0.9738 0.7712 1.69 3.1654 0.8599 0.6431
1.34 1.9282 0.9718 0.7664 1.70 3.2050 0.8557 0.6405
1.35 1.9596 0.9697 0.7618 1.71 3.2448 0.8516 0.6380

I - (1) Shock relationships

Useful tables 407

p2 (p0 )2 p2 (p0 )2
(M1 )n M2 (M1 )n M2
p1 (p0 )1 p1 (p0 )1
(1D shock) (1D shock)
1.72 3.2848 0.8474 0.6355 2.08 4.8808 0.6835 0.5643
1.73 3.3250 0.8431 0.6330 2.09 4.9295 0.6789 0.5628
1.74 3.3655 0.8389 0.6305 2.10 4.9783 0.6742 0.5613
1.75 3.4062 0.8346 0.6281 2.11 5.0275 0.6696 0.5598
1.76 3.4472 0.8302 0.6257 2.12 5.0768 0.6649 0.5583
1.77 3.4884 0.8259 0.6234 2.13 5.1264 0.6603 0.5568
1.78 3.5298 0.8215 0.6210 2.14 5.1762 0.6557 0.5554
1.79 3.5714 0.8171 0.6188 2.15 5.2263 0.6511 0.5540
1.80 3.6133 0.8127 0.6165 2.16 5.2765 0.6464 0.5525
1.81 3.6554 0.8082 0.6143 2.17 5.3271 0.6419 0.5511
1.82 3.6978 0.8038 0.6121 2.18 5.3778 0.6373 0.5498
1.83 3.7404 0.7993 0.6099 2.19 5.4288 0.6327 0.5484
1.84 3.7832 0.7948 0.6078 2.20 5.4800 0.6281 0.5471
1.85 3.8262 0.7902 0.6057 2.21 5.5314 0.6236 0.5457
1.86 3.8695 0.7857 0.6036 2.22 5.5831 0.6191 0.5444
1.87 3.9130 0.7811 0.6016 2.23 5.6350 0.6145 0.5431
1.88 3.9568 0.7765 0.5996 2.24 5.6872 0.6100 0.5418
1.89 4.0008 0.7720 0.5976 2.25 5.7396 0.6055 0.5406
1.90 4.0450 0.7674 0.5956 2.26 5.7922 0.6011 0.5393
1.91 4.0894 0.7627 0.5937 2.27 5.8450 0.5966 0.5381
1.92 4.1341 0.7581 0.5918 2.28 5.8981 0.5921 0.5368
1.93 4.1791 0.7535 0.5899 2.29 5.9514 0.5877 0.5356
1.94 4.2242 0.7488 0.5880 2.30 6.0050 0.5833 0.5344
1.95 4.2696 0.7442 0.5862 2.31 6.0588 0.5789 0.5332
1.96 4.3152 0.7395 0.5844 2.32 6.1128 0.5745 0.5321
1.97 4.3611 0.7349 0.5826 2.33 6.1670 0.5702 0.5309
1.98 4.4071 0.7302 0.5808 2.34 6.2215 0.5658 0.5297
1.99 4.4534 0.7255 0.5791 2.35 6.2762 0.5615 0.5286
2.00 4.5000 0.7209 0.5774 2.36 6.3312 0.5572 0.5275
2.01 4.5468 0.7162 0.5757 2.37 6.3864 0.5529 0.5264
2.02 4.5938 0.7115 0.5740 2.38 6.4418 0.5486 0.5253
2.03 4.6410 0.7069 0.5723 2.39 6.4974 0.5444 0.5242
2.04 4.6885 0.7022 0.5707 2.40 6.5533 0.5401 0.5231
2.05 4.7362 0.6975 0.5691 2.41 6.6094 0.5359 0.5221
2.06 4.7842 0.6928 0.5675 2.42 6.6658 0.5317 0.5210
2.07 4.8324 0.6882 0.5659 2.43 6.7224 0.5276 0.5200

I - (2) Shock relationships

 408 Useful tables

p2 (p0 )2 p2 (p0 )2
(M1 )n M2 (M1 )n M2
p1 (p0 )1 p1 (p0 )1
(1D shock) (1D shock)
2.44 6.7792 0.5234 0.5189 2.80 8.9800 0.3895 0.4882
2.45 6.8362 0.5193 0.5179 2.81 9.0454 0.3862 0.4875
2.46 6.8935 0.5152 0.5169 2.82 9.1111 0.3829 0.4868
2.47 6.9510 0.5111 0.5159 2.83 9.1770 0.3797 0.4861
2.48 7.0088 0.5071 0.5149 2.84 9.2432 0.3765 0.4854
2.49 7.0668 0.5030 0.5140 2.85 9.3096 0.3733 0.4847
2.50 7.1250 0.4990 0.5130 2.86 9.3762 0.3701 0.4840
2.51 7.1834 0.4950 0.5120 2.87 9.4430 0.3670 0.4833
2.52 7.2421 0.4911 0.5111 2.88 9.5101 0.3639 0.4827
2.53 7.3010 0.4871 0.5102 2.89 9.5774 0.3608 0.4820
2.54 7.3602 0.4832 0.5092 2.90 9.6450 0.3577 0.4814
2.55 7.4196 0.4793 0.5083 2.91 9.7128 0.3547 0.4807
2.56 7.4792 0.4754 0.5074 2.92 9.7808 0.3517 0.4801
2.57 7.5390 0.4715 0.5065 2.93 9.8490 0.3487 0.4795
2.58 7.5991 0.4677 0.5056 2.94 9.9175 0.3457 0.4788
2.59 7.6594 0.4639 0.5047 2.95 9.9862 0.3428 0.4782
2.60 7.7200 0.4601 0.5039 2.96 10.0552 0.3398 0.4776
2.61 7.7808 0.4564 0.5030 2.97 10.1244 0.3369 0.4770
2.62 7.8418 0.4526 0.5022 2.98 10.1938 0.3340 0.4764
2.63 7.9031 0.4489 0.5013 2.99 10.2634 0.3312 0.4758
2.64 7.9645 0.4452 0.5005 3.00 10.3333 0.3283 0.4752
2.65 8.0262 0.4416 0.4996 3.10 11.0450 0.3012 0.4695
2.66 8.0882 0.4379 0.4988 3.20 11.7800 0.2762 0.4643
2.67 8.1504 0.4343 0.4980 3.30 12.5383 0.2533 0.4596
2.68 8.2128 0.4307 0.4972 3.40 13.3200 0.2322 0.4552
2.69 8.2754 0.4271 0.4964 3.50 14.1250 0.2129 0.4512
2.70 8.3383 0.4236 0.4956 3.60 14.9533 0.1953 0.4474
2.71 8.4014 0.4201 0.4949 3.70 15.8050 0.1792 0.4439
2.72 8.4648 0.4166 0.4941 3.80 16.6800 0.1645 0.4407
2.73 8.5284 0.4131 0.4933 3.90 17.5783 0.1510 0.4377
2.74 8.5922 0.4097 0.4926 4.00 18.5000 0.1388 0.4350
2.75 8.6562 0.4062 0.4918 5.00 29.0000 0.0617 0.4152
2.76 8.7205 0.4028 0.4911
2.77 8.7850 0.3994 0.4903
2.78 8.8498 0.3961 0.4896
2.79 8.9148 0.3928 0.4889

I - (3) Shock relationships

Useful tables 409

M p/p0 A/A∗ M p/p0 A/A∗

0.01 0.999930 57.873840 0.41 0.890711 1.558673
0.02 0.999720 28.942123 0.42 0.885722 1.528905
0.03 0.999370 19.300543 0.43 0.880651 1.500718
0.04 0.998881 14.481482 0.44 0.875498 1.474005
0.05 0.998252 11.591443 0.45 0.870267 1.448672
0.06 0.997484 9.665911 0.46 0.864960 1.424629
0.07 0.996577 8.291526 0.47 0.859580 1.401795
0.08 0.995533 7.261608 0.48 0.854128 1.380097
0.09 0.994350 6.461342 0.49 0.848607 1.359468
0.10 0.993032 5.821829 0.50 0.843019 1.339843
0.11 0.991576 5.299228 0.51 0.837367 1.321168
0.12 0.989985 4.864317 0.52 0.831654 1.303388
0.13 0.988260 4.496858 0.53 0.825881 1.286454
0.14 0.986400 4.182399 0.54 0.820050 1.270321
0.15 0.984408 3.910343 0.55 0.814165 1.254948
0.16 0.982284 3.672738 0.56 0.808228 1.240295
0.17 0.980030 3.463508 0.57 0.802241 1.226326
0.18 0.977647 3.277926 0.58 0.796206 1.213007
0.19 0.975135 3.112258 0.59 0.790127 1.200308
0.20 0.972497 2.963520 0.60 0.784004 1.188200
0.21 0.969733 2.829294 0.61 0.777841 1.176654
0.22 0.966845 2.707602 0.62 0.771640 1.165645
0.23 0.963835 2.596811 0.63 0.765402 1.155151
0.24 0.960703 2.495563 0.64 0.759131 1.145148
0.25 0.957453 2.402710 0.65 0.752829 1.135616
0.26 0.954085 2.317287 0.66 0.746498 1.126535
0.27 0.950600 2.238470 0.67 0.740140 1.117887
0.28 0.947002 2.165553 0.68 0.733758 1.109654
0.29 0.943291 2.097927 0.69 0.727353 1.101822
0.30 0.939470 2.035065 0.70 0.720928 1.094373
0.31 0.935540 1.976507 0.71 0.714485 1.087293
0.32 0.931503 1.921851 0.72 0.708025 1.080571
0.33 0.927361 1.870745 0.73 0.701552 1.074192
0.34 0.923117 1.822876 0.74 0.695068 1.068144
0.35 0.918773 1.777969 0.75 0.688573 1.062417
0.36 0.914329 1.735778 0.76 0.682070 1.056999
0.37 0.909790 1.696086 0.77 0.675562 1.051881
0.38 0.905156 1.658696 0.78 0.669050 1.047053
0.39 0.900430 1.623433 0.79 0.662536 1.042505
0.40 0.895615 1.590140 0.80 0.656022 1.038230

II - (1) Isentropic flow with varying section (γ = 1.4)

410 Useful tables

M p/p0 A/A∗ M p/p0 A/A∗

0.81 0.649509 1.034219 1.21 0.407021 1.033439
0.82 0.643000 1.030464 1.22 0.401711 1.036572
0.83 0.636496 1.026959 1.23 0.396446 1.039835
0.84 0.630000 1.023696 1.24 0.391229 1.043229
0.85 0.623512 1.020669 1.25 0.386058 1.046753
0.86 0.617034 1.017871 1.26 0.380934 1.050406
0.87 0.610569 1.015297 1.27 0.375857 1.054189
0.88 0.604117 1.012941 1.28 0.370828 1.058100
0.89 0.597680 1.010798 1.29 0.365847 1.062138
0.90 0.591260 1.008863 1.30 0.360914 1.066304
0.91 0.584858 1.007131 1.31 0.356029 1.070598
0.92 0.578476 1.005597 1.32 0.351192 1.075018
0.93 0.572114 1.004257 1.33 0.346403 1.079565
0.94 0.565775 1.003108 1.34 0.341663 1.084238
0.95 0.559460 1.002145 1.35 0.336971 1.089038
0.96 0.553170 1.001364 1.36 0.332328 1.093964
0.97 0.546905 1.000763 1.37 0.327733 1.099015
0.98 0.540668 1.000337 1.38 0.323187 1.104193
0.99 0.534460 1.000084 1.39 0.318690 1.109496
1.00 0.528282 1.000000 1.40 0.314241 1.114926
1.01 0.522134 1.000083 1.41 0.309840 1.120481
1.02 0.516018 1.000330 1.42 0.305488 1.126162
1.03 0.509935 1.000738 1.43 0.301185 1.131969
1.04 0.503886 1.001305 1.44 0.296929 1.137903
1.05 0.497872 1.002029 1.45 0.292722 1.143963
1.06 0.491894 1.002907 1.46 0.288563 1.150149
1.07 0.485952 1.003938 1.47 0.284452 1.156463
1.08 0.480047 1.005119 1.48 0.280388 1.162903
1.09 0.474181 1.006449 1.49 0.276372 1.169471
1.10 0.468354 1.007925 1.50 0.272403 1.176167
1.11 0.462567 1.009547 1.51 0.268481 1.182991
1.12 0.456820 1.011312 1.52 0.264607 1.189943
1.13 0.451115 1.013219 1.53 0.260779 1.197023
1.14 0.445451 1.015267 1.54 0.256997 1.204234
1.15 0.439829 1.017454 1.55 0.253262 1.211574
1.16 0.434251 1.019779 1.56 0.249573 1.219044
1.17 0.428716 1.022242 1.57 0.245930 1.226644
1.18 0.423225 1.024840 1.58 0.242332 1.234376
1.19 0.417779 1.027573 1.59 0.238779 1.242239
1.20 0.412377 1.030439 1.60 0.235271 1.250235

II - (2) Isentropic flow with varying section (γ = 1.4)

Useful tables 411

M p/p0 A/A∗ M p/p0 A/A∗

1.61 0.231808 1.258364 2.02 0.123888 1.715971
1.62 0.228389 1.266625 2.04 0.120087 1.745139
1.63 0.225014 1.275021 2.06 0.116399 1.775016
1.64 0.221683 1.283553 2.08 0.112823 1.805614
1.65 0.218395 1.292219 2.10 0.109353 1.836943
1.66 0.215150 1.301021 2.12 0.105988 1.869015
1.67 0.211948 1.309960 2.14 0.102726 1.901843
1.68 0.208788 1.319037 2.16 0.099562 1.935437
1.69 0.205670 1.328252 2.18 0.096495 1.969810
1.70 0.202593 1.337606 2.20 0.093522 2.004974
1.71 0.199559 1.347100 2.22 0.090640 2.040944
1.72 0.196564 1.356735 2.24 0.087846 2.077731
1.73 0.193611 1.366512 2.26 0.085139 2.115348
1.74 0.190697 1.376430 2.28 0.082515 2.153811
1.75 0.187824 1.386492 2.30 0.079973 2.193130
1.76 0.184990 1.396698 2.32 0.077509 2.233323
1.77 0.182195 1.407049 2.34 0.075122 2.274402
1.78 0.179438 1.417546 2.36 0.072810 2.316380
1.79 0.176720 1.428190 2.38 0.070570 2.359275
1.80 0.174040 1.438982 2.40 0.068399 2.403100
1.81 0.171398 1.449923 2.42 0.066297 2.447870
1.82 0.168792 1.461013 2.44 0.064261 2.493602
1.83 0.166224 1.472254 2.46 0.062288 2.540309
1.84 0.163691 1.483648 2.48 0.060378 2.588010
1.85 0.161195 1.495194 2.50 0.058528 2.636719
1.86 0.158734 1.506894 2.52 0.056736 2.686453
1.87 0.156309 1.518749 2.54 0.055000 2.737228
1.88 0.153918 1.530761 2.56 0.053319 2.789063
1.89 0.151562 1.542929 2.58 0.051692 2.841972
1.90 0.149240 1.555256 2.60 0.050115 2.895975
1.91 0.146951 1.567744 2.62 0.048589 2.951088
1.92 0.144696 1.580390 2.64 0.047110 3.007330
1.93 0.142473 1.593200 2.66 0.045679 3.064719
1.94 0.140283 1.606172 2.68 0.044292 3.123274
1.95 0.138126 1.619309 2.70 0.042950 3.183011
1.96 0.135999 1.632611 2.72 0.041650 3.243952
1.97 0.133905 1.646080 2.74 0.040391 3.306113
1.98 0.131841 1.659717 2.76 0.039172 3.369515
1.99 0.129808 1.673523 2.78 0.037992 3.434179
2.00 0.127805 1.687500 2.80 0.036848 3.500123

II - (3) Isentropic flow with varying section (γ = 1.4)

412 Useful tables

M p/p0 A/A∗
2.82 0.035741 3.567368
2.84 0.034669 3.635934
2.86 0.033631 3.705841
2.88 0.032625 3.777113
2.90 0.031651 3.849770
2.92 0.030708 3.923830
2.94 0.029795 3.999319
2.96 0.028910 4.076255
2.98 0.028054 4.154664
3.00 0.027224 4.234569
3.10 0.023449 4.657311
3.20 0.020228 5.120959
3.30 0.017477 5.628648
3.40 0.015125 6.183700
3.50 0.013111 6.789621
3.60 0.011385 7.450110
3.70 0.009903 8.169066
3.80 0.008629 8.950586
3.90 0.007532 9.798974
4.00 0.006586 10.718751
4.10 0.005769 11.714651
4.20 0.005062 12.791640
4.30 0.004449 13.954907
4.40 0.003918 15.209865
4.50 0.003455 16.562197
4.60 0.003053 18.017792
4.70 0.002701 19.582827
4.80 0.002394 21.263721
4.90 0.002126 23.067127
5.00 0.001890 25.000004
5.10 0.001683 27.069584
5.20 0.001501 29.283327
5.30 0.001341 31.649059
5.40 0.001200 34.174816
5.50 0.001075 36.868961
5.60 0.000964 39.740196
5.70 0.000866 42.797436
5.80 0.000779 46.050018
5.90 0.000702 49.507496
6.00 0.000633 53.179802

II - (4) Isentropic flow with varying section (γ = 1.4)

 Useful tables 413

90

80 M1=5.0

70 1.1 1.3 1.5 1.7 2.0 2.5 3.0 4.0

Inclinaison du choc β

60

50

40

30

20

10

0
0 4 8 12 16 20 24 28 32 36 40 44 48
Angle de deflection θ
III - Evolution of the shock angle β (shock angle = angle of inclination for the oblique shock / inclinaison du choc)
as a function of the flow deviation angle θ (angle de déflection)
for various values of the upstream Mach number M1
 414 Useful tables

ν(deg) M µ(deg) ν(deg) M µ(deg) ν(deg) M µ(deg)

0.0 1.000 90.000 17.5 1.689 36.293 35.0 2.329 25.430
0.5 1.051 72.099 18.0 1.706 35.874 35.5 2.349 25.196
1.0 1.082 67.574 18.5 1.724 35.465 36.0 2.369 24.965
1.5 1.108 64.451 19.0 1.741 35.065 36.5 2.390 24.736
2.0 1.133 61.997 19.5 1.758 34.673 37.0 2.410 24.510

2.5 1.155 59.950 20.0 1.775 34.290 37.5 2.431 24.287

3.0 1.177 58.180 20.5 1.792 33.915 38.0 2.452 24.066
3.5 1.198 56.614 21.0 1.810 33.548 38.5 2.473 23.847
4.0 1.218 55.205 21.5 1.827 33.188 39.0 2.495 23.631
4.5 1.237 53.920 22.0 1.844 32.834 39.5 2.516 23.418

5.0 1.256 52.738 22.5 1.862 32.488 40.0 2.538 23.206

5.5 1.275 51.642 23.0 1.879 32.148 40.5 2.560 22.997
6.0 1.294 50.619 23.5 1.897 31.814 41.0 2.582 22.790
6.5 1.312 49.658 24.0 1.915 31.486 41.5 2.604 22.585
7.0 1.330 48.753 24.5 1.932 31.164 42.0 2.626 22.382

7.5 1.348 47.896 25.0 1.950 30.847 42.5 2.649 22.182

8.0 1.366 47.082 25.5 1.968 30.536 43.0 2.671 21.983
8.5 1.383 46.306 26.0 1.986 30.229 43.5 2.694 21.786
9.0 1.400 45.566 26.5 2.004 29.928 44.0 2.718 21.591
9.5 1.418 44.857 27.0 2.023 29.632 44.5 2.741 21.398

10.0 1.435 44.177 27.5 2.041 29.340 45.0 2.764 21.207

10.5 1.452 43.523 28.0 2.059 29.052 45.5 2.788 21.017
11.0 1.469 42.894 28.5 2.078 28.769 46.0 2.812 20.830
11.5 1.486 42.287 29.0 2.096 28.491 46.5 2.836 20.644
12.0 1.503 41.701 29.5 2.115 28.216 47.0 2.861 20.459

12.5 1.520 41.134 30.0 2.134 27.945 47.5 2.886 20.277

13.0 1.537 40.585 30.5 2.153 27.678 48.0 2.910 20.096
13.5 1.554 40.053 31.0 2.172 27.415 48.5 2.936 19.916
14.0 1.571 39.537 31.5 2.191 27.155 49.0 2.961 19.738
14.5 1.588 39.035 32.0 2.210 26.899 49.5 2.987 19.561

15.0 1.605 38.547 32.5 2.230 26.646 50.0 3.013 19.386

15.5 1.622 38.073 33.0 2.249 26.397 50.5 3.039 19.213
16.0 1.639 37.611 33.5 2.269 26.151 51.0 3.065 19.041
16.5 1.655 37.160 34.0 2.289 25.908 51.5 3.092 18.870
17.0 1.672 36.721 34.5 2.309 25.668 52.0 3.119 18.701

IV - (1) Prandlt-Meyer function ν(M ), Mach number M and Mach angle µ

Useful tables 415

ν(deg) M µ(deg) ν(deg) M µ(deg) ν(deg) M µ(deg)

52.5 3.146 18.532 70.0 4.339 13.325 87.5 6.390 9.003
53.0 3.174 18.366 70.5 4.382 13.191 88.0 6.472 8.888
53.5 3.202 18.200 71.0 4.426 13.059 88.5 6.556 8.774
54.0 3.230 18.036 71.5 4.470 12.927 89.0 6.642 8.660
54.5 3.258 17.873 72.0 4.515 12.795 89.5 6.729 8.546

55.0 3.287 17.711 72.5 4.561 12.665 90.0 6.819 8.433

55.5 3.316 17.551 73.0 4.608 12.535 90.5 6.911 8.320
56.0 3.346 17.391 73.5 4.655 12.406 91.0 7.005 8.207
56.5 3.375 17.233 74.0 4.703 12.277 91.5 7.102 8.095
57.0 3.406 17.076 74.5 4.752 12.149 92.0 7.201 7.983

57.5 3.436 16.920 75.0 4.801 12.021 92.5 7.302 7.871

58.0 3.467 16.765 75.5 4.852 11.894 93.0 7.406 7.760
58.5 3.498 16.611 76.0 4.903 11.768 93.5 7.513 7.649
59.0 3.530 16.458 76.5 4.955 11.642 94.0 7.623 7.538
59.5 3.562 16.306 77.0 5.009 11.517 94.5 7.735 7.428

60.0 3.594 16.155 77.5 5.063 11.392 95.0 7.851 7.318

60.5 3.627 16.005 78.0 5.118 11.268 95.5 7.970 7.208
61.0 3.660 15.856 78.5 5.174 11.145 96.0 8.092 7.099
61.5 3.694 15.708 79.0 5.231 11.022 96.5 8.218 6.989
62.0 3.728 15.561 79.5 5.289 10.899 97.0 8.347 6.881

62.5 3.762 15.415 80.0 5.348 10.777 97.5 8.480 6.772

63.0 3.797 15.270 80.5 5.408 10.656 98.0 8.618 6.664
63.5 3.832 15.126 81.0 5.470 10.535 98.5 8.759 6.556
64.0 3.868 14.983 81.5 5.532 10.414 99.0 8.905 6.448
64.5 3.904 14.840 82.0 5.596 10.294 99.5 9.055 6.340

65.0 3.941 14.698 82.5 5.661 10.175 100.0 9.210 6.233

65.5 3.979 14.557 83.0 5.727 10.056 100.5 9.371 6.126
66.0 4.016 14.417 83.5 5.795 9.937 101.0 9.536 6.019
66.5 4.055 14.278 84.0 5.864 9.819 101.5 9.708 5.913
67.0 4.094 14.140 84.5 5.935 9.701 102.0 9.885 5.806

67.5 4.133 14.002 85.0 6.006 9.584

68.0 4.173 13.865 85.5 6.080 9.467
68.5 4.214 13.729 86.0 6.155 9.350
69.0 4.255 13.593 86.5 6.232 9.234
69.5 4.297 13.459 87.0 6.310 9.119

VI - (2) Prandlt-Meyer function ν(M ), Mach number M and Mach angle µ