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

Анализ напряженно-деформированного состояния в окрестности выработки в массиве горных пород

Г.Е. Уцын, А.А. Беспалько1, Б.А. Люкшин, Н.Ю. Матолыгина

Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия 1 Томский политехнический университет, Томск, 634050, Россия

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

Analysis of the stress-strain state in the vicinity of a mine opening in the rock massif

G.E. Utsyn, A.A. Bespalko1, B.A. Lyukshin, and N.Yu. Matolygina

Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia 1Tomsk Polytechnic University, Tomsk, 634050, Russia

The paper is devoted to mathematical simulation and analysis of the static stress-strain state in the vicinity of a mine opening in the natural state with account for overlaying mass loading or mass forces near the opening. We also discuss problems of simulating wave processes occurring near the opening under pulse loading on a part of its inner surface.

1. Введение

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

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

2. Физическая постановка задачи

Задача анализа напряженно-деформированного состояния в окрестности выработки в массиве горных пород в общем случае является трехмерной в силу отсутствия точечной, осевой или плоской симметрии. В то же время, известно [2, 3], что анализ распространения волн в двух измерениях дает возможность обнаружить те же эффекты, что и в случае трехмерной постановки задачи. В частности, это относится к наличию двух типов волн (продольных и поперечных), распространяющихся с одинаковыми характерными скоростями в двумерном и трехмерном случаях. С физической точки зрения это связано с наличием двух упругих характеристик однородной упругой изотропной среды.

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

© Уцын Г.Е., Беспалько А.А., Люкшин Б.А., Матолыгина Н.Ю., 2005

Анализируемая область представляет собой прямоугольник ABCD, линия АВ соответствует поверхности земли. На границах AD, ВС необходимо задавать условия, известные как условия скольжения вдоль жесткой стенки, на линии АВ — условия на свободной поверхности, и на линии CD можно принять условия опирания на жесткую недеформируемую поверхность.

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

В связи с этим при решении задачи определения напряженно-деформированного состояния в статической постановке предлагается подвергать анализу не всю область ABCD, а часть ее, показанную на рис. 1 в виде прямоугольника abcd. В этом случае в дополнение к массовым силам, приложенным к каждому из элементов массива, необходимо учесть нагрузку от «отброшенной» верхней части. Эта часть нагрузки показана на рис. 1 в виде эпюры равномерного давления Р. Такое давление задается как одно из граничных условий на линии аЬ. В качестве второго условия можно использовать либо условие отсутствия касательных напряжений, либо условие отсутствия перемещений вдоль этой линии. На линиях ad и Ьс по-прежнему задаются условия сколь-

1

Я

\У \|/ \|/ \У У\У\|/\|/ФуУ>УУ

і

жения вдоль жесткой стенки, а на линии cd — условия опирания на жесткую недеформируемую поверхность.

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

3. Математическая постановка и расчетные методы

Напряженно-деформированное состояние области abcd (рис. 1) под действием статической нагрузки исследуется с помощью метода конечных элементов в сочетании с процедурой последовательных нагружений. В соответствии с этой процедурой анализ напряженно-деформированного состояния расчетной области проводится как многошаговый процесс, каждый из шагов настолько мал, что допускает линеаризацию всех соотношений — геометрических, физических и статических — при переходе от текущего состояния к смежному [2].

Формулировку статической задачи можно записать в виде принципа виртуальной работы для конечных приращений. Этот принцип лежит в основе конечно-элементной формулировки задачи:

ст9. АикЪАик ,г + ст Аик,г 8Ди,.} )<Т =

У///////////////////////////

Рис. 1. Расчетная схема

= I АТ, ЗДи,^ -1 ст® 8Ди,-} +| Г/^Ди^. (1)

3 V 3

Здесь ст у — компоненты тензора напряжений; индексО относится к исходному (текущему) состоянию, соответствующие величины без индексов отвечают смежному состоянию; Т, — вектор внешних нагрузок; V — объем (в двумерном случае — площадь) расчетной области; 5—площадь поверхности (длина контура), ограничивающего расчетную область; Аык — приращения перемещений; 8Д% — виртуальные изменения приращений перемещений; С,]к[ — коэффициенты, связывающие между собой приращения напряжений Дстгу и приращения деформаций Дву:

= СуЫ АвЫ . (2)

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

Р

пользуются уже при вариационной формулировке задачи («заложены» в функционал).

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

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

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

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

В качестве начальных условий принимается отсутствие начальных смещений и их скоростей, напряжений и деформаций внутри расчетной области.

На линии у = 0, совпадающей с линией контура выработки, приложена нормальная нагрузка, которая задается в виде

р(х, о = х(х) щ. (3)

Первая из функций в правой части определяет закон изменения нагрузки по пространственной переменной, вторая — во времени.

В расчетах принимался закон изменения нагрузки вдоль оси х в виде так называемой колоколообразной функции [4]. В плоской задаче эта функция принимается в виде

Т (?) = в ^ бш—

где Т/2 определяет длительность действия импульса.

На остальных границах расчетной области нужно ставить неотражающие граничные условия.

Варианты распространенных постановок граничных условий на контуре расчетной области отвечают двум крайним случаям: 1) свободная от напряжений поверхность; 2) опирание на жесткую стенку.

В первом случае нормальная и касательная компоненты вектора напряжений обращаются в нуль:

ст = 0. (5)

Во втором случае требуется равенство нулю компонент вектора перемещений или вектора скорости перемещений:

2Ш

( X-х0)

X (х) = е 2й?

(4)

Величина d определяет форму «колокола». Чем больше значение d, тем медленнее меняется нагрузка по мере удаления от центра ее приложения.

Закон изменения нагрузки во времени принимается в виде:

^ = 0, дг

(6)

где ст — вектор напряжении; и — вектор перемещении (или скоростеИ перемещении) на границе.

С физическоИ точки зрения на искусственно вводимой границе не реализуется ни один из этих случаев.

В работе [5] для волн Рэлея (плоская задача) предлагаются и обсуждаются так называемые неотражающие классические граничные условия вида

. дм

а-п + М— = 0,

Ы

(7)

где

а = (ст„

( Л

}22 >

М =

Л2 - П

Т12),

0

п1

П2

Л2 - П

0 V,

ч / ч / ч /

В этих уравнениях п(пь п2) — направляющий вектор внешней нормали к границе; пь п2 — косинусы углов направляющего вектора с осями х, у в случае использования декартовой системы координат; стп, а22, ст12 — компоненты тензора напряжений; р — плотность среды; V Vs — скорости продольной и поперечной упругих волн соответственно.

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

Исследования авторов [5] показали, что такие условия порождают отражение на границе расчетной области. Интенсивность отраженной волны имеет порядок 0(02) для продольных волн и порядок 0(9) для волн сдвига, где 0 — угол падения волны, отсчитываемый от нормали к поверхности.

С формальной точки зрения условия (7) можно рассматривать как линейную комбинацию условий (5) и (6). Исследования показывают, что подбором коэффициентов в этой комбинации можно значительно снизить эффект отражения.

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

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

- алгоритм дает монотонное решение;

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

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

- алгоритм обобщается на случай неупругой среды;

- используются неотражающие граничные условия на границах расчетной области.

Численный анализ нестационарных задач проводится с использованием нецентральной разностной схемы [6]. Эта схема относится к схемам сквозного счета. Для моделирования волновых процессов привлекается метод, который имеет второй порядок точности относительно шагов по пространственным переменным и по времени.

4. Результаты расчетов

При решении плоской задачи в каждой ячейке конечно-элементной или конечно-разностной сетки определяются 8 величин: две компоненты вектора перемещений и по три компоненты тензоров напряжений и деформаций. Наиболее информативным представляется распределение по расчетной области интенсивности деформаций.

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

Распределение интенсивности деформаций при отсутствии каких-либо включений (прожилков) (рис. 2, а) симметрично относительно вертикальной срединной оси расчетной области. Учет горизонтально ориентированного прожилка справа от проходки эту симметрию нарушает, что отражается в форме поверхности, иллюстрирующей распределение интенсивности деформаций (рис. 2, д), и в изменении изолиний (рис. 2, б). При проведении данного расчета принималось, что модуль упругости материала включения составляет 2/3 от соответствующей характеристики основной породы, для последней принималось значение модуля упругости Е = 90 ГПа. Введение вертикально расположенного про-

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

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

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

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

Анализу подвергалась расчетная область прямоугольной формы (рис. 3), размер этой области 0.05x0.10 м2 соответствовал размерам 0.05x0.05x0.10 м3 экспериментально исследуемых лабораторных образцов в виде параллелепипедов. На границе прямоугольной расчетной области у = 0 задается возмущение в виде (3) и анализируется влияние присутствующих в области неоднородностей на распространение этого возмущения. Основная порода имитируется цементным камнем (модуль упругости 21.6 ГПа, коэффициент Пуассона 0.25, плотность 1.78 -103 кг/м3), компактное включение — медная вставка (соответствующие характеристики — 114 ГПа, 0.34, 8.8-103 кг/м3).

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

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

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

Рис. 3. Форма расчетной области при наличии в ней компактного включения: включение (1); основная порода (2); слева показана форма приложенного импульса как функция координаты х

большей скоростью, что приводит к наличию изломов на изолиниях распределения интенсивности деформаций на оси (рис. 4, в).

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

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

Єі

0.00010

0.00005

0.00000

Сі

0.0004

0.0002

0.0000

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

Рис. 5. Поверхности (а, б) и изолинии (в) интенсивности деформаций при наличии жесткого компактного включения

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

5. Выводы

СтатическиИ анализ позволяет оценить параметры напряженно-деформированного состояния в окрестности горноИ выработки, в том числе и при наличии неоднородностей. Эта информация является исходноИ и должна использоваться при задании начальных условиИ в формулировке динамическоИ задачи.

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

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

Работа выполнена при финансовой поддержке РФФИ, грант № 05-01-98005.

Литература

1. Беспалько А. А., Гольд P.M., Яворович Л. В., Дацко Д. И. Возбуждение

электромагнитной эмиссии в слоистых горных породах при акустическом возбуждении // ФТПРПИ. - 2003. - № 2. - С. 8-14.

2. Люкшин Б.А., Герасимов А.В., Кректулева Р.А., Люкшин П.А. Моделирование физико-механических процессов в неоднородных конструкциях. - Новосибирск: Изд-во СО РАН, 2001. - 272 с.

3. Иванов Г.В., Волчков Ю.М., Богульский И.О. и др. Численное решение динамических задач упругопластического деформирования твердых тел. - Новосибирск: Сиб. унив. изд-во, 2002. - 352 с.

4. Немирович-Данченко М.М. Численное моделирование трехмерных

динамических задач сейсмологии // Физ. мезомех. - 2002. - Т. 5.-№ 5. - С. 99-106.

5. Bamberger A., Chalinder B., Joly P., Roberts J.E., Teron J.L. Absorbing

boundary conditions for Rayleigh waves // SIAM, J. Scientific and Statistical Computing. - 1988. - V. 9. - No. 6. - P. 1016-1049.

6. Уорминг П.Ф., Кутлер П., Ломакс Г. Нецентральные разностные схемы II и III порядков точности для решения нелинейных гиперболических уравнений // Рак. техн. и косм. - 1973. - Т. 11. - № 2. -С. 76-85.