Poroelasticity
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee. This
dissertation does not include proprietary or classifled information.
Chadia Afiane Aji
Certiflcate of Approval:
Greg Harris
Associate Professor
Mathematics
Amnon J. Meir, Chair
Professor
Mathematics
Georg Hetzer
Professor
Mathematics
Frank Uhlig
Professor
Mathematics
Richard A. Zalik
Professor
Mathematics
Joe F. Pittman
Interim Dean
Graduate School
Poroelasticity
Chadia Afiane Aji
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulflllment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
August 4, 2007
Poroelasticity
Chadia Afiane Aji
Permission is granted to Auburn University to make copies of this dissertation at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Chadia Afiane Aji, daughter of Abbes Afiane Aji and Nouflssa Hassan, was born Oc-
tober 2, 1968, in Fes, Morocco. She is married and has three children. In 1991, she got an
associate degree in Chemical Engineering from Ecole Superieure de Technologie, E.S.T., in
Fes, Morocco. She entered Texas A&M University in 1997 and graduated with the degree
of Bachelor of Chemical engineering in 1999. Then she went back to her country where
she taught Mathematics for two years at Ifrane School in Morocco. She came to Auburn
in December 2001 when her husband joined the Computer Science faculty. She felt drawn
to learn more mathematics of an applied kind thus she became a graduate student in fall
2002. She flnished her Master of Applied Mathematics in December 2003. Since then she
has been in the Ph. D. program in Mathematics.
iv
Dissertation Abstract
Poroelasticity
Chadia Afiane Aji
Doctor of Philosophy, August 4, 2007
(M.A., Auburn University, 2003)
(B.S., Texas A&M University, 1999)
151 Typed Pages
Directed by Amnon J. Meir
Poroelasticity is the study of elastic deformation of porous materials saturated with a
 uid and the coupling between the  uid pressure and the solid deformation.
Considerable progress has been made in formulating analytical and numerical models of
subsurface  uid  ow, but only few models explain the interrelations between  uid- ow
pressure changes and seismicity.
In this work, we describe the quasi-static poroelasticity system of partial difierential
equations consisting of the equilibrium equation for momentum conservation and the difiu-
sion equation for Darcy  ow. We prove existence and uniqueness of weak solutions to the
equations of the quasi-static poroelasticity system and derive error estimates. We describe
a coupled numerical algorithm that accounts for the interrelations between the  uid pres-
sure changes and the deformation of the porous elastic material based on the flnite element
method using MATLAB.
v
Acknowledgments
I would like to thank my major professor Prof. A.J. Meir for his scientiflc guidance. He
has shared his excellent knowledge with me through his creativity. He was a very patient
and understanding advisor throughout this work. Without his support, this thesis would
not have been completed. I also wish to thank Dr. Uhlig for his support. He encouraged
me to engage in the Ph.D. program. His constant source of motivation kept me focused
throughout my graduate work at Auburn University. My sincere appreciation and gratitude
goes to Dr. Wolf and Dr. Lee from the Geology Department for their kindness, support,
and scientiflc guidance. I owe Dr. Wolf, Dr. Lee, and Dr. Meir the topic and the flnancial
support for this study through the U.S. Geological Survey (National Earthquake Hazards
Reduction Program to Wolf and Lee, under grant number 05HQ6R0081).
I also express my gratitude to my committee members Dr. Harris, Dr. Hetzer, and Dr.
Zalik.
This work is dedicated to my lovely father Abbes Afiane Aji (1927-1992). A special
thanks to my lovely mother Nouflssa Hassan, my brother Si Mohammed Afiane, my sisters,
and all my family for their ongoing encouragement and love.
Of all I am most thankful to my great kids Youssef, Rim, Yasmine, and I am very
grateful to my husband Saad Biaz for the unconditioned love and care. I would not be
where I am today without his invaluable support throughout my graduate studies.
vi
Style manual or journal used Journal of Approximation Theory (together with the style
known as \aums"). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package TEX (speciflcally LATEX)
together with the departmental style-flle aums.sty.
vii
Table of Contents
List of Figures x
1 Introduction 1
2 Poroelasticity model 5
2.1 The elasticity equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.2 Strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.3 Stress-strain relationship . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 The pore pressure equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 Boundary and initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . 23
3 Analysis 26
3.1 Existence and uniqueness of solutions (Showalter) . . . . . . . . . . . . . . . 27
3.2 Existence and uniqueness of weak solutions . . . . . . . . . . . . . . . . . . 30
3.3 Rothe?s method of lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3.1 Existence and uniqueness of weak solutions for
homogeneous initial and boundary conditions . . . . . . . . . . . . . 43
3.3.2 Energy norm estimate for
homogeneous initial and boundary conditions . . . . . . . . . . . . . 72
3.3.3 Existence and uniqueness of weak solutions for
nonhomogeneous initial conditions . . . . . . . . . . . . . . . . . . . 77
3.3.4 Energy norm estimate for nonhomogeneous initial conditions . . . . 86
3.3.5 Existence and uniqueness of weak solutions for
nonhomogeneous boundary conditions . . . . . . . . . . . . . . . . . 88
3.3.6 Energy norm estimate for
nonhomogeneous boundary condition . . . . . . . . . . . . . . . . . . 95
3.4 Error estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
3.4.1 Error estimates for the semi-discrete problem . . . . . . . . . . . . . 99
3.4.2 Error estimates for the fully discrete problem . . . . . . . . . . . . . 110
4 Numerical methods 117
4.1 2-D algorithm for the difiusion equation: 2dp ow . . . . . . . . . . . . . . . 117
4.2 3-D algorithm for the difiusion equation: 3dp ow . . . . . . . . . . . . . . . 122
4.3 3-D algorithm for the elasticity equation: 3dfem . . . . . . . . . . . . . . . . 123
4.4 Segregated algorithm: it3dupfem . . . . . . . . . . . . . . . . . . . . . . . . 125
4.5 Coupled algorithm: c3dupfem . . . . . . . . . . . . . . . . . . . . . . . . . . 126
viii
5 Conclusion and future work 129
Bibliography 132
Appendices 134
A Notation and Inequalities 135
A.1 Notation for derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
A.2 Spaces of continuous and difierentiable functions . . . . . . . . . . . . . . . 136
B Preliminaries 138
ix
List of Figures
2.1 Body force Fidv (dv = dx1dx2dx3, where dx1;dx2;dx3 are lengths of the
edges of the elements in x, y, and z-direction respectively) (see [6]) . . . . . 6
2.2 Stresses acting on surfaces and body force Fidx1dx2dx3 . . . . . . . . . . . 7
2.3 Stress components in x1-direction . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Normal strain in x1-direction . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Shear strain in x1, x2-direction . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.1 Comparing approximate solution pn and exact solution p from 2dp ow at the
last time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
4.2 Comparing approximate solution pn and exact solution p from 3dp ow at the
last time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
4.3 Comparing approximate solution un and exact solution u from 3dfem . . . 124
4.4 Comparing approximate solution pn and exact solution p from c3dupfem at
the last time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4.5 Comparing approximate solution un and exact solution u from c3dupfem at
the last time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
x
Chapter 1
Introduction
A porous medium, such as rock, sediment, or artiflcial porous material, is a material
with empty cavities called pores. These cavities may be fllled with liquids or gases. Peat
and clay are porous materials; about half of their volume consists of empty cavities. A
kitchen sponge is an artiflcial porous medium. Due to its nature, a porous medium is
usually elastic: when subjected to a force, a porous material may change its form but it
will often return to its original shape when the force is removed. The notion and study of
porous material in geology was flrst introduced in 1943 by Karl von Terzaghi (see [17]), the
father of soil mechanics. When a saturated porous medium is deformed, the volume of the
cavities changes, producing a change in the gas or liquid pressure. The relationship between
the deformation and the pressure changes is of interest in many geologic and engineering
applications.
Two mechanisms play a role in the interaction between  uid pressure changes and
deformation of the porous elastic material: (1) dilation of the medium results in a decrease
of pore pressure and, (2) compression of the material causes a rise of pore pressure, if
the compression is faster than the  uid  ow rate. For example, the water level in a well
changes when a train passes nearby. In 1892, F. H. King (see [18]) noticed that the water
level in a well near the train station at Whitewater, Wisconsin, went up when the train
approached the station and it went down when the train left the station. The change of
the water level depends on the weight of the train, that is, the water level increases more
for a heavy loaded freight train than for a passenger train. Another example of the coupled
1
pressure-deformation is a sponge whose pores are saturated with water. By compressing the
sponge, its form changes. The decrease of the volume of the pores creates an overpressure.
Therefore, the  uid is pressed out of the material and  ows away because of the increase of
the pore pressure. When releasing the sponge, i.e., reducing the pore pressure, the sponge
returns to its original form. This is explained by the elastic behavior of the material. This
coupled mechanism { namely the coupling of the stress in the solid with the pressure of the
 uid { plays an essential role in poroelasticity.
Poroelasticity is the study of elastic deformation of porous materials saturated with
a  uid and the coupling between the  uid pressure and the solid deformation. Anthony
Biot was the flrst to develop a model for such a relationship. His seminal paper [1] in 1941
describes a linear theory of poroelasticity which relates the evolution of  uid pressure p (a
scalar fleld) and the solid displacement u (a vector fleld). This system of equations is given
as follows
?utt ???u?(?+?)r(r?u)+firp = f; in ??(0;T); (1.1)
@
@t(c0p+fir?u)?r?krp = h; in ??(0;T); (1.2)
where ? is an open bounded non-empty set in R3 and T is a positive time.
This system consists of the momentum balance equations for the displacement of the
medium (1.1) and the mass balance equation for the pressure distribution (1.2). The co-
e?cient ? represents the local density of the porous medium. The constants, ?, called the
Lame constant, and ?, the shear modulus, are a measure of the strength of the material,
and are determined from the elasticity of the medium. The constant fi > 0 is the Biot-
Willis constant and accounts for the mechanical coupling of the porous media and the  uid
2
pressure. The coe?cient c0 > 0, called speciflc storage, is the amount of  uid which can
be forced into or out the medium by a unit pressure increment under constant volume.
The parameter k involves the permeability of the medium and the viscosity of the  uid
in Darcy?s law. The functions f and h are suitably given functions. Note that utt is the
second order partial derivative of u with respect to time, ? denotes the Laplace operator
in R3, r is the gradient, and r? is the divergence. These mathematical terms are deflned
in Appendix A.
Biot?s poroelasticity model is very general as it is independent of the application domain.
This model was extended by Coussy (see [4]) to take into account the heat-convection
phenomena. Other researchers have reflned the model for speciflc engineering flelds (see
[11], [12], [13], and [20]) such as geomechanics or petroleum engineering. The Biot?s system
model is complex and, in general, does not have closed form solutions. In 1993, Gomberg
and Ellis (see [7]) provided an algorithm dubbed 3D-DEF (a three-dimensional boundary
element program) that approximates Biot?s system for the displacement from which the
strain ? and the stress  can be calculated. In order to calculate pore pressure changes,
Lee and Wolf proposed in 1998 an algorithm dubbed 3P-Flow (see [9]) that uses the above
calculated strain ? and stress  . The algorithm 3D-DEF approximates solutions of the
quasi-static case of the elasticity equation (1.1) for the vector displacement u. Using these
results, 3P- ow can approximate the pressure in the difiusion equation (1.2). Thus the two
algorithms together do not treat the fully coupled system of the two partial difierential
equations. Furthermore, at the time there was no guarantee that a solution for the system
exists.
3
In 2000, using abstract theory (non constructive), Showalter (see [14]) showed existence
and uniqueness of strong solutions and weak solutions to Biot?s system in the quasi-static
case. A summary of his results is described in Section 3.1.
In this work, using a constructive approach (based on Babuska-Brezzi theory and
Rothe?s method of lines), we proved existence and uniqueness of weak solutions to the
equations of quasi-static poroelasticity (1.1)-(1.2).
This approach (a constructive approach) suggests numerical approximation methods and
allows derivation of error estimates.
Developing solvers for the coupled system (1.1)-(1.2) is an area of current research. To our
knowledge, there is no rigorous 3-dimensional error analysis for this coupled system.
The main contribution of this work is the construction of two algorithms for approx-
imating solutions of the quasi-static poroelasticity system of partial difierential equations:
a segregated algorithm where the solution is approximated by an iterative method and
a coupled algorithm where the system is concurrently approximated for both the vector
displacement u and the scalar pore pressure p.
This work is organized as follows. In the next chapter, we describe the mathematical model.
That is, we describe the quasi-static poroelasticity system of partial difierential equations.
Existence and uniqueness of weak solutions are proved and error estimates are derived in
Chapter 3. Chapter 4 describes numerical experiments. Finally, conclusion and proposed
future work are given in Chapter 5.
4
Chapter 2
Poroelasticity model
A quasi-static poroelastic problem is described by the following basic variables: stress
( ); the normalized force, strain (?); the symmetric part of the deformation gradient, the
vector displacement (u), the scalar pore pressure (p), and the increment of  uid content
(?).
In this section, we formulate the equations describing the coupling of elastic deformation
and pore  uid pressure in a porous medium. We flrst consider poroelastic constitutive
response { i.e., the dependence of strain and  uid content on stress and pore pressure { and
the Darcy law for pore  uid transport. Then we formulate the governing fleld equations
using considerations of stress equilibrium and mass conservation.
2.1 The elasticity equation
The equilibrium equation for momentum conservation will be formulated based on the
force equilibrium equation and the linear constitutive equation, the dependence of strain
and  uid content on stress. Therefore, we flrst need to deflne stress and strain to derive the
elasticity equation.
2.1.1 Stress
Consider a volume, an inflnitesimal cube with faces pointing in the coordinate direc-
tions. There are two types of external forces acting on the material body:
1. The force acting on volume elements of the body, called body force.
5
F  dv3
F  dv
2
F  dv
1
dx
1
dx
2
dx
3
x 3
x 2
x 1
Figure 2.1: Body force Fidv (dv = dx1dx2dx3, where dx1;dx2;dx3 are lengths of the edges
of the elements in x, y, and z-direction respectively) (see [6])
The force vector F = [F1;F2;F3] is called body force per unit volume. Examples
of body forces that are due to the action at a distance, are gravitational forces and
electromagnetic forces.
2. The forces acting on surface elements called stresses. Stresses can be deflned with
reference to an inflnitesimal cube with faces pointing in the coordinate directions:  ji
is the force in the xj direction, per unit area, acting on a face of the cube whose
normal points in the xi direction (see [16]).
Examples of stresses are aerodynamic pressure acting on a body and pressure due to
mechanical contact of two bodies.
6
F  
x 2
x 1
x 3
11
12
13
21
22
23
12
x1
dx112+ 
x1
dx1+ 11 11
x1
dx1+ 13 13
+ dx
x22
22
2
2
+ dx
x2 223
23
+ dx
x2 221
21
1
2
3
F  
F  
Figure 2.2: Stresses acting on surfaces and body force Fidx1dx2dx3
For example, as shown in Figure 2.2 (see [6]), the force  11dx2dx3 acts on the left hand face
of the cube, the force ( 11 + @ 11@x1 dx1)dx2dx3 acts on the right hand face of the cube, and
so forth.
Every stress component is a function of position, that is,  11 is a function of (x1;x2;x3). The
value of the stress  11 at a point slightly to the right of (x1;x2;x3), namely (x1+dx1;x2;x3),
is  11(x1+dx1;x2;x3). Now, since  11 is a continuously difierentiable function of (x1;x2;x3),
then according to Taylor?s theorem we have
 11(x1 +dx1;x2;x3) =  11(x1;x2;x3) + dx1@ 11@x
1
(x1;x2;x3)
+ dx21 12 @
2 11
@x21 (x1 +fidx1;x2;x3); (2.1)
7
where 0 ? fi ? 1. The last term can be made arbitrarily small by choosing dx1 su?ciently
small.
Neglecting the last term in equation (2.1), we get
 11(x1 +dx1;x2;x3) =  11(x1;x2;x3) + dx1@ 11@x
1
(x1;x2;x3):
Let us consider balance of forces in the x1-direction:
1
21
+ dxx
2 221
21
11
x1 dx1+ 11 1131
+ 31 31x dx3
3
x 2
x 1x 
3
F  
Figure 2.3: Stress components in x1-direction
At equilibrium, the sum of forces on the body vanishes, then we have
( 11 + @ 11@x1 dx1)dx2dx3 ?  11dx2dx3 + ( 21 + @ 21@x2 dx2)dx1dx3 ?  21dx3dx1
+ ( 31 + @ 31@x3 dx3)dx1dx2 ?  31dx1dx2 + F1dx1dx2dx3 = 0.
Simplifying and dividing by dx1dx2dx3, we obtain
@ 11
@x1 +
@ 21
@x2 +
@ 31
@x3 +F1 = 0: (2.2)
8
Repeating the same process (used in the x1-direction), in x2-direction and in x3-direction,
we get
@ 12
@x1 +
@ 22
@x2 +
@ 32
@x3 +F2 = 0; (2.3)
@ 13
@x1 +
@ 23
@x2 +
@ 33
@x3 +F3 = 0: (2.4)
Equations (2.2){(2.3) can be expressed in index notation as
3X
j=1
@ ji
@xj + Fi = 0 i = 1;2;3: (2.5)
Equation (2.5) expresses the force equilibrium where  ji is the total stress, per unit area
(in the j-direction acting on the surface with normal in the i-direction) and Fi is the body
force per unit volume.
The components  11; 22; and  33 are the normal stresses, they are perpendicular to the
face. The other three stresses are the shear stresses where the force is tangent to the face.
Rotational equilibrium on all such inflnitesimal elements of material requires that shear
stresses be equal on adjoining faces, which is concisely expressed by requiring that  ji =  ij
for all i and j (see [18]) (the symmetry of the stress tensor).
9
2.1.2 Strain
Stresses cause solids to deform. The quantities describing deformations of the body are
called strains (denoted by ?). Strains can be deflned most simply in the case of extremely
small deformations, in which case the coordinates directionsx1, x2, and x3 of material points
are virtually the same before and after deformation. For normal strain in x1-direction, let
us consider two points A and B, a small distance dx1 apart (see [6]),
B
1 u1 1xd
u1
1x
+ 
1xdA
u
Figure 2.4: Normal strain in x1-direction
and let u1 be the x1-displacement at A and u1 + @u1@x1dx1 be the x1-displacement at B.
The unit displacement in the x1-direction @u1@x1 deflnes the normal strain denoted by ?11,
?11 = @u1@x
1
:
Similarly, the unit displacement in the x2-direction and x3-direction respectively are
?22 = @u2@x
2
;
and
?33 = @u3@x
3
:
The shear strains ?12;?13;?23 are the small changes of angle between the line segments in
the x1 and x2-directions, x1 and x3-directions, and x2 and x3-directions respectively (see
10
[16]).
To illustrate this for the shear strain ?12, consider the line segements AB and AC, initially
making a right angle with dx1 a small distance between A and B and dx2 a small distance
between A and C. After deformation the points are at A?, B?, and C? and the lines A?B?
and A?C? no longer meet at a right angle at A?.
B?
1 u1+ x dx22
1xd1x+ u u22
u1
u2
A B
C
A?
C?
u
Figure 2.5: Shear strain in x1, x2-direction
The shear strain is deflned as the average of the angles @u2@x1 and @u1@x2 that A?B? and A?C?
make with the x1 and x2 directions respectively (see [6]),
?12 = 12
 @u
1
@x2 +
@u2
@x1
?
:
Similarly,
?13 = 12
 @u
1
@x3 +
@u3
@x1
?
;
11
and
?23 = 12
 @u
2
@x3 +
@u3
@x2
?
:
The normal and shear strains can be compactly written using Einstein notation, in which
repeated indices are summed, as follows
?ij = 12
 @u
i
@xj +
@uj
@xi
?
i;j = 1;2;3: (2.6)
Note that ?ij = ?ji and that strains are dimensionless since they are the ratio of lengths.
2.1.3 Stress-strain relationship
The relationship between stress and strain can be derived using the following consti-
tutive equation (see [18])
trace(?) = 13Ktrace( )+ 1Hp; (2.7)
where 1K is the compressibility of the material (K is bulk modulus) measured under drained
conditions. Drained conditions correspond to the deformation at flxed pressure p, with the
 uid being allowed to  ow in or out of the deforming element. The coe?cient 1K is obtained
( 1K = ??? jp=0) by measuring the change in volumetric strain due to changes in applied stress
while holding the pressure constant. The coe?cient 1H represents the poroelastic expansion
coe?cient. It describes how much the bulk volume changes due to a pore pressure change
( 1H = ???pj =0) while holding the applied stress constant (see [18]).
Equation (2.7), says that the fractional volume change is the result of change in applied
stress and pore pressure.
In equation (2.7), trace( )3 is the average of normal stresses and trace(?) is the volumetric
12
strain, i.e.,
trace( )
3 =
 11 + 22 + 33
3 ;
trace(?) = ?11 +?22 +?33:
Then
?11 +?22 +?33 = 1K ( 11 + 22 + 33)3 + pH:
The coe?cient K can be expressed in terms of Young?s modulus E (see [18]) by
K = E3(1?2?):
Young?s modulus is the measure of the stifiness of an elastic material and is deflned as
the ratio of the rate of change of stress with strain. The constant ? represents Poisson?s
ratio. When an elastic material is stretched or compressed in one direction, it deforms in
perpendicular directions (becoming thicker or thinner), the measure of this deformation is
given by the Poisson?s ratio ? (see [16]).
From the above we get
?11 +?22 +?33 = 1E( 11 ?? 11 ?? 11)+ 1E( 22 ?? 22 ?? 22)
+ 1E( 33 ?? 33 ?? 33)+ pH:
That is,
?11 = 1E 11 ? ?E 22 ? ?E 33 + p3H;
?22 = ??E 11 + 1E 22 ? ?E 33 + p3H;
13
?33 = ??E 11 ? ?E 22 + 1E 33 + p3H:
This form is chosen (see [8]) to express the fact that one constant, 1E, connects strain
and stress in the same direction. The other constant, ?E, relates strain and stress in two
perpendicular directions.
We use the following expressions (see [18])
E = 2G(1+?) and 1H = fiK:
The coe?cient G is the shear modulus and is a quantity measuring the strength of the
material deflned as a ratio of shear stress to the shear strain. The positive constant fi is the
Biot-Willis coe?cient, the ratio of volume of  uid that is added to storage and the change
in bulk volume under the constraint that the pore pressure remains constant. Note that
the constant  uid pressure condition means that the volume of  uid that goes into or out
of storage is equal to the change in pore volume (see [16]).
Substituting the previous two expressions into the three normal strain equations, we obtain
after simpliflcation (using that 11+? = 1? ?1+?):
?11 = 12G[ 11 ? ?1+? kk]+ fi3Kp:
Similarly,
?22 = 12G[ 22 ? ?1+? kk]+ fi3Kp;
and
?33 = 12G[ 33 ? ?1+? kk]+ fi3Kp:
14
Because changes in pore pressure are assumed not to induce shear strain, the following shear
strain and shear stress relationships are independent of pore pressure (see [18]).
?12 = 12G 12;
?23 = 12G 23;
and
?13 = 12G 13:
Expressing the previous six strain-stress equations using Einstein summation convention,
we get
?ij = 12G[ ij ? ?1+? kk?ij]+ fi3Kp?ij i;j = 1;2;3 (2.8)
where ?ij is the Kronecker delta.
?ij =
8
><
>:
1 if i = j
0 if i 6= j.
The above equation of strain in terms of stress and pore pressure may be inverted to solve
for stress, that is,
 ij = 2G?ij + ?1+? kk?ij ?2G fi3Kp?ij: (2.9)
We have
?11 = 12G[ 11 ? ?1+? kk]+ fi3Kp;
?22 = 12G[ 22 ? ?1+? kk]+ fi3Kp;
15
and
?33 = 12G[ 33 ? ?1+? kk]+ fi3Kp:
Adding these three equations yields
?kk = 12G(1?2?)(1+?)  kk + fiKp;
which implies that
 kk = 2G (1+?)(1?2?)?kk ?2G (1+?)(1?2?) fiKp:
Substituting  kk into (2.9) and simplifying, we get
 ij = 2G?ij +2G ?1?2??kk?ij ?fip?ij: (2.10)
Writing equation (2.10) explicitly for the normal stresses yields
 11 = 2G?11 +2G ?1?2??kk ?fip; (2.11)
 22 = 2G?22 +2G ?1?2??kk ?fip; (2.12)
 33 = 2G?33 +2G ?1?2??kk ?fip; (2.13)
 12 = 2G?12; (2.14)
 13 = 2G?13; (2.15)
and
 23 = 2G?23: (2.16)
16
We also have the force equilibrium equations
@ 11
@x1 +
@ 21
@x2 +
@ 31
@x3 +F1 = 0; (2.17)
@ 12
@x1 +
@ 22
@x2 +
@ 32
@x3 +F2 = 0; (2.18)
and
@ 13
@x1 +
@ 23
@x2 +
@ 33
@x3 +F3 = 0: (2.19)
Substituting the six stress equations (2.11){(2.16) into the force equilibrium equations
(2.17){(2.19), we obtain
2G@?11@x
1
+2G ?(1?2?)@?kk@x
1
+2G@?12@x
2
+2G@?13@x
3
?fi @p@x
1
+F1 = 0; (2.20)
2G@?12@x
1
+2G ?(1?2?)@?kk@x
2
+2G@?22@x
2
+2G@?23@x
3
?fi @p@x
2
+F2 = 0; (2.21)
and
2G@?13@x
1
+2G ?(1?2?)@?kk@x
3
+2G@?23@x
2
+2G@?33@x
3
?fi @p@x
3
+F3 = 0: (2.22)
We write explicitly the strain in terms of the displacement
?11 = @u1@x
1
; (2.23)
?22 = @u2@x
2
; (2.24)
?33 = @u3@x
3
; (2.25)
17
?12 = 12
 @u
1
@x2 +
@u2
@x1
?
; (2.26)
?13 = 12
 @u
1
@x3 +
@u3
@x1
?
; (2.27)
and
?23 = 12
 @u
2
@x3 +
@u3
@x2
?
; (2.28)
where the displacementsu1, u2, u3 are the displacements inx1, x2, x3 directions respectively.
We flrst substitute equations (2.23){(2.28) into equations (2.20){(2.22) and simplify to get
G
 @2u
1
@x21 +
@2u1
@x22 +
@2u1
@x23
?
+ G(1?2?)
 @2u
1
@x21 +
@2u2
@x1@x2 +
@2u3
@x1@x3
?
?fi @p@x
1
+F1 = 0;
G
 @2u
1
@x21 +
@2u1
@x22 +
@2u1
@x23
?
+ G(1?2?)
 @2u
1
@x1@x2 +
@2u2
@x22 +
@2x3
@x2@x3
?
?fi @p@x
2
+F2 = 0;
and
G
 @2u
1
@x21 +
@2u1
@x22 +
@2u1
@x23
?
+ G(1?2?)
 @2u
1
@x1@x3 +
@2u2
@x2@x3 +
@2u3
@x23
?
?fi @p@x
3
+F3 = 0:
These last three equations are expressed as
?G4ui ? G(1?2?) @
2uk
@xi@xk +fi
@p
@xi = Fi i;k = 1;2;3:
This is just the conservation of momentum and can be written in vector form as
?G4u? G(1?2?)r(r?u)+firp = F: (2.29)
18
An equivalent expression to the conservation of momentum equation (2.29) is
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F: (2.30)
Equation (2.30) is obtained by using
4u = r?(ru+ruT)?r(r?u);
since r?(ruT) = r(r?u).
Then equation (2.29) becomes
?Gr?(ru+ruT)+Gr(r?u)? G(1?2?)r(r?u)+firp = F:
By simplifying we get the equivalent equation for momentum conservation (2.30).
19
2.2 The pore pressure equation
The simplest mathematical description of the coupling of the pressure and deformation is
the constitutive equation
? = fi3Ktrace( )+ 1Rp; (2.31)
(see [18]) where ? is increment of  uid which is positive for  uid added to the control
volume and negative for  uid withdrawn from the control volume. The coe?cient 1K is
the compressibility of the material as deflned previously. The coe?cient 1R (???pj =0) is the
speciflc storage coe?cient measured under conditions of constant applied stress.
The Skempton?s coe?cient B = RH is the ratio of the induced pore pressure to the change
in applied stress for undrained condition, that is, no  uid is allowed to move into or out of
the control volume (see [18]).
Using the Biot-Willis coe?cient fi = KH and Skempton?s coe?cient B = RH, we get 1R = fiKB.
Equation (2.31) can be expressed as
? = fi3Ktrace( )+ fiKBp: (2.32)
The average velocity, v = q`, is interpreted as the relative velocity between the  uid and
solid, that is,
v = 1`q = r?(Uf ?Us); (2.33)
where Uf is the average displacement of the  uid, Us is the average displacement of the
solid, q is the  uid  ux, and ` is the porosity (see [18]).
20
The increment of  uid is expressed by Biot and Willis (see [18]) (1957) in terms of Uf and
Us as
? = ?`r?(Uf ?Us): (2.34)
Taking derivative of equation (2.34) with respect to time and substituting equation (2.33)
into it yields
@?
@t = ?r?q:
Now, substituting q from Darcy?s law: q = ?k?rp into this last equation (here k is the
permeability of the rock and ? is the viscosity) gives
@?
@t = r?(
k
?rp):
Accounting for quantity of  uid from an external source Q, we have
@?
@t ?r?(
k
?rp) = Q:
Finally, substituting equation (2.32) into the previous equation yields
@
@t
h fi
KBp+
fi
3K kk
i
?r?(k?rp) = Q:
If in this last equation, displacement is chosen as the mechanical variable instead of mean
normal stress then we get the general difiusion equation for Darcy  ow
@
@t(Se p+fir?u)?r?(
k
?rp) = Q: (2.35)
21
Equation (2.35) was derived using the following steps:
? write the normal stress  kk3 in terms of strains using equation (2.10),
? replace strain by displacement using equation (2.6),
? use K = E3(1?2?), E = 2G(1+?), and the speciflc storage Se = fiKB.
In summary, we derived the following system of partial difierential equations:
?G4u? G1?2?r(r?u)+firp = F; in ??(0;T);
@
@t(Se p+fir?u)?r?(
k
?rp) = Q; in ??(0;T):
Equivalently,
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F; in ??(0;T);
@
@t(Se p+fir?u)?r?(
k
?rp) = Q; in ??(0;T):
The flrst equation represents the equilibrium equation for conservation of momentum and
the second equation is the general difiusion equation.
22
2.3 Boundary and initial conditions
In this section, we discuss the boundary and initial conditions for the quasi-static
poroelasticity system. Let ? be a bounded open connected subset of R3 with Lipschitz
boundary. Denote the boundary by ? = @?.
The boundary ? is divided into two disjoint parts; the clamped boundary denoted by ?c
with strictly positive measure and the traction boundary denoted by ?t. The boundary, ?,
can further be divided into drained boundary ?d and the  ux boundary ?f.
Under certain geological conditions the boundary condition can belong to both ?t and ?f.
Let us denote this boundary by ?tf (?tf = ?t \?f). The boundaries ?c and ?t correspond
to the momentum equation, the boundaries ?d and ?f correspond to the  uid equation,
and there is a coupling between the two equations on ?tf.
The initial boundary value problem (IBVP) becomes (see [15]):
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F in ??(0;T); (2.36)
@
@t(Se p+fir?u)?r?(
k
?rp) = Q in ??(0;T); (2.37)
u = uc on ?c ?(0;T);(2.38)
[G(ru+(ru)T)+G 2?1?2?r?uI]^n?flfip^n?tf =  t on ?t ?(0;T); (2.39)
p = pd on ?d ?(0;T);(2.40)
? @@t((1?fl)fiu? ^n)?tf + k?rp? ^n = h1?tf on ?f ?(0;T);(2.41)
Se p+fir?u = v0 on ??f0g; (2.42)
(1?fl)fiu? ^n = v1 on ?tf ?f0g: (2.43)
23
Equations (2.36) and (2.37) are the partial difierential equations for the quasi static poroe-
lasticity system. Equation (2.36) represents the general force equilibrium equation and
equation (2.37) is the general difiusion equation.
Boundary conditions (2.38) and (2.40) correspond to the clamped boundary ?c and the
drained boundary ?d. The boundary conditions (2.39) and (2.41) consist of a balance
forces on the traction boundary ?t and a balance of  uid mass on the  ux boundary ?f.
Motivated by the geological application and for simplicity, the boundary functions uc; t;
and pd are set equal to zero.
Here I is the identity tensor and ^n is the unit outward pointing normal vector on the bound-
ary. The fraction 0 ? fl ? 1 deflned on the boundary ?tf, the portion of the boundary
which is neither clamped nor drained, denotes the surface fraction of the matrix pores which
are sealed along ?tf. The remaining portion (1?fl) is exposed along the  ux boundary ?f
and contributes to the  ux.
Here ?tf denotes the characteristic function of ?tf, that is, ?tf = 1 on ?tf and 0 otherwise.
The transverse  ow on the  ux boundary ?f is h1. More speciflcally h1 = ?(1?fl)v(t)?^n,
where v(t) is the  uid velocity on the boundary ?f.
Finally, equations (2.42) and (2.43) represent the initial conditions where v0 and v1 are the
given initial data.
In [14] Showalter showed that the system (2.36)-(2.43) has a unique strong solution under
mild (smoothness) requirements on the data, these will be clarifled later. He also proved
existence and uniqueness of weak solutions for the system.
We will use a constructive approach to prove the existence and uniqueness of weak solutions
24
for the system. Our approach immediately suggests a numerical algorithm which can be
used to approximate solutions of the quasi-static poroelasticity system.
25
Chapter 3
Analysis
In the previous Chapter, we introduced the system of partial difierential equations that
consists of the equilibrium equation for momentum conservation and the difiusion equation
for Darcy  ows. We also discussed the boundary and initial conditions for the system. In
this Chapter, we brie y recall existence and uniqueness of strong and weak solutions derived
by Showalter in [14]. We give an alternative constructive proof of existence and uniqueness
of weak solutions of the quasi-static poroelasticity system. We describe a discretization of
the problem and derive error estimates.
26
3.1 Existence and uniqueness of solutions (Showalter)
We start this section by brie y recalling existence and uniqueness results for strong
and weak solutions proved by Showalter (see [14]) for the system (3.1){(3.8).
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F in ??(0;T); (3.1)
@
@t(Se p+fir?u)?r?(
k
?rp) = Q in ??(0;T); (3.2)
u = 0 on ?c; (3.3)
[G(ru+(ru)T)+G 2?1?2?r?uI]^n?flfip^n?tf = 0 on ?t; (3.4)
p = 0 on ?d; (3.5)
? @@t((1?fl)fiu? ^n)?tf + k?rp? ^n = h1?tf on ?f; (3.6)
Se p+fir?u = v0 on ??f0g; (3.7)
(1?fl)fiu? ^n = v1 on ?tf ?f0g: (3.8)
Throughout the course of this work, we will use the following Hilbert spaces: L2(?)
which is the space of square integrable functions on ? and H1(?) which is the space of
functions in L2(?) whose flrst distribution derivatives are square integrable. The L2 inner
product and L2 and H1 norms are given by
(f;g) =
Z
?
f(x)g(x)dx;
jjfjjL2(?) =
 Z
?
jf(x)j2dx
?1
2 = ((f;f))1
2 ;
27
and
jjfjjH1(?) =
?
jjfjj2L2(?) +jjrfjj2L2(?)
?1
2
respectively.
The space H1?(?) is the closure of fv 2 C1(?)3 : v(x) = 0 for x 2 ?g with respect to the
jj?jj1-norm.
Notation: we denote the L2 inner product by (?;?), the L2 norm by jj?jj, and the H1 norm
by jj?jj1.
In addition to the above spaces, we will need the subspaces
V = fv 2 (H1(?))3 : v = 0 on ?cg
and
M = fq 2 H1(?) : q = 0 on ?dg:
For the strong solution and in the special case where Se 6= 0, Showalter?s theorem (see
theorem 3.1 [14]) becomes
Theorem 3.1 Let T > 0, v0 2 L2(?), v1 2 L2(?tf), and the Holder continuous functions
F;Q 2 Cfi([0;T];L2(?)), h1(:) 2 Cfi([0;T];L2(?tf)) be given, then there exists a pair of
functions p : (0;T] ! M and u : (0;T] ! V for which (Sep +r?u) 2 C0([0;T];L2(?))\
C1([0;T];L2(?)) and u 2 C0([0;T];V)\C1([0;T];V). The system (3.1){(3.8) is satisfled
and the function u is unique. Furthermore, if the measure of ?c is strictly positive then p
is unique.
28
In [14] Showalter also proved existence and uniqueness of weak solutions to the system
under weaker assumptions on the data compared to the strong solution. His results is given
in the following theorem (see theorem 4.1 [14]).
Theorem 3.2 Let T > 0, v0 2 M0, and Q 2 Cfi([0;T];M0) be given. Then there exists a
unique pair of functions p : (0;T] ! M and u : (0;T] ! V for which the system (3.1){(3.8)
is satisfled in a weak sense. The function u is unique. Furthermore, if the measure of ?c is
strictly positive then p is unique.
29
3.2 Existence and uniqueness of weak solutions
We will prove existence and uniqueness of weak solutions of the quasi-static poroelas-
ticity system (3.1){(3.8). We will use Rothe?s method of lines, and at each time step we
will use Babuska-Brezzi theory to show that the elliptic system has a unique solution.
We flrst derive a weak formulation of the following system of partial difierential equations:
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F; in ??(0;T);
@
@t(Se p+fir?u)?r?(
k
?rp) = Q; in ??(0;T):
Let the Hilbert spaces V and M such that:
V = fv 2 (H1(?))3 : v = 0 on ?cg;
M = fq 2 H1(?) : q = 0 on ?dg:
Multiply the previous two partial difierential equations by the test functions v 2 V and
q 2 M respectively and integrate over ?, to get
Z
?
?Gr?(ru + (ru)T)v ?
Z
?
G 2?1?2?r(r?u)v +
Z
?
firp v =
Z
?
F v 8 v 2 V;
Z
?
Septq +
Z
?
fir?utq ?
Z
?
k
??p q =
Z
?
Q q 8 q 2 M:
Applying Green?s formula:
?
Z
?
r?(ru + (ru)T)v =
Z
?
(ru + (ru)T) : rv ?
Z
?t
(ru + (ru)T)?^nv; for u;v 2 V;
30
and
?
Z
?
r(r?u)v =
Z
?
(r?u)(r?v) ?
Z
?t
(r?u)^n?v; for u;v 2 V:
Apply Green?s formula
?
Z
?
?p?q =
Z
?
rp?rq ?
Z
?f
rp? ^nq; for p;q 2 M:
Therefore, the system becomes:
Z
?
[G(ru+(ru)T) : rv +G 2?1?2?(r?u)(r?v)] +
Z
?
firp v =
Z
?
F v +
Z
?t
G(ru + (ru)T)? ^nv +
Z
?t
G 2?1?2?(r?u)^n?v;
Z
?
Se pt q +
Z
?
fi(r?u)tq +
Z
?
k
?rp?rq =
Z
?
Q q +
Z
?f
k
?rp? ^nq:
Discretizing in time using  -scheme, we get
Z
?
[G(run+1 +(run+1)T) : rv +G 2?1?2?(r?un+1)(r?v)] +
Z
?
firpn+1 v =
Z
?
Fn+1v +
Z
?t
G(run+1 + (run+1)T)? ^nv +
Z
?t
G 2?1?2?(r?un+1)^n?v;
31
Z
?
Se p
n+1 ?pn
? q +
Z
?
fir?u
n+1 ?r?un
? ?q +
Z
?
k
?( rp
n+1 +(1? )rpn)?rq =
Z
?
( Qn+1 +(1? )Qn) q +
Z
?f
k
?( rp
n+1 +(1? )rpn+1)? ^nq:
Here 0 ?  ? 1, the superscript n denotes the discrete time level at which the functions are
evaluated, and ? is the time step. That is, ? = TN where N is the number of time steps.
Hence un = u(tn) where tn = n??.
Using the divergence theorem: R?r?uq = R?f u? ^nq ? R? u?rq, and rearranging the
second equation of the previous system, we get
Z
?
[G(run+1 +(run+1)T) : rv +G 2?1?2?(r?un+1)(r?v)]+
Z
?
firpn+1 v =
Z
?
Fn+1v +
Z
?t
G(run+1 + (run+1)T)? ^nv +
Z
?t
G 2?1?2?(r?un+1)^n?v;
?
Z
?
fiun+1 ?rq +
Z
?
[Se pn+1q + k??  rpn+1 ?rq] =
Z
?
[?( Qn+1 +(1? )Qn) + fir?un + Se pn]q ?
Z
?
k?
? (1? )rp
nrq
?
Z
?f
fiun+1 ? ^nq +
Z
?f
k
?( rp
n+1 +(1? )rpn+1)? ^nq:
32
We introduce the bilinear forms a, b, and c as follows:
a(u;v) :=
Z
?
[G(ru + (ru)T) : rv + G 2?1?2?(r?u)(r?v)];
b(v;p) :=
Z
?
firp v;
c(p;q) :=
Z
?
[Se pq + k??  rp?rq];
and
l1(F;g;v) :=
Z
?
Fv +
Z
?t
G(rg + (rg)T)? ^nv +
Z
?t
2G ?1?2?(r?g)^n?v;
l2(Q1;r1;s1;q) :=
Z
?
[?(1? )Q1 + fir?r1 + Se s1]q ?
Z
?
k?
? (1? )?rs1rq;
l3(Q2;r2;s2;q) :=
Z
?
? Q2q ?
Z
?f
fir2 ? ^nq +
Z
?f
k
?( rs2 +(1? )rs2)? ^nq:
The weak formulation of this problem is: flnd (un+1;pn+1) 2 V ?M such that:
a(un+1;v) + b(v;pn+1) = l1(Fn+1;un+1;v) 8 v 2 V;
?b(un+1;q) + c(pn+1;q) = l2(Qn;un;pn;q)
+ l3(Qn+1;un+1;pn+1;q) 8 q 2 M:
That is,
a(un+1;v) + b(v;pn+1) = l1(Fn+1;un+1;v) 8 v 2 V; (3.9)
b(un+1;q) ? c(pn+1;q) = ?l2(Qn;un;pn;q)
? l3(Qn+1;un+1;pn+1;q) 8 q 2 M: (3.10)
Using Babuska-Brezzi theory (see [3]), we will show that the system (3.9){(3.10) has a
unique solution.
33
Deflnition 1 : The bilinear form a(?;?) : V?V !Rand the linear form b(?;?) : V?M !R
are continuous provided that positive constants fl and  exist such that:
ja(u;v)j? fljjujjVjjvjjV 8 u;v 2 V:
and
jb(u;v)j?  jjujjVjjvjjM 8u 2 V; 8v 2 M:
Deflnition 2 : The bilinear form a(?;?) : V ?V ! R is coercive (or V-elliptic) provided
that a positive constant fi exists such that:
a(v;v) ? fijjvjj2V 8 v 2 V:
Theorem 3.3 : If the bilinear forms a(?;?) : V ? V ! R is continuous and coercive,
c(?;?) : M ?M !R is continuous and coercive, and b(?;?) : V ?M !R is continuous then
for every f 2 V0 and g 2 M0
a(u;v) + b(v;p) = (f;v)
b(u;q) ? c(p;q) = (g;q)
has a unique solution (u;p).
We now show that the bilinear forms a(?;?) and c(?;?) are continuous and coercive and the
bilinear form b(?;?) is continuous on the respective spaces, hence there exist unique solutions
u and p of the semi-discrete problem.
34
Recall that:
ru : rv =
X
i;j
@ui
@xj ?
@vi
@xj;
r?u =
X
i
@ui
@xi;
and
V = fv 2 (H1(?))3 : v = 0 on ?cg;
M = fq 2 H1(?) : q = 0 on ?dg:
Continuity and coercivity of a(?;?)
ja(u;v)j =
flfl
flfl
Z
?
Gru : rv + GruT : rv + G 2?1?2?(r?u)(r?v)
flfl
flfl
Using the triangle inequality
ja(u;v)j?
flfl
flfl
flfl
Z
?
G
X
i;j
 @u
i
@xj
? @v
i
@xj
?flflfl
flfl
fl
+
flfl
flfl
flfl
Z
?
G
X
i;j
 @u
j
@xi
? @v
i
@xj
?flflfl
flfl
fl
+
flfl
flfl
fl
Z
?
G 2?1?2?
X
i
@ui
@xi
X
i
@vi
@xi
flfl
flfl
fl
The inner product
(ru;rv) =
Z
?
X
i;j
 @u
i
@xj
? @v
i
@xj
?
?
Z
?
jrujjrvj
? jjrujj jjrvjj (By H?older?s inequality (Appendix A))
? jjujj1jjvjj1 (Since jjujj1 = jjrujj+jjujj so jjrujj?jjujj1):
35
Thus
j a(u;v) j ? Gjjrujj jjrvjj + Gjjrujj jjrvjj + G 2?1?2?jjr?ujj jjr?vjj:
Since jjr?ujj?p3jjrujj (is shown below),
j a(u;v) j ? 2Gjjrujj jjrvjj + 6G?1?2?jjrujj jjrvjj;
hence,
j a(u;v) j? max
 
2G; 6G?1?2?
?
jj u jj1jj v jj1 8 u;v 2 (H1(?))3: (3.11)
Hence a(u;v) is continuous.
The inequality jjr?ujj?p3jjrujj
We have
r?u =
X
i
@ui
@xi;
thus,
(r?u)(r?v) =
 @u
1
@x1 +
@u2
@x2 +
@u3
@x3
? @v
1
@x1 +
@v2
@x2 +
@v3
@x3
?
;
(r?u)(r?v) = @u1@x
1
@v1
@x1 +
@u1
@x1
@v2
@x2 +
@u1
@x1
@v3
@x3 +
@u2
@x2
@v1
@x1 +
@u2
@x2
@v2
@x2 +
@u2
@x2
@v3
@x3
+ @u3@x
3
@v1
@x1 +
@u3
@x3
@v2
@x2 +
@u3
@x3
@v3
@x3:
Using the fact that ab ? 12a2 + 12b2
(r?u)(r?v) ? 32
 @u
1
@x1
?2
+32
 @u
2
@x2
?2
+32
 @u
3
@x3
?2
+32
 @v
1
@x1
?2
+32
 @v
2
@x2
?2
+32
 @v
3
@x3
?2
:
36
Therefore,
(r?u)(r?v) ? 32
X
i
 @u
i
@xi
?2
+ 32
X
i
 @v
i
@xi
?2
:
That is,
(r?u)(r?u) ? 3
X
i
 @u
i
@xi
?2
;
thus,
jjr?ujj2 ? 3
Z
?
X
i
 @u
i
@xi
?2
? 3jjrujj2:
Hence
jjr?ujj ?p3jjrujj: (3.12)
Coercivity of a(?;?),
a(u;v) =
Z
?
G(ru + ruT) : rv + G 2?1?2?(r?u)(r?v);
and
(ru + ruT) : rv = ru : rv + ruT : rv
= 12[ru : rv + ruT : rvT + ruT : rv + ru : rvT]
= 12(ru + ruT) : (rv + rvT): (3.13)
Note that the strain
?(u) = 12(ru + ruT);
37
and
a(u;v) =
Z
?
2G(12(ru + ruT))(12(rv + rvT))+G 2?1?2?(r?u)(r?v);
where again we use Einsteins? summation convention.
Therefore,
a(v;v) ?
Z
?
2G(?(v))2:
Korn?s inequality (see [2]): Let ? ? R3 be an open bounded set with piecewise smooth
boundary. In addition, suppose ?0 ? @? has positive two-dimentional measure. Then there
exists a positive constant c = c(?;?0) such that:
Z
?
?(v) : ?(v) ? cjjvjj21 for all v 2 H1?(?):
Since we assumed that the measure of the clamped boundary ?c is positive, then from
Korn?s inequality a(?;?) is coercive.
a(v;v) ? 2G jj v jj21 8 v 2 (H1(?))3; (3.14)
where G > 0 is the shear modulus.
Furthermore, the bilinear form a(:;:) is symmetric. Recall that
a(u;v) =
Z
?
G(ru + ruT) : rv + G 2?1?2?(r?u)(r?v);
38
and from (3.13)
(ru + ruT) : rv = 12(ru + ruT) : (rv + rvT):
Therefore,
a(u;v) =
Z
?
1
2G(ru + ru
T) : (rv + rvT) + G 2?
1?2?(r?u)(r?v)
=
Z
?
1
2G(ru : rv +ru : rv
T +ruT : rv +ruT : rvT) + G 2?
1?2?(r?v)(r?u)
=
Z
?
1
2G(rv : ru+rv
T : ru+rv : ruT +rvT : ruT) + G 2?
1?2?(r?v)(r?u)
=
Z
?
1
2G(rv + rv
T) : (ru + ruT) + G 2?
1?2?(r?v)(r?u)
=
Z
?
G(rv + rvT) : ru + G 2?1?2?(r?v)(r?u)
= a(v;u): (3.15)
Continuity of b(?;?)
Recall that
b(v;p) :=
Z
?
fi(rp)v;
hence,
j b(v;p) j ? fijjrpjj jjvjj
? fi jj p jj1jj v jj1 8 p 2 H1(?) and v 2 (H1(?))3; (3.16)
where the constant fi > 0 is the Biot-Willis coe?cient.
39
Continuity and coercivity of c(?;?)
We have
c(p;q) :=
Z
?
Se p q + k??  (rp)(rq);
thus,
j c(p;q) j ? Sejjpjj jjqjj+ k??  jjrpjj jjrqjj
? max(Se; k??  ) jj p jj1jj q jj1 8 p;q 2 H1(?): (3.17)
Thus c is continuous.
Now,
c(q;q) =
Z
?
Se q2 + k??  (rq)2
? Sejjqjj2 + k??  jjrqjj2
? min(Se; k?? ) jj q jj21 8 q 2 H1(?): (3.18)
Thus c is coercive.
corollary 1 : The bilinear form a(?;?) : V ?V ! R is continuous. If the measure of the
clamped boundary ?c is positive, then a(?;?) is coercive. The bilinear forms c(?;?) : M?M !
Ris continuous and coercive, and b(?;?) : V ?M !Ris continuous. Hence for every F 2 V0
and Q 2 M0 the semi-discrete system (3.9){(3.10) has a unique weak solution (u;p).
40
3.3 Rothe?s method of lines
Using the previous result (Corollary 1), we can now use Rothe?s method to prove
existence and uniqueness of weak solutions of the equations of poroelasticity:
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F in ??(0;T);(3.19)
@
@t(Se p+fir?u)?r?(
k
?rp) = Q in ??(0;T);(3.20)
u = 0 on ?c; (3.21)
[G(ru+(ru)T)+G 2?1?2?r?uI]^n?flfip^n?tf = 0 on ?t; (3.22)
p = 0 on ?d; (3.23)
? @@t((1?fl)fiu? ^n)?tf + k?rp? ^n = h1?tf on ?f; (3.24)
Se p+fir?u = v0 on ??f0g; (3.25)
(1?fl)fiu? ^n = v1 on ?tf ?f0g:(3.26)
As deflned above, ? is a bounded open connected subset of R3 with Lipschitz boundary
and T is a positive time. Given functions F, Q 2 C0;1(0;T;L2(?)), h1 2 C0;1(0;T;L2(?f)),
v0 2 L2(0;T;H1(?)) and v1 2 L2(0;T;L2(?tf)).
Since the problem (3.19){(3.26) is linear its solution can be written as the sum of the
solutions of the following three problems.
41
Probelem I:
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F in ??(0;T); (3.27)
@
@t(Se p+fir?u)?r?(
k
?rp) = Q in ??(0;T); (3.28)
u = 0 on ?c; (3.29)
[G(ru+(ru)T)+G 2?1?2?r?uI]^n?flfip^n?tf = 0 on ?t; (3.30)
p = 0 on ?d; (3.31)
? @@t((1?fl)fiu? ^n)?tf + k?rp? ^n = 0 on ?f; (3.32)
Se p+fir?u = 0 on ??f0g; (3.33)
(1?fl)fiu? ^n = 0 on ?tf ?f0g; (3.34)
problem II:
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = 0 in ??(0;T); (3.35)
@
@t(Se p+fir?u)?r?(
k
?rp) = 0 in ??(0;T); (3.36)
u = 0 on ?c; (3.37)
[G(ru+(ru)T)+G 2?1?2?r?uI]^n?flfip^n?tf = 0 on ?t; (3.38)
p = 0 on ?d; (3.39)
? @@t((1?fl)fiu? ^n)?tf + k?rp? ^n = 0 on ?f; (3.40)
Se p+fir?u = v0 on ??f0g; (3.41)
(1?fl)fiu? ^n = v1 on ?tf ?f0g;(3.42)
42
and problem III:
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = 0 in ??(0;T); (3.43)
@
@t(Se p+fir?u)?r?(
k
?rp) = 0 in ??(0;T); (3.44)
u = 0 on ?c; (3.45)
[G(ru+(ru)T)+G 2?1?2?r?uI]^n?flfip^n?tf = 0 on ?t; (3.46)
p = 0 on ?d; (3.47)
? @@t((1?fl)fiu? ^n)?tf + k?rp? ^n = h1?tf on ?f; (3.48)
Se p+fir?u = 0 on ??f0g; (3.49)
(1?fl)fiu? ^n = 0 on ?tf ?f0g: (3.50)
3.3.1 Existence and uniqueness of weak solutions for
homogeneous initial and boundary conditions
We flrst consider problem (3.27){(3.34) which has homogeneous boundary and initial con-
ditions.
?Gr?(ru+(ru)T)?G 2?(1?2?)r(r?u)+firp = F; in ??(0;T);(3.51)
@
@t(Se p+fir?u)?r?(
k
?rp) = Q; in ??(0;T):(3.52)
Construct a mesh d1 on the interval I = [0;T]; divide I into m subintervals Ij := [tj?1;tj]
each of length h = Tm and tj = jh, j = 1;:::;m.
43
Using flnite difierence backward time discretization, we get
?Gr?(ruj +(ruj)T)?G 2?(1?2?)r(r?uj)+firpj = Fj; in ??(0;T);
Se
h (pj ?pj?1)+
fi
h(r?uj ?r?uj?1)?r?
k
?rpj = Qj; in ??(0;T):
Let wj 2 L2(0;t;V) and zj 2 L2(0;t;M) be the approximates solutions of the system, i.e.,
wj = uj and zj = pj, for j = 1;:::;m. Here wj = w(tj) and zj = z(tj), so the system (in
terms of wj and zj) is
?Gr?(rwj +(rwj)T)?G 2?(1?2?)r(r?wj)+firzj = Fj; in ??(0;T);
Se
h (zj ?zj?1)+
fi
h(r?wj ?r?wj?1)?r?
k
?rzj = Qj; in ??(0;T):
The weak formulation of this problem is: flnd wj 2 L2(0;t;V) and zj 2 L2(0;t;M) such
that:
Z
?
[G(rwj + (rwj)T) : rv + G 2?1?2?(r?wj)(r?v)] +
Z
?
firzj ?v =
Z
?
F v +
Z
?t
[G(rwj + (rwj)T)+G 2?1?2?r?wjI]^n?v; 8v 2 V;
Z
?
Se
h (zj ?zj?1)q +
Z
?
fi
h(r?wj ?r?wj?1)q +
Z
?
k
?rzj ?rq =
Z
?
Q q +
Z
?f
k
?rzj ? ^nq; 8q 2 M:
44
We have (using the divergence theorem R? zj r?v = R?t zj ? ^nv ? R?rzj ?v)
Z
?
[G(rwj + (rwj)T) : rv + G 2?1?2?(r?wj)(r?v)] ?
Z
?
fizjr?v =
Z
?
Fjv +
Z
?t
flfizj^n?tf ?v ?
Z
?t
fizj^n?v; (3.53)
Z
?
fi(r?wj ?r?wj?1)q +
Z
?
Se(zj ?zj?1)q +h
Z
?
k
?rzj ?rq =
h
Z
?
Qjq + h
Z
?f
h1?tfq +
Z
?f
(1?fl)fi(wj ?wj?1)?tf ?q: (3.54)
Let
a(wj;v) =
Z
?
[G(rwj + (rwj)T) : rv + G 2?1?2?(r?wj)(r?v)];
b(v;zj) = ?
Z
?
fizjr?v;
c(zj;q) =
Z
?
[Se zjq + k?hrzj ?rq]:
Hence the system (3.51){(3.52) is in the form
a(wj;v) + b(v;zj) =
Z
?
Fjv ?
Z
?tf
(1?fl)fizj^nv; (3.55)
b(wj;q) ? c(zj;q) = h
Z
?
Qjq + h
Z
?tf
h1q +
Z
?tf
(1?fl)fi(wj ?wj?1)q
+
Z
?
(fir?wj?1 +Se zj?1)q; (3.56)
It was shown in section 3.2 that the bilinear forms a(:;:) and c(:;:) are continuous and
coercive and the linear form b(:;:) is continuous. Therefore, from Corollary 1 (3.55)-(3.56)
has a unique solution (zj; wj) 2 V ?M.
45
The functions zj 2 M and wj 2 V, j = 1;:::;m, are the approximates to the functions p
and u.
Deflne the Rothe functions for the mesh d1 (recall that the mesh d1 corresponds to the
division of the interval I into m subintervals Ij) by
p1(x;t) = zj?1 + zj ?zj?1h (t?tj?1);
u1(x;t) = wj?1 + wj ?wj?1h (t?tj?1);
for all t 2 Ij = [tj?1;tj] and j = 1;:::;m.
Instead of the mesh d1 (division of the interval I into m subintervals Ij of lengths h = Tm),
consider the mesh dn, n = 2;3;:::, which consists of m2n?1 subintervals Inj := [tnj?1;tnj],
j = 1;:::;m2n?1, each of length hn = Tm2n?1. (Note that the superscript n corresponds to
the mesh dn).
The Rothe functions pn and un which correspond to the mesh dn are deflned as follows
pn(x;t) = znj?1 + z
nj ?znj?1
hn (t?t
n
j?1);
un(x;t) = wnj?1 + w
nj ?wnj?1
hn (t?t
n
j?1);
for all t 2 Inj , j = 1;:::;m2n?1.
We constructed the sequences fpn(x;t)g and fun(x;t)g, we will show that these sequences
converge to the solution p(x;t) and u(x;t) of problem I.
46
Consider the system of equations (3.53)-(3.54)
Z
?
[G(rwj + (rwj)T) : rv + G 2?1?2?(r?wj)(r?v)] ?
Z
?
fizjr?v =
Z
?
Fjv +
Z
?t
flfizj^n?tfv ?
Z
?t
fizj^nv
Z
?
Se(zj ?zj?1)q +
Z
?
fi(r?wj ?r?wj?1)q +h
Z
?
k
?rzj ?rq =
h
Z
?
Qjq + h
Z
?f
h1?tfq +
Z
?f
(1?fl)fi(wj ?wj?1)?tfq
We flrst consider the case that F = 0 and then return to the case of an arbitrary F.
Recall the bilinear form
a(u;v) = G
??
ru+ruT : rvrvT?+ 2?1?2? (r?u;r?v)
?
;
it induces a norm
[[u]] =
?
Gjuj21 +GjuTj21 +G 2?1?2?jjr?ujj20
?1
2 :
Thus
a(u;v) ? [[u]] [[v]] and a(u;u) = [[u]]2:
Recall (3.11) and (3.14)
a(u;v) ? max
 
2G; 6G?1?2?
?
jjujj1jjvjj1 the continuity of a;
a(u;u) = [[u]]2 ? 2Gjjujj21 the coercivity of a:
47
Let us denote the boundary integrals by < ?;? >, that is,
< f;g >=
Z
?tf
f^n?g:
Then the system (3.51)-(3.52) (with F = 0) can be written at time j, (j = 1;:::;m) as
a(wj;v) ? fi(zj;r?v) = ?(1?fl)fi < zj;v > (3.57)
Se(zj ?zj?1;q) + fi(r?wj ?r?wj?1;q) + hk?(rzj;rq)
= h(Qj;q)+ (1?fl)fi < wj ?wj?1;q > (3.58)
Recall that in problem I, we have homogeneous boundary conditions, i.e., h1(t) = 0 and
homogeneous initial conditions, i.e., v0 = 0 and v1 = 0.
Equation (3.57) can be written at time (j ?1) as
a(wj?1;v) ? fi(zj?1;r?v) = ?(1?fl)fi < zj?1;v > (3.59)
Subtracting equation (3.59) from equation (3.57) and using equation (3.58), we get
a(wj ?wj?1;v) ? fi(zj ?zj?1;r?v) = ?(1?fl)fi < zj ?zj?1;v >; (3.60)
Se(zj ?zj?1;q) + fi(r?wj ?r?wj?1;q) + hk?(rzj;rq)
= h(Qj;q)+ (1?fl)fi < wj ?wj?1;q >; (3.61)
48
for j = 1;:::;m. Assuming that the the system (3.51)-(3.52) holds at time zero.
Setting j = 1 in (3.60) with v = w1 ?w0 (this is possible since v 2 V), we obtain
a(w1?w0;w1?w0) ? fi(z1?z0;r?w1?r?w0) = ?(1?fl)fi < z1?z0;w1?w0 >; (3.62)
Substituting q = z1 ? z0 and q = z0 (again we can do this since q 2 M) into (3.61)
respectively, we get
Se(z1 ?z0;z1 ?z0) + fi(r?w1 ?r?w0;z1 ?z0) + hk?(rz1;rz1 ?rz0)
= h(Q1;z1 ?z0) + (1?fl)fi < w1 ?w0;z1 ?z0 >; (3.63)
and
Se(z1 ?z0;z0) + fi(r?w1 ?r?w0;z0) + hk?(rz1;rz0)
= h(Q1;z0) + (1?fl)fi < w1 ?w0;z0 > : (3.64)
Using equation (3.59) with v = w1 ?w0, we have
a(w0;w1 ?w0) ? fi(z0;r?w1 ?r?w0) = ?(1?fl)fi < w1 ?w0;z0 > : (3.65)
Adding (3.62){(3.65) and simplifying, we get
[[w1?w0]]2 + Sejjz1?z0jj2 + a(w0;w1?w0) + Se(z1?z0;z0) + hk?jjrz1jj2 = h(Q1;z1):
49
Since hk?jjrz1jj2 ? 0,
[[w1 ?w0]]2 + Sejjz1 ?z0jj2 ? [[w0]] [[w1 ?w0]] + Sejjz1 ?z0jj jjz0jj + hjjQ1jj jjz1jj;
hence,
[[w1 ?w0]]2 +Sejjz1 ?z0jj2 ?
?
[[w1 ?w0]]2 +Sejjz1 ?z0jj2
?1
2
?
[[w0]]2 +Sejjz0jj2
?1
2
+ hjjQ1jj jjz1jj: (3.66)
Taking (3.57){(3.58) with j = 1 and v = w1, q = z1 we get
a(w1;w1) ? fi(z1;r?w1) = ?(1?fl)fi < z1;w1 >;
Se(z1 ?z0;z1) + fi(r?w1 ?r?w0;z1) + hk?(rz1;rz1)
= h(Q1;z1) + (1?fl)fi < w1 ?w0;z1 > :
Adding these two equations, we have
a(w1;w1) ? fi(r?w0;z1) ? Se(z0;z1) + Se(z1;z1) + hk?(rz1;rz1)
= h(Q1;z1)?(1?fl)fi < w0;z1 >,
that is,
a(w1;w1) ? (Se z0 +fir?w0;z1) + Se(z1;z1) + hk?(rz1;rz1)
= h(Q1;z1)? < (1?fl)fiw0;z1 > :
Using the initial conditions Se z0 +fir?w0 = 0 and (1?fl)fiw0 ? ^n = 0, we get
[[w1]]2 + Sejjz1jj2 ? hjjQ1jj jjz1jj;
which implies
jjz1jj ? hjjQ1jjSe : (3.67)
50
The initial conditions Se z0 +fir?w0 = v0 = 0 obviously implies that
? fi(r?w0;z0) = Se(z0;z0): (3.68)
From (3.59) with j = 1 and v = w0, we have
a(w0;w0) ? fi(z0;r?w0) = ?(1?fl)fi < z0;w0 > :
Using (3.68) and homogeneous initial condition ((1?fl)fiw0 ? ^n = 0) implies that
a(w0;w0) + Se(z0;z0) = 0:
Therefore,
[[w0]]2 + Sejjz0jj2 = 0; (3.69)
since each term on the left hand side of equation (3.69) is positive
w0 = z0 = 0: (3.70)
Hence (3.66) becomes
[[w1 ?w0]]2 + Sejjz1 ?z0jj2 ? hjjQ1jj jjz1jj;
and using (3.67)
[[w1 ?w0]]2 + Sejjz1 ?z0jj2 ? h2jjQ1jj
2
Se : (3.71)
51
At time j, (j = 2;:::;m),
a(wj;v) ? fi(zj;r?v) = ?(1?fl)fi < zj;v >; (3.72)
Se(zj ?zj?1;q)+fi(r?wj ?r?wj?1;q)+hk?(rzj;rq)
= h(Qj;q) + (1?fl)fi < wj ?wj?1;q >; (3.73)
and at time (j ?1), (j = 2;:::;m),
a(wj?1;v) ? fi(zj?1;r?v) = ?(1?fl)fi < zj?1;v >; (3.74)
Se(zj?1 ?zj?2;q) + fi(r?wj?1 ?r?wj?2;q) + hk?(rzj?1;rq)
= h(Qj?1;q) + (1?fl)fi < wj?1 ?wj?2;q > : (3.75)
Subtract (3.74) from (3.72) and (3.75) from (3.73) to get
a(wj ?wj?1;v) ? fi(zj ?zj?1;r?v) = ?(1?fl)fi < zj ?zj?1;v >; (3.76)
Se(zj ?zj?1;q)?Se(zj?1 ?zj?2;q)+fi(r?wj ?r?wj?1;q)
?fi(r?wj?1 ?r?wj?2;q)+hk?(rzj ?rzj?1;rq)
= h(Qj;q)?h(Qj?1;q)+(1?fl)fi < wj ?wj?1;q >
?(1?fl)fi < wj?1 ?wj?2;q > : (3.77)
Adding equations (3.76) and (3.77) with v = wj ? wj?1, q = zj ? zj?1, subtracting the
following equation ((3.76) with v = wj?1 ?wj?2)
a(wj ?wj?1;wj?1 ?wj?2) ? fi(zj ?zj?1;r?wj?1 ?r?wj?2) =
?(1?fl)fi < zj ?zj?1;wj?1 ?wj?2 >;
52
and simplifying, we obtain
[[wj ?wj?1]]2 + Sejjzj ?zj?1jj2
? [[wj ?wj?1]] [[wj?1 ?wj?2]] + Sejjzj ?zj?1jj jjzj?1 ?zj?2jj
+ h(jjQjjj?jjQj?1jj)jjzj ?zj?1jj
?
h
[[wj ?wj?1]]2 + Sejjzj ?zj?1jj2
i1
2
h
[[wj?1 ?wj?2]]2 + Se
?
jjzj?1 ?zj?2jj + hjjQjjj ? jjQj?1jjSe
?2i1
2.
Squaring both sides and simplifying, we get
[[wj ?wj?1]]2 + Sejjzj ?zj?1jj2
? [[wj?1 ?wj?2]]2 + Se
?
jjzj?1 ?zj?2jj + hjjQjjj ? jjQj?1jjSe
?2
[[wj ?wj?1]]2 +Sejjzj ?zj?1jj2 ? [[wj?1 ?wj?2]]2 +Sejjzj?1 ?zj?2jj2
+h2
?jjQ
jjj?jjQj?1jj
?2
Se
+2h(jjQjjj?jjQj?1jj)jjzj?1 ?zj?2jj: (3.78)
Recalling inequality (3.71)
[[w1 ?w0]]2 + Sejjz1 ?z0jj2 ? h2jjQ1jj
2
Se ;
and setting j = 2 in (3.78)
[[w2 ?w1]]2 + Sejjz2 ?z1jj2
? [[w1 ?w0]]2 +Sejjz1 ?z0jj2 +h2
?
jjQ2jj?jjQ1jj
?2
Se
+2h?jjQ2jj?jjQ1jj?jjz1 ?z0jj
? h2jjQ1jj2Se +h2
?
jjQ2jj?jjQ1jj
?2
Se +2h
?jjQ
2jj?jjQ1jj
?jjz
1 ?z0jj;
53
also from inequality (3.71): jjz1 ?z0jj ? hjjQ1jjSe , hence,
[[w2 ?w1]]2 + Sejjz2 ?z1jj2
? h2jjQ1jj2Se + h2
?
jjQ2jj?jjQ1jj
?2
Se + 2h
2?jjQ2jj?jjQ1jj?jjQ1jj
Se
= h2Se
?
jjQ1jj + ?jjQ2jj?jjQ1jj?
?2
= h2SejjQ2jj2:
Repeating the same process, we get for j = 2;:::;m
[[wj ?wj?1]]2 + Sejjzj ?zj?1jj2 ? h
2
SejjQjjj
2: (3.79)
Let us now deflne the norm:
jjj(wj;zj)?(wj?1;zj?1)jjj =
?
[[wj ?wj?1]]2 + Sejjzj ?zj?1jj2
?1
2
So,
jjj(wj;zj)?(w0;z0)jjj = jjj(wj;zj)?(wj?1;zj?1)+ ??? +(w1;z1)?(w0;z0)jjj
? jjj(wj;zj)?(wj?1;zj?1)jjj+ ??? +jjj(w1;z1)?(w0;z0)jjj;
by the triangle inequality. Then using (3.71) and (3.79), we obtain
jjj(wj;zj)?(w0;z0)jjj ? hpSe
?
jjQ1jj + jjQ2jj + ??? + jjQjjj
?
: (3.80)
Since Q(t) 2 C0;1(0;T;L2(?)), then for all t in I, there exists a constant d such that
jjQ(t+h)?Q(t)h jj? d, for all t; t+h 2 I, see (see [10]). Then jjQ(t)jj is a continuous function
54
on I and so jjQ(t)jj attains a maximum on I, say jjQjj, i.e.,
max
t 2 I
jjQ(t)jj = jjQjj:
From (3.80),
jjj(wj;zj)?(w0;z0)jjj ? jhjjQjjpSe;
and
jjj(wj;zj)?(w0;z0)jjj2 ? j2h2jjQjj
2
Se :
Therefore
[[wj ?w0]]2 + Sejjzj ?z0jj2 ? j2h2jjQjj
2
Se :
Using the fact that z0 = w0 = 0 (from (3.70))
jjzjjj ? jhjjQjjSe and jjwjjj1 ? jh jjQjjp2GSe:
Since h = Tm,
jjzjjj ? TjjQjjSe jjwjjj1 ? T jjQjjp2GSe: (3.81)
The estimates in (3.81) are obviously independent of h, thus remain valid for an arbitrary
mesh dn. Thus for every positive integer n and j = 1;:::;m2n?1, we have
jjznj jj ? TjjQjjSe ; jjwnjjj1 ? T jjQjjp2GSe: (3.82)
55
Let Zj = zj?zj?1h and Wj = wj?wj?1h , j = 1;??? ;m. From (3.79),
[[wj ?wj?1]]2 + Sejjzj ?zj?1jj2 ? h2jjQjjj
2
Se ;
we get,
[[Wj]]2 + SejjZjjj2 ? jjQjjj
2
Se :
Using (3.14), [[u]]2 > 2Gjjujj21, we have
2GjjWjjj21 + SejjZjjj2 ? jjQjjj
2
Se :
Therefore,
jjZjjj? jjQjjjSe ? jjQjjSe and jjWjjj1 ? jjQjjjp2GSe ? jjQjjp2GSe: (3.83)
Again, the estimates in (3.83) are independent of h, and so remain valid for an arbitrary
mesh dn. Thus
jjZnj jj? jjQjjSe and jjWnj jj1 ? jjQjjp2GSe: (3.84)
From (3.82) and (3.84), we see that the norms (in L2(?) and H1(?)) of the functions znj ,
Znj and wnj , Wnj are uniformly bounded with respect to j and n, thus independently of the
mesh dn. Hence
jjwnjjj1 ? c1; 8j = 0;1;??? ;m2n?1; n = 1;2;??? ; (3.85)
jjznj jj? c2; 8j = 0;1;??? ;m2n?1; n = 1;2;??? : (3.86)
56
We now examine the Rothe sequences fun(x;t)g and fpn(x;t)g in the spaces L2(I;V) and
L2?I;L2(?)? of abstract functions which are square integrable in the Bochner sense. See
Appendix B for the deflnitions of abstract function, Bochner integral and square integra-
bility in the Bochner sense.
The Rothe functions are
un(x;t) = wnj?1 + (t?tnj?1)w
nj ?wnj?1
hn ; in I
n
j = [t
n
j?1;t
n
j];
pn(x;t) = znj?1 + (t?tnj?1)z
nj ?znj?1
hn ; in I
n
j = [t
n
j?1;t
n
j]:
Since 0 ? t?t
n
j?1
hn ? 1 in I
nj , for arbitrary t 2 I,
jjun(t)jjV =
flfl
flfl
flfl
flfl
 
1? t?t
nj?1
hn
?
wnj?1 + t?t
nj?1
hn w
n
j
flfl
flfl
flfl
flfl
V
?
flfl
flfl
flfl
flfl
 
1? t?t
nj?1
hn
?
wnj?1
flfl
flfl
flfl
flfl
V
+
flfl
flfl
flfl
flflt?t
nj?1
hn w
n
j
flfl
flfl
flfl
flfl
V
? c1
 
1? t?t
nj?1
hn
?
+ c1t?t
nj?1
hn
jjun(t)jjV ? c1:
Similarly, we get that jjpn(t)jj? c2.
Hence from (B.1),
jjun(t)jj2L2(I;V) =
Z T
0
jjun(t)jj2V dt ? c21T;
and
jjpn(t)jj2L2(I;L2(?)) =
Z T
0
jjpn(t)jj2dt ? c22T:
57
Therefore the Rothe sequences fung and fpng are bounded in the spaces L2(I;V) and
L2?I;L2(?)? respectively. Since these spaces are Hilbert spaces (see [10]), there exist sub-
sequences funkg and fpnkg which converge weakly to abstract functions u in L2(I;V) and
p in L2?I;L2(?)? respectively.
Let Wnj (t) = w
n
j (t)?w
n
j?1(t)
hn and Z
nj (t) = znj (t)?znj?1(t)
hn , we have
un(t) = wnj?1 + (t?tnj?1)Wnj ; in Inj ; (3.87)
pn(t) = znj?1 + (t?tnj?1)Znj ; in Inj : (3.88)
Deflne the abstract functions Un(t) : I ! H1(?) and Pn(t) : I ! L2(?) by
Un(0) = Wn1
Un(t) = Wnj for t 2 ~Inj := (tnj?1;tnj]; j = 1;??? ;m2n?1; (3.89)
and
Pn(0) = Zn1
Pn(t) = Znj for t 2 ~Inj =: (tnj?1;tnj]; j = 1;??? ;m2n?1: (3.90)
Since jjWnj jj1 ? jjQjjp2GSe and jjZnj jj? jjQjjSe , the sequences fUng and fPng are bounded in the
spaces L2?I;H1(?)? and L2?I;L2(?)? respectively. Since these spaces are Hilbert spaces,
there exist subsequences fUnkg and fPnkg which converge weakly to U in L2?I;H1(?)? and
to P in L2?I;L2(?)? respectively (see [10]).
58
Hence the integrals
Z t
0
U(?)d? = r(t) and
Z t
0
P(?)d? = s(t) (3.91)
exist. From (3.87), (3.88) and (3.89), (3.90), we get that
Z t
0
Unk(?)d? = unk(t) and
Z t
0
Pnk(?)d? = pnk(t): (3.92)
We now show that
r = u in L2(I;H1(?)) and s = p in L2(I;L2(?)): (3.93)
It is su?cient to show that unk * r in L2?I;H1(?)?. We show that
limn
k!1
Z T
0
?u
nk(t);v(t)
?dt ? Z T
0
?r(t);v(t)?dt = 0 8v 2 L2?I;H1(?)?:
Let v(t) be the constant v in H1(?) for all t in I, then
limn
k!1
Z T
0
?u
nk(t)?r(t);v
? = lim
nk!1
?Z T
0
?U
nk(?)?U(?)
?d?;v? by (3.92) and (3.91)
= limn
k!1
Z T
0
?
Unk(?)?U(?);v
?
d?
= 0 (since Unk * U)
We can now apply the Lebesgue theorem that
0 = limn
k!1
Z T
0
?
unk(t)?r(t);v
?
dt =
Z T
0
h
limn
k!1
?u
nk(t)?r(t);v
?idt:
59
implies
unk * r:
The same proof also applies for piecewise constant functions v(t), t 2 I and so for every func-
tion v 2 L2(I;H1(?)) (since the piecewise constant functions are dense in L2(I;H1(?))).
Hence r = u.
In the same manner, it can be shown that s = p.
From (3.91) and (3.93), we get that Rt0 U(?)d? = u and Rt0 P(?)d? = p, hence
u 2 AC(I;H1(?)); p 2 AC(I;L2(?))
and
ut(t) = U(t); pt(t) = P(t);
in H1(?) and L2(?), respectively, for almost all t 2 I.
Since u(t) = Rt0 U(?)d? and p(t) = Rt0 P(?)d? then u(0) = 0 and p(0) = 0 in C(I;H1(?))
and C(I;L2(?)). Thus the initial conditions of the problem are satisfled. Since u 2 L2(I;V)
and p 2 L2(I;L2(?)), then for almost all t 2 I, u 2 V and p 2 M, which imply that the
boundary conditions are satisfled in the sense of traces.
We now have to show that the functions u and p satisfy the system of partial difierential
equations.
Deflne sequences f~unk(t)g and f~pnk(t)g by
~unk(0) = wn1
~unk(t) = wnj for t 2 ~Inj = (tnj?1;tnj]; j = 1;??? ;m2n?1;
60
and
~pnk(0) = zn1
~pnk(t) = znj for t 2 ~Inj = (tnj?1;tnj]; j = 1;??? ;m2n?1;
We show that if unk * u in L2(I;V) and pnk * p in L2(I;L2(?)), then ~unk * u in L2(I;V)
and ~pnk * p in L2(I;L2(?)).
We show that
limn
k!1
Z T
0
?u(t)? ~u
nk(t);v(t)
?
V dt = 0 8v 2 L
2(I;V)
and
limn
k!1
Z T
0
?p(t)? ~p
nk(t);q(t)
?
Mdt = 0 8q 2 L
2(I;M)
limn
k!1
Z T
0
?u(t)? ~u
nk(t);v(t)
?
V dt = limnk!1
Z T
0
?
u(t)?unk(t)+unk(t)? ~unk(t);v(t)
?
V
dt
= limn
k!1
Z T
0
?u(t)?u
nk(t);v(t)
?
V dt
+ limn
k!1
Z T
0
?u
nk(t)? ~unk(t);v(t)
?
V dt:
The limit of the flrst term on the right hand side is zero since unk * u, then we have to
show that the limit of the second term is equal to zero.
Let K be a set of abstract functions v 2 L2(I;V) such that v = g where g 2 V is a certain
function on an interval [fi;fl] ? I and v = 0 on I n[fi;fl].
61
Assume that for su?ciently large n
fi = ~fihn; fl = ~flhn; where 0 ? ~fi ? ~fl ? m2n?1:
Let X be the set of all linear combinations of the functions from K. The set X is dense in
L2(I;V).
To show that
limn
k!1
Z T
0
?u
nk(t)? ~unk(t);v(t)
?
V dt = 0 8v 2 L
2(I;V);
it is su?cient to show that
limn
k!1
Z T
0
?u
nk(t)? ~unk(t);v(t)
?
V dt = 0 8v 2 K;
since each of the functions from X is a linear combination of functions from K.
Fix a function v(t) from K and assume that nk is su?ciently large, so
fi = ~fihnk; fl = ~flhnk; where 0 ? ~fi ? ~fl ? m2nk?1:
Then 8v 2 K
Z T
0
?u
nk(t)? ~unk(t);v
?
V dt =
Z fl
fi
?u
nk(t)? ~unk(t);v
?
V dt
=
Z ~flhn
k
~fihnk
?u
nk(t)? ~unk(t);v
?
V dt:
62
Recall that
unk(t) = wnkj?1 +(wnkj ?wnkj?1)t?t
nk
j?1
hnk ; t 2
~Inkj = (tnkj?1;tnkj ];
and
~unk(t) = wnkj ; t 2 ~Inkj = (tnkj?1;tnkj ]:
This implies that
unk(t)? ~unk(t) = (wnkj ?wnkj?1)
"t?tn
kj?1
hnk ?1
#
= (wnkj ?wnkj?1)t?t
nk
j
hnk (since hnk = t
nk
j ?t
nk
j?1):
Then we have
Z tnk
j
tnkj?1
?
(wnkj ?wnkj?1)t?t
nk
j
hnk ;v
?
V
dt =
Z tnk
j
tnkj?1
?
wnkj ?wnkj?1;v
?
V
t?tnkj
hnk dt
=
?
wnkj ?wnkj?1;v
?
V
1
hnk
"(t?tn
kj )2
2
#tnkj
tnkj?1
:
Since hnk = tnkj ?tnkj?1,
Z tnk
j
tnkj?1
?
(wnkj ?wnkj?1)t?t
nk
j
hnk ;v
?
V
dt = ?wnkj?1 ?wnkj ;v?V hnk2 : (3.94)
63
From which it follows that
Z fl
fi
?
unk(t)? ~unk(t);v
?
V
dt
= hnk2
?
(wnk~fi ?wnk~fi+1)+(wnk~fi+1 ?wnk~fi+2)+???+(wnk~fl?1 ?wnk~fl );v
?
V
= hnk2
?
wnk~fi ?wnk~fl ;v
?
V
:
Recall (3.85), jjwnjjj1 ? c1, then
flfl
fl(wnk~fi ?wnk~fl ;v)
flfl
fl?jjvjjVjjwnk~fi ?wnk~fl jj?jjvjjV
?
jjwnk~fi jj+jjwnk~fl jj
?
? 2c1jjvjjV
Now as nk !1, hnk ! 0, therefore since v is flxed
limn
k!1
Z T
0
?
unk(t)? ~unk(t);v
?
V
dt = 0:
Similarly, we can show that this limit is zero when v(t) is a piecewise constant function of
t 2 I. Since the piecewise constant functions are dense in X, this proof is also valid for
every function v 2 L2(I;V). Therefore, ~unk * u in L2(I;V).
Similarly, using the same approach ~pnk * p in L2(I;L2(?)).
We now consider the question in which sense the functions u(t) and p(t) satisfy the given
system of partial difierential equations. We have by (3.57) and (3.58) the system for j =
1;??? ;m2nk?1
a(wnkj ;v) ? fi(znkj ;r?v) = ?(1?fl)fi < znkj ;v > 8v 2 V
64
Se 1h
nk
(znkj ?znkj?1;q) + fi 1h
nk
(r?wnkj ?r?wnkj?1;q) + k?(rznkj ;rq)
= (Qj;q) + (1?fl)fi 1h
nk
< wnkj ?wnkj?1;q > 8q 2 M
Deflne the abstract function Q(t) to be the constant Q for all t 2 [0;T] and let v(t) and
q(t) be arbitrary functions in L2(I;V) and L2(I;M) respectively. With Wnkj = w
nk
j ?w
nk
j?1
hnk
and Pnkj = p
nk
j ?p
nk
j?1
hnk , we get
a(wnkj ;v) ? fi(znkj ;r?v) = ?(1?fl)fi < znkj ;v >; (3.95)
Se(Znkj ;q) + fi(r?Wnkj ;q) + k?(rznkj ;rq)
= (Q;q) + (1?fl)fi < Wnkj ;q > : (3.96)
Recall that
~unk(t) = wnkj ; Unk(t) = Wnkj ; for t 2 ~Inj = (tnkj?1;tnkj ]; j = 2;??? ;m2nk?1;
~pnk(t) = znkj ; Pnk(t) = Znkj ; for t 2 ~Inj = (tnkj?1;tnkj ]; j = 2;??? ;m2nk?1:
Integrate from 0 to T to obtain
Z T
0
a
?
~unk(t);v(t)
?
dt ?
Z T
0
fi
?
~pnk(t);r?v(t)
?
dt = ?
Z T
0
(1?fl)fi < ~pnk(t);v(t) > dt;
Z T
0
Se
?
Pnk(t);q(t)
?
dt +
Z T
0
fi
?
r?Unk(t);q(t)
?
dt +
Z T
0
k
?
?
r~pnk(t);rq(t)
?
dt
=
Z T
0
?Q;q(t)?dt + Z T
0
(1?fl)fi < Unk(t);q(t) > dt:
65
Each of these integrals exists since v 2 L2(I;V), q 2 L2(I;L2(?)) (consequently, v 2
L2(I;H1(?)) and q 2 L2(I;H1(?))), ~pnk 2 L2(I;L2(?)), Unk 2 L2(I;H1(?)), Pnk 2
L2(I;H1(?)), and Q 2 L2(I;L2(?)).
The integral RT0 a(~unk;v)dt deflnes a bounded linear functional on L2(I;V) since a(:;:) is a
bounded bilinear form on V. For a flxed v 2 L2(I;V),
 Z T
0
a(~unk(t);v(t))dt
?2
? C2
 Z T
0
jj~unk(t)jjVjjv(t)jjV dt
?2
? C2
Z T
0
jj~unk(t)jj2V dt
Z T
0
jjv(t)jj2V dt
? C2jj~unkjj2L2(I;V)jjvjj2L2(I;V):
Thus for a flxed v, the integral RT0 a(~unk(t);v(t))dt ? Cjj~unkjjL2(I;V).
For nk !1, ~unk * u in L2(I;V), thus
Z T
0
a
?
~unk(t);v(t)
?
dt !
Z T
0
?
u(t);v(t)
?
dt; for nk !1:
Furthermore, for nk !1 we have that
~pnk * p in L2(I;L2(?));
hence, Z
T
0
fi
?
~pnk(t);r?v(t)
?
dt !
Z T
0
fi
?
p(t);r?v(t)
?
dt;
and Z
T
0
(1?fl) < ~pnk(t);v(t) > dt !
Z T
0
(1?fl) < p(t);v(t) > dt;
66
and Z
T
0
k
?
?
r~pnk(t);rq(t)
?
dt !
Z T
0
k
?
?
rp(t);rq(t)
?
dt;
as nk !1.
For nk !1, Pnk * pt in L2(I;L2(?)), thus
Z T
0
?
Pnk(t);r?v(t)
?
dt !
Z T
0
?
pt(t);r?v(t)
?
dt; as nk !1:
Furthermore, as nk !1, we have Unk * ut in L2(I;H1(?)), hence
Z T
0
fi
?
r?Unk(t);q(t)
?
dt !
Z T
0
fi
?
r?ut(t);q(t)
?
dt
and Z
T
0
(1?fl)fi < Unk(t);q(t) > dt !
Z T
0
(1?fl)fi < ut(t);q(t) > dt;
for nk !1.
Since v(t) and q(t) were arbitrary functions from L2(I;V) and L2(I;M), we have that
Z T
0
a
?
u(t);v(t)
?
dt?
Z T
0
fi
?
p(t);r?v(t)
?
dt
= ?
Z T
0
(1?fl) < p(t);v(t) > dt 8v 2 L2(I;V);
and
Z T
0
fi
?
r?ut(t);q(t)
?
dt+
Z T
0
Se
?
pt(t);q(t)
?
dt+
Z T
0
k
?
?
rp(t);rq(t)
?
dt
=
Z T
0
?
Q;q(t)
?
dt+
Z T
0
(1?fl)fi < ut(t);q(t) > dt 8q 2 L2(I;M):
67
Therefore, u(t) and p(t) satisfy the given system of partial difierential equations weakly and
have the following properties
u 2 L2(I;V); p 2 L2(I;L2(?));
u 2 AC(I;H1(?)); p 2 AC(I;H1(?));
ut 2 AC(I;H1(?)); pt 2 AC(I;H1(?));
u(0) = 0 in C(I;H1(?)); p(0) = 0 in C(I;H1(?));
Z T
0
a(u;v)dt?
Z T
0
fi(p;r?v)dt
= ?
Z T
0
(1?fl) < p;v > dt 8v 2 L2(I;V);
Z T
0
Se(pt;q)dt+
Z T
0
fi(r?ut;q)dt+
Z T
0
k
?(rp;rq)dt
=
Z T
0
(Q;q)dt+
Z T
0
(1?fl)fi < ut;q > dt 8q 2 L2(I;M):
In conclusion, we just proved existence of weak solutions u(t) and p(t) for problem (3.27){
(3.34).
68
Uniqueness of weak solutions
Let (~u; ~p) and (^u; ^p) be two solutions of problem (3.27){(3.34). Then u = ~u?^u and p = ~p?^p
are also solutions of this problem.
We have
a(~u;v) ? fi?~p;r?v? = (F;v)?(1?fl)fi < ~p;v > (3.97)
Se(~pt;q) + fi?(r? ~u)t;q? + k?(r~p;rq) = (Q;q)
+ < h1;q > +(1?fl)fi < ~ut;q > (3.98)
and
a(^u;v) ? fi?^p;r?v? = (F;v)?(1?fl)fi < ^p;v > (3.99)
Se(^pt;q) + fi?(r? ^u)t;q? + k?(r^p;rq) = (Q;q)
+ < h1;q > +(1?fl)fi < ^ut;q > (3.100)
Setting u = ~u? ^u and p = ~p? ^p, we subtract (3.99) from (3.97) and (3.100) from (3.98)
then integrate over I to get
Z T
0
a(u;v)dt ?
Z T
0
fi?p;r?v?dt = ?
Z T
0
(1?fl)fi < p;v > dt (3.101)
Z T
0
Se(pt;q)dt +
Z T
0
fi?(r?u)t;q?dt +
Z T
0
k
?(rp;rq)dt =
Z T
0
(1?fl)fi < ut;q > dt: (3.102)
69
Choose an arbitrary a 2 I and let
v(t) =
8
><
>:
ut(t) if 0 ? t ? a;
0 if a < t ? T;
and
q(t) =
8
><
>:
p(t) if 0 ? t ? a;
0 if a < t ? T:
Hence
Z a
0
a(u;ut)dt ?
Z a
0
fi?p;(r?u)t?dt = ?
Z a
0
(1?fl)fi < p;ut > dt (3.103)
Z a
0
Se(pt;p)dt +
Z a
0
fi?(r?u)t;p?dt +
Z a
0
k
?(rp;rp)dt =
Z a
0
(1?fl)fi < ut;p > dt:
(3.104)
Adding (3.103) and (3.104), we get
Z a
0
Se(pt;p)dt +
Z a
0
a(u;ut)dt +
Z a
0
k
?(rp;rp)dt = 0 (3.105)
Obviously, Z
a
0
k
?(rp;rp)dt =
Z a
0
k
?jjrpjj
2 ? 0; (3.106)
Z a
0
Se(pt;p)dt = 12Sejjp(a)jj2 ? 12Sejjp(0)jj2 = 12Sejjp(a)jj2 ? 0; (3.107)
and Z
a
0
a(u;ut)dt = 12[[u(a)]]2 ? 12[[u(0)]]2 = 12[[u(a)]]2 ? 0: (3.108)
70
Hence we conclude that
jju(a)jj = 0 and jjp(a)jj = 0:
And since a was arbitrary then
jju(t)jj = 0 in I;
jjp(t)jj = 0 in I.
Therefore, ~u(t) = ^u(t) and ~p(t) = ^p(t) and we conclude that the system has a unique
solution (u(t);p(t)).
71
3.3.2 Energy norm estimate for
homogeneous initial and boundary conditions
Given the quasi-static poroelasticity system of partial difierential equations with ho-
mogeneous initial and boundary conditions then the weak formulation yields
a(u;v) ? fi(p;r?v) = ?(1?fl)fi < p;v >; 8v 2 V;
Se(pt;q) + fi((r?u)t;q) + k?(rp;rq) = (Q(t);q)
+ (1?fl)fi < ut;q >; 8q 2 M;
for almost every t 2 I.
Let v = ut and q = p and add the two equations to obtain
a(u;ut) + Se(pt;p) + k?(rp;rp) = (Q;p);
and since k?jjrpjj2 ? 0
1
2
d
dt[[u]]
2 + Se
2
d
dtjjpjj
2 ?jjQjj jjpjj; (3.109)
That is,
d
dt
?[[u]]2 +jjpjj2? ? 2
min(1;Se)jjQjj jjpjj;
? jjQjj
2
(min(1;Se))2 +jjpjj
2; (using 2ab ? a2 +b2)
? jjQjj
2
(min(1;Se))2 +jjpjj
2 +jjujj2
1: (3.110)
72
Lemma 1 (Gronwall?s Inequality (see [5])): Let ? be a nonnegative, absolutely continuous
function on [0;T], which satisfles for a.e. t the difierential inequality
?0(t) ? `(t)?(t)+?(t);
where `(t) and ?(t) are nonnegative functions on [0;T]. Then
?(t) ? e
Rt
0 `(s)ds
?
?(0)+
Z t
0
?(s)ds
?
;
for all 0 ? t ? T.
Therefore applying Gronwall?s inequality to (3.109), we obtain
[[u]]2 +jjp(t)jj2 ? e
Rt
0 ds
?
[[u]]2 +jjp(0)jj2 +
Z t
0
jjQjj2
(min(1;Se))2ds
?
: (3.111)
Using the fact that u(0) = p(0) = 0 (from (3.70)), we get
[[u]]2 +jjp(t)jj2 ? e
Rt
0 ds
?Z t
0
jjQjj2
(min(1;Se))2ds
?
:
Then
jjp(t)jj2 ? e
Rt
0 ds
?Z t
0
jjQjj2
(min(1;Se))2ds
?
;
and
2Gjju(t)jj21 ? e
Rt
0 ds
?Z t
0
jjQjj2
(min(1;Se))2ds
?
:
Therefore,
max
0?t?T
jjp(t)jj2 ? e
T
(min(1;Se))2jjQjj
2
L2(0;T;L2(?));
73
and
max
0?t?T
jju(t)jj21 ? e
T
2G(min(1;Se))2jjQjj
2
L2(0;T;L2(?));
Hence
jjp(t)jj ?
peT
min(1;Se)jjQjjL2(0;T;L2(?)); (3.112)
and
jju(t)jj1 ?
peT
p2Gmin(1;Se)jjQjjL2(0;T;L2(?)): (3.113)
74
The right hand side of the elasticity equation
Let ?u be the solution of the stationary elasticity problem
?Gr?(r?u+(r?u)T)?G 2?(1?2?)r(r? ?u) = F:
The weak formulation of this problem is: flnd ?u 2 V such that
a(?u;v) = (F;v)+
Z
?t
?
G(r?u+(r?u)T)+G 2?1?2?(r? ?u)
?
? ^nv 8v 2 V; (3.114)
where a(?u;v) = R?[G(r?u + (r?u)T) : rv + G 2?1?2?(r? ?u)(r?v)].
And let ~u be the solution of the quasi-static poroelasticity system with F = 0:
?Gr?(r~u+(r~u)T)?G 2?(1?2?)r(r? ~u)+firp = 0;
@
@t(Se p+fir? ~u)?
k
??p = Q:
The weak formulation of this problem is: flnd ~u 2 V and p 2 M such that
a(~u;v) + fi(rp;v) =
Z
?t
?
G(r~u+(r~u)T)+G 2?1?2?(r? ~u)
?
? ^nv 8v 2 V;(3.115)
Se(pt;q)+fi
?
(r? ~u)t;q
?
+ k?(rp;rq) = (Q;q)+
Z
?f
k
?rp? ^nq 8q 2 M;(3.116)
75
The above equation holds for almost every t 2 I. We will show that u = ?u+ ~u satisfles the
poroelasticity system with F 6= 0
?Gr?(r~u+(r~u)T)?G 2?(1?2?)r(r? ~u)+firp = F;
@
@t(Se p+fir? ~u)?
k
??p = Q:
Adding equations (3.114), (3.115), and (fi(r??u)t;q) = 0 (which is zero since ?u is the solution
to a stationary problem) to (3.116), we get
a(?u;v) + a(~u;v) + fi(rp;v) = (F;v)+
Z
?t
?
G(r?u+(r?u)T)+G 2?1?2?(r? ?u)
?
? ^nv
+
Z
?t
?
G(r~u+(r~u)T)+G 2?1?2?(r? ~u)
?
? ^nv;
Se(pt;q)+fi
?
(r? ?u)t;q
?
+fi
?
(r? ~u)t;q
?
+ k?(rp;rq) = (Q;q)+
Z
?f
k
?rp? ^nq:
That is,
a(?u+ ~u;v)+fi(rp;v) = R?t ?Gr(?u+ ~u)+(r(?u+ ~u))T
+G 2?1?2?r?(?u+ ~u)
i
? ^nv +(F;v)
Se(pt;q)+fi
?
(r?(?u+ ~u))t;q
?
+ k?(rp;rq) = (Q;q)+R?f k?rp? ^nq:
Therefore,
a(u;v) + fi(rp;v) = (F;v)+
Z
?t
?
G(ru+(ru)T)+G 2?1?2?(r?u)
?
? ^nv 8v 2 V
Se(pt;q)+fi
?
(r?u)t;q
?
+ k?(rp;rq) = (Q;q)+
Z
?f
k
?rp? ^nq; 8q 2 M:
Again the above equation holds for almost every t 2 I. Hence we got the weak formulation
of the poroelasticity system with F 6= 0.
76
3.3.3 Existence and uniqueness of weak solutions for
nonhomogeneous initial conditions
Consider problem II (3.35){(3.42) (with nonhomogeneous initial conditions).
From equations (3.57) and (3.58), we have for j = 1;??? ;m
a(wj;v) ? fi(zj;r?v) = ?(1?fl)fi < zj;v > 8v 2 V; (3.117)
Se(zj?zj?1;q) + fi(r?wj?r?wj?1;q) + hk?(rzj;rq) =
(1?fl)fi < wj ?wj?1;q > 8q 2 M; (3.118)
with initial conditions
Se z0 +fir?w0 = v0; on ?;
(1?fl)fiw0 ? ^n = v1; on ?tf:
Let
~zj = zj ? v0se
and
~wj = wj;
then the system becomes
a(~wj;v) ? fi
?
~zj + v0se;r?v
?
= ?(1?fl)fi < ~zj + v0se;v >;
77
Se(~zj ? ~zj?1;q)+fi(r? ~wj ?r? ~wj?1;q)+hk?
?
r~zj +rv0Se;rq
?
=
(1?fl)fi < ~wj ? ~wj?1;q > :
That is,
a(~wj;v) ? fi(~zj;r?v) = fi
?v0
Se;r?v
?
?(1?fl)fi < v0Se;v > ?(1?fl)fi < ~zj;v >; (3.119)
and
Se(~zj ? ~zj?1;q) + fi(r? ~wj ?r? ~wj?1;q) + hk?(r~zj;rq) =
?hk?
?
rv0Se;rq
?
+(1?fl)fi < ~wj ? ~wj?1;q > : (3.120)
With initial conditions
Se
?
~z0 + v0Se
?
+fir? ~w0 = v0;
(1?fl)fi~w0 ? ^n = v1;
which implies that
Se ~z0 +fir? ~w0 = 0; (3.121)
(1?fl)fi~w0 ? ^n = v1: (3.122)
78
Therefore, we transformed the given problem (3.35){(3.42) to an equivalent one with ho-
mogeneous initial condition (3.121). Furthermore, each of these integrals ?v0Se;r?v?, <
v0
Se;v >, and
?rv
0Se;rq
? exists since v 2 V, q 2 M (consequently, v 2 L2(I;H1(?)),
q 2 L2(I;H1(?))), and v0 2 L2(I;H1(?)), (fi;fl;Se are positive constants).
At time j ?1, j = 2;??? ;m, (3.119) and (3.120) become
a(~wj?1;v)?fi(~zj?1;r?v) =
fi
?v0
Se;r?v
?
?(1?fl)fi < v0Se;v > ?(1?fl)fi < ~zj?1;v >; (3.123)
Se(~zj?1?~zj?2;q) + fi(r?~wj?1?r?~wj?2;q) + hk?(r~zj?1;rq) =
?hk?
?
rv0Se;rq
?
+(1?fl)fi < ~wj?1 ? ~wj?2;q > : (3.124)
Subtracting (3.123) from (3.119) and (3.124) from (3.120), we obtain
a(~wj ? ~wj?1;v) ? fi(~zj ? ~zj?1;r?v) = ?(1?fl)fi < ~zj ? ~zj?1;v >; (3.125)
Se
?
~zj?~zj?1?(~zj?1?~zj?2);q
?
+ fi
?
r? ~wj?r? ~wj?1?(r? ~wj?1?r? ~wj?2);q
?
+ hk?(r~zj ?r~zj?1;rq) = (1?fl)fi < ~wj ? ~wj?1 ?(~wj?1 ? ~wj?2);q > : (3.126)
Adding (3.125) and (3.126) with v = ~wj ? ~wj?1 and q = ~zj ? ~zj?1, we get
79
a(~wj?~wj?1; ~wj?~wj?1)?fi(r?~wj?1?r?~wj?2;~zj?~zj?1)+Se(~zj?~zj?1;~zj?~zj?1)
?Se(~zj?1 ? ~zj?2;~zj ? ~zj?1)+hk?(r~zj ?r~zj?1;r~zj ?r~zj?1) =
?(1?fl)fi < ~wj?1 ? ~wj?2;~zj ? ~zj?1 > : (3.127)
Equation (3.125) with v = ~wj?1 ? ~wj?2, j = 2;??? ;m, is
a(~wj ? ~wj?1; ~wj?1 ? ~wj?2)?fi(~zj ? ~zj?1;r? ~wj?1 ?r? ~wj?2) =
?(1?fl)fi < ~zj ? ~zj?1; ~wj?1 ? ~wj?2 > : (3.128)
Subtracting (3.128) from (3.127), we get
a(~wj ? ~wj?1; ~wj ? ~wj?1)+Se(~zj ? ~zj?1;~zj ? ~zj?1)+hk?(r~zj ?r~zj?1;r~zj ?r~zj?1) =
a(~wj ? ~wj?1; ~wj?1 ? ~wj?2)+Se(~zj ? ~zj?1;~zj?1 ? ~zj?2);
which implies that,
2Gjj~wj?~wj?1jj21+Sejj~zj?~zj?1jj2
? C1jj~wj ? ~wj?1jj1jj~wj?1 ? ~wj?2jj1 +Sejj~zj ? ~zj?1jj jj~zj?1 ? ~zj?2jj;
80
where C1 =
?
2G; 6G?1?2?
?
is the continuity constant of the bilinear form a(?;?) (3.11). Then,
jj~wj?~wj?1jj21+jj~zj?~zj?1jj2
? max(C1;Se)min(2G;Se)
?
jj~wj ? ~wj?1jj21 +jj~zj ? ~zj?1jj2
?1
2
?
jj~wj?1 ? ~wj?2jj21 + jj~zj?1 ? ~zj?2jj2
?1
2:
Squaring both sides and simplifying, we obtain
jj~wj?~wj?1jj21+jj~zj?~zj?1jj2
?
 max(C
1;Se)
min(2G;Se)
?2h
jj~wj?1 ? ~wj?2jj21 + jj~zj?1 ? ~zj?2jj2
i
: (3.129)
Setting j = 1 in (3.125) and letting v = ~w1 ? ~w0, we get
a(~w1? ~w0; ~w1? ~w0) ? fi(~z1?~z0;r? ~w1?r? ~w0) = ?(1?fl)fi < ~z1?~z0:~w1? ~w0 >; (3.130)
Setting j = 1 in (3.120) and letting q = ~z1 ? ~z0, we have
Se(~z1?~z0;~z1?~z0)+fi(r?~w1?r?~w0;~z1?~z0)+hk?(r~z1;r~z1?r~z0) =
?hk?
?
rv0Se;r~z1 ?r~z0
?
+(1?fl)fi < ~w1 ? ~w0;~z1 ? ~z0 > : (3.131)
Using (r~z1;r~z1 ?r~z0) = (r~z1 ?r~z0 + r~z0;r~z1 ?r~z0) = (r~z1 ?r~z0;r~z1 ?r~z0) +
(r~z0;r~z1 ?r~z0); then adding (3.130) and (3.131), we obtain
81
2Gjj~w1 ? ~w0jj21 +Sejj~z1 ? ~z0jj2 +hk?jjr~z1 ?r~z0jj2
? hk?jjrv0jjSe jjr~z1 ?r~z0jj+hk?jjr~z0jj jjr~z1 ?r~z0jj
? hk?
?jjrv
0jj
Se +jjr~z0jj
?
jjr~z1 ?r~z0jj:
Sincezj solvesthehomogeneousinitialconditionsporoelasticityproblem, thenfromproblem
I (z0 = w0 = 0) we have
~z0 = zo ? v0Se = ?v0Se; and ~w0 = w0 = 0: (3.132)
Therefore,
2Gjj~w1 ? ~w0jj21 +Sejj~z1 ? ~z0jj2 +hk?jjr~z1 ?r~z0jj2
? 12?
?
hk?
 2jjrv
0jj
Se
??2
+ ?2jjr~z1 ?r~z0jj2:
Choose ? small enough such that hk? ? ?2 > 0, then (hk? ? ?2)jjr~z1 ?r~z0jj2 ? 0. Hence
jj~w1 ? ~w0jj21 +jj~z1 ? ~z0jj2 ? h2 2? min(2G;Se)
?k
?
jjrv0jj
Se
?2
: (3.133)
Using (3.129) and (3.133), we obtain (j = 1;??? ;m)
jj~wj ? ~wj?1jj21 +jj~zj ? ~zj?1jj2 ? h2
 max(C
1;Se)
min(2G;Se)
?2 2
? min(2G;Se)
?k
?
jjrv0jj
Se
?2
: (3.134)
82
Let
C =
s
2
? min(2G;Se)
max(C1;Se)
min(2G;Se)
k
?
jjrv0jj
Se ;
therefore
jj~wj ? ~wj?1jj21 +jj~zj ? ~zj?1jj2 ? h2C2: (3.135)
Recall the norm:
jjj(~wj;~zj)?(~wj?1;~zj?1)jjj =
?
jj~wj ? ~wj?1jj21 + jj~zj ? ~zj?1jj2
?1
2;
then,
jjj(~wj;~zj)?(~w0;~z0)jjj = jjj(~wj;~zj)?(~wj?1;~zj?1)+ ??? +(~w1;~z1)?(~w0;~z0)jjj
? jjj(~wj;~zj)?(~wj?1;~zj?1)jjj+ ??? +jjj(~w1;~z1)?(~w0;~z0)jjj;
? jhC
From which it follows that
jj~wj ? ~w0jj21 +jj~zj ? ~z0jj2 ? j2h2C2 (3.136)
Since h = Tm, then
jj~wj ? ~w0jj21 ? T2C2
and
jj~zj ? ~z0jj2 ? T2C2;
that is,
jj~wjjj1 ? TC +jj~w0jj1
83
and
jj~zjjj? TC +jj~z0jj:
Hence using the fact that ~z0 = ?v0Se and ~w0 = 0 (3.132), we get
jj~wjjj1 ? TC and jj~zjjj? TC + jjv0jjSe : (3.137)
Using (3.135) with ~Wj = ~wj?~wj?1h and ~Zj = ~zj?~zj?1h , we obtain
jj ~Wjjj21 +jj~Zjjj2 ? C2:
Hence
jj ~Wjjj1 ? C and jj~Zjjj? C: (3.138)
The estimates obtained in (3.137) and (3.138) are independent of h, thus remain valid
for an arbitrary mesh dn. That is, for every positive integer n and j = 1;??? ;m2n?1
jj~wnjjj1 ? TC jj~znj jj ? TC + jjv0jjSe ; (3.139)
jj ~Wnj jj1 ? C; and jj~Znj jj ? C: (3.140)
Recall that ~wj = wj, ~zj = zj ? v0Se, and the constant
C =
s
2
? min(2G;Se)
max(C1;Se)
min(2G;Se)
k
?
jjrv0jj
Se ;
84
thus the estimates (3.139) and (3.140) become
jjwnjjj1 ? 2Tkp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se ; (3.141)
jjznj jj ? 2Tkp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se +2
jjv0jj
Se ; (3.142)
jjWnj jj1 ? 2kp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se ; (3.143)
and
jjZnj jj ? 2kp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se : (3.144)
These basic estimates can then be used to show existence and uniqueness of weak
solutions as done for problem I (3.27){(3.34).
85
3.3.4 Energy norm estimate for nonhomogeneous initial conditions
For the given quasi-static poroelasticity system with nonhomogeneous initial condi-
tions, we transformed the problem to an equivalent one with nonhomogeneous right hand
side and homogeneous initial condition (3.121). Then the weak formulation of the trans-
formed problem is: flnd u 2 V and p 2 M such that
a(u;v) ? fi
?
p+ v0Se;r?v
?
= ?(1?fl)fi < p+ v0Se;v >; 8 v 2 V; (3.145)
Se(pt;q) + fi(r?ut;q) + k?(rp+rv0Se;rq) = (1?fl)fi < ut;q >; 8 q 2 M; (3.146)
for almost every t 2 I.
Letting v = ut and q = p+ v0Se and adding (3.145) and (3.146), we obtain
a(u;ut) + Se
?
pt + v0tSe;p+ v0Se
?
+ k?
?
rp+rv0Se;rp+rv0Se
?
= Se
?v0
t
Se;p+
v0
Se
?
;
therefore
1
2
d
dt[[u]]
2 + Se
2
d
dt
flfl
fl
flfl
flp+ v0Se
flfl
fl
flfl
fl
2 + k
?
flfl
fl
flfl
flrp+rv0Se
flfl
fl
flfl
fl
2 ?jjv
0tjj
flfl
fl
flfl
flp+ v0Se
flfl
fl
flfl
fl:
That is,
d
dt[[u]]
2 + d
dt
flfl
fl
flfl
flp+ v0Se
flfl
fl
flfl
fl
2 ? 2 jjv0tjj
min(1;Se)
flfl
fl
flfl
flp+ v0Se
flfl
fl
flfl
fl:
Then
d
dt
 
[[u]]2 +
flfl
fl
flfl
flp+ v0Se
flfl
fl
flfl
fl
2?? jjv0tjj
min(1;Se)
?2
+
flfl
fl
flfl
flp+ v0Se
flfl
fl
flfl
fl
2 +[[u]]2:
86
Integrating from 0 to t, where t 2 [0;T], using Gronwall?s inequality, and using the initial
conditions for the transformed problem u(0) = 0 and p(0) = ?v0Se, we obtain
[[u]]2 +
flfl
fl
flfl
flp+ v0Se
flfl
fl
flfl
fl
2 ? eT
(min(1;Se))2jjv0tjj
2
L2(0;T;H1(?)):
Hence
jjp(t)jj ?
peT
min(1;Se)jjv0tjjL2(0;T;H1(?)) +
jjv0jjL2(0;T;H1(?))
Se ; (3.147)
and
jju(t)jj1 ?
peT
p2Gmin(1;Se)jjv0tjjL2(0;T;H1(?)): (3.148)
87
3.3.5 Existence and uniqueness of weak solutions for
nonhomogeneous boundary conditions
For problem III (3.43){(3.50) (with nonhomogeneous boundary conditions, i.e., h1 6= 0),
we assume that h1(t) 2 C0;1(0;T;L2(?f)) and use the Rothe method to approximate the
solution as done for problem I (3.27){(3.34). Then for j = 2;??? ;m we have
a(wj;v) ? fi(zj;r?v) = ?(1?fl)fi < zj;v > 8v 2 V; (3.149)
and
Se(zj ?zj?1;q) + fi(r?wj ?r?wj?1;q) + hk?(rzj;rq) =
h < h1j;q > + (1?fl)fi < wj ?wj?1;q > 8q 2 M: (3.150)
At time (j ?1), j = 2;??? ;m
a(wj?1;v) ? fi(zj?1;r?v) = ?(1?fl)fi < zj?1;v > (3.151)
Subtracting (3.151) from (3.149), we get
a(wj ?wj?1;v) ? fi(zj ?zj?1;r?v) = ?(1?fl)fi < zj ?zj?1;v >; (3.152)
and
Se(zj ?zj?1;q) + fi(r?wj ?r?wj?1;q) + hk?(rzj;rq) =
h < h1j;q > + (1?fl)fi < wj ?wj?1;q > : (3.153)
88
Setting j = 1, assuming that (3.149) and (3.150) hold at j = 0, and using the previous
equations (3.151){(3.153), we obtain
a(w1?w0;w1?w0) ? fi(z1?z0;r?w1?r?w0) = ?(1?fl)fi < z1?z0;w1?w0 >; (3.154)
Se(z1 ?z0;z1 ?z0) + fi(r?w1 ?r?w0;z1 ?z0) + hk?(rz1;rz1 ?rz0) =
h < h11;z1 ?z0 > + (1?fl)fi < z1 ?z0;w1 ?w0 >; (3.155)
Se(z1 ?z0;z0) + fi(r?w1 ?r?w0;z0) + hk?(rz1;rz0) =
h < h11;z0 > + (1?fl)fi < z0;w1 ?w0 >; (3.156)
a(w0;w1 ?w0) ? fi(z0;r?w1 ?r?w0) = ?(1?fl)fi < z0;w1 ?w0 > : (3.157)
Adding (3.154){(3.157) and simplifying, we get
Sejjz1?z0jj2 +[[w1?w0]]2 +hk?jjrz1jj2 ? [[w1?w0]][[w0]]+Sejjz1?z0jjjjz0jj+h < h11;z1 >
Using the homogeneous initial conditions (w0 = z0 = 0 (3.70)), we obtain
Sejjz1 ?z0jj2 + [[w1 ?w0]]2 + hk?jjrz1jj2 ? hjjh11jjL2(?f)jjz1jjL2(?f): (3.158)
Given q 2 M and applying theorem B.1, we have
jjqjjL2(?f) ? cjjqjjM; for some c > 0;and 8q 2 M: (3.159)
89
We have by Poincare?s inequality (Appendix A)
(rq;rq) = jjrqjj2 ? ?jjqjj2M; for some ? > 0; 8q 2 M: (3.160)
Therefore, using (3.159) and (3.160) with q = z1, inequality (3.158) becomes
Sejjz1 ?z0jj2 +[[w1 ?w0]]2 +hk??jjz1jj2M ? hcjjh11jjL2(?f)jjz1jjM
? hc?2jjh11jj2L2(?f) + hc2?jjz1jj2M
? h?jjh11jj2L2(?f) +h?1jjz1jj2M: (3.161)
Where ? and ?1 are positive constants that can be found with ?1 < k??, then
hk??jjz1jj2M ? h?1jjz1jj2M ? 0; (3.162)
and (3.161) yields
Sejjz1 ?z0jj2 + [[w1 ?w0]]2 ? h?jjh11jj2L2(?f): (3.163)
And (w0 = z0 = 0 (3.70))
Sejjz1jj2 + [[w1]]2 ? h?jjh11jj2L2(?f): (3.164)
Setting now j = 2 we have
a(w2 ?w1;w2 ?w1) ? fi(z2 ?z1;r?w2 ?r?w1) = ?(1?fl)fi < z2 ?z1;w2 ?w1 >;
90
Se(z2 ?z1;z2 ?z1) + fi(r?w2 ?r?w1;z2 ?z1) + hk?(rz2;rz2 ?rz1) =
h < h11;z2 ?z1 > + (1?fl)fi < w2 ?w1;z2 ?z1 >;
a(w1;w2 ?w1) ? fi(z1;r?w2 ?r?w1) = ?(1?fl)fi < z1;w2 ?w1 >;
Se(z2 ?z1;z1) + fi(r?w2 ?r?w1;z1) + hk?(rz2;rz1) =
h < h11;z1 > + (1?fl)fi < w2 ?w1;z1 >;
a(w2;w2 ?w1) ? fi(z2;r?w2 ?r?w1) = ?(1?fl)fi < z2;w2 ?w1 >;
and
Se(z2 ?z1;z2) + fi(r?w2 ?r?w1;z2) + hk?(rz2;rz2) =
h < h11;z2 > + (1?fl)fi < w2 ?w1;z2 > :
If we add the previous six equations and use (3.159) and (3.160), we get
[[w2 ?w1]]2 + [[w2]]2 + Sejjz2 ?z1jj2 + Sejjz2jj2 + 2hk?jjrz2jj2M
? [[w1]]2 + Sejjz1jj2 + 2h?jjh12jj2 + 2h?1jjz2jj2M:
Note that we are again using here the symmetry of the bilinear form a(:;:) (3.15). Since
?1 < k??, then (by (3.162))
hk??jjz2jj2M ? h?1jjz2jj2M ? 0;
91
therefore,
[[w2 ?w1]]2 +[[w2]]2 +Sejjz2 ?z1jj2 +Sejjz2jj2 ? [[w1]]2 + Sejjz1jj2 + 2h?jjh12jj2:
The flrst and third terms in this inequality are positive, so we can write
[[w2]]2 + Sejjz2jj2 ? [[w1]]2 + Sejjz1jj2 + 2h?jjh12jj2: (3.165)
Repeating the same process, we get
[[wj?wj?1]]2+[[wj]]2+Sejjzj?zj?1jj2+Sejjzjjj2
? [[wj?1]]2 +Sejjzj?1jj2 +2h?jjh1jjj2; (3.166)
or
[[wj]]2 + Sejjzjjj2 ? [[wj?1]]2 +Sejjzj?1jj2 +2h?jjh1jjj2: (3.167)
That is we obtained
[[w1]]2 + Sejjz1jj2 ? h?jjh11jj2
[[w2]]2 + Sejjz2jj2 ? [[w1]]2 +Sejjz1jj2 +2h?jjh12jj2
...
[[wj]]2 + Sejjzjjj2 ? [[wj?1]]2 +Sejjzj?1jj2 +2h?jjh1jjj2
Adding these to obtain
[[wj]]2 + Sejjzjjj2 ? h?jjh11jj2 +2h?
h
jjh12jj2 +???+jjh1jjj2
i
? 2h?
h
jjh11jj2 +jjh12jj2 +???+jjh1jjj2
i
:
92
Since h1(t) 2 C0;1(0;T;L2(?f)), then there exists a constant d such that jjh1(t+h)?h1(t)h jj? d
for all t, t + h 2 I, (see [10]). Then, jjh1(t)jj is a continuous function on I and so jjh1(t)jj
attains a maximum on I, say jj ~Hjj i.e.,
max
t2I
jjh1(t)jj = jj ~Hjj:
Consequently,
[[wj]]2 + Sejjzjjj2 ? 2h?jjj ~Hjj2
? 2?Tjj ~Hjj2 (since h = Tm);
and
jjwjjj1 ? jj ~Hjj
r
?T
G and jjzjjj ? jj
~Hjj
r
2?T
Se : (3.168)
We have h1(t) 2 C0;1(0;T;L2(?f)) and assume that h1(0) = 0, then there exists a constant
C > 0 such that jjh1(?2)?h1(?1)jj ? C2j?2 ??1j2:
Then, jjh1(h)jj2 = jjh1(h)?h1(0)jj ? C2h2:
Inequalities (3.164){(3.167) with jjh1jjj2 ? C2h2 yield
[[w1]]2 + Sejjz1jj2 ? h?jjh11jj2 ? h3?C2
[[w2]]2 + Sejjz2jj2 ? [[w1]]2 +Sejjz1jj2 +2h3?C2
...
[[wj?1]]2 +Sejjzj?1jj2 ? [[wj?2]]2 +Sejjzj?2jj2 +2h3?C2
[[wj ?wj?1]]2 +Sejjzj ?zj?1jj2 ? [[wj?1]]2 +Sejjzj?1jj2 +2h3?C2
93
Adding these inequalities and simplifying, we get
[[wj ?wj?1]]2 +Sejjzj ?zj?1jj2 ? h3?C2 + 2(j ?1)h3?C2
Using Wj = wj?wj?1h and Zj = zj?zj?1h , the previous inequality becomes
[[Wj]]2 +SejjZjjj2 ? h?C2 + 2(j ?1)h?C2
? 2jh?C2
? 2T?C2 (since h = Tm);
therefore,
jjWjjj1 ? C
r
?T
G and jjZjjj ? C
r
2?T
Se : (3.169)
The estimates obtained in (3.168) and (3.169) are independent of h, thus remain valid for
an arbitrary mesh dn. That is, for every positive integer n and j = 1;??? ;m2n?1
jjwnjjj1 ? jj ~Hjj
r
?T
G ; jjz
n
j jj ? jj ~Hjj
r
2?T
Se ; (3.170)
jjWnj jj1 ? C
r
?T
G ; and jjZ
n
j jj ? C
r
2?T
Se : (3.171)
These basic estimates can then be used to show existence and uniqueness of weak solutions
as done for problem I (3.27){(3.34).
94
3.3.6 Energy norm estimate for
nonhomogeneous boundary condition
To flnd the energy norm for the poroelasticity problem with homogeneous initial con-
dition and nonhomogeneous boundary condition, we follow the same steps done in Section
3.3.2 for homogeneous initial and boundary condition. Then we get
a(u;ut) + Se(pt;p) + k?(rp;rp) = < h1;p >;
that is,
1
2
d
dt[[u]]
2 + Se
2
d
dtjjpjj
2 + k
?jjrpjj
2 ? ?1jjh1jj2
L2(?tf) +?2jjpjj
2
M (3.172)
with ?2 < k?c. By Poincare?s inequality, jjrpjj2 ? cjjpjj2, we have
k
?cjjpjj
2
M ??2jjpjj
2
M ? 0;
and so,
d
dt[[u]]
2 +Se d
dtjjpjj
2 ? 2?1jjh1jj2
L2(?tf):
Using the norm: jjj(u;p)jjj =
?
[[u]]2 +Sejjpjj2
?1
2, we get
d
dt
?
jjj(u;p)jjj2
?
? 2?1jjh1jj2L2(?tf):
Integrating from 0 to t, where t 2 [0;T], and using the homogeneous initial conditions
(u(0) = p(0) = 0), we obtain
jjj(u(t);p(t))jjj2 ? 2?1Tjjh1jj2L2(0;T;L2(?tf)):
Therefore,
[[u]]2 +Sejjpjj2 ? 2?1Tjjh1jj2L2(0;T;L2(?tf));
95
from which
jjp(t)jj ?
r
2?1T
Se jjh1jjL2(0;T;L2(?tf)); (3.173)
and
jju(t)jj1 ?
r
?1T
G jjh1jjL2(0;T;L2(?tf)): (3.174)
We can now obtain the energy norm estimates for the fully coupled poroelasticity system.
We obtained the energy norm estimates for homogeneous boundary and initial conditions
(3.112) and (3.113)
jju(t)jj1 ?
peT
p2Gmin(1;Se)jjQjjL2(0;T;L2(?)); jjp(t)jj?
peT
min(1;Se)jjQjjL2(0;T;L2(?));(3.175)
the energy norm estimates for nonhomogeneous initial conditions (3.147) and (3.148)
jju(t)jj1 ?
s
jjv0jj2
2GSe +
T
2?G
 k
?
jjrv0jj
Se
?2
;
jjp(t)jj ?
s
jjv0jj2
Se2 +
T
?Se
 k
?
jjrv0jj
Se
?2
; (3.176)
andtheenergynormestimatesfornonhomogeneousboundaryconditions(3.173)and(3.174)
jju(t)jj1 ?
r
?1T
G jjh1jjL2(0;T;L2(?tf)); jjp(t)jj?
r
2?1T
Se jjh1jjL2(0;T;L2(?tf)): (3.177)
Denote by
C1 =
peT
p2Gmin(1;Se)jjQjjL2(0;T;L2(?));
C2 =
peT
min(1;Se)jjQjjL2(0;T;L2(?));
96
C3 =
s
jjv0jj2
2GSe +
T
2?G
 k
?
jjrv0jj
Se
?2
;
C4 =
s
jjv0jj2
Se2 +
T
?Se
 k
?
jjrv0jj
Se
?2
;
C5 =
r
?1T
G jjh1jjL2(0;T;L2(?tf));
and
C6 =
r
2?1T
Se jjh1jjL2(0;T;L2(?tf)):
Hence the energy norm estimates for the fully coupled quasi-static poroelasticity problem
is
jju(t)jj1 ? C1 +C3 +C5; and jjp(t)jj? C2 +C4 +C6: (3.178)
Energy norm estimates for the discrete poroelasticity problem
Let V h ? V and Mh ? M be flnite dimensional spaces. The weak formulation for the
discrete problem with homogeneous initial and boundary conditions is: flnd uh 2 V h and
ph 2 Mh such that
a(uh;vh)?fi(ph;r?vh) = ?(1?fl)fi < ph;vh >; 8vh 2 V h;
Se(pht ;qh)+fi(r?uht ;qh)+ k?(rph;rqh) =
(Q;qh)+(1?fl)fi < uht ;qh >; 8qh 2 Mh:
97
Letting vh = uht and qh = ph and adding the two previous equation, we get
a(uh;uht ) + Se(pht ;ph) + k?(rph;rph) = (Q;ph);
that is,
1
2
d
dt[[u
h]]2 + Se
2
d
dtjjp
hjj2 ?jjQjj jjphjj:
From which it follows that
d
dt
?
[[uh]]2 +jjphjj2
?
? 2jjQjjmin(1;Se)jjphjj;
? jjQjj
2
(min(1;Se))2 +jjp
hjj2 +[[u]]2:
Applying Gronwall?s inequality and the fact that uh(0) = ph(0) = 0 (from (3.70)), we obtain
jjph(t)jj ?
peT
min(1;Se)jjQjjL2(0;T;L2(?)); (3.179)
jjuh(t)jj1 ?
peT
p2Gmin(1;Se)jjQjjL2(0;T;L2(?)): (3.180)
The estimates (3.179) and (3.179) are the same as (3.112) and (3.113) obtained for the semi-
discrete problem with homogeneous initial and boundary conditions. Similarly, we can get
the same energy estimates for the nonhomogeneous initial conditions and for the nonhomo-
geneous boundary conditions problems as done in sections 3.3.4 and 3.3.6. Therefore, we
can obtain
jjuh(t)jj1 ? C1 +C3 +C5; and jjph(t)jj? C2 +C4 +C6; (3.181)
with Ci;i = 1;??? ;6 are as deflned above.
98
3.4 Error estimates
In this section, we will obtain error estimates for the semi-discrete and for the fully
discrete poroelasticity problem (3.19){(3.26).
3.4.1 Error estimates for the semi-discrete problem
We will derive estimates of the difierence between the solutions u(t) and p(t) and the
approximations un(t) and pn(t) by the method of discretization in time.
Since the poroelasticity problem (3.19){(3.26) is linear, its error estimates can be writ-
ten as the sum of error estimates for homogeneous initial conditions and error estimates for
nonhomogeneous initial conditions.
Error estimates for homogeneous initial conditions
Consider the poroelasticity problem (3.19){(3.26) with homogeneous initial conditions. Re-
call the deflnitions of the functions un(t), pn(t), ~un(t), ~pn(t), Un(t), and Pn(t):
un(t) = wnj?1 + (t?tnj?1)Wnj ; t 2 Inj = [tnj?1;tnj];
pn(t) = znj?1 + (t?tnj?1)Znj ; t 2 Inj = [tnj?1;tnj];
~un(0) = wn1
~un(t) = wnj for t 2 ~Inj = (tnj?1;tnj]; j = 1;??? ;m2n?1;
99
~pn(0) = zn1
~pn(t) = znj for t 2 ~Inj = (tnj?1;tnj]; j = 1;??? ;m2n?1;
Un(0) = Wn1
Un(t) = Wnj for t 2 ~Inj := (tnj?1;tnj]; j = 1;??? ;m2n?1;
and,
Pn(0) = Zn1
Pn(t) = Znj for t 2 ~Inj =: (tnj?1;tnj]; j = 1;??? ;m2n?1:
In section 3.3.1 (homogeneous initial and boundary conditions), we obtained (3.84)
jjWnj jj1 ? jjQjjp2GSe and jjZnj jj ? jjQjjSe ; (3.182)
and in section 3.3.5 (nonhomogeneous boundary conditions), we obtained (3.171)
jjWnj jj1 ? C
r
?T
G ; and jjZ
n
j jj ? C
r
2?T
Se : (3.183)
Therefore since the poroelasticity system is linear, we have
jjWnj jj1 ? jjQjjp2GSe +C
r
?T
G and jjZ
nj jj ? jjQjj
Se +C
r
2?T
Se : (3.184)
100
From which it follows that
jjUn(t)jj1 ? jjQjjp2GSe +C
r
?T
G ; jjPn(t)jj ?
jjQjj
Se +C
r
2?T
Se 8t 2 I = [0;T];(3.185)
and
jj~un(0)?un(0)jj1 = jjwn1 ?wn0jj1 = jjhnWn1 jj1 ?
?
jjQjjp
2GSe +C
r
?T
G
!
hn;
~un(t)?un(t) = wnj ?wnj?1 ?(t?tnj?1)Wnj
= ?hn ?(t?tnj?1)?Wnj 8t 2 ~Inj :
Similarly,
jj~pn(0)?pn(0)jj = jjzn1 ?zn0jj = jjZn1hnjj?
?
jjQjj
Se +C
r
2?T
Se
!
hn;
~pn(t)?pn(t) = znj ?znj?1 ?(t?tnj?1)Znj
= ?hn ?(t?tnj?1)?Znj 8t 2 ~Inj :
We have in ~Inj : 0 < t?tnj?1 ? hn, then
jj~un(t)?un(t)jj1 ?
?
jjQjjp
2GSe +C
r
?T
G
!
hn; 8 t 2 I = [0;T]; (3.186)
and
jj~pn(t)?pn(t)jj ?
?
jjQjj
Se +C
r
2?T
Se
!
hn; 8 t 2 I = [0;T]: (3.187)
101
Rewriting the system (3.95) and (3.96) for n (instead of for nk), then for almost every t 2 I
we have
a(~un;v) ? fi(~pn;r?v) = ?(1?fl)fi < ~pn;v >; (3.188)
and
Se(Pn;q) + fi(r?Un;q) + k?(r~pn;rq) =
(Q;q) + < h1;q > + (1?fl)fi < Un;q > : (3.189)
Rewrite (3.188) and (3.189) for m
a(~um;v) ? fi(~pm;r?v) = ?(1?fl)fi < ~pm;v >; (3.190)
and
Se(Pm;q) + fi(r?Um;q) + k?(r~pm;rq) =
(Q;q) + < h1;q > + (1?fl)fi < Um;q > : (3.191)
Subtracting (3.190) from (3.188) and (3.191) from (3.189), we get
a(~un ? ~um;v) ? fi(~pn ? ~pm;r?v) = ?(1?fl)fi < ~pn ? ~pm;v >; (3.192)
Se(Pn ?Pm;q) + fi(r?Un ?r?Um;q) + k?(r~pn ?r~pm;rq) =
+ (1?fl)fi < Un ?Um;q > : (3.193)
102
Letting v = Un ?Um, q = ~pn ? ~pm, and adding (3.192) and (3.193), we obtain
a(~un?~um;Un?Um)+Se(Pn?Pm; ~pn? ~pm)+ k?(r~pn?r~pm;r~pn?r~pm) = 0: (3.194)
That is
a(~un?~um?(un?um);Un?Um)+Se(Pn?Pm; ~pn?~pm?(pn?pm))
+k?(r~pn?r~pm;r~pn?r~pm) = ?a(un?um;Un?Um)?Se(Pn?Pm;pn?pm): (3.195)
Since Un ?Um and Pn ?Pm are the derivatives of un ?um and pn ?pm, respectively, then
a(un ?um;Un ?Um) = 12 ddt[[un ?um]]2
and
(pn ?pm;Pn ?Pm) = 12 ddtjjpn ?pmjj2:
Furthermore, using the continuity of the bilinear form a(?;?), (3.11) and denoting by C1 =
max
?
2G; 6G?1?2?
?
, (3.195) becomes
1
2
d
dt[[un ?um]]
2 + Se
2
d
dtjjpn ?pmjj
2
? C1
?
jj~un ? ~um ?(un ?um)jj1jjUn ?Umjj1
?
+Sejj~pn ? ~pm ?(pn ?pm)jj jjPn ?Pmjj
? C1
?
jj~un ?unjj1 +jj~um ?umjj1
??
jjUnjj1 +jjUmjj1
?
+Se
?
jj~pn ?pnjj+jj~pm ?pmjj
? ?
jjPnjj+jjPmjj
?
(3.196)
103
From (3.185){(3.187), (3.196) becomes
1
2
d
dt[[un ?um]]
2 + Se
2
d
dtjjpn ?pmjj
2
? 2C1
?
jjQjjp
2GSe +C
r
?T
G
!2
(hn +hm)
+2Se
?
jjQjj
Se +C
r
2?T
Se
!2
(hn +hm): (3.197)
Let us denote by
C2 = 2
2
4C1
?
jjQjjp
2GSe +C
r
?T
G
!2
+Se
?
jjQjj
Se +C
r
2?T
Se
!23
5;
then (3.197) becomes
d
dt
?[[u
n ?um]]2 +Sejjpn ?pmjj2
?? C
2(hn +hm): (3.198)
Integrating (3.198) from 0 to T, we obtain
[[un ?um]]2 +Sejjpn ?pmjj2 ?
Z T
0
C2(hn +hm)dt; 8t 2 I:
Thus
2Gjjun(t)?um(t)jj21 +Sejjpn(t)?pm(t)jj2 ? C2T(hn +hm): (3.199)
We have um(t) ! u(t) in H1(?) for almost every t 2 I, pm(t) ! p(t) in L2(?) for almost
every t 2 I, and hm ! 0 for m !1, thus
2Gjjun(t)?u(t)jj21 ? C2Thn; 8 t 2 I;
104
and
Sejjpn(t)?p(t)jj2 ? C2Thn; 8 t 2 I:
Hence
jjun(t)?u(t)jj1 ?
r
C2Thn
2G ; 8 t 2 I; (3.200)
and
jjpn(t)?p(t)jj?
r
C2Thn
Se ; 8 t 2 I; (3.201)
where
C2 = 2
2
4C1
?
jjQjjp
2GSe +C
r
?T
G
!2
+Se
?
jjQjj
Se +C
r
2?T
Se
!23
5;
and
C1 = max
 
2G; 6G?1?2?
?
:
105
Error estimates for nonhomogeneous initial conditions
Consider now the poroelasticity system (3.19){(3.26) with nonhomogeneous initial condi-
tions. Then in section 3.3.3, we obtained (3.143) and (3.144)
jjWnj jj1 ? 2kp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se
and
jjZnj jj ? 2kp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se :
Denote by
C = 2kp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se ;
then
jjWnj jj1 ? C and jjZnj jj ? C:
The functions un(t), pn(t), ~un(t), ~pn(t), Un(t), and Pn(t) are as deflned above for the
homogeneous initial condition case. Therefore,
jjUn(t)jj1 ? C; jjPnjj ? C; (3.202)
jj~un(0)?un(0)jj1 = jjwn1 ?wn0jj1 = jjhnWn1 jj1 ? Chn;
and
~un(t)?un(t) = wnj ?wnj?1 ?(t?tnj?1)Wnj
= ?hn ?(t?tnj?1)?Wnj 8t 2 ~Inj :
106
Similarly,
jj~pn(0)?pn(0)jj = jjzn1 ?zn0jj = jjhnZn1jj? Chn;
and
~pn(t)?pn(t) = znj ?znj?1 ?(t?tnj?1)Znj
= ?hn ?(t?tnj?1)?Znj 8t 2 ~Inj :
Since in ~Inj : 0 < t?tnj?1 ? hn, then
jj~un(t)?un(t)jj1 ? Chn; 8 t 2 I = [0;T]; (3.203)
and
jj~pn(t)?pn(t)jj ? Chn; 8 t 2 I = [0;T]: (3.204)
For almost every t 2 I, the system (3.117) and (3.118), corresponding to mesh dn, can be
written as
a(~un;v) ? fi(~pn;r?v) = ?(1?fl)fi < ~pn;v >; (3.205)
Se(Pn;q) + fi(r?Un;q) + k?(r~pn;rq) = (1?fl)fi < Un;q > : (3.206)
Rewriting (3.205) and (3.206) for m, we get
a(~um;v) ? fi(~pm;r?v) = ?(1?fl)fi < ~pm;v >; (3.207)
Se(Pm;q) + fi(r?Um;q) + k?(r~pm;rq) = (1?fl)fi < Um;q > : (3.208)
107
Subtracting now (3.207) from (3.205) and (3.208) from (3.206), we get
a(~un ? ~um;v) ? fi(~pn ? ~pm;r?v) = ?(1?fl)fi < ~pn ? ~pm;v >; (3.209)
Se(Pn ?Pm;q) + fi(r?Un ?r?Um;q) + k?(r~pn ?r~pm;rq) =
(1?fl)fi < Un ?Um;q > : (3.210)
Adding (3.209) and (3.210) with v = Un ?Um and q = ~pn ? ~pm, we get
a(~un ? ~um;Un ?Um)+Se(Pn ?Pm; ~pn ? ~pm)+ k?(r~pn ?r~pm;r~pn ?r~pm) = 0;
which is exactly (3.194). Hence we can obtain (3.196)
1
2
d
dt[[un ?um]]
2 + Se
2
d
dtjjpn ?pmjj
2
? C1
?
jj~un ?unjj1 +jj~um ?umjj1
??
jjUnjj1 +jjUmjj1
?
+Se
?
jj~pn ?pnjj+jj~pm ?pmjj
? ?
jjPnjj+jjPmjj
?
; (3.211)
where C1 = max
?
2G; 6G?1?2?
?
is the continuity constant for the bilinear form a(?;?).
Using now (3.202){(3.204), (3.211) becomes
d
dt
?
[[un ?um]]2 +jjpn ?pmjj2
?
? 2 C
2
min(1;Se)(C1 +Se)(hn +hm): (3.212)
Integrating (3.212) from 0 to T, we get
2Gjjun(t)?um(t)jj21 +jjpn(t)?pm(t)jj2 ? 2 C
2
min(1;Se)(C1 +Se)T(hn +hm):
108
We have um(t) ! u(t) in H1(?) for almost every t 2 I, pm(t) ! p(t) in L2(?) for almost
every t 2 I, and hm ! 0 for m !1, hence
jjun(t)?u(t)jj1 ? C
s
(C1 +Se)
Gmin(1;Se)Thn; 8 t 2 I; (3.213)
and
jjpn(t)?p(t)jj? C
s
2(C1 +Se)
min(1;Se) Thn; 8 t 2 I; (3.214)
with
C = 2kp?? max(C1;Se)
(min(2G;Se))32
jjrv0jj
Se and C1 = max
 
2G; 6G?1?2?
?
:
109
3.4.2 Error estimates for the fully discrete problem
The weak formulation of the quasi-static poroelasticity problem is: flnd u 2 V and
p 2 M such that
a(u;v)?fi(p;r?v) = (F;v)?(1?fl)fi < p;v >; 8v 2 V; (3.215)
Se(pt;q)+fi(r?ut;q)+ k?(rp;rq) =
(Q;q)+ < h1;q > +(1?fl)fi < ut;q >; 8q 2 M; (3.216)
for almost every t 2 I. Using backward time discretization, we get
a(ui;v)?fi(pi;r?v) = (Fi;v)?(1?fl)fi < pi;v >; 8v 2 V; (3.217)
Se(pi ?pi?1;q)+fi(r?ui ?r?ui?1;q)+hk?(rpi;rq) =
h(Qi;q)+h < h1i;q > +(1?fl)fi < ui ?ui?1;q >; 8q 2 M: (3.218)
Let V h ? V and Mh ? M be flnite dimensional spaces. The weak formulation for the
discrete problem is: flnd uh 2 V h and ph 2 Mh such that
a(uhi ;vh)?fi(phi ;r?vh) = (Fi;vh)?(1?fl)fi < phi ;vh >; 8vh 2 V h; (3.219)
Se(phi ?phi?1;qh)+fi(r?uhi ?r?uhi?1;qh)+hk?(rphi ;rqh) =
h(Qi;qh)+h < h1i;qh > +(1?fl)fi < uhi ?uhi?1;qh >; 8qh 2 Mh: (3.220)
110
Since V h ? V and Mh ? M, we have
a(ui;vh)?fi(pi;r?vh) = (Fi;vh)?(1?fl)fi < pi;vh >; 8vh 2 V h; (3.221)
Se(pi ?pi?1;qh)+fi(r?ui ?r?ui?1;qh)+hk?(rpi;rqh) =
h(Qi;qh)+h < h1i;qh > +(1?fl)fi < ui ?ui?1;qh >; 8qh 2 Mh: (3.222)
Subtracting (3.221) from (3.219) and (3.222) from (3.220), we obtain
a(uhi ?ui;vh)?fi(phi ?pi;r?vh) = ?(1?fl)fi < phi ?pi;vh >; (3.223)
Se(phi ?pi;qh)+fi(r?uhi ?r?ui;qh)+hk?(rphi ?rpi;rqh) = Se(phi?1?pi?1;qh)
+fi(r?uhi?1?r?ui?1;qh)+(1?fl)fi < uhi ?ui;qh > ?(1?fl)fi < uhi?1?ui?1;qh > : (3.224)
Letting vh = uhi ?wh and qh = phi ?zh and using
a(uhi ?ui;uhi ?wh) = a(uhi ?ui;uhi ?ui)+a(uhi ?ui;ui ?wh), we get
a(uhi ?ui;uhi ?ui)?fi(phi ?pi;r?uhi ?r?wh) =
?(1?fl)fi < phi ?pi;uhi ?wh > ?a(uhi ?ui;ui ?wh); (3.225)
and
111
Se(phi ?pi;phi ?pi)+fi(r?uhi ?r?ui;phi ?zh)+hk?(rphi ?rpi;rphi ?rpi) =
Se(phi?1 ?pi?1;phi ?zh)+fi(r?uhi?1 ?r?ui?1;phi ?zh)+(1?fl)fi < uhi ?ui;phi ?zh >
?(1?fl)fi < uhi?1?ui?1;phi ?zh > ?Se(phi ?pi;pi?zh)?hk?(rphi ?rpi;rpi?rzh): (3.226)
Adding (3.225) and (3.226) and using the coercivity and the continuity of the bilinear form
a(?;?) (with C is the continuity constant), we have
2Gjjui?uhijj21+Sejjpi?phijj2+hk?jjrpi?rphijj2
? fijjpi?phijjjjr?uhi?r?whjj+fijjr?ui?r?uhijjjjphi?zhjj+Cjjui?uhijj1jjui?whjj1
+Sejjpi?phijjjjpi?zhjj+hk?jjrpi?rphijjjjrpi?rzhjj+Sejjpi?1?phi?1jjjjphi?zhjj
+fijjr?ui?1 ?r?uhi?1jj jjphi ?zhjj+(1?fl)fijjpi ?phijj?tfjjuhi ?whjj?tf
+(1?fl)fijjui?uhijj?tfjjphi ?zhjj?tf +(1?fl)fijjuhi?1?ui?1jj?tfjjphi ?zhjj?tf: (3.227)
Using now jjr?ujj ? p3jjrujj ? p3jjujj1, jjpi ?phijj?tf ? jjpi ?phijjM ? jjpi ?phijj1, and
Young?s inequality (ab ? 12?a2 + ?2b2), (3.227) becomes
112
2Gjjui?uhijj21+Sejjpi?phijj2+hk?jjrpi?rphijj2
? 12?
1
?
fip3jjuhi ?whjj1 +Sejjpi ?zhjj
?2
+ ?12 jjpi?phijj2
+ 12?
2
?
fip3jjphi ?zhjj+Cjjui ?whjj1
?2
+?22 jjui?uhijj21+hk??32 jjrpi?rphijj2
+hk? 12?
3
jjrpi?rzhjj2+Sejjpi?1?phi?1jjjjphi?zhjj+fip3jjui?1?uhi?1jj1jjphi?zhjj
+(1?fl)fi?42 jjpi ?phijj21 +(1?fl)fi 12?
4
jjui ?whjj21 +(1?fl)fi?52 jjui ?uhijj21
+(1?fl)fi 12?
5
jjphi ?zhjj21 +(1?fl)fijjuhi?1 ?ui?1jj1jjphi ?zhjj1: (3.228)
By Poincare?s inequality jjrpi ?rphijj2 ? ?jjpi ?phijj21, for some ? > 0. Choose ?3 and ?4
su?ciently small such that 2hk???(hk??3 ?(1?fl)fi?4) > 0, choose ?2 and ?5 small enough
so that 4G?(?2 +(1?fl)fi?5) > 0, and choose ?1 su?ciently small so that (2Se??1) > 0.
Furthermore, jjpi ?zhjj? hjjpijjH2(?) and jjphi ?zhjj? hjjpijjH2(?). Hence (3.228) becomes
(4G?(?2+(1?fl)fi?5))jjui?uhijj21+(2Se??1)jjpi?phijj2
? h
2
?1
?
fip3jjuijjH2(?) +SejjpijjH2(?)
?2
+h
2
?2
?
fip3jjpijjH2(?) +CjjuijjH2(?)
?2
+h2 k??
3
jjpijj2H2(?) +2hSejjpi?1 ?phi?1jj jjpijjH2(?)
+2hfip3jjui?1?uhi?1jj1jjpijjH2(?)+h2(1?fl)fi?
4
jjuijj2H2(?)
+h2(1?fl)fi?
5
jjpijj2H2(?) +2h(1?fl)fijjuhi?1 ?ui?1jj1jjpijjH2(?): (3.229)
113
Using again Young?s inequality, we obtain
(4G?(?2+(1?fl)fi?5))jjui?uhijj21+(2Se??1)jjpi?phijj2
? h
2
?1
?
fip3jjuijjH2(?) +SejjpijjH2(?)
?2
+h
2
?2
?
fip3jjpijjH2(?) +CjjuijjH2(?)
?2
+h2jjpijj2H2(?)
? k
??3 +Se
2jjpi?1 ?ph
i?1jj
2 +fi2(3+(1?fl)2)jjui?1 ?uh
i?1jj
2
+(1?fl)fi?
5
?
+h2(1?fl)fi?
4
jjuijj2H2(?): (3.230)
In section 3.3 we derived energy norm estimates (3.178) and (3.181)
jju(t)jj1 ? C1 +C3 +C5; jjp(t)jj? C2 +C4 +C6;
and
jjuh(t)jj1 ? C1 +C3 +C5; jjph(t)jj? C2 +C4 +C6;
where
C1 =
peT
p2Gmin(1;Se)jjQjjL2(0;T;L2(?));
C2 =
peT
min(1;Se)jjQjjL2(0;T;L2(?));
C3 =
s
jjv0jj2
2GSe +
T
2?G
 k
?
jjrv0jj
Se
?2
;
C4 =
s
jjv0jj2
Se2 +
T
?Se
 k
?
jjrv0jj
Se
?2
;
114
C5 =
r
?1T
G jjh1jjL2(0;T;L2(?tf));
and
C6 =
r
2?1T
Se jjh1jjL2(0;T;L2(?tf)):
Therefore,
jjui?1 ?uhi?1jj2 ?
?
jjui?1jj+jjuhi?1jj
?2
? 4(C1 +C3 +C5)2 (3.231)
and
jjpi?1 ?phi?1jj2 ?
?
jjpi?1jj+jjphi?1jj
?2
? 4(C2 +C4 +C6)2 : (3.232)
Hence (3.230) becomes
(4G?(?2+(1?fl)fi?5))jjui?uhijj21+(2Se??1)jjpi?phijj2
? h
2
?1
?
fip3jjuijjH2(?) +SejjpijjH2(?)
?2
+h
2
?2
?
fip3jjpijjH2(?) +CjjuijjH2(?)
?2
+h2jjpijj2H2(?)
? k
??3 +4Se
2 (C2 +C4 +C6)2 +4fi2(3+(1?fl)2)(C1 +C3 +C5)2
+(1?fl)fi?
5
?
+h2(1?fl)fi?
4
jjuijj2H2(?): (3.233)
Denoting by K the right hand side of (3.233), we get
115
jjui ?uhijj1 ?
s
K
4G?(?2 +(1?fl)fi?5) (3.234)
and
jjpi ?phijj1 ?
r K
2Se??1 (3.235)
116
Chapter 4
Numerical methods
Our objective is to approximate concurrently solutions of the system of partial difier-
ential equations
?Gr2u? G1?2?r(r?u)+firp = F; in ??(0;T); (4.1)
@
@t(Se p+fir?u)?r?(
k
?rp) = Q; in ??(0;T); (4.2)
for both the solid displacement u (a vector fleld) and the  uid pressure p (a scalar fleld).
To this end, we developed several algorithms: 2dp ow, 3dp ow, and 3dupfem. These
algorithms approximate the solution of each equation separately, approximating the dis-
placement u in equation (4.1) assuming that the pressure p is given or approximating p in
equation (4.2) assuming that u is given. For the fully coupled system, we flrst developed a
segregated algorithm (it3dupfem) then a coupled algorithm (c3dupfem).
4.1 2-D algorithm for the difiusion equation: 2dp ow
A 2-dimension flnite element method (2dp ow) was used to approximate the solution
of the difiusion equation:
( @@t(Se?p+fir?u)?r?(k?rp) = Q
117
for the  uid pressure p, assuming that the vector displacement u is known. The domain
considered was a box with Dirichlet boundary condition on the top and homogeneous Neu-
mann boundary condition on the left/right side and bottom.
A brief description of the numerical method:
Let V = fv : rv is a piecewise continuous on ? and vj? = 0g:
We start with flnite element discretization in space: we multiply the difiusion equation by
a test function v and integrate over the domain ? to obtain
Z
?
Se pt ?v d??
Z
?
k
?(4p)?v d? =
Z
?
f ?v d?:
Here f is our right hand side consisting of the two terms: the source term Q and the term
containing the known vector displacement u.
Applying Green?s formula and using the boundary condition, we get
Z
?
Sept ?v d?+
Z
?
k
?(rp)?(rv) d? =
Z
?
f ?v d?:
To approximate a solution on ??(0;T), divide (0;T) into n subintervals, each of length
? = Tn, and p(x;n?) ? pn(x).
Using flnite difierence backward time discretization, we obtain
Z
?
Sep
n+1 ?pn
? ?v d?+
Z
?
k
?(rp
n)?(rv) d? =
Z
?
fn ?v d?:
The superscript n denotes the discrete time level at which the function is evaluated and ?
is the time step.
118
Rearranging the previous equation and assuming that Se, k, and ? are constants, then we
have
Z
?
pn+1 ?v d? =
Z
?
pn ?v d?? ?Se k?
Z
?
(rpn)?(rv) d?+ ?Se
Z
?
fn ?v d?: (4.3)
Construct V h ? V, where V h is a flnite dimensional space (the set of all functions which are
linear on each subinterval and continuous on ?). Construct a basis for V h, choose ?j 2 V h,
1 ? j ? n, with
?j(xi) =
8
><
>:
1 if i = j, 1 ? i;j ? n;
0 if i 6= j, 1 ? i;j ? n:
Equation (4.3) is a system of linear algebraic equations of the form
M pn+1 = M pn ?C1 A pn +C2 b;
M pn+1 = (M ?C1 A)pn +C2 b; (4.4)
where C1 = ?Se k?, C2 = ?Se, M = R? ?(i)??(j), A = R?(r?(i)?r?(j)), and b = R? f?(j).
Now, instead of expressing the right hand side of (4.4) entirely at time n, it is averaged at
n and n+1. This is called the Crank-Nicolson method, the result is as follows
M pn+1 ?M pn = ?C12 A pn+1 ? C12 A pn +C2 b:
Or equivalently,
(M + C12 A)pn+1 = (M ? C12 A)pn +C2 b:
We then approximate this system of equations for the scalar pore pressure p.
119
To test this program, data for which the exact solution is known are generated and compared
to the approximate solution obtained from the developed algorithm.
The graph of 2Dp ow from MATLAB comparing the exact solution and the approximate
solution is shown below:
120
0 0.2 0.4 0.6 0.80
0.2
0.4
0.6
0.8
The approximate sol.
x
y
?1.5
?1
?0.5
0
0.5
1
1.5
0 0.2 0.4 0.6 0.80
0.2
0.4
0.6
0.8
the exact sol. (mesh)
x
y
?1.5
?1
?0.5
0
0.5
1
1.5
0
0.5
1
0
0.5
1?2
0
2
x
The approximate sol.
y
?1.5
?1
?0.5
0
0.5
1
1.5
0
0.5
1
0
0.5
1?2
0
2
x
the exact sol.
y
?1.5
?1
?0.5
0
0.5
1
1.5
2
Figure 4.1: Comparing approximate solution pn and exact solution p from 2dp ow at the
last time step
Figure 4.1 depicts the  uid pressure in the square (0, 1)x(0, 1) at the last time step (n = 1)
for the approximate solution pn and the exact solution p. As we can see from the graph, the
approximate solution pn on the left hand side looks exactly the same as the exact solution
p on the right hand side.
121
4.2 3-D algorithm for the difiusion equation: 3dp ow
The same equation - the difiusion equation - is solved for the pressure p assuming
u is given using a 3-dimensional flnite element discretization in space and second order
Crank Nicolson discretization in time. We consider the equation posed on a cube with
homogeneous Dirichlet boundary conditions. The (MATLAB) code is 3dp ow.
Again the approximate solution and the exact solution (we solve the equation for data for
which the exact solution is known) are compared to test and validate the program.
The plot for the approximate solution and the exact solution at the last time step is shown
below.
Figure 4.2: Comparing approximate solution pn and exact solution p from 3dp ow at the
last time step
122
From Figure 4.2, we clearly see that the approximate solution pn and the exact solution
p are similar which is evidence for the validity of our code 3dp ow.
4.3 3-D algorithm for the elasticity equation: 3dfem
The program 3dfem uses a 3-dimensional flnite element method to approximate the
displacements u in the elasticity equation:
?Gr2u? G1?2?r(r?u)+firp = F;
assuming that the pore pressure p is given. The equation was approximated in a box with
homogeneous Dirichlet boundary conditions.
This program was tested the same way by approximating the displacement un and com-
paring it to the exact solution u. So, the exact solution and the approximate solution are
compared in the following plot.
123
Figure 4.3: Comparing approximate solution un and exact solution u from 3dfem
Here we are plotting the vector displacement u in the box (0,1)x(0,1)x(0,1). The flrst,
second, and third row correspond to the x-component, y-component, and z-component of
the displacement u respectively. The left column is for the approximate solution un, and the
right column is for the exact solution u providing evidence for the validity of 3dfem (since
the graphs for the exact solution look exactly the same as the ones for the approximate
solution as shown in Figure 4.3).
124
4.4 Segregated algorithm: it3dupfem
This program approximates solutions of the system of the two partial difierential equa-
tions: the elasticity equation and the difiusion equation with an iterative method using
3-dimensional flnite element method. That is, approximating the elasticity equation
?Gr2u? G1?2?r(r?u)+firp = F;
for the displacements u. The body force per unit bulk volume F is set to be the gravity
force, and the pressure p is initialized using the program 3D-DEF (Gomberg and Ellis [7]).
Then the displacement is used in the difiusion equation:
@
@t(Sep+fir?u)?r?(
k
?rp) = Q;
to solve for the pore pressure p.
In this equation - the difiusion equation - the initial pressure is set as before using 3D-DEF.
The system solved at each time step is:
?G4(un+1(i) )? G1?2?r(r?(un+1(i) ))+fir(pn+1(i) ) = Fn+1; (4.5)
S
pn+1(i+1) ?pn(i)
? +fi
r?(un+1(i) )?r?(un(i))
? ?
k
?4p
n+1
(i+1) = Q
n: (4.6)
The superscript n denotes the discrete time level at which the function is evaluated and the
subscript i denoted the inner iteration (counter).
125
Giving un and pn and guessing pn+1(0) , we flrst calculate un+1(0) using equation (4.5) then we
substitute un+1(0) into equation (4.6) to flnd pn+1(1) .
The iteration yields un+1 = un+1(i) and pn+1 = pn+1(i+1). The process is repeated several times
until convergence, then the solution at the next time step is computed in a similar manner.
This algorithm did converge, i.e., the difierence between the previous calculated displace-
ment and the next calculated displacement is less than or equal to some tolerance, similarly
the difierence between the previous calculated pressure p and the next calculated pressure
p is less than or equal to some tolerance.
4.5 Coupled algorithm: c3dupfem
A coupled algorithm (with 3-D flnite element method) is used to approximate the
solution of the system of the two coupled partial difierential equations.
?Gr2u? G1?2?r(r?u)+firp = F;
@
@t(Sp+fir?u)?r?(
k
?rp) = Q:
In other words, after discretization using flnite elements in space and second order Crank-
Nicolson in time, the system of linear algebraic equations to be solved has the form:
2
64 A B
BT ?C
3
75
2
64U
P
3
75 =
2
64 eF
eQ
3
75
126
The system is solved for the vector displacement u and pore pressure p.
The(MATLAB)codec3dupfemapproximatedsolutionsofthesysteminthebox(0,1)x(0,1)x(0,1)
with homogeneous Dirichlet boundary condition. In this program, data for which the exact
solution is known are generated and compared to the approximate solution obtained from
the developed program. The following graph compares the approximate solution and the
exact solution for the pore pressure at the last time step.
0 0.2
0.4 0.6
0.8 1
0
0.5
1?1
0
1
2
xy
z
P: the approximate solution
?14
?12
?10
?8
?6
?4
?2
0x 10
?3
0 0.2
0.4 0.6
0.8 1
0
0.5
1?1
0
1
2
xy
z
P: the exact solution
?0.015
?0.01
?0.005
0
Figure 4.4: Comparing approximate solution pn and exact solution p from c3dupfem at the
last time step
127
The graph below compares the approximate solution and the exact solution for the
vector displacement at the last time step.
Figure 4.5: Comparing approximate solution un and exact solution u from c3dupfem at the
last time step
Figure 4.4 clearly shows that the approximate solution is the same as the exact so-
lution for the pressure, and in Figure 4.5 the approximate solution for the x-component,
y-component, and z-component of the displacement on the left hand side look exactly the
same as the ones of the exact solutions on the right hand side which is evidence for the
validity of c3dupfem code.
128
Chapter 5
Conclusion and future work
In this work, we considered the interaction between  uid pressure changes and the
deformation of a porous elastic material. Starting from the force equilibrium equation
and the linear constitutive equations, we formulated the equations describing the coupled
processes of elastic deformation and the pore  uid pressure in a porous medium. The
fully coupled system of equations does not in general yield closed form solutions. The
algorithm 3D-DEF (Gomberg and Ellis [7]) approximates Biot?s system for the displacement
from which the strain ? and the stress  can be calculated. In order to calculate pore
pressure changes, the 3P-Flow (see [9]) algorithm uses the above calculated strain ? and
stress  . In other words, 3D-DEF approximates the quasi-static elasticity equation for the
vector displacement u. Using these results, 3P-Flow then approximates the pressure in the
difiusion equation. Thus the two algorithms together do not approximate the fully coupled
system of the two partial difierential equations. Our main objective in this work was to
derive numerical algorithms for approximating solutions to the fully coupled system by
concurrently approximating solutions for the vector displacement u and the scalar pressure
p. This objective was attained. Our numerical algorithms were extensively tested. After
numerically approximating the fully coupled system, we considered the problem of existence
and uniqueness of solutions.
In [14] Showalter showed existence and uniqueness of strong and weak solutions using
abstract theory. This work proposed a constructive approach based on Babuska-Brezzi
theory and Rothe?s method to show existence and uniqueness of weak solutions for the
129
quasi-static poroelasticity system. Our approach suggested numerical methods which were
used to approximate solutions of the quasi-static poroelasticity system. Moreover, error
estimates were derived.
In the numerical experiment for the fully coupled system (c3dupfem), all the coe?cients
in the equilibrium equation for momentum conservation and the difiusion equation for Darcy
 ow were set to one except Poisson?s ratio ? that was set to 1=3. If we use the physical
coe?cients, then the matrix
M =
2
64 A B
BT ?C
3
75
has high condition number since this matrix M is \close" to being singular. Our future
work is to construct and solve the system with approximate Schur complement. In other
words, we compute the Schur complement of the matrix M and precondition it with its
diagonal. That is, we solve instead the following problem
D?1
2
64A B
0 ?C ?BTA?1B
3
75
2
64U
P
3
75 = D?1
2
64 eF
eQ?BTA?1eF
3
75
Here D is the diagonal matrix whose diagonal is
diag
2
64A B
0 ?C ?BTA?1B
3
75
The condition number of the matrix
D?1
2
64A B
0 ?C ?BTA?1B
3
75
130
is of order 1. It seems now that we can obtain accurate approximate solutions since the
matrix is far from being singular (this will be our future work).
131
Bibliography
[1] Biot, M. A. General theory of three dimensional consolidation. Journal of Applied
Physics 12 (1941), 155{164.
[2] Braess, D. Finite elements, 2 ed. Cambridge University Press, New York, 1997.
[3] Brezzi, F., and Fortin, M. Mixed and Hybrid Finite Element Methods. Springer-
Verlag, New York, 1991.
[4] Coussy, O. A general theory of thermoporoelasticity for saturated porous materials.
Transport in Porous Media, 4 (1989), 281{293.
[5] Evans, L. C. Partial Difierential Equations. American Mathematical Society, Provi-
dence, Rhode Island, 1998.
[6] Fung, Y. C. A First Course In Continuum Mechanics. Prentic-Hall, Inc., Englewood
Clifis, New Jersey, 1969.
[7] Gomberg, and Ellis. Crustal deformation model.
[8] Jaeger, J. C., and Cook, N. G. W. Fundamentals of Rock Mechanics. John Wiley
& Sons, Inc., New York, 1976.
[9] Lee, M.-K., and Wolf, L. W. Analysis of  uid pressure propagation in heteroge-
neousrocks: Implicationsforhydrologically-inducedearthquakes. Geophysical Research
Letters 25, 13 (July 1998), 2329{2332.
[10] Rektorys, K. The Method of Discretization in Time. D. Reidel Publishing Company,
Dordrecht, Holland, 1982.
[11] Rice, J., and Cleary, M. Some basic stress difiusion solutions for  uid-saturated
elastic porous media with compressible constituents. Reviews in Geophysics and Space
Physics 14 (1976), 227{241.
[12] Schrefler, B., Shiomi, T., Chan, A., Zienkiewicz, O., and Pastor, M. Com-
putational Geomechanics. Wiley, Chichester, 1999.
[13] Settari, A., and Mourits, F. Coupling of geomechanics and reservoir simulation
models. Computer Meth. and Adv. in Geomechanics (1994).
[14] Showalter, R. Difiusion in poro-elastic media. Jour. Math. Anal. Appl. 251 (2000),
310{340.
132
[15] Showalter, R. Difiusion in deforming porous media. Jour. Math. Anal. Appl. 25
(Oct. 2002), 115{139.
[16] Turcotte, D. L., and Schubert, G. Geodynamics Applications of Continuum
Physics to Geological Problems. John Wiley & Sons, Inc., New York, 1982.
[17] von Terzaghi, K. Theoretical Soil Mechanics. John Wiley & Sons, New York, 1943.
[18] Wang, H. F. Theory of Linear Poroelasticity with applications to Geomechanics and
Hydrogeology. Princeton University Press, Princeton, New Jersey, 2000.
[19] Wloka, J. Partial difierential equations. Cambridge University Press, New York,
1987.
[20] Zienkiewicz, O., Chang, C., and Battess, P. Drained, undrained, consolidating,
and dynamic behaviour assumptions in soils, limits of validity, vol. 30. Geotechnique,
1980.
133
Appendices
134
Appendix A
Notation and Inequalities
A.1 Notation for derivatives
Let U be a subset of Rn. Assume that u : U !R, x 2 U.
? Gradient vector: ru = ( @u@x1; @u@x2;??? ; @u@xn).
? Laplacian of u: 4u = Pni=1 @2u@x2
i
.
Vector-valued function
If m > 1 and u : U !Rm;u = (u1;u2;??? ;um), then
? Gradient matrix:
ru =
0
BB
BB
B@
@u1
@x1 ???
@u1
@xn
... ... ...
@um
@x1 ???
@um
@xn
1
CC
CC
CA
? If m = n, then divergence of u is
r?u =
nX
i=1
@ui
@xi
Multi-index notation
? @iu = @u@xi, i = 1;??? ;n.
? @mi u = @i???@i| {z }
m times
u, i = 1;??? ;n, m 2Z+, (@0i u = u).
135
? Let fi 2Zn+, fi = (fi1;??? ;fin)
Dfiu = @fi11 ???@finn u. The order of this derivative is the order of fi, i.e. jfij := Pni=1 fii.
A.2 Spaces of continuous and difierentiable functions
? C(U): the set of all continuous functions u : U !R.
? Ck(U), k 2N: the set of all continuous functions u : U !R with continuous partial
derivatives up to and including k.
? C0(U) = C(U) and C1(U) = \m2Z+Ck(U).
? Lp(U): the set of all functions u : U ! R such that u is Lebesgue measurable,
jjujjLp(U) < 1, where jjujjLp(U) = ?RU jfjpdx?1p (1 ? p < 1).
? L1(U) the set of all functions u : U ! R such that u is Lebesgue measurable,
jjujjL1(U) < 1:
? H1(U): space of all functions u 2 L2(U) whose flrst derivatives are square integrable.
? H2(U): space of all functions u 2 L2(U) whose flrst and second derivatives are square
integrable.
? H10(U): space of all functions u 2 H1(U) such that uj@U = 0.
? Wm;p(U): the set of all functions u 2 Lp(U) that have weak derivatives Dfiu 2 Lp(U)
for all fi 2ZN+ with jfij? m. Also,
jjujjWm;p(U) :=
 P
fi2ZN+jfij?m jjD
fiujjp
Lp(U)
?1
p if p < 1, and
jjujjWm;1(U) := maxjjDfiujjLp(U)jfi 2ZN+;jfij? m.
136
? Wm;p0 (U): we say that a function u is in Wm;p0 (U) if u is the limit in Wm;p(U), of a
sequence of Cm-functions with compact support in U.
137
Appendix B
Preliminaries
Deflnition 3 (see [10]): Let I = [0;T] and let H be a Hilbert space. A mapping y(t) : I ! H
is called an abstract function from I into H.
The set of all abstract functions continuous in I, equipped with the norm
jjyjjC(I;H) = max
t2I
jjy(t)jjH
is called the space C(I;H).
Deflnition 4 (see [10]): A simple function is an abstract function which attains, on I, only
a flnite number of "values" f1;??? ;fm 2 H, on (Lebesgue) measurable sets N1;??? ;Nm with
measures ?1;??? ;?m, respectively.
The Bochner integral of a simple function is deflned by
Z
I
y(t)dt =
mX
i=1
fi?i:
Measurable functions in the Bochner sense (see [10]) are functions which can be approxi-
mated, to arbitrary accuracy, by simple functions.
The space L2(I;H) is the space of functions which are square integrable in the Bochner
sense, i.e. Bochner integrable and satisfying
Z
I
jjy(t)jj2Hdt < 1
138
with the scalar product
(y1;y2)L2(I;H) =
Z
I
(y1(t);y2(t))Hdt
and the norm
jjyjj2L2(I;H) =
Z
I
jjy(t)jj2Hdt: (B.1)
Convergence yn ! y in L2(I;H) means that
limn!1
Z
I
jjy ?ynjj2Hdt = 0:
By the Riesz theorem, a primitive function Y(t) is deflned by
(Y(t);f)H =
Z t
0
(y(?);f)Hd? 8f 2 H:
Then
Y 2 C(I;H)
(that is, Y(t) is continuous abstract function in the interval I (see [10])) and
Y 2 AC(I;H)
(that is, Y(t) is absolutely continuous (see [10])) for every y 2 L2(I;H).
The derivative which is Y0(t) = y(t) in L2(I;H), exists almost everywhere.
139
Theorem B.1 There exists a constant c > 0 depending only on the domain G, such that
for every function u 2 W(1)2 (G) we have
jjujjL2(?) ? cjjujjW(1)
2 (G)
Consider the boundary value problem (bvp):
?4u = f in D;
u = 0 on @D;
where D ?RN is a bounded domain, and f : D !R is given.
Strong solution of (bvp): Given p 2 (1;1), then u 2 W2;p(D)\W2;p0 (D) satisfying the
partial difierential equation ?4u = f in the sense of weak derivatives is the strong solution
of (bvp).
Weak solution of (bvp): Given p 2 (1;1), then u 2 W1;p0 (D) satisfying RD(ru)?(rv) =
R
D f v for all v 2 W
1;p0
0 (D) is the weak solution of (bvp).
Poincare inequality (see [19]): Let ? be bounded and l = 1;2;???. Then there exists a
constant c dependent only on the diameter of ?, such that for all ` 2 W2;10 (?)
jj`jj2l ? c
X
jsj=l
Z
?
jDs`(x)j2dx:
140
Holder?s inequality (see [5]): Assume 1 ? p;q ? 1, 1p + 1q = 1. Then if u 2 Lp(U),
v 2 Lq(U), we have Z
U
juvjdx ? jjujjLp(U)jjvjjLq(U)
.
141