Mathematical modeling of seismic explosion waves impact on rock mass with a working Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

^AleksandrP. Gospodarikov, YaroslavN. Vykhodtsev, Mikhail A. Zatsepin

Mathematical Modelling of Seismic Explosion Waves Impact on Rock Mass ...

Mining

UDC 622.235.535.2

MATHEMATICAL MODELING OF SEISMIC EXPLOSION WAVES IMPACT ON ROCK

MASS WITH A WORKING

Aleksandr P. GOSPODARIKOV, Yaroslav N. VYKHODTSEV, Mikhail A. ZATSEPIN

Saint-Petersburg Mining University, Saint-Petersburg, Russia

In the article, within the framework of the dynamic theory of elasticity, a mathematical model of the impact of seismic blast waves on rock mass is presented, including a working. The increase in the volume of mining operations in complex mining and geological conditions, taking into account the influence of the explosion energy, is closely connected with the analysis of the main parameters of the stress-strain state of the rock massif including a working. The latter leads to the need to determine the safe parameters of drilling and blasting operations that ensure the operational state of mining. The main danger in detonation of an explosive charge near an active working is a seismic explosive wave which characteristics are determined by the properties of soil and parameters of drilling and blasting operations. The determination of stress fields and displacement velocities in rock mass requires the use of a modern mathematical apparatus for its solution. For numerical solution of the given boundary value problem by the method of finite differences, an original calculation-difference scheme is constructed. The application of the splitting method for solving a two-dimensional boundary value problem is reduced to the solution of spatially one-dimensional differential equations. For the obtained numerical algorithm, an effective computational software has been developed. Numerical solutions of the model problem are given for the case when the shape of the working has a form of an ellipse.

Key words: mine working, mathematical model, seismic explosion wave, difference scheme, numerical algorithm

How to cite this article: Gospodarikov A.P., Vykhodtsev Y.N., Zatsepin M.A. Mathematical Modeling of Seismic Explosion Waves Impact on Rock Mass with a Working. Zapiski Gornogo instituta. 2017. Vol. 226. P. 405-411. DOI: 10.25515/PMI.2017.4.405

t = 0

Introduction. Stresses that arise in the working area during the seismic explosion waves impact on it depend on many factors, but, first of all, on the explosive agent (EA) charge power, the distance to the working, the detonation velocity, the mixture of explosive charge, etc. [5, 6, 10]. To determine the impact of the impact of a seismic explosion wave on a working, it is necessary to conduct a large number of full-scale tests, which is not always possible from the economic and technical point of view. Therefore, in order to evaluate the explosive effect on rock mass, the article uses a numerical simulation of the interaction of a longitudinal wave propagating in an elastic medium with a working [12].

Statement of the problem. To create a mathematical model of seismic explosion waves impact on a mine working the paper uses the equations of dynamic theory of elasticity of Mises [3, 9] written in the form of curvilinear coordinates (Fig. 1).

The figure shows two coordinate systems with the following notations being used:

1) O1Zn - rectangular coordinate system; the origin of coordinates O1 - mass center of a mine working; axis O^ is parallel to displacement velocity of undisturbed wave front; Fig. 1. Introduced coordinate systems

O

^AleksandrP. Gospodarikov, YaroslavN. Vykhodtsev, Mikhail A. Zatsepin

Mathematical Modelling of Seismic Explosion Waves Impact on Rock Mass ...

2) M1xy - curvilinear coordinate system; x - distance (MM1) from point M to a neat line of mining working; y - a length of G curve, calculated from point O to point M1; О - meeting point of front C1 of an incident blast wave with a boundary surface at the initial instant (t = 0). Taking into account the introduced coordinate system, we have a relation of the following form

r = R+xn,

where r - radius vector of point M; R - radius vector of point M1; n - unit normal vector n, x -unit vector of tangent line x in point M1. Then the following is correct

dr = dR + xdn + dxn = dyx + dxn + k(y) xdyx = (1 + k(y) x)dyx + dxn ,

where k(y) - line G curvature in point M1. Finally, we have

(dr )2 = dx2 + (1 + k ( y ) x)

2 7 2

= dx2 + H

22

where H = 1 + k ( y ) x - Lame coefficient.

Let us introduce the following notations: v1, v2 - components of velocity vector of rock mass particles in coordinate axes M1x and M1y accordingly; and through o11, a12, a22 - components of stress. Then the equations motion and Hooke's law differentiated by time, written in dimensionless form will have the following form [1, 7, 8]:

ÔV1 ÔCn 1 ÔC12 1 ÔH .

—1 = —11 +--— +--(c„ -C22);

Ôt Ôx H ôy H Ôx

ÔV2 _ ÔC12

+

1 ÔC

22

2 ÔH

+--

Ôt Ôx H Ôy H Ôx

C

12

ÔC

ÔV1 (1 - 2b) ÔV2 (1 - 2b) ÔH

Ôt Ôx

ÔC

+

H Ôy

ÔV

+

H Ôx

1

(1)

= (1 - 2b)^ + — Ôt Ôx H

1 ( Ôv2 ÔH \

— +-Vj

Ôy Ôx 1

ÔC12 ( 1 ÔV, ÔV.

Ôt

__1+ ÔV^ - ÔH

\

H Ôy Ôx H Ôx

Here b =

1 - 2v 2(1 -v)

; v - Poisson ratio; material coordinates x, y are related to characteristic size of

working L = S , m; S - working cross-sectional area, m2; components of velocity vector v1 and v2 -V %

for longitudinal waves propagation velocity in rock mass c =.

E (1 -v)

p(1 + v)(1 - 2v)

, m/s; E - Young's

modulus, Pa; p - medium density, kg/m3; stress components are related to a variable pc2; velocity vector components - to a variable c; time t - to a variable L/c.

The system of first-order differential equations in partial derivatives (1) can be written in the matrix form:

ÔU _ ÔU nÔU TJ

-+ B-+ pC-+ qQU = 0,

Ôt Ôx Ôy

(2)

r „ , 1 1 ÔH

where U = {v1,v2,c11,c22,c12}r - vector of unknowns; p = — ; q =--

H H Ôx

1 Ôp p Ôx

2

406 -

Journal of Mining Institute. 2017. Vol. 226. P. 405-411 • Mining

Aleksandr P. Gospodarikov, Yaroslav N. Vykhodtsev, Mikhail A. Zatsepin

Mathematical Modelling of Seismic Explosion Waves Impact on Rock Mass ...

The fifth-order constant matrixes present in matrix equation (2), have the form

( 0 0 1 0 01 0 0 0 0 1

B = -

1

0000

1 - 2b 0 0 0 0

0

b000

C = -

(0 0 0 0 0 0 0 0 1 0 0 1 - 2b 0 0 0

1

0

000 000

Q = -

(0 0 1 -1 01

0 0 0 0 2

1 - 2b 0 0 0 0

1 0 0 0 0

V 0 -b 0 0 0,

To close the system (1) with boundary conditions, the following conditions are considered in the paper:

1. Boundary conditions on the cavity surface (a mine working):

or in the matrix form

where rectangular matrix S has a form

Periodicity condition

c,, = c,J = 0

11 lx=0 12 \x=0

SU\ = 0.

x=0

(3)

S =

(0 0 1 0 01 0 0 0 0 1

U = U , •

ly=0 ly=l '

(4)

where l - line G curve length.

2. The initial conditions at the initial time t = 0 determine the stress and velocity fields by for mulas

o0j =a°(c 0 -o(sin2 e- b cos2 e) a02 =-a°(C0 -C)(1 -b)sinecose 0 -C)(cos2 e + b sin2 e)

c^-, = c

(5)

v? = -c0(^0 -Ç)sin0; v2° =c0(^0 -Ç)cos0,

f (s), s > 0,

where a0(s) = \J vu/' "; function f(s) sets the diagram of the incident seismic blast wave; s - dis-

[ 0, s < 0

tance of the medium point to the wave front at time t = 0.

Solution of the problem. For the numerical solution of the boundary value problem (2) - (5) (differential matrix equation with the corresponding boundary and initial conditions) the finite difference method is constructed. The calculated difference scheme is constructed. In this case, we write equation (2) in a divergent form [4]:

U I d(BU) | C(pCU) + TU = 0 dt dx Cy

(6)

where T = qQ -CpC .

Cy

Further, the domain of varuables variation x, y (computational domain) is divided into rectangles by straight lines x = xj ( j = 1,2,...,J) and y = yn (n = 1,2,...,N), and in space (x, y, t), allocating

an elementary parallelepiped V, bounded by the planes x = xj-1, x = xj, y = yk-1, y = yk, t = t,

t = t'+r, we obtain a finite-difference scheme for the set boundary-value problem (Fig.2) for determining main parameters of stress condition of rock mass, including a working. Integrating the matrix equation (6) by volume (parallelepiped V)

CU 1 CBU 1 CpCU cy

\

- +

Ct Cx

+ -

dV = - TUdV

V \ S J V

and transforming the left-hand side of the last expression according to Gauss's formula: Journal of Mining Institute. 2017. Vol. 226. P. 405-411 • Mining

^AleksandrP. Gospodarikov, YaroslavN. Vykhodtsev, Mikhail A. Zatsepin

Mathematical Modelling of Seismic Explosion Waves Impact on Rock Mass ...

t = t' + x

yn-i V

hX, ^ Sjn n-i t = t'

hy yn

Fig.2. Finite-difference scheme

t' + T

v* > 0/ /V . /V ■ / ' x j, n < 0 Vx,j + i,

x:-i

xj+i

Fig.3. Riemann invariants

U+B SU=0,

St 8x f U,n, x < x,

I [U cos(n, t) + BU cos(n, X) +

S

+ pCU cos(n, * )]dS = -J TUdV,

V

where S - surface of the examined parallelepiped V, i.e. S = S: u Sjn u

u Sjn u Sj-i,n u S:n U S:n-i; n - directi°n

of the outer normal line to it, as well as assuming that on each face a vector of unknowns U keeps a constant value to a precision of first-order small quantities, we have

jn

(U: - Un )hxhy + B(U,n - Uj-i,n)hyx +

+ C(P,nU,n - P,,n-iU,,n-1 )hxX =

= -(TU):nhxhyx,

(7)

U:, U n - values of vector of unknowns

n

: S]n ac-

U at upper and lower faces S cordingly; U:n, U:-i n - values of vector of unknowns U at lateral faces Sjn, Sj-hn, perpendicular to axis Mix;

Ujn, Ujn n-i - values of vector of unknowns U at lateral faces Sjn, Sjn-i, perpendicular to axis Miy.

Variables Ujn, U Jn are defined by

splitting method [3, ii], in accordance with which the values of the unknown vector U at lateral faces are calculated by solving spatial one-dimension equations [2, 3]. Then for determining the unknown vector U we consider the boundary-value problem in the following form

(8)

t < t < t'+x .

Denoting through Ax a matrix of left eigenvectors of the matrix B, corresponding to its eigenvalues (k = i, 2, ..., 5), and assuming U = A-.Vx, from equation (8) we have

SVx wSVx n —x + M—x = 0 St Sx

(9)

where M = AxBAx1 = diag(^i, ^2,..., .

From equation (9) it follows, that components V^ (Riemann invariants) of vector Vx keep constant values at straight lines ^kt-x = const (Fig.3).

t

0

t

t'

0

x

x

^AleksandrP. Gospodarikov, YaroslavN. Vykhodtsev, Mikhail A. Zatsepin

Mathematical Modelling of Seismic Explosion Waves Impact on Rock Mass ...

Thus, we have

**) = { ^ * 0;

xj [j, ^ < 0. The last expression we write in the matrix form:

K, jn = P+K, jn + P~VjU, (10)

where P+ = 1(|P| ± P); P = diag^ign^)}; P = diag{|Ptt |}.

Finally, the solution of the difference equation (7) takes the form

Ujn =

X X

E--O--(p O++ p )

h h i y 1 ,n~i y'

hx hy

X

Ujn + — (O+Uj_i,n + O"Uj+i,n) +

x

X

+ Y (phn_O+yUhn_i + Pjn O_Uj,n+i) _ (TU) in X , (11)

y

where

O± = AxM±A_x1; O± = AyM±A_/; M+ = !(|m| + M);

O = 0+ + O- = A Ml A-1 • O = 0+ + O- = A Ml A-1

x ^x^^x -"-x^ | x J y y y yr I y '

Formula (11) allows to obtain an expression for the values of the unknown vector U at the inner nodes of the upper layer corresponding to the time instant (t' + x), through the values of this vector at the nodes of the lower layer for the time instant t'.

Then, with the help of the obtained formula (11), a transition from the time layer corresponding to the instant of time t = t is made to the next time layer t = t'+x. We note that formula (11) allows us to find new values of the vector of unknowns U only in the inner cells of the computational domain, i.e. in cells that do not share points with the cavity boundary. For this reason, in order to find a solution in boundary cells, it is necessary to involve the boundary conditions (3) - (4).

Denoting by S1 p,p e 1: N cells adjacent to a cavity boundary, the boundary conditions (3) on the surface will have the following form

Sf>0,p = 0, (12)

or shifting to (12) to Reimann invariants:

SA-Vx ,0, p = 0. (13)

The matrix equality (12) have two conditions for determining the component of the unknowns vector f/0 . We obtain three more necessary relations from the condition of saving the Reimann invariants, corresponding to the nonpositive eigenvalues of matrix B:

Vx,0,p = Vx(,i)p ; k = 3,4,5, (14)

or in matrix form we can write it in the form

GVx,0,p = GVx,1,p , (15)

(0 0 1 0 0^ 0 0 0 1 0

0 0 0 0 1

V /

Combining the matrix equations (13) and (15), we obtain the required system of algebraic equations for the determination of all components of the vector of unknowns Vx,0, p , the solution of which is written in the form

where matrix G =

Aleksandr P. Gospodarikov, Yaroslav N. Vykhodtsev, Mikhail A. Zatsepin

Mathematical Modelling of Seismic Explosion Waves Impact on Rock Mass ...

V

'SA"1 V1 ( 0 ^

x.0. p

V G V

V GVx .1. p y

= Y .

The last three components of the vector Vx,0,p are defined. as before. by the relations (14) Thus. the vector Vx 0 can be represented in the form

V = P~V + P+Y

y x.0.p r y x.1.p ^ r 1 •

Then we introduce ghost cells S0,p and assume Vx,0, = Y . We have

0. p

Vx .0. p = P+VX .0. p + P'Vx .1. p I Up =A+xU0.p +A"xU!.p .

Where A±x = A-1P ±A x.

U 0 p =A"1Y = A x

-1

v G y V

0 ^

GU , p

^SA-1 Ï V1 ( S Y1 ( 0 >

A

VV G y y

V GA x y

GU

1.p

(16)

(17)

f n \

where matrix G A x = eigenvalues.

Ve5 y

; e3. e4. e5 - eigenvectors of matrix B. corresponding to its nonpositive

Fig.4. Calculation configuration

Fig.5. Visualization of a process with defined points

Fig.6. Visualization of a scenario: defined points and stress diagrams (left), velocity diagrams (right)

3

e

4

J\ AleksandrP. Gospodarikov, YaroslavN. Vykhodtsev, Mikhail A. Zatsepin DOI: 10.25515IPMI.2017.4.405

Mathematical Modelling of Seismic Explosion Waves Impact on Rock Mass ...

In addition, we introduce ghost cells Sj0, SjN ( j e 1: J ) and assume U0,j = Un-1j, UN,j = Uu. We

note that formulas (16) completely coincide with formulas (10), and, hence, the relation (11) is valid for the entire calculated area.

To solve the boundary value problem, we use a step-by-step algorithm. If on the time layer corresponding to the instant of time t = t', the state of the medium is already known, then to shift to the next time layer t = t' + x it is necessary with formulas like (17) to fill in the ghost cells and later use the relation (11). Moreover, as it follows from (10), during each transition the length of the calculated area decreases by one cell in the direction of axis Ox (see Fig.2). Therefore, if it is required to

find a solution in the arbitrary rectangle \x e [0, x ]; for all solutions t e [0, T], then the initial calcul y e [0, l]

lated area should at least occupy a rectangle \ x e [0, x + T ];

l y e [0,l].

Figures 4-6 show the results of the solution of the model boundary-value problem with the following parameters: working in the form of an ellipse with semi-axes 2 and 4 m is located in rocky ground (Young's modulus E = 57.9 GPa, Poisson's ratio v = 0.35, velocity of longitudinal waves v = 5800 m/s), the diagram of the incident wave is shown in Fig. 4.

Conclusions

1. A mathematical model of the seismic explosion wave impact on a mine working has been obtained.

2. An original calculation-difference scheme for the solution of boundary-value problem under consideration has been developed, which describes the impact of a seismic explosion wave on rock mass, including mine working.

3. An algorithm for the numerical solution of a boundary-value problem that realizes the obtained mathematical model has been developed.

4. An effective software product has been developed in the JavaScript language, which was tested for the solution of a typical model problem.

REFERENCES

1. Vallander S.V. Lectures on hydroaeromechanics. Leningrad: Nauka, 1978, p.296 (in Russian).

2. Godunov S.K., Ryaben'kii V.S. Differential schemes. Moscow: Nauka, 1977, p.440 (in Russian).

3. Gospodarikov A.P., Kolton G.A., Buldakov E.L. One approach to mathematical modeling of the effect of blast waves on an underground oil pipeline. Zapiski Gornogo instituta. 2014. Vol. 210, p. 37-42 (in Russian).

4. Gospodarikov A.P., Gorokhov N.L. Dynamic calculation of pipelines for seismic impacts. Zapiski Gornogo instituta. 2011. Vol. 193, p. 318-321 (in Russian).

5. Lyakhov G.M. Fundamentals of the dynamics of blast waves in soils and rocks. Moscow: Nedra, 1974, p. 192 (in Russian).

6. Novozhilov V.V. Theory of elasticity. Leningrad: Sudpromgiz, 1958, p. 371 (in Russian).

7. Filonenko-Borodich M.M. Theory of elasticity. Moscow: Fizmatlit, 1959, p. 361 (in Russian).

8. Godunov S.K., Zabrodin A.V., Ivanov M.Ya., Kraiko A.N., Prokopov G.N. Numerical solution of multidimensional problems of gas dynamics. Moscow: Nauka, 1976, p. 400 (in Russian).

9. Shemyakin E.I. Dynamic problems of the theory of elasticity and plasticity. Novosibirsk: Izd-vo NGU, 1968, p.337 (in Russian).

10. Bormann P., Engdahl E.R., Kind R. Seismic wave propagation and Earth models. Ed. Bormann. New Manual of Seismol-ogical Observatory Practice. Potsdam: German Research Center for Geosciences, 2012, p. 1-105. DOI: 10.2312/GFZ.NMS0P-2_ch2.

11. Yan Bo, Xinwu Zeng, Yuan Li. Subsection forward modeling method of blasting stress wave underground. Mathematical problems in engineering. 2015. Vol. 215, p. 9. DOI: 10.1155/2015/678468.

12. Ziaran S., Musil M., Cekan M., Chlebo O. Analysis of seismic waves generated by blasting operations and their response on buildings. International Journal of Environmental, Chemical, Ecological, Geological and Geophysical Engineering. 2013. Vol.7. N 11, p. 769-774.

Authors: Aleksandr P. Gospodarikov, Doctor of Engineering Sciences, Professor, kafmatem@spmi.ru (Saint-Petersburg Mining University, Saint-Petersburg, Russia), Yaroslav N. Vykhodtsev, Postgraduate Student, mm1000@list.ru (Saint-Petersburg Mining University, Saint-Petersburg, Russia), Mikhail AZatsepin, Candidate of Pshysics and Mathematics, Associate Professor, zatsepin@spmi.ru (Saint-Petersburg Mining University, Saint-Petersburg, Russia). The paper was accepted for publication on 17 April, 2017.