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

УДК 539.3

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

А.В. Круподеров, М.А. Журавков

Белорусский государственный университет, Беларусь, Минск

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

Проблема моделирования и прогнозирования воздействия крупномасштабных подземных горных работ (например, извлечение полезных ископаемых, подземное инженерное строительство и др.) на земную поверхность является актуальной и важной. Так, например, неправильная оценка величин и характера смещений земной поверхности может привести к значительному ущербу всей наземной инфраструктуры, находящейся в рассматриваемом регионе.

К настоящему времени существует достаточно большое количество подходов и методов к расчету перемещений (сдвижений) земной поверхности от ведения подземных горных работ [2].

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

Сначала опишем первую технологию. В ее основе лежит следующая феноменологическую модель: поскольку

любая искусственная полость (выработка) в массиве вносит возмущение в исходное НДС, то ее воздействие на земную поверхность и приповерхностные области можно моделировать воздействием некоторых усилий, действующих в массиве (рис. 1).

Рис. 1. Переход от реальной геомеханической задачи к модельной

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

Аналог решения Миндлина для трансверсально-изотропной среды можно записать в таком виде [4]:

«=Е i=1

2

«2 =Е

АХ1

V 2 (

sign(Х3 - Ъ) аах-— + £-

Ъ (Ъ + 5 I Хз - ъ I) £1 Г, (Г, + 2, ) у

sign(х3 - Ъ)

Е-

Л..х„

, 2

V

2

Ъ ,(Ъ,,+ 5 , | Х3 - Ъ I) г..(г. + 2 ,, )

,, ^ 11 ,1^1/ ,=1 у V у у /у

=1

«з=Е«, ^+ЕА

л

где

5 =

(С - С1з)(С + С13 + 2с44)

4сз3С44

+ (-1Г

(С + С1з)(С - С13 - 2С44)

' 4С33С44

= 1 2 53 , С = у]СПС33 .

Здесь, в свою очередь, (с),, ,, у = 1,6 - матрица упругих констант, связывающая напряжения с деформациями, т.е.

С1 = ^ С2 = ^22, С3 = ^33, С4 = ^23 , С5 = ^ С6 = С12 с 81 = Б11, 82 = 822,83 = 833, 8 4 = 2823, 85 = 2В13, 86 = 2В12 .

В общем случае эффективные упругие параметры массива таковы, что 5 Ф 52. Тогда, с учетом данного обстоятельства, величины, входящие в формулы (2), имеют следующий вид (частный случай 5 = 52 приведен в [4]):

, Г, =\1 х12 + х22 + 2,у,

2 , = 5 ,Х3 - Ну , 2 , = 5 ,Х3 + Н у ,

И, — 5 Ъ, а —

с13 + с44

1 4л СззС44(522 - 52)' 2

_ + v2 w1 _ v2

— а , — 2а0

А

А

а Х2

Л = 2«,^, Л, = Л,,, А

= kIc66, V,. = сзза& - c1з,

2

, ai + si с,, -с^,.

К = ~ , а =■

^ (с13 + c44)sl

Используя указанные фундаментальные решения (1) и (2), можно построить решение для нагрузки, распределенной по произвольной площади. Применительно к рассматриваемому классу прикладных задач достаточно рассмотреть нагрузку, распределенную по площадке, параллельной поверхности полупространства. В этом случае выражение для компонент перемещений, вызванных такой нагрузкой, можно записать в следующем виде:

V, (X, у, z) = Л и.((х -%, у - ц, z) f (£, ц^ц, (3)

я

где через V. обозначены компоненты перемещений, вызванные распределенной нагрузкой, а и{ - компоненты, являющиеся следствием воздействия сосредоточенной силой и определяемые (1) или (2), А(ху) - плотность нагрузки, распределенной по площади Я.

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

Исходя из механического смысла прикладной задачи, установлено, что нагрузку можно прикладывать на естественной глубине залегания выработки. При этом функция плотности нагрузки коррелирует с реальным распределением напряжений возле выработки. Предлагается следующий вид для функции Ах,у):

/(x,,x2) = d <

1, if | Xj |< width,

aXj2 + bx1 + c, if width <| x |< mp, (4)

2

J -(x, -mp)/с . r , .

md e 1 , if mp <| x |

Здесь width - ширина выработанного пространства, mp - координата расположения точки максимума опорного давления, md - отношение максимального опорного давления к давлению над выработкой. Параметры a,b,c в (4) выбираются из условия непрерывности и гладкости нагрузки в точке максимума опорного давления, а также из условия ее непрерывности на контуре выработки. Параметры md, d, с зависят от глубины и ширины выработанного пространства. Параметр mp зависит только от ширины выработки.

На рис. 2-5 отражены сравнения кривых оседаний в главном поперечном сечении мульды сдвижения с кривыми, полученными по [3] (пунктирная кривая соответствует модели) при различных значениях ширины выработки и ее ширины.

Далее приведем описание второй технологии. Поскольку реальный массив горных пород представляет собой гетерогенную среду, то для выбора модельной среды, поведение которой наиболее адекватно состоянию реального массива, рассматривалась модель трансверсально-изотропной упругой деформируемой среды.

При рассмотрении массива как трансверсально-изотропной среды (в таком массиве существует направление, в нашем случае вертикальное, поворот вокруг которого на любой угол не влечет за собой изменение упругих свойств) ее состояние определяется заданием пяти констант: Ep, v р, Epz, v pZ, GZp . Здесь Ep, v p - модуль Юнга и

коэффициент Пуассона в горизонтальной плоскости, а Ер2, v р2, Gzp - модуль Юнга, коэффициент Пуассона и модуль сдвига в вертикальном направлении.

Рис. 2. D = 265, Н = 830 Рис. 3. D = 120, Н = 680

Рис. 4. D = 120, Н = 680 Рис. 5. Отработка двух блоков

Если вертикальную ось обозначить как ось Оz, то в этом случае закон Гука запишется в следующей форме:

Г1 -л

ЕЛ А

ЕГА А

елА

ЕЛ А

1 -у р? г

ЕЛ А

ЕЕ А

ЕЕ А

ЕГЕ А

1 -ур

ЕрА 0

0 0

0 0 ^ 0 0

0 0

2Gp 0 0

0 0

0

ер

1 + у„

(5)

Е.

(1 + ув)(1 -V -2у„уи)

где у^ = ^у_., А- р ^ ^

Е

Е 2 Е

у +у у

р рг г

у +у у

гр рг г

0

Га Л

Гк Л

у +у у

р рг гр

X

0

у +у у

гр рг гр

V + у у

гр рг гр

0

0

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

Таким образом, полная система разрешающих уравнений включает в себя:

1) уравнения равновесия

^хх = 0,

дх ду

да да

—— +--- — -pg.

дх ду

2) уравнение совместности

д2в д2 в д2в и __уУ. = 2- ху

ду2 дх2 дудх '

3) закон Гука;

4) граничные условия соответствующего типа на всех частях контура области.

В работе [1] показано, что основное влияние на характеристики формирующейся мульды сдвижения для массивов горных пород с пластовым залеганием слоев оказывает модуль сдвига в вертикальном направлении. В случае отсутствия установленных в лабораторных условиях значений модуля сдвига в вертикальном направлении рекомендовано брать его величину 10-15 % от значений модуля сдвига в горизонтальном направлении.

В массиве было выделено три слоя с эффективными упругими свойствами. Учет действия обрушенных пород достигался введением так называемой эффективной мощности блока, которая зависит от ширины блока и реальной вынимаемой мощности. Схема расчета выглядит следующим образом:

- вычисление НДС ненарушенного тяжелого массива;

- удаление элементов, моделируемых выработанное пространство. Расчет НДС массива с выработкой и определение зон нарушений в соответствии с принятым критерием;

- определение алгоритма моделирования «схождения кровли и подошвы выработки, вследствие формирования области разрушения». Алгоритм состоит в склейке узлов, расстояние между которыми при деформировании не превышает некоторой заданной величины.

На рис. 6-7 отражены сравнения кривых оседаний в главном поперечном сечении мульды сдвижения с кривыми, полученными по [3] (красная кривая соответствует [3]) при различных значениях ширины выработки и ее глубины.

Рис. 6. D = 150, Н = 530

Рис. 7. D = 300, Н = 430

Список литературы

1. Гавриленко Ю.Н. Численное моделирование процессов сдвижения земной поверхности методом конечных элементов в объемной постановке / Ю.Н. Гавриленко // Фи-

зико-технические проблемы горного производства. -Вып. 3. - 2001. - С. 12-25.

2. Журавков М.А. Геомеханический мониторинг горных массивов / М.А. Журавков, О.В. Стагурова. - Минск: Юникап, 2002.

3. Указания по охране сооружений и природных объектов от вредного влияния подземных горных разработок в условиях Старобинского месторождения калийных солей / Степанов К.А., Журавков М.А., Невельсон И.С. [и др.]. -Солигорск: Белгорхимпром, 2002. - 341 с.

4. Ding H. Elasticity of transversely isotropic materials / H. Ding, W. Chen, L. Zhang. - Springer, 2006.