CC BY
48
УДК 51-72:539.4.013:539.411.4
РАЗВИТИЕ ФОРМАЛИЗМА МЕТОДА ЧАСТИЦ ДЛЯ МОДЕЛИРОВАНИЯ ОТКЛИКА ФЛЮИДОНАСЫЩЕННЫХ ПОРИСТЫХ ГЕОЛОГИЧЕСКИХ МАТЕРИАЛОВ
Димаки Андрей Викторович,
канд. техн. наук, науч. сотр. лаборатории компьютерного конструирования материалов Института физики прочности и материаловедения СО РАН, Россия, 634021, г. Томск, пр. Академический, 2/4. E-mail: dav@ispms.tsc.ru
Шилько Евгений Викторович,
д-р физ.-мат. наук, ведущ. науч. сотр. лаборатории компьютерного конструирования материалов Института физики прочности и материаловедения СО РАН, Россия, 634021, г. Томск, пр. Академический, 2/4. E-mail: shilko@ispms.tsc.ru
Астафуров Сергей Владимирович,
канд. физ.-мат. наук, науч. сотр. лаборатории компьютерного конструирования материалов Института физики прочности и материаловедения СО РАН, Россия, 634021, г. Томск, пр. Академический, 2/4. E-mail: svastafurov@gmail.com
Псахье Сергей Григорьевич,
член-корреспондент РАН, директор Института физики прочности и материаловедения СО РАН, Россия, 634021, г. Томск, пр. Академический, 2/4; зав. кафедрой физики высоких технологий в машиностроении Института физики высоких технологий ТПУ Россия, 634050, г. Томск, пр. Ленина, д. 30. E-mail: sp@ispms.tsc.ru
Проведено обобщение математического формализма гибридного метода частиц (метода гибридных клеточных автоматов) для численного моделирования процессов деформирования и разрушения флюидонасыщенных пористых геологических материалов на мезо- и макроскопическом масштабных уровнях. Развитый формализм учитывает возможность неупругого деформирования, дилатансии и разрушения твердофазного скелета, а также влияние порового давления на напряженное состояние скелета и перераспределение флюида в трещинно-поровом пространстве геосреды. Сформулирована связанная модель газонасыщенного геоматериала, в рамках которой механический отклик скелета описывается на основе модели пластичности геосред с неассоциированным законом течения, а перераспределение газовой фазы моделируется путем численного решения уравнений фильтрации и диффузии. Адекватность развитого формализма иллюстрируется результатами теоретического изучения влияния порового давления газа на прочность образцов угля в условиях сжатия. Показано, что степень влияния порового давления на интегральные (макроскопические) прочностные характеристики газонасыщенных материалов в значительной степени определяется соотношением параметров критерия прочности, отражающим содержание в твердофазном скелете дефектов и повреждений высокого ранга.
Ключевые слова:
Трещиновато-пористая геологическая среда, флюид, численное моделирование, метод клеточных автоматов, метод частиц, пластичность, дилатансия, поровое давление, разрушение, прочность.
Введение
Известно, что наличие флюидов в трещинно-по-ровом пространстве геологических материалов оказывает существенное влияние на их механический отклик, включая разрушение. В частности, наличие адсорбированных газов может служить причиной внезапных выбросов газопылевой смеси, эпизодически происходящих в ходе выработки угольных пластов [1-6]. При этом, как свидетельствуют результаты многочисленных экспериментальных исследований, степень влияния флюида на характер деформирования и разрушения геоматериалов может быть существенно различной и определяется целым комплексом факторов. К ним относятся как
факторы, непосредственно связанные с пространственным распределением и состоянием флюида в трещинно-поровом пространстве материала (величина открытой/связанной и закрытой пористости, особенности геометрии и пространственного распределения несплошностей, адсорбция на стенках пор и т. д.), так и механические (включая реологические) характеристики самого скелета горной породы (предел упругости, коэффициент деформационного упрочнения, коэффициенты внутреннего трения и дилатансии и другие).
Корректное описание деформирования и разрушения заполненных флюидом пород невозможно без учета процессов перераспределения флюида в
трещинно-поровом пространстве материала. В свою очередь, протекание процессов перераспределения флюида взаимосвязано с деформированием твердофазного скелета, а также с изменением объема и просветности связанного трещинно-поро-вого пространства вследствие дилатансии либо уплотнения породы. Таким образом, теоретическое описание механического отклика флюидонасыщенных трещиновато-пористых сред должно осуществляться в связанной постановке, учитывающей взаимовлияние процессов деформирования твердофазного скелета (приводящего к изменению трещинно-порового объема) и пространственного перераспределения флюида вследствие неоднородного изменения порового давления.
В связи со сказанным актуальной задачей механики деформируемого твердого тела является развитие методов численного моделирования процессов деформирования и разрушения флюидонасыщенных геоматериалов и геосред на различных масштабных уровнях. Перспективным численным методом, созданным для решения задач такого рода, является «гибридный» метод клеточных автоматов (ГКА). Метод ГКА объединяет формализмы метода частиц в формулировке дискретных элементов (МДЭ) для описания механического отклика твердофазного скелета и метода классических клеточных автоматов для моделирования процессов фильтрации и диффузии флюида в трещинно-поро-вом пространстве [7]. Отметим, что преимуществом использования метода частиц при решении связанных задач является возможность непосредственного моделирования явлений разрушения, в том числе множественного, сопровождающегося эффектами перемешивания и локализованного массопереноса.
В настоящей работе предложена и реализована в рамках метода ГКА связанная модель газонасыщенного трещиновато-пористого материала. Данная модель эффективно учитывает взаимовлияние процессов деградации твердофазного скелета при нагружении и перераспределения флюида. Показана адекватность развитого формализма вычислительного метода ГКА на примере изучения влияния порового давления газа на прочность газонасыщенных углей.
Формализм метода гибридных клеточных автоматов
для моделирования флюидонасыщенных
трещиновато-пористых материалов
Основные положения метода
В рамках метода ГКА исследуемая среда рассматривается как суперпозиция двух взаимопроникающих слоев (рис. 1, а), один из которых представлен ансамблем взаимодействующих подвижных клеточных автоматов (ПКА) [8, 9], а другой -сеткой классических клеточных автоматов (ККА) [10]. В такой комбинированной схеме ансамбль подвижных клеточных автоматов моделирует твердофазный скелет, а сетка классических клеточных автоматов используется для моделирования процессов перераспределения флюида в объеме вмещающего твердого тела (иными словами, в объеме подвижных клеточных автоматов) [7].
В рамках метода ГКА шаг интегрирования уравнений динамики моделируемой многофазной системы состоит из двух этапов. Первый из них условно называется «механическим» и представляет собой шаг интегрирования уравнений движения подвижных клеточных автоматов. На этом этапе моделируются процессы деформирования и разрушения материала твердого каркаса под действием механических нагрузок. Второй этап связан с описанием эволюции (переноса) газовой/жидкой фаз в трещинно-поровом пространстве подвижных клеточных автоматов и окружающей моделируемое твердое тело среде (эта среда также может быть заполнена жидкостью или газом). Данные процессы моделируются на сетке классических клеточных автоматов. На основе рассчитанных значений концентраций флюида в трещинно-поровом пространстве подвижных клеточных автоматов вычисляется вклад порового давления в величину локальных внутренних напряжений в твердофазном каркасе подвижных клеточных автоматов. Таким образом, метод ГКА позволяет связать решения динамических задач механики твердого тела и газа или жидкости и, тем самым, описывать поведение многофазной гетерогенной среды. Отметим, что точность решения такой связанной задачи зависит от соотношения
а
б
в
Рис. 1. а) слои подвижных и классических клеточных автоматов; б) проекция пары подвижных клеточных автоматов на сетку классических клеточных автоматов; в) определение сорта классического клеточного автомата
между размерами классического и подвижного автоматов (с его уменьшением повышается точность «проецирования» свойств со слоя ПКА на слой ККА и обратно).
Формализм метода ГКА обладает достаточно высокой «гибкостью», позволяющей моделировать отклик многофазных систем на различных структурных уровнях и реализовывать различные (в том числе многомасштабные) реологические модели твердофазного скелета и модели динамики флюидов. Ниже предлагается реализация в рамках метода ГКА связанной модели трещиноватопористых газонасыщенных геосред на мезо- и макроскопическом масштабах. В ее рамках отклик упругохрупкого твердофазного скелета описывается с использованием модели пластичности Николаевского, перераспределение газовой фазы моделируется путем численного решения уравнений фильтрации и диффузии газа.
Модель механического отклика
твердофазного скелета
Метод подвижных клеточных автоматов, использующийся для моделирования деформирования и разрушения твердофазного каркаса, является представителем метода частиц в формулировке дискретных элементов [8, 9]. Основным преимуществом метода ПКА в сравнении с традиционными реализациями МДЭ [11, 12] является обобщенная (многочастичная) форма записи соотношений для сил центрального взаимодействия в парах элементов, аналогичная записи потенциала погруженного атома. Как показано в работе [9], это позволяет корректно реализовывать различные реологические модели твердых тел.
Для моделирования отклика трещиновато-пористых геоматериалов в рамках метода ПКА реализована модель пластичности горных пород с неассоциированным законом течения и критерием Мизеса-Шлейхера достижения предела упругости (модель Николаевского) [13, 14]. Выбор данной модели связан с тем, что она адекватно описывает отклик геоматериалов на различных масштабных уровнях с учетом вкладов «нижележащих» структурных масштабов. Особенностью модели Николаевского является постулируемая линейная связь между скоростями объемной и сдвиговой составляющих пластической деформации с коэффициентом пропорциональности Л, называемым коэффициентом дилатансии.
Адаптация модели Николаевского к методу частиц осуществлена с использованием алгоритма Уилкинса [15], в рамках которого решение упругопластической задачи осуществляется путем решения упругой задачи и последующей корректировки потенциальных сил взаимодействия частиц с соблюдением необходимых требований модели Николаевского, предъявляемых к величине локального давления и критерия Мизеса-Шлейхера.
Решение упругой задачи на ансамбле частиц сводится к определению приращений центральной
и тангенциальной компонент силы упругого отклика подвижного клеточного автомата Ь на воздействие со стороны соседнего автомата у. Соответствующие выражения записываются на основе обобщенного закона Гука в гипоупругой форме [9]:
;г/х, =
= 20,
1Ае'<, >“КЛ1 - К)
К, >” ‘/в, = 20 Аг„
|АстТ
(1)
;( Л ;(, > / V ' ' 'О >
где символ А означает приращение величины соответствующего параметра за шаг по времени А1 численной схемы интегрирования; ат и тад - удельные значения центральной (¥су,г) и тангенциальной (^“0 составляющих потенциальной силы реакции подвижного клеточного автомата Ь на воздействие со стороны соседа у - площадь поверхности контакта пары; ^ и К - модули сдвига и всестороннего сжатия материала автомата Ь; Ает и Аут - приращения нормальной и сдвиговой деформации автомата Ь в паре Ь-у, а“““ - среднее напряжение в объеме автомата Ь; Р?ш - вклад порового давления флюида в величину среднего напряжения в объеме подвижного клеточного автомата. Среднее напряжение а“““ вычисляется с использованием соотношения Лява, связывающего компоненты тензора усредненных напряжений (о^) в некотором замкнутом объеме с силами, действующими на поверхность данного объема [9, 12].
Известно, что напряженное состояние твердофазного каркаса является сложным и обусловлено не только величиной удельного объема несплошно-стей, но и особенностями их геометрии и пространственного распределения [16, 17]. В случае отсутствия ярко выраженной геометрии и ориентации (текстуры) трещинно-порового пространства адекватным является приближение, в соответствии с которым на мезо- или макроскопическом масштабном уровне учитывается только вклад давления газов в величину гидростатического напряжения в скелете (гидростатическое растяжение). Это приближение и использовано в соотношении (1).
После решения упругой задачи для автомата Ь проверяется условие достижения предельного состояния (критерий Мизеса-Шлейхера), учитывающее вклад порового давления флюида:
ф, = р (а,теа" + Р‘ш > + 5;.,П773 > X, (2)
где У1 - предел упругости материала автомата Ь при сдвиге, р - коэффициент внутреннего трения, а“ - интенсивность усредненных напряжений в объеме автомата Ь.
Приведение сил межчастичного взаимодействия к предельной поверхности (релаксация) осуществляется в соответствии с выражениями, полученными в работе [9]:
Н,, > = К, > - ёТ >М + ёГ - N > Ь'(, > = т,.( ,>М,
где (c'ilj),z'ilj)) - приведенные значения удельных сил реакции автомата; Mr1-(^/-r%3G№rY)/(KA3+3G)) -коэффициент приведения девиаторов напряжений; N=KA (Ф-Yi)/(KAify+Gi) - поправка к величине среднего напряжения, вычисленного в упругом приближении; A - коэффициент дилатансии материала автомата i.
Развитый подход к построению многочастичного взаимодействия дает возможность использовать многопараметрические критерии разрушения (Кулона-Мора, Друккера-Прагера и т. д.) в качестве критерия разрыва связи (потери сцепления) в парах взаимодействующих элементов [9]. В развитой модели применяется модифицированный критерий разрушения Друккера-Прагера, учитывающий вклад локального порового давления флюида:
a dp = 0,5(a + 1)ст1п4 +1,5(a - l)(o-mean + Pas ) > cc, (3)
где ac - пороговое значение для рассматриваемой пары частиц (величина, характеризующая прочность химической связи), a=cjct - отношение прочностей сцепления пары при одноосном сжатии (cc) и растяжении (ct). Отметим, что, хотя критерии наступления неупругой деформации (2) и разрушения (3) имеют одинаковую структурную форму, смысл коэффициентов при слагаемом cmean в них совершенно различный. Так, величина коэффициента 3 определяет степень влияния давления на зарождение и развитие сети микроскопических повреждений и трещин в объеме материала. В то же время коэффициент 1,5(a-1) определяет влияние давления на условие локализованного объединения отдельных микроповреждений в единый мезо- или макроразрыв. Сказанное свидетельствует о том, что значения этих двух коэффициентов не обязательно являются согласованными, поскольку роль локального давления в эволюции повреждения, как правило, существенно возрастает с ростом пространственного масштаба повреждения.
Фильтрационно-диффузионная модель переноса
газовой фазы
При моделировании процессов переноса газов в твердофазном скелете необходимо учитывать, что в общем случае такие системы представляют собой среду с многомасштабной структурой. В порядке убывания размеров можно выделить следующие структурные уровни: среда в целом как сплошное тело (макроуровень); фрагменты и блоки; открытые поры и фильтрационные каналы различной протяженности и диаметров; повреждения и микротрещины; закрытые поры; отдельные кристаллиты [18-20]. Очевидно, что при описании эволюции фрагментов контрастной среды мезо- и макроскопического масштаба (то есть, масштабного уровня, характерный размер которого на порядки величины превышает среднее расстояние между порами и каналами) характеристики трещинно-поровой структуры более низких масштабов должны учитываться интегрально. Среди них прин-
ципиально важными являются следующие: 1) величина открытой пористости (или фильтрационного объема), образованной сообщающейся системой пор, каналов и трещин, имеющих выход на внешнюю поверхность; 2) величина закрытой пористости, представленной ансамблем изолированных микропор и микроповреждений, не имеющих связи с фильтрационным объемом материала.
Математический формализм, применяемый для описания переноса газа в такой контрастной среде, определяется наличием открытой/закрытой пористости и состоянием заключенного в объеме материала газа. В настоящей модели для каждого компонента смеси газов используются приближения идеального газа и изотермических условий (Г=со^) [7]. Перенос газовой смеси в материале описывается системой уравнений, включающей уравнение фильтрации [20, 21]:
To,
dp
д_( kP dP
dt dx( /л dx
(4)
и уравнение диффузии газа из закрытых пор (и твердофазного каркаса) в фильтрационный объем (и обратно) [20]:
др д2р
T closed
dt
■ = Df
dx2
(5)
где р - плотность газа, k - коэффициент проницаемости пористого каркаса, л - динамическая вязкость газа, P - парциальное давление газа в тре-щинно-поровом пространстве, yoipa - открытая пористость, yd()Sed - закрытая пористость, Deff - эффективный коэффициент диффузии, зависящий от коэффициента диффузии материала твердого каркаса, параметров пористости, среднего размера монолитных блоков среды, а также от растворимости газа в твердофазном каркасе [20, 22].
Связь между концентрацией газа в твердофазном каркасе (в том числе в закрытых микропорах) и его плотностью в фильтрационном объеме устанавливается законом Генри, справедливым в широком интервале значений гидростатического давления [23]:
7closedP=vP, V<< !• (6)
где v - растворимость газа в материале твердого каркаса. В качестве уравнения состояния газа в модели используется уравнение Менделеева-Клапейрона:
PV = NRgasT, (7)
где N - количество молей газа; Rgas - универсальная газовая постоянная, T - температура газа (в рамках используемого приближения температура каркаса и газа полагаются постоянными и равными между собой). Отметим, что в общем случае давление газа в фильтрационном объеме и закрытых (изолированных) порах и трещинах является различным.
Поскольку компоненты газовой смеси (отдельные газы) рассматриваются как идеальные, их
парциальные давления могут рассчитываться независимо. Коэффициент проницаемости твердого каркаса удовлетворительно описывается формулой [20]:
k = yd2,
где й - характерный диаметр каналов просачивания.
Уравнения фильтрации (4) и диффузии (5) решаются на сетке классических клеточных автоматов с использованием неявной разностной схемы расщепления по пространственным переменным с соответствующими начальными и граничными условиями, конкретный вид которых определяются физической постановкой решаемой задачи. Результатом их решения на каждом временном шаге являются новые значения давления газа в классическом клеточном автомате.
В случае присутствия в моделируемой системе смеси нескольких газов уравнения (4)-(7) решаются независимо для каждого компонента. Полное давление смеси газов в фильтрационном объеме или закрытых порах классического клеточного автомата представляет собой сумму парциальных давлений.
Следует отметить, что в рамках описанной модели некоторые свойства газов учтены лишь приближенно. В частности, применяется уравнение состояния идеального газа (в порах и трещинах на-носкопического масштаба оно может существенно отличаться от идеального и характеризуется сильным влиянием адсорбированных на поверхности пор молекулярных слоев газа и наличием квантовых эффектов [19, 24]). Кроме того, не учитывается возможность физико-химического взаимодействия различных компонентов в многокомпонентной газовой смеси (в этом случае используемое приближение независимого определения характеристик компонентов является неточным). Тем не менее, потенциал развитого подхода обеспечивает широкие возможности реализации в рамках метода ГКА различных уравнений состояния флюидов и моделей переноса многокомпонентных смесей газов и жидкостей в проницаемых трещиноватопористых системах [25].
Взаимосвязь моделей вмещающего
твердофазного скелета и газовой фазы
Математические модели механического отклика твердофазного скелета и заполняющего его нес-плошности флюида образуют взаимосвязанную систему, моделируемую посредством наложения двух слоев клеточных автоматов (ансамбля ПКА и неподвижной сетки ККА). Влияние слоя классических клеточных автоматов на ансамбль частиц (подвижных автоматов) осуществляется посредством вклада Рв соотношениях (1), (2) и (3). В настоящей модели вычисление давления газа осуществляется в приближении линейной зависимости от величины удельного объема фильтрационного и изолированного трещинно-порового пространства:
п gas r¡closed closed . -rppen open
Pi =P Y +P h • (8)
Здесь Yi°p‘" и Yilosed - фильтрационная и закрытая пористость в объеме подвижного клеточного автомата i, P°p’" и P^ed - полное давление газовых смесей, содержащихся в фильтрационном объеме и изолированных порах/микротрещинах подвижного автомата. Давления P¡°р’а и Pf°”d определяются усреднением соответствующих значений по классическим клеточным автоматам, относящимся к подвижному клеточному автомату i (рис. 1, б). Выражение (8) следует рассматривать как первое приближение, предполагающее однородное распределение пор, трещин и фильтрационных каналов в объеме подвижного клеточного автомата и не учитывающее особенности их геометрии и пространственной ориентации. Кроме того, развитая модель не учитывает влияния молекул газа, абсорбированных в кристаллической решетке, на величину упругой энергии твердофазного каркаса. Тем не менее, как показано в [23], такое приближение является вполне допустимым.
Отметим также, что наличие газов в окружающем моделируемое твердое тело пространстве (в «макропорах», рис. 1, в) приводит к возникновению дополнительных сил, действующих на поверхность твердого тела и вызывающих внутренние напряжения в скелете. В рассматриваемой модели влияние давления «внешних» газов на напряженное состояние подвижных клеточных автоматов, моделирующих поверхностный слой твердофазного каркаса, вычисляется в соответствии с ранее полученными соотношениями [7]. Полагается, что давление окружающего газа на свободную (не занятую соседями) часть поверхности подвижного клеточного автомата i дает вклад в полную силу, приложенную к центру масс автомата, а также приводит к изменению величины компонентов усредненных напряжений а^.
Изучение влияния порового давления газа
на распределение напряжений
и характер разрушения образцов угля
Корректность развитого формализма метода гибридных клеточных автоматов подтверждена тестовыми расчетами фильтрации и диффузии газов в газоносных геоматериалах [7], а также моделированием механического отклика таких систем при механическом нагружении. В качестве примера ниже приведены результаты моделирования одноосного нагружения образца молодого бурого угля (лигнита) [26], трещинно-поровое пространство которого заполнено газом (использовался режим нагружения с постоянной скоростью 7=0,02 м/с, рис. 2, а). Давление газа в окружающем образец пространстве полагалось равным 0 бар. Образец помещен между пуансонами, нижний из которых оставался неподвижным, а верхний перемещался вниз с постоянной скоростью 7. Механические свойства материала образца соответствовали высококачественной литотипной ка-
Пуансон (Уу = -0.02 м/сек)
а
Деформация
б
Рис. 2. Схема моделируемого образца (а), линеаризованная диаграмма одноосного сжатия дробного детрита (б)
тегории лигнита (дробному детриту), упрощенная (линеаризованная) диаграмма одноосного сжатия которого приведена на рис. 2, б. Коэффициент ди-латансии материала полагался равным Л=0,173, коэффициент внутреннего трения ^=0,49. Макроструктура моделируемого образца предполагалась однородной, макроскопические поры и включения отсутствовали, что является типичным для дробного детрита. Параметры пористости лигнита, (таблица) были получены на основании экспериментальных данных по сорбционным свойствам лигнита, полученных Й. Пездичем и Я. Жулой на образцах из Веленской угольной шахты (Веленье, Словения) [27, 28].
Таблица. Оценки параметров пористости и сорбционных свойств лигнита (Веленская угольная шахта, Веленье, Словения)
Открытая пористость 0,042
Коэффициент диффузии, см2/с 3,4-10-7
Закрытая пористость 0,55
Коэффициент растворимости 0,05
Средний диаметр фильтрационного канала, м 2,5-10-9
Средний размер «монолитного» блока, м 1-10-4
Результаты моделирования показали, что в условиях одноосного сжатия образцы дробного детрита демонстрируют хрупкий отклик, хорошо согласующийся с поведением реальных образцов цилиндрической формы. Разрушение образцов происходит путем формирования одной или нескольких диагональных трещин. Расчеты, проводимые при различных значениях порового давления газа в интервале от 0 до 5,5 бар, показали, что прочность образцов снижается с ростом порового давления по линейному закону (рис. 3, а). Это соответствует экспериментальным данным для различ-
ных углей и свидетельствует об адекватной формулировке критерия разрушения (3) и корректной реализации связанной модели газонасыщенных геоматериалов.
Хорошо известно, что хрупкие материалы, в том числе большинство горных пород, характеризуются низкими значениями прочности на растяжение (а() при достаточно высоких значениях прочности на сжатие (а). Отношение а/а=а может достигать 10 и более и определяется особенностями гетерогенной структуры материала, в том числе иерархией исходных повреждений. Так, с увеличением концентрации «крупных» дефектов и повреждений (много больших характерных размеров пор и каналов фильтрационного объема) величина прочности материала на растяжение снижается намного быстрее величины прочности на сжатие. Это связано с различиями в функционировании повреждений плоскостного типа в условиях растяжения (разрыв) и сжатия (сухое/вязкое трение сжатых поверхностей повреждения). Отметим, что наличие повреждений достаточно высокого ранга может приводить к кратному уменьшению величины а, даже если их концентрация (и вклад в фильтрационную пористость) незначительна. В связи со сказанным актуальным представляется анализ степени влияния порового давления газа на прочность трещинно-пористых геоматериалов, характеризующихся различным отношением а/ а=а. В настоящей работе такое теоретическое исследование проведено для образцов дробного детрита (рис. 2, а) при одноосном сжатии.
Рассматривались образцы дробного детрита с различным соотношением параметров ас и а1 критерия разрушения Друккера-Прагера (3). При этом значение в проводимых расчетах полагалось постоянным (12,5 МПа). Значение второго па-
раметра (а) варьировалось в пределах от 4 до 10 МПа. На рис. 3, а приведены зависимости прочностей моделируемых образцов с различными значениями а от величины вклада порового давления газа в величину среднего напряжения в твердофазном скелете (Рт), определяемого выражением (8). Приведенные кривые имеют вид, близкий к линейному. При этом влияние внутрипорового газа на прочность образца уменьшается с ростом отношения а/а и в пределе стремится к нулю. Это проявляется, в частности, в нелинейном характере снижения тангенса угла наклона кривых, приведенных на рис. 3, а (рис. 3, б).
Результаты численного моделирования, представленные на рис. 3, б, могут быть верифицированы аналитической оценкой. Применяемый в рамках изложенной модели модифицированный критерий Друккера-Прагера (3) может быть переписан в следующей форме:
:0,5( а +1) <7ш, +1,5( а -1)( Кєш
+ )'
<*Р„п
- = 1,5(а -1) (Кгш + Р^) =
а gаs
= 1,5(а -1) = 1,5 ^^ -1 ^.
(9)
В рамках используемого приближения поровое давление газа не оказывает интегрального влияния на касательные напряжения во вмещающем твердом теле. Поэтому, пренебрегая влиянием объемной деформации на величину Рт, можно записать:
сте с результатами численного моделирования. Можно видеть как качественное совпадение характера построенных зависимостей, так и хорошее количественное согласие.
Различие между аналитической зависимостью
а ® ЭР / ч
—-----(а) и результатами численного моделирова-
ар
gаs
ния составляет менее 20 % и приблизительно одинаково во всем рассмотренном интервале прочностей. Данное различие обусловлено главным образом эффектами локализации неупругой деформации и дилатансионного расширения в образцах моделируемого газонасыщенного материала при сжатии. Действительно, аналогичные расчеты, проведенные на образцах с упругохрупкими механическими характеристиками (линейно-упругий отклик вплоть до разрушения) и теми же значениями прочностных параметров, показали более чем двукратное уменьшение различия между аналитическими и численными оценками зависимости
а&п
ара
-(а).
(10)
Выражение (10) является аналитической оценкой скорости снижения прочности образцов с увеличением порового давления газа. Аналитическая
^ а 0р , ч .
зависимость —р— (а) приведена на рис. 3, б вме-
Заключение
В работе проведено развитие формализма метода частиц для моделирования связанных процессов деформации и разрушения трещиновато-пористых материалов и перераспределения флюидов в их трещинно-поровом пространстве. Механический отклик твердофазного скелета описывается на основе модели пластичности горных пород с неассоциированным законом течения и критерием Мизеса-Шлейхера достижения предела упругости (модель Николаевского). Перенос газа (или смеси газов) в трещинно-поровом пространстве материала описывается уравнениями фильтрации и диффузии. Взаимосвязь двух систем уравнений осу-
Давление газа, МПа а
Прочность материала на растяжение, МПа б
Рис. 3.
Зависимость прочности образца модельного газонасыщенного материала от давления газа в порах при различных значениях прочности материала на растяжение (а), зависимость угла наклона кривых (рис. а) от значения прочности материала на растяжение (б)
ществляется посредством учета влияния порового давления газов на пределы упругости и прочности твердофазного скелета. Предложенная связанная модель реализована в рамках «гибридного» метода клеточных автоматов, объединяющего вычислительные методы частиц и классических клеточных автоматов.
Адекватность развитого формализма продемонстрирована на примере моделирования одноосного сжатия газонасыщенных образцов бурого угля. Показано, что с увеличением давления газа в поро-вом пространстве прочность образцов при одноосном сжатии линейно снижается. Вид и наклон данной зависимости соответствуют экспериментальным данным и теоретическим представлениям о механизмах влияния порового давления газа на прочностные свойства геоматериалов. Полученные результаты подтвердили применимость
СПИСОК ЛИТЕРАТУРЫ
1. Петухов И.М., Линьков А.М. Механика горных ударов и выбросов. - М.: Наука, 1983. - 280 с.
2. Полевщиков Г.Я. Динамические газопроявления при проведении подготовительных и вскрывающих выработок в угольных шахтах. - Кемерово: Институт угля и углехимии СО РАН, 2003. - 317 с.
3. Flores R.M. Coalbed methane: from hazard to resource // Int. J. Coal Geol. - 1998. - V. 35. - P. 3-26.
4. Beamish B.B., Crosdale P.J. Instantaneous outbursts in underground coal mines: an overview and association with coal type // Int. J. Coal Geol. - 1998. - V. 35. - P. 27-55.
5. Cao Y., He D., Glick D.C. Coal and gas outbursts in footwalls of reverse faults // Int. J. Coal Geol. - 2001. - V. 48. - P. 47-63.
6. Larsen J.W. The effects of dissolved CO2 on coal structure and properties // Int. J. Coal Geol. - 2004. - V. 57. - P. 63-70.
7. Псахье С.Г., Шилько Е.В., Смолин А.Ю. и др. Развитие подхода к моделированию деформирования и разрушения иерархически организованных гетерогенных, в том числе контрастных, сред // Физическая мезомеханика. - 2011. - Т. 14. -№ 3. - С. 27-54.
8. Псахье С.Г., Хори Я., Коростелев С.Ю. и др. Метод подвижных клеточных автоматов как инструмент для моделирования в рамках физической мезомеханики // Известия вузов. Физика. - 1995. - Т. 38. - № 11. - С. 58-69.
9. Psakhie S., Shilko E., Smolin A. et al. Development of a formalism of movable cellular automaton method for numerical modeling of fracture of heterogeneous elastic-plastic materials // Fracture and Structural Integrity. - 2013. - V. 24. - P. 59-91.
10. Лоскутов А.Ю., Михайлов А.С. Введение в синергетику. - М.: Наука, 1990. - 272 с.
11. Jing L., Stephansson O. Fundamentals of discrete element method for rock engineering: theory and applications. - London: Elsevier, 2007. - 562 p.
12. Potyondy D.O., Cundall P.A. A bonded-particle model for rock // Int. J. Rock Mech. Min. Sci. - 2004. - V. 41. - P. 1329-1364.
13. Гарагаш И.А., Николаевский В.Н. Неассоциированные законы течения и локализации пластической деформации // Успехи механики. - 1989. - Т. 12. - № 1. - С. 131-183.
14. Стефанов Ю.П. Локализация деформации и разрушение в геоматериалах. Численное моделирование // Физическая мезомеханика. - 2002. - Т. 5. - № 5. - С. 107-118.
метода ГКА для изучения механического отклика (включая разрушение) газоносных геоматериалов в сложных условиях нагружения.
Формализм метода гибридных клеточных автоматов является достаточно гибким и обеспечивает возможность реализации в его рамках различных связанных моделей контрастных сред, содержащих компоненты в твердой, жидкой и газообразной фазах. Это дает потенциальную возможность исследовать большой комплекс проблем, связанных с влиянием флюидов на состояние геоматериалов и участков геосреды (например, горного массива в окрестности выработки или разломных зон) вблизи критического состояния.
Работа выполнена в рамках проекта 111.23.1.4 Программы фундаментальных исследований СО РАН (2013-2016 гг.) и при поддержке проекта № 4 Программы Президиума РАН № 25.
15. Wilkins M.L. Computer simulation of dynamic phenomena. -Heidelberg: Springer-Verlag, 1999. - 246 p.
16. Kushch V.I., Shmegera S.V., Sevostianov I. SIF statistics in micro cracked solid: effect of crack density, orientation and clustering // Int. J. Eng. Sci. - 2009. - V. 47. - P. 192-208.
17. Черных В.А. Модель изменения проницаемости неоднородной горной породы в процессе деформации // ФТПРПИ. - 1990. -№ 1. - С. 126-128.
18. Cushman J.H. The physics of fluids in hierarchical porous media: angstroms to miles. - Heidelberg: Springer, 1997. - 469 p.
19. Хавкин А.Я. Наноявления и нанотехнологии в добыче нефти и газа. - М.; Ижевск: НИЦ «Регулярная и хаотическая динамика», 2010. - 692 с.
20. Алексеев А.Д., Василенко Т.А., Гуменник К.В. и др. Диффузионно-фильтрационная модель выхода метана из угольного пласта // ЖТФ. - 2007. - Т. 77. - № 4. - С. 65-74.
21. Христианович С.А. Об основах теории фильтрации // ФТПРПИ. - 1989. - № 5. - С. 3-18.
22. Эттингер И.Л., Ковалева И.Б. Напряжения набухания и свободная энергия в системе газ-уголь // ДАН СССР. - 1979. -Т. 244. - № 3. - С. 659-663.
23. Алексеев А.Д., Синолицкий В.В. Теория абсорбции газа пористыми веществами // Физика и техника высоких давлений. -1983. - Вып. 12. - С. 103-106.
24. Гриняев Ю.В., Псахье С.Г. О фазовых переходах газа в 2D на-нопоре // Физическая мезомеханика. - 2013. - Т. 16. - № 4. -С. 105-108.
25. Нигматулин Р.И. Динамика многофазных сред. Ч. II. - М.: Наука, 1987. - 360 с.
26. Markic M., Sachsenhofer R.F. Petrographic composition and de-positional environments of the Pliocene Velenje lignite seam (Slovenia) // Int. J. Coal Geol. - 1997. - V. 33. - P. 229-254.
27. Pezdic J., Markic M., Letic M. et al. Laboratory simulation of adsorption-desorption processes on different lignite lithotypes from the Velenje lignite mine // RMZ - Materials and Geoenvironment. - 1999. - V. 46. - № 3. - P. 555-568.
28. Zula J., Pezdic J., Zavsek S. et al. Adsorption capacity of the Velenje lignite: methodology and equipment // RMZ - Materials and Geoenvironment. - 2011. - V. 58. - № 2. - P. 193-216.
Поступила 19.11.2013 г.
UDC 51-72:539.4.013:539.411.4
DEVELOPMENT OF PARTICLE METHOD FORMALISM FOR NUMERICAL MODELING OF THE RESPONSE OF FLUID-SATURATED POROUS GEOLOGICAL MATERIALS
Andrey V. Dimaki ,
Cand. Sc., Institute of Strength Physics and Materials Science, Siberian Branch of Russian Academy of Science, Russia, 634021, Tomsk, Akademicheskii avenue, 2/4. E-mail: dav@ispms.tsc.ru
Evgeny V. Shilko,
Dr. Sc., Institute of Strength Physics and Materials Science, Siberian Branch of Russian Academy of Science, Russia, 634021, Tomsk, Akademicheskii avenue, 2/4. E-mail: shilko@ispms.tsc.ru
Sergey V. Astafurov,
Cand. Sc., Institute of Strength Physics and Materials Science, Siberian Branch of Russian Academy of Science, Russia, 634021, Tomsk, Akademicheskii avenue, 2/4. E-mail: svastafurov@gmail.com
Sergey G. Psakhie,
Dr. Sc., Institute of Strength Physics and Materials Science, Tomsk Polytechnic University, Russia, 634050, Tomsk, Lenin Avenue, 30. E-mail: sp@ispms.tsc.ru
The authors have generalized the mathematical formalism of hybrid particle method (method of hybrid cellular automata) for the numerical simulation of deformation and fracture of fluid-saturated porous geological materials at the meso- and macroscopic scales. The formalism allows taking into account inelastic deformation, dilatancy and fracture of solid skeleton as well as the influence of pore pressure on skeleton stress state and the fluid redistribution in pore space of the fractured geological medium. In the framework of the model the mechanical response of the skeleton is described within the model of plasticity of geological medium with unassociated flow law. The redistribution of the gas phase is simulated by numerical solution of equations of filtration and diffusion. The adequacy of the formalism is illustrated with the results of the theoretical study of gas pore pressure influence on coal sample strength under uniaxial compression. It is shown that the degree of pore pressure influence on macroscopic strength of gas-saturated material is strongly determined by the ratio of parameters of strength criterion, which reflects the content of defects and damages of high rank in the solid skeleton.
Key words:
Fractured porous geological medium, fluid, numerical modeling, cellular automaton method, particle method, plasticity, dilatancy, pore pressure, fracture, strength.
REFERENCES
1. Petukhov I.M., Linkov A.M. Mekhanika gornykh udarov i vybro-sov [Mechanics of Rockbursts and Outbursts]. Moscow, Nedra Publ. House, 1983. 280 p.
2. Polevshchikov G.Ya. Dinamicheskie gazoproyavleniya pri prove-denii podgotovitelnykh i vskryvayushchikh vyrabotok v ugolnykh shakhtah [Dynamic gas effects in the preparatory and revealing workings in coal mines]. Kemerovo, Institute of Coal and Coal Chemistry SB RAS, 2003. 317 p.
3. Flores R.M. Coalbed methane: from hazard to resource. Int. J. Coal Geol., 1998, vol. 35. pp. 3-26.
4. Beamish B.B., Crosdale P.J. Instantaneous outbursts in underground coal mines: an overview and association with coal type. Int. J. Coal Geol., 1998, vol. 35, pp. 27-55.
5. Cao Y., He D., Glick D.C. Coal and gas outbursts in footwalls of reverse faults. Int. J. Coal Geol., 2001, vol. 48, pp. 47-63.
6. Larsen J.W. The effects of dissolved CO2 on coal structure and properties. Int. J. Coal Geol., 2004, vol. 57, pp. 63-70.
7. Psakhie S.G., Shilko E.V., Smolin A.Yu. Razvitie podkhoda k mo-delirovaniyu deformirovaniya i razrusheniya ierarkhicheskikh or-ganizovannykh geterogennykh v tom chisel kontrastnykh sred [Approach to simulation of deformation and fracture of hierarchically organized heterogeneous media, including contrast media]. Fi-sicheskaya mezomekhanika, 2011, vol. 14, no. 5-6, pp. 224-248.
8. Psakhie S.G., Horie Y., Korostelev S.Yu. Metod podvizhnykh kle-tochnykh avtomatov kak instrument dlya modelirovaniya v ram-kakh fizicheskoy mezomekhaniki [Method of movable cellular automata as a tool for simulation within the framework of mesomecha-nics]. Izvestiya vuzov. Fizika, 1995, vol. 38, no. 11, pp. 1157-1168.
9. Psakhie S., Shilko E., Smolin A. Development of a formalism of movable cellular automaton method for numerical modeling of fracture of heterogeneous elastic-plastic materials. Fracture and Structural Integrity, 2013, vol. 24, pp. 59-91.
10. Loskutov A.Yu., Mikhailov A.S. Vvedenie v sinergetiku [Foundations of Synergetics]. Moscow, Nauka, 1990, 272 p.
11. Jing L., Stephansson O. Fundamentals of discrete element method for rock engineering: theory and applications. London, Elsevier, 2007. 562 p.
12. Potyondy D.O., Cundall P.A. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci., 2004, vol. 41, pp. 1329-1364.
13. Garagash I.A., Nikolaevskiy V.N. Neassotsiirovannye zakony techeniya i lokalizatsii plasticheskoy deformatsii [Non-associated laws of plastic flow and localization of deformation]. Uspekhi mekhaniki - Advanced Mech., 1989, vol. 12, no. 1, pp. 131-183.
14. Stefanov Yu.P. Lokalizatsiya deformatsii i razrushenie v geoma-terialakh. Chislennoe modelirovanie [Deformation localization and fracture in geomaterials. Numerical simulation]. Fisicheska-ya mezomekhanika, 2002, no. 5-6, pp. 67-77.
15. Wilkins M.L. Computer simulation of dynamic phenomena. Heidelberg, Springer-Verlag, 1999. 246 p.
16. Kushch V.I., Shmegera S.V., Sevostianov I. SIF statistics in micro cracked solid: effect of crack density, orientation and clustering. Int. J. Eng. Sci., 2009, vol. 47, pp. 192-208.
17. Chernykh V.A. Model izmeneniya pronitsaemosti neodnorodnoy gornoy porody v protsesse deformatsii [The model of variation of permeability of inhomogeneous rock during deformation process]. FTPRPI - Journal of Mining Science, 1990, no. 1, pp. 126-128.
18. Cushman J.H. The physics of fluids in hierarchical porous media: angstroms to miles. Heidelberg, Springer, 1997. 469 p.
19. Khavkin A.Ya. Nanoyavleniya i nanotekhnologii v dobyche nefti i gaza [Nano-effects and nano-technologies in oil and gas mining]. Moscow; Izhevsk, HITS «Regular and Chaotic Dynamics», 2010. 692 p.
20. Alexeev A.D., Vasilenko T.A., Gumennik K.V., Kalugina N.A., Feldman E.P. Diffuzionno-filtratsionnaya model vykhoda meta-na iz ugolnogo plasta [Diffusion-filtration model of methane escape from a coal seam]. Zhurnal Tekhnicheskoy Fiziki, 2007, vol. 52, no. 4, pp. 456-465.
21. Khristianovich S.A. Ob osnovakh teorii filtratsii [Fundamentals of seepage theory]. Journal of Mining Science, 1989, vol. 25, no. 5, pp. 397-412.
22. Ettinger I.L., Kovaleva I.B. Napryazheniya nabukhaniya i svo-bodnaya energiya v sisteme gaz-ugol [Swelling stress and the free
energy of the gas-coal system]. Dokladi Akademii Nauk SSSR -Proceedings of the Russian Academy of Science, 1979, vol. 244, no. 3, pp. 659-663.
23. Alekseev A.D., Sinolitskiy V.V. Teoriya absorbtsii gaza poristymi veshchestvami [The theory of gas absorption by porous matters]. Fizika i tekhnika visokikh davleniy - Physics and Techniques of High Pressure, 1983, no. 12, pp. 103-106.
24. Grinyaev Yu.V., Psakhie S.G. O fazovykh perekhodakh gaza v 2D nanopore [On phase transitions of a gas in 2D nano-pore]. Fi-zicheskaya mesomekhanika - Physical Mesomechanics, 2013, vol. 16, no. 4, pp. 105-108.
25. Nigmatulin R.I. Dinamika mnogofaznykh sred [Dynamics of Multiphase Media]. New York, Hemisphere, 1990. Vol. 2, 371 p.
26. Markic M., Sachsenhofer R.F. Petrographic composition and de-positional environments of the Pliocene Velenje lignite seam (Slovenia). Int. J. Coal Geol., 1997, vol. 33, pp. 229-254.
27. Pezdic J., Markic M., Letic M. Laboratory simulation of adsorption-desorption processes on different lignite lithotypes from the Velenje lignite mine. RMZ - Materials and Geoenvironment, 1999, vol. 46, no. 3, pp. 555-568.
28. Zula J., Pezdic J., Zavsek S. Adsorption capacity of the Velenje lignite: methodology and equipment. RMZ - Materials and Geoenvironment, 2011, vol. 58, no. 2, pp. 193-216.
