Пакет прикладных программ для анализа напряженно-деформированного состояния многоуровневой шахты с целиками с учетом воздействия гравитационных и тектонических сил Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 11, № 3, 2006

ПАКЕТ ПРИКЛАДНЫХ ПРОГРАММ ДЛЯ АНАЛИЗА НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ МНОГОУРОВНЕВОЙ ШАХТЫ С ЦЕЛИКАМИ С УЧЕТОМ ВОЗДЕЙСТВИЯ ГРАВИТАЦИОННЫХ И ТЕКТОНИЧЕСКИХ СИЛ

М.М. МУЗДАКБАЕВ НИИ математики и механики при КазНУ им. аль-Фараби,

Алматы, Казахстан e-mail: danaev@kazsu.kz

An algorithm for analysis of tensional-deformed conditions of the complicated underground construction accounting for the gravitational and tectonic forces is developed. The package of applied programs based on this algorithm is implemented.

1. Постановка задачи

Рассматривается известная постановка задачи теории упругости, где вектор перемещения и и тензор напряжений а^ удовлетворяют уравнениям равновесия Коши [1]:

divjaij} + pF = 0

и закону Гука

Gij = Xeôij +

г]-,

-г]

U ди. + ди.

2 \ дх,

дхг

i,j = 1, 2, 3.

(1)

(2) (3)

Здесь Г — вектор объемных сил; е^ — компоненты тензора деформаций; А Е

Ev

(1 + V)(1 - 2v)

и ^ = 2(1+—) — параметры Ламе; Е и V — модуль Юнга и коэффициент Пуассона; р —

плотность среды; е = ец = ец + е22 + езз, г = 1, 2, 3; 8^ — символ Кронекера. Для системы уравнений (1)-(3) в рассматриваемой области П ставим следующие краевые условия на границе Б:

1) задача Дирихле

ul

g (х)

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

S

2) задача Неймана

nj | £> ¿г •

(5)

Для единственности решения задачи (1)-(5) необходимо, чтобы компоненты тензора деформаций (е) удовлетворяли уравнениям совместности.

В случае плоской деформации уравнения равновесия (1) имеют вид

Л. д2и д2и . д2и ^

<2" + А) + + <А + " дХдУ + * = 0'

д 2и д 2и д 2и

<А + " аХду + + <2" + А) ду? + ^ = 0

Пусть имеет место вторая краевая задача

(6)

, ди л ди

<А + 2") дХ+

ди ди

сов <П,Х) + " ( дХ + ду ^ 008 <П'У) =

/ ди ди \ , .

"(ду + 8^) 008 <п'х) +

ди ди

АдХ + <а + 2") ву.

008 <П'У) = -/2'

(7)

где Х'У € £, п — внутренняя нормаль к границе £.

Как известно из литературы [2], если решение задачи (6), (7) существует, то оно обеспечивает минимум следующему функционалу:

I <и) = W <и) + / (Ди + ¿Ц ¿П - [ </1и + Ди) ¿в'

(8)

где

W <й) = ^ ^2"

ди \ / ди дх / \ ду

/ ди ди

ди ди \'

+ А ^ дх + д^ + " ^ дх + д^ ^ (9)

есть энергия упругой деформации тела.

Будем называть решением задачи (6), (7) вектор и = ^и), и, и € W21(П), который минимизирует функционал (8). Приближенным решением задачи (6), (7) называется вектор и €' <М''у) € V, который минимизирует функционал (8) на подпространстве V пространства W21(П).

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

2

2

2. Алгоритм вычисления

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

алгоритму метода конечных элементов (МКЭ) составляем систему уравнений равновесия. Базисные функции для четырехугольника сирендипового типа имеют вид [3, 4]

к = -0.25(1 - 0(1 - п)(£ + П + 1), кз = 0.25(1+ £)(1 - п)(£ - п - 1), к5 = 0.25(1+ £)(1 + п)(£ + п - 1), к7 = -0.25(1 - £)(1 + п)(£ - П + 1),

к2 = 0.5(1 - £2)(1 - п),

4 = 0.5(1 - £2)(1 + п), 2

к

ко = 0.5(1 - £2)(1 + п), к8 = 0.5(1 - £)(1 - п2).

(10)

Нумерация узлов таких конечных элементов начинается с левого нижнего углового узла против часовой стрелки. Итак, для восьми узлов имеем восемь соотношений по Зенкевичу [3].

Каждое из этих соотношений равно единице только в тех индексах, которые соответствуют номерам точек, и равно нулю в остальных узлах, и координату х, г любой точки можем аппроксимировать с помощью декартовых координат вершин в виде

х

хг

г=1

8

У] кггг.

г=1

(11)

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

дкг дкг дх дкг дг дкг ' ' + ' '

дкг дх дкг дг

д£ дх д£ дг д£ дп дх дп дг дп

или

дкг

Ж

дкг

дп

Отсюда матрица Якоби имеет вид

дхдг

д£д£ дх дг

дп дп

дкг

дх

дкг

дг

(12)

(13)

V ]

дх дг дх дг

дп дп

Обратный переход осуществляется обращением (13):

дкг ^ ( дкг

дх

дкг

дг

[V ]

-1

д£ дкг

дп

(14)

(15)

После всевозможных подстановок для перехода в систему координат £Оп формула (15) принимает вид

дкг

дх

дкг

дг

[V]

8 дк - 8 дк1

/ .. о гг п а гг

г=1 дп г=1 д£

8 дЫ 8 дк1

г=1 дп г=1 д£

дкг дкг

дп

г

1

После окончательного перехода в систему £0п, который подробно изложен в [8], получим формулу для построения матрицы жесткости произвольного четырехугольника с учетом внутренних точек интегрирования:

[К]е = £ ^ ау [В]т (17)

i,j

где aij — гауссовы точки интегрирования [3, 4], а ^ — толщина элемента в точке ¿']; [В], [Д — матрицы базисных функций и упругости элемента е; det [У] — якобиан. Матрица жесткости системы вычисляется суммированием матриц жесткости отдельных элементов, рассчитанных по (17). Координаты пО, ^^ пО и перемещения пг) запишем в виде

(ж) = Шж„}' (*) =[д](гс)' (и) = [д](ЦЬ)' (18)

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

(Р ) = ДОГ (/)* det[J]' (19)

где (/ — вектор объемной силы для элемента. Горизонтальная составляющая тектонической силы имеет вид [7]

(Т) = Х(Р)' (20)

где X — коэффициент тектонического сжатия; V < х < 1, V — коэффициент Пуассона. Уравнения равновесия всей системы запишем в виде

[К](и„) = (Р) + (Т). (21)

Векторы деформаций и напряжений вычисляются в точках , :

(ек}„• = [Вк,т]i>j■ (и„т)' к = 1' 2' 3' т = 1' 2...8' г = 3 = 1' 2'...' 9; (22)

(°к)i>j■ = [Ок,и] (ек)' п =1' 2' 3'к =1' 2' 3' (23)

где (ек)i)j = (ек,Х'ек,г 'Ак,хг )i)j, (ak)ij = (^к , х' 0к, г' Тк,х:х

)ij — деформации и напряжения;

и 0

(Щ

0 и

I = 1' 2'...' 8, — перемещения узловых точек.

3. Расчетная область задачи

Как известно, по мере разработки верхних слоев полезных ископаемых приходится добывать руду из глубоких нижних пластов. В целях предотвращения обрушения кровли выработки на определенных расстояниях друг от друга оставляются целики из породы, а пустые отработанные забои закладываются пустыми породами или быстро твердеющими растворами. На рис. 1 показана схема расчетной области наклонно-слоистого массива вокруг шахты с многоуровневыми выработками. В данном случае рассматриваются три выработки с оставленными целиками и уровнями залегания выработок Н1 = 250, Н2 = 310 и = 400 м. В каждой расчетной точке области действуют геостатическая сила, направленная вниз, и горизонтальная составляющая тектонической силы, направленная в сторону выработки либо от нее. Расчетные и измеренные данные тектонической силы показывают, что существенное значение имеет ее горизонтальная составляющая, которая по

с

Рис. 1. Схема расположения выработок трехуровневой шахты: L — общая длина расчетной области ABCD; Li — L3 — длины выработок, соответственно верхней, средней и нижней; H1 — H3 — их глубины залегания; hB — высота выработок.

величине значительно превосходит силы тяжести с удалением в глубь массива от земной поверхности [5]. Размеры расчетной области для МКЭ выбираются из условия соблюдения гипотезы Динника. Численные эксперименты показали, что для этого ширину и высоту внешней границы расчетной области независимо от типа неоднородного строения массива

необходимо выбирать в соотношении [6]

Ь = 8Я' 2.2 < 8 < 2.3. (24)

4. Исходные данные

Длина исследуемой области АБСБ Ь = 2000 м, высота Я = 800 м. Выработки находятся на глубинах Я1 = 250, Я2 = 310, Я3 = 400 м и все имеют высоту, равную к = 4 м. Длины выработок Ь1 = 160, Ь2 = 250 и Ь3 = 300 м, причем 1/5 правой части третьей нижней выработки имеет наклон под углом ^ = 15

Физико-механические характеристики анизотропного транстропного массива, состоящего из горного массива, пласта угля, закладочного материала и геологического разлома, взяты из работы [1]: Е1 = 2' 0 ■ 104 МПа, Е2 = 1' 5 ■ 103 МПа, Е3 = 3 ■ 102 МПа, Е4 = 2 ■ 103 МПа, v1 = 0' 28, v2 = 0' 3, ^ = 0' 32, ^ = 0' 26, 71 = 2' 46 т/м3, 72 = 2.4 т/м3, 73 = 2' 0 т/м3, 74 = 2' 3 т/м3.

Сформулированы следующие граничные условия: поверхность земли СД свободна от напряжений ох = тгх = 0 на ВС, а на АД и = 0, и = 0.

Геостатические напряжения в массиве вычисляются по известным формулам

ах = -7Я' ах = ^' П = v/(1 - V). (25)

Область численного решения задачи показана на рис. 1. Рассматривается шахта с трехуровневой выработкой, учитываются массовые и тектонические силы. В каждой из выработок оставлены целики: в первой — два, во второй — три, в третьей — четыре, как показано на рис. 1. Для удобства пронумеруем целики слева направо, например, в первой выработке 1 и 2; во второй — 1, 2 и 3 и т. д.

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

Рассматриваемая область разбивается на 430 изопараметрических четырехугольных элементов с общим количеством узлов 550. Вокруг выработок для ожидаемой концентрации напряжений сетка конечных элементов сгущалась.

Система основных уравнений МКЭ (21) решалась итерационным методом Гаусса — Зейделя.

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

На стороне а1с1 имеют место сжимающие тангенциальные напряжения ав, наибольшие значения которых приходятся над целиками 1 и 2, причем в случае воздействия геостатических и тектонических (+) сил напряжения ав над вторым целиком на 10-30% больше, чем над первым, тогда как в случае тектонических (—) они почти равны. По мере удаления от первого целика в сторону угловой точки а1 напряжения а в резко уменьшаются в десятки раз и становятся близкими к нулю, а в угловой точке меняют знак и становятся растягивающими. На левом торце а1к1 первой выработки имеются слабые растягивающие напряжения ах. Близкие к нулю растягивающие напряжения а в появляются и на стороне

Рис. 2. Эпюры тангенциальных напряжений вдоль контуров выработок, целиков и стволов шахт: 1 — тонкие сплошные линии соответствуют геостатическому (—7Н) случаю, 2 — тонкие штриховые линии — воздействию тектонических (+) сил (в направлении оси Ох), 3 — жирные штрих-пунктирные линии — воздействию тектонических (—) сил в направлении, противоположном оси Ох.

к\Г\. Сжимающие напряжения о® на стороне g\f\ растут на контакте с первым и вторым целиками, их значения в случае геостатических и тектонических (+) сил увеличиваются до значения, равного 50-60 % их верхнего максимального значения на контакте с целиками, тогда как при тектонических (—) силах нижние растягивающие напряжения о® почти равны верхним максимальным около целиков. Такое напряженное состояние в массиве имеет место в случае геостатического воздействия. Сжимающие напряжения на сторонах с2а2, а2к2, к2г2, достигают максимального значения на контакте с первым целиком второй выработки, причем его нижнее значение в 2.7 раза превосходит верхнее.

Тангенциальное напряжение на стороне ¿2е2 ведет себя очень сложным образом, при-

нимая растягивающее значение на контактах со вторым и третьим целиками, а также между ними, принимая затем сжимающие значения. Следует отметить, что все эти значения малы.

Такая картина наблюдается вдоль контура третьей выработки. Стволы шахты находятся в слабосжатом либо слаборастянутом состоянии. Если рассматривать напряженное состояние целиков, то у двух в первой выработке, у первого во второй и первого в третьей выработках они находятся примерно в одинаковых условиях, т. е. на контакте с контурами выработок они испытывают сильное сжатие, а в срединной части — растяжение, но значение последнего меньше сжатия. Исключение составляет четвертый целик в третьей выработке: в случае воздействия тектонических (+) и геостатических сил он находится под воздействием сил растяжения, тогда как в случае воздействия тектонических (—) сил он находится в слаборастянутом состоянии. Также сложным образом ведут себя второй и третий целики второй выработки, например, второй в случае воздействия тектонических (+) и геостатических сил находится в слабосжатом и слаборастянутом состоянии в срединной части, тогда как третий в случае тектонических (—) сил снизу на контакте с контуром испытывает сжатие. Следует отметить также появление зон слаборастягиваю-щих напряжений ав на нижней части контура третьей выработки перед наклоном вниз на <15 о и близкие к нулю значения а в напротив верхней части контура в случае воздействия геостатических и тектонических (+) сил.

В заключение отметим, что величины тангенциальных напряжений почти во всей области в случае воздействия тектонических (+) сил наибольшие, а в случае геостатических сил напряжения становятся меньше. Самые наименьшие напряжения возникают при воздействии тектонических (—) сил. Их сравнение со случаями воздействия геостатических и тектонических (+) сил показывает, что они близки друг к другу.

Точность проведенных расчетов сравнивалась с результатами, полученными в [8]. Установлено хорошее совпадение с ними, т. е. по мере удаления от выработок в сплошном массиве везде и всюду выполняются условия равенства геостатистического давления ах = 7Я, где ах получен из результатов численных расчетов. Поскольку для трехуровневой шахты в литературе отсутствуют какие-либо оценки точности с помощью прямых или косвенных замеров напряженных состояний в забоях, так как ранее такие задачи в комплексе не рассматривались, в предельном переходе за основу приняты оценки точности по результатам счета для одиночной выработки. Проведено сравнение аналитических решений этой задачи с численными результатами, что показало их точное совпадение [8].

Следует сказать, что модель, предложенная Ж.С. Ержановым, Ш.М. Айталиевым и Ж.К. Масановым для расчета напряженно-деформированного состояния подземных сооружений в наклонно-слоистом анизотропном массиве, применима для анализа напряженно-деформированного состояния целого класса задач, как утверждалось многочисленными авторами, и нашла еще одно подтверждение в приведенных расчетах.

Таким образом, установлены основные закономерности распределения концентрации напряжений на контурах сложной системы многоуровневых подземных сооружений.

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

[1] Тимошенко С.П., Гудвер Дж. Теория упругости. М.: Наука, 1975.

[2] Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970.

[3] Зенкевич О. Метод конечных элементов в технике. М.: Мир, 1975. 541 с.

[4] Сегерлинд Л. Применение метода конечных элементов. М.: Мир, 1979.

[5] Ержанов Ж.С., Айталиев Ш.М., Масанов Ж.К. Сейсмонапряженное состояние подземных сооружений в анизотропном слоистом массиве. Алма-Ата: Наука, 1980. 211 с.

[6] Булычев Н.С. Механика подземных сооружений. М., 1986.

[7] Тектонические напряжения в земной коре и устойчивость горных выработок / И.А. Турчанинов, Г.А. Марков и др. Л.: Наука, 1978. 256 с.

[8] Баймахан Р.Б. Расчет сейсмонапряженного состояния подземных сооружений в неоднородной толще методом конечных элементов. Алматы: Дауир, 2002. 230 с.

Поступила в редакцию 9 декабря 2005 г., в переработанном виде — 20 февраля 2006 г.