CC BY
21
УДК 519.63
DOI 10.25205/1818-7900-2018-16-1-61-73
Б. М. Глинский 1 2, С. И. Кабанихин 1 2, С. А. Калинин 3, И. М. Куликов 1 2 Р. З. Курмангалиев 4, Н. С. Новиков 1 2, А. Э. Рязанцев 5, А. Ф. Сапетина 1 2 И. Г. Черных \ М. А. Шишленин 1 2
1 Институт вычислительной математики и математической геофизики СО РАН пр. Академика Лаврентьева, 6, Новосибирск, 630090, Россия
2 Новосибирский государственный университет ул. Пирогова, 1, Новосибирск, 630090, Россия
3АО «Геологика» ул. Мусы Джалиля, 23, Новосибирск, 630055, Россия
4 Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН ул. Институтская, 4/1, Новосибирск, 630090, Россия
5 ООО «Исследовательский комплекс центра технологического обеспечения» ул. Инженерная, 20, Новосибирск, 630090, Россия
gbm@sscc.ru, kabanikhin@sscc.ru, sergey.kalinin@list.ru, kulikov@ssd.sscc.ru rkurmangaliev1@mail.ru, novikov-1989@yandex.ru, Ryazantsev1@gmail.com afsapetina@gmail.com, chernykh@parbz.sscc.ru, mshishlenin@ngs.ru
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ УПРУГОПЛАСТИЧЕСКИХ ДЕФОРМАЦИЙ ДЛЯ ЗАДАЧ ОБРАЗОВАНИЯ И ЭВОЛЮЦИИ ГЕОЛОГИЧЕСКИХ ТРЕЩИН *
Представлена математическая модель упругопластических сред с учетом фазовых переходов между твердым телом, жидкостью и газом. Также учтен процесс откола вещества от газовой фазы. Сформулирован один из подходов к записи уравнения состояния упругопластических сред, выписаны основные инварианты, используемые для построения уравнения состояния. Представлены результаты математического моделирования геологических деформаций, приводящих к образованию трещин.
Ключевые слова: гидроразрыв пласта, упругопластичная среда, математическая модель, метод Годунова, адаптивные и подвижные сетки.
Введение
Задача определения геометрии трещины гидроразрыва пласта (ГРП) является ключевой при разработке дизайна ГРП. При проведении ГРП осуществляется закачка жидкости гидро-
* Работа выполнена в рамках гранта РФФИ № 16-29-15120 и государственного задания ФАНО России № 03152016-009 и 0315-2016-0005.
Глинский Б. М., Кабанихин С. И., Калинин С. А., Куликов И. М., Курмангалиев Р. З., Новиков Н. С., Рязанцев А. Э., Сапетина А. Ф., Черных И. Г., Шишленин М. А. Математическое моделирование упругопластических деформаций для задач образования и эволюции геологических трещин // Вестн. НГУ. Серия: Информационные технологии. 2018. Т. 16, № 1. С. 61-73.
ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2018. Том 16, № 1 © Б. М. Глинский, С. И. Кабанихин, С. А. Калинин, И. М. Куликов, Р. 3. Курмангалиев, Н. С. Новиков, А. Э. Рязанцев, А. Ф. Сапетина, И. Г. Черных, М. А. Шишленин, 2018
разрыва в продуктивный пласт, ограниченный малопроницаемыми пластами. Трещина ГРП по мере закачки жидкости гидроразрыва растет в размерах, при этом наблюдается несколько качественно различных этапов в поведении трещины ГРП. На первом этапе проведения процесса гидравлического разрыва пласта происходит формирование зародыша будущей трещины, для чего в скважину спускается так называемый перфоратор; далее производится отстрел зарядов, формирующих отверстия заданного диаметра в стенке скважины и в прискважинной области. После проведения перфорации производится серия тестов, часто объединяемых общим названием «мини-ГРП». На этапе мини-ГРП в продуктивном пласте создается трещина некоторых размеров и формы путем нагнетания жидкости ГРП в перфорационные отверстия, при этом ведется тщательная запись поведения графика давления (и ряда других параметров) приборами, установленными в соответствующих местах наземного оборудования. В некоторых случаях используются забойные приборы для измерения параметров мини-ГРП (например, давления и температуры). По характеру изменения давления во время проведения мини-ГРП делается оценка параметров пласта (определение коэффициента утечек, минимальное напряжение закрытия трещины и ряд других). Если после некоторого времени закачки жидкости выключить насосы и предоставить систему самой себе, то будет наблюдаться постепенная утечка жидкости в пласт. Этот процесс может проходить с различной интенсивностью, зависящей от добавок, содержащихся в жидкости ГРП, и от проницаемости пласта. После того как трещина мини-ГРП закроется, а все необходимые параметры теста мини-ГРП будут записаны, производится подготовка дизайна основного ГРП с учетом измеренных параметров. Целью дизайна ГРП является формирование так называемого расписания закачки (таблицы с концентрациями закачиваемых агентов, расходами и рядом других контролируемых параметров). По сформированному расписанию начинается закачка основного ГРП с повторным открытием уже сформированной трещины. На последних этапах вместе с жидкостью ГРП производится закачка расклинивающего агента (ргоррап1;), который предотвращает смыкание стенок трещины после остановки закачки. В технологии ГРП используется основное соотношение теории фильтрации - закон Дарси. Данное соотношение связывает расход жидкости, фильтрующейся через проницаемый пласт, с градиентом давления, проницаемостью пород, вязкостью флюидов и площадью притока флюидов. Целью проведения операции ГРП является создание формы и размеров пропированной трещины ГРП, которые позволили бы увеличить площадь притока пластовых флюидов и объем добываемого флюида в единицу времени.
Из приведенного описания технологии ГРП можно понять, что размеры трещины контролируются сразу несколькими параметрами: прочностью (трещиностойкостью) пород, вязкостью жидкости ГРП и скоростью фильтрации жидкости в пласт. При определенной скорости утечек давление, создаваемое жидкостью ГРП, может становиться недостаточным для увеличения размеров трещины. Кроме того, непроницаемые пласты, ограничивающие продуктивный пласт, обладают прочностными характеристиками, отличными от аналогичных характеристик продуктивного пласта. В то же время сжимающие напряжения, действующие в покрышке и подошве продуктивного пласта, как правило, превосходят напряжения в продуктивном пласте: трещина не может неограниченно расти в высоту из-за действия напряжений в покрышке и подошве продуктивного пласта.
Из анализа системы напряжений, действующих на пласт, можно понять, что характер распределения этих напряжений зависит от принятых допущений о свойствах материалов, а также качества экспериментальных данных, описывающих геометрию слагающих пластов (сейсмические исследования). При допущениях об изотропном (слоисто-анизотропном) материале и параллельном горизонтальном залегании пластов постоянной толщины, распределение напряжений будет таким, что трещина ГРП будет располагаться преимущественно в вертикальной плоскости. Однако в случае геологического строения, осложненного наличием пластов с переменной мощностью, разломов и большими углами падения пластов, может формироваться такое напряженное состояние, которое будет приводить к значительному искривлению плоскости трещины.
Решение плоской задачи о деформации прямоугольной области с трещиной, нагруженной давлением на ее берегах, было получено в работах Ю. П. Желтова и С. А. Христиановича [1], Г. И. Баренблатта [2] и др. Последовательный обзор первоначальных решений с ясной систе-
мой обозначений и подробным объяснением физических допущений для рассматриваемой плоской задачи приведен в сборнике А. Ю. Ишлинского [3]. Среди работ последнего времени, которые следует выделить в контексте задачи, решаемой в данной работе, можно отметить решения, полученные в [4]. Нельзя не упомянуть многочисленные решения зарубежных авторов [5], которые были получены независимо от работ отечественных ученых. Общим местом для всех решений является допущение о плоской деформации горизонтального сечения вертикальной трещины. В литературе принято обозначать модели такого класса аббревиатурой KGD по первым буквам фамилий авторов (Khristianovich - Geerstma -de Klerk). Альтернативным взглядом на решение плоской задачи о росте трещины ГРП является класс моделей, обозначаемых аббревиатурой PKN (Perkins - Kern - Nordgren). В данном классе моделей рассматривается вертикальное сечение вертикальной же трещины ГРП. Интересной особенностью данного класса моделей является допущение о том, что трещина останавливает свой рост в высоту, когда подвижный фронт трещины достигает покрышки и подошвы продуктивного пласта, в котором и формируется трещина ГРП. Дальнейший рост объема трещины типа PKN происходит за счет роста в длину. Более сложные модели (P3D, Planar-3D) предполагают, что фронт трещины может проникать на некоторую глубину в ограничивающие пласты. Целью данной работы является решение плоской задачи о росте трещины типа KGD, однако некоторые выводы являются более общими, чем это предполагается в исходной модели KGD. Следует отметить, что в первоначальных моделях использовался целый ряд очень сильных допущений: трещина остается плоской все время процесса ГРП, утечки жидкости ГРП в пласт отсутствуют, пласт считается однородным и изотропным, напряжения на границах прямоугольной области являются однородными и т. п. Допущение о том, что трещина остается плоской, прямо следует из условий симметрии задачи. Общей особенностью всех упомянутых решений является использование так называемых интегральных уравнений. В частности, решение о перемещении точек, расположенных на берегах трещины, получают путем совместного решения уравнения смазки (уравнение Рейнольдса) и интегрального уравнения связи между перемещениями и давлением на стенки трещины. Упомянутое интегральное уравнение получают путем того или иного обобщения классической задачи теории упругости о точечной силе, приложенной к упругому полупространству (задача Буссинеска). Теоретические решения плоской задачи теории упругости для областей с трещинами приведены в работах Мусхелишвили [6], где представлен подробный вывод соответствующих интегральных уравнений с привлечением ТФКП. Поскольку задача о росте трещины ГРП предполагает рассмотрение физических процессов из различных разделов физики (теория трещин, теория фильтрации, гидродинамика и т. д.), то получение даже простейшего аналитического решения является очень полезным в теоретических работах по ГРП. Тем не менее, в работах последнего времени основной упор делается на численные решения задачи о росте трещины ГРП. Наиболее удачное численное решение для трещин типа KGD было получено Жаном Дероше (Jean Desroches) и Марком Терцилином (Marc Thiercelin) после разработки кода Loramec [7]. Многие решения плоской задачи сравнивались именно с решением Loramec. Отметим, что код Loramec позволил смоделировать характер деформации берегов трещины вблизи ее кончика. Вообще говоря, распределение деформаций и давления в кончике трещины ГРП является определяющим для корректной количественной оценки. Более того, известно, что подвижный фронт жидкости ГРП отстает от подвижного фронта растущей трещины, при этом образуется так называемый «fluid lag» - сухой кончик трещины, преимущественно заполненной жидкостью ГРП. На существование лага указывалось ещё в работе [1], однако изучение влияния лага было выполнено в работе [8], а первые численные решения получены с помощью кода Loramec. Отдельно следует упомянуть, что для лабораторной количественной проверки решений Loramec была изготовлена специальная лабораторная установка. Этот факт является особенно интересным, поскольку большинство моделей ГРП не могут быть проверены в условиях лабораторного эксперимента из-за чрезвычайно сложной физики решаемой связанной задачи. В частности, прогиб стенок трещины типа KGD, а также распределение давления по длине трещины подробно рассмотрены в [9; 10]. Изучение этих работ показывает, что распределение прогибов стенок трещины вблизи кончика трещины подчиняется очень сложным физическим закономерно-
стям, поэтому привлекаются различные асимптотические приближения решения либо используется достаточно мелкая расчетная сетка вблизи кончика трещины ГРП.
Теоретический интерес представляет расчет геометрии трещины ГРП, распространяющейся в слоисто-неоднородном пласте при наличии такой системы напряжений на границах прямоугольной области, при которой трещина может искривляться. Кроме того, исследование поведения кончика трещины ГРП численными методами требует привлечения достаточно подробной расчетной сетки, как это было показано в работах [7-10]. Для того чтобы избежать адаптивного перестроения сетки в подвижном кончике трещины, была предложена автомодельная формулировка соответствующих уравнений в подвижной системе координат, расположенной в кончике прямолинейно распространяющейся трещины. При этом все переменные были представлены в безразмерной форме.
К сожалению, применение интегральных уравнений для областей, заполненных неоднородным материалом, вызывает большие теоретические сложности [11], поэтому для приемлемой точности решения требуется привлечение конечно-элементной формулировки. Наиболее полная кончено-элементная постановка задачи о росте трещины типа KGD при наличии эффектов фильтрации жидкости ГРП в пласт выполнена в работе [12]. В качестве условия на кончике трещины рассматривается специальная формулировка для действующих в кончике сил когезионного сцепления материала, удовлетворяющая условию Г. И. Баренблатта о плавном смыкании берегов трещины. В теоретических расчетах использование такой формулировки предполагает либо выбор изначально мелкой сетки, либо привлечение адаптивной сетки, либо специальной автомодельной формулировки уравнений по аналогии с [7; 8].
Целью настоящей работы является привлечение механизмов адаптивной подвижной сетки для моделирования роста трещин ГРП. Рассмотрим упрощенную постановку задачи, которая раскрывает основные сложности моделирования с применением адаптивных сеток.
Система координат, граничные условия и исходная сетка показаны на (рис. 1). Квадратная область Q с размерами 10 х 10 м находится в условиях плоской деформации. На границе v = ±5 м перемещения в направлении v отсутствуют. На границе u = ±5 м задано давление p = 25 МПа. Вертикальный зародыш трещины длиной 1 м расположен в начале координат u = v = w = 0. Две равные и противоположные силы 12,5 МН приложены в начале координат к противоположным берегам зародыша будущей трещины. По всей области задано начальное напряжение с^ = 25 МПа. Остальные силы и напряжения в начальный момент времени
1
равны нулю. Материал среды обладает модельной трещиностойкостью KIC = 500H х м2. С момента t = 0 силы, приложенные в начале координат, постепенно возрастают, вызывая рост трещины. При достижении определенной величины нагрузки происходит прирост длины трещины вследствие превышения интенсивности напряжений в кончике K¡ трещино-стойкости материала KIC, при этом осуществляется адаптивное перестроение сетки с ее измельчением в области кончика трещины. Для анализа используется упрощенный алгоритм, в основе которого лежит VCCT (virtual crack closure technique), который хорошо описан в [13]. Предложенные в литературе методы адаптивного перестроения сеток в применении к задачам гидроразрыва пласта не позволяют моделировать трещины, способные искривляться в поле действующих напряжений и сил. Учитывая относительно сложную физику процесса ГРП, требуются достаточно эффективные формулировки системы уравнений и алгоритмы перестроения адаптивных подвижных сеток, которые позволили бы стабилизировать расчет и сократить расчетное время. Например, из анализа рис. 2 видно, что даже в условиях симметричного нагружения (закрепления) может происходить искривление траектории трещины, что, очевидно, вызвано ошибками округления при адаптивном перестроении сетки, а не физическими эффектами перераспределения напряжений.
Таким образом, можно констатировать, что даже простейшие расчеты динамики трещин с применением современных универсальных симуляторов напряженно-деформированного состояния необходимо проводить с большой осторожностью, поскольку алгоритмы адаптивного перестроения расчетных сеток могут давать ошибочные решения. Проблема выбора эффективного алгоритма генерации адаптивных подвижных сеток особенно актуальна для
Модуль Юнга: Е=40 ГПа
Коэффициент Пуассона: у=о,25
»
ОООе+(Ю7 Э.200*+006 •1.360е«007 •3.0406+007
-«.тгое+оо?
-6.400е+007 •8.0в0б+007 -9.760+00? -1.Ж6+009 -1.3126+008
-мэомоов 1,№«008 -1.016е+<И0 •1.98+64008 -2.1526+0«
-г.згос+оов
-2.4886+ОО0 •г.б5бб+оо0 _ -г.ег+моое
_ 2.992с4008
Щ *3.160б4С№
_ -з.з^е+сое
Н -3,4№с+008 Н -3.6646*0« Н -3 832е+008 Н .00064000
Шаг 1. Напряжения акл
Шаг 1. УССТ Energv Яе1еа5в
Шаг 25, Напряжения ахк
Шаг 25. Х/ССТ Епегёу Ие1еа$е
& +
Рис. 1. Постановка задачи. Эволюция напряжений и положение зародыша трещины
задач гидроразрыва пласта, поскольку в мире практически не существует лабораторных установок для проверки численных решений по всему спектру рассматриваемых физических эффектов ГРП.
Для эффективного моделирования роста трещин ГРП необходимо так подобрать физико-математическую формулировку основных законов сохранения, уравнений состояния, начальных и граничных условий и уравнений, отвечающих за генерацию расчетных сеток, чтобы обусловить перестроение подвижных сеток на каждом шаге вычислительного процесса. Этому условию удовлетворяет подход, предложенный в [14], который применялся для решения задач о сварке взрывом. Отличительной особенностью этого подхода является применение динамических уравнений равновесия многофазной упругой среды, которые лучше подходят для численного решения задач с применением параллельных вычислений. Следует отметить, что динамические уравнения равновесия пороупругопластической среды применялись в работах В. Н. Николаевского [15], однако уравнения равновесия были записаны в смешанной форме (Эйлер - Лагранж), что сразу определяло характер решения и ограничивало возможности метода.
Рис. 2. Эволюция адаптивной сетки (повернуто на 90°). Деформации показаны в увеличенном масштабе
Настоящая работа посвящена подбору модельного уравнения состояния для описания процесса сварки взрывом с учетом трехфазности модели и является обобщением формулировок уравнения состояния и условий их корректности, которые выработались в результате оригинальных вычислительных экспериментов и их сравнения с экспериментом на задаче о сварке взрывом металлических пластин [14; 16], а также на задачах о разрушении материала [17]. Оказалось, что созданная математическая модель естественным образом может быть расширена на задачи об образовании трещин в геологической среде, что также исследуется в настоящей работе.
В общем виде условие корректности уравнения состояния можно сформулировать в виде интервала допустимых значений для термодинамического потенциала - функции энтропии. Достижение нижней границы такого интервала свидетельствует о достижении «больших» отрицательных величин нормальных напряжений, что фактически является условием разру-
шения материала или отколом. На верхней границе такого интервала достигаются «большие» касательные напряжения, что свидетельствует о необходимости перестройки структуры материала, и для этого вводятся пластические деформации, которые уменьшают девиатор несимметричного тензора напряжений, таким образом, происходит «разгрузка». При дальнейшем уменьшении роли девиатора тензора напряжений осуществляется переход к гидродинамической модели и корректировка эффективной скорости звука, возникающей благодаря росту энтропии.
1. Математическая модель
уравнения состояния упругопластической среды
Напряженное состояние упругой среды зависит от деформации, которой эта среда подвергалась и от ее теплового состояния, описываемого термодинамической переменной (температурой или энтропией S). Для определения напряжений по заданным деформациям используется уравнение состояния среды с энергетическими затратами, потребными для осуществления этих деформаций. Деформация описывается при помощи отображения координат
х' = х' (е1,Е2,Е3,X),
составляющих среду материальных точек из «начального» состояния в деформированное. Начальные координаты Е1,Е2,Е3 носят название лагранжевых. Возникшее при деформации искажение координатной сетки задается так называемой матрицей дисторсии С :
С =
ж (е1, е2, е3)
ае 1
т. е. якобианом отображения х' (Е1,Е2,Е3,X). Для вычислений нам удобно рассматривать сингулярный вид матрицы дисторсии. Выполним в какой-либо точке нашей среды в рассматриваемый момент времени сингулярное разложение для матрицы С :
С =
(и1
V
м^ м-.
Ук о
о V V
0 к
3 /
2 2 v3 V?
3
= иКУт (иит = УУт = 1).
Мы будем использовать такой вид далее.
Уравнение состояния среды - это выражение зависимости внутренней энергии Е (на единицу объема и р0Е на единицу массы), от тензора деформации С и от термодинамического
параметра (энтропии) Е (С, S) = Е (С1, S). При адиабатическом (ёБ = 0) изменении полный
дифференциал матрицы дисторсии есть
ёЕ (С, Б ) = Ес,ёС1 = ЕС,ё
( 5х' ^
аЕ1
из которого следует формула Мурнагана для тензора напряжений Пиола - Кирхгофа:
^ = ЕС,.
Если среда изотропная, а именно такие среды мы и будем рассматривать, то энергетические затраты на деформацию можно выразить через инварианты симметрической матрицы
ССТ, которой присваивается название «тензор деформации». Это позволяет сформулировать общий вид уравнения состояния для моделирования широкого круга упругопластических сред. Для формулировки и исследования уравнения состояния нам оказалось удобнее использовать вместо матрицы ССТ ее квадратный корень \1ССТ. В этом случае уравнение состояния мы можем записать в виде
Е (С, 5 ) = Е (^СС7,5). Здесь матрица ^ССТ может быть представлена в виде
VССТ =^1 и куту к иТ =^1 и к киТ = ^ик (и к )Т = и киТ.
Тензор напряжений вычисляется по формуле
% = ЕС =Е (ССТ, 5)С = ЕссТ (ССТ )С = и ЕкК-1 иТ и К УТ =и ЕКУТ.
Следует отметить, что при изучении вращений системы координат тензор удобно
представить в виде шарового тензора О и девиатора (тензора с нулевым следом) V :
4ССТ = ^гл/СС^ + 4ССТ - 1ггу[ССТ13 = О + V, 3 3 3 3
где
О = ^гл/СС7^,
33
V = 4ССТ - 1ггу[ССТ13.
33
В наших математических моделях мы использовали более простой вид уравнения состояния вида и сделали три допущения.
1. Будем предполагать аддитивность по инвариантам уравнения состояния. Это допущение вполне обоснованно, так как каждый инвариант формулируется для описания определенного свойства упругопластического материала, а уравнение состояния в этом случае формулируется в виде суперпозиции рассматриваемых свойств.
2. Термодинамическую переменную 5 мы будем рассматривать отдельно и в качестве множителя к одному из инвариантов. Такая форма записи не умаляет сферу применения уравнения состояния, так как энтропийная переменная характеризует влияние от теплового состояния и связана, прежде всего, с изменением плотности упругопластической среды, которая описывается отдельным инвариантом.
3. Будем рассматривать изотропную среду, а именно такие среды и представляют интерес. В этом случае энергетические затраты на деформацию можно выразить через инварианты симметричной матрицы ССТ , которой присваивается название «тензор деформации».
В наших расчетах мы будем использовать следующие виды уравнения состояния для каждой фазы.
Уравнение состояния упругой среды:
Е = -(-^(5№У-1)+ - + J2.
у(У-1) У 3
Уравнение состояния жидкой среды:
У(У_1) У
Уравнение состояния газовой фазы:
Е а(5 )^(у-1)
У (у- 1)
Решаемые уравнения имеют следующий вид:
дыi дК _о дС1 _ди_0 дт дЕ * _ ' дт дЕ' _ '
_ Х0ЧЕ, Ес Ео1
дт ~ Е5 '
дс, _ _£,_ дс, _ с,
дт х0' дт х1 ' дт х,
Домножив уравнения на ы', Е(,], Ех, ЕСо, Ес,, Ес соответственно, получим закон сохранения полной энергии, который записывается в следующей форме:
дт
ыыг ^ V -- + Е
\
(и' к* )_ 0. дЕл ';
Метод решения уравнений теории упругости и его анализ подробно описан в книге [14].
2. Метод Годунова для решения уравнений теории упругости в лагранжевых координатах
Для разрешения модели упругопластической среды помимо решения уравнений теории упругости необходимо использовать процедуру перемещения расчетной сетки. Таким образом, вычислительный алгоритм можно представить в виде следующих шагов, которые выполняются для каждой ячейки расчетной области:
1) сингулярное разложение матрицы дисторсии С ;
2) при достижении необходимых условий учет фазового перехода и формулировка уравнения состояния;
3) контроль корректности уравнения состояния;
4) решение задачи Римана для сформулированного уравнения состояния;
5) расчет законов сохранения;
6) пересчет функции энтропии;
7) движение расчетной сетки.
Будем считать, что на начало временного шага нам известны следующие величины: скоп П 1 П 1 о п
рость ы1 , элементы тензора дисторсии ц. и напряжения л. , энтропия Л , полная энергия
1 ' П П Г"7 п п п и " тл
—ы' ыi + Е , скорость звука с01, и плотность р0 в каждой ячейке. Кроме этого, нам известны геометрические характеристики ячейки X. Схему Годунова для такой системы уравнений (законов сохранения момента импульса и полной энергии) можно записать в виде
и" - и" П] + -П]
= 0,
С I
ии.
^ + Е
V+1 ( , 1 ии
У 1тк
^ + Е
и + П1 + - и- П1 -
У1тк 1,1тк 1,1тк 1,1тк 1,1тк
= 0.
к
т
к
х
Величины потоков скорости и^1тк и напряжения Пбудут находиться из решения задачи
Римана. Затем происходит пересчет энтропии и движение расчетной сетки. После чего на основании нового положения ячейки пересчитывается матрица дисторсии С (подробности метода см. в работе [14]).
3. Моделирование упругопластических деформаций при движении фронта напряжений
Будем моделировать пластину размером 40 х 20 см с плотностью р = 3,0 г/см3, продольной и поперечной скоростями звука с1 = 5,2 км/с и с = 3,4 км/с соответственно. Критическое давление л:сг = 60 ГПа. По верхней границе пластины идет постоянный профиль давления, движущийся с постоянной скоростью V = 1 км/с. Результаты моделирования показаны на рис. 3, из которого видно, что в области фронта давления образуется довольно сложная геометрия сетки.
в
г
Рис. 3. Эволюция сетки на моменты времени t = 100 мкс (а), Г = 200 мкс (б), t = 300 мкс (в), / = 400 мкс (г)
Заключение
В работе представлена математическая модель упругопластических сред с учетом фазовых переходов между твердым телом, жидкостью и газом. Также учтен процесс откола вещества от газовой фазы. Сформулирован один из подходов к записи уравнения состояния упру-
гопластических сред, выписаны основные инварианты, используемые для построения уравнения состояния. Представлены результаты математического моделирования геологических деформаций, приводящих к образованию трещин.
Настоящая модель упругопластической среды была построена под руководством академика РАН Сергея Константиновича Годунова [16], которому благодарны все авторы данной статьи. В настоящем изложении исправлены неточности записи полной системы уравнений, которые имели место в книге [14]. В дальнейшем построенную упругопластическую модель планируется активно использовать для решения задач образования и эволюции трещин в геологической среде.
Список литературы
1. Желтов Ю. П., Христианович С. А. О гидравлическом разрыве нефтеносного пласта // Изв. АН СССР, Отдел технических наук. 1955. № 5. С. 3-41.
2. Баренблатт Г. И. О равновесных трещинах, образующихся при хрупком разрушении. Общие представления и гипотезы. Осесимметричные трещины // Прикладная математика и механика. 1959. Т. 23, № З. С. 434-444.
3. Ишлинский А. Ю. Прикладные задачи механики. М.: Наука, 1986.
4. Алексеенко О. П., Есипов Д. В., Куранаков Д. С., Лапин В. Н., Черный С. Г. Двумерная пошаговая модель распространения трещины гидроразрыва // Вестн. НГУ. Серия: Математика, механика, информатика. 2011. Т. 11, № 3. C. 36-59.
5. EconomidesM. J., Nolte K. G. Reservoir Stimulation. John Wiley & Sons, 2000.
6. Мусхелишвили Н. И. Сингулярные интегральные уравнения. М.: Наука, 1968.
7. Desroches J., Thiercelin M. Modeling the propagation and closure of micro-hydraulic fractures // International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts. 1993. Vol. 30, iss. 7. P. 1231-1234.
8. Garagash D., Detournay E. Similarity solution of a semi-infinite fluid-driven fracture in a linear elastic solid // Comptes Rendus de l'Académie des Sciences. Series IIB. Mechanics-Physics-Chemistry-Astronomy. 1998. Vol. 326, iss. 5. P. 285-292.
9. Desroches J., Lenoach B., Papanastasiou P., Thiercelin M. On the Modelling of Near Tip Processes in Hydraulic Fractures // International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts. 1993. Vol. 30, iss. 7. P. 1127-1134.
10. Desroches J., Detournay E., Lenoach B., Papanastasiou P., Pearson J. R. A., Thiercelin M., Cheng A. The Crack Tip Region in Hydraulic Fracturing // Proc. of the Royal Society A. 1994. Vol. 447, iss. 1929. P. 39-48
11. Katsikadelis J. T. The Boundary Element Method for Engineers and Scientists. Theory and Applications. Elsevier, 2016.
12. Байкин А. Н. Динамика трещины гидроразрыва пласта в неоднородной пороупругой среде: Автореф. дис. ... канд. физ.-мат. наук: 01.02.05. Новосибирск, 2016.
13. Krueger. R. Virtual Crack Closure Technique: History, Approach and Applications // Appl. Mech. Rev. 2004. Vol. 57: 2. P. 109-143.
14. Годунов С. К., Киселев С. П., Куликов И. М., Мали В. И. Моделирование ударно-волновых процессов в упругопластических материалах на различных (атомный, мезо и термодинамический) структурных уровнях / Ижевск. ин-т компьютерных исследований. Ижевск, 2014. 296 с.
15. Aidagulov G. R., Thiercelin M. J., Nikolaevsly V. N., Kapustyanskiy S. M., Zhilenkov A. G. Prediction of Long-Term Proppant Flowback on Weak Rocks // SPE Hydraulic Fracturing Technology Conference. 2007.
16. Годунов С. К., Киселев С. П., Куликов И. М., Мали В. И. Численное и экспериментальное моделирование образования волн при сварке взрывом // Тр. Математического института им. В. А. Стеклова. 2013. Т. 281. С. 16-31.
17. Favrie N., Gavrilyuk S., Ndanou S. A thermodynamically compatible splitting procedure in hyperelasticity // Journal of Computational Physics. 2014. Vol. 270. P. 300-324.
Материал поступил в редколлегию 01.03.2018
B. M. Glinskiy 1 2, S. I. Kabanikhin 1 2, S. A. Kalinin 3, I. M. Kulikov 1 2 R. Z. Kurmangaliev 4, N. S. Novikov 1 2, A. E. Ryazantsev 5, A. F. Sapetina 1 2 I. G. Chernykh \ M. A. Shishlenin 1 2
1 Institute of Computational Mathematics and Mathematical Geophysics SB RAS 6 Academician Lavrentiev Ave., Novosibirsk, 6630090, Russian Federation
2 Novosibirsk State University 1 Pirogov Str., Novosibirsk, 6630090, Russian Federation
3 JSC Geologika 23 Musa Dzhalil Str., Novosibirsk, 630055, Russian Federation
4 Khristianovich Institute of Theoretical and Applied Mechanics SB RAS 4/1 Institutskaya Str., Novosibirsk, 630090, Russian Federation
5 Research Complex of Technological Support Center Ltd. 20 Inzhenernaya Str., Novosibirsk, 630090, Russian Federation
gbm@sscc.ru, kabanikhin@sscc.ru, sergey.kalinin@list.ru, kulikov@ssd.sscc.ru rkurmangalievl@mail.ru, novikov-1989@yandex.ru, Ryazantsev1@gmail.com afsapetina@gmail.com, chernykh@parbz.sscc.ru, mshishlenin@ngs.ru
MATHEMATICAL MODELING OF ELASTOPLASTIC DEFORMATIONS
FOR PROBLEMS OF GEOLOGICAL CRACKS FORMATION AND EVALUATION
In this paper, the mathematical model of elastoplastic media with allowance for phase transitions between a solid, a liquid and a gas is presented. The process of splitting off the substance from the gas phase is also taken into account. An approach to the writing the state equation of elastoplastic media is formulated, basic invariants used for constructing the state equation are received. The results of mathematical modeling of geological deformations leading to the formation of cracks are presented.
Keywords: hydraulic fracture, elastoplastic medium, mathematical model, Godunov method, adaptive and moving mesh.
References
1. Zheltov Yu. P., Khristianovich S. A. O gidravlicheskom razryve neftenosnogo plasta [About hydraulic fracturing of oil reservoir]. Proc. of the AS USSR, Department of Technical Sciences, 1955, no. 5, p. 3-41. (In Russ.)
2. Barenblatt G. I. O ravnovesnyh treshchinah, obrazuyushchihsya pri hrupkom razrushenii. Obshchie predstavleniya i gipotezy. Osesimmetrichnye treshchiny [On the equilibrium cracks formed during brittle fracture. General representations and hypotheses. Axisymmetric cracks]. Journal of Applied Mathematics and Mechanics, 1959, vol. 23, no. 3, p. 434-444. (In Russ.)
3. Ishlinsky A. Yu. Prikladnye zadachi mekhaniki [Applied problems of mechanics]. Moscow, Nauka, 1986. (In Russ.)
4. Alekseenko O. P., Esipov D. V., Kuranakov D. S., Lapin V. N., Cherny S. G. Dvumernaya poshagovaya model' rasprostraneniya treshchiny gidrorazryva [Two-dimensional step-by-step model of propagation of a hydraulic fracture]. VestnikNSU. Series: Mathematics, Mechanics, Informatics, 2011, vol. 11, no. 3, p. 36-59. (In Russ.)
5. Economides M. J., Nolte K. G. Reservoir Stimulation. John Wiley & Sons, 2000.
6. Muskhelishvili N. I. Singulyarnye integral'nye uravneniya [Singular integral equations]. Moscow, Nauka, 1968. (In Russ.)
7. Desroches J., Thiercelin M. Modelling the propagation and closure of micro-hydraulic fractures. International Journal of Rock Mechanics and Mining Sciences & Geomechanics. Abstracts, 1993, vol. 30, iss. 7, p. 1231-1234.
8. Garagash D., Detournay E. Similarity solution of a semi-infinite fluid-driven fracture in a linear elastic solid. Comptes Rendus de l'Académie des Sciences. Series IIB. Mechanics-Physics-Chemistry-Astronomy, 1998, vol. 326, iss. 5, p. 285-292.
9. Desroches J., Lenoach B., Papanastasiou P., Thiercelin M. On the Modelling of Near Tip Processes in Hydraulic Fractures. International Journal of Rock Mechanics and Mining Sciences & Geomechanics. Abstracts, 1993, vol. 30, iss. 7, p. 1127-1134.
10. Desroches J., Detournay E., Lenoach B., Papanastasiou P., Pearson J. R. A., Thiercelin M., Cheng A. The Crack Tip Region in Hydraulic Fracturing. Proc. of the Royal Society A, 1994, vol. 447, iss. 1929, p. 39-48.
11. Katsikadelis J. T. The Boundary Element Method for Engineers and Scientists. Theory and Applications. Elsevier, 2016.
12. Baikin A. N. Dinamika treshchiny gidrorazryva plasta v neodnorodnoj porouprugoj srede [Dynamics of hydraulic fracturing in an inhomogeneous poroelastic medium]. PhD Thesis. Novosibirsk, 2016. (In Russ.)
13. Krueger. R., Virtual Crack Closure Technique: History, Approach and Applications. Appl. Mech. Rev, 2004, vol. 57: 2, p. 109-143.
14. Godunov S. K., Kiselev S. P., Kulikov I. M., Mali V. I. Modelirovanie udarno-volnovyh processov v uprugoplasticheskih materialah na razlichnyh (atomnyj, mezo i termodinamicheskij) strukturnyh urovnyah [Simulation of shock wave propagation in elastoplastic materials at different (atomic, meso and thermodynamic) structural levels]. Institute of Computer Science. Izhevsk, 2014, 296 p. (In Russ.)
15. Aidagulov G. R., Thiercelin M. J., Nikolaevsly V. N., Kapustyanskiy S. M., Zhilenkov A. G. Prediction of Long-Term Proppant Flowback on Weak Rocks. SPE Hydraulic Fracturing Technology Conference, 2007.
16. Godunov S. K., Kiselev S. P., Kulikov I. M., Mali V. I. Numerical and experimental simulation of wave formation during explosion welding. Proc. of the Steklov Institute of Mathematics, 2013, vol. 281, p. 12-26.
17. Favrie N., Gavrilyuk S., Ndanou S. A thermodynamically compatible splitting procedure in hyperelasticity. Journal of Computational Physics, 2014, vol. 270, p. 300-324.
For citation :
Glinskiy B. M., Kabanikhin S. I., Kalinin S. A., Kulikov I. M., Kurmangaliev R. Z., Novi-kov N. S., Ryazantsev A. E., Sapetina A. F., Chernykh I. G., Shishlenin M. A. Mathematical Modeling of Elastoplastic Deformations for Problems of Geological Cracks Formation and Evaluation.
VestnikNSU. Series: Information Technologies, 2018, vol. 16, no. 1, p. 61-73. (In Russ.)
