Generalized Multiscale Finite Element Method for the poroelasticity problem in multicontinuum media

TL;DR

This paper develops a multiscale finite element method for simulating poroelasticity in complex heterogeneous multicontinuum media, enabling accurate solutions with fewer computational resources.

Contribution

It introduces a generalized multiscale finite element method tailored for poroelasticity in multicontinuum media, including fracture networks, with demonstrated accuracy in 2D and 3D models.

Findings

Accurate coarse grid solutions with few basis functions.

Effective modeling of fractured porous media.

Reduced computational cost for complex simulations.

Abstract

In this paper, we consider a poroelasticity problem in heterogeneous multicontinuum media that is widely used in simulations of the unconventional hydrocarbon reservoirs and geothermal fields. Mathematical model contains a coupled system of equations for pressures in each continuum and effective equation for displacement with volume force sources that are proportional to the sum of the pressure gradients for each continuum. To illustrate the idea of our approach, we consider a dual continuum background model with discrete fracture networks that can be generalized to a multicontinuum model for poroelasticity problem in complex heterogeneous media. We present a fine grid approximation based on the finite element method and Discrete Fracture Model (DFM) approach for two and three-dimensional formulations. The coarse grid approximation is constructed using the Generalized Multiscale Finite Element Method (GMsFEM), where we solve local spectral problems for construction of the multiscale basis functions for displacement and pressures in multicontinuum media. We present numerical results for the two and three dimensional model problems in heterogeneous fractured porous media. We investigate relative errors between reference fine grid solution and presented coarse grid approximation using GMsFEM with different numbers of multiscale basis functions. Our results indicate that the proposed method is able to give accurate solutions with few degrees of freedoms.

Figures

Tables

Table 1. Table 1: Relative errors (%) for displacement and pressures with different number of multiscale basis functions for two - dimensional problem. Coarse grids 5 × 5 5 5 5\times 5 and 10 × 10 10 10 10\times 10 . D ​ O ​ F h = 57504 𝐷 𝑂 subscript 𝐹 ℎ 57504 DOF_{h}=57504

$M$	$DO{F}_H$	$p_1$	$p_2$	$u$
		$e_{L^2}^{p_1}$ (%)	$e_{L^2}^{p_2}$ (%)	$e_{L^2}^u$ (%)
Coarse mesh $5\times 5$
1	144	28.501	9.514	93.582
2	288	24.552	7.933	90.760
4	576	10.228	3.258	53.017
8	1152	5.695	2.080	27.283
12	1728	2.526	0.887	12.646
16	2304	1.755	0.582	9.076
Coarse mesh $10\times 10$
1	484	24.385	7.825	111.247
2	968	11.947	3.906	40.595
4	1936	5.046	1.766	21.797
8	3872	1.172	0.384	4.406
12	5808	0.270	0.149	3.166
16	7744	0.154	0.083	2.410

Table 2. Table 2: Relative errors (%) for displacement and pressures with different number of multiscale basis functions for two - dimensional problem. Coarse grid 10 × 10 10 10 10\times 10 . D ​ O ​ F h = 57504 𝐷 𝑂 subscript 𝐹 ℎ 57504 DOF_{h}=57504

$M_u$	$DO{F}_H$	$p_1$	$p_2$	$u$
		$e_{L^2}^{p_1}$ (%)	$e_{L^2}^{p_2}$ (%)	$e_{L^2}^u$ (%)
$M_p=1$
1	484	24.385	7.825	111.247
$M_p=2$
1	726	11.905	3.741	72.299
2	968	11.947	3.906	40.595
$M_p=4$
1	1210	5.029	1.723	51.183
2	1452	5.046	1.772	23.402
4	1936	5.046	1.766	21.797
$M_p=8$
1	2178	1.169	0.585	38.968
2	2420	1.170	0.426	9.469
4	2904	1.170	0.406	8.285
8	3872	1.172	0.384	4.406
$M_p=12$
1	3146	0.302	0.504	38.068
2	3388	0.272	0.237	8.940
4	3872	0.271	0.194	7.705
8	4840	0.271	0.154	3.814
12	5808	0.270	0.149	3.166
$M_p=16$
1	4114	0.212	0.485	37.974
2	4356	0.156	0.196	9.044
4	4840	0.155	0.148	7.775
8	5808	0.155	0.095	3.836
12	6776	0.155	0.087	3.211
16	7744	0.154	0.083	2.410

Table 3. Table 3: Relative errors for displacement and pressures with different number of multiscale basis functions for three - dimensional problem. D ​ O ​ F h = 108045 𝐷 𝑂 subscript 𝐹 ℎ 108045 DOF_{h}=108045

$M_u$	$DO{F}_H^u$	$p_1$	$p_2$	$u$
		$e_{L^2}^{p_1}$ (%)	$e_{L^2}^{p_2}$ (%)	$e_{L^2}^u$ (%)
$M_p=1$
1	1080	25.713	29.035	92.385
$M_p=2$
1	1512	20.461	22.644	130.73
2	2106	20.496	22.648	82.179
$M_p=4$
1	2376	6.453	6.991	86.899
2	3024	6.557	7.114	59.549
4	4320	6.626	7.190	22.580
$M_p=8$
1	4104	3.563	3.869	78.466
2	4752	3.647	3.974	44.919
4	6048	3.711	4.050	12.562
8	8640	3.709	4.052	11.773
$M_p=12$
1	5832	2.444	2.710	77.902
2	6480	2.481	2.755	43.286
4	7776	2.530	2.813	10.389
8	10368	2.529	2.813	8.518
12	12960	2.530	2.814	8.518
$M_p=16$
1	7560	1.898	2.080	76.300
2	8208	1.897	2.077	41.473
4	9504	1.889	2.064	9.687
8	12096	1.888	2.063	6.908
12	14688	1.888	2.064	6.170
16	17280	1.887	2.063	5.881

Equations

$$\begin{array}{r}{\alpha }_1\frac{\partial \mathrm{div}u}{\partial t}+\frac{1}{M_1}\frac{\partial p_1}{\partial t}+\mathrm{div}{q}_1+q_{12}+q_{1f}=f_1,x\in \Omega ,\\ q_1=-k_1\mathrm{grad}{p}_1,x\in \Omega ,\\ {\alpha }_2\frac{\partial \mathrm{div}u}{\partial t}+\frac{1}{M_2}\frac{\partial p_2}{\partial t}+\mathrm{div}{q}_2-q_{21}+q_{2f}=f_2,x\in \Omega ,\\ q_2=-k_2\mathrm{grad}{p}_2,x\in \Omega ,\\ {\alpha }_f\frac{\partial \mathrm{div}u}{\partial t}+\frac{1}{M_f}\frac{\partial p_f}{\partial t}+\mathrm{div}{q}_f-q_{f1}-q_{f2}=f_f,x\in \gamma ,\\ q_f=-k_f\mathrm{grad}{p}_f,x\in \gamma ,\end{array}$$

$$q_{12}=q_{21}=r_{12}(p_1-p_2),r_{if}={\eta }_i{r}_{if}(p_i-p_f),r_{fi}={\eta }_f{r}_{if}(p_i-p_f),$$

$$\begin{array}{r}-\mathrm{div}{\sigma }_T(u,p_1,p_2,p_f)=0,x\in \Omega ,\\ {\sigma }_T(u,p_1,p_2,p_f)=\sigma (u)-\sum _{i=1,2,f}{\alpha }_i{p}_i\mathcal{I},x\in \Omega ,\end{array}$$

$$\sigma (u)=2\mu \epsilon (u)+\lambda \mathrm{div}u\mathcal{I},\epsilon (u)=\frac{1}{2}(\mathrm{grad}u+\mathrm{grad}{u}^T),$$

$$\begin{array}{r}{\alpha }_1\frac{\partial \mathrm{div}u}{\partial t}+\frac{1}{M_1}\frac{\partial p_1}{\partial t}-\mathrm{div}\cdot (k_1\mathrm{grad}{p}_1)+q_{12}+q_{1f}=f_1,x\in \Omega ,\\ {\alpha }_2\frac{\partial \mathrm{div}u}{\partial t}+\frac{1}{M_2}\frac{\partial p_2}{\partial t}-\mathrm{div}\cdot (k_2\mathrm{grad}{p}_2)-q_{21}+q_{2f}=f_2,x\in \Omega ,\\ {\alpha }_f\frac{\partial \mathrm{div}u}{\partial t}+\frac{1}{M_f}\frac{\partial p_f}{\partial t}-\mathrm{div}\cdot (k_f\mathrm{grad}{p}_f)-q_{f1}-q_{f2}=f_f,x\in \gamma ,\\ -\mathrm{div}\sigma (u)+{\alpha }_1\mathrm{grad}{p}_1+{\alpha }_2\mathrm{grad}{p}_2+{\alpha }_f\mathrm{grad}{p}_f=0,x\in \Omega .\end{array}$$

$$p_1=p_2=p_f=p^0,u=0,x\in \Omega ,$$

$$\begin{array}{r}{p}_f=g,x\in {\Gamma }_{\gamma }-k_f\frac{\partial p_f}{\partial n}=0,x\in \partial \gamma /{\Gamma }_{\gamma },\\ -k_1\frac{\partial p_1}{\partial n}=0,x\in \Omega ,-k_2\frac{\partial p_2}{\partial n}=0,x\in \Omega ,\\ u=0,x\in {\Gamma }_D\sigma =0,x\in {\Gamma }_N,\end{array}$$

$$\begin{array}{r}{\alpha }_i\frac{\partial \mathrm{div}u}{\partial t}+\frac{1}{M_i}\frac{\partial p_i}{\partial t}-\mathrm{div}\cdot (k_i\mathrm{grad}{p}_i)+\sum _{j\ne i}{r}_{ij}(p_i-p_j)=f_i,x\in \Omega ,\\ -\mathrm{div}\sigma (u)+\sum _j{\alpha }_j\mathrm{grad}{p}_j=0,x\in \Omega .\end{array}$$

$$V=\{v\in H^1(\Omega ):v=0\text{ on }{\Gamma }_D\},W_1=W_2=H^1(\Omega ),$$

$$W_f=\{w\in H^1(\gamma ):w=g\text{ on }{\Gamma }_{\gamma }\},{\hat{W}}_f=\{w\in H^1(\gamma ):w=0\text{ on }{\Gamma }_{\gamma }\}.$$

$$\begin{array}{rl}& d_i\left(\frac{\partial u}{\partial t},w_i\right)+c_i\left(\frac{\partial p_i}{\partial t},w_i\right)+b_i(p_i,w_i)+\sum _j{q}_{ij}(p_i-p_j,w_i)=l(w_i),\forall w_i\in W_i,\\ & a(u,v)+\sum _j{g}_i(p_i,v)=0,\forall v\in V,\end{array}$$

$$b_i(p_i,w_i)={\int }_{{\Omega }_i}{k}_i\nabla p_i\cdot \nabla w_{i}dx,l(w_i)={\int }_{{\Omega }_i}{f}_i{w}_{i}dx,$$

$$c_i(p_i,w_i)={\int }_{{\Omega }_i}\frac{1}{M_i}{p}_i{w}_{i}dx,d_i(u,w_i)={\int }_{{\Omega }_i}{\alpha }_i\mathrm{div}u{w}_{i}dx,$$

$$q_{ij}(p_i-p_j,w_i)={\int }_{{\Omega }_i}{r}_{ij}(p_i-p_j)w_{i}dx,$$

$$a(u,v)={\int }_{\Omega }\sigma (u)\cdot \epsilon (v)dx,g_i(p_i,v)={\int }_{{\Omega }_i}{\alpha }_i\mathrm{grad}{p}_{i}vdx,$$

$$C\frac{y^h-{\check{y}}^h}{\tau }+A{y}^h=F,$$

$$C=\left(\begin{array}{cccc}{C}_1& 0& 0& D_1\\ 0& C_2& 0& D_2\\ 0& 0& C_f& D_f\\ 0& 0& 0& 0\end{array}\right),A=\left(\begin{array}{cccc}{A}_1+Q_{12}+Q_{1f}& -Q_{12}& -Q_{1f}& 0\\ -Q_{12}& A_2+Q_{12}+Q_{2f}& -Q_{2f}& 0\\ -Q_{1f}& -Q_{2f}& A_f+Q_{1f}+Q_{2f}& 0\\ D_1^T& D_2^T& D_f^T& A_u\end{array}\right)$$

$$A_i=[a_{i,ln}],a_{i,ln}={\int }_{{\Omega }_i}{k}_i\mathrm{grad}{\varphi }_l^i\cdot \mathrm{grad}{\varphi }_n^{i}dx,A_u=[a_{u,ln}],a_{u,ln}={\int }_{\Omega }\sigma ({\Phi }_l)\cdot \epsilon ({\Phi }_n)dx,$$

$$Q_{ij}=[q_{ij,ln}],q_{ij,ln}={\int }_{{\Omega }_i}{r}_{ij}{\varphi }_l^i{\varphi }_n^{j}dx,C_i=[c_{i,ln}],c_{i,ln}={\int }_{{\Omega }_i}\frac{1}{M_i}{\varphi }_l^i{\varphi }_n^{i}dx,$$

$$D_i=[d_{i,ln}],d_{i,ln}={\int }_{\Omega }{\alpha }_i\mathrm{div}{\Phi }_l{\varphi }_n^{i}dx,F_i=[f_{i,l}],f_{i,l}={\int }_{{\Omega }_i}{f}_i{\varphi }_l^{i}dx,$$

$$\tilde{C}\frac{y^h-{\check{y}}^h}{\tau }+\tilde{A}{y}^h=\tilde{F},$$

$$\tilde{C}=\left(\begin{array}{ccc}{C}_1+C_f& 0& D_1\\ 0& C_2& D_2\\ 0& 0& 0\end{array}\right),\tilde{A}=\left(\begin{array}{ccc}{A}_1+A_f+Q_{12}& -Q_{12}& 0\\ -Q_{12}& A_2+Q_{12}& 0\\ D_1^T& D_2^T& A_u\end{array}\right),\tilde{F}=\left(\begin{array}{c}{F}_1+F_f\\ F_2\\ 0\end{array}\right).$$

$${\mathcal{T}}^H=\bigcup _j{K}_j,$$

$$b_i({\psi }_i^{l,j},w_i)+\sum _j{q}_{ij}({\psi }_i^{l,j}-{\psi }_j^{l,j},w_i)=0,\forall w_i\in {\hat{W}}_i^h$$

$$W_i^h=\{w\in H^1({\omega }_l):w={\delta }_i^j\text{ on }\partial {\omega }_l\},{\hat{W}}_i^h=\{w\in H^1({\omega }_l):w=0\text{ on }\partial {\omega }_l\},$$

$$W_{snap}({\omega }_l)=\text{span}\{{\psi }^{l,j},l=1,...,N_v^H,j=1,...,L_p^{{\omega }_l}\}.$$

$${\tilde{A}}_p{\tilde{\varphi }}^l={\lambda }_p{\tilde{S}}_p{\tilde{\varphi }}^l,$$

$${\tilde{A}}_p=R_{snap}^p{A}_p(R_{snap}^p{)}^T,{\tilde{S}}_p=R_{snap}^p{S}_p(R_{snap}^p{)}^T,R_{snap}^p=({\psi }^{l,1},...,{\psi }^{l,L_p^{{\omega }_l}}{)}^T$$

$$S_p=\left(\begin{array}{ccc}{S}_1& 0& 0\\ 0& S_2& 0\\ 0& 0& S_f\end{array}\right),A_p=\left(\begin{array}{ccc}{A}_1+Q_{12}+Q_{1f}& -Q_{12}& -Q_{1f}\\ -Q_{12}& A_2+Q_{12}+Q_{2f}& -Q_{2f}\\ -Q_{1f}& -Q_{2f}& A_f+Q_{1f}+Q_{2f}\end{array}\right)$$

$$A_i=[a_{i,mn}],a_{i,mn}={\int }_{{\omega }_l^i}{k}_i\mathrm{grad}{\varphi }_m^i\cdot \mathrm{grad}{\varphi }_n^{i}dx,S_i=[s_{i,mn}],s_{i,mn}={\int }_{{\omega }_l^i}{k}_i{\varphi }_m^i{\varphi }_n^{i}dx.$$

Full text

1 Introduction

Mathematical models of the poroelasticity problems in multicontinuum media for heterogeneous fractured media are used for simulation of the unconventional hydrocarbon reservoirs, geothermal fields, underground disposal of radioactive waste in subsurface collectors, etc [1, 2, 3, 4, 5]. Mathematical model is described by a system of equations for pressures in each continuum and equation for displacements. The most important feature of the mathematical model is that the equations are coupled. We can highlight two main coupling types: (1) multicontinuum pressures coupling via mass transfer term [6, 7, 8], and (2) pressure and displacement coupling via term that describes the compressibility of the medium and volume force which is proportional to the pressure gradient [9, 10].

Mathematical models for flow in the multicontinuum media is used for describing of the complex flow processes in multiscale fractured heterogeneous porous media [11, 12, 13, 14]. The flow in fractures has a significant impact on filtration processes and requires careful consideration [15, 16]. Moreover, since the fractures are characterized by high permeability and their thickness is significantly smaller than the size of the simulated field, this leads to the need to build special mathematical models of multicontinuum, where independent variables are distinguished to describe the flow in a porous medium and in the network of fractures taking into account the special flow function [17, 18, 19]. At the same time, the size of fractures should be separated, since they can exist on different scales and can differ in the nature of their occurrence. In the case of naturally fractured porous media, the fracture system is basically connected and dual porosity models are traditionally used [6, 7]. The interaction of the continua is described by specifying the flow functions between continua (mass transfer) [20, 18, 3].

For numerical simulation of the problems in fractured and heterogeneous porous media, we should use a sufficiently fine grid that resolve all small scale features in the level of mesh construction. Therefore, the discrete formulation of such problems leads to the large system of equations that is computationally expensive. To reduce the size of the discrete system, various multiscale methods have been developed [21, 22, 23, 24, 25, 26]. In [27, 28], the multiscale finite volume method was presented for solution of the flow problems in fractured porous media. Multiscale finite volume method for solution of the poroelasticity problem is presented in [29]. Generalized multiscale finite element (GMsFEM) for solution of the flow problems in fractured porous media is considered in [30, 31]. In [32, 33], we considered construction of the coarse grid problem for poroelasticity problem in heterogeneous media. Extension of the GMsFEM method for solution of the poroelasticity problem in fractured media with discontinuous Galerkin method for displacements and continuous Galerkin method for pressure is presented in our previous work [17]. Recently a new method was presented for solution of the flow and poroelasticity problems in fractured and heterogeneous porous media [34, 35, 36]. Generalization of the GMsFEM and NLMC approach for solution of the flow problems in multicontinuum media is considered in [3, 37].

In this paper, we consider the Generalized Multiscale Finite Element method for solution of the poroelasticity problems in multicontinuum heterogeneous media. We construct a multiscale basis functions for pressures and displacement. Construction of the basis functions for flow problem in multicontinuum media is based on the solution of the coupled system of equations in each local domains. Method automatically identify each continuum flow features via solution of the local spectral problems. For displacement, we use similar approach where main features are captured by the local spectral problem. We numerically investigate presented method for two and three-dimensional model problems in fractured and heterogeneous porous media.

The work is organized as follows. In Section 2, we present the mathematical model of the poroelasticity problem in multicontinuum medium. Then in Section 3, a fine grid approximation is constructed using the finite element method. In Section 4, we present a coarse grid approximation using Generalized Multiscale Finite Element method, where we describe construction of the multiscale basis functions and coarse grid system construction. Numerical results for two and three-dimensional model poroelasticity problems are presented in Section 5. Finally, we present Conclusions.

2 Mathematical model

Let $\Omega\subset R^{d}$ is the computational domain for background medium with dual continuum approach, where $d=2$ for two-dimensional problems and $d=3$ for three-dimensional problems. For example, the first continuum can describe a flow in the matrix of the porous media, and the second continuum belongs to the network of small highly connected fracture network. Furthermore, we define $\gamma\in R^{d-1}$ as a computational domain for low dimensional fracture networks model that describes the flow in the large-scale fractures. For the fluid flow in the poroelastic medium, we have the mass balance equation and the Darcy’s law

[TABLE]

where $u$ is the displacement, $q_{i}$ is the velocity, $p_{i}$ is the pressure, $\kappa_{i}$ is the permeability ($k_{i}=\kappa_{i}/\mu$, $\mu$ is the fluid viscosity), $f_{i}$ refer to source and sink terms, $\alpha_{i}$ is Biot coefficients, $M_{i}$ is the Biot modulus and $i$ is continuum index ($i=1,2,f)$.

Here $q_{12}$, $q_{1f}$ and $L_{2f}$ are the transfer term between first and second continua; fist continuum and fractures; and second continuum and fractures, respectively. We have

[TABLE]

and $\int_{\Omega}q_{12}dx=\int_{\Omega}q_{21}dx$, $\int_{\Omega}q_{if}dx=\int_{\gamma}q_{fi}ds$, $\eta_{i}$ is the geometric factor and $r_{ij}$ is the mass transfer term that proportional to the continuum permeabilities.

For the mechanics of the poroelastic multicontinuum media, we use an effective equation for displacement with volume force sources that proportional to the sum of the pressure gradients for each continuum.

[TABLE]

where $\sigma_{T}$ is the total stress tensor, $\sigma$ is the stress tensor [17, 36, 10]. In the case of a linear elastic stress-strain constitutive relation, we have

[TABLE]

where $\varepsilon$ is the strain tensor, $\lambda$ and $\mu$ are the Lame’s coefficients.

Then, we have the following coupled system of equations

[TABLE]

We consider a system of equations (3) with the following initial conditions

[TABLE]

and boundary conditions

[TABLE]

where boundaries $\Gamma_{D}\cup\Gamma_{N}=\partial\Omega$ and $\Gamma_{\gamma}=\partial\Omega\cap\partial\gamma$.

We can generalize presented model as poroelasticity model for multicontinuum media

[TABLE]

where $i=1,...,M$ and $M$ is the number of continua.

3 Fine grid finite element approximation

For the approximation of the system of equations (3) with boundary conditions (5), we use a finite element method. To do so, we define the following functional spaces

[TABLE]

[TABLE]

The variational formulation of the poroelasticity problem in multicontinuum media can be written as follows: find $(p_{1},p_{2},p_{f},u)\in W_{1}\times W_{2}\times W_{f}\times V$ such that

[TABLE]

where bilinear and linear forms are following

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and $\Omega_{1}=\Omega_{2}=\Omega$, $\Omega_{f}=\gamma$, $i=1,2,f$.

Let $\mathcal{T}^{h}$ denote a finite element partition of the domain $\Omega$. For the fracture continuum, we use a discrete fracture model and use an unstructured fine grid $\mathcal{T}^{h}$ that explicitly resolve fracture geometry. We assume that $\cup_{j}\gamma_{j}$ is the subset of faces for $\mathcal{T}^{h}$ that represent fractures, where $j=1,...,N_{frac}$ and $N_{frac}$ is the number of discrete fractures. For approximation by time, we use an implicit finite difference scheme with time step $\tau$. Therefore, we have following discrete system in matrix form on the fine grid for the triple-continuum media and $y^{h}=(p_{1}^{h},p_{2}^{h},p_{f}^{h},u^{h})^{T}$

[TABLE]

where

[TABLE]

and $F=(F_{1},F_{2},F_{f},0)^{T}$, $\check{y}^{h}$ is the solution from the previous time step. Here

[TABLE]

[TABLE]

[TABLE]

and $p_{i}=\sum_{l}p_{i,l}^{h}\phi^{i}_{l}$, $u=\sum_{l}u_{l}^{h}\Phi_{l}$, where $\Phi_{l}$ is the linear basis functions for displacements, $\phi^{1}_{l}=\phi^{2}_{l}$ is the $d$ - dimensional linear basis functions for pressure, $\phi^{f}_{l}$ is the $(d-1)$ - dimensional linear basis functions for pressure.

In this paper, for simplification of the matrix construction, we use a modified DFM approach and consider the case when $\alpha_{f}=0$, $\sigma_{2f}=0$. We assume that $p^{h}_{1}=p^{h}_{f}$ and using superposition principle [30, 31], we eliminate $p^{h}_{f}$ from equations and obtain following coupled system of equations for $y^{h}=(p_{1}^{h},p_{2}^{h},u^{h})^{T}$

[TABLE]

where

[TABLE]

with matrices of the size $N_{h}=4N_{v}$, $N_{v}$ is the number of vertices in $\mathcal{T}^{h}$.

4 Coarse grid approximation using GMsFEM

For construction of the coarse grid approximation of the poroelasticity problems in fractured and heterogeneous media, we use a Generalized Multiscale Finite Element Method (GMsFEM). GMsFEM contains following steps:

coarse grid and local domains construction; 2. 2.

solution of the local problems with different boundary conditions to construct a snapshot space in each local domain; 3. 3.

multiscale basis functions construction via solution of the local spectral problems on the snapshot space; 4. 4.

generation of the projection matrix using local multiscale basis functions; 5. 5.

construction of the coarse grid system using projection matrix; 6. 6.

solution of the unsteady problem on the coarse grid and reconstruction of the fine grid solution.

In this computational algorithm, first four steps are offline (preprocessing) steps for a given fracture geometry and heterogeneity. Fifth step is also offline for linear problems and time - independent right-hand side, but should be online step for nonlinear problems, where fine grid system is change on each nonlinear or/and time iteration. After that on the sixth (online) step, we can perform fast and accurate solution of the reduced order model on the coarse grid.

In this work, we construct multiscale basis functions for displacements and pressures separately, but basis functions for multicontinuum pressure equations are constructed in the coupled way. We start with definition of the snapshot space, after that we define a local eigenvalue problems for pressures and displacements. Finally, we define projection matrix and present construction of the coarse grid poroelasticity system on the multiscale space. Let $\mathcal{T}^{H}$ is the coarse grid partitioning of the domain

[TABLE]

where $K_{j}$ is the coarse grid cell. We will use a continuous Galerkin approximation on the coarse grid, and define local domain $\omega_{l}$ for multiscale basis functions as combination of the several coarse grid cells that share same coarse grid node ($l=1,...,N_{v}^{H}$, $N_{v}^{H}$ is the number of coarse grid vertices).

Multiscale basis functions for pressures in multicontinuum media. To construct a snapshot space, we solve following local problem in domain $\omega_{l}$: find $\psi^{l,j}=(\psi^{l,j}_{1},...,\psi^{l,j}_{M})\in W^{h}_{1}\times...\times W^{h}_{M}$ such that

[TABLE]

where

[TABLE]

and $\delta^{j}_{i}$ is the piecewise constant function (delta function) for $j=1,..,N_{v}^{\omega_{l}}$ ($,N_{v}^{\omega_{l}}$ is the number of nodes on the computation mesh for $\omega_{l}$), $i$ is the index of continuum ($i=1,...,M$). Therefore, we solve $L_{p}^{\omega_{l}}=N_{v}^{\omega_{l}}\cdot M$ local problems.

We define snapshot space for pressures in multicontinuum media as follows.

[TABLE]

Next, we solve following local spectral problem on the snapshot space

[TABLE]

where $\hat{\phi}^{l}=(R^{p}_{snap})^{T}\tilde{\phi}^{l}$ and

[TABLE]

Here for matrices in triple continuum case, we have

[TABLE]

where

[TABLE]

We choose an eigenvector $\hat{\phi}_{j}$ ($j=1,..,M^{l,p}$) corresponding to the first smallest $M^{l,p}$ eigenvalues and multiply to the linear partition of unity functions $\chi^{l}$ for obtaining conforming basis functions

[TABLE]

where $\phi^{l,j}=\chi^{l}\hat{\phi}^{l,j}$.

Multiscale basis functions for displacements. We construct the multiscale basis functions by solution following problem in local domain $\omega_{l}$: find $\Psi^{l,j}\in V^{h}$ such that

[TABLE]

where

[TABLE]

and $\bar{\delta}^{j}_{i}$ is the vector for each component for $d$ - dimensional problem ($d=2,3$) i.e. $\bar{\delta}^{j}_{i}=({\delta}^{j}_{i},0,0)$ or $\bar{\delta}^{j}_{i}=(0,{\delta}^{j}_{i},0)$ or $\bar{\delta}^{j}_{i}=(0,0,{\delta}^{j}_{i})$ for $d=3$. We solve $L_{u}^{\omega_{l}}=d\cdot N_{v}^{\omega_{l}}$ local problems.

We define snapshot space for pressures in multicontinuum media as follows

[TABLE]

For construction of the multiscale basis, we solve following local spectral problem on the snapshot space

[TABLE]

where $\hat{\Phi}^{l}=(R^{u}_{snap})^{T}\tilde{\Phi}^{l}$,

[TABLE]

and

[TABLE]

We choose an eigenvector $\hat{\Phi}_{j}$ ($j=1,..,M^{l,u}$) corresponding to the first smallest $M^{l,u}$ eigenvalues and multiply to the linear partition of unity functions for obtaining conforming basis functions

[TABLE]

where $\Phi^{l,j}=\chi^{l}\hat{\Phi}^{l,j}$.

Coarse grid system. Using constructed multiscale basis functions, we define projection matrix

[TABLE]

where

[TABLE]

[TABLE]

Finally, we obtain following reduced order model

[TABLE]

where $C^{H}=RCR^{T}$, $A^{H}=RAR^{T}$ and $F^{H}=RF$. with matrices of the size $N_{f}=4N_{v}$, $N_{v}$ is the number of vertices in $\mathcal{T}^{h}$. Size of the system is $N_{H}=\sum_{l=1}^{N^{H}_{v}}(M^{l,p}+M^{l,u})$ or $N_{H}=(M^{p}+M^{u})\cdot N^{H}_{v}$ for $M^{p}=M^{l,p}$ and $M^{u}=M^{l,u}$ ($\forall l=1,...,N^{H}_{v}$), where $N^{H}_{v}$ is the number of vertices of coarse grid $\mathcal{T}^{H}$.

After obtaining of a coarse-scale solution, we reconstruct fine-scale solution

[TABLE]

We note that, in general, multiscale basis functions for the displacements and pressures can be calculated by solution coupled poroelasticity problem similarly to the multicontinuum pressures.

5 Numerical results

In this section, we consider poroelasticity problem in fractured and heterogeneous media. We consider two problems in $\Omega=[0,10]^{d}$ ($d=2,3$):

• •

Two - dimensional model problem. Fine grid contains 14376 vertices and 28350 cells. Coarse grid contains 121 vertices and 100 cells (Figure 1);

• •

Three - dimensional model problem. Fine grid contains 21609 vertices and 118500 cells. Coarse grid contains 216 vertices and 125 cells (Figure 2).

We set parameters of model problem as follows: $\alpha_{1}=1$, $\alpha_{2}=1$, $\alpha_{f}=0$, $M_{1}=M_{2}=M_{f}=10^{8}$. The calculation is performed by $T_{max}=100$ with time step $\tau=10$ for two - dimensional problem. For three - dimensional problem, we set $T_{max}=6000$ with time step $\tau=300$. Initial condition is $p^{0}=0.1*10^{7}$ and $g=0.2*10^{7}$ for pressure boundary condition.

To compare the results, we use relative errors $L_{2}$ between multiscale solution and fine-scale solution

[TABLE]

where $i$ is the index of the continuum ($i=1,2$), $y^{ms}=(p^{ms}_{1},p^{ms}_{2},u^{ms})$ is the multiscale solution using GMsFEM and $y=(p_{1},p_{2},u)$ is the fine grid solution.

We use $DOF_{h}$ (Degree of Freedom) to denote fine grid system size and $DOF_{H}$ to denote problem size of the coarse scale system using GMsFEM. On GMsFEM, we choose a same number of multiscale basis functions in each local domains $M=M_{p}=M_{u}$, where $M_{p}$ and $M_{u}$ are the number of basis functions. We use GMSH software to construct computational domains and grids [38]. The implementation is based on the open-source library FEniCS [39].

5.1 Two-dimensional problem

We simulate a two-dimensional model problem in heterogeneous fractured porous media. We set $k_{f}=10^{-4}$ and $r_{12}=1000\cdot k_{2}$. Heterogeneous coefficients for elasticity modulus and heterogeneous permeability for first and second continuum are presented in Figure 3.

In Figures 4 distribution of pressure for first continuum and second continuum, displacement along $X$ and $Y$ directions at final time are presented. On the first row, we depict a fine scale solution and multiscale solution with 16 multiscale basis functions for GMsFEM is presented on second row. Comparing the fine-scale solution with the multiscale solution with 16 basis functions in Figure 4 for displacement along $X$ and $Y$ pressure, we can observe good accuracy. In general, solutions look good without visible oscillations.

In Table 1, we present an errors for $5\times 5$ and $10\times 10$ coarse grid with $M=M_{p}=M_{u}$. The results show that 8 multiscale basis functions are enough to achieve good results with $4.406\%$ of $L^{2}$ error for displacement, $1.172\%$ and $0.384\%$ of $L_{2}$ errors for first continuum and second continuum pressures. When we increase number of the multiscale basis functions, the relative $L^{2}$ errors are decrease two times when we take two multiscale basis functions instead of one. We have similar improvements for further increasing of the multiscale basis functions.

Table 2 shows a comparison of the difference between multiscale and fine grid solutions, where we present relative $L_{2}$ for different number of the multiscale basis functions on the coarse grid $10\times 10$. From the Table 2, we observe that number of basis functions for pressure $M_{p}$ highly impact to the displacements errors. For example in the case $M_{u}=2$, we have $40\%$ of $L^{2}$ displacements error when we have 2 multiscale basis functions for pressure, and reduce error to $9\%$ for $M_{p}=8$. In conclusion, according to the results of the comparison, we observe a good convergence of the method when we increase number of the multiscale basis functions.

5.2 Three-dimensional problem

We simulate a three-dimensional model problem in heterogeneous fractured porous media. We set $k_{f}=10^{-4}$ and $r_{12}=25\cdot k_{2}$. Heterogeneous coefficients for elasticity modulus and heterogeneous permeability for first continuum and second continuum are presented in Figure 5.

In Figures 6-7 distribution of pressure for first continuum and second continuum, displacement along $X$,$Y$ and $Z$ directions at the final time are presented. On the first row, we depict a fine scale solution and multiscale solution with 16 multiscale basis functions for GMsFEM is presented on second row. We observe a good results of the multiscale method compared with the fine grid solution.

Table 3 shows a comparison of the difference between multiscale and fine grid solutions. We present relative $L^{2}$ for displacements and pressures. For three - dimensional model problem using 8 multiscale basis functions, we have $11.773\%$ of $L^{2}$ error for displacement, $3.709\%$ and $4.052\%$ of $L_{2}$ errors for first continuum and second continuum pressures. When we take 16 multiscale basis functions, we obtain two times better results with $5.881\%$ of $L^{2}$ error for displacement, $1.887\%$ and $2.063\%$ of $L^{2}$ errors for first continuum and second continuum pressures. For the three-dimensional problem, we also observe good convergences for the poroelasticity problem in heterogeneous and fractured media.

6 Conclusion

In this work, we considered the poroelasticity problem in heterogeneous and fractured medium. We presented a mathematical model and fine grid approximation using finite element method for general multicontinuum poroelasticity problem. We developed a generalized multiscale finite element framework and explained construction of the multiscale basis functions. We presented results of the numerical investigation for model problems in two and three-dimensional formulations. We compared a relative error between multiscale and fine-scale solutions for different number of multiscale basis functions. The proposed multiscale method provides a good accuracy for two and three-dimensional problems in heterogeneous fractures media with few degrees of freedoms.

7 Acknowledgements

AT’s, MV’s and DS’s works are supported by the mega-grant of the Russian Federation Government N14.Y26.31.0013 and RSF N17-71-20055. The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund (Project numbers 14304217 and 14302018) and CUHK Faculty of Science Direct Grant 2018-19.