CC BY
52
Моделирование обрушения кровли над выработанным пространством
П.В. Макаров, И.Ю. Смолин, Е.П. Евтушенко, А.А. Трубицын1,
Н.В. Трубицына1, С.П. Ворошилов1
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия 1 Научный центр по безопасности работ в угольной промышленности ВостНИИ, Кемерово, 650002, Россия
Предложен и развит подход к численному решению прикладных задач разрушения прочных сред и твердых тел, основанный на идеях эволюции и самоорганизации нелинейных систем в процессе нагружения. Предложенная методика позволяет рассматривать катастрофические явления как развитие событий в режимах с обострением.
Нагружаемая полем сил тяжести геосреда, в том числе содержащая различные полости антропогенного происхождения (скважины, штреки и лавы в шахтах, различные подземные сооружения), рассматривается как нелинейная динамическая блочная система, эволюционирующая в поле сил по законам синергетики. Механическое поведение таких нелинейных динамических систем исследуется с континуальной макроскопической точки зрения, что позволяет эффективно использовать подходы и численные методы механики сплошных сред для моделирования процессов деформации и разрушения.
На высокопроизводительном вычислительном кластере выполнены расчеты эволюции состояния и обрушения кровли над выработанным пространством в зависимости от скорости подвигания забоя. Показано, что при высоких скоростях подвигания забоя реализуются существенно неравновесные режимы эволюции нагруженного массива, в результате чего зависают протяженные участки пород покрывающей толщи, что может привести к их катастрофическому обрушению.
Simulation of roof collapse of the worked-out area
P.V. Makarov, I.Yu. Smolin, E.P. Evtushenko, A.A. Trubitsyn1,
N.V. Trubitsyna1, and S.P. Voroshilov1
Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia 1 Research Center for Work Safety in Coal Industry VostNII, Kemerovo, 650002, Russia
The paper puts forward an approach to the numerical solution of applied problems on the failure of media and solids, which is based on the notions of evolution and self-organization of nonlinear systems under loading. The method allows treating catastrophic events as the development of events in blow-up regimes.
A geomedium under gravity, including the one with various man-made hollows (boreholes, galleries and mine faces, different underground structures), is considered as a nonlinear dynamic block system that evolves in the gravity field by the synergetic rules. We study the mechanical behavior of such nonlinear dynamic systems from a continuous macroscopic standpoint. This permits efficient application of continuum mechanics methods, approaches and numerical methods of simulating the deformation and fracture processes.
With a high-performance computing cluster we calculate the roof evolution and collapse above the worked-out area depending on the rate of face advance. It is shown that at high advance rates the evolution regimes of the loaded mass are highly nonequilibrium. As a result, extended hanging roof zones in this area are prone to catastrophic collapse.
1. Введение. Формирование неравновесных состояний горного массива над выработкой
За последние 10-15 лет производительность труда в шахтах многократно возросла. В разы увеличилась скорость подвигания забоя. С точки зрения безопасности работ стала складываться чрезвычайно тревожная ситуация. Состояние горных массивов в кровле, в приза-
бойной зоне сместилось в области существенно неравновесных состояний и размеры зависающей кровли достигают катастрофически опасных величин, часто превышающих 100 м до первой посадки (обрушения). Некоторые данные о призабойной зоне и временах первой и последующих посадок кровли приведены в табл. 1.
© Макаров П.В., Смолин И.Ю., Евтушенко Е.П., Трубицын А.А., Трубицына Н.В., Ворошилов С.П., 2008
Тенденцию заметного увеличения шагов начального и последующих обрушений кровли при увеличении скорости подвигания забоя при прочих равных условиях (глубине выработки, длине лавы, физико-механических параметрах горного массива) мы связываем с тем, что столь быстрое увеличение выработанного пространства приводит окружающие выработку породы в существенно неравновесное состояние. Это обусловлено тем, что скорость накопления опасных повреждений в разных массивах, определяемая геологией массива, глубиной выработки и т.д., «отстает» от скорости подвигания забоя. Забой «ушел» далеко вперед, а повреждения в зависающей кровле еще очень далеки от критических значений и деградация прочностных свойств геосреды незначительна.
Условно можно выделить три режима:
1. Равновесный, реализуемый при малых скоростях или слабых породах, в которых скорости накопления повреждений велики на данных глубинах, или при значительной глубине проводимых горных работ, что также
приводит к более быстрому накоплению повреждений. В этом случае скорость накопления повреждений оказывается столь быстрой, что обрушения кровли следуют сразу за продвижением забоя (табл. 1, № 1, 2).
2. Квазиравновесный режим. В этом случае заметные повреждения успевают накапливаться вместе с по-двиганием забоя, т.е. «волна» повреждений в кровле движется приблизительно со скоростью подвигания забоя, значительных и опасных зависаний кровли не наблюдается. Первая посадка кровли составляет первые десятки метров, шаги последующих обрушений кровли также невелики (табл. 1, № 3, 4).
3. Существенно неравновесный режим. Это наиболее опасная ситуация. Забой успевает уйти достаточно далеко, прежде чем в кровле накопятся заметные повреждения. Зависают очень протяженные участки кровли, шаги последующих обрушений также велики (табл. 1, № 5-7). Катастрофическое обрушение такого огромного массива могут спровоцировать различные причины.
Таблица 1
Сведения о временах посадок кровли над очистными забоями
№ п/п Шахта, пласт, марка Наиме- нование выра- ботки Скорость подвигания забоя, м/мес Расстояние, пройденное до первой посадки кровли, м Последующие шаги обрушения кровли, м Кол-во дней, пройденное до первой посадки кровли Мощность непосредственной кровли, м Мощность основной кровли, м Глубина залегания Н, м Длина лавы, м Мощность угольного пласта h, м
1 Томусинс-кая 5-6, пласт V, марка КО Лава 0-1-5 120 Сразу при выходе из монтажной камеры Обрушение вслед за по-двиганием комплекса 4-5.5 30 310-420 2.5
2 Есауль- ская, пласт 26а, марка Ж Лава 26-25 290 Сразу при выходе из монтажной камеры Обрушение вслед за по-двиганием комплекса 7.5-8 20-26 340-526
3 Алар-динская, пласт 21, марка ГЖ Лава 21-1-52 «восток» 150 25 5 7 0.6-10 17-20 160-260 220 4.5
4 Ш/У Анжерское, пласт Румянцевский, марка ТС Лава № 31 30 20 6 17 4.3 2.2-4 105-150 190 1.5
5 Полосу-хинская, пласт 29а, марка ГЖ Лава 29-320 180 82 28 14 9.4-13.8 20-23 250-350 200 3.0
6 Ульяновская, пласт 50, марка ГЖ Лава 50-11 270-300 73 26 8 6 20 160-230 263 2.55
7 Ульяновская, пласт 50, марка ГЖ Лава 50-09 188 60 23 7 6 30 150 300 2.5
Понятно, что на величину зависающего участка кровли (фактически на скорость накопления повреждений) влияют многие факторы. Это и глубина выработки, и геометрия и размеры забойной зоны, и физикомеханические характеристики массива пород, и особенности геологии, обводненность, наличие крупных трещин и разломов и т.д. Упорядочить действие столь разнообразных факторов и учесть их влияние на деградацию прочностных свойств массива — очень непростая задача. В настоящем сообщении мы сделаем акцент на качественной стороне проблемы.
Ряд особенностей массива задается его исходными прочностными характеристиками и исходной поврежден-ностью. Изменения геометрии выработанного пространства, обусловленные подвиганием забоя, и силы тяжести формируют в массиве специфическое напряженно-деформированное состояние.
Скорость накопления повреждений в каждой точке нагружаемой среды зависит как от величины локальных напряжений в этой точке, так и от накопленных неупругих деформаций. Другой важнейший фактор — это деградация локальных прочностных свойств среды в зависимости от накопленных повреждений.
Таким образом, по известным данным величин первой посадки кровли и последующих обрушений для известных исходных данных глубины выработки и ее размеров тестируется модель, в которой устанавливаются две важнейшие характеристики: скорость накопления повреждений и функция деградации прочностных характеристик массива. Такова самая общая схема решения задачи.
2. Постановка задачи об обрушении кровли как задачи эволюции геосреды в поле сил тяжести
Ранее в работах [1-5] было показано, что уравнения механики сплошных сред описывают развитие деформационных процессов в нагружаемых твердых телах по законам синергетики. Более того, было показано, что прочные тела и среды являются идеальными нелинейными самоорганизующимися системами, так как в них имеются объемные нелинейные источники и стоки, которые обеспечивают накопление и перераспределение под внешними воздействиями деформационных дефектов и повреждений [2, 4]. Нелинейность среды и, прежде всего, нелинейные источники дефектов и повреждений приводят к локализации неупругих деформаций и повреждений, формируя в среде нестационарные эволюционирующие структуры (в геосреде это блочные структуры). Следовательно, под нестационарными структурами мы понимаем локализованные на определенных участках процессы деформирования, повреждаемости, деструкции. Определяющие уравнения реализуют также положительные и отрицательные обратные связи, причем отрицательные обратные связи обеспечивают
релаксацию напряжений, стабилизируя процесс. Роль положительных обратных связей иная: она двояка. С одной стороны, они, ускоряя процессы локализации, способствуют образованию диссипативных нестационарных структур, а с другой — ускоряя процессы деградации прочностных характеристик среды в областях локализации, приводят к сверхбыстрым режимам эволюции системы — режимам с обострениями [6], т.е. в данном случае к развитию катастрофического сверхбыстрого обрушения определенных участков кровли.
В табл. 2 для типичной геометрии призабойной зоны приведены некоторые данные о прочностных характеристиках почвы и кровли реального пласта. В расчетах забой от некоторого начального состояния (заштрихованная область в табл. 2) продвигался вдоль угольного пласта, наращивая выработанное пространство со скоростями проходки от 1 до 10 м/сут, что соответствует реальным скоростям подвигания забоя (табл. 1).
Полная система уравнений вместе с определяющими эволюционными уравнениями первой и второй группы подробно обсуждена в работах [2, 4]. Она включает в себя уравнения, выражающие законы сохранения:
Ау- да, і + р Аіу у = 0, р—- = —- +pFi
ф Аі =_ ді р *
Аі дх]
ді
(1)
Таблица 2
Стратиграфическая колонка пласта Владимировский II
Мощность, м Характеристики пласта и вмещающих пород
Основная кровля 4.5 Переслаивание песчаников с алевролитами
50.0 Песчаник средне-мелкозернистый, крепкий, устойчивый f = 6-8, р = 2.52 г/см3, К = 75-109 МПа
Непосредственная кровля 5.0-6.0 Алевролит с тонкими прослойками песчаника средней устойчивости f = 4-6, р = 2.41 г/см3, К = 70-86 МПа
Пласт угля 2.1—4.2 Щщ Уголь средней-слабой крепости, марка КО f = 0.6-1.0, р = 1.39 г/см3
ЙШ
Ложная почва 0.1-2.5 Алевролит серый средней крепо сти f = 4-6, р = 2.42 г/см3
Непосредственная почва 4.0-10.0
Основная почва 10.0 Песчаник средне-мелкозернистый, крепкий, устойчивый f = 6-8, р = 2.5 г/см3, К = 50-95 МПа
f — коэффициент крепости, р — плотность, К — временное сопротивление сжатию
эволюционные уравнения первой группы:
ст- =Х(0‘ -0р)8- + 2ц(ё- -ёр), (2)
где ст- =-р8у + Щ Щ, -Р = У3сти, р = /(р, £), эволюционные уравнения второй группы:
0р = А----В---0р + С (0),
I I
дх дх
(3)
например,
3 ёр
&р = 3 ёе® ЪУ
"е®
2 а,
где
& р.
р. = &.. + &
/ / /
= ё£
ё/ 2
Эу, Эу-
—Ц + —-
Эх1 Эх1
Здесь ёр® и ае1Г — вторые инварианты скоростей пластических деформаций и напряжений соответственно; А, В, С—некоторые функции, определяемые при выборе конкретной кинетики.
Скорость неупругой деформации ё р определяется заданной «поверхностью текучести» и «законом течения». Для хрупкопластичных геосред и геоматериалов под этими терминами мы понимаем законы ограничения напряжений.
Другой принципиально важной частью развиваемой эволюционной методологии является наличие так называемого характерного реального времени процесса Т. Это время задается условиями нагружения. Оно зависит от масштаба явления и может составлять секунды, минуты, часы, годы, миллионы лет. Время t в приведенной выше системе уравнений — это счетное время, по которому проводятся численные расчеты и для которого шаги по времени удовлетворяют известным условиям устойчивости разностной схемы. Характерное время процесса Т — это реальное время эволюции, определяемое, например, скоростью движения структурного элемента (блока выбранного масштаба) и скоростью накопления повреждений в геосреде на меньших масштабах.
Фактически задача моделирования эволюции системы на больших временах (вплоть до геологических времен во многие миллионы лет) решается методом установления. Необходимо, чтобы за некий интервал АТ, задаваемый условиями нагружения, волны напряжений, распространяющиеся со скоростью звука в нагружаемой среде, несколько раз пробегали по нагружаемой области, что обеспечивает квазистационарное на данном шаге нагружения распределение всех параметров. Следующий шаг по АТ приводит как к перераспределению напряжений и деформаций в среде, так и к дополнительному накоплению в ней повреждений и деградации прочностных характеристик, а также к локализации неупругих деформаций и повреждений.
Таким образом, характерное время Т зависит от рассматриваемой задачи, в том числе, и от пространственного масштаба (чем большие пространственные масштабы мы рассматриваем, тем больше времена Т: часы, дни, месяцы или, например, миллионы лет для плитной тектоники). Такой подход в геодинамике и механике позволяет решать задачи геодинамики как эволюционные задачи на любых характерных для изучаемого процесса временах с помощью динамических уравнений механики сплошных сред.
Развиваемая методология позволяет также автоматически выявлять сверхбыстрые катастрофические режимы и выполнять соответствующие расчеты развития режима с обострением. Если к некоторому шагу по времени процесса АТ, система оказалась подготовленной к сверхбыстрой эволюции, то последующий шаг АТ+1 позволит подробно изучить этот сверхбыстрый режим, даже если его скорости будут близки к ударноволновым, так как внутренние шаги А{1 << АТ позволяют детально изучать любые динамические процессы, в том числе и ударноволновые.
Подобная методология при численном моделировании требует наличия высокопроизводительных ЭВМ, поэтому для численного решения поставленных задач был использован вычислительный кластер СКИФ СуЬепа Томского государственного университета.
3. Модель среды
Рассматривается модель упруго-хрупкопластичной среды с учетом внутреннего трения, сдвиговой дилатан-сии (т.е. повреждений, появившихся вследствие сдвигов), накопления повреждений и деградации прочностных характеристик.
Наиболее простой формой уравнений, задающих предельную поверхность для чувствительных к давлению хрупких материалов, служат зависимости, включающие первый инвариант напряжений J1. В данной работе использовалась модель Николаевского [7], особенности которой подробно были изучены в работах [8, 9]. Для этой модели уравнение предельной поверхности, взятое в форме Мизеса-Шлейхера, имеет вид:
| Jl + Jf = У, (4)
а связь между сдвиговой и объемной неупругой деформацией задается выражением [7, 8]:
/р = 2Лр1/2, (5)
где а — коэффициент внутреннего трения; Л — в общем случае функция дилатансии.
В выражении для предельной величины напряжений Y деградация зависит как от накопленных неупругих деформаций, так и от локального напряженного состояния:
Рис. 1. Распределение повреждений и области разрушений в кровле при подвигании забоя со скоростью 1 м/сут на 20 (а), 27 (б), 36 сутки (в)
У = Уо(1 + А(е) - De(e))(1 - D(ст)). (6)
Здесь А(е) — функция, описывающая упрочнение материала (если оно есть), а также залечивание повреждений и частичное восстановление прочности. Функция De(e) описывает разупрочнение вследствие накопления повреждений в ходе неупругой деформации. Для упрочнения была принята линейная зависимость:
А(е) = &, (7)
е
а для разупрочнения — квадратичная:
De(e) = 2^2^^ . (8)
Здесь к1, к2 — коэффициенты модели; е* — критическая деформация, после достижения которой преобладает деградация материала.
Функция повреждаемости зависит от напряженного состояния и времени и задается выражением:
= |(СТ *ая°! & для ст>ст0. (9)
•’ (ст*)nt
где ст — эффективное напряжение; ст0, ст*, t* и п = 2 — параметры модели, определяющие пороговое напряжение, начиная с которого начинают накапливаться повреждения, предельное напряжение и скорость процесса накопления повреждений.
Это позволяет накапливать повреждения и учитывать деградацию свойств нагружаемого материала в упругой области, когда напряжения существенно меньше предела прочности среды (ст0 мало).
Усечение предельной поверхности в области растяжения также определяется с учетом повреждений. Этот критерий также отражает различие прочностных свойств геосреды при сжатии и растяжении. Фактически среда всегда разрушается в области растягивающих напряжений, абсолютная величина которых невелика и уменьшается с ростом поврежденности среды:
Р* = р*(1 - Да)). (10)
Таким образом, принимается, что разрушение и раскрытие трещин происходит при наличии растягивающих напряжений, которые присутствуют в локальных
Рис. 2. Распределение повреждений и области разрушений в кровле при подвигании забоя со скоростью 2 м/сут на 12 (а), 27 (б), 33 (в), 39 (г), 42 (д), 48 сутки (е). На 48 сутки непосредственная кровля сильно повреждена и рост повреждений продолжается
Рис. 3. Распределение повреждений и области разрушений в кровле при подвигании забоя со скоростью 4 м/сут на 9 (а), 12 (б), 15 (в), 20 (г), 27 (д), 28 сутки (е). На 27 сутки достигнут конец забоя, на 28 сутки развивается в катастрофическом режиме обвал основной кровли позади забоя
Рис. 4. Распределение повреждений и области разрушений в кровле при подвигании забоя со скоростью 8 м/сут на 9 (а), 18 (б), 20 (в), 24 (г), 27 (д), 28 сутки (е). На 18 сутки достигнут конец забоя, забой не удлиняется; на 20 сутки идут накопления повреждений; на 27 сутки происходит катастрофа, обвал основной кровли позади забоя; на 28 сутки развивается обвал основной кровли позади забоя
областях неоднородной среды практически при любом виде нагрузок. Определение параметров модели в выражениях (4)-(10) — отдельная задача. Особенно сложно их подбирать для геосред, так как непосредственных экспериментальных данных нет, а есть только ограниченные наблюдения за скоростями деградации. В настоящей работе эти параметры выбирались как на основании известных данных для геоматериалов, так и по данным о временах посадок и обрушений кровли (см. табл. 1).
4. Результаты расчетов
Расчеты проведены в плоской двухмерной постановке с применением программы, использующей параллельный многопроцессорный счет, по явной разностной схеме второго порядка точности. Представленные расчеты следует рассматривать как тестовые, хотя они хорошо качественно отражают эволюцию горного массива над выработанным пространством. Все приведенные расчеты выполнены для горизонтально залегающего пласта (H = const) и средней глубины H = 200 м, что об-
легчает сравнение результатов для разных скоростей по-двигания забоя. Скорость подвигания забоя варьировалась приблизительно на порядок — от 1 до 8 м/сут.
На рис. 1, 2 представлены результаты для малой скорости подвигания забоя (1 и 2 м/сут). Длина расчетной области — 256 м. Выработанное пространство на рисунках не изображено, оно находится непосредственно снизу и растет слева направо от 10 до 200 м. Все деформации и повреждения показаны как градации серого цвета, области разрушений выделены черным цветом. Сопряженные системы полос локализации неупругих сдвигов и повреждений выявляют иерархию формирующихся блоков. В обоих случаях фронт повреждений двигается практически вместе с забоем и накопленные повреждения далеки от критических только в начальные моменты времени. В нагружаемой геосреде как неупругие деформации, так и повреждения формируют своеобразный купол. Такие режимы с малыми скоростями подвигания забоя можно отнести к квазиравновесным, когда критические значения повреждений достигаются при пролетах зависающей кровли не более 20-40 м (глубина Н сравнительно небольшая и повреждения копятся достаточно медленно, так как массив пород прочный (табл. 2)).
Существенно неравновесные режимы накопления повреждений и катастрофические обвалы зависшей кровли (порядка 100 м длиной) показаны на рис. 3, 4. При скорости подвигания забоя 4 м/сут обвал произошел на 28 сутки, при скорости 8 м/сут — несколько раньше, на 27 сутки. Так как все прочие условия, кроме скорости подвигания забоя, одинаковы, то скорости накопления повреждений почти равны. При большей скорости движения забоя больше длина зависшего пролета забоя и в оставшееся время до катастрофы на массив действуют несколько большие напряжения, поэтому катастрофа происходит на сутки раньше.
5. Заключительные замечания
В работе обсуждена модель эволюции горного массива в полях силы тяжести, основанная на решении нелинейных динамических уравнений механики деформируемого твердого тела. Выполнены расчеты развития повреждений в кровле над очистным забоем вплоть до ее обрушения. Показано, что нестационарные существенно неравновесные режимы эволюции кровли над вы-
работанным пространством реализуются при высоких скоростях подвигания забоя (более 4 м/сут в приведенных примерах) и приводят к зависанию протяженных участков кровли и ее катастрофическим обвалам. Этими расчетами показано, что эволюционные процессы и самоорганизация в нелинейных сплошных средах, в данном случае в геосредах и горных массивах, описываются нелинейными уравнениями в частных производных математической физики (механики сплошных сред), демонстрируя характерные особенности эволюции нелинейных динамических систем — локализацию деформаций и повреждений, образование субструктур и сверхбыстрые (катастрофические) режимы эволюции. Развиваемая методика неявно содержит времена реальных процессов, что позволяет решать задачи механики деформируемых твердых тел как задачи эволюции на любых больших временных отрезках вплоть до геологических времен для геосред.
Работа выполнена при поддержке СО РАН (проекты фундаментальных исследований на 2007-2009 гг. №№ 7.11.1.6 и 3.6.2.3) и программы Президиума РАН (проект 16.3).
Литература
1. Макаров П.В., Трубицын А.А., Ворошилов С.П. Самоподобие разру-
шения углей и эволюция нагружаемых твердых тел // Уголь. -2006. - № 10. - С. 55-58.
2. Макаров П.В. Эволюционная природа деструкции твердых тел и сред // Физ. мезомех. - 2007. - Т. 10. - № 3. - С. 23-38.
3. Макаров П.В. Нагружаемый материал как нелинейная динамическая система. Проблемы моделирования // Физ. мезомех. -2005.- Т. 8. - № 6. - С. 39-56.
4. Макаров П.В. Эволюционная природа блочной организации геоматериалов и геосред. Универсальный критерий фрактальной делимости // Геология и геофизика. - 2007. - Т. 48. - № 7. - С. 724-746.
5. Макаров П.В., Трубицын А.А., Трубицына Н.В., Кузнецов П.В., Смолин И.Ю., Стефанов Ю.П., Ворошилов С.П., Ищук И.Г., Голосков С.И. Численное изучение разрушения угля на мезо- и макроуровнях // Уголь. - 2005. - № 2. - С. 33-36.
6. Князева Е.Н., Курдюмов С.П. Законы эволюции и самоорганизации
сложных систем. - М.: Наука, 1994. - 236 с.
7. Гарагаш И.А., Николаевский В.Н. Неассоциированные законы течения и локализации пластической деформации // Успехи механики. - 1989. - Т. 12. - № 1. - С. 131-183.
8. Стефанов Ю.П. Локализация деформации и разрушение в геоматериалах. Численное моделирование // Физ. мезомех. - 2002. -Т.5. - №5. - С. 107-118.
9. Стефанов Ю.П. Некоторые особенности численного моделирования поведения упруго-хрупкопластичных материалов // Физ. ме-зомех. - 2005. - Т. 8. - № 3. - С. 129-142.
Поступила в редакцию 21.10.2007 г.
