Abstract

Particle scale diffusion is implemented in the discrete element code, Esys-Particle. We focus on the question of how to calibrate the particle scale diffusion coefficient. For the regular 2D packing, theoretical relation between micro- and macrodiffusion coefficients is derived. This relation is then verified in several numerical tests where the macroscopic diffusion coefficient is determined numerically based on the half-time of a desorption scheme. To further test the coupled model, we simulate the diffusion and desorption in the circular sample. The numerical results match the analytical solution very well. An example of gas diffusion and desorption during sample crushing and fragmenting is given at the last. The current approach is the first step towards a realistic and comprehensive modelling of coal and gas outbursts.

1. Introduction

The migration of gas (i.e., methane or CO2) in coal plays an important role in the process of coalbed methane recovery or gas drainage prior to coal mining. It is generally accepted that gas transport in coal mainly occurs in two stages: gas flow within the coal matrix and flow in the cleat system which forms the natural fractures in coal [14]. Flow through the cleat system (or the fracture system) is believed to be pressure-driven laminar flow and can be described by Darcy’s law. Flow through the matrix is a diffusion process, which can be described by Fick’s diffusion law where the concentration gradient in the matrix is the driving force. From a microscopic point of view, the two kinds of flow patterns start from the desorption of gas which was initially adsorbed in the coal matrix, followed by the diffusion through the matrix into the cleat and then the flow through the cleat system into a production well or a drainage borehole. Hence, the gas production rate is mainly controlled by the gas diffusivity in the matrix and gas permeability in the cleat system.

Meanwhile, gas diffusion and gas flow in the cleat system are also found to be major contributing factors to the start of coal and gas outbursts. An outburst of coal and gas is the sudden release of a large amount of gas in conjunction with the ejection of coal and possibly associated with rocks. Previous studies have suggested that the major contributing factors of outbursts include stress condition, gassiness of coal seams, geological structures, and mechanical and physical properties of coal [510]. These important mechanisms should be implemented in a coal and gas outburst model which will permit certain parameters to be varied so that their effects on the outburst can be quantitatively studied. However most of the current models do not model solid fracture and fragmentation explicitly. Free flow of fluid and two-way interactions between the solid and fluid are also missing in most of the existing outburst models. To overcome these difficulties, we have developed a new outburst model which couples two well-developed numerical approaches: the discrete element method (DEM) and the lattice Boltzmann method (LBM) [11, 12]. The DEM is used to model the deformation and fracture of solid, while LBM simulates fluid flow, including free flow and Darcy flow. These two methods are coupled in a two-way process: the solid part provides a moving boundary condition and transfers momentum to the fluid, and the fluid exerts a dragging force to the solid. The new model includes the most recognized factors of outbursts, including deformation, fracture and fragmentation of solids, free flow of fluid, Darcy flow, diffusion, desorption of gas, and two-way coupling of solid and fluid. The preliminary results with small scale simulations suggest that the new model has potentials to numerically investigate the underlying mechanism and interaction of contributing factors of outbursts.

As the first step towards a realistic and comprehensive modelling for coal and gas outbursts, we mainly focus on the coupling of diffusion mechanism with the DEM model in this paper. Since the DEM models require some input parameters at particle scale which are not directly linked to the macroscopic parameters, calibration of the input parameters has to be carefully carried out. In this paper, we first introduce the implementation of diffusion mechanism into the DEM code and then discuss how to determine the particle scale diffusion coefficient so as to reproduce the macroscopic diffusion coefficient. Several numerical tests are carried out to verify the results by comparing with analytical solutions.

2. Implementation of Fluid Diffusion in DEM

The DEM is a widely used numerical tool to model the behaviour of rock and granular materials [13]. In DEM simulations, the specimen to be modelled is represented as an assembly of indivisible particles interacting with their nearest neighbours. By applying boundary conditions and solving Newton’s laws for each particle, the complex behaviour of the material may be simulated. The Esys-Particle, an open source DEM code, is used in this study. Detailed information about the model and the code can be found in the literature [14].

In this study, diffusion is implemented into the DEM code in the following way. Solid particles are treated as porous materials. It is assumed that the voids inside particle are much smaller than the particle sizes. Therefore the porosity is just an average concept for each particle. There is an average and uniform concentration for each particle . For two contacted particles and , the fluid exchange at each time step is determined by Fick’s first law of diffusion: where is the diffusion coefficient of the link at particle scale. For the regular packing in 2D case, we can derive theoretically the relation between and the macroscopic diffusion coefficient [15]:

3. Numerical Verifications of (2) with 1D Diffusion Simulation

The method described in the book of Crank [16] is adopted in this study to determine the macroscopic diffusion coefficient. In this approach, we assume a constant diffusion coefficient and consider the case of one-dimensional diffusion in a medium bounded by two parallel planes, for example, the planes at and . For diffusion through a plane sheet or membrane of thickness with diffusion coefficient , the diffusion equation is

If the boundary and initial conditions are the analytical solution of this problem is [16]

Let be the total amount of gas desorbed by the sheet at time and the corresponding amount theoretically after infinite time; we have

The value of for which , written as , is given by [16]

The macroscopic diffusion coefficient is found to be

Numerical tests have been carried out to calculate the macroscopic diffusion coefficients based on (8). In these tests, we use the regular packing which consists of 2006 particles with the same size of 1 unit and 5824 bonds (Figure 1). The width of the sheet in vertical direction is . We keep the initial concentration and the boundary condition (top), (bottom). Desorption can only occur on the top and bottom boundaries but is not allowed from the left and right boundaries. Therefore this is exactly a one-dimensional problem described in [16].

The fractional loss versus time for five different values, = 2, 6, 10, 50, and 100, are plotted in Figure 2. Note that larger input values of are used in this study since the small requires huge amount of computer time. After reading values from Figure 2, the macroscopic diffusion coefficient can be calculated based on (8). The relation between and is plotted in Figure 3, where the straight line is the analytical solution of (2), and dots are the simulation results. A very good agreement is found, which verifies the method described in (8).

4. Further Verifications of Circle (or Cylinder) Diffusion

To further test the coupled model, we simulate diffusion and desorption in a circular (or cylinder) sample. We assume that the cylinder length is infinite, diffusion coefficient () is constant, and the gas concentration () depends on the radial coordinate of the cylinder only; then the diffusion equation is given as

The initial and boundary conditions can be expressed as where is the radius of the cylinder, is the initial uniform concentration, and is the constant concentration at the surface of the cylinder. The analytical solution of (9) with the boundary and initial conditions above is [16] where is the Bessel function of the first kind of order zero and is the Bessel function of the first order, and , are roots of .

Similarly, the analytical solution of the fractional loss can be obtained [16]:

Several numerical tests have been carried out to further test the coupled model based on the calculated macroscopic diffusion coefficients in Section 3. Figure 4 shows the snapshots of concentration distribution at three different times for particle scale diffusion coefficient . Figure 5 presents the fractional loss versus . The solid line is the analytical solution of (12), and the dots are simulation results. Good agreements are observed between the numerical results and the analytical solutions.

5. Diffusion and Desorption during Fracturing Process

In coal mines, the coal sample may be crushed to quickly measure the desorbed gas content and residual gas. This process is modelled by moving two rigid loading walls towards each other in vertical direction. In this case neighboring particles are bonded and the bonds can break if the stress condition is reached, explicitly modelling microscopic fracturing events. Figure 6 shows the distribution of concentration during crushing process at = 10, 20, and 25. It can be seen that once fractures occur and new surfaces are generated, the desorption process speeds up due to the increase of new surfaces, which further enhances the diffusion process and greatly reduces the total concentration in the sample. Therefore in the relative intact areas, the concentration keeps higher, while in the fractured areas in the central part close to the major crack, the concentration is very low.

Figure 7 shows the fluid velocity at the fixed grids during the crushing process. In this test, the fluid flow outside solid particles is modelled using lattice Boltzmann method (LBM). Gas desorption and moving boundary conditions are implemented at the surface of particles which are exposed to fluid grids. The detailed coupling of DEM and LBM is not discussed here but can be found in our published papers [11, 12]. We can clearly see that the gas is desorbed from the sample surface in the beginning with high speed (left plot). With the generation of cracks and loss of gas, the desorption rate is reduced. From the middle and right plot in Figure 7, we can still observe higher fluid velocity at the bottom part of the major crack. This is because the desorbed gas, although not so strong now, is trying to escape through the narrow passage. To our knowledge, this is the first attempt to model this process.

6. Discussions and Conclusions

In this paper, particle scale diffusion was implemented in the open source discrete element code, the Esys-Particle. Diffusion of gas occurs between particles since each particle was regarded as a porous particle. We specifically focused on how to calibrate the particle scale diffusion coefficient. For the regular 2D packing, we derived theoretically the relation between micro- and macrodiffusion coefficients. This relation was verified in several numerical tests where the macroscopic diffusion coefficient can be determined numerically based on the half-time of a desorption scheme. To further test the coupled model, we simulated diffusion and desorption in the circular sample, and the numerical results match the analytical solution very well. We also simulated the process of diffusion and desorption while the sample is loaded and crushed into small pieces. Desorption occurs at the solid surface when the particle is exposed to the fluid. In this case the analytical solution is not available. Since diffusion process can be described by Fick’s diffusion law or diffusion equation at macroscopic scale, this process can easily be modelled using continuum approaches (i.e., finite difference, finite element method). However, once fractures occur, the continuum based methods are not suitable to model the discontinuity caused by large deformation of solids. The discrete element method does not suffer from this limitation. This is the major advantage of the approach presented in this paper.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.