Finite element method with the total stress variable for Biot's consolidation model

In this work, semi-discrete and fully-discrete error estimates are derived for the Biot's consolidation model described using a three-field finite element formulation. The fields include displacements, total stress and pressure. The model is implemented using a backward Euler discretization in time for the fully-discrete scheme and validated for benchmark examples. Computational experiments presented verifies the convergence orders for the lowest order finite elements with discontinuous and continuous finite element appropriation for the total stress.

Authors

• 2 publications
• 1 publication
• 14 publications
12/24/2021

Error Estimates of a Fully Discrete Multiphysics Finite Element Method for a Nonlinear Poroelasticity Model

In this paper, we propose a multiphysics finite element method for a non...
10/13/2019

A Conservative Finite Element Method for the Incompressible Euler Equations with Variable Density

We construct a finite element discretization and time-stepping scheme fo...
02/08/2018

A mixed finite element for weakly-symmetric elasticity

We develop a finite element discretization for the weakly symmetric equa...
06/24/2021

A micropolar isotropic plasticity formulation for non associated flow rule and softening featuring multiple classical yield criteria. Part II – FE integration and application

A Finite Element procedure based on a full implicit backward Euler predi...
08/04/2017

A Clinical and Finite Elements Study of Stress Urinary Incontinence in Women Using Fluid-Structure Interactions

Stress Urinary Incontinence (SUI) or urine leakage from urethra occurs d...
06/04/2021

Finite element stress analysis of a combined stacker-reclaimer machine: A design audit report

Design audit or design verification is an important step in engineering ...
05/03/2020

A mixed finite element method with reduced symmetry for the standard model in linear viscoelasticity

We introduce and analyze a new mixed finite element method with reduced ...
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

Biot’s fundamental equations for the soil consolidation process that describes the gradual adaptation of the soil to a load variation have been established in [1]. The mechanism of consolidation phenomenon described using a linear isotropic model is identical with the process of squeezing water out of an elastic porous medium in many cases. This solid-fluid coupling was extended to general anisotropy case in [2]. Such poroelasticity models have a lot applications in many areas including geomechanics [3], medicine [4], biomechanics [5], reservoir engineering [6]. The Biot’s consolidation models have also been used to combine transvascular and interstitial fluid movement with the mechanics of soft tissue in [7], which can be applied to improve drug delivery in solid tumors. Existence, uniqueness, and regularity theory were developed in [8] for poroelasticity and quasi-static problem in thermoelasticity. In [9], experiments with finite element method for the model in [7] have demonstrated the effects of fluid flow on the spatio-temporal patterns of soft-tissue elastic strain under a variety of physiological condition, which emphasized simulations relevant to a quasistatic elasticity imaging for the characterization of fluid flow in solid tumors.

Biot’s consolidation model has been considered by many researchers using finite element methods. In [10], a variational principle and the finite element method for a model with applications to a nonhomogeneous, anisotropic soil were developed. The fully discretization with backward Euler time discrete finite element method has been carried out and the existence and uniqueness were proved in [11]. Moreover, the simplest Taylor-Hood finite elements were employed. The stability and error estimates for the semi-discrete and fully-discrete finite element approximations were derived in [12] and [13]. Decay functions and post processing technique also were employed to improve the pore pressure accuracy. With the displacement and pore pressure fields as unknowns, the short and long time behavior of spatially discrete finite element solutions have been studied in [14]. Asymptotic error estimates have been derived for both stable and unstable combinations of the finite element spaces. For the coupled problem, a least squares mixed finite element method was presented and analyzed in [15]

with the unknowns fluid displacement, stress tensor, flux, and pressure. In

[16], coupling of mixed and continuous Galerkin finite element methods for pressure and displacements have been formulated deriving the convergence error estimates in time continuous setting. Methods for coupling mixed and discontinuous Galerkin have been presented in [17]. The error estimates for a fully-discrete stabilized discontinuous Galerkin method were obtained with the unknowns pressure and displacement in [18]. In [19], a discretization method in irregular domains with general grids for discontinuous full tensor permeabilities was developed. A new mixed finite element method for Biot’s consolidation problem in four variables was proposed in [20] and later, a three field mixed finite element which was free of pressure oscillations and Poisson locking has been proposed [21]. The priori error estimates that were robust for material parameters were provided in [22] with a four-field mixed method formulation. A three field finite element formulation with nonconforming finite element space for the displacements was considered in [23]. Based on the parameter dependent norms, the parameter-robust stability was established in [24], and parameter robust inf-sup stability and strong mass conservation were derived for three field mixed discontinuous Galerkin discretizations. A stabilized finite element method with equal order elements for the unknowns pressure and displacement was proposed in [25] to reduce the effects of non-physical oscillations. Combining the mixed method with symmetric interior penalty discontinuous Galerkin method obtained a H(div) conforming finite element method in [26]. The method achieved strong mass conservation.

Moreover, in [27], the total stress (or the soil pressure) expressed as a combination of the divergence of the velocity and pressure has been introduced coupling the solid and fluid robustly for Biot’s consolidation problem with the unknown displacement, pressure, and volumetric stress. Using a Fredholm argument for a static model, error estimates were derived independently of the Lam constants for both continuous and discrete formulations. A three field formulation of Biot’s model with the total stress variable has also been proposed in [28], which developed a parameter-robust block diagonal preconditioner for the associated discrete systems. Then, in [29], a priori error estimates for semi-discrete scheme has been presented by introducing the total pressure variable for quasi-static multiple-network poroelasticity equations. Our goal in this paper is to emphasize the time dependence of field variables in error estimates, so we choose to consider the fully-discrete scheme. Thus, using the total stress as a new variable for the three field formulation, we give the error estimates for semi-discrete and fully-discrete with backward Euler time discretization schemes.

Let be an open bounded polygonal or polyhedral domain in

. Biot’s consolidation model is described as follows: Find the displacement vector

and the fluid pressure in

 (1) −∇⋅(2με(u)+λ∇⋅uI−pI) =f,          in Ω×(0,¯T], (2) ∇⋅(Dtu)−∇⋅(k∇p) =g,          in Ω×(0,¯T], (3) u=0,  p =0,          on ΓD×(0,¯T], (4) (2με(u)+λ∇⋅uI−pI)⋅n =β,          on ΓN×(0,¯T], (5) (κ∇p)⋅n =γ,          on ΓN×(0,¯T],

where and are disjoint closed subsets with and the Dirichlet boundary . Here, is the strain tensor expressed in terms of symmetrized gradient of displacements, is the permeability of the porous solid and parameters , are the elastic Lam constants. The right hand side term in (1) represents the density of the applied body forces, and the source term in (2) represents a forced fluid extraction or injection. The outward unit normal vector on is denoted by .

Next, we introduce the coupling between the solid and fluid using , where is the total stress, is the displacement and is the fluid pressure [27, 28]. This can be rewritten as

 (6) λ−1(q−p)+∇⋅u=0.

For , the initial conditions are given by

 (7) u(x,0) =φ, p(x,0) =ϕ.

Hence, the initial condition for the total stress is . Denote the spaces and .

Multiplying (1), (6) and (2) by test functions, integrating by parts and applying boundary conditions yield the following weak formulation: Find , , such that

 (8) 2μ(ε(u),ε(v))−(q,∇⋅v) =(f,v)+⟨β,v⟩ΓN,      ∀v∈[H10,D(Ω)]d, (λ−1(q−p),wq)+(∇⋅u,wq) =0,                         ∀wq∈L2(Ω), −(λ−1(qt−pt),wp)+κ(∇p,∇wp) =(g,wp)+⟨γ,wp⟩ΓN,  ∀wp∈H10,D(Ω).

In Section 2, we present the semi-discrete and fully-discrete finite element formulation for system (8) with the unknown displacement, total stress and pressure and prove uniqueness of the solution for each of these schemes. Section 3 presents the derivation and analysis of the error estimates for both the semi-discrete and fully-discrete schemes. In Section 4 we present computational experiments on benchmark problems that validate the theoretical convergence rates with respect to mesh size and time step .

In the paper, we denote the arbitrary constants by , where is positive integer and is a constant which is independent of time step and mesh size . Let be the space of polynomials of degree less than or equal to in all variables. Moreover, let be the norm in space and be the norm in space. Denote the space-time space by for a Banach space (see details in [34]).

2 Finite Element Discretization and Uniqueness

Let be a regular and quasi-uniform triangulation of domain into triangular or tetrahedron elements [30, 31]. For each element , is its diameter and is the mesh size of triangulation . We consider the following finite element spaces on ,

 Uh:= {v∈[H10,D(Ω)]d∩[C0(Ω)]d:v∣T∈[Pk(T)]d,∀T∈Th}, Zh:= {q∈L2(Ω):q∣T∈Pl(T),∀T∈Th}, Ph:= {p∈H10,D(Ω)∩C0(Ω):p∣T∈Pk−1(T),∀T∈Th},

where In order to describe the initial conditions of discretization schemes, we define the following two projection operators. Let us define the Stokes projection and by

 (9) 2μ(ε(Quhu),ε(v))−(Qqhq,∇⋅v) =2μ(ε(u),ε(v))−(q,∇⋅v),      ∀v∈Uh, (wq,∇⋅Quhu) =(wq,∇⋅u),                          ∀wq∈Zh.

Also the elliptic projection is defined with the following properties,

 (10) (∇Qphp,∇wp) =(∇p,∇wp),  ∀wp∈Ph.

Hence, given a suitable approximation of initial conditions , and , the semi-discrete scheme corresponding to a three field formulation (8) is for all , to seek , , such that

 (11) 2μ(ε(uh),ε(v))−(qh,∇⋅v) =(f,v)+⟨β,v⟩ΓN,     ∀v∈Uh, (λ−1(qh−ph),wq)+(∇⋅uh,wq) =0,                         ∀wq∈Zh, −(λ−1(qht−pht),wp)+κ(∇ph,∇wp) =(g,wp)+⟨γ,wp⟩ΓN,  ∀wp∈Ph.

To obtain a fully-discrete formulation, we denote time step by , and , where is non-negative integer. Thus, given a suitable approximation of initial conditions , and , the fully-discrete scheme with backward Euler time discretization corresponding to the three field formulation (8) is to find , , such that

 (12) 2μ(ε(un),ε(v))−(qn,∇⋅v) =(fn,v)+⟨βn,v⟩ΓN,       ∀v∈Uh, (λ−1(qn−pn),wq)+(∇⋅un,wq) =0,                             ∀wq∈Zh, −(λ−1(¯∂tqn−¯∂tpn),wp)+κ(∇pn,∇wp) =(gn,wp)+⟨γn,wp⟩ΓN,   ∀wp∈Ph,

where and .

Then, we present an inequality, which will be useful in the uniqueness and convergence analysis. (Korn’s inequality) [32, 33] For each , there exists a positive constant such that

 (13) |u|1≤C(∥ε(u)∥+∥u∥).

Next we use Lemma 2 to prove the uniqueness for semi-discrete (11) and fully-discrete (12) schemes. For each , the semi-discrete scheme (11) has a unique solution. We need to prove that the homogeneous problem of (11) has only the trivial solution. Taking the time derivative of the second equation in (11), the homogeneous problem is rewritten as seeking , , with , , such that

 (14) 2μ(ε(uh),ε(v))−(qh,∇⋅v) =0,     ∀v∈Uh, (λ−1(qht−pht),wq)+(∇⋅uht,wq) =0,    ∀wq∈Zh, −(λ−1(qht−pht),wp)+κ(∇ph,∇wp) =0,    ∀wp∈Ph.

Using the test functions , and in (14) and simplifying, we can derive the following identity

 μddt∥ε(uh)∥2+λ−12ddt∥qh−ph∥2+κ∥∇ph∥2=0.

Integrating the above identity over , one finds

 μ∥ε(uh(t))∥2+λ−12∥qh(t)−ph(t)∥2+κ∫t0∥∇ph∥2ds=0.

Therefore, with the conditions , , , we have

 ∥ε(uh(t))∥=0, ∥qh(t)−ph(t)∥=0 and ∥∇ph∥=0.

Then, using Korn’s inequality from Lemma 2 when leads to . For , the fully-discrete scheme (12) has a unique solution. Similar to the semi-discrete, with , and , we rewrite the second equation of (12) and consider the homogeneous problem for fully-discrete scheme (12) is to find , , such that

 (15) 2μ(ε(un),ε(v))−(qn,∇⋅v) =0,       ∀v∈Uh, (λ−1(¯∂tqn−¯∂tpn),wq)+(∇⋅¯∂tun,wq) =0,      ∀wq∈Zh, −(λ−1(¯∂tqn−¯∂tpn),wp)+κ(∇pn,∇wp) =0,      ∀wp∈Ph.

Choosing , and , equation (15) can be simplified to

 μ∥ε(un)∥2−μ∥ε(un−1)∥2+λ−12∥qn−pn∥2−λ−12∥qn−1−pn−1∥2+κτ∥∇pn∥2≤0.

Here, we have used the inequality

 (ε(un),ε(un)−ε(un−1))≥12(∥ε(un)∥2−∥ε(un−1)∥2).

Summing over from to , it follows that

 μ∥ε(uN)∥2+λ−12∥qN−pN∥2+κτN∑n=1∥∇pn∥2≤0.

Note that the assumptions , , , and using Korn’s inequality (13) from Lemma 2, we have , and .

3 Error Estimates

In order to derive the error estimates, we first show the inf-sup condition for space pair with . Denote . For each , there exists satisfying and . Thus, we can deduce that for [35, 27]

Since not all the discrete finite element spaces meet the inf-sup condition, we assume that the space pair satisfy the inf-sup condition, i.e. there exists a positive constant such that for

 (16) sup0≠vh∈Uhb(vh,wqh)|vh|1≥C∥wqh∥.

3.1 Error estimate for the semi-discrete scheme

For the semi-discrete scheme, we denote the error in displacement by , where and . Similarly, we decompose the errors of the total stress and pressure into two parts, respectively, i.e. and . Assume that the inf-sup condition (16) is satisfied for the space pair . Let and be the solutions of (8) and (11) respectively. Then there exists a constant such that for each

 (17) μ∥ε(euh(t))∥2+∥eqh(t)∥2≤C(μ∥ε(ηuh(t))∥2+∥ηqh(t)∥2+∫t0∥ηqht−ηpht∥2ds),

and

 (18) κ∥∇eph(t)∥2≤C(κ∥∇ηph(t)∥2+∫t0∥ηqht−ηpht∥2ds).

Subtracting (11) from (8), and taking the derivative of the second equation with respect to time , we get

 2μ(ε(euh),ε(v))−(eqh,∇⋅v) =0,     ∀v∈Uh, (λ−1(eqht−epht),wq)+(∇⋅euht,wq) =0,    ∀wq∈Zh, −(λ−1(eqht−epht),wp)+κ(∇eph,∇wp) =0,    ∀wp∈Ph.

With the use of the definitions of projections (9) and (10), it follows

 (19) 2μ(ε(ξuh),ε(v))−(ξqh,∇⋅v) =0, (λ−1(ξqht−ξpht),wq)+(∇⋅ξuht,wq) =−(λ−1(ηqht−ηpht),wq), −(λ−1(ξqht−ξpht),wp)+κ(∇ξph,∇wp) =(λ−1(ηqht−ηpht),wp).

Taking , and , then we can deduce that

 μddt∥ε(ξuh)∥2+λ−12ddt∥ξqh−ξph∥2+κ∥∇ξph∥2 =−λ−1(ηqht−ηpht,ξqh−ξph) ≤λ−12∥ηqht−ηpht∥2+λ−12∥ξqh−ξph∥2.

Next, integrating over , since and , we have

 μ∥ε(ξuh(t))∥2+λ−12∥ξqh(t)−ξph(t)∥2+κ∫t0∥∇ξph∥2ds ≤ λ−12∫t0∥ηqht−ηpht∥2ds+λ−12∫t0∥ξqh−ξph∥2ds.

Using the Gronwall Lemma [34], one finds

 (20) μ∥ε(ξuh(t))∥2+λ−12∥ξqh(t)−ξph(t)∥2+κ∫t0∥∇ξph∥2ds≤C∫t0∥ηqht−ηpht∥2ds.

Now, to estimate , we deduce from the first equation of (19), for each

 b(v,ξqh)=(ξqh,∇⋅v)=2μ(ε(ξuh),ε(v)).

Thus, using the inf-sup condition (16), it follows that

 (21) ∥ξqh∥≤Csup|v|1≠0(ξqh,∇⋅v)|v|1=Csup|v|1≠02μ(ε(ξuh),ε(v))|v|1≤C∥ε(ξuh)∥.

On the other hand, differentiating the first equation of (19) with respect to , it follows

 2μ(ε(ξuht),ε(v))−(ξqht,∇⋅v) =0,     ∀v∈Uh.

Taking in the above identity and , in (19), we have

 2μ∥ε(ξuht)∥2+λ−1∥ξqht−ξpht∥2+κ2ddt∥∇ξph∥2 =−λ−1(ηqht−ηpht,ξqht−ξpht) ≤λ−12∥ηqht−ηpht∥2+λ−12∥ξqht−ξpht∥2.

So we rewrite the above inequality as following

 2μ∥ε(ξuht)∥2+λ−12∥ξqht−ξpht∥2+κ2ddt∥∇ξph∥2≤λ−12∥ηqht−ηpht∥2.

Integrating on , we have

 (22) 4μ∫t0∥ε(ξuht)∥2ds+λ−1∫t0∥ξqht−ξpht∥2ds+κ∥∇ξph(t)∥2≤C∫t0∥ηqht−ηpht∥2ds.

The proof is completed combining (20), (21), (22) with the definitions of errors. Moreover, as an immediate application of (20), (21), and (22), we have the following corollary. Assume that the inf-sup condition (16) is satisfied for the space pair . Let and be the solutions of (8) and (11), then there exists a constant such that for each

 ∫t0(μ∥ε(euht)∥2+∥eqht∥2)ds≤C∫t0(μ∥ε(ηuht)∥2+∥ηqht∥2+∥ηqht−ηpht∥2)ds,

and

 κ∫t0∥∇eph∥2≤C∫t0(κ∥∇ηph∥2+∫t0∥ηqht−ηpht∥2)ds.

3.2 Error estimate for the fully-discrete scheme

Denote the error in displacement at time by . We define and in a similar fashion as . Assume that the inf-sup condition (16) is satisfied for the space pair . Let and be the solutions of (8) and (12), then there exists a constant such that for

 (23) μ∥ε(eNu)∥2+∥eNq∥2≤ C(μ∥ε(ηNu)∥2+∥ηNq∥2+τ∫tN0(∥ηqht∥2+∥ηpht∥2)ds +τ3∫tN0(∥qtt∥2+∥ptt∥2)ds+τ2∫tN0∥utt∥21ds),

and

 (24) κ∥∇eNp∥2 ≤κ∥∇ηNp∥2+C(τ2∫tN0(∥qtt