Computational Multiscale Methods for Linear Poroelasticity with High Contrast

TL;DR

This paper presents a multiscale finite element method for efficiently solving high-contrast linear poroelasticity problems, demonstrating convergence and improved computational performance through numerical tests.

Contribution

It introduces a novel CEM-GMsFEM approach with energy minimization and oversampling for high-contrast poroelasticity, enhancing efficiency and accuracy.

Findings

Method achieves first-order convergence.

Numerical tests confirm computational efficiency.

Effective handling of high contrast in coefficients.

Abstract

In this work, we employ the Constraint Energy Minimizing Generalized Multiscale Finite Element Method (CEM-GMsFEM) to solve the problem of linear heterogeneous poroelasticity with coefficients of high contrast. The proposed method makes use of the idea of energy minimization with suitable constraints in order to generate efficient basis functions for the displacement and the pressure. These basis functions are constructed by solving a class of local auxiliary optimization problems based on eigenfunctions containing local information on the heterogeneity. Techniques of oversampling are adapted to enhance the computational performance. Convergence of first order is shown and illustrated by a number of numerical tests.

Tables

Table 1. Table 1: Numerical results with varying coarse grid size H 𝐻 H for test model 1 1 1 .

$J$	$H$	$m$	$e_{L^2}^u$	$e_a^u$	$e_{L^2}^p$	$e_b^p$
4	$\sqrt{2}$/10	4	9.41e-03	1.14e-01	6.05e-03	5.79e-02
4	$\sqrt{2}$/20	5	1.22e-03	7.39e-02	8.75e-04	2.29e-02
4	$\sqrt{2}$/40	6	2.08e-04	2.08e-02	1.58e-04	9.64e-03

Table 2. Table 2: Numerical results with varying coarse grid size H 𝐻 H for test model 2.

$J$	$H$	$m$	$e_{L^2}^u$	$e_a^u$	$e_{L^2}^p$	$e_b^p$
4	$\sqrt{2}$/10	4	2.22e-02	5.14e-01	9.64e-05	3.59e-02
4	$\sqrt{2}$/20	5	3.95e-03	2.06e-01	2.77e-05	1.49e-02
4	$\sqrt{2}$/40	6	4.94e-04	5.60e-02	7.81e-06	4.50e-03

Equations

$$-\nabla \cdot \sigma (u)+\nabla (\alpha p)$$

$$\displaystyle\partial_{t}\bigg{(}\alpha\nabla\cdot u+\frac{1}{M}p\bigg{)}-\nabla\cdot\bigg{(}\frac{\kappa}{\nu}\nabla p\bigg{)}$$

$$u$$

$$p$$

$$p(\cdot ,0)$$

$${\sigma(u)}=2\mu{\epsilon(u)}+\lambda(\nabla\cdot{u})\,{\mathcal{I}}\quad\text{and}\quad{\epsilon}({u})=\frac{1}{2}\Big{(}\nabla{u}+(\nabla{u})^{T}\Big{)},$$

$$\lambda =\frac{{\nu }_p}{(1-2{\nu }_p)(1+{\nu }_p)}E,\mu =\frac{1}{2(1+{\nu }_p)}E.$$

$$\|v\|_{L^{r}(0,T;X)}:=\bigg{(}\int_{0}^{T}\|v\|_{X}^{r}\,\text{d}t\bigg{)}^{1/r},\quad 1\leq r<\infty,$$

$$\Vert v{\Vert }_{L^{\infty }(0,T;X)}:=\underset{0\le t\le T}{\sup }\Vert v{\Vert }_X,$$

$$V_0:=[H_0^1(\Omega ){]}^d,Q_0:=H_0^1(\Omega ).$$

$$a(u,v)-d(v,p)$$

$$d({\partial }_{t}u,q)+c({\partial }_{t}p,q)+b(p,q)$$

$$p(\cdot ,0)=p^0\in Q_0.$$

$$a(u,v)={\int }_{\Omega }\sigma (u):ϵ(v)\text{d}x,b(p,q)={\int }_{\Omega }\frac{\kappa }{\nu }\nabla p\cdot \nabla q\text{d}x,$$

$$c(p,q)={\int }_{\Omega }\frac{1}{M}pq\text{d}x,d(u,q)={\int }_{\Omega }\alpha (\nabla \cdot u)q\text{d}x.$$

$$c_{\sigma }\Vert v{\Vert }_1^2\le a(v,v)=:\Vert v{\Vert }_a^2\le C_{\sigma }\Vert v{\Vert }_1^2$$

$$c_{\kappa }\Vert q{\Vert }_1^2\le b(q,q)=:\Vert q{\Vert }_b^2\le C_{\kappa }\Vert q{\Vert }_1^2$$

$$V_h:=\{v\in V_0:v{|}_K\text{ is a polynomial of degree }\le 1\text{ for all }K\in {\mathcal{T}}^h\},$$

$$Q_h:=\{q\in Q_0:q{|}_K\text{ is a polynomial of degree }\le 1\text{ for all }K\in {\mathcal{T}}^h\}.$$

$$a(u_h^n,v)-d(v,p_h^n)$$

$$d(D_{\tau }{u}_h^n,q)+c(D_{\tau }{p}_h^n,q)+b(p_h^n,q)$$

$$a(u_h^0,v)=d(v,p_h^0)$$

$$a_i(v_j^i,v)={\lambda }_j^i{s}_i^1(v_j^i,v)$$

$$b_i(q_j^i,q)={\zeta }_j^i{s}_i^2(q_j^i,q)$$

$$\tilde{\sigma }:=\sum _{i=1}^{N_{\mathrm{v}}}(\lambda +2\mu )|\nabla {\chi }_i^1{|}^2,\tilde{\kappa }:=\sum _{i=1}^{N_{\mathrm{v}}}\frac{\kappa }{\nu }|\nabla {\chi }_i^2{|}^2.$$

$$V_{\text{aux}}:=\underset{i=1}{\overset{N_{\mathrm{e}}}{⨁}}{V}_{\text{aux}}(K_i)\subseteq V_0\text{and}{Q}_{\text{aux}}:=\underset{i=1}{\overset{N_{\mathrm{e}}}{⨁}}{Q}_{\text{aux}}(K_i)\subseteq Q_0.$$

$$s^1(u,v):=\sum _{i=1}^{N_{\mathrm{e}}}{s}_i^1(u,v),s^2(p,q):=\sum _{i=1}^{N_{\mathrm{e}}}{s}_i^2(p,q)$$

$${\pi }^1(v):=\sum _{i=1}^{N_{\mathrm{e}}}\sum _{j=1}^{J_i^1}{s}_i^1(v,v_j^i)v_j^i,{\pi }^2(q):=\sum _{i=1}^{N_{\mathrm{e}}}\sum _{j=1}^{J_i^2}{s}_i^2(q,q_j^i)q_j^i.$$

$$K_{i,0}:=K_i,K_{i,m}:=\bigcup \left\{K\in {\mathcal{T}}^H:K\cap \overline{K_{i,m-1}}\ne \varnothing \right\},m=1,2,\dots ,$$

$$\psi_{j,m}^{i}=\text{argmin}\Big{\{}a(\psi,\psi)+s^{1}\big{(}\pi^{1}(\psi)-v_{j}^{i},\pi^{1}(\psi)-v_{j}^{i}\big{)}:\,\psi\in V_{0}(K_{i,m})\Big{\}}$$

$$\phi_{j,m}^{i}=\text{argmin}\Big{\{}b(\phi,\phi)+s^{2}\big{(}\pi^{2}(\phi)-q_{j}^{i},\pi^{2}(\phi)-q_{j}^{i}\big{)}:\,\phi\in Q_{0}(K_{i,m})\Big{\}}.$$

$$a(\psi_{j,m}^{i},v)+s^{1}\big{(}\pi^{1}(\psi_{j,m}^{i}),\pi^{1}(v)\big{)}=s^{1}\big{(}v_{j}^{i},\pi^{1}(v)\big{)}$$

$$b(\phi_{j,m}^{i},q)+s^{2}\big{(}\pi^{2}(\phi_{j,m}^{i}),\pi^{2}(v)\big{)}=s^{2}\big{(}q_{j}^{i},\pi^{2}(q)\big{)}$$

$$V_{\text{ms}}:=\text{span}\{{\psi }_{j,m}^i:1\le j\le J_i^1,1\le i\le N_{\mathrm{e}}\}\text{and}{Q}_{\text{ms}}:=\text{span}\{{\varphi }_{j,m}^i:1\le j\le J_i^2,1\le i\le N_{\mathrm{e}}\},$$

Full text

Keywords: Linear poroelasticity, High contrast values, Generalized multiscale finite element method, Constraint energy minimization.

1 Introduction

Modeling and simulating the deformation of porous media saturated by an incompressible viscous fluid is important for a wide range of applications such as reservoir engineering in the field of geomechanics [30] as well as environmental safety due to overburden subsidence and compaction [27]. A reasonable model should couple the flow of the fluid with the behavior of the surrounding solid. Biot [4] proposed a poroelasticity model that couples a Darcy flow with the linear elastic behavior of the porous medium. This model consists of coupled equations for the pressure and the displacements. The stress equation represents quasi-static elasticity coupled to the pressure gradients as a forcing term. On the other hand, the pressure equation is a Darcy-type parabolic equation with a time-dependent coupling to the volumetric strain.

Standard numerical methods such as the finite element method can be used to solve the poroelastic system in case that the medium is homogeneous [14], i.e., present parameters are constant. If the medium is strongly heterogeneous, however, the discretization of the domain needs to be sufficiently fine to obtain accurate results. Such a method is costly in the sense that the dimension of the resulting linear system is huge and therefore not feasible for practical computations.

To alleviate the computational burden, model reduction techniques such as upscaling and multiscale methods can be applied. In upscaling methods [10, 15, 24, 26, 29], one typically upscales the properties of the medium based on the theory of homogenization. The resulting systems can be solved on a coarse grid with standard techniques and the dimensions of the corresponding linear systems are much smaller. In multiscale methods [2, 6, 12, 18, 19], however, one still solves the problem on a coarse grid but with precomputed multiscale basis functions using local information of the medium. For the problem of linear poroelasticity in highly heterogeneous media, this was done, e.g., in  [5]. Therein, a set of multiscale basis functions is constructed under the framework of the Generalized Multiscale Finite Element Method (GMsFEM), see [7, 11]. Another technique to solve multiscale problems is the Localized Orthogonal Decomposition method (LOD), cf. [13, 16, 17, 23, 25]. This method was applied in [22] to linear heterogeneous thermoelasticity and optimal first-order convergence of the fully discretized system based on LOD and an implicit Euler discretization in time was proven. This approach was transferred to the present poroelastic setting in [1].

In the present work, we combine the ideas from LOD and GMsFEM as recently proposed in [8]. For this, the basis functions are constructed by the principle of constraint energy minimization. Based on the GMsFEM, we first create so-called snapshot spaces and afterwards perform model reduction within those spaces by locally solving a class of well-designed spectral problems. The convergence of this Constraint Energy Minimizing GMsFEM (CEM-GMsFEM) is analyzed in [8]. Therein, it is shown that this method has a convergence rate proportional to the coarse grid size, which remains valid even in the presence of high contrast provided that sufficiently many basis functions are selected. The approach makes use of the ideas of localization [20, 21, 23, 25] and oversampling to compute multiscale basis functions in some oversampled subregions with the aim to obtain an appropriate orthogonality condition.

Adopting the idea of CEM-GMsFEM, we propose a multiscale method for the problem of linear poroelasticity with high contrast and construct multiscale spaces for both, the pressure and the displacement. Based on the previous work [1], we prove the first-order convergence of the implicit Euler scheme combined with CEM-GMsFEM for the spatial discretization. Numerical results are provided to demonstrate the efficiency of the proposed method. While in the LOD approach the number of basis functions is limited by the number of coarse nodes, the present CEM-GMsFEM setting allows to add additional basis functions in a flexible way based on spectral properties of the differential operators. This improves the accuracy of the method in the presence of high contrast. It is shown that if enough basis functions are selected, the convergence of the method can be shown independently of the contrast. Unfortunately, a high number of basis functions directly influences the computational complexity of the method. The direct influence of the contrast on the needed number of basis functions is not known but numerical results indicate that a moderate number of basis functions, depending logarithmically on the contrast, seems sufficient.

The paper is organized as follows. In Section 2, we introduce the model problem and the functional analytical setting. The framework of CEM-GMsFEM including the construction of the basis functions and the resulting fully discrete method are presented in Section 3. In Section 4, we analyze the method and provide the corresponding convergence results. Numerical experiments, proving the expected performance of the proposed method, are shown in Section 5.

2 Preliminaries

2.1 Model problem

Let $\Omega\subset\mathbb{R}^{d}$ ($d=2,3$) be a bounded and polyhedral Lipschitz domain and $T>0$ a fixed time. We consider the problem of linear poroelasticity where we are interested in finding the pressure $p\colon[0,T]\times\Omega\to\mathbb{R}$ and the displacement field $u\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ satisfying

[TABLE]

with boundary and initial conditions

[TABLE]

For the sake of simplicity, we only consider homogeneous Dirichlet boundary here. The extension to other types of boundary conditions is straightforward. In this model, the primary sources of the heterogeneity are the stress tensor $\sigma$, the permeability $\kappa$, and the Biot-Willis fluid-solid coupling coefficient $\alpha$. We denote by $M$ the Biot modulus and by $\nu$ the fluid viscosity. Both are assumed to be constant. Moreover, $f$ is a source term representing injection or production processes. Body forces, such as gravity, are neglected. In the case of a linear elastic stress-strain constitutive relation, the stress and strain tensors may be expressed as

[TABLE]

where $\mathcal{I}$ is the identity tensor and $\lambda,\,\mu>0$ are the Lamé coefficients, which can also be expressed in terms of the Young’s modulus $E>0$ and the Poisson ratio $\nu_{p}\in(-1,1/2)$,

[TABLE]

In the considered case of heterogeneous media, the coefficients $\mu$, $\lambda$, $\kappa$, and $\alpha$ may be highly oscillatory.

2.2 Function spaces

In this subsection, we clarify the notation used throughout the article. We write $(\cdot,\cdot)$ to denote the inner product in $L^{2}(\Omega)$ and $\|\cdot\|$ for the corresponding norm. Let $H^{1}(\Omega)$ be the classical Sobolev space with norm \|v\|_{1}:=\big{(}\|v\|^{2}+\|\nabla v\|^{2}\big{)}^{1/2} and $H_{0}^{1}(\Omega)$ the subspace of functions having a vanishing trace. We denote the corresponding dual space by $H^{-1}(\Omega)$. Moreover, we write $L^{r}(0,T;X)$ for the Bochner space with the norm

[TABLE]

[TABLE]

where $(X,\|\cdot\|_{X})$ is a Banach space. Also, we define $H^{1}(0,T;X):=\{v\in L^{2}(0,T;X):\partial_{t}v\in L^{2}(0,T;X)\}$. To shorten notation, we define the spaces for the displacement $u$ and the pressure $p$ by

[TABLE]

2.3 Variational formulation and discretization

In this subsection, we provide the variational formulation corresponding to the system (1). We first multiply the equations (1a) and (1b) with test functions from $V_{0}$ and $Q_{0}$, respectively. Then, applying Green’s formula and making use of the boundary conditions (2a) and (2b), we obtain the following variational problem: find $u(\cdot,t)\in V_{0}$ and $p(\cdot,t)\in Q_{0}$ such that

[TABLE]

The bilinear forms are defined by

[TABLE]

Note that (3a) can be used to define a consistent initial value $u^{0}:=u(\cdot,0)\in V_{0}$. Using Korn’s inequality [9], we get

[TABLE]

for all $v\in V_{0}$, where $c_{\sigma}$ and $C_{\sigma}$ are positive constants. Similarly, there exist two positive constants $c_{\kappa}$ and $C_{\kappa}$ such that

[TABLE]

for all $q\in Q_{0}$. The proof of existence and uniqueness of solutions $u$ and $p$ to (3) can be found in [28].

To discretize the variational problem (3), let $\mathcal{T}^{h}$ be a conforming partition for the computational domain $\Omega$ with (local) grid sizes $h_{K}:=\text{diam}(K)$ for $K\in\mathcal{T}^{h}$ and $h:=\text{max}_{K\in\mathcal{T}^{h}}h_{K}$. We remark that $\mathcal{T}^{h}$ is referred to as the fine grid. Next, let $V_{h}$ and $Q_{h}$ be the standard finite element spaces of first order with respect to the fine grid $\mathcal{T}^{h}$, i.e.,

[TABLE]

[TABLE]

For the time discretization, let $\tau$ be a uniform time step and define $t_{n}=n\tau$ for $n=0,1,\dots,N$ and $T=N\tau$. We use the backward Euler method, i.e., for $n=1,\dots,N$ and given $p_{h}^{0},\ u_{h}^{0}$, we aim to find $u_{h}^{n}\in V_{h}$ and $p_{h}^{n}\in Q_{h}$ such that

[TABLE]

for all $v\in V_{h}$ and $q\in Q_{h}$. Here, $D_{\tau}$ denotes the discrete time derivative, i.e., $D_{\tau}u_{h}^{n}:=(u_{h}^{n}-u_{h}^{n-1})/\tau$, and $f^{n}:=f(t_{n})$. The initial value $p_{h}^{0}\in Q_{h}$ is set to be the $L^{2}$ projection of $p^{0}\in Q_{0}$. The initial value $u_{h}^{0}$ for the displacement can be obtained by solving

[TABLE]

for all $v\in V_{h}$. We remark that this classical approach will serve as a reference solution. The aim of this research is to construct a reduced system based on (4). To this end, we introduce finite-dimensional multiscale spaces $V_{\text{ms}}\subseteq V_{0}$ and $Q_{\text{ms}}\subseteq Q_{0}$, whose dimensions are much smaller, for approximating the solution on some feasible coarse grid.

3 Construction of the multiscale spaces

In this section, we construct multiscale spaces on a coarse grid.

Let $\mathcal{T}^{H}$ be a conforming partition of the computational domain $\Omega$ such that $\mathcal{T}^{h}$ is a refinement of $\mathcal{T}^{H}$. We call $\mathcal{T}^{H}$ the coarse grid and each element of $\mathcal{T}^{H}$ is a coarse block. We denote with $H:=\text{max}_{K\in\mathcal{T}^{H}}\text{diam}(K)$ the coarse grid size. Let $N_{\mathrm{v}}$ be the total number of (interior) vertices of $\mathcal{T}^{H}$ and $N_{\mathrm{e}}$ be the total number of coarse elements. Let $\{x_{i}\}_{i=1}^{N_{\mathrm{v}}}$ be the set of nodes in $\mathcal{T}^{H}$. Figure 1 illustrates the fine grid and a coarse element $K_{i}$, and the oversampling domain $K_{i,1}$. The construction of the multiscale spaces consists of two steps. The first step is to construct auxiliary multiscale spaces using the concept of GMsFEM. Based on the auxiliary spaces, we can then construct multiscale spaces containing basis functions whose energy are minimized in some subregions of the domain. These energy-minimized basis functions can then be used to construct a multiscale solution.

3.1 Auxiliary spaces

Here, we construct auxiliary multiscale basis functions by solving spectral problems on each coarse element $K_{i}$ using the spaces $V(K_{i}):={V_{0}|}_{K_{i}}$ and $Q(K_{i}):={Q_{0}|}_{K_{i}}$. More precisely, we consider the local eigenvalue problems: find $(\lambda_{j}^{i},v_{j}^{i})\in\mathbb{R}\times V(K_{i})$ such that

[TABLE]

for all $v\in V(K_{i})$ and find $(\zeta_{j}^{i},q_{j}^{i})\in\mathbb{R}\times Q(K_{i})$ such that

[TABLE]

for all $q\in Q(K_{i})$, where $a_{i}(u,v):=\int_{K_{i}}\sigma(u):\epsilon(v)\,\text{d}x$, $b_{i}(p,q):=\int_{K_{i}}\frac{\kappa}{\nu}\nabla p\cdot\nabla q\,\text{d}x$, $s^{1}_{i}(u,v):=\int_{K_{i}}\tilde{\sigma}u\cdot v\,\text{d}x$, and $s^{2}_{i}(p,q):=\int_{K_{i}}\tilde{\kappa}^{2}pq\,\text{d}x$ with

[TABLE]

The functions ${\chi_{i}^{1}}$ and ${\chi_{i}^{2}}$ are neighborhood-wise defined partition of unity functions [3] on the coarse grid. To be more precise, for $k=1,2$ the function $\chi_{i}^{k}$ satisfies $H\lvert\nabla\chi_{i}^{k}\rvert=O(1)$, $0\leq\chi_{i}^{k}\leq 1$, and $\sum_{i=1}^{N_{\mathrm{v}}}\chi_{i}^{k}=1$. Assume that the eigenvalues $\{\lambda_{j}^{i}\}$ (resp. $\{\zeta_{j}^{i}\}$) are arranged in ascending order and that the eigenfunctions satisfy the normalization condition $s_{i}^{1}(v_{j}^{i},v_{j}^{i})=1$ as well as $s_{i}^{2}(q_{j}^{i},q_{j}^{i})=1$. Next, choose $J_{i}^{1}\in\mathbb{N}^{+}$ and define the local auxiliary space $V_{\text{aux}}(K_{i}):=\text{span}\{v_{j}^{i}:1\leq j\leq J_{i}^{1}\}$. Similarly, we choose $J_{i}^{2}\in\mathbb{N}^{+}$ and define $Q_{\text{aux}}(K_{i}):=\text{span}\{q_{j}^{i}:1\leq j\leq J_{i}^{2}\}$. Based on these local spaces, we define the global auxiliary spaces $V_{\text{aux}}$ and $Q_{\text{aux}}$ by

[TABLE]

The corresponding inner products for the global auxiliary multiscale spaces are defined by

[TABLE]

for all $u,v\in V_{\text{aux}}$ and $p,q\in Q_{\text{aux}}$. Further, we define projection operators $\pi^{1}:V_{0}\to V_{\text{aux}}$ and $\pi^{2}:Q_{0}\to Q_{\text{aux}}$ such that for all $v\in V_{0},\ q\in Q_{0}$ it holds that

[TABLE]

3.2 Multiscale spaces

In this subsection, we construct the multiscale spaces for the practical computations. For each coarse element $K_{i}$, we define the oversampled region $K_{i,m}\subset\Omega$ obtained by enlarging $K_{i}$ by $m$ layers, i.e.,

[TABLE]

see Figure 1 for an illustration of $K_{i,1}$. We define $V_{0}(K_{i,m}):=[H_{0}^{1}(K_{i,m})]^{d}$ and $Q_{0}(K_{i,m}):=H^{1}_{0}(K_{i,m})$. Then, for each pair of auxiliary functions $v_{j}^{i}\in V_{\text{aux}}$ and $q_{j}^{i}\in Q_{\text{aux}}$, we solve the following minimization problems: find $\psi_{j,m}^{i}\in V_{0}(K_{i,m})$ such that

[TABLE]

and find $\phi_{j,m}^{i}\in Q_{0}(K_{i,m})$ such that

[TABLE]

Note that problem (8) is equivalent to the local problem

[TABLE]

for all $v\in V_{0}(K_{i,m})$, whereas problem (9) is equivalent to

[TABLE]

for all $q\in Q_{0}(K_{i,m})$. Finally, for fixed parameters $m$, $J_{i}^{1}$, and $J_{i}^{2}$, the multiscale spaces $V_{\text{ms}}$ and $Q_{\text{ms}}$ are defined by

[TABLE]

see also Figure 2 for an illustration of such a multiscale basis function.

The multiscale basis functions can be interpreted as approximations to global multiscale basis functions $\psi_{j}^{i}\in V_{0}$ and $\phi_{j}^{i}\in Q_{0}$, similarly defined by

[TABLE]

These basis functions have global support in the domain $\Omega$ but, as shown in [8], decay exponentially outside some local (oversampled) region. This property plays a vital role in the convergence analysis of CEM-GMsFEM and justifies the use of local basis functions in $V_{\text{ms}}$ and $Q_{\text{ms}}$.

3.3 The multiscale method

In order to make the multiscale spaces $V_{\text{ms}}$ and $Q_{\text{ms}}$ suitable for computations, we need finite-dimensional analogons. For this, we follow the construction of the previous subsections, restricted the finite element space based on the fine grid $\mathcal{T}^{h}$. This then yields the following fully discrete scheme: for $n=1,2,\dots,N$ find $(u_{\text{ms}}^{n},p_{\text{ms}}^{n})\in V_{\text{ms}}\times Q_{\text{ms}}$ that solve

[TABLE]

for all $(v,q)\in V_{\text{ms}}\times Q_{\text{ms}}$ with initial condition $p_{\text{ms}}^{0}\in Q_{\text{ms}}$ defined by

[TABLE]

for all $q\in Q_{\text{ms}}$.

4 Convergence analysis

In this section, we analyze the proposed multiscale method (10). First, we recall some theoretical results related to the discretization of the problem of linear poroelasticity with finite elements. Throughout this section, $C$ denotes a generic constant which is independent of spatial discretization parameters and the time step size. Further, the notation $a\lesssim b$ is used equivalently to $a\leq Cb$.

Lemma 4.1** ([1, Lem. 3.1]).**

Given initial data $p_{h}^{0}\in Q_{h}$ and $u_{h}^{0}\in V_{h}$ defined in (5), system (4) is well-posed. That is, there exists a unique solution, which can be bounded in terms of the initial values and the source function.

Lemma 4.2** ([22, Thm. 3.3]).**

Assume $f\in L^{\infty}(0,T;L^{2}(\Omega))\cap H^{1}(0,T;H^{-1}(\Omega))$. Then, for all $n=1,\dots,N$, the fully discrete solution $(u_{h}^{n},p_{h}^{n})$ of (4) satisfies the stability bound

[TABLE]

Further, if $p_{h}^{0}=0$, we have

[TABLE]

and for $f=0$ it holds that

[TABLE]

Lemma 4.3** (cf. [14, Thm. 3.1]).**

Assume that the coefficients satisfy $\mu,\lambda,\kappa,\alpha\in W^{1,\infty}(\Omega)$. Further, let the solution $(u,p)$ of (3) be sufficiently smooth. Then, for each $n=1,\dots,N$ the fully discrete solution $(u_{h}^{n},p_{h}^{n})$ of (4) satisfies the estimate

[TABLE]

where the constant $C_{\text{osc}}$ scales with $\max\{\|\mu\|_{W^{1,\infty}(\Omega)},\|\lambda\|_{W^{1,\infty}(\Omega)},\|\kappa\|_{W^{1,\infty}(\Omega)},\|\alpha\|_{W^{1,\infty}(\Omega)}\}$.

Remark**.**

The previous lemma states the first-order convergence of the finite element method but also reveals the dependence of the involved constant on possible oscillations in the coefficients describing the media.

Next, for any $u\in V_{h}$, $p\in Q_{h}$, we define $\widehat{u}\in V_{\text{ms}}$, $\widehat{p}\in Q_{\text{ms}}$ as the elliptic projection of $u$ and $p$ with respect to $a$ and $b$, i.e.,

[TABLE]

for all $v\in V_{\text{ms}}$ and $q\in Q_{\text{ms}}$. Moreover, let $\mathcal{A}\colon V_{h}\to V_{h}$ and $\mathcal{B}\colon Q_{h}\to Q_{h}$ be the Riesz’s projections defined by

[TABLE]

for all $v\in V_{h}$ and $q\in Q_{h}$. The elliptic projections $\mathcal{A}$ and $\mathcal{B}$ fulfill the following estimates.

Lemma 4.4**.**

Let $V_{\mathrm{ms}},\,Q_{\mathrm{ms}}$ be the multiscale spaces defined in Section 3 with parameters $m,\ J^{1}_{i},\ J^{2}_{i}$ sufficiently large. Then, for all $v\in V_{h}$ and $q\in Q_{h}$, it holds that

[TABLE]

and

[TABLE]

Proof.

For $v\in V_{h}$, consider the variational problem

[TABLE]

for all $w\in V_{h}$. It was shown in [8, Thm. 2] that for sufficiently large $m$ and $J^{1}_{i}$ the ellipticity of $a(\cdot,\cdot)$ implies $\|z-\widehat{z}\|_{1}\lesssim H\|v-\widehat{v}\|_{1}$. Therefore, we have

[TABLE]

and, thus, $\|v-\widehat{v}\|\lesssim H\|v-\widehat{v}\|_{1}$, which proves (11). Further, we have

[TABLE]

With (11) and the inequality above, we get

[TABLE]

This proves (13). Similarly, one may prove (12) and (14) using the ellipticity of $b(\cdot,\cdot)$. ∎

Remark**.**

In order to achieve the desired estimates in Lemma 4.4, it is necessary to choose the number of oversampling layers m\approx O\big{(}\log(HC_{p})\big{)}, where $C_{p}$ is a constant depending on the Lamé coefficients, the permeability $\kappa$, and $\alpha$, see [8, Sect. 6] for more details. Further, the parameters $J^{1}_{i}$ and $J^{2}_{i}$ account for high contrast within the media and numerical results suggest a logarithmic dependence on the contrast.

The main theoretical result reads as follows.

Theorem 4.1**.**

Assume sufficiently large parameters $m$, $J^{1}_{i}$, $J^{2}_{i}$, a source function $f\in L^{\infty}(0,T;L^{2}(\Omega))\cap H^{1}(0,T;H^{-1}(\Omega))$ as well as initial data $p_{h}^{0}\in Q_{h}$ and $u_{h}^{0}\in V_{h}$ defined in (5). Then, the error between the multiscale solution $(u_{\mathrm{ms}}^{n},p_{\text{ms}}^{n})\in V_{\mathrm{ms}}\times Q_{\mathrm{ms}}$ of (10) and the fine-scale solution $(u_{h}^{n},p_{h}^{n})\in V_{h}\times Q_{h}$ of (4) satisfies

[TABLE]

for $n=1,2,\dots,N$. Here, $\mathcal{S}_{n}$ only depends on the data and is defined by

[TABLE]

Proof.

The proof mainly follows the lines of the proof of [1, Thm. 3.7] and makes use of the results of Lemma 4.4. Thus, we only give the ideas of the proof.

Due to the linearity of the problem, we can decompose the fine-scale solutions into $u_{h}^{n}=\bar{u}_{h}^{n}+\tilde{u}_{h}^{n}$ and $p_{h}^{n}=\bar{p}_{h}^{n}+\tilde{p}_{h}^{n}$, where $(\bar{u}_{h}^{n},\bar{p}_{h}^{n})$ solves (4) with $f=0$ and $(\tilde{u}_{h}^{n},\tilde{p}_{h}^{n})$ solves (4) with $p_{h}^{0}=0$. With this, the stability estimates in Lemma 4.2 can be applied. As a second step, we decompose the difference of the fine-scale and multiscale solution with the help of the introduced projection. The differences $u^{n}_{h}-{\widehat{u}}_{h}^{n}$ and $p^{n}_{h}-{\widehat{p}}_{h}^{n}$ can then be estimated by Lemma 4.4. For the remaining parts ${\widehat{u}}_{h}^{n}-u^{n}_{\mathrm{ms}}$ and ${\widehat{p}}_{h}^{n}-p^{n}_{\mathrm{ms}}$ one needs to consider system (10) with appropriate test functions in combination with properties of the projection and Lemma 4.2. ∎

Remark**.**

Combining Lemma 4.3 and Theorem 4.1, it follows that the proposed multiscale method converges as $H+\tau$ with respect to the $L^{\infty}(0,T;V_{0})$ norm for the displacement and the $L^{2}(0,T;Q_{0})\cap L^{\infty}(0,T;L^{2}(\Omega))$ norm for the pressure.

5 Numerical results

In this section, we present some numerical experiments to demonstrate the performance of our method. The computational domain is set to $\Omega:=(0,1)^{2}$, $T:=100$, and the time step size is chosen as $\tau:=5$. We consider homogeneous Dirichlet boundary conditions for both the pressure and the displacement. The initial pressure is $p^{0}(x):=x_{1}(1-x_{1})x_{2}(1-x_{2})$ and the force $f$ is set to be constant $1$. We assume heterogeneous coefficients that have different values in two subdomains. Considering Figure 3, we call the blue region $\Omega_{1}$ and the red region $\Omega_{2}$.

We choose $\kappa=1$ in $\Omega_{1}$ and $\kappa=10^{4}$ in $\Omega_{2}$, which implied high contrast in the parameters. The Young’s modulus $E$ is set to be equal to the permeability $\kappa$ and the Biot-Willis fluid-solid coupling coefficient $\alpha$ is chosen uniformly distributed on each element, i.e., $\alpha|_{K}\sim U[0.5,1]$ for any $K\in\mathcal{T}^{H}$. We further set $\nu=1$, $M=1$, and $\nu_{p}=0.2$. As shown in Figure 3, the first model includes isolated short and long channels, while the second model is a fracture-type media including randomly distributed fractures. The fine grid size is chosen as $h=\sqrt{2}/200$ in both models.

Recall that $m$ is the number of oversampling layers used to compute the multiscale basis and that $H$ denotes the coarse grid size. Further, we set $J:=J_{i}^{1}=J_{i}^{2}$ as the number of basis functions used in the auxiliary spaces $V_{\text{aux}}$ and $Q_{\text{aux}}$. We use the same parameter settings for the construction of $V_{\text{ms}}$ and $Q_{\text{ms}}$. To quantify the accuracy of CEM-GMsFEM, we define relative weighted $L^{2}$ errors and energy errors for the displacement and the pressure at time $T$,

[TABLE]

[TABLE]

where $(u_{h}(\cdot,T),p_{h}(\cdot,T))$ is the reference solution computed on the fine grid and $(u_{\text{ms}}(\cdot,T),p_{\text{ms}}(\cdot,T))$ is the multiscale solution obtained by the proposed method (10). Based on the above theory, we expect that the multiscale solution converges linearly in $H$ with respect to the energy norms and that the results may be improved by either using more basis functions (increase $J$) or enlarging the oversampling region (increase $m$).

5.1 Test model 1

In this subsection, the numerical results for test model $1$ are presented. First, we investigate the convergence behavior of the CEM-GMsFEM solution with respect to the coarse grid size. As suggested in Remark Remark, we set the number of oversampling layers to $m=4\lfloor\text{log}(H)/\text{log}(\sqrt{2}/10)\rfloor$ and $J=4$ to form the auxiliary spaces. The results are presented in Table 1. One may find that the error decreases as $H$ decreases and that the CEM-GMsFEM solution converges linearly in the energy norms and quadratically in $L^{2}$. In Figure 4 (left), we present the results when using different numbers of oversampling layers with fixed coarse grid size $H=\sqrt{2}/40$.

It is clearly visible that the error decreases faster if more oversampling layers are added. Further, once the number of oversampling layers exceeds a certain number, the error stagnates. To investigate the influence of using different numbers of basis functions, we use a fixed coarse grid with grid size $H=\sqrt{2}/40$ and a fixed number of oversampling layers. The errors in the energy norms are presented in Figure 4 (right). For $m=6$ there is no observable decay for $e_{b}^{p}$ by using more basis functions. One of the reasons for that is the fact that the error for the pressure is already very small even if only $1$ basis is used. On the other hand, due to the high contrast, the error of the displacement $e_{a}^{u}$ decreases as the number of local basis functions increases. Once the number of local basis functions exceeds a certain level, the error decays slower. This happens when also the decay of the eigenvalues slows down. Figures 5 and 6 show the reference solution and the multiscale solution at $T=100$, respectively. One can see that the proposed method captures most of the details of the reference solution.

5.2 Test model 2

In this subsection, we apply the proposed method to test model $2$ with the same choices for $H$, $J$, and $m$ as in the previous subsection. Table 2 shows the results for different values of the coarse grid size $H$ (and thus $m$) with a fixed number of local basis functions. The reference and the multiscale solutions are sketched in Figures 8 and 9, respectively. As before, the impact of the number of oversampling layers and number of basis functions are investigated and depicted in Figure 7. One can observe that the multiscale approximation converges to the reference solution as $H$ decreases. Note that Figure 7 (left) suggests to choose $m\geq 5$, which is in line with the theoretical findings in [8] that the multiscale features are captured when using a sufficiently large number of oversampling layers. As before, also the number of eigenfunctions in the auxiliary spaces can enhance the performance of the multiscale method, see Figure 7 (right).

6 Conclusion

In this work, we have proposed a Generalized Multiscale Finite Element Method based on the idea of constraint energy minimization [8] for solving the problem of linear heterogeneous poroelasticity. The spatial discretization is based on CEM-GMsFEM which provides a framework to systematically construct multiscale basis functions for both displacement and pressure. The multiscale basis functions with locally minimal energy are constructed by employing the techniques of oversampling, which leads to an improved accuracy in the simulations. Combined with the implicit Euler scheme for the time discretization, we have shown that the fully discrete method has optimal convergence rates despite the heterogeneities of the media. Numerical results have been presented to illustrate the performance of the proposed method.