Математическое моделирование воздействия сейсмовзрывных волн на горный массив, включающий выработку Текст научной статьи по специальности «Математика»

А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин DOI: 10.25515/РМ1.2017.4.405

Математическое моделирование воздействия сейсмовзрывных волн на горный массив...

Горное дело

УДК 622.235.535.2

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВОЗДЕЙСТВИЯ СЕЙСМОВЗРЫВНЫХ ВОЛН НА ГОРНЫЙ МАССИВ, ВКЛЮЧАЮЩИЙ ВЫРАБОТКУ

А.П.ГОСПОДАРИКОВ, Я.Н.ВЫХОДЦЕВ, М. А. ЗАЦЕПИН

Санкт-Петербургский горный университет, Санкт-Петербург, Россия

В статье в рамках динамической теории упругости представлена математическая модель воздействия сейсмовзрывных волн на массив горных пород, включающий выработку. Увеличение объемов добычи полезных ископаемых в сложных горно-геологических условиях с учетом влияния энергии взрыва тесно связано с анализом основных параметров напряженно-деформированного состояния массива горных пород, включающего выработку. Последнее приводит к необходимости определения безопасных параметров буровзрывных работ, обеспечивающих эксплуатационное состояния горной выработки. Основную опасность при взрыве заряда взрывчатого вещества вблизи действующей выработки представляет сейсмовзрывная волна, характеристики которой определяются свойствами грунта и параметрами буровзрывных работ. Определение полей напряжений и скоростей смещений в массиве горных пород требует привлечения для своего решения современного математического аппарата. Для численного решения поставленной краевой задачи методом конечных разностей авторами построена оригинальная расчетно-разностная схема. Применение метода расщепления для решения двухмерной краевой задачи сводится к решению пространственно одномерных дифференциальных уравнений. Для полученного численного алгоритма разработана эффективная вычислительная программа. Приведены численные решения модельной задачи для случая, когда форма горной выработки представлена эллипсом.

Ключевые слова: горная выработка, математическая модель, сейсмовзрывная волна, разностная схема, численный алгоритм

Как цитировать эту статью: Господариков А.П. Математическое моделирование воздействия сейсмовзрывных волн на горный массив, включающий выработку / А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин // Записки Горного института. 2017. Т. 226. С. 405-411. DOI: 10.25515/РМ1.2017.4.405

t = 0

Введение. Напряжения, возникающие в окрестности горной выработки при сейсмическом воздействии взрыва на нее, зависят от многих причин, но, в первую очередь, от мощности заряда взрывчатого вещества (ВВ), расстояния до выработки, скорости детонации, смеси заряда ВВ и т.д. [5, 6, 10]. Для определения влияния воздействия взрывной волны на горную выработку необходимо провести большое количество натурных испытаний, что экономически и технически не всегда возможно. Поэтому для оценки взрывного воздействия на горную выработку в работе используется численное моделирование взаимодействия продольной волны, распространяющейся в упругой среде, с выработкой [12].

Постановка задачи. В качестве математической модели воздействия взрывных волн на горную выработку в работе используются уравнения динамической теории упругости Мизеса [3, 9], записанные в криволинейных координатах (рис.1).

На рисунке представлены две системы координат, где введены следующие обозначения:

1) O\zn - прямоугольная система координат; начало координат Ol -центр масс горной выработки; ось Оц параллельна скорости перемещения

фронта невозмущенной волны; Рис. 1. Введенные системы координат

О

А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин

Математическое моделирование воздействия сейсмовзрывных волн на горный массив.

DOI: 10.25515/РМ1.2017.4.000

2) Ы\ху - криволинейная система координат; x - расстояние (ЫЫ\) от точки Ы до границы горной выработки; у - длина дуги Г, отсчитываемая от точки О до точки Ы\, О - точка встречи фронта С\ падающей волны с граничной поверхностью выработки в начальный момент времени (7 = 0).

С учетом введенных систем координат имеем соотношение вида

Г = Я+хп,

где г - радиус-вектор точки Ы; Я - радиус-вектор точки Ы\, п - орт нормали п, т - орт касательной т в точке Ы\. Тогда справедливо

ёг = ёЯ + хёп + ёхп = ёут + ёхп + k(у) хёут = (1 + k(у) х)ёут + ёхп ,

где к(у) - кривизна линии Г в точке Ы1. Окончательно получим

.,2 , тт2 1,2

(ёг )2 = ёх2 + (1 + k (у) х)2 ёу2 = ёх2 + Н

где Н = 1 + к (у) х - коэффициент Ламе.

Введем обозначения: v1, VI - компоненты вектора скорости частицы массива горных пород по координатным осям Ы1х и Ы\у соответственно; а через ац, о12, а22 - компоненты тензора напряжений. Тогда уравнения движения и продифференцированный по времени закон Гука, записанные в безразмерной форме, примут вид [1, 7, 8]:

ду да1

+

1 да

12

1 дН . .

(аи-а22);

д7 дх Н ду Н дх

дУт да„ 1 да„ 2 дН

—2 = —— +--22 +--а12;

д7 дх Н ду Н дх

дап _ ду + (1 - 2Ь) ду2 + (1 - 2Ь) дН у ;

д7 да

дх

Н ду ду

Н дх

(1)

= (1 - 2Ь)—1 + — д7 дх Н

1 (ду дН ^ — +-у

ду дх

да12 ( 1 ду, ду.

д7

__1 +ду2 - дН

\

Н ду дх Н дх

Здесь Ь =

1 - 2у 2(1 -V)

; у - коэффициент Пуассона; материальные координаты х, у отнесены к ха-

5

рактерному размеру горной выработки L = Л — , м; 5 - площадь выработки, м ; компоненты век-

%

тора скорости VI и у2

к скорости распространения продольных волн в массиве

с = -Е-——-, м/с; Е - модуль Юнга, Па; р - плотность среды, кг/м3; компоненты тензора

р(1 +—)(1 - 2—)

2.

напряжений отнесены к величине рс ; компоненты вектора скорости - к величине с; время 7 - к величине Ь/с.

Систему дифференциальных уравнений первого порядка в частных производных (1) можно записать в матричной форме:

ди ди ди тт

-+ В-+ рС-+ qQU = 0,

д7 дх ду

(2)

т 1

где и = {у,у2,а11,а22,а12} - вектор неизвестных; р = —; q =

1 дН 1 др

Н Н дх р дх

2

А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин

Математическое моделирование воздействия сейсмовзрывных волн на горный массив...

DOI: 10.25515/PMI. 2017.4.405

Постоянные матрицы пятого порядка, присутствующие в матричном уравнении (2), имеют вид

( 0 0 1 0 0 0 0 0 0 1

B = -

1

0000

1 - 2b 0 0 0 0

0

b000

C = -

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

1 0 0 0 0 0 0 0

(0 0

Q = -

0 1 -10 ^ 0 0 0 2

1 - 2b 0 0 0 0

1

0 0 0 0 -b 0 0 0

Для замыкания системы (1) краевыми условиями в работе рассматриваются следующие условия: 1. Граничные условия на поверхности полости (горной выработки):

с,, = с, J = 0

11lx=0 12 lx=0

или в матричном виде

где прямоугольная матрица £ имеет вид

S =

Условие периодичности

SU\ = 0.

x=0

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

и = U , ■

ly=0 ly=l '

(3)

(4)

где I - длина дуги линии Г.

2. Начальные условия в начальный момент времени t = 0 определяют поля напряжений и скоростей по формулам

c0j = с0(С 0-C)(sin2 е- b cos2 е); с02 = -с0(С0 -с)(1 -b)sinеcosе;

Cv, = с

'(С0 -C)(cos2 е + bsin2 е);

(5)

v? = -с0(С0 -с)sinе; v0 = с0(С0 -с )C0S е,

f (s), s > 0

; функция fs) задает эпюру падающей сейсмовзрывной волны; s - рас-

где а0(5) =,

0, 5 < 0

стояние точки среды до фронта волны в момент времени t = 0.

Решение задачи. Для численного решения краевой задачи (2)-(5) (дифференциальное матричное уравнение с соответствующими граничными и начальными условиями) методом конечных разностей в работе построена расчетная разностная схема. В этом случае уравнение (2) запишем в дивергентной форме [4]:

ди д(ВП) д(рСП)

dt

- + -

dx

- + -

dy

- + TU = 0,

(6)

где Т = qQ-дрС.

ду

Далее область изменения переменных х, у (расчетная область) разбивается на прямоугольники прямыми линиями х = ху (у = 1,2,...,J) и у = уп (п = 1,2,...,N), а в пространстве (х, у, 0, выделяя элементарный параллелепипед V, ограниченный плоскостями х = Ху-1, х = Ху, у = ук-1, у = ук, t = t', t = , получим конечно-разностную схему поставленной краевой задачи (рис.2) для определения основных параметров напряженного состояния массива горных пород, включающего выработку.

Интегрируя матричное уравнение (6) по объему (параллелепипед V)

dU t dBU , dpCU dy

- +

dt dx

+ -

dV = - J TUdV

V

и преобразовав левую часть последнего выражения по формуле Гаусса:

0

А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин

Математическое моделирование воздействия сейсмовзрывных волн на горный массив...

DOI: 10.25515/PMI. 2017.4.000

+

t' + т

t = t' + х

Уп-1 - — J V

hS: Г Snn п-1 t = t'

hy Уп

Рис.2. Конечно-разностная схема

и* > о/ /V . /V ■ / ' x j, п < 0 Vx,j + 1,

xj-1 Xj Xj+J

Рис.3. Инварианты Римана

U (x, у,, С) = ■

SU+B S-U=0

St 8x ' U,n, X < Xj

| [и cos(", ^ + Ви cos(", X) + реи шз(й, .у = TUdV,

V

где 5 - поверхность рассматриваемого параллелепипеда V, т.е. 5 = 5]п и 5]п и

и 51"и 51-1,пи 51"и 51,п-1; " - направление внешней нормали к ней, а также считая, что на каждой грани вектор неизвестных и сохраняет постоянное значение с точностью до малых первого порядка получим

(и1 - и]п + Вф"п - и1 -!," ят +

+ С (Р1пи1п - Р1,п-1и1,п-1 Ж Т =

= -(Ти) т, (7)

и1, ип - значения вектора неизвестных и на верхней и нижней гранях Б1, 5п соответственно; и п , и,-1 п - значения вектора неизвестных и на боковых гранях 5п, 5п-1п, перпендикулярных

оси М1Х; и 1п, ип, п-1 - значения вектора неизвестных и на боковых гранях

, 5у,п-1, перпендикулярных оси М\у.

Величины ип , и 1п определяются

методом расщепления [3, 11], в соответствии с которым значения вектора неизвестных и на боковых гранях находятся из решения пространственно одномерных уравнений [2, 3]. Тогда для определения вектора неизвестных и рассматривается краевая задача вида

(8)

Uj+1,n , x > Xj

; t < t < t'+х .

Обозначая через Лх матрицу левых собственных векторов матрицы В, отвечающих ее собственным значениям ^ (к = 1, 2, ..., 5), и полагая и = ЛС^¥Х, из уравнения (8) получим

dFX wSVx -

—x + M—- = 0 St Sx

(9)

где M = ЛxBA;: = diag^ ц2,..., .

Из уравнения (9) следует, что компоненты VX (инварианты Римана) вектора Vx сохраняют постоянные значения на прямых линиях - x = const (рис.3).

t

о

t

t'

о

x

А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин

Математическое моделирование воздействия сейсмовзрывных волн на горный массив.

DOI: 10.25515/РМ1.2017.4.405

Следовательно, имеем

*к) = 1 ^ Цк *0;

[У, Цк < 0.

Последнее выражение запишем в матричном виде:

х, уп х, уп х, у+1,п

где р+ = 2(|Р| ± Р); Р = ^^(Цк)}; |р| = ^{ркк |}.

Окончательно решение разностного уравнения (7) примет вид

иуп =

Vх, уп = пп + Р^х, ,.+1_ , (10)

X X

Е--Ф--(р Ф++ р ,Ф-)

и и } У у ,п-1 У'

Пх Пу

х

иуп + — (Ф+и;-1,п + Ф-и;+1,п ) +

х

+ т~ (Р;,п-1Ф+и;,п-1 + р]п Ф-и;,п+1) - (Ти) уп X, (11)

у

где

ф^ = а м ±л;1; ф;у=л уМ ±л-;; м+ = 2(|м|+м);

Ф х = Ф+ + Ф- = л х|М| л;1; фу = Ф++ + Ф- = л у М1 А1.

Формула (11) позволяет получить выражение для значений вектора неизвестных и во внутренних узлах верхнего слоя, отвечающего моменту времени (V + х), через значения этого вектора в узлах нижнего слоя для момента времени ^.

Далее с помощью полученной формулы (11) осуществляется переход с временного слоя, отвечающего моменту времени t = t', на следующий временной слой t = ^+х. При этом отметим, что формула (11) позволяет найти новые значения вектора неизвестных и только во внутренних ячейках расчетной области, т.е. в ячейках, не имеющих общих точек с границей полости. Поэтому для нахождения решения в граничных ячейках необходимо привлечь краевые условия (3)-(4).

Обозначая через 51 р, р е 1: N ячейки, примыкающие к границе полости, краевые условия (3)

^0,р = 0, (12)

, р

на поверхности полости примут вид

или перейдя в (12) к инвариантам Римана:

ял-^хдр = 0. (13)

Матричное равенство (12) имеет два условия для определения компонент вектора неизвестных и0 . Еще три необходимых соотношения получим из условия сохранения инвариантов Римана, отвечающих неположительным собственным значениям матрицы В:

Vx,0,p = V®р ; к = 3,4,5, (14)

или в матричной форме данное равенство запишем в виде

^,0,р =Г^,1,р, (15)

'0 0 1 0 0^ где матрица Г = 0 0 0 1 0 ч0 0 0 0 1 у

Объединяя матричные уравнения (13) и (15), получим искомую систему алгебраических уравнений для определения всех компонент вектора неизвестных Ух0 р , решение которой записывается в виде

А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин

Математическое моделирование воздействия сейсмовзрывных волн на горный массив...

DOl: 10.25515/РМ1.2017.4.000

/„.-Л-1

V

X ,0, р

5Л

V Г ,

0

VГVX ,1, р ,

= Y .

При этом последние три компоненты вектора V. 0 по-прежнему определяются соотношениями (14). Таким образом, вектор Ух 0 представим в виде

V = р ~у + р+у

у х,0,р Г у х,1,р ^ г 1 •

Далее введем фиктивные ячейки 50 р и положим V. 0 р = Y . Тогда получим

V ,0, р = р X ,0, р + р V д, р; I Ц>р =Л+.и0,р + Л;и1,р .

Здесь Л* = Л.Р ±Л

и 0 р =Л-х1У = Л;1

^Г1 Г 0 >

V Г V

V" р

5Л-Г

-1

Л

^ у у

5

VГЛ х ,

-1

ги1

(16)

(17)

, р у

где матрица ГЛ х =

V е5 у

; ез, е4, е5 - собственные векторы матрицы В, отвечающие ее неположи-

тельным собственным значениям.

Рис.4. Конфигурация расчета

Рис.5. Визуализация процесса с реперными точками

Рис.6. Визуализация сценария: реперные точки и графики напряжений (слева), графики скоростей (справа)

X

X

0

е

4

Ж А.П.Господариков, Я.Н.Выходцев, М.А.Зацепин 001: 10.25515/РМ1.2017.4.405

чи Математическое моделирование воздействия сейсмовзрывных волн на горный массив...

Кроме того, вводим фиктивные ячейки Sj0, SjN (у е 1: J) и положим и0 у = ип-1 у, UN у = и1 у.

Отметим, что формулы (16) полностью совпадают с формулами (10), и, следовательно, соотношение (11) справедливо для всей расчетной области.

Для решения поставленной краевой задачи используем пошаговый алгоритм. Если на временном слое, отвечающем моменту времени t = t', состояние среды уже известно, то для перехода на следующий временной слой t = ^ + х необходимо с помощью формул типа (17) заполнить все фиктивные ячейки, а затем воспользоваться соотношением (11). При этом, как следует из (10), при каждом таком переходе протяженность расчетной области уменьшается на одну ячейку в направлении оси Ох (см. рис.2). Поэтому, если требуется построить решение в произвольном

прямоугольнике \ х е [0, х ]; для всех значений t е [0, Т], то начальная расчетная область должна,

[ у е [0,1]

„ „ |хе [0,х1 + Т];

по крайней мере, занимать прямоугольник <

1 у е [0,/].

На рис.4-6 представлены результаты решения модельной краевой задачи со следующими параметрами: горная выработка в форме эллипса с полуосями 2 и 4 м расположена в скальном грунте (модуль Юнга Е = 57,9 ГПа; коэффициент Пуассона V = 0,35; скорость продольных волн у = 5800 м/с), эпюра падающей волны изображена на рис.4.

Выводы

1. Получена математическая модель воздействия сейсмовзрывной волны на горную выработку.

2. Разработана оригинальная расчетно-разностная схема решения рассматриваемой краевой задачи, описывающая воздействие сейсмовзрывной волны на массив пород, включающих горную выработку.

3. Построен алгоритм численного решения краевой задачи, реализующий полученную математическую модель.

4. Разработан эффективный программный продукт на языке JavaScript, апробированный на решении характерной модельной задачи.

ЛИТЕРАТУРА

1. Валландер С.В. Лекции по гидроаэромеханике. Л.: Наука, 1978. 296 с.

2. Годунов С.К. Разностные схемы / С.К.Годунов, В.С.Рябенький. М.: Наука, 1977. 440 с.

3. Господариков А.П. Об одном подходе к математическому моделированию воздействия взрывных волн на подземный нефтепровод / А.П.Господариков, Г.А.Колтон, Е.Л.Булдаков // Записки Горного института. 2014. Т. 210. С. 37-42.

4. Господариков А.П. Динамический расчет трубопроводов на сейсмические воздействия / А.П.Господариков, Н.Л.Горохов // Записки Горного института. 2011. Т. 193. С. 318-321.

5. ЛяховГ.М. Основы динамики взрывных волн в грунтах и горных породах. М.: Недра, 1974. 192 с.

6. НовожиловВ.В. Теория упругости. Л.: Судпромгиз, 1958. 371 с.

7. Филоненко-БородичМ.М. Теория упругости. М.: Физматлит, 1959. 361 с.

8. Численное решение многомерных задач газовой динамики / С.К. Годунов, А.В. Забродин, М.Я. Иванов, А.Н. Край-ко, Г.Н. Прокопов. М.: Наука, 1976. 400 с.

9. Шемякин Е.И. Динамические задачи теории упругости и пластичности. Новосибирск: Изд-во НГУ, 1968. 337 с.

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

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

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

Авторы: А.П.Господариков, д-р. техн. наук, профессор, kafmatem@spmi.ru (Санкт-Петербургский горный университет, Санкт-Петербург, Россия), Я.Н.Выходцев, аспирант, mm1000@list.ru (Санкт-Петербургский горный университет, Санкт-Петербург, Россия), М.АЗацепин, канд. физ.-мат. наук, доцент, zatsepin@spmi.ru (Санкт-Петербургский горный университет, Санкт-Петербург, Россия).

Статья принята к публикации 17.04.2017.