Iterative solvers for Biot model under small and large deformation

TL;DR

This paper develops and analyzes iterative solvers, including Newton and L-scheme methods, for the Biot model under both small and large deformations, incorporating nonlinear fluid compressibility and Saint Venant-Kirchoff law.

Contribution

It introduces a unified framework for linearization schemes applicable to any spatial discretization of the Biot model with nonlinear fluid and large deformation considerations.

Findings

Convergence of schemes proven analytically for small deformation cases.

Numerical examples demonstrate effectiveness for large deformation scenarios.

Schemes are flexible for different discretization methods.

Abstract

We consider L-scheme and Newton based solvers for Biot model under small or large deformation. The mechanical deformation follows the Saint Venant-Kirchoff constitutive law. Further, the fluid compressibility is assumed to be nonlinear. A Lagrangian frame of reference is used to keep track of the deformation. We perform an implicit discretization in time (backward Euler) and propose two linearization schemes for solving the nonlinear problems appearing within each time step: Newton's method and L-scheme. The linearizations are used monolithically or in combination with a splitting algorithm. The resulting schemes can be applied for any spatial discretization. The convergences of all schemes are shown analytically for cases under small deformation. Illustrative numerical examples are presented to confirm the applicability of the schemes, in particular, for large deformation.

Tables

Table 1. Table 1 : The coefficient functions 𝔟 ​ ( ⋅ ) , 𝔠 ​ ( ⋅ ) 𝔟 ⋅ 𝔠 ⋅ \mathfrak{b}(\cdot),\mathfrak{c}(\cdot) for test problem 1.

Case	$𝔟\left(p\right)$	$𝔠\left(\nabla \cdot 𝐮\right)$
1	$e^p$	${\left(\nabla \cdot 𝐮\right)}^3+\nabla \cdot 𝐮$
2	$e^p$	${\left(\nabla \cdot 𝐮\right)}^3$
3	$e^p$	$\sqrt[3]{\nabla \cdot {𝐮}^5}+\nabla \cdot 𝐮$
4	$p^2$	$\nabla \cdot {𝐮}^2$

Table 2. Table 2 : Boundary conditions for Traction and Rotation case respectively.

Face	Flow	Mechanics
Top	$p=0$	$𝐮=\left(𝐑\left(\theta \left(t\right)\right)-𝐈\right)X_0$
Bottom	$p=0$	$𝐮\cdot \overrightarrow{n}=0$
Lateral	$p=0$	$\overrightarrow{\Pi }\cdot \overrightarrow{n}=0$

Equations

$$-\nabla \cdot \mathbf{\Pi }={\rho }_b\mathbf{g},$$

$$\mathbf{\Sigma }={\mathbf{\Sigma }}^{eff}-J{\mathbf{F}}^{-1}{\mathbf{F}}^{\top }p,$$

$$\dot{\Gamma }+\nabla \cdot \mathbf{q}=S_f.$$

$$\mathbf{q}=-\mathbf{K}(\mathbf{u})\left(\nabla p-{\rho }_f{\mathbf{g}}_0\right),$$

$$-\nabla \cdot \mathbf{\Pi }\left(\nabla \mathbf{u},p\right)$$

$$\mathbf{q}$$

$$\dot{\Gamma }(\mathbf{u},p)+\nabla \cdot \mathbf{q}$$

$$\displaystyle\begin{array}[]{ll}\mathcal{F}_{\rm{mech}}(\mathbf{u},p)&=-\nabla\cdot{\bf\Pi}\left(\nabla\mathbf{u},p\right)-\rho_{b}\mathbf{g},\\[5.59721pt] \mathcal{F}_{\rm{darcy}}(\mathbf{u},p)&=\mathbf{q}+{\mathbf{K}}(\mathbf{u})\left(\nabla p-\rho_{f}\Upsilon\right),\\[5.59721pt] \mathcal{F}_{\rm{mass}}(\mathbf{u},p)&=\dot{\Gamma}(\mathbf{u},p)+\nabla\cdot\mathbf{q}-S_{f}.\end{array}$$

$$\displaystyle\begin{array}[]{ll}-\nabla\cdot\left(\mathcal{\partial}_{\mathbf{u}}{\bf\Pi}\left(\nabla\mathbf{u}^{i-1},p^{i-1}\right)\nabla\delta\mathbf{u}^{i}-\mathcal{\partial}_{p}{\bf\Pi}\left(\nabla\mathbf{u}^{i-1},p^{i-1}\right)\delta p^{i}\right)&=-\mathcal{F}_{mech}(\mathbf{u}^{i-1},p^{i-1}),\\[5.59721pt] \delta\mathbf{q}^{i}+{\mathbf{K}}(\mathbf{u}^{i-1})\nabla\delta p^{i-1}+\mathcal{\partial}_{\mathbf{u}}{\mathbf{K}}(\mathbf{u}^{i-1})\nabla p^{i-1}\delta\mathbf{u}^{i}&=-\mathcal{F}_{darcy}(\mathbf{u}^{i-1},p^{i-1}),\\[5.59721pt] \mathcal{\partial}_{p}\dot{\Gamma}(\mathbf{u}^{i-1},p^{i-1})\delta p^{i}+\mathcal{\partial}_{\mathbf{u}}\dot{\Gamma}(\mathbf{u}^{i-1},p^{i-1})\delta{\mathbf{u}}^{i}+\nabla\cdot\delta\mathbf{q}^{i}&=-\mathcal{F}_{mass}(\mathbf{u}^{i-1},p^{i-1}),\end{array}$$

$$\left({\mathbf{u}}^i,{\mathbf{q}}^i,p^i\right)=\left({\mathbf{u}}^{i-1},{\mathbf{q}}^{i-1},p^{i-1}\right)+\left(\delta {\mathbf{u}}^i,\delta {\mathbf{q}}^i,\delta p^i\right).$$

$$\displaystyle\begin{array}[]{ll}\delta\mathbf{q}^{i}+{\mathbf{K}}(\mathbf{u}^{i-1})\nabla\delta p^{i-1}&=-\mathcal{F}_{darcy}(\mathbf{u}^{i-1},p^{i-1}),\\[5.59721pt] \mathcal{\partial}_{p}\dot{\Gamma}(\mathbf{u}^{i-1},p^{i-1})\delta p^{i}+\nabla\cdot\delta\mathbf{q}^{i}&=-\mathcal{F}_{mass}(\mathbf{u}^{i-1},p^{i-1}),\end{array}$$

$$\left({\mathbf{q}}^i,p^i\right)=\left({\mathbf{q}}^{i-1},p^{i-1}\right)+\left(\delta {\mathbf{q}}^i,\delta p^i\right).$$

$$-\nabla \cdot \left({\partial }_{\mathbf{u}}\mathbf{\Pi }\left(\nabla {\mathbf{u}}^{i-1},p^i\right)\nabla \delta {\mathbf{u}}^i-L_s(\nabla \cdot \delta {\mathbf{u}}^i)\mathbf{I}\right)$$

$${\mathbf{u}}^i={\mathbf{u}}^{i-1}+\delta {\mathbf{u}}^i.$$

$$\displaystyle\begin{array}[]{ll}-\nabla\cdot{\bf{L}}_{\mathbf{u}}\nabla\delta\mathbf{u}^{i}-\nabla\cdot{\bf{L}}_{p}\delta p^{i}&=-\mathcal{F}_{mech}(\mathbf{u}^{i-1},p^{i-1}),\\[5.59721pt] \delta\mathbf{q}^{i}+{\mathbf{K}}(\mathbf{u}^{i-1})\nabla\delta p^{i}+{\bf{L}}_{\mathbf{q}}\delta\mathbf{u}^{i}&=-\mathcal{F}_{darcy}(\mathbf{u}^{i-1},p^{i-1}),\\[5.59721pt] L_{p}\delta p^{i}+L_{\mathbf{u}}\delta{\mathbf{u}}^{i}+\nabla\cdot\delta\mathbf{q}^{i}&=-\mathcal{F}_{mass}(\mathbf{u}^{i-1},p^{i-1}),\end{array}$$

$$\left({\mathbf{u}}^i,{\mathbf{q}}^i,p^i\right)=\left({\mathbf{u}}^{i-1},{\mathbf{q}}^{i-1},p^{i-1}\right)+\left(\delta {\mathbf{u}}^i,\delta {\mathbf{q}}^i,\delta p^i\right).$$

$$\displaystyle\begin{array}[]{ll}\delta\mathbf{q}^{i}+{\mathbf{K}}(\mathbf{u}^{i-1})\nabla\delta p^{i}&=-\mathcal{F}_{darcy}(\mathbf{u}^{i-1},p^{i-1}),\\[5.59721pt] L_{p}\delta p^{i}+\nabla\cdot\delta\mathbf{q}^{i}&=-\mathcal{F}_{mass}(\mathbf{u}^{i-1},p^{i-1}),\end{array}$$

$$\left({\mathbf{q}}^i,p^i\right)=\left({\mathbf{q}}^{i-1},p^{i-1}\right)+\left(\delta {\mathbf{q}}^i,\delta p^i\right).$$

$$-\nabla \cdot \left({\mathbf{L}}_{\mathbf{u}}\nabla \delta {\mathbf{u}}^i+L_s(\nabla \delta \cdot {\mathbf{u}}^i)\mathbf{I}\right)$$

$${\mathbf{u}}^i={\mathbf{u}}^{i-1}+\delta {\mathbf{u}}^i.$$

$$\mathbf{\Pi }(\mathbf{u},p)=\sigma (\mathbf{u},p)=2\mu \epsilon (\nabla \mathbf{u})+\mathfrak{c}(\text{tr}(\epsilon (\nabla \mathbf{u})))-\alpha p\mathbf{I},$$

$$\Gamma (\mathbf{u},p)$$

$$\displaystyle\begin{array}[]{rl}\left(\varepsilon({\mathbf{u}}),\varepsilon({\mathbf{v}})\right)+\left(\mathfrak{c}(\nabla\cdot{\mathbf{u}})-\alpha p,\nabla\cdot{\mathbf{v}}\right)&=(\rho_{b}\mathbf{g},{\mathbf{v}}),\ \ \forall{\mathbf{v}}\in\left(H(\Omega)\right)^{d},\\[5.59721pt] \left({\mathbf{K}}^{-1}\mathbf{q},{\bf z}\right)-\left(p,\nabla\cdot\bf z\right)&=\left(\rho_{f}\mathbf{g},{\bf z}\right),\ \ \forall{\bf z}\in H^{1}(\text{div},\Omega),\\[5.59721pt] \left(\dot{\mathfrak{b}}(p))+\alpha\nabla\cdot\dot{{\mathbf{u}}},w\right)+\tau\left(\nabla\cdot\mathbf{q},w\right)&=\tau\left(S_{f},w\right),\ \forall w\in L^{2}(\Omega),\end{array}$$

$$\left(\mathfrak{b}(p_0))+\alpha \nabla \cdot {\mathbf{u}}_0,w\right)=0,\forall w\in L^2(\Omega ).$$

$${\mathbf{V}}_h$$

$$W_h$$

$${\mathbf{Z}}_h$$

$$\displaystyle\begin{array}[]{rl}\left(\varepsilon({\mathbf{u}}_{h}^{n}),\varepsilon({\mathbf{v}}_{h})\right)+\left(\mathfrak{c}(\nabla\cdot{\mathbf{u}}_{h}^{n}),\nabla\cdot{\mathbf{v}}_{h}\right)-\alpha\left(p_{h}^{n},\nabla\cdot{\mathbf{v}}_{h}\right)&=(\rho_{b}\mathbf{g},{\mathbf{v}}_{h}),\\[5.59721pt] \left({\mathbf{K}}^{-1}\mathbf{q}_{h}^{n},{\bf z}_{h}\right)-\left(p_{h}^{n},\nabla\cdot\bf z_{h}\right)&=\left(\rho_{f}\mathbf{g},{\bf z}_{h}\right),\\[5.59721pt] \left(\mathfrak{b}(p_{h}^{n})-\mathfrak{b}(p_{h}^{n-1}),w_{h}\right)+\alpha\left(\nabla\cdot({\mathbf{u}}_{h}^{n}-{\mathbf{u}}_{h}^{n-1}),w_{h}\right)\\[5.59721pt] +\tau\left(\nabla\cdot\mathbf{q}_{h}^{n},\nabla w_{h}\right)&=\tau(S_{f},w_{h}),\end{array}$$

$${\mathbf{e}}_{\mathbf{u}}^{n,i}$$

$${\mathbf{e}}_{\mathbf{q}}^{n,i}$$

$$e_p^{n,i}$$

$$\displaystyle\begin{array}[]{rl}\left(\varepsilon({\mathbf{u}}_{h}^{n,i}),\varepsilon({\mathbf{v}}_{h})\right)+\left(\mathfrak{c}(\nabla\cdot{\mathbf{u}}_{h}^{n,i-1})+\mathfrak{c}^{\prime}(\nabla\cdot{\mathbf{u}}_{h}^{n,i-1})\nabla\cdot\delta{\mathbf{u}}_{h}^{n,i},\nabla\cdot{\mathbf{v}}_{h}\right)\\[5.59721pt] -\left(\alpha p_{h}^{n,i},\nabla\cdot{\mathbf{v}}_{h}\right)=(\rho_{b}\mathbf{g},{\mathbf{v}}_{h}),\\[5.59721pt] \left({\mathbf{K}}^{-1}\mathbf{q}_{h}^{n,i},{\bf z}_{h}\right)-\left(p_{h}^{n,i},\nabla\cdot\bf z_{h}\right)=\left(\rho_{f}\mathbf{g},{\bf z}_{h}\right),\\[5.59721pt] \left(\mathfrak{b}(p_{h}^{n,i-1})+\mathfrak{b}^{\prime}(p_{h}^{n,i-1})\delta p_{h}^{n,i}-\mathfrak{b}(p_{h}^{n-1}),w_{h}\right)+\left(\alpha\nabla\cdot({\mathbf{u}}_{h}^{n,i}-{\mathbf{u}}_{h}^{n-1}),w_{h}\right)\\[5.59721pt] +\tau\left(\nabla\cdot\mathbf{q}_{h}^{n,i},\nabla w_{h}\right)=\tau(S_{f},w_{h}),\end{array}$$

$$x_n\le a{x}_{n-1}^2+b{x}_{n-1}\forall n\ge 1,$$

Full text

1 Introduction

The coupling of flow and mechanics in a porous medium, typically referred to as poromechanics, plays a crucial role in many socially relevant applications. These include geothermal energy extraction, energy storage in the subsurface, $CO_{2}$ sequestration, and understanding of biological tissues. The increased role played by computing in the development and optimisation of (industrial) technologies for these applications implies the need for improved mathematical models in poromechanics and robust numerical solvers for them.

The most common mathematical model for coupled flow and mechanics in porous media is the linear, quasi-stationary Biot model [8, 9, 10, 52]. The model consists of two coupled partial differential equations, representing balance of forces for the mechanics and conservation of mass and momentum for (single-phase) flow in porous media.

In terms of modelling, Biot’s model has been extended to unsaturated flow [14, 37], multiphase flow [27, 28, 34, 36, 48], thermo-poro-elasticity [19], and reactive transport in porous media [33, 49], where nonlinearities arise in the flow model, specifically in the diffusion term, the time derivative term and/or in Biot’s coupling term. The mechanics model can also be extended to the elasto-plastic [3, 56], the fracture propagation [35] and the hyperelasticity [20, 21], where the nonlinearities appear in the constitutive law of the material, in the compatibility condition and/or the conservation of momentum equation. Furthermore, elastodynamics or non-stationary Biot, i.e. Biot-Allard model [38], includes a convolution in the coupling term of both mechanics and flow equations. In this paper, we are going to explore a general case that allows large deformations. The mechanical deformation follows the Saint Venant-Kirchoff constitutive law and the fluid compressibility in the fluid equation is assumed to be nonlinear. This model formulation is needed to later consider extensions of Biot’s model to plasticity, more general hyperelastic materials, and elastodynamics.

Finding closed-form solutions for coupled problems is very difficult, and commonly based on various simplifications. We, therefore, resort to numerical approximations. In general, there are two approaches to solve such problems, the fully coupled and the weakly coupled scheme. In general the fully coupled schemes for fluid potential and mechanical deformation are stable, have excellent convergence properties, and ensure that the numerical solution is consistent with the underlying continuous differential equations [29, 55]. Despite obvious advantages, the monolithic solver for the fully coupled problem are more difficult to implement, and have difficulties solving the resulting linear system, particularly in the context of existing legacy codes for separate physics. In the weakly coupled approach, while marching in time, we time-lag the flow problem (or the mechanics), thereby fully decoupling the two problems. Due to the complexities associated with the fully coupled scheme, the industry standard remains to use weakly coupled or iteratively coupled approaches [18, 42, 51, 59]. An iteratively coupled approach takes somewhat of a middle path; at each time step, it decouples the flow and mechanics, but iterates so that the convergence is achieved. Weakly coupled schemes, wherein there are not iterations within time step, have in particular been questioned in previous works [17, 22, 42, 45]; they have been shown to lack robustness and even convergence, if not properly designed. In order to ensure the robustness and accuracy of the resulting computations, it is therefore essential to understand the efficiency, stability, and convergence of iterative coupling schemes, in particular in the presence of nonlinearities.

In this work, we present monolithic and splitting approaches for solving this nonlinear system, that is, nonlinear compressibility and the Saint Venant-Kirchoff constitutive law for stress-strain. Moreover, we rigorously study the convergence of our schemes, including the Newton based ones, under the assumption of small deformations. As for splitting approach, we use the undrained split method, see [31, 39]. We use linear conformal Galerkin elements for the discretization of the mechanics equation and mixed finite elements for the flow equation [7, 23, 30, 43, 58]. Precisely, the lowest order Raviart-Thomas elements are used [16]. We expect, however, that the solution strategy discussed herein will be applicable to other combinations of spatial discretizations such as those discussed in [40, 50] and the references therein. Backward Euler is used for the temporal discretization.

To summarise, the new contributions of this paper are

• $\bullet$

We propose Newton and $L$-scheme based monolithic and splitting schemes for solving the Biot model under small or large deformation.

• $\bullet$

The convergence analysis of all schemes is shown rigorously under the assumption of small deformations.

• $\bullet$

We provide a benchmark for the convergence of splitting algorithms for a general nonlinear Biot model that includes large deformations.

We mention some relevant works in this direction. For the convergence analysis of the undrained split method applied to the linear Biot model, we refer to [5, 6, 12, 24, 25, 39]. For a discussion on the stabilization/tuning parameter used in the undrained split approach, we refer to [12, 15]. A theoretical investigation on the optimal choice for this parameter is performed in [53]. The linearization is based on either Newton’s method, or the $L$-scheme [37, 44, 48] or a combination of them [14, 37]. For monolithic and splitting schemes based solely on $L$-scheme, we refer to [11]. Multirate time discretizations or higher order space-time Galerkin method has also been proposed for the linear Biot model in [1] and [6], respectively.

The paper is structured as follows. In the next section, we present the mathematical model. In Section 3, we propose four iterative schemes. Section 4 shows the analysis of iterative schemes under the assumption of small deformations. Numerical results are presented in Section 5 followed by the conclusion in Section 6.

2 Governing equations

We consider a fluid flow problem in a poroelastic bounded reference domain $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$ under large deformation. A Lagrangian frame of reference is used to keep track of the invertible transformation $x:=\{x(X,t)=X+\mathbf{u}(X,t):X\in\Omega\rightarrow x\in\Omega_{t}\}$, where $\Omega_{t}$ is the deformed domain at time $t$ and $\mathbf{u}$ represents the deformation field. The gradient of the transformation and its determinant are given by ${\mathbf{F}}=\nabla\,x(X,t)$ and $J=\det({\mathbf{F}})$. All differentials are with respect to the undeformed coordinates $X$, unless otherwise stated.

We will now write the conservation of momentum and mass equation in $\Omega$. The conservation of momentum represents the balance between the first Piola-Kirchhoff poroelastic stress ${\bf\Pi}$ in $\Omega$ and the forces acting on $\Omega_{t}$, and is given by

[TABLE]

where $\rho_{b}=J\varrho_{b}$ is the bulk density in $\Omega$, $\varrho_{b}$ is the bulk density in $\Omega_{t}$ and $\mathbf{g}$ is gravity.

We exploit the relation ${\bf\Pi}={\mathbf{F}}{\bf\Sigma}$ since the constitutive laws are developed for the second Piola-Kirchhoff poroelastic stress $\bf\Sigma$. This stress tensor is composed of the effective mechanical stress ${\bf\Sigma}^{eff}$ and the pore pressure $p$ by the following relation

[TABLE]

where $J{\mathbf{F}}^{-1}{\mathbf{F}}^{\top}$ ensures that pressure $p$ exerts an isotropic stress in $\Omega_{t}$. We assume an isotropic poroelastic material with constant shear modulus $\mu$ and a nonlinear function of the volumetric strain $\mathfrak{c}(\cdot)$ [11, 54]. The effective stress is given by Saint Venant-Kirchhoff constitutive law: $\bf\Sigma^{eff}=2\mu{\bf E}+\mathfrak{c}\left(\text{tr}({\bf E})\right),$ where the Green strain tensor ${\bf E}$ is defined by

${\bf E}=\frac{1}{2}\left(\nabla\mathbf{u}+\nabla^{\top}\mathbf{u}+(\nabla\mathbf{u})^{\top}\nabla\mathbf{u}\right)$.

The conservation of fluid mass is given by

[TABLE]

We consider a fluid mass ${\Gamma}=J\rho_{f}\phi$ of a slightly compressible fluid, where $\phi$ is the porosity and $\rho_{f}$ the fluid density and $S_{f}$ the source term in $\Omega$ respectively. The time derivative of the fluid content $\dot{\Gamma}=\dot{\Gamma}(\mathbf{u},p)$ is considered to be a function of the pressure and the pore volume change due to the deformation field. We consider Darcy’s law

[TABLE]

where the flux variable $\mathbf{q}$ is the first Piola transform of the corresponding flux variable in $\Omega_{t}$, ${\mathbf{K}}=J{\mathbf{F}}^{-1}{\bf{k}}{\mathbf{F}}^{-\top}$ is the corresponding transformation of the mobility tensor ${\bf{k}}$ in $\Omega_{t}$ and $\Upsilon={\mathbf{F}}^{\top}\mathbf{g}$. Finally, the general nonlinear Biot model considered in this paper reads as:

Find $\left(\mathbf{u},\mathbf{q},p\right)$ such that

[TABLE]

To complete the model we consider Dirichlet boundary conditions (BC) and initial conditions given by $(\mathbf{u}_{0},p_{0})$ such that $\Gamma(\mathbf{u}_{0},p_{0})=\Gamma_{0}$ and ${\bf\Pi}(\mathbf{u}_{0},p_{0})={\bf\Pi}_{0}$ at time $t=0$. The functions $\Gamma_{0}$ and ${\bf\Pi}_{0}$ are supposed to be given (and to be sufficiently regular).

In practice, the initial data $\mathbf{u}_{0}$ and $p_{0}$ are not independent and can be obtained by solving the flow equation for $p_{0}$ and then solving the mechanics equation for getting $\mathbf{u}_{0}$.

3 Iterative schemes

In this section, we present several monolithic and splitting iterative schemes for solving Eqs. (4). First, we propose the Newton method which is well known for having quadratic convergence. Secondly, we combine the Newton method with a stabilized splitting scheme based on the undrained split method. Finally, for the third and fourth schemes, we propose monolithic and splitting $L$-schemes. The iterative schemes will be written using an incremental formulation. In this regard, we introduce naturally defined residuals for the nonlinear Eqs. (4).

[TABLE]

We will denote by $\delta(\cdot)^{i}=(\cdot)^{i}-(\cdot)^{i-1}$ the incremental operator, $i$ the incremental counter, $\partial_{(\cdot)}$ the partial derivative operator respect to $(\cdot)$.

3.1 A monolithic Newton solver

The Newton method is usually the first choice of the linearization methods due to its quadratic convergence. However, the convergence is local and it requires relatively small time steps to ensure the quadratic convergence [47]. The method starts by using initial solution $\left(\mathbf{u}^{0},\mathbf{q}^{0},p^{0}\right)$, solves for $\left(\delta\mathbf{u}^{i},\delta\mathbf{q}^{i},\delta p^{i}\right)$ satisfying

[TABLE]

and finally updates the variables

[TABLE]

3.2 A splitting Newton solver

The splitting Newton method combines a splitting method with the Newton linearization. We introduce a stabilization parameter $L_{s}\geq 0$ to stabilize the mechanics equation. The precise condition on $L_{s}$ to ensure convergence is shown in Theorem 2. The method consists on two steps: starting with the initial condition $\left(\mathbf{u}^{0},\mathbf{q}^{0},p^{0}\right)$:

Step 1: solve for $\left(\delta\mathbf{q}^{i},\delta p^{i}\right)$

[TABLE]

and update the variables

[TABLE]

Step 2: solve for $\delta\mathbf{u}^{i}$ satisfying

[TABLE]

and update the variable

[TABLE]

The stability of the scheme is controlled by $L_{s}$ as it is shown in [47].

3.3 A monolithic $L$-scheme

The $L$-scheme can be interpreted as either a stabilized Picard method or a quasi-Newton method. This scheme is robust but only linearly convergent. Moreover, it can be applied to non-smooth but monotonically increasing nonlinearities. For example, for the case of Hölder continuous (not Lipschitz) nonlinearities we refer to [13]. As it is a fixed point scheme, it can be speeded up by using the Anderson acceleration [2, 15]. To summarize, the main advantages of the $L$-scheme are:

• •

It does not involve computation of derivatives.

• •

The arising linear systems are well-conditioned.

• •

It can be applied to non-smooth nonlinearities.

• •

It is easy to understand and implement.

A monolithic $L$-scheme requires three constant tensors ${\bf{L}}_{\mathbf{u}},\ {\bf{L}}_{p},\ {\bf{L}}_{\mathbf{q}}\in\mathbb{R}^{d\times d}$ and two positive constants $L_{p}$ and $L_{\mathbf{u}}$ as linearization parameters. A practical choice of the linearization parameters will be discussed in the numerical section. We refer to [11, 22] for a discussion regarding the best choice for the linearization parameters $L_{p}$ and $L_{\mathbf{u}}.$

The method starts with the given initial solution $\left(\mathbf{u}^{0},\mathbf{q}^{0},p^{0}\right)$ and solve for $\left(\delta\mathbf{u}^{i},\delta\mathbf{q}^{i},\delta p^{i}\right)$

[TABLE]

and then update the variables

[TABLE]

3.4 A splitting $L$-scheme

The splitting scheme requires less linearization terms: two constants $\mathbf{L}_{\mathbf{u}}\in\mathbb{R}^{d\times d}$, $L_{p}\geq 0$ and a positive stabilisation term $L_{s}$. This makes it suitable for quick implementation since there is no need to calculate any Jacobian. The method is split in two steps, given initial solution $\left(\mathbf{u}^{0},\mathbf{q}^{0},p^{0}\right)$:

Step 1: solve for $\left(\delta\mathbf{q}^{i},\delta p^{i}\right)$

[TABLE]

update the variables

[TABLE]

Step 2: solve for $\delta\mathbf{u}^{i}$

[TABLE]

and then update the variables

[TABLE]

4 The Biot model under small deformations

The convergence analysis of the iterative schemes proposed cannot be addressed with standard techniques [11, 15, 14, 37, 39]. This is due to the nonlinearities being non-monotone. Nevertheless, a rigorous analysis can be performed for the case of small deformations. Accordingly, we assume the porous medium to be under small deformation and present the convergence of the iterative schemes proposed in the previous section.

Under small deformation, the different between $\Omega_{t}$ and $\Omega$ can be neglected. The gradient of the transformation is approximated by ${\mathbf{F}}\approx{\mathbf{I}}$ and the determinant of the transformation by $J\approx 1$. Additionally, the Green strain tensor $\bf E$ can be approximated by the infinitesimal strain tensor ${\bf E}\approx\varepsilon=\frac{1}{2}\left(\nabla\mathbf{u}+(\nabla\mathbf{u})^{\top}\right)$. Then, the poroelastic stress tensor can be expressed by

[TABLE]

where $\alpha$ is the Biot constant. The mobility tensor is considered isotropic ${\bf K}(\mathbf{u},p)=k{\mathbf{I}}$, but the results of the convergence analysis can be extended without difficulties to a more general anisotropic case. Additionally, the time derivative of the volumetric deformation is approximated by $\dot{J}\approx\nabla\cdot\dot{\mathbf{u}}$. In this regard the fluid mass can be expressed as

[TABLE]

where the relative density $\mathfrak{b}(\cdot)$ is a nonlinear function of the pressure $p$. The variational formulation for the Biot model, under small deformation, reads as follows:

For each $t\in(0,T]$, find ${\mathbf{u}}(t)\in\left(H_{0}^{1}(\Omega)\right)^{d}$, $\mathbf{q}\in H^{1}(\rm{div},\Omega)$ and $p(t)\in L^{2}(\Omega)$ such that there holds

[TABLE]

with the initial condition

[TABLE]

In the above, we have used the standard notations. We denote by $L^{2}(\Omega)$ the space of square integrable functions and by $H^{1}(\Omega)$ the Sobolev space $H^{1}(\Omega)=\{v\in L^{2}(\Omega)\,;\,\nabla\,v\in L^{2}(\Omega)^{d}\}.$ Furthermore, $H^{1}_{0}(\Omega)$ will be the space of functions in $H^{1}(\Omega)$ vanishing on $\partial\Omega$ and $H({\rm div};\Omega)$ the space of vector valued function having all the components and the divergence in $L^{2}(\Omega)$. As usual we denote by $(\cdot,\cdot)$ the inner product in $L^{2}(\Omega)$, and by $||\cdot||$ its associated norm.

Next, we make structural assumptions on the nonlinearities:

• (A1)

$\mathfrak{c},\,\mathfrak{b}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ differentiable with $\mathfrak{c}^{\prime}$ and $\mathfrak{b}^{\prime}$ Lipschitz continuous.

• (A2)

There exists a constant $\alpha_{\mathfrak{c}}$ such that $\mathfrak{c}^{\prime}(\xi)>~{}\alpha_{\mathfrak{c}}$, $\forall\,\xi\in{\mathbb{R}}$.

• (A3)

There exists a constant $\alpha_{\mathfrak{b}}$ such that $\mathfrak{b}^{\prime}(\xi)>\alpha_{\mathfrak{b}}$, $\forall\,\xi\in{\mathbb{R}}$.

• (A4)

There exists constant $k_{m}>0$ and $k_{M}$ such that $k_{m}\leq k(\vec{\xi})\leq k_{M},$ $\forall{\vec{\xi}}\in\Omega$.

For the discretization of problem (30) we use conformal Galerkin finite elements for the displacement variable and mixed finite elements for the flow [23, 43]. More precisely, we use linear elements for the displacement and lowest order Raviart-Thomas elements [16] for the flow. Backward Euler is used for the temporal discretization.

Let $\Omega=\cup_{K\in\mathcal{T}_{h}}K$ be a regular decomposition of $\Omega$ into $d$-simplices. We denote by $h$ the mesh size. The discrete spaces are given by

[TABLE]

where ${\mathbb{P}}_{0},{\mathbb{P}}_{1}$ denote the spaces of constant functions and of linear polynomials, respectively. For $N\in\mathbb{N}$, we discretize the time interval uniformly and define the time step $\tau=\frac{T}{N}$ and $t_{n}=n\tau$. We use the index $n$ for the primary variable $\mathbf{u}^{n}$, $\mathbf{q}^{n}$ and $p^{n}$ at corresponding time step $t_{n}$. In this way, the fully discrete weak problem reads:

For $n\geq 1$ and given $\left(\mathbf{u}_{h}^{n-1},\mathbf{q}_{h}^{n-1},p_{h}^{n-1}\right)$ find $\left(\mathbf{u}_{h}^{n},\mathbf{q}_{h}^{n},p_{h}^{n}\right)\in\left({\mathbf{V}}_{h},{\mathbf{Z}}_{h},W_{h}\right)$, such that

[TABLE]

$\text{for all }\left({\bf v}_{h},{\bf z}_{h},w_{h}\right)\in\left({\mathbf{V}}_{h},{\mathbf{Z}}_{h},W_{h}\right)$.

Following the notation previously introduced, we denote by $n$ the time level, whereas $i$ will refer to the iteration number of the Newton method. We further denote the approximate solution of the linearized problem (36) by $({{\mathbf{u}}}_{h}^{n,i},\mathbf{q}_{h}^{n,i},p_{h}^{n,i})$. At this stage we can introduce the notations

[TABLE]

These will be used subsequently in the convergence analysis of the monolithic Newton method and the alternate version. For the monolithic and splitting $L$-scheme the convergence analysis can be found in [11].

4.1 Convergence analysis of the monolithic Newton method

In this section, we analyse the monolithic Newton method introduced in Section 3 used for solving the simplified nonlinear Biot model given in (36). As we have previously stated, we perform the analysis for the case of small deformation. Here we present a variational formulation of the scheme and demonstrate its quadratic convergence in a rigorous manner. The Newton scheme reads as follows:

For $i=1,2,\ldots$ solve

[TABLE]

$\forall\left({\mathbf{v}}_{h},{\mathbf{z}}_{h},w_{h}\right)\in\left({\mathbf{V}}_{h},{\mathbf{Z}}_{h},W_{h}\right)$, where the initial approximation $({\mathbf{u}}_{h}^{n,0}\mathbf{q}_{h}^{n,0},p_{h}^{n,0})$ is taken as the solution at the previous time step, that is $({\mathbf{u}}_{h}^{n-1},\mathbf{q}_{h}^{n-1},p_{h}^{n-1})$.

In order to prove the convergence of the considered Newton method, the following lemmas will be used.

Lemma 1**.**

Let $\left\{x_{n}\right\}_{n\geq 0}$ be a sequence of real positive number satisfying

[TABLE]

where $a,b\geq 0$. Assuming that

[TABLE]

holds, then the sequence $\left\{x_{n}\right\}_{n\geq 0}$ converges to zero.

Proof.

The result can be shown by induction, see page 52 in [46] for more details. ∎

Lemma 2**.**

If $f:{\mathbb{R}}\rightarrow{\mathbb{R}}$ is differentiable and $f^{\prime}$ is Lipschitz continuous, then there holds

[TABLE]

Proof.

See page $350$ in [32], for example. ∎

Next, the following result provides the quadratic convergence of the Newton method (42) for $\tau$ sufficiently small.

Theorem 1**.**

Assuming (A1)-(A4), the Newton method in (42) converges quadratically if $\tau=O(h^{d})$.

Proof.

By subtracting equations (36) from (42), taking as test functions ${\mathbf{e}}_{\mathbf{u}}^{n,i}$, ${{\mathbf{e}}}_{\mathbf{q}}^{n,i}$ and $e_{p}^{n,i}$ and rearranging some terms to the right hand side we obtain,

[TABLE]

where we have rewritten,

[TABLE]

We obtain an analogous expression for the term with $\mathfrak{b}^{\prime}(\cdot)$. From (A1), $\mathfrak{c}(\cdot)$ is differentiable with $\mathfrak{c}^{\prime}(\cdot)$ Lipschitz continuous, then from Lemma 2 we have,

[TABLE]

where $L_{\mathfrak{c}^{\prime}}$ represents the Lipshitz constant of $\mathfrak{c}^{\prime}(\cdot)$. Then, by using Young’s inequality $(a,b)\leq\displaystyle\frac{||a||^{2}}{2\gamma}+\frac{\gamma||b||^{2}}{2}$, for $\gamma\geq 0$, and by choosing $x=\nabla\cdot{\mathbf{u}}_{h}^{n}$ and $y=\nabla\cdot{\mathbf{u}}_{h}^{n,i-1}$ in (50), from (44) we obtain the following bound, for any $\gamma\geq 0$

[TABLE]

Next, by using the inverse inequality for discrete spaces $||\cdot||_{L^{4}(\Omega)}\leq Ch^{-d/4}||\cdot||$ [41], (pg. 111) the latter reads,

[TABLE]

Finally, by using (A2) and choosing $\gamma=\alpha_{\mathfrak{c}}$, we obtain the following inequality,

[TABLE]

In a similar way, we obtain the following expression from (46),

[TABLE]

Adding  (57), (58), and (45) multiplied by $\tau$ yields,

[TABLE]

By defining $\alpha_{\mathfrak{c},\mathfrak{b}}=\min\left(\alpha_{\mathfrak{c}},\alpha_{\mathfrak{b}},\frac{\tau}{k_{M}}\right)$ and $C_{\mathfrak{c},\mathfrak{b}}=\max\left(\frac{C_{1}L_{\mathfrak{c}^{\prime}}^{2}}{\alpha_{\mathfrak{c}}},\frac{C_{2}L_{\mathfrak{b}^{\prime}}^{2}}{\alpha_{\mathfrak{b}}}\right)$ we can rewrite (61) as

[TABLE]

Using $\lVert\nabla\cdot{\mathbf{e}}_{\mathbf{u}}^{n,0}\rVert\leq C\tau$, $\lVert{e}_{p}^{n,0}\rVert\leq C\tau$ (which can be proven) and Lemma 1, the quadratic convergence of Newton’s method is ensured if

[TABLE]

which holds true for $\tau=O(h^{\frac{d}{2}})$. ∎

4.2 Convergence analysis of the alternate splitting Newton scheme

In this section we present the splitting Newton scheme for solving the nonlinear Biot model given in (36). We present the solver in a variational form and demonstrate its linear convergence.

Let $i\geq 1$, $L_{s}\geq 0$ and $({\mathbf{u}}_{h}^{n,i-1},\mathbf{q}_{h}^{n,i-1},p_{h}^{n,i-1})\in\left({\mathbf{V}}_{h},{\mathbf{Z}}_{h},W_{h}\right)$ be given.

Step 1: find $(\mathbf{q}_{h}^{n,i},p_{h}^{n,i})\in\left({\mathbf{Z}}_{h},W_{h}\right)$ such that

[TABLE]

$\forall\left({\mathbf{z}}_{h},w_{h}\right)\in\left({\mathbf{Z}}_{h},W_{h}\right)$

Step 2: find ${\mathbf{u}}_{h}^{n,i}\in{\mathbf{V}}_{h}$ such that

[TABLE]

$\forall{\mathbf{v}}_{h}\in{\mathbf{V}}_{h}$.

Theorem 2**.**

Assuming (A1)-(A4) and $L_{s}\geq\frac{\alpha^{2}}{\alpha_{\mathfrak{b}}}$, the alternate Newton splitting method in (66)-(69) converges linearly if $\tau$ is small enough.

Proof.

The proof is similar to that of Theorem 1. Nevertheless, for the sake of completion we give it in Appendix A. ∎

5 Numerical examples

In this section, we present numerical experiments that illustrate the performance of the proposed iterative schemes. We study two test problems: a 2D academic problem with a manufactured analytical solution, and a 3D large deformation case on a unit cube. All numerical experiments were implemented using the open-source finite element library Deal II [4]. For all numerical experiments, a Backward Euler scheme has been used for the time discretization. We consider continuous linear Galerkin FE for $\mathbf{u}$, lowest order of Raviart-Thomas FE and discontinuous Galerkin FE for $\mathbf{q}$ and $p$. However, we would like to mention that any stable discretization can be considered instead. For all cases, as stopping criterion for the schemes, we use

[TABLE]

Test problem 1: an academic example for Biot’s model under small deformation

We solve the nonlinear Biot problem under small deformation in the unit-square $\Omega=(0,1)^{2}$ and until final time $T=1$. This test case was proposed in [11] to study the performance of the monolithic and splitting $L$-scheme. We extend the Newton method and the alternate Newton method described in Section 4.

Here, we introduce a manufactured right hand side such that the problem admits the following analytical solution

[TABLE]

which has homogeneous boundary values for $p$ and $\mathbf{u}$.

For infinitesimal deformations and rotations, there is no distinction between the reference and the deformed domains. In this regard, we solve problem (36) using the iterative schemes proposed in Section 4. The mesh size and the time step are set as $h=\tau=0.1$. For this case, all initial conditions are zero. The linearization parameters $L_{p}$ and $L_{\mathbf{u}}$ are equal to the Lipschitz constant $L_{\mathfrak{b}}$ and $L_{\mathfrak{c}}$ corresponding to the nonlinearities $\mathfrak{b}(\cdot)$ and $\mathfrak{c}(\cdot)$ [11].

In order to study the performance of the considered schemes, we propose four coefficient functions for $\mathfrak{b(\cdot)}$ and two for $\mathfrak{c(\cdot)}$, and define four test cases as given in Table 1. Figure 1 shows the performance of the numerical methods at the last time step $T=1$. The monolithic Newton method shows quadratic convergence in all cases. Nevertheless, the alternate Newton and the $L$-scheme methods show linear convergence as predicted in Section 4.

Figure 2 shows the performance of the considered schemes for different time steps. The Newton method has better convergence for smaller time steps while the $L$-scheme has it for larger time steps; all this is in agreement with the Theorems 1 and 2. The performance of the considered schemes are independent of the mesh discretization.

Test problem 2: a unit cube under large deformation

We now solve a large deformation problem on the unit-cube $\Omega=(0,1)^{3}$. A Lagrangian frame of reference is necessary to keep track of the deformed domain $\Omega_{t}$ at time $t$. We study the performance of the iterative schemes presented in Section 3 for solving Eqs. (4). The material is supposed to be isotropic and with constant Lamé parameters $\mu$ and $\mathfrak{c}(\cdot)$. We consider a Lagrangian fluid mass $m_{f}=\rho_{f}J\phi$ of a slightly compressible fluid, where $\phi$ is the porosity. Under this assumption, the time derivative of the fluid content reads as

[TABLE]

where the compressibility $c_{p}$ and Biot’s coefficient $c_{\alpha}=J\frac{\partial\phi}{\partial J}+\phi\approx 1$ for simplicity. We will compare the iterative schemes for a torsion case on a unit cube. On the top face, we apply the rotation tensor $\mathbf{R}(\theta)$ of a time dependent angle $\theta(t)=\pi/4\ t$, which gives a rotation of $\pi/4$ at $T=1$. We set homogeneous initial condition for $(\mathbf{q}_{0},p_{0})$ and $\nabla\mathbf{u}_{0}=\left(\mathbf{R}(\theta)-\mathbf{I}\right)$. In the alternate Newton method, the stabilization parameter is set to $L_{s}=1.$ In the L-scheme method, the linearisation tensor parameters are set as follows: ${\bf{L}}_{\mathbf{u}}=\mathcal{\partial}_{\mathbf{u}}{\bf\Pi}\left(\nabla\mathbf{u}_{0},p_{0}\right),$ ${\bf{L}}_{p}=\mathcal{\partial}_{p}{\bf\Pi}\left(\nabla\mathbf{u}_{0},p_{0}\right),$ $\ {\bf{L}}_{\mathbf{q}}=\mathcal{\partial}_{p}{\bf K}\left(\nabla\mathbf{u}_{0}\right),$ $L_{p}=\mathcal{\partial}_{p}{\bf\Gamma}\left(\nabla\mathbf{u}_{0},p_{0}\right)$ and $L_{\mathbf{u}}=\mathcal{\partial}_{\mathbf{u}}{\bf\Gamma}\left(\nabla\mathbf{u}_{0},p_{0}\right)$. The mesh size and the time step are set as $h=\tau=2^{-3}$. We denote by top face of the unit-cube the region $z=1$, the bottom face $z=0$ and the lateral faces are $x=0$, $x=1$, $y=0$ and $y=1$. The boundary conditions are listed in Table 2 and the displacement and pressure field are shown in Figure 4.

We compare the performance of the schemes proposed in Section 3 and we observe that the numerical convergence is in accordance with the theory developed in Section 4, even though the analysis is done for small deformation. Newton’s method has quadratic convergence for the smaller time steps and linear convergence for the larger time steps. In contrast, the monolithic $L$-scheme has the same rate of convergence regardless of size of the time step (see Figure 4). All splitting schemes have better convergence when the stability term is used (we use $L_{s}=1.0$).

6 Conclusions

We considered Biot’s model under small and large deformation. Different nonlinear solvers based on the $L$-scheme, Newton’s method, and the undrained splitting method were presented. The only quadratic convergent scheme is the monolithic Newton method. The splitting Newton method also requires a stabilization parameter, otherwise the (linear) convergence cannot be guaranteed. The analysis of the schemes and illustrative numerical experiments were presented.

We tested the performance of the schemes on two test problems: a unit square under small deformation and a unit cube under large deformation. To summarise, we make the following remarks:

• $\bullet$

Monolithic and splitting $L$-schemes are robust with respect to the choice of the linearization parameter, the mesh size, and time step size.

• $\bullet$

The stabilization parameter $L_{s}$ has a strong influence on the speed of the convergence of the splitting Newton scheme.

• $\bullet$

The splitting $L$-scheme can be used both as a robust solver or even as a preconditioner (as it is established in [26, 57]) to improve the performance of the monolithic Newton method and the $L$-scheme.

Appendix A Convergence proof of the alternate Newton method

The following result provides the linear convergence of the alternate Newton method in (66)- (69) for $\tau$ sufficiently small.

Theorem 3**.**

Assuming (A1)-(A4) and $L_{s}\geq\frac{\alpha^{2}}{\alpha_{\mathfrak{b}}}$, the alternate Newton splitting method in (66)-(69) converges linearly if $\tau$ is small enough.

Proof.

By subtracting problems (66)-(69) and (36), taking as test functions ${\mathbf{e}}_{\mathbf{q}}^{n,i}$, $e_{p}^{n,i}$ and ${\mathbf{e}}_{\mathbf{u}}^{n,i}$, and rearranging some elements to the right hand side we obtain,

[TABLE]

The mechanics equation then gives,

[TABLE]

By using similar steps as in Theorem 1, we obtain the following

[TABLE]

Next, by using the inverse inequality $||\cdot||_{L^{4}(\Omega)}\leq Ch^{-d/4}||\cdot||$ [41], and by using the following formula $(x-y,x)=\displaystyle\frac{||x||^{2}}{2}+\frac{||x-y||^{2}}{2}-\frac{||y||^{2}}{2}$, by choosing $x=\nabla\cdot{{\mathbf{e}}_{\mathbf{u}}}^{n,i}$ and $y=\nabla\cdot{\mathbf{e}}_{\mathbf{u}}^{n,i-1},$ we obtain from (82)

[TABLE]

Finally, by reorganizing (86), using (A2) and choosing $\gamma_{1}=\alpha_{\mathfrak{c}}$, we obtain the following inequality,

[TABLE]

In a similar way, we obtain the following expression from (46),

[TABLE]

Adding equations (89) and (90) yields,

[TABLE]

By using Young’s inequality $(a,b)\leq\displaystyle\frac{||a||^{2}}{2\gamma}+\frac{\gamma||b||^{2}}{2}$, for $\gamma>0$ and choosing $b=e_{p}^{n,i}$ and $a=\nabla\cdot\delta{\mathbf{e}}_{\mathbf{u}}^{n,i}$ we bound the coupling term (for $\gamma_{2}>0$),

[TABLE]

Then by using (95) and choosing $\gamma_{2}=\frac{\alpha_{\mathfrak{b}}}{2}$ we obtain from (94)

[TABLE]

Since $L_{s}\geq\frac{\alpha^{2}}{\alpha_{\mathfrak{b}}}$, we obtain

[TABLE]

By using $\lVert\nabla\cdot{\mathbf{e}}_{\mathbf{u}}^{n,0}\rVert\leq C\tau$, $\lVert{e}_{p}^{n,0}\rVert\leq C\tau$ wich can be proven and the estimate in Lemma 1, the convergence is ensured if $\tau=O(h^{\frac{d}{2}})$.

∎

Acknowledgement

The research was supported by the University of Bergen in cooperation with the FME-SUCCESS center (grant 193825/S60) funded by the Research Council of Norway. The work has also been partly supported by: the NFR-DAADppp grant 255715, the NFR-Toppforsk project 250223, the NRC-CHI grant 255510 and the NRC-IMMENS grant 255426.