Chapter 14

The Navier-Stokes equations

14.1 Conservation laws

A general conservation law
We consider the flow of a quantity with density (x, t) at x 2 ⌦ ⇢ Rn ,
with n = 2, 3. For a time t > 0, the total flow of the quantity through the
boundary @⌦, is given by Z
u · n ds, (14.1)
@⌦

where n is the outward unit normal of @⌦, and u = u(x, t) is the velocity
of the flow.
For an arbitrary subdomain ! ⇢ ⌦, the change of the integral of is
equal to the volume source or sink s = s(x, t) minus the total flow of the
quantity through the boundary @!,
Z Z Z
d
(x, t) dx = u · n ds + s(x, t) dx, (14.2)
dt ! @! !

which by Gauss’ theorem leads to

Z ✓ ◆
@
(x, t) + r · ( u) s dx = 0, (14.3)
! @t

for any ! ⇢ ⌦, and thus we get the general continuity equation

˙ + r · ( u) s = 0, (14.4)

for any x 2 ⌦, and t > 0.

129
130 CHAPTER 14. THE NAVIER-STOKES EQUATIONS

Mass conservation
We now consider the flow of mass of a continuum, with ⇢ = ⇢(x, t) the mass
density of the continuum. The general continuity equation with = ⇢, and
zero sink s = 0, gives the equation for conservation of mass

⇢˙ + r · (⇢u) = 0. (14.5)

We say that a flow is incompressible if

r · u = 0, (14.6)

or equivalently if the material derivative is zero,

D⇢
= ⇢˙ + u · r⇢ = 0, (14.7)
Dt
since
D⇢
0 = ⇢˙ + r · (⇢u) = + ⇢r · u. (14.8)
Dt

Conservation of momentum
Newton’s 2nd Law states that the change of momentum ⇢u, is equal to the
sum of all forces, including volume forces,
Z
⇢f dx, (14.9)
!

for a force density f = f (x, t) = (f1 (x, t), ..., fn (x, t)), and surface forces,
Z
n · ds, (14.10)
@!

with the Cauchy stress tensor = ( ij ), where ij = ij (x, t), and we define
n · = ( ji nj ). Gauss’ theorem gives the total force as
Z Z Z
⇢f dx + n · ds = (⇢f + r · ) dx. (14.11)
! @! !

The general continuity equation with = ⇢u, and the sink given by the
sum of all forces, gives the equation for conservation of momentum

@
(⇢u) + r · (⇢u ⌦ u) = ⇢f + r · , (14.12)
@t
14.2. THE NAVIER-STOKES EQUATIONS 131

with u ⌦ u = uuT , the tensor product of the velocity vector field u. With
the help of conservation of mass, we can rewrite the left hand side as

@
(⇢u) + r · (⇢u ⌦ u) = u(⇢˙ + r · (⇢u)) + ⇢(u̇ + (u · r)u) = ⇢(u̇ + (u · r)u),
@t

so that we get
⇢(u̇ + (u · r)u) = ⇢f + r · . (14.13)

The Cauchy stress tensor consists of normal stresses on the diagonal,

and shear stresses on the o↵-diagonal. We can decompose into a dynamic
pressure
1
pd = tr( ), (14.14)
3
and a deviatoric stress tensor ⌧ = + pd I, with I the identity matrix,

= pd I + ⌧, (14.15)

so that
⇢(u̇ + (u · r)u) = ⇢f rpd + r · ⌧. (14.16)

14.2 The Navier-Stokes equations

We now consider incompressible flow, so that the velocity is divergence free,
and we assume the density to be constant. To determine the deviatoric
stress we need a constitutive model of the fluid.
For a Newtonian fluid, the deviatoric stress depends linearly on the
strain rate tensor
✓ ◆
1 T 1 @ui @uj
✏ = (ru + (ru) ) = + , (14.17)
2 2 @xj @xi

with ⌧ = 2µ✏, where µ is the dynamic viscosity.

The incompressible Navier-Stokes equations takes the form,

u̇ + (u · r)u + rp ⌫ u = f, (14.18)
r · u = 0, (14.19)

with the kinematic viscosity ⌫ = µ/⇢, and the kinematic pressure p = pd /⇢.
132 CHAPTER 14. THE NAVIER-STOKES EQUATIONS

Non-dimensionalization
Solutions to the Navier-Stokes equations may take quite di↵erent forms,
depending on the balance of the inertial and dissipative terms of the equa-
tions. To exhibit this balance, we express the Navier-Stokes equations in
terms of the non-dimensional variables u⇤ , p⇤ , f⇤ , x⇤ , t⇤ ,

u = U u⇤ , p = P p⇤ , x = Lx⇤ , f = F f⇤ , t = T t⇤ , (14.20)

where U, P, L, T are characteristic scales of the velocity, pressure, force,

length and time, respectively. The resulting non-dimensionalized di↵eren-
tial operators are scaled as,
@ 1 @ 1 1
= , r= r⇤ , = ⇤, (14.21)
@t T @t⇤ L L2
which gives
U @ U2 P ⌫U
u⇤ + (u⇤ · r⇤ )u⇤ + r⇤ p⇤ ⇤ u⇤ = F f⇤ , (14.22)
T @t⇤ L L L2
U
r · u⇤ = 0, (14.23)
L
or,
1
u̇ + (u · r)u + rp Re u = f, (14.24)
r · u = 0. (14.25)

Here we have dropped the non-dimensional notation for simplicity, with

U2 UL
T = L/U, P = U 2, F = , Re = , (14.26)
L ⌫
where the Reynolds number Re determines the balance between inertial
and viscous characteristics in the flow. For low Re linear viscous e↵ects
dominate, whereas for high Re we have a flow dominated by nonlinear
inertial e↵ect, and turbulence for sufficiently high Reynolds number.
Formally, in the limit Re ! 1, the viscous term vanishes and we are
left with the inviscid Euler equations,

u̇ + (u · r)u + rp = f, (14.27)
r · u = 0, (14.28)

traditionally seen as a model for flow at high Reynolds numbers. Although,

this simple analysis is too naive and relies on strong assumptions on the
14.3. STOKES FLOW 133

regularity of solutions to the equations, an open problem posed as one of

the Clay $1 million Prize problems.
In the limit Re ! 0, we obtain the Stokes equations as a model of
viscous flow,

u + rp = f, (14.29)
r · u = 0, (14.30)

with now a di↵erent scaling of the pressure and the force,

⌫U ⌫U
P = , F = . (14.31)
L L2

14.3 Stokes flow

The Stokes equations
The Stokes equations for a domain ⌦ ⇢ R3 with boundary = D [ N,
and associated normal n, takes the form,

u + rp = f, x 2 ⌦, (14.32)
r · u = 0, x 2 ⌦, (14.33)
u = gD , x 2 D, (14.34)
ru · n + pn = gN , x 2 N. (14.35)

Homogeneous Dirichlet boundary conditions

First assume that we have = D and that gD = 0, that is homogeneous
Dirichlet boundary conditions for the velocity. We then seek a weak solution
to the Stokes equations in the following spaces,

V = [H01 (⌦)]3 (14.36)

Z
Q = {q 2 L2 (⌦) : q dx = 0}, (14.37)
⌦

where the extra condition on Q is needed to assure uniqueness of the pres-

sure, which otherwise is undetermined up to a constant.
We derive the variational formulation by taking the inner product of the
momentum equation with a test function v 2 V , and the inner product of
the continuity equation with a test function q 2 Q. By Green’s formula and
134 CHAPTER 14. THE NAVIER-STOKES EQUATIONS

the homogeneous Dirichlet boundary condition, we obtain the variational

formulation as: find (u, p) 2 V ⇥ Q, such that,
a(u, v) + b(v, p) = (f, v), (14.38)
b(u, q) = 0, (14.39)
for all (v, q) 2 V ⇥ Q, with
Z
a(v, w) = (rv, rw) = rv : rw dx, (14.40)
⌦Z

b(v, q) = (r · v, q) = (r · v)q dx, (14.41)

⌦
and
X 3
@vi @wi
rv : rw = . (14.42)
i,j=1
@x j @x j

The saddle-point problem

The solution (u, p) to the Stokes equations (14.38-14.39), is also the solution
to the constrained minimization problem,
1
min J(v) = a(v, v) (f, v) (14.43)
2
under the constraint
b(v, q) = 0, (14.44)
for which we can formulate the Lagrangian
L(v, q) = J(v) + b(v, q), (14.45)
so that p 2 Q represents a Lagrange multiplier for the constraint r · u = 0.
The Stokes equations thus represent a saddle-point problem, since
L(u, q)  L(u, p)  L(v, p), 8(v, q) 2 V ⇥ Q. (14.46)
Theorem 20. The saddle-point point problem (14.38-14.39) has a unique
solution, if
(i) the bilinear form a is coercive, i.e. that exists an ↵ > 0, such that
a(v, v) ↵kvkV , (14.47)
for all v 2 Z = {v 2 V : b(v, q) = 0, 8q 2 Q},
(ii) the bilinear form b satisfies the inf-sup condition, i.e. there exists a
> 0, such that
b(v, q)
inf sup . (14.48)
q2Q v2V kvkV kqkQ
14.3. STOKES FLOW 135

Mixed finite element approximation

We now formulate a finite element method for solving Stokes equations.
Since we use di↵erent approximation spaces for the velocity and the pres-
sure, we refer to the method as a mixed finite element method. We seek an
approximation (U, P ) 2 Vh ⇥ Qh , such that,

a(U, v) + b(v, P ) = (f, v), (14.49)

b(U, q) = 0, (14.50)

for all (v, q) 2 Vh ⇥ Qh , where Vh and Qh are finite element approxima-

tion spaces. There exists a unique solution to (14.49-14.50), under certain
conditions on the approximation spaces Vh and Qh .

Theorem 21. The mixed finite element problem (14.49-14.50) has a unique
solution (U, P ) 2 Vh ⇥ Qh , if

(i) the bilinear form a is coercive, i.e. that exists an ↵h > 0, such that

a(v, v) ↵h kvkV , (14.51)

for all v 2 Zh = {v 2 Vh : b(v, q) = 0, 8q 2 Qh },

(ii) the bilinear form b satisfies the inf-sup condition, i.e. there exists a
h > 0, such that
b(v, q)
inf sup h, (14.52)
q2Qh v2Vh kvkV kqkQ

and this unique solution satisfies the following error estimate,

✓ ◆
ku U kV + kp P kQ  C inf ku vk + inf kp qk , (14.53)
v2vh q2Qh

for a constant C > 0.

The pair of approximation spaces must be chosen to satisfy the inf-sup

condition, with the velocity space sufficiently rich compared to the pres-
sure space. For example, continuous piecewise quadratic approximation of
the velocity and continuous piecewise linear approximation of the pressure,
referred to as the Taylor-Hood elements. On the other hand, continuous
piecewise linear approximation of both velocity and pressure is not inf-sup
stable.
136 CHAPTER 14. THE NAVIER-STOKES EQUATIONS

Schur complement methods

We seek finite element approximations in the following spaces,
N
X
Vh = {v = (v1 , v2 , v3 ) : vk (x) = vkj j (x), k = 1, 2, 3} (14.54)
j=1

and
M
X
Qh = {q : q(x) = qj j (x)}, (14.55)
j=1

which leads to the following discrete system in matrix form,

  
A B u f
= , (14.56)
BT 0 p 0

with u and p vectors holding the coordinates of U and P in the respective

bases of Vh and Qh .
The matrix A is symmetric positive definite and thus invertible, so we
can express
u = A 1 (f Bp), (14.57)
and since B T u = 0,
B T A 1 Bp = B T A 1 f, (14.58)
which is the Schur complement equation. If null(B) = {0}, then the matrix
S = B T A 1 B is symmetric positive definite and can also be inverted.
Schur complement methods take the form
1
pk = pk 1 C (B T A 1 Bpk 1 B T A 1 f ), (14.59)

where C 1 is a preconditioner for S = B T A 1 B. The Usawa algorithm is

based on C 1 as a scaled identity matrix, which gives
1. Solve Auk = f Bpk 1 ,
2. Set pk = pk 1 + ↵B T uk .

Stabilized methods
Approximation spaces of equal order is possible, by stabilization of the
standard Galerkin finite element method: find (U, P ) 2 Vh ⇥ Qh , such that,

a(U, v) + b(v, P ) = (f, v), (14.60)

b(U, q) + s(P, q) = 0, (14.61)
14.4. THE TRANSIENT NAVIER-STOKES EQUATIONS 137

for all (v, q) 2 Vh ⇥ Qh , where s(P, q) is a pressure stabilization term. The

resulting discrete system takes the form,
  
A B u f
= , (14.62)
BT S p 0
where the stabilization term is chosen so that the matrix S is invertible.
For example, the Brezzi-Pitkäranta stabilization takes the form,
Z
s(P, q) = h2 rP · rq dx. (14.63)
⌦

14.4 The transient Navier-Stokes equations

Semi-discretization
We now formulate a finite element method for solving the unsteady Navier-
Stokes equations (14.18)-(14.19) by semi-discretization.
For each t > 0, we seek approximations (U (t), P (t)) 2 Vh ⇥ Qh , with
U (t) = (U1 (t), U2 (t), U3 (t)), of the form,
N
X M
X
Uk (x, t) = Ukj (t) j (x), k = 1, 2, 3, P (x, t) = P j (t) j (x), (14.64)
j=1 j=1

such that
(U̇ , v) + c(U ; U, v) + a(U, v) + b(v, P ) b(U, q) = (f, v), (14.65)
for all (v, q) 2 Vh ⇥ Qh , where the bilinear forms are defined by (14.40)-
(14.41), with the trilinear form,
Z
c(u; v, w) = ((u · r)v, w) = (u · r)v · w dx. (14.66)
⌦

The semi-discretization (14.65) is a system of ODEs, which we can solve

by a suitable time-stepping method. For r · u = 0, we have that
c(u; v, w) = c̄(u; v, w), (14.67)
with
1 1
c̄(u; v, w) = ((u · r)v, w) (v, (u · r)w). (14.68)
2 2
We may alternatively use this form in (14.65), where in particular we note
that c(u; v, v) = 0.
138 CHAPTER 14. THE NAVIER-STOKES EQUATIONS

The ✓-method
Semi-discretizaton by the ✓-method takes the form: for each time interval
In = (tn 1 , tn ), with the time step length kn = tn tn 1 , find (Un , Pn ) =
(U (tn ), P (tn )) 2 Vh ⇥ Qh , such that

1
((Un , v) (Un 1 , v))+c(U✓ ; U✓ , v)+⌫a(U✓ , v)+b(v, P✓ ) b(U✓ , q) = (f, v),
kn

for all (v, q) 2 Vh ⇥ Qh , with

U✓ = (1 ✓)Un + ✓Un 1 , P✓ = (1 ✓)Pn + ✓Pn 1 . (14.69)

Here e.g. ✓ = 0 corresponds to the Implicit Euler method, and ✓ = 0.5

corresponds to the Trapezoidal method.

14.5 Stabilized finite element methods

Stabilization techniques
Previously we have found that the inf-sup condition of the Stokes equations
could be circumvented by a pressure stabilization technique to allow for
equal order approximation spaces for the velocity and the pressure. The idea
of stabilization through a small regularizing perturbation is fundamental for
finite element methods, not only in the context of saddle-point problems.
Whereas low Reynolds number flow is dominated by viscosity, for high
Reynolds numbers the dominating phenomenon is nonlinear transport. Stan-
dard Galerkin finite element methods are not optimal for discretization of
transport dominated equations, instead we will use stabilization techniques
to formulate suitable finite element methods.

Linear transport
To understand the basic mechanisms, we first study a linear transport equa-
tion for a scalar quantity u = u(x, t), convected by a divergence-free vector
field = (x, t),

u̇ + ( · r)u ✏ u = f, (x, t) 2 ⌦ ⇥ I, (14.70)

r · = 0, (x, t) 2 ⌦ ⇥ I, (14.71)

with suitable initial and boundary conditions, and ✏ > 0 a small di↵usion
coefficient.
14.5. STABILIZED FINITE ELEMENT METHODS 139

Model problem
We study the simple model problem in one space,
✏u00 + u0 = 0, x 2 (0, 1), (14.72)
u(0) = 1, u(1) = 0, (14.73)
for which we formulate a standard Galerkin finite element method: find
U 2 Vh such that, Z 1 Z 1
0 0
✏u v dx + u0 v dx = 0, (14.74)
0 0
for all test functions v 2 Vh0 , with
Vh = {v 2 H 1 (0, 1) : v(0) = 1, v(1) = 0}, (14.75)
Vh0 = {v 2 H 1 (0, 1) : v(0) = 0, v(1) = 0}. (14.76)
Divide the interval (0, 1) into M uniform subintervals Ii = (xi 1 , xi ) of
length h = xi xi 1 , with nodes {xi }M +1
i=0 and associated piecewise linear
basis functions i = i (x).
Then we can write the finite element approximation as
M
X
U (x) = uj j (x) + u0 0 (x) + uM +1 M +1 (x), (14.77)
j=1

with uj = u(xj ) (since we have a nodal basis), and from the boundary
conditions we have that
M
X
U (x) = uj j (x) + 0 (x). (14.78)
j=1

The discrete system takes the form Ax = b, with A = (aij ), b = (bi ) and
x = (xj ),
Z 1 Z 1
0 0 0
aij = ✏ j (x) i (x) dx + j (x) i (x) dx, (14.79)
0 0
Z 1 Z 1
0 0 0
bi = ✏ 0 (x) i (x) dx + 0 (x) i (x) dx. (14.80)
0 0
(14.81)
Equation i takes the form
X M ✓ ◆ ✓ ◆
✏ 1 2✏ ✏ 1
aij xj = xi 1 + xi + xi+1 + = 0. (14.82)
j=1
h 2 h h 2
140 CHAPTER 14. THE NAVIER-STOKES EQUATIONS

We observe two di↵erent regimes,

✏ 1
>> ) xi 1 + 2xi xi+1 = 0, (14.83)
h 2
✏ 1
<< ) xi 1 + xi+1 = 0, (14.84)
h 2
with a combination of the two when ✏ ⇡ h. In the convection dominated
case, the boundary conditions lead to two cases depending on if M is an odd
or even number; either no solution exists, or the solution oscillates between
0 and 1.
To obtain a finite element approximation that is close to the exact so-
lution in the convection dominated case, we stabilize the method by an
artificial di↵usion ✏ = h/2. We also refer to this as an upwind method, since
the resulting equation takes the form
xi 1 + xi = 0, (14.85)
where information is propagated from the upwind direction.

Streamline di↵usion stabilization

For the Navier-Stokes equations we can use artificial viscosity to stabilize
the finite element method for high Reynolds numbers. But there are also
more accurate, less di↵usive, stabilization methods.
For example, we may use streamline di↵usion stabilization, where we
add artificial viscosity in the streamline direction only, which is enough
to remove the spurious oscillations. Combined with pressure stabilization to
allow for equal order approximation spaces, the method takes the following
form.
For each t > 0, find (U (t), P (t)) 2 Vh ⇥ Qh , such that
(U̇ , v)+c̄(U ; U, v)+a(U, v)+b(v, P ) b(U, q)+s1 (U ; U, v)+s2 (P, q) = (f, v),
for all (v, q) 2 Vh ⇥ Qh , with the stabilization terms
s1 (U ; U, v) = ( 1 (U · r)U, (U · r)v), (14.86)
s2 (P, q) = ( 2 rP, rq), (14.87)
with stabilization parameters 1 ⇠ h/Un 1 and 2 ⇠ h.
By choosing (v, q) = (U, P ), we obtain a stability estimate of the method,
d1 p p p
kU k2 + k ⌫rU k2 + k 1 (U · r)U k2 + k 2 rP k2 = 0, (14.88)
dt 2
where we can observe the regularizing e↵ect of the stabilization terms.
14.6. A POSTERIORI ERROR ESTIMATION 141

Least squares stabilization of the residual

The Galerkin Least Squares (GLS) method is based on a combination of
Galerkin’s method with a least squares minimization of the residual of the
Navier-Stokes equations. GLS is a consistent method in the sense that all
terms in the method are based on the residual of the equations, no artificial
stabilization terms are added.
For each t > 0, find (U (t), P (t)) 2 Vh ⇥ Qh , such that

(U̇ , v)+c̄(U ; U, v)+a(U, v)+b(v, P ) b(U, q)+s1 (U ; U, v)+s2 (U, v) = (f, v),

for all (v, q) 2 Vh ⇥ Qh , with the stabilization terms

s1 (w; U, v) = ( 1 (U̇ + (w · r)U + rP ), v̇ + (w · r)v + rq)

s2 (U, v) = ( 2 r · U, r · v),

with stabilization parameters 1 ⇠ h/Un 1 and 2 ⇠ hUn 1 .

By choosing (v, q) = (U, P ), we obtain a stability estimate of the method,

d1 p p p
kU k2 +k ⌫rU k2 +k 1 (U̇ +(U ·r)U +rP ⌫ U )k2 +k 2 r·U k2 = 0,
dt 2
where we can observe the regularizing e↵ect of the stabilization terms.

14.6 A posteriori error estimation

The model problem
Consider the linear model problem: find u 2 V such that,

a(u, v) = L(v), (14.89)

for all v 2 V . The Galerkin finite element method takes form: find U 2 Vh
such that,
a(U, v) = L(v), (14.90)
for all v 2 Vh , with Vh ⇢ V . Galerkin orthogonality is expressed as

a(e, v) = 0, 8v 2 Vh , (14.91)

for the error e = u U.

142 CHAPTER 14. THE NAVIER-STOKES EQUATIONS

The adjoint equation

We define an output functional of the form
M (u) = (u, ), (14.92)
and to estimate the output error
M (u) M (U ), (14.93)
we introduce the adjoint problem: find ' 2 V , such that
a(v, ') = M (v), 8v 2 V. (14.94)
We have that
M (u) M (U ) = a(u, ') a(U, ') = L(') a(U, ') = r(U, '), (14.95)
with the weak residual
r(U, ') = L(') a(U, '). (14.96)
We can split the integral over the elements K in the mesh T h , so that
X X
M (u) M (U ) = r(U, ') = rK (U, ') = EK , (14.97)
K2T h K2T h

with the error indicator

EK = rK (U, '). (14.98)
To approximate the error indicator we can compute an approximation
to the adjoint problem, ⇡ ', so that
EK ⇡ rK (U, ). (14.99)

The Navier-Stokes equations

The weak form of the Navier-Stokes equations can be formulated in terms
of the weak residual, as the problem to find û = (u, p) 2 V ⇥ Q, such that,
r(û, v̂) = 0, (14.100)
for all v̂ = (v, q) 2 V ⇥ Q. The Galerkin finite element method then takes
the form: find Û = (U, P ) 2 Vh ⇥ Qh , such that,
r(Û , v̂) = 0, (14.101)
for all v̂ = (v, q) 2 Vh ⇥ Qh .
The a posteriori error analysis can then be extended also to the Navier-
Stokes equations, with modifications for the nonlinearity of the equations,
and possible stabilization terms in the finite element method.