Francoise J.-P., Naber G.L., Tsun T.S. (editors) Encyclopedia of Mathematical Physics
Подождите немного. Документ загружается.
adequate hypotheses on time and space scales,
various models are obtained. An important class of
fluid model is one for which the stress tensor is a
linear function of the strain tensor, yielding the so-
called Newtonian fluid model:
s ¼ðp þ ruÞI þ 2 e ½5
Here p is the pressure, I is the identity matrix, and
and are viscosity coefficients. Still assuming
linearity, common models for heat and solute fluxes
consist of assuming
j
T
¼rT; j
¼Dr ½6
where T is the temperature. These are the so-called
Fourier’s law and Fick’s law, respectively.
Having introduced two new quantities, namely
the pressure p and the temperature T, two new
scalar relations are needed to close the system. These
are the state equations. One admissible assumption
consists of setting = (p, T). Another usual addi-
tional hypothesis consists of assuming that the
variations in the internal energy are proportional
to those in the temperature, that is, @e
i
= c
P
@T.
Let us now simplify the above models by
assuming that is constant. Then, mass conserva-
tion implies that the flow is incompressible, that is,
ru = 0. Let us further assume that neither , ,
nor p depend on e
i
. Then, upon abusing the
notation and still denoting by p the ratio p=, the
above set of assumptions yields the so-called
incompressible Navier–Stokes equations:
ru ¼ 0 ½7
@
t
u þ u ru u þrp ¼ f ½8
As a result, the mass and momentum conservation
equations are independent of that of the energy and
those of the solutes:
c
P
ð@
t
T þ u rTÞrðrTÞ¼2 e : e þ q
T
½9
@
t
þ u r
1
rðDrÞ¼
1
q
½10
Another model allowing for a weak dependency of
on the temperature, while still enforcing incom-
pressibility, consists of setting =
0
(1 (T T
0
)).
If buoyancy effects induced by gravity are important,
it is then possible to account for them by setting
f =
0
g(1 (T T
0
)), where g is the gravitational
acceleration, yielding the so-called Boussinesq model.
Variations on these themes are numerous and a
wide range of fluids can be modeled by using
nonlinear constitutive laws and nonlinear state
laws. For the purpose of numerical simulations,
however, it is important to focus on simplified
models.
The Building Blocks
From the above considerations we now extract a
small set of elementary problems which constitute
the building blocks of most numerical methods in
fluid mechanics.
Elliptic Equations
By taking the divergence of the momentum equation
[8] and assuming u to be known and renaming p to
, one obtains the Poisson equation
¼ f ½11
where f is a given source term. This equation plays a
key role in the computation of the pressure when
solving the Navier–Stokes equations; see [54b].
Assuming that adequate boundary conditions are
enforced, this model equation is the prototype for
the class of the so-called elliptic equations. A simple
generalization of the Poisson equation consists of the
advection–diffusion equation
u r rðrÞ¼f ½12
where >0. Admissible boundary conditions are
(@
n
þ r)
j@
= a, r 0, or
j@
= a. This type of
equation is obtained by neglecting the time deriva-
tive in the heat equation [9] or in the solute
conservation equation [10]. Mathematically speak-
ing, [12] is also elliptic since its properties (in
particular, the way the boundary conditions must
be enforced) are controlled by the second-order
derivatives. For the sake of simplicity, assume that
u = 0 in the above equation and that the boundary
condition is
j@
= 0, then it is possible to show that
solves [12] if and only if minimizes the
functional
Jð Þ¼
Z
ðjr j
2
f Þdx
where jjis the Euclidean norm and spans
H ¼ ;
Z
jr j
2
dx < 1;
j@
¼ 0
½13
Writing the first-order optimality condition for this
optimization problem yields
Z
r r ¼
Z
f
for all 2 H. This is the so-called variational
formulation of [12]. When u is not zero, no
variational principle holds but a similar way to
366 Fluid Mechanics: Numerical Methods
reformulate [12] consists of multiplying the equation
by arbitrary functions in H and integrating by parts
the second-order term to give
Z
ðu rÞ þ r r ¼
Z
f ; 8 2 H ½14
This is the so-called weak formulation of [12]. Weak
and variational formulations are the starting point
for finite-element approximations.
Stokes Equations
Another elementary building block is deduced from
[8] by assuming that the time derivative and the
nonlinear term are both small. The corresponding
model is the so-called Stokes equations,
u þrp ¼ f ½15
ru ¼ 0 ½16
Assume for the sake of simplicity that the no-slip
boundary condition is enforced: u
j@
= 0. Introduce
the Lagrangian functional
Lðv; qÞ¼
Z
ðru: rv qrv f vÞdx
Set
X ¼ v;
Z
jrvj
2
dx < 1; v
j@
¼ 0
M ¼ q;
Z
q
2
dx < 1
Then, the pair (u, p) 2 X M solves the Stokes
equations if and only if it is a saddle point of L,thatis,
Lðu; qÞLðu; pÞLðv; pÞ; 8ðv; qÞ2X M ½17
In other words, the pressure p is the Lagrange multi-
plier of the incompressibility constraint ru = 0.
Realizing this fact helps to understand the nature of
the Stokes equations, specially when it comes to
constructing discrete approximations. A variational
formulation of the Stokes equations is obtained by
writing the first-order optimality condition, namely:
Z
ðru : rv prv f vÞdx ¼ 0 8v 2 X
Z
qru dx ¼ 0 8q 2 M
When the nonlinear term is not zero in the
momentum equation, or when this term is linear-
ized, there is no saddle point, but a weak formula-
tion is obtained by multiplying the momentum
equation by arbitrary functions v in X and integrat-
ing by parts the Laplacian, and by multiplying the
mass equation by arbitrary functions q in M:
Z
ððu ruÞv þru:rv prvÞdx ¼
Z
f v ½ 18
Z
qru ¼ 0 ½19
Parabolic Equations
The class of elliptic equations generalizes to that of
the parabolic equations when time is accounted for:
@
t
þ u r rðrÞ¼f ;
jt¼0
¼
0
½20
Fundamentally, this equation has many similarities
with the elliptic equation
þ u r rðrÞ¼f ½21
where >0. In particular, the set of boundary
conditions that are admissible for [20] and [21] are
identical, that is, it is legitimate to enforce (@
n
þ
r)
j@
= a, r 0, or
j@
= a. Moreover, solving [21]
is always a building block of any algorithm solving
[20]. The important fact to remember here is that if
a good approximation technique for solving [21] is
at hand, then extending it to solve [20] is usually
straightforward.
Hyperbolic Equations
When =UL ! 0, where U is the reference velocity
scale and L is the reference length scale, [20]
degenerates into the so-called transport equation
@
t
þ u r ¼ f ½22
This is the prototypical example for the class of
hyperbolic equations. For this equation to be well-
posed, it is necessary to enforce an initial condition
jt = 0
=
0
and an inflow boundary condition, that
is,
j@
= a, where @
= {x 2 @;(u n)(x) < 0} is
the so-called inflow boundary of the domain. To
better understand the nature of this equation,
introduce the characteristic lines X(x, s; t)ofu(x, t)
defined as follows:
d
t
Xðx; s; tÞ¼uðXðx; s; tÞ; t Þ
Xðx; s; sÞ¼x
½23
If u is continuous with respect to t and Lipschitz
with respect to x, this ordinary differential equation
has a unique solution. Furthermore, [22] becomes
d
t
ðXðx; s; tÞ; tÞ½¼f ðXðx; s; tÞ; tÞ½24
Fluid Mechanics: Numerical Methods 367
Then
ðx; tÞ¼
0
ðXðx; t ; 0ÞÞ þ
Z
t
0
f ðXðx; t; Þ;Þd
provided X(x, t ; ) 2 for all 2 [0, t]. This shows
that the concept of characteristic curves is important
to construct an approximation to [22].
Meshes
The starting point of every approximation technique
for solving any of the above model problems consists
of defining a mesh of on which the approximate
solution is defined. To avoid having to account for
curved boundaries, let us assume that the domain is
a two-dimensional polygon (resp. three-dimensional
polyhedron). A mesh of ,sayT
h
, is a partition of
into small cells, hereafter assumed to be simple
convex polygons in two dimensions (resp. polyhe-
drons in three dimensions), say triangles or quad-
rangles (resp. tetrahedrons or cuboids). Moreover,
this partition is usually assumed to be such that if
two different cells have a nonempty intersection, then
the intersection is a vertex, or an entire edge, or an
entire face. The left panel of Figure 1 shows a mesh
satisfying the above requirement. The mesh in the
right panel is not admissible.
Finite Elements: Interpolation
The finite-element method is foremost an interpola-
tion technique. The goal of this section is to
illustrate this idea by giving examples.
Let T
h
= {K
m
}
1mN
el
be a mesh composed of N
el
simplices, that is, triangles in two dimensions or
tetrahedrons in three dimensions. Consider the
following vector spaces of functions:
V
h
¼fv
h
2C
0
ð
Þ; v
hjK
m
2 P
k
; 1 m N
el
g½25
where P
k
denotes the space of polynomials of global
degree at most k. V
h
is called a finite-element
approximation space. We now construct a basis for V
h
.
Given a simplex K
m
in R
d
, let v
n
be a vertex of
K
m
, let F
n
be the face of K
m
opposite to v
n
, and
define n
n
to be the outward normal to F
n
,1 n
d þ 1. Define the barycentric coordinates
n
ðxÞ¼1
ðx v
n
Þn
n
ðv
l
v
n
Þn
n
; 1 n d þ 1 ½26
where v
l
is an arbitrary vertex in F
n
(the definition
of
n
is clearly independent of v
l
provided v
l
belongs
to F
n
). The barycentric coordinate
n
is an affine
function; it is equal to 1 at v
n
and vanishes on F
n
; its
level sets are hyperplanes parallel to F
n
. The
barycenter of K
m
has barycentric coordinates
1
d þ 1
; ...;
1
d þ 1
The barycentric coordinates satisfy the following
properties: for all x 2 K
m
,0
n
(x) 1, and for all
x 2 R
d
,
X
dþ1
n¼1
n
ðxÞ¼1 and
X
dþ1
n¼1
n
ðxÞðx v
n
Þ¼0
Consider the set of nodes {a
n,m
}
1nn
sh
of K
m
with
barycentric coordinates
i
0
k
; ...;
i
d
k
; 0 i
0
; ...; i
d
k; i
0
þþi
d
¼ k
These points are called the Lagrange nodes of K
m
.It
is clear that there are n
sh
= (1=2)(k þ 1)(k þ 2)
of these points in two dimensions and n
sh
= (1=6)
(k þ 1)(k þ 2)(k þ 3) in three dimensions. It is
remarkable that n
sh
= dim P
k
.
Let {b
1
, ..., b
N
} =
S
K
m
2T
h
{a
1,m
, ..., a
n
sh
,m
}bethe
set of all the Lagrange nodes in the mesh. For K
m
2T
h
and n 2 {1, ..., n
sh
}, let j(n, m) 2 {1, ..., N}bethe
integer such that a
n, m
= b
j(n, m)
; j(n, m) is the global
index of the Lagrange node a
n, m
. Let {’
1
, ..., ’
N
}be
the set of functions in V
h
defined by ’
i
(b
j
) =
ij
,thenit
can be shown that
f’
1
; ...;’
N
g is a basis for V
h
½27
The functions ’
i
are called global shape functions.
An important property of global shape functions is
that their supports are small sets of cells. More
precisely, let i 2 {1, ..., N} and let V
i
= {m; 9n;
i = j(n, m)} be the set of cell indices to which the
node b
i
belongs, then the support of ’
i
is
S
m2V
i
K
m
.
For k = 1, it is clear that ’
i jK
m
=
n
for all m 2V
i
and all n such that i = j(n, m), and ’
i jK
m
= 0
otherwise. The graph of such a shape function in
two dimensions is shown in the left panel of Figure 2.
For k = 2, enumerate from 1 to d þ 1 the vertices of
K
m
, and enumerate from d þ 2ton
sh
the Lagrange
nodes located at the midedges. For a midedge node
of index d þ 2 n n
sh
, let b(n), e(n) 2 {1, ...,
d þ 1} be the two indices of the two Lagrange
Figure 1 Admissible (left) and nonadmissible (right) meshes.
368 Fluid Mechanics: Numerical Methods
nodes at the extremities of the edge in question. Then,
the restriction to K
m
of a P
2
shape function ’
i
is
’
ijK
m
¼
n
ð2
n
1Þ; if 1 n d þ 1
4
bðnÞ
eðnÞ
; if d þ 2 n n
sh
½28
Figure 2 shows the graph of two P
2
shape functions
in two dimensions.
Once the space V
h
is introduced, it is natural to
define the interpolation operator
h
: C
0
ð
Þ3v 7!
X
N
i¼1
vðb
i
Þ’
i
2 V
h
½29
This operator is such that for all continuous
functions v, the restriction of
h
(v) to each mesh
cell is a polynomial in P
k
and
h
(v) takes the same
values as v at the Lagrange nodes. Moreover, setting
h = max
K
m
2T
h
diam(K
m
), and defining
krk
L
p
¼
Z
jrj
p
dx
1=p
for 1 p < 1
the following approximation holds:
kv
h
ðvÞk
L
p
þ hkrðv
h
ðvÞÞk
L
p
ch
kþ1
kvk
C
kþ1
ð
Þ
½30
where c is a constant that depends on the quality of
the mesh. More precisely, for K
m
2T
h
, let
K
m
be
the diameter of the largest ball that can be inscribed
into K
m
and let h
K
m
be the diameter of K
m
. Then, c
depends on = max
K
m
2T
h
h
K
m
=
K
m
. Hence, for the
mesh to have good interpolation properties, it is
recommended that the cells be not too flat. Families
of meshes for which is bounded uniformly with
respect to h as h ! 0 are said to be shape-regular
families.
The above example of finite-element approxima-
tion space generalizes easily to meshes composed of
quadrangles or cuboids. In this case, the shape
functions are piecewise polynomials of partial
degree at most k. These spaces are usually referred
to as Q
k
approximation spaces.
Finite Elements: Approximation
We show in this section how finite-element approx-
imation spaces can be used to approximate some
modelproblemsexhibitedinthesection‘‘Building
blocks.’’
Advection–Diffusion
Consider the model problem [21] supplemented
with the boundary condition (@
n
þ r)
j@
= g.
Assume >0, þ (1=2)ru 0, and r 0. Define
að; Þ¼
Z
ð þ rÞ þ r r ðÞdx
þ
Z
@
r ds
Then, the weak formulation of [21] is: seek 2 H
(H defined in [13]) such that for all 2 H
að; Þ¼
Z
f dx þ
Z
@
g ds ½31
Using the approximation space V
h
defined in [25]
together with the basis defined in [27], we seek an
approximate solution to the above problem in the
form
h
=
P
N
i = 1
U
i
’
i
2 V
h
. Then, a simple way of
approximating [31] consists of seeking U =
(U
1
, ..., U
N
)
T
2 R
N
such that for all 1 i N
að
h
;’
i
Þ¼
Z
f ’
i
dx þ
Z
@
g’
i
ds ½32
This problem finally amounts to solving the follow-
ing linear system:
AU ¼ F ½33
Figure 2 Two-dimensional Lagrange shape functions: piecewise P
1
(left) and piecewise P
2
(center and right).
Fluid Mechanics: Numerical Methods 369
where A
ij
= a(’
j
, ’
i
) and
F
i
¼
Z
f ’
i
dx þ
Z
@
g’
i
ds
The above approximation technique is usually
referred to as the Galerkin method. The following
error estimate can be proved:
k
h
k
L
p
þ hkrð
h
Þk
L
p
ch
kþ1
kk
C
kþ1
ð
Þ
½34
where, in addition to depending on the shape
regularity of the mesh, the constant c also depends
on , , and .
Stokes Equations
The line of thought developed above can be used to
approximate the Navier–Stokes problem [15]–[16].
Let us assume that the nonlinear term u ru is
linearized in the form v ru, where v is known. Let
T
h
be a mesh of , and assume that finite-element
approximation spaces have been constructed to
approximate the velocity and the pressure, say X
h
and M
h
. Assume for the sake of simplicity that X
h
X and M
h
M. Assume that bases for X
h
and M
h
are at hand, say {’
1
, ..., ’
N
u
} and {
1
, ...,
N
p
},
respectively. Set
aðu; jÞ¼
Z
ððv ru Þj þ ru : rjdx
and
bðv; Þ¼
Z
rvdx
Then, we seek an approximate velocity u
h
=
P
N
u
i = 1
U
i
j
i
and an approximate pressure p
h
=
P
N
p
k = 1
P
k
k
such that for all i 2 {1, ..., N
u
} and all
k 2 {1, ..., N
p
} the following holds:
aðu
h
; j
i
Þþbðj
i
; p
h
Þ¼
Z
f j
i
dx ½35
bðu
h
;
k
Þ¼0 ½36
Define the matrix A2R
N
u
, N
u
such that
A
ij
= a(j
j
, j
i
). Define the matrix B2R
N
p
, N
u
such
that B
ki
= b(j
i
,
k
). Then, the above problem can be
recast into the following partitioned linear system:
AB
T
B 0
U
P
¼
F
0
½37
where the vector F 2 R
N
u
is such that F
i
=
R
f j
i
.
An important aspect of the above approximation
technique is that, for the linear system to be
invertible, the matrix B
T
must have full row rank
(i.e., B has full column rank). This amounts to
9
h
> 0; inf
q
h
2M
h
sup
v
h
2X
h
R
q
h
rv
h
dx
kv
h
k
X
kq
h
k
M
h
½38
where
kv
h
k
2
X
¼
Z
jrv
h
j
2
dx; kq
h
k
2
M
¼
Z
q
2
h
dx
This nontrivial condition is called the Ladyzˇenskaja–
Babusˇka–Brezzi condition (LBB) in the literature.
For instance, if P
1
finite elements are used to approx-
imate both the velocity and the pressure, the above
condition does not hold, since there are nonzero
pressure fields q
h
in M
h
such that
R
q
h
rv
h
dx = 0
for all v
h
in X
h
. Such fields are called spurious
pressure modes. An example is shown in Figure 3.
The spurious function alternatively takes the values
1, 0, and þ1 at the vertices of the mesh so that its
mean value on each cell is zero.
Couples of finite-element spaces satisfying the
LBB condition are numerous. For instance, assuming
k 2, using P
k
finite elements to approximate the
velocity and P
k1
finite elements to approximate the
pressure is acceptable. Likewise, using Q
k
elements
for the velocity and Q
k1
elements for the pressure
on meshes composed of quadrangles or cuboids is
admissible.
Approximation techniques for which the pressure
and the velocity degrees of freedom are not
associated with the same nodes are usually called
staggered approximations. Staggering pressure and
velocity unknowns is common in solution methods
for the incompressible Stokes and Navier–Stokes
equations;seealsothesubsection‘‘Stokes
equations.’’
Finite Volumes: Principles
The finite-volume method is an approximation
technique whose primary goal is to approximate
conservation equations, whether time dependent or
+1 +1
+1 +1
+1
+1
+1 +1
+1
+1
–1
–1
–1
–1
–1
–1
–1
–1
–100
00
0
0
0
0
00–1+1 +1
0
0
–1
–1
Figure 3 The P
1
=P
1
finite element: the mesh (left); one
pressure spurious mode (right).
370 Fluid Mechanics: Numerical Methods
not. Given a mesh, say T
h
= {K
m
}
1mN
el
,
and a
conservation equation
@
t
þrFð; r; x; tÞ¼f ½39
( = 0 if the problem is time independent and = 1
otherwise), the main idea underlying every finite-
volume method is to represent the approximate
solution by its mean values over the mesh cells
(
K
1
, ...,
K
N
el
)
T
2 R
N
el
and to test the conservation
equation by the characteristic functions of the mesh
cells {1
K
1
, ...,1
K
N
el
}. For each cell K
m
2T
h
,denoteby
n
K
m
the outward unit normal vector and denote by F
m
the set of the faces of K
m
. The finite-volume approx-
imation to [39] consists of seeking (
K
1
, ...,
K
N
el
)
T
2
R
N
el
such that the function
h
=
P
N
el
m = 1
K
m
1
K
1
m
satisfies the following: for all 1 m N
el
jK
m
jd
t
K
m
ðtÞþ
X
2F
m
F
m;
h
ð
h
;r
h
h
;tÞ¼
Z
K
f dx ½40
where
jK
m
j¼
Z
K
dx
r
h
h
is an approximation of r, and F
m,
h
is an
approximation of
Z
Fð; r; x; t Þn
K
m
ds
The precise definition of the so-called approximate
flux F
m,
h
depends on the nature of the problem
(e.g., elliptic, parabolic, hyperbolic, saddle point)
and the desired accuracy. In general, the approx-
imate fluxes are required to satisfy the following
two important properties:
1. Conservativity: for K
m
, K
l
2T
h
such that
= K
m
\ K
l
, F
m,
h
= F
l,
h
.
2. Consistency: let be the solution to [39], and set
h
¼
1
K
1
jK
1
j
Z
K
1
dx þþ
1
K
N
el
jK
N
el
j
Z
K
N
el
dx
then
F
m;
h
ð
h
;r
h
h
;tÞ!
Z
Fð ;r ;x;tÞnds as h !0
The quantity
F
m;
h
ð
h
; r
h
h
; tÞ
Z
Fð ; r ; x; tÞn ds
is called the consistency error.
Note that [40] is a system of ordinary differential
equations. This system is usually discretized in time
by using standard time-marching techniques such as
explicit Euler, Runge–Kutta, etc.
The discretization technique described above is
sometimes referred to as cell-centered finite-volume
method. Another method, called vertex-centered
finite volume method, consists of using the char-
acteristic functions associated with the vertices of
the mesh instead of those associated with the cells.
Finite Volumes: Examples
In this section we illustrate the ideas introduced
above. Three examples are developed: the Poisson
equation, the transport equation, and the Stokes
equations.
Poisson Problem
Consider the Poisson equation [11] equipped with
the boundary condition @
n
j@
= a. To avoid techni-
cal details, assume that =[0, 1]
d
. Let K
h
be a mesh
of composed of rectangles (or cuboids in three
dimensions).
The flux function is F(, r, x) = r; hence,
F
m,
h
must be a consistent conservative approxima-
tion of
R
n
K
m
r ds.Let be an interior face of
the mesh and let K
m
, K
l
be the two cells such that
= K
m
\ K
l
. Let x
K
m
, x
K
l
be the barycenters of K
m
and K
l
, respectively. Then, an admissible formula
for the approximate flux is
F
m;
h
¼
jj
jx
K
m
x
K
l
j
ð
K
l
K
m
Þ½41
where jj=
R
ds. The consistency error is O(h)in
general, and is O(h
2
) if the mesh is composed of
identical cuboids. The conservativity is evident. If
is part of @, an admissible formula for the
approximate flux is F
m,
h
=
R
a ds. Then, upon
defining F
i
K
m
= F
K
m
n@ and F
@
K
m
= F
K
m
\ @, the
finite-volume approximation of the Poisson problem
is: seek
h
2 R
N
el
such that for all 1 m N
el
X
2F
i
K
m
F
m;
h
¼
Z
K
m
f dx þ
X
2F
@
K
m
Z
a ds ½42
Transport Equation
Consider the transport equation
@
t
þrðuÞ¼f ½43
jt¼0
¼
0
;
j@
¼ a ½44
where u(x, t) is a given field in C
1
(
[0, T]). Let T
h
be a mesh of . For the sake of simplicity, let us use
the explicit Euler time-stepping to approximate [40].
Fluid Mechanics: Numerical Methods 371
Let N be positive integer, set t = T=N, set t
n
= nt
for 0 n N, and partition [0, T] as follows:
½0; T¼
[
N1
n¼0
½t
n
; t
nþ1
Denote by
n
h
2 R
N
el
the finite-volume approxima-
tion of
h
(t
n
). Then, [40] is approximated as
follows:
jK
m
j
t
ð
nþ1
K
m
n
K
m
Þþ
X
2F
m
F
m;
h
ð
h
; r
h
h
; t
n
Þ
¼
Z
K
f ðx; t
n
Þdx ½45
where
0
K
m
=
R
K
m
0
dx. The approximate flux F
m,
h
must be a consistent conservative approximation of
R
(u n
K
m
) ds. Let be a face of the mesh and let
K
m
, K
l
be the two cells such that = K
m
\ K
l
(note
that if is on @, belongs to one cell only and we
set K
m
= K
l
). If is on @
, set
F
m;
h
¼
Z
ðu n
K
m
Þa ds ½46
If is not on @
, set u
n
m,
=
R
(u n
K
m
)ds and
define
F
m;
h
¼
n
K
m
u
n
m;
if u
n
m;
0
n
K
l
u
n
m;
if u
n
m;
< 0
(
½47
The above choice for the approximate flux is usually
called the upwind flux. It is consistent with the analysis
that has been done for [22], that is, information flows
along the characteristic lines of the field u;see[24].In
other words, the updating of
nþ1
k
m
must be done by
using the approximate values
n
h
coming from the cells
that are upstream the flow field.
An important feature of the above approximation
technique is that it is L
1
-stable, in the sense that
max
0nN;1mN
el
n
K
m
cðu
0
; f Þ
if the two mesh parameters t and h satisfy
the so-called Courant–Friedrichs–Levy (CFL)
condition kuk
L
1
t=h c(), where c() is a con-
stant that depends on the mesh regularity parameter
= max
K
m
2T
h
h
K
m
=
K
m
. In one dimension, c() = 1.
Stokes Equations
To finish this short review of finite-volume methods,
we turn our attention to the Stokes problem (15)–(16)
equipped with the homogeneous Dirichlet boundary
condition u
j@
= 0.
Let T
h
be a mesh of composed of triangles (or
tetrahedrons). All the angles in the triangulation are
assumed to be acute so that, for all K 2T
h
, the
intersection of the orthogonal bisectors of the sides
of K, say x
K
,isinK. We propose a finite-volume
approximation for the velocity and a finite-element
approximation for the pressure. Let {e
1
, ..., e
d
}bea
Cartesian basis for R
d
. Set 1
k
K
m
= 1
K
m
e
k
for all 1
m N
el
and 1 k d ; then define
X
h
¼ span 1
1
K
1
; ...; 1
d
K
1
; ...; 1
1
K
N
el
; ...; 1
d
K
N
el
no
Let {b
1
, ..., b
N
v
} be the vertices of the mesh, and let
{’
1
, ..., ’
N
v
} be the associated piecewise linear
global shape functions. Then, set (see the section
‘‘Finiteelements:interpolation’’)
N
h
¼ spanf’
1
; ...;’
N
v
g
M
h
¼fq 2 N
h
;
Z
q dx ¼ 0g
The approximate problem consists of seeking
(u
K
1
, ..., u
K
N
el
) 2 R
dN
el
and p
h
2 M
h
such that for
all 1 m N
el
,1 k d, and all 1 i N
v
,
X
2F
m
1
k
K
m
F
m;
h
þ c 1
k
K
m
; p
h
¼
Z
K
m
1
k
K
m
f dx ½48
cðu
K
m
;’
i
Þ¼0 ½49
where
cðv
K
m
; p
h
Þ¼
Z
K
m
v
K
m
rp
h
dx
Moreover,
F
m;
h
¼
jj
jx
m
x
l
j
ðu
K
m
u
K
l
Þ if ¼ K
m
\ K
l
jj
dðx
m
;Þ
u
K
m
if ¼ K
m
\ @
8
>
>
>
<
>
>
>
:
where d(x
K
m
, ) is the Euclidean distance between
x
K
m
and . This formulation yields a linear system
with the same structure as in [37]. Note in particular
that
sup
v
h
2X
h
cðv
h
; p
h
Þ
kv
h
k
L
1
¼krp
h
k
L
1
½50
Since the mean value of p
h
is zero, krp
h
k
L
1
is a norm
on M
h
. As a result, an inequality similar to [38] holds.
This inequality is a key step to proving that the linear
system is wellposed and the approximate solution
converges to the exact solution of (15)–(16).
372 Fluid Mechanics: Numerical Methods
Projection Methods for Navier–Stokes
In this section we focus on the time approximation
of the Navier–Stokes problem:
@
t
u u þ u ru þrp ¼ f ½51a
ru ¼ 0 ½51b
u
j@
¼ 0 ½51c
u
jt¼0
¼ u
0
½51d
where f is a body force and u
0
is a solenoidal
velocity field. There are numerous ways to discretize
this problem in time, but, undoubtedly, one of the
most popular strategies is to use projection methods,
sometimes also referred to as Chorin–Temam
methods.
A projection method is a fractional-step time-
marching technique. It is a predictor–corrector
strategy aiming at uncoupling viscous diffusion and
incompressibility effects. One time step is composed
of three substeps: in the first substep, the pressure is
made explicit and a provisional velocity field is
computed using the momentum equation; in the
second substep, the provisional velocity field is
projected onto the space of incompressible (solenoi-
dal) vector fields; in the third substep, the pressure is
updated.
Let q > 0 be an integer and approximate the time
derivative of u using a backward difference formula of
order q.Tothisend,introduceapositiveintegerN,set
t = T=N,sett
n
= nt for 0 n N, and consider a
partitioning of the time interval in the form
½0; T¼
[
N1
n¼0
½t
n
; t
nþ1
For all sequences v
t
= (v
0
, v
1
, ..., v
N
), set
D
ðqÞ
v
nþ1
¼
q
v
nþ1
X
q1
j¼0
j
v
nj
½52
where q 1 n N 1. The coefficients
j
are
such that
1
t
ð
q
uðt
nþ1
Þ
X
q1
j¼0
j
uðt
nj
ÞÞ
is a qth-order backward difference formula approx-
imating @
t
u(t
nþ1
). For instance,
D
ð1Þ
v
nþ1
¼ v
nþ1
v
n
D
ð2Þ
v
nþ1
¼
3
2
v
nþ1
2v
n
þ
1
2
v
n1
Furthermore, for all sequences
t
= (
0
,
1
, ...,
N
),
define
?;nþ1
¼
X
q1
j¼0
j
nj
½53
so that
P
q1
j = 0
j
p(t
nj
)isa(q 1)th-order extrapola-
tion of p(t
nþ1
). For instance, p
?, nþ1
= 0 for
q = 1, p
?, nþ1
= p
n
for q = 2, and p
?, nþ1
= 2p
n
p
n1
for q = 3. Finally, denote by (u ru )
?, nþ1
a qth-
order extrapolation of (u ru)(t
nþ1
). For instance,
ðu ruÞ
?;nþ1
¼
u
n
ru
n
for q ¼1
2u
n
ru
n
u
n1
ru
n1
if q ¼2
A general projection algorithm is as follows. Set
~
u
0
= u
0
and
l
= 0 for 0 l q 1. If q > 1;
assume that
~
u
1
, ...,
~
u
q1
, p
?, q
and (u ru)
?, q
have
been initialized properly. For n q 1, seek
~
u
nþ1
such that
~
u
nþ1
j@
= 0and
D
ðqÞ
t
~
u
nþ1
~
u
nþ1
þr p
?;nþ1
þ
X
q1
j¼0
j
t
nj
!
¼ S
nþ1
½54a
where S
nþ1
= f (t
nþ1
) (u ru)
?, nþ1
. Then solve
nþ1
¼r
~
u
nþ1
;@
n
nþ1
j@
¼ 0 ½54b
Finally, update the pressure as follows:
p
nþ1
¼
q
t
nþ1
þ p
?;nþ1
r
~
u
nþ1
½54c
The algorithm [54a–c] is known in the literature as
the rotational form of the pressure-correction
method. Upon denoting u
t
= (u(t
0
), ..., u(t
N
))
and p
t
= (p(t
0
), ..., p(t
N
)), the above algorithm
has been proved to yield the following error
estimates:
ku
t
~
u
t
k
‘
2
ðL
2
Þ
ct
2
krðu
t
~
u
t
Þk
‘
2
ðL
2
Þ
þkp
t
p
t
k
‘
2
ðL
2
Þ
ct
3=2
where k
t
k
2
‘
2
(L
2
)
=t
P
N
n = 0
R
j
n
j
2
dx.
A simple strategy to initialize the algorithm
consists of using D
(1)
u
1
at the first step in [54a];
then using D
(2)
u
2
at the second step, and proceed-
ing likewise until
~
u
1
, ...,
~
u
q1
have all been
computed.
At the present time, projection methods count among
the few methods that are capable of solving the time-
dependent incompressible Navier–Stokes equations in
three dimensions on fine meshes within reasonable
Fluid Mechanics: Numerical Methods 373
computation times. The reason for this success is that
the unsplit strategy, which consists of solving
D
ðqÞ
t
u
nþ1
u
nþ1
þrp
nþ1
¼ S
nþ1
½55a
ru
nþ1
¼ 0; u
nþ1
j@
¼ 0 ½ 55b
yields a linear system similar to [37], which usually takes
far more time to solve than sequentially solving [54a]
and [54b]. It is commonly reported in the literature that
the ratio of the CPU time for solving [55a]–[55b] to that
for solving [54a–c] rangesbetween10to30.
See also: Compressible Flows: Mathematical Theory;
Computational Methods in General Relativity: The Theory;
Geophysical Dynamics; Image Processing: Mathematics;
Incompressible Euler Equations: Mathematical Theory;
Interfaces and Multicomponent Fluids;
Magnetohydrodynamics; Newtonian Fluids and
Thermohydraulics; Non-Newtonian Fluids; Partial
Differential Equations: Some Examples; Variational
Methods in Turbulence.
Further Reading
Doering CR and Gibbon JD (1995) Applied Analysis of the
Navier–Stokes Equations, Cambridge Texts in Applied Mathe-
matics. Cambridge: Cambridge University Press.
Ern A and Guermond J-L (2004) Theory and Practice of Finite
Elements. Springer Series in Applied Mathematical Sciences,
vol. 159. New York: Springer-Verlag.
Eymard R, Galloue¨t T, and Herbin R (2000) Finite volume methods.
In: Ciarlet PG and Lions JL (eds.) Handbook of Numerical
Analysis, vol. VII , pp. 713–1020. Amsterdam: North-Holland.
Karniadakis GE and Sherwin SJ (1999) Spectral/hp Element
Methods for CFD, Numerical Mathematics and Scientific
Computation. New York: Oxford University Press.
Rappaz M, Bellet M, and Deville M (2003) Numerical Modeling in
Material Science and Engineering. Springer Series in Computa-
tional Mathematics, vol. 32. Berlin: Springer Verlag.
Temam R (1984) Navier–Stokes Equations. Theory and Numer-
ical Analysis. Studies in Mathematics and its Applications,
vol. 2. Amsterdam: North-Holland.
Toro EF (1997) Riemann Solvers and Numerical Methods for
Fluid Dynamics. A Practical Introduction. Berlin: Springer.
Wesseling P (2001) Principles of Computational Fluid Dynamics.
Springer Series in Computational Mathematics, vol. 29. Berlin:
Springer.
Fourier Law
F Bonetto, Georgia Institute of Technology, Atlanta,
GA, USA
L Rey-Bellet, University of Massachusetts, Amherst,
MA, USA
ª 2006 Elsevier Ltd. All rights reserved.
Introduction
In the famous 1822 treatise by Jean Baptiste Joseph
Fourier, The´orie analytique de la chaleur,theDiscours
pre´liminaire opens with: ‘‘Primary causes are
unknown to us; but are subject to simple and constant
laws, which may be discovered by observation, the
study of them being the subject of natural philosophy.
Heat, like gravity, penetrates every substance of the
universe, its rays occupy all parts of space. The object
of our work is to set forth the mathematical laws
which this element obeys. The theory of heat will
hereafter form one of the most important branches of
general physics.’’ After a brief discussion of rational
mechanics, he continues with the sentence: ‘‘But
whatever may be the range of mechanical theories,
they do not apply to the effects of heat. These make up
a special order of phenomena, which cannot be
explained by the principles of motion and equilibria.’’
Fourier goes on with a thorough description of the
phenomenology of heat transport and the derivation of
the partial differential equation describing heat trans-
port: the heat equation. A large part of the treatise is
then devoted to solving the heat equation for various
geometries and boundary conditions. Fourier’s treatise
marks the birth of Fourier analysis. After Boltzmann,
Gibbs, and Maxwell and the invention of statistical
mechanics in the decades after Fourier’s work, we
believe that Fourier was wrong and that, in principle,
heat transport can and should be explained ‘‘by the
principles of motion and equilibria,’’ that is, within the
formalism of statistical mechanics. But well over a
century after the foundations of statistical mechanics
were laid down, we still lack a mathematically
reasonable derivation of Fourier’s law from first
principles. Fourier’s law describes the macroscopic
transport properties of heat, that is, energy, in none-
quilibrium systems. Similar laws are valid for the
transport of other locally conserved quantities, for
example, charge, particle density, momentum, etc. We
will not discuss these laws here, except to point out
that in none of these cases macroscopic transport laws
have been derived from microscopic dynamics. As
Peierls once put it: ‘‘It seems there is no problem in
modern physics for which there are on record as many
false starts, and as many theories which overlook some
essential feature, as in the problem of the thermal
conductivity of [electrically] non-conducting crystals.’’
Macroscopic Law
Consider a macroscopic system characterized at
some initial time, say t = 0, by a nonuniform
374 Fourier Law
temperature profi le T
0
(r). This temperature profile
will generate a heat, that is, energy current J(r).
Due to energy conservation and basic
thermodynamics:
c
v
ðTÞ
@
@t
Tðr; tÞ ¼ r J ½1
where c
v
(T) is the specific heat per unit volume. On the
other hand, we know that if the temperature profile is
uniform, that is, if T
0
(r) T
0
, there is no current in
the system. It is then natural to assume that, for small
temperature gradients, the current is given by
JðrÞ¼ðTðrÞÞrTðrÞ½2
where (T) is the conductivity. Here we have
assumed that t here is no mass flow or other mode
of energy transport besides heat conduction (we
also ignore, for simplicity, any variations in density
or pressure). Equation [2] is normally called as
Fourier’s law. Putting together e qns [1] and [2],we
get the heat equation:
c
v
ðTÞ
@
@t
Tðr; tÞ¼rðTÞrT½½3
This equation must be completed with suitable
boundary conditions. Let us consider two distinct
situations in which the heat equation is observed to
hold experimentally with high precision:
1. An isolated macroscopic system, for example,
a fluid or solid in a domain surrounded
by effectively adiabatic walls. In this case,
eqn [3] is to be solved subject to the initial
condition T(r,0)= T
0
(r) and no heat flux
across the boundary of (denoted by @), that
is, n(r) rT(r) = 0ifr 2 @ with n the normal
vector to @ at r.Ast !1, the system reaches a
stationary state characterized by a uniform
temperature
T determined by the constancy of
the total energy.
2. A system in contact with heat reservoirs. Each
reservoir fixes the temperature of some portion
(@)
of the boundary @. The rest of the
boundary is insulated. When the system reaches
a stationary state (again assuming no matter
flow), its temperature will be given by the
solution of eqn [3] with the left-hand side set
equal to zero,
r
~
JðrÞ¼rðr
~
TðrÞÞ ¼ 0 ½4
subject to the boundary condition
~
T(r) = T
for
r 2 (@)
and no flux across the rest of the
boundary.
The simplest geometry for a conducting system is
that of a cylindrical slab of height h and cross-
sectional area A. It can be either a cylindrical
container filled with a fluid or a piece of crystalline
solid. In both cases, one keeps the lateral surface of
the cylinder insulated. If the top and the bottom of
the cylinder are also insulated we are in case (1). If
one keeps the top and the bottom in contact with
thermostats at temperatures T
h
and T
b
, respectively,
this is (for a fluid) the usual setup for a Benard
experiment. To avoid convection, one has to make
T
h
> T
b
or keep jT
h
T
b
j small. Assuming unifor-
mity in the direction perpendicular to the vertical
x-axis one has, in the stationary state, a tempera-
ture profile
~
T(x) with
~
T(0) = T
b
,
~
T(h) = T
h
and
(
~
T)d
~
T=dx = const. for x 2 (0, h).
In deriving the heat equation, we have implicitly
assumed that the system is described fully by specifying
its temperature T(r, t)everywherein.Whatthis
means on the microscopic level is that we imagine the
system to be in local thermal equilibrium (LTE).
Heuristically, we might think of the system as being
divided up (mentally) into many little cubes, each large
enough to contain very many atoms yet small enough
on the macroscopic scale to be accurately described, at
a specified time t, as a system in equilibrium at
temperature T(r
i
, t), where r
i
is the center of the ith
cube. For slow variation in space and time, we can
then use a continuous description T(r, t). The theory
of the heat equation is very developed and, together
with its generalizations, plays a central role in modern
analysis. In particular, one can consider more general
boundary conditions. Here we are interested in the
derivation of eqn [2] from first principles. This clearly
presupposes, as a first fundamental step, a precise
definition of the concept of LTE and its justification
within the law of mechanics.
Empirical Argument
A theory of heat conduction has as a goal the
computation of the conductivity (T) for realistic
models, or, at the very least, the derivation of
behavior of (T) as a function of T. The early
analysis was based on ‘‘kinetic theory.’’ Its applica-
tion to heat conduction goes back to the works of
Clausius, Maxwell, and Boltzmann, who obtained a
theoretical expression for the heat conductivity of
gases,
ffiffiffiffi
T
p
, independent of the gas density. This
agrees with experiment (when the density is not too
high) and was a major early achievement of the
atomic theory of matter.
Heat Conduction in Gases
Clausius and Maxwell used the concept of a ‘‘mean
free path’’ : the average distance a particle (atom or
Fourier Law 375