<< . .

. 36
( : 59)

. . >>

m CMl .
C ’1
Here CMk is the number of operations for an iteration with the iteration
method ¦k .
In the model problem of the Friedrichs“Keller triangulation with uniform
re¬nement we have C/(C ’ 1) = 4/3 and C3 = 2. For · = · 0 as basic
norm, ± = 2 is a typical case according to Theorem 3.37. The existence of
the constant C2 will hereby ¬nally be ensured consistently by the condition
(5.103), observing (5.100). Assuming also that the constants C1 , C2 , Pl
are “small” and the iteration method has a “small” contraction number ,
only a small number of iterations m is necessary, in the ideal case m = 1.
At least in this situation we can count on only a small increase of the
necessary e¬ort through the process of nested iterations, which provides an
“appropriate” approximation xk on all levels k of discretization.
Finally, it is to be observed that the sequence of the discrete problems
has to be de¬ned only during the process of the nested iteration. This o¬ers
the possibility to combine it with a posteriori error estimators as discussed
in Section 4.2, in order to develop a grid Tk+1 on which the discrete problem
of level k + 1 is determined, on the basis of xk as a re¬nement of Tk .
The Finite Volume Method

Finite volume methods are widely applied when di¬erential equations in
divergence form (cf. Section 0.5) or di¬erential equations involving such
di¬erential expressions (for example, parabolic di¬erential equations) are to
be solved numerically. In the class of second-order linear elliptic di¬erential
equations, expressions of the form
Lu := ’∇ · (K ∇u ’ c u) + r u = f (6.1)
are typical (cf. (0.33)), where
K : „¦ ’ Rd,d , c : „¦ ’ Rd , r, f : „¦ ’ R .
The corresponding “parabolic version” is
+ Lu = f
and will be treated in Chapter 7.
First-order partial di¬erential equations such as the classical conservation
∇ · q(u) = 0 ,
where q : R ’ Rd is a nonlinear vector ¬eld depending on u, or higher-order
partial di¬erential equations (such as the biharmonic equation (3.36)), or
even systems of partial di¬erential equations can be successfully discretized
by the ¬nite volume method.
In correspondence to the comparatively large class of problems that can
be treated by the ¬nite volume method, there are rather di¬erent sources
256 6. Finite Volume Method

1960 Forsythe and Wasow computation of neutron di¬usion
1961 Marˇuk
c computation of nuclear reactors
1971 McDonald ¬‚uid mechanics
1972 MacCormack and Paullay ¬‚uid mechanics
1973 Rizzi and Inouye ¬‚uid mechanics in 3D
1977 Samarski integro-interpolation method,
balance method
1979 Jameson ¬nite volume method
1984 Heinrich integro-balance method,
generalized ¬nite di¬erence method
1987 Bank and Rose box method
Table 6.1. Some sources of the ¬nite volume method.

originating mainly from practical applications. Some of these sources are
listed in Table 6.1. In contrast to ¬nite di¬erence or ¬nite element methods,
the theoretical understanding of the ¬nite volume method remained at an
early stage for a long time; only in recent years has essential progress been
The ¬nite volume method can be viewed as a discretization method of
its own right. It includes ideas from both ¬nite di¬erence and ¬nite element
methods. So in the literature approaches can be found that interpret it as
a “generalized ¬nite di¬erence method” or rather as a variant of the ¬nite
element method. In this chapter, we will consider only equations of the
type (6.1).

6.1 The Basic Idea of the Finite Volume Method
Now we will describe the fundamental steps in the derivation of the ¬nite
volume method. For simplicity, we restrict ourselves to the case d = 2 and
r = 0. Furthermore, we set q(u) := ’K ∇u + c u. Then equation (6.1)
∇ · q(u) = f . (6.2)
In order to obtain a ¬nite volume discretization, the domain „¦ will be
subdivided into M subdomains „¦i such that the collection of all those
subdomains forms a partition of „¦, that is:
(1) each „¦i is an open, simply connected, and polygonally bounded set
without slits,
6.1. Basics 257

(2) „¦i © „¦j = … (i = j),
(3) ∪M „¦i = „¦ .

These subdomains „¦i are called control volumes or control domains.
Without going into more detail we mention that there also exist ¬nite
volume methods with a well-de¬ned overlapping of the control volumes
(that is, condition 2 is violated).
The next step, which is in common with all ¬nite volume methods, con-
sists in integrating equation (6.2) over each control volume „¦i . After that,
Gauss™s divergence theorem is applied:

ν · q(u) dσ = i ∈ {1, . . . , M } ,
f dx ,
‚„¦i „¦i

where ν denotes the outer unit normal to ‚„¦i . By the ¬rst condition of the
partition, the boundary ‚„¦i is formed by straight-line segments “ij (j =
1, . . . , ni ), along which the normal ν|“ij =: νij is constant (see Figure 6.1).
So the line integral can be decomposed into a sum of line integrals from
which the following equation results:
νij · q(u) dσ = f dx . (6.3)
“ij „¦i

νi 2 νi 1

νi 3
νi 5

νi 4

Figure 6.1. A control volume.

Now the integrals occurring in (6.3) have to be approximated. This can
be done in very di¬erent ways, and so di¬erent ¬nal discretizations are
In general, ¬nite volume methods can be distinguished by the following
(1) the geometric shape of the control volumes „¦i ,
(2) the position of the unknowns (“problem variables”) with respect to
the control volumes,
258 6. Finite Volume Method

(3) the approximation of the boundary (line (d = 2) or surface (d = 3))
Especially the second criterion divides the ¬nite volume methods into two
large classes: the cell-centred and the cell-vertex ¬nite volume methods. In
the cell-centred methods, the unknowns are associated with the control
volumes (for example, any control volume corresponds to a function value
at some interior point (e.g., at the barycentre)). In the cell-vertex methods,
the unknowns are located at the vertices of the control volumes. Sometimes,
instead of the ¬rst-mentioned class a subdivision into two classes, the so-
called cell-centred and node-centred methods, is considered. The di¬erence
is whether the problem variables are assigned to the control volumes or,
given the problem variables, associated control volumes are de¬ned.
Example 6.1 Consider the homogeneous Dirichlet problem for the Pois-
son equation on the unit square:
’∆u = in „¦ = (0, 1)2 ,
u= 0 on ‚„¦ .

Problem variables:
Function values at the nodes ai
of a square grid with mesh width
aj3 ai aj1
aj4 Control volumes:
„¦i := {x ∈ „¦ : |x ’ ai |∞ < h }

Figure 6.2. Problem variables and control volumes in a cell-centred ¬nite volume

For an inner control volume „¦i (i.e., ai ∈ „¦), equation (6.3) takes the
’ νijk · ∇u dσ = f dx ,
“ijk „¦i

where “ijk := ‚„¦i © ‚„¦jk . A closer look at the directional derivatives shows
νij1 · ∇u = ‚1 u , νij2 · ∇u = ‚2 u ,
νij3 · ∇u = ’‚1 u , νij4 · ∇u = ’‚2 u .
i.e. they are just partial derivatives with respect to the ¬rst or the second
variable on the corresponding parts of the boundary.
6.1. Basics 259

Approximating the integrals on “ijk by means of the midpoint rule and
replacing the derivatives by di¬erence quotients, we have

4 4
ai + ajk
’ νijk · ∇u dσ ≈ ’ νijk · ∇u h
k=1 k=1
u(aj1 )’u(ai ) u(aj2 )’u(ai ) u(ai )’u(aj3 ) u(ai )’u(aj4 )
≈’ ’ ’
+ h
h h h h
= 4 u(ai ) ’ u(ajk ) .

Thus, we obtain exactly the expression that results from the application of
a ¬nite element method with continuous, piecewise linear ansatz and test
functions on a Friedrichs“Keller triangulation (cf. Figure 2.9).
Furthermore, if we approximate the integral „¦i f dx by f (ai )h2 , we see
that this term coincides with the trapeziodal rule applied to the right-hand
side of the mentioned ¬nite element formulation (cf. Lemma 2.13).
Actually, it is no accident that both discretization methods lead to the
same algebraic system. Later on we will prove a more general result to
con¬rm the above observation.
The boundary control volumes are treated as follows:
If ai ∈ ‚„¦, then parts of the boundary ‚„¦i lie on ‚„¦. At these nodes,
the Dirichlet boundary conditions already prescribe values of the unknown
function, and so there is no need to include the boundary control volumes
into the balance equations (6.3).
A detailed description for the case of ¬‚ux boundary conditions will be
given later, in Section 6.2.4; see (6.23).

Example 6.2 We consider the same boundary value problem as in
Example 6.1.

Problem variables:
Function values at the nodes ai
of a square grid with mesh width
Control volumes:
Subsquares of the grid

Figure 6.3. Problem variables and control volumes in a cell-vertex ¬nite volume
260 6. Finite Volume Method

In the interior of „¦, the resulting discretization yields a 12-point stencil
(in the terminology of ¬nite di¬erence methods).
Remark 6.3 In the ¬nite volume discretization of systems of partial dif-
ferential equations (resulting from ¬‚uid mechanics, for example), both
methods are used simultaneously for di¬erent variables; see Figure 6.4.

. . . ..

. . . . : problem variable of type 1

. . . . : problem variable of type 2

. . . .


Figure 6.4. Finite volume discretization of systems of partial di¬erential

Assets and Drawbacks of the Finite Volume Method
• Flexibility with respect to the geometry of the domain „¦ (as in ¬nite
element methods).
• Admissibility of unstructured grids (as in ¬nite element methods,
important for adaptive methods).
• Simple assembling.
• Conservation of certain laws valid for the continuous problem
(for example, conservation laws or maximum principles). This
property is important in the numerical solution of di¬erential equa-
tions with discontinuous coe¬cients or of convection-dominated
di¬usion-convection equations (see Section 6.2.4).
• Easy linearization of nonlinear problems (simpler than in ¬nite
element methods (Newton™s method)).
• Simple discretization of boundary conditions (as in ¬nite element
methods, especially a “natural” treatment of Neumann or mixed
boundary conditions).
• In principle, no restriction of the spatial dimension d of the domain
6.1. Basics 261


• Smaller ¬eld of applications in comparison with ¬nite element or ¬nite
di¬erence methods.

• Di¬culties in the design of higher order methods (no so-called p-
version available as in the ¬nite element method).

• In higher spatial dimensions (d ≥ 3), the construction of some classes
or types of control volumes may be a complex task and thus may lead
to a time-consuming assembling.

• Di¬cult mathematical analysis (stability, convergence, . . . ).

6.1 Given the boundary value problem

’(au ) = 0 in (0, 1) , u(0) = 1 , u(1) = 0 ,

with piecewise constant coe¬cients

κ± , x ∈ (0, ξ) ,
a(x) :=
x ∈ (ξ, 1) ,

where ±, κ are positive constants and ξ ∈ (0, 1) \ Q :

(a) What is the weak solution u ∈ H 1 (0, 1) of this problem?
(b) For general “smooth” coe¬cients a, the di¬erential equation is
obviously equivalent to

’au ’ a u = 0 .

Therefore, the following discretization is suggested:
ui’1 ’ 2ui + ui+1 ai+1 ’ ai’1 ui+1 ’ ui’1
’ai ’ = 0,
h 2h 2h
where an equidistant grid with the nodes xi = ih (i = 0, . . . , N + 1)
and ai := a(xi ), ui :≈ u(xi ) is used.
This discretization is also formally correct in the given situation of
discontinuous coe¬cients. Find the discrete solution (ui )N in this
(c) Under what conditions do the values ui converge to u(xi ) for h ’ 0?
262 6. Finite Volume Method

6.2 The Finite Volume Method for Linear Elliptic
Di¬erential Equations of Second Order on
Triangular Grids
In this section we will explain the development and the analysis of a ¬nite
volume method of “cell-centred” type for a model problem. Here, „¦ ‚ R2
is a bounded, simply connected domain with a polygonal boundary, but
without slits.

6.2.1 Admissible Control Volumes

The Voronoi Diagram
By {ai }i∈Λ ‚ „¦ we denote a consecutively numbered point set that includes
all vertices of „¦, where Λ is the corresponding set of indices. Typically, the
points ai are placed at those positions where the values u(ai ) of the exact
solution u are to be approximated. The convex set
„¦i := x ∈ R2 |x ’ ai | < |x ’ aj |
˜ for all j = i
is called the Voronoi polygon (or Dirichlet domain, Wigner“Seitz cell,
Thiessen polygon, . . . ). The family „¦i i∈Λ is called the Voronoi diagram
of the point set {ai }i∈Λ .

. . .
.. boundary of „¦

. ∼
.. boundary of „¦i

Figure 6.5. Voronoi diagram.

The Voronoi polygons are convex, but not necessarily bounded, sets (con-
sider the situation near the boundary in Figure 6.5). Their boundaries are
polygons. The vertices of these polygons are called Voronoi vertices.
It can be shown that at any Voronoi vertex at least three Voronoi poly-
gons meet. According to this property, Voronoi vertices are classi¬ed into
regular and degenerate Voronoi vertices: In a regular Voronoi vertex, the
boundaries of exactly three Voronoi polygons meet, whereas a degenerate
Voronoi vertex is shared by at least four Voronoi polygons. In the latter
case, all the corresponding nodes ai are located at some circle (they are

<< . .

. 36
( : 59)

. . >>