Advanced modeling and simulation of flow in subsurface reservoirs with fractures and wells for a sustainable industry
Open Access
Issue
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 75, 2020
Advanced modeling and simulation of flow in subsurface reservoirs with fractures and wells for a sustainable industry
Article Number 33
Number of page(s) 8
DOI https://doi.org/10.2516/ogst/2020025
Published online 05 June 2020

© M.F. El-Amin et al., published by IFP Energies nouvelles, 2020

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Nomenclature

bd : Factor of slippage

k0,d: Intrinsic permeability

ϕ m $ {\phi }_m$ : Porosity of the matrix blocks

ϕ f $ {\phi }_f$ : Porosity of the fractures

ρ m $ {\rho }_m$ : Mass density of the gas in the matrix blocks

ρ f $ {\rho }_f$ : Mass density of the gas in the fractures

ρ s $ {\rho }_s$ : Mass density of the core

M w $ {M}_w$ : Gas molar weight

P L $ {P}_L$ : Langmuir pressure

p m $ {p}_m$ : Gas pressure in matrix blocks

p f $ {p}_f$ : Gas pressure in fractures

p c $ {p}_c$ : Critical pressure

p fe $ {p}_{\mathrm{fe}}$ : Average fractures pressure around the well

V L $ {V}_L$ : Langmuir volume

V std $ {V}_{\mathrm{std}}$ : Mole volume under standard conditions

R $ R$ : General constant of gases

T $ T$ : Temperature

T c $ {T}_c$ : Critical temperature

Z $ Z$ : Gas compressibility factor

r w $ {r}_w$ : Well-radius

r e $ {r}_e$ : Drainage-radius

n $ n$ : Set of normal fractures

σ m $ {\sigma }_m$ : Mean stress in matrix blocks

σ f $ {\sigma }_f$ : Mean stress in fractures

σ m $ {\sigma }_m^\mathrm{\prime}$ : Effective stress in matrix blocks

σ f $ {\sigma }_f^\mathrm{\prime}$ : Effective stress in fractures

α d $ {\alpha }_d$ : Biot’s effective parameter

L x $ {L}_x$ : Fracture spacing in x$ x$

L y $ {L}_y$ : Fracture spacing in y$ y$

L z $ {L}_z$ : Fracture spacing in z$ z$

Superscripts and subscripts

f $ f$ : Fractures

m $ m$ : Matrix

d $ d$ : m$ m$ or f$ f$

1 Introduction

Finite Element Methods (FEMs) are effective numerical techniques for solving the complex engineering problems. FEMs appeared in the middle of the last century and were widely used in solid mechanical applications [1]. The FEM has been extensively investigated and developed in the last five decades [2], and now is used for highly nonlinear problems [3]. However, the standard FEM faced certain instabilities for solving elliptic problem that describe the flow in porous media [4, 5]. On the other hand, the Mixed Finite Element Methods (MFEMs) have succeeded in eliminating such instabilities [6, 7] as it may be extended to higher-order approximations as well as it is a locally conservative method [8]. Different physical variables can be treated differently in the MFEM from other variables. For example, the enhanced linear Brezzi-Douglas-Marini (BDM1) space is used for speed approximation, which requires more precision, while the piece-wise constant space is used for pressure approximation. The models of shale-gas transport were mainly developed by adapting the traditional models of flow in fractured porous media. Warren and Root [9] have developed an idealized geometric model to investigate the behavioral flow in fractured porous media. The dual porosity model including stress of fracture was used to investigate the flow of gas in the shale [10, 11], while, the dual-continuum model was used to describe gas flow Kerogen shale [12, 13]. Several versions of the dual-continuum model have been developed to include, for instance, adsorption, heat, deformation and Knudsen diffusion [1417]. Several authors presented research work in shale gas reservoirs with geomechanical effects, such as Wang et al. [18], Lin et al. [19], and Yang et al. [20]. Wang et al. [18] developed a general model including embedded discrete fractures, multiple interacting continua and geomechanics in shale gas reservoirs with multi-scale fractures. In [19], authors provide insights on fracture propagation using reservoir flow simulation integrated with geomechanics in unconventional tight reservoirs. Effects of multicomponent adsorption and geomechanics are investigated in a shale condensate reservoir [20]. Girault et al. [21] have introduced a priori error estimates for a discretized poro-elastic-elastic system, in which the flow pressure equation is discretized by either a continuous Galerkin scheme or a mixed scheme, while, elastic displacement equations are discretized by a continuous Galerkin scheme. El-Amin et al. [22] have used the MFEM with stability analysis to simulate the problem of natural gas transport in a low-permeability reservoir without considering fractures. The authors [23] extended their work to cover fractured porous media and the rock stress-sensitivity with considered stability analysis of the MFEM. In this work, we present a theoretical basis with proofs of the stability analysis of the MFEM (in Ref. [23]) including the necessary lemmas and theorem.

2 Modeling and formulation

In this section, the mathematical model of the problem under consideration is developed. The Dual Porosity Dual Permeability (DPDP) model is employed to describe the gas transport in fractured porous media which consists of matrix blocks and fractures. The absorbed gas and the free gas coexist in the matrix blocks, however, the free gas exists in fractures only. The mass of free gas accumulation is ϕmρm$ {\phi }_m{\rho }_m$, and the adsorbed gas accumulation (Langmuir isotherm model) is [24, 25],(1-ϕm)ρsMwVLpmVstd(PL+pm).$$ \frac{(1-{\phi }_m){\rho }_s{M}_w{V}_L{p}_m}{{V}_{\mathrm{std}}({P}_L+{p}_m)}. $$(1)

The mass density of the gas can be written as,ρg,d=pdMwZRT, d=m,f.$$ {\rho }_{g,d}=\frac{{p}_d{M}_w}{\mathrm{ZRT}},\enspace \hspace{1em}d=m,f. $$(2)

The gas compressibility factor, Z$ Z$ which is given by Peng-Robinson equation of state [26],Z3-(1-B)Z2+(A-3B3-2B)Z-(AB-B2-B3)=0,$$ {Z}^3-(1-B){Z}^2+(A-3{B}^3-2B)Z-({AB}-{B}^2-{B}^3)=0, $$(3)A=aTpR2T2, B=bTpRT,$$ A=\frac{{a}_Tp}{{R}^2{T}^2},\hspace{1em}\enspace B=\frac{{b}_Tp}{{RT}}, $$(4)aT=0.45724R2Tc2pc, bT=0.0778RTcpc.$$ {a}_T=0.45724\frac{{R}^2{T}_c^2}{{p}_c},\enspace \hspace{1em}{b}_T=0.0778\frac{R{T}_c}{{p}_c}. $$(5)

In low-permeability formations, Klinkenberg effect, which is described by the apparent permeability, takes place and written as,kapp,d=k0,d(1+bdpd), d=m,f.$$ {k}_{\mathrm{app},d}={k}_{0,d}\left(1+\frac{{b}_d}{{p}_d}\right),\hspace{1em}\enspace d=m,f. $$(6)

The DPDP model of gas transport in fracture shale strata is represented as,f1(pm)pmt-[ρg,m(pm)μk0,m(1+bmpm)pm]=-S(pm,pf),$$ {f}_1({p}_m)\frac{\mathrm{\partial }{p}_m}{\mathrm{\partial }t}-\nabla \cdot \left[\frac{{\rho }_{g,m}({p}_m)}{\mu }{k}_{0,m}\left(1+\frac{{b}_m}{{p}_m}\right)\nabla {p}_m\right]=-S({p}_m,{p}_f), $$(7)andf2(pf)pft-[ρg,f(pf)μk0,f(1+bfpf)pf]=S(pm,pf)-Q(pf),$$ {f}_2({p}_f)\frac{\mathrm{\partial }{p}_f}{\mathrm{\partial }t}-\nabla \cdot \left[\frac{{\rho }_{g,f}({p}_f)}{\mu }{k}_{0,f}\left(1+\frac{{b}_f}{{p}_f}\right)\nabla {p}_f\right]=S({p}_m,{p}_f)-Q({p}_f), $$(8)wheref1(pm)=MwZRT[ϕm(pm)+ϕm'(pm)pm]+MwVLρsVstd(PL(1-ϕm(pm))(PL-pm)2-ϕm'(pm)PL+pm),$$ {f}_1\left({p}_m\right)=\frac{{M}_w}{\mathrm{ZRT}}\left[{\phi }_m\left({p}_m\right)+{\phi }_m^{\prime}\left({p}_m\right){p}_m\right]+\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}\left(\frac{{P}_L\left(1-{\phi }_m\left({p}_m\right)\right)}{{\left({P}_L-{p}_m\right)}^2}-\frac{{\phi }_m^{\prime}\left({p}_m\right)}{{P}_L+{p}_m}\right), $$(9)andf2(pf)=MwZRT[ϕf(pf)+ϕ'f(pf)pf],$$ {f}_2({p}_f)=\frac{{M}_w}{\mathrm{ZRT}}\left[{\phi }_f({p}_f)+{{\phi \prime}_f({p}_f){p}_f\right], $$(10)such that ϕm$ {\phi }_m$ and ϕf$ {\phi }_f$ is the matrix and fractures porosities, respectively,ϕ'd(pd)=dϕddpd, d=m,f.$$ {{\phi \prime}_d\left({p}_d\right)=\frac{\mathrm{d}{\phi }_d}{\mathrm{d}{p}_d},\enspace \hspace{1em}d=m,f. $$(11)

Moreover, based on the following definitions, bm$ {b}_m$ and bf$ {b}_f$ are treated as constants [17],bm=8πRTMw1r(2α-0.995)μg,$$ {b}_m=\sqrt{\frac{8{\pi RT}}{{M}_w}}\frac{1}{r}\left(\frac{2}{\alpha }-0.995\right){\mu }_g, $$(12)andbf=πRTϕ0,fMwk0,f μg.$$ {b}_f=\sqrt{\frac{{\pi RT}{\phi }_{0,f}}{{M}_w{k}_{0,f}}}\enspace {\mu }_g. $$(13)

The Knudsen diffusion coefficient is defined as,Dkf=πRTk0,fϕ0,fMw.$$ {D}_{\mathrm{kf}}=\sqrt{\frac{{\pi RT}{k}_{0,f}{\phi }_{0,f}}{{M}_w}}. $$(14)Here, S(pm,pf)$ S({p}_m,{p}_f)$ is the transfer parameter which is connecting of the matrix-fracture domains, and defined as,S(pm,pf)=σρg,mkmμ(pm-pf).$$ S\left({p}_m,{p}_f\right)=\frac{\sigma {\rho }_{g,m}{k}_m}{\mu }\left({p}_m-{p}_f\right). $$(15)

One may provide the definition of the source of the production well by,Q(pf)=θρg,fkfμlnrerw(pfe-pw),$$ Q({p}_f)=\frac{\theta {\rho }_{g,f}{k}_f}{\mu \mathrm{ln}\frac{{r}_e}{{r}_w}}\left({p}_{\mathrm{fe}}-{p}_w\right), $$(16)θ={θ=2πif the production well at the center of the reservoirθ=π2if the production well at the corner of the reservoir.$$ \theta =\left\{\begin{array}{l}\theta =2\pi \hspace{1em}\mathrm{if}\enspace \mathrm{the}\enspace \mathrm{production}\enspace \mathrm{well}\enspace \mathrm{at}\enspace \mathrm{the}\enspace \mathrm{center}\enspace \mathrm{of}\enspace \mathrm{the}\enspace \mathrm{reservoir}\\ \theta =\frac{\pi }{2}\hspace{1em}\mathrm{if}\enspace \mathrm{the}\enspace \mathrm{production}\enspace \mathrm{well}\enspace \mathrm{at}\enspace \mathrm{the}\enspace \mathrm{corner}\enspace \mathrm{of}\enspace \mathrm{the}\enspace \mathrm{reservoir}\end{array}.\right. $$

Given the constant, rc$ {r}_c$, the drainage-radius is represented as,re=0.14(Δx)2+(Δy)2.$$ {r}_e=0.14\sqrt{(\Delta x{)}^2+(\Delta y{)}^2}. $$(17)

The coefficient of crossflow between the matrix and fracture domains is defined as [9],σ=2n(n+1)l2,$$ \sigma =\frac{2n(n+1)}{{l}^2}, $$(18)such that l$ l$ is given by,l=3LxLyLzLxLy+LyLz+LxLz.$$ l=\frac{3{L}_x{L}_y{L}_z}{{L}_x{L}_y+{L}_y{L}_z+{L}_x{L}_z}. $$(19)

2.1 Geomechanical effects

Invoking the rock stress-sensitivity [27], the porosity can be expressed as a function of the mean effective-stress,ϕd(σ'd(pd))=ϕr,d+(ϕ0,d-ϕr,d)exp(-'d), d=m,f,$$ {{{\phi }_d({\sigma \prime}_d({p}_d))={\phi }_{r,d}+({\phi }_{0,d}-{\phi }_{r,d})\mathrm{exp}(-{a\sigma \prime}_d),\enspace \hspace{1em}d=m,f, $$(20)σ'd=σd+αdpd, d=m,f.$$ {{\sigma \prime}_d={\sigma }_d+{\alpha }_d{p}_d,\enspace \hspace{1em}d=m,f. $$(21)From (20) and (21), we can obtain,ϕd(pd)=ϕr,d+C1,dexp(-apd), C1,d=(ϕ0,d-ϕr,d)exp(-aσd), d=m,f.$$ {\phi }_d({p}_d)={\phi }_{r,d}+{C}_{1,d}\mathrm{exp}(-a{p}_d),\enspace {C}_{1,d}=({\phi }_{0,d}-{\phi }_{r,d})\mathrm{exp}(-a{\sigma }_d),\enspace \hspace{1em}d=m,f. $$(22)

As the porosity increases the difference (ϕ0,d-ϕr,d)$ ({\phi }_{0,d}-{\phi }_{r,d})$ becomes positive and vice versa. Also, the coefficient C1,d$ {C}_{1,d}$ becomes small based on the change in the porosity, and it becomes positive if porosity increases and negative if porosity decreases.

2.2 Initial and boundary conditions

The initial conditions can be represented by,pm(,0)=pf(,0)=p0in ΩmΩf.$$ {p}_m(\cdot,0)={p}_f(\cdot,0)={p}_0\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}\enspace {\mathrm{\Omega }}_m\cup {\mathrm{\Omega }}_f. $$(23)

The pressure on the production well (boundary condition) is given as,pf(,t)=pwonΓDf×(0,T).$$ {p}_f\left(\cdot,t\right)={p}_w\hspace{1em}\mathrm{o}\mathrm{n}\hspace{1em}{\mathrm{\Gamma }}_D^f\times \left(0,T\right). $$(24)

The no-flow boundary conditions on matrix are given as,umn=0onΓNmΓNf×(0,T),$$ {u}_m\cdot n=0\hspace{1em}\mathrm{o}\mathrm{n}\hspace{1em}{\mathrm{\Gamma }}_N^m\cup {\mathrm{\Gamma }}_N^f\times (0,T), $$(25)while the no-flow boundary conditions for fracture domain are,ufn=0onΓNf×(0,T).$$ {u}_f\cdot n=0\hspace{1em}\mathrm{o}\mathrm{n}\hspace{1em}{\mathrm{\Gamma }}_N^f\times \left(0,T\right). $$(26)

3 MFEM spaces

The MFEM contains two spaces for a scalar variable and it’s flux. The MFEM was developed to approximate both variables simultaneously and to give a higher order approximation for the flux. A compatibility condition must hold for the two spaces to insure stability, consistency and convergence of the mixed method and by considering more constraints to the numerical discretization.

Now, let us define the inner product in Ω$ \mathrm{\Omega }$ as,(f,g)Ω=Ωf(x)g(x)dV f,g:ΩR,$$ (f,g{)}_{\mathrm{\Omega }}={\int }_{\mathrm{\Omega }} f(\mathbf{x})g(\mathbf{x})\mathrm{d}\mathbf{V}\enspace \hspace{1em}\forall f,g:\mathrm{\Omega }\to \mathbb{R}, $$(27)and the inner product on Ω$ \mathrm{\partial \Omega }$ as,f,gΩ=ΩfgdS f,g:ΩR.$$ \left\langle f,g{\right\rangle_{\mathrm{\partial \Omega }}={\int }_{\mathrm{\partial \Omega }} {fg}\mathrm{d}S\enspace \hspace{1em}\forall f,g:\mathrm{\partial \Omega }\to \mathbb{R}. $$(28)Given,L2(Ω)={f:ΩR:Ωf2dx<+},$$ {\mathbb{L}}^2(\mathrm{\Omega })=\{f:\mathrm{\Omega }\to \mathbb{R}:{\int }_{\mathrm{\Omega }} {f}^2\mathrm{d}x < +\infty \}, $$(29)is the largest Hilbert space, such that (D-1u,w)$ ({D}^{-1}\mathbf{u},w)$ and (p,w)$ (p,\nabla \cdot w)$ are well defined, therefore, p,wL2(Ω)$ p,\nabla \cdot w\in {\mathbb{L}}^2(\mathrm{\Omega })$. It is needed that p,φL2(Ω)$ p,\phi \in {\mathbb{L}}^2(\mathrm{\Omega })$ and u,wH(div,Ω)$ \mathbf{u},w\in \mathbb{H}(\mathrm{div},\mathrm{\Omega })$. The Hilbert space with norm given by,H(div,Ω)={u:uL2(Ω)},$$ \mathbb{H}(\mathrm{div},\mathrm{\Omega })=\{\mathbf{u}:\nabla \cdot \mathbf{u}\in {\mathbb{L}}^2(\mathrm{\Omega })\}, $$(30)uH(div,Ω)2=u2 + u2.$$ {\Vert \mathbf{u}\Vert }_{\mathbf{H}(\mathrm{d}\mathrm{iv},\mathrm{\Omega })}^2={\Vert \mathbf{u}\Vert }^{\mathbf{2}}\enspace +\enspace {\Vert \nabla \cdot \mathbf{u}\Vert }^2. $$(31)

The seeking solution is,(p,u)L2(Ω)×H(div,Ω),$$ (p,\mathbf{u})\in {\mathbb{L}}^2(\mathrm{\Omega })\times \mathbb{H}(\mathrm{div},\mathrm{\Omega }), $$(32)such that p$ p$ and u$ \mathbf{u}$ are smooth. For the approximate numerical solution, the two spaces become,WhL2(Ω), dim(Wh)<+$$ {W}_h\subset {\mathbb{L}}^2(\mathrm{\Omega }),\enspace \hspace{1em}\mathrm{dim}({W}_h) < +\infty $$(33)and,VhH(div,Ω),dim(Vh)<+.$$ {V}_h\subset \mathbb{H}\left(\mathrm{div},\mathrm{\Omega }\right),\hspace{1em}\mathrm{dim}\left({V}_h\right) < +\infty. $$(34)

Thus, the normal component un$ \mathbf{u}\cdot \mathbf{n}$ is continuous across the inter-element boundaries.

Moreover, it is useful to present the RTr$ \mathbf{R}{\mathbf{T}}_r$ elements that are designed to approximate H(div,Ω)$ \mathbb{H}(\mathrm{div},\mathrm{\Omega })$ [7], which satisfies,Wh={wL2(Ω):w|EP0(E),EEh},$$ {W}_h=\left\{w\in {\mathbb{L}}^2\left(\mathrm{\Omega }\right):w{|}_E\in {\mathbb{P}}_0(E),E\in {\mathcal{E}}_h\right\}, $$(35)and,Vh={uH(div,Ω):u|EP1(E),EEh},$$ {V}_h=\left\{\mathbf{u}\in \mathbb{H}\left(\mathrm{div},\mathrm{\Omega }\right):\mathbf{u}{|}_E\in {\mathbb{P}}_1(E),E\in {\mathcal{E}}_h\right\}, $$(36)where w$ w$ is discontinuous piecewise constant and u$ \mathbf{u}$ is piecewise linear.

4 Mixed finite element approximation

Let Ωm$ {\mathrm{\Omega }}_m$ and Ωf$ {\mathrm{\Omega }}_f$ are, respectively, the matrix and fracture in a polygonal/polyhedral Lipschitz domain ΩRd,d{1,2,3}$ \mathrm{\Omega }\subset {\mathbb{R}}^d,\hspace{1em}d\in \{\mathrm{1,2},3\}$ (on which we define the standard L2(Ω)$ {L}^2(\mathrm{\Omega })$ space such that L2(Ω)(L2(Ω))d$ {\mathbb{L}}^2(\mathrm{\Omega })\equiv {\left({L}^2(\mathrm{\Omega })\right)}^d$), with the boundaries, Ωm=ΓDmΓNm$ \mathrm{\partial }{\mathrm{\Omega }}_m={\mathrm{\Gamma }}_D^m\cup {\mathrm{\Gamma }}_N^m$ and Ωf=ΓDfΓNf$ \mathrm{\partial }{\mathrm{\Omega }}_f={\mathrm{\Gamma }}_D^f\cup {\mathrm{\Gamma }}_N^f$. One may write the above dual-porosity model in the following general form,f1(pm)pmt+um=-S(pm,pf)inΩm×(0,T),$$ {f}_1\left({p}_m\right)\frac{\mathrm{\partial }{p}_m}{\mathrm{\partial }t}+\nabla \cdot {u}_m=-S\left({p}_m,{p}_f\right)\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_m\times (0,T), $$(37)Dm(pm)-1um=-pminΩm×(0,T),$$ {\mathbf{D}}_m({p}_m{)}^{-1}{u}_m=-\nabla {p}_m\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_m\times (0,T), $$(38)f2(pf)pft+uf=S(pm,pf)-Q(pf)inΩf×(0,T),$$ {f}_2({p}_f)\frac{\mathrm{\partial }{p}_f}{\mathrm{\partial }t}+\nabla \cdot {u}_f=S({p}_m,{p}_f)-Q({p}_f)\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_f\times (0,T), $$(39)Df(pf)-1uf=-pfinΩf×(0,T),$$ {\mathbf{D}}_f({p}_f{)}^{-1}{u}_f=-\nabla {p}_f\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_f\times (0,T), $$(40)whereDm(pm)=ρ(pm)μk0,m(1+bmpm),$$ {\mathbf{D}}_m({p}_m)=\frac{\rho ({p}_m)}{\mu }{k}_{0,m}\left(1+\frac{{b}_m}{{p}_m}\right), $$(41)andDf(pf)=ρ(pf)μk0,f(1+bfpf).$$ {\mathbf{D}}_f({p}_f)=\frac{\rho ({p}_f)}{\mu }{k}_{0,f}\left(1+\frac{{b}_f}{{p}_f}\right). $$(42)

The functions Dm(pm)-1$ {\mathbf{D}}_m({p}_m{)}^{-1}$ and Df(pf)-1$ {\mathbf{D}}_f({p}_f{)}^{-1}$ are moved to the left hand side to avoid discontinuity when we integrate pm$ \nabla {p}_m$ and pf$ \nabla {p}_f$ by parts. Selecting any φWh$ \phi \in {W}_h$ and ωVh$ \omega \in {V}_h$, the mixed finite element weak formulation can be written in the following form,(f1(pm)pmt,φ)+(um,φ)+(S(pm,pf),φ)=0,$$ \left({f}_1({p}_m)\frac{\mathrm{\partial }{p}_m}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_m,\phi )+(\mathcal{S}({p}_m,{p}_f),\phi )=0, $$(43)(Dm(pm)-1um,ω)=(pm,ω),$$ ({\mathbf{D}}_m({p}_m{)}^{-1}{\mathbf{u}}_m,\omega )=({p}_m,\nabla \cdot \omega ), $$(44)(f2(pf)pft,φ)+(uf,φ)-(S(pm,pf),φ)=-(Q(pf),φ),$$ \left({f}_2({p}_f)\frac{\mathrm{\partial }{p}_f}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_f,\phi )-(\mathcal{S}({p}_m,{p}_f),\phi )=-(Q({p}_f),\phi ), $$(45)(Df(pf)-1uf,ω)=(pf,ω)-pw,ωΓfD.$$ ({\mathbf{D}}_f({p}_f{)}^{-1}{\mathbf{u}}_f,\omega )=({p}_f,\nabla \cdot \omega )-\left\langle {p}_w,\omega \right\rangle_{{\mathrm{\Gamma }}_f^D}. $$(46)

Now, let the approximating subspace duality VhH(Ω;div)$ {V}_h\subset H(\mathrm{\Omega };\mathrm{div})$ and WhL2(Ω)$ {W}_h\subset {L}^2(\mathrm{\Omega })$ be the r$ r$-th order (r0$ r\ge 0$) Raviart–Thomas space (RTr) on the partition Th$ {\mathcal{T}}_h$. The mixed finite element formulations are stated as below: find pmh,pfhWh$ {p}_m^h,{p}_f^h\in {W}_h$ and umh,ufhVh$ {u}_m^h,{u}_f^h\in {V}_h$ such that,(f1(pmh)pmht,φ)+(umh,φ)+(S(pmh,pfh),φ)=0,$$ \left({f}_1({p}_m^h)\frac{\mathrm{\partial }{p}_m^h}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_m^h,\phi )+(\mathcal{S}({p}_m^h,{p}_f^h),\phi )=0, $$(47)(Dm(pmh)-1umh,ω)=(pmh,ω),$$ ({\mathbf{D}}_m({p}_m^h{)}^{-1}{\mathbf{u}}_m^h,\omega )=({p}_m^h,\nabla \cdot \omega ), $$(48)(f2(pfh)pfht,φ)+(ufh,φ)-(S(pmh,pfh),φ)=-(Q(pfh),φ),$$ \left({f}_2({p}_f^h)\frac{\mathrm{\partial }{p}_f^h}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_f^h,\phi )-(\mathcal{S}({p}_m^h,{p}_f^h),\phi )=-(Q({p}_f^h),\phi ), $$(49)(Df(pfh)-1ufh,ω)=(pfh,ω)-pw,ωΓfD,$$ ({\mathbf{D}}_f({p}_f^h{)}^{-1}{\mathbf{u}}_f^h,\omega )=({p}_f^h,\nabla \cdot \omega )-\left\langle {p}_w,\omega \right\rangle_{{\mathrm{\Gamma }}_f^D}, $$(50)for any φWh$ \phi \in {W}_h$ and ωVh$ \omega \in {V}_h$.

In order to obtain an explicit formulation for the flux, we employ a quadrature rule along with the MFE method to [8]. The total time interval [0,T]$ [0,T]$ is divided into NT$ {N}_T$ time steps with length Δtn=tn+1-tn$ \Delta {t}^n={t}^{n+1}-{t}^n$. The superscript n+1$ n+1$ denotes for the current time step, while n$ n$ denotes for the previous one. We use a backward Euler semi-implicit discretization for the time derivative terms. The following scheme has been developed,(f1(pmh,n)pmh,n+1-pmh,nΔt,φ)+(umh,n+1,φ)+(S(pmh,n+1,pfh,n),φ)=0,$$ \left({f}_1({p}_m^{h,n})\frac{{p}_m^{h,n+1}-{p}_m^{h,n}}{\Delta t},\phi \right)+(\nabla \cdot {\mathbf{u}}_m^{h,n+1},\phi )+(\mathcal{S}({p}_m^{h,n+1},{p}_f^{h,n}),\phi )=0, $$(51)(Dm(pmh,n)-1umh,n+1,ω)=(pmh,n+1,ω),$$ ({\mathbf{D}}_m({p}_m^{h,n}{)}^{-1}{\mathbf{u}}_m^{h,n+1},\omega )=({p}_m^{h,n+1},\nabla \cdot \omega ), $$(52)(f2(pfh,n)pfh,n+1-pfh,nΔt,φ)+(ufh,n+1,φ)-(S(pmh,n+1,pfh,n+1),φ)=-(Q(pfh,n+1),φ),$$ \left({f}_2({p}_f^{h,n})\frac{{p}_f^{h,n+1}-{p}_f^{h,n}}{\Delta t},\phi \right)+(\nabla \cdot {\mathbf{u}}_f^{h,n+1},\phi )-(\mathcal{S}({p}_m^{h,n+1},{p}_f^{h,n+1}),\phi )=-(Q({p}_f^{h,n+1}),\phi ), $$(53)(Df(pfh,n)-1ufh,n+1,ω)=(pfh,n+1,ω)-pw,ωΓfD.$$ ({\mathbf{D}}_f({p}_f^{h,n}{)}^{-1}{\mathbf{u}}_f^{h,n+1},\omega )=({p}_f^{h,n+1},\nabla \cdot \omega )-\left\langle {p}_w,\omega \right\rangle_{{\mathrm{\Gamma }}_f^D}. $$(54)Given pmh,n$ {p}_m^{h,n}$ and pfh,n$ {p}_f^{h,n}$, the numerical procedure to calculate pressure and velocity is presented here,

  1. Update the thermodynamical variables explicitly.

  2. Solve the equations (51) and (52) to obtain pmh,n+1$ {p}_m^{h,n+1}$ and umh,n+1$ {\mathbf{u}}_m^{h,n+1}$.

  3. Solve the equations (53) and (54) to obtain pfh,n+1$ {p}_f^{h,n+1}$ and ufh,n+1$ {\mathbf{u}}_f^{h,n+1}$.

5 Stability analysis

In this section, we carry out the stability analysis of the proposed MFE method, which ensures that the discrete solutions are bounded in a physically reasonable range. The key issue encountered in the stability analysis is the nonlinearity of the matrix and fracture pressures. In order to resolve this issue, we need to define some auxiliary functions and analyze their boundedness. Now let us define the following functions,F1(pdh)=0pdhf1(pdh)dpdh, G1(pdh)=0pdhF1(pdh)dpdh,$$ {F}_1({p}_d^h)={\int }_0^{{p}_d^h} {f}_1({p}_d^h)\mathrm{d}{p}_d^h,\enspace {G}_1({p}_d^h)={\int }_0^{{p}_d^h} {F}_1({p}_d^h)\mathrm{d}{p}_d^h, $$(55)F2(pdh)=0pdhf2(pdh)dpdh, G2(pdh)=0pdhF2(pdh)dpdh.$$ {F}_2({p}_d^h)={\int }_0^{{p}_d^h} {f}_2({p}_d^h)\mathrm{d}{p}_d^h,\enspace {G}_2({p}_d^h)={\int }_0^{{p}_d^h} {F}_2({p}_d^h)\mathrm{d}{p}_d^h. $$(56)

Therefore, one mayF1(pmh)=(ϕr,m+C1,m)MwZRTpmhe-apmh+MwVLρsVstd[C1,m(aeaPL(1+PL)Ei(-x)|aPLa(PL+pmh)+PLe-apmhPL+pmh-1)+(1-ϕr,m)pmhPL+pmh],$$ \begin{array}{l}{F}_1({p}_m^h)=\left({\phi }_{r,m}+{C}_{1,m}\right)\frac{{M}_w}{\mathrm{ZRT}}{p}_m^h{e}^{-a{p}_m^h}+\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}\left[{C}_{1,m}\left(a{e}^{a{P}_L}\left(1+{P}_L\right){Ei}(-x){|}_{a{P}_L}^{a({P}_L+{p}_m^h)}+\frac{{P}_L{e}^{-a{p}_m^h}}{{P}_L+{p}_m^h}-1\right)\right.\left.+(1-{\phi }_{r,m})\frac{{p}_m^h}{{P}_L+{p}_m^h}\right],\\ \end{array} $$(57)G1(pmh)=(ϕr,m+C1,m)Mwa2ZRT[1-(apmh+1)e-apmh]+MwVLρsVstd[C1,m(aeaPL(1+PL)((PL+pmh)Ei(-x)|aPLa(PL+pmh)-1a(1-e-apmh))PLeaPL+Ei(-x)|aPLa(PL+pmh))+(1-ϕr,m-C1,m)pmh+(1-ϕr,m)PLln(PLPL+pmh)],$$ \begin{array}{l}{G}_1({p}_m^h)=\left({\phi }_{r,m}+{C}_{1,m}\right)\frac{{M}_w}{{a}^2\mathrm{ZRT}}\left[1-(a{p}_m^h+1){e}^{-a{p}_m^h}\right]+\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}\left[{C}_{1,m}\left(a{e}^{a{P}_L}\left(1+{P}_L\right)\left(({P}_L+{p}_m^h){Ei}(-x){|}_{a{P}_L}^{a({P}_L+{p}_m^h)}-\frac{1}{a}(1-{e}^{-a{p}_m^h})\right){P}_L{e}^{a{P}_L}+{Ei}(-x){|}_{a{P}_L}^{a({P}_L+{p}_m^h)}\right)+(1-{\phi }_{r,m}-{C}_{1,m}){p}_m^h+(1-{\phi }_{r,m}){P}_L\mathrm{ln}\left(\frac{{P}_L}{{P}_L+{p}_m^h}\right)\right],\\ \\ \end{array} $$(58)F2(pfh)=(ϕr,f+C1,f)MwZRTpfhe-apfh,$$ {F}_2({p}_f^h)=\left({\mathrm{\phi }}_{r,f}+{C}_{1,f}\right)\frac{{M}_w}{\mathrm{ZRT}}{p}_f^h{e}^{-a{p}_f^h}, $$(59)andG2(pfh)=(ϕr,f+C1,f)Mwa2ZRT[1-(apfh+1)e-apfh].$$ {G}_2({p}_f^h)=\left({\phi }_{r,f}+{C}_{1,f}\right)\frac{{M}_w}{{a}^2\mathrm{ZRT}}\left[1-(a{p}_f^h+1){e}^{-a{p}_f^h}\right]. $$(60)

Ei(x)$ {Ei}(x)$ is a special function called the exponent integral function. As stated above in the model formulation that the coefficient C1,d,d=m,f$ {C}_{1,d},\hspace{1em}d=m,f$ is positive in the case of increasing porosity and negative in the case of decreasing porosity.

Lemma 5.1

For the case of increasing matrix porosity, C1,m>0$ {C}_{1,m}>0$, and sufficiently large pmh$ {p}_m^h$, there exists a positive constant,γ1=(ϕr,m+C1,m)C4lMwZRT+MwVLρsVstd(1-ϕr,m)C3l,$$ {\gamma }_1=\left({\phi }_{r,m}+{C}_{1,m}\right){C}_4^l\frac{{M}_w}{\mathrm{ZRT}}+\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}(1-{\phi }_{r,m}){C}_3^l, $$(61)such that,(F1(pmh),pmh)γ1pmh2.$$ ({F}_1({p}_m^h),{p}_m^h)\ge {\gamma }_1{\Vert {p}_m^h\Vert }^2. $$(62)

Proof. Note that the value of the exponent integral function Ei(-x)$ {Ei}(-x)$ is negative and (small for big values of x$ x$), therefore the quantity, Ei(-x)|aPLa(PL+pmh)$ {Ei}(-x){|}_{a{P}_L}^{a({P}_L+{p}_m^h)}$ has always a positive value which has lower and upper bounds, namely,C2,muEi(-x)|aPLa(PL+pmh)C2,ml,$$ {C}_{2,m}^u\ge {Ei}\left(-x\right){|}_{a{P}_L}^{a\left({P}_L+{p}_m^h\right)}\ge {C}_{2,m}^l, $$(63)where C2,ml,C2,mu>0$ {C}_{2,m}^l,{C}_{2,m}^u>0$. On the other hand, in shale reservoir, it is well known that the initial pressure has the maximum value, i.e., max(pmh)=P0$ \mathrm{max}({p}_m^h)={P}_0$, while, the pressure of the well has the minimum value, i.e., min(pmh)=Pw$ \mathrm{min}({p}_m^h)={P}_w$. Therefore, we have,C3u=1PL+pw1PL+pmh1PL+p0=C3l,$$ {C}_3^u=\frac{1}{{P}_L+{p}_w}\ge \frac{1}{{P}_L+{p}_m^h}\ge \frac{1}{{P}_L+{p}_0}={C}_3^l, $$(64)C4u=e-apwe-apmhe-ap0=C4l,$$ {C}_4^u={e}^{-a{p}_w}\ge {e}^{-a{p}_m^h}\ge {e}^{-a{p}_0}={C}_4^l, $$(65)and,C5u=e-apwPL+pwe-apmhPL+pmhe-ap0PL+p0=C5l,$$ {C}_5^u=\frac{{e}^{-a{p}_w}}{{P}_L+{p}_w}\ge \frac{{e}^{-a{p}_m^h}}{{P}_L+{p}_m^h}\ge \frac{{e}^{-a{p}_0}}{{P}_L+{p}_0}={C}_5^l, $$(66)where C3l,C3u,C4l,C4u,C5l,C5u0$ {C}_3^l,{C}_3^u,{C}_4^l,{C}_4^u,{C}_5^l,{C}_5^u\ge 0$. Therefore, for the case of increasing matrix porosity, C1,m>0$ {C}_{1,m}>0$, and sufficiently large pmh$ {p}_m^h$, and holding (63)(66), the following coefficient is positive, i.e.,aeaPL(1+PL)Ei(-x)|aPLa(PL+pmh)+PLe-apmhPL+pmh-1>,$$ a{e}^{a{P}_L}\left(1+{P}_L\right){Ei}\left(-x\right){|}_{a{P}_L}^{a\left({P}_L+{p}_m^h\right)}+\frac{{P}_L{e}^{-a{p}_m^h}}{{P}_L+{p}_m^h}-1>, $$(67)aeaPL(1+PL)C2,ml+PLC5l-1=C3,m>0.$$ a{e}^{a{P}_L}\left(1+{P}_L\right){C}_{2,m}^l+{P}_L{C}_5^l-1={C}_{3,m}>0. $$(68)Therefore,(F1(pmh),pmh)(ϕr,m+C1,m)C4lMwZRT(pmh,pmh)+MwVLρsVstd[C1,mC3,m+(1-ϕr,m)C3l(pmh,pmh)](ϕr,m+C1,m)C4lMwZRT(pmh,pmh)+MwVLρsVstd(1-ϕr,m)C3l(pmh,pmh)=[(ϕr,m+C1,m)C4lMwZRT+MwVLρsVstd(1-ϕr,m)C3l]pmh2=γ1pmh2,$$ \begin{array}{cc}({F}_1({p}_m^h),{p}_m^h)& \ge ({\phi }_{r,m}+{C}_{1,m}){C}_4^l\frac{{M}_w}{\mathrm{ZRT}}({p}_m^h,{p}_m^h)+\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}\left[{C}_{1,m}{C}_{3,m}+(1-{\phi }_{r,m}){C}_3^l({p}_m^h,{p}_m^h)\right]\\ \begin{array}{c}\\ \end{array}& \begin{array}{c}\ge \left({\phi }_{r,m}+{C}_{1,m}\right){C}_4^l\frac{{M}_w}{\mathrm{ZRT}}({p}_m^h,{p}_m^h)+\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}(1-{\phi }_{r,m}){C}_3^l({p}_m^h,{p}_m^h)\\ =\left[\left({\phi }_{r,m}+{C}_{1,m}\right){C}_4^l\frac{{M}_w}{\mathrm{ZRT}}+\frac{{M}_w{V}_L{\rho }_s}{{V}_{{std}}}(1-{\mathrm{\phi }}_{r,m}){C}_3^l\right]{\Vert {p}_m^h\Vert }^2={\gamma }_1{\Vert {p}_m^h\Vert }^2,\end{array}\end{array} $$(69)and this completes Lemma 5.1 proof.

Lemma 5.2

For the case of decreasing matrix porosity, C1,m<0$ {C}_{1,m} < 0$, assuming,aeaPL(1+PL)Ei(-x)|aPLa(PL+pmh)+PLe-apmhPL+pmh-1<0,$$ a{e}^{a{P}_L}\left(1+{P}_L\right){Ei}(-x){|}_{a{P}_L}^{a\left({P}_L+{p}_m^h\right)}+\frac{{P}_L{e}^{-a{p}_m^h}}{{P}_L+{p}_m^h}-1 < 0, $$(70)there exists a positive constant γ1$ {\gamma }_1$, such that,(F1(pmh),pmh)γ1pmh2.$$ ({F}_1({p}_m^h),{p}_m^h)\ge {\gamma }_1{\Vert {p}_m^h\Vert }^2. $$(71)

Proof. It is clear that, C1,m<ϕr,m$ {C}_{1,m} < {\phi }_{r,m}$, so, we always have, ϕr,m-C1,m>0$ {\phi }_{r,m}-{C}_{1,m}>0$. Holding the assumption (70), one can easily prove this lemma in a similar way to Lemma 5.1.

Lemma 5.3

For the case of increasing matrix porosity, C1,m>0$ {C}_{1,m}>0$, and for sufficiently large pmh$ {p}_m^h$, there exist two positive constants,γ2=Mwa2ZRT(ϕr,m+C1,m)(1-C4,mu)+MwVLρsVstd[C1,m(aeaPL(1+PL)PLC2,mu+PLeaPLC2,mu)+(1-ϕr,m)PLln(PLC3u)],$$ \begin{array}{l}{\gamma }_2=\frac{{M}_w}{{a}^2\mathrm{ZRT}}\left({\phi }_{r,m}+{C}_{1,m}\right)\left(1-{C}_{4,m}^u\right)+\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}\left[{C}_{1,m}\left(a{e}^{a{P}_L}\left(1+{P}_L\right){P}_L{C}_{2,m}^u+\left.\left.{P}_L{e}^{a{P}_L}{C}_{2,m}^u\right)+\left(1-{\phi }_{r,m}\right){P}_L\mathrm{ln}\left({P}_L{C}_3^u\right)\right]\right.\right.,\\ \end{array} $$(72)andγ'2=MwVLρsVstd[C1,maeaPL(1+PL)C2,mu+(1-ϕr,m-C1,m)],$$ {\gamma \prime}_2^{}=\frac{{M}_w{V}_L{\rho }_s}{{V}_{\mathrm{std}}}\left[{C}_{1,m}a{e}^{a{P}_L}(1+{P}_L){C}_{2,m}^u+(1-{\phi }_{r,m}-{C}_{1,m})\right], $$(73)such that,(G1(pmh),1)γ2|Ω|+γ'2pmh.$$ \left({G}_1\left({p}_m^h\right),1\right)\le {\gamma }_2\left|\mathrm{\Omega }\right|+{{\gamma }^{\prime}_2^{}\Vert {p}_m^h\Vert. $$(74)

Proof. Holding (63)(66) with C4,mu=apmh+1eapmh<1$ {C}_{4,m}^u=\frac{a{p}_m^h+1}{{e}^{a{p}_m^h}} < 1$ for sufficiently large pmh$ {p}_m^h$, and integrate (58) over Ωm$ {\mathrm{\Omega }}_m$ for the case of increasing porosity, C1,m>0$ {C}_{1,m}>0$, it is easy to prove (74), which completes the lemma proof.

Lemma 5.4

For the case of decreasing matrix porosity, C1,m<0$ {C}_{1,m} < 0$, with assuming that,C1,m(aeaPL(1+PL)PLC2,mu+PLeaPLC2,mu)+(1-ϕr,m)PLln(PLC3u)0,$$ {C}_{1,m}\left(a{e}^{a{P}_L}(1+{P}_L){P}_L{C}_{2,m}^u+{P}_L{e}^{a{P}_L}{C}_{2,m}^u\right)+(1-{\phi }_{r,m}){P}_L\mathrm{ln}({P}_L{C}_3^u)\ge 0, $$(75)andC1,maeaPL(1+PL)C2,mu+(1-ϕr,m-C1,m)0,$$ {C}_{1,m}a{e}^{a{P}_L}(1+{P}_L){C}_{2,m}^u+(1-{\phi }_{r,m}-{C}_{1,m})\ge 0, $$(76)one can prove,(G1(pmh),1)γ2|Ω|+γ'2pmh.$$ ({G}_1({p}_m^h),1)\le {\gamma }_2|\mathrm{\Omega }|+{\gamma \prime}_2^{}\Vert {p}_m^h\Vert. $$(77)

Proof. It is clear that, C1,m<ϕr,m$ {C}_{1,m} < {\phi }_{r,m}$, so,we always have, ϕr,m-C1,m>0$ {\phi }_{r,m}-{C}_{1,m}>0$. Holding the assumptions (75) and (76), therefore we find, γ2>0,γ'2>0$ {\gamma }_2>0,\hspace{1em}{\gamma \prime}_2^{}>0$, and integrate (58) over Ωm$ {\mathrm{\Omega }}_m$ for the case of decreasing porosity, C1,m<0$ {C}_{1,m} < 0$, which proves (74), and this completes the proof.

Lemma 5.5

For both cases of increasing or decreasing fracture porosity, namely, C1,f>0$ {C}_{1,f}>0$ or C1,f<0$ {C}_{1,f} < 0$, and sufficiently large pfh$ {p}_f^h$, there exists a positive constant,γ3=(ϕr,f+C1,f)MwZRTC4l,$$ {\gamma }_3=\left({\phi }_{r,f}+{C}_{1,f}\right)\frac{{M}_w}{\mathrm{ZRT}}{C}_4^l, $$(78)such that,(F2(pfh),pfh)γ3pfh2.$$ ({F}_2({p}_f^h),{p}_f^h)\ge {\gamma }_3{\Vert {p}_f^h\Vert }^2. $$(79)

Proof. Again, we have,C4u=e-apwe-apfhe-ap0=C4l.$$ {C}_4^u={e}^{-a{p}_w}\ge {e}^{-a{p}_f^h}\ge {e}^{-a{p}_0}={C}_4^l. $$(80)where C4l,C4u>0$ {C}_4^l,{C}_4^u>0$. Therefore, for the case of increasing fracture porosity, C1,f>0$ {C}_{1,f}>0$, and sufficiently large pfh$ {p}_f^h$, Therefore,(F2(pfh),pfh)(ϕr,f+C1,f)MwZRTC4l(pfh,pfh)=γ3pfh2.$$ ({F}_2({p}_f^h),{p}_f^h)\ge \left({\phi }_{r,f}+{C}_{1,f}\right)\frac{{M}_w}{\mathrm{ZRT}}{C}_4^l({p}_f^h,{p}_f^h)={\gamma }_3{\Vert {p}_f^h\Vert }^2. $$(81)

As stated above, for the case of the decreasing porosity, C1,f<0$ {C}_{1,f} < 0$, the coefficient (ϕr,f+C1,f)$ ({\phi }_{r,f}+{C}_{1,f})$ remains positive, then the above inequality holds. This completes the lemma proof.

Lemma 5.6

For both cases of increasing or decreasing fracture porosity, namely, C1,f>0$ {C}_{1,f}>0$ or C1,f<0$ {C}_{1,f} < 0$, and sufficiently large pfh$ {p}_f^h$, there exists a positive constant,γ4=Mwa2ZRT(ϕr,m+C1,m)(1-C4,mu),$$ {\gamma }_4=\frac{{M}_w}{{a}^2\mathrm{ZRT}}({\phi }_{r,m}+{C}_{1,m})(1-{C}_{4,m}^u), $$(82)such that,(G2(pfh),1)γ4|Ω|.$$ ({G}_2({p}_f^h),1)\le {\gamma }_4|\mathrm{\Omega }|. $$(83)

Proof. Similar to the previous proofs.

Lemma 5.7

Assuming that,0<ρg*ρgρg*,0<kf*kfkf*,0<μg*μgμg*,$$ 0<{\rho }_{g\mathrm{*}}\le {\rho }_g\le {\rho }_g^{\mathrm{*}},\hspace{1em}0<{k}_{f\mathrm{*}}\le {k}_f\le {k}_f^{\mathrm{*}},\hspace{1em}0 < {\mu }_{g\mathrm{*}}\le {\mu }_g\le {\mu }_g^{\mathrm{*}}, $$(84)then, Q(pfh)$ Q({p}_f^h)$ is bounded.

Proof. Since pfe$ {p}_{\mathrm{fe}}$ is bounded by pw$ {p}_w$ and p0$ {p}_0$, i.e., pwpfep0$ {p}_w\le {p}_{\mathrm{fe}}\le {p}_0$, then, 0<pfe-pwp0-pw=Cp$ 0 < {p}_{\mathrm{fe}}-{p}_w\le {p}_0-{p}_w={C}_p$, such that Cp$ {C}_p$ is positive number. Holding the conditions (84), one may find that,0<kf(pfh)ρg(pfh)[pfe-pw]μgkf*ρg*Cpμg*=CQ>0.$$ 0 < \frac{{k}_f({p}_f^h){\rho }_g({p}_f^h)[{p}_{\mathrm{fe}}-{p}_w]}{{\mu }_g}\le \frac{{k}_f^{\mathrm{*}}{\rho }_g^{\mathrm{*}}{C}_p}{{\mu }_g^{\mathrm{*}}}={C}_Q>0. $$(85)

Theorem 8

Holding the results in the above lemmas ( 5.1 5.6 ), namely, ( F 1 ( p m h ) , p m h ) γ 1 p m h 2 , ( G 1 ( p m h ) , 1 ) γ 2 | Ω | + γ ' 2 p m h , $$ \left({F}_1\left({p}_m^h\right),{p}_m^h\right)\ge {\gamma }_1{\Vert {p}_m^h\Vert }^2,\hspace{1em}\left({G}_1\left({p}_m^h\right),1\right)\le {\gamma }_2\left|\mathrm{\Omega }\right|+{\gamma \prime}_2^{}\Vert {p}_m^h\Vert, $$(86) and ( F 2 ( p f h ) , p f h ) γ 3 p f h 2 , ( G 2 ( p f h ) , 1 ) γ 4 1 , $$ \left({F}_2\left({p}_f^h\right),{p}_f^h\right)\ge {\gamma }_3{\Vert {p}_f^h\Vert }^2,\hspace{1em}\left({G}_2\left({p}_f^h\right),1\right)\le {\gamma }_4\Vert 1\Vert, $$(87) with assuming that, γ 1 p m h 2 γ 2 | Ω | + γ ' 2 p m h , $$ {\gamma }_1{\Vert {p}_m^h\Vert }^2\ge {\gamma }_2|\mathrm{\Omega }|+{\gamma \prime}_2^{}\Vert {p}_m^h\Vert, $$(88) and γ 3 p f h 2 γ 4 | Ω | . $$ {\gamma }_3{\Vert {p}_f^h\Vert }^2\ge {\gamma }_4|\mathrm{\Omega }|. $$(89)

Moreover, suppose that pfe,pwL2(0,T;L2(Ω))$ {p}_{\mathrm{fe}},{p}_w\in {L}^2(0,T;{L}^2(\Omega ))$. Then,pmhL(0,T;L2(Ω))2+Dm(pmh)-1/2umhL2(0,T;L2(Ω))2+pfhL(0,T;L2(Ω))2+Df(pfh)-1/2ufhL2(0,T;L2(Ω))2+pmh-pfhL2(0,T;L2(Ω))2C(F1(p0),p0)+C(F2(p0),p0)+C(G1(p0),1)+C(G2(p0),1)+C(pfeL2(0,T;L2(Ω))2+pwL2(0,T;L2(Ω))2).$$ \begin{array}{l}{\Vert {p}_m^h\Vert }_{{L}^{\mathrm{\infty }}(0,T;{L}^2(\mathrm{\Omega }))}^2+{\Vert {\mathbf{D}}_m({p}_m^h{)}^{-1/2}{u}_m^h\Vert }_{{L}^2(0,T;{L}^2(\mathrm{\Omega }))}^2+{\Vert {p}_f^h\Vert }_{{L}^{\mathrm{\infty }}(0,T;{L}^2(\mathrm{\Omega }))}^2+\\ {\Vert {\mathbf{D}}_f({p}_f^h{)}^{-1/2}{u}_f^h\Vert }_{{L}^2(0,T;{L}^2(\mathrm{\Omega }))}^2+{\Vert {p}_m^h-{p}_f^h\Vert }_{{L}^2(0,T;{L}^2(\mathrm{\Omega }))}^2\le C({F}_1({p}_0),{p}_0)+C({F}_2({p}_0),{p}_0)\\ +C({G}_1({p}_0),1)+C({G}_2({p}_0),1)+C({\Vert {p}_{\mathrm{fe}}\Vert }_{{L}^2(0,T;{L}^2(\mathrm{\Omega }))}^2+{\Vert {p}_w\Vert }_{{L}^2(0,T;{L}^2(\mathrm{\Omega }))}^2).\end{array} $$(90)

Proof. Let φ=pmh$ \phi ={p}_m^h$ in (47) and ω=umh$ \omega ={\mathbf{u}}_m^h$ in (48), and sum the two equations, we get,(f1(pmh)pmht,pmh)+Dm(pmh)-1/2umh2+(S(pmh,pfh),pmh)=0.$$ \left({f}_1({p}_m^h)\frac{\mathrm{\partial }{p}_m^h}{\mathrm{\partial }t},{p}_m^h\right)+{\Vert {D}_m({p}_m^h{)}^{-1/2}{\mathbf{u}}_m^h\Vert }^2+(\mathcal{S}({p}_m^h,{p}_f^h),{p}_m^h)=0. $$(91)

It follows from the definitions of F1$ {F}_1$ and G1$ {G}_1$ that,(f1(pmh)pmht,pmh)=t(F1(pmh),pmh)-(F1(pmh),pmht)=t(F1(pmh),pmh)-(G1(pmh)t,1).$$ \left({f}_1({p}_m^h)\frac{\mathrm{\partial }{p}_m^h}{\mathrm{\partial }t},{p}_m^h\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }t}({F}_1({p}_m^h),{p}_m^h)-\left({F}_1({p}_m^h),\frac{\mathrm{\partial }{p}_m^h}{\mathrm{\partial }t}\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }t}({F}_1({p}_m^h),{p}_m^h)-\left(\frac{\mathrm{\partial }{G}_1({p}_m^h)}{\mathrm{\partial }t},1\right). $$(92)

Integrating (92) from 0 to t$ t$ (0<tT$ 0 < t\le T$) yields,γ1pmh2(t)-γ2|Ω|(t)-γ'2pmh(t)+0tDm(pmh)-1/2umh2+0t(S(pmh,pfh),pmh)(F1(p0),p0)-(G1(p0),1).$$ {\gamma }_1{\Vert {p}_m^h\Vert }^2(t)-{\gamma }_2|\mathrm{\Omega }|(t)-{\gamma \prime}_2^{}\Vert {p}_m^h\Vert (t)+{\int }_0^t {\Vert {\mathbf{D}}_m({p}_m^h{)}^{-1/2}{\mathbf{u}}_m^h\Vert }^2+{\int }_0^t (\mathcal{S}({p}_m^h,{p}_f^h),{p}_m^h)\le ({F}_1({p}_0),{p}_0)-({G}_1({p}_0),1). $$(93)

It is similar to obtain,γ3pfh2(t)-γ4|Ω|(t)+0tDf(pfh)-1/2ufh2+0tpw,uwΓfD-0t(S(pmh,pfh),pfh)(F2(p0),p0)-(G2(p0),1)-0t(Q(pfh),pfh).$$ {\gamma }_3{\Vert {p}_f^h\Vert }^2(t)-{\gamma }_4|\mathrm{\Omega }|(t)+{\int }_0^t {\Vert {\mathbf{D}}_f({p}_f^h{)}^{-1/2}{\mathbf{u}}_f^h\Vert }^2+{\int }_0^t \left\langle {p}_w,{\mathbf{u}}_w\right\rangle_{{\mathrm{\Gamma }}_f^D}-{\int }_0^t (\mathcal{S}({p}_m^h,{p}_f^h),{p}_f^h)\le ({F}_2({p}_0),{p}_0)-({G}_2({p}_0),1)-{\int }_0^t (Q({p}_f^h),{p}_f^h). $$(94)

There exists a positive constant α0$ {\alpha }_0$ such that,(S(pmh,pfh),pmh)-(S(pmh,pfh),pfh)α0pmh-pfh2.$$ \left(\mathcal{S}\left({p}_m^h,{p}_f^h\right),{p}_m^h\right)-\left(\mathcal{S}\left({p}_m^h,{p}_f^h\right),{p}_f^h\right)\ge {\alpha }_0{\Vert {p}_m^h-{p}_f^h\Vert }^2. $$(95)

Thus, it is obtained from, (93)(95) that,γ1pmh2(t)-γ2|Ω|(t)-γ'2pmh(t)+0tDm(pmh)-1/2umh2+γ3pfh2(t)-γ4|Ω|(t)+0tDf(pfh)-1/2ufh2+α00tpmh-pfh2(F1(p0),p0)+(F2(p0),p0)-(G1(p0),1)-(G2(p0),1)+C0t(pfe2+pw2+pfh2).$$ {{\gamma }_1{\Vert {p}_m^h\Vert }^2(t)-{\gamma }_2|\mathrm{\Omega }|(t)-\gamma \prime}_2\Vert {p}_m^h\Vert (t)+{\int }_0^t {\Vert {\mathbf{D}}_m({p}_m^h{)}^{-1/2}{u}_m^h\Vert }^2+{\gamma }_3{\Vert {p}_f^h\Vert }^2(t)-{\gamma }_4|\mathrm{\Omega }|(t)+{\int }_0^t {\Vert {\mathbf{D}}_f({p}_f^h{)}^{-1/2}{u}_f^h\Vert }^2+{\alpha }_0{\int }_0^t{ \Vert {p}_m^h-{p}_f^h\Vert }^2\le ({F}_1({p}_0),{p}_0)+({F}_2({p}_0),{p}_0)-({G}_1({p}_0),1)-({G}_2({p}_0),1)+C{\int }_0^t ({\Vert {p}_{\mathrm{fe}}\Vert }^2+{\Vert {p}_w\Vert }^2+{\Vert {p}_f^h\Vert }^2). $$(96)Finally, (90) is obtained by 1Gronwall’s lemma.

Remark 5.1. In Lemma 5.6, if we do not use the coefficient, C4,mu=apmh+1eapmh<1$ {C}_{4,m}^u=\frac{a{p}_m^h+1}{{e}^{a{p}_m^h}} < 1$, and just use the coefficient, Cmu=e-apw$ {C}_m^u={e}^{-a{p}_w}$, then the coefficient of the term, apmh$ a{p}_m^h$ in C4,mu$ {C}_{4,m}^u$ will be moved from γ2$ {\gamma }_2$ into γ'2$ {\gamma \prime}_2^{}$. Thus, the definitions of both γ2$ {\gamma }_2$ and γ'2$ {\gamma \prime}_2^{}$ will be slightly changed without loss of the stability concept.

6 Conclusion

In this work, stability analysis of the MFE solution of the gas transport in a low-permeability fractured reservoir with stress-sensitivity effect has been considered. The dual-porosity dual permeability model with the slippage effect and the apparent permeability have been used to describe the flow in a low-permeability fractured porous media. Stability analysis of the MFEM is presented theoretically and numerically. We proved a seven lemmas and a theorem on the stability of the MFEM. Stability conditions are stated and estimated. Variation in the porosity is correlated to the stress-sensitivity effect and depends on the values of the corresponding physical parameters. For example, the coefficient C1,d$ {C}_{1,d}$ is relatively small depending on the variation of the porosity. It has a positive value in the case of porosity increasing while it has a negative value in the case of porosity decreasing. Lemmas 5.1 and 5.3 discuss the positive case, while Lemmas 5.2 and 5.4 discuss the positive case. Lemmas 5.5 and 5.6 discuss the two cases.

References

  1. Zienkiewicz O.C., Taylor R.L. (2000) The finite element method (5th edn.), vol. 1 – The basis, Butterworth-Heinemann, Oxford.
  2. Bathe K.J. (1996) Finite element procedures, Prentice Hall, Englewood Cliffs, NJ.
  3. Zienkiewicz O.C., Taylor R.L. (2000) The finite element method (5th edn), vol. 3 – fluid dynamics, Butterworth-Heinemann, Oxford.
  4. Brezzi F., Fortin V. (1991) Mixed and hybrid finite element methods, Springer-Verlag, New York.
  5. Chen Z. (2010) Finite element methods and their applications, Ch. 3. Springer-Verlag Berlin and Heidelberg GmbH & Co. KG.
  6. Brezzi F., Douglas J., Marini L.D. (1985) Two families of mixed finite elements for second order elliptic problems, Numer. Math. 47, 217–235.
  7. Raviart P.A., Thomas J.M. (1977) A mixed finite element method for 2nd order elliptic problems. Mathematical aspects of finite element methods, in: Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975, Lect. Notes Math., Vol. 606, Springer, Berlin, pp. 292–315.
  8. Arbogast T., Wheeler M.F., Yotov I. (1997) Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences, SIAM J. Num. Anal. 34, 2, 828–852.
  9. Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, Soc. Petrol. Eng. J. 3, 245–255.
  10. Ozkan E., Raghavan R., Apaydin O. (2010) Modeling of fluid transfer from shale matrix to fracture network, in: Annual Technical Conference and Exhibition, Lima, Peru, 1–3 December. SPE-134830-MS.
  11. Ertekin T., King G.R., Schwerer F.C. (1986) Dynamic gas slippage: A unique dual-mechanism approach to the flow of gas in tight formations, SPE Form. Evalu. 1,1, 43–52.
  12. Bustin A., Bustin R., Cui X. (2008) Importance of fabric on the production of gas shales, in: Unconventional Reservoirs Conference in Colorado, USA, 10–12 February. SPE-114167-MS.
  13. Moridis G., Blasingame T., Freeman C. (2010) Analysis of mechanisms of flow in fractured tight-gas and shale-gas reservoirs, in: Latin American and Caribbean Petroleum Engineering Conference, 1–3 December. SPE-139250-MS.
  14. Shu W.Y., Fakcharoenphol P. (2011) A unified mathematical model for unconventional reservoir simulation, in: EUROPEC/EAGE Annual Conference and Exhibition in Vienna, Austria, 23–26 May. SPE-142884-MS.
  15. Kazemi H. (1969) Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution, Soc. Pet. Eng. 9, 4, 451–462.
  16. Guo C, Wei M, Chen H, He X, Bai B (2014) Improved numerical simulation for shale gas reservoirs, in: Offshore Technology Conference-Asia, 25–28 March 2014, Kuala Lumpur, Malaysia.
  17. Javadpour F. (2009) Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone), J. Can. Pet. Technol. 48, 8, 16–21.
  18. Wang K., Liu H., Luo J., Wu K., Chen Z. (2017) A comprehensive model coupling embedded discrete fractures, multiple interacting continua and geomechanics in shale gas reservoirs with multi-scale fractures, Energy Fuel 31, 7758–7776. doi: 10.1021/acs.energyfuels.7b00394.
  19. Lin M., Chen S., Mbia S., Chen Z. (2018) Application of reservoir flow simulation integrated with geomechanics in unconventional tight play, Rock Mech. Rock Eng. 51, 315–328. doi: 10.1007/s00603-017-1304-1.
  20. Yang S., Wu K., Xu J., Chen Z. (2018) Roles of multicomponent adsorption and geomechanics in the development of an Eagle Ford shale condensate reservoir, Fuel 242, 710–718. doi: 10.1016/j.fuel.2019.01.016.
  21. Girault V., Wheeler M.F., Almani T., Dana S. (2019) A priori error estimates for a discretized poro-elastic–elastic system solved by a fixed-stress algorithm, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 74, 24. doi: 10.2516/ogst/2018071.
  22. El-Amin M.F., Kou J., Sun S. (2018) Mixed finite element simulation with stability analysis for gas transport in low-permeability reservoirs, Energies 111, 208–226.
  23. El-Amin M.F., Kou J., Sun S. (2018) Numerical modeling and simulation of shale-gas transport with geomechanical effect, Trans. Porous Med. 126, 779–806. doi: 10.1007/s11242-018-1206-z.
  24. Cui X., Bustin A.M.M., Bustin R.M. (2009) Measurements of gas permeability and diffusivity of tight reservoir rocks: Different approaches and their applications, Geofluids 9, 3, 208–223.
  25. Civan F., Rai C.S., Sondergeld C.H. (2011) Shale-gas permeability and diffusivity inferred by improved formulation of relevant retention and transport mechanisms, Trans. Porous Med. 86, 3, 925–944.
  26. Firoozabadi A. (2015) Thermodynamics and applications in hydrocarbon reservoirs and production. McGraw-Hill Education – Europe, United States.
  27. Terzaghi K. (1936) The shearing resistance of saturated soils and the angle between the planes of shear, in: Proceedings of International Conference on Soil Mechanics and Foundation Engineering, Harvard University Press, Cambridge, MA., Vol. 1, pp. 54–56.

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.