# A Mixed Element Method for the Desorption-Diffusion-Seepage Model of Gas Flow in Deformable Coalbed Methane Reservoirs

Authors: Lei Yang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We present a desorption-diffusion-seepage model for the gas flow problem in deformable coalbed methane reservoirs. Effects of fracture systems deformation on permeability have been considered in the proposed model. A mixed finite element method is introduced to solve the gas flow model, in which the coalbed gas pressure and velocity can be approximated simultaneously. Numerical experiments using the lowest order Raviart-Thomas (RT0) mixed element are carried out to describe the dynamic characteristics of gas pressure, velocity, and concentration. Error estimate results indicate that approximation solutions
could achieve first-order convergence rates.


---

## Body

## 1. Introduction

Coalbed methane (CBM) is an abundant, low cost energy source, which has become a viable alternative source to conventional fuel. CBM reservoirs are dual porosity systems consisting of coal matrix and fracture network [1]. In CBM reservoirs, most of coalbed gas is stored in the coal matrix as adsorbed gas, and only a small amount is stored in fracture systems as free gas. Comparing with conventional natural gas reservoirs, the effect of fracture systems deformation on permeability is significant in CBM reservoirs [2–4]. Meanwhile, the gas flow in fracture systems is driven not only by pressure field but also by gas concentration field [5, 6], so we need to consider Darcy seepage and Fick diffusion simultaneously in multifield. As coalbed gas is released from CBM reservoirs, fracture system pressure reduces, resulting in gas diffusion in the coal matrix and gas desorption from coal matrix surface to fracture systems. Therefore, the gas flow process in CBM reservoirs includes gas desorption-diffusion in coal matrix and gas seepage in fracture systems [7–10]. Usually, the fracture systems are considered as the gas flow passage and the coal matrix is treated as the source for fracture systems.Based on the above understandings about coalbed gas flow mechanism, a series of conventional mathematical and numerical models [5, 6, 11, 12] have been developed, obtaining some useful computational and simulation results. Young [11] presented a computer model for CBM reservoirs to determine key data and to describe the variability in CBM reservoir properties. Unsal et al. [12] proposed a numerical model for multiphase flow in CBM reservoirs using a fracture-only model. Thararoop et al. [5, 6] developed a compositional dual porosity, dual permeability CBM numerical model to describe the pressure and concentration dynamic behaviors. Nevertheless, the existing CBM numerical simulation methods are based on finite difference scheme and mainly focus on gas pressure and concentration, ignoring gas flow velocity.Mixed finite element method was initially introduced by engineers [13] in the 1960s and has been applied to many areas such as solid and fluid mechanics, which could approximate both vector variable (e.g., the fluid velocity) and scalar variable (e.g., the pressure) simultaneously and give a high order approximation of both variables. Comparing with the standard finite element method only employing a single finite element space, mixed finite element method need construct two different finite element spaces. Raviart and Thomas [14] introduced the first family of mixed finite element spaces for second-order elliptic problems. Nedelec [15] extended these spaces to three-dimensional problems. Then, Brezzi et al. [16, 17] and Chen and Douglas [18] presented many mixed finite element spaces.In addition, the study about numerical method and numerical analysis for fluid flow problems in porous media has been a research hotspot for the past decades [19–24], due to its wide applications in various engineering areas. Douglas et al. [19–21] have done a lot of useful work for the fluid flow problem in porous media. For the incompressible miscible displacement problem in porous media, Wang [22] introduced ELLAM-MFEM to solve equation of the problem. A component-based Eulerian-Lagrangian method [23] was used to treat the multicomponent, multiphase flow problem in porous media. Li and Sun [24] presented an unconditional convergence Galerkin-mixed element approximation for the incompressible miscible flow problem. However, until now, we have not found any paper considering mixed element approximation scheme for the gas flow problem in deformable CBM reservoirs.The aim of this paper is to study the gas flow problem in deformable CBM reservoirs using mixed element method. For this purpose, we introduce a mixed finite element method to approximate the coalbed gas pressure and velocity simultaneously. Another objective of this paper is to carry out some numerical experiments using the lowest order Raviart-Thomas mixed element, so as to show the convergence rate of the mixed element approximation method, and to analyze coalbed gas flow characteristics.This paper is organized as follows. In Section2, considering the effects of fracture system deformation, we present a desorption-diffusion-seepage model of gas flow in deformable CBM reservoirs for the two-dimensional problem. In Section 3 we introduce a mixed finite element approximation scheme for the coalbed gas flow model. In Section 4 numerical experiments using the lowest order Raviart-Thomas (RT0) mixed element are carried out.


## 2. Coalbed Gas Flow Problem in Deformable CBM Reservoirs

In this section, considering the pressure dependence of permeability and porosity, we derive a mathematical model for the gas flow problem in deformable CBM reservoirs.


### 2.1. Basic Assumptions and Gas State Descriptions

(i)
CBM reservoirs are treated as dual porosity systems consisting of coal matrix and fracture network.
(ii)
CBM reservoirs are isothermal. When pressure variation is not significant, we can assume the gas viscosity is constant under isothermal conditions [25].
(iii)
The gas diffusion process in the coal matrix is pseudostatic and described by Fick’s first diffusion law.
(iv)
Gas absorption/adsorption is described by Langmuir adsorption isotherm.The equation of real gas state is(1)pV=nZRT=mMZRT,
where p is the gas pressure, V is the gas volume, n is the amount of substance, T is the gas temperature, m is the gas mass, M is the gas molar mass, Z is the gas deviation factor, and R is the universal gas constant. From (1) the gas density ρ can be written as
(2)ρ=mV=pMZRT.According to the Langmuir isotherm and the state equation of real gas, the adsorbed gas concentrationCm in coal matrix and free gas concentration Cf in fracture systems can be written, respectively, as
(3)Cm(pm)=CLpmpL+pm,Cf=ϕfρ=ϕfpfMZRT,
where pm is the gas pressure in coal matrix, CL is the maximum adsorption concentration of the coal matrix, pL is the Langmuir pressure, pf is the gas pressure in fracture systems, and ϕf is the porosity of fracture systems.


### 2.2. Gas Flow Model

Considering the effects of fracture system deformation, fracture system permeability tensorKf(x,pf) can be written as
(4)Kf(x,pf)=K0eβ(pf-p0),
where p0 is the initial gas pressure in fracture systems, K0 is the fracture system permeability tensor under p0, and β is the permeability modulus with respect to gas pressure.Coalbed gas flow in fracture systems is driven by pressure field and concentration field simultaneously. According to Darcy’s law and Fick’s diffusion law, the motion equation of gas flow in fracture systems can be written as(5)u=-(Kf(x,pf)μg+Dfϕfpf)∇pf,
where u is the total volume velocity of gas flow driven by multifield, μg is the gas viscosity, and Df is the gas diffusion coefficient in fracture systems.Then, according to the mass conservation law, the continuity equation of gas flow in fracture systems is(6)∂ρϕf∂t+∇·(ρu)=ρqm,
where qm is the gas volume of interporosity flow from coal matrix to fracture systems. Using the state equation of real gas, the divergence term of (6) could be unfolded as the following form:
(7)∇·(ρu)=ρ∇·u+u·∇ρ=ρ∇·u+ρpfu·∇pf.
Similarly, the capacity term could be unfolded as the following form:
(8)∂ρϕf∂t=ϕf∂ρ∂t+ρ∂ϕf∂t=ϕfρpf∂pf∂t-cZϕfρ∂pf∂t+cϕϕfρ∂pf∂t,
where cϕ is the compressibility factor of ϕf and cz is the compressibility factor of Z, defined as
(9)cϕ=1ϕf∂ϕf∂pf,cz=1Z∂Z∂pf.
For the expression brevity, define
(10)K(x,pf)=Kf(x,pf)μg+Dfϕfpf,c(x,pf)=ϕfpf-czϕf+cϕϕf,
and substitute motion equation (5) into (7); the continuity equation (6) can be rewritten as follows:
(11)c(x,pf)∂pf∂t-∇·(K(x,pf)∇pf)-K(x,pf)-1pfu2=qm.
Since the seepage velocity u is very low, the quadratic term u2 usually can be omitted in engineering application. Therefore, the more common form of continuity equation is as follows:
(12)c(x,pf)∂pf∂t-∇·(K(x,pf)∇pf)=qm.
In addition, there are some other studies [26, 27] considering the influence of quadratic term u2. For example, Ranjbar et al. [26] have considered the gas density variation in space and solved the nonlinear equations using semianalytical methods. However, they did not consider the effect of gas desorption and deformable media. But their method may form a basis for application of semianalytical methods for problem like the one studied here and it will be useful to the readers.According to Fick’s first diffusion law, the pseudosteady state equation of gas diffusion in coal matrix can be written as(13)dCamdt=Dmσ(CE(pf)-Cam),
where Cam is the average gas concentration of a coal matrix block, CE(pf) is the adsorbed concentration on the coal matrix surface and is defined as CE(pf)=CLpf/(pL+pf), Dm is the gas diffusion coefficient in coal matrix block, and σ is the geometrical factor.Combining (11) and (13), with corresponding initial and boundary conditions, we have the coalbed gas flow model as follows:
(14)(i)c(x,pf)∂pf∂t-∇·(K(x,pf)∇pf)=qm,x∈Ω,t∈[0,T],(ii)dCamdt=Dmσ(CE(pf)-Cam),t∈[0,T],(iii)pf=pd,x∈∂Ω,t∈[0,T],(iv)pf=p0,x∈Ω,t=0,(v)Cam=C0,x∈Vm,t=0,
where Ω is the objective region for the two-dimensional problem with boundary ∂Ω, [0,T] is a time interval with T∈(0,∞], and Vm is the region of a coal matrix block.


## 2.1. Basic Assumptions and Gas State Descriptions

(i)
CBM reservoirs are treated as dual porosity systems consisting of coal matrix and fracture network.
(ii)
CBM reservoirs are isothermal. When pressure variation is not significant, we can assume the gas viscosity is constant under isothermal conditions [25].
(iii)
The gas diffusion process in the coal matrix is pseudostatic and described by Fick’s first diffusion law.
(iv)
Gas absorption/adsorption is described by Langmuir adsorption isotherm.The equation of real gas state is(1)pV=nZRT=mMZRT,
where p is the gas pressure, V is the gas volume, n is the amount of substance, T is the gas temperature, m is the gas mass, M is the gas molar mass, Z is the gas deviation factor, and R is the universal gas constant. From (1) the gas density ρ can be written as
(2)ρ=mV=pMZRT.According to the Langmuir isotherm and the state equation of real gas, the adsorbed gas concentrationCm in coal matrix and free gas concentration Cf in fracture systems can be written, respectively, as
(3)Cm(pm)=CLpmpL+pm,Cf=ϕfρ=ϕfpfMZRT,
where pm is the gas pressure in coal matrix, CL is the maximum adsorption concentration of the coal matrix, pL is the Langmuir pressure, pf is the gas pressure in fracture systems, and ϕf is the porosity of fracture systems.


## 2.2. Gas Flow Model

Considering the effects of fracture system deformation, fracture system permeability tensorKf(x,pf) can be written as
(4)Kf(x,pf)=K0eβ(pf-p0),
where p0 is the initial gas pressure in fracture systems, K0 is the fracture system permeability tensor under p0, and β is the permeability modulus with respect to gas pressure.Coalbed gas flow in fracture systems is driven by pressure field and concentration field simultaneously. According to Darcy’s law and Fick’s diffusion law, the motion equation of gas flow in fracture systems can be written as(5)u=-(Kf(x,pf)μg+Dfϕfpf)∇pf,
where u is the total volume velocity of gas flow driven by multifield, μg is the gas viscosity, and Df is the gas diffusion coefficient in fracture systems.Then, according to the mass conservation law, the continuity equation of gas flow in fracture systems is(6)∂ρϕf∂t+∇·(ρu)=ρqm,
where qm is the gas volume of interporosity flow from coal matrix to fracture systems. Using the state equation of real gas, the divergence term of (6) could be unfolded as the following form:
(7)∇·(ρu)=ρ∇·u+u·∇ρ=ρ∇·u+ρpfu·∇pf.
Similarly, the capacity term could be unfolded as the following form:
(8)∂ρϕf∂t=ϕf∂ρ∂t+ρ∂ϕf∂t=ϕfρpf∂pf∂t-cZϕfρ∂pf∂t+cϕϕfρ∂pf∂t,
where cϕ is the compressibility factor of ϕf and cz is the compressibility factor of Z, defined as
(9)cϕ=1ϕf∂ϕf∂pf,cz=1Z∂Z∂pf.
For the expression brevity, define
(10)K(x,pf)=Kf(x,pf)μg+Dfϕfpf,c(x,pf)=ϕfpf-czϕf+cϕϕf,
and substitute motion equation (5) into (7); the continuity equation (6) can be rewritten as follows:
(11)c(x,pf)∂pf∂t-∇·(K(x,pf)∇pf)-K(x,pf)-1pfu2=qm.
Since the seepage velocity u is very low, the quadratic term u2 usually can be omitted in engineering application. Therefore, the more common form of continuity equation is as follows:
(12)c(x,pf)∂pf∂t-∇·(K(x,pf)∇pf)=qm.
In addition, there are some other studies [26, 27] considering the influence of quadratic term u2. For example, Ranjbar et al. [26] have considered the gas density variation in space and solved the nonlinear equations using semianalytical methods. However, they did not consider the effect of gas desorption and deformable media. But their method may form a basis for application of semianalytical methods for problem like the one studied here and it will be useful to the readers.According to Fick’s first diffusion law, the pseudosteady state equation of gas diffusion in coal matrix can be written as(13)dCamdt=Dmσ(CE(pf)-Cam),
where Cam is the average gas concentration of a coal matrix block, CE(pf) is the adsorbed concentration on the coal matrix surface and is defined as CE(pf)=CLpf/(pL+pf), Dm is the gas diffusion coefficient in coal matrix block, and σ is the geometrical factor.Combining (11) and (13), with corresponding initial and boundary conditions, we have the coalbed gas flow model as follows:
(14)(i)c(x,pf)∂pf∂t-∇·(K(x,pf)∇pf)=qm,x∈Ω,t∈[0,T],(ii)dCamdt=Dmσ(CE(pf)-Cam),t∈[0,T],(iii)pf=pd,x∈∂Ω,t∈[0,T],(iv)pf=p0,x∈Ω,t=0,(v)Cam=C0,x∈Vm,t=0,
where Ω is the objective region for the two-dimensional problem with boundary ∂Ω, [0,T] is a time interval with T∈(0,∞], and Vm is the region of a coal matrix block.


## 3. Notations and Mixed Element Approximation Scheme

In this section, we will present the mixed element approximation scheme for the coalbed gas flow model and give some notations.


### 3.1. Some Notations and Basic Approximation Results

For the convenience of subsequent analysis, we first give notations and basic approximation results. Define the inner product onΩ and its norm:
(15)(u,v)=∫Ω‍u(x)v(x)dx,∥u∥L2(Ω)=(∫Ω‍u2dx)1/2.
Recall that the Sobolev space Hk(Ω) is the closure of C∞(Ω) in the norm
(16)∥u∥k2=∑|α|≤k‍∥Dαu∥L2(Ω)2.
Define the function spaces W, V and their norms as follows:
(17)W=L2(Ω),V=H(div⁡,Ω)={v=(v1,v2)∈(L2(Ω))2∣∇·v∈L2(Ω)},∥ω∥W=(∫Ω‍ω2dx)1/2,∥v∥V={∥v∥2+∥∇·v∥2}1/2.
Let Wh×Vh⊂W×V be the mixed element function spaces such as RTk with index k and discretization parameter h. Now we define the RT projection Πh:V→Vh, satisfying
(18)(∇·(v-Πhv),ωh)=0,v∈V,ωh∈Wh,∥v-Πhv∥0≤C∥v∥shs,1≤s≤k+1,∀v∈V∩(Hs(Ω))2,∥∇·(v-Πhv)∥0≤C∥∇·v∥shs,0≤s≤k+1,∀v∈V∩Hs(div⁡,Ω).
Define standard L2 projection Ph:W→Wh, satisfying
(19)(Phz-z,ωh)=0,z∈W,ωh∈Wh,∥z-Phz∥0≤C∥z∥shs,0≤s≤k+1,∀z∈W∩Hs(Ω).


### 3.2. Mixed Element Approximation Scheme

We consider the case that the seepage-diffusion tensorK(x,pf) is a diagonal matrix and ki,i>0(i=1,2) for the two-dimensional problem. That is to say, there exist positive constants km and kM such that, for all ξ∈R2,
(20)kmξ′ξ≤ξ′K(x,pf)ξ≤kMξ′ξ.
For the practical situation, there exist cm and cM such that the scalar coefficient c(x,pf) satisfies
(21)0<cm≤c(x,pf)≤cM<∞.
Also, we can easily verify that the coefficients ki,i(pf), c(pf), and qm(pf) satisfy the Lipschitz continuity. Since K(x,pf) is a diagonal and invertible matrix, from (5) we have that
(22)K(x,pf)-1u+∇pf=0.Using the integration by parts and applying the boundary condition, we can define the mixed weak formation, which is to findpf∈W, u∈V such that
(23)(K-1(x,pf)u,v)-(∇·v,pf)=-(pd,v·n)∂Ω,∀v∈V,(c(x,pf)∂pf∂t,ω)+(∇·u,ω)=(qm(pf),ω),∀ω∈W,
where n is the unit exterior normal vector to the ∂Ω.Replacing the original pressurepf and velocity u by their approximations, we get the semidiscrete mixed element approximate problem, which is to find ph∈Wh, uh∈Vh such that
(24)(K-1(x,ph)uh,vh)-(∇·vh,ph)=-(pd,vh·n)∂Ω,∀vh∈Vh,(c(x,ph)∂ph∂t,ωh)+(∇·uh,ωh)=(qm(t,ph),ωh),∀ωh∈Wh.LetΔt>0, M=T/Δt, an integer, and tn=nΔt, n=0,1,…,M. We can define the full discrete mixed element approximate scheme with backward Euler time-discretization as follows:
(25)(K-1(phn-1)uhn,vh)-(∇·vh,phn)=-(pd,vh·n)∂Ω,∀vh∈Vh,(c(phn-1)phn-phn-1Δt,ωh)+(∇·uhn,ωh)=(qm(phn-1),ωh),∀ωh∈Wh.
Since gas diffusion equation (13) is an ordinary differential equation, when phn is known, we can get Cam(tn) as follows:
(26)Cam(tn)=CE(phn)-(CE(phn)-C0)e-Dmσtn.
Then, the gas volume of interporosity flow from coal matrix to fracture systems could be calculated as follows:
(27)qm(phn)=-1ρdCamdt|t=tn=-1ρDmσ(CE(phn)-C0)e-Dmσtn.
In practical calculation, alternative and iterative method is used to solve (25)–(27), and specific calculation procedure is listed as follows. (i)
Forn=1,2,…,M, when phn-1 is known, qm(phn-1) can be calculated using (27).(ii)
Then, substitutingqm(phn-1) into (25), (phn,uhn) can be approximated using (25).(iii)
Lastly, whenphn is known, Cam(tn) can be calculated using (26).Letτh be a quasiregular triangulation for the two-dimensional rectangular region Ω with mesh size h. For a triangular element T(A1A2A3), let xi=(xi,yi) be the coordinate of vertex Ai and γi the edge opposite to vertex Ai. Let ni be the outward unit vector on the edge γi and hi the length of the perpendicular dropped from the vertex Ai onto the edge γi, i=1,2,3. We denote by Ϝh the set of all edges of the triangulation τh.The pressure function spaceWh is the space of piecewise constant functions:
(28)Wh={ωh∣ωh=ωk≡constantkinTk,Tk∈τh}
and the dimension of the space Wh is equal to m1—the total number of elements Tk in the triangulation τh. Each function ωh∈Wh can be represented as a linear combination:
(29)ωh=∑k=1m1‍ωkψk(x),
where ψk(x) is the basis function associated with the element Tk, satisfying
(30)ψk(x)={1,x∈Tk0,x∉Tkk=1⋯m1.The velocity function spaceVh is chosen to be the RT0(Ω,τh), defined as
(31)Vh=RT0(Ω,τh)={vh∣vh=(ax+bxay+by)inTk,Tk∈τh}
and the dimension of the space Vh is equal to m2—the total number of edges γj∈Ϝh. Each function vh∈Wh can be represented as a linear combination
(32)vh=∑j=1m2‍vjϕj(x),
where vj=vh·nj and ϕj is the basis function associated with the edge γj∈Ϝ, satisfying
(33)ϕj·nj′={1,j=j′0,j≠j′.
To be specific, for a triangle element Tk, the basis functions ϕi(i=1,2,3) are determined by the following formulas:
(34)ϕi(x)=1hi(x-xi)≡|γi|2|T|(x-xi),
where xi=(xi,yi) are the coordinates of the vertices Ai(i=1,2,3) and |T| is the area of triangle element Tk.Now we use the basis functionsψk(x)(k=1⋯m1) and ϕj(x)(j=1⋯m2) to represent phn and uhn:
(35)phn=∑k=1m1‍pknψk(x),uhn=∑j=1m2‍ujnϕj(x).
Substituting (35) into (25) and setting ωh=ψk(k=1⋯m1) and vh=ϕj(x)(j=1⋯m2), full discrete mixed element approximation scheme (25) can be written as follows:
(36)∑j=1m2‍(K-1(phn-1)ϕj,ϕt)ujn-∑k=1m1‍(∇·ϕt,ψk)pkn=-(pd,ϕt·n)∂Ω,t=1,2⋯m2,Δt∑j=1m2‍(∇·ϕj,ψl)ujn+∑k=1m1‍(c(phn-1)ψk,ψl)pkn=Δt(qm(phn-1),ψl)+(c(phn-1)phn-1,ψl)l=1,2⋯m1.We introduce matrices and vectors as follows:(37)A=(Aij)m2×m2,B=(Bij)m2×m1,C=(Cij)m1×m1,Un=(ujn)m2×1,Pn=(pkn)m1×1,b1=(b1,i)m2×1,b2=(b2,i)m1×1,
where
(38)Aij=∫Ω‍K-1(phn-1)ϕjϕjdx,Bij=∫Ω‍∇·ϕiψjdx,Cij=∫Ω‍c(phn-1)ψiψjdx,b1,i=-∫∂Ω‍pd(ϕi·n)dx,b2,i=Δt∫Ω‍qm(phn-1)ψidx+∫Ω‍c(phn-1)phn-1ψidx.
Therefore, (36) can be written in a matrix form as follows:
(39)(A-BΔtBTC)(UnPn)=(b1b2),
where BT is the transpose of B. Solving linear algebra problem (39), we can get approximation solutions Pn and Un at nth time step, simultaneously.


## 3.1. Some Notations and Basic Approximation Results

For the convenience of subsequent analysis, we first give notations and basic approximation results. Define the inner product onΩ and its norm:
(15)(u,v)=∫Ω‍u(x)v(x)dx,∥u∥L2(Ω)=(∫Ω‍u2dx)1/2.
Recall that the Sobolev space Hk(Ω) is the closure of C∞(Ω) in the norm
(16)∥u∥k2=∑|α|≤k‍∥Dαu∥L2(Ω)2.
Define the function spaces W, V and their norms as follows:
(17)W=L2(Ω),V=H(div⁡,Ω)={v=(v1,v2)∈(L2(Ω))2∣∇·v∈L2(Ω)},∥ω∥W=(∫Ω‍ω2dx)1/2,∥v∥V={∥v∥2+∥∇·v∥2}1/2.
Let Wh×Vh⊂W×V be the mixed element function spaces such as RTk with index k and discretization parameter h. Now we define the RT projection Πh:V→Vh, satisfying
(18)(∇·(v-Πhv),ωh)=0,v∈V,ωh∈Wh,∥v-Πhv∥0≤C∥v∥shs,1≤s≤k+1,∀v∈V∩(Hs(Ω))2,∥∇·(v-Πhv)∥0≤C∥∇·v∥shs,0≤s≤k+1,∀v∈V∩Hs(div⁡,Ω).
Define standard L2 projection Ph:W→Wh, satisfying
(19)(Phz-z,ωh)=0,z∈W,ωh∈Wh,∥z-Phz∥0≤C∥z∥shs,0≤s≤k+1,∀z∈W∩Hs(Ω).


## 3.2. Mixed Element Approximation Scheme

We consider the case that the seepage-diffusion tensorK(x,pf) is a diagonal matrix and ki,i>0(i=1,2) for the two-dimensional problem. That is to say, there exist positive constants km and kM such that, for all ξ∈R2,
(20)kmξ′ξ≤ξ′K(x,pf)ξ≤kMξ′ξ.
For the practical situation, there exist cm and cM such that the scalar coefficient c(x,pf) satisfies
(21)0<cm≤c(x,pf)≤cM<∞.
Also, we can easily verify that the coefficients ki,i(pf), c(pf), and qm(pf) satisfy the Lipschitz continuity. Since K(x,pf) is a diagonal and invertible matrix, from (5) we have that
(22)K(x,pf)-1u+∇pf=0.Using the integration by parts and applying the boundary condition, we can define the mixed weak formation, which is to findpf∈W, u∈V such that
(23)(K-1(x,pf)u,v)-(∇·v,pf)=-(pd,v·n)∂Ω,∀v∈V,(c(x,pf)∂pf∂t,ω)+(∇·u,ω)=(qm(pf),ω),∀ω∈W,
where n is the unit exterior normal vector to the ∂Ω.Replacing the original pressurepf and velocity u by their approximations, we get the semidiscrete mixed element approximate problem, which is to find ph∈Wh, uh∈Vh such that
(24)(K-1(x,ph)uh,vh)-(∇·vh,ph)=-(pd,vh·n)∂Ω,∀vh∈Vh,(c(x,ph)∂ph∂t,ωh)+(∇·uh,ωh)=(qm(t,ph),ωh),∀ωh∈Wh.LetΔt>0, M=T/Δt, an integer, and tn=nΔt, n=0,1,…,M. We can define the full discrete mixed element approximate scheme with backward Euler time-discretization as follows:
(25)(K-1(phn-1)uhn,vh)-(∇·vh,phn)=-(pd,vh·n)∂Ω,∀vh∈Vh,(c(phn-1)phn-phn-1Δt,ωh)+(∇·uhn,ωh)=(qm(phn-1),ωh),∀ωh∈Wh.
Since gas diffusion equation (13) is an ordinary differential equation, when phn is known, we can get Cam(tn) as follows:
(26)Cam(tn)=CE(phn)-(CE(phn)-C0)e-Dmσtn.
Then, the gas volume of interporosity flow from coal matrix to fracture systems could be calculated as follows:
(27)qm(phn)=-1ρdCamdt|t=tn=-1ρDmσ(CE(phn)-C0)e-Dmσtn.
In practical calculation, alternative and iterative method is used to solve (25)–(27), and specific calculation procedure is listed as follows. (i)
Forn=1,2,…,M, when phn-1 is known, qm(phn-1) can be calculated using (27).(ii)
Then, substitutingqm(phn-1) into (25), (phn,uhn) can be approximated using (25).(iii)
Lastly, whenphn is known, Cam(tn) can be calculated using (26).Letτh be a quasiregular triangulation for the two-dimensional rectangular region Ω with mesh size h. For a triangular element T(A1A2A3), let xi=(xi,yi) be the coordinate of vertex Ai and γi the edge opposite to vertex Ai. Let ni be the outward unit vector on the edge γi and hi the length of the perpendicular dropped from the vertex Ai onto the edge γi, i=1,2,3. We denote by Ϝh the set of all edges of the triangulation τh.The pressure function spaceWh is the space of piecewise constant functions:
(28)Wh={ωh∣ωh=ωk≡constantkinTk,Tk∈τh}
and the dimension of the space Wh is equal to m1—the total number of elements Tk in the triangulation τh. Each function ωh∈Wh can be represented as a linear combination:
(29)ωh=∑k=1m1‍ωkψk(x),
where ψk(x) is the basis function associated with the element Tk, satisfying
(30)ψk(x)={1,x∈Tk0,x∉Tkk=1⋯m1.The velocity function spaceVh is chosen to be the RT0(Ω,τh), defined as
(31)Vh=RT0(Ω,τh)={vh∣vh=(ax+bxay+by)inTk,Tk∈τh}
and the dimension of the space Vh is equal to m2—the total number of edges γj∈Ϝh. Each function vh∈Wh can be represented as a linear combination
(32)vh=∑j=1m2‍vjϕj(x),
where vj=vh·nj and ϕj is the basis function associated with the edge γj∈Ϝ, satisfying
(33)ϕj·nj′={1,j=j′0,j≠j′.
To be specific, for a triangle element Tk, the basis functions ϕi(i=1,2,3) are determined by the following formulas:
(34)ϕi(x)=1hi(x-xi)≡|γi|2|T|(x-xi),
where xi=(xi,yi) are the coordinates of the vertices Ai(i=1,2,3) and |T| is the area of triangle element Tk.Now we use the basis functionsψk(x)(k=1⋯m1) and ϕj(x)(j=1⋯m2) to represent phn and uhn:
(35)phn=∑k=1m1‍pknψk(x),uhn=∑j=1m2‍ujnϕj(x).
Substituting (35) into (25) and setting ωh=ψk(k=1⋯m1) and vh=ϕj(x)(j=1⋯m2), full discrete mixed element approximation scheme (25) can be written as follows:
(36)∑j=1m2‍(K-1(phn-1)ϕj,ϕt)ujn-∑k=1m1‍(∇·ϕt,ψk)pkn=-(pd,ϕt·n)∂Ω,t=1,2⋯m2,Δt∑j=1m2‍(∇·ϕj,ψl)ujn+∑k=1m1‍(c(phn-1)ψk,ψl)pkn=Δt(qm(phn-1),ψl)+(c(phn-1)phn-1,ψl)l=1,2⋯m1.We introduce matrices and vectors as follows:(37)A=(Aij)m2×m2,B=(Bij)m2×m1,C=(Cij)m1×m1,Un=(ujn)m2×1,Pn=(pkn)m1×1,b1=(b1,i)m2×1,b2=(b2,i)m1×1,
where
(38)Aij=∫Ω‍K-1(phn-1)ϕjϕjdx,Bij=∫Ω‍∇·ϕiψjdx,Cij=∫Ω‍c(phn-1)ψiψjdx,b1,i=-∫∂Ω‍pd(ϕi·n)dx,b2,i=Δt∫Ω‍qm(phn-1)ψidx+∫Ω‍c(phn-1)phn-1ψidx.
Therefore, (36) can be written in a matrix form as follows:
(39)(A-BΔtBTC)(UnPn)=(b1b2),
where BT is the transpose of B. Solving linear algebra problem (39), we can get approximation solutions Pn and Un at nth time step, simultaneously.


## 4. Numerical Examples

In this section, we carry out numerical experiments using the lowest order Raviart-Thomas mixed finite elementsRT0 for the gas flow problem in deformable CBM reservoirs.Example 1.
In order to verify the convergence, the test region is selected as unit square; that is,Ω=[0,1]×[0,1], and 6 levels for h=1/4, 1/8, 1/16, 1/32, 1/64, 1/128 are computed to estimate the convergence rate. Since we cannot get the analytical solution pf, u, and Cam for the problem (14), we use ∥phN-p1/128N∥L2/∥p1/128N∥L2, ∥uhN-u1/128N∥L2/∥u1/128N∥L2, and ∥Cam,hN-Cam,1/128N∥L2/∥Cam,1/128N∥L2 as the criterion of convergence for the pressure, velocity, and concentration. Error estimates in L2-norm and convergence rate estimates for ph, uh, and Cam,h are listed in Tables 1, 2, and 3, respectively. We can find that the approximation solutions ph, uh, and Cam,h could achieve first-order convergence rate (see Figure 1).Table 1
Relative errors and convergence rates forph.
Resolution
Error(∥phN-p1/128N∥L2∥p1/128N∥L2)
Rate
h
=
1
4
6.67
E
-
2
—
h
=
1
8
3.60
E
-
2
-
0.8903
h
=
1
16
1.91
E
-
2
-
0.9131
h
=
1
32
9.63
E
-
3
-
0.9880
h
=
1
64
4.81
E
-
3
-
1.0015Table 2
Relative errors and convergence rates foruh.
Resolution
Error(∥uhN-u1/128N∥L2∥u1/128N∥L2)
Rate
h
=
1
4
1.17
E
-
2
—
h
=
1
8
6.23
E
-
3
-
0.9092
h
=
1
16
3.24
E
-
3
-
0.9432
h
=
1
32
1.64
E
-
3
-
0.9823
h
=
1
64
8.17
E
-
4
-
1.0053Table 3
Relative errors and convergence rates forCam,h.
Resolution
Error(∥Cam,hN-Cam,1/128N∥L2∥Cam,1/128N∥L2)
Rate
h
=
1
4
4.77
E
-
2
—
h
=
1
8
2.54
E
-
2
-
0.9092
h
=
1
16
1.31
E
-
2
-
0.9553
h
=
1
32
6.59
E
-
3
-
0.9912
h
=
1
64
3.29
E
-
3
-
1.0022Figure 1
Convergence rates for∥phN-p1/128N∥L2/∥p1/128N∥L2, ∥uhN-u1/128N∥L2/∥u1/128N∥L2, and ∥Cam,hN-Cam,1/128N∥L2/∥Cam,1/128N∥L2.Example 2.
For a simple application, the scale of objective regionΩ is 100m×100m and the mesh size h=2.5m. There is a production well at the center of region Ω with the production rate qg. Simulation parameters are listed in Table 4.
Substituting the parameter value into calculation procedure, we can get the numerical solutions, so as to describe the dynamic characteristic of gas flow in deformable CBM reservoirs.
Since average pressure could reflect the decrement of reservoir driving energy during production process, Figures2 and 3 illustrate the gas average pressure in fracture systems varying with time t and show the effects of Dm and β on average pressure, respectively. In Figure 2, the model has the initial reservoir pressure of 10 MPa, and diffusion coefficient of coal matrix Dm is taken as 1×10-7, 5×10-8, 1×10-8, and 5×10-9 m3/s, respectively, so as to analyze the effects of Dm on average pressure in fracture systems. In Figure 2, we have found that pressure variation curves could be divided into 3 stages. Firstly, at the initial stage, since coalbed gas could not desorb and diffuse from matrix into fracture systems instantly, average pressure in fracture systems declines rapidly. At this stage, we have found when Dm is increasing, the initial stage became shorter. It is because that larger Dm could make fracture systems get pressure compensation from coal matrix more quickly. Secondly, at the stable stage, with gas desorption increasing from the coal matrix into fracture systems, fracture systems could get more pressure compensation; thus, the average pressure drop speed slows down. We have also found, with a bigger Dm, fracture systems could maintain a higher average pressure, while the stable stage would be shorter. Lastly, at the last stage, due to the decline of gas concentration in matrix gas desorption rate decreases, resulting in the pressure drop speed becoming rapid once again. Meanwhile, when Dm is larger, the last stage appears at an earlier time. The double porosity characteristics of CBM reservoirs have been just shown through the pressure variation process presented above.
In Figure3, the permeability modulus β is set as 1×10-1, 1×10-3, and 1×10-5, respectively, so as to analyze the effects of β on average pressure. We have found that pressure drop curves with different β are almost coincident at early time, which indicates that fracture deformation has no evident effects on permeability and pressure at this stage. Essentially, with the fracture pressure dropping, the effects of fracture deformation on permeability become significant, and, conversely, the permeability could also affect fracture pressure. That is to say, the effects of permeability and pore pressure are mutual. Subsequently, we have found pressure drop curves become apart in Figure 3. Also, we have found that the larger the β is, the lower the average pressure will be. From the definition of β, we know that β reflects the effects of fracture system deformation on the permeability Kf. Thus, when β is larger, Kf decreases more significantly with pf dropping, such that CBM reservoirs can only get less pressure compensation from outer boundary, resulting in a lower average pressure.
Figure4 shows the bottom hole pressure varying with time t, and the diffusion coefficient Dm of coal matrix is taken as 1×10-7, 5×10-8, 1×10-8, and 5×10-9 m3/s, respectively, so as to analyze the effects of Dm on bottom hole pressure. Since the bottom hole pressure is the direct reflection of production effects, the bottom hole pressure drops sharply from 10 MPa to about 4.25 MPa after starting production, and then there is a pressure recovery due to the gas desorption from coal matrix to fracture systems. Thus, this causes an early fluctuation before reaching a stable state. In the subsequent stable stage, we have found that a bigger diffusion coefficient Dm could make more coalbed gas desorb into fracture systems and result in a higher bottom hole pressure. Then, at the last stage, due to the gas desorption decreasing, the effects of Dm become weak and bottom hole pressure curves tend to be close. Put simply, the effects of Dm are evident at stable stage, while they are not so evident at the early and last stage.
Coalbed gas flow in fracture systems is driven by pressure field and concentration field simultaneously; thus, the effects of gas diffusion in fracture systems need to be investigated. In Figure5, we set fracture diffusion coefficient Df as 0,1×10-2, 1×10-3, and 1×10-5 m3/s, respectively, so as to analyze the effects of gas diffusion in fracture systems. We have found that when Df is larger, the bottom hole pressure is higher. Specifically, the pressure difference between Df=1×10-2 and Df=0 is about 0.35 MPa all through the stable stage in Figure 5. That is to say, existence of gas diffusion in fracture systems could make fracture systems get more pressure compensation from outer boundary, resulting in a higher bottom hole pressure. Thus, it is very meaningful to consider gas diffusion in fracture system.
In order to illustrate the gas desorption dynamic process from coal matrix into fracture systems, we show the gas average desorption rate varying with timet in Figures 6 and 7. In Figure 6, the diffusion coefficient Dm is taken as 1×10-7, 5×10-8, 1×10-8, and 5×10-9 m3/s, respectively, so as to investigate the effects of matrix diffusion coefficient on gas desorption process. According to equilibrium desorption model, we know that gas desorption rate is affected by gas average concentration Cam in matrix, fracture pressure pf, and diffusion coefficient Dm. During production process, as fracture pressure pf and gas average concentration Cam decline, there are two stages in gas desorption process. At first half stage, due to quick drawdown of fracture pressure, gas concentration CE(pf) on matrix surface declines sharply, and then the gas average desorption rate increased to the maximum. Then, at the second half stage, the desorption rate declines due to the reduction of gas average concentration Cam. In Figure 6, we have also found when Dm is larger, average desorption rate will be higher and peak desorption rate appears at an earlier moment. In Figure 7, the fracture diffusion coefficient Df is taken as 0,1×10-2, 1×10-3, and 1×10-5 m3/s, respectively. Since the effects of Df on desorption rate are indirect and reflected by means of fracture pressure pf, the effects of Df are not evident at the starting and ending time. However, we have still found that the larger the Df is, the higher the peak desorption rate will be.
In order to reflect the remaining reserves of coalbed gas in coal matrix, we show the gas average concentration in coal matrix varying with timet in Figures 8 and 9. In Figure 8, Dm is taken as 1×10-7, 5×10-8, 1×10-8, and 5×10-9 m3/s, respectively. We have found that the gas average concentration declines slowly at the early stage on account of the existence of desorption delay phenomenon. Then, due to gas desorption increasing from coal matrix into fracture systems, the average concentration in coal matrix declines quickly at the later stage. At last stage, due to gas desorption rate decreasing, the concentration drop speed becomes slow once again. Also, in Figure 8, we have found concentration curves are almost coincident at early time. At subsequent stages, when Dm is larger, the average concentration is lower and the concentration drop speed is quicker. It is because that larger Dm could accelerate gas desorption rate. In Figure 9, we have found that the effects of initial permeability k0 are not evident at the early stage, due to the existence of desorption delay phenomenon. And then, at the later half stage, the average concentration increases with k0 increasing. It indicates that the CBM reservoir with higher permeability could get more pressure compensation from outer boundary, resulting in higher pore pressure and matrix concentration.Table 4
Parameters for simulation.
Coefficient
Value
Ω
100 m × 100 m
ϕ
f
0.15
p
0
10 MPa
C
0
22.5 kg/m3
D
m
1
0
-
8 m3/s
P
L
6 MPa
c
Z
1
0
-
3
q
g
180 m3/h
h
2.5 m
k
0
0.01 md
μ
g
0.65 mPa·s
D
f
1
0
-
6 m3/s
C
L
36 kg/m3
c
ϕ
5
×
1
0
-
2
β
1.25
×
1
0
-
3Figure 2
Gas average pressure of fracture systems with differentDm.Figure 3
Gas average pressure of fracture systems with differentβ.Figure 4
Bottom hole pressure of CBM reservoirs with differentDm.Figure 5
Bottom hole pressure of CBM reservoirs with differentDf.Figure 6
Gas average desorption rate with differentDm.Figure 7
Gas average desorption rate with differentDf.Figure 8
Gas average concentration in coal matrix with differentDm.Figure 9
Gas average concentration in coal matrix with differentk0.


## 5. Conclusions

(1)
In this paper, we present a desorption-diffusion-seepage model of gas flow in deformable CBM reservoirs, in which the effects of fracture system deformation have been considered.
(2)
We introduce a mixed finite element method to approximate the coalbed gas pressure and velocity, simultaneously, and establish error estimates for the numerical solutions.
(3)
We carry out numerical experiments using the lowest order Raviart-Thomas mixed finite elementsRT0. Numerical results in Example 1 indicate that the approximation solutions ph, uh, and Cam,h could achieve first-order convergence rate. Numerical results in Example 2 have shown that the existence of fracture system deformation phenomenon will result in a lower gas average pressure.



---

*Source: 735931-2014-09-28.xml*