# Numerical Simulation of Unsteady-State Flow in Dual Porous Coalbed Methane Horizontal Wells with Complex Boundary Conditions

Authors: Cheng-yong Li; Jun Zhou; Xiang-yi Yi; Yi Luo; Ping-zhi Gong
Journal: Journal of Chemistry
 (2015)

---

## Abstract

The bottom-hole pressure response which can reflect the gas flow characteristics is important to study. A mathematical model for description of gas from porous coalbed methane (CBM) reservoirs with complex boundary conditions flowing into horizontal wells has been developed. Meanwhile, basic solution of boundary elements has been acquired by combination of Lord Kelvin point source solution, the integral of Bessel function, and Poisson superimpose formula for CBM horizontal wells with complex boundary conditions. Using this model, type curves of dimensionless pressure and pressure derivative are obtained, and flow characteristics of horizontal wells in complex boundary reservoirs and relevant factors are accordingly analyzed.


---

## Body

## 1. Introduction

Coalbed methane (CBM) is a kind of green and clean energy. The development and utilization of coalbed methane could not only relieve the tense situation of conventional oil and gas in supply, but also reduce the atmospheric environment pollution.Different from conventional gas reservoirs, the migration mechanism of gas in coal is more complex and diverse [1]. Coal is dual porous media reservoir, matrix is main storage space of CBM adsorption, and fractures are main transport routes of CBM diffusion-seepage. Analyzing bottom-hole pressure helps to figure out CBM production status. Ertekin and Sung [2] established a nonequilibrium adsorption nonsteady seepage flow model in dual porosity media with Henry’s law. Applying Fick’s law, Anbarci and Ertekin [3] suggested a single-phase CBM seepage flow mathematical model considering pseudosteady and unsteady diffusion phenomenon. Clarkson and Bustin [4, 5] put forward a new double diffusion model, which assumes that adsorption occurs only in micropores and conforms to the law of nonlinear adsorption. The macropores accumulate the free gas or provide channels of gas migration between micropores and fractures. Reeve [6] proposed a new gas-water two-phase triple-porosity dual-seepage flow mathematical model, and this model can increase the third matrix pore system; however, the double-permeability model is very complex and difficult to describe and calculate. Tong et al. [7–9] introduced permeability modulus considering the deformation of coal and developed a pseudosteady diffusion nonequilibrium adsorption nonsteady seepage mathematical model. Hu et al. [10] established a gas-water two-phase percolation mathematical model of well test interpretation for CBM reservoirs, and its correctness has been verified by simulation of CBM seepage flow characteristics. Clarkson et al. [11] introduced pseudopressure function to analyze production of gas and water flow performance using the numerical simulation method. Aminian and Ameri [12] set up a function to predict gas production based on storage and transport mechanisms in CBM reservoirs. Recently, Cai and Yu [13] introduced the fractal theory to study the enhancing recovery mechanism in natural gas-saturated porous media by spontaneous imbibition effect.In order to improve the production, a large number of horizontal wells have been used in the CBM reservoirs. Sung and Ertekin [14–16] established a two-dimensional, two-phase, multiwell gas-water flow model, and this model has the ability to simulate multiple horizontal wells. Engler and Rajtar [17–19] established a mathematical model of single-phase gas flow in the horizontal wells, and the analytical solutions are given for the horizontal well pressure drop and pressure recovery. Wang et al. [20] established a mathematical model, in which anisotropy formation heterogeneity, permeability stress sensitivity, and influence of wellbore pressure drop on directional pinnate horizontal wells in CBM reservoirs are considered. Nie et al. [21] deduced CBM flow equations based on Langmuir adsorption in matrix and Darcy flow in fracture and analyzed the transient transport characteristics of gas from CBM reservoirs to horizontal wells.Generally, the theories of calculating reservoir and bottom-hole pressures which can reflect the gas flow characteristics are mostly based on homogeneous reservoirs and regular geometry, such as infinite boundary or circular boundary. Outer boundary conditions of a reservoir have also been simplified; it is generally regarded as simple situations as constant pressure or closed boundary. However, influenced by characteristics of geological structures, the true reservoirs usually have complex and diversiform boundaries. In this case, the conventional flow theories and solving methods could do nothing to calculate reservoir or bottom-hole pressures with mixed boundary conditions.The boundary element method (BEM) is a numerical computational method of solving linear partial differential equations developed after the better-known finite element method (FEM) and finite difference method (FDM). The BEM could be able to reduce dimension and save computer memory and running time. Coupling with the BEM [22], bottom-hole pressures and complicated flow characteristics under the condition of irregularly shaped area with mixed boundary could be calculated and analyzed. Numbere and Tiab [23] developed a streamline simulation with the BEM for homogeneous or partly homogeneous reservoirs with irregular boundary to make the simulation match physical model better. Kikani and Horne [24] employed the BEM to analyze the transient pressure response in reservoirs with arbitrary boundary and developed two formulas, namely, convolution formula and Laplace domain space formula, to solve transient fluid flow through porous medium in homogeneous reservoirs. Hou et al. [25] used the BEM to simulate flow line map in homogeneous reservoirs with irregular boundary. Chaiyo et al. [26] used the BEM to solve free boundary saturated seepage problem. Rafiezadeh and Ataie-Ashtiani [27] studied the flow mechanism in anisotropic media by three-dimensional boundary elements.In this paper, a mathematical model is developed to describe gas flow in horizontal wells in CBM reservoirs based on the theory of fluid flow through porous. The type curves of pressure derivative, characteristics of gas flow with complex external boundary, and relevant affecting factors are analyzed.


## 2. Physical Model of Gas Flow

Coal reservoir is a dual porous medium composed of matrix and fracture. The matrix is the main reservoir space of coalbed methane (Figure1(a)), and the fracture is the flow channel of fluid (Figure 1(b)). The average pore size of CBM reservoirs is much smaller than those of conventional reservoirs. Pore size can be divided into three categories [28]: macropore (aperture > 20 nm), mesopore (2 nm < aperture < 20 nm), and micropore (aperture < 2 nm). The small porosity of coal leads to the great specific surface area. Hence, methane could be strongly adsorbed, resulting in the fact that the content of CBM is far more than its pore volume.Figure 1
Pore characteristics of coal ((a) micropore; (b) fracture).
(a)
(b)According to the reservoir dual pore structure of coal, a physical model is set up for gas flow. The migration of gas in coal is shown in Figure2. In the production of CBM, formation pressure keeps declining. When formation pressure drops below the critical desorption pressure, CBM starts to desorb from the coal matrix surface. Meanwhile, the original state of equilibrium is broken, causing the flow of gas in fracture. This process is very similar to spontaneous imbibition in fractured porous media [29].Figure 2
Desorption and diffusion in the CBM reservoirs.To facilitate the derivation, it is assumed that the length of a horizontal well isL, the center location of the horizontal well is xw,yw,zw, and the thickness of a CBM reservoir is h. In this paper, gas flow into horizontal wells with closed boundary (Figure 3(a)), constant pressure boundary (Figure 3(b)), and mixed boundary (Figure 3(c)) is considered.Figure 3
Horizontal wells with different boundaries ((a) closed boundary; (b) constant pressure boundary; (c) mixed boundary).
(a)
(b)
(c)The fundamental assumptions are as follows:(1)
CBM diffuses directly from matrix to fracture, and the process of diffusion is unsteady.(2)
Gas flow in fracture is radial laminar flow, in agreement with Darcy’s law.(3)
Only single-phase gas flow exists in coal.(4)
The effects of gravity and capillary force are negligible.(5)
The effects of temperature are negligible.(6)
CBM isothermal adsorption process is in line with the Langmuir isotherm adsorption law, and the initial state conforms to the isothermal adsorption curve.(7)
Radius of the gas well is regarded infinitesimal, and gas well production is constant.


## 3. Mathematical Model

### 3.1. Gas Diffusion Model in Matrix

In combination with the mass conservation equation with the second Fick’s diffusion law, the change of gas concentration with time is given by(1)ϕ∂c∂t=∇D′∇c.The planar radial flow equation is as follows:(2)∂c∂t=DriI∂∂ririI∂c∂ri.The dimensionless equation of gas diffusion in matrix is(3)1riD2∂∂riDriD2∂cD∂riD=λ1∂cD∂tD,where riD is the dimensionless radius defined by riD=ri/R; cD is the dimensionless diffusion concentration defined by cD=c-ci; tD is the dimensionless time defined by tD=3.6kt/θrw2; λ1 is interporosity flow coefficient defined by λ1=3.6kR2/θDrw2; θ is comprehensive storage coefficient defined by θ=ϕfctμ+6pscT/qDTscψi; qD is the dimensionless production defined by qD=1.842×10-3qscpscT/khTscψi; ci is the initial concentration of gas, kg/m3; k is the permeability, μm2; and ψi is the pseudopressure.


### 3.2. Gas Flow Model in Matrix and Fracture

In the process of diffusion of CBM, the concentration often changes. Therefore, unsteady diffusion equation which accords with actual situation is used.Free gas concentration is as follows:(4)c1=ρ1ϕ=MpϕRTZ.Concentration of the adsorbed gas is(5)c2=VLPPL+P.Gas concentration in coal is(6)c=MpϕRTZ+VLPPL+P.Volume of gas desorption from coal is(7)qd=-∂Md∂t=-ρsc∂Vd∂t,where Md is the gas mass under standard condition, kg/m3; Vd is the volume of gas adsorption in coal under standard condition, m3/m3. ρsc is the density of gas under standard condition, kg/m3.Density of gas is defined as(8)ρsc=MpsczRTsc.Spherical matrix has the following relationship:(9)∂Vd∂t=-AVJ=-AVD∂c∂ri=3DR∂c∂ri,where V is the volume of coal matrix, m3; A is the surface area of coal matrix, m2; D is the mass diffusion coefficient, m2/s; J is the diffusion flux, g/m2s.Combining (7), (8), and (9), the volume of gas desorption is given by(10)qd=MpscRTsc3DR∂c∂ri.Based on the basic principle of material balance, the governing equations of gas flow in fracture system are given by(11)1r∂∂rrρkμg∂p∂r+qd=ρϕCt∂p∂t.Equation (11) is simplified into the following equation: (12)1rD∂∂rDrD∂ψD∂rD+1-ωλ∂cD∂riD=ω∂ψD∂tD,where ω is the fracture storage ratio given by ω=ϕctμ/θ; λ is the diffusion coefficient given by λ=3.6kτ/θrw2; ψD is the dimensionless pseudopressure given by ψD=(ψi-ψ)/ψiqD; τ is the adsorption time given by τ=R2/D.


### 3.3. Mathematical Model Solution

Diffusion and governing equations of gas flow in fracture are combined as follows:(13)1riD2∂∂riDriD2∂cD∂riD=λ∂cD∂tD,1riD2∂∂ridiffusioninmatrix,1rD∂∂rDrD∂ψD∂rD+1-ωλ∂cD∂riD=ω∂ψD∂tD,1rD∂∂rDrD∂ψD∂rD+fgigflowinfracture.Dimensionless initial and boundary conditions are given as follows.The initial condition is(14)ψDrD,tD=0=0.The boundary condition is(15)∂ψD∂rDrD=1,tD=-1.The infinite outer boundary condition is(16)ψDrD⟶∞,tD=0.The constant pressure boundary condition is(17)ψDreD,tD=0.The closed boundary condition is(18)∂ψD∂rDreD,tD=0.WhenM=riDcD, (13) can be transformed into(19)∂2M∂riD2=λ∂M∂tD.The Laplace transform of (19) is(20)∂2M-∂riD2=λuM-.The general solution of (20) is (21)M-=Asinh⁡λuriD+Bcosh⁡λuriD,where sinh⁡x=(ex-e-x)/2 and cosh⁡x=(ex+e-x)/2.The coefficientsA and B could be obtained by the initial and boundary conditions:(22)B=0,A=M-asinh⁡λu.Hence, (23)M-=M-asinh⁡λuriDsinh⁡λu.Therefore, the diffusion equation is transformed into (24)∂c-D∂riDriD=1=c-Dλucoth⁡λu-1.Combining the dimensionless definitioncD and Langmuir isothermal adsorption equation V=VLPm/(PL+Pm), (24) is transformed into (25)∂c-D∂rDrD=1=LVLPPL+P-VLPiPL+Piλucoth⁡λu-1,where L is Laplace transform: (26)LVLPPL+P-VLPiPL+Pi=βP-D.So,(27)∂c-D∂rDrD=1=-βP-Dλucoth⁡λu-1.Laplace transform of (20) is (28)1rD∂∂rDrD∂ψ-D∂rD+1-ωλ∂c-D∂riD=ωuψ-D.The initial conditions are(29)ψ-Du→∞=0.The inner boundary conditions are(30)rD∂ψ-D∂rDrD=1=-1u.The outer boundary conditions are(31)ψ-DrD→∞=0infinite,ψ-DrD=reD=0constantpressure,∂ψ-D∂rDrD=reD=0closed.The solution of diffusion equation (27) is(32)∂c-D∂riDriD=1=c-Dλucoth⁡λu-1.Substituting the dimensionless definitioncD and Langmuir isothermal adsorption equation V=VLPm/(PL+Pm) in (32), it can be written as follows:(33)∂c-D∂riDrD=1=LVLPPL+P-VLPiPL+Piλucoth⁡λu-1.SubstitutingψD=(ψi-ψ)/ψiqD in (33), it can be expressed as follows:(34)VLPPL+P-VLPiPL+Pi=-ψLVLψiqDψL+ψψL+ψiψ+ψiψ-D.Definingβ=ψLVLψiqD/ψL+ψψL+ψiψ+ψi,(35)∂c-D∂riDriD=1=-βψ-Dλucoth⁡λu-1.Substituting (35) in (28),(36)1rD∂∂rDrD∂ψ-D∂rD=fuψ-D,fu=ωu+1-ωλβλucoth⁡λu-1.


## 3.1. Gas Diffusion Model in Matrix

In combination with the mass conservation equation with the second Fick’s diffusion law, the change of gas concentration with time is given by(1)ϕ∂c∂t=∇D′∇c.The planar radial flow equation is as follows:(2)∂c∂t=DriI∂∂ririI∂c∂ri.The dimensionless equation of gas diffusion in matrix is(3)1riD2∂∂riDriD2∂cD∂riD=λ1∂cD∂tD,where riD is the dimensionless radius defined by riD=ri/R; cD is the dimensionless diffusion concentration defined by cD=c-ci; tD is the dimensionless time defined by tD=3.6kt/θrw2; λ1 is interporosity flow coefficient defined by λ1=3.6kR2/θDrw2; θ is comprehensive storage coefficient defined by θ=ϕfctμ+6pscT/qDTscψi; qD is the dimensionless production defined by qD=1.842×10-3qscpscT/khTscψi; ci is the initial concentration of gas, kg/m3; k is the permeability, μm2; and ψi is the pseudopressure.


## 3.2. Gas Flow Model in Matrix and Fracture

In the process of diffusion of CBM, the concentration often changes. Therefore, unsteady diffusion equation which accords with actual situation is used.Free gas concentration is as follows:(4)c1=ρ1ϕ=MpϕRTZ.Concentration of the adsorbed gas is(5)c2=VLPPL+P.Gas concentration in coal is(6)c=MpϕRTZ+VLPPL+P.Volume of gas desorption from coal is(7)qd=-∂Md∂t=-ρsc∂Vd∂t,where Md is the gas mass under standard condition, kg/m3; Vd is the volume of gas adsorption in coal under standard condition, m3/m3. ρsc is the density of gas under standard condition, kg/m3.Density of gas is defined as(8)ρsc=MpsczRTsc.Spherical matrix has the following relationship:(9)∂Vd∂t=-AVJ=-AVD∂c∂ri=3DR∂c∂ri,where V is the volume of coal matrix, m3; A is the surface area of coal matrix, m2; D is the mass diffusion coefficient, m2/s; J is the diffusion flux, g/m2s.Combining (7), (8), and (9), the volume of gas desorption is given by(10)qd=MpscRTsc3DR∂c∂ri.Based on the basic principle of material balance, the governing equations of gas flow in fracture system are given by(11)1r∂∂rrρkμg∂p∂r+qd=ρϕCt∂p∂t.Equation (11) is simplified into the following equation: (12)1rD∂∂rDrD∂ψD∂rD+1-ωλ∂cD∂riD=ω∂ψD∂tD,where ω is the fracture storage ratio given by ω=ϕctμ/θ; λ is the diffusion coefficient given by λ=3.6kτ/θrw2; ψD is the dimensionless pseudopressure given by ψD=(ψi-ψ)/ψiqD; τ is the adsorption time given by τ=R2/D.


## 3.3. Mathematical Model Solution

Diffusion and governing equations of gas flow in fracture are combined as follows:(13)1riD2∂∂riDriD2∂cD∂riD=λ∂cD∂tD,1riD2∂∂ridiffusioninmatrix,1rD∂∂rDrD∂ψD∂rD+1-ωλ∂cD∂riD=ω∂ψD∂tD,1rD∂∂rDrD∂ψD∂rD+fgigflowinfracture.Dimensionless initial and boundary conditions are given as follows.The initial condition is(14)ψDrD,tD=0=0.The boundary condition is(15)∂ψD∂rDrD=1,tD=-1.The infinite outer boundary condition is(16)ψDrD⟶∞,tD=0.The constant pressure boundary condition is(17)ψDreD,tD=0.The closed boundary condition is(18)∂ψD∂rDreD,tD=0.WhenM=riDcD, (13) can be transformed into(19)∂2M∂riD2=λ∂M∂tD.The Laplace transform of (19) is(20)∂2M-∂riD2=λuM-.The general solution of (20) is (21)M-=Asinh⁡λuriD+Bcosh⁡λuriD,where sinh⁡x=(ex-e-x)/2 and cosh⁡x=(ex+e-x)/2.The coefficientsA and B could be obtained by the initial and boundary conditions:(22)B=0,A=M-asinh⁡λu.Hence, (23)M-=M-asinh⁡λuriDsinh⁡λu.Therefore, the diffusion equation is transformed into (24)∂c-D∂riDriD=1=c-Dλucoth⁡λu-1.Combining the dimensionless definitioncD and Langmuir isothermal adsorption equation V=VLPm/(PL+Pm), (24) is transformed into (25)∂c-D∂rDrD=1=LVLPPL+P-VLPiPL+Piλucoth⁡λu-1,where L is Laplace transform: (26)LVLPPL+P-VLPiPL+Pi=βP-D.So,(27)∂c-D∂rDrD=1=-βP-Dλucoth⁡λu-1.Laplace transform of (20) is (28)1rD∂∂rDrD∂ψ-D∂rD+1-ωλ∂c-D∂riD=ωuψ-D.The initial conditions are(29)ψ-Du→∞=0.The inner boundary conditions are(30)rD∂ψ-D∂rDrD=1=-1u.The outer boundary conditions are(31)ψ-DrD→∞=0infinite,ψ-DrD=reD=0constantpressure,∂ψ-D∂rDrD=reD=0closed.The solution of diffusion equation (27) is(32)∂c-D∂riDriD=1=c-Dλucoth⁡λu-1.Substituting the dimensionless definitioncD and Langmuir isothermal adsorption equation V=VLPm/(PL+Pm) in (32), it can be written as follows:(33)∂c-D∂riDrD=1=LVLPPL+P-VLPiPL+Piλucoth⁡λu-1.SubstitutingψD=(ψi-ψ)/ψiqD in (33), it can be expressed as follows:(34)VLPPL+P-VLPiPL+Pi=-ψLVLψiqDψL+ψψL+ψiψ+ψiψ-D.Definingβ=ψLVLψiqD/ψL+ψψL+ψiψ+ψi,(35)∂c-D∂riDriD=1=-βψ-Dλucoth⁡λu-1.Substituting (35) in (28),(36)1rD∂∂rDrD∂ψ-D∂rD=fuψ-D,fu=ωu+1-ωλβλucoth⁡λu-1.


## 4. Equation of Boundary Condition

### 4.1. Boundary Integral Equation

Applying the theory of boundary element, previous equation can be solved from the integral transformation of the governing equation based on the expression of the fundamental solutions and differential equation of fluid flow through porous medium:(37)∫ΩψD(P,u)∇2G(P,Q,u)-G(P,Q,u)∇2ψD(P,u)∑i=1NwqDiδ(xD-xDi,yD-yDi)+δP,QψDP,u-1u∑i=1NwqDiδ(xD-xDi,yD-yDi)∑i=1NwqDiδxD-xDi,yD-yDi·GP,Q,udΩ=0,where G(P,Q,u) is the fundamental solution of horizontal wells in complex boundary reservoirs.According to the properties ofδ function and the second order of Green formula, it can be simplified to the boundary integral equation:(38)ψDQk,u=∫ΓGP,Qk,u∂ψDP′,u∂n-ψDP,u∂GP′,Qk,u∂ndΓP′+1u∑i=1Nw‍qDiGP,Qi,u.The boundaryΓ is divided into Nb cells, which are located at the end point and are taken as the nodes of boundary elements. Cell properties are assumed as linear distribution. Meanwhile, the boundary sections near nodes are assumed as arcs with nodes as their centers; the resulting boundary integral equation is(39)θkψD(Qk,u)=∑i=1Nb∫ΓiGP′,Qk,u∂ψDP′,u∂n-ψDP′,u∂GP′,Qk,u∂ndΓi(P′)+1u∑i=1NwqDiGP,Qi,u,where θk represents interior angles between any two adjacent boundary elements. Consider(40)θk=1,thepointindomainθi=2π,0.5,thepointatsmoothboundaryθi=π,θi2π,thepointatsmoothlessboundary.Using the linear interpolation in boundary element, the boundary integral formula is deformed as follows:(41)θkψD(Qk,u)=∑i=1Nbli2∫-11GP′,Qk,uφ1ξ∂ψDi∂n+φ2ξ∂ψDi+1∂n∂G(P′,Qk,u)∂n-(φ1(ξ)ψDi+φ2(ξ)ψDi+1)∂G(P′,Qk,u)∂ndξ+1u∑i=1NwqDiGP′,Qi,u,where ϕ1(ξ)=(1-ξ)/2 and ϕ2(ξ)=(1+ξ)/2 are linear interpolation formula, li=(xi+1-xi)2+(yi+1-yi)2, -1<ξ<1, and Γi is the length of linearity cell.


### 4.2. Fundamental Solution of Boundary Element Integral Equation

It is crucial to find its fundamental solution when the horizontal flow problem in complex reservoirs is resolved using the boundary element method. According to the properties of the boundary element and the mathematic equation governing pressure transmission in porous medium, the fundamental solution must satisfy the modified Helmholtz operator. The equation is given by(42)1r∂∂rr∂G∂r-fuG=-2πδMD,MD′.With Lord Kelvin point source solution, the fundamental solution of (42) can be derived as follows:(43)γ-=exp⁡-ρDfu4πρD.With mirror image method of fluid mechanics in porous medium, the transient point source fundamental solution of closed boundary is(44)γ-=14π∑-∞+∞exp⁡(-f(u)RD2+ZD+ZD′-2nZeD2)RD2+ZD+ZD′-2nZeD2+exp⁡(-f(u)RD2+ZD-ZD′-2nZeD2)RD2+ZD-ZD′-2nZeD2.With Poisson superposition formula, (44) can be simplified, and the transient point source fundamental solution of sealed boundary at Z=0 and Z=Ze is(45)γ-=12πZeDK0RDfuRDfu+n2π2ZeD2+2∑n=1n=∞K0RDfu+n2π2ZeD212πZeDRDfu+n2π2ZeD2·cos⁡⁡nπZDZeDcos⁡⁡nπZD′ZeD.The fundamental boundary element solution of horizontal wells in a reservoir with closed top and bottom boundaries is(46)GP′,Q,u=12∫-11K0RDfudα+∑n=1n=∞cos⁡(nπzD)cos⁡(nπzwD)·∫-11K0xD-α2+yD2fu+n2π2ZeD2dα,∂GP′,Q,u∂n=-12∫-11fuK1rD-αfu∂rD∂ndα-∑n=1n=∞cos⁡(nπzD)·cos⁡⁡nπzwD∫-11fu+n2π2ZeD2mmmmmmmmmm·K1fu+n2π2ZeD2(rD-α)mmmmmmmmmmmmm·fu+n2π2ZeD2∂rD∂ndα,where(47)∂rD∂n=±(xξ-xyi-yi+1-yξ-yxi-xi+1)(xi-xi+1)2+(yi-yi+1)2-1·(xi-xi+1)2+(yi-yi+1)2-1·(xξ-xi)2+(yξ-yi)2-1.Ifn and P′ are in the same direction, then ∂rD/∂n is positive; otherwise it is negative.


### 4.3. Solving of Integral Equation

The integral boundary equation (37) can be simplified as(48)θkψDQk,u=∑i=1NbHk1′∂ψDi∂n+Hk2′∂ψDi+1∂n∂ψDi∂n+Hk3′ψDi+Hk4′ψDi+1+1u∑i=1Nw‍qDiGP′,Q,u,where (49)H1′=li2∫-11GP′,Qk,uϕ1ξdξ;H2′=li2∫-11GP′,Qk,uϕ2ξdξ;H3′=li2∫-11-∂GP′,Qk,u∂nϕ1ξdξ;H4′=li2∫-11-∂GP′,Qk,u∂nϕ2ξdξ.The unknown variable in the integral boundary equation is∂ψDi/∂n or ψDi. The number of the nodes at the boundary is Nb, so Nb equation with form of (48) can be established. When the boundary properties are known, there are just Nb unknown variables. So we can solve the set of equations, whose matrix expression is(50)H11H12⋯H1NbH21H22⋯H2Nb⋮⋮⋮⋮HNb1HNb2⋯HNbNbx1x2⋮xNb=F1F2⋮FNb,where xi is ∂ψDi/∂n or ψDi and Fi is (1/u)∑i=1NwqDiG(P′,Q,u).Once the unknown variables are acquired, we can solvePD of arbitrary point in the research domain using the boundary integral equation (50):(51)ψDQ,u=∑i=1NbHk1′∂ψDi∂n+Hk2′∂ψDi+1∂n+Hk3′ψDi∂ψDi∂n+Hk4′ψDi+1+1u∑i=1NwqDiGP′,Q,u.


## 4.1. Boundary Integral Equation

Applying the theory of boundary element, previous equation can be solved from the integral transformation of the governing equation based on the expression of the fundamental solutions and differential equation of fluid flow through porous medium:(37)∫ΩψD(P,u)∇2G(P,Q,u)-G(P,Q,u)∇2ψD(P,u)∑i=1NwqDiδ(xD-xDi,yD-yDi)+δP,QψDP,u-1u∑i=1NwqDiδ(xD-xDi,yD-yDi)∑i=1NwqDiδxD-xDi,yD-yDi·GP,Q,udΩ=0,where G(P,Q,u) is the fundamental solution of horizontal wells in complex boundary reservoirs.According to the properties ofδ function and the second order of Green formula, it can be simplified to the boundary integral equation:(38)ψDQk,u=∫ΓGP,Qk,u∂ψDP′,u∂n-ψDP,u∂GP′,Qk,u∂ndΓP′+1u∑i=1Nw‍qDiGP,Qi,u.The boundaryΓ is divided into Nb cells, which are located at the end point and are taken as the nodes of boundary elements. Cell properties are assumed as linear distribution. Meanwhile, the boundary sections near nodes are assumed as arcs with nodes as their centers; the resulting boundary integral equation is(39)θkψD(Qk,u)=∑i=1Nb∫ΓiGP′,Qk,u∂ψDP′,u∂n-ψDP′,u∂GP′,Qk,u∂ndΓi(P′)+1u∑i=1NwqDiGP,Qi,u,where θk represents interior angles between any two adjacent boundary elements. Consider(40)θk=1,thepointindomainθi=2π,0.5,thepointatsmoothboundaryθi=π,θi2π,thepointatsmoothlessboundary.Using the linear interpolation in boundary element, the boundary integral formula is deformed as follows:(41)θkψD(Qk,u)=∑i=1Nbli2∫-11GP′,Qk,uφ1ξ∂ψDi∂n+φ2ξ∂ψDi+1∂n∂G(P′,Qk,u)∂n-(φ1(ξ)ψDi+φ2(ξ)ψDi+1)∂G(P′,Qk,u)∂ndξ+1u∑i=1NwqDiGP′,Qi,u,where ϕ1(ξ)=(1-ξ)/2 and ϕ2(ξ)=(1+ξ)/2 are linear interpolation formula, li=(xi+1-xi)2+(yi+1-yi)2, -1<ξ<1, and Γi is the length of linearity cell.


## 4.2. Fundamental Solution of Boundary Element Integral Equation

It is crucial to find its fundamental solution when the horizontal flow problem in complex reservoirs is resolved using the boundary element method. According to the properties of the boundary element and the mathematic equation governing pressure transmission in porous medium, the fundamental solution must satisfy the modified Helmholtz operator. The equation is given by(42)1r∂∂rr∂G∂r-fuG=-2πδMD,MD′.With Lord Kelvin point source solution, the fundamental solution of (42) can be derived as follows:(43)γ-=exp⁡-ρDfu4πρD.With mirror image method of fluid mechanics in porous medium, the transient point source fundamental solution of closed boundary is(44)γ-=14π∑-∞+∞exp⁡(-f(u)RD2+ZD+ZD′-2nZeD2)RD2+ZD+ZD′-2nZeD2+exp⁡(-f(u)RD2+ZD-ZD′-2nZeD2)RD2+ZD-ZD′-2nZeD2.With Poisson superposition formula, (44) can be simplified, and the transient point source fundamental solution of sealed boundary at Z=0 and Z=Ze is(45)γ-=12πZeDK0RDfuRDfu+n2π2ZeD2+2∑n=1n=∞K0RDfu+n2π2ZeD212πZeDRDfu+n2π2ZeD2·cos⁡⁡nπZDZeDcos⁡⁡nπZD′ZeD.The fundamental boundary element solution of horizontal wells in a reservoir with closed top and bottom boundaries is(46)GP′,Q,u=12∫-11K0RDfudα+∑n=1n=∞cos⁡(nπzD)cos⁡(nπzwD)·∫-11K0xD-α2+yD2fu+n2π2ZeD2dα,∂GP′,Q,u∂n=-12∫-11fuK1rD-αfu∂rD∂ndα-∑n=1n=∞cos⁡(nπzD)·cos⁡⁡nπzwD∫-11fu+n2π2ZeD2mmmmmmmmmm·K1fu+n2π2ZeD2(rD-α)mmmmmmmmmmmmm·fu+n2π2ZeD2∂rD∂ndα,where(47)∂rD∂n=±(xξ-xyi-yi+1-yξ-yxi-xi+1)(xi-xi+1)2+(yi-yi+1)2-1·(xi-xi+1)2+(yi-yi+1)2-1·(xξ-xi)2+(yξ-yi)2-1.Ifn and P′ are in the same direction, then ∂rD/∂n is positive; otherwise it is negative.


## 4.3. Solving of Integral Equation

The integral boundary equation (37) can be simplified as(48)θkψDQk,u=∑i=1NbHk1′∂ψDi∂n+Hk2′∂ψDi+1∂n∂ψDi∂n+Hk3′ψDi+Hk4′ψDi+1+1u∑i=1Nw‍qDiGP′,Q,u,where (49)H1′=li2∫-11GP′,Qk,uϕ1ξdξ;H2′=li2∫-11GP′,Qk,uϕ2ξdξ;H3′=li2∫-11-∂GP′,Qk,u∂nϕ1ξdξ;H4′=li2∫-11-∂GP′,Qk,u∂nϕ2ξdξ.The unknown variable in the integral boundary equation is∂ψDi/∂n or ψDi. The number of the nodes at the boundary is Nb, so Nb equation with form of (48) can be established. When the boundary properties are known, there are just Nb unknown variables. So we can solve the set of equations, whose matrix expression is(50)H11H12⋯H1NbH21H22⋯H2Nb⋮⋮⋮⋮HNb1HNb2⋯HNbNbx1x2⋮xNb=F1F2⋮FNb,where xi is ∂ψDi/∂n or ψDi and Fi is (1/u)∑i=1NwqDiG(P′,Q,u).Once the unknown variables are acquired, we can solvePD of arbitrary point in the research domain using the boundary integral equation (50):(51)ψDQ,u=∑i=1NbHk1′∂ψDi∂n+Hk2′∂ψDi+1∂n+Hk3′ψDi∂ψDi∂n+Hk4′ψDi+1+1u∑i=1NwqDiGP′,Q,u.


## 5. Validation of the Model

In order to validate the gas flow model proposed in this paper, we test the adsorption and desorption data of coal. Then, we compare the boundary element model with those published works.


### 5.1. Validation of the Gas Flow Model in CBM Reservoirs

The adsorption and desorption data of coal are tested in the laboratory; the pressure continues to drop as the coal continues to adsorb methane. We deal with the experimental data in log-log coordinates as shown in Figure4. It can be seen from the diagram that the results of the numerical simulation are in good match with the experimental data and are consistent with the gas flow law in CBM reservoirs.Figure 4
Double logarithmic curve of pressure drop of test data.


### 5.2. Validation of Gas Flow in a Horizontal Well with Complex Boundary

A pressure buildup test example of a horizontal gas well in tight gas reservoir has been presented by Han [30]; the production and pressure data fitting results show that the horizontal gas well has a closed boundary.Figure5 is the comparison of pressure curves between horizontal wells in CBM and conventional gas reservoirs, and the values of parameters are from Han’s paper [30]. According to the analysis of flow characteristics, there is an additional radial flow which reflects gas flow in fracture. After this stage, the pressure transmits from fracture to matrix and gas starts to desorb. The flow characteristics of this stage would be influenced by desorption velocity and desorption quantity. The last stage of pressure derivative curves in different reservoirs is in good match, which reflects the type of boundary.Figure 5
Comparison of typical curves between horizontal wells in CBM and conventional gas reservoirs.


## 5.1. Validation of the Gas Flow Model in CBM Reservoirs

The adsorption and desorption data of coal are tested in the laboratory; the pressure continues to drop as the coal continues to adsorb methane. We deal with the experimental data in log-log coordinates as shown in Figure4. It can be seen from the diagram that the results of the numerical simulation are in good match with the experimental data and are consistent with the gas flow law in CBM reservoirs.Figure 4
Double logarithmic curve of pressure drop of test data.


## 5.2. Validation of Gas Flow in a Horizontal Well with Complex Boundary

A pressure buildup test example of a horizontal gas well in tight gas reservoir has been presented by Han [30]; the production and pressure data fitting results show that the horizontal gas well has a closed boundary.Figure5 is the comparison of pressure curves between horizontal wells in CBM and conventional gas reservoirs, and the values of parameters are from Han’s paper [30]. According to the analysis of flow characteristics, there is an additional radial flow which reflects gas flow in fracture. After this stage, the pressure transmits from fracture to matrix and gas starts to desorb. The flow characteristics of this stage would be influenced by desorption velocity and desorption quantity. The last stage of pressure derivative curves in different reservoirs is in good match, which reflects the type of boundary.Figure 5
Comparison of typical curves between horizontal wells in CBM and conventional gas reservoirs.


## 6. Analysis of Flow Characteristics and Field Application

### 6.1. Flow Regime Recognition

Type curves for horizontal wells in a CBM reservoir with complex boundary calculated by the BEM in this paper are shown in Figure6. The figure shows that flow characteristics can be divided into seven flow periods: (1) wellbore storage period; influenced by early wellbore storage effect, the curves of pressure and pressure derivative are straight lines with the same slope in this period; (2) the first transition flow section, which reflects the degree of pollution near the bottom and is mainly affected by skin factor S; (3) the first radial flow section, which is perpendicular to the horizontal wellbore; the pressure derivative curve is shown as a horizontal line, which reflects the early radial flow perpendicular to the horizontal wellbore; (4) the second transition section, which is the transition stage of radial flow from wellbore to fracture; (5) the second radial flow stage, which reflects the free gas flow in fracture; this period is affected by the flow of free gas and adsorbed gas, diffusion coefficient λ, and Langmuir adsorption parameter β; (6) the third radial flow section, which reflects radial flow of whole system after pressure balance, and pressure derivative curve appears as a flat behavior; (7) boundary response section. The pressure and pressure derivative curves with closed boundary would be upward after pressure transmitting to the boundary. The pressure derivative curve with constant pressure boundary would be in decline after pressure transmitting to the boundary. Compared with the constant pressure boundary, the falling range of pressure derivative curve with mixed boundary is less.Figure 6
Typical curve of CBM horizontal well with complex boundary conditions (ω=0.01, λ=300, β=50, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).


### 6.2. Parameter Sensitivity Analysis

Figure7 shows the influence of diffusion coefficient λ on the bottom-hole pressure with mixed boundary. There is relevance between diffusion coefficient λ=αkτ/θL2 and adsorption time τ=R2/D. The smaller the τ, the shorter the time of desorption-diffusion and the shorter the time of pressure balance between fracture and matrix. As shown in the diagram, diffusion coefficient λ mainly affects the duration of the second radial flow period and appearance of third radial flow. The greater the diffusion coefficient λ, the longer the duration of second radial flow period and the later the t appearance of third radial flow, and vice versa. If λ is small enough, the second radial flow period would be covered.Figure 7
Influence of diffusion coefficientλ on bottom-hole pressure with mixed boundary (ω=0.01, β=50, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).Figure8 shows the influence of Langmuir parameter β on the bottom-hole pressure with mixed boundary. β=ψLVLψiqD/ψL+ψψL+ψiψ+ψi is related to Langmuir adsorption pressure pL and Langmuir adsorption volume VL. The larger the β, the greater the adsorption ability of coal. As shown in the graph, β only affects the second radial flow section. The greater the β, the smaller the pressure derivative values of second radial flow period.Figure 8
Influence ofβ on the bottom-hole pressure with mixed boundary (ω=0.01, λ=300, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).Figure9 shows the influence of fracture storage ratio ω on the bottom-hole pressure with mixed boundary. Its expression ω=ϕfctμ/(ϕfctμ+6pscT/qDTscψi) indicates that the smaller the ϕfct, the smaller the ω and the more the radial flow in fracture. From Figure 8, almost all of the flow periods would be affected by ω, except for wellbore storage period. The greater the ω, the smaller the hump value of pressure derivative, and the longer the duration of first radial flow period, the greater the pressure derivative value in second radial flow period.Figure 9
Influence ofω on the bottom-hole pressure with mixed boundary (λ=300, β=50, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).Figure10 shows the influence of eccentricity of horizontal well Zwd on the bottom-hole pressure with mixed boundary. The streamline shape and pressure distribution of a horizontal well would be affected by the size of Zwd. The smaller the Zwd, the larger the pressure derivative value in first radial flow period.Figure 10
Influence ofZwd on the bottom-hole pressure with mixed boundary (ω=0.01, λ=300, β=50, LD=2.5, Cd=0.0001, s=1, and reD=50).Figure11 shows the influence of length of horizontal well LD on the bottom-hole pressure with mixed boundary. It shows large influence of LD on the pressure derivative value in first radial flow period. The smaller the LD, the larger the pressure derivative value in first radial flow period.Figure 11
Influence ofLD on the bottom-hole pressure with mixed boundary (ω=0.01, λ=300, β=50, Zwd=0.5, Cd=0.0001, s=1, and reD=50).


### 6.3. Field Example

To demonstrate the application of the model proposed in this paper, a horizontal well ZX-12 in a CBM reservoir is studied. The vertical depth and horizontal length of this well are 689 m and 536 m, respectively, with well radiusrw of 0.1 m, coal thickness of 4.5 m, and initial pressure of 5.5 MPa. According to test results, the Langmuir volume VL is 32.75 m3/t and the Langmuir pressure PL is 2.49 MPa. Figure 12 is the fitted curves which were obtained by calculating the actual production and pressure data. As shown in Figure 12, the field data are in good match with the theoretical results calculated by the new model in this paper.Figure 12
Fitted curves of test data from an actual well.According to the actual field data, some of reservoir parameters are fitted as follows: the porosity is 0.04 and the permeability is 3 mD. The stage of desorption and boundary response can be clearly seen in Figure12. Because of the falling of the pressure derivative curve in the last stage, well ZX-12 has a constant pressure boundary, implying that the gas production of well ZX-12 is affected by the production of adjacent well.


## 6.1. Flow Regime Recognition

Type curves for horizontal wells in a CBM reservoir with complex boundary calculated by the BEM in this paper are shown in Figure6. The figure shows that flow characteristics can be divided into seven flow periods: (1) wellbore storage period; influenced by early wellbore storage effect, the curves of pressure and pressure derivative are straight lines with the same slope in this period; (2) the first transition flow section, which reflects the degree of pollution near the bottom and is mainly affected by skin factor S; (3) the first radial flow section, which is perpendicular to the horizontal wellbore; the pressure derivative curve is shown as a horizontal line, which reflects the early radial flow perpendicular to the horizontal wellbore; (4) the second transition section, which is the transition stage of radial flow from wellbore to fracture; (5) the second radial flow stage, which reflects the free gas flow in fracture; this period is affected by the flow of free gas and adsorbed gas, diffusion coefficient λ, and Langmuir adsorption parameter β; (6) the third radial flow section, which reflects radial flow of whole system after pressure balance, and pressure derivative curve appears as a flat behavior; (7) boundary response section. The pressure and pressure derivative curves with closed boundary would be upward after pressure transmitting to the boundary. The pressure derivative curve with constant pressure boundary would be in decline after pressure transmitting to the boundary. Compared with the constant pressure boundary, the falling range of pressure derivative curve with mixed boundary is less.Figure 6
Typical curve of CBM horizontal well with complex boundary conditions (ω=0.01, λ=300, β=50, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).


## 6.2. Parameter Sensitivity Analysis

Figure7 shows the influence of diffusion coefficient λ on the bottom-hole pressure with mixed boundary. There is relevance between diffusion coefficient λ=αkτ/θL2 and adsorption time τ=R2/D. The smaller the τ, the shorter the time of desorption-diffusion and the shorter the time of pressure balance between fracture and matrix. As shown in the diagram, diffusion coefficient λ mainly affects the duration of the second radial flow period and appearance of third radial flow. The greater the diffusion coefficient λ, the longer the duration of second radial flow period and the later the t appearance of third radial flow, and vice versa. If λ is small enough, the second radial flow period would be covered.Figure 7
Influence of diffusion coefficientλ on bottom-hole pressure with mixed boundary (ω=0.01, β=50, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).Figure8 shows the influence of Langmuir parameter β on the bottom-hole pressure with mixed boundary. β=ψLVLψiqD/ψL+ψψL+ψiψ+ψi is related to Langmuir adsorption pressure pL and Langmuir adsorption volume VL. The larger the β, the greater the adsorption ability of coal. As shown in the graph, β only affects the second radial flow section. The greater the β, the smaller the pressure derivative values of second radial flow period.Figure 8
Influence ofβ on the bottom-hole pressure with mixed boundary (ω=0.01, λ=300, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).Figure9 shows the influence of fracture storage ratio ω on the bottom-hole pressure with mixed boundary. Its expression ω=ϕfctμ/(ϕfctμ+6pscT/qDTscψi) indicates that the smaller the ϕfct, the smaller the ω and the more the radial flow in fracture. From Figure 8, almost all of the flow periods would be affected by ω, except for wellbore storage period. The greater the ω, the smaller the hump value of pressure derivative, and the longer the duration of first radial flow period, the greater the pressure derivative value in second radial flow period.Figure 9
Influence ofω on the bottom-hole pressure with mixed boundary (λ=300, β=50, LD=2.5, Zwd=0.5, Cd=0.0001, s=1, and reD=50).Figure10 shows the influence of eccentricity of horizontal well Zwd on the bottom-hole pressure with mixed boundary. The streamline shape and pressure distribution of a horizontal well would be affected by the size of Zwd. The smaller the Zwd, the larger the pressure derivative value in first radial flow period.Figure 10
Influence ofZwd on the bottom-hole pressure with mixed boundary (ω=0.01, λ=300, β=50, LD=2.5, Cd=0.0001, s=1, and reD=50).Figure11 shows the influence of length of horizontal well LD on the bottom-hole pressure with mixed boundary. It shows large influence of LD on the pressure derivative value in first radial flow period. The smaller the LD, the larger the pressure derivative value in first radial flow period.Figure 11
Influence ofLD on the bottom-hole pressure with mixed boundary (ω=0.01, λ=300, β=50, Zwd=0.5, Cd=0.0001, s=1, and reD=50).


## 6.3. Field Example

To demonstrate the application of the model proposed in this paper, a horizontal well ZX-12 in a CBM reservoir is studied. The vertical depth and horizontal length of this well are 689 m and 536 m, respectively, with well radiusrw of 0.1 m, coal thickness of 4.5 m, and initial pressure of 5.5 MPa. According to test results, the Langmuir volume VL is 32.75 m3/t and the Langmuir pressure PL is 2.49 MPa. Figure 12 is the fitted curves which were obtained by calculating the actual production and pressure data. As shown in Figure 12, the field data are in good match with the theoretical results calculated by the new model in this paper.Figure 12
Fitted curves of test data from an actual well.According to the actual field data, some of reservoir parameters are fitted as follows: the porosity is 0.04 and the permeability is 3 mD. The stage of desorption and boundary response can be clearly seen in Figure12. Because of the falling of the pressure derivative curve in the last stage, well ZX-12 has a constant pressure boundary, implying that the gas production of well ZX-12 is affected by the production of adjacent well.


## 7. Conclusions

The mathematic model of gas flowing into horizontal wells in CBM reservoirs with complex boundary conditions is derived based on the percolation theory. The curves of bottom-hole pressure and pressure derivative are obtained by using boundary element method and Laplace transform. The conclusions are as follows.(1)
The fundamental boundary element solutions for transient pressure response of horizontal wells in CBM reservoirs with complex boundary conditions could be obtained by using Lord Kelvin’s point source solution, point source function theory, and Poisson’s summation formula.(2)
Comparison of the typical curves of flow characteristics between horizontal wells in CBM and conventional gas reservoirs shows that there is an additional radial flow which reflects gas flow in fracture.(3)
Comparing the typical curves of flow characteristics with complex boundary conditions, the pressure and pressure derivative curves with closed boundary would be upward after pressure transmitting to the boundary. The pressure derivative curves with constant pressure boundary and mixed boundary would be fallen, but the falling range of pressure derivative curve with mixed boundary is less.



---

*Source: 173975-2015-09-27.xml*