Analysis and  nite element approximation for nonlinear problems in
poroelasticity and bioconvection
by
Song Chen
A dissertation submitted to the Graduate Faculty of
Auburn University
in partial ful llment of the
requirements for the Degree of
Doctor of Philosophy
Auburn, Alabama
May 5, 2013
Keywords: Nemytskii operator, poroelasticity, bioconvection, Rothe?s method
Copyright 2013 by Song Chen
Approved by
Yanzhao Cao, Chair, Professor of Mathematics and Statistics
A. J. Meir, Co-chair, Professor of Mathematics and Statistics
Wenxian Shen, Professor of Mathematics and Statistics
Georg Hetzer, Professor of Mathematics and Statistics
Abstract
This dissertation is concerned with nonlinear systems of partial di erential equation
with solution dependent physical coe cients satisfying the Nemytskii assumption. Such
equations arise from two important application  elds: poroelasticity and bioconvection.
First, we consider a quasi-static poroelasticity model with dilatation dependent hy-
draulic conductivity and an implicit time derivative. We derive the existence and unique-
ness of solutions using the modi ed Rothe?s method, Brouwer?s  xed point Theorem and the
Sobolev embedding Theorem. Next we construct a  nite element approximation with linear
elements and establish the optimal error estimate. We then conduct numerical examples to
verify the convergence and simulate the di usion in a  uid saturated sponge.
Second, we study the bioconvection model, a coupled Navier-Stokes type equation, with
concentration dependent viscosity. We combine the theory of the Navier-Stokes equation
and the modi ed Rothe?s method to establish existence and uniqueness of solutions of both
steady and time dependent bioconvection. After that we perform  nite element analysis with
Taylor-Hood elements and prove the convergence theorem. Finally numerical examples are
constructed using lab data to verify the convergence of the numerical scheme and simulate
the convection pattern formed by micro-organisms inside a culture  uid.
ii
Acknowledgments
It would have been impossible to  nish my dissertation without the guidance of my
committee members and support from my family and friends. First of all, I would like to
gratefully and sincerely thank my adviser, Dr. Yanzhao Cao, for his research guidance,
passion, caring and patience. Not only did Dr. Cao provide all the advices and directions
as a professional mathematician, he also shared with me his valuable life experiences to
inspire my life and career, as a wise elder and a close friend. I also want to thank my co-
adviser, Dr. A. J. Meir, for his speciality in scienti c computing, my committee members Dr.
Wenxian Shen and Dr. George Hetzer for guiding my research for the past several years and
especially Dr. Xing Fang who was kind enough to be the university reader of my dissertation
and provided many valuable advices on my defence. Finally, I want to thank my parents
and my lovely wife, Yilanna Hu, who have been there supporting me, now and always.
iii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Poroelasticity with dilatation dependent hydraulic conductivity . . . . . . . 2
1.2 Bioconvection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Plan of dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4 Notations and Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Poroelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1 Steady poroelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Existence and uniqueness of a weak solution of quasi-static poroelasticity . . 16
2.3 Numerical approximation with the  nite element method . . . . . . . . . . . 33
2.4 Numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3 Steady bioconvection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1 Existence and uniqueness of a weak solution . . . . . . . . . . . . . . . . . . 51
3.2 Numerical approximation using the  nite element method . . . . . . . . . . . 60
3.3 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4 Time dependent bioconvection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.1 Existence and uniqueness of a solution . . . . . . . . . . . . . . . . . . . . . 73
4.2 Numerical approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.3 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
iv
4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5 Conclusion and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
v
List of Figures
2.1 The plot of errors for a  xed time step . . . . . . . . . . . . . . . . . . . . . . . 46
2.2 The plot of errors for time step k = h2 . . . . . . . . . . . . . . . . . . . . . . . 47
2.3 Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.4 Case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.5 Case 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1 Concentration and velocity  eld for  = 1% . . . . . . . . . . . . . . . . . . . . 70
3.2 Concentration and velocity  eld for  = 20% . . . . . . . . . . . . . . . . . . . 71
3.3 Concentration and velocity  eld for  = 20% with constant viscosity . . . . . . 71
3.4 Concentration and velocity  eld for  = 30% . . . . . . . . . . . . . . . . . . . 72
vi
List of Tables
2.1 The converge rate for a  xed time step in the L2(I;L2( )) norm . . . . . . . . . 45
2.2 The converge rate for a  xed time step in the L2(I;H10 ( )) norm . . . . . . . . 46
2.3 The converge rate for time step k = h2 in the L2(I;L2( )) norm . . . . . . . . . 46
2.4 The converge rate for time step k = h2 in the L2(I;H1( )) norm . . . . . . . . 46
3.1 Convergence rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.2 Parameter values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.1 Convergence rate in L2(I;L2( )) . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.2 Convergence rate in L2(I;H1( )) . . . . . . . . . . . . . . . . . . . . . . . . . . 104
vii
Chapter 1
Introduction
This dissertation is devoted to the study of the well posedness and numerical approxi-
mations of systems of partial di erential equation (PDE) with solution dependent physical
parameters arising from two application  elds: poroelasticity and bioconvection.
In mathematical modeling, we usually consider a linear system under certain ideal ho-
mogeneous assumptions about the physical parameters. In practice, however, the physical
parameters are generally inhomogeneous and related to the solution of the system of PDE.
A generic PDE with solution dependent physical coe cient can be written as
L
 
x;t;u(t;x); (t;x;u)
 
= f; x2 ; t> 0; (1.1)
where  is the physical domain, t refers to time, L denotes a general di erential operator
and  (t;x;u) is the physical parameter depending on the solution u. When L is a linear
di erential operator, (1.1) is often reduced to a quasi-linear PDE. Theoretical analysis of
this problem concerning the well posedness and regularity of solutions can be found in [1],
[2], and [3]. Finite element approximations of some simple quasi-linear elliptic PDEs and
parabolic PDEs (1.1) are presented in [4], [5] and [6].
In this dissertation, we consider two mathematical models taking the form of coupled
PDEs with solution dependent coe cients. The  rst model arises in poroelasitcity with di-
latation dependent hydraulic conductivity and the second model is about bioconvections with
concentration dependent viscosity. Though the physical backgrounds of the two problems
are quite di erent, we will use similar approaches to study the well posedness and numerical
approximations of the two underlying PDE systems. For instance we will use the modi ed
1
Rothe?s method to study the existence of weak solutions of both of the two evolution PDE
systems. We will use the fact that the physical parameters in both of the two models give
rise to Nemiskii operators to derive the well-posedness results under minimum regularity
assumptions on the input data.
1.1 Poroelasticity with dilatation dependent hydraulic conductivity
Motivation. We are surrounded by porous elastic solid materials: natural (e.g., rocks,
soils, shale, living tissue, the brain, the heart) and man-made (e.g., cement, concrete,  lters,
foams, ceramics, gels, clays). Porous material has a solid matrix structure with small pores
inside which contain air or  uid. Because of their ubiquity and unique properties, porous
materials are of great interest to natural scientists and engineers ([7, 8, 9, 10] ). Porous media
 nds applications in diverse areas include reservoir engineering [11], biomechanics [12, 13, 14]
and environmental engineering [15, 16, 11].
Mathematical model. Poromechanics is a branch of physics and speci cally contin-
uum mechanics and acoustics that studies the behaviour of  uid-saturated porous media.
The particular mathematical model we are interested in describes the swelling and shrinking
of an elastic deforming porous medium coupled with the  uid.
Due to the elastic nature of the porous medium, the study of  uid saturated porous
media is called poroelasticity. Terzaghi [17]  rst derived a one dimensional model in soil
mechanics to involve the in uence of the  uid inside a solid body. The model was later
extended to three dimension by Rendulik [18]. A mathematical formulation is derived in
Biot?s work [19] and studies in his other work between 1955 and 1962 [20, 21, 22, 23, 24],
which are later considered to be the foundation of modern poroelasticity theory.
2
The following equations are the essence of Biot?s poroelasticity theory [19] (details of
the theory can be found in [25]) in which the coupled constitutive equations takes the form
8
><
>:
@
@t(momentum) + stress = external force;
@
@t( uid content) +  ux = external  uid source:
(1.2)
Next we convert (1.2) into a coupled system of partial di erential equations for the  uid
pressure and the displacement of the poroelastic medium. To this end, we denote the pressure
by p and the displacement by u. According to Hooke?s law in 3-D the e ective stress caused
by the deformation of the solid matrix is given by
 e = 2 "+ tr(")I: (1.3)
Here " = 12(ru +ruT) is the strain, tr(") is the trace of ", I is the identity matrix, and
 ; are the Lam e constants, corresponding to the dilatation and shear modulus respectively,
given by
 = E (1 + )(1 2 ) ;  = E2(1 + ) ;
where the Young?s modulus E captures the elastic sti ness in the direction of a load and the
Poisson?s ratio  describes the stretch (compression) in the direction perpendicular to the
load. By introducing  and  , the e ective stress is written in a symmetric form. To include
the e ect of the  uid pressure inside the porous body, we introduce the addition stress due
to the pore pressure pp:
pp =  p; (1.4)
where the Biot-Willis constant  [23] satis es  = 1 K=Ks with K being the bulk modulus
of the porous matrix and Ks is the bulk modulus of the solid material. In most situations
  1, corresponding to an incompressible solid matrix. In a certain models of secondary
3
consolidation in clays, we also add an additional term
 s =   @@tr u; (1.5)
where   is the coe cient of the secondary consolidation [26]. Together, the total stress is
given by
 =  e +pp + s: (1.6)
In the second equation of (1.2), the  uid content is given by
 :=  f + d = c0p+ r u: (1.7)
Here  f = c0p measures the  uid content that can be forced into the medium by pressure
increments within constant volume with the constant c0  0 combining the porosity and the
compressibility, and  d =  r u denotes the  uid content due to the change of the pores size,
i.e., the dilatation of the void volume, which is proportional to r u, the total dilatation of
the body. If  = 1, the solid matrix is incompressible. Thus the total dilatation is the same
as the dilatation of the pores. According to Darcy?s law,
q =  rp (1.8)
represents the linear relationship between  uid  ux and the pressure drop, where the hy-
draulic conductivity  > 0 measures the permeability and the viscosity of the  uid. Now
letting the function f denote the volume-distributed external force, g the  uid source density,
and  the density of the medium, we substitute (1.3), (1.4), (1.5), (1.6), (1.7) and (1.8) into
(1.2) and apply divergence theorem to obtain the fully dynamic system
8>
<
>:
 @
2
@t2 u   r(
@
@tr u) ( + )r(r u)   u + rp = f ;
@
@t(c0p+ r u) r ( rp) = g;
(1.9)
4
subject to initial and boundary conditions. The  rst equation is a hyperbolic linear elasticity
equation for the displacement of the poroelastic medium while the second equation is an
implicit parabolic equation for the  uid?s pressure.
In general, complicated boundary conditions must be considered. In particular, the
boundary may have to be partitioned into disjoint, regular open set, on which various bound-
ary conditions are imposed. For example, we may partition the boundary  as  =   c[  t
with  c\ t = ;. We say  c is the clamped boundary on which the Dirichlet condition
uj c = ub is given, and  t is the traction boundary on which the normal stress  j t n = s is
prescribed. For the pressure, we may partition  as  =   d[  f with  d\ f =;. We say
 d is the drained boundary on which the Dirichlet condition pj d = pb is given, and  f is the
 ux boundary on which the  uid  ux ( (r u)rp)j f  n = r is prescribed. On  t\ f, we
have the balance of force and  ux at the same time. To handle this, we introduce  to be
the surface fraction of the sealed portion of  t\ f. Then 1  corresponds to the exposed
portion. Next we de ne the characteristic function
 tf =
8
><
>:
1; x2 t\ f ;
0; otherwise:
On the sealed part, the pore pressure contributes to the total force. Thus
h
( + )r uI + ru
i
 n   p n tf = s on  t;
On the exposed part, the changing rate of the  uid content caused by the dilatation must
be included in the balance of  ux. As a result we have
 @@t
 
(1  ) u n
 
 tf + rp n = r on  f :
Quasi-static poroelasticity. In this dissertation we restrict our considerations to the linear
quasi-static  ow in a saturated deformable poroelastic medium, i.e., we neglect the e ects
5
of  @
2
@t2 u and the secondary consolidation (1.5). Then the system takes the form
8>
<
>:
 ( + )r(r u)   u + rp = f ;
@
@t(c0p+ r u) r ( rp) = g:
(1.10)
In this case, we focus on a quasi-static problem with a coupling of an elliptic equation and
a parabolic equation.
In a homogeneous and isotropic medium, where the permeability and the viscosity are
constants,  is also a constant. In this case, the well posedness and regularity of system
(1.10) was studied in [25] as an application of the semi-group theory to linear degenerate
evolution equations in Hilbert spaces. Galerkin and discontinuous Galerkin approximations
for this linear system can be found in [27] and [28].
In this dissertation, we consider the case when the hydraulic conductivity  depends on
the dilatationr u, i.e.,  =  (r u). A typical example for dilatation dependent hydraulic
conductivity can be found in [29], which is used in calculating the pressure drop of a  uid
 owing through a packed bed of solids. Di culties in studying (1.10) arise from
1. Nonlinearity of the system due to the dependence of  on r u.
2. Implicit evolution in the second equation of (1.10) due to the term @@t( r u).
A proof of the existence and uniqueness of a solution of similar type of equations using
Rothe?s method can be found in Chapter 5 of [1] under the assumption that  : R!R is
continuous and satis es
0 <    (x)   8x2R; (1.11)
for some constants   ,   . In general, the uniqueness of the solution usually requires Lip-
schitz continuity condition of  and additional coe cient assumptions. The  nite element
approximations of similar problems can be found in Chapter 8.7 of [4] and in Chapter 13 of
[5]. However, none of these studies involves implicit evolution equations.
6
Linear implicit evolution equations have been studies in [30], [31], and [2] using the semi-
group theory. The quasi-static poroelasticity model with constant hydraulic conductivity was
studies in [25], also using the semi-group theory on Hilbert spaces. Formally, note that the
 rst equation in (1.10) is of second order for displacement u and  rst order for pressure p.
Therefore c0p and  r u should have the same regularity.
1.2 Bioconvection
Motivation. Bio-convection occurs due to on average upwardly swimming micro-
organisms which are slightly denser than water. The micro-organisms swim upward to meet
the sunlight. When the micro-organisms on the surface become too dense, they sink under
the e ect of gravity. Repeating this process, a convection pattern is formed.
Mathematical model. A  uid dynamical model treating the micro-organisms as col-
lections of particles was derived by M. Levandowsky, W. S. Hunter and E. A. Spiegel [32]
and independently by Y. Moribe [33]. We describe the model as follows. Let   R3 be a
bounded domain with boundary @ . At point x2 , let u(x) =fuj(x)g3j=1 and p(x) denote
the velocity and pressure of the culture  uid while c(x) refers to the concentration of the
micro-organism. Note that we use the volume concentration
c = nv0;
where n is the number of organisms per unit volume and v0 is the volume of an individual
organism. Let  0 be the density of the micro-organism and  m be the density of the culture
 uid. Then the density of the suspension is the sum
 =  0c+ m(1 c) =  m(1 + c); (1.12)
where  =  0= m 1. Assume that the organisms a ect the  uid dynamics only through
their in uence on its density and that the suspension is nearly incompressible. Then the
7
 uid satis es the Navier-Stokes equation
8
><
>:
 m
 @u
@t + (u r)u
 
 r ( ru) +rp = g i3 + f ;
r u = 0:
(1.13)
Here  is the kinematic viscosity of the culture  uid, g is the acceleration of gravity, f is
the volume-distributed external force, and i3 = (0;0;1) is the vertical unit vector. For the
concentration c, mass conservation gives
D
Dt c+r q = 0: (1.14)
Here DDt = @@t + u r is the material derivative and q is the  ux of micro-organisms which
is given by
q =  rc+Uci3; (1.15)
where  and U are the di usion rate and mean upward swimming velocity of the micro-
organism, respectively.
Combining (1.12), (1.13), (1.14), and (1.15), we derive the fully dynamic system, in  
8
>>>
>><
>>>
>>:
 m
 @u
@t + (u r)u
 
 r ( (c)ru) +rp = g m(1 + c)i3 + f ;
r u = 0;
@c
@t   c+ u rc+U
@c
@x3 = 0:
(1.16)
We assume the following boundary conditions for u and c, on  
8
><
>:
u = 0;
 @c@n Ucn3 = 0:
(1.17)
The second equation of (1.17) refers to a zero  ux condition on the boundary and n =
(n1;n2;n3) is the outward pointing unit normal vector on @ . We further assume the  xed
8
total mass for the micro-organisms, i.e.,
1
j j
Z
 
c(x)dx =  ; (1.18)
for some constant  . This means that no micro-organisms are allowed to leave or enter the
container. Finally the complete system describing the motion of micro-organisms takes the
form, in  
8
>>>
>>>
>>>>
>>>
><
>>>
>>>>
>>>
>>>
>:
 m
 @u
@t + (u r)u
 
 r ( (c)ru) +rp = g m(1 + c)i3 + f ;
r u = 0;
@c
@t   c+ u rc+U
@c
@x3 = 0;
u = 0;  @c@n Ucn3 = 0; on @ 
1
j j
Z
 
c(x)dx =  :
(1.19)
Note that the bioconvection model (1.19) is a special case of a more general equation
describing the di usion of an admixture in a region [34].
In an ideal Newtonian  uid, the viscosity  is a constant. In this case, the existence
of the solution as well as the positivity of the concentration are proved in [35] where the
author considered both the stationary and evolutionary cases. The evolutionary case of
system (1.16) with constant viscosity  is studied numerically in [36]. The numerical study
of slightly di erent bioconvection models can be found in [37], [38], [39], [40] and [41].
In practice, the viscosity is related to the concentration of the solute. Albert Einstein
 rst showed in his Ph.D thesis [42] that
 
 0 = 1 + c (1.20)
9
where  is the viscosity of the suspension,  0 is the viscosity of the pure solution, c is
the volume fraction of the particle spheres and  is a proportionality coe cient chosen
(experimentally) to be 2:5. Equation (1.20) is valid only for low concentration cases 0 <c<
10%. Therefore the result was later extended to add the c2 term by Batchelor [43] for larger
concentration (c 10%). When the concentration is much higher, the relative viscosity   0
varies as an exponential function of the concentration c ([44], [45] and [46]).
A recent work [47] showed the existence and uniqueness of a periodic solution of (1.19)
under the assumption that  ( ) is a C1 function and that for some positive constants   and
  
  < (x) <  8x2R and sup
x2R
 0(x) <1:
In this dissertation, we relax the above condition to assume that  : R!R is continuous
and
0 <    (x)   8x2R; (1.21)
for some constants   and   .
1.3 Plan of dissertation
In the rest of Chapter 1 we introduce notations and assumptions that will be used
throughout the thesis. In chapter 2 the existence of a weak solution of (1.10) with homo-
geneous boundary conditions will be proved using the modi ed Rothe?s approach ([48], [49]
and [1]) and analysis of implicit evolutionary equations [25] based on existing results for
the steady poroelasticity [50]. The uniqueness will also be proved under certain regularity
assumptions on the exact solution. A fully discrete  nite element approximation using linear
elements and backward Euler scheme will be studied and an optimal priori error estimate
will be proved using approximation theory and functional analysis ([5], Chapter 13). Nu-
merical experiment will be constructed to demonstrate the e ciency and the accuracy of
the numerical method. Then a numerical experiment with dilatation dependent hydraulic
10
conductivity  (r u) will be used to simulate a  uid saturated sponge. In chapter 3 and
4, the existence and uniqueness of a solution of (1.19) will be proved for both steady and
the fully dynamic case using the modi ed Rothe?s method combined with the theory of the
Navier-Stokes equations [51, 52]. A complete  nite element method will be constructed us-
ing Taylor-Hood elements [53] and a convergence theorem will be proved for both steady
and evolutionary cases. The error estimate will be established for the steady bioconvection
using theories in [52, 54, 55]. Numerical examples will be constructed to illustrate the con-
vergence rate and several practical simulations will be studied for the steady bioconvection
 ow. We conclude this dissertation in Chapter 5 with remarks and a plan for the future work.
1.4 Notations and Assumptions
Throughout the paper, we consider system (1.10) and (1.19) on an open bounded region
 in R2 or R3 with a smooth boundary  and on a time interval I = (0;T]. We denote by
C10 ( ) the space of in nitely di erentiable functions with compact support in  and by
L2( ) the space of square integrable functions on  . Let Wk;p( ) be the Sobolev space
consisting of functions in Lp( ) with each of their partial derivatives through order k also
in Lp( ). Speci cally, Hk( ) denotes the Hilbert space Wk;2( ). The space H10 ( ) is the
closure of C10 ( ) in the H1( ) norm. As usual, H 1( ) and (H1( )) denote the dual
of H10 ( ) and H1( ), respectively. Let Hk( ) =
 
Hk( )
 d
, H10( ) =
 
H10 ( )
 d
, and
L2( ) =
 
L2( )
 d
, d = 2;3, with k kk and k k denoting the respective norms of the two
spaces. Let H1=2( ) and H1=2( ) denote the trace space of H1( ) and H1( ) while H 1=2( )
is the dual space of H1=2( ). We shall use ( ; ) to denote both the L2 and L2 inner product
and denote by h ; i the duality pairing. The Poincar e inequality implies that there exists a
constant Cp such that
kwk1  Cpkrwk; 8w2H10( ) or H10 ( ): (1.22)
11
Throughout the paper, we will use C as a generic constant whose value may vary from
one occurence to the next.
12
Chapter 2
Poroelasticity
2.1 Steady poroelasticity
We  rst introduce some existing results for steady poroelasticty which is described by
the system (1.10) without the time derivative @@t(c0p+ r u), that is
8
><
>:
 ( + )r(r u)   u + rp = f ;
 r 
 
 (r u)rp
 
= g;
(2.1)
with boundary conditions
uj = ub;
and
 r 
 
 (r u)rp
 
j  n = r:
We assume that the body force f 2H 1( ), the  uid source g2(H1( )) , ub2H1=2( ), and
that r2H 1=2( ). Furthermore we assume that the data satis es the following compatibility
condition
hr;1i =hg;1i :
We denote by Q the quotient space H1( )=R and de ne the bilinear forms
8
>><
>>:
e(u;v) : =
 
( + )(r u);r v
 
+ ( ru;rv) 8u;v2H1( );
b(p;v) : =
Z
 
pr v 8p2Q; v2H1( ):
(2.2)
13
Multiply the  rst equation and second equation in (2.1) by test functions v 2 H10( ) and
q2Q respectively, the weak formulation of (2.1) takes the form
De nition 2.1. Given f 2 H 1( ), g2 (H1( )) , ub 2 H1=2( ), and r2H 1=2( ),  nd
(u;p)2H1( ) Q with uj = ub such that
8>
<
>:
e(u;v) b(p;v) =hf;vi 8v2H10( );
 
 (r u)rp;rq
 
=hg;qi +hl;qi 8q2Q:
(2.3)
Assumption (1.11) assures that  (r u(x)) is a so called Nemytskii operator whose
de nition and properties are stated in the following Lemma.
Lemma 2.2 ([56]). Assume that a function f :   Rm ! R satis es the Carath eodory
conditions:
(i) f(x;u) is a continuous function of u for almost all x2 ;
(ii) f(x;u) is a measurable function of x for all u2Rm. Furthermore for some constant
C and g2Lq( )
jf(x;u)j Cjujp 1 +g(x) x2 ; u2Rm
where 1 <q<1 and 1p + 1q = 1. Then the Nemytskii operator F(u) :  !R de ned by
F(u)(x) = f(x;u(x))
is a bounded and continuous map from Lp( ;Rm) into Lq( ;R).
The existence of a weak solution of (2.3) is given in the following Theorem.
Theorem 2.3 (Y. Cao, S. Chen and A. J. Meir). Given f 2 H 1( ), g 2 (H1( )) ,
ub2H1=2( ), and r2H 1=2( ), System (2.3) admits at least one weak solution satisfying
uj = ub.
Remark 2.4. The uniqueness of (2.3) can be obtained under additional coe cient assump-
tions and regularity assumptions if we assume that  is Lipschitz continuous (see [50]).
14
Finite element approximation. We construct the Galerkin  nite element approx-
imation fort he weak solution (u;p) of (2.3) on a convex polygonal, or polyhedral domain
 . Let  h be a family of quasi-uniform triangulations  h satisfying max 2 
h
diam   h. Then
we de ne the  nite dimensional subspace Qh  Q and Vh  H10( ) with the following
approximation properties
lim
h!0
inf
vh2Vh
kv vhk1 = 0 8v2H10( );
and
lim
h!0
inf
qh2Qh
kq qhk1 = 0 8q2Q:
The numerical scheme is to seek (uh;ph)2Vh Qh such that
8
><
>:
e(uh;vh) b(ph;vh) =hf;vhi 8vh2Vh;
 
 (r uh)rph;rqh
 
=hg;qhi +hl;qhi 8qh2Qh:
(2.4)
Following an argument similar to the proof of Theorem 2.3, we can show that the solution
(uh;ph) of (2.4) exists and the convergence (uh;ph) to the exact solution (u;p) of (2.3) is
established in what follows.
Theorem 2.5 (Y. Cao, S. Chen and A. J. Meir). Assume that the weak solution (u;p) 2
H1( ) Q of (2.3) is unique. In addition, assume that p2W1;1( ). Then
lim
h!0
(ku uhk1 +kp phk1) = 0:
15
2.2 Existence and uniqueness of a weak solution of quasi-static poroelasticity
We consider the quasi-static poroelasticity (1.10). For brevity, we assume that both u
and p satisfy homogeneous boundary conditions. Also for brevity, we set  = 1 (correspond-
ing to an incompressible solid matrix). For  6= 1, one may convert it to the  = 1 case by
rescaling the problem. We also assume that the external force f = 0. The nonzero case can
be handled through a simple transformation (see [25]).
The weak formulation. By de nition (2.2), the weak formulation of system (1.10) is
given in the following de nition
De nition 2.6. Given g in L2(I;L2( )) and l in L2( ), a pair (u;p) in H10( ) H10 ( ) is
said to be a weak solution of system (1.10) if it satis es for all t2I
8
>>>
>><
>>>
>>:
e(u;v) = (rp;v) 8v2H10( );
h@@t(c0p+r u);qi+
 
 (r u)rp;rq
 
= (g;q) 8q2H10 ( );
c0p( ;0) +r u( ;0) = l:
(2.5)
It is easy to verify that the bilinear form e( ; ) satis ed the hypothesis of Lax-Milgram
Theorem. Thus for a  xed t2I and p( ;t) in L2( ), the  rst equation of (2.5) can be solved
for u. De ne the operator B : L2( )!L2( ) such that for p2L2( ), Bp =r u where u
satis es 8
><
>:
e(u;v) = (rp;v) 8v2H10( );
u = 0 on @ :
The following Lemma can be found in [25].
Lemma 2.7. The operator B : L2( )!L2( ) de ned above is linear, continuous, monotone
and self-adjoint with Ker(B) = Ker(r) and Rg(B) = Ker(r)?.
16
Remark 2.8. Due to the assumption that the boundary is smooth, B is also continuous from
H10 ( ) into itself. Therefore there exist constants CB0andCB1 > 0 such that
kBqk CB0kqk 8q2L2( );
kBqk1  CB1kqk1 8q2H10 ( ):
(2.6)
Moreover, from Lemma 2.7 and the homogeneous boundary condition, Ker(B) = Ker(r) =
f0g. As a result, the mapping B : L2( ) ! L2( ) is one to one thus B is a continuous
bijection from L2( ) into itself. According to the bounded inverse Theorem, B and c0 + B
have bounded inverses.
Substitutingr u = Bp in the second equation of (2.5), we obtain the decoupled initial
value problem for p.
8
><
>:
h@@t(c0 +B)p;qi+
 
 (Bp)rp;rq
 
= (g;q) 8q2H10 ( );
(c0 +B)p( ;0) = l:
(2.7)
To prove the existence of a weak solution, we use the modi ed Rothe?s method to
construct a convergent sequence of approximate solutions of (2.7) using the backward Euler
approximation of the time derivative in (2.7). Let k = T=n for some positive integer n.
Partition I uniformly with time step k and denote nodal points by ti = tni = ik, for i =
1;2;:::;n. Let p0n2L2( ) be such that (c0 +B)p0n = l and de ne
8>
><
>>:
gin := 1=k
Z ti
ti 1
g(t)dt;
 (c0 +B)pin := (c0 +B)(pin pi 1n )=k; i = 1;:::;n:
(2.8)
We apply the following scheme inductively to obtain a sequence pin, i = 1;   ;n.
 
 (c0 +B)pin;q
 
+
 
 (Bpin)rpin;rq
 
= (gin;q) 8q2H10 ( ): (2.9)
17
Multiplying the above equation by k, we have
 
(c0 +B)pin;q
 
+k
 
 (Bpin)rpin;rq
 
= k(gin;q) +
 
(c0 +B)pi 1n ;q
 
8q2H10 ( ): (2.10)
To prove the existence of a solution of (2.10), we need the following direct corollary of
Brouwer?s  xed point Theorem.
Lemma 2.9. Let H be a  nite-dimensional Hilbert space with scalar product ( ; ) and the
corresponding norm j j. Let  be a continuous mapping from H to H and assume that there
exists  > 0 such that:
( (u);u) 0; 8u2H with juj=  :
Then there exists an elment u2H such that:
 (u) = 0; with juj  :
The following Lemma shows that (2.10) is well posed.
Lemma 2.10. Given pi 1n in L2( ), and gin 2L2( ), equation (2.10) has a weak solution
pin in H10 ( ), i = 1;   ;n.
Proof. For notational simplicity, we write  g = kgin + (c0 + B)pi 1n and  p = pin. Given  g in
L2( ), we show that there exists a  p in H10 ( ) such that
 
(c0 +B) p;q
 
+k
 
 (B p)r p;rq
 
= ( g;q) 8q2H10 ( ): (2.11)
For this purpose we let fqig1i=1 be an orthonormal basis of H10 ( ). Denote by Vm the  nite
dimensional space spanned by fq1;q2;:::;qmg and de ne the mapping  m : Vm!Vm by
( mq;w) =
 
(c0 +B)q;w
 
+k
 
 (Bq)rq;rw
 
 ( g;w) 8w2Vm:
18
From (1.11), (1.22), and the monotonicity of B,
( mq;q) =
 
(c0 +B)q;q
 
+k
 
 (Bq)rq;rq
 
 ( g;q)
 c0kqk2 + (Bq;q) +kk krqk2 k gkkqk
 (kk kqk1=C2p k gk)kqk1:
Thus ( mq;q)  0, for all q with kqk1 = C2pk gk=(kk ). Because Vm is  nite dimensional
and  de ned above is continuous, from Lemma 2.9, there exists  pm in Vm, such that
k pmk1  C2pk gk=(kk ) and  pm satis es  m pm = 0, i.e.
 
(c0 +B) pm;q
 
+k
 
 (B pm)r pm;rq
 
= ( g;q); 8q2Vj; j m: (2.12)
Since
k pmk1  C2pk gk=(kk ); (2.13)
i.e., f pmg1m=1 is a uniformly bounded sequence in H10 ( ), there exsits a subsequence of
f pmg1m=1, still denoted by f pmg1m=1, and a function  p2H10 ( ) such that
 pm *  p in H10 ( ): (2.14)
Due to the compact embedding of H10 ( ) into L2( ),
 pm!  p in L2( ): (2.15)
We now show that the weak limit  p is a solution of (3.9). Choose a test function
q2W1;1( )\H10 ( ): (2.16)
19
Using (2.6), (2.15), (2.16), and applying the Cauchy-Schwarz inequality we have
  
 
 
(c0 +B) pm;q
 
 
 
(c0 +B) p;q
   
 
 k(c0 +B)( pm  p)kkqk
 (c0 +CB0)kqkk pm  pk!0 as m!1:
The boundedness of B, the Nemytskii property of  and the fact  pm!  p imply that
 (B pm)! (B p) in L2( ) as m!1: (2.17)
From (1.11), (2.13), (2.14), (2.15), (2.16) and (2.17)
  
 
 
 (B p)r p;rq
 
 
 
 (B pm)r pm;rq
   
 
 
  
 
 
 (B p)r( p  pm);rq
   
 +
  
 
  
 (B p)  (B pm) r pm;rq
   
 
 k 
  
 
 
r( p  pm);rq
   
 +k (B p)  (B pm)kkr pmk
!0 m!1:
Because W1;1( ) is dense in H10 ( ), for all q in H10 ( ) we have that
lim
m!1
 
(c0 +B) pm;q
 
+
 
 (B pm)r pm;rq
 
=
 
(c0 +B) p;q
 
+
 
 (B p)r p;rq
 
:
Combining with (2.12), we obtain
 
(c0 +B) p;q
 
+
 
 (B p)r p;rq
 
= ( g;q) 8q2Vj: (2.18)
Since  nite linear combinations offqjg1j=1 are dense in H10 ( ), we conclude that  p is a solutin
of (3.9).
20
Lemma 2.11 (Discrete energy estimate). There exists a constant C independent of n such
that
k
nX
i=1
kpink21  C:
Proof. Writing (2.9) with the test function q = pin, multiplying by k and summing from i = 2
to any n0  n, we obtain
n0X
i=2
 
(c0 +B)(pin pi 1n );pin
 
+k
n0X
i=2
 
 (Bpin)rpin;rpin
 
= k
n0X
i=2
(gin;pin): (2.19)
Using summation by parts and the fact that B is self-adjoint, we have that
n0X
i=2
 
(c0 +B)(pin pi 1n );pin
 
= 12
n0X
i=2
h 
(c0 +B)(pin pi 1n );pin +pi 1n
 
+
 
(c0 +B)(pin pi 1n );pin pi 1n
 i
= 12
n0X
i=2
h 
(c0 +B)pin;pin
 
 
 
(c0 +B)pi 1n ;pi 1n
 
+
 
(c0 +B)(pin pi 1n );pin pi 1n
 i
= 12
n0X
i=2
 
(c0 +B)(pin pi 1n );pin pi 1n
 
+ 12
 
(c0 +B)pn0n ;pn0n
 
 12
 
(c0 +B)p1n;p1n
 
:
It follows from (1.11) and (2.19) that
1
2
n0X
i=2
 
(c0 +B)(pin pi 1n );pin pi 1n
 
+ 12
 
(c0 +B)pn0n ;pn0n
 
+k
n0X
i=2
k krpink2
 k
n0X
i=2
(gin;pin) + 12
 
(c0 +B)p1n;p1n
 
;
(2.20)
and the monotonicity of B leads to
k
n0X
i=2
k krpink2  k
n0X
i=2
(gin;pin) + 12
 
(c0 +B)p1n;p1n
 
: (2.21)
21
To estimate the second term on the right hand side of the above inequality, we write (2.9)
with i = 1, multiply by k and set the test function q = p1n to obtain
 
(c0 +B)p1n;p1n
 
+kk krp1nk2
 
 
(c0 +B)p1n;p1n
 
+k
 
 (Bp1n)rp1n;rp1n
 
=k(g1n;p1n) + (l;p1n):
The monotonicity of B, Young?s inequality, and inequality (2.21) yield
k
n0X
i=1
krpink2  Ck
n0X
i=1
(gin;pin) + (l;p1n)
 k
n0X
i=1
(C(")kgink2 +"kpink2) +C(")klk2;
(2.22)
where C(") is a constant depending on ". Choosing a su ciently small " > 0 such that
"C2p < 1 where Cp is given by (1.22) and noticing that for some constant C independent of
n
k
n0X
i=1
kgink2  C
nX
i=1
Z ti
ti 1
kg(t)k2dt = Ckgk2L2(I;L2( )); (2.23)
we obtain that for any n0  n
k
n0X
i=1
kpink21  C1 "C2
p
kgkL2(I;L2( )) :
Lemma 2.12 (see [1], p. 327). For ? in C1(I), de ne two piecewise constant functions ?n
and ~?n such that
8
><
>:
?n(t) = ?(ti); t2(ti 1;ti];
~?n(t) =
 
?(ti+1) ?(ti)
 
=k; t2(ti 1;ti];
22
for i = 1;2;:::;n, with
8>
<
>:
?n(0) = ?(t1);
~?n(0) = ~?(t1); ~?n(tn+1) = ~?(tn):
Then
k?n ?kL2(I)  C(?)k; k~?n ?0kL2(I)  C(?)k1=2
where the constant C(?) depends on ? only.
Lemma 2.13 (Aubin-Lions, see [57]). Let X0;X;X1 be three Banach spaces with X0  X 
X1. Suppose that X0 is compactly embedded in X and X is continuously embedded in X1.
Furthermore assume that X0 andd X2 are re exive spaces. Let 1 <p<1 and 1 <q <1.
De ne
W =fu2Lp([0;T];X0); u02Lq([0;T];X1g:
Then the embedding of W into Lp([0;T];X) is compact.
Remark 2.14. Speci cally, let X0 = H10 ( ), X = L2( ), X1 = H 1( ) the dual space of
H10 ( ) and p = q = 2. Then
W =fu2L2(I;H10 ( )); u02Lq(I;H 1( ))g
,!L2(I;L2( )) compactly:
(2.24)
Armed with the above Lemmas, we are now ready to show the existence of a solution
of equation (2.7).
Theorem 2.15. Given l in L2( ) and g in L2(I;L2( )), equation (2.7) has at least one
solution p in L2(I;H10 ( )).
Proof. For each positive integer n, de ne the piecewise constant function Pn2L2(I;H10 ( ))
as
Pn(t) = pin; t2(ti 1;ti]; (2.25)
23
for i = 1;2;:::;n with Pn(0) = p1n. From Lemma 2.11
kPnk2L2(I;H1
0( ))
= k
nX
i=1
kpink21  C: (2.26)
Hence there exists a weak limit of the sequence fPng1n=1 which is the candidate solution
of (2.7). However, the derivative of Pn is zero for a.e. t. To approximate the weak time
derivative we de ne ~Pn2L2(I;H10 ( )) as
~Pn(t) = pi 1n + (t ti 1)(pin pi 1n )=k; t2(ti 1;ti]; i = 1;2;:::;n: (2.27)
Obviously ~Pn is piecewise linear in t. Moreover, the time derivative of ~Pn is a piecewise
constant in t satisfying
~P0n =  pin = (pin pi 1n )=k; t2(ti 1;ti]; i = 1;2;:::;n: (2.28)
It is easy to verify that bothfPng1n=1 andf~Png1n=1 satisfy the same energy estimate and
have the same weak and strong limit. Thereforek~Pn(t)kL2(I;H10( )) is also uniformly bounded.
Due to the boundedness of B on H10 ( ), there exists a constant C such that
k(c0 +B) ~Pn(t)kL2(I;H10( ))  C: (2.29)
Meanwhile, it follows from (1.11), (2.9) and the Cauchy-Schwarz inequality that
  
 
 
(c0 +B) pin;q
   
 =j(gin;q) ( (Bpin)rpin;rq)j (kgink+k kpink1)kqk1 8q2H10 ( ):
Hence
k(c0 +B) pink 1  kgink+k kpink1:
24
It follows from Lemma 2.11 and (2.23) that
k(c0 +B) ~P0nk2L2(I;H 1( )) = k
nX
i=1
k(c0 +B) pink2 1
 k
nX
i=1
kgink2 +k
nX
i=1
k kpink21
 C:
(2.30)
Thus we can  nd a subsequence of f(c0 + B) ~Png1n=1, still denoted as f(c0 + B) ~Png1n=1, and
a function r in L2(I;H10 ( )) with r0 in L2(I;H 1( )) ([58], p. 356) such that
8
><
>:
(c0 +B) ~Pn *r in L2(I;H10 ( ));
(c0 +B) ~P0n *r0 in L2(I;H 1( )):
From the embedding (2.24)
(c0 +B) ~Pn!r in L2(I;L2( )):
Due to the relationship between the piecewise constant Pn and the piecewise linear functions
~Pn we have that
(c0 +B)Pn!r in L2(I;L2( )): (2.31)
Since c0+B is invertible on L2( ) (see Remark 2.8), there exists a p in L2(I;L2( )) satisfying
(c0 +B)p = q, such that
Pn!p = (c0 +B) 1r in L2(I;L2( )): (2.32)
Recall that fPng1n=1 is bounded in L2(I;H10 ( )) (see (2.26)). Hence there exists  p in
L2(I;H10 ( )) such that
Pn *  p in L2(I;H10 ( )): (2.33)
25
Then p =  p because the weak limit is unique.
Next we prove that p satis es (2.7). De ne
gn(t) = gin; t2(ti 1;ti]; i = 1;2;:::;n:
Proceeding as in Lemma 2.12, we obtain that
kgn gkL2(I;L2( )) !0; as n!1: (2.34)
Let  q be a test function of the form  q = q? satisfying
q2W1;1( )\H10 ( ) and ?(t)2C10 (I): (2.35)
It follows from summation by parts that
nX
i=1
 
(c0+B)(pin pi 1n );q
 
?(ti)
=
 
(c0 +B)pnn;q
 
?(tn) 
 
(c0 +B)p0n;q
 
?(t1)
+k
n 1X
i=1
 
(c0 +B)pin;q
  
?(ti+1) ?(ti)
 
=k:
26
Multiplying (2.9) with k?(ti), summing up from i = 1 to n, and using the summation above,
we obtain
 
(c0 +B)pnn;q
 
?(T) 
 
(c0 +B)p0n;q
 
?(t1)
+
n 1X
i=1
 
(c0 +B)pin;q
  
?(ti+1) ?(ti)
 
+k
nX
i=1
 
 (Bpin)rpin;rq
 
?(ti)
=k
nX
i=1
(gin;q)?(ti):
Recall that ?(T) = 0. The above equation takes the form
 (l;q)?(k) 
Z T
0
 
(c0 +B)Pn;q
 
~?n dt
+
Z T
0
 
 (BPn)rPn;rq
 
?n dt
=
Z T
0
(gn;q)?n:
(2.36)
Letting n!1 in (2.36) and using the continuity of ?(t), we have that
 (l;q)?(k)! (l;q)?(0) = 0: (2.37)
Notice that
Z T
0
 
(c0 +B)Pn;q
 
~?n dt =
Z T
0
 
(c0 +B)Pn;q
 
( ~?n ?0) dt
+
Z T
0
 
(c0 +B)Pn;q
 
?0 dt:
27
It follows from Lemma 2.12, (2.6), and (2.26) that
  
 
Z T
0
 
(c0 +B)Pn;q
 
( ~?n ?0) dt
  
 
 kqk
Z T
0
k(c0 +B)Pnkj~?n ?0jdt
 CkqkkPnkL2(I;L2( ))k~?n ?0kL2(I) !0 as n!1:
Since q?02L2(I;H10 ( )), the above and (2.33) yields that
Z T
0
 
(c0 +B)Pn;q
 
?0 dt
!
Z T
0
 
(c0 +B)p;q
 
?0 dt
= 
Z T
0
 
(c0 +B)p0;q
 
? dt:
Hence
Z T
0
 
(c0 +B)Pn;q
 
~?n dt! 
Z T
0
 
(c0 +B)p0;q
 
? dt; as n!1: (2.38)
Next we prove the convergence of the third term
Z T
0
 
 (BPn)rPn;rq
 
?n dt on the left
hand side of (2.36) to
Z T
0
 
 (Bp)rp;rq
 
? dt. Write
Z T
0
 
 (BPn)rPn;rq
 
?n dt 
Z T
0
 
 (Bp)rp;rq
 
? dt
=
Z T
0
 
 (BPn)rPn;r
 
(?n ?) dt
+
Z T
0
 
 (Bp)r(Pn p);rq
 
? dt
+
Z T
0
  
 (BPn)  (Bp) rPn;rq
 
? dt:
:= I + II + III
28
Using (1.11), Lemma 2.12, (2.26) and (4.11), we have that
I =
  
 
Z T
0
 
 (BPn)rPn;rq
  
?n(t) ?(t)
 
dt
  
 
 k krqk1kPnkL2(I;H10( ))k?n(t) ?(t)kL2(I) !0; as n!1:
The test function rq? belongs to L2(I;L2( )). Therefore the weak convergence of Pn to p
in L2(I;L2( ) and (1.11) imply that
II =
Z T
0
 
 (Bp)r(Pn p);rq
 
?(t) dt!0; as n!1:
For III, we  rst use (1.11) and (2.32) to deduce that
 (BPn)! (Bp) in L2(I;L2( )):
Thus
III =
  
 
Z T
0
  
 (BPn)  (Bp) rPn;rq
 
?(t) dt
  
 
 C
Z T
0
k (BPn)  (Bp)kkrPnkdt
 Ck (BPn)  (Bp)kL2(I;L2( ))kPnkL2(I;H10( ))
!0 as n!1:
Combining the above estimates we obtain
Z T
0
 
 (BPn)rPn;rq
 
?n dt!
Z T
0
 
 (Bp)rp;rq
 
? dt; as n!1: (2.39)
From a similar argument, it is straightforward to show that
Z T
0
(gn;q)?n dt =
Z T
0
(g;q)?n(t) dt!
Z T
0
(g;q)?(t) dt: (2.40)
29
Combing (2.37), (2.38), (2.39) and (2.40), we conclude that for any test function  q of the
form q?(t) satisfying (2.35),
 
Z T
0
 
(c0 +B)p0;  q
 
dt+
Z T
0
 
 (Bp)rp;r q
 
dt =
Z T
0
(g;  q) dt: (2.41)
The test functions  q de ned in (2.35) is dense in L2(I;H10 ( )). Therefore (2.41) holds for
all the test function q in L2(I;H10 ( )), i.e., p satis es (2.7) for a.e. t2I.
Finally we consider the initial condition. The facts
(c0 +B)p2L2(I;H10 ( )) and (c0 +B)p02L2(I;H 1( ))
imply that (c0 +B)p belongs to C(I;L2( )), the space of all continuous functions that value
in L2( ) (see [58], p. 288). Hence the initial condition (c0 +B)p( ;0) = l should be given in
L2( ).
In the next Theorem we give some conditions that guarantee the uniqueness of solutions
of equation (2.7).
Theorem 2.16. Assume
(H1)  is Lipschitz continuous with Lipschitz constant kl, i.e.,
j (x)  (y)j<kljx yj; 8x;y2R;
(H2) rp2L1( ) and there exists a constant C such that krpk1 C;
(H3) c0k =(C2pk CB1) > 1.
Then the solution p of equation (2.7) is unique.
30
Proof. Suppose p1;p2 both satisfy (2.7), i.e.,
h@@t(c0 +B)p1;qi+
 
 (Bp1)rp1;rq
 
= (g;q); 8q2H10 ( );
h@@t(c0 +B)p2;qi+
 
 (Bp2)rp2;rq
 
= (g;q); 8q2H10 ( ):
Taking the substraction and set q = (c0 +B)(p1 p2) we have
h(c0 +B)(p1 p2)0;(c0 +B)(p1 p2)i
+
 
 (Bp1)rp1  (Bp2)rp2;r(c0 +B)(p1 p2)
 
=h(c0 +B)(p1 p2)0;(c0 +B)(p1 p2)i
+
  
 (Bp1)  (Bp2) rp1;r(c0 +B)(p1 p2)
 
+
 
 (Bp2)r(p1 p2);rc0(p1 p2)
 
+
 
 (Bp2)r(p1 p2);rB(p1 p2)
 
= 0:
(2.42)
Recall
h(c0 +B)(p1 p2)0;(c0 +B)(p1 p2)i= 12 ddtk(c0 +B)(p1 p2)k2: (2.43)
Then (H1), (H2), (2.6), and Young?s inequality yields that
j(( (Bp1)  (Bp2))rp1;r(c0 +B)(p1 p2))j
 "k( (Bp1)  (Bp2)k2 +C(")k(c0 +B)(p1 p2)k21
 "kp1 p2k2 +C(")kp1 p2k21:
(2.44)
It follows from (1.11), (1.22) and (2.6) that
 
 (Bp2)r(p1 p2);rc0(p1 p2)
 
 (c0k =C2p)kp1 p2k21 (2.45)
31
and
j( (Bp2)r(p1 p2);rB(p1 p2))j k CB1kp1 p2k21: (2.46)
(2.43), (2.44), (2.45), (2.46), and (2.43) yield
1
2
d
dtk(c0 +B)(p1 p2)k
2 C(")kp1 p2k2
+ (c0k =C2p k CB1  ")kp1 p2k21
 0:
(2.47)
Hypothesis (H3) allows us to choose small " > 0 such that c0k =C2p k CB1  " > 0. The
boundedness of the inverse of B and (2.47) give that
1
2
d
dtk(c0 +B)(p1 p2)k
2  Ckp1 p2k2  Ck(c0 +B)(p1 p2)k2:
It follows from the Gronwall?s inequality that (c0 + B)(p1  p2) = 0. Therefore p1 = p2
because c0 +B has a bounded inverse.
Remark 2.17. After obtaining a solution p of equation (2.7), we can solve the  rst equation
of (2.5) for u with p substituted to the right hand side. According to the inverse estimate
of the elliptic equation
ku(t)k1  Ckp(t)k; 8t2I;
u belongs to L2(I; H10( )) and the pair (u;p) is the solution of system (2.5). Furthermore,
the linear dependence of u on p guarantees the uniqueness of u as long as p is unique.
32
2.3 Numerical approximation with the  nite element method
In this section, we consider the numerical approximation to the solutions of (2.7). In
particular we will derive error estimates for a fully discretized numerical scheme using back-
ward Euler method for the temporal discretization and  nite element method on the spatial
dimension. Throughout this section, we assume that the weak solution (u;p) of (1.10) is
unique, i.e., we assume that hypothesis (H1), (H2) and (H3) are satis ed.
We start by constructing the  nite element spaces as follows. Let  h be a family of
quasi-uniform triangulations (see [5], p. 2-3) of a convex polygonal, or polyhedral, domain
 satisfying max 2 
h
diam   h (where  is a geometrical element, e.g., a triangle, or a tetra-
hedron). Let Qh be the space consisting of continuous functions on  which are linear on
each triangle, or tetrahedron, and vanish on @ . Let fPjgNhj=1 be the set of all the interior
vertices of the triangulation. Assume that  pj is the pyramid function which value 1 at Pj
and vanishes at all the other vertices. It is easy to see thatf pjgNhj=1 forms a basis of Qh. Let
Vh = (Qh)d, d = 2;3. Then f( pj;0);(0; pj)gNhj=1 ( f( pj;0;0);(0; pj;0);(0;0; pj)gNhj=1 ) forms
a basis of Vh for d=2 or d=3, respectively.
As in the previous section, we denote the time stepsize by k, that is, k = TN, for
some positive integer N, and tn = nk, n = 0   ;N. Let uh 2 Vh and ph 2 Qh be
the approximation solution of (2.5). We write Pn = ph(tn) and Un = uh(tn) to be the
approximation of u, p at tn. De ne  @Pn := (Pn Pn 1)=k, n = 1   ;N. The fully
discretized  nite element approximation with backward Euler method is to  nd Un 2 Vh,
Pn2Qh, for n = 1;   ;N, such that
8
><
>:
e(Un;v) = (rPn;v); 8v2Vh;
(c0  @Pn;q) +  @(r Un;q) +
 
 (r Un)rPn;rq
 
=
 
g(tn);q
 
; 8q2Qh:
(2.48)
33
Here we set c0P0 +r U0 = l0, where l0 is the approximation to l in Qh.
Remark 2.18. We can prove the existence of a weak solution of (2.48) following an argument
similar to the one used for the continuous problem.
Next we consider the error estimates of the approximate solutions given by (2.48). We
need the Ritz-Galerkin type projections of the steady state problem corresponding to (2.5).
Given r2H10 ( ), de ne the projeciton Rhr of r onto Qh by
 
 (r u)r(r Rhr);rq
 
=
 
 (Bp)r(r Rhr);rq
 
= 0; 8q2Qh: (2.49)
We  rst introduce the following Lemmas which can be found in [27] and [5].
Lemma 2.19 ([27], Lemma 3.4.2). Let B be a continuous linear operator on a Banach space
X and let f : [0;T]!X be continuously di erentiable with respect to t. Then
B@f@t = @@tBf:
Lemma 2.20 ([5], p. 3). Let Qh be given as above. De ne the interpolation operator
Ih : H2( )\H10 ( )!Qh such that for any q2H2( )\H10 ( ). Then we have the estimate
kq Ihqk+hkr(q Ihq)k Ch2kqk2:
Lemma 2.21 ([5], Lemma 13.1). Assume that q2H2( )\H10 ( ). Then
kr(q Rhq)k C1hkqk2;
and
kq Rhqk C2h2kqk2;
34
where C1 and C2 depend on u and p.
To estimate the error between between pn and Pn, we de ne
 n := pn Rhpn;  n := Rhpn Pn: (2.50)
With these we can write the di erence between pn and Pn as
pn Pn =  n + n: (2.51)
Lemma 2.22. Assume that  is di erentiable, p2C1(I;H2( )\H10 ( )) and u2C1(I;W1;1( )).
Then
k (t)k+hkr (t)k C(u;p)h2; t2(0;T];
and
k t(t)k+hkr t(t)k C(u;p;pt)h2; t2(0;T]:
Proof. The  rst estimate follows directly from Lemma 2.21. Di erentiating (2.49) with
respect to t and setting r = p we obtain
 
 (Bp)r t;rq
 
+
  
 (Bp) tr ;rq
 
= 0; 8q2Qh: (2.52)
Hypothesis (H1) guarantees that  0 is uniformly bounded. Thus
  (Bp) 
t =
  (r u) 
t =  
0(r u)(r u)t
35
is also uniformly bounded due to the assumptions on the regularity of u. It follows from
(1.11), (2.49), (H1), and (2.52) that
k kr tk2  
 
 (Bp)r t;r t
 
=
 
 (Bp)r t;r(pt q)
 
+
 
 (Bp)r t;r(q Rhpt)
 
=
 
 (Bp)r t;r(pt q)
 
 
  
 (Bp) tr ;r(q Rhpt)
 
 C(kr tkkr(pt q)k+Ckr kkr(q Rhpt)k):
Taking w = Ihpt, applying Lemma 2.20 and Young?s inequality twice we obtain
k kr tk2  Chkptk2kr tk+Ckr k
 
kr(Ihpt pt)k+kr(pt Rhpt)k
 
 Chkptk2kr tk+Chkr kkptk2 +kr kkr tk
 "kr tk2 +h2C(")kptk22 +Ch2kptk22 +"kr tk2 +C(")kr k2:
Recall the  rst estimate that kr k Ch2. Choose "> 0 such that "<k . Then
kr tk2  C(h2kptk22 +kr k2) Ch2:
Next we estimate k tk using a standard duality argument. For ?2L2( ), let  2H10 ( )
be the solution of
 r 
 
 (Bp)r 
 
=  (Bp)   r (Bp) r = ?:
The existence of ? is guaranteed by the Lax-Milgram Lemma. Furthermore there exists a
constant C such that
k k2  Ck?k:
36
Using integration by parts and (2.52) we have that for any q2Qh
( t;?) =
 
 (Bp)r t;r 
 
=
 
 (Bp)r t;r(  q)
 
+
 
 (Bp)r t;rq
 
=
 
 (Bp)r t;r(  q)
 
 
 
( (Bp))tr ;rq
 
=
 
 (Bp)r t;r(  q)
 
+
  
 (Bp) tr ;r(  q)
 
 
 
r ;  (Bp) tr 
 
:
Integration by parts yields
 
 
r ;  (Bp) tr 
 
=
 
 ;  (Bp) t  
 
+
 
 ;r  ( (Bp))t r 
 
:
Recall that both  (Bp) and   (Bp) t are uniformly bounded. Choosing w = Ih and using
the previous estimates, we obtain
j( t;?)j C
 
kr tjkr(  Ih k+kr kkr(  Ih )k+k kk  k
 
 C
 
kr tkhk k2 +kr khk k2 +k kk  k
 
 C(p)h2k k2  C(p)h2k?k:
The statement above holds for any ?2L2( ). Therefore k t(t)k Ch2.
Lemma 2.23. krRhp(t)kL1( )  C(p) for any p(t)2H2( )\H10 ( ), 8t2[0;T].
Proof. For q2H2( )\H10 ( ), we have the inverse estimate(since rRhp(t) is constant on
teach triangle)
krqkL1( )  Ch 1krqk; for q2Qh;
37
which implies that
kr(Rhp Ihp)kL1( )  Ch 1kr(Rhp Ihp)k
 Ch 1
 
kr(Rhp p)k+kr(p Ihp)k
 
 Ch 1
 
C(p)h+Chkpk2
 
 C(p):
Using the fact that krIhpkL1( )  CkrpkL1( ) we obtain
krRhpkL1( )  kr(Rhp Ihp)kL1( ) +krIhpkL1( )  C(p):
Writing un = u(tn), pn = p(tn) and pnt = pt(tn), we are now ready to derive the error
estimates for the  nite element approximation (2.48).
Theorem 2.24. Assume
1. The hypothesis of Theorem 2.16 holds;
2.  is di erentiable;
3. p2C2(I;H2( )\H10 ( )) and u2C1(I;W1;1( )).
Then there exists k0 > 0 such that for k k0 , there exists constants C1 and C2 depending
on u, p, pt, ptt, and l, such that
ku(tn) Unk+kpn Pnk C1kl0 lk+C2(h2 +k):
Proof. We  rst estimate kp(tn) Pnk. In light of (2.51) and Lemma 2.22, it su ces to
estimate  n de ned in (2.50).
38
From (2.49)
 
 (Bpn)rpn  (BPn)rPn;rq
 
=
 
 (Bpn)rRhpn;rq
 
 
 
 (BPn)rRhpn;rq
 
+
 
 (BPn)rRhpn;rq
 
 
 
 (BPn)rPn;rq
 
; 8q2Qh:
Substracting (2.48) from the continuous equation (2.9) and using the above equation we have
that
0 =
 
(c0 +B)(pnt   @Pn);q
 
+
 
 (Bpn)rpn  (BPn)rPn;q
 
=
 
(c0 +B)(pnt   @pn);q
 
+
 
(c0 +B) @ n;q
 
+
 
(c0 +B) @ n;q
 
+
  
 (Bpn)  (BPn) rRhpn;rq
 
+
 
 (BPn)r n;rq
 
:
Letting q =  n in the above equation we obtain
 
(c0 +B) @ n; n
 
+
 
 (BPn)r n;r n
 
=
 
(c0 +B)(pnt   @pn); n
 
+
 
(c0 +B) @ n; n
 
+
  
 (Bpn)  (BPn) rRhpn;r n
 
:
(2.53)
The fact that B is self-adjoint and monotone leads to
 
(c0 +B) @ n; n
 
= 12
h 
(c0 +B) @ n; n  n 1
 
+
 
(c0 +B) @ n; n + n 1
 i
= k2
 
(c0 +B) @ n;  @ n
 
+ 12
h 
(c0 +B) n; n
 
 
 
(c0 +B) n 1; n 1
 i
 12  @
 
(c0 +B) n; n
 
:
39
It follows from (1.11), (2.53), (H1), Remark 2.8, Lemma 2.23, and Young?s inequality that
1
2
 @ (c0 +B) n; n +k kr nk2
 C(")k(c0 +B)(pnt   @pn)k2 +"k nk2
+C(")k(c0 +B) @ nk2 +"k nk2
+C(")k (BPn)  (Bpn)k2 +"kr nk2
 C(")kpnt   @pnk2 +C(")k @ nk2
+C(")kPn pnk2 + (2C2p + 1)"kr k2;
(2.54)
where Cp is given by (1.22). Choose "> 0 such that (2Cp + 1)" k . Then
1
2
 @ (c0 +B) n; n  C kpnt   @pnk2 +k @ nk2 +kPn pnk2 :
Multiplying the above inequality by k, we have that
 
(c0 +B) n; n
 
 
 
(c0 +B) n 1; n 1
 
 Ck
 
kpnt   @pnk2 +k @ nk2 +kPn pnk2
 
 Ck
 
kpnt   @pnk2 +k @ nk2 +k nk2 +k nk2 + (B n; n)
 
 Ck
 
kpnt   @pnk2 +k @ nk2 +k nk2
 
+Ck
 
(c0 +B) n; n
 
:
Equivalently ([5], p. 238)
(1 Ck)
 
(c0 +B) n; n
 
 
 
(c0 +B) n 1; n 1
 
+CkRn;
where the remainder Rn is given by
Rn =kpnt   @pnk2 +k @ nk2 +k nk2:
40
For su ciently small k, the above equation is equivalent to
 
(c0 +B) n; n
 
 (1 +Ck)
 
(c0 +B) n 1; n 1
 
+CkRn:
Inductively for n 2 we have that
 
(c0 +B) n; n
 
 (1 +Ck)n 1
 
(c0 +B) 1; 1
 
+Ck
nX
j=2
(1 +Ck)n jRj
 C
 
(c0 +B) 1; 1
 
+Ck
nX
j=2
Rj:
(2.55)
It follows from Lemma 2.22 that
k jk Ch2 (2.56)
and
k @ jk= 1kk
Z tj
tj 1
 t(s) dsk Ch2: (2.57)
From (H3)
kpjt  @pjk= 1kkkpjt p(tj) +p(tj 1)k
= 1kk
Z tj
tj 1
(s tj 1)ptt dsk
 k
Z tn
0
kpttkds Ck:
(2.58)
Combining (2.56), (2.57), and (2.58), we obatin
Rj C(h2 +k)2; 2 j n: (2.59)
41
For n = 1, denote  (l) = l Rhl, (c0 +B) 0 =  (l) = Rhl l0 Then (2.54) yields
1
k
 
(c0+B) 1  (l); 1
 
+C1kr 1k2
 C2
 
k(c0 +B)(p1t  @p1)k2 +k(c0 +B) @ 1k2 +k (BP1)  (Bp1)k2
 
:
Applying the Young?s inequality we have
 
(c0+B) 1; 1
 
+C1kr 1k2
 C2k
 
k(c0 +B)(p1t  @p1)k2 +k(c0 +B) @ 1k2
+kP1 p1k2
 
+ ( (l); 1):
 C2
 
k(c0 +B)(p1t  @p1)k2 +k(c0 +B) @ 1k2
+k 1k2 +k 1k2
 
+C(")k (l)k2 +"k 1k2:
Choose su ciently small "> 0 such that "<C1. Then the above inequality leads to
(1 Ck)
 
(c0 +B) 1; 1
 
 C
 
k (l)k2 +k(c0 +B)(p1t  @p1)k2 +k(c0 +B) @ 1k2 +k 1k2
 (2.60)
Proceeding as in (2.56), (2.57), (2.58), using Lemma 2.22, we have
k 1k Ch2:
As a consequence of the boundedness of B and Lemma 2.22
k @(c0 +B) 1k= 1kk
Z t1
0
(c0 +B) t(s) dsk Ch2
42
and
k(c0+B)(p1t  @p1)k
= 1kkk(c0 +B)p1t (c0 +B)p(k) + (c0 +B)p(0)k
= 1kk
Z t1
0
s(c0 +B)ptt dsk
 k
Z t1
0
k(c0 +B)pttkds = Ck:
Thus
((c0 +B) 1; 1) Ck (l)k2 +C(h2 +k): (2.61)
It follows from Lemma 2.22 that
k (l)k kl0 lk+k (l)k kl l0k+Ch2:
Hence (2.59) and (2.61) yield
k nk Ckl l0k+C(h2 +k):
Hence
kpn Pnk Ckl l0k+C(h2 +k):
To obtain the estimate for u(tn) Un, we de ne the projection Rh : H10( )!Vh by
e(u Rhu;v) = 0; 8v2Vh; u2H10( ): (2.62)
We split the error un Un as
un Un = un Rhun + Rhun Un =  nu + nu;
43
where  nu = un Rhun and  nu = Rhun Un. Similarly to Lemma 2.22, it is easy to verify
that the projection error  u satis es
k uk C(u)h2:
Substracting the  rst equation in (2.48) from (2.5), using (2.62), we have that
e( nu;v) = a
  
Un un ;v
 
= 
 
r(Pn pn);v
 
; 8v2Vh:
With v =  u(tn) in the above equation, Young?s inequality and the inverse estimate of elliptic
equation yield
k nuk1  kpn Pnk Ckl l0k+C(h2 +k)
from which the assertion in Theorem follows.
2.4 Numerical example
We  rst test the convergence rate in two dimension for brevity. Let  = [ 1;1] [ 1;1],
I = [0;1]. We choose Kozeny-Carmen-type ([59], [60], [61]) hydraulic conductivity  , which
is de ned by
 (s) =
8
>>>>
<
>>>>
:
k0
 
 3(s)
(1  (s))2 ;
 0
 0 1 <s <s<s
 < 1;
k ; s s ;
k ; s  s:
(2.63)
Here
 (s) =  0 + (1  0)s
and
k = k0 ( 0 + (1  0)s )
3
(1 s )2(1  0)2 ; k
 = k0
 
( 0 + (1  0)s )3
(1 s )2(1  0)2 ;
where k0,  ,  0, s , and s are given constants.
44
Table 2.1: The converge rate for a  xed time step in the L2(I;L2( )) norm
h kp phk ku uhk conv. rate p conv. rate u
1/4 0.1023 0.0295
1/8 0.0554 0.0072 0.88 2.03
1/16 0.0145 0.0018 1.93 2.00
1/32 0.0037 4.61e-4 1.97 1.97
1/64 9.20e-4 1.16e-4 2.01 1.99
1/128 2.30e-4 2.89e-5 2.00 2.00
It is easy to see that the continuous function  given by (2.63) satis es (1.11) and (H1).
We set  =  = c0 =  = 1 in system (1.10) and  0 =  = 1,  0 = 0:5, s = 0:75, s = 0:75
in our numerical experiment. Also we choose the forcing term f, g in (1.10) and the initial
value l in (2.5) such that (u;p) de ned by
u = e t(sin( x) sin( y);sin( x) sin( y))T
and
p = tsin( x) sin( y)
is the exact solution of (1.10). Notice that the solution (u;p) de ned above is analytic.
Therefore (H2) is satis ed.
For each time step, we solve the nonlinear problem using the Picard iteration. We set
the tolerance of this iteration to be 10 10 .
To verify the convergence rate in spatial dimension, we let the time step k be a  xed
value k = 10 4. Table 2.1 and Table 2.2 list the errors and rates of convergence in spatial
dimension. Figure 2.1 plots the errors in both L2(I;L2( )) norm and L2(I;H1( )) norm.
Both the tables and the  gure verify the second order convergence rate for the  nite element
approximation. Moreover, the H1 norms of u and p both have  rst order convergence, which
is natural for linear  nite elements with pyramid basis.
45
Table 2.2: The converge rate for a  xed time step in the L2(I;H10 ( )) norm
h kp phk1 ku uhk1 conv. rate p conv. rate u
1/4 0.7685 0.3061
1/8 0.5556 0.1646 0.47 0.90
1/16 0.2767 0.0838 1.01 0.97
1/32 0.1385 0.0421 1.00 0.99
1/64 0.0693 0.0211 1.00 1.00
1/128 0.0346 0.0105 1.00 1.00
Table 2.3: The converge rate for time step k = h2 in the L2(I;L2( )) norm
h k kp phk ku uhk conv. rate p conv. rate u
1/4 1/16 0.0980 0.0290
1/8 1/64 0.0519 0.0070 0.92 2.05
1/16 1/256 0.0135 0.0018 1.94 1.96
1/32 1/1024 0.0034 4.52e-4 1.99 1.99
Table 2.4: The converge rate for time step k = h2 in the L2(I;H1( )) norm
h k kp phk1 ku uhk1 conv. rate p conv. rate u
1/4 1/16 0.7530 0.2041
1/8 1/64 0.5347 0.1629 0.49 0.33
1/16 1/256 0.2662 0.0829 1.07 0.97
1/32 1/1024 0.1332 0.0416 1.00 0.99
Figure 2.1: The plot of errors for a  xed time step
?5 ?4.5?4 ?3.5?3 ?2.5?2 ?1.5?1?12
?10
?8
?6
?4
?2
0
ln(h)
ln(error)
convergence rate
pL2(L2)
uL2(L2)
pL2(H1)
uL2(H1)
46
Figure 2.2: The plot of errors for time step k = h2
?3.5 ?3 ?2.5 ?2 ?1.5 ?1?8
?7
?6
?5
?4
?3
?2
?1
0
ln(h)
ln(error)
convergence rate
pL2(L2)
uL2(L2)
pL2(H1)
uL2(H1)
To see the convergence rate with respect to the time step k, we set k = h2. The errors
are presented in table 2.3, table 2.4 and  gure 2.2, from which we know that when k is small,
the convergence rate with respect to k is the same as h2.
We next show some simulation of a real life example. Our intention is to show the e ect
of a nonlinear hydraulic conductivity  thus we choose a model in three dimensions and
we only consider the numerical approximation of the steady system to save computations.
Suppose we have a 2 2 2 inches cube sponge on  = [ 1;1] [ 1;1] [ 1;1]. We assume
that there is no external force of  uid resources, i.e., f = 0 and g = 0. We assume zero
boundary conditions of both u and p except that p = 1 on x =  1 (the left surface) and
p = 0 on x = 1 (the right surface). In this case, there is a pressure drop from the left surface
to the right surface. If  is a constant, we expect a linear pressure with constant pressure
drop. As a result, we have uniform  ux pointing to the right, parallel to the x axis.
Case 1. We want to show the improvements of simulation when we introduce the
nonlinear hydraulic conductivity  in (2.63). The  ux and hydraulic conductivity are then
given in  gure 2.3. From the graph we can see that the  ux is not uniform due to the
nonlinear  and the value of  varies from 0:3 to 0:8 at di erent part of the elastic body.
The highest  occurs close to the left surface since the pore structure has been pushed to
the right by the high pressure. As a result, the size of the pores become bigger which leads
47
Figure 2.3: Case 1
?1 ?0.5
0 0.5
1
?1?0.5
00.5
1?1
?0.5
0
0.5
1
xy
z
(a) Flux (b) Hydraulic conductivity
to a higher  . The lower  near the right surface where most of solid matrix stays is due to
the small size of the pores.
Case 2. We apply a displacement of compression in the middle of the front and back
surface, that is
u =
8
>>>>
<
>>>
>:
(0;0:1;0); y = 1; 0:25 x;z 0:25;
(0; 0:1;0); y = 1; 0:25 x;z 0:25;
0; otherwise;
and we assume the same boundary condition of p as in case 1. The  ux and  are shown
in  gure 2.4. We observe two facts. First, we have a very low hydraulic conductivity in
the middle and the  ux are avoiding where we pressed the body. Only very small  ux can
pass through the middle and this e ect becomes weaker as the  ux are away from the front
surface into the heart of the object. In real life case, this is because the pores close to the
front and back surface become small when we press the sponge. Another fact is that the
highest  ow occurs right around where we pressed the sponge. This is because of the elastic
nature of the sponge that the water originally kept in the pores near the middle of the front
and back surface was forced to come out due to pressing.
48
Figure 2.4: Case 2
?1 ?0.5 0 0.5 1 1.5?1
?0.8
?0.6
?0.4
?0.2
0
0.2
0.4
0.6
0.8
1
(a) Flux on y = 1 (b) Hydraulic conductivity
Figure 2.5: Case 3
?1 ?0.5 0 0.5 1 1.5?1
?0.8
?0.6
?0.4
?0.2
0
0.2
0.4
0.6
0.8
1
(a) Flux on y = 1 (b) Hydraulic conductivity
Case 3 We observe another case similar to case 2. We squeeze the lower half of the
sponge, that is
u =
8
>>>
><
>>>
>:
(0;0:1;0); y = 1; 1 z 0;
(0; 0:1;0); y = 1; 1 z 0;
0; otherwise:
The result is shown in  gure 2.5. Again, we observe a low hydraulic conductivity in the
lower half of the sponge and a high hydraulic conductivity around where we squeeze it.
49
2.5 Conclusion
We studied the mathematical model of  ows through elastic porous media. The PDE
system under consideration is comprised of the linear elasticity equation for the deforma-
tion of the matrix and an implicit parabolic equation for the  uid pressure with hydraulic
conductivity dependent on the displacement. We showed the existence of a weak solution
under minimum regularity assumption on the input data and its uniqueness under stronger
conditions. We then focused our attention on the  nite element approximations of the weak
solution and derived rigorous error estimates. We used our model to simulate the distribu-
tion of  ux and hydraulic conductivity inside a sponge. In the simulation we showed that
the nonlinear hydraulic conductivity  de ned in (2.63) captures the distribution of the hy-
draulic conductivity at di erent part of the elastic body, which cannot be obtained by the
linear system with a constant hydraulic conductivity assumption. The di erence is signi -
cant and has a huge e ect on the  ow through the porous media. We established a better
approximation to simulate the  ow through a porous media, which matches the common
sense.
50
Chapter 3
Steady bioconvection
3.1 Existence and uniqueness of a weak solution
The weak formulation. we consider a steady case of system (1.19). Choose the
pressure space L20( ) = L2( )=R to be the quotient subspace of L2( ), the subspace of
functions which are orthogonal to constants. It is easy to see that ([62], p. 23)
L20( ) =fp2L2( );
Z
 
p dx = 0g:
De ne the following bilinear and trilinear forms
8
>>>
>>>
>>>>
<
>>>
>>>
>>>
>:
a(c;r) = (rc;rr) 8c;r2H1( );
B2(u;v;w) =
Z
 
u rv w dx 8u;v;w2H10( );
B3(u;c;r) =
Z
 
u rc rdx; 8u2H10( ); c;r2H1( );
b(q;v) = (q;r v) 8q2L20( ); v2H10( )
(3.1)
and set
~H = H1( )\L20( ) =fc2H1( ) :
Z
 
c dx = 0g:
Condition (1.18) is equivalent to requiring c  j j2 ~H. For brevity, we write  =  j j and
we assume  m = 1, otherwise we rescale the problem. Also for brevity, de ne an auxiliary
concentration c = c  and f = f  g m  i3. Multiplying the equations of (1.19) by
51
test functions and integrating by parts we obtain the following weak formulation for system
(1.19) (without confussion, we write c = c and f = f ).
De nition 3.1. Given f in L2( ), a triple (u;p;c) 2 H10( ) L20( ) ~H is said to be a
weak solution of the steady bioconvection if it satis es
8
>>>
>>>
>>><
>>>
>>>>
>>:
 
 (c+ )ru;rv
 
+B2(u;u;v) +b(p;v)
= 
 
g(1 + c)i3;v
 
+ (f;v) 8v2H10( );
b(q;u) = 0 8q2L20( );
 a(c;r) +B3(u;c;r) U(c; @r@x
3
) = U ( @r@x
3
;1) 8r2 ~H:
(3.2)
We adopt the theory for solving Navier-Stokes type equations (see Chapter 4 of [62] and
Chapter 2-3 of [52]). De ne
V =fu2H10( ) : r u = 0 in  g:
Then to solve system (3.2), it su ces to solve the following auxiliary problem ([62]).
Find a pair (u;c)2V ~H such that
8>
><
>>:
 
 (c+ )ru;rv
 
+B2(u;u;v) = 
 
g(1 + c)i3 + f;v
 
8v2V;
 a(c;r) +B3(u;c;r) U(c; @r@x
3
) = U ( @r@x
3
;1) 8r2 ~H:
(3.3)
Remark 3.2. Obviously, if (u;c;p) is a solution of system (3.2), then (u;c) must be a solution
of (3.3). The converse is also true since the bilinear form b( ; ) de ned above saties es the
inf-sup condition (see [62], p. 113), i.e., for some  > 0
sup
v2H10( )
(q;r v)
kvk1   kqk 8q2L
2
0( ):
52
Existence. To prove the existence of a weak solution of (3.3), we construct a sequence
of approximate weak solutions using the Galerkin method. This will also be helpful in our
later discussion and analysis of the  nite element method. The following inequalities can be
found in [32]: 8
><
>:
kvk1  C krvk 8v2H10( );
krk1  C krrk 8r2 ~H;
(3.4)
for some constant C independent of v and r. The  rst inequality is Poincar e?s inequality
while the second is due to the fact that R rdx = 0. We can obviously use the same constant
C in both inequalities. From Lemma 2.2, it is straightforward that  (x) satisfying (1.21) is
a Nemytskii operator.
Next we study the properties of the trilinear forms B2 and B3.
Lemma 3.3. The trilinear form B2( ; ; ) and B3( ; ; ) are continuous on H1( ).
Proof. From Holder?s inequality and the Sobolev imbedding theorem in three dimensions we
have
B2(u;v;w) CkukL4( )krvkL2( )kwkL4( )  CB2kuk1kvk1kwk1: (3.5)
Similarly
B(u;c;r) CB3kuk1kck1krk1: (3.6)
Lemma 3.4. Assume that v;w2H10( ), c;r2 ~H, and u2V, then
8
><
>:
B2(u;v;w) +B2(u;w;v) = 0;
B3(u;c;r) +B3(u;r;c) = 0;
(3.7)
and
B2(u;v;v) = 0; B3(u;r;r) = 0: (3.8)
53
Proof. Notice that properties (3.7) and (3.8) are equivalent. Thus it su ces to prove (3.8).
According to Green?s formula, if u2V
B2(u;v;v) = 12
3X
i;j=1
Z
 
uj@(v
2
i)
@xj = 
1
2
3X
i=1
Z
 
r u v2i dx = 0
and
B3(u;r;r) = 12
3X
j=1
Z
 
uj@(r
2)
@xj = 
1
2
Z
 
r u r2 dx = 0:
Since V and ~H are both separable Hilbert spaces, there exist sequences fvjg1j=1 and
frjg1j=1 such that fvjg1j=1 and frjg1j=1 are orthonormal bases of V and ~H, respectively. Let
Vm and Hm be the spaces spanned byfv1;v2;:::;vmgandfr1;r2;:::;rmgrespectively. We
seek (um;cm)2Vm Hm such that
8
>><
>>:
 
 (cm + )rum;rv
 
+B2(um;um;v) = 
 
(g + cm)i3;v
 
+ (f;v) 8v2Vm;
 a(cm;r) +B3(um;cm;r) U(cm; @r@x
3
) = U ( @r@x
3
;1) 8r2Hm:
(3.9)
Proceeding as in Lemma 2.10, using properties (3.5) and (3.6), we can prove the existence of
a weak solution of (3.9) for any integar m> 0 either by a direct corollary of Brouwer  xed
point theorem or using Riesz? theorem.
Next we show that fumg1m=1 and fcmg1m=1 are uniformly bounded in V and ~H respec-
tively.
Lemma 3.5. Assume that
 
C2 >U: (3.10)
Then there exists a constant C independent of m such that
kcmk1 +kumk1 <C:
54
Proof. Let v = um, r = cm in (3.9). Equation (3.8) implies that
8
>><
>>:
 
 (cm + )rum;rum
 
= 
 
(g + cm)i3;um
 
+ (f;um);
 a(cm;cm) U(cm;@c
m
@x3 ) = U (
@cm
@x3 ;1):
Thus it follows from (1.21) and Young?s inequality that
  krumk2  
 
 (cm + )rum;rum
 
 j 
 
(g + cm)i3;um
 
+ (f;um)j
 
 
 kcmk+kf gi3k
 
kumk
(3.11)
and
 krcmk2   a(cm;cm)
 
  
 U(cm;@c
m
@x3 ) +U (
@cm
@x3 ;1)
  
 
 Ukcmk21 +U j j12kcmk1:
(3.12)
Using the above inequality, (3.4), and assumption (3.10) we obtain
kcmk1  
  
C2  U
  1
U j j12 : (3.13)
Substituting (3.13) into (3.11) gives
kumk1  C
2
 
  
 
(  C2
 
 U) 1 U j j12 ) +kf gi3k
 
:
We are now ready to show the existence of a solution of (3.3).
Theorem 3.6. Assume that (1.21) and (3.10) hold, and f 2L2( ), then system (3.3) has
a weak solution.
55
Proof. Consider sequences of solutions fumg1m=1 and fcmg1m=1 de ned by (3.9). According
to Lemma 3.5, both sequences are bounded therefore there exist u 2 V and c 2 ~H (via
subsequences if necessary) such that
um * u in V and cm *c in ~H as m!1: (3.14)
Then the Sobolev embedding theorem implies that
um!u in L2( ) and cm!c in L2( ) as m!1: (3.15)
We now show that the weak limit (u;c) is a solution of (3.3). Let v and r be test functions
such that
v2V\(C10 ( ))3; r2C1( )\ ~H: (3.16)
First notice that
 
 (c+ )ru;rv
 
 
 
 (cm + )rum;rv
 
=
 
 (c+ )r(u um);rv
 
+
  
 (c+ )  (cm + ) rum;rv
 
:= I + II:
From (1.21) and (3.14) it follows that
jIj   j(r(um u);rv)j!0 as m!1:
The limits (3.15) and the fact that  is a Nemytskii operator imply that
 (cm + )! (c+ ) in L2( ) as m!1: (3.17)
56
Hence from (3.16), Lemma 3.5, and Holder?s inequality
jIIj Ck (cm + )  (c+ )kkumk1 !0 as m!1:
Combing the above estimates we obtain
 
 (cm + )rum;rv
 
!
 
 (c+ )ru;rv
 
as m!1: (3.18)
Next using Green?s formula
B2(um;um;v) =
3X
i;j=1
Z
 
umj (@u
m
i
@xj vi) dx = 
3X
i;j=1
Z
 
umi umj (@vi@x
j
) dx:
By assumption (3.16), @v@xj is uniformly bounded. The fact um ! u in L2( ) implies that
umi umj !uiuj in L1( ) as m!1. Therefore
limm!1B2(um;um;v) = 
3X
i;j=1
Z
 
uiuj(@vi@x
j
) dx = B2(u;v;u) = B2(u;u;v):
Following the same argument, we have
B2(um;um;v)!B2(u;u;v); B3(um;cm;r)!B3(u;c;r) as m!1: (3.19)
It follows from the weak convergence of fcmg1n=1 to c that
8
>>>
>><
>>>>
>:
 
g(1 + cm)i3;v
 
!
 
g(1 + c)i3;v
 
as m!1;
 a(cm;r)! a(c;r) as m!1;
U(cm; @r@x
3
)!U(c; @r@x
3
) as m!1:
(3.20)
Because the test functions v, r de ned in (3.16) are dense in V and ~H respectively, conclu-
sions (3.18), (3.19), and (3.20) hold for 8v 2V and 8r2 ~H. Letting m!1 in (3.9) and
57
using the above results we obtain
8>
><
>>:
 
 (c+ )ru;rv
 
+B2(u;u;v) = 
 
(g + c)i3;v
 
+ (f;v) 8v2Vm;
 a(c;r) +B3(u;c;r) U(c; @r@x
3
) = U ( @r@x
3
;1) 8r2Hm:
Since the basis fvjg1j=1 and frjg1j=1 are dense in V and ~H, we conclude that (u;c) is a
solution of (3.3).
Uniqueness. First notice that the bilinear form b( ; ) satis es the inf-sup condition
(see Remark 3.2). Therefore for each solution (u;c)2V ~H of system (3.3), there exists a
unique p2L20( ) satisfying system (3.2) (see [62], p. 113). To prove the uniqueness of the
solution of (3.2), it su ces to prove that system (3.3) has a unique solution.
Following the proof of Lemma 3.5 we can obtain a priori estimates for u and c.
kuk1  C3 and kck1  C4: (3.21)
where
C3 = C
2
 
  ( C4 +kf gi3k); C4 =
U j j12
 
C2  U
:
Theorem 3.7. Assume that
(H4) The hypothesis of Theorem 3.6 holds;
(H5) The viscosity  ( ) is Lipschitz continuous, i.e., there exists a constant  L > 0 such that
j (x1)  (x2)j  Ljx1 x2j 8x1;x2 2R;
(H6) There exists a constant C0 such that kruk1 C0;
58
(H7) The following inequality
  
C2  
 
CB3C4
 
C2  U
( LC0 +g ) +CB2C3
!
> 0
holds.
Then the solution (u;c) of system (3.3) is unique.
Proof. Let (u;c) and ( u; c) be two di erent solutions of (3.3). Substituting both solutions
into (3.3) with v = u  u and r = c  c, taking the the di erence of the two equations, we
have
 
 (c+ )ru;r(u  u)
 
 
 
 ( c+ )r u;r(u  u)
 
+B2(u;u;u  u) B2( u;  u;u  u)
= g 
 
(c  c)i3;u  u
 
(3.22)
and
 a(c  c;c  c) +B3(u;c;c  c) B3( u; c;c  c) U(c  c;@(c  c)@x
3
) = 0: (3.23)
The skew symmetry (3.8) leads to the identity
8
><
>:
B2(u;u;u  u) B2( u;  u;u  u) = B2(u  u;u;u  u);
B3(u;c;c  c) B3( u; c;c  c) = B3(u  u;c;c  c):
(3.24)
Thus it follows from (3.6), (3.4), (3.21), (3.23), and (3.24) that
 
C2 kc  ck
2
1  jB3(u  u;c;c  c)j+U(c  c;
@(c  c)
@x3 )
 CB3C4jkc  ck1ku  uk1 +Ukc  ck21:
59
Then from (3.10)
kc  ck1  CB3C4 
C2  U
ku  uk1: (3.25)
Substituting the above estimate into (3.22) and combining (1.21), (3.4), (3.5), (3.21), (H6)
and (3.24), we obatin
  
C2 ku  uk
2
1  
 
 (c+ )r(u  u);r(u  u)
 
 
  
 
  
 (c+ )  ( c+ ) r u ;r(u  u)
   
 +jB2(u  u;u;u  u)j+g 
   (c  c)i
3;u  u
   
  LC0kc  ck1ku  uk1 +CB2ku  uk21kuk1 +g kc  ck1ku  uk1
 
 
CB3C4
 
C2  U
( LC0 +g ) +CB2C3
!
ku  uk21:
From assumption (H7) we conclude that
ku  uk1 =kc  ck1 = 0:
Remark 3.8. In practice, it is necessary to verify condition (3.10) and (H6). Since the micro-
organism is slightly denser than water,  =  0= m 1 is small. Therefore to ful ll (3.10)
and (H6),   and  must be su ciently large while U and C are su ciently small. In other
words, the model is valid for a suspension containing culture  uid with large viscosity, large
di usion rate, slowly upswimming micro-organisms in a relatively small container.
3.2 Numerical approximation using the  nite element method
In this section, we construct and analyze a  nite element method for approximating
weak solutions of (3.2). Throughout this section, we assume that the hypotheses of Theorem
3.7 hold.
60
Let  h be a regular triangulation of  and Xh, Mh, and Sh be  nite element subspaces
of H10( ), L20( ), and ~H, respectively. Assume that the following discrete inf-sup condition
holds.
sup
v2Xh
b(v;q)
kvkXh   kqkMh 8q2Mh; (3.26)
where  > 0 is a  xed constant. Furthermore we assume that Xh, Mh, and Sh satisfy the
following approximation properties.
inf
vh2Xh
kv vhk1  Chskvks+1 8v2Hs+1( ); 0 <s k; (3.27)
inf
qh2Mh
kq qhk Chskqks 8q2Hs( ); 0 <s k; (3.28)
inf
th2Sh
kt thk1  Chsktks+1 8t2Hs+1( ); 0 <s k; (3.29)
for k = 2, 3. For the construction of these spaces, see [63, 62, 52]. Next we de ne the discrete
divergence free space
Vh =fv2Xh : (r v;qh) = 0 8qh2Mhg:
Notice that in general, Vh is not a subspace of V. Thus in general the identity (3.8) does
not hold. To obtain a skew symmetry similar to (3.8) on Vh, we de ne auxiliary forms ^B2
and ^B3 by 8
><
>:
^B2(u;v;w) = 1
2B2(u;v;w) 
1
2B2(u;w;v);
^B3(u;c;r) = 1
2B3(u;c;r) 
1
2B3(u;r;c):
It is easy to verify that
^B2(u;v;w) = B2(u;w;v) 8u2V; 8v;w2H10( );
^B3(u;c;r) = B3(u;c;r) 8u2V; 8c;r2 ~H:
(3.30)
61
In addition, we have the identity
^B2(u;v;v) = 0; ^B3(u;c;c) = 0 8u;v2H10( ); 8c2 ~H; (3.31)
and the tricontinuous property
8
><
>:
^B2(u;v;w) CB
2kuk1kvk1kwk1 8u;v;w2H
1
0( );
^B3(u;c;r) CB
3kuk1kck1krk1 8u2H
1
0( ); 8c;r2 ~H;
(3.32)
where CB2 and CB3 are the same as in (3.5) and (3.6).
We de ne the  nite element approximation of (3.3) as follows.
De nition 3.9. Find (uh;ph;ch)2Xh Mh Sh, such that
8
>>>
>>>
>>><
>>>
>>>>
>>:
 
 (ch + )ruh;rv
 
+ ^B2(uh;uh;v) (ph;r v)
= 
 
g(1 + ch
 
i3;v) + (f;v) 8v2Xh;
(r uh;qh) = 0 8qh2Mh;
 a(ch;r) + ^B3(uh;ch;r) U(ch; @r@x
3
) = U ( @r@x
3
;1) 8r2Sh:
(3.33)
Analogous to the continuous case, we  rst solve an auxiliary discrete system.
Find (uh;ch)2Vh Sh such that
8
>>>
>><
>>>
>>:
 
 (ch + )ruh;rv
 
+ ^B2(uh;uh;v)
= 
 
g(1 + ch)i3;v
 
+ (f;v); 8v2Vh;
 a(ch;r) + ^B3(uh;ch;r) U(ch; @r@x
3
) = U ( @r@x
3
;1); 8r2Sh:
(3.34)
The skew symmetry (3.31) and inequality (3.32) guarantee the existence of a weak solution
of (3.34) by using anargument similar to the one used in the continuous case. A solution
62
(uh;ph;ch) of (3.33) can be completed by solving (see [62], p. 59, Theorem I.4.1)
(ph;r v) =
 
 (ch + )ruh;rv
 
+ ^B2(uh;unh;v)
+
 
g(1 + ch)i2;v) (fn;v) 8v2Xh:
(3.35)
The right hand side de nes a functional on Xh which vanishes on Vh. Due to (4.55) and the
property of Lagrange multipliers, equation (4.66) is always solvable and the solution ph2Mh
is unique in the quotient space Mh=Nh where
Nh =fqh2Sh : (qh;r v) = 0; 8v2Xhg:
Following the same approach as in Lemma 3.5, we can show that kuhk1 and kchk1 are
uniformly bounded, i.e., there exist constants C3 and C4 independent of h such that
kuhk1  C3; kchk1  C4: (3.36)
To carry out the error estimate, we introduce the Ritz Galerkin projectionsrh : H10( )!
Vh, sh : ~H!Sh , and the L2 projection  h : L20( )!Mh. In this way we split the errors
into two parts 8
>>>
>><
>>>>
>:
u uh = u rhu +rhu uh :=  hu + hu;
p ph = p  hp+ hp ph :=  hp + hp ;
c ch = c shc+shc ch :=  hc + hc :
(3.37)
From the approximation properties (3.27) { (3.29) we known that
8
>>>
><
>>>>
:
krhuk1  C(u); k huk1  Chskvks+1; u2Hs+1( ); 0 <s k;
kshck1  C(c); k hck1  Chskcks+1; c2Hs+1( ); 0 <s k;
k hpk C(p); k hpk Chskpks; p2Hs( ); 0 <s k:
(3.38)
63
Theorem 3.10. Assume that the hypotheses of Theorem 3.6 and Theorem 3.7 hold. Then
for u2Hs+1( ), p2Hs( ), and c2Hs+1( ), there exists a constant C independent of h
such that
ku uhk1 +kc chk1 +kp phk Chs; 0 <s k: (3.39)
Proof. Due to (3.38), it su ces to estimate  hu,  hp, and  hc. Subtracting (3.33) from (3.2)
with v =  hu, r =  hc and using (3.30) we have that
 
 (c+ )ru;r hu
 
 
 
 (ch + )ruh;r hu
 
+ ^B2(u;u; hu) ^B2(uh;uh; hu)
+b(p ph; hu)
= g 
 
(c ch)i3; hu
 
(3.40)
and
 a(c ch; hc) + ^B3(u;c; hc) ^B3(uh;ch; hc) U(c ch;@ 
h
c
@x3 ) = 0: (3.41)
It follows from (1.21), (3.31), and (3.32) that
 
C2 k 
h
ck
2
1   a( 
h
c; 
h
c)
=  a( hc; hc) ^B3( hu;ch; hc) ^B3(rhu; hc; hc)
 ^B3( hu;c; hc) +U( hc;@ 
h
c
@x3 ) +U( 
h
c;
@ hc
@x3 )
  k hck1k hck1 +Uk hck21 +Uk hck1k hck1
+CB3k hck1
 
k huk1kchk1 +krhuk1k hck1 +k huk1kck1
 
:
64
Due to assumption (3.10), moving the term Uk hck21 to the left and dividing by k hck1, we
have from (3.21) and (3.36) that
k hck1  1 
C2  U
 
( +U +CB3krhuk1)k hck1 +CB3C4(k huk1 +k huk1)
 
: (3.42)
Notice that  hu2Vh and  hp 2Mh. The de nition of Vh implies that b( hp; hu) = 0. Therefore
b(p ph; hu) = b( hp; hu):
Then according to (1.21), (H5), (H6), (3.21), (3.32), and (3.36), equation (3.40) yields
  
C2 k 
h
uk
2
1  
 
 (ch + )r hu;r hu
 
=
 
 (ch + )r hu;r hu
 
+
  
 (ch + )  (c+ ) ru;r hu
 
 ^B2( hu;uh; hu) B2(rhu; hu; hu) ^B2( hu;u; hu)
+b( hp; hu) +g 
 
(c ch)i3; hu
 
   k huk1k huk1 + ( LC0 +g )kch ck1k huk1
+CB2k huk1
 
C3k huk1 +krhuk1k huk1 +C3k huk1
 
+k hpkk huk1:
(3.43)
From (3.42)
kch ck1  k hck1 +k hck1
 k hck1 + 1 
C2  U
 
( +U +CB3krhuk1)k hck1 +CB3C4(k huk1 +k huk1)
 
:
65
Dividing (3.43) by k huk1 we obtain
(   C2
 
 CB2C3 ( LC0 +g ) CB3C4 
C2  U
)k huk1
 
 
  +CB2krhuk1 +CB2C3 + CB3C4( LC0 +g ) 
C2  U
 
k huk1
+ ( LC0 +g )
 
1 +  +U +CB3krhuk1 
C2  U
!
k hck1 +k hpk:
From assumption (H7) and (3.38)
k huk1  C
 
k huk1 +k hck1 +k hpk
 
; (3.44)
and from (3.42) and (3.38)
k hck1  C
 
k huk1 +k hck1 +k hpk
 
: (3.45)
Combining (4.67), (3.44) and (3.45) we have that
ku uhk1 +kc chk1  C
 
k huk1 +k hck1 +k hpk
 
: (3.46)
It remains to estimate kp phk. Subtracting (3.33) from (3.2) gives
 b(v; hp) =
 
 (c+ )ru;rv
 
 
 
 (ch + )ruh;rv
 
+ ^B2(u;u;v) ^B2(uh;uh;v) +b( hp;v) +g 
 
(c ch)i3;v
 
:
66
From (3.26), (3.32), (3.21), (3.36), and (3.46) it follows that
k hpk 1 sup
v2Xh
1
kvk1
 
( (ch + )r(u uh);rv
 
+
  
 (c+ )  (ch + ) ru;rv
 
+ ^B2(u uh;u;v) + ^B2(uh;u uh;v)
+b( hp;v) +g  (c ch)i3;v 
 
 1 sup
v2Xh
1
kvk1
 
(  + 2CB2C3)ku uhk1kvk1
+ ( LC0 + )kc chk1kvk1 +k hpkkvk1
 
 C
 
ku uhk1 +kc chk1 +k hpk
 
:
Combing the above estimate with (3.44) and (3.45), we obtain
kp phk1  C
 
k huk1 +k hck1 +k hpk
 
: (3.47)
Then the estimate of the pressure error follows from (3.38).
3.3 Numerical experiments
In this section we describe  ve numerical experiments that were conducted. The  rst
one uses arti cial data to verify the error estimates while the other four use data obtained
from lab experiments. We used Taylor-Hood  nite element spaces ([53]) for Vh and Mh,
and continuous piecewise quadratic function spaces for Sh. In particular, the theoretical
convergence rate is given by
k huk= O(hm); k huk1 = O(hm 1); k hpk= O(hm 1); k hck= O(hm);
67
Table 3.1: Convergence rate
h kp phk ku uhk kc chk kp phk1 ku uhk1 kc chk1
1/2 0.2520 0.0078 0.0049 0.9846 0.0854 0.0460
1/4 0.0323 0.0010 6.8E-04 0.3847 0.0207 0.0118
1/8 0.0055 1.31E-04 8.88E-05 0.1786 0.0050 0.030
1/16 0.0011 1.65E-05 1.13E-05 0.0877 0.0012 7.45E-04
1/32 2.28E-04 2.07E-06 1.42E-06 0.0436 3.09E-04 1.86E-04
1/64 4.90E-05 2.60E-07 1.80E-07 0.0216 7.68E-05 4.7E-05
conv. rate 2.22 2.99 2.97 1.02 2.01 2.00
where
u2Hm( )\H10( ); p2Hm 1( )\L20( ); c2Hm( )\L20( ); m = 2;3:
In the numerical experiment we used the following parameters
 = 0:1; U = 0:1;  = 1;
and
 (x) = sin2x+ 1; x2R:
The forcing term f was chosen so that the exact solution is
8
>>>>
><
>>>
>>:
u = (sin xsin y;sin xsin y)T ;
p = sin xsin y;
c = sin xsin y:
The numerical errors for di erent mesh sizes are listed in Table 3.1. The convergence
rates listed in the table are consistent with our theoretical result.
Example 2. In this example we consider a 10 cm 10 cm container  lled with microor-
ganisms in a suspension under zero external force, i.e., f  0. For computational simplicity,
68
we study the domain in two dimensional horizontal { vertical plan . The parameters of the
model, obtained from lab experiments (see [36]) are given in Table 3.2.
Table 3.2: Parameter values
 0 g   U
cm2=sec m=sec2 cm2=sec cm/sec
0.01 9.81 0.1 0.0025 0.01
De ne
 (c) =
8
>>>>
>>>
<
>>>
>>>
>:
 0; c< 0;
 0(1 + 2:5 c+ 5:3 c2); 0 <c< 10%;
 0 exp( 2:5 c1 1:4 c); 10% <c< 60%;
 0 exp(9:375); c> 60%:
(3.48)
Here  0 is the viscosity of the culture  uid. The viscosity (3.48) combines the work of Batche-
lor [43] for low concentration and Mooney [44] for high concentration. Note that exp( 2:5 c1 1:4 c)
is bounded below by   =  0 but tends to in nity when the maximum concentration ?m = 11:4
is reached since the suspension behaves as a solid. In this case no movement of neighboring
particles are allowed. Therefore we set the upper bound   =  0 exp(9:375) so that the
viscosity de ned in (3.48) satis es property (1.21).
We  rst chose  = 1%. The velocity and concentration are given in Figure 3.1. We can
see that a bioconvection pattern can not be formed and the concentration has a homogeneous
horizontal distribution. This is because the right hand side of the  rst equation in (1.16) is
almost equal to  g. As a result, u 0 while p is almost linear with rp  g and @c@x  0
because of the nearly zero velocity u. In this case, the micro-organisms do not move and
the concentration stays linear in the vertical direction with zero horizontal gradient. From
observed experiments, for a shallow container with low concentration of micro-organisms, the
micro-organisms will stay at the surface of the suspension due to the upswimming. Actually,
the e ect of gravity is due to high density of the micro-organisms but the micro-organisms
are almost isolated thus gravity can be neglected. In fact, bioconvection only occurs for
69
Figure 3.1: Concentration and velocity  eld for  = 1%
su ciently deep container. The higher the concentration is, the shallower the container can
be. We conclude that a %1 concentration is not enough to form a bioconvection pattern in
a 10 cm deep container.
Example 3. In this example, we assign the same parameter as in Example 2 except
that  = 20%. Figure 3.2 shows the concentration distribution and the velocity  led using
streamlines. Here the color denotes the magnitude of the velocity. The  gure shows that
a bioconvection pattern can be formed for su ciently large concentration. Our simulation
result is consistent with the results obtained in [36]. From the  gure we can see that two
convections, separated from the center,  ow steadily in opposite directions. The highest
velocity occurs in the center, where the concentration is low, due to the upswimming under
small e ect of gravity. Another high speed motion is observed on the left and right hand
sides of the container, which is caused mostly by gravity due to high concentration in the
upper left and right corners. Only a few micro-organisms remains at the bottom while most
of the micro-organisms stay close to the surface.
Example 4. In this example,  = 20% and constant viscosity  (c) 0:01. The result
is shown in Figure 3.3. From the graph, we can see that both example 3 and example 4
capture the motion of the bioconvection but the concentration distribution and velocity  eld
are slightly di erent. The velocity in the nonlinear case are slower and smoother due to a
relatively higher viscosity. The di erence is more notable in regions where the concentration
is high. Example 3 re ects higher concentrated micro-organisms at the top corners since more
70
Figure 3.2: Concentration and velocity  eld for  = 20%
Figure 3.3: Concentration and velocity  eld for  = 20% with constant viscosity
micro-organisms are washed up by the drag force and stay there due to the high viscosity.
One can see that the involvement of a nonhomogeneous viscosity improves the accuracy of
the simulation.
Example 5. In the last numerical experiment, all data are the same as in Example 3
except that  = 30%. The velocity  eld and concentration distribution are given in  gure 3.4.
As the concentration increases, the e ect of the gravity becomes more signi cant. However,
once the pattern is formed, the distribution of the concentration stays the same.
71
Figure 3.4: Concentration and velocity  eld for  = 30%
72
Chapter 4
Time dependent bioconvection
4.1 Existence and uniqueness of a solution
The weak formulation. In this chapter, we consider system (1.19) in a two dimen-
sional domain  . Similarly to the steady bioconvection case we chooseL20( ) as the functional
space for the pressure. Using the same auxiliary concentration c = c , force f = f and the
same bilinear and trilinear forms de ned in (3.1), we de ne the weak solution of system
(1.19) as follows.
De nition 4.1. Given f in L2(I; L2( )), u0 2 L2( ), c0 2 L2( ), a triple (u;p;c) 2
L2(I; H10( )) L2(I;L20( )) L2(I; ~H) is said to be a weak solution of system (1.19) if for
any t2(0;T]
8
>>>
>>>>
>>>
>>><
>>>
>>>>
>>>
>>>
:
hu0;vi+
 
 (c+ )ru;rv
 
+B2(u;u;v) +b(p;v)
= 
 
g(1 + c
 
i2;v) + (f;v) 8v2H10( );
b(q;u) = 0 8q2L20( );
hc0;ri+ a(c;r) +B3(u;c;r) U(c; @r@x
2
) = 0 8r2 ~H;
u(0) = u0; c(0) = c0:
(4.1)
Analogous to the steady state case, to solve system (4.1), it su ces to solve the following
associated problem.
73
Given f in L2(I; L2( )), u0 2 L2( ), and c0 2L2( ),  nd a pair (u;c) 2L2(I; V) 
L2(I; ~H) such that for any t2(0;T],
8
>>>
>><
>>>
>>:
hu0;vi+
 
 (c+ )ru;rv
 
+B2(u;u;v) = 
 
g(1 + c)i2;v
 
+ (f;v) 8v2V;
(c0;r) + a(c;r) +B3(u;c;r) U(c; @r@x
2
) = U ( @r@x
2
;1) 8r2 ~H;
u(0) = u0; c(0) = c0  :
(4.2)
Remark 4.2. Similarly to the steady case, if (u;c) is a solution of system (4.2), then (u;p;c),
the solution of (4.1), can be recovered because the bilinear form b( ; ) de ned above saties es
the inf-sup condition (3.2), i.e., for some  > 0
sup
v2H10( )
(q;r v)
kvk1   kqk; 8q2L
2
0( ):
For more details see [62], p. 59, Theorem I.4.1.
Existence. To establish the existence of a weak solution, we use the same modi ed
Rothe?s method developed in chapter 2 to construct a convergent sequence of approximate
solutions of (4.2) using the backward Euler approximation of the time derivative u0 and c0.
To this end we let k = T=n for some positive integer n, partition I uniformly with time step
k, and denote nodal points ti = tni = ik, for i = 1;2;:::;n. Let u0n = u0 and c0n = c0 and
de ne 8
>>>
>>>
<
>>>
>>>
:
fin := 1=k
Z ti
ti 1
f(t) dt;
 uin := (uin ui 1n )=k;
 cin := (cin ci 1n )=k:
For each integer n > 0, we apply the following scheme inductively to  nd uin 2 V and
cin 2 ~H, the approximations of the solution u and c at ti, i = 1;2;:::;n respectively, with
74
u0n = u0 and c0n = c0.
8
>>>
>><
>>>
>>:
( uin;v) +
 
 (cin + )ruin;rv
 
+B2(uin;uin;v)
= 
 
g(1 + cin)i2;v
 
+ (fin;v) 8v2V;
( cin;r) + a(cin;r) +B3(uin;cin;r) U(cin; @r@x
2
) = U ( @r@x
2
;1) 8r2 ~H:
(4.3)
Multiplying both sides of (4.3) by k, we obtain the following scheme.
Given ui 1n 2V, ci 1n 2 ~H, seek (uin;cin)2V ~H such that
8
>>>>
>>>
>>><
>>>
>>>>
>>>
:
(uin;v) +k
 
 (cin + )ruin;rv
 
+kB2(uin;uin;v)
= k
 
g(1 + cin)i2;v
 
+ (kfin + ui 1n ;v) 8v2V;
(cin;r) +k a(cin;r) +kB3(uin;cin;r) kU(cin; @r@x
2
)
= kU ( @r@x
2
;1) + (ci 1n ;r) 8r2 ~H:
(4.4)
Using (3.5), (3.6), (3.7), and (3.8), we establish the existence of a weak solution of (4.4)
as follows.
Lemma 4.3. Given ui 1n 2V, ci 1n 2 ~H, assume that
 
C2 >U: (4.5)
Then system (4.4) has a weak solution (uin;cin)2V ~H.
Proof. Analogous to the proof of Theorem 3.6, since V and ~H are both separable, there
existfvjg1j=1 andfrjg1j=1 orthonormal bases of V and ~H respectively. Let Vm, Hm be  nite
subspaces of V, ~H spanned by fv1;v2;:::;vmg and fr1;r2;:::;rmg. We consider the  nite
dimensional approximations of u and c in Vm and Hm, i.e., we seek um 2 Vm, cm 2Hm
75
such that 8
>>>
>>>>
>>>
<
>>>
>>>
>>>
>:
(um;v) +k
 
 (cm + )rum;rv
 
+kB2(um;um;v)
= k
 
g(1 + cm)i2;v
 
+ (kfin + ui 1n ;v) 8v2Vm;
(cm;r) +k a(cm;r) +kB3(um;cm;r) kU(cm; @r@x
2
)
= kU ( @r@x
2
;1) + (ci 1n ;r) 8r2Hm:
(4.6)
For any integar m> 0, the existence of a solution (um;cm) of (4.6) is guaranteed by Lemma
2.9. To show thatfumg1m=1 andfcmg1m=1 are bounded sequence in Vm and Hm respectively,
we take v = um and r = cm in (4.6), and use (1.21) and (3.8) to  nd
kumk2 + k  C2
 
kumk21  (um;um) +k
 
 (cm + )rum;rum
 
 
  
  k
 
g(1 + cm)i2;um
 
+ (kfin + ui 1n ;um)
  
 
 "kumk2 +C(")kcmk2 +C(")kkfin + ui 1n  kgi2k2
(4.7)
and
kcmk2 + k C2
 
kcmk21  (cm;cm) +k a(cm;cm)
=
  
 kU(cm;@c
m
@x2 ) +kU (
@cm
@x2 ;1) + (c
i 1
n ;c
m)
  
 
 kUkcmk21 +"kcmk21 +C(")(1 +kci 1n k2):
By assumption (4.5), we may choose a su ciently small "> 0 such that "< k c2
 
 kU. Then
the above inequality gives
kcmk1  C: (4.8)
Furthermore substituting (4.8) into (4.7) and choosing 0 <"< 1 we have that
kumk1  C:
76
This means that the sequencesfumg1m=1 andfcmg1m=1 are bounded in V and ~H respectively.
Therefore there exist uin2V and cin2 ~H such that
um * uin in V and cm *cin in ~H as m!1: (4.9)
Due to Sobolev embedding theorem, this yields
um!uin in L2( ) and cm!cin in L2( ) as m!1: (4.10)
Next we show that the weak limit (u;c) is a solution of (4.4). Choose test functions
v2V\(C10 ( ))3; r2C1( )\ ~H: (4.11)
Similarly to the strategy used in the proof of Theorem 3.6, we let m tend zero and use the
weak and strong convergence of fumg1n=1 and fcmg1n=1 to conclude that
8
>>>
>>>>
>>>
>>>
<
>>>
>>>
>>>>
>>>
:
 
 (cm + )rum;rv
 
!
 
 (c+ )ruin;rv
 
as m!1;
B2(um;um;v)!B2(uin;uin;v); B3(um;cm;r)!B3(uin;cin;r) as m!1;
(um;v)!(uin;v); (cm;r)!(cin;r) as m!1;
(g(1 + cm)i2;v)!(g(1 + cin)i2;v);  a(cm;r)! a(cin;r) as m!1;
U(cm; @r@x
2
)!U(cin; @r@x
2
) as m!1:
Since v and r de ned in (4.11) are dense in V and ~H, the limits uin and cin satisfy (4.6), i.e.,
8
>>>
>>>>
>>>
<
>>>
>>>
>>>
>:
(uin;v) +k
 
 (cin + )ruin;rv
 
+kB2(uin;uin;v)
= k
 
g(1 + cin)i2;v
 
+ (kfin + ui 1n ;v) 8v2Vm;
(cin;r) +k a(cin;r) +kB3(uin;cin;r) kU(cin; @r@x
2
)
= kU ( @r@x
2
;1) + (ci 1n ;r) 8r2Hm:
(4.12)
77
Because such test functions v and r are dense in V and ~H, the above system holds for all
v2V and r2 ~H. This means that (uin;cin), is a solution of (4.4).
By the modi ed Rothe?s method, we now construct a sequence to approximate the
solution of (4.2). To this end, for each positive integer n, we de ne two piecewise constant
functions
Un(t) = uin; Cn(t) = cin; t2(ti 1;ti); i = 1;2;:::;n; (4.13)
with Un(0) = u0 and Cn(0) = c0. Notice that piecewise constant function does not have
derivative in L2(I;L2( )). Therefore to approximate the time derivative in (3.3), we de ne
the following two piecewise linear functions
~Un(t) = uin + (t ti 1)(uin ui 1n )=k;
~Cn(t) = cin + (t ti 1)(cin ci 1n )=k;
t2(ti 1;ti]; i = 1;2;:::;n;
(4.14)
with ~Un(0) = u0 and ~Cn(0) = c0. The following estimate is true for any integer n> 0.
Lemma 4.4. Assume (4.5) is satis ed, then for any positive integer n0  n, we have the
estimate
kun0n k+kcn0n k+k
n0X
i=1
(kuink21 +kcink21 +k uink2V0 +kcink2~H0) C; (4.15)
where C is a constant independent of n and n0.
Proof. Taking v = uin and r = cin in (4.4), we deduce from (3.8) that
8
>><
>>:
( uin;uin) +
 
 (cin + )ruin;ruin
 
= 
 
g 1 + cin)i2;uin
 
+ (fin;uin);
( cin;cin) + a(cin;cin) U(cin;@c
i
n
@x2 ) =
U 
j j(
@cin
@x2;1):
(4.16)
78
Notice the identities
8
><
>:
(uin ui 1n ;uin) = 1=2(kuin ui 1n k2 +kuink2 kui 1n k2);
(cin ci 1n ;cin) = 1=2(kcin ci 1n k2 +kcink2 kci 1n k2):
(4.17)
Similarly to the way of obtaining (2.20), we multiply (4.16) and sum from i = 1 to n0.
Applying (1.21) and Young?s inequality we  nd that
1
2
n0X
i=1
kuin ui 1n k2 + 12kun0n k2 + k  C2
 
n0X
i=1
kuink21
 k
n0X
i=1
  
  
 
g(1 + cin)i2;uin
 
+ (fin;uin)
  
 + 12ku0nk2
 k
n0X
i=1
 
"kuink2 +C(") kfink2 +kcink2 + 1 
 
+ 12ku0nk2
(4.18)
and
1
2
n0X
i=1
kcin ci 1n k2 + 12kcn0n k2 + k C2
 
n0X
i=1
kcink21
 k
n0X
i=1
jU(cin;@c
i
n
@x2 ) +U (
@cin
@x2;1)j+
1
2kc
0
nk
2
 k
n0X
i=1
 
Ukcink21 +"kcink2 +C(")
 
+ 12kc0nk2:
(4.19)
By (4.5), we may choose "> 0 to be su ciently small such that "< k C2
 
 kU. Then (4.19)
gives
kcn0n k2 +k
n0X
i=1
kcink21  C: (4.20)
In addition we have
kfink2  1=k2
 Z
ti
ti 1
f(t) dt
!
 1=k
Z ti
ti 1
kf(t)k2 dt:
79
Hence
k
n0X
i=1
kfink2  
n0X
i=1
Z ti
ti 1
kf(t)k2 dt kfkL2(I;L2( )): (4.21)
Substituting (4.20) and (4.21) into (4.18) we conclude that
kun0n k2 +k
n0X
i=1
kuink21  Ck
n0X
i=1
 
kcink21 +kfink2 + 1
 
+ 12ku0nk2  C: (4.22)
To obtain an estimate for  uin and  cin, we  rst recall the following estimate ([52], p. 292,
Lemma 3.3 and Lemma 3.4)
krkL4( )  214krk12krrk12 ; 8r2H1( ): (4.23)
Applying Holder?s inequality to (3.7) we  nd
jB3(u;c;r)j=jB3(u;r;c)j kukL4( )krrkkckL4( )
 
 2X
i=1
kuik2L4( )
!1
2
kckL4( )krk1
 C
 2X
i=1
kuikkruik
!1
2
kck12krck12krk1
 C(kukkuk1kckkck1)12krk1 8u2V; 8c;r2H1( ):
(4.24)
Similarly
jB2(u;v;w)j C(kukkuk1kvkkvk1)12kwk1 8u2V; 8v;w2H10( ): (4.25)
Recall, from (4.22), thatkuinkandkcinkare uniformly bounded with respect to i. Combining
(4.20), (4.22), (1.21), and (4.25), we apply Cauchy-Schwarz inequality to deduce from (4.3)
80
that
( uin;v) = 
 
 (cin + )ruin;rv
 
 B2(uin;uin;v) 
 
g(1 + cin)i2;v
 
+ (fin;v)
   kuink1kvk1 +Ckuinkkuink1kvk1 +Ckcink1kvk1 +kfin gi2kkvk1;
 C
 
kuink1 +kcink1 +kfin gi2k
 
kvk1 8v2H10( )
(4.26)
and
( cin;r) =  a(cin;r) B(uin;cin;r) +U(cin; @r@x
2
) + U j j( @r@x
2
;1)
 Ckcink1krk1 +C(kuinkkuink1kcinkkcink1)12krk1 +Ckcink1krk1 +Ckrk1
 C
 
kcink1 +kuink1 + 1
 
krk1 8r2 ~H:
(4.27)
Consequently 8
><
>:
k uinkV0 C
 
kcink1 +kuink1 +kfink
 
;
k cink~H0 C
 
kcink1 +kuink1 + 1
 
:
Using (4.20), (4.22), and (4.21)
k
n0X
i=1
 
k uink2V0 +k cink2~H0
 
 C:
In view of Lemma 4.4 there exists a constant C such that
kUnk2L2(I;V) +kCnk2L2(I; ~H) = k
nX
i=1
(kuink21 +kcink21) C: (4.28)
Thus there exist u2V, c2 ~H such that
Un * u in L2(I; V) and Cn *c in L2(I; ~H): (4.29)
81
It is straight forward to verify that the piecewise constant function de ned in (4.13) share
the same weak and strong limits with the piecewise linear function in (4.14). As a result
~Un * u in L2(I; H10( )) and ~Cn *c in L2(I; ~H):
Moreover Lemma 4.4 implies that
k~U0nk2L2(I;V0) +k~C0nk2L2(I; ~H0 = k
nX
i=1
(k uink2V +k cink2~H0) C:
Hence there exist  u2L2(I; V0) and  c2L2(I; ~H0) such that
~U0n *  u in L2(I; V0) and ~C0n *c in L2(I; ~H0): (4.30)
It is easy to verify that ([58], p. 356)
 u = u0 and  c = c0: (4.31)
Then in view of Lemma 2.24, we set X0 = V, X = L2( )\V, and X1 = V0 to use the
compact embedding
fu2L2(I; V); u02Lq(I; V0)g,!L2(I; L2( )\V);
and we set X0 = ~H, X = L20( ), and X1 = ~H0 to obtain that
fc2L2(I; ~H); u02Lq(I; ~H0)g,!L2(I;L20( )):
Consequently
~Un!u in L2(I; L2( )) and ~Cn!c in L2(I;L2( )):
82
Thus
Un!u in L2(I; L2( )) and Cn!c in L2(I;L2( )): (4.32)
Armed with the above Lemmas and Lemma 2.12, we are now ready to show that (u;c)
is a solution of (4.2) with weak derivative u0, c0 as de ned in (4.30) and (4.31) .
Theorem 4.5 (Existence). Suppose (4.5) is satis ed. Given f 2L2(I; L2( )), u0 2L2( ),
and c0 2L2( ), system (4.2) has a weak solution (u;c) 2L2(I; V) L2(I; ~H) with u0 2
L2(I; V0) and c02L2(I; ~H0).
Proof. De ne
fn(t) = fin; t2(ti 1;ti]; i = 1;2;:::;n:
Proceeding as in Lemma 2.12, we have
kfn fkL2(I;L2( )) !0 as n!1: (4.33)
Construct the test functions of the form ~v = v? and ~r = r? with
v2(C10 ( ))3\V; r2C10 ( )\ ~H; and ?(t)2C10 (I): (4.34)
We have the following identities
nX
i=1
(uin ui 1n ;v)?(ti)
= (unn;v)?(tn) (u0n;v)?(t1) k
n 1X
i=1
(uin;v)
 
?(ti+1) ?(ti)
 
=k
83
and
nX
i=1
(cin ci 1n ;r)?(ti)
= (cnn;r)?(tn) (c0n;r)?(t1) k
n 1X
i=1
(cin;r)
 
?(ti+1) ?(ti)
 
=k:
Multiplying (4.3) with k?(ti) and summing from i = 1 to n, we obtain
(unn;v)?(T) (u0n;v)?(t1) k
n 1X
i=1
(uin;v)
 
?(ti+1) ?(ti)
 
=k
+k
nX
i=1
 
 (cin + )ruin);rv
 
?(ti) +B2(uin;uin;v)?(ti)
=k
nX
i=1
 
  g(1 +cin)i2;v + (fin;v)
 
?(ti)
and
(cnn;r)?(T) (c0n;r)?(t1) k
n 1X
i=1
(cin;r)
 
?(ti+1) ?(ti)
 
=k
+k
nX
i=1
 
 a(cin;r) +kB3(uin;cin;r) U(cin; @r@x
2
)
 
?(ti)
=k
nX
i=1
U 
j j(
@r
@x2;1)?(ti):
84
From the de nition of Un, Cn, ~Un, and ~Cn the above equations yield
8
>>>
>>>>
>>>
>>>
>>>>
>>>
><
>>>>
>>>
>>>
>>>>
>>>
>>>>
:
 (u0;v)?(k) 
Z T
0
(Un;v) ~?n dt
+
Z T
0
 
 (Cn + )rUn;rv
 
?n dt+
Z T
0
B2(Un;Un;v)?n dt
= 
Z T
0
(g(1 +Cn)i2;v)?n dt+
Z T
0
(fn;v)?n dt;
 (c0;r)?(k) +
Z T
0
(Cn;r) ~?n dt
+
Z T
0
 a(Cn;r)?n dt+
Z T
0
B3(Un;Cn;r)?n dt
 
Z T
0
U(Cn; @r@x
2
)?n dt =
Z T
0
U 
j j(
@r
@x2;1)?n dt:
(4.35)
The continuity of ? implies (notice that ?(0) = 0) that
(u0;v)?(k)!0; and (c0;v)?(k)!0 as n!1: (4.36)
From Lemma 2.12, (4.28), and Holder?s inequality it follows that
  
 
Z T
0
(Un;v)( ~?n ?0) dt
  
  kvkkUnkL2(I;L2( ))k~?n ?0kL2(I) !0 as n!1
and
  
 
Z T
0
(Cn;r)( ~?n ?0) dt
  
  krkkCnkL2(I;L2( ))k~?n ?0kL2(I) !0 as n!1:
85
Together with (4.29) we have
8>
>>>
>>>>
>>>
><
>>>
>>>>
>>>
>>:
Z T
0
(Un;v) ~?n dt =
Z T
0
(Un;v)( ~?n ?0) dt+
Z T
0
(Un;v)?0 dt
!
Z T
0
(u;v)?0 dt;
Z T
0
(Cn;r) ~?n dt =
Z T
0
(Cn;r)( ~?n ?0) dt+
Z T
0
(Cn;r)?0 dt
!
Z T
0
(c;r)?0 dt
(4.37)
as n !1. Next the Nemytskii property (1.21) leads to
 (Cn + )! (c+ ) in L2(I;L2( )) as n!1:
Thus (1.21), (4.29), (4.28), Lemma 2.12, and Holder?s inequality yield
  
 
Z T
0
 
 (Cn + )rUn;rv
 
?n dt 
Z T
0
 
 (C + )ru;rv
 
?n dt
  
 
=
  
 
Z T
0
 
 (Cn + )rUn;rv
 
(?n ?) dt
+
Z T
0
 
 (C + )(rUn ru);rv)
 
? dt
+
Z T
0
  
 (Cn + )  (C + ) rUn;rv
 
? dt
  
 
   kvk1kUnkL2(I;V)k?n ?kL2(I)
+  
  
 
Z T
0
(rUn ru;rv)? dt
  
 
+kvk1k (Cn + )  (C + )kL2(I;L2( ))kUnkL2(I;V)
!0; as n!1:
(4.38)
86
Combing (4.28), Lemma 2.12, and Holder?s inequality we let n tend 1 to  nd
8
>>>>
>>>
>>>>
><
>>>
>>>
>>>>
>>:
  
 
Z T
0
 
g(1 +Cn)i2;v
 
(?n ?) dt
  
  
 
kvk1kCnkL2(I;L2( )) +C
 
k?n ?kL2(I) !0;
  
 
Z T
0
 a(Cn;r)(?n ?) dt
  
   krk1kCnkL2(I;V)k?n ?kL2(I) !0;
  
 
Z T
0
U(Cn; @r@x
2
)(?n ?) dt
  
  Ukrk1kCnkL2(I;L2( ))k?n ?kL2(I) !0;
  
 
Z T
0
U 
j j(
@r
@x2;1)(?n ?) dt
  
  Ckrk1k?n ?kL2(I) !0:
Hence by (4.29) we conclude that as n!1
8
>>>
>>>
>>>>
>>>
>>>>
>>>
>>>>
>>>
>>>
<
>>>
>>>
>>>>
>>>
>>>>
>>>
>>>>
>>>
>>>
:
Z T
0
 
g(1 +Cn)i2;v
 
?n dt =
Z T
0
 
g(1 +Cn)i2;v
 
(?n ?) dt
+
Z T
0
 
g(1 +Cn)i2;v
 
? dt
!
Z T
0
 
g(1 +C)i2;v
 
? dt;
Z T
0
 a(Cn;r)?n dt =
Z T
0
 a(Cn;r)(?n ?) dt+
Z T
0
 a(Cn;r)? dt
!
Z T
0
 a(C;r)? dt;
Z T
0
U(Cn; @r@x
2
)?n dt =
Z T
0
U(Cn; @r@x
2
)(?n ?) dt+
Z T
0
U(Cn; @r@x
2
)? dt
!
Z T
0
U(C; @r@x
2
)? dt;
Z T
0
U 
j j(
@r
@x2;1)?n dt!
Z T
0
U 
j j(
@r
@x2;1)? dt:
(4.39)
87
Note that
Z T
0
(fn;v)?n dt =
nX
i=1
Z ti
ti 1
 
1=k
Z ti
ti 1
f( )d ;v
 
?(ti) dt
=
nX
i=1
 Z ti
ti 1
f( )d ;v
 
?(ti)
=
nX
i=1
 Z ti
ti 1
f( );v
 
?n( ) d 
=
Z T
0
 
f(t);v
 
?n(t) dt:
Together with (4.33) and Lemma 2.12
  
 
Z T
0
(fn;v)?n dt 
Z T
0
(f;v)? dt
  
 
 
  
 
Z T
0
(f;v)(?n ?) dt
  
 
 kvk1kfkL2(I;L2( ))k?n ?kL2(I) !0 as n!1:
(4.40)
Next we show that
Z T
0
(B2(Un;Un;v)?n dt!
Z T
0
B2(u;u;v)? dt as n!1: (4.41)
Write
 Z T
0
(B2(Un;Un;v)?n dt 
Z T
0
B2(u;u;v)? dt
 
=
Z T
0
B2(Un;Un;v)(?n ?) dt+
Z T
0
 
B2(Un;Un;v) B2(u;u;v)
 
? dt:
88
By (4.25), Lemma 2.12, Lemma 4.4, and the uniform boundedness of kuink we have that
j
Z T
0
(B2(Un;Un;v)(?n ?) dtj
 Cj
Z T
0
kUnkkUnk1kvk1(?n ?) dtj
 Ckvk1kUnkL2(I;V)k?n ?kL2(I) !0 as n!1:
By (3.7),(4.32), and (4.34), we use integration by parts to  nd that as n!1
Z T
0
B2(Un;Un;v)? dt = 
Z T
0
B2(Un;v;Un)? dt
=
2X
i;j=1
Z T
0
Z
 
(Un)i(@vj@x
i
)(Un)j dx ? dt
!
2X
i;j=1
Z T
0
Z
 
ui(@vj@x
i
)uj dx ? dt
= 
Z T
0
B2(u;v;u)? dt
=
Z T
0
B2(u;u;v)? dt;
from which (4.41) follows. Analogously, (4.24), Lemma (4.4), and Holder?s inequality yield
j
Z T
0
B3(Un;Cn;r)(?n ?) dtj
 C
Z T
0
kUnkkUnk1kCnkkCnk1)12 )kvk1j?n ?jdt
 C
Z T
0
(kUnk1 +kCnk1)kvk1j?n ?jdt
 Ckvk1
 
kUnkL2(I;V) +kCnkL2(I; ~H)
 
k?n ?kL2(I) !0 as n!1
and Z
T
0
B3(Un;Un;v)? dt!
Z T
0
B3(u;u;v)? dt:
89
Together the above estimates give
Z T
0
B(Un;Cn;r)?n dt!
Z T
0
B(u;c;r)? dt; as n!1: (4.42)
Combing (4.36), (4.37), (4.38), (4.39), (4.40), (4.41), and (4.42), we deduce from system
(4.35) that
8
>>>>
>>>
>>>>
><
>>>
>>>
>>>>
>>:
 
Z T
0
(u;~v0)dt+
Z T
0
 
 (c+ )ru;r~v
 
dt+
Z T
0
B2(u;u;~v) dt
= 
Z T
0
 
g(1 + c)i2;~v) dt+
Z T
0
(f;~v) dt;
 
Z T
0
(c;~r0) dt+
Z T
0
 a(c;~r) dt+
Z T
0
B3(u;c;~r) dt 
Z T
0
U(c; @~r@x
2
) dt
=
Z T
0
U 
j j(
@~r
@x2;1) dt;
(4.43)
for any ~v and ~r of the form in (4.11). Because such ~v and ~r are dense in L2(I; V) and
L2(I; ~H) respectively, system (4.43) holds for all ~v 2 L2(I; V) and ~r 2 L2(I; ~H). Hence
(u;c)2V ~H satis es (4.2) in sense of distribution. Since the weak derivative u02L2(I; V0)
and c02L2(I; ~H0) both exist, the pointwise version (4.2) is true for a.e.t2I, i.e., (u;c) is
the solution of (4.2).
Uniqueness. Analogous to the steady case, the bilinear form b( ; ) satis es the inf-
sup condition (3.2). Thus for each solution (u;c) of system (4.2), there exists a unique
p 2 L2(I;L20( )) satisfying system (4.1) (see [62], p. 59, Theorem I.4.1). As a result, to
prove the uniqueness of solution to system (4.1), it su ces to prove that system (4.2) has a
unique solution (u;c).
Theorem 4.6. Suppose
(H8) The hypotheses of Theorem 4.5 hold;
90
(H9) The viscosity  ( ) is Lipschitz continuous, i.e., there exists  L > 0 such that
j (x1)  (x2)j  Ljx1 x2j 8x1;x2 2R;
(H10) There exists a constant C0 such that
kru)k1 C0 8t2I:
Then the solution (u;c) of system (4.2) is unique.
Proof. Let (u1;c2) and (u2;c2) be two di erent solutions of (4.2). Substituting (u;c) with
(u1;c1) and (u2;c2) in (4.2) we have
8
>>>
>>>>
>><
>>>>
>>>
>>:
hu01;vi+
 
 (c1 + )ru1;rv
 
+B2(u1;u1;v)
= 
 
g(1 + c1)i2;v
 
+ (f;v) 8v2V;
(c01;r) + a(c1;r) +B3(u1;c1;r) U(c1; @r@x
2
) = U ( @r@x
2
;1) 8r2 ~H;
u1(0) = u0; c1(0) = c0  
(4.44)
and
8
>>>
>>>>
>><
>>>>
>>>
>>:
hu02;vi+
 
 (c2 + )ru2;rv
 
+B2(u2;u2;v)
= 
 
g(1 + c2)i2;v
 
+ (f;v) 8v2V;
(c02;r) + a(c2;r) +B3(u2;c2;r) U(c2; @r@x
2
) = U ( @r@x
2
;1) 8r2 ~H;
u2(0) = u0; c2(0) = c0  :
(4.45)
91
Subtracting (4.45) from (4.44) with v = u1 u2 and r = c1 c2, we  nd that
h(u1 u2)0;u1 u2i+
 
 (c1 + )ru1  (c2 + )ru2;r(u1 u2)
 
+B2(u1;u1;u1 u2) B2(u2;u2;u1 u2)
= g 
 
(c1 c2)i2;u1 u2
 
(4.46)
and
h(c1 c2)0;c1 c2i+ a(c1 c2;c1 c2)
+B3(u1;c1;c1 c2) B3(u2;c2;c1 c2)
 U
 
c1 c2;@(c1 c2)@x
2
 
= 0:
(4.47)
From (3.8)
8
><
>:
B2(u1;u1;u1 u2) B2(u2;u2;u1 u2) = B2(u1 u2;u2;u1 u2);
B3(u1;c1;c1 c2) B3(u2;c2;c1 c2) = B3(u1 u2;c2;c1 c2):
Using (3.4) and(4.24), we apply Young?s inequality to deduce from (4.47) that
1
2
d
dtkc1 c2k
2 +  
C2 kc1 c2k
2
1
 h(c1 c2)0;c1 c2i+ a(c1 c2;c1 c2)
 jB3(u1 u2;c2;c1 c2)j+Ukc1 c2kkc1 c2k1
 C
 
ku1 u2kku1 u2k1kc2kkc2k1
 1
2kc
1 c2k1 +Ukc1 c2k21
 ("+U)kc1 c2k21 +C(")
 
ku1 u2kku1 u2k1kc2kkc2k1
 
:
(4.48)
92
We rewrite (4.46) in the following form
h(u1 u2)0;u1 u2i+
  
 (c1 + )  (c2 + ) ru1;r(u1 u2)
 
+
 
 (c2 + )r(u1 u2);r(u1 u2)
 
+B2(u1 u2;u2;u1 u2)
= g 
 
(c1 c2)i3;u1 u2
 
:
Then using (1.21), (3.4), (4.25), (H9), and Young?s inequality yield that
1
2
d
dtku1 u2k
2 +   
C2 ku1 u2k
2
1
 h(u1 u2)0;u1 u2i+
 
 (c2 + )r(u1 u2);r(u1 u2)
 
 
  
 
  
 (c1 + )  (c2 + ) ru1;r(u1 u2)
 
+B2(u1 u2;u2;u1 u2) +g 
 
(c1 c2)i3;u1 u2
   
 
 C0 Lkc1 c2kku1 u2k1 +ku1 u2kkc1 c2k
+C
 
ku1 u2kku1 u2k1ku2kku2k1
 1
2ku
1 u2k1
 "ku1 u2k21 +C(")
 
ku1 u2kku1 u2k1ku2kku2k1 +kc1 c2k2
 
:
(4.49)
Hypothesis (4.5) implies that we may choose a su ciently small "> 0 such that
"< minf  C2
 
;  C2
 
 Ug:
93
Taking the sum of (4.48) and (4.49) and applying Young?s inequality we  nd that
d
dtku1 u2k
2 +ku1 u2k2
1 +
d
dtkc1 c2k
2 +kc1 c2k2
1
 C
 
ku1 u2kku1 u2k1kc2kkc2k1
+ku1 u2kku1 u2k1ku2kku2k1
+kc1 c2k2
 
 "ku1 u2k21
 
kc2k2 +ku2k2
 
+C(")ku1 u2k2
 
kc2k21 +ku2k21
 
+Ckc1 c2k2:
Proceeding as in Lemma 4.4, we can show that ku2k and kc2k are uniformly bounded with
respect to t. Hence choosing "> 0 such that
"< 1maxfkc
2k2;ku2k2g
;
we obtain that
d
dtku1 u2k
2 + d
dtkc1 c2k
2
 Cku1 u2k2
 
kc2k21 +ku2k21
 
+Ckc1 c2k2:
Using Gronwall?s inequality we conclude that for any t2I (recall that bothku2k21 andkc2k21
are integrable on I)
ku1 u2k2(t)+kc1 c2k2(t)
 ku1 u2k2(0) exp(
Z t
0
 
ku2k21 +kc2k21
 
ds
+Ckc1 c2k2(0):
94
Thus
ku1 u2k=kc1 c2k= 0; a:e:t2I: (4.50)
Remark 4.7. Comparing Theorem 4.6 with Theorem 3.7, we can see condition (H7) is not
needed to prove the uniqueness of the solution, i.e., the uniqueness of solution of time
dependent bioconvection can be obtained under weaker conditions than in the steady case.
This is general for partial di erential equations with solution dependent coe cient. We also
note that we restrict the equation to a two dimensional model.
4.2 Numerical approximation
In this section, we consider the semi-discrete  nite element approximation to solutions
of (4.1). Throughout this section, we assume that the weak solution (u;p;c) of (3.2) exists
and is unique.
We use the same  nite element spaces as for the steady bioconvection. Let  h be a
family of quasi-uniform triangulations of the convex polygonal domain  satisfying
max 2 
h
diam   h;
where h refers to the mesh size. Let Xh, Mh, and Sh be the corresponding  nite dimensional
subspaces of H10( ), L20( ), and ~H, respectively, that satisfy the following approximation
properties
inf
vh2Xh
kv vhk1  Chskvks+1 8v2Hs+1( ); 0 <s k; (4.51)
inf
qh2Mh
kq qhk Chskqks 8q2Hs( ); 0 <s k; (4.52)
inf
th2Sh
kt thk1  Chsktks+1 8t2Hs+1( ); 0 <s k; (4.53)
95
for k = 2, 3. For the construction of these spaces, see [63, 62, 52]. Let Uh and Ch be the
approximation of the initial conditions u0 and c0 in Xh and Sh respectively. The semi-discrete
 nite dimensional approximation of (4.1) is de ned as follows.
Given Uh2Xh, Ch2Sh,  nd (uh;ph;ch)2Xh Mh Sh such that for any t2I
8
>>>
>>>>
>>>
>>><
>>>
>>>>
>>>
>>>:
hu0h;vi+
 
 (ch + )ruh;rv
 
+B2(uh;uh;v) +b(ph;v)
= 
 
g(1 + ch
 
i2;v) + (f;v) 8v2Xh;
b(q;uh) = 0 8q2Mh;
hc0h;ri+ a(ch;r) +B3(uh;ch;r) U(ch; @r@x
2
) = 0 8r2Sh;
uh(0) = Uh; ch(0) = Ch:
(4.54)
As for the steady biocovnection, we assume that these  nite spaces have explicit bases and
that a discrete version of the inf-sup condition (3.26) is satis ed, i.e., for some  > 0,
sup
v2Xh
b(v;q)
kvkXh   kqkMh 8q2Mh: (4.55)
For the construction of these spaces, see [63, 62, 52]. De ne the discrete divergence free
space
Vh =fv2Xh : (r v;qh) = 0; 8qh2Mhg
and the auxiliary forms ^B2 and ^B3
8
><
>:
^B2(u;v;w) = 1
2B2(u;v;w) 
1
2B2(u;w;v);
^B3(u;c;r) = 1
2B3(u;c;r) 
1
2B3(u;r;c):
We recall the following properties
^B2(u;v;w) = B2(u;w;v) 8v;w2H10( );
^B3(u;c;r) = B3(u;c;r) 8u2V; 8c;r2 ~H;
(4.56)
96
and the identities
^B2(u;v;v) = 0; ^B3(u;c;c) = 0 8u;v2H10( ); 8c2 ~H; (4.57)
and the tricontinuous properties
8
><
>:
^B2(u;v;w) CB
2kuk1kvk1kwk1 8u;v;w2H
1
0( );
^B3(u;c;r) CB
3kuk1kck1krk1 8u2H
1
0( ); 8c;r2H
1( ):
(4.58)
It is straight forward to verify that the discrete versions of estimates (4.25) and (4.24) hold,
that is,
jB3(u;c;r)j C(kukkuk1kckkck1)12krk1 8u2V; 8c;r2H1( ); (4.59)
and
jB2(u;v;w)j C(kukkuk1kvkkvk1)12kwk1 8u2V; v;w2H10( ): (4.60)
We  rst consider the discrete version of (4.2).
Find a pair (uh;ch) such that for 8t2I,
8>
>>>
><
>>>>
>:
(u0h;v) +
 
 (ch + )ruh;rv
 
+ ^B2(uh;uh;v)
= (g(1 + ch)i2;v) + (f;v) 8v2Vh;
(c0h;r) + a(ch;r) + ^B3(uh;ch;r) U(ch; @r@x
2
) = U ( @r@x
2
;1) 8r2Sh:
(4.61)
The existence of the above scheme is guaranteed by general ordinary di erential equation
theory. Taking v = uh and r = ch in (4.61), using (1.21), (1.22), and (3.31), and applying
97
Young?s inequality we obtain
1
2
d
dtkuhk
2 +   
C2 kuhk
2
1
 (u0h;uh) +
 
 (ch + )ruh;ruh
 
 j 
 
g(1 + ch)i2;uh
 
+ (f;uh)j "kuhk21 +C
 
kf gi2k2 +kchk2
 
(4.62)
and
1
2
d
dtkchk
2 +  
C2 kchk
2
1  (c
0
h;ch) + a(ch;ch)
 jU(ch;@ch@x
2
) +U (@ch@x
2
;1)j
 ("+U)kchk21 +C:
(4.63)
By (4.5) we can choose small "> 0 such that "<  C2
 
 U. Then (4.63) gives
d
dtkchk
2 +kchk2
1  C:
Thus (4.62) yields
d
dtkuhk
2 +kuhk2
1  kf gi2k
2 +kchk2  C:
Integrating (4.62) and (4.63) with respect to t and taking the sum, we obtain the estimate
kuhk2 +kchk2 +
Z t
0
 
kuhk21 +kchk21
 
d  C 8t2I: (4.64)
Following the same argument, we establish a similar estimate for the exact solution (u;c),
that is
kuk2 +kck2 +
Z t
0
 
kuk21 +kck21
 
ds C 8t2I: (4.65)
98
After (uh;ch) are computed, we compute ph2Mh by solving
(ph;r v) =(u0h;v) +
 
 (ch + )ruh;rv
 
+B2(uh;uh;v)
+
 
g(1 + ch)i2;v
 
 (f;v) 8v2Xh:
(4.66)
As in the steady case, by the property of Lagrange multiplier and (4.55), the above equation
is always solvable and the solution ph2Mh is unique in the quotient space Mh=Nh where
Nh =fqh2Mh : (qh;r v) = 0 8v2Xhg:
In this way, the pressure ph depends continuously on the discrete solution uh.
To obtain the error estimates, we use the Ritz Galerkin projections ([62], p. 132-139)
rh : H10( )!Vh, sh : ~H!Sh , and the L2 projection  h : L20( )!Mh to split the errors
into two parts: 8
>>>
>><
>>>>
>:
u uh = u rhu +rhu uh :=  hu + hu;
p ph = p  hp+ hp ph :=  hp + hp ;
c ch = c shc+shc ch :=  hc + hc :
(4.67)
The convergence of the projection error is guaranteed by (4.51), that is
k huk1 !0; k hpk1 !0; k hck1 !0 as h!0: (4.68)
We also have the following estimate for the projection
krhuk1  C(kuk1); kshck1  C(kck1); k hpk C(kpk1): (4.69)
Then the main convergence theorem of the numerical approximation is stated as follows.
Theorem 4.8. Assume that
(H11) The assumptions of Theorem 4.6 hold;
99
(H12) u2C1(I; H10( ))\V and c2C1(I;H10 ( ))\ ~H (see [52], p. 211).
Then the solution (uh;ch) converges to the exact solution (u;c), i.e.,
ku uhk+kc chk!0 as h!0:
Proof. In light of (4.68), it su ces to estimate  hu,  hp and  hc. Subtracting (4.1) from (4.54)
with v =  hu, r =  c we have that
(u0h u0; hu)+
 
 (ch + )ruh;r hu
 
 
 
 (c+ )ru;r hu
 
+ ^B2(uh;uh; hu) ^B2(u;u; hu)
+b(p ph; hu)
= g 
 
(ch c))i2; hu
 
(4.70)
and
(c0h c0; hc) + a(ch c; hc) + ^B3(uh;ch; hc) ^B3(u;c; hc) U(ch c;@ 
h
c
@x2 ) = 0: (4.71)
Notice that  hu2Vh and  hp 2Mh. Therefore the de nition of Vh implies that b( hp; hu) = 0,
i.e.,
b(p ph; hu) = b( hp; hu):
The following identities are guaranteed by (4.57) .
^B2(uh;uh; hu) ^B2(u;u; hu) = ^B2( hu;uh; hu) + ^B2(rhu; hu; hu) + ^B2( hu;u; hu)
and
^B3(uh;ch; hc) ^B3(u;c; hc) = ^B3( hu;ch; hc) + ^B3(rhu; hc; hc) + ^B3( hu;c; hc):
100
Then combing (1.21), (3.4), (4.60), (4.59), (4.64), (4.65), and (H12), we apply Young?s
inequality to deduce from (4.70) and (4.71) that
1
2
d
dtk 
h
uk
2 +   
C2 k 
h
uk
2
1
 (( hu)0; hu) +
 
 (ch + )r hu;r hu
 
 j(( hu)0; hu) +
 
 (c+ )r hu;r hu
 
+
  
 (ch + )  (c+ ) ru;r hu
 
+ ^B2( hu;uh; hu) + ^B2(rhu; hu; hu) + ^B2( hu;u; hu)
+b( hp; hu) +g 
 
(ch c)i2; hu
 
j
 k( hu)0kk huk+  k huk1k huk1 +C0 Lkc chkk huk1
+k hpkk huk1 +g kch ckk huk
+Ck uk1
 
(k hukk huk1kuhkkuhk1)12 +krhuk1k huk1 +k huk1kuk1
 
 "k uk21 +C(kuk1;")
 
k( hu)0k2 +k huk21 +k ukk uk1kuhkkuhk1
+k hpk2 +k hck2 +k hck2
 
(4.72)
and
1
2
d
dtk 
h
ck
2 +  
C2 k 
h
ck
2
1
 (( hc)0; hc) + a( hc; hc)
 
  
 (( hc)0; hc) + a( hc; hc) +U(ch c;@ 
h
c
@x2 )
+ ^B3( hu;ch; hc) + ^B3(rhu; hc; hc) + ^B3( hu;c; hc)
  
 
 k( hc)0kk hck+ k hck1k hck1 +Ukch ckk hck1
+Ck hck1
 
(k hukk huk1kchkkchk1)12 +krhuk1k hck1 +k huk1kck1
 
 "k hck21 +C(u;c)
 
k( hc)0k2 +k hck21 +k huk21
+k hukk huk1kchkkchk1 +k hck2 +k hck2
 
:
(4.73)
101
Taking the sum of (4.72) and (4.73), applying Young?s inequality again we obtain
d
dtk 
h
uk
2 +k h
uk
2
1 +
d
dtk 
h
ck
2 +k h
ck
2
1
 "(k huk21 +k hck21)
+C
 
k( hu)0k2 +k( hc)0k2 +k huk21 +k hck21 +k hck2 +k hck2 +k hpk2
+k hukk huk1kuhkkuhk1 +k hukk huk1kchkkchk1
 
 "
 
1 +kchk2 +kuhk2
 
k huk21 +"k hck21
+C
 
k( hu)0k2 +k( hc)0k2 +k huk21 +k hck21
+k hpk2 +k hck2 +k ck2 + kuhk21 +kchk21 k huk2
 
:
(4.74)
Proceeding as in Lemma 4.4, we can prove the uniform boundedness of kuhk and kchk with
respect to h. Choose "> 0 such that
"< minf 11 +ku
hk2 +kchk2
;1g:
Then (4.74) leads to
d
dt
 
k uk2 +k ck2
 
 C
 
k( hu)0k2 +k( hc)0k2 +k huk21 +k hck21 +k hck2 +k hpk2
 
+C
 
1 +kuhk21 +kchk21
  
k uk2 +k ck2
 
:
Integrating with respect to s from 0 to t gives
k huk2 +k hck2  k huk2(0) +k hck2(0)
+C
Z t
0
 
k( hu)0k2 +k( hc)0k2 +k huk21 +k hck21 +k hck2 +k hpk2
 
ds
+C
Z t
0
 
1 +kuhk21 +kchk21
  
k huk2 +k hck2
 
ds:
102
Applying Gronwall?s inequality we have
k huk2 +k hck2  
 
k huk2(0) +k hck2(0)
 
exp
 
C
Z t
0
(1 +kuhk21 +kchk21) ds
 
+C
Z t
0
 
k( hu)0k2 +k( hc)0k2 +k huk21 +k hck21 +k hck2 +k hpk2
 
ds
 C
 
kUh u0k2 +kCh c0k2
 
+
Z t
0
 
k 0uk2 +k 0ck2 +k uk21 +k ck21 +k ck2 +k hpk2
 
ds);
(4.75)
since kuhk21 and kchk21 are integrable on I according to a similar argument as in Lemma
4.4.
Remark 4.9. The proof of the convergence of pressure ph remains open.
4.3 Numerical experiments
In this section we describe a numerical experiment in two dimensions to verify the
convergence of numerical scheme (4.54). In this experiment, the parameters used were as in
Section 3.3. Construct the Taylor-Hood element and assign the parameters as
 = 0:1; U = 0:1;  = 1;  = 10%;
and
 (x) = sin2x+ 1; x2R:
The forcing term f was chosen so that the exact solution is
8
>>>
>><
>>>>
>:
u =pt(sin xsin y;sin xsin y)T ;
p =ptsin xsin y;
c =ptsin xsin y:
103
Table 4.1: Convergence rate in L2(I;L2( ))
h ku uhk kp phk kc chk
1/2 0.0067 0.2386 0.0128
1/4 7.72e-04 0.0531 0.0014
1/8 8.84e-05 0.0128 1.75e-04
1/16 1.16e-05 0.0030 2.20e-05
1/32 1.51e-06 7.23e-04 3.06e-06
conv. rate 2.94 2.05 2.85
Table 4.2: Convergence rate in L2(I;H1( ))
h ku uhk1 kp phk1 kc chk1
1/2 0.1065 0.5044 0.0415
1/4 0.0264 0.2298 0.0103
1/8 0.0066 0.1070 0.0027
1/16 0.0016 0.0474 6.85e-04
1/32 0.0004 0.0233 1.71e-04
conv. rate 2.00 1.02 2.00
Since our focus is on the convergence of the numerical solution with respect to the mesh size
h, we used the time step k = 10 5. The numerical errors for di erent mesh sizes are shown in
table 4.1 and 4.2, from which we can see that the error tends zero as h become smaller just as
stated in Theorem 4.8. Furthermore the numerical method achieves the optimal convergence
order although we did not prove it.
4.4 Conclusion
In chapter 3 and 4, we studied the mathematical model of bioconvection caused by
average upswimming micro-organisms. The PDE system consists of a Navier-Stokes type
equation for the velocity and pressure coupled with a parabolic equation for the concentra-
tion. The viscosity is assumed to be dependent on the concentration. We established the
existence and uniqueness of a weak solution of both steady and evolutionary bioconvection.
We then constructed  nite element approximation of the weak solutions and proved the
104
convergence of the numerical approximation to the exact solution. We used our numerical
method to simulate the velocity and concentration distribution inside a small container. The
uniqueness result and the convergence theorem for the evolutionary case are only valid in
the two dimensional case. The convergence of the discrete pressure remains open for the
time dependent bioconvection.
105
Chapter 5
Conclusion and future work
Conclusion. We studied two systems of partial di erential equations with coe cients
that depend on the solution. The  rst, a quasi-static poroelasticity is a system comprised
of the equation of linear elasticity and and a nonlinear di usion equation, and the second,
a Navier-Stokes type system, both systems were studied using the modi ed Rothe?s method
to handle the solution dependent coe cients. We established the existence and uniqueness
of solutions, and conducted numerical experiments approximating weak solutions of both
systems using the  nite element method . Numerical simulations were constructed to show
the accuracy of the models.
Future work. Many equations remain unresolved and various extensions to the present
work are possible.
1. More e cient numerical methods may be considered, including  nite volume method
and discontinuous Galerkin methods.
2. The equation of quasi-static poroelasticty can be considered subject to general bound-
ary conditions. Introducing general boundary condition may change the null space of oper-
ator B de ned in chapter 2. The energy estimate may need to be modi ed.
3. We may consider the fully dynamic model (1.9) with secondary consolidation, that
is, with the two terms  @
2
@t2 u and   r
d
dt(r u). Both cases result in coupled systems of
hyperbolic and parabolic equations.
4. Many poroelasticity problems have multi-scale features. For instance, we consider a
rock with small pores. On the macro-scale we consider a deformation equation of the rock
while on the micro-scale we solve a di usion equation inside the pores. In Biot?s model
the two equations are coupled through the Biot-Willis constant  . However, in view of
106
multi-scale  nite element method, the two equations can be coupled through mathematical
homogenization by constructing multi-scale  nite elements. More work and details can be
found in [64].
5. System (1.19) can be extended to equations that model the convection caused by the
admixture in the atmosphere and ocean which is important in the study of earth ecology
and which involves large scale computation and multi-species creatures [34].
107
Bibliography
[1] A.  Zen   sek. Nonlinear elliptic and evolution problems and their  nite element approxi-
mations. Computational Mathematics and Applications. Academic Press Inc. [Harcourt
Brace Jovanovich Publishers], London, 1990. With a foreword by P.-A. Raviart.
[2] R. E. Showalter. Monotone operators in Banach space and nonlinear partial di erential
equations, volume 49 of Mathematical Surveys and Monographs. American Mathematical
Society, Providence, RI, 1997.
[3] A. Pazy. Semigroups of linear operators and applications to partial di erential equations.
Springer-Verlag, New York, 1983.
[4] Susanne C. Brenner and L. Ridgway Scott. The Mathematical Theory of Finite Element
Methods, volume 15 of Texts in Applied Mathematics. Springer-Verlag, second edition,
2002.
[5] Vidar Thom ee. Galerkin  nite element methods for parabolic problems, volume 25 of
Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition,
2006.
[6] Phillip Joseph Phillips and Mary F. Wheeler. A coupling of mixed and continuous
Galerkin  nite element methods for poroelasticity. I. The continuous in time case. Com-
put. Geosci., 11(2):131{144, 2007.
[7] Jean-Louis Auriault, Claude Boutin, and Christian Geindreau. Homogenization of Cou-
pled Phenomena in Hetrogenous Media. ISTE Ltd, London, 2009.
[8] Olivier Coussy. Poromechanics. John Wiley & Sons Ltd, Chichester, 2004.
[9] Olivier Coussy. Mechanics and Physics of Porous Solids. John Wiley & Sons Ltd,
Chichester, 2010.
[10] Herbert F. Wang. Theory of Linear Poroelasticity with Applications to Geomechanics
and Hydrology. Princeton Series in Geophysics. Princeton University Press, Princeton,
2000.
[11] N. Lubick. Modeling complex, multiphase porous media systems. SIAM News, 5(3),
2002.
[12] H. Bryne and L. Preziosi. Modeling tumour growth using the theory of mixtures.
Mathematical Medicine and Biology, 2003.
108
[13] T. Roose, P.A. Netti, L. Munn, Y. Boucher, and R. Jain. Solid stress generated by
spheroid growth estimated using a linear poroelastic model. Microvascular Research,
66(204{212), 2003.
[14] A. Smillie, I. Sobey, and Z. Molnar. A hydro-elastic model of hydrocephalus. Technical
report, Oxford University Computing Laboratory: Numerical Analysis Group, 2004.
[15] J. Hudson, O. Stephansson, J. Andersson, C.-F. Tsang, and L. Ling. Coupled t-h-m
issues related to radioactive waste repository design and performance. International
Journal of Rock Mechanics and Mining Sciences, 38(143{161), 2001.
[16] J.-M. Kim and R. Parizek. Numerical simulation of the noordbergum e ect resulting
from groundwater pumping in a layered aquifer system. Journal of Hydrology, 202:231{
243, 1997.
[17] K. Terzaghi. Principle of soil mechanics. Eng. New Record, A Series of Articles, 1925.
[18] L. Rendulic. Prenzi er und porenwasserdruck in tonen. Der Bauingenieur, 17:559{564,
1936.
[19] M. A. Biot. General theory of three-dimensional consolidation. J. Appl. Phys.,
12(2):155{164, 1941.
[20] M. A. Biot. Theory of elasticity and consolidation for a porous anisotropic solid. J.
Appl. Phys., 26(2):182{185, 1955.
[21] M. A. Biot. General solutions of the equations of elasticity and consolidation for a
porous maaterial. J. Appl. Mech., Trans. ASME, 23:91{96, 1956.
[22] M. A. Biot. Theory of propagation of elastic waves in a  uid-saturated porous solid. I.
Low-frequency range. J. Acoust. Soc. Am., 28(2):168{178, 1956.
[23] M. A. Biot and D. G. Willis. The elastic coe cients of the theory of consolidation. J.
Appl. Mech., Trans. ASME, 24:594{601, 1957.
[24] M. A. Biot. Mechanics of deformation and acoustic propagation on porous media. J.
Appl. Phys., 33(4):1482{1498, 1962.
[25] R. E. Showalter. Di usion in poro-elastic media. J. Math. Anal. Appl., 251(1):310{340,
2000.
[26] M. A. Murad and J. H. Cushman. Multiscale  ow and deformation in hydrophilic
swelling porous media. Internat. J. Engrg. Sci., 34:M. A. Murad and J. H. Cushman,
Multiscale  ow and deformation in hydrophilic swelling 313, 1996.
[27] Phillip Joseph Phillips. Finite element methods in linear poroelasticity: Theoretical
and computational results. ProQuest LLC, Ann Arbor, MI, 2005. Thesis (Ph.D.){The
University of Texas at Austin.
109
[28] Phillip Joseph Phillips and Mary F. Wheeler. A coupling of mixed and discontinuous
Galerkin  nite-element methods for poroelasticity. Comput. Geosci., 12(4):417{435,
2008.
[29] P.C. Carman. Permeability of saturated sands, soils and clays. Journal of Agricultural
Science, 29:263{273, 1939.
[30] Robert Wayne Carroll and Ralph E. Showalter. Singular and degenerate Cauchy prob-
lems. Academic Press [Harcourt Brace Jovanovich Publishers], New York, 1976. Math-
ematics in Science and Engineering, Vol. 127.
[31] Kenneth L. Kuttler, Jr. Time-dependent implicit evolution equations. Nonlinear Anal.,
10(5):447{463, 1986.
[32] Spiegel E.A. Levandowsky M., Hunter W.S. A mathematical model of pattern formation
by swimming microorganisms. J. Protozoology 22, pages 296{306, 1975.
[33] Moribe Y. On the bioconvection of tetrahymena pyriformis. Mastersthesis(inJapanese),
OsakaUniversity., 1973.
[34] G.V. Alekseev and E.A. Adomavichus. Theoretical analysis of inverse extremal prob-
lems of admixture di usion in viscous  uids. Journal of Inverse & Ill-Posed Problems,
9(5):435, 2001.
[35] Yukio Kan-on, Kimiaki Narukawa, and Yoshiaki Teramoto. On the equations of bio-
convective  ow. J. Math. Kyoto Univ., 32(1):135{153, 1992.
[36] Akira Harashima, Masataka Watanabe, and Issei Fujishiro. Evolution of bioconvection
patterns in a culture of motile  agellates. Physics of Fluids, 31(4):764{775, 1988.
[37] R. Childress, S. & Peyret. A numerical study of two-dimensional convection by motile
particles. J. M?ec, 15:753{779, 1976.
[38] S. Ghorai and N. A. Hill. Development and stability of gyrotactic plumes in bioconvec-
tion. Journal of Fluid Mechanics, 400:1{31, December 1999.
[39] S. Ghorai and N. A. Hill. Periodic arrays of gyrotactic plumes in bioconvection. Physics
of Fluids, 12(1):5{22, 2000.
[40] S. Ghorai and N. Hill. Wavelengths of gyrotactic plumes in bioconvection. Bulletin of
Mathematical Biology, 62:429{450, 2000. 10.1006/bulm.1999.0160.
[41] MATTHEW M. HOPKINS and LISA J. FAUCI. A computational model of the collec-
tive  uid dynamics of motile micro-organisms. Journal of Fluid Mechanics, 455:149{174,
2002.
[42] Albert Einstein. A new determination of molecular dimensions. Annalen der Physik,
19:289{306, 1906.
110
[43] G. K. Batchelor and J. T. Green. The determination of the bulk stress in a suspension
of spherical particles to order c2. Journal of Fluid Mechanics, 56(03):401{427, 1972.
[44] M. Mooney. The viscosity of a concentrated suspension of spherical particles. Journal
of Colloid Science, 6(2):162 { 170, 1951.
[45] Irvin M. Krieger and Thomas J. Dougherty. A mechanism for non-newtonian  ow in
suspensions of rigid spheres. Transactions of the Society of Rheology, 3(1):137{152,
1959.
[46] John F. Brady. The rheological behavior of concentrated colloidal dispersions. The
Journal of Chemical Physics, 99(1):567{581, 1993.
[47] B. Climent-Ezquerra, L. Friz, and M. Rojas-Medar. Time-reproductive solutions
for a bioconvective  ow. Annali di Matematica Pura ed Applicata, pages 1{20.
10.1007/s10231-011-0245-7.
[48] Jozef Ka cur. Method of Rothe in evolution equations, volume 80 of Teubner-Texte zur
Mathematik [Teubner Texts in Mathematics]. BSB B. G. Teubner Verlagsgesellschaft,
Leipzig, 1985. With German, French and Russian summaries.
[49] Karel Rektorys. The method of discretization in time and partial di erential equations,
volume 4 of Mathematics and Its Applications (East European Series). D. Reidel Pub-
lishing Co., Dordrecht, 1982. Translated from the Czech by the author.
[50] A. J. Meir Y. Cao, S. Chen. Steady  ow in a deformable porous medium. Mathematical
Methods in Applied Sciences, 2013, Coming Soon.
[51] J.-L. Lions. Optimal control of systems governed by partial di erential equations. Trans-
lated from the French by S. K. Mitter. Die Grundlehren der mathematischen Wis-
senschaften, Band 170. Springer-Verlag, New York, 1971.
[52] Roger Temam. Navier-Stokes equations. AMS Chelsea Publishing, Providence, RI,
2001. Theory and numerical analysis, Reprint of the 1984 edition.
[53] C. Taylor and P. Hood. A numerical solution of the Navier-Stokes equations using the
 nite element technique. Internat. J. Comput. & Fluids, 1(1):73{100, 1973.
[54] John G. Heywood and Rolf Rannacher. Finite-element approximation of the nonsta-
tionary navier-stokes problem part ii: stability of solutions and error estimates uniform
in time. SIAM Journal, 23:750{777, 1986.
[55] John G. Heywood and Rolf Rannacher. Finite-element approximation of the nonsta-
tionary navier-stokes problem part iv:error analysis for second-order time discretization.
SIAM Journal, 27:353{384, 1990.
[56] Michael Renardy and Robert C. Rogers. An Introduction to Partial Di erential Equa-
tions, volume 13 of Texts in Applied Mathematics. Springer-Verlag, New York, second
edition, 2004.
111
[57] J.-L. Lions. Quelques m ethodes de r esolution des probl emes aux limites non lin eaires.
Dunod, 1969.
[58] Lawrence C. Evans. Partial di erential equations, volume 19 of Graduate Studies in
Mathematics. American Mathematical Society, Providence, RI, second edition, 2010.
[59] J. Kozeny. Ueber kapillare leitung des wassers im boden. Sitzungsber Akad. Wiss.,
Wien, 136(2a):271{306, 1927.
[60] P.C. Carman. Fluid  ow through granular beds. Transactions, Institution of Chemical
Engineers, London, 15:150{166, 1937.
[61] P.C. Carman. Flow of gases through porous media. Butterworths, London, 1956.
[62] Vivette Girault and Pierre-Arnaud Raviart. Finite element methods for Navier-Stokes
equations, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag,
Berlin, 1986. Theory and algorithms.
[63] Max D. Gunzburger. Finite element methods for viscous incompressible  ows. Computer
Science and Scienti c Computing. Academic Press Inc., Boston, MA, 1989. A guide to
theory, practice, and algorithms.
[64] Weinan E. and Bjorn Engquist. The heterogeneous multiscale methods. COMM. MATH.
SCI, 1:87{132, 2003.
112