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

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 CBMhorizontal 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.


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][8][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 2 Journal of Chemistry 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][15][16] established a two-dimensional, twophase, multiwell gas-water flow model, and this model has the ability to simulate multiple horizontal wells.Engler and Rajtar [17][18][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 threedimensional 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.

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 (Figure 1(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.
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 Figure 2. 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].
To facilitate the derivation, it is assumed that the length of a horizontal well is , the center location of the horizontal well is (  ,   ,   ), and the thickness of a CBM reservoir is ℎ.In this paper, gas flow into horizontal wells with closed boundary (Figure 3 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.

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   The planar radial flow equation is as follows: The dimensionless equation of gas diffusion in matrix is where   is the dimensionless radius defined by   =   /;   is the dimensionless diffusion concentration defined by   =  −   ;   is the dimensionless time defined by   = 3.6/  2 ;  1 is interporosity flow coefficient defined by  1 = 3.6 2 /  2 ;  is comprehensive storage coefficient defined by  =     +6  /      ;   is the dimensionless production defined by   = 1.842×10 −3     /ℎ    ;   is the initial concentration of gas, kg/m 3 ;  is the permeability, m 2 ; and   is the pseudopressure.

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: Concentration of the adsorbed gas is Gas concentration in coal is Volume of gas desorption from coal is where   is the gas mass under standard condition, kg/m 3 ;   is the volume of gas adsorption in coal under standard condition, m 3 /m 3 .  is the density of gas under standard condition, kg/m 3 .Density of gas is defined as Spherical matrix has the following relationship: where  is the volume of coal matrix, m 3 ;  is the surface area of coal matrix, m 2 ;  is the mass diffusion coefficient, m 2 /s;  is the diffusion flux, g/m 2 s.Combining ( 7), (8), and ( 9), the volume of gas desorption is given by Based on the basic principle of material balance, the governing equations of gas flow in fracture system are given by Equation ( 11) is simplified into the following equation: where  is the fracture storage ratio given by  = (  )/;  is the diffusion coefficient given by  = 3.6/  2 ;   is the dimensionless pseudopressure given by   = (  − )/    ;  is the adsorption time given by  =  2 /.

Mathematical Model Solution.
Diffusion and governing equations of gas flow in fracture are combined as follows: Dimensionless initial and boundary conditions are given as follows.
The initial condition is The boundary condition is The infinite outer boundary condition is The constant pressure boundary condition is The closed boundary condition is When  =     , ( 13) can be transformed into The Laplace transform of ( 19) is The general solution of ( 20) is where sinh() = (  −  − )/2 and cosh() = (  +  − )/2.The coefficients  and  could be obtained by the initial and boundary conditions: Hence, Therefore, the diffusion equation is transformed into Combining the dimensionless definition   and Langmuir isothermal adsorption equation  =     /(  +   ), (24) is transformed into where  is Laplace transform: So, Laplace transform of ( 20) is The initial conditions are The inner boundary conditions are The outer boundary conditions are The solution of diffusion equation ( 27) is Substituting the dimensionless definition   and Langmuir isothermal adsorption equation  =     /(  +   ) in (32), it can be written as follows: Substituting   = (  − )/    in (33), it can be expressed as follows: Substituting ( 35) in (28),

Equation of Boundary Condition
where (, , ) 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: The boundary Γ is divided into   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 where   represents interior angles between any two adjacent boundary elements.Consider the point in domain   = 2, 0.5, the point at smooth boundary   = ,   2 , the point at smoothless boundary. (40) Using the linear interpolation in boundary element, the boundary integral formula is deformed as follows: where , and Γ  is the length of linearity cell.

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 With Lord Kelvin point source solution, the fundamental solution of (42) can be derived as follows: With mirror image method of fluid mechanics in porous medium, the transient point source fundamental solution of closed boundary is With Poisson superposition formula, (44) can be simplified, and the transient point source fundamental solution of sealed boundary at  = 0 and  =  is The fundamental boundary element solution of horizontal wells in a reservoir with closed top and bottom boundaries is where If  and   are in the same direction, then   / is positive; otherwise it is negative..The integral boundary equation (37) can be simplified as

Solving of Integral Equation
where The unknown variable in the integral boundary equation is   / or   .The number of the nodes at the boundary is   , so   equation with form of (48) can be established.When the boundary properties are known, there are just   unknown variables.So we can solve the set of equations, whose matrix expression is where   is   / or   and   is (1/) ∑  =1   (  , , ).Once the unknown variables are acquired, we can solve   of arbitrary point in the research domain using the boundary integral equation ( 50): (  , , ) . (51)

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.

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 Figure 4.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.

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.Figure 5 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.

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 Figure 6.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 ; (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.

Parameter Sensitivity Analysis.
Figure 7 shows the influence of diffusion coefficient  on the bottom-hole pressure with mixed boundary.There is relevance between diffusion coefficient  = / 2 and adsorption time  =  2 /.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  appearance of third radial flow, and vice versa.If  is small enough, the second radial flow period would be covered.Figure 8 shows the influence of Langmuir parameter  on the bottom-hole pressure with mixed boundary. =         /(  + )(  +   )( +   ) is related to Langmuir adsorption pressure   and Langmuir adsorption volume   .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 9 shows the influence of fracture storage ratio  on the bottom-hole pressure with mixed boundary.Its expression  =     /(     + 6  /      ) indicates that the smaller the     , 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 10 shows the influence of eccentricity of horizontal well   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   .The smaller the   , the larger the pressure derivative value in first radial flow period.
Figure 11 shows the influence of length of horizontal well   on the bottom-hole pressure with mixed boundary.It shows large influence of   on the pressure derivative value in first radial flow period.The smaller the   , the larger the pressure derivative value in first radial flow period.

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 radius   of 0.1 m, coal thickness of 4.5 m, and initial pressure of 5.5 MPa.According to test results, the Langmuir volume   is 32.75 m 3 /t and the Langmuir pressure   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.
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 Figure 12.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.

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 bottomhole pressure and pressure derivative are obtained by using boundary element method and Laplace transform.The conclusions are as follows.

Figure 2 :
Figure 2: Desorption and diffusion in the CBM reservoirs.

Figure 4 :
Figure 4: Double logarithmic curve of pressure drop of test data.

Figure 5 :Figure 6 :
Figure 5: Comparison of typical curves between horizontal wells in CBM and conventional gas reservoirs.

Figure 12 :
Figure 12: Fitted curves of test data from an actual well.

( 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.Nomenclature   : Radius, m : Permeability, mD : Porosity, fraction : Diffusion flux, g/m 2 s : Laplace transform variable ℎ: Thickness, m : Volume factor, fraction : Influx into the wellbore, m 3 /d : Half length of horizontal well, m : Viscosity, mPa⋅s : Concentration of coalbed methane, kg/m 3   : Initial pressure, MPa : Production of well, m 3 /d   : Langmuir volume constant, m 3 /ton : Volume of coal matrix, m 3   : Pseudopressure : Fracture storage ratio, fraction.