Finite volume method

Integral conservation law

Z Z Z
d
u dx + f · nds = q dx
dt Vi Si Vi

Application to control volumes

N
[
V̄ = V̄i , Vi ∩ Vj = ∅, ∀i 6= j
i=1

Discrete solution values

Z
1
ui (t) = u(x, t) dx ≈ u(xi , t)
|Vi | Vi
Control volumes

Vertex-centered FVM Cell-centered FVM

1D ui 1D ui
ui+1 ui+1
ui 1 ui 1

Vi Vi

xi 1 xi xi+1 xi 1 xi 1=2 xi
xi 1=2 xi+1=2

2D 2D

Vi

Vi Vi
Vi
Discretization of local problems

Integral conservation law for a single control volume

Z
∂ui 1 X
+ f · nk dS = qi
∂t |Vi | Sk
k
Z Z
1 1
ui = u dV, qi = q dV
|Vi | Vi |Vi | Vi

Approximation of integrals in terms of the mean values ui

One algebraic equation or ODE for each control volume

Numerical integration

Approximation of integrals using numerical quadrature

Z p
X
Φ(x)dx ≈ |K| ωr Φ(xr )
K r=1

where xr ∈ K̄ are nodes and ωr are weights satisfying

p
X
ωr = 1
r=1

Derivation via exact integration of interpolation polynomials

p Z
X 1
Φ(x) ≈ Φ(xr )lr (x) ⇒ ωr = lr (x)dx
r=1
|K| K
Newton-Cotes quadrature in 1D

Midpoint rule
R
K
Φ(x) dx ≈ |K|Φ12 exact for Φ ∈ P1 (K)

Trapezoidal rule
Φ(x) dx ≈ |K| Φ1 +Φ
R
K 2
2
exact for Φ ∈ P1 (K)

Simpson’s rule
Φ(x) dx ≈ |K| Φ1 +4Φ612 +Φ2
R
K
exact for Φ ∈ P3 (K)

x1 + x2
K = (x1 , x2 ), |K| = x2 − x1 , x12 =
2
Gaussian quadrature in 1D

Optimal distribution of nodes: p-point Gaussian quadrature is exact

for polynomials of degree 2p − 1

Higher accuracy than Newton-Cotes approximations using the same

number of (equidistant) nodes

Simpson’s rule (3 points) is exact for Φ ∈ P3 (K)

2-point Gaussian formula is exact for Φ ∈ P3 (K)
3-point Gaussian formula is exact for Φ ∈ P5 (K)

Gaussian quadrature points are roots of Legendre polynomials

Gaussian quadrature in 1D

Nodes and weights for K = (0, 1)

p i xi ωi exact for
1 1 1/2 1 P1 (K)
√
1 (3 − 3)/6 1/2
√
2 2 (3 + 3)/6 1/2 P3 (K)
√
1 (5 − 15)/6 5/18
1
3 2 2 4/9 P5 (K)
√
3 (5 + 15)/6 5/18
Interpolation techniques

Volume integrals
Z
1
ui = u dV ≈ u(x̄i ) midpoint rule
|Vi | Vi

Cell averages represent second-order approximations to midpoint values

Surface integrals
Z
1 X
f = vu − d∇u ⇒ f · nk dS = Ic + Id
|Vi | Sk
k
Z Z
1 X 1 X
Ic = (v · nk )u dS, Id = d(nk · ∇u) dS
|Vi | Sk |Vi | Sk
k k

Interpolation is used to obtain function values at quadrature points

Convective fluxes

One-dimensional case: v = const, Ω = (0, 1)

1
∆x = , xi = i∆x, i = 0, . . . , N
N

Vertex-centered finite volume approximation

ui
ui+1
ui 1
Vi = (xi−1/2 , xi+1/2 ), |Vi | = ∆x
Vi

ui+1/2 − ui−1/2 xi 1 xi xi+1

Ic = v xi 1=2 xi+1=2
∆x

Approximation of ui±1/2 in terms of uj , |i − j| ≤ s

Upwind difference scheme (UDS)

Vertex-centered approximation
ui+1/2 − ui−1/2
Ic = v , v = const
∆x

P0 reconstruction
v>0 ui−1/2 ≈ ui−1 , ui+1/2 ≈ ui
ui 1=2 ui+1=2
ui+1 ui − ui−1
ui
Ic ≈ v BD
∆x
ui 1

x
v<0 ui−1/2 ≈ ui , ui+1/2 ≈ ui+1
xi 1 xi xi+1
xi 1=2 xi+1=2 ui+1 − ui
Ic ≈ v FD
∆x
Central difference scheme (CDS)

Vertex-centered approximation
ui+1/2 − ui−1/2
Ic = v , v = const
∆x

P1 reconstruction i −x
pi,1 (x) = ui−1 xix−x i−1
+ ui xx−x i−1
i −xi−1

ui 1=2 ui+1=2
ui−1 + ui
ui+1 ui−1/2 ≈ pi,1 (xi−1/2 ) =
ui 2
ui 1
ui + ui+1
ui+1/2 ≈ pi+1,1 (xi+1/2 ) =

x 2
xi 1 xi xi+1 ui+1 − ui−1
Ic ≈ v CD
xi 1=2 xi+1=2 2∆x
Linear upwind scheme (LUDS)

Vertex-centered approximation
ui+1/2 − ui−1/2
Ic = v , v = const
∆x

3ui−1 −ui−2
P1 reconstruction ui−1/2 ≈ 2
v>0 3ui −ui−1
ui+1/2 ≈ 2
ui 1=2 ui+1=2
v ui+1 3ui − 4ui−1 + ui−2
ui Ic ≈ v
2∆x
ui 1
ui 2 3ui −ui+1
ui−1/2 ≈

x v<0 2
3ui+1 −ui+2
ui+1/2 ≈ 2
xi 2 xi 1 xi xi+1
xi 1=2 xi+1=2
3ui − 4ui+1 + ui+2
Ic ≈ −v
2∆x
Quadratic upwind scheme (QUICK)

Vertex-centered approximation
ui+1/2 − ui−1/2
Ic = v , v = const
∆x

Quadratic Upwind Interpolation for Convective Kinematics

x − xi xi+1 − x x − xi−1 xi+1 − x

pi,2 (x) = ui−1 + ui
xi−1 − xi xi+1 − xi−1 xi − xi−1 xi+1 − xi
x − xi−1 x − xi
+ ui+1
xi+1 − xi−1 xi+1 − xi

(fitting a parabola to the data at three control volume centers)

Quadratic upwind scheme (QUICK)

3ui +6ui−1 −ui−2

ui−1/2 ≈ 8
v>0 3ui+1 +6ui −ui−1
P2 reconstruction ui+1/2 ≈ 8

ui 1=2 ui+1=2 3ui+1 + 3ui − 7ui−1 + ui−2

v ui+1 Ic ≈ v
ui 8∆x
ui 1 3ui−1 +6ui −ui+1
ui 2 ui−1/2 ≈
v<0 8

x ui+1/2 ≈ 3ui +6ui+1 −ui+2
8
xi 2 xi 1 xi xi+1
xi 1=2 xi+1=2 3ui−1 + 3ui − 7ui+1 + ui+2
Ic ≈ −v
8∆x

just marginal improvement in accuracy compared to LUDS

formal order can be determined using Taylor expansions
Error analysis

Taylor series expansions

 
∆x ∂u

(∆x)2 ∂2u
ui+1 = ui+1/2 + 2 ∂x i+1/2 + 8 ∂x2 + ...
2
 2
i+1/2
∆x ∂u

(∆x) ∂ u
ui = ui+1/2 − 2 ∂x i+1/2 + 8 ∂x2 − ...
i+1/2

Local truncation errors

 
v(∆x)2 ∂2u
vui = vui+1/2 − v∆x ∂u


2 ∂x i+1/2 + 8 ∂x2 + ...
  i+1/2
2
∂2u
v ui +u
2
i+1
= vui+1/2 + v(∆x)
8 ∂x2 + ...
i+1/2

LUDS is equivalent to the second-order one-sided difference scheme

QUICK is also second-order despite third-order flux approximation

Convective fluxes

Any second-order scheme is a linear combination of CDS and LUDS

3 1
IQU ICK = ICDS + ILU DS
4 4

High-order flux approximations can be derived using additional

points for construction of interpolation polynomials

Overall order also depends on the accuracy of the quadrature rule

for approximation of volume integrals (as in QUICK)

A high-order scheme is guaranteed to produce better results than a

low-order one only asymptotically, i.e., on sufficiently fine meshes
Diffusive fluxes

Vertex-centered approximation
d ∂u ∂u



∂x i+1/2 − d ∂x i−1/2
Id = − , d = const
∆x

Derivatives at the interfaces of control volumes

   
∂u ui − ui−1 ∂u ui+1 − ui
≈ , ≈
∂x i−1/2 ∆x ∂x i+1/2 ∆x

Second-order central difference

ui+1 − 2ui + ui−1
Id ≈ −d
(∆x)2
Finite element methods

Stationary partial differential equation

Lu = f in Ω ⊂ Rd , d ∈ {1, 2, 3}

Variational form: Find u ∈ V such that

a(u, v) = b(v) ∀v ∈ V

where V is a suitably chosen functional space, a : V × V → R is a

bilinear form and b : V → R is a linear functional, e.g.,
Z Z
a(u, v) = ∇u · ∇v dx, b(v) = f v dx
Ω Ω
Variational problems

Variational forms of partial differential equations can be derived and

discretized (i.e., converted into linear systems) in two ways:

1 by using a weighted residual formulation (Galerkin method)

2 by using an equivalent minimization problem (Ritz method)

Both approaches lead to the same discrete problem

Au = f

for the coefficients of the finite element approximation

X
uh = uj ϕj ∈ V h ⊂ V
j
Galerkin method

Let u ∈ V be the classical solution of the PDE

Lu = f in Ω

The weighted residual formulation is given by

Z
(Lu − f )v dx = 0 ∀v ∈ V
Ω

Integration by parts is usually performed to shift partial derivatives

from u to v and impose boundary conditions in a weak sense

The variational form is valid for piecewise-differentiable functions

Weak solutions

For the PDE form of a boundary value problem to be well-posed, all

functions that appear in the governing equation must be smooth

A classical solution of an n-th order PDE model must be n times

continuously differentiable and satisfy the boundary conditions

The concept of weak derivatives makes it possible to consider

generalized, piecewise-differentiable solutions

This generalization is based on integration by parts

Weak derivatives in 1D

Let u : [a, b] 7→ R and ϕ : [a, b] 7→ R be differentiable on (a, b). If

ϕ(a) = ϕ(b) = 0, then integration by parts yields
Z b Z b
0
u (x)ϕ(x)dx = − u(x)ϕ0 (x)dx
a a

A function g : (a, b) 7→ R is called weak derivative of u : (a, b) 7→ R if

Z b Z b
g(x)ϕ(x)dx = − u(x)ϕ0 (x)dx
a a

for all differentiable test functions ϕ : [a, b] 7→ R s.t. ϕ(a) = ϕ(b) = 0

Example


x if x ≥ 0
The function u(x) = |x| =
−x if x ≤ 0

is differentiable everywhere except at x = 0


1 if x > 0
The function g(x) =
−1 if x < 0

is a weak derivative of u. We have

g(x) = u0 (x) ∀x ∈ R\{0}

Integrals involving g are independent of g(0)

Integration by parts in 2D/3D

Divergence theorem
Z Z
∇ · f dx = f · nds
V S

for any differentiable vector field f

For f = ϕv we have ∇ · (ϕv) = ϕ∇ · v + v · ∇ϕ

Integration by parts in multidimensions

Z Z Z
ϕ∇ · vdx = ϕv · nds − v · ∇ϕdx
V S V
Integration by parts in 2D/3D

Let ϕ be a test function differentiable in Ω with ϕ = 0 on ∂Ω

Z Z
∂u ∂ϕ
v = ui ⇒ ϕ dx = − u dx
Ω ∂x Ω ∂x
Z Z
∂u ∂ϕ
v = uj ⇒ ϕ dx = − u dx
Ω ∂y Ω ∂y
Z Z
∂u ∂ϕ
v = uk ⇒ ϕ dx = − u dx
Ω ∂z Ω ∂z

Z Z
Vector notation ϕ∇udx = − u∇ϕdx
Ω Ω
Higher-order weak derivatives

m
P
Multiindex notation: α = (α1 , . . . , αm ), |α| := αi
i=1

∂ |α| ϕ
Dα ϕ := , ϕ : Ω ⊂ Rm → R
∂ α1 x α
1 . . . ∂ m xm

Let u : Ω ⊂ Rm 7→ R be a square integrable function

Z
|u(x)|2 dx < ∞
Ω

A function g is called weak derivative of u w.r.t. α if

Z Z
ϕ(x)g(x)dx = (−1)|α| u(x)Dα ϕ(x)dx
Ω Ω

for all sufficiently smooth test functions ϕ s.t. ϕ = 0 on ∂Ω

Strong vs. weak solutions

The classical solution of a second-order PDE must be (at least)

twice continuously differentiable

The space of functions satisfying this requirement is too small to

guarantee existence for many problems of practical interest

The method of weighted residuals yields variational formulations

which are valid under much weaker regularity assumptions

Weak solutions of PDE models reside in Sobolev spaces

Functional spaces

Notation used in calculus of variations

C 0 (Ω) space of continuous functions
C 1 (Ω) space of continuously differentiable functions
2
C (Ω) space of twice continuously differentiable functions
L2 (Ω) space of square integrable functions
1
H (Ω) Sobolev space of functions u ∈ L2 (Ω) with first
weak derivatives Dα u ∈ L2 (Ω), |α| = 1
H 2 (Ω) Sobolev space of functions u ∈ H 1 (Ω) with second
weak derivatives Dα u ∈ L2 (Ω), |α| = 2
Weighted residuals formulation

Generic transport equation

∂u
+ ∇ · (vu) − ∇ · (d∇u) − q = 0 in Ω
∂t

The left-hand side of this formula is the residual

∂u
R(u) = + ∇ · (vu) − ∇ · (d∇u) − q
∂t

The classical solution u ∈ C 2 (Ω) satisfies

Z
wR(u)dx = 0 ∀w ∈ L2 (Ω)
Ω
Initial and boundary conditions

Initial condition (unsteady problems)

u(x, 0) = u0 (x), x∈Ω

Dirichlet and flux boundary conditions

u|Γ1 = g1 , (vu − d∇u) · n|Γ2 = g2

Weak imposition of Dirichlet BC

(vu − d∇u) · n|Γ1 = (v · n)g1

Weak formulation

Let u ∈ V = {v ∈ H 2 (Ω) : v|Γ1 = g1 } satisfy

Z  
∂u
w + ∇ · (vu) − ∇ · (d∇u) − q dx = 0
Ω ∂t

for all test functions w ∈ W = {v ∈ L2 (Ω) : v|Γ1 = 0}

This is a weak formulation of the generic transport equation

The Sobolev space H 2 (Ω) is larger than the space C 2 (Ω)

It can be shown that u is the classical solution if u ∈ C 2 (Ω)

Weak formulation

Integration by parts formula for the diffusive term

Z Z Z
w∇ · (d∇u) dx = w(d∇u) · nds − ∇w · (d∇u) dx
Ω Γ Ω

Integration by parts formula for the convective term

Z Z Z
w∇ · (vu) dx = w(vu) · nds − ∇w · (vu) dx
Ω Γ Ω

Substitution of the Dirichlet and flux boundary conditions

w|Γ1 = 0, (vu − d∇u) · n|Γ2 = g2

Weak formulation

Variational problem: find u ∈ V such that

Z   Z
∂u
w − ∇w · (vu − d∇u) − wq dx + g2 ds = 0 ∀w ∈ W
Ω ∂t Γ2

Functional spaces for u and w

V = {v ∈ H 1 (Ω) : v|Γ1 = g1 }

W = {v ∈ H 1 (Ω) : v|Γ1 = 0}

Ideally suited for problems with flux boundary conditions

Conservation property

In the case Γ2 = Γ, the function w ≡ 1 belongs to the space

W = H 1 (Ω) of admissible test functions

Substitution into the weak form with integration by parts yields

Z Z Z
d
u dx = q dx − g ds
dt Ω Ω Γ

where
g = (vu − d∇u) · n on Γ

This is the integral form of the underlying conservation law