CC BY
120
Модель гипоупругой хрупкой среды: применение к расчету деформирования и разрушения горных пород
М.М. Немирович-Данченко
Томский филиал Института геологии нефти и газа СО РАН, Томск, 634055, Россия
В работе рассмотрено конечное деформирование и разрушение горных пород. Изучение проводится в рамках численного моделирования конечно-разностным методом в переменных Лагранжа. Для описания конечного деформирования предложено использовать модель гипоупругой среды, а специальная методика разделения узлов расчетной сетки позволяет рассчитывать разрушение породы. Изучено деформирование среды с имеющейся трещиной под действием растяжения и сдвига. Оценивается излучение сейсмических волн при раскрытии трещин.
1. Введение
Хорошо известно, что обычно в горных породах имеют место деформации, не всегда достаточно малые для того, чтобы проводить их расчет в рамках классической теории упругости. Кроме того, геологические и техногенные процессы, сопровождающиеся конечными деформациями среды, часто влекут за собой разрушение и/или относительное скольжение в горных породах [1]. Типичные примеры — подготовка землетрясения и поведение горных пород непосредственно в зоне очага, обрушение кровли выработки и т.п. Следствием и зачастую единственным внешним проявлением разрушения горных пород является высвобождение энергии и распространение ее в виде упругих (сейсмических) волн.
В последнее время стали появляться работы, в которых моделирование деформирования и разрушения среды используется для более адекватного описания реальных процессов, влекущих за собой излучение сейсмических волн [2, 3, 4]. Численные методы обычно применяются по следующей схеме [5]: строится подробная дискретная модель среды, в гипоцентре задаются функции источника, затем рассчитывается система волн. При таком подходе, который имеет и свои преимущества, нет собственно расчетов подготовки (нуклеации), рассматривается просто свершившееся событие, т.е. действие наперед заданного источника в известной среде. При использовании аналитических и полуанали-тических методов рассмотрение проводится в рамках линейной механики разрушения, что позволяет получать фундаментальные оценки и качественно описывать реальный процесс. Здесь один из наиболее интересных
результатов приведен в [4], где в полной постановке решена вся задача до акта излучения, само излучение не рассмотрено.
Вычислительные методы механики сплошных сред, например [6], позволяют проводить полновесные численные эксперименты в прямых задачах сейсмики. В настоящей статье излагается наш подход к моделированию деформирования и разрушения горных пород. Для описания деформирования используется модель гипоупругой среды, определяющие соотношения которой приводятся в разделе 1.1. Кроме того, применяется специальная методика для описания процесса разрушения среды и для расчета движения разрушенной среды.
2. Особенности математического моделирования для больших деформаций и разрушения среды
1.1. Модель гипоупругой среды
У Жермена [7] дается следующее понятие гипоупругой среды: среда является гипоупругой, если в каждой точке и в любой момент времени тензор скоростей изменений напряжений есть линейная функция тензора скоростей деформаций, причем эта функция может, в свою очередь, зависеть от тензора напряжений как от параметра.
Скорость изменения тензора напряжения может задаваться различным образом. Так, если заданы поле скоростей V(r, t) и поле напряжений ак (г, t), то можно записать следующие выражения для обычной полной производной изменения компонент аЛ относи-
© Немирович-Данченко М.М., 1998
тельно пространственной системы координат:
_ да»_V
/0 0 0 ,ч
хі = & (х1> х2 >х3 >* )’
dt
дt дх,
(1)
где в правой части первый член — локальная часть (локальная производная), а второй член — конвективная часть.
Определенная таким образом скорость изменения напряжений будет зависеть от собственного вращения элемента среды и не является поэтому объективной (материальной) величиной. Для построения закона поведения необходимо использовать производную тензора, которая будет обращаться в нуль при вращении тела как твердого относительно реальной системы отсчета. Яуманном было предложено для этого кроме локального и конвективного членов в (1) рассматривать еще вращательный член [8]
(2)
В последней формуле использованы компоненты тензора скоростей вращения. Вычитание такого члена из полной производной учтет то обстоятельство, что напряженное состояние частицы относительно неподвижной системы координат не должно меняться, если это состояние не меняется относительно координатной системы, которая вращается вместе с окрестностью частицы (так называемой вмороженной системой координат). Итак, производная относительно собственного вращения или производная Яуманна имеет вид
йа"’ (3)
А =-
----- Ю;.СТ..
dt 4 Д * 4
Тензор скоростей деформаций записывается обычным образом
11 ди, ди,
2 I дхдх, I ’
(4)
а закон поведения гипоупругои среды выглядит так:
(5)
В общем случае сгктп есть функция а гк; если же это не так, то для изотропной среды, как обычно,
Сгктп _ *8 тп + М-(8 ш8 кп + § т8 кт ) •
Важнейшей характеристикой (5) является то, что это не есть закон поведения в обычном смысле, позволяющий точно сказать, чему будет равна деформация при данном значении напряжения. Закон (5) по сути своей инкрементальный, что делает его наиболее удобным для численного моделирования в нестационарных задачах.
1.2. Описание разрушения при численном моделировании
В механике сплошных сред движение считается заданным, если известны три функции
отображающие координаты х° в х1 при перемещении материальной точки М0 в положение М. Система всех частиц {М} , составляющих сплошную среду, может быть охарактеризована некоторой начальной конфигурацией {М}0 и семейством множеств {М}‘, описывающих движение среды. Взаимнооднозначное соответствие между {МУ1 и {М}2 определяется функциями gi из (6).
Как правило, предполагается, что gi непрерывны и непрерывно-дифференцируемы по своим аргументам столько раз, сколько это окажется необходимым для рассуждений. При изучении процессов разрушения очевидна необходимость того, чтобы относительно переменных х°, х20, х3 функции gi были бы кусочнонепрерывными и непрерывно-дифференцируемыми.
При деформировании в среде происходят разрушения различного характера и масштаба. Для численного моделирования существенным является отношение характерного размера расчетной ячейки Дг к некоему характерному размеру разрушения I. При Дг/1 >> 1 разрушение можно трактовать как появление небольших пор и микротрещин. Такое разрушение не требует изменения расчетной сетки и учитывается в осредненных значениях модулей упругости. В этом случае последующий расчет с новыми константами проводится как для сплошной среды. Очевидно, для такого разрушения сохраняется взаимнооднозначное отображение (6).
Однако на любом масштабном уровне при нагружении геологической среды наиболее характерна ситуация Дг ~ I (или даже Дг/1 << 1). Это — разрушение типа трещинообразования с высвобождением значительной энергии и с нарушением диффеоморфизма (6) (из одной математической частицы среды получаются, к примеру, две — одна отходит к правому берегу трещины, другая — к левому).
Для численного моделирования такого разрушения автором используется следующий метод, в простейшем случае упомянутый в [6, 9].
Предположим, что при аппроксимации уравнений движения координаты и скорости (ускорения) определяются в узлах лагранжевой расчетной сетки. Предположим далее для простоты, что в расчетах используются обычные прямоугольные в начальный момент времени ячейки. Тогда в каждом внутреннем расчетном узле сходятся четыре угла соседних ячеек. Будем такой узел мысленно считать состоящим из 4 точек, рис. 1. Зададим с самого начала вычислений в каждом узле не один набор координат и скоростей, как это делается традиционно, а четыре. Первый набор координат и скоростей (точка 1) отнесем к верхней правой ячейке, второй набор (точку 2) — к верхней левой и т.д. Пока тело сплошное, между четырьмя точками имеются четыре
2(1, J) 1 (I, J)
щ 4 % 4
3(1, J) _ 'iU 4(l, J)
1+1
Рис. 1. Часть модифицированной расчетной сетки с мысленными связями между узлами
связи и четыре точки сливаются в одну (т.е. для сплошных участков тела эти наборы совпадают). Если разрывающее усилие в любой связи достигло предела или выполнился иной критерий разрушения, связь рвется. В этом случае для четырех точек разрушенного узла записываются граничные условия с учетом вновь образованных свободных поверхностей. Для точек, соседних с берегами трещины, будут иначе рассчитываться как компоненты скорости, так и компоненты тензора скоростей деформаций. При таком описании каждый расчетный узел состоит из отдельных лагранжевых точек, имеющих при разрушении различные траектории, причем сохраняется взаимнооднозначное соответствие, определяемое в (6).
3. Концентрация напряжений, разрушение и излучение сейсмических волн для трещин моды I и моды II
Все расчеты, результаты которых приводятся в данной работе, проводились по конечно-разностной схеме, приведенной в [10], с определяющими соотношениями
(5) и описанным в п. 1.2 алгоритмом расчета разрушения. Рассчитан один вариант нагружения по моде I (для трещины отрыва), рис. 2, а, и два варианта нагружения по моде II (для трещины сдвига), рис. 2, б, в.
2.1. Трещина моды I (трещина отрыва)
Рассмотрим однородный образец модельного материала с трещиной, рис. 2, а. Скорость продольной волны в среде принимаем 3 000 м/с. Размеры расчетной области 800 X 400 м. К верхней и нижней поверхностям образца приложено растягивающее напряжение, на боковых гранях ставится условие свободной поверхности (это трещина моды I или трещина отрыва). В терминах п. 1.2 такая трещина описывается следующим образом: по всей линии трещины разорваны связи 1-4 и 2-3, а в вершине — 2-3. (В силу осевой симметрии рассмотрена только правая часть образца).
Основное физическое явление, связанное с трещинами и влияющее на деформирование реальных тел, — это концентрация напряжений. На рис. 3. приведена динамика изменения компоненты напряжений о yy (x), коллинеарной нагрузке, вдоль линии трещины для нескольких моментов времени.
Для безразмерного момента времени Т4 на рис. 4 приведены изолинии максимальных касательных напряжений. Для этой задачи хорошо известно качественное поведение поля напряжений в окрестностях трещины, поэтому рис. 3 и рис. 4 приводятся здесь как результаты тестирования. Можно сделать вывод, что полученные эпюры концентраций напряжений и изолинии максимальных касательных напряжений находятся в хорошем соответствии с инструментальными данными.
На рис. 5 приведен фрагмент векторного поля скоростей смещений расчетных точек в окрестности вершины трещины для начальной стадии взаимодействия трещины с волной нагрузки. По фрагменту векторного
Рис. 2. Исходная геометрия задач: а — трещина отрыва; б, в— трещины сдвига; 1 — трещина; 2 — направление действия прикладываемой нагрузки
Рис. 3. Динамика концентрации напряжений вдоль линии трещины. Эпюры показаны для пяти безразмерных моментов времени Т= ^/Ь, где V— скорость продольной волны в среде; t—текущее время с начала растяжения; Ь — длина расчетной модели. Т1 = 1.25, Т2 = 1.87, Т3 = 2.50, Т4 = 3.12, Т5 = 3.75. В скобках показаны максимальные значения эпюры в килобарах (1 кбар = 108 Па). 1 — положение трещины, 2 — боковая поверхность, свободная от напряжений
поля можно судить о характере потока энергии в вершине трещины. Аткинсоном и Эшелби [11] было показано, что для статических (квазистатических) трещин вершина является энергетическим стоком. Для трещины моды I это хорошо видно из приведенного на рис. 5 векторного поля. Вихревой характер поля скоростей для динамической постановки задачи также естественным образом вытекает из наличия трещины.
Таким образом, изложенная методика качественно хорошо описывает явление концентрации напряжений при деформировании тела с трещиной по моде I.
Далее рассмотрено разрушение в вершине трещины. Нами при расчетах использовался интегральный критерий накопления повреждений [12], выполнение которого проверялось в каждой расчетной ячейке. Собственно критерий и соответствующие константы горных пород здесь не обсуждаются. При выполнении критерия обра-
/о.срУ— -<1Ш5
г'*™“—' \
зуются новые свободные поверхности, как это описано в п. 1.2. Применение специальной методики позволяет выделить из напряженно-деформированного состояния ту часть, которая порождена высвобождением упругой энергии при росте трещины.
Обозначим время выполнения критерия Т*. На рис. 6 приводится векторное поле скоростей смещений в момент Т* + 0.008 с, а на рис. 7 — в момент Т* + 0.04 с. Хорошо видно, что разрушение с образованием свободных поверхностей приводит к излучению волн всех типов. Так, вдоль линии трещины бежит продольная волна, порождающая головную коническую волну, рис. 7. Амплитуда этих волн на поверхности трещины отлична от нуля. Наиболее мощной волной является волна Рэлея, а в направлении, нормальном к разрыву (т.е. направлении нагружения), значителен вклад продольной волны.
По результатам численных расчетов в полярных координатах построены интенсивности излучения, зна-
R
чения вычислялись с шагом 4.5° по формуле Iv dr ,
0
т.е. для данного момента времени просуммирована вся излученная в данном направлении энергия. Можно сказать, что эта функция характеризует направленность источников, возникающих при разрыве среды. На рис. 8 построены графики интенсивности полей скоростей для моментов времени Т* + 0.008 с и Т* + 0.04 с. Видно, что наибольшая часть энергии вначале излучается в направлении приложения нагрузки; затем происходит перераспределение энергии, после чего основная часть перемещается в виде головной и рэлеевской волн вдоль линии разрыва назад от направления разрыва.
2.2. Трещина моды II (трещина сдвига)
На рис. 2, в показан один из возможных способов нагружения сдвигом. На нижней и верхней гранях заданы небольшие х- компоненты скорости смещения, у-
Рис. 4. Изолинии максимальных касательных напряжений. Числовые значения даны в ГПа. Жирной линией показана трещина
Рис. 5. Векторное поле скоростей смещений (локальных скоростей) в окрестности вершины трещины. Жирной линией показана трещина
Рис. 6. Векторное поле скоростей смещений для трещины моды I в момент времени Т* + 0.008 с
компонента скорости зафиксирована. На рис. 9 приведены кадры поля скоростей вокруг трещины для восьми последовательных моментов времени. Кроме того на каждом кадре по координатам построены берега трещины; критерий разрушения не проверялся. На первом кадре нагрузка еще не дошла до линии трещины, течение можно считать одномерным. На втором кадре волна нагружения достигла линии трещины, все поле скоростей однородно, за исключением небольшой области высотой в 4 расчетных ячейки с каждой стороны вокруг трещины. Образно говоря, расчетное тело еще «не знает», что имеется трещина; информация об этом распространяется с конечной скоростью продольной волны. На кадре 3 линия трещины уже искажена, появились зоны интерференции. Кадр 4 весьма характерен. Трещина является, как и при отрыве, источником вих-реобразования. Светлая разреженная область, охватывающая центр трещины, содержит два небольших вихря. На кадре 5 область разрежения расширилась, образовав на концах трещины две полосовые структуры, вихрь вдоль верхнего берега распространяется вправо, вдоль нижнего берега — влево, что отчетливо видно из кадра 6. На кадрах 6 и 7 видно не только отрывание вихрей от вершин трещины, но и заметно, как изменение напряженного состояния приводит к раскрытию трещины. Векторы скоростей смещений направлены под значительным углом к линии трещины. На кадре 8 хорошо видно равномерное раскрытие трещины и перераспределение векторного поля скоростей смещений.
Если при нагружении такого типа не фиксировать скорость границ в направлении оси Y, то расчетная область начинает разворачиваться, что также приводит к плавному переходу от чистого сдвига к сдвигу с отрывом. На рис. 10 показано векторное поле скоростей после такого деформирования, причем приведена почти вся расчетная область (200x100 ячеек), а на рис. 11 дана увеличенная картина вблизи трещины. Хорошо видно,
Рис. 7. Векторное поле скоростей смещений для разрушения по моде I в момент времени Т* + 0.04 с: 1 — сферический фронт продольной волны; 2 — фронт поперечной волны; 3 — коническая волна; 4 — волна Рэлея. Трещина обозначена жирной линией
что первоначально однородное тело из-за наличия трещины стало в напряженном состоянии существенно неоднородным. Попытка аккомодации к нагрузке приводит в таком теле к образованию особенностей поля скоростей смещений, имеющих внешний облик полосовых, гиперболических, вихревых и прочих структур. Перейдем к задаче рис. 2, б. Здесь нагрузка задана в верхней части левого торца среды, нижняя часть левого торца при этом жестко закреплена. Для этой трещины
Рис. 8. Интенсивности излучения полей скоростей для трещины моды I: 1 — трещина; 2 — направление скачка трещины; 3 — график для момента времени Т* + 0.008 с; 4 — график для момента времени Т„ + 0.04 с
Рис. 9. Численные снимки процесса сдвигового нагружения тела с трещиной. На первом кадре стрелками показано направление нагрузки
приведен фрагмент расчетной сетки до разрушения в вершине трещины, рис. 12.
Хорошо виден сдвиг верхнего берега трещины относительно нижнего; относительное сжатие ячейки достигает 10 %. Необходимо подчеркнуть, что именно применение модели гипоупругого тела, по сути описываемой инкрементальным законом поведения, позволяет рассчитать конечные обратимые деформации горной породы. Выделение в тензоре напряжений шаровой части и девиатора позволяет легко дополнить алгоритм условием пластичности и проводить расчеты в рамках модели гипопластического тела.
Немалую роль в волновой картине для трещины моды II играет характер единичного скачка. Для задачи, геометрия которой показана на рис. 2, б, единичный скачок приведен на рис. 13. Здесь также в каждой точке проверялось выполнение интегрального критерия [12]. На рис. 14 показана развитая картина разрушения при продолжении действия нагрузки после первого скачка.
Рис. 10. Численный снимок деформирования тела с трещиной сдвиговой нагрузкой с незафиксированными границами. В центре видна раскрывшаяся трещина
Рис. 11. Поле скоростей вблизи трещины. Отчетливо видны два вихря, расположенные антисимметрично линии трещины
Волновые картины для развитого разрушения не приводятся.
Естественно, нас интересовал характер сейсмического источника, возникшего при разрушении по моде II. На рис. 15 приведено поле скоростей для момента времени Т* + 0.008 с, а на рис. 16 — для момента времени Т* + 0.04 с. Соответствующие функции направленности
Рис. 12. Фрагмент расчетной сетки до скачка трещины
}
Рис. 13. Фрагмент расчетной сетки после скачка трещины
■
Рис. 14. Развитая картина разрушения
Рис. 15. Векторное поле скоростей смещений для трещины моды II в момент времени Т* + 0.008 с
(интенсивности излучения) приводятся на рис. 17. Итак, для трещины моды II рассчитанное поле принципиально отличается от трещины моды I. Прежде всего, из приведенного на рис. 16 численного снимка видно появление моментного источника. Свободная поверхность берегов трещины по-прежнему является зоной распространения головной конической волны, но в случае трещины моды II головная коническая волна несет энергию, сопоставимую с энергией волны Рэлея. В интегральном виде это сказывается на интенсивности излучения, рис. 17. На рис. 13 и 14 можно видеть, что основная часть энергии для рассчитанного акта разрушения распространяется в нижнем фрагменте в на-
Рис. 16. Векторное поле скоростей смещений для трещины моды II
Рис. 17. Интенсивности излучения полей скоростей для трещины моды II: 1 — трещина; 2 — направление скачка трещины; 3 — график для момента времени Т* + 0.008 с; 4 — график для момента времени Т„ + 0.04 с
правлении, обратном направлению разрыва. Головная волна также, как и для моды I, бежит вдоль берега трещины, принадлежащего нижнему фрагменту.
Литература
1. McGarr A., Spottiswoode S.M., Gray N.C., Ortleep W.D. Observation relevant to seismic driving stress, stress drop and efficiency // J. Geophys. Res. - 1979. - V. 84. - P. 2251-2261.
2. Rice J.R. The Mechanics of Earthquake Rupture // Physics of the Earth’s Interior. - Amsterdam, Holland, 1980.
3. Kame N., Yamashita T. Dynamic nucleation process of shallow earthquake faulting in a fault zone // Geophys. J. Int. - 1997. - V. 128. -No. 1.- P. 204-216.
4. Дядъков П.Г., Назаров Л.А., Назарова Л.А. Численное моделирование напряженного состояния земной коры и условий возникновения динамической неустойчивости сейсмоактивных разломов при рифтогенезе // Геология и геофизика. - 1997. - Т. 38. - № 12. -С. 2001-2010.
5. OlsenK.B., ArchuletaR.J. andMatarese J.R. Three-Dimensional Simulation of a Magnitude 7.75 Earthquake on the San Andreas Fault// Science. - 1995. - V. 270. - С. 1628-1632.
6. Wilkins M.L. Calculation of elastic-plastic flow // Methods in Computational Physics, Vol. 3, edited by B. Alder, S. Fernbach and M. Ro-tenberg. - New York: Academic Press, 1964. - P. 211-263.
7. Жермен П. Курс механики сплошных сред. Общая теория. - М.: Высшая школа, 1983. - 399 с.
8. Прагер В. Введение в механику сплошных сред. - М.: Иностр. лит-ра, 1963. - 312 с.
9. Гриднева В.А., Немирович-Данченко М.М. Метод раздвоения точек
сетки для численного расчета разрушения твердых тел. - Томск: ТГУ, 1983. - 12 с. - Деп. в ВИНИТИ 14.06.83, № 3258.
10. Немирович-Данченко М.М., Стефанов Ю.П. Применение конечно-разностного метода в переменных Лагранжа для численного расчета волновых полей в сложнопостроенных средах // Геология и геофизика. - 1995. - № 11. - C. 94-105.
11. Atkinson C., Echelby J.D. The flow of energy into the tip of a moving cracks // Int. J. Fracture Mech. - 1968. - V. 4 (1). - P. 3-8.
12. Гриднева В.А., Корнеев А.И., Трушков В.Г. Численный расчет напряженного состояния и разрушения плиты конечной толщины при ударе бойками различной формы // Изв. АН СССР. МТТ. -1977.- №. 1.- С. 146-157.
