ńņš. 1 |

Elliptic and Parabolic

Partial Differential

Equations

Peter Knabner

Lutz Angermann

Springer

44

Texts in Applied Mathematics

Editors

J.E. Marsden

L. Sirovich

S.S. Antman

Advisors

G. Iooss

P. Holmes

D. Barkley

M. Dellnitz

P. Newton

This page intentionally left blank

Peter Knabner Lutz Angermann

Numerical Methods for

Elliptic and Parabolic Partial

Differential Equations

With 67 Figures

Peter Knabner Lutz Angermann

Institute for Applied Mathematics Institute for Mathematics

University of Erlangen University of Clausthal

Martensstrasse 3 Erzstrasse 1

D-91058 Erlangen D-38678 Clausthal-Zellerfeld

Germany Germany

knabner@am.uni-erlangen.de angermann@math.tu-clausthal.de

Series Editors

J.E. Marsden L. Sirovich

Control and Dynamical Systems, 107ā“81 Division of Applied Mathematics

California Institute of Technology Brown University

Pasadena, CA 91125 Providence, RI 02912

USA USA

marsden@cds.caltech.edu chico@camelot.mssm.edu

S.S. Antman

Department of Mathematics

and

Institute for Physical Science

and Technology

University of Maryland

College Park, MD 20742-4015

USA

ssa@math.umd.edu

Mathematics Subject Classification (2000): 65Nxx, 65Mxx, 65F10, 65H10

Library of Congress Cataloging-in-Publication Data

Knabner, Peter.

[Numerik partieller Differentialgleichungen. English]

Numerical methods for elliptic and parabolic partial differential equations /

Peter Knabner, Lutz Angermann.

p. cm. ā” (Texts in applied mathematics ; 44)

Include bibliographical references and index.

ISBN 0-387-95449-X (alk. paper)

1. Differential equations, Partialā”Numerical solutions. I. Angermann, Lutz. II. Title.

III. Series.

QA377.K575 2003

515ā².353ā”dc21 2002044522

ISBN 0-387-95449-X Printed on acid-free paper.

ļ™ 2003 Springer-Verlag New York, Inc.

All rights reserved. This work may not be translated or copied in whole or in part without the written

permission of the publisher (Springer-Verlag New York, Inc., 175 Fifth Avenue, New York, NY 10010,

USA), except for brief excerpts in connection with reviews or scholarly analysis. Use in connection with

any form of information storage and retrieval, electronic adaptation, computer software, or by similar or

dissimilar methodology now known or hereafter developed is forbidden.

The use in this publication of trade names, trademarks, service marks, and similar terms, even if they are

not identified as such, is not to be taken as an expression of opinion as to whether or not they are subject

to proprietary rights.

Printed in the United States of America.

987654321 SPIN 10867187

Typesetting: Pages created by the authors in 2e using Springerā™s svsing6.cls macro.

www.springer-ny.com

Springer-Verlag New York Berlin Heidelberg

A member of BertelsmannSpringer Science+Business Media GmbH

Series Preface

Mathematics is playing an ever more important role in the physical and

biological sciences, provoking a blurring of boundaries between scientific

disciplines and a resurgence of interest in the modern as well as the classical

techniques of applied mathematics. This renewal of interest, both in re-

search and teaching, has led to the establishment of the series Texts in

Applied Mathematics (TAM).

The development of new courses is a natural consequence of a high level

of excitement on the research frontier as newer techniques, such as numeri-

cal and symbolic computer systems, dynamical systems, and chaos, mix

with and reinforce the traditional methods of applied mathematics. Thus,

the purpose of this textbook series is to meet the current and future needs

of these advances and to encourage the teaching of new courses.

TAM will publish textbooks suitable for use in advanced undergraduate

and beginning graduate courses, and will complement the Applied Mathe-

matical Sciences (AMS) series, which will focus on advanced textbooks and

research-level monographs.

Pasadena, California J.E. Marsden

Providence, Rhode Island L. Sirovich

College Park, Maryland S.S. Antman

This page intentionally left blank

Preface to the English Edition

Shortly after the appearance of the German edition we were asked by

Springer to create an English version of our book, and we gratefully ac-

cepted. We took this opportunity not only to correct some misprints and

mistakes that have come to our knowledge1 but also to extend the text at

various places. This mainly concerns the role of the ļ¬nite diļ¬erence and

the ļ¬nite volume methods, which have gained more attention by a slight

extension of Chapters 1 and 6 and by a considerable extension of Chapter

7. Time-dependent problems are now treated with all three approaches (ļ¬-

nite diļ¬erences, ļ¬nite elements, and ļ¬nite volumes), doing this in a uniform

way as far as possible. This also made a reordering of Chapters 6ā“8 nec-

essary. Also, the index has been enlarged. To improve the direct usability

in courses, exercises now follow each section and should provide enough

material for homework.

This new version of the book would not have come into existence without

our already mentioned team of helpers, who also carried out ļ¬rst versions

of translations of parts of the book. Beyond those already mentioned, the

team was enforced by Cecilia David, Basca Jadamba, Dr. Serge KrĀØutle, a

Dr. Wilhelm Merz, and Peter Mirsch. Alexander Prechtel now took charge

of the diļ¬cult modiļ¬cation process. Prof. Paul DuChateau suggested im-

provements. We want to extend our gratitude to all of them. Finally, we

1 Users

of the German edition may consult

http:/

/www.math.tu-clausthal.de/Ėmala/publications/errata.pdf

viii Preface to the English Edition

thank senior editor Achi Dosanjh, from Springer-Verlag New York, Inc., for

her constant encouragement.

Remarks for the Reader and the Use in Lectures

The size of the text corresponds roughly to four hours of lectures per week

over two terms. If the course lasts only one term, then a selection is nec-

essary, which should be orientated to the audience. We recommend the

following ācutsā:

Chapter 0 may be skipped if the partial diļ¬erential equations treated

therein are familiar. Section 0.5 should be consulted because of the notation

collected there. The same is true for Chapter 1; possibly Section 1.4 may

be integrated into Chapter 3 if one wants to deal with Section 3.9 or with

Section 7.5.

Chapters 2 and 3 are the core of the book. The inductive presenta-

tion that we preferred for some theoretical aspects may be shortened for

students of mathematics. To the lecturerā™s taste and depending on the

knowledge of the audience in numerical mathematics Section 2.5 may be

skipped. This might impede the treatment of the ILU preconditioning in

Section 5.3. Observe that in Sections 2.1ā“2.3 the treatment of the model

problem is merged with basic abstract statements. Skipping the treatment

of the model problem, in turn, requires an integration of these statements

into Chapter 3. In doing so Section 2.4 may be easily combined with Sec-

tion 3.5. In Chapter 3 the theoretical kernel consists of Sections 3.1, 3.2.1,

3.3ā“3.4.

Chapter 4 presents an overview of its subject, not a detailed development,

and is an extension of the classical subjects, as are Chapters 6 and 9 and

the related parts of Chapter 7.

In the extensive Chapter 5 one might focus on special subjects or just con-

sider Sections 5.2, 5.3 (and 5.4) in order to present at least one practically

relevant and modern iterative method.

Section 8.1 and the ļ¬rst part of Section 8.2 contain basic knowledge of

numerical mathematics and, depending on the audience, may be omitted.

The appendices are meant only for consultation and may complete

the basic lectures, such as in analysis, linear algebra, and advanced

mathematics for engineers.

Concerning related textbooks for supplementary use, to the best of our

knowledge there is none covering approximately the same topics. Quite a

few deal with ļ¬nite element methods, and the closest one in spirit probably

is [21], but also [6] or [7] have a certain overlap, and also oļ¬er additional

material not covered here. From the books specialised in ļ¬nite diļ¬erence

methods, we mention [32] as an example. The (node-oriented) ļ¬nite volume

method is popular in engineering, in particular in ļ¬‚uid dynamics, but to

the best of our knowledge there is no presentation similar to ours in a

Preface to the English Edition ix

mathematical textbook. References to textbooks specialised in the topics

of Chapters 4, 5 and 8 are given there.

Remarks on the Notation

Printing in italics emphasizes deļ¬nitions of notation, even if this is not

carried out as a numbered deļ¬nition.

Vectors appear in diļ¬erent forms: Besides the āshortā space vectors

x ā Rd there are ālongā representation vectors u ā Rm , which describe

in general the degrees of freedom of a ļ¬nite element (or volume) approxi-

mation or represent the values on grid points of a ļ¬nite diļ¬erence method.

Here we choose bold type, also in order to have a distinctive feature from

the generated functions, which frequently have the same notation, or from

the grid functions.

Deviations can be found in Chapter 0, where vectorial quantities belong-

ing to Rd are boldly typed, and in Chapters 5 and 8, where the unknowns

of linear and nonlinear systems of equations, which are treated in a general

manner there, are denoted by x ā Rm .

Components of vectors will be designated by a subindex, creating a

double index for indexed quantities. Sequences of vectors will be supplied

with a superindex (in parentheses); only in an abstract setting do we use

subindices.

Erlangen, Germany Peter Knabner

Clausthal-Zellerfeld, Germany Lutz Angermann

January 2002

This page intentionally left blank

Preface to the German Edition

This book resulted from lectures given at the University of Erlangenā“

Nuremberg and at the University of Magdeburg. On these occasions we

often had to deal with the problem of a heterogeneous audience composed

of students of mathematics and of diļ¬erent natural or engineering sciences.

Thus the expectations of the students concerning the mathematical accu-

racy and the applicability of the results were widely spread. On the other

hand, neither relevant models of partial diļ¬erential equations nor some

knowledge of the (modern) theory of partial diļ¬erential equations could be

assumed among the whole audience. Consequently, in order to overcome the

given situation, we have chosen a selection of models and methods relevant

for applications (which might be extended) and attempted to illuminate the

whole spectrum, extending from the theory to the implementation, with-

out assuming advanced mathematical background. Most of the theoretical

obstacles, diļ¬cult for nonmathematicians, will be treated in an āinduc-

tiveā manner. In general, we use an explanatory style without (hopefully)

compromising the mathematical accuracy.

We hope to supply especially students of mathematics with the in-

formation necessary for the comprehension and implementation of ļ¬nite

element/ļ¬nite volume methods. For students of the various natural or

engineering sciences the text oļ¬ers, beyond the possibly already existing

knowledge concerning the application of the methods in special ļ¬elds, an

introduction into the mathematical foundations, which should facilitate the

transformation of speciļ¬c knowledge to other ļ¬elds of applications.

We want to express our gratitude for the valuable help that we received

during the writing of this book: Dr. Markus Bause, Sandro Bitterlich,

xii Preface to the German Edition

Dr. Christof Eck, Alexander Prechtel, Joachim Rang, and Dr. Eckhard

Schneid did the proofreading and suggested important improvements. From

the anonymous referees we received useful comments. Very special thanks

go to Mrs. Magdalena Ihle and Dr. Gerhard Summ. Mrs. Ihle transposed

the text quickly and precisely into TEX. Dr. Summ not only worked on the

original script and on the TEX-form, he also organized the complex and

distributed rewriting and extension procedure. The elimination of many

inconsistencies is due to him. Additionally he inļ¬‚uenced parts of Sec-

tions 3.4 and 3.8 by his outstanding diploma thesis. We also want to thank

Dr. Chistoph Tapp for the preparation of the graphic of the title and for

providing other graphics from his doctoral thesis [70].

Of course, hints concerning (typing) mistakes and general improvements

are always welcome.

We thank Springer-Verlag for their constructive collaboration.

Last, but not least, we want to express our gratitude to our families for

their understanding and forbearance, which were necessary for us especially

during the last months of writing.

Erlangen, Germany Peter Knabner

Magdeburg, Germany Lutz Angermann

February 2000

Contents

Series Preface v

Preface to the English Edition vii

Preface to the German Edition xi

0 For Example: Modelling Processes in Porous Media

with Diļ¬erential Equations 1

0.1 The Basic Partial Diļ¬erential Equation Models . . . . . 1

0.2 Reactions and Transport in Porous Media . . . . . . . . 5

0.3 Fluid Flow in Porous Media . . . . . . . . . . . . . . . . 7

0.4 Reactive Solute Transport in Porous Media . . . . . . . . 11

0.5 Boundary and Initial Value Problems . . . . . . . . . . . 14

1 For the Beginning: The Finite Diļ¬erence Method

for the Poisson Equation 19

1.1 The Dirichlet Problem for the Poisson Equation . . ... 19

1.2 The Finite Diļ¬erence Method . . . . . . . . . . . . ... 21

1.3 Generalizations and Limitations

of the Finite Diļ¬erence Method . . . . . . . . . . . ... 29

1.4 Maximum Principles and Stability . . . . . . . . . . ... 36

2 The Finite Element Method for the Poisson Equation 46

2.1 Variational Formulation for the Model Problem . . . . . 46

xiv Contents

2.2 The Finite Element Method with Linear Elements . ... 55

2.3 Stability and Convergence of the

Finite Element Method . . . . . . . . . . . . . . . . ... 68

2.4 The Implementation of the Finite Element Method:

Part 1 . . . . . . . . . . . . . . . . . . . . . . . . . ... 74

2.5 Solving Sparse Systems of Linear Equations

by Direct Methods . . . . . . . . . . . . . . . . . . ... 82

3 The Finite Element Method for Linear Elliptic

Boundary Value Problems of Second Order 92

3.1 Variational Equations and Sobolev Spaces . . . . . . .. 92

3.2 Elliptic Boundary Value Problems of Second Order . .. 100

3.3 Element Types and Aļ¬ne

Equivalent Triangulations . . . . . . . . . . . . . . . .. 114

3.4 Convergence Rate Estimates . . . . . . . . . . . . . . .. 131

3.5 The Implementation of the Finite Element Method:

Part 2 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 148

3.6 Convergence Rate Results in Case of

Quadrature and Interpolation . . . . . . . . . . . . . . . 155

3.7 The Condition Number of Finite Element Matrices . . . 163

3.8 General Domains and Isoparametric Elements . . . . . . 167

3.9 The Maximum Principle for Finite Element Methods . . 171

4 Grid Generation and A Posteriori Error Estimation 176

4.1 Grid Generation . . . . . . . . . . . . . . . . . . . . . . . 176

4.2 A Posteriori Error Estimates and Grid Adaptation . . . 185

5 Iterative Methods for Systems of Linear Equations 198

5.1 Linear Stationary Iterative Methods . . . . . . . . . ... 200

5.2 Gradient and Conjugate Gradient Methods . . . . . ... 217

5.3 Preconditioned Conjugate Gradient Method . . . . ... 227

5.4 Krylov Subspace Methods

for Nonsymmetric Systems of Equations . . . . . . ... 233

5.5 The Multigrid Method . . . . . . . . . . . . . . . . ... 238

5.6 Nested Iterations . . . . . . . . . . . . . . . . . . . ... 251

6 The Finite Volume Method 255

6.1 The Basic Idea of the Finite Volume Method . . . . . . . 256

6.2 The Finite Volume Method for Linear Elliptic Diļ¬eren-

tial Equations of Second Order on Triangular Grids . . . 262

7 Discretization Methods for Parabolic Initial Boundary

Value Problems 283

7.1 Problem Setting and Solution Concept . . . . . . . . . . 283

7.2 Semidiscretization by the Vertical Method of Lines . . . 293

Contents xv

7.3 Fully Discrete Schemes . . . . . .............. 311

7.4 Stability . . . . . . . . . . . . . .............. 315

7.5 The Maximum Principle for the

One-Step-Theta Method . . . . .............. 323

7.6 Order of Convergence Estimates .............. 330

8 Iterative Methods for Nonlinear Equations 342

8.1 Fixed-Point Iterations . . . . . . . . . . . . . . . . . . . . 344

8.2 Newtonā™s Method and Its Variants . . . . . . . . . . . . 348

8.3 Semilinear Boundary Value Problems for Elliptic

and Parabolic Equations . . . . . . . . . . . . . . . . . . 360

9 Discretization Methods

for Convection-Dominated Problems 368

9.1 Standard Methods and

Convection-Dominated Problems . . . . . . . . . . . . . 368

9.2 The Streamline-Diļ¬usion Method . . . . . . . . . . . . . 375

9.3 Finite Volume Methods . . . . . . . . . . . . . . . . . . . 383

9.4 The Lagrangeā“Galerkin Method . . . . . . . . . . . . . . 387

A Appendices 390

A.1 Notation . . . . . . . . . . . . . . . . . . . ........ 390

A.2 Basic Concepts of Analysis . . . . . . . . . ........ 393

A.3 Basic Concepts of Linear Algebra . . . . . ........ 394

A.4 Some Deļ¬nitions and Arguments of Linear

Functional Analysis . . . . . . . . . . . . . ........ 399

A.5 Function Spaces . . . . . . . . . . . . . . . ........ 404

References: Textbooks and Monographs 409

References: Journal Papers 412

Index 415

This page intentionally left blank

0

For Example:

Modelling Processes in Porous

Media with Diļ¬erential Equations

This chapter illustrates the scientiļ¬c context in which diļ¬erential equation

models may occur, in general, and also in a speciļ¬c example. Section 0.1

reviews the fundamental equations, for some of them discretization tech-

niques will be developed and investigated in this book. In Sections 0.2 ā“

0.4 we focus on reaction and transport processes in porous media. These

sections are independent of the remaining parts and may be skipped by

the reader. Section 0.5, however, should be consulted because it ļ¬xes some

notation to be used later on.

0.1 The Basic Partial Diļ¬erential Equation Models

Partial diļ¬erential equations are equations involving some partial deriva-

tives of an unknown function u in several independent variables. Partial

diļ¬erential equations which arise from the modelling of spatial (and tempo-

ral) processes in nature or technology are of particular interest. Therefore,

we assume that the variables of u are x = (x1 , . . . , xd )T ā Rd for d ā„ 1,

representing a spatial point, and possibly t ā R, representing time. Thus

the minimal set of variables is (x1 , x2 ) or (x1 , t), otherwise we have ordinary

diļ¬erential equations. We will assume that x ā ā„¦, where ā„¦ is a bounded

domain, e.g., a metal workpiece, or a groundwater aquifer, and t ā (0, T ] for

some (time horizon) T > 0. Nevertheless also processes acting in the whole

Rd Ć— R, or in unbounded subsets of it, are of interest. One may consult the

Appendix for notations from analysis etc. used here. Often the function u

2 0. Modelling Processes in Porous Media with Diļ¬erential Equations

represents, or is related to, the volume density of an extensive quantity like

mass, energy, or momentum, which is conserved. In their original form all

quantities have dimensions that we denote in accordance with the Inter-

national System of Units (SI) and write in square brackets [ ]. Let a be

a symbol for the unit of the extensive quantity, then its volume density

is assumed to have the form S = S(u), i.e., the unit of S(u) is a/m3 . For

example, for mass conservation a = kg, and S(u) is a concentration. For

describing the conservation we consider an arbitrary ānot too badā sub-

set ā„¦ ā‚ ā„¦, the control volume. The time variation of the total extensive

Ė

Ė

quantity in ā„¦ is then

ā‚t S(u(x, t))dx . (0.1)

Ė

ā„¦

If this function does not vanish, only two reasons are possible due to con-

servation:

ā” There is an internally distributed source density Q = Q(x, t, u) [a/m3 /s],

being positive if S(u) is produced, and negative if it is destroyed, i.e., one

term to balance (0.1) is ā„¦ Q(x, t, u(x, t))dx.

Ė

Ė

ā” There is a net ļ¬‚ux of the extensive quantity over the boundary ā‚ ā„¦ of

ā„¦. Let J = J (x, t) [a/m2 /s] denote the ļ¬‚ux density, i.e., J i is the amount,

Ė

that passes a unit square perpendicular to the ith axis in one second in

the direction of the ith axis (if positive), and in the opposite direction

otherwise. Then another term to balance (0.1) is given by

ā’ J (x, t) Ā· Ī½(x)dĻ ,

ā‚ā„¦

where Ī½ denotes the outer unit normal on ā‚ā„¦. Summarizing the conserva-

tion reads

S(u(x, t))dx = ā’ J (x, t) Ā· Ī½(x)dĻ +

ā‚t Q(x, t, u(x, t))dx . (0.2)

Ė Ė Ė

ā„¦ ā‚ā„¦ ā„¦

The integral theorem of Gauss (see (2.3)) and an exchange of time

derivative and integral leads to

[ā‚t S(u(x, t)) + ā Ā· J (x, t) ā’ Q(x, t, u(x, t))]dx = 0 ,

Ė

ā„¦

Ė

and, as ā„¦ is arbitrary, also to

ā‚t S(u(x, t)) + ā Ā· J(x, t) = Q(x, t, u(x, t)) for x ā ā„¦, t ā (0, T ] . (0.3)

All manipulations here are formal assuming that the functions involved

have the necessary properties. The partial diļ¬erential equation (0.3) is the

basic pointwise conservation equation, (0.2) its corresponding integral form.

Equation (0.3) is one requirement for the two unknowns u and J , thus it

0.1. The Basic Partial Diļ¬erential Equation Models 3

has to be closed by a (phenomenological) constitutive law, postulating a

relation between J and u.

Assume ā„¦ is a container ļ¬lled with a ļ¬‚uid in which a substance is dis-

solved. If u is the concentration of this substance, then S(u) = u and a

= kg. The description of J depends on the processes involved. If the ļ¬‚uid

is at rest, then ļ¬‚ux is only possible due to molecular diļ¬usion, i.e., a ļ¬‚ux

from high to low concentrations due to random motion of the dissolved

particles. Experimental evidence leads to

J (1) = ā’Kāu (0.4)

with a parameter K > 0 [m2 /s], the molecular diļ¬usivity. Equation (0.4)

is called Fickā™s law.

In other situations, like heat conduction in a solid, a similar model occurs.

Here, u represents the temperature, and the underlying principle is energy

conservation. The constitutive law is Fourierā™s law, which also has the form

(0.4), but as K is a material parameter, it may vary with space or, for

anisotropic materials, be a matrix instead of a scalar.

Thus we obtain the diļ¬usion equation

ā‚t u ā’ ā Ā· (Kāu) = Q . (0.5)

If K is scalar and constant ā” let K = 1 by scaling ā”, and f := Q is

independent of u, the equation simpliļ¬es further to

ā‚t u ā’ āu = f ,

where āu := āĀ·(āu) . We mentioned already that this equation also occurs

in the modelling of heat conduction, therefore this equation or (0.5) is also

called the heat equation.

If the ļ¬‚uid is in motion with a (given) velocity c then (forced) convection

of the particles takes place, being described by

J (2) = uc , (0.6)

i.e., taking both processes into account, the model takes the form of the

convection-diļ¬usion equation

ā‚t u ā’ ā Ā· (Kāu ā’ cu) = Q . (0.7)

The relative strength of the two processes is measured by the PĀ“clete

number (deļ¬ned in Section 0.4). If convection is dominating one may ignore

diļ¬usion and only consider the transport equation

ā‚t u + ā Ā· (cu) = Q . (0.8)

The diļ¬erent nature of the two processes has to be reļ¬‚ected in the models,

therefore, adapted discretization techniques will be necessary. In this book

we will consider models like (0.7), usually with a signiļ¬cant contribution

of diļ¬usion, and the case of dominating convection is studied in Chapter

9. The pure convective case like (0.8) will not be treated.

4 0. Modelling Processes in Porous Media with Diļ¬erential Equations

In more general versions of (0.7) ā‚t u is replaced by ā‚t S(u), where S

depends linearly or nonlinearly on u. In the case of heat conduction S is

the internal energy density, which is related to the temperature u via the

factors mass density and speciļ¬c heat. For some materials the speciļ¬c heat

depends on the temperature, then S is a nonlinear function of u.

Further aspects come into play by the source term Q if it depends linearly

or nonlinearly on u, in particular due to (chemical) reactions. Examples for

these cases will be developed in the following sections. Since equation (0.3)

and its examples describe conservation in general, it still has to be adapted

to a concrete situation to ensure a unique solution u. This is done by the

speciļ¬cation of an initial condition

for x ā ā„¦ ,

S(u(x, 0)) = S0 (x)

and by boundary conditions. In the example of the water ļ¬lled container

no mass ļ¬‚ux will occur across its walls, therefore, the following boundary

condition

J Ā· Ī½(x, t) = 0 for x ā ā‚ā„¦, t ā (0, T ) (0.9)

is appropriate, which ā” depending on the deļ¬nition of J ā” prescribes the

normal derivative of u, or a linear combination of it and u. In Section 0.5

additional situations are depicted.

If a process is stationary, i.e. time-independent, then equation (0.3)

reduces to

ā Ā· J (x) = Q(x, u(x)) for x ā ā„¦ ,

which in the case of diļ¬usion and convection is speciļ¬ed to

ā’ā Ā· (Kāu ā’ cu) = Q .

For constant K ā” let K = 1 by scaling ā”, c = 0, and f := Q, being

independent of u, this equation reduces to

ā’āu = f in ā„¦ ,

the Poisson equation.

Instead of the boundary condition (0.9), one can prescribe the values of

the function u at the boundary:

for x ā ā‚ā„¦ .

u(x) = g(x)

For models , where u is a concentration or temperature, the physical reali-

sation of such a boundary condition may raise questions, but in mechanical

models, where u is to interpreted as a displacement, such a boundary con-

dition seems reasonable. The last boundary value problem will be the ļ¬rst

model, whose discretization will be discussed in Chapters 1 and 2.

Finally it should be noted that it is advisable to non-dimensionalise the

ļ¬nal model before numerical methods are applied. This means that both

the independent variables xi (and t), and the dependent one u, are replaced

0.2. Reactions and Transport in Porous Media 5

by xi /xi,ref , t/tref , and u/uref , where xi,ref , tref , and uref are ļ¬xed reference

values of the same dimension as xi , t, and u, respectively. These reference

values are considered to be of typical size for the problems under investiga-

tion. This procedure has two advantages: On the one hand, the typical size

is now 1, such that there is an absolute scale for (an error in) a quantity

to be small or large. On the other hand, if the reference values are chosen

appropriately a reduction in the number of equation parameters like K

and c in (0.7) might be possible, having only fewer algebraic expressions of

the original material parameters in the equation. This facilitates numerical

parameter studies.

0.2 Reactions and Transport in Porous Media

A porous medium is a heterogeneous material consisting of a solid matrix

and a pore space contained therein. We consider the pore space (of the

porous medium) as connected; otherwise, the transport of ļ¬‚uids in the

pore space would not be possible. Porous media occur in nature and man-

ufactured materials. Soils and aquifers are examples in geosciences; porous

catalysts, chromatographic columns, and ceramic foams play important

roles in chemical engineering. Even the human skin can be considered a

porous medium. In the following we focus on applications in the geosciences.

Thus we use a terminology referring to the natural soil as a porous medium.

On the micro or pore scale of a single grain or pore, i.e., in a range of Āµm

to mm, the ļ¬‚uids constitute diļ¬erent phases in the thermodynamic sense.

Thus we name this system in the case of k ļ¬‚uids including the solid matrix

as (k + 1)-phase system or we speak of k-phase ļ¬‚ow.

We distinguish three classes of ļ¬‚uids with diļ¬erent aļ¬nities to the solid

matrix. These are an aqueous phase, marked with the index āwā for water,

a nonaqueous phase liquid (like oil or gasoline as natural resources or con-

taminants), marked with the index āo,ā and a gaseous phase, marked with

the index āgā (e.g., soil air). Locally, at least one of these phases has al-

ways to be present; during a transient process phases can locally disappear

or be generated. These ļ¬‚uid phases are in turn mixtures of several com-

ponents. In applications of the earth sciences, for example, we do not deal

with pure water but encounter diļ¬erent species in true or colloidal solu-

tion in the solvent water. The wide range of chemical components includes

plant nutrients, mineral nutrients from salt domes, organic decomposition

products, and various organic and inorganic chemicals. These substances

are normally not inert, but are subject to reactions and transformation

processes. Along with diļ¬usion, forced convection induced by the motion

of the ļ¬‚uid is the essential driving mechanism for the transport of solutes.

But we also encounter natural convection by the coupling of the dynamics

of the substance to the ļ¬‚uid ļ¬‚ow. The description level at the microscale

6 0. Modelling Processes in Porous Media with Diļ¬erential Equations

that we have used so far is not suitable for processes at the laboratory or

technical scale, which take place in ranges of cm to m, or even for processes

in a catchment area with units of km. For those macroscales new models

have to be developed, which emerge from averaging procedures of the mod-

els on the microscale. There may also exist principal diļ¬erences among the

various macroscales that let us expect diļ¬erent models, which arise from

each other by upscaling. But this aspect will not be investigated here fur-

ther. For the transition of micro to macro scales the engineering sciences

provide the heuristic method of volume averaging, and mathematics the

rigorous (but of only limited use) approach of homogenization (see [36] or

[19]). None of the two possibilities can be depicted here completely. Where

necessary we will refer to volume averaging for (heuristic) motivation.

Let ā„¦ ā‚ Rd be the domain of interest. All subsequent considerations are

formal in the sense that the admissibility of the analytic manipulations is

supposed. This can be achieved by the assumption of suļ¬cient smoothness

for the corresponding functions and domains.

Let V ā‚ ā„¦ be an admissible representative elementary volume in the

sense of volume averaging around a point x ā ā„¦. Typically the shape and

the size of a representative elementary volume are selected in such a manner

that the averaged values of all geometric characteristics of the microstruc-

ture of the pore space are independent of the size of V but depend on

the location of the point x. Then we obtain for a given variable ĻĪ± in the

phase Ī± (after continuation of ĻĪ± with 0 outside of Ī±) the corresponding

macroscopic quantities, assigned to the location x, as the extrinsic phase

average

1

ĻĪ± := ĻĪ±

|V | V

or as the intrinsic phase average

1

Ī±

ĻĪ± := ĻĪ± .

|VĪ± | VĪ±

Here VĪ± denotes the subset of V corresponding to Ī±. Let t ā (0, T ) be

the time at which the process is observed. The notation x ā ā„¦ means the

vector in Cartesian coordinates, whose coordinates are referred to by x,

y, and z ā R. Despite this ambiguity the meaning can always be clearly

derived from the context.

Let the index āsā (for solid) stand for the solid phase; then

Ļ(x) := |V \ Vs | |V | > 0

denotes the porosity, and for every liquid phase Ī±,

SĪ± (x, t) := |VĪ± | |V \ Vs | ā„ 0

0.3. Fluid Flow in Porous Media 7

is the saturation of the phase Ī±. Here we suppose that the solid phase is

stable and immobile. Thus

Ī±

ĻĪ± = ĻSĪ± ĻĪ±

for a ļ¬‚uid phase Ī± and

SĪ± = 1 . (0.10)

Ī±:ļ¬‚uid

So if the ļ¬‚uid phases are immiscible on the micro scale, they may be miscible

on the macro scale, and the immiscibility on the macro scale is an additional

assumption for the model.

As in other disciplines the diļ¬erential equation models are derived here

from conservation laws for the extensive quantities mass, impulse, and en-

ergy, supplemented by constitutive relationships, where we want to focus

on the mass.

0.3 Fluid Flow in Porous Media

Consider a liquid phase Ī± on the micro scale. In this chapter, for clarity, we

write āshortā vectors in Rd also in bold with the exception of the coordinate

vector x. Let ĖĪ± [kg/m3 ] be the (microscopic) density, q Ī± :=

Ė Ė

Ī· ĖĪ· v Ī· ĖĪ±

Ė

[m/s] the mass average mixture velocity based on the particle velocity v Ī· of

3

a component Ī· and its concentration in solution ĖĪ· [kg/m ]. The transport

theorem of Reynolds (see, for example, [10]) leads to the mass conservation

law

Ė

ā‚t ĖĪ± + ā Ā· (ĖĪ± q Ī± ) = fĪ±

Ė (0.11)

Ė

with a distributed mass source density fĪ± . By averaging we obtain from

here the mass conservation law

+āĀ·(

ā‚t (ĻSĪ± Ī±) Ī±qĪ±) = fĪ± (0.12)

with Ī± , the density of phase Ī±, as the intrinsic phase average of ĖĪ± and

q Ī± , the volumetric ļ¬‚uid velocity or Darcy velocity of the phase Ī±, as the

Ė

extrinsic phase average of q Ī± . Correspondingly, fĪ± is an average mass source

density.

Before we proceed in the general discussion, we want to consider some

speciļ¬c situations: The area between the groundwater table and the imper-

meable body of an aquifer is characterized by the fact that the whole pore

space is occupied by a ļ¬‚uid phase, the soil water. The corresponding satu-

ration thus equals 1 everywhere, and with omission of the index equation

(0.12) takes the form

ā‚t (Ļ ) + ā Ā· ( q) = f . (0.13)

8 0. Modelling Processes in Porous Media with Diļ¬erential Equations

If the density of water is assumed to be constant, due to neglecting

the mass of solutes and compressibility of water, equation (0.13) simpliļ¬es

further to the stationary equation

āĀ·q =f , (0.14)

where f has been replaced by the volume source density f / , keeping the

same notation. This equation will be completed by a relationship that

can be interpreted as the macroscopic analogue of the conservation of mo-

mentum, but should be accounted here only as an experimentally derived

constitutive relationship. This relationship is called Darcyā™s law, which

reads as

q = ā’K (āp + gez ) (0.15)

and can be applied in the range of laminar ļ¬‚ow. Here p [N/m2 ] is the intrinsic

average of the water pressure, g [m/s2 ] the gravitational acceleration, ez the

unit vector in the z-direction oriented against the gravitation,

K = k/Āµ , (0.16)

a quantity, which is given by the permeability k determined by the solid

phase, and the viscosity Āµ determined by the ļ¬‚uid phase. For an anisotropic

solid, the matrix k = k(x) is a symmetric positive deļ¬nite matrix.

Inserting (0.15) in (0.14) and replacing K by K g, known as hydraulic

conductivity in the literature, and keeping the same notation gives the

following linear equation for

1

h(x, t) := p(x, t) + z ,

g

the piezometric head h [m]:

ā’ā Ā· (Kāh) = f . (0.17)

The resulting equation is stationary and linear. We call a diļ¬erential equa-

tion model stationary if it depends only on the location x and not on the

time t, and instationary otherwise. A diļ¬erential equation and correspond-

ing boundary conditions (cf. Section 0.5) are called linear if the sum or a

scalar multiple of a solution again forms a solution for the sum, respectively

the scalar multiple, of the sources.

If we deal with an isotropic solid matrix, we have K = KI with the dĆ— d

unit matrix I and a scalar function K. Equation (0.17) in this case reads

ā’ā Ā· (Kāh) = f . (0.18)

Finally if the solid matrix is homogeneous, i.e., K is constant, we get from

division by K and maintaining the notation f the Poisson equation

ā’āh = f , (0.19)

0.3. Fluid Flow in Porous Media 9

which is termed the Laplace equation for f = 0. This model and its more

general formulations occur in various contexts. If, contrary to the above as-

sumption, the solid matrix is compressible under the pressure of the water,

and if we suppose (0.13) to be valid, then we can establish a relationship

Ļ = Ļ(x, t) = Ļ0 (x)Ļf (p)

with Ļ0 (x) > 0 and a monotone increasing Ļf such that with S(p) := Ļf (p)

we get the equation

Ļ0 S(p) ā‚t p + ā Ā· q = f

and the instationary equations corresponding to (0.17)ā“(0.19), respectively.

For constant S(p) > 0 this yields the following linear equation:

Ļ0 S ā‚t h ā’ ā Ā· (Kāh) = f , (0.20)

which also represents a common model in many contexts and is known from

corresponding ļ¬elds of application as the heat conduction equation.

We consider single phase ļ¬‚ow further, but now we will consider gas as

ļ¬‚uid phase. Because of the compressibility, the density is a function of the

pressure, which is invertible due to its strict monotonicity to

p = P( ) .

Together with (0.13) and (0.15) we get a nonlinear variant of the heat

conduction equation in the unknown :

ā‚t (Ļ ) ā’ ā Ā· K( āP ( ) + 2

gez ) = f , (0.21)

which also contains derivatives of ļ¬rst order in space. If P ( ) = ln(Ī± ) holds

for a constant Ī± > 0, then āP ( ) simpliļ¬es to Ī±ā . Thus for horizontal

ļ¬‚ow we again encounter the heat conduction equation. For the relationship

P ( ) = Ī± suggested by the universal gas law, Ī± ā = 1 Ī±ā 2 remains

2

nonlinear. The choice of the variable u := 2 would result in u1/2 in the

time derivative as the only nonlinearity. Thus in the formulation in the

coeļ¬cient of ā disappears in the divergence of = 0. Correspondingly,

the coeļ¬cient S(u) = 1 Ļuā’1/2 of ā‚t u in the formulation in u becomes

2

unbounded for u = 0. In both versions the equations are degenerate, whose

treatment is beyond the scope of this book. A variant of this equation has

gained much attention as the porous medium equation (with convection) in

the ļ¬eld of analysis (see, for example, [42]).

Returning to the general framework, the following generalization of

Darcyā™s law can be justiļ¬ed experimentally for several liquid phases:

krĪ±

qĪ± = ā’ k (āpĪ± + Ī± gez ) .

ĀµĪ±

Here the relative permeability krĪ± of the phase Ī± depends upon the

saturations of the present phases and takes values in [0, 1].

10 0. Modelling Processes in Porous Media with Diļ¬erential Equations

At the interface of two liquid phases Ī±1 and Ī±2 we observe a diļ¬erence of

the pressures, the so-called capillary pressure, that turns out experimentally

to be a function of the saturations:

pc Ī±1 Ī±2 := pĪ±1 ā’ pĪ±2 = FĪ±1 Ī±2 (Sw , So , Sg ) . (0.22)

A general model for multiphase ļ¬‚ow, formulated for the moment in terms

of the variables pĪ± , SĪ± , is thus given by the equations

ā’āĀ·(

ā‚t (ĻSĪ± Ī±) Ī± Ī»Ī± k(āpĪ± + Ī± gez )) = fĪ± (0.23)

with the mobilities Ī»Ī± := krĪ± /ĀµĪ± , and the equations (0.22) and (0.10),

where one of the SĪ± ā™s can be eliminated. For two liquid phases w and g,

e.g., water and air, equations (0.22) and (0.10) for Ī± = w, g read pc =

pg ā’ pw = F (Sw ) and Sg = 1 ā’ Sw . Apparently, this is a time-dependent,

nonlinear model in the variables pw , pg , Sw , where one of the variables can

be eliminated. Assuming constant densities Ī± , further formulations based

on

ā Ā· q w + q g = fw / + fg / (0.24)

w g

can be given as consequences of (0.10). These equations consist of a sta-

tionary equation for a new quantity, the global pressure, based on (0.24),

and a time-dependent equation for one of the saturations (see Exercise 0.2).

In many situations it is justiļ¬ed to assume a gaseous phase with constant

pressure in the whole domain and to scale this pressure to pg = 0. Thus

for Ļ := pw = ā’pc we have

Ļā‚t S(Ļ) ā’ ā Ā· (Ī»(Ļ)k(āĻ + gez )) = fw / (0.25)

w

:= w , and S(Ļ) := F ā’1 (ā’Ļ) as a strictly

with constant pressure

monotone increasing nonlinearity as well as Ī».

With the convention to set the value of the air pressure to 0, the pressure

in the aqueous phase is in the unsaturated state, where the gaseous phase is

also present, and represented by negative values. The water pressure Ļ = 0

marks the transition from the unsaturated to the saturated zone. Thus

in the unsaturated zone, equation (0.25) represents a nonlinear variant

of the heat conduction equation for Ļ < 0, the Richards equation. As

most functional relationships have the property S (0) = 0, the equation

degenerates in the absence of a gaseous phase, namely to a stationary

equation in a way that is diļ¬erent from above.

Equation (0.25) with S(Ļ) := 1 and Ī»(Ļ) := Ī»(0) can be continued in a

consistent way with (0.14) and (0.15) also for Ļ ā„ 0, i.e., for the case of a

sole aqueous phase. The resulting equation is also called Richards equation

or a model of saturated-unsaturated ļ¬‚ow.

0.4. Reactive Solute Transport in Porous Media 11

0.4 Reactive Solute Transport in Porous Media

In this chapter we will discuss the transport of a single component in a

liquid phase and some selected reactions. We will always refer to water

as liquid phase explicitly. Although we treat inhomogeneous reactions in

terms of surface reactions with the solid phase, we want to ignore exchange

processes between the ļ¬‚uid phases. On the microscopic scale the mass con-

servation law for a single component Ī· is, in the notation of (0.11) by

omitting the phase index w,

ā‚t ĖĪ· + ā Ā· (ĖĪ· q ) + ā Ā· J Ī· = QĪ· ,

Ė

Ė

where

J Ī· := ĖĪ· (Ė Ī· ā’ q ) [kg/m2 /s]

Ė

v (0.26)

Ė

represents the diļ¬usive mass ļ¬‚ux of the component Ī· and QĪ· [kg/m3 /s] is

its volumetric production rate. For a description of reactions via the mass

action law it is appropriate to choose the mole as the unit of mass. The

diļ¬usive mass ļ¬‚ux requires a phenomenological description. The assump-

tion that solely binary molecular diļ¬usion, described by Fickā™s law, acts

between the component Ī· and the solvent, means that

J Ī· = ā’ ĖDĪ· ā (ĖĪ· / Ė) (0.27)

with a molecular diļ¬usivity DĪ· > 0 [m2 /s]. The averaging procedure applied

on (0.26), (0.27) leads to

ā‚t (Ī˜cĪ· ) + ā Ā· (qcĪ· ) + ā Ā· J (1) + ā Ā· J (2) = Q(1) + Q(2)

Ī· Ī·

for the solute concentration of the component Ī·, cĪ· [kg/m3 ], as intrinsic

phase average of ĖĪ· . Here, we have J (1) as the average of J Ī· and J (2) ,

the mass ļ¬‚ux due to mechanical dispersion, a newly emerging term at the

(1) Ė

macroscopic scale. Analogously, QĪ· is the intrinsic phase average of QĪ· ,

(2)

and QĪ· is a newly emerging term describing the exchange between the

liquid and solid phases.

The volumetric water content is given by Ī˜ := ĻSw with the water

saturation Sw . Experimentally, the following phenomenological descriptions

are suggested:

J (1) = ā’Ī˜Ļ„ DĪ· ācĪ·

with a tortuosity factor Ļ„ ā (0, 1],

J (2) = ā’Ī˜Dmech ācĪ· , (0.28)

and a symmetric positive deļ¬nite matrix of mechanical dispersion D mech ,

which depends on q/Ī˜. Consequently, the resulting diļ¬erential equation

reads

ā‚t (Ī˜cĪ· ) + ā Ā· (qcĪ· ā’ Ī˜DācĪ· ) = QĪ· (0.29)

12 0. Modelling Processes in Porous Media with Diļ¬erential Equations

(1) (2)

with D := Ļ„ D Ī· + D mech, QĪ· := QĪ· + QĪ· .

Because the mass ļ¬‚ux consists of qcĪ· , a part due to forced convection, and

of J (1) + J (2) , a part that corresponds to a generalized Fickā™s law, an equa-

tion like (0.29) is called a convection-diļ¬usion equation. Accordingly, for

the part with ļ¬rst spatial derivatives like ā Ā· (qcĪ· ) the term convective part

is used, and for the part with second spatial derivatives like ā’ā Ā· (Ī˜DācĪ· )

the term diļ¬usive part is used. If the ļ¬rst term determines the character of

the solution, the equation is called convection-dominated. The occurrence

of such a situation is measured by the quantity Pe, the global PĀ“clet num-

e

ber, that has the form Pe = q L/ Ī˜D [ - ]. Here L is a characteristic

length of the domain ā„¦. The extreme case of purely convective transport

results in a conservation equation of ļ¬rst order. Since the common mod-

els for the dispersion matrix lead to a bound for Pe, the reduction to the

purely convective transport is not reasonable. However, we have to take

convection-dominated problems into consideration.

Likewise, we speak of diļ¬usive parts in (0.17) and (0.20) and of (nonlin-

ear) diļ¬usive and convective parts in (0.21) and (0.25). Also, the multiphase

transport equation can be formulated as a nonlinear convection-diļ¬usion

equation by use of (0.24) (see Exercise 0.2), where convection often dom-

inates. If the production rate QĪ· is independent of cĪ· , equation (0.29) is

linear.

In general, in case of a surface reaction of the component Ī·, the kinetics of

the reaction have to be described . If this component is not in competition

with the other components, one speaks of adsorption. The kinetic equation

thus takes the general form

ā‚t sĪ· (x, t) = kĪ· fĪ· (x, cĪ· (x, t), sĪ· (x, t)) (0.30)

with a rate parameter kĪ· for the sorbed concentration sĪ· [kg/kg], which is

given in reference to the mass of the solid matrix. Here, the components

in sorbed form are considered spatially immobile. The conservation of the

total mass of the component undergoing sorption gives

Q(2) = ā’ b ā‚t sĪ· (0.31)

Ī·

with the bulk density b = s (1ā’Ļ), where s denotes the density of the solid

phase. With (0.30), (0.31) we have a system consisting of an instationary

partial and an ordinary diļ¬erential equation (with x ā ā„¦ as parameter). A

widespread model by Langmuir reads

fĪ· = ka cĪ· (sĪ· ā’ sĪ· ) ā’ kd sĪ·

with constants ka , kd that depend upon the temperature (among other

factors), and a saturation concentration sĪ· (cf. for example [24]). If we

assume fĪ· = fĪ· (x, cĪ· ) for simplicity, we get a scalar nonlinear equation in

cĪ· ,

ā‚t (Ī˜cĪ· ) + ā Ā· (qcĪ· ā’ Ī˜DācĪ· ) + = Q(1) ,

b kĪ· fĪ· (Ā·, cĪ· ) (0.32)

Ī·

0.4. Reactive Solute Transport in Porous Media 13

and sĪ· is decoupled and extracted from (0.30). If the time scales of transport

and reaction diļ¬er greatly, and the limit case kĪ· ā’ ā is reasonable, then

(0.30) is replaced by

fĪ· (x, cĪ· (x, t), sĪ· (x, t)) = 0 .

If this equation is solvable for sĪ· , i.e.,

sĪ· (x, t) = Ļ•Ī· (x, cĪ· (x, t)) ,

the following scalar equation for cĪ· with a nonlinearity in the time

derivative emerges:

+ ā Ā· (qcĪ· ā’ Ī˜DācĪ· ) = Q(1) .

ā‚t (Ī˜cĪ· + b Ļ•Ī· (Ā·, cĪ· )) Ī·

If the component Ī· is in competition with other components in the sur-

face reaction, as, e.g., in ion exchange, then fĪ· has to be replaced by a

nonlinearity that depends on the concentrations of all involved components

c1 , . . . , cN , s1 , . . . , sN . Thus we obtain a coupled system in these variables.

Finally, if we encounter homogeneous reactions that take place solely in the

(1)

ļ¬‚uid phase, an analogous statement is true for the source term QĪ· .

Exercises

0.1 Give a geometric interpretation for the matrix condition of k in (0.16)

and D mech in (0.28).

Ī± ā {w, g})

0.2 Consider the two-phase ļ¬‚ow (with constant Ī±,

ā‚t (ĻSĪ± ) + ā Ā· q Ī± = fĪ± ,

ā’Ī»Ī± k (āpĪ± +

qĪ± = Ī± gez ) ,

Sw + Sg = 1,

pg ā’ pw = pc

with coeļ¬cient functions

Ī± ā {w, g}.

pc = pc (Sw ) , Ī»Ī± = Ī»Ī± (Sw ) ,

Starting from equation (0.23), perform a transformation to the new

variables

q = qw + qg , ātotal ļ¬‚ow,ā

Ī»g ā’ Ī»w dpc

S

1 1

p= (pw + pg ) + dĪ¾ , āglobal pressure,ā

2 2 Ī»g + Ī»w dĪ¾

Sc

and the water saturation Sw . Derive a representation of the phase ļ¬‚ows in

the new variables.

14 0. Modelling Processes in Porous Media with Diļ¬erential Equations

0.3 A frequently employed model for mechanical dispersion is

Dmech = Ī»L |v|2 Pv + Ī»T |v|2 (I ā’ Pv )

with parameters Ī»L > Ī»T , where v = q/Ī˜ and Pv = vv T /|v|2 . Here2

Ī»L and Ī»T are the longitudinal and transversal dispersion lengths. Give a

geometrical interpretation.

0.5 Boundary and Initial Value Problems

The diļ¬erential equations that we derived in Sections 0.3 and 0.4 have the

common form

ā‚t S(u) + ā Ā· (C(u) ā’ K(āu)) = Q(u) (0.33)

with a source term S, a convective part C, a diļ¬usive part K, i.e., a total

ļ¬‚ux C ā’ K and a source term Q, which depend linearly or nonlinearly

on the unknown u. For simpliļ¬cation, we assume u to be a scalar. The

nonlinearities S, C, K, and Q may also depend on x and t, which shall be

suppressed in the notation in the following. Such an equation is said to be

in divergence form or in conservative form; a more general formulation is

obtained by diļ¬erentiating ā Ā· C(u) = ā‚u C(u) Ā· āu + (ā Ā· C)(u) or by

ā‚

introducing a generalized āsource termā Q = Q(u, āu). Up to now we have

considered diļ¬erential equations pointwise in x ā ā„¦ (and t ā (0, T )) under

the assumption that all occurring functions are well-deļ¬ned. Due to the

applicability of the integral theorem of Gauss on ā„¦ ā‚ ā„¦ (cf. (3.10)), the

Ė

integral form of the conservation equation follows straightforwardly from

the above:

(C(u) ā’ K(āu)) Ā· Ī½ dĻ = Q(u, āu) dx

ā‚t S(u) dx + (0.34)

Ė

Ė Ė

ā‚ā„¦

ā„¦ ā„¦

with the outer unit normal Ī½ (see Theorem 3.8) for a ļ¬xed time t or also

in t integrated over (0, T ). Indeed, this equation (on the microscopic scale)

is the primary description of the conservation of an extensive quantity:

Ė

Changes in time through storage and sources in ā„¦ are compensated by the

normal ļ¬‚ux over ā‚ ā„¦. Moreover, for ā‚t S, ā Ā· (C ā’ K), and Q continuous

Ė

Ė

on the closure of ā„¦, (0.33) follows from (0.34). If, on the other hand, F is

Ė

a hyperplane in ā„¦ where the material properties may rapidly change, the

jump condition

[(C(u) ā’ K(āu)) Ā· Ī½] = 0 (0.35)

for a ļ¬xed unit normal Ī½ on F follows from (0.34), where [ Ā· ] denotes the

diļ¬erence of the one-sided limits (see Exercise 0.4).

Since the diļ¬erential equation describes conservation only in general,

it has to be supplemented by initial and boundary conditions in order to

0.5. Boundary and Initial Value Problems 15

specify a particular situation where a unique solution is expected. Boundary

conditions are speciļ¬cations on ā‚ā„¦, where Ī½ denotes the outer unit normal

ā¢ of the normal component of the ļ¬‚ux (inwards):

ā’ (C(u) ā’ K(āu)) Ā· Ī½ = g1 on Ī“1 (0.36)

(ļ¬‚ux boundary condition),

ā¢ of a linear combination of the normal ļ¬‚ux and the unknown itself:

ā’ (C(u) ā’ K(āu)) Ā· Ī½ + Ī±u = g2 on Ī“2 (0.37)

(mixed boundary condition),

ā¢ of the unknown itself:

u = g3 on Ī“3 (0.38)

(Dirichlet boundary condition).

Here Ī“1 , Ī“2 , Ī“3 form a disjoint decomposition of ā‚ā„¦:

ā‚ā„¦ = Ī“1 āŖ Ī“2 āŖ Ī“3 , (0.39)

where Ī“3 is supposed to be a closed subset of ā‚ā„¦. The inhomogeneities

gi and the factor Ī± in general depend on x ā ā„¦, and for nonstationary

problems (where S(u) = 0 holds) on t ā (0, T ). The boundary conditions

are linear if the gi do not depend (nonlinearly) on u (see below). If the gi

are zero, we speak of homogeneous, otherwise of inhomogeneous, boundary

conditions.

Thus the pointwise formulation of a nonstationary equation (where S

does not vanish) requires the validity of the equation in the space-time

cylinder

QT := ā„¦ Ć— (0, T )

and the boundary conditions on the lateral surface of the space-time

cylinder

ST := ā‚ā„¦ Ć— (0, T ) .

Diļ¬erent types of boundary conditions are possible with decompositions

of the type (0.39). Additionally, an initial condition on the bottom of the

space-time cylinder is necessary:

for x ā ā„¦ .

S(u(x, 0)) = S0 (x) (0.40)

These are so-called initial-boundary value problems; for stationary prob-

lems we speak of boundary value problems. As shown in (0.34) and (0.35)

ļ¬‚ux boundary conditions have a natural relationship with the diļ¬erential

equation (0.33). For a linear diļ¬usive part K(āu) = Kāu alternatively

we may require

ā‚Ī½K u := Kāu Ā· Ī½ = g1 on Ī“1 , (0.41)

16 0. Modelling Processes in Porous Media with Diļ¬erential Equations

and an analogous mixed boundary condition. This boundary condition is

the so-called Neumann boundary condition. Since K is symmetric, ā‚Ī½K u =

āu Ā· KĪ½ holds; i.e., ā‚Ī½K u is the derivative in direction of the conormal KĪ½.

For the special case K = I the normal derivative is given.

In contrast to ordinary diļ¬erential equations, there is hardly any general

theory of partial diļ¬erential equations. In fact, we have to distinguish dif-

ferent types of diļ¬erential equations according to the various described

physical phenomena. These determine, as discussed, diļ¬erent (initial-)

boundary value speciļ¬cations to render the problem well-posed. Well-

posedness means that the problem possesses a unique solution (with certain

properties yet to be deļ¬ned) that depends continuously (in appropriate

norms) on the data of the problem, in particular on the (initial and)

boundary values. There exist also ill-posed boundary value problems for

partial diļ¬erential equations, which correspond to physical and technical

applications. They require special techniques and shall not be treated here.

The classiļ¬cation into diļ¬erent types is simple if the problem is lin-

ear and the diļ¬erential equation is of second order as in (0.33). By order

we mean the highest order of the derivative with respect to the variables

(x1 , . . . , xd , t) that appears, where the time derivative is considered to be

like a spatial derivative. Almost all diļ¬erential equations treated in this

book will be of second order, although important models in elasticity the-

ory are of fourth order or certain transport phenomena are modelled by

systems of ļ¬rst order.

The diļ¬erential equation (0.33) is generally nonlinear due to the nonlin-

ear relationships S, C, K, and Q. Such an equation is called quasilinear if

all derivatives of the highest order are linear, i.e., we have

K(āu) = Kāu (0.42)

with a matrix K, which may also depend (nonlinearly) on x, t, and u.

Furthermore, (0.33) is called semilinear if nonlinearities are present only

in u, but not in the derivatives, i.e., if in addition to (0.42) with K being

independent of u, we have

S(u) = Su , C(u) = uc (0.43)

with scalar and vectorial functions S and c, respectively, which may depend

on x and t. Such variable factors standing before u or diļ¬erential terms are

called coeļ¬cients in general.

Finally, the diļ¬erential equation is linear if we have, in addition to the

above requirements,

Q(u) = ā’ru + f

with functions r and f of x and t.

In the case f = 0 the linear diļ¬erential equation is termed homoge-

neous, otherwise inhomogeneous. A linear diļ¬erential equation obeys the

superposition principle: Suppose u1 and u2 are solutions of (0.33) with the

0.5. Boundary and Initial Value Problems 17

source terms f1 and f2 and otherwise identical coeļ¬cient functions. Then

u1 + Ī³u2 is a solution of the same diļ¬erential equation with the source

term f1 + Ī³f2 for arbitrary Ī³ ā R. The same holds for linear boundary

conditions. The term solution of an (initial-) boundary value problem is

used here in a classical sense, yet to be speciļ¬ed, where all the quantities

occurring should satisfy pointwise certain regularity conditions (see Deļ¬ni-

tion 1.1 for the Poisson equation). However, for variational solutions (see

Deļ¬nition 2.2), which are appropriate in the framework of ļ¬nite element

methods, the above statements are also valid.

Linear diļ¬erential equations of second order in two variables (x, y) (in-

cluding possibly the time variable) can be classiļ¬ed in diļ¬erent types as

follows:

To the homogeneous diļ¬erential equation

ā‚2 ā‚2 ā‚2

Lu = a(x, y) 2 u + b(x, y) u + c(x, y) 2 u

ā‚x ā‚xā‚y ā‚y (0.44)

ā‚ ā‚

+ d(x, y) u + e(x, y) u + f (x, y)u = 0

ā‚x ā‚y

the following quadratic form is assigned:

(Ī¾, Ī·) ā’ a(x, y)Ī¾ 2 + b(x, y)Ī¾Ī· + c(x, y)Ī· 2 . (0.45)

According to its eigenvalues, i.e., the eigenvalues of the matrix

1

a(x, y) 2 b(x, y)

, (0.46)

1

2 b(x, y) c(x, y)

we classify the types. In analogy with the classiļ¬cation of conic sections,

which are described by (0.45) (for ļ¬xed (x, y)), the diļ¬erential equation

(0.44) is called at the point (x, y)

ā¢ elliptic if the eigenvalues of (0.46) are not 0 and have the same sign,

ā¢ hyperbolic if one eigenvalue is positive and the other is negative,

ā¢ parabolic if exactly one eigenvalue is equal to 0.

For the corresponding generalization of the terms for d + 1 variables and

arbitrary order, the stationary boundary value problems we treat in this

book will be elliptic, of second order, and ā” except in Chapter 8 ā” also

linear; the nonstationary initial-boundary value problems will be parabolic.

Systems of hyperbolic diļ¬erential equations of ļ¬rst order require partic-

ular approaches, which are beyond the scope of this book. Nevertheless,

we dedicate Chapter 9 to convection-dominated problems, i.e., elliptic or

parabolic problems close to the hyperbolic limit case.

The diļ¬erent discretization strategies are based on various formulations

of the (initial-) boundary value problems: The ļ¬nite diļ¬erence method,

which is presented in Section 1, and further outlined for nonstationary prob-

lems in Chapter 7, has the pointwise formulation of (0.33), (0.36)ā“(0.38)

18 0. Modelling Processes in Porous Media with Diļ¬erential Equations

(and (0.40)) as a starting point. The ļ¬nite element method, , which lies in

the focus of our book (Chapters 2, 3, and 7), is based on an integral formu-

lation of (0.33) (which we still have to depict) that incorporates (0.36) and

(0.37). The conditions (0.38) and (0.40) have to be enforced additionally.

Finally, the ļ¬nite volume method (Chapters 6 and 7) will be derived from

the integral formulation (0.34), where also initial and boundary conditions

come along as in the ļ¬nite element approach.

Exercises

0.4 Derive (formally) (0.35) from (0.34).

0.5 Derive the orders of the given diļ¬erential operators and diļ¬er-

ential equations, and decide in every case whether the operator is

linear or nonlinear, and whether the linear equation is homogeneous or

inhomogeneous:

(a) Lu := uxx + xuy ,

(b) Lu := ux + uuy ,

ā

(c) Lu := 1 + x2 (cos y)ux + uyxy ā’ arctan x u = ln(x2 + y 2 ) ,

y

ā

(d) Lu := ut + uxxxx + 1 + u = 0 ,

(e) utt ā’ uxx + x2 = 0 .

0.6 (a) Determine the type of the given diļ¬erential operator:

(i) Lu := uxx ā’ uxy + 2uy + uyy ā’ 3uyx + 4u ,

(ii) Lu = 9uxx + 6uxy + uyy + ux .

(b) Determine the parts of the plane where the diļ¬erential operator Lu :=

yuxx ā’ 2uxy + xuyy is elliptic, hyperbolic, or parabolic.

(c) (i) Determine the type of Lu := 3uy + uxy .

(ii) Compute the general solution of Lu = 0.

0.7 Consider the equation Lu = f with a linear diļ¬erential operator of

second order, deļ¬ned for functions in d variables (d ā N) in x ā ā„¦ ā‚ Rd .

The transformation Ī¦ : ā„¦ ā’ ā„¦ ā‚ Rd has a continuously diļ¬erentiable,

nonsingular Jacobi matrix DĪ¦ := ā‚Ī¦ .

ā‚x

Show that the partial diļ¬erential equation does not change its type if it

is written in the new coordinates Ī¾ = Ī¦(x).

1

For the Beginning:

The Finite Diļ¬erence Method for the

Poisson Equation

1.1 The Dirichlet Problem for the Poisson

Equation

In this section we want to introduce the ļ¬nite diļ¬erence method using

the Poisson equation on a rectangle as an example. By means of this ex-

ample and generalizations of the problem, advantages and limitations of

the approach will be elucidated. Also, in the following section the Poisson

equation will be the main topic, but then on an arbitrary domain. For the

spatial basic set of the diļ¬erential equation ā„¦ ā‚ Rd we assume as minimal

requirement that ā„¦ is a domain, where a domain is a nonempty, open, and

connected set. The boundary of this domain will be denoted by ā‚ā„¦, the

closure ā„¦ āŖ ā‚ā„¦ by ā„¦ (see Appendix A.2). The Dirichlet problem for the

Poisson equation is then deļ¬ned as follows: Given functions g : ā‚ā„¦ ā’ R

and f : ā„¦ ā’ R, we are looking for a function u : ā„¦ ā’ R such that

d

ā‚2

ā’ u = f in ā„¦ , (1.1)

ā‚x2

i

i=1

u = g on ā‚ā„¦ . (1.2)

This diļ¬erential equation model has already appeared in (0.19) and

(0.38) and beyond this application has an importance in a wide spectrum of

disciplines. The unknown function u can be interpreted as an electromag-

netic potential, a displacement of an elastic membrane, or a temperature.

Similar to the multi-index notation to be introduced in (2.16) (but with

20 1. Finite Diļ¬erence Method for the Poisson Equation

indices at the top) from now on for partial derivatives we use the following

notation.

Notation: For u : ā„¦ ā‚ Rd ā’ R we set

ā‚

ā‚i u := u for i = 1, . . . , d ,

ā‚xi

ā‚2

for i, j = 1, Ā· Ā· Ā· , d ,

ā‚ij u := u

ā‚xi ā‚xj

āu := (ā‚11 + . . . + ā‚dd ) u .

The expression āu is called the Laplace operator. By means of this, (1.1)

can be written in abbreviated form as

ā’āu = f in ā„¦ . (1.3)

We could also deļ¬ne the Laplace operator by

āu = ā Ā· (āu) ,

T

where āu = (ā‚1 u, . . . , ā‚d u) denotes the gradient of a function u, and

ā Ā· v = ā‚1 v1 + Ā· Ā· Ā· + ā‚d vd the divergence of a vector ļ¬eld v. Therefore,

an alternative notation exists, which will not be used in the following:

āu = ā2 u. The incorporation of the minus sign in the left-hand side of

(1.3), which looks strange at ļ¬rst glance, is related to the monotonicity and

deļ¬niteness properties of ā’ā (see Sections 1.4 and 2.1, respectively).

The notion of a solution for (1.1), (1.2) still has to speciļ¬ed more pre-

cisely. Considering the equations in a pointwise sense, which will be pursued

in this chapter, the functions in (1.1), (1.2) have to exist, and the equations

have to be satisļ¬ed pointwise. Since (1.1) is an equation on an open set ā„¦,

there are no implications for the behaviour of u up to the boundary ā‚ā„¦. To

have a real requirement due to the boundary condition, u has to be at least

continuous up to the boundary, that is, on ā„¦. These requirements can be

formulated in a compact way by means of corresponding function spaces.

The function spaces are introduced more precisely in Appendix A.5. Some

examples are

u : ā„¦ ā’ R u continuous in ā„¦ ,

C(ā„¦) :=

u : ā„¦ ā’ R u ā C(ā„¦) , ā‚i u exists in ā„¦ ,

C 1 (ā„¦) :=

ā‚i u ā C(ā„¦) for all i = 1, . . . , d .

The spaces C k (ā„¦) for k ā N, C(ā„¦), and C k (ā„¦), as well as C(ā‚ā„¦), are

deļ¬ned analogously. In general, the requirements related to the (contin-

uous) existence of derivatives are called, a little bit vaguely, smoothness

requirements.

In the following, in view of the ļ¬nite diļ¬erence method, f and g will also

be assumed continuous in ā„¦ and ā‚ā„¦, respectively.

1.2. The Finite Diļ¬erence Method 21

Deļ¬nition 1.1 Assume f ā C(ā„¦) and g ā C(ā‚ā„¦). A function u is called

a (classical) solution of (1.1), (1.2) if u ā C 2 (ā„¦) ā© C(ā„¦), (1.1) holds for all

x ā ā„¦, and (1.2) holds for all x ā ā‚ā„¦.

1.2 The Finite Diļ¬erence Method

The ļ¬nite diļ¬erence method is based on the following approach: We are

looking for an approximation to the solution of a boundary value problem

at a ļ¬nite number of points in ā„¦ (the grid points). For this reason we

substitute the derivatives in (1.1) by diļ¬erence quotients, which involve

only function values at grid points in ā„¦ and require (1.2) only at grid

points. By this we obtain algebraic equations for the approximating values

at grid points. In general, such a procedure is called the discretization of the

boundary value problem. Since the boundary value problem is linear, the

system of equations for the approximate values is also linear. In general, for

other (diļ¬erential equation) problems and other discretization approaches

we also speak of the discrete problem as an approximation of the continuous

problem. The aim of further investigations will be to estimate the resulting

error and thus to judge the quality of the approximative solution.

Generation of Grid Points

In the following, for the beginning, we will restrict our attention to problems

in two space dimensions (d = 2). For simpliļ¬cation we consider the case

of a constant step size (or mesh width) h > 0 in both space directions.

The quantity h here is the discretization parameter, which in particular

determines the dimension of the discrete problem.

ā—¦ā—¦ā—¦ ā—¦ ā—¦ ā—¦ ā—¦ ā—¦ ā—¦ ā¢ : ā„¦h

ā—¦ ā¢ā¢ ā¢ ā¢ ā¢ ā¢ ā¢ ā—¦

3

ā—¦ : ā‚ā„¦h

ā—¦ ā¢ā¢ ā¢ ā¢ ā¢ ā¢ ā¢ ā—¦

l=8 2

ā—¦ ā¢ā¢ ā¢ ā¢ ā¢ ā¢ ā¢ ā—¦

2 : far from boundary

m=5

ā—¦ ā¢ā¢ ā¢ ā¢ ā¢ ā¢ ā¢ ā—¦

ā—¦ā—¦ā—¦ ā—¦ ā—¦ ā—¦ ā—¦ ā—¦ ā—¦ 3 : close to boundary

Figure 1.1. Grid points in a square domain.

For the time being, let ā„¦ be a rectangle, which represents the simplest

case for the ļ¬nite diļ¬erence method (see Figure 1.1). By translation of the

coordinate system the situation can be reduced to ā„¦ = (0, a) Ć— (0, b) with

a, b > 0 . We assume that the lengths a, b, and h are such that

for certain l, m ā N.

a = lh, b = mh (1.4)

22 1. Finite Diļ¬erence Method for the Poisson Equation

We deļ¬ne

(ih, jh) i = 1, . . . , l ā’ 1 , j = 1, . . . , m ā’ 1

ā„¦h :=

(1.5)

(x, y) ā ā„¦ x = ih , y = jh with i, j ā Z

=

as a set of grid points in ā„¦ in which an approximation of the diļ¬erential

equation has to be satisļ¬ed. In the same way,

(ih, jh) i ā {0, l} , j ā {0, . . . , m} or i ā {0, . . . , l} , j ā {0, m}

ā‚ā„¦h :=

(x, y) ā ā‚ā„¦ x = ih , y = jh with i, j ā Z

=

deļ¬nes the grid points on ā‚ā„¦ in which an approximation of the boundary

condition has to be satisļ¬ed. The union of grid points will be denoted by

ā„¦h := ā„¦h āŖ ā‚ā„¦h .

Setup of the System of Equations

Lemma 1.2 Let ā„¦ := (x ā’ h, x + h) for x ā R, h > 0. Then there exists

a quantity R, depending on u and h, the absolute value of which can be

bounded independently of h and such that

(1) for u ā C 2 (ā„¦):

u(x + h) ā’ u(x) 1

|R| ā¤

u (x) = + hR and u ,

ā

h 2

(2) for u ā C 2 (ā„¦):

u(x) ā’ u(x ā’ h) 1

|R| ā¤

u (x) = + hR and u ,

ā

h 2

(3) for u ā C 3 (ā„¦):

u(x + h) ā’ u(x ā’ h) 1

|R| ā¤

+ h2 R

u (x) = and u ,

ā

2h 6

(4) for u ā C 4 (ā„¦):

u(x + h) ā’ 2u(x) + u(x ā’ h) 1 (4)

+ h2 R and |R| ā¤

u (x) = u .

ā

2

h 12

Here the maximum norm Ā· ā (see Appendix A.5) has to be taken over

the interval of the involved points (x, x + h), (x ā’ h, x), or (x ā’ h, x + h).

Proof: The proof follows immediately by Taylor expansion. As an example

we consider statement 3: From

h2 h3

u(x Ā± h) = u(x) Ā± hu (x) + u (x) Ā± u (x Ā± Ī¾Ā± ) for certain Ī¾Ā± ā (0, h)

2 6

2

the assertion follows by linear combination.

1.2. Derivation and Properties 23

Notation: The quotient in statement 1 is called the forward diļ¬erence

quotient, and it is denoted by ā‚ + u(x). The quotient in statement 2 is

called the backward diļ¬erence quotient (ā‚ ā’ u(x)), and the one in statement

3 the symmetric diļ¬erence quotient (ā‚ 0 u(x)). The quotient appearing in

statement 4 can be written as ā‚ ā’ ā‚ + u(x) by means of the above notation.

In order to use statement 4 in every space direction for the approximation

of ā‚11 u and ā‚22 u in a grid point (ih, jh), in addition to the conditions of

Deļ¬nition 1.1, the further smoothness properties ā‚ (3,0) u, ā‚ (4,0) u ā C(ā„¦)

and analogously for the second coordinate are necessary. Here we use, e.g.,

the notation ā‚ (3,0) u := ā‚ 3 u/ā‚x3 (see (2.16)).

1

Using these approximations for the boundary value problem (1.1), (1.2),

at each grid point (ih, jh) ā ā„¦h we get

u ((i + 1)h, jh) ā’ 2u(ih, jh) + u ((i ā’ 1)h, jh)

ā’

h2

u (ih, (j + 1)h) ā’ 2u(ih, jh) + u (ih, (j ā’ 1)h)

+ = (1.6)

h2

f (ih, jh) + R(ih, jh)h2 .

=

Here R is as described in statement 4 of Lemma 1.2, a function depending

on the solution u and on the step size h, but the absolute value of which can

be bounded independently of h. In cases where we have less smoothness of

the solution u, we can nevertheless formulate the approximation (1.6) for

ā’āu, but the size of the error in the equation is unclear at the moment.

For the grid points (ih, jh) ā ā‚ā„¦h no approximation of the boundary

condition is necessary:

u(ih, jh) = g(ih, jh) .

If we neglect the term Rh2 in (1.6), we get a system of linear equations

for the approximating values uij for u(x, y) at points (x, y) = (ih, jh) ā ā„¦h .

They have the form

1

ā’ ui,jā’1 ā’ uiā’1,j + 4uij ā’ ui+1,j ā’ ui,j+1 = fij (1.7)

h2

for i = 1, . . . , l ā’ 1 , j = 1, . . . , m ā’ 1 ,

if i ā {0, l}, j = 0, . . . , m or j ā {0, m}, i = 0, . . . , l .

uij = gij (1.8)

Here we used the abbreviations

fij := f (ih, jh), gij := g(ih, jh) . (1.9)

Therefore, for each unknown grid value uij we get an equation. The grid

points (ih, jh) and the approximating values uij located at these have a

natural two-dimensional indexing.

In equation (1.7) for a grid point (i, j) only the neighbours at the four

cardinal points of the compass appear, as it is displayed in Figure 1.2. This

24 1. Finite Diļ¬erence Method for the Poisson Equation

interconnection is also called the ļ¬ve-point stencil of the diļ¬erence method

and the method the ļ¬ve-point stencil discretization.

y

T (i,j+1)

ā¢

(iā’1,j) (i,j) (i+1,j)

ā¢ ā¢ ā¢

(i,jā’1)

ā¢

Ex

Figure 1.2. Five-point stencil.

At the interior grid points (x, y) = (ih, jh) ā ā„¦h , two cases can be

distinguished:

(1) (i, j) has a position such that its all neighbouring grid points lie in

ā„¦h (far from the boundary).

(2) (i, j) has a position such that at least one neighbouring grid point

(r, s) lies on ā‚ā„¦h (close to the boundary). Then in equation (1.7) the

value urs is known due to (1.8) (urs = grs ), and (1.7) can be modiļ¬ed

in the following way:

Remove the values urs with (rh, sh) ā ā‚ā„¦h in the equations for (i, j)

close to the boundary and add the value grs /h2 to the right-hand

side of (1.7). The set of equations that arises by this elimination of

boundary unknowns by means of Dirichlet boundary conditions we

call (1.7)ā— ; it is equivalent to (1.7), (1.8).

Instead of considering the values uij , i = 1, . . . , l ā’ 1, j = 1, . . . , m ā’ 1,

ńņš. 1 |