Физика Земли, 2023, № 3, стр. 158-167

ВВЕДЕНИЕ

За последние 20–30 лет проведено большое число экспериментов, в которых исследовалось развитие разрушения в горных породах [Lockner et al., 1991; Botvina, 2011; Petružálek et al., 2013; Hamiel et al., 2006; Kuksenko et al., 1996; Smirnov et al., 2010; Xinglin Lei, Shengli Ma, 2014; Tal et al., 2020]. Основной задачей этих исследований был поиск признаков приближающегося катастрофического разрушения (в случае лабораторных образцов – образование магистральной трещины). Однако до сих пор физические механизмы, управляющие переходом от “неопасного” делокализованного образования трещин к ускоренному локализованному процессу, остаются невыясненными [Botvina, 2011; Carpinteri et al., 2011; Ponomarev et al., 1997; Naimark, 2003; Panteleev et al., 2012; Lei, 2019]. В связи с этим, мы решили обратиться к компьютерному моделированию, предполагая, что численные эксперименты позволят не только выявить факторы, существенным образом влияющие на развитие процесса, но и предложить новые схемы лабораторных экспериментов.

Большинство ранее предложенных компьютерных моделей (например, [Krajcinovic, 1989]), используемых для описания механического поведения горных пород, основаны на представлениях механики сплошных сред, которые не позволяют явно учесть образование дефекта (нарушение сплошности).

Однако для понимания физических закономерностей образования и эволюции дефектов (трещин), наиболее адекватными являются модели, основанные на методе дискретных элементов (DEM) [Lisjak, Grasselli, 2014; Cundall, 1971]. Модель дискретных элементов позволяет явно учитывать локальные нарушения сплошности в процессе деформирования, естественным образом имитируя образование и развитие трещин.

Рис. 1.

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

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

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

ОПИСАНИЕ КОМПЬЮТЕРНОГО ЭКСПЕРИМЕНТА

В работе использовалась модель связных частиц (bonded particle model – BPM), которая подробно описана в работе [Potyondy, Cundall, 2004]. Различные модификации этой модели применяются для изучения поведения материалов при механическом разрушении. Модель материала (горной породы) – сферические частицы одного или разных размеров, которые имитируют зерна, и связи между частицами, которые имитируют межзеренные границы. Физико-механические параметры материалов, из которых состоят частицы и связи, приведены в таблице.

Модельные эксперименты выполнялись в свободно распространяемом пакете программ MUSEN [Dosta, Skorych, 2020].

Связи формировались в местах контактов частиц. Материал связей выбирался из таблицы. Частицы из одного материала соединялись связью из того же материала, а частицы из разных материалов соединялись либо низкопрочными хрупкими стеклянными связями, либо низкомодульными связями (набор связей типа 6, 7, 8 из таблицы). Диаметр связей (d) выбирался генератором связей автоматически равным меньшему диаметру пары соединяемых частиц 1 и 2: $d = \min \left\{ {{{d}_{1}},{{d}_{2}}} \right\}$ (d₁, d₂ – диаметры частиц 1 и 2) [Dosta, Skorych, 2020]. Максимальная длина связей выбиралась из того условия, чтобы между парой соединяемых частиц не могла поместиться еще одна частица. Минимальная длина ${{L}_{{{\text{min}}}}}$ обычно была нулевой.

Для описания трения на контактах между частицами, а также между частицами и поверхностью плит, использовалась модель Герца–Миндлина [Hertz, 1882; Tsuji et al., 1992]. Параметры трения на контактах не представляются важными для данной задачи, так как этот вид взаимодействия вступает в силу после разрыва связей. Он будет определять, например, проскальзывание двух частей образца относительно друг друга, после его разрушения. Поэтому выбирались значения параметров по умолчанию: коэффициент восстановления = 0.6; коэффициент трения скольжения = 0.45, коэффициент трения качения = 0.05. При этих значениях параметров проскальзывание в зажимах отсутствовало.

Моделировались образцы цилиндрической формы диаметром 10 мм и высотой 20 мм. Цилиндры заполнялись сферическими частицами одного или разных размеров и упаковывались генератором упаковки MUSEN до достижения пористости 0.35−0.37. Образец помещался в виртуальный пресс. Нижняя плита была неподвижна, а верхняя перемещалась вниз с постоянной скоростью v = 0.02 м/с. Тем самым имитировалось одноосное сжатие.

В настоящей работе рассмотрено изотермическое нагружение, поэтому вопросы, связанные с характерным временем, определяемым теплопроводностью, не рассматриваются, предполагая, что тепло передается мгновенно. Характерное время релаксации упругих напряжений определяется скоростью звука ${{{{\tau }}}_{e}} = \frac{H}{V} \approx 4 \times {{10}^{{ - 6}}}~~{\text{с}}$ , что намного меньше времени разрушения образца, а скорость звука в гранитах 5400 м/с намного превосходит скорость нагружения. С этой точки зрения нагружение медленное. Поскольку расчет делается с шагом ~5e–9 с, за время жизни образца ~0.02 с система проходит около 4 000 000 состояний, что сводит динамические эффекты к минимуму [Brown, 2013]. Наконец, процессы, связанные с конечной вязкостью, в данной работе тоже не рассматриваются.

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

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

В гомогенном образце все частицы и связи состоят из одного материала. В гетерогенном образце (как и в реальных материалах) присутствуют частицы разных минералов и связи (межзеренные границы) также имеют свойства разных материалов: гранит, ортоклаз, олигоклаз, кварц, стекло. Изменялась степень гетерогенности модельного материала: варьировались размеры частиц, длина и толщина связей, материал частиц и связей. Число частиц гомогенного образца – 28 125, гетерогенного – 48 695.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

На рис. 2 приведены диаграммы нагружения гомогенного и гетерогенного и образцов. Для вычисления деформации использовалась формула $\varepsilon = {{\left( {vt} \right)} \mathord{\left/ {\vphantom {{\left( {vt} \right)} H}} \right. \kern-0em} H}$, где H – высота образца в начальный момент времени. Напряжения рассчитывались исходя из сил, действующих на нагружающие плиты. Поскольку в численном эксперименте невозможно в общем случае выдержать равенство сил, действующих на плиты, напряжение рассчитывалось по формуле $\sigma = {{0.5\left( {{{F}_{t}} + {{F}_{b}}} \right)} \mathord{\left/ {\vphantom {{0.5\left( {{{F}_{t}} + {{F}_{b}}} \right)} S}} \right. \kern-0em} S}$, где индексы t и b означают верх и низ, S – сечение образца (использовалось сечение в начальный момент времени) [Brown, 2013].

Рис. 2.

Диаграммы нагружения (напряжение (σ) – деформация (ε) ) при деформировании гомогенного (a) и гетерогенного (б) образцов.

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

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

Модель позволяет изучать “акустическую эмиссию”. В первом приближении можно считать, что каждый разрыв связи (разрушение межзеренной границы) сопровождается излучением сигнала акустической эмиссии. На рис. 3 показаны графики зависимости числа разорванных связей за определенный промежуток времени – аналог активности акустической эмиссии (красные кривые). Хорошо видна разница в поведении “акустической активности” (AA) при деформировании гомогенного и гетерогенного материала.

Рис. 3.

Изменение напряжения (синие кривые) и активность акустической эмиссии (красные кривые) при деформировании гомогенного (a) и гетерогенного (б) образцов.

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

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

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

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

В гомогенном образце наблюдается хрупкое разрушение: незначительное уменьшение числа целых связей в течение длительного времени (рис. 4а). При этом коэффициент вариации сохраняет низкие значения (рис. 4а). Это говорит об однородности разрушения по всему объему. Затем наблюдается резкое уменьшение числа целых связей. При этом происходит взрывоподобный скачок коэффициента вариации. Это указывает на локализацию разрушения и прорастание трещины.

Рис. 4.

Зависимости от времени числа целых связей усредненных по слоям (N) и коэффициента вариации (CV): (а) – гомогенный образец; (б) – гетерогенный образец.

На рис. 4б показано изменение этих же параметров при разрушении гетерогенного образца. Видно, что в гетерогенном образце происходит накопление повреждений при значительно меньших временах. Однако коэффициент вариации на этой стадии также небольшой, что говорит о том, что повреждения накапливаются более или менее однородно по объему образца. Это подтверждает справедливость модели двухстадийного разрушения гетерогенных материалов [Kuksenko et al., 1996].

Далее был проведен анализ растягивающих напряжений на связях.

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

Рис. 5.

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

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

Для проверки адекватности построенной компьютерной модели реальным процессам, происходящим при деформировании образцов горных пород, необходимо выполнить сопоставление с лабораторным экспериментом. У нас имеются обширные экспериментальные данные, полученные при деформировании образцов различных горных пород [Damaskinskaya et al., 2018; 2021]. Обратимся к данным экспериментов, в которых проводилось одноосное сжатие цилиндрических образцов (d = 10 мм, h = 20 мм) песчаника Berea и гранита Westerly.

Для качественного сопоставления с компьютерными экспериментами выбраны именно эти горные породы, поскольку песчаник Berea является более однородным материалом по сравнению с существенно гетерогенным гранитом Westerly. Это демонстрируют томографические срезы данных материалов, показанные на рис. 6.

Рис. 6.

Примеры томографических срезов: (а) – песчаник Berea; (б) – гранит Westerly. Градации серого соответствуют различной плотности материала: наиболее плотным фракциям соответствует белый цвет, пустотам (поры, трещины и т.д.) – черный цвет.

Сигналы акустической эмиссии (AE) регистрировалась двумя пьезопреобразователями, прикрепленными к торцам образца. Каждый сигнал акустической эмиссии характеризовался временем излучения, координатой источника по высоте образца и энергией. Точность определения координат источника сигнала AE составляла около 2 мм. Постановка экспериментов и результаты подробно описаны ранее в работах [Damaskinskaya et al., 2018; 2021; 2022].

На рис. 7 и рис. 8 приведены примеры распределений по координате (высоте образца) суммарной энергии сигналов акустической эмиссии. Распределения построены для последовательных временных отрезков. Данное рассмотрение позволяет проследить изменения распределений энергии по координате в процессе нагружения. Предполагается, что наиболее интенсивное дефектообразование происходит в области, в которой наблюдается наибольшее энерговыделение. Этот факт был проверен экспериментально в работах [Damaskinskaya et al., 2021; 2022].

Рис. 7.

Изменение суммарной энергии сигналов AE по координате (высоте образца) при деформировании образца песчаника Berea.

Рис. 8.

Изменение суммарной энергии сигналов AE по координате (высоте образца) при деформировании образца гранита Westerly.

На рис. 7 показаны такие зависимости, полученные при деформировании образца песчаника Berea. Видно, что на каждом временном отрезке наблюдается одна область повышенного энерговыделения. Это указывает на то, что имеется один очаг, координаты которого незначительно изменяются в процессе деформирования.

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

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

ЗАКЛЮЧЕНИЕ

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

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

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

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