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

УДК 519.6

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

Основным способом разрушения скального массива при открытой разработке месторождений полезных ископаемых является буровзрывной. Реализация математический модели развала буровзрывного блока процесс трудоёмкий, требующий большого количества машинного времени. В работе предложены алгоритмы и численные методы компьютерного моделирования взрыва, адекватно описывающие математическую модель, являющиеся простыми в реализации, в том числе с применением методов параллельного программирования на e-mail: kabelko@mail.ru многопроцессорных серверах и графических станциях.

Ключевые слова: алгоритмы моделирования взрыва, численные методы моделирования взрыва, распределение полезного компонента в развале.

С.Г. КАБЕЛКО

ФГУП ВИОГЕМ

В настоящее время доминирующим способом подготовки скальных пород и рудных залежей к выемке на карьерах является буровзрывной. Математическая модель процесса развала буровзрывного блока описывает зависимости количественных показателей природно-технологических факторов влияющих на формирование развала и распределение полезного компонента во взорванной горной массе. Для построения численной модели используется расчетная сетка, которой буровзрывной блок разбивается взаимно перпендикулярными плоскостями на элементарные ячейки (блочная модель). На практике, блочная модель буровзрывного блока с размером элементарной ячейки 1 м3 представляет собой совокупность нескольких миллионов ячеек, на которые оказывают влияние природно-технологические факторы взрыва и которые взаимодействуют между собой в различные моменты времени. Очевидно, что реализация математический модели процесс трудоёмкий, требующий большого количества машинного времени. А с учетом того, что на производстве специалисты получают исходную информацию за сутки до взрыва, для практического применения данной разработки необходимо при выборе алгоритмов и численных методов расчёта руководствоваться критерием быстродействия.

При короткозамедленном взрывании моделирование осуществляется последовательно для каждого ряда буровзрывной сети скважин с выделением трёх этапов:

- расчет распределения потенциалов и начальных скоростей, определение области разрушения;

- дезинтеграция разрушенной части буровзрывного блока;

- баллистический разлет и формирование развала.

На первом этапе для нахождения потенциального поля скоростей необходимо решить уравнение Лапласа[1]:

Фкс(X У.2) + фуу (х. У.2)+ Н (X У.2) = 0, (1)

с граничными условиями на открытой поверхности:

<р\0 = 0, (2)

граничные условия на поверхности п-ой скважины Sn:

Н^ = Яя, (3)

граничные условия на поверхности удалённой более 50 м от взрываемых скважин Sm:

Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1

НS* = 0 . (4)

Для нахождения потенциала на поверхности n-ой скважины фэп для каждой скважины требуется вычисление интеграла[2,з]:

Н _ =- Т 2p<t)<V/V0>' 25 dt, (5)

I (' + P„D/(P,DC)PC

где p(t) - мгновенное давление в продуктах взрыва, определяемое по формуле[з]:

p(t) =1 рВВ D8/3r2/3t2/3eUDt/L)2/3,

8

где Тв - время до вылета забойки, t(z) - момент времени, при котором нижний срез забойки проходит расстояние z от своего первоначального положения, D - скорость детонации, V,Vo - конечный и начальный объемы продуктов взрыва, Dc - скорость упругой волны в среде, L - длина заряда, рвв - плотность заряда, рс - плотность среды. Уравнение для определения функции t(z) имеет вид[з]:

3 2

d2t/dz2 = Pmaxh3 g (dt/dz)—, (6)

(hзаб -Z)(hз + Z) Yзаб

где, hзаб,hз - соответственно длина забойки и заряда; Pmax = %D2p; D - скорость детонации заряда; рвв - плотность заряда; g - ускорение свободного падения; узаб - удельный вес забойки.

Уравнение (6) для каждой n-ой скважины решается на отрезке [о, Ьзабп-Az] (нижняя кромка забойки смещена в начало координат) при начальных условиях t(o) = 0,6 мс, dt/dz(o) = i/Укр. Здесь Az - шаг разбиения интервала [о, Ьзабп] на большое число промежутков, Укр - скорость нижней части забойки при ее срезе. Это значение достигается через о,2-о,6 мс после детонации заряда ВВ и величина Укр для различных видов забоек лежит в диапазоне (20-200) м/с. Для различных видов забоек значения Укр приведены в таблице [3].

Для численного решения уравнения (6) применен алгоритм метода Рунге-Кутта [4]. Введем следующие обозначения:

z0 = zk = kAz,t (zk ) = tk,t\zk ) = tk'

Pmax h3 g (tk ')

(Кабп - Zk )(h3 + Zk Y Узаб

з = f (zk, h ')■

t = t + t Az +—1------------- -3 Az *

lk+1 k k ;

Тогда вычислительные формулы алгоритма имеют следующий вид:

а + а2 + а 6

. . а + 2а2 + 2а3 + а4

г к+1 = г к + -1-2-----3—4,

6

где аі = /(^к, гк )Лг;

а2 = /(zk + Лг /2, гк + а1 / 2)Лг; а3 = /(гк + Лг /2, гк + а2 / 2)Лг ;

а4 = /(гк + Лг,гк + а3)Лг , к = 0,1,2,3,....

Применяя этот алгоритм, находим с некоторой погрешностью время до вылета забойки для каждой скважины. Для достижения требуемой точности вычисляем значения Тві при Д7=Ьз/М, Д7=Ьз/(2Ы), Д7=Ьз/(4М) и.т.д. и, сравнивая при этом значения Тві, Тв2, Твз ,... Счет прекращается при выполнении условия | Тві- Тві+і | /Тві< 0,01.

Интеграл нахождение величины ф8п (5) - функция координат ячеек блока, через который проходит та или иная скважина. При разбиении блока на элементарные

ячейки граничные условия можно поставить в центрах таких ячеек. Вычисляемый интеграл зависят от Тв и функции 1(7), которые определяются при численном решении уравнения (6). Разобьем интервал интегрирования Тв-Ь^) на четное (п=2ш) число частей с шагом Ь. Применим для вычисления интеграла обобщенную формулу Симпсона [4] с последующим методом половинного деления. Тогда, имеем ниже описанный алгоритм расчета.

На первом шаге приближенное значение интеграла вычисляем по формуле к

*/1 = 3(^0 + У2т + 4(у1 + У3 + ••• + У2т-1) + 2(у2 + У4 + - + У2т-2)) ,

- * (^)

где h =

2m

= 2p(tt )(V / Vo)

y t

1.25

1 + РввD /(PcDc )

Для уточнения значения Jibh применим удвоение разбиения интервала [о, Тв-t(z)]. Погрешность определим как | Jk- Jk+i | / Jk, где Jk- значение интеграла на к-ом шаге удвоения. Счет прекращаем при значении погрешности < 0,01. Pc,Dc,p - величины, зависящие от координат и, в конечном итоге, от номера ячейки. Можно считать их постоянными в пределах ячейки. Эти распределения должны быть заданы.

Численное решение уравнения Лапласа (1) с граничными условиями (2,3,4) при зависимости потенциала от двух переменных (плоская задача теории потенциала) используют методы конечных разностей [5,6], конечных элементов [7], граничных элементов [8].

Отличительной особенностью этой задачи в настоящей работе является то, что искомый потенциал ф(х,у^) зависит от трех переменных и задача ставится в многосвязной области. При этом возникают большие сложности в построении расчетной схемы и вычислительного алгоритма. Следует также учесть большое число получающихся расчетных узлов при дискретизации рассматриваемого объема. Поэтому, для решения задачи Дирихле применим простейший алгоритм метода установления, используемый при исследовании явлений установившейся диффузии [9].

Вычисления проводятся методом итераций. В случае пространственной задачи расчетные формулы для внутренних узлов области имеют вид:

.. m | m | m . m . m . m

m+1 ^t+1 jk + 9t-1 jk + Vij+1k + Vij-1k + Vijk+1 + Vijk -1

Vyk =-----------------------6--------------------,

где 9ijkm+1 - потенциалы в центрах ячеек блока на (т+1)-ой итерации;

Фi+iikm,Фi-iikm,Фij+ikm, 9ij-ikm, 9ijk+im, 9ijk-im - потенциалы в центрах ячеек блока на т-ой итерации.

Смысл метода итераций состоит в то, что мы заменяем стационарный процесс, процессом установления, приняв нулевое приближение за начальное распределение потенциала. Выражение в правой части является приближенным представлением в объеме теоремы о среднем для гармонической функции. Удобство этого метода расчета обусловлено тем, что расчетные формулы предельно просты. Процесс итераций сходящийся [9]. Для повышения скорости счета в задаче Дирихле используются неявные вычислительные схемы, упрощения вычислительных процедур достигается применением метода дробных шагов или переменных направлений [5,6]. Для поставленной задачи ускорения сходимости можно добиться, применяя метод Зейделя. После произвольной последовательной нумерации узлов к первому узлу применяется формула (1) с m=0. После этого формула применяется ко второй точке. При этом для первой точки в правой части формулы используется уже вычисленное значение для точки предыдущей. Этот простой механизм легко программируется, ускоряет сходимость, экономит ресурс памяти.

Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1

Эффективность решения задачи в сильной степени зависит от качественного выбора начального приближения. В качестве начального приближения приняты следующие значения потенциалов фук0:

• в ячейках открытой поверхности фук°=°;

• в ячейках, через которые проходит поверхность с граничными условиями вида (4) фук°=°;

• , а в ячейках, через которые проходят скважины, фук° определяются по формуле 5;

• во всех внутренних ячейках области рассчитываемого распределения потенциалов фук° =0.

После численного решения вышеописанной задачи Дирихле компоненты скорости в узлах элементарных ячеек Ух(хукУчк,7|к)= Ушк, Уу(хук,Уук,7ук)= Ууик ,У7(хук,Уук,7ук)= у^к можно получить, применяя правила численного дифференцирования [4] к формуле:

(ф+1]к - ф<-1 ]к )

vyjk = -

h

(ф у+1k j 1- k

h

((Pijk+1 - ф ijk -1)

h '

(7)

где h - расстояние между узлами кубической сетки узлов.

По О.Е. Власову критерий разрушения имеет вид:

а < Vkp/л/£/3 ,

где а - размер области, которая должна сохраниться целой, D - коэффициент дроби-мости, который равен:

D =

+

гдф Л 2 чдУ2 у

+

^д 2уЛ 2 vdz 2 у

+ 2

^д 2уЛ 2 дхду

+ 2

^д2 уЛ 2 дyдz

+ 2

^д2 уЛ 2 дzдx

(8)

Критическая скорость Укр вычисляется по формуле:

Укр _ °кр /л[рСе

где Окр - предельное напряжение на сжатие, отрыв или сдвиг для материала, заполняющего ячейку, рс и Е - соответственно плотность среды и и модуль упругости среды в этой ячейке.

Для вычисления коэффициента дробимости (8) требуется приближенное вычисление вторых производных фхх,фуу,ф77,фху,фу7,ф7х. Для этого используются следующие приближения в узлах[1°]:

г ] _ Ф+1 У - 2Фик + Ф-1 У

Гфхх лук 7 2 5

к

г Л _ фу+1к - 2фук + фу-1к

[фуу ]ук _ к 2 ’

\фху ] ijk =

\<Pzz ] ijk = ф+1 j+1k - ф

ф ijk+1 - 2ф ijk + ф ijk -1

h2

i+1 j-1k — ф-1 j+1k

-1 j-1k

4h2

[ ] фу+1k+1 фу+1k-1 фу-1k+1 + фу-1k-1

\ фyz ] j = 4h^

Г Фхг ]ук ли2

ф1+1 jk+1 ф1+1 jk-1 фI-1 jk+1 + фI-1 jk-1 ^к2

На этапе дезинтеграции формула для вычисления коэффициентов разрыхления имеет вид:

V V2 ,дкх д^ дvz

Kp _ — _ 1 + тp (—х у + —-),

p V p дx дy дг

где У1, У2 - соответственно значения элементарных объемов до и после разрыхления, Ух,Уу, Vz - проекции вектора скорости на координатные оси, Тр - время дезинтеграции.

Для нахождения коэффициента разрыхления применяется формула численного дифференцирования[1о]:

1 / Vi+1,у,к Vi-1,у ,к Vi, у+1,к Vi,у-1,к Vi,у,к+1 Vi, у,к-1 \

Крук _ 1 + Тр (----- ---------------------------------------------------+- -+- -) .

у,к Х'-1, У ,к уг, У+1,к У1', у-1,к г1, у ,к+1 гг, у ,к-1

На третьем этапе происходит перемещение разрушенной среды по баллистической траектории. В проекциях на оси декартовой системы координат уравнения баллистики будут иметь вид [11]:

= -Кт

а х

■ _ -Ыу | V | , (9)

_ ~Ь^ 1 V 1 - £

ш

где У - модуль вектора скорости, Ух _ йх / а; Уу _ аУ / аг; / а - проекции

вектора скорости на оси декартовой системы координат, § - значение ускорения свободного падения , 1 - время; Ь - коэффициент сопротивления.

Начальные условия для решения системы уравнений баллистики принимают вид х(°)=х°; у(°)=у°; z(o)=Zo; Ух(°)=Ух°; Уу(°)=Уу°; Vz(o)=Vzo. Значения Ух°, Уу°^° для каждой ячейки определяют при решении второго этапа задачи.

Коэффициент Ь находится по формуле[11,12]:

Ь _-------—-----.

ХсрРс (Х» У» г)

Средний размер кусков раздробленного массива хср приближенно определяется выражением хср= (dскв)0,95 [13] (¿скв - диаметр скважины).

Для численного решения системы уравнений применяется алгоритм метода Рунге-Кутта общего вида [4]. Введем в соответствии с (9) следующие обозначения:

гк _ г0 + Ш,

vxk _ vx (гк), vyk _ vy (гк), ^к _ V (гк),

Г2 , 2 , 2~~

^ _ V vл + vyk + Vzk ,

а _ Аг,

Ь1 _ -Ь^ЛАг

С1 _ (-Чл - £)Аг

а2 _ -Ь(vxk + а2\1 ^хк + ^ + (vyk + Ь) + (^к + ^2 )^ :

Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1

Ь2 _ _Ь(^к + Ь^^! (ухк + “21) + (vyk + Ь^) + (vzk + ^2 )Аг ,

С2 _ ~Ь(Угк + 2)(]1 ^хк + а1) + (Vyk + Ь1) + (Угк + ^ - £)А , аз _ -Ь(vxk + аг\1 (vxk + аг) + (^к + ь2) + (vZk + Сг)Аг,

Ь3 _ -Ь(vyk + Ь^х](Ухк + ^ + (vyk + Ь^ + (Лк + ^ А ,

С3 _ _Ь(^к + С2)(} ^хк + ^ + (VУk + + (Угк + С^) - £)Аг ,

а4 _ -Ь(vxk+азз^ (vxk+^+(vyk+ь23)+(^+^)Аг,

Ь4 _ -Ь(^к + Ь3^(Vxk + ^ + (Vyk + Ьз) + (Лк + ^А ,

С4 _ -Ь(vzk + С3^(Vxk + ^ + (VУk + Ь3) + (^к + С^) - £)А '

Тогда вычислительные формулы алгоритма будут иметь следующий вид:

а, + 2а2 + 2а3 + а4 .

Vxk+1 _ Vxk + -1--3---------4 Аг,

6

Ь + 2Ь9 + 2Ьз + Ь4 .

Vyk+1 _ Vyk + Ь1-2~6Г-Ь3-4 Аг,

с + 2с + 2сз + С4 .

vzk+1 _ ^к + ---------—3-4 Аг, где к = 0,1,2,3,-~.

6

Перераспределение объемов между ячейками на баллистической стадии происходит в соответствии с формулой[10]:

А+АА 1 (А-Ао - № )2

Р(А, г) _ [ -^= е М аА, (10)

А-АА^лМ

где Р(А, 1) - вероятность перехода среды из положения А° в положение Л за время 1; А={х,у^}; ?= {Ух,Уу^}-проекции скорости в декартовой системе координат; N - коэффициент диффузии; ДА - половина размера ячейки.

Если в формуле (10) сделать замену переменных:

т _ (А - Ао ~ Иг)

_ 4Мг ’

Тогда формула (1°) примет вид:

аи _ ^ а а .

л/М

их

■ ,-и2

Р(А, г) _^ I е -и аи.

лП и

ио

Полученный интеграл приводится к функции ошибки[4], который является табличным интегралом.

Выбранные численные методы компьютерного моделирования адекватно описывают математическую модель развала, являются простыми в реализации и широко используются на практике.

Реализация алгоритмов и численных методов математической модели развала осуществляется на базе комплекса ГИС ГЕОМИКС[14]. Это позволяет использовать исходную геологическую, маркшейдерскую и буровзрывную информацию по взрываемому блоку путём импорта из соответствующих модулей, а также использовать инструмент для работы с блочными моделями и поверхностями.

В качестве исходных данных для проведения численных экспериментов использовались результаты промышленных взрывов по 12 буровзрывным блокам на Лебединском карьере. Машинное время компьютерного моделирования развала составило от 45 минут (125 буровзрывных скважин, 3 500 000 ячеек блочной модели) до 2 часов 10 минут (520 буровзрывных скважин, 14 300 000 ячеек блочной модели).

В настоящее время ведётся работа по оптимизации программного кода с применением методов параллельного программирования на многопроцессорных серверах и графических станциях.

1. Боровиков, В.А. Моделирование действия взрыва при разрушении горных пород [Текст] / В.А. Боровиков, И.Ф. Ванягин. - М.: Недра, 1990, - 231 с.

2. Кузнецов, В.М. Математические модели взрывного дела [Текст] / В.М. Кузнецов. Изд-во “Наука”, Новосибирское отд., Новосибирск, 1977. - 762 с.

3. Друкованый, М.Ф. Управление действием взрыва скважинных зарядов на карьерах [Текст] / М.Ф. Друкованый, В.С.Куц, В.И.Ильин. - М.: Недра,198о, - 223 с.

4. Корн, Г. Т. Справочник по математике для научных работников и инженеров [Текст] / Г. Корн, Т.Корн. - М.: Наука, 1974, - 832 с.

5. Самарский А. А. Теория разностных схем. Учебное пособие. М.: Наука. 1983, - 616 с.

6. Годунов, С.К. Разностные схемы [Текст] / Годунов С.К., Рябенький В.С. М.: Наука, 1977. - 439 с.

7. Стренч Г. Теория методов конечных элементов. /Г. Стренч - М.: Мир, • 1977 - 325 с.

8. Алексидзе М.А. Фундаментальные функции в приближенных решениях граничных задач [Текст] : научное издание / Алексидзе М.А. - М. : Наука, 1991. - 352с.

9. Зельдович, Я.Б. Элементы математической физики [Текст] / Зельдович, Я.Б.

A.Д. Мышкис. - М.: Наука,1973, - 351 с.

10. Тихонов, В.И. Марковские процессы [Текст] / В.И. Тихонов, М.А.Миронов. - М.: Радио, 1977, - 488 с.

11. Черниговский, А. А. Метод плоских систем зарядов в горном деле и строительстве. [Текст] / А.А. Черниговский. - М.: Недра, 1977. - 244 с.

12. Гальянов, А.В. Трансформация структуры горных массивов при взрывных работах на карьерах [Текст] / А.В. Гальянов, В.Н. Рождественский, Ф.Н. Блинов. Екатеринбург, 1999. - 140 с.

13. Мосинец, В.Н. Дробящее и сейсмическое действие взрыва в горных породах [Текст] / В.Н. Мосинец. - М.: Недра, 1976. - 271 с.

14. Серый, С.С. ГИС ГЕОМИКС - интегрированная информационная система геологомаркшейдерского обеспечения открытых и подземных горных работ [Текст] / С.С. Серый,

B.А. Дунаев, А.В.Герасимов.// Сборник докладов Международного семинара "Передовые технологии проектирования буровзрывных работ на карьерах". Бишкек, 2006, с.87-89.

Литература

Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1

ALGORITHMS AND NUMERICAL METHODS FOR COMPUTER SIMULATION OF THE COLLAPSE OF DRILLING AND BLASTING UNIT AND DISTRIBUTION OF MINERAL COMPONENTS IN THE MINING EXPLODE AMIS

S.G. KABELKO

FSUE VIOGEM e-mail: kabelko@mail.ru

The main way to destroy rock massif in the open mining is blasting. Implementation of mathematical models of the collapse of explosive unit labor intensive process requiring large quantities of computer time. The algorithm and numerical methods for computer simulation of the explosion that adequately describe the mathematical model, are simple in implementation, including the use of methods of parallel programming for multiprocessor servers and workstations.

Key words: algorithms for modeling the explosion, numerical methods for modeling the explosion, the distribution of useful component in the collapse.