An outburst of coal and gas is the rapid release of a large quantity of gas in conjunction with the ejection of coal and possibly associated rock, into the working face in underground coal mines. With an increase in depth of coal mining, outburst intensity and frequency tend to increase, posing more threats to lives of miners and facilities. Therefore in recent decades, the subject of outbursts has been a focus of interest in major coal-producing countries, particularly in China where some coal mines extract coal seams of greater than 1000 m depth. A number of hypothesis and theoretical models have been proposed to explain the outburst process, such as cavity theory (Briggs 1920), pocket theory (Farmer and Pooley 1967), dynamic theory (Shepherd et al. 1981) and spherical shell destabilization theory (Jiang 1998). However there is still no single theory or hypothesis that can explain the entire outburst process due to the wide variety of conditions under which outbursts occur. Field observations and laboratory studies reveal that the occurrence and development of an outburst is the result of combined effects of stress redistribution, coal gassiness and physical–mechanical properties of coal. It is generally recognized that for an outburst to occur coal must be deformed and failed under an effective stress and gas in coal must be able to desorb rapidly from the coal and eject the failed coal into a mining opening instantaneously (Lama and Bodziony 1998; Li 2001; Cao et al. 2003; Aguado and Nicieza 2007; Yuan et al. 2011; Torano et al. 2012).

Current knowledge about outbursts is drawn largely from field observations and laboratory studies (Yuan 2004, 2008). The rapid advance of computer technology has enabled the use of numerical simulations to gain useful insights on outbursts. A reliable numerical model permits certain factors to be varied such that their effect on outbursts can be studied. Along this strand of research, some attempts have been made to numerically model the process of an outburst, including mainly a phase transformation model (Litwiniszyn 1985), a gas desorption and flow model (Paterson 1986), a boundary element model (Barron and Kullmann 1990), an airway gas flow model (Otuonye and Sheng, 1994), a fracture mechanics model (Odintsev 1997), a finite element model (Xu et al. 2006), a plasticity model, and a coupled solid–fluid model (Xue et al. 2011). Despite these great efforts there is no single numerical model that can accurately simulate the entire process of an outburst because an outburst includes several interacting processes, including coal deformation and failure, coal fracture and fragmentation, gas desorption, mass transfer between adsorbed gas and free gas, flow of gas and water within coal cleats, gas dynamics and transport of failed and fragmented coal (Xue et al. 2014). While most of these models include deformation of the solid, desorption of gas, and Darcy flow of the fluid, however the fragmentation of the solid, which is a key process of an outburst, is either ignored or modeled with continuum damage mechanics which cannot naturally model the discrete nature of solid fragmentation and movement.

In this study, a new model is developed to simulate the entire process of an outburst. The new model couples two well-developed numerical approaches: the discrete element method (DEM) and the lattice Boltzmann method (LBM). The former explicitly models the deformation, fracture and fracture of the solid, while the later models fluid flow. This paper describes the basic principles and data interaction in the coupled DEM and LBM model. A simple example to simulate the process of an outburst is also presented to demonstrate the potential capability of the coupled model.

The DEM is based on the concept that the material to be modeled can be represented as a collection of discrete solid particles interacting with one another at their contacts (Cundall and Strack 1979). At each time step, the calculations performed in DEM alternate between integrating equations of motion for each particle, and applying the force–displacement law at each contact, through which the contact forces are updated based on the relative motions between two particles and their relevant contact stiffness. A number of DEM codes are available depending on the nature of interest and the details of simulation. In this study we use ESyS-Particle which permits particles to be bonded so that tensile forces can be transmitted. Fracturing is represented explicitly as broken bonds, which form and coalesce into macroscopic fracture of intact materials such as rocks.

ESyS-Particle is an open source DEM code developed by the Australian Computational Earth Systems Simulator (ACeESS) and designed for execution on parallel supercomputers (Mora and Place 1993; Abe et al. 2004; Wang 2009). The major features that distinguish ESyS-Particle from other DEM codes are the explicit representation of particle orientations using unit quaternion, complete interactions, and a new way of decomposing relative rotations between two rigid bodies (Wang et al. 2006). ESyS-Particle has been successfully used in the study of rock fracture and earthquake dynamics (Place et al. 2002; Abe et al. 2006; Wang and Mora 2008; Wang and Alonson-Marriquin 2009).

In ESyS-Particle, the solid particle motion is decomposed into translational motion of the centre of mass and rotation about the centre of mass. The translational motion is governed by the Newtonian equation (Eq. (1)). The particle rotation is governed by the Euler’s equation in the body-fixed frame (Eq. (2)). The orientation of each solid particle is explicitly described by the unit quaternion. A quaternion for each particle satisfies the following equation (Eq. (3)) (Evans 1977). The relationship between interactions and relative displacements between two bonded particles can be written in the linear form (Eq. (4)).

where \(\ddot{\varvec{r}}(t)\) and M are the position of a particle and its particle mass respectively. \(\varvec{f}(t)\) is the total forces acting on the particle, which may include the spring forces by neighboring particles, the forces by the walls, viscous force, and gravitational force.

where \(\tau_{x}^{b}\), \(\tau_{y}^{b}\), and \(\tau_{z}^{b}\) are the components of total torque \(\varvec{\tau}^{b}\) expressed in body-fixed frame, \(\omega_{x}^{b}\), \(\omega_{y}^{b}\) and \(\omega_{z}^{b}\) are components of angular velocities \(\omega^{b}\) measured in body-fixed frame, and \(I_{xx}\), \(I_{yy}\) and \(I_{zz}\) are the three principle moments of inertia in body-fixed frame in which the inertia tensor is diagonal.

where:

where \(\Delta \varvec{r}\), \(\Delta \varvec{s}_{{\mathbf{\rm 1}}}\), \(\Delta \varvec{s}_{{\mathbf{\rm 2}}}\) are the relative displacements in normal and tangent directions, \(\Delta\varvec{\alpha}_{\varvec{t}}\), \(\Delta\varvec{\alpha}_{{_{{\varvec{b}{\rm 1} }} }}\) and \(\Delta\varvec{\alpha}_{{_{{\varvec{b}{\rm 2} }} }}\) are the relative angular displacements caused by torsion and rolling, \(\varvec{f}_{\varvec{r}} ,\, \varvec{f}_{\varvec{{s}{{\rm 1}} }} ,\,\varvec{f}_{{\varvec{s}{{\mathbf{\rm 2}}} }} ,\varvec{\tau}_{\varvec{t}} ,\,\varvec{\tau}_{{\varvec{b}{\varvec{\rm 1}} }} ,\;{\text{and}}\;\varvec{\tau}_{{\varvec{b}{\varvec{\rm 2}} }}\) are forces and torques, and \(K_{r} ,K_{{{{s}}{\rm 1} }} ,K_{{{{s}}{\rm 2} }} ,K_{t} ,K_{{{{b}}{\rm 1} }} ,K_{{{{b}}{\rm 2} }}\) are relevant stiffness. Assuming that the bonds are identical in every direction, \(K_{s} = K_{{s{\rm 1} }} = K_{{s{\rm 2} }}\) and \(K_{{b}} = K_{{b{\rm 1} }} = K_{{b{\rm 2} }} .\)

The LBM is built on a mesoscopic scale in which the fluid is described by a group of discrete particles that propagate along a regular lattice and collide with each other. The LBM solves the particle velocity distribution function f. The completely discretized equation, with the time step and space step, is given by BGK model (Chen and Doolen 1998).

where \(\tau\) denotes the lattice relaxation time, \(\varvec{e}_{\alpha }\) is the discrete lattice velocity in direction \(\alpha\), \(\varvec{x}_{\varvec{i}}\) is a point in the discretized physical space, and \(f_{\alpha }^{eq}\) is the equilibrium distribution function. Equation (5) is usually solved in the following two steps: a collision step and a streaming step:

where \(\tilde{f}_{\alpha }\) represents the post-collision state. The equilibrium distribution function \(f_{\alpha }^{eq}\) in a square and nine-velocity lattice that is typically referred to as the D2Q9 model (Fig. 1) is of the form:

A 2-D 9-velocity lattice model (D2Q9)

where \(c = {{\Delta x} \mathord{\left/ {\vphantom {{\Delta x} {\Delta t}}} \right. \kern-0pt} {\Delta t}}\) is the lattice speed, \(\rho {\kern 1pt}\) is the lattice fluid density, \(\varvec{u}\) is the macroscopic velocity, \(w_{\alpha }\) is the weighting factor given by

Then the macroscopic quantities such as density, momentum and the fluid pressure can be obtained

where \(c_{s} = {c \mathord{\left/ {\vphantom {c {\sqrt 3 }}} \right. \kern-0pt} {\sqrt 3 }}\) is the speed of sound in this model.

A number of LBM codes are available. In this study we use OpenLB which is a C++ package for the implementation of Lattice Boltzmann simulations, addresses a vast range of problems in computational fluid dynamics, and is publicly available. The core of the implementation is a rectangular grid which can be used to build higher-level structures such as locally refined grids or complex geometries by connecting a set of them together. Efficient parallelization is achieved through the message passage interface (MPI) and OpenMP extensions. It supports advanced data structures that take into account complex geometries and parallel program executions. The programming concepts strongly rely on dynamic genericity via the use of object oriented interfaces as well as static genericity by means of templates. This design allows an efficient, straightforward and intuitive implementation of LBM.

A number of issues are considered in coupling of DEM (the Esys_Particle code) and LBM (OpenLB), including moving boundary conditions for a curved shape, momentum transfer between a solid particle and the fluid, and the force from the fluid to solid particles. A curved wall separates the solid region from the fluid region. The lattice node on the fluid side of the boundary is denoted as \(\varvec{x}_{\varvec{f}}\) and that on the solid side is denoted as \(\varvec{x}_{\varvec{b}}\) (Fig. 2). The particle momentum moving from \(\varvec{x}_{\varvec{f}}\) to \(\varvec{x}_{\varvec{b}}\) is \(\varvec{e}_{\alpha }\) and the reversed one from \(\varvec{x}_{\varvec{b}}\) to \(\varvec{x}_{\varvec{f}}\) is \(\varvec{e}_{{\tilde{\alpha }}} = - \varvec{e}_{\alpha }\). The intersection of the wall with the lattice link is denoted by \(\varvec{x}_{\varvec{w}}\). The particle surface can intersect the link between two nodes at arbitrary distance and the fraction of an intersected link in the fluid region is

The moving curved wall boundary condition

The reflected distribution function at node \(\varvec{x}_{\varvec{f}}\) can be calculated using an interpolation scheme (Yu et al. 2003)

where \(w_{\alpha }\) is the weight factor, \(\rho_{w}\) is the fluid density at node \(\varvec{x}_{\varvec{f}}\), and \(\varvec{u}_{\varvec{w}}\) is the velocity of the solid particle. The last term in Eq. (12) represents the momentum transferred from the solid particle to the fluid. The fluid force acted on the particle surface can be obtained using

where the first summation is taken over all fluid nodes at \(\varvec{x}_{\varvec{b}}\) adjacent to the particle boundary and the second is taken over all possible lattice directions pointing towards a particle cell. This force is added to the particle force in DEM code.

The coupling procedure is as follows: at each time step, the equations of motion for DEM particles are first solved by obtaining the particle positions and velocities; this enables the state (or flag) of each LBM cell to be determined (i.e., “Solid”, “Fluid”, “Solid Boundary” or “Darcy”); the LBM calculation is then carried out to yield the fluid velocity and pressure fields, from which the drag forces acting on particles can then be computed (Fig. 3); the forces acting on particles and each fluid cell are used to update the positions and velocities of particles at the next time step. This procedure is repeated until a specified time step is reached.

The coupling procedure of DEM-LBM

Darcy flow is modeled by LBM code. To implement the LBM for Darcy flow, consider the translational collision step as a second intermediate step after streaming, denoted by \(\it \it {\varvec f^{*}}\)

Then the porous medium step has the form

where \(n_{s}\) is a damping parameter. Note that for \(n_{s} = 0\) we recover the normal free-fluid node and for \(n_{s} = 1\) we have a bounce-back-like condition that effectively makes the medium impermeable. For values of \(n_{s}\) between 0 and 1 we have partial bounce-back. Dardis and McCloskey (1998a, b) indicate the permeability k of a medium with damping \(n_{s}\) could be computed as

where \(\nu\) is the kinematic viscosity of the fluid.

Diffusion is managed by DEM part in the coupled scheme. In this case, the solid particles represent porous coal. It is assumed that the voids inside particles are much smaller than the size of particles and the porosity is just an average concept for each particle. There is an average and uniform pore pressure \(p_{i}\) and concentration \(c_{i}\) for each particle i. The fluid exchange between two contacted particles i and j is determined by the Fick’s first law of diffusion:

where D is the diffusion coefficient of the link.

The hydro-mechanical coupling is implemented based on Biot’s linear pore-elastic theory. According to this theory, the constitutive equations of a porous medium can be written as (Detournay and Cheng 1993).

where \(p\) is pore pressure, \(P = - {{\sigma_{kk} } \mathord{\left/ {\vphantom {{\sigma_{kk} } 3}} \right. \kern-0pt} 3}\) is the mean or total mechanical pressure (isotropic compressive stress), \(\varepsilon = \varepsilon_{kk} = {{\Delta V} \mathord{\left/ {\vphantom {{\Delta V} V}} \right. \kern-0pt} V}\) is the volumetric strain (positive for extension), \(\varsigma = {{V_{f} } \mathord{\left/ {\vphantom {{V_{f} } V}} \right. \kern-0pt} V}\) is the variation of fluid content (positive corresponds to a “gain” fluid), \(\alpha\) is Biot coefficient, \(B\) is the Skempton pore pressure coefficient and \(K_{m}\) is the drained bulk modulus of the material \(V\) and \(V_{f}\) are the volume of the material and fluid respectively. From Eq. (19), the following equation is obtained

The pore pressure for particle i is updated according to

where the summation j goes through all the neighboring particles of particle I, and

Gas desorption is implemented in both DEM and LBM codes. The reduction of concentration in a solid particle is described by

where c is average matrix gas concentration; p is gas pressure (Eq. (18)); and \(\tau_{d}\) is sorption time. \(c(p) = {{V_{c} } \mathord{\left/ {\vphantom {{V_{c} } V}} \right. \kern-0pt} V}\), where V is the volume of the particle and \(V_{c}\) is the volume of adsorbed gas in the coal matrix at pressure \(p\) that is described by the Langmuir adsorption isotherms:

where V _L and P _L are Langmuir volume and Langmuir pressure.

Then in the LBM code, the loss of concentration in a solid particle calculated by Eq. (23) is distributed to its fluid neighboring nodes. This is done by first distributing it to the boundary nodes (which have direct link with the fluid grid) in the solid particle, and then steaming to the fluid nodes.

To verify the coupled DEM-LBM model for the outburst simulation, two of the key processes in outbursting, namely particle motion in the fluid and gas desorption of a particle aggregate, are simulated and the results are presented here. Figure 4 shows particle motion in the fluid. Three particles bonded as a rigid body are driven by the fluid flow (Poiseuille flow). The motion and rotation of the rigid body, as well as vortex, are clearly observed. In this simulation Reynolds Number R _n  = 100, and particle size is five times of the grid size of the fluid. This example verifies the algorithm of the two-way coupling of solid and fluid: moving boundary condition and dragging forces.

Three bonded particle driving by fluid flow

The simulation of gas desorption of a particle aggregate is shown in Fig. 5. Seven bonded particles having the same initial gas content are exposed to the fluid. In this simulation gas desorbs from particle interfaces, flows into the fluid grids, and continues to move. Larger velocity (red color) can be seen close to the boundary where desorption occurs.

Simulation of desorption of gas from a bulk of coal consisting of several particles

Finally, a small outburst simulation is carried out. The model consists of 1003 equal-sized particles and a fluid grid. Constant confining pressure is applied at left, right, top and bottom boundaries. The soft and fragile middle part (shown in dark and light grey in Fig. 6) sandwiched by harder and stronger rock models a coal seam. Excavation starts from left side with constant speed by removing the dark grey part. Two supporting walls (horizontally dashed lines) are placed after the excavation, growing, and following the moving excavation wall (vertical dashed line) at the right end (Fig. 6).

Schematic diagram of excavation and outburst simulation

Figure 7 shows the results of particle motion and fluid velocity at several time steps. In the early stage of excavation, there are no fracture events, and desorption of gas from the exposed surface is clearly observed (Fig. 7a). At the time step 7000 (Fig. 7b) fractures occur at the excavation surface, and the excavation is stopped. After this the whole system of particle and fluid evolves by itself (Fig. 7c, d), and eventually fractured coal is propelled by the fluid and ejected from the face. A cave with a small mouth and big inside void is formed, which is widely observed after the occurrence of an outburst in coal mines.

Particle motion and fluid velocity during the outburst development

Figure 8 shows the concentration of gas adsorbed in the solid coal at the same time steps as Fig. 7. Lower concentration at the surface is observed due to desorption (Fig. 8a, b). The concentration changes in particles are caused by diffusion. The particles from the burst coal have the lowest gas concentration (Fig. 8c, d) because they are totally exposed and have larger surface areas, therefore continue gas desorption while moving away from the face. It is observed that lower initial gas concentration will delay outbursts when all other parameters remain the same. This clearly suggests the important role of initial gas content played in outbursts.

Gas concentration at different time steps

A coupled DEM-LBM model is developed to simulate the entire process of an outburst by coupling two well developed numerical approaches: DEM for solid and LBM for fluid. ESyS-Particle is used as the DEM code because it permits particles to be bonded so that tensile forces can be transmitted, and fracturing is represented explicitly as broken bonds, which form and coalesce into macroscopic fracture of intact materials such as coal. OpenLB is used as the LBM code as it supports advanced data structures that take into account complex geometries and parallel program executions. Both of ESyS-Particle and OpenLB are open source codes. The coupled model includes the most important factors of outbursts, including deformation, fracture and fragmentation of solids, free flow of fluid, Darcy flow, diffusion, desorption of gas, and the coupling of these factors.

A simple 2-dimensional model of coal excavation in an underground coal mine is studied with the coupled simulator. The preliminary results with small scale simulations with the coupled model are encouraging as the entire process of an outburst is well reproduced. It should be noted that the current version of the coupled DEM-LBM model requires fine-tuning such as calibration against field data to simulate a real-case outburst scenario. This is currently being undertaken and will be reported later. Nevertheless, the new model has the potential to numerically investigate the full mechanism and interaction of contributing factors of outbursts.