Rock Deformation In Parallel Coal Seams

Paper 2A 21 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

A NUMERICAL APPROACH FOR COUPLED GAS LEAK FLOW AND COAL/ROCK DEFORMATION IN PARALLEL COAL SEAMS P.D. Sun Department of Environmental Engineering, Hangzhou University of Commerce, Hangzhou, 310035, China. [email protected]

Abstract: It is a well-established fact that the description of the interaction between coal/rock and gas in the parallel coal seams is difficult. Nonetheless, the problem of gas leak flow in multi-coal-seam must be studied for the sakes of coal mining safety and coal gas draining.From the viewpoint of interaction mechanics for solid and gas, a coupled mathematical model is presented for solid coal/rock deformation and gas leak flow in parallel deformable coal seams. Numerical solutions using the SIP (Strong Implicit Procedure) method to the coupled mathematical model for double parallel coal seams are also developed in detail. Numerical simulations for the prediction of the safety range using protection layer mining are performed with experimental data from a mine with potential danger of coal/gas outbursts. Analyses show that the numerical simulation results are consistent with the measured data on the spot. The coupled model shows a positive future for applications in a wide range of gas-leak-flow-related problems in mining engineering, gas drainage engineering and mining safety engineering. Keywords: coupled models, rock deformation, gas leak flow, strong implicit procedure, numerical simulation, parallel coal seams

1. INTRODUCTION Underground coal seams are mostly distributed in adjacent multi-layers in most countries of the world, especially in China. Many problems in mining engineering and rock engineering remain unsolved such as prediction and control of methane gas emission from the coalface during mining, safety mining area prediction issues related to coal and gas outbursts in protective layer mining, prediction of methane gas drainage rate near to the mining seam and from original coal seam, as well as prediction of gas seepage rate in multi-coalseam. These problems, in essence, are linked with rock mass deformation and gas leak flow by Sun(1998, 1999, 2000, 2002).. Under the action of a pore gas pressure gradient in the multi-coal-seam system, gas flow emission occurs near to the mining seam, which in turn leaks through the layer of intercalation with lower permeability, and flows into the mining coal faces or gas drainage bores. This is referred to as gas leak flow by Sun(1998). As gas flows during coal mining and gas drainage, the effective stress of the solid skeleton of coal/rock mass changes, and the pore pressure in the coal seam or gas drainage bores also changes, which lead to a gas flow through a network made of micro-pores and micro-

cracks. Moreover, the gas flowing through coal/rock cracks also alters with the variance of the gas permeability in Zhao(1994) and Sun(1990,1998). Therefore, the interaction between the coal/rock mass and gas must be taken seriously into account for simulating the action of methane gas. It is a well-established fact that the description of the interaction between coal/rock mass and gas in the gas leak flow system is difficult. Nonetheless, the problem of gas leak flow must be studied for the sake of coal mining safety and gas drainage in the multi-coal-seam. In this study the numerical approach for coupled coal/rock mass deformation and gas seepage in double coal seams are studied.

2. SOLID-GAS COUPLED MODELS 2.1 Equations for gas leak flow According to the continuity equations for the gas leak flow by Sun(1998) , we have Eq.1, which is expressed in the tensor notation in the Cartesian coordinate system. T ∂P 2 ∂em (1) ( T P 2 , ) ± 3 ( P2 − P 2 ) = Sun( P , t) m − 2 P mi m i

,i

bm* b3*

1

2

m

∂t

m

∂t

(m=1, 2; i, j=1, 2, 3)

1

Paper 2A 21 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

nm a m b m (1 − e − λ m t ) + Pm P m (1 + b m P m ) 2 *

Sun ( P m , t ) =

λ*m =

12 dm

D n Pn πb m Pm

(2) (3)

Where, Sun( P, t ) is named as Sun function which is a state one describing the methane gas flow. In general, the gas flow grows fast as the value of the Sun function increases, and vice versa. The full deduction detail of the Sun function is given by Sun (1998). Where, the subscript m denotes the coal seam or intercalation level with m=1 representing the bottom coal seam, m=2 standing for the top one and m=3 symbolizing the intercalation level; Tm is the gas permeability; The Pm is the pore gas pressure; b*m is the average thickness; n m is the porosity; a m is the maximum gas adsorption coefficient; bm is the gas adsorption coefficient; Dm is the gas diffusion coefficient; e m is the deformation rate; d m is the average coal micro-grain size; the Pn equals one atmosphere; and the Dn is the gas diffusion coefficient when the pore pressure equals Pn .

The equations for the coal/rock mass deformations are obtained by Sun(1998) and using the equilibrium Eq.4 for steady-state-equation of motion with the overall stress in Bear(1972), Jaeger and Cook(1979).

(4)

( m = 1, 2 )

Taking into account the action of the pore gas pressure, we have the Eq.5 by Sun(1998) and Zhao(1994). Umj, ji (λm + Gm ) + Umi, jjGm + Fmi + (αm Pm ),i = 0

(5)

Where, the U mx , xx is the displacement function; λm is the Lame-constant, δ ij is the Kronecker function, Gm is the shear modulus and the term for the stress.

∂e1   bb ∂t ∂t  P2 − P2 ∂P 2 ∂e  (T2i P22,i ), i −T3 1 * * 2 = Sun( P2 , t) 2 − 2 P2 2  b2 b3 ∂t ∂t   U mj, ji (λm + Gm ) + U mi, jj Gm + Fmi + (αm Pm ), i = 0   em = U mi,i ,   Pm ( x, y, z,0) = Pmo( x, y, z );  Pm ( xo , y o , z o , t) t∈(0,+∞ ) = Pa ;  r  r Tmi∂Pm / ∂n t∈(0,+∞ ) = Qm ( x, y, z , t);   em ( x, y, z, t) t =0 = emo ;   U mi ( x, y, z, t ) t= 0 = 0;   U mi ( x, y, z, t ) t∈(0, +∞) = U mi .  (T1i P12,i ), i +T3

P12 − P22 * * 1 3

= Sun(P1 , t )

∂P12

− 2 P1

(6)

2.2 Equations for coal/rock deformation

σ ijm, j + F mi = 0

described with the boundary and initial conditions (Eq.6) by Sun(1998). Where, m=1,2; j=1,2,3; Pm0 is the initial pore gas pressure; Pa is the pressure in free space; e m0 is the initial deformation rate; U mi is the boundary value for displacement; and Qm is the quality vector for gas seepage rate from boundary coal-mass surface.

Fmi

is the free

2.3 Solid-gas coupled model for the gas leak flow system With the continuity equations for gas leak flow and the equilibrium equations for coal/rock mass defor- mation provided, the coupled mathematical model for the interaction between coal/rock mass and gas in the gas leak flow system can be

3. NUMERICAL SOLUTION OF THE COUPLED MATHEMATICAL MODEL For the task of solving the solid-gas coupled mathematical model, several numerical approaches, such as finite differential method (FDM), finite element method (FEM), and discrete element method (DEM), can be adopted. But the most effective method is clearly dependent on the practical conditions. Regarding the approximately regular geometric shapes of the coal seams as well as the non-linearity of the equations, we adopt the Strong Implicit Procedure (SIP) method by Stone(1986) for numerically solving our coupled model. The SIP method is one of the FDM methods and is characterized by: 1. the convergence speed for iteration of the method is weakly dependent on the coefficient matrix; 2. it is insensitive to the selection of the iteration parameters; 3. it is suitable for solving both linear equations and non-linear ones; and 4. it is often more effective than other approaches regarding the numerical solution for extremely complex equations.

2

Paper 2A 21 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

Let Pm2 = P~m ; Qm* = 0; ( m = 1,2) in Eq.6, then the equations for the solid-gas coupled model can be simplified into: Let Pm2 = P~m ; Qm* = 0; ( m = 1,2) in Eq.6, then the equations for the solid-gas coupled model can be simplified into: ~ ~ ~ ~ ∂ ∂P1 ∂ ∂P1 ∂ ∂P1 ~ ∂P ~ ∂e T ~ ~ [T1x ] + [T1y ] + [T1z ] = Sun(P1,t) 1 −2 P1 − *3 * ( P1 − P2 ) ∂x ∂x ∂y ∂y ∂z ∂z ∂t ∂t b1 b3

~ ~ ~ ~ ∂ ∂P2 ∂ ∂P2 ∂ ∂P2 ~ ∂P ~ ∂e T ~ ~ [T2x ] + [T2y ]+ [T2z ] =Sun(P2,t) 2 −2 P2 + *3* (P1 −P2) ∂x ∂x ∂y ∂y ∂z ∂z ∂t ∂t b2b3 G(

G(

G(

∂2u ∂x

2

∂2v ∂x

2

∂ 2w ∂x 2

+

+

+

∂ 2u ∂y

2

∂ 2v ∂y

2

∂2w ∂y 2

e=

∂ 2u

+

+

) + (λ + G )

∂z 2 ∂ 2v

+

∂z 2

) + (λ + G )

∂2 w ∂z 2

∂e

( α1 P1 ) = 0

(10)

∂e ∂ + f z + (α1P1 ) = 0 ∂z ∂z

(11)

∂y

+ fx +

+ fy +

∂x ∂ ∂y

∂u ∂v ∂w + + ∂x ∂y ∂z

N A B (1 − exp(− λm * t )) ~ Sun( Pm , t) = ~m + m m~ ~ Pm Pm (1 + Bm Pm )

Bi , j = Ti , j ,k Ti , j −1 ,k /(∆Y ) 2

(16)

Di , j = Ti , j , k Ti −1, j ,k /(∆X ) 2

(17)

Fi , j = Ti , j , k Ti +1 , j ,,k /(∆X ) 2

(18)

H i , j = Ti , j , k Ti , j +1,k /(∆Y ) 2

(19)

(8)

(9)

∂e

) + (λ + G )

∂

(7)

(α1 P1 ) = 0

∂x

Where, the i and j are the indices for the grid points with M and N the size of the grid, the k is the index for the discrete points in time domain, and the coefficients Bi , j , Di , j , Fi , j , Hi , j , Ei , j and Ai , j , k are represented respectively from Eq.16 to Eq.21.

~ E i, j = −( Bi , j + Di , j + Fi , j + H i , j + Sun[(P1 )i , j, k∆t ] / ∆t) (12) (13)

Eqs. 7 & 8 are respectively the equations for describing the gas leak flow of each of the double parallel coal seams, Eqs. 9~12 are the equations for describing the coal/rock deformations with the displacements, and Eq.13 is Sun function which is a state one describing the methane gas flow. From Eqs.7~13, we can clearly see that the six unknown variables P1 , P2 , e , u , v and w can be uniquely determined by the six independent equations.

(20) ~ ~ Ai, j,k = −Sun[(P1 )i , j,k∆t ](P1 )i, j ,k−1 / ∆t ~ ~ ~ − 2 (P1)i , j,k −1 (ei, j ,k −ei, j,k−1) / ∆t − T3[(P1 )i, j,k−1 − (P2 )i , j ,k−1 ]/ b1*b*3

(21) i-1

j+1

j

i

j-1

Ei,j

Fi,j

¦ x ¤i

¦ ¤ xi+1

(a) coefficients

3.1 Discretization of the differential equations for gas leak flow For simplicity of description, we investigate Eq. 7 in two-dimensional case: ~ ~ ~ ∂P1 ∂P ∂ ∂ ~ ∂P ~ ∂e T 3 ~ ~ [T ] + [T1 y 1 ] = Sun(P1, t) 1 − 2 P1 − ( P − P2 ) ∂x 1x ∂x ∂y ∂y ∂t ∂t b1*b*3 1

(14) As shown in Figure 1, the gas leak flow equation (Eq.14) can be discretized into Eq. 15 with the FDM method. Note that the geometric mean is used to approximately score the distribution of gas permeability within the heterogeneous gas field. ~ ~ ~ ~ ~ Bi, j (P1)i, j−1,k +Di, j (P1)i−1, j,k +Ei, j(P1)i, j,k +Fi, j(P1)i+1, j,k +Hi, j(P1)i, j+1,k = Ai, j,k (i =1,2,...,M; j =1,2,...N; k=1,2,...,NT)

(15)

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

j=3

.

.

.

.

.

.

.

.

.

.

j=2

.

.

.

.

.

.

.

.

.

.

j=1

.

.

.

.

.

.

.

.

.

i=

1

2

3

¦ ¤ yj

¦ y ¤j-1

Bi,j ¦ ¤ xi- 1

.

. ¦ ¤ yj+1

Hi,j Di,j

.

j=N

i+1

…

. M

(b) node arrangements

Figure 1. Schematic view of discretization with FDM By applying Eq.15 to each grid point, together with the boundary conditions, we can obtain a matrix equation:

[A]

[P~]

=

[F ] (22)

Where [A ] is an L × L (with L = N × M ) square ~ matrix named coefficient matrix, [P ] is a column vector with L unknown variables of pore pressure, and [F ] is a column vector with L free items. To describe Eq.22 more concretely, without losing generality, let k denotes a time step in the time

3

Paper 2A 21 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

domain, and N = M = 3 , then when j = 1 , i will be 1, 2, or 3. Similarly, when j is 2 or 3, i will also be 1, 2, or 3. Therefore totally 9 grid points will be acquired from discretizing the space and 9 equations will be derived from Eq.15.  E1,1 F1,1 0 H1,1 0 D  2,1 E2,1 F2,1 0 H2,1  0 D3,1 E3,1 0 0  0 E1,2 F1,2  B1,2 0  0 B2,2 0 D2,2 E2,2  0 B3,2 0 D3,2  0  0 0 0 B1,3 0  0 0 0 B2,3  0  0 0 0 0 0 

~ ~ 0 0 0 0   P1,1   f1.1   ~   ~  0 0 0 0  P2,1   f 2,1  ~ ~ H3,1 0 0 0   P3,1   f3,1   ~   ~  0 H1,2 0 0  P1,2   f1,2  ~ ~ F2,2 0 H2,2 0  × P2,2  =  f2,2   ~   ~  E3,2 0 0 H3,2  P3,2   f 3,2  ~  ~  0 E1,3 F1,3 0   P1,3   f1,3   ~ 0 D2,3 E2,3 F2,3  P2,3   ~ f   ~   ~2,3   B3,3 0 D3,3 E3,3  P3,3   f 3,3     

(23) The elements in the column vector [F ] at the right side of Eq.23 are free items. They can be calculated from Eq.15 with boundary conditions: ~ ~ ~ f1,1 = A1,1,k − B1,1 ( P1 )1,0,k − D1,1 ( P1 ) 0,1,k (24) ~ ~ f 2,1 = A2,1,k − B2,1 ( P1 ) 2,0,k

(25)

~ ~ f 3,1 = A3,1,k − B3,1 ( P1 )3,0, k

(26)

~ ~ f1,2 = A1,2,k − D1,2 ( P1 ) 0,2,k

(27)

~ f 2,2 = A2,2,k

(28)

~ ~ f 3,2 = A3,2,k − F3,2 ( P1 ) 4,2,k ~ ~ ~ f1,3 = A1,3,k − D1,3 ( P1 ) 0,3,k − H1,3 ( P1 )1,4, k

(29)

solve the linear equation system. The LU decomposition is performed on the coefficient matrix [A ] . Thus we have

[A]

=

[A]

[P~]

= ([L] [U ])

[]

=

[L] ([U ]

[P~])

[]

~ ([A ] + [M ]) P

=

[F ] (34)

[]

~ ~ = ( [A ] + [M ]) P − ([A ] P − [F ])

(37)

Its iterative representation can be expressed as:

[]

[]

[]

~ n +1 ~n ~n ([A ] + [M ]) P = ([A ] + [M ]) P − ([A ] P − [F ])

(38)

Where, n is the iterative number. To reduce the rounding error, let

~ ~ f 2,3 = A2,3,k − H 2,3 ( P1 ) 2,4,k

(31)

[γ ]

~ ~ ~ f 3,3 = A3,3,k − F3,3 ( P1 ) 4,3,k − H 3,3 ( P1 ) 3,4,k

(32)

~

[P~]

By first solving for the column vector [V ] such that [L] [V ] = [F ] (35) And then solving [U ] [P~] = [V ] (36) ~ We can get the required [P ] .To facilitate the LU decomposition, the coefficient matrix [A ] is approximated with ([A] + [M ]) .The process of solving the linear equation system (Eq.22) is described briefly as follows. The solving of the Eq.22 is equivalent to the solving of the following equation:

[δ ]n+ 1

Eq.16~20. Hence the 9 unknown variables of Pi, j (hence the pore pressure Pi, j (i = 1,2,3; j = 1,2,3) ) can be acquired by solving the linear equation system (Eq.23). For the more general case of Eq.22, the similar result will also be obtained. Let’s take a closer look at Eq.22. From Eq.23, we can find that the coefficient matrix [A ] is a 5diagonal matrix with non-zero elements near the main diagonal. Therefore we adopt LU decomposition rather than Gaussian elimination to

(33)

Where [L] is lower triangular with elements only on the diagonal and below, and [U ] is upper triangular with elements only on the diagonal and above. We can use the decomposition to solve the linear equation system (Eq.22):

(30)

In addition, the 9×9 elements of the coefficient matrix [A ] for the discrete matrix equation (Eq.23) are either constants or can be calculated with

[L] [U ]

=

[P~]n +1 − [P~]n [F ] − [A] [P~]n

(39)

(40) With Eqs.39~40, Eq.38 can be expressed as: ([A ]+ [M ]) [δ ]n + 1 = [γ ]n (41) After the LU decomposition, we have, [L ] [U ] [δ ]n +1 = [L ] [V ] = [γ ]n (42) Thus, we have Vi n, +j 1 = γ in, j − bi , jVi n, j+−11 − d i , jVi n−+1,1j / ei , j (43) n

=

With the variables

Vi n, +j 1

calculated, we can get

from Eq.44 the values of δ in, +j 1 .

[U ] [δ ]n +1 δ in, +j 1

= [V ]

= Vi n, +j 1 − f i , j δ in++11, j

(44) − hi , j δ in, +j +1 1

(45)

where bi , j = Bi , j − ωc i , j d i , j = Di , j − ωg i , j e i , j = E i , j + ω (ci , j + gi , j ) − bi , j hi , j − 1 − d i , j f i −1 , j

(46) (47) (48)

4

Paper 2A 21 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

(49) hi , j = ( H i , j − ωgi , j ) / ei , j (50) c i , j = f i , j Bi , j /(1 + ωf i , j −1 ) (51) g i , j = hi , j Di , j /(1 + ωhi −1 , j ) (52) Where, the ω in Eqs.46~52 is an iterative parameter given by the user. The iterative procedure is initiated by first calculating the [γ ]n with Eq. 40, then calculating the Vi n, +j 1 with Eq. 43. f i , j = Fi , j − ωci , j / ei , j

Following this, the

δ in, +j 1

is calculated with Eq.45

and substituted into Eq. 39 to calculate the [P ] . The procedure iterates till the following in equation is satisfied: ~ n +1

~ ~ Max Pi ,nj+ 1 − Pi n, j < ε

of the mechanical and physical properties of the two coal seams are given respectively in Reference by Sun(1998,2002). The numerical solution for the testing case is performed with the SIP method. In Figure 2a~2c, 3D contour for distributions of the pore gas pressure in the coal seam 8 when the upper-protect seam (the coal seam 7) is being exploited along the direction of the run of the coal seam at respectively the 100th , the 140th and the 194th day. In the figures, the horizontal axis shows the length along the direction of the coalface and the value is 2.4 * Number in meter, and the vertical axis shows the length of the coal seam along the direction of the run of the coal seam and the value is 4.0 * Number in meter. The contour line is distribution of pore gas pressure with unit MPa .

(53) Where ε is a small enough positive real number given by the user. In a similar way, another equation for gas leak flow (Eq.8) and deformation equations can also be solved with the SIP approach.

4. VISUAL SIMULATIONS OF THE COUPLED MODEL

Figure 2a. 100 th day

4.1 Implementations The SIP methods by Stone (1986) for the gas leak flow equations and the rock mass deformation equations for double parallel coal seams have been developed by Sun(1998) with Microsoft Visual C++ 6.0 under Windows2000 on a Pentium IV PC. The program is suitable for isotropic heterogeneous coal seams with irregular shapes as well as anisotropic heterogeneous coal seams. The systems with the first or the second boundary conditions as well as the hybrid boundary conditions can all be solved. The overview of the numerical implementations of the SIP methods for the solidgas coupled mathematical model for double parallel coal seams is described as follows.

Figure 2b. 140 th day

4.2 Testing case The proposed solid-gas coupled mathematical model and its numerical implementations have been tested with experimental data measured from the coalface N1709 in the Datong mine 2, Songzao mining district, China. The coal seams 7 and 8 are mostly horizontal and the average distance between the two seams is 7.2 meters. The coalface N1709 in coal seam 7 has a length of 120 meters. The details

Figure 2c. 194 th day Figure 2. 3D contour line for distributions of pore pressure in the coal seam 8.

5

Paper 2A 21 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

The 3-dimensional contour line in Figure 2 clearly shows that the pore pressure distributions in the coal seam 8 are time-dependent with the coal seam 7 (as the protective-seam) exploited firstly. If the minimum value of pore gas pressure were chosen as 0.5 MPa, the zone with the pore pressure less than 0.5 MPa in the seam 8 in Figure 2 would be the safety protective one. That is to say, it would be safe to work in this zone with the pore gas pressure less than 0.5 MPa in the coal seam 8, and coal/gas outbursts would not appear. The above deduction has been proved to be right by practice on the spot in Sun(1998). It can also be seen from Figure 3 that the pore pressure of methane gas drops fast in the direction of gob within 100 meters from the coalface in the coal seam 8. The pore pressures in the coal seam 8 drops dramatically from 0.9Mpa to 0.4 MPa in the direction of gob within 50 meters to 100 meters away from the coalface.

leak flow and coal/rock mass elastic -deformation is also described by the coupled model. The interaction between the elastic deformation of coal/rock mass and gas leak flow has been thoroughly studied with the built equations. Numerical simulations have been performed with the SIP approach with robustness and efficiency. Furthermore, the simulation results have shown a promising consistency with the measured data on the spot.

REFERENCES Bear, J. 1972. Dynamics of Fluids in Porous Media . American Elsevier Publishing Co.Inc. New York. Jaeger, J. C. and Cook, N. G. W. 1979. Fundamentals of Rock Mechanics. Third Edition, Chapman and Hall, London. Stone, H.L. 1986. Iteration solution of implicit approximations of multi-dimensional partial differential equations. J. Soc. Ind. Appl. Math. 5(6): pp.530-558. Sun, P.D. 1990. A new method for calculating the gas permeability of a coal seam. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 27(4): pp.325-327.

Figure 3. Curves for pressure down angle along the run of the coal seam 8.

5. CONCLUSIONS In this paper, a mathematical model for coupled coal/rock mass deformations and gas leak flow in parallel coal seams is presented from a new viewpoint of solid-gas interaction. It consists of equations for the gas leak flow, equations for coal/rock mass deformation and the boundary/initial conditions. The model is characterized by: (1) The main consolidation of gas leak flow through the coal/rock mass due to gas drainage or coal mining is described by the linearelastic deformation in the coal/rock mass; (2) The mechanism for gas leak flow with both seepage and diffusion is given by the equations of gas leak flow, and the gas desorption due to the pore pressure decreasing is also given by equations for the gas leak flow; and (3) The interaction between the gas

Sun, P.D. 1998. A study on the interaction mechanics for coal seam deformation and gas leak flow and its numerical simulations. University of Chongqing, Ph.D. thesis, (in Chinese with English abstract). Sun, P.D. and Xian, X.F. 1999. A study on coupled models for coal seam deformation, gas leak flow and its applications. Journal of China Coal Society. 24(1): pp.60-64 (in Chinese with English abstract). Sun, P.D. 2000. Interaction modeling for mining safety range and numerical simulations. Journal of Coal Science & Engineering. 6(2): pp.41-46. Sun, P.D. 2002. Sun model and its applications. Zhejiang University Press, Hangzhou (in Chinese with English summary). Sun, P.D. 2002. SIP analysis on coupled models for coal seam deformation and gas leakage flow. Journal of China Coal Society. 27(5): pp.276280 (in Chinese with English abstract). Zhao, Y.S.1994. Rock fluid mechanics in mine. China Coal Industry Press, Beijing (in Chinese).

6