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 () 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.
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 Section 2, 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 () 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 where is the gas pressure, is the gas volume, is the amount of substance, is the gas temperature, is the gas mass, is the gas molar mass, is the gas deviation factor, and is the universal gas constant. From (1) the gas density can be written as
According to the Langmuir isotherm and the state equation of real gas, the adsorbed gas concentration in coal matrix and free gas concentration in fracture systems can be written, respectively, as where is the gas pressure in coal matrix, is the maximum adsorption concentration of the coal matrix, is the Langmuir pressure, is the gas pressure in fracture systems, and is the porosity of fracture systems.
2.2. Gas Flow Model
Considering the effects of fracture system deformation, fracture system permeability tensor can be written as where is the initial gas pressure in fracture systems, is the fracture system permeability tensor under , 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 where is the total volume velocity of gas flow driven by multifield, is the gas viscosity, and 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 where 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: Similarly, the capacity term could be unfolded as the following form: where is the compressibility factor of and is the compressibility factor of , defined as For the expression brevity, define and substitute motion equation (5) into (7); the continuity equation (6) can be rewritten as follows: Since the seepage velocity is very low, the quadratic term usually can be omitted in engineering application. Therefore, the more common form of continuity equation is as follows: In addition, there are some other studies [26, 27] considering the influence of quadratic term . 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 where is the average gas concentration of a coal matrix block, is the adsorbed concentration on the coal matrix surface and is defined as , 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: where is the objective region for the two-dimensional problem with boundary , is a time interval with , and 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: Recall that the Sobolev space is the closure of in the norm Define the function spaces , and their norms as follows: Let be the mixed element function spaces such as with index and discretization parameter . Now we define the projection , satisfying Define standard projection , satisfying
3.2. Mixed Element Approximation Scheme
We consider the case that the seepage-diffusion tensor is a diagonal matrix and for the two-dimensional problem. That is to say, there exist positive constants and such that, for all , For the practical situation, there exist and such that the scalar coefficient satisfies Also, we can easily verify that the coefficients , , and satisfy the Lipschitz continuity. Since is a diagonal and invertible matrix, from (5) we have that
Using the integration by parts and applying the boundary condition, we can define the mixed weak formation, which is to find , such that where is the unit exterior normal vector to the .
Replacing the original pressure and velocity by their approximations, we get the semidiscrete mixed element approximate problem, which is to find , such that
Let , , an integer, and , . We can define the full discrete mixed element approximate scheme with backward Euler time-discretization as follows: Since gas diffusion equation (13) is an ordinary differential equation, when is known, we can get as follows: Then, the gas volume of interporosity flow from coal matrix to fracture systems could be calculated as follows: In practical calculation, alternative and iterative method is used to solve (25)–(27), and specific calculation procedure is listed as follows. (i)For , when is known, can be calculated using (27).(ii)Then, substituting into (25), can be approximated using (25).(iii)Lastly, when is known, can be calculated using (26).
Let be a quasiregular triangulation for the two-dimensional rectangular region with mesh size . For a triangular element , let be the coordinate of vertex and the edge opposite to vertex . Let be the outward unit vector on the edge and the length of the perpendicular dropped from the vertex onto the edge , . We denote by the set of all edges of the triangulation .
The pressure function space is the space of piecewise constant functions: and the dimension of the space is equal to —the total number of elements in the triangulation . Each function can be represented as a linear combination: where is the basis function associated with the element , satisfying
The velocity function space is chosen to be the , defined as and the dimension of the space is equal to —the total number of edges . Each function can be represented as a linear combination where and is the basis function associated with the edge , satisfying To be specific, for a triangle element , the basis functions are determined by the following formulas: where are the coordinates of the vertices and is the area of triangle element .
Now we use the basis functions and to represent and : Substituting (35) into (25) and setting and , full discrete mixed element approximation scheme (25) can be written as follows:
We introduce matrices and vectors as follows: where Therefore, (36) can be written in a matrix form as follows: where is the transpose of . Solving linear algebra problem (39), we can get approximation solutions and at th time step, simultaneously.
4. Numerical Examples
In this section, we carry out numerical experiments using the lowest order Raviart-Thomas mixed finite elements 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, , and 6 levels for , , , , , are computed to estimate the convergence rate. Since we cannot get the analytical solution , , and for the problem (14), we use , , and as the criterion of convergence for the pressure, velocity, and concentration. Error estimates in -norm and convergence rate estimates for , , and are listed in Tables 1, 2, and 3, respectively. We can find that the approximation solutions , , and could achieve first-order convergence rate (see Figure 1).
Example 2. For a simple application, the scale of objective region is m and the mesh size m. There is a production well at the center of region with the production rate . 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, Figures 2 and 3 illustrate the gas average pressure in fracture systems varying with time and show the effects of and on average pressure, respectively. In Figure 2, the model has the initial reservoir pressure of 10 MPa, and diffusion coefficient of coal matrix is taken as , , , and m3/s, respectively, so as to analyze the effects of 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 is increasing, the initial stage became shorter. It is because that larger 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 , 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 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 Figure 3, the permeability modulus is set as , , and , 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 . Thus, when is larger, decreases more significantly with dropping, such that CBM reservoirs can only get less pressure compensation from outer boundary, resulting in a lower average pressure.
Figure 4 shows the bottom hole pressure varying with time , and the diffusion coefficient of coal matrix is taken as , , , and /s, respectively, so as to analyze the effects of 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 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 become weak and bottom hole pressure curves tend to be close. Put simply, the effects of 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 Figure 5, we set fracture diffusion coefficient as , , and m3/s, respectively, so as to analyze the effects of gas diffusion in fracture systems. We have found that when is larger, the bottom hole pressure is higher. Specifically, the pressure difference between and 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 time in Figures 6 and 7. In Figure 6, the diffusion coefficient is taken as , , , and 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 in matrix, fracture pressure , and diffusion coefficient . During production process, as fracture pressure and gas average concentration decline, there are two stages in gas desorption process. At first half stage, due to quick drawdown of fracture pressure, gas concentration 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 . In Figure 6, we have also found when is larger, average desorption rate will be higher and peak desorption rate appears at an earlier moment. In Figure 7, the fracture diffusion coefficient is taken as , , and m3/s, respectively. Since the effects of on desorption rate are indirect and reflected by means of fracture pressure , the effects of are not evident at the starting and ending time. However, we have still found that the larger the 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 time in Figures 8 and 9. In Figure 8, is taken as , , , and 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 is larger, the average concentration is lower and the concentration drop speed is quicker. It is because that larger could accelerate gas desorption rate. In Figure 9, we have found that the effects of initial permeability 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 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.
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 elements . Numerical results in Example 1 indicate that the approximation solutions , , and 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.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The work described in this paper was a part of research work of project supported by National Natural Science Foundation of China (Grant nos. 91330106 and 11171190).