Unconditionally stable second order convergent partitioned methods for multiple-network poroelasticity

TL;DR

This paper introduces unconditionally stable, second-order convergent partitioned numerical methods for solving quasi-static multiple-network poroelasticity equations, enabling efficient and accurate simulations without iterative coupling.

Contribution

The paper presents novel partitioned methods that split the equations into subproblems, achieving unconditional stability and high-order convergence without iterative coupling.

Findings

Methods are unconditionally stable.

Achieve second-order convergence in time.

Numerical results confirm good performance.

Abstract

In this paper, we consider partitioned numerical methods for quasi-static multiple-network poroelasticity (MPET) equations, generalizations of the Biot model in poroelasticity for multiple pore networks. Two partitioned numerical methods are presented for the equations which split time discretization into solving two subequations, a Lame equation and a system of heat equations, alternatively. In contrast to the iterative coupling methods which require multiple iterations at each time step, our numerical methods solve these smaller equations only once at each time step. We prove their unconditional stability and high order convergence in time with a novel error analysis. A number of numerical results are presented to illustrate good performances of these partitioned methods.

Tables

Table 1. Table 1 . The errors and convergence rates given by the diffusion-then-elasticity scheme with k = 1 𝑘 1 k=1

$M$	${\Vert 𝒖-{𝒖}_h\Vert }_1$		${\Vert p_t-p_{t,h}\Vert }_0$		${\Vert p_1-p_{1,h}\Vert }_1$		${\Vert p_2-p_{2,h}\Vert }_1$
	error	rate	error	rate	error	rate	error	rate
8	1.290e+0	-	2.146e-1	-	2.661e-1	-	5.323e-1	-
16	3.195e-1	2.01	3.898e-2	2.46	1.865e-1	0.51	3.729e-1	0.51
32	7.700e-2	2.05	8.856e-3	2.14	1.059e-1	0.82	2.118e-1	0.82
64	1.872e-2	2.04	2.154e-3	2.04	5.599e-2	0.92	1.120e-1	0.92
128	4.603e-3	2.02	5.333e-4	2.01	2.873e-2	0.96	5.747e-2	0.96

Table 2. Table 2 . The errors and convergence rates given by the elsticity-then-diffusion scheme with k = 1 𝑘 1 k=1

$M$	${\Vert 𝒖-{𝒖}_h\Vert }_1$		${\Vert p_t-p_{t,h}\Vert }_0$		${\Vert p_1-p_{1,h}\Vert }_1$		${\Vert p_2-p_{2,h}\Vert }_1$
	error	rate	error	rate	error	rate	error	rate
8	1.290e+0	-	2.146e-1	-	2.661e-1	-	5.323e-1	-
16	3.195e-1	2.01	3.898e-2	2.46	1.865e-1	0.51	3.729e-1	0.51
32	7.700e-2	2.05	8.856e-3	2.14	1.059e-1	0.82	2.118e-1	0.82
64	1.872e-2	2.04	2.154e-3	2.04	5.599e-2	0.92	1.120e-1	0.92
128	4.603e-3	2.02	5.333e-4	2.01	2.873e-2	0.96	5.747e-2	0.96

Table 3. Table 3 . The errors and convergence rates given by the diffusion-then-elasticity scheme with k = 2 𝑘 2 k=2

$M$	${\Vert 𝒖-{𝒖}_h\Vert }_1$		${\Vert p_t-p_{t,h}\Vert }_0$		${\Vert p_1-p_{1,h}\Vert }_1$		${\Vert p_2-p_{2,h}\Vert }_1$
	error	rate	error	rate	error	rate	error	rate
8	2.682e-1	-	3.405e-2	-	4.082e-2	-	8.165e-2	-
16	3.153e-2	3.09	3.615e-3	3.24	1.440e-2	1.50	2.880e-2	1.50
32	3.698e-3	3.09	4.082e-4	3.15	4.098e-3	1.81	8.196e-3	1.81
64	4.451e-4	3.05	4.865e-5	3.07	1.084e-3	1.92	2.168e-3	1.92
128	5.454e-5	3.03	5.943e-6	3.03	2.781e-4	1.96	5.563e-4	1.96

Table 4. Table 4 . The errors and convergence rates given by the elsticity-then-diffusion scheme with k = 2 𝑘 2 k=2

$M$	${\Vert 𝒖-{𝒖}_h\Vert }_1$		${\Vert p_t-p_{t,h}\Vert }_0$		${\Vert p_1-p_{1,h}\Vert }_1$		${\Vert p_2-p_{2,h}\Vert }_1$
	error	rate	error	rate	error	rate	error	rate
8	2.682e-1	-	3.405e-2	-	4.082e-2	-	8.165e-2	-
16	3.153e-2	3.09	3.615e-3	3.24	1.440e-2	1.50	2.880e-2	1.50
32	3.698e-3	3.09	4.082e-4	3.15	4.098e-3	1.81	8.196e-3	1.81
64	4.451e-4	3.05	4.865e-5	3.07	1.084e-3	1.92	2.168e-3	1.92
128	5.454e-5	3.03	5.943e-6	3.03	2.781e-4	1.96	5.563e-4	1.96

Table 5. Table 5 . The errors and convergence rates given by the diffusion-then-elasticity scheme with k = 3 𝑘 3 k=3

$M$	${\Vert 𝒖-{𝒖}_h\Vert }_1$		${\Vert p_t-p_{t,h}\Vert }_0$		${\Vert p_1-p_{1,h}\Vert }_1$		${\Vert p_2-p_{2,h}\Vert }_1$
	error	rate	error	rate	error	rate	error	rate
8	4.942e-2	-	8.388e-3	-	4.240e-3	-	8.479e-3	-
16	3.108e-3	3.99	4.581e-4	4.19	7.292e-4	2.54	1.458e-3	2.54
32	1.888e-4	4.04	2.626e-5	4.12	1.058e-4	2.78	2.114e-4	2.79
64	1.150e-5	4.04	1.559e-6	4.07	1.556e-5	2.77	3.092e-5	2.77
128	7.069e-7	4.02	9.467e-8	4.04	2.719e-6	2.52	5.280e-6	2.55

Table 6. Table 6 . The errors and convergence rates given by the elsticity-then-diffusion scheme with k = 3 𝑘 3 k=3

$M$	${\Vert 𝒖-{𝒖}_h\Vert }_1$		${\Vert p_t-p_{t,h}\Vert }_0$		${\Vert p_1-p_{1,h}\Vert }_1$		${\Vert p_2-p_{2,h}\Vert }_1$
	error	rate	error	rate	error	rate	error	rate
8	4.942e-2	-	8.388e-3	-	4.240e-3	-	8.479e-3	-
16	3.108e-3	3.99	4.581e-4	4.19	7.292e-4	2.54	1.458e-3	2.54
32	1.888e-4	4.04	2.626e-5	4.12	1.058e-4	2.78	2.114e-4	2.79
64	1.150e-5	4.04	1.559e-6	4.07	1.556e-5	2.77	3.092e-5	2.77
128	7.069e-7	4.02	9.467e-8	4.04	2.719e-6	2.52	5.280e-6	2.55

Equations

$$-\mathrm{div}\mathcal{C}\epsilon (\mathbf{u})+\sum _{i=1}^N{\alpha }_i\nabla p_i$$

$$s_i{\dot{p}}_i+{\alpha }_i\mathrm{div}\dot{\mathbf{u}}-\mathrm{div}\left(K_i\nabla p_i\right)+{\xi }_i(\mathbf{p})$$

$$\mathcal{C}\epsilon (\mathbf{u})=2\mu \epsilon (\mathbf{u})+\lambda \mathrm{div}\mathbf{u}\mathbb{I}$$

$${\xi }_i(\mathbf{p})=\sum _{i=1}^N{\xi }_{j\leftarrow i}(p_j-p_i),$$

$$C^m([0,T];\mathcal{X})=\{f|{\partial }^{i}f/\partial t^i\in C^0([0,T];\mathcal{X}),0\le i\le m\},$$

$$\Vert f{\Vert }_{L^r((0,T);\mathcal{X})}=\{\begin{array}{l}{\left({\int }_0^T\Vert f(s){\Vert }_{\mathcal{X}}^r\mathrm{d}s\right)}^{1/r},1\le r<\infty ,\\ {\mathrm{esssup}}_{t\in (0,T)}\Vert f(t){\Vert }_{\mathcal{X}},r=\infty .\end{array}$$

$$K_i^{-}\le K_i(x)\le K_i^{+}\text{a.e.}x\in \Omega .$$

$$\mathbf{u}$$

$$\left(\mathcal{C}\epsilon (\mathbf{u})-\sum _{i=1}^N{\alpha }_i{p}_i\mathbb{I}\right)\mathbf{n}$$

$$p_i$$

$$p_i(x,0)=p_i^0(x),x\in \Omega ,i=1,\dots ,N.$$

$$\mathrm{div}\mathbf{u}={\lambda }^{-1}\left(-p_t+\mathbf{\alpha }\cdot \mathbf{p}\right).$$

$$-\mathrm{div}\left(2\mu \epsilon (\mathbf{u})-p_t\mathbb{I}\right)$$

$$-\mathrm{div}\mathbf{u}-{\lambda }^{-1}{p}_t+{\lambda }^{-1}\mathbf{\alpha }\cdot \mathbf{p}$$

$$s_i{\dot{p}}_i+{\alpha }_i{\lambda }^{-1}(-{\dot{p}}_t+\mathbf{\alpha }\cdot \dot{\mathbf{p}})-\mathrm{div}(K_i\nabla p_i)+{\xi }_i(\mathbf{p})$$

$$\mathbf{V}=H_{0,{\Gamma }_d}^1(\Omega ;{\mathbb{R}}^d),Q_t=L^2(\Omega ),Q_i=H_0^1(\Omega )i=1,\dots ,N.$$

$${\Vert \mathbf{v}\Vert }_{\mathbf{V}}^2=⟨2\mu \epsilon (\mathbf{v}),\epsilon (\mathbf{v})⟩,{\Vert q_t\Vert }_{Q_t}^2=⟨(2\mu {)}^{-1}{q}_t,q_t⟩.$$

$$⟨2\mu \epsilon (\mathbf{u}),\epsilon (\mathbf{v})⟩-⟨p_t,\mathrm{div}\mathbf{v}⟩$$

$$-⟨\mathrm{div}\mathbf{u},q_t⟩-⟨{\lambda }^{-1}{p}_t,q_t⟩+⟨{\lambda }^{-1}\mathbf{\alpha }\cdot \mathbf{p},q_t⟩$$

$$⟨s_i{\dot{p}}_i+{\alpha }_i{\lambda }^{-1}(-{\dot{p}}_t+\mathbf{\alpha }\cdot \dot{\mathbf{p}}+{\xi }_i(\mathbf{p}),q_i⟩+⟨K_i\nabla p_i,\nabla q_i⟩$$

$${\sigma }^{n+\frac{1}{2}}=\frac{{\sigma }^n+{\sigma }^{n+1}}{2}.$$

$$\underset{\mathbf{v}\in {\mathbf{V}}_h}{\sup }\frac{⟨\mathrm{div}\mathbf{v},q⟩}{{\Vert \mathbf{v}\Vert }_{\mathbf{V}}}\ge C{\Vert q\Vert }_{Q_t}\forall q\in Q_{t,h}$$

$$\Vert q_i{\Vert }_{a_{h,i}}^2=a_{h,i}(q_i,q_i).$$

$$\underset{\mathbf{v}\in {\mathbf{V}}_h}{\inf }\Vert \mathbf{u}-\mathbf{v}{\Vert }_{\mathbf{V}}\lesssim h^{k_{\mathbf{u}}}\Vert \mathbf{u}{\Vert }_{k_{\mathbf{u}}+1},\underset{q_t\in Q_t}{\inf }\Vert p_t-q_t{\Vert }_{Q_t}\lesssim h^{k_{\mathbf{u}}}\Vert p_t{\Vert }_{k_{\mathbf{u}}}$$

$$\underset{q_i\in Q_{i,h}}{\inf }\Vert p_i-q_i{\Vert }_{a_{h,i}}\lesssim h^{k_p}\Vert p_i{\Vert }_{k_p+1}1\le i\le N$$

$$S(\mathbf{p},\mathbf{q}):=\sum _{i=1}^N⟨s_i{p}_i,q_i⟩+⟨{\lambda }^{-1}\mathbf{\alpha }\cdot \mathbf{p},\mathbf{\alpha }\cdot \mathbf{q}⟩,$$

$$A(\mathbf{p},\mathbf{q}):=\sum _{i=1}^N\left(⟨{\xi }_i(\mathbf{p}),q_i⟩+a_{h,i}(p_i,q_i)\right),$$

$$⟨2\mu \epsilon ({\mathbf{u}}^n),\epsilon (\mathbf{v})⟩-⟨p_t^n,\mathrm{div}\mathbf{v}⟩=⟨{\mathbf{f}}^n,\mathbf{v}⟩,$$

$$-⟨\mathrm{div}{\mathbf{u}}^n,q_t⟩-⟨{\lambda }^{-1}{p}_t^n,q_t⟩=-⟨{\lambda }^{-1}\mathbf{\alpha }\cdot {\mathbf{p}}^n,q_t⟩,$$

$$-⟨s_i{\dot{p}}_i^{n+\frac{1}{2}},q_i⟩-⟨{\alpha }_i{\lambda }^{-1}\mathbf{\alpha }\cdot {\dot{\mathbf{p}}}^{n+\frac{1}{2}}+{\xi }_i\left({\mathbf{p}}^{n+\frac{1}{2}}\right),q_i⟩$$

$$-a_{h,i}\left(p_i^{n+\frac{1}{2}},q_i\right)=-⟨g^{n+\frac{1}{2}},q_i⟩-⟨{\alpha }_i{\lambda }^{-1}{\dot{p}}_t^{n+\frac{1}{2}},q_i⟩,1\le i\le N,$$

$$⟨2\mu \epsilon ({\Pi }_h^{\mathbf{V}}\mathbf{u}),\epsilon (\mathbf{v})⟩-⟨{\Pi }_h^{Q_t}{p}_t,\mathrm{div}\mathbf{v}⟩$$

$$-⟨\mathrm{div}{\Pi }_h^{\mathbf{V}}\mathbf{u},q_t⟩-⟨{\lambda }^{-1}{\Pi }_h^{Q_t}{p}_t,q_t⟩$$

$$\sum _{i=1}^N\left[a_{h,i}({\Pi }_h^{Q_i}{p}_i,q_i)+⟨{\xi }_i({\Pi }_h^{\mathbf{Q}}\mathbf{p}),q_i⟩\right]=\sum _{i=1}^N\left[a_{h,i}\left(p_i,q_i\right)+⟨{\xi }_i(\mathbf{p}),q_i⟩\right],q_i\in Q_{i,h}$$

$$\sum _{1\le i\le N}⟨{\xi }_i(\mathbf{p}),q_i⟩=\sum _{1\le i,j\le N}\frac{1}{2}⟨{\xi }_{j\leftarrow i}(p_i-p_j),q_i-q_j⟩,$$

$$\sum _{1\le i\le N}⟨{\xi }_i(\mathbf{p}),p_i⟩=\sum _{1\le i,j\le N}\frac{1}{2}⟨{\xi }_{j\leftarrow i}(p_i-p_j),p_i-p_j⟩\ge 0.$$

Full text

1. Introduction

In this paper, we consider partitioned numerical methods for a family of quasi-static multiple-network poroelasticity (MPET111The abbreviation MPET is from the term multiple-network poroelastic theory in literature e.g. [23]. Here, we instead refer to the multiple-network poroelasticity equations but keep the abbreviation for the sake of convenience.) equations reading as follows: for a given number of networks $N\in\mathbb{N}$, find the displacement $\boldsymbol{u}$ and the network pressures $p_{i}$ for $i=1,\dots,N$ such that

[TABLE]

where $\boldsymbol{p}$ is the $N$-tuple $(p_{1},\cdots,p_{N})$. Here $\boldsymbol{u}=\boldsymbol{u}(x,t)$, $p_{i}=p_{i}(x,t)$, $1\leq i\leq N$ are time-dependent functions for $x\in\Omega\subset{\mathbb{R}}^{d}$ ($d=2,3$) and for $t\in[0,T]$. The operators and parameters are as follows: $\mathcal{C}$ is the elastic stiffness tensor, each network $i$ is associated with a Biot-Willis coefficient $\alpha_{i}\in(0,1]$, storage coefficient $s_{i}\geq 0$, and hydraulic conductivity tensor $K_{i}=\kappa_{i}/\mu_{i}>0$ (where $\kappa_{i}$ and $\mu_{i}$ represent the network permeability and the network fluid viscosity, respectively). In (1.1a), $\operatorname{\nabla}$, $\varepsilon$, and $\operatorname{div}$ denote the (column-wise) gradient, the symmetric (row-wise) gradient, and the row-wise divergence, respectively. In (1.1b), $\operatorname{\nabla}$ and $\operatorname{div}$ are the standard gradient and divergence operators, and the superposed dot denotes the time derivative. Further, $\boldsymbol{f}$ represents a body force and $g_{i}$ represents sources in network $i$ for $i=1,\dots,N$, while $\xi_{i}$ represents transfer terms from network $i$ to other networks. From a physical point of view, (1.1) represents the balance of linear momentum and the mass conservation in a porous, linearly elastic medium permeated by $N$ segregated viscous fluid networks.

In this paper, we consider the case of an isotropic stiffness tensor for which

[TABLE]

where $\mu,\lambda$ are the standard non-negative Lamé parameters and $\mathbb{I}$ denotes the identity tensor. Moreover, we will consider the case where the transfer terms $\xi_{i}$, quantifying the transfer out of network $i$ into the other fluid networks, are proportional to pressure differences between the networks. More precisely, we assume that $\xi_{i}$ takes the form:

[TABLE]

where $\xi_{j\leftarrow i}$ are non-negative transfer coefficients for $1\leq i,j\leq N$. We will also assume that these transfer coefficients are symmetric in the sense that $\xi_{j\leftarrow i}=\xi_{i\leftarrow j}$, and note that $\xi_{j\leftarrow j}$ is arbitrary.

The MPET equations have an abundance of both geophysical and biological applications. For example, the case with $N=2$ is the Biot-Barenblatt model which models dual porosity property of poroelastic media. In biomechanics, Tully and Ventikos [23] considered (1.1) with four different networks ($N=4$) as a macroscopic model for the dynamics of fluid flows in brain tissue. The fluid networks represent the arteries, the arterioles/capillaries, the veins and the interstitial fluid-filled extracellular space, and each network may have a different permeability $\kappa_{i}$ and different transfer coefficients $\xi_{j\leftarrow i}$ for $1\leq j\leq N$, $j\not=i$.

While the Biot model has been throughly studied, see e.g. [1, 6, 15, 19, 20, 21, 24], only a few numerical methods for the MPET equations are available. To the best of our knowledge, the first numerical method was proposed in [17], and another method is proposed more recently in [14]. Robust preconditioners of numerical methods for MPET equations can be constructed by extending the block preconditioners of the Biot model (see e.g., [13, 16]). Nevertheless, the problem sizes of the MPET equations are intrinsically large, so partitioned numerical methods for time discretization are valuable approaches to reduce computational costs.

There are only a few of previous studies on partitioned numerical methods for poroelasticity equations. The conditional stability of various partitioned methods for a dynamic poroelasticity model was studied in [10]. In [11] a partitioned numerical method was developed for quasi-static poroelasticity equations using the discontinuous Galerkin method with a stabilization technique. However, the error analysis of the method gives $1/2$ convergence rate of time discretization errors which is not regarded as the optimal order of time discretization errors.

The objective of this paper is to develop and analyze new partitioned numerical schemes for the MPET equations. We propose two partitioned numerical schemes which are unconditionally stable without any stabilization terms and have second or higher order convergence in time. A novel error analysis will be presented to show the stability and convergence of the partitioned schemes, and numerical results will be given to illustrate their performances. We use the formulation proposed in [17] with the total pressure, so the implicit constants in the error estimates are uniformly bounded for arbitrarily large $\lambda>0$ and for small storage coefficients $s_{i}$’s. In contrast to monolithic numerical methods, our partitioned numerical schemes solve two subproblems sequentially, one is a linear elasticity equation and the other is a system of parabolic equations, so required computational resources at each solve can be significantly reduced. We point out that our partitioned schemes are intrinsically different from the iterative coupling methods, which are equivalent to monolithic methods with block triangular preconditioners which use exact LU solvers in diagonal blocks of the preconditioners.

This paper is organized as follows. Section 2 presents notation and general preliminaries. In Section 3, we introduce a variational formulation with the total pressure for the quasi-static MPET equations (1.1). In Section 4 we present two partitioned numerical methods and prove their convergence with the a priori error estimates. We present numerical results in Section 5 to illustrate these theoretical results. Conclusions and future research directions are highlighted in Section 6.

2. Notation and preliminaries

In this paper we use $X\lesssim Y$ to denote an inequality $X\leq CY$ with a generic constant $C>0$ which is independent of mesh sizes. If needed, we will write $C$ explicitly in inequalities but it can vary across expressions.

2.1. Sobolev spaces

Let $\Omega$ be a bounded polyhedral domain in ${\mathbb{R}}^{d}$ ($d=2,3$) with boundary $\partial\Omega$. We let $L^{2}(\Omega)$ be the set of square-integrable real-valued functions on $\Omega$. The inner product of $L^{2}(\Omega)$ and the induced norm are denoted by $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$, respectively. For a finite-dimensional inner product space ${\mathbb{X}}$, typically ${\mathbb{X}}={\mathbb{R}}^{d}$, let $L^{2}(\Omega;{\mathbb{X}})$ be the space of ${\mathbb{X}}$-valued functions such that each component is in $L^{2}(\Omega)$. The inner product of $L^{2}(\Omega;{\mathbb{X}})$ is naturally defined by the inner product of ${\mathbb{X}}$ and $L^{2}(\Omega)$, so we use the same notation $\langle\cdot,\cdot\rangle$ and $\|\cdot\|_{0}$ to denote the inner product and norm on $L^{2}(\Omega;{\mathbb{X}})$.

For a non-negative integer $m$, $H^{m}(\Omega)$ denotes the standard Sobolev spaces of real-valued functions based on the $L^{2}$-norm, and $H^{m}(\Omega;{\mathbb{X}})$ is defined similarly based on $L^{2}(\Omega;{\mathbb{X}})$. To avoid confusion with the weighted $L^{2}$-norms we use $\|\cdot\|_{m}$ to denote the $H^{m}$-norm (both for $H^{m}(\Omega)$ and $H^{m}(\Omega;{\mathbb{X}})$). For $m\geq 1$, we use $H^{m}_{0,\Gamma}(\Omega)$ to denote the subspace of $H^{m}(\Omega)$ with vanishing trace on $\Gamma\subset\partial\Omega$, and $H^{m}_{0,\Gamma}(\Omega;{\mathbb{X}})$ is defined similarly [12]. For $\Gamma=\partial\Omega$, we write $H^{m}_{0}(\Omega)$ and analogously $H^{m}_{0}(\Omega;{\mathbb{X}})$.

2.2. Spaces involving time

For $T>0$ and a reflexive Banach space $\mathcal{X}$, let $C^{0}([0,T];\mathcal{X})$ denote the set of functions $f:[0,T]\rightarrow\mathcal{X}$ that are continuous in $t\in[0,T]$. For an integer $m\geq 1$, we define

[TABLE]

where $\partial^{i}f/\partial t^{i}$ is the $i$-th time derivative in the sense of the Fréchet derivative in $\mathcal{X}$ (see e.g. [25]).

For a function $f:[0,T]\rightarrow\mathcal{X}$, the Bochner norm is defined as

[TABLE]

We define $W^{k,r}((0,T);\mathcal{X})$ for a non-negative integer $k$ and $1\leq r\leq\infty$ as the closure of $C^{k}([0,T];\mathcal{X})$ with the norm $\|f\|_{W^{k,r}((0,T);\mathcal{X})}=\sum_{i=0}^{k}\|\partial^{i}f/\partial t^{i}\|_{L^{r}((0,T);\mathcal{X})}$.

2.3. Finite element spaces

Let $\mathcal{T}_{h}$ be a shape-regular triangulation of $\Omega$. For any integer $k\geq 1$, we let $\mathcal{P}_{k}(\mathcal{T}_{h})$ denote the space of continuous piecewise polynomials of order $k$ associated to $\mathcal{T}_{h}$, and $\mathcal{P}_{k}(\mathcal{T}_{h};{\mathbb{R}}^{d})$ as the space of $d$-tuples with components in $\mathcal{P}_{k}(\mathcal{T}_{h})$. We will omit $\mathcal{T}_{h}$ when it is clear in context.

2.4. Parameter values

We will make the following assumptions on the material parameter values. First, we assume that the Biot-Willis coefficients $\alpha_{i}\in(0,1]$, $i=1,\dots,N$, and the storage coefficients $s_{i}>0$ are constant in time for $i=1,\dots,N$. In the analysis, we will pay particular attention to robustness of estimates with respect to arbitrarily large $\lambda$ because large $\lambda$ correspond to nearly incompressible materials which are common in biomechanical modelling.

We will assume that the hydraulic conductivities $K_{i}$ are constant in time, but possibly spatially-varying and that these satisfy standard ellipticity constraints: i.e. there exist positive constants $K_{i}^{-}$ and $K_{i}^{+}$ such that

[TABLE]

We assume that the transfer coefficients $\xi_{j\leftarrow i}$ are constant in time and non-negative, i.e., $\xi_{j\leftarrow i}(x)\geq 0$ for $1\leq i,j\leq N$, $x\in\Omega$.

2.5. Boundary conditions

We will consider (1.1) augmented by the following standard boundary conditions. First, we assume that the boundary decomposes into two parts: $\partial\Omega=\Gamma_{d}\cup\Gamma_{t}$ with $|\Gamma_{d}\cap\Gamma_{t}|=0$ and $|\Gamma_{d}|,|\Gamma_{t}|>0$ where $|\Gamma|$ is the $(d-1)$-dimensional Lebesgue measure of $\Gamma\subset\partial\Omega$. We use $\boldsymbol{n}$ to denote the outward unit normal vector field on $\partial\Omega$. Relative to this partition, we consider the homogeneous boundary conditions

[TABLE]

The subsequent formulation and analysis can easily be extended to cover inhomogeneous and other types of boundary conditions with suitable modifications.

2.6. Initial conditions

The MPET equations (1.1) need appropriate initial conditions. In particular, in agreement with the assumption that $c_{i}>0$ for $i=1,\dots,N$, we assume that initial conditions are given for all $p_{i}$:

[TABLE]

Given such $p_{i}^{0}$, we note that we may compute $\boldsymbol{u}(x,0)=\boldsymbol{u}^{0}(x)$ from (1.1a), which in particular yields a $\operatorname{div}\boldsymbol{u}(x,0)=\operatorname{div}\boldsymbol{u}^{0}(x)$ for $x\in\Omega$. In the following, we will assume that any initial conditions given are compatible in the sense described here.

3. The formulation with total pressure

In this section, we review the formulation for the quasi-static multiple-network poroelasticity equations with the total pressure. In order to be consistent with the pressure definitions in solid mechanics total pressure definition in this paper and the one in [17] have different sign. In the subsequent subsections, we present the augmented governing equations and introduce a corresponding variational formulation.

3.1. Governing equations introducing the total pressure

Let $\boldsymbol{u}$ and $p_{i}$ for $i=1,\dots,N$ be solutions of (1.1) with boundary conditions given by (2.1), initial conditions given by (2.2) and recall the isotropic stiffness tensor assumption, cf. (1.2).

For simplicity, we denote $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{N})$ and $\boldsymbol{p}=(p_{1},\dots,p_{N})$, and we write $\boldsymbol{\alpha}\cdot\boldsymbol{p}=\operatorname{\textstyle\sum}_{i=1}^{N}\alpha_{i}p_{i}$. Introducing the total pressure $p_{t}$ defined as $-p_{t}=\lambda\operatorname{div}\boldsymbol{u}-\boldsymbol{\alpha}\cdot\boldsymbol{p}$, we have

[TABLE]

Inserting (3.1) and its time-derivative into (1.1b), we obtain an augmented system of quasi-static multiple-network poroelasticity equations: for $t\in(0,T]$, find the displacement vector field $\boldsymbol{u}$ and the pressure scalar fields $p_{t}$ and $\boldsymbol{p}$ such that

[TABLE]

We note that $p_{t}(x,0)$ can be computed from (2.2) and the definition of $p_{t}$.

3.2. Variational formulation

For $\Gamma_{d}\subset\partial\Omega$ satisfying $0<|\Gamma_{d}|<|\partial\Omega|$, using the notations introduced in Section 2, let

[TABLE]

If $\Gamma_{d}=\partial\Omega$, then $Q_{t}=L_{0}^{2}(\Omega)$, the space of all mean-value zero functions in $L^{2}(\Omega)$. If $\Gamma_{d}=\emptyset$, then $\boldsymbol{V}$ must be the subspace of $H^{1}(\Omega;{\mathbb{R}}^{d})$ orthogonal to the space of rigid motions on $\Omega$. For simplicity of presentation we will assume that $0<|\Gamma_{d}|<|\partial\Omega|$ in the rest of this paper. We also use $\boldsymbol{Q}=Q_{1}\times\dots\times Q_{N}$ for simplicity. Here we define the norms of $\left\|\cdot\right\|_{\boldsymbol{V}}$, $\left\|\cdot\right\|_{Q_{t}}$ as

[TABLE]

Throughout the paper we assume that $\mu\lesssim\lambda$, therefore $\|q_{t}\|_{\lambda^{-1}}\leq C\|q_{t}\|_{Q_{t}}$ holds with a constant $C$ which is uniformly bounded above for arbitrarily large $\lambda$.

Multiplying (3.2) by test functions and integrating by parts with boundary conditions given by (2.1) and initial conditions given by (2.2) yield the following variational formulation: given a compatible initial data $(\boldsymbol{u}^{0},p_{t}^{0},\boldsymbol{p}^{0})$ satisfying (3.2a), (3.2b), and given $\boldsymbol{f}$ and $g_{i}$’s for $i=1,\dots,N$, find $\boldsymbol{u}\in C^{0}([0,T];\boldsymbol{V})$, $p_{t}\in C^{0}([0,T];Q_{t})\cap C^{1}((0,T],Q_{t})$, and $\boldsymbol{p}\in C^{0}([0,T];\boldsymbol{Q})\cap C^{1}((0,T],\boldsymbol{Q})$ such that

[TABLE]

for $i=1,\dots,N$ and such that $\boldsymbol{u}(\cdot,0)=\boldsymbol{u}^{0}(\cdot)$, $p_{t}(\cdot,0)=p_{t}^{0}(\cdot)$, and $p_{i}(\cdot,0)=p_{i}^{0}(\cdot)$ for $i=1,\dots,N$.

4. Partitioned numerical methods with error analysis

In this section, we present two partitioned time discretization algorithms and their a priori error analyses. In Subsection 4.1, we introduce notations and definitions for the algorithms and the associated error analyses. In addition, we prove estimates of initial errors which are commonly necessary for error analyses of the two algorithms. In Subsection 4.2, the first one, called “elasticity-then-diffusion” algorithm, will be defined and the error analysis will be presented with the full details. In Subsection 4.3, we will present another algorithm, so called “diffusion-then-elasticity” algorithm. Since the error analysis of the diffusion-then-elasticity method is similar to the one in Subsection 4.2, we will state main intermediate estimates in the analysis and details will be omitted.

In our error analyses the first method (elasticity-then-diffusion) methods has the second order convergence in time whereas the second method (diffusion-then-elasticity) has the third order convergence in time. In spite of this lower order convergence rate, the first method can be advantageous when local mass conservation is significant because then numerical solutions at the same time step are used in the mass conservation equations, so a locally mass conservative numerical flux can be easily recovered by post-processing. This is not the case in the second method because the numerical solutions of $\boldsymbol{p}$ and $p_{t}$ are at different time steps as we will see in Subsection 4.3.

4.1. Preliminaries

In our numerical algorithms, we take time steps $t_{n}=n\Delta t$ for $n=0,1,\ldots$ and given time step size $\Delta t>0$. For a function $\sigma$ in $C^{0}([0,T];\mathcal{X})$, we use $\sigma^{n}$ to denote $\sigma(t_{n})$. $\sigma_{h}$ is a discrete solution of $\sigma$ if $\sigma$ is a variable in equations and $\sigma_{h}^{n}$ is the $n$-th time step solution of $\sigma_{h}$. For error terms we define $e_{\sigma}^{n}=\sigma^{n}-\sigma_{h}^{n}$. For all quantities with superscript indices, we will use the convention

[TABLE]

For the finite element discretization of $\boldsymbol{u}$ and $p_{t}$ we use a pair of finite element spaces $\boldsymbol{V}_{h}\subset\boldsymbol{V}$, $Q_{t,h}\subset Q_{t}$ which satisfy the inf-sup condition

[TABLE]

with $C>0$ independent of the mesh sizes. For discretization of $p_{i}$’s for $i=1,\ldots,N$, we use a finite element method for the primal form of the Poisson equation yielding a symmetric bilinear form for $\left<K_{i}\nabla p_{i},\nabla q_{i}\right>$. Such methods include the continuous Galerkin (CG) methods, the discontinuous or enriched Galerkin methods (DG or EG) with symmetric bilinear forms (see e.g., [3, 4, 18]), and the finite element space is denoted by $Q_{i,h}$ for $i=1,\ldots,N$. In order to keep this generality, we use $a_{h,i}(\cdot,\cdot)$ to denote the discrete bilinear form corresponding to $\left<K_{i}\nabla p_{i},\nabla q_{i}\right>$. The corresponding discrete norm $\|\cdot\|_{a_{h,i}}$ is defined by

[TABLE]

The convergence orders of numerical solutions depend on the approximation properties of $\boldsymbol{V}_{h}$, $Q_{t,h}$, $Q_{i,h}$ for $i=1,\ldots,N$. For simplicity of presentation we assume that the approximation properties of $\boldsymbol{V}_{h}$ and $Q_{t,h}$ satisfy

[TABLE]

for a positive integer $k_{\boldsymbol{u}}$ if $\boldsymbol{u}$ and $p_{t}$ are sufficiently regular. This assumption holds for the family of Taylor–Hood elements [7, 8, 9, 22] ($\mathcal{P}_{r+1}({\mathbb{R}}^{d})-\mathcal{P}_{r}$, $r\geq 1$) with $k_{\boldsymbol{u}}=r+1$ and for the MINI element [5] with $k_{\boldsymbol{u}}=1$. Similarly, assuming $Q_{1,h}=\cdots=Q_{N,h}$ for simplicity, we assume that

[TABLE]

for $k_{p}\geq 1$ and for sufficiently regular $\boldsymbol{p}$. This holds for CG methods with $Q_{i,h}=\mathcal{P}_{k_{p}}$, or DG/EG methods with piecewise polynomials of degree $k_{p}$. For convenience we define two additional bilinear forms

[TABLE]

with the corresponding norms $\|\cdot\|_{S}$ and $\|\cdot\|_{A}$. By a discrete Poincaré inequality $\|\boldsymbol{q}\|_{S}\lesssim\|\boldsymbol{q}\|_{A}$ holds for all $\boldsymbol{q}\in\boldsymbol{Q}_{h}$ with an implicit constant which is uniform for small $s_{i}$’s and arbitrarily large $\lambda$.

The continuous solutions (with a regularity assumption $p_{i}\in H^{s}(\Omega)$ for $s>3/2$ if DG or EG is used for the discretization of $Q_{i}$) satisfy the variational equations:

[TABLE]

for $\boldsymbol{v}\in\boldsymbol{V}_{h}$, $q_{t}\in Q_{t,h}$, $q_{i}\in Q_{i,h}$. For the error analysis we split error terms into two parts using appropriate interpolation operators. The interpolation operators for the variables $\boldsymbol{u}$, $p_{t}$, $p_{i}$ for $i=1,\ldots,N$ will be denoted by $\Pi_{h}^{\boldsymbol{V}}$, $\Pi_{h}^{Q_{t}}$, $\Pi_{h}^{Q_{i}}$ for $i=1,\ldots,N$. Specifically, we define $(\Pi_{h}^{\boldsymbol{V}}\boldsymbol{u},\Pi_{h}^{Q_{t}}p_{t})$ as the solution of the Lamé equation

[TABLE]

For simplicity, we use $\Pi_{h}^{\boldsymbol{Q}}\boldsymbol{p}$ to denote the $N$-tuple $(\Pi_{h}^{Q_{1}}p_{1},\ldots,\Pi_{h}^{Q_{N}}p_{N})$. Then $\Pi_{h}^{\boldsymbol{Q}}\boldsymbol{p}$ is defined as the solution of the elliptic system

[TABLE]

for $1\leq i\leq N$. Well-posedness of this problem is not difficult to show from the property of $\xi_{i}$

[TABLE]

which, in particular, gives

[TABLE]

By a standard error analysis argument, we can show that these interpolation operators have optimal approximation properties in $H^{1}$ norm for $\boldsymbol{V}_{h}$, in $L^{2}$ norm for $Q_{t,h}$, and in (discrete) $H^{1}$ norm for $Q_{i,h}$, i.e.,

[TABLE]

for sufficiently regular $\boldsymbol{u}$, $p_{t}$, $\boldsymbol{p}$.

For the error analysis we use the notation

[TABLE]

We also define

[TABLE]

for variable $\sigma$.

Throughout the paper we assume that $(\boldsymbol{u}_{h}^{0},p_{t,h}^{0},\boldsymbol{p}_{h}^{0})$, our discrete initial data, satisfies the following:

[TABLE]

As a consequence of (4.10) and the inequality $\|\cdot\|_{S}\lesssim\|\cdot\|_{A}$, we have

[TABLE]

in which the implicit constants of these inequalities depend on the norms of $\boldsymbol{u}^{0}$, $p_{t}^{0}$, $\boldsymbol{p}^{0}$.

To define our partitioned algorithms we need $(\boldsymbol{u}_{h}^{1},p_{t,h}^{1},\boldsymbol{p}_{h}^{1})$ as well. To obtain this, we use the monolithic method combining the backward Euler and Crank–Nicolson methods, i.e., $(\boldsymbol{u}_{h}^{1},p_{t,h}^{1},\boldsymbol{p}_{h}^{1})$ satisfies

[TABLE]

For the error estimates, we need estimates of some error terms in the beginning time steps.

Theorem 1**.**

For given compatible initial data $(\boldsymbol{u}^{0},p_{t}^{0},\boldsymbol{p}^{0})$ and given $\boldsymbol{f}$ and $\{g_{i}\}_{i=1}^{N}$, suppose that $(\boldsymbol{u},p_{t},\boldsymbol{p})$ is the solution of (3.2). For numerical initial data $(\boldsymbol{u}_{h}^{0},p_{t,h}^{0},\boldsymbol{p}_{h}^{0})$ satisfying (4.10) suppose that $(\boldsymbol{u}_{h}^{1},p_{t,h}^{1},\boldsymbol{p}_{h}^{1})$ is obtained by (4.12). Then we have

[TABLE]

and the implicit constants in these inequalities depend on the norms of the exact solutions and are uniformly bounded above for small $s_{i}$’s and arbitrarily large $\lambda$.

Proof.

From $-\operatorname{div}\boldsymbol{u}^{0}-\lambda^{-1}p_{t}^{0}=-\lambda^{-1}\boldsymbol{\alpha}\cdot\boldsymbol{p}^{0}$ we see that

[TABLE]

where $R_{h}^{0}:=\operatorname{div}\boldsymbol{u}_{h}^{0}+\lambda^{-1}p_{t,h}^{0}-\lambda^{-1}\boldsymbol{\alpha}\cdot\boldsymbol{p}_{h}^{0}$. Using this and the interpolation operators defined above we have

[TABLE]

Moreover, the differences of (4.12a), (4.12b) and (4.4a), (4.4b) with $n=1$ with the interpolation operators $\Pi_{h}^{\boldsymbol{V}}$ and $\Pi_{h}^{Q_{t}}$, give

[TABLE]

The differences of these equations are

[TABLE]

The difference of (4.4c) and (4.12c) with $\Delta t$ multiple in consideration of the interpolation $\Pi_{h}^{\boldsymbol{Q}}$ is

[TABLE]

in which $\boldsymbol{I}_{3}^{n}=(\boldsymbol{I}_{3,1}^{n},\ldots,\boldsymbol{I}_{3,N}^{n})$, $I_{4}^{n}$ are defined by

[TABLE]

By the inf-sup condition (4.1), there exists $\boldsymbol{w}^{0}\in\boldsymbol{V}_{h}$ such that

[TABLE]

With a sufficiently small $\delta>0$ independent of the mesh sizes, we can get

[TABLE]

with $C_{0}>0$ independent of the mesh sizes.

If we take $\boldsymbol{v}=D_{\boldsymbol{u}}^{0}+\delta\boldsymbol{w}^{0}$, $q_{t}=-D_{p_{t}}^{0}$, $q_{i}=-D_{p_{i}}^{0}$ in (4.19), (4.20), (4.21), and add them altogether, then we get

[TABLE]

using (4.24). We first remark that

[TABLE]

holds due to the equality $-\operatorname{div}\boldsymbol{u}^{0}-\lambda^{-1}p_{t}^{0}=-\lambda^{-1}\boldsymbol{\alpha}\cdot\boldsymbol{p}^{0}$. Moreover,

[TABLE]

hold. We also remark that

[TABLE]

holds by the triangle inequality $\|\boldsymbol{\alpha}\cdot D_{\boldsymbol{p}}^{0}\|_{\lambda^{-1}}\leq\|D_{p_{t}}^{0}\|_{\lambda^{-1}}+\|D_{p_{t}}^{0}-\boldsymbol{\alpha}\cdot D_{\boldsymbol{p}}^{0}\|_{\lambda^{-1}}$. Then applying Young’s inequality to (4.25), we have

[TABLE]

so (4.13) is proved.

To estimate $\|e_{\boldsymbol{u}}^{h,1}\|_{\boldsymbol{V}}$, $\|e_{p_{t}}^{h,1}\|_{Q_{t}}$, $\|e_{\boldsymbol{p}}^{h,1}\|_{S}$, we need another equation which is the sum of (4.15) and (4.17). If we take $\boldsymbol{v}=e_{\boldsymbol{u}}^{h,0}+e_{\boldsymbol{u}}^{h,1}$ in this equation, $q_{t}=-e_{p_{t}}^{h,0}-e_{p_{t}}^{h,1}$ in (4.19), $q_{i}=-e_{p_{i}}^{h,0}-e_{p_{i}}^{h,1}$ in (4.21) for $1\leq i\leq N$, and add these equations altogether, then we get

[TABLE]

If

[TABLE]

is true, then

[TABLE]

so (4.14) for $\|e_{\boldsymbol{u}}^{h,1}\|_{\boldsymbol{V}}$ follows from (4.11), the Cauchy–Schwarz inquality, and the estimate of $\|D_{p_{t}}^{0}\|_{\lambda^{-1}}$. If (4.28) is not true, then we have

[TABLE]

If we use the Cauchy–Schwarz inquality and divide both sides by

[TABLE]

we can obtain

[TABLE]

Then (4.14) for $\|e_{\boldsymbol{u}}^{h,1}\|_{\boldsymbol{V}}$ is proved. The estimate (4.14) for $\|e_{p_{t}}^{h,1}\|_{Q_{t}}$ can be obtained from (4.17) by taking $\boldsymbol{v}\in\boldsymbol{V}_{h}$ such that $\left<e_{p_{t}}^{h,1},\operatorname{div}\boldsymbol{v}\right>=\|e_{p_{t}}^{h,1}\|_{Q_{t}}^{2}$ and $\|\boldsymbol{v}\|_{\boldsymbol{V}}\lesssim\|e_{p_{t}}^{h,1}\|_{Q_{t}}$. Finally, the same estimate for $\|e_{\boldsymbol{p}}^{h,1}\|_{S}$ follows from the same argument using the triangle inequality as before. ∎

4.2. The first partitioned method (elasticity-then-diffusion)

We here present the first partitioned method inspired by the differential form of (3.2b), i.e.,

[TABLE]

Method 1

Suppose that $(\boldsymbol{u}_{h}^{0},p_{t,h}^{0},\boldsymbol{p}_{h}^{0})$ and $(\boldsymbol{u}_{h}^{1},p_{t,h}^{1},\boldsymbol{p}_{h}^{1})$ are provided by the monolithic numerical method described in the subsection 4.1.

• Step 1

For $n\geq 1$, given $\boldsymbol{u}_{h}^{n}$, $p_{t,h}^{n}$, $\boldsymbol{p}_{h}^{n}$, $\boldsymbol{p}_{h}^{n-1}$, compute $(\boldsymbol{u}_{h}^{n+1},p_{t,h}^{n+1})$ by

[TABLE]

• Step 2

Compute $\boldsymbol{p}_{h}^{n+1}$ with $p_{t,h}^{n}$ and $p_{t,h}^{n+1}$ by

[TABLE]

• Step 3

Repeat Step 1 for $n\leftarrow n+1$

Theorem 2**.**

Let $(\boldsymbol{u},p_{t},\boldsymbol{p})$ be the exact solution of (3.4) for compatible initial data $(\boldsymbol{u}^{0},p_{t}^{0},\boldsymbol{p}^{0})$. For numerical initial data $(\boldsymbol{u}_{h}^{0},p_{t,h}^{0},\boldsymbol{p}_{h}^{0})$ satisfying (4.10), suppose that $\{(\boldsymbol{u}_{h}^{n},p_{t,h}^{n},\boldsymbol{p}_{h}^{n})\}$ is a numerical solution obtained by Method 1. Assuming that the exact solution is sufficiently regular,

[TABLE]

hold with implicit constants depending on the norms of the exact solution.

Proof.

By the triangle inequality and the optimal approximation properties (4.5) and (4.6), it suffices to estimate $\|e_{\boldsymbol{u}}^{h,n}\|_{\boldsymbol{V}}$, $\|e_{p_{t}}^{h,n}\|_{Q_{t}}$, $\|e_{\boldsymbol{p}}^{h,n}\|_{S}$, and $\|e_{\boldsymbol{p}}^{h,n}\|_{A}$. To estimate these terms, we consider the differences of continuous equations and discrete equations. The three continuous equations are (4.4a) with $n=l$, the difference of (4.4b) with $n=l+1$ and $n=l$, and (4.4c) with $n=l$. If we subtract the discrete equations (4.29), (4.30), (4.31) with $n=l$ from these equations, assuming that $\boldsymbol{p}$ is sufficiently regular in case for the DG or EG methods, the differences of the continuous equations and the discrete equations are

[TABLE]

for $l\geq 1$. These equations with the interpolations $\Pi_{h}^{\boldsymbol{V}}$, $\Pi_{h}^{Q_{t}}$, $\Pi_{h}^{\boldsymbol{Q}}$, lead to the error equations

[TABLE]

for $l\geq 1$, $i=1,\ldots,N$ where $\boldsymbol{I}_{j}^{l}=(I_{j,1}^{l},\ldots,I_{j,N}^{l})$ are defined as

[TABLE]

and in (4.22). Considering the differences of two consecutive time steps of (4.32a), we get

[TABLE]

for $l\geq 1$.

We are ready to prove the error estimates. Since the proof is long, we split the proof into two steps. In the first step, we prove

[TABLE]

In the second step, we prove

[TABLE]

It is easy to see that these two estimates complete the proof.

For the first step let $\boldsymbol{w}^{l}$ be an element in $\boldsymbol{V}_{h}$ satisfying the condition (4.23) for $D_{p_{t}}^{l}$. If we take $\boldsymbol{v}=D_{\boldsymbol{u}}^{l}+\delta\boldsymbol{w}^{l}$ in (4.33), $q_{t}=-D_{p_{t}}^{l}$ with $\delta>0$ as above in (4.32b), $q_{i}=-D_{p_{i}}^{l}$ in (4.32c), and add them altogether, then the sum gives

[TABLE]

for $l\geq 1$. By Young’s inequality,

[TABLE]

with any $\epsilon_{1},\epsilon_{2},\epsilon_{3}>0$. From (4.24) and the above two inequalities we can obtain

[TABLE]

We choose sufficiently small $\epsilon_{1}$ and $\epsilon_{2}$ such that $\left(\frac{\epsilon_{1}}{2}+\frac{\epsilon_{2}}{4}\right)\|q_{t}\|_{\lambda^{-1}}^{2}\leq\frac{C_{0}}{2}\|q_{t}\|_{Q_{t}}^{2}$ for all $q_{t}\in Q_{t,h}$ and for the $C_{0}$ in (4.24). Then the above inequality gives

[TABLE]

We take $\epsilon_{3}$ to satisfy $1-\epsilon_{3}-\frac{1}{2+\epsilon_{2}}=\epsilon_{3}$. Taking the summation of the above inequality over the index $l$ from 1 to $n-1$ and ignoring some nonnegative terms, we can get

[TABLE]

Since $\epsilon_{1}$, $\epsilon_{2}$, $\epsilon_{3}$ depend only on $C_{0}$, they are independent of the mesh sizes. A standard argument gives

[TABLE]

From these, (4.26), (4.27), (4.13), and the assumption (4.10), we can obtain

[TABLE]

Now we begin the second step and prove an estimate of $\|e_{\boldsymbol{u}}^{h,n}\|_{\boldsymbol{V}}+\|e_{\boldsymbol{p}}^{h,n}\|_{Q_{t}}+\|e_{\boldsymbol{p}}^{h,n}\|_{S}$. To do it we add the following three equations:

the sum of (4.32a) of indices $l+1$ and $l$ with $\boldsymbol{v}=D_{\boldsymbol{u}}^{l}$

• 2.

(4.32b) with $q_{t}=-2e_{p_{t}}^{h,l+\frac{1}{2}}$

• 3.

the equations (4.32c) with $q_{i}=-2e_{p_{i}}^{h,l+\frac{1}{2}}$ for $1\leq i\leq N$.

From the sum of these equations we get

[TABLE]

Defining $X_{l},Y_{l}\geq 0$ as

[TABLE]

the left-hand side of the above equality is $X_{l+1}^{2}-X_{l}^{2}+Y_{l}^{2}$. Noting that we can obtain $\|e_{p_{t}}^{h,l}\|_{Q_{t}}\lesssim\|e_{\boldsymbol{u}}^{h,l}\|_{\boldsymbol{V}}$ from (4.32a), we can also show that

[TABLE]

by the triangle inequality and the definition of $\|\cdot\|_{S}$ with some $C_{2}>0$ independent of $h$. If we use the formula

[TABLE]

obtained by a summation by parts argument, taking the summation of (4.39) for $l$ over $1$ to $n-1$ gives

[TABLE]

If $X_{n}<\max_{1\leq l\leq n}X_{l}$, then it suffices to estimate $X_{n_{0}}$ for the smallest $1\leq n_{0}<n$ such that $X_{n_{0}}=\max_{1\leq l\leq n_{0}}X_{l}$ because the estimate of $X_{n_{0}}$ ($>X_{n}$) will also give an estimate of $X_{n}$ which can be used to obtain (4.35). Therefore we will show an estimate of $X_{n}$ below with the assumption $X_{n}=\max_{1\leq l\leq n}X_{l}$.

By the Cauchy–Schwarz inequality and (4.42),

[TABLE]

This inequality has a form of $A^{2}+B^{2}\leq CA+D^{2}$ with $A=X_{n}$, $B=\sum_{l=1}^{n-1}Y_{l}$,

[TABLE]

We may assume $A,B>0$ without loss of generality. Then it is easy to show that the above inequality implies

[TABLE]

Therefore either

[TABLE]

or

[TABLE]

holds. Recall that (4.14) gives an estimate of $X_{1}$. Then from the previous estimates (4.36), (4.37), (4.26), (4.27), and (4.38), we can conclude that

[TABLE]

This proves (4.35) for $\|e_{\boldsymbol{u}}^{h,n}\|_{\boldsymbol{V}}$. The estimate for other two terms easily follows because $\|e_{\boldsymbol{p}}^{h,n}\|_{S}\leq C_{2}X_{n}$ and $\|e_{p_{t}}^{h,n}\|_{Q_{t}}\lesssim\|e_{\boldsymbol{u}}^{h,n}\|_{\boldsymbol{V}}$ hold.

∎

4.3. The second partitioned algorithm (diffusion-then-elasticity)

We now present the second partitioned method.

Method 2

Suppose that $(\boldsymbol{u}_{h}^{0},p_{t,h}^{0},\boldsymbol{p}_{h}^{0})$ and $(\boldsymbol{u}_{h}^{1},p_{t,h}^{1},\boldsymbol{p}_{h}^{1})$ are provided by the monolithic numerical method described in the subsection 4.1.

• Step 1

For given $p_{t,h}^{n}$, $p_{t,h}^{n-1}$, compute $\boldsymbol{p}_{h}^{n+1}$ from

[TABLE]

• Step 2

For given $\boldsymbol{p}_{h}^{n+1}$, compute $(\boldsymbol{u}_{h}^{n+1},p_{t,h}^{n+1})$ with

[TABLE]

for any $\boldsymbol{v}\in\boldsymbol{V}$ and $q_{t}\in Q_{t,h}$

• Step 3

Repeat Step 1 with $n\leftarrow n+1$

Theorem 3**.**

Let $(\boldsymbol{u},p_{t},\boldsymbol{p})$ be an exact solution of (3.4) for compatible initial data $(\boldsymbol{u}^{0},p_{t}^{0},\boldsymbol{p}^{0})$. For numerical initial data $(\boldsymbol{u}_{h}^{0},p_{t,h}^{0},\boldsymbol{p}_{h}^{0})$ satisfying (4.10), suppose that $\{(\boldsymbol{u}_{h}^{n},p_{t,h}^{n},\boldsymbol{p}_{h}^{n})\}$ is a numerical solution obtained by Method 2. Assuming that the exact solutions are sufficiently regular, the following hold

[TABLE]

Remark 1**.**

Compared to the elasticity-then-diffusion method, this method has one higher order convergence of time discretization errors.

Proof.

Since the proof is similar to the one of Theorem 2, we sketch the proof without full details. As before, interpolation errors denoted by $e_{\sigma}^{I,n}$ for unknown $\sigma$ are of optimal order, so it is enough to show the estimates for the errors of $e_{\sigma}^{h,n}$.

The differences of (4.4) and (4.44), (4.45), (4.43) with index $l$ are

[TABLE]

By the interpolation operators $\Pi_{h}^{\boldsymbol{V}}$, $\Pi_{h}^{Q_{t}}$, $\Pi_{h}^{\boldsymbol{Q}}$, and by (4.7), we can obtain error equations from (4.46) which are

[TABLE]

for $l\geq 1$ where

[TABLE]

For the difference of the equations (4.47a) and (4.47b) with indices $l+1$ and $l$, we take $\boldsymbol{v}=D_{\boldsymbol{u}}^{l}+\delta\boldsymbol{w}^{l}$, $q_{t}=-D_{p_{t}}^{l}$, and take $q_{i}=-D_{p_{i}}^{l}$ in (4.47c). The sum of these equations yield

[TABLE]

We can use an argument completely analogous to the proof of (4.34) and get

[TABLE]

with implicit constants independent of $h$. This gives the assertion on $\|\boldsymbol{p}^{n}-\boldsymbol{p}_{h}^{n}\|_{A}$.

We take $\boldsymbol{v}=D_{\boldsymbol{u}}^{l-1}$ for the sum of (4.47a) with indices $l+1$, $l$, and $q_{t}=-2e_{p_{t}}^{h,l-\frac{1}{2}}$ in the difference of the equation (4.47b) with indices $l+1$, $l$, and $q_{i}=-2e_{p_{i}}^{h,l+\frac{1}{2}}$ in (4.47c). If we add these equations altogether, then we get

[TABLE]

With $X_{l}$, $Y_{l}$ in (4.40), (4.41), we have

[TABLE]

Using

[TABLE]

by a similar argument for the estimate of $X_{n}$ in the proof of Theorem 2, we can estimate ${X}_{n}+\left(\sum_{l=1}^{n-1}Y_{l}^{2}\right)^{\frac{1}{2}}$ by a linear combination of $X_{1}$,

[TABLE]

These terms are bounded by $(\Delta t)^{3}+h^{k_{\boldsymbol{u}}}+h^{k_{p}}$ by the estimates of interpolation error terms, (4.48), and Theorem 1. As in the proof of Theorem 2 we can prove the asserted error estimates from the estimate of $X_{n}$. ∎

5. Numerical convergence experiments

In this section, we present a set of numerical results to illustrate the performances of the new partitioned methods. For this we construct a manufactured smooth solution and compute convergence of errors of the numerical solutions. All numerical simulations in this section were run using the FEniCS finite element software [2] (version 2017.2.0). For simplicity we only consider the two-network case, i.e., $N=2$ in our numerical experiments. The exact solutions in our experiments are

[TABLE]

The physical parameters are given as $E=1.0$, $\nu=0.49999$, $\alpha_{1}=\alpha_{2}=1$, $s_{1}=s_{2}=1$, $K_{1}=K_{2}=1$, and corresponding $\mu$ and $\lambda$ are computed by $\mu=E/(2(1+\nu))$ and $\lambda=\nu E/((1-2\nu)(1+\nu))$. The domain $\Omega$ is $[0,1]\times[0,1]$. To check convergence of errors we use nested structured triangular meshes obtained by dividing $\Omega$ into $M\times M$ rectangles ($M=8,16,32,64,128$), i.e., the relation $h\sim 1/M$ holds with uniform implicit constants. For simplicity we impose Dirichlet boundary conditions of $\boldsymbol{u}$ on the vertical sides of $\Omega$ whereas we impose Dirichlet boundary conditions of $p_{1}$ and $p_{2}$ on the whole boundary of $\Omega$. The time step size $\Delta t$ is taken as $\Delta t=1/M$ in all experiments, so $\Delta t\sim h$ holds. In all the numerical experiments we use the Taylor–Hood elements with polynomial degrees of $k+1$ and $k$ for $\boldsymbol{V}_{h}$ and $Q_{t,h}$, and the Lagrange finite elements of degree $k$ for $Q_{1,h}$ and $Q_{2,h}$.

In Table 1 and Table 2, we present the errors of variables computed at $t=1$ with $k=1$ for the diffusion-then-elasticity (DTE) and the elasticity-then-diffusion (ETD) schemes, respectively. The convergence rates of the errors of $p_{1}$ and $p_{2}$ in the $H^{1}$ norm are bounded by 1 because of the best approximation property of linear polynomials. Other errors show convergence rates higher than the rates expected by our error analysis.

To confirm the second order convergence in time, we test the numerical methods with $k=2$, and the results for DTE and ETD schemes are given in Table 3 and Table 4. In the results we observe that some errors show the convergence rates which are higher than the rates expected by the error analysis. The expected convergence rates of $\|p_{1}-p_{1,h}\|_{1}$ and $\|p_{2}-p_{2,h}\|_{1}$ are $\frac{3}{2}$ in the error analysis but the numerical results show second order convergence. We conjecture that there is an improved way to analyze the errors with second order convergence of these errors but we leave it as a future research work. The convergence rates of the errors $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}$ and $\|p_{t}-p_{t,h}\|_{0}$ are also not covered by the present error analysis. Note that these errors have $\frac{1}{2}$ higher order convergence rates in the current error analysis and note also that some time error terms are of orders higher than 2 (see (4.26) and (4.27)). Thus we believe that the improved error analysis in our conjecture can be used to obtain these higher convergence rates of $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}$ and $\|p_{t}-p_{t,h}\|_{0}$.

Finally, we run experiments with $k=3$ and the results are presented in Table 5 and Table 6. Interestingly, the convergence rates of $\|p_{1}-p_{1,h}\|_{1}$ and $\|p_{2}-p_{2,h}\|_{1}$ are approximately $\frac{5}{2}$ but those of $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}$ and $\|p_{t}-p_{t,h}\|_{0}$ are asymptotically 4. The analysis in Subsection 4.3 explains the $\frac{5}{2}$ convergence of $\|p_{1}-p_{1,h}\|_{1}$ and $\|p_{2}-p_{2,h}\|_{1}$ in the DTE scheme but it still cannot explain the superconvergence of $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}$ and $\|p_{t}-p_{t,h}\|_{0}$ of order 4. In the ETD scheme, the superconvergent errors are observed and they are beyond the scope of the analysis in Subsection 4.2. At present we have no theory to explain this superconvergence, so it will be investigated in our future research.

6. Conclusions

In this paper, we have presented two partitioned time discretization schemes for the total-pressure-based formulation of quasi-static multiple-network poroelasticity. We proved that the partitioned schemes are unconditionally stable and have second and third order convergence in time by a novel error analysis. The analyses also show that the numerical schemes are still robust in the limit of incompressibility and other parameter variations such as small storage coefficients. We also presented a number of numerical experiments to illustrate the validaty of our theoretical analysis but the numerical results often show superconvergence results in time discretization errors which are not completely understood by the current error analysis. Thus, an improved error analysis for the schemes will be investigated in the future work.