Finite Difference Schemes for Poro-elastic ProblemS

COMPUTATIONAL METHODS IN APPLIED MATHEMATICS, Vol.2(2002), No.2, pp.132–142 c Institute of Mathematics of the National Academy of Sciences of Belarus °

FINITE DIFFERENCE SCHEMES FOR PORO-ELASTIC PROBLEMS 1 FRANCISCO J. GASPAR Departamento de Matem´atica Aplicada, Universidad de Zaragoza 50009 Zaragoza, Spain E-mail: [email protected]

FRANCISCO J. LISBONA Departamento de Matem´atica Aplicada, Universidad de Zaragoza 50009 Zaragoza, Spain E-mail: [email protected]

PETR N. VABISHCHEVICH Institute for Mathematical Modelling RAS 4-A, Miusskaya Sq., 124047 Moscow, Russia E-mail: [email protected] Abstract — In this paper, we present a finite difference analysis of the consolidation problem for saturated porous media. In the classical model, the behaviour of the porous environment – fluid system is described by a set of equations for the unknown vector displacements of the matrix skeleton and the fluid pressure. For simplicity we consider a model problem with constant coefficients in a rectangular domain. A priori estimates for the difference solution of the problem are obtained and on their basis the convergence of two-level difference schemes is investigated. 2000 Mathematics Subject Classification: 65M06; 65M12. Keywords: Biot’s model, finite difference method, stability estimates.

Introduction Filtration and consolidation problems are concerned with the response of a saturated porous material under some modification of the equilibrium conditions. For instance, pumping water from an aquifer produces a reduction of the water pressure and, therefore, a growth in the effective stress in the solid skeleton, which results in a consolidation process. Also, an increase in the loads on the ground surface leads to an increase in the stress, which results in a deformation and time-dependent consolidation process associated with the drainage of the pore fluid [1, 14]. 1

This research has been partially supported by the Spanish project MCYT-FEDER BFM 2001-2521 and The Russian Foundation for Basic Research RFBR 99–01–00958 Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

Finite Difference Schemes in Poro-Elasticity

133

A simple mathematical model was first proposed in [18] with a one-dimensional dissipation equation for the pore pressure without simultaneous consideration of the strain in the solid matrix. Another important approach is the one proposed by Biot ( [2–4]). Here the strain and the flow are considered simultaneously in a three-dimensional domain and the resulting evolutionary problem consists of a set of equations for the unknown vector of displacements of the solid skeleton and the pore pressure of the fluid. Models for the description of the poro-elastic problems under more general conditions were proposed, for example, in [7], where the problems with the nonlinear filtration law and rheological properties are considered. The existence and uniqueness of the solution for the initial-boundary value problem for the system of equations given by Biot’s model have been investigated in several papers [17,22]. The problem of the well-posedness for nonlinear models is considered, for example, in [5]. The numerical solution of the poro-elastic problems by the finite element method is considered for the incompressible Biot’s model in several papers, see for example [9, 21]. For a more comprehensive analysis see [11–13]; these papers discuss the error estimates for Euler-Galerkin discretizations using stable and unstable combinations of finite element spaces of displacement and pore pressure fields and the short- and long-time behaviour of such approximations. Finite difference methods for the flow and deformation problems are used in [10], without analyzing the convergence. For the most common nonlinear problems, with incompressible fluids, the difference schemes are investigated in [6]. In this paper, the construction of difference schemes for the numerical solution of the poroelastic problems for a slightly compressible fluid is considered. For technical simplicity in the finite difference analysis, the proposed model considers only constant coefficient equations in a rectangular domain. A priori estimates (stability estimates) for the semi-discrete and fully discrete difference solution of the problem are obtained using the appropriate Hilbert norms. These estimates and the approximation error bounds of the proposed schemes lead to convergence.

1. Continuous problem Let us consider a domain Ω ⊂ R2 (a two-dimensional skeleton matrix) with a piecewise smooth boundary ∂Ω and unit outward normal n. We use the notation x = (x1 , x2 ) for the generic point in Ω, v = (v1 , v2 ) for the displacement vector, and p for the pore pressure of the fluid. We shall consider the elementary consolidation problem for a saturated, homogeneous and isotropic porous material. The classical quasi-static Biot model is based on the assumption of incompressibility for the filtering fluid. Here, we consider a different situation where the fluid is slightly compressible. Under these conditions, the process of filtration and consolidation is described in [1, 17]. Neglecting the body forces, the problem is governed by the set of equations ˜ − (λ + µ)grad div v + grad p = 0, −µ∆v

(1)

∂ (a p + div v) − χ∆p = f (x, t), x ∈ Ω, 0 < t 6 T , (2) ∂t where λ and µ are the Lam´e coefficients; a = nf β 6= 0, with nf the porosity and β the compressibility coefficient of the fluid; χ = k/µf , being k the permeability of the porous Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

134 F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich environment, µf the viscosity of the fluid and µ ¶ ∆ 0 ˜ ∆= . 0 ∆ The source term f (x, t) is used, for example, when a forced filtration process is present. For simplicity, we assume that ∂Ω is rigid and permeable (free drainage condition) so that we can impose the homogeneous Dirichlet boundary conditions v(x, t) = 0 ,

p(x, t) = 0,

x ∈ ∂Ω.

(3)

The initial condition can be given by a p(x, 0) + div v (x, 0) = s(x),

x ∈ Ω,

(4)

that is just the fluid content in the system at the initial time. To formulate the operator for (1)–(4), we first introduce appropriate functional spaces and operators. Let H = L2 (Ω) be the set of square-integrable scalar valued functions defined on Ω, with the scalar product and the corresponding norm Z (u, v) = u(x)v(x)dx, kuk = (u, u)1/2 . Ω

Let V denote the usual Sobolev space H01 (Ω) of functions vanishing at the boundary ∂Ω equipped with the inner product and norm 2 Z X ∂u ∂v (u, v)1 = dx , ∂xα ∂xα α=1

kuk1 = (u, u)1/2 .

Ω

e = For the two-dimensional vector-valued functions u, v we shall set the Hilbert space H H ⊕ H with the inner product and norm (u, v) = (u1 , v1 ) + (u2 , v2 ),

kuk = (ku1 k2 + ku2 k2 )1/2 ,

e = V ⊕ V, with and the Hilbert space V (u, v)1 = (u1 , v1 )1 + (u2 , v2 )1 ,

kuk1 = (ku1 k21 + ku2 k21 )1/2 .

e we consider the elasticity operator A by In H e − (λ + µ)grad div v, Av = −µ∆v

(5)

e | Av ∈ H} e . The operator A is positive and self-adjoint with the domain D(A) = {v ∈ V (A > 0, A = A∗ ). Also, the following energetic equivalence holds: ˜ v) 6 (Av, v) 6 −(λ + 2µ)(∆v, ˜ v). −µ(∆v,

(6)

e where Ee is the identity operator on V e and δ0 is the minimum In view of this A = A∗ > µδ0 E, eigenvalue of −∆. Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

Finite Difference Schemes in Poro-Elasticity

135

Similarly, in H we define Bp = −χ∆p,

(7)

with the domain D(B) = {p ∈ V | Bp ∈ H} . The operator B is positive and satisfies B = B ∗ > χδ0 E, where E now denotes the identity operator on H . In constructing discrete analogs for A and B, we will be oriented to the fulfillment of the same properties. Problem (1)–(4) can be written in differential operator form as the abstract initial value e×H problem in H Av + grad p = 0,

(8)

d (a p + div v) + Bp = f (t), dt a p(0) + div v(0) = s,

(9) (10)

with s ∈ H and f (t) ∈ H, ∀ t 6 T . We assume T 1 that problem (8)-(10) has a solution (v, p), 0 for which ap(·) + divv(·) ∈ C ([0, T ], H) C ((0, T ], H). Moreover in the analysis of the numerical approximation below, we take the classical solution as regular as we need. We associate to the operators A and B the inner products e (y, w)A = (Ay, w), y, w ∈ H, (y, w)B = (By, w), y, w ∈ H, and the corresponding norms kykA and kykB . We use the same association for the selfadjoint and definite positive operator B −1 . By the properties of A and B it is easy to prove the following a priori estimate for the solution of (8), (10): kv(t)k2A

2

+ akp(t)k 6

kv(0)k2A

1 + akp(0)k + 2

Zt

2

kf (θ)k2B−1 dθ,

(11)

0

for t 6 T , which ensures stability with respect to the initial data and the right-hand side. To obtain estimates for the pressure gradient field, we use the equations (Aut , ut ) + (Gpt , ut ) = 0, a(pt , pt ) + (Dut , pt ) + (Bp, pt ) = (f, pt ),

(12) (13)

from which it follows that 1d 1 α kpk2B 6 kf k2 + kpt k, 2 dt 2α 2 As kut kA > 0, taking α = 2a, we have the stability estimate kut k2A + akpt k2 +

k p(t)

k2B 6k

p(0)

k2B

1 + 4a

∀α > 0.

Zt kf (x, θ)k2 dθ.

(14)

0

At the initial stage, the displacements and the pressure satisfy the system of equations ˜ 0 − (λ + µ)grad div v0 + grad p0 = 0, −µ∆v a p0 + div v0 = s(x), x ∈ Ω,

(15) (16)

with the boundary condition v0 (x) = 0 on ∂Ω. The solution of this problem gives the values for the displacements and pressure at t = 0 required in estimate (11). Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

136 F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich

2. Space Discretization In considering of difference schemes for the solution of problem (1) – (4), we begin with making space approximation. We consider the problem on the rectangle Ω = {x | x = (x1 , x2 ), 0 < xα < lα , α = 1, 2} discretized by a uniform rectangular grid with mesh steps hα , α = 1, 2. Let ω be the set of internal nodes of the grid ω ={x | x = (x1 , x2 ), xα = iα hα , Nα hα = lα , α = 1, 2},

iα = 1, 2..., Nα − 1,

and ∂ω – the set of boundary nodes. The finite difference solution of problem (1) - (4) will be denoted by vh (x, t), ph (x, t) , x ∈ ω ∪ ∂ω, 0 < t 6 T . Using the standard index-free notation of the theory of difference schemes [16], for the right and left difference derivatives we write wx =

w(x + h) − w(x) , h

wx¯ =

w(x) − w(x − h) , h

and the second difference derivative is given by the expression wx¯x =

1 w(x + h) − 2w(x) + w(x − h) (wx − wx¯ ) = . h h2

For the grid functions vanishing at ∂ω, we define the Hilbert space H = L2 (ω), with the scalar product and norm given by X (y, w) = ywh1 h2 , kyk = (y, y)1/2 . x∈ω

e = H ⊕ H for vector-valued functions that Similarly to the differential case, we introduce H are null at ∂ω, with (y, w) = (y1 , w1 ) + (y2 , w2 ),

kyk = (y, y)1/2 .

Given a self-adjoint and definite positive operator C , HC denotes the space H provided by the scalar product (y, w)C = (Cy, w) and norm kyk = (Cy, y)1/2 . We approximate the differential operator A by the difference operator ¶ µ A11 A12 , (17) A= A21 A22 A11 y = −µ∆h y − (λ + µ) yx¯1 x1 , λ+µ (yx¯1 x2 + yx1 x¯2 ), A12 y = A21 y = − 2 A22 y = −µ∆h y − (λ + µ) yx¯2 x2 , where the usual five-point stencil approximation of the Laplace operator ∆h y = yx¯1 x1 + yx¯2 x2 Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

Finite Difference Schemes in Poro-Elasticity

137

is used. e we have For y(x), w(x) ∈ H (Ay, w) = (y, Aw), e In addition, we have i.e., A is self-adjoint in H. ˜ h 6 A 6 −(λ + 2µ)∆ ˜ h, −µ∆ where

µ ˜h = ∆

∆h 0 0 ∆h

(18)

¶ .

Relation (18) is a discrete analog of (6) given for the differential operator A. Due to this e where δh > 0 is the minimum eigenvalue of operator ∆h and E e is the A = A∗ > µδh E, e i.e. the difference operator A is, like the differential operator A, a identity operator in H, self-adjoint and positive definite operator. To approximate the diffusion operator B, we use By = −χ∆h y,

(19)

so that B = B ∗ > χδh E where E is the identity operator in H. The approximation of the coupling terms deserves separate consideration. For the funce and p ∈ H, tions v ∈ H (v, grad p) = −(p, div v), i.e. the divergence operator coincides with the opposite conjugate of the gradient operator. Therefore, we will approximate these operators so that the same property takes place at the difference level (vh , gradh ph ) = −(ph , divh vh ). (20) Several types of this kind of approximations for the gradient and divergence are widely used in the numerical simulation of viscous incompressible flows (see, for example [8, 15, 19, 20]). We now report the main possibilities to be used when the same grid is used for components of the displacement vector and pressure. The first variant is connected with one-sided differences. For example, if forward differences are used for the gradient gradh y = (yx1 , yx2 ),

x ∈ ω,

for y ∈ H,

the discrete divergence is divh w = (w1 )x¯1 + (w2 )x¯2 ,

x ∈ ω,

e for w ∈ H.

Also, it is possible to use second-order approximations for the gradient and divergence by taking gradh y = (yx◦ , yx◦ ), x ∈ ω, for y ∈ H, 1

2

divh w = (w1 )x◦ + (w2 )x◦ , x ∈ ω, 1

2

e for w ∈ H.

After the space approximation, from (8), (9) we arrive at the Cauchy problem for the system of differential-difference equations Avh + gradh ph = 0, Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

(21)

138 F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich d (aph + divh vh ) + Bph = fh (x, t), dt with the initial condition a ph (0) + divh vh (0) = sh (x),

x ∈ ω,

(22)

x ∈ ω.

(23)

We now construct a simple difference scheme for the approximation of the solution {vh , ph } of the Cauchy problem for system (21), (22).

3. Fully Discrete Approximation We use a uniform grid for time discretization with a step-size τ > 0. Let y n (x) = y(x, tn ), where tn = nτ, n = 0, 1, . . . , N, N τ = T . We consider the two-level scheme with weight σ Avhn+1 + gradh pn+1 = 0, h a

n = 0, 1, . . . , N − 1,

pn+1 − pnh divh vhn+1 − divh vhn h + + B(σpn+1 + (1 − σ)pnh ) h τ τ = fh (x, σtn+1 + (1 − σ)tn ), n = 0, 1, 2, . . . , N − 1,

(24) (25) x ∈ ω.

Under standard restrictions for σ, stability of the difference scheme (24), (25) can be established. More precisely, the following result holds. Proposition 1. For σ > 0.5 the solution of the difference scheme (24), (25) satisfies the a priori estimate τ 2 n 2 n 2 2 kvhn+1 k2A + akpn+1 h k 6 kvh kA + akph k + kfh (x, σtn+1 + (1 − σ)tn )kB −1 . 2

(26)

Proof. Let us denote n+1 pn+1 + (1 − σ)pnh h,σ = σph

and n+1 vh,σ = σvhn+1 + (1 − σ)vhn .

From (24), ( 25 ) we have for 0 6 n < N µ ¶ µ ¶ n+1 n+1 − vhn − vhn n+1 vh n+1 vh + gradh ph,σ , = 0, Avh,σ , τ τ and

µ a

¶ µ ¶ − pnh n+1 pn+1 div vhn+1 − div vhn n+1 h , ph,σ + , ph,σ τ τ n+1 n+1 n+1 + (Bpn+1 h,σ , ph,σ ) = (fh,σ , ph,σ ).

The addition of (27) and (28) yields µ ¶ ¶ µ n+1 n+1 − vhn ph − pnh n+1 n+1 vh n+1 n+1 n+1 Avh,σ , , ph,σ + (Bpn+1 +a h,σ , ph,σ ) = (fh,σ , ph,σ ) τ τ and by

1 n+1 n+1 n+1 , ph,σ ) 6 kpn+1 (fh,σ h,σ kB + kfh,σ kB −1 , 4 Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

(27)

(28)

Finite Difference Schemes in Poro-Elasticity

139

we obtain

τ n+1 kf kB −1 . (29) 4 h,σ Now, using the identity σξ + (1 − σ)ζ = 1/2(ξ + ζ) + (σ − 1/2)(ξ − ζ) in the expression of n+1 vh,σ and of pn+1 h,σ , we have the inequality n+1 n+1 (Avh,σ , vhn+1 − vhn ) + a(pn+1 − pnh ) 6 h,σ , ph

a 1 2 n 2 (kvhn+1 k2A − kvhn k2A ) + (kpn+1 h k − kph k ) 2 2 µ ¶ 1 τ n+1 + σ− (kvhn+1 − vhn k2A + akpn+1 − pnh k2 ) 6 kfh,σ kB −1 , h 2 2

(30)

and if σ > 0.5, estimate (26) follows. Estimate (26) provides a difference analogue of (11) for the solution of problem (21), (22). To derive the error bounds for the difference scheme, we apply the methodology of analyzing the of finite difference methods. We consider the difference equations for the error in the displacements and pressure δvhn (x) = vhn (x) − v(x, tn ), δpnh (x) = pn (x) − p(x, tn ), x ∈ ω , Aδvhn+1 + gradh δpn+1 = ψ1n+1 (x), h a

δpn+1 h

− τ

δpnh

divh δvhn+1

n = 0, 1, . . . , N − 1,

− τ + B(σδpn+1 + (1 − σ)δpnh ) = ψ2n+1 (x), h +

(31)

divh δvhn

(32) n = 0, 1, . . . , N − 1,

x ∈ ω.

The right-hand side terms in equations (31) and (32) are respectively the approximation errors for equations (1) and (2). For smooth solutions, we have ψ1n+1 (x) = −Av(x, tn+1 ) − gradh p(x, tn+1 ) = O(|h|α )

(33)

and p(x, tn+1 ) − p(x, tn ) τ divh v(x, tn+1 ) − divh v(x, tn ) − τ − B(σp(x, tn+1 ) + (1 − σ)p(x, tn+1 )) = O(τ ν + |h|α ),

ψ2n+1 (x) =f (x, σtn+1 + (1 − σ)tn ) − a

(34)

p where |h| = h21 + h22 , α = 1 if the gradient and divergence operators are approximated by the one-sided differences or α = 2 if the central differences are used, and ν = 2 if σ = 0.5 or ν = 1 if σ 6= 0.5. Attempting to apply (25) in order to estimate errors δvhn (x) and δpn (x), 1 6 n 6 N , n n 0 n we split the displacement error δvhn (x) = w1,h (x) + w2,h (x), where w1,h (x) = 0 and w1,h (x), 1 6 n 6 N is the solution of n (35) (x) = ψ1n (x). Aw1,h This part of the error satisfies n kw1,h kA 6 kψ1n kA−1 ,

n = 0, 1, . . . , N.

Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

(36)

140 F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich n In this situation, w2,h (x) and δpn (x) are the solutions of problem n+1 Aw2,h (x) + gradh δpn+1 = 0, h

n = 0, 1, . . . , N − 1,

(37)

n+1 n (x) δpn+1 − δpnh divh w2,h (x) − divh w2,h h a + + B(σδpn+1 + (1 − σ)δpnh ) h τ τ n+1 n div h w1,h (x) − divh w1,h (x) n+1 = ψ2 (x) − , τ n = 0, 1, . . . , N − 1, x ∈ ω,

(38)

where the right-hand side is null in equation (37). If σ > 0.5, from (26) we obtain n+1 2 0 2 2 2 0 kw2,h kA + akδpn+1 h k 6 kw2,h kA + akph k  ° °2  n k+1 k ° ° τ X  k+1 2 ° divh w1,h − divh w1,h °  + kψ2 kB −1 + ° , ° ° ° −1 2 τ k=0

(39)

B

n = 0, 1, . . . , N − 1,

x ∈ ω.

From (34) one can straightforward by obtain kψ2k+1 kB −1 = O(τ ν + |h|α ). The estimate for ° ° ° div wk+1 − div wk ° h 1,h ° ° h 1,h Dk+1 = ° ° ° −1 ° τ B

is given by the following lemma. Lemma 1. With the previous notation we have r ° n+1 ° n° 2 ° 1 ψ − ψ 1 1 n+1 ° ° D 6 ° −1 , χδh µ ° τ A

n = 0, . . . N − 1

(40)

and consequently Dn = O(τ ν + |h|α ), where α = 1 if the gradient and divergence operators are approximated by one-sided differences or α = 2 if the central differences are used and ν = 2 if σ = 0.5 or ν = 1 if σ 6= 0.5. (s1n , s2n )

n+1 n w1,h − w1,h . From (19) the difference derivative τ

Proof. Let us denote by sn = we have 1 kdivh sn k. Dn+1 = kdivh sn kB −1 6 χ δh Now, taking into account that

kdivh sn k2 6 2ks1n k2−∆h + 2ks2n k2−∆h = 2ksn k2−∆e and using (18), we obtain kdivh sn k2 6 Dn+1 By definition, ksn kA

2 ksn k2A . Then µ r 1 2 6 ksn kA . χδ0,h µ ° n+1 ° ° ψ1 − ψ1n ° ° 6° ° ° τ

A−1

and (40) follows. Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

h

Finite Difference Schemes in Poro-Elasticity

141

Proposition 2. Let vh0 (x) and p0h (x), O(|h|α )be approximations of v(x, 0) and p(x, 0). For σ > 0.5 the difference solution of the Cauchy problem (21)–(23) converges to smooth enough solutions of problem (1)–(4).More precisely k δvhn kA +a k δpnh k= O(τ ν + hα ) , where the parameters ν and α are specified in the previous Lemma . Proof. Note that 1 n+1 2 n+1 2 n+1 2 2 kδvhn+1 k2A + akpn+1 h k 6 kw1,h kA + kw2,h kA + akδph k . 2 If σ > 0.5, using (36) and (39), we have 1 2 0 2 0 2 n+1 2 kδvhn+1 k2A + akδpn+1 kA−1 h k 6 kw2,h kA + akδph k + kψ1 2 Ã ° ! n n ° n+1 n °2 X X ° 2τ ψ − ψ 1° ° 1 +τ kψ2k+1 k2B −1 + 2 2 ° ° −1 , χ δ µ τ h k=0 A k=0 and now it is sufficient to use the order estimates for the approximation errors ψ1n , ψ2n and for initial error, to obtain the result. From the previous proposition, the convergence properties for the displacements in the A– norm and for the pressure in the discrete L2 –norm follows. We now consider the convergence of the pressure gradients or equivalently convergence of pressure in the B–norm. From (24)– (25) we have ° n+1 ° ° ° n+1 ° vh − vhn °2 ° ph − pnh °2 ° ° ° + (Bpn+1 , pn+1 − pnh ) = (f n+1 , pn+1 − pnh ), ° τ° h,σ h h,σ h ° + aτ ° ° τ τ

(41)

A

which provides the inequality ° n+1 ° n °2 v − v 1° h h ° ° + a k pn+1 − pnh k2 + 1 (k pn+1 k2B − k pnh k2B ) h h ° ° τ τ τ 2 A µ ¶ 1 α 1 n+1 2 k + k pn+1 − pnh k2 . − pnh k2B 6 k fh,σ + σ− k pn+1 h h 2 2α 2 If σ > 0.5, taking α =

(42)

2a we have τ k2B 6k pnh k2B + k pn+1 h

τ n+1 2 k . k fh,σ 2a

Proposition 3. Let us assume that the same conditions as in proposition 2 hold. For σ > 0.5 we have the estimate k δpnh kB = O(τ ν + hα ), where ν = 2 if σ = 0.5 or ν = 1 if σ > 0.5 and α = 1 if the gradient and divergence operators are approximated by the one-sided differences or α = 2 if the central differences are used. Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM

142 F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich

References [1] J. Bear and Y. Bachmat, Introduction to Modelling Phenomena of Transport in Porous Media, Kluwer Academic, Dordrecht, 1991. [2] M. A. Biot, General theory of three dimensional consolidation, J. Appl. Phys., 12 (1941), pp. 155–169. [3] M. A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid , J. Appl. Phys., 26 (1955), pp. 182–185. [4] M. A. Biot, General solutions of the equations of elasticity and consolidation for a porous material , J. Appl. Phys., 78 (1956), pp. 91–96. [5] R. Z. Dautov, M. I. Drobotenko, and A. D. Lyashko, Study on well-posedness of the generalized solution of the problem of filtration consolidation, Differents. Uravnenia, 33 (1997), pp. 515–521, in Russian. [6] M. I. Drobotenko and A. D. Lyashko, Approximate solution of the problem of filtration consolidation, Izv. Vuzov. Matematika, 358 (1992), No. 3, pp. 3–6, in Russian. [7] A. G. Egorov, A. V. Kosterin, and E. V. Skvortsov, Consolidation and Acoustic Waves in Porous Media, Kazan, 1990, in Russian. [8] F. H. Harlow and J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids, 8 (1965), pp. 2182–2189, in Russian. [9] A. V. Kosterin, M. F. Pavlova, and Y. Shemuranova, Numerical investigation of filtration consolidation, Matem. Modelirovanie, 13 (2001), No. 9, pp. 63–70, in Russian. [10] G. M. Mercer and S. Barry, Flow and deformation in poroelasticity. II: Numerical method , Matem. Modelirovanie, 30 (1999), pp. 31–38. [11] M. A. Murad and A. F. D. Loula, Improved accuracy in finite element analysis of Biot’s consolidation problem, Comput. Methods Appl. Mech. Eng., 95 (1992), pp. 359–382. [12] M. A. Murad and A. F. D. Loula, On stability and convergence of finite element approximations of Biot’s consolidation problem, Int. J. Numer. Methods Eng., 37 (1994), pp. 645–667. [13] M. A. Murad, V. Thom´ee, and A. F. D. Loula, Asymptotic behaviour of semidiscrete finite-element approximations of Biot’s consolidation problem, SIAM J. Numer. Anal., 33 (1996), pp. 1065–1083. [14] V. N. Nikolaevskii, Mechanics of Porous and Cracked Media, Nedra, Moscow, 1984, in Russian. [15] S. V. Patankar, Numerical Heat Transfer and Fluid Flow , Hemispher, Washington, 1980. [16] A. A. Samarskii, Theory of Difference Schemes, Marcel Dekker, New York, 2001. [17] R. E. Showalter, Diffusion in poro-elastic media, J. Math. Anal. Appl., 251 (2000), pp. 310–340. [18] K. Terzaghi, Theoretical Soil Mechanics, Wiley, New York, 1942. [19] P. N. Vabishchevich, A. N. Pavlov, and A. G. Churbanov, Numerical methods for solving time-dependent Navier-Stokes equations in the primitive variables on partially staggered grids, Matem. Modelirovanie, 9 (1997), No. 4, pp. 85–115, in Russian. [20] P. N. Vabishchevich, A. N. Pavlov, and A. G. Churbanov, Solution methods for unsteady incompressible flows in the primitive variables on non-staggered grids, Matem. Modelirovanie, 8 (1999), No. 7, pp. 81–108, in Russian. [21] Y. Yokoo, K. Yamagata, and H. Nagaoka, Finite element method applied to Biot’s consolidation theory, Soils and Foundations, 11 (1971), pp. 9–46. ˇ sek, Finite element method applied to Biot’s consolidation theory, Appl. Math., 29 (1988), pp. [22] A. Zeniˇ 194–210.

Received 26 Feb. 2002 Revised 4 Jul. 2002 Brought to you by | University of Notre Dame Authenticated Download Date | 11/17/14 8:26 PM