On the structure preserving high-order approximation of quasistatic poroelasticity

12/30/2019 ∙ by Herbert Egger, et al. ∙ Technische Universität Darmstadt 0

We consider the systematic numerical approximation of Biot's quasistatic model for the consolidation of a poroelastic medium. Various discretization schemes have been analysed for this problem and inf-sup stable finite elements have been found suitable to avoid spurios pressure oscillations in the initial phase of the evolution. In this paper, we first clarify the role of the inf-sup condition for the well-posedness of the continuous problem and discuss the choice of appropriate initial conditions. We then develop an abstract error analysis that allows us to analyse some approximation schemes discussed in the literature in a unified manner. In addition, we propose and analyse the high-order time discretization by a scheme that can be interpreted as a variant of continuous-Galerkin or particular Runge-Kutta methods applied to a modified system. The scheme is designed to preserve both, the underlying differential-algebraic structure and energy-dissipation property of the problem. In summary, we obtain high-order Galerkin approximations with respect to space and time and derive order-optimal convergence rates. The numerical analysis is carried out in detail for the discretization of the two-field formulation by Taylor-Hood elements and a variant of a Runge-Kutta time discretization. Our arguments can however be extended to three- and four field formulations and other time discretization strategies.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

In linear quasistatic theory [3] the consolidation of a poroelastic solid, which is fully saturated by an incompressible fluid, is usually described by Biot’s equations


together with appropriate initial and boundary conditions. Here is the solid displacement of the porous medium and is the pressure of the residing fluid whereas denotes the density of external forces and is the fluid source density. The Biot parameter is usually close to one and the hydraulic conductivity is assumed strictly positive.

Due to its many applications, e.g., in geosciences or biomathematics, the theoretical and numerical analysis of (1)–(2) has attracted significant interest in the mathematical literature. The existence of unique solutions for the Biot system has been established by Zenisek [24, 25]

under some regularity conditions on the data via discretization with finite elements and the implicit Euler method using a-priori estimates and compactness arguments. Showalter

[21] established well-posedness for different formulations of the Biot model by the method of semi-groups. Murad and Loula [15, 16] investigated the Galerkin approximation by stable and unstable finite element pairs and established decay estimates for the discrete error via improved energy estimates. Mixed finite element approximations of the three-field formulation in which the seapage velocity is introduced as a new variable, were considered by Phillips and Wheeler [18, 19, 20]. Yi [27] investigated the discretization of a four-field formulation in which also the elastic stress field is introduced additionally. In [12] Kanschat and Riviere considered the approximation of the three-field formulation with a non-conforming approximation of the elastic deformation by finite elements. While most of the previous papers only utilized low order approximations in time, Bause et al. [2, 13] considered the efficient implementation of high order time approximations for poroelasticity, but no convergence analysis was conducted. Various available results concern the discretization of the static Biot systems that arise after time discretization, see e.g., [11, 17] in which stability with respect to model and discretization parameters has been studied. It is well-accepted by now [16, 20, 28], that inf-sup stable finite element pairs should be used to avoid spurious pressure oscillations that might appear in the initial phase of the simulations.

In this paper, we consider the systematic construction of high-order approximations for the quasistatic system (1)–(2) by Galerkin methods in space and time. We give a short and concise proof of well-posedness of the continuous problem which motivates our functional analytic setting and provides guidelines for the choice of initial conditions. The regularity requirements for the data are based on the physically relevant energy-dissipation structure of the evolution problem and compatibility conditions for the initial values are derived from the differential-algebraic structure. We then consider the Galerkin approximation in space and discuss the properties of the resulting differential-algebraic equations. In particular, we show that the index of the differential-algebraic system [5, 14] is one, independent of the approximation spaces, while a discrete inf-sup condition, i.e., the surjectivity of the discrete divergence operator, is required to guarantee the well-posedness of the semi-discretization under natural compatibility conditions for the initial conditions. Otherwise, additional non-physical constraints for the initial conditions arise. Afterwards, we show that the discrete error, i.e., the difference between the Galerkin semi-discretization and an appropriate elliptic projection, only depends on the approximation error in the displacement which implies a certain super-convergence for the discrete error. We then consider the time-discretization by a strategy which is capable of preserving both, the particular differential-algebraic structure and the energy-dissipation property of the underlying problem. The resulting scheme can be interpreted as a continuous Galerkin method or a variant of certain Runge-Kutta methods applied to a reformulation of the problem in which the algebraic equation is differentiated. The latter interpretation also allows us for an efficient implementation.

The remainder of the manuscript is organized as follows: In Section 2, we introduce an abstract evolution problem which covers the weak formulation of the Biot system (1)–(2) as a special case, and we establish its well-posedness under mild assumptions on the problem structure and the data. A natural compatibility condition for the initial conditions and the underlying energy-dissipation structure are presented together with the resulting a-priori bounds. In Section 3, we discuss the Galerkin approximation of the abstract problem and investigate the effect of the discrete inf-sup stability on the index of the resulting differential-algebraic system. We then establish an abstract convergence result for inf-sup stable approximations. In Section 4, we propose a time discretization scheme which can be interpreted as an inexact continuous Galerkin approximation or a variant of certain Runge-Kutta collocation methods. In the spirit of [1], this method can be interpreted in a pointwise sense which allows us to show that it preserves the underlying energy-dissipation structure and to conduct the error analysis with similar arguments as on the continuous level. In Section 5, we apply the results to the discretization of the Biot system by Taylor-Hood finite elements and derive order optimal error estimates in space and time. Numerical tests are presented in Section 6 for illustration of our theoretical results, and some directions for possible extensions and further investigations are discussed in the last section.

2. An abstract model problem

We now introduce an abstract model problem which covers the Biot system as a special case, and then establish its well-posedness under mild structural assumptions on the data. In addition, we briefly discuss the regularity of the solution as well as the underlying energy-dissipation identity which us later used to establish a-priori estimates.

2.1. Notation and summary of results

Let , , be real Hilbert spaces with compact embedding . By identifying with its dual , we obtain a Gelfand-triple ; see [7, Chapter XVIII]. The symbol will be used to denote the duality product on , , and .

In the following, we consider an abstract evolution problem stated in weak form as


Here , , and are given bilinear forms and , are functions of time with values in and , respectively. To establish the well-posedness for the corresponding Cauchy problem, we make the following assumptions.

Assumption 1.
  • is bounded, elliptic, and symmetric.

  • is bounded, elliptic, and symmetric.

  • is bounded and inf-sup stable; see below.

We further use and as scalar products on and , and thus obtain


The norm of is chosen such that and (A3) yields


for some . All estimates hold uniformly with respect to their arguments.

The following theorem summarizes the basic results about well-posedness and regularity of the solution for the above system with appropriate initial conditions.

Theorem 2.

Let Assumption 1 hold. Then for any , , and , there exists a unique weak solution

which satisfies (3)–(4) for a.e. and the initial conditions , where is defined by the compatibility conditions


If, in addition, and , then

In both cases, the solution can be bounded by the data in the natural norms.

The remainder of this section is devoted to a proof of these assertions.

2.2. Well-posedness

As usual, we associate via , , and to any of the bilinear forms a linear operator , , and its adjoint . This allows us to rewrite (3)–(4) as equivalent operator equations


which have to be understood as equations in the sense of and , respectively. From the properties of the bilinear forms, we immediately conclude the following result.

Lemma 3.

Let Assumption 1 hold. Then

  • and are symmetric, bounded, and boundedly invertible
    with , and , .

  • is bounded, surjective with closed range, and .

  • is bounded and injective with .

Using property (i) of the previous lemma, we can eliminate from (9) leading to


Inserting this into equation (10) yields the Schur complement problem


From the properties of the operators in Lemma 3, we deduce the following assertions.

Lemma 4.

Let Assumption 1 and the conditions on the data in Theorem 2 hold. Then

  • The operator is symmetric, bounded, and elliptic.

  • .

Note that, according to property (v), the operator induces a symmetric, bounded, and elliptic bilinear form , and the operator equation (12) can hence be written in an equivalent weak form as


This is an abstract parabolic equation whose well-posedness can be proven via Galerkin approximation; see [9, Chapter 7] or [7, Chapter XVIII]. We thus obtain

Lemma 5.

Let Assumption 1 hold. Then for any and , the reduced evolution problem (13) has a unique weak solution


with initial value . Moreover, by embedding.

This proves existence and uniqueness of a solution for problem (13) as well as a-priori estimates in the corresponding norms. Inserting into (11) and using the mapping properties of the operators and yields . The compatibility condition (8) then follows from continuity. Furthermore from (10), one can see that . By equivalence of (9)–(10) with the variational form, we obtain existence and uniqueness of a weak solution for the system (3)–(4) with the given initial conditions. Moreover, we have established the a-priori bounds in the first part of Theorem 2.

2.3. Regularity

We now turn to the improved a-priori estimates stated in Theorem 2. Additional regularity of the solution to the reduced problem (13) can be obtained with similar arguments as in [9, Chapter 7]. Assume that and with and . Then from the properties of the operators and , we deduce that with and . By formally testing (13) with , we obtain

and by integrating with respect to time, we further get

Here we used and , which follow directly from the properties of the bilinear forms and norms. The last two terms in the above estimate can be absorbed into the left hand side, which yields the required estimates for . Inserting this into (11) leads to the improved regularity for . This concludes the proof of Theorem 2.

2.4. Choice of initial conditions

Before we proceed, let us discuss in more detail the choice of initial conditions. Under Assumptions 1, the variational problem


is uniquely solvable for all and , which is a direct consequence of Brezzi’s splitting lemma [6]. In the above arguments, we simply chose and specified . By assumption (A1), we can determine from (15), and inserting into the equation (16) determines . Alternatively, one could set , choose freely, and then determine , by solving the coupled system (15)–(16). This again provides initial values and satisfying (8), where can be used as initial value for (12). The following choices of initial conditions are therefore equivalently possible:

  • choose and determine by (15);

  • choose and determine , by (15)–(16).

While the first choice is the natural one for the reduced problem (13), the second choice seems more natural for the coupled system (3)–(4). As indicated above, both choices are possible and actually equivalent in the considered functional analytic setting.

2.5. Energy-dissipation

The following property of the dynamical system will serve as the basic tool for the stability analysis of approximation schemes in later sections.

Lemma 6.

Let Assumption 1 hold and let denote a regular solution of (9)–(10) in the sense of Theorem 2. Then


Formal differentiation in time yields

which after rearrangement of the terms already yields the assertion of the lemma. ∎

Remark 7.

By integration in time, the validity of the stability estimate can be extended to less regular weak solutions, and these estimates provide an alternative route for proving uniqueness and a-priori estimates for weak solutions.

3. Galerkin approximation

For discretization of the variational problem (3)–(4), we now consider Galerkin approximations in finite dimensional sub-spaces and , i.e., we search for semi-discrete functions and satisfying


together with appropriate initial conditions to be specified below. Using similar arguments as for the analysis on the continuous level, we will show the following result.

Theorem 8.

Let Assumption 1 hold. Then (18)–(19) is a regular system of differential-algebraic equations of index . If is inf-sup stable on , i.e.,

  • with some constant ,

then a unique solution to the system (18)–(19) exists for any choice of initial values and satisfying


If (A3h) is not valid, additional compatibility conditions for and are required.

Remark 9.

A quick comparison with Theorem 2 shows that (20) are the natural compatibility conditions for the initial values of the problem under consideration, while the additional conditions required when (A3h) is not valid are artificial; see below. Inf-sup stable approximation spaces , are therefore required to guarantee the well-posedness of the semi-discrete problem without artificial conditions on the initial values.

In the following, we give a detailed proof of the above theorem and then turn to the error analysis of the Galerkin approximations defined in the beginning of the section.

3.1. Proof of Theorem 8

Choice of a basis for the spaces , allows to convert the semi-discrete problem into an equivalent system of differential-algebraic equations


with and

denoting the coordinate vectors of the discrete solutions and data, and

, , being matrices of appropriate size. From (A1)–(A2), one can see that the matrices and are symmetric and positive definite, and therefore, the matrix

is positive definite, e.g., for . As a consequence, the matrix pencil is regular and so the system (21) is a regular differential-algebraic equation; see [5, 14] for details. Differentiating the first equation in (21) leads to


Note that the matrix in front of the time derivatives is regular if, and only if, is surjective. In that case, (22

) is an (implicit) ordinary differential equation, and the existence of a unique solution

follows for any choice of initial conditions and . By integration of the first equation, one obtains

with . Hence solves the original system (21) if, and only if, , i.e., when the compatibility condition


is satisfied. If, on the other hand, is not surjective, then (22) is still a differential-algebraic equation. By change of basis, we may transform the system into

with being surjective and , both being positive definite. Differentiation of the second equation then leads to the system


Using the fact that is surjective and and are positive definite, one can deduce that the matrix in front of the time derivatives is regular, hence (24) is an ordinary differential equation. An inspection of the right hand side shows that none of the equations has been differentiated more than once, and hence the index of (21) is also one in this case. To obtain equivalence of (24) with (21), not only the compatibility condition (23) is required, but the additional artificial condition


has to be enforced, which is only caused by the inappropriate numerical approximation.

The assertions of Theorem 8 follow immediately from the above results by equivalence of the differential-algebraic system (21) with the weak formulation (18)–(19).

3.2. Abstract error analysis

We next turn to the a-priori analysis of Galerkin discretizations introduced in the beginning of this section. As usual, see e.g. [22, 23, 26], we decompose the error between the continuous and the semi-discrete solution via

into approximation errors and discrete error components. Following [15, 16], we utilize the variational problem corresponding to the stationary system associated with (3)–(4) to define the approximation and , i.e.,


Note that and can be computed by solving elliptic variational problems and error estimates for the elliptic projection can therefore be obtained by standard arguments; see Section 5 below. From the definition of the projections and the discrete solution, we immediately obtain the following discrete error equation.

Lemma 10.

Let Assumption 1 hold and denote a solution of (18)–(19) with initial values and . Then the discrete error components and satisfy , , and


With similar arguments as on the continuous level, we further obtain the following discrete energy–dissipation estimate for the discrete error.

Lemma 11.

Let Assumption 1 hold and be a solution of (18)–(19). Then

If also (A3h) holds, then additionally


From the discrete error equations (28)–(29), one can deduce that

Using , the embedding estimate , the definition of the norm , and Young’s inequality allow us to estimate

Then integrating with respect to time and using yield

The first estimate then follows by noting that and . From the discrete error equation (28) and continuity of the bilinear form , we deduce that

The discrete inf-sup condition (A3h) then leads to the second estimate of the lemma. ∎

Remark 12.

Let us emphasize that the particular choice of the elliptic projection in the error decomposition allows us to bound the discrete error by the approximation error in the component alone; this leads to improved error estimates and allows the use of post-processing techniques to obtain approximations for the pressure in polynomial spaces of higher order; see [4, 16] for details.

Using the previous bounds for the discrete error components, we now obtain the following estimates for semi-discrete approximation (18)–(19).

Theorem 13.

Let Assumption 1 hold. Then


If also the discrete inf-sup stability condition (A3h) holds, then additionally


After choosing the approximation spaces and , this estimate allows to derive quantitative error bounds via corresponding estimates for the elliptic projection. Details for a particular discretization will be given in Section 5.

4. Time discretization

We now turn to the discretization of the semi-discrete variational problem (18)–(19) in time. The guiding principle will be to preserve the underlying differential-algebraic structure and energy–dissipation identity elaborated in Section 3 as good as possible.

4.1. Preliminaries

Let be a partition of the time interval with time steps and . We denote by the space of piecewise polynomial functions of with values in . Following the notation of [22], we utilize capital letters to denote piecewise polynomial functions of time in the following.

In the spirit of [1] and to highlight the preservation of the problem structure, we first give a pointwise definition of our method, which requires the following two projection operators in time: We denote by the -orthogonal projection and by , , the -conforming projection, which is defined by the relations


We will refer to (33) as the commuting-diagram property. By integration of this relation and use of (34), one can see that and consequently is continuous on ; thus is an -conforming projection in time. We now consider the following time-discretization of the semi-discrete problem.

Problem 14 (Fully discrete scheme).

Find and such that , , and