CC BY
111
Савенков Е.Б., Борисов В.Е. Математическая модель развития трещины гидроразрыва пласта в трехмерной пороупругой среде // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2018. - № 1. - С. 5-17. БОТ: 10.15593/регт.шесЬ/2018.1.01
Savenkov E.B., Borisov V.E. A mathematical model for hydraulic fracture propagation in three dimensional poroelastic medium. PNRPU Mechanics Bulletin, 2018, no. 1, pp. 5-17. DOI: 10.15593/perm.mech/2018.1.01
ВЕСТНИК ПНИПУ. МЕХАНИКА № 1,2018 PNRPU MECHANICS BULLETIN
http ://vestnik.pstu. ru/mechanics/about/inf/
БО1: 10.15593/регт.тесИ/2018.1.01 УДК 519.6
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАЗВИТИЯ ТРЕЩИНЫ ГИДРОРАЗРЫВА ПЛАСТА В ТРЕХМЕРНОЙ ПОРОУПРУГОЙ СРЕДЕ
Е.Б. Савенков, В.Е. Борисов
Институт прикладной математики им. М.В. Келдыша РАН, Москва, Россия
АННОТАЦИЯ
В настоящее время гидроразрыв пласта (ГРП) является одним из самых распространенных методов увеличения нефтеотдачи, применяемым при промышленной разработке нефтегазовых месторождений. В то же время используемые модели развития трещины ГРП являются во многом упрощенными: геометрия трещины предполагается плоской и заранее заданной, практически всегда отсутствует полноценный учет геомеханических эффектов, эволюция трещины зачастую описывается эмпирическими критериями. Несмотря на успешность применения таких моделей, их возможностей недостаточно для решения ряда важных задач разработки.
В работе представлена полностью трехмерная самосогласованная физико-математическая модель развития крупномасштабной трещины гидроразрыва пласта. Модель включает в себя несколько групп уравнений, включая неизотермическую пороупругую модель Био для описания поведения вмещающей трещину среды, двумерные уравнения смазочного слоя, описывающие течение в трещине, соответствующие условия согласования на границе «трещина-среда». Геометрическая модель трещины предполагает, что она является произвольной гладкой поверхностью с краем. Эволюция поверхности трещины в ходе ее роста описывается физически обоснованными критериями на основе векторной формы ^интеграла Райса-Черепанова. Модель пригодна как для анализа непосредственно процесса развития трещины ГРП, так и для анализа вызванных наличием трещины фильтрационных и геомеханических эффектов в ходе эксплуатации скважины с ГРП.
Основное назначение предложенной модели - согласованное описание процесса развития трещины ГРП в достаточно общей постановке с минимумом априорных допущений о характере протекания процесса и вместе с тем пригодное для использования на практике с применением современных вычислительных подходов. В связи с этим в работе дан краткий обзор вычислительных алгоритмов, пригодных для практической реализации модели.
©ПНИПУ
О СТАТЬЕ
Получена: 28 октября 2017 г. Принята: 26 февраля 2018 г. Опубликована: 30 марта 2018 г.
Ключевые слова:
геомеханика, пороупругая среда, трещина гидроразрыва пласта.
© Савенков Евгений Борисович - кандидат физико-математических наук, ведущий научный сотрудник, e-mail: e.savenkov@googlemail.com
Борисов Виталий Евгеньевич - кандидат физико-математических наук, научный сотрудник, e-mail: narelen@gmail.com
Evgeny B. Savenkov - CSc in Physical and Mathematical Sciences, Leading Researcher, e-mail: e.savenkov@googlemail.com Vitaly E. Borisov - CSc in Physical and Mathematical Sciences, Researcher, e-mail: narelen@gmail.com
Эта статья доступна в соответствии с условиями лицензии Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)
A MATHEMATICAL MODEL FOR HYDRAULIC FRACTURE PROPAGATION IN THREE DIMENSIONAL POROELASTIC MEDIUM
E.B. Savenkov, V.E. Borisov
Keldysh Institute of Applied Mathematics RAS, Moscow, Russian Federation
ARTICLE INFO ABSTRACT
Currently hydraulic fracturing (HF) is a stimulation technique which is most widely used during industrial development of gas and oil reservoirs. At the same time widely used mathematical models of hydraulic fracturing are often oversimplified, as fracture geometry is assumed to be planar and predefined, a comprehensive treatment of geomechanical effects is seldom considered, and the fracture growth is often assumed using empirical criteria. Despite their successful applications, their possibilities are not sufficient for solving a number of important problems of reservoir development.
This paper considers a complete three dimensional self-consistent mathematical model for large scale hydraulic fracture development. The model consist of several groups of equations including non-isothermal Biot poroelastic model to describe reservoir behavior, Reynold's lubrication equations to describe flow inside fracture and corresponding "reservoir"/"fracture" interface conditions. The fracture's geometric model assumes that it is an arbitrary smooth surface with a boundary. The fracture's surface evolution is governed by the physically-sound criteria based on J-integral of Rice and Cherepanov in the vector from. The model is suitable for describing hydraulic fracture development as well as for the analysis of flow and geomechanical effects induced by normal operations of a fractured well.
The main purpose of the suggested model is a consistent description of hydraulic fracture development in a general setting with a minimal number of a-priory assumptions and, at the same time, useful for solution of applied reservoir development problems using advanced numerical simulation techniques. Thus, we have also given a short review of the computational algorithms suitable for the model's implementation.
©PNRPU
Received: 28 October 2017 Accepted: 26 February 2018 Published: 30 March 2018
Keywords:
geomechanics, poroelastic medium, hydraulic fracture.
Введение
В настоящей работе предложена связанная математическая модель развития крупномасштабной трещины гидроразрыва пласта в пороупругой среде. Гидравлический разрыв пласта (ГРП) является одним из самых распространенных методов увеличения нефтеотдачи, используемым при промышленной разработке нефтегазовых месторождений. Сущность технологии ГРП заключается в закачке в нефтеносный пласт специальной жидкости с целью создания искусственной (техногенной) трещины длиной ~ 100 м, высотой ~10 м и средним раскрытием ~5—10 мм. В результате создается соединенный со скважиной искусственный канал с большой площадью притока, который имеет высокую (на порядки превышающую пластовую) проницаемость. Это обеспечивает значительное увеличение притока пластового флюида к скважине. Инженерные аспекты технологии рассматриваются, например, в [1].
Используемые в настоящее время математические модели развития трещины ГРП в большинстве своем имеют «инженерный» характер: в них используются двумерные постановки, априорно заданные геометрические модели трещин, практически никогда не производится корректный учет геомеханических эффектов и т.д. (см., например, [1-3]. Использование таких математических моделей давно стало стандартной инженерной практикой. Однако ряд задач, связанных с динами-
кой развития трещин ГРП, требует использования полных корректных (самосогласованных) пространственно-трехмерных постановок. Прежде всего к таким задачам относится анализ развития трещин в существенно неоднородном поле напряжений, индуцированном наличием соседних техногенных (например, множественный ГРП на горизонтальных скважинах) или естественных трещин, неоднородных полей давления, порожденных закачкой или отбором флюида в пласт и т. д.
Одновременно с полнотой постановки математические модели должны быть пригодны для использования на практике.
Целью настоящей работы является попытка построения такой модели. В связи с прикладной направленностью предлагаемой модели она дополнена разделом, который носит характер краткого обзора и посвящен описанию методов вычислительной математики, которые могут быть использованы при реализации модели в рамках программы-симулятора для математического моделирования ГРП в «расширенной» постановке.
Физико-математическое описание динамики трещины ГРП в ходе ее развития сводится к решению сложной связанной задачи, включающей в себя следующие системы уравнений:
• система уравнений пороупругости, описывающих эволюцию напряженно-деформированного состояния среды и полей давления флюида в ней в ходе развития трещины;
• уравнения течения жидкости в трещине;
• механические условия развития трещины, определяющие направление ее развития в каждой точке ее фронта;
• заданные на боковых поверхностях трещины условия согласования между полями давления в трещине и в среде, потоками массы, импульса и энергии, а также кинематические условия, связывающие раскрытие трещины и перемещение точек пласта.
Построенная в работе математическая модель позволяет корректно описывать согласованное распределение возмущений начального поля напряжений регионального масштаба с учетом эффектов, индуцированных процессами фильтрации при наличии добывающих и нагнетательных скважин, наличия трещины и ее развития с учетом как региональных, так и индуцированных полей напряжений, течения флюида в трещине. Модель пригодна для анализа эффектов самопроизвольного развития трещины ГРП («автоГРП») в пространственных масштабах элемента (ячейки) системы заводнения месторождения.
В качестве модели вмещающей среды используется связанная модель Био для описания возмущений начального поля напряжений регионального масштаба, вызванных присутствием и развитием трещины и эффектами фильтрации. Все основные соотношения модели приведены для неизотермического случая. Описание группы уравнений пороупругости модели основано на [4]. Отметим, что, вообще говоря, геологические среды на временных и пространственных масштабах, связанных с развитием трещины ГРП, могут проявлять пластические и/или вязкоупругие свойства. Учет таких эффектов в рамках рассматриваемой модели потребует изменения «механической» части модели Био и является скорее техническим усложнением, не меняющим принципиально структуру математической модели и способ связывания соответствующих групп уравнений. Вместе с тем при моделировании процессов развития трещины ГРП на практике зачастую используются именно линейно-упругие среды. Это связано с простотой указанных моделей и сравнительно высокой степенью надежности определения параметров среды (упругих модулей, модулей Био и критериев разрушения) в случае линейной среды.
Для описания течения флюида в трещине использованы уравнения гидродинамики в приближении смазочного слоя. Гидродинамическая модель течения дополнена уравнениями переноса энергии и позволяет производить учет неизотермических эффектов.
При описании процесса развития трещины использованы физически обоснованные критерии разрушения с учетом комбинированного режима нагружения в окрестности фронта распространяющейся трещины.
1. Модель вмещающей среды
В настоящем разделе описана согласованная математическая модель развития крупномасштабной трещины ГРП в напряженно-деформированном насыщенном
пласте. Основные допущения модели: напряженно-деформированное состояние пласта основано на физически и геометрически линейной модели; задача рассматривается в неизотермическом приближении, при этом отклонения температуры от начальной предполагаются малыми; пластовый флюид является однофазным и слабо сжимаемым; динамическими эффектами (силами инерции) при описании напряженно-деформированного состояния вмещающей среды можно пренебречь (то есть процесс является квазистационарным).
Уравнения механического равновесия (уравнения закона сохранения импульса) имеют вид
V-Т + р% = 0,
(1)
где р = р(х) - плотность насыщенной флюидом вмещающей среды; % - ускорение свободного падения; Т -тензор напряжений.
Описанные уравнения должны быть дополнены граничными условиями. При моделировании развития трещины ГРП граничные условия для давления можно задать как однородные граничные условия Неймана, обеспечивающие нулевой поток массы флюида, либо как граничное условие Дирихле, выражающее постоянство давления как функции времени (граница расчетной области при этом рассматривается как «контур питания»). Перемещения точек среды на внешней границе области в пластовых условиях неизвестны, по этой причине «механические» граничные условия обычно имеют вид граничных условий Неймана, которые выражают заданное значение нормальных напряжений на границе области.
В обоих случаях предполагается, что расстояние от трещины до границы области, занятой вмещающей средой, достаточно велико по сравнению с длиной трещины.
В общем случае анизотропного пороупругого насыщенного тела в линейном приближении определяющие соотношения имеют вид [4]
ДТ = С: Е - БД» - С: АД©, Дф = В : Е +—Др - 3а Д©,
N ф
где р - давление флюида; Т - тензор напряжений; Е = + ^®и)г]/2 - тензор деформаций; и -
поле перемещений точек вмещающей среды; С = [Сун ] -симметричный тензор упругих коэффициентов 4-го ранга; Б = [Бу ^ - симметричный тензор Био; А = |ау ^ -
симметричный тензор коэффициентов термического расширения; © - температура; ф - пористость; аф -
объемный коэффициент термического расширения поро-вого пространства; N - модуль Био, связывающий изменение пористости с изменением давления флюида; С - коэффициент теплоемкости скелета при постоянном объеме, Д/ = / - /0 для какой-либо величины /
индексом «0» обозначены значения величины в исходном состоянии.
В случае изотропной среды имеем
C = XI ® I + 201, B = Ы, A = а!,
где X, О - коэффициенты Ламе скелета; е = Е : I = Еи -объемная деформация; X = К - (2/3)0; К - объемный модуль упругости.
Плотность р в (1) зависит от плотности скелета и флюида и пористости вмещающей среды: р = (1-ф)ру +фр/, где ру - плотность скелета; р / - плотность флюида.
При деформации породы происходит изменение объема, занятого флюидом. Это изменение связано как с изменением объема скелета, так и с изменением объема порового пространства за счет деформации скелета. Уравнения закона сохранения массы флюида имеют вид
дт, №
= д, V = = К-(-У/7 + р/8), (2)
д р/ где - масса флюида в единице объема насыщенной среды; № т = р/ V - вектор плотности потока массы; V -
скорость фильтрации, определяемая законом Дарси; р / - плотность флюида; К - симметричный тензор коэффициентов проницаемости; д - массовая плотность источников.
Считая, что в элементарном объеме насыщенной среды скелет и флюид находятся в состоянии локального термодинамического равновесия, уравнение закона сохранения энергии (при отсутствии внешних источников энергии) запишется в виде
де ~8t
+ Vh w т + q®) = 0, q® = -AV®,
где к/ - энтальпия флюида; q© - вектор плотности теплового потока за счет эффектов теплопроводности; е -полная внутренняя энергия элемента объема насыщенной среды; е = ру (1 - ф)еу + р/фе/ ; Л - симметричный тензор
коэффициентов теплопроводности; еу / и ру / - внутренние энергии и плотности твердой и подвижной фаз.
Приведенные выше уравнения должны быть дополнены определяющими соотношениями. В дальнейшем будем считать, что отклонения от начального состояния малы и флюид слабо сжимаем, то есть
Ар, Др А© „ Ар ^ А© _
--3а,-, Ду, = -3а + Ср-, (3)
Kf
®„
f РУ~' ®0
Р/ 'V
где p - давление флюида; рf - его плотность; af -коэффициент температурного расширения флюида; ® -его температура; sf - энтропия; Cp - теплоемкость при постоянном давлении; Kf - объемный модуль сжатия флюида. Далее будем считать, что Kf = const, C = const, a . = const - заданные параметры, завися-
щие от точки пространства, но не зависящие от времени и состояния среды. В частном случае несжимаемой жидкости (К/ , а=0) и малых изменений температуры имеем р/ = р/, Ду/ = СрА© / ©0. Зависимости внутренней энергии от давления и температуры (уравнения состояния) могут быть получены непосредственно из уравнений (3).
С учетом условий малых изменений свойств насыщенной среды, можно показать, что
Amf = P/vf,
vf = B : AE + — Ap - 3a A®
f M m
(4)
где vf - относительное изменение объема жидкости в элементарном объеме насыщенной среды.
2. Геометрическая модель трещины
Толщина (раскрытие) трещины ГРП имеет среднюю величину порядка единиц миллиметров, максимальную -порядка единиц сантиметров. В то же время характерный размер области □, в которой решается задача (масштаб элемента системы заводнения), составляет десятки-сотни метров, характерная полудлина трещины - десятки-сотни метров, высота - до единиц-десятков метров. По этой причине не имеет смысла рассматривать трещину как трехмерный объект при изучении пороупругих процессов во вмещающей трещину среде. Поэтому при построении модели насыщенного напряженно-деформированного пласта будем считать трещину разрезом нулевой толщины, который соответствует срединной поверхности Т трещины. Срединной поверхностью назовем условную поверхность, равноудаленную от берегов трещины.
Будем считать, что срединная поверхность трещины является гладкой, односвязной поверхностью
Т: DT ^ ]
(5)
с краем дТ , который соответствует фронту трещины.
Обозначим через Т, рассматриваемую поверхность в момент времени ,, при этом для > ,2 выполняется условие Т с Т2 .
В дальнейшем будем считать, что в каждый момент времени на линии у(,) = дТ, задано векторное поле
скорости V = v(x,,), х еу(,), описывающее ее эволюцию. Таким образом, движение точки границы описывается уравнением
dx
— = V(xtX x |t=0= x0 ey(0X dt
(6)
где поле скорости V = v(x,,) является гладкой функцией точки х е у(,) в каждый фиксированный момент времени , и гладкой функцией времени для каждой фиксированной (лагранжевой) точки границы.
3. Модель течения в трещине
Будем считать, что флюид распространяется по трещине с заданной геометрией, стремясь заполнить всю трещину. Таким образом, части трещины, занятой жидкостью, соответствует лишь часть срединной поверхности Т трещины. Раскрытие трещины является функцией точки поверхности Т , ^ = ^(х), х е Т .
Относительно срединной поверхности будем предполагать, что ее кривизна мала.
Поверхность Т является срединной поверхностью трещины в том смысле, что боковые поверхности Т± реальной, «физической» трещины находятся на расстоянии /2 от срединной поверхности Т , Т± = х + п(х)^(х)/2, х е Т . Здесь п(х) - нормаль к поверхности Т в точке х, направленная в сторону Т-
(таким образом, п+ = п, п- = -п). Естественными ограничениями на вид функции ^ является ее гладкость, неотрицательность и равенство нулю на крае дТ, ^(х) > 0 , х е Т , ^(х) = 0 , х е дТ .
В приближении смазочного слоя уравнение течения сжимаемой жидкости, связывающее раскрытие w и давление р/ в трещине, имеет вид [5]
др ^ ~дГ
+ Vт(vт) = , ут = -Б(/-Р/%Т) (7)
где ут - средняя по раскрытию трещины скорость течения флюида в трещине, коэффициент Б является эффективной проницаемостью трещины и в общем случае определяется свойствами жидкости и состоянием среды; % -
проекция вектора ускорения свободного падения на срединную поверхность трещины, %т =(1 - п ® п) . В простейшем случае ньютоновской жидкости Б = w2 / (12ц/), где ц/ - динамическая вязкость флюида.
В общем случае неньютоновской жидкости Б = Б(w, VТр/) (см., например, [6]). Символом VТ обозначен поверхностный оператор Гамильтона [7].
Величина Qf в правой части уравнения (7) описывает мощность источников. Для трещины ГРП она может быть представлена в виде Qf = + Q0" , где первое слагаемое описывает мощность потока флюида из скважины в трещину (и задается как точечный источник, мощность которого определяется параметрами закачки и моделью течения в скважине), второе - утечку флюида в пороупругий пласт (определяется условиями согласования, см. п. 4).
Уравнение (7) рассматривается в ограниченной области О/ = О/ (/) с Т , занятой флюидом. Границу области О/ (/) будем обозначать как дО/ (/). Будем считать, что О/ - односвязная область, а дО/ - замкнутая
гладкая кривая (рис. 2). На границе дО/ этой области
должны быть заданы соответствующие граничные условия, которые будут рассмотрены в п. 5.
В ряде практически важных случаев уравнение (7) может быть дополнено уравнением, описывающим перенос флюидом в трещине того или иного пассивного компонента (например, «проппанта» - специального калиброванного «песка», предназначенного для того, чтобы трещина не «закрылась» в случае, если давление флюида в ней упадет достаточно сильно). Соответствующие модели широко известны (см., например, [6, 8]) и практически не меняют структуру рассматриваемой модели и вид ее уравнений. По этой причине здесь они не рассматриваются.
В общем случае эволюция фронта жидкости описывается так же, как и фронта трещины (см. п.п. 2 и 5).
В случае необходимости учета неизотермических эффектов при описании течения в трещине уравнения сохранения массы флюида дополняются уравнением сохранения энергии. В случае слабосжимаемой жидкости в трещине (трехмерное) уравнение закона сохранения энергии
имеет вид (без учета внешних источников тепла)
д
— (р У е) + V ((©V© + р С © У Т ) - Б : Т = 0, (8)
где р/Де = р/СД© - объемная плотность внутренней
энергии флюида; X© - коэффициент теплопроводности
флюида; р£ - объемная теплоемкость флюида; Б -
тензор скоростей деформации. Последнее слагаемое в (8) описывает диссипацию энергии за счет работы вязких сил во флюиде.
Интегрирование уравнения (8) по толщине трещины в интервале от ^ /2 до +w /2 дает
д
- (р/С©Т )
д/
х (-wX Т V Т ©Т + wp уС©Т ут ) - V + д© =0, (9)
где ©Т - средняя по толщине трещины температура; V - диссипация за счет работы вязких сил; д© - приток тепла через боковые поверхности трещины.
С учетом допущений о структуре течения в тонком смазочном слое (локально пуазейлевский профиль, малые числа Рейнольдса, малость толщины слоя) для величины диссипации V может быть получено выражение V = (12цТ) / w уТ . Для течений рассматриваемого
класса эта величина обычно мала и ей можно пренебречь, то есть считать V = 0 .
4.Условия согласования «вмещающая среда»/«трещина»
Величины, описывающие напряженно-деформированное и насыщенное состояние пласта (поле перемещений, тензоры напряжений и деформаций, давление
флюида в пласте, скорость фильтрации) терпят при переходе через ее срединную поверхность Т скачок, величина которого определяется специальными условиями согласования. Эти условия связывают рассмотренные выше независимые группы уравнений, описывающие состояние пороупругой вмещающей среды и течение в трещине. Условия согласования включают в себя кинематические условия и условия, обеспечивающие непрерывность потоков импульса, массы и энергии. Рассмотрим их последовательно.
Пусть далее х е Т - точка на срединной поверхности
трещины, х - соответствующие ей (в смысле проекции кратчайшего расстояния) точки на боковых поверхностях («берегах») Т± трещины, п± - вектор единичной внешней нормали к вмещающей среде на боковых поверхностях трещины (рис. 1). Как и ранее (см. п. 3), считаем
п+ = п, п- = -п, (10)
где п - единичная нормаль к срединной поверхности Т , направленная в сторону от Т + к Т- (см. рис. 1).
Рис. 1. Ориентация боковых поверхностей (берегов) трещины Fig. l. Orientation of lateral surfaces (shores) of the crack
К первой группе условий согласования относятся кинематические условия, связывающие перемещения берегов трещины и ее раскрытие. Эти условия имеют вид
w(x) = (u + - u- ) • n = [u] • n = [un ],
где u± = u± • n; [un ] = u,+ - u- ; u± = u(x± ). В рассматриваемой постановке задачи нельзя сделать естественных априорных предположений о скачке тангенциальных компонент поля перемещений при переходе через поверхность Т : величина этого скачка неизвестна и определяется решением задачи.
Рассмотрим теперь динамические условия согласования, которые имеют вид условий непрерывности нормальной компоненты тензора напряжений на берегах трещины:
T+ • n+ = -p*n +, T- • n- = -p*n-, откуда с учетом (1G) имеем
T+- n = -p*n, T- • n = -p*n, где T± = T(x± ).
Эти соотношения можно записать в эквивалентном виде как соотношения для скачка и среднего значения нормальных напряжений при переходе через срединную поверхность Т :
[Т] • п = 0, {Т} • п = -р* • п, (11)
где (здесь и далее) [а] = а + - а-; {а} = 1/2 (а + + а-) для какой-либо величины а. В последнем уравнении величина р* имеет смысл давления, приложенного к берегам трещины. Конкретный вид р* и его связь со значением давления во вмещающей среде и трещине будут рассмотрены ниже. Условия (11) предполагают, что силы вязкого трения, действующие на берега трещины со стороны заполняющего трещину флюида, пренебрежимо малы.
Рассмотрим теперь условия согласования для «гидродинамической» части задачи, включающей в себя уравнения фильтрации в пористой вмещающей трещину среде и уравнение течения непосредственно в трещине. Эти условия связывают давление в трещине, вмещающей среде и соответствующие потоки массы. Они будут иметь различный вид в зависимости от допущений, накладываемых на условия течения.
В первом случае будем считать, что трещина и вмещающая среда заполнены одним и тем же флюидом, который заполняет весь объем трещины, причем гидродинамический контакт между трещиной и вмещающей средой - идеальный. В этом случае в точке х е Т
Рис. 2. Взаимное расположение срединной поверхности трещины и области, занятой флюидом Fig. 2. Relative position of fracture mid-surface and fluid domain
срединной поверхности имеем
p + = p- = Pf, (12)
где, как и ранее, p± = p(x±), или, что эквивалентно
[p] = 0 , {p}= pf .
Условия непрерывности потока массы с учетом соответствующих выражений для течения во вмещающей
lG
среде (2) и вида уравнения (7), выражающего закон сохранения массы в трещине, имеют вид
Qf = Q + + Q-, Q± = 'т • п±| ±. (13)
Таким образом, величина потока из среды в трещину определяет правую часть закона сохранения массы флюида в трещине, а давление в трещине играет роль граничного условия для давления на берегах трещины.
В рассматриваемом случае естественно считать, что область срединной поверхности, соответствующая объему трещины, занятому флюидом, совпадает со всей срединной поверхностью, то есть О/ = Т . При этом
р* в (11) определено как р*(х) = р/(х).
Во втором случае будем считать, что флюид в трещине не проникает во вмещающую среду. Это может быть связано как механической изоляцией объема вмещающей среды, так и существенным отличием свойств флюида в трещине от флюида во вмещающей среде. В этом случае будем считать, что часть срединной поверхности, соответствующая объему, занятому флюидом, не совпадает со всей срединной поверхностью трещины, то есть О/ с Т, дО/ п дТ = 0 (см. рис. 2).
Тогда условия согласования для потока массы флюида в области О / имеют вид
Qf =0, 0± =0, х еО/. (14)
Другими словами, условия согласования вырождаются в граничные условия Неймана для уравнения фильтрации во вмещающей среде. В этом случае условия для давления во вмещающей среде не могут быть сформулированы, а давление в трещине определяется уравнением (7) с заданной (нулевой) правой частью Qf .
Величина р* в (11) задается как р (х) = pf (х).
Аналогичные условия согласования в области Т \ О / (области «лага») могут быть заданы по-разному
в зависимости от того, являются ли берега трещины проницаемыми или нет для флюида во вмещающей среде. В первом случае в точках х е Т \ О/ имеем Q± =0
и р*(х) = 0 в (11). Во втором случае естественно считать, что в рассматриваемой области трещины, соответствующей части срединной поверхности Т \ О/, находится поровый флюид, давление которого непрерывно при переходе через берега трещины, то есть [ р] = 0, [О] = 0.
Наконец, рассмотрим практически важный случай, когда флюид в трещине не проникает во вмещающую среду, но содержит в своем составе компонент, который может проникать во вмещающую среду. Тогда обычно считается, что потоки на берегах трещины определены как функции свойств флюида, фильтрационно-емкостных свойств среды и полей давления в трещине и ее окрестности. В рамках рассматриваемой модели это приводит к условиям
Qf = Q + + Q-, Q± = Q±(р/,р±) = '±т-п±.
Последняя зависимость может быть задана, например, в соответствии с моделью утечки Картера или ее обобщений [9, 10].
Рассмотрим теперь баланс потока энергии. Будем считать, что в точках границы Т± насыщенная пороуп-ругая среда (скелет и флюид) и флюид в трещине находятся в состоянии локального теплового равновесия. В этом случае имеем ©± = ©Т , или, эквивалентно, [©] = 0 , {©} = ©Т . В случае если утечка массы из трещины во вмещающую среду не происходит, условия непрерывности потока энергии будут иметь вид
д© = 0+ + <2©, <2© = я©-п±.
При этом первое уравнение, по существу, является определением для мощности источников энергии в уравнении (9), а температура в трещине играет роль граничного условия для уравнения закона сохранения энергии в среде.
В том случае если утечка массы происходит и, так или иначе, задана значениями Q±, условия непрерывности потока энергии будут иметь вид
д© = 0+ + О©, <2© = Ь/О* = (Ь/'т + я© )|х=х±- п±.
Тогда заданным можно считать поток из трещины во вмещающую среду. Значение температуры при этом определятся полным решением задачи.
5. Граничные условия на фронте трещины
Под действием внешних нагрузок (объемных сил и «механических» граничных условий) и в процессе закачки флюида в трещину она может развиваться. С математической точки зрения это соответствует развитию срединной поверхности трещины. В дальнейшем будем считать, что развитие трещины происходит квазиравновесным образом. Это означает, что в любой момент времени среда находится в состоянии механического равновесия, которое определяет стационарное положение фронта трещины. Условия, определяющие рост трещины, по существу, являются уравнениями, определяющими положение фронта трещины в каждый конкретный момент времени.
Изучение условий развития трещин является предметом механики разрушения. Более полный обзор этих вопросов в настоящей работе не предусматривается, для знакомства с ними см., например, [11, 12].
В рамках настоящей работы будем считать, что разрушение среды является хрупким и рассматривается в рамках линейно-упругой модели среды. Для таких задач вопрос о возможности развития трещины обычно решается путем анализа значений коэффициентов интенсивности напряжений в точках фронта трещины. Коэффициенты интенсивности напряжений, в свою
очередь, определяются как коэффициенты при главных сингулярных асимптотиках поля напряжений в окрестности точки фронта трещины. В общем случае можно определить три таких коэффициента, соответствующих трем режимам нагружения трещины (которые определяются относительной ориентацией касательным к фронту трещины вектором и направлением нагрузки, приложенной к берегам трещины).
Альтернативным подходом является использование энергетического критерия разрушения. В этом случае с элементом площади срединной поверхности трещины связывается некоторая энергия, необходимая для его образования. Соответственно, развитие поверхности трещины связано с потоком энергии, который относится к соответствующей точке фронта трещины. Формализация этого рассуждения приводит к понятию J -интеграла Черепанова-Райса [11]. В настоящей модели используется именно такой подход. При этом используется так называемая векторная форма J -интеграла [11]. Вариант векторного J -интеграла для задачи теории пороупругости рассмотрен в [13]. Конкретный вид и способ вычисления векторного J -интеграла в настоящей работе не рассматривается. Укажем лишь его «функциональную» форму, достаточную для замыкания уравнений рассматриваемой модели.
А именно, в соответствии с [13] будем считать, что при заданном положении фронта 8Tt трещины в момент времени t в каждой его точке x е 8Tt можно определить вектор J = J(x). С физической точки зрения вектор J определяет количество энергии J , высвобождаемой при локальном развитии трещины в заданном направлении l (|| 11|=1), J = J• l. В соответствии с энергетическим критерием развития трещины ее рост происходит в направлении вектора 1, при котором величина высвобождаемой энергии максимальна,
J = Jmax = J(1 max). В случае если Jmax < J*, ^ J* - некоторое критическое значение энергии, которое обычно считают свойством материала, роста трещины не происходит. При Jmax = J, происходит устойчивый рост трещины, при Jmax > J, - неустойчивый.
Таким образом, если рост трещины в заданной точке происходит, локальное направление роста совпадает с направлением вектора J , то есть в уравнении (6) имеем v = a J , a > 0 .
Вместе с тем теория квазиравновесной трещины в упругой среде не позволяет непосредственно определить скорость развития трещины, то есть значение коэффициента a в случае, если он положителен. Поэтому уравнение (6) является «кинематическим», то есть дает удобную параметризацию поверхности трещины. Оно должно быть дополнено непосредственно уравнением для определения положения фронта трещины, которое уже было сформулировано выше и имеет вид Jmax < J,. Это уравнение является нелокальным и зависит от ис-
тории развития трещины, текущей формы ее поверхности, решения задачи пороупругости во всей области и распределения давления внутри трещины.
В совокупности приведенные уравнения означают, что в ходе развития трещины в каждой точке ее фронта вектор, касательный к поверхности трещины и нормальный к ее фронту, соосен вектору Л в заданной точке, при этом в этой точке выполняется соотношение Jmax < J*. При этом а = 0 в точках, где Jmax < J*. В точках, где происходит рост трещины, справедливо уравнение Jmax = J*. Это условие, по существу, является уравнением относительно параметра а, который определяет скорость, с которой точка фронта трещины «продвигается» в направлении вектора Л.
Приведенные выше рассуждения касаются исключительно вопроса о развитии трещины с точки зрения модели поведения вмещающей среды. Рассмотрим теперь случай трещины, заполненной флюидом.
Будем считать сначала, что область, заполненная флюидом, «меньше», чем область трещины. Другими словами, в любой момент времени , имеем О/ (,) с Т
и дО/ (,) п дТ = 0 . Это условие предполагает (если не
указано обратное), что в каждой точке фронта раскрытие положительно, то есть ^(х) > > 0, х е Г(,).
В этом случае флюид в трещине «растекается» по уже существующей трещине, пытаясь заполнить весь ее объем. Часть О/ срединной поверхности, соответствующая области трещины, занятой флюидом, в общем случае изменяется за счет притока флюида в трещину и изменения раскрытия трещины за счет деформации вмещающей среды.
Пусть п - единичная внешняя нормаль к дО/ , касательная к срединной поверхности трещины, V - скорость движения границы Г(,) = дО/ (,), ^ - как и ранее,
раскрытие трещины, vТ - средняя по раскрытию скорость движения частиц жидкости, р / - плотность жидкости, Qf - мощность массового источника жидкости, Ч/ = Qf / р/ - мощность объемного источника.
Интегрируя уравнение (7) по области, занятой флюидом, имеем
, др м , ] +.|ру =
О г (,) Г(,)
dt
= j Qfds,
а. (t)
vn = v^ • nlr(t)■
(15)
Это уравнение будет справедливо, если поток массы через границу области флюида связан исключительно с движением границы, то есть отсутствуют потоки массы, связанные с испарением через фронт флюида и так далее, и источник в правой части не является сингулярным на фронте. В таком случае нормальная компонента
скорости распространения фронта жидкости ¥п совпадает с нормальной компонентной уп средней скорости движения частиц жидкости на фронте.
dx f
lim vn = Vn = -f- • n, xf еГ(/).
x^xz dt
(16)
Несмотря на свою очевидность, соотношение (16) носит фундаментальный характер и связывает скорость движения точек границы со скоростью жидкости. Его называют уравнением для скорости фронта («speed equation») [14] или условием Стефана («Stefan condition») [15]. Как правило, на фронте тангенциальная компонента V частиц жидкости мала по сравнению с vn. В таком случае V = v , то есть скорости фронта и частиц жидкости совпадают (в общем случае это касается лишь их нормальных компонент).
В терминах потока уравнение для скорости (16) можно записать в виде
.ln _ q^ • n Vn = lim —1- = lim-,
йл.
w
х/ еГ(/); Ят = р.^Ут. (17)
х/ PfW
Отметим, что уравнения (16) и (17) верны и в том случае, если на фронте жидкости раскрытие равно нулю (именно поэтому в этих уравнениях используется обозначение х ^ х/, а не х = х/).
Теперь перейдем к обсуждению возможных граничных условий, которые могут быть наложены на уравнение (7) с учетом (16) и (17).
Первый вариант предполагает задание на фронте жидкости граничного условия (Дирихле) для давления,
р/ (х) = р/ (x), х ^(Д
где pf - заданная функция. Скорость движения фронта
жидкости в этом случае определяется уравнением (17). Эволюция области, занятой флюидом, определяется уравнением, аналогичным (6), где скорость определяется уравнением (17). Значение р* в точке фронта можно
положить равным нулю в том случае, если пластовый флюид не проникает в область трещины, не занятую жидкостью, и р(х/), если область «лага» заполнена
пластовым флюидом.
Во втором варианте предполагается задание в каждой точке границы фронта флюида значения скорости движения точки фронта, или, что эквивалентно в силу (17), потока массы на фронте. Эволюция фронта трещины определяется указанным выше способом. При этом значение скорости движения фронта не может быть задано произвольно - она должна быть согласована с потоком массы в трещине и текущим значением раскрытия трещины соотношением (15).
Отметим, что в рассматриваемом случае естественным является первый способ задания граничных условий. При этом в задаче присутствуют два фронта - трещины и флюида, эволюция которых определяется независимо.
Наконец, рассмотрим случай, когда флюид может заполнять весь объем трещины. С математической точки зрения этот случай является наиболее сложным. Это связано с тем, что в точках фронта трещины раскрытие равно нулю, и поток в уравнении (7) формально вырождается. Однако в силу того, что фронт трещины и совпадающий с ним фронт флюида распространяются с конечной скоростью, скорость флюида на фронте, определенная в соответствии с (17), конечна. В силу этого значение давления в точках фронта имеет особенность и является неограниченным. Корректная постановка граничных условий для уравнения течения в трещине в этом случае является нетривиальной задачей. «Естественным» граничным условием в этом случае является задание скорости движения фронта, которая, в свою очередь, определяется решением задачи о развитии трещины. При этом содержательным становится использование соотношения (17), которое в данном случае рассматривается не как определение скорости движения фронта, а как уравнение, связывающее скорость движения фронта флюида и трещины. Различные способы корректной постановки граничных условий в этом случае рассмотрены, например, в [14, 15].
Сингулярность поля давления в трещине может быть устранена допущением о том, что в ходе развития трещины фронт флюида всегда немного отстает от фронта трещины на заданную величину. В этом случае говорят, что рассматривается модель трещины с «лагом». Величина лага заранее неизвестна. В общем случае она может быть определена рассмотрением постановки задачи с одновременным учетом эволюции двух фронтов - флюида и трещины.
Вместе с тем известно, что а) для практических расчетов важно существование лага как такового, а не его конкретная величина и б) учет лага существенен только на первоначальном этапе развития трещины ГРП (в дальнейшем величина лага становится несущественной и ей можно пренебречь). Поэтому зачастую величина лага определяется из тех или иных дополнительных соображений (см., например [16-18]) или считается малой фиксированной величиной.
6. Алгоритмические аспекты реализации модели
Описанная в предыдущих разделах математическая модель мало пригодна для теоретического анализа процесса развития трещины в силу своей сложности. Ее основное назначение - использование в виде математической модели для анализа средствами математического моделирования и вычислительного эксперимента. По этой причине ниже кратко представлены базовые алгоритмы вычислительной математики, которые могут быть использованы для реализации модели в описанном выше виде. Наличие таких методов, вообще говоря, обосновывает сложность модели и демонстрирует возможность ее реализации на практике.
Для решения уравнений, описывающих состояние вмещающей среды, то есть уравнений пороупругости, в настоящее время предложено значительное количество подходов. Последние включают в себя метод конечных разностей, метод конечных элементов и метод граничных интегральных уравнений. В силу того что модель предполагает произвольное пространственное распределение свойств вмещающей среды, наиболее приемлемым из указанных является метод конечных элементов [19]. Для аппроксимаций по времени чаще всего используются полностью неявные схемы (например, неявный метод Эйлера). Для решения конечномерной задачи могут использоваться как «монолитные» постановки (в рамках которых одновременно определяется как поле перемещений, так и поле давления), так и методы на основе итераций между группами уравнений упругости и фильтрации [20].
Характерным свойством системы уравнений поро-упругости в рассматриваемом приближении является то, что после аппроксимации по времени она имеет вид задачи о седловой точке [21]. В этом случае для устойчивости решения как континуальной, так и конечномерной задач необходимо выполнение так называемых inf-sup условий (условий Ладыженской-Бабушки-Бреззи) [21]. При нарушении этих условий в задачах пороупругости могут возникать численные неустойчивости конечномерного решения, особенно в предельном случае несжимаемой жидкости). Простейшим примером конечного элемента, удовлетворяющего этим условиям, является элемент Тэйлора-Худа, в котором для аппроксимации перемещений используются конечные элементы второго порядка, а для аппроксимации давления -первого [22-24]. Альтернативным подходом является регуляризация исходной задачи (см., например, [2527]). В случае если для конкретной задачи inf-sup устойчивые пары пространств неизвестны, необходим контроль за особенностями решения и проверка выполнения численных inf-sup условий [28].
Наличие в постановке задачи трещины приводит к тому, что решение связанной задачи термопороупруго-сти терпит разрыв при переходе через ее срединную поверхность. При этом компоненты полей перемещений, тензора деформаций и напряжений терпят сильный разрыв; поля давления и температуры являются непрерывными, соответствующие потоки (массы флюида и энергии) терпят, вообще говоря, сильный разрыв. В рамках традиционных подходов считается, что разрыв проходит по ребрам и/или граням сетки. Использование таких подходов требует применения дорогостоящих процедур перестройки расчетной сетки на каждом шаге развития трещины и корректной (консервативной) интерполяции расчетных величин со «старой» сетки на «новую».
Более оправданным является использование методов, предполагающих подсеточное разрешение трещины. В методах этого класса считается, что геометрия срединной поверхности трещины не согласована с геометрией сетки во вмещающем трещину пространстве. Для учета
разрыва решения в ячейках сетки, через которые проходит трещина, используются специальные базисные функции, допускающие наличие сильного или слабого разрыва конечномерного решения внутри этих ячеек. К методам этого класса относятся методы типа Assumed Strain Elements [29, 30], X-FEM (extended Finite Elements), метод Нитше (Nitsche method) и др. При использовании методов этого класса перестройка сетки не является необходимой для корректного описания развития трещины процедурой: если она и проводится, то с целью повышения фактической точности решения. Отметим, что в контексте метода X-FEM обычно используется неявное представление срединной поверхности трещины, когда она описывается как поверхность уровня ноль заданной в объеме функции. Для описания эволюции таких поверхностей разработан целый ряд эффективных методов на основе решений уравнений типа Га-мильтона-Якоби (см., например, [31, 32]).
В случае если срединная поверхность трещины не согласована с геометрией расчетной сетки во вмещающей среде, возникает вопрос о конечномерной аппроксимации условий согласования на границе «трещина-пласт». В частности, условия равенства двух функций, определенных в объеме и на срединной поверхности, не могут быть выполнены поточечно. В этом случае стандартным способом является аппроксимация условий равенства в слабом (в смысле вариационной постановки задачи) смысле. Типичными примерами таких подходов являются методы штрафа, классический и регуляризо-ванный метод множителей Лагранжа и метод Нитше.
Построение inf-sup устойчивых пространств для множителей Лагранжа для конечно-элементных аппроксимаций метода X-FEM рассмотрено в работах [33-36]. Метод Нитше изначально был предложен в работе [37] и в дальнейшем развит, в том числе в контексте метода X-FEM в работе [38].
Наконец, рассмотрим вопросы, связанные с решением уравнения смазочного слоя, определенного на срединной поверхности трещины. В том случае, если срединная поверхность трещины образована гранями расчетной сетки во вмещающей среде, решение этой задачи не представляет принципиальных проблем (см., например, [39] и ссылки там). В том случае, если поверхность трещины не согласованна с сеткой, а сама поверхность задана неявно, задача становится значительно сложнее. Можно выделить целый ряд методов, разработанных для ее решения.
В первом классе методов (условно - «методы на основе погружения поверхности») поверхность, на которой решается уравнение, считается погруженной в трехмерное пространство. Геометрия поверхности задается неявно, как поверхность уровня ноль заданной функции, определенной в трехмерном пространстве. Суть метода заключается в том, что заданное на поверхности уровня ноль уравнение продолжается во все пространство (то есть на все поверхности уровня) так, что решение этого уравнения совпадает с решением
исходного уравнения на поверхности. Далее «продолженное» уравнение аппроксимируется во всей трехмерной области «стандартным» методом конечных элементов. След трехмерного решения на исходной поверхности дает решение исходной задачи. Метод этого типа изначально был предложен в [40, 41].
Второй класс методов относится к классу методов типа «диффузной границы», «diffuse interface approach» [42]. Этот класс методов близок к рассмотренному ранее. Общая идея сохраняется, однако способ продолжения уравнения на поверхности в трехмерное пространство отличается от рассмотренного ранее. Метод основан на том, что коэффициенты уравнения рассматриваются как пространственные обобщенные дельта-функции (меры Дирака) с носителем на поверхности. Замена этих распределений на сходящиеся к ним элементы дельта-образных последовательностей позволяет построить задачу во вмещающей среде, след которой сходится к решению задачи на поверхности.
Еще один подход рассмотрен в работе ([43], см. также ссылки там). В соответствии с ним уравнение на поверхности не модифицируется, но для аппроксимации его решения используются конечно-элементные базисные функции, заданные на трехмерной сетке во вмещающей среде.
Наконец, удобной альтернативой указанным классам методов можно считать метод на основе проекции ближайшей точки, предложенный в [44]. Его принципиальной особенностью является способ описания поверхности, который естественно позволяет рассматривать поверхности с краем.
Библиографический список
1. Гидравлический разрыв карбонатных пластов / В.Г. Са-лимов, Н.Г. Ибрагимов, А.В. Насыбуллин, О.В. Салимов. - М.: Нефтяное хозяйство, 2013. - 471 с.
2. Экономидес М., Олини Р., Валько П. Унифицированный дизайн гидроразрыва пласта. От теории к практике / Институт компьютерных исследований. - М., 2007. - 236 с.
3. Michael J. Economides, Kenneth G. Nolte (Eds.) Reservoir Stimulation. - Wiley, 2000. - 856 p.
4. Coussy O. Poromechanics, John Wiley and Sons, 2004.
5. Chipot M., Luskin M. The compressible Reynolds lubrication equation // IMA Preprint Series. - 1986. - No. 232.
6. Computer simulation of hydraulic fractures / J. Adachi, E. Siebrits, A. Peirce, J. Desroches // Int. J. Rock Mech. Min. Sci. -2007. - Vol. 44. - P. 739-757.
7. Дубровин Б.А., Новиков С.П., Фоменко А.Т. Современная геометрия. Методы и приложения. - М.: Наука, 1986. -760 c.
8. Eskin D., Miller M.J. A model of non-Newtonian slurry flow in a fracture // Powder Technology. - 2008. - Vol. 182. - P. 313-322.
9. Clifton R.J., Wang J.-J. Multiple fluids, proppant transport, and thermal effects in threedimensional simulation of hydraulic fracturing // SPE Paper. - 1988. - 18198.
10. Clifton R.J., Brown U., Wang J.-J. Modeling of Poroelastic Effects in Hydraulic Fracturing // SPE Paper. - 1991. - 21871.
11. Черепанов Г.П. Механика хрупкого разрушения. -М.: Наука, 1974. - 640 с.
Заключение
В настоящей работе описана математическая модель развития трещины гидроразрыва пласта в полной связанной трехмерной поставке. Модель включает в себя группу уравнений пороупругости, описывающих связанные механические и фильтрационные процессы во вмещающей трещину среде, уравнение смазочного слоя, описывающее течение в трещине, группу условий согласования на границе «трещина-пласт» и граничные условия на фронте трещины и флюида в ней. С геометрической точки зрения моделью трещины является произвольная гладкая поверхность с краем. Подробно описаны основные допущения, использованные при выводе уравнений модели. Основное назначение предложенной модели - согласованное описание развития трещины ГРП в достаточно общей постановке, пригодное для использования с применением современных вычислительных подходов. В связи с этим в работе дан краткий обзор вычислительных алгоритмов, пригодных для практической реализации модели.
Благодарности
Работа выполнена при поддержке Российского научного фонда, проект № 15-11-00021.
Acknowledgments
The work has been carried out with the financial support of the Russian Science Foundation, Project No. 15-11-00021.
12. Качанов Л.М. Основы механики разрушения. - М.: Наука, 1974. - 312 с.
13. Рамазанов М.М., Критский Б.В., Савенков Е.Б. Формулировка J-интеграла для модели пороупругой среды Био // ИФЖ. - 2018. - T. 91, № 6.
14. Линьков А.М. Уравнение скорости и его применение для решения некорректных задач о гидроразрыве // Докл. Акад. наук. - 2011. - Т. 439, № 4. - С. 473-475.
15. Detournay E., Peirce A. On the moving boundary conditions for a hydraulic fracture // International Journal of Engineering Science. - 2014. - Vol. 84. - P. 147-155.
16. Garagash D.I., Detournay E. The tip region of a fluid-driven fracture in an elastic medium // ASME J. Appl. Mech. -2000. - Vol. 67. - No. 1. - P. 183-192.
17. Garagash D.I. Propagation of a plane-strain fluid-driven fracture with a fluid lag: early-time solution // Int J Solids Struct. -2006. - Vol. 43. - P. 5811-5835.
18. Linkov A.M. Particle velocity, speed equation and universal asymptotics for efficient modelling of hydraulic fracturing // Journal of Applied Mathematics and Mechanics. -2015. - Vol. 79. - No. 1. - P. 54-63.
19. Lewis R.W., Schrefler B.A. The Finite Element Method in the Static and Dynamic Deformation and Consolidation in Porous Media. - J. Wiley, Chichester, 1998.
20. Kim J., Tchelepi H.A., Juanes R. Stability, Accuracy and Efficiency of Sequential Methods for Coupled Flow and Geomechanics // SPE Paper. - 2009. - 119084.
21. Brezzi F., Fortin M. Mixed and Hybrid Finite Elements Methods. - Springer, 1991.
22. Ern A., Meunier S. A-posteriori error analysis of Euler-Galerkin approximations to coupled elliptic-parabolic problems // ESAIM: M2AN. - 2009. - Vol. 43. - No. 2. - P. 353-375.
23. Murad M.A., Thomee V., Loula A.F.D. Asymptotic Behavior of Semidiscrete Finite-Element Approximations of Biot's Consolidation Problem // SIAM Journal on Numerical Analysis. - 1996. - Vol. 33. - No. 3. - P. 1065-1083.
24. Showalter R.E. Diffusion in Poro-Elastic Media // Journal of Mathematical Analysis and Applications. - 2000. -Vol. 251. - No. 1. - P. 310-340.
25. White J.A., Borja R.I. Stabilized low-order finite elements for coupled solid-deformation/fluid-diffusion and their application to fault zone transients // Comput. Methods Appl. Mech. Engrg. - 2008. - Vol. 197. - No. 49-50. - P. 4353-4366.
26. Commend S., Truty A., Zimmermann Th. Stabilized finite elements applied to elastoplasticity: I. Mixed displacement-pressure formulation // Comput. Methods Appl. Mech. Engrg. -2004. - Vol. 193. - P. 3559-3586.
27. Xia K., Masud A.A stabilized finite element formulation for finite deformation elastoplasticity in geomechanics // Computers and Geotechnics. - 2009. - Vol. 36. - P. 396-405.
28. Chappele D., Bathe K.J. The inf- sup Test // Computers
& Structures. - 1993. - Vol. 47. - No. 4-5. - P. 537-545.
29. Simo J.C., Rifai M.S. A class of mixed assumed strain methods and the method of incompatible modes // Int. J. Numer. Methods Engrg. - 1990. - Vol. 29. - P. 1595-1638.
30. Oliver J., Cervera M., Manzoli O. Strong discontinuities and continuum plasticity models: the strong discontinuity approach // Int. J. Plasticity. - 1999. - Vol. 15. - P. 319-351.
31. A PDE-based fast local level set method / D. Peng, B. Merriman, S. Osher, H. Zhao, M. Kang // Journal of Computational Physics. - 1999. - Vol. 155. - P. 410-438.
32. Chopp D.L. Computing minimal surfaces via level-set curvature flow // Journal of Computational Physics. - 1993. -Vol. 106. - P. 77-91.
33. Mourad H.M., Dolbow J., Harari I. A bubble-stabilized finite element method for Dirichlet constraints on embedded
References
1. Salimov V.G., Ibragimov N.G., Nasybullin A.V., Salimov O.V. Gidravlicheskii razryv karbonatnykh plastov [Hydraulic fracturing of carbonate formations]. Moscow, Neftyanoe khozyaistvo, 2013, 471 p.
2. Economides M.J., Oligney R.E., Valko P. Unified fracture design: bridging the gap between theory and practice alvin. Texas, Orsa Press, 2001, 200 p.
3. Michael J. Economides, Kenneth G. Nolte (Eds.) Reservoir Stimulation. Wiley; 3 edition, 2000, 856 p.
4. Coussy O. Poromechanics. John Wiley and Sons, 2004.
5. Chipot M., Luskin M. The compressible Reynolds lubrication equation. IMA Preprint Series no. 232, 1986.
6. Adachi J., Siebrits E., Peirce A., Desroches J. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 2007. vol. 44. pp. 739-757.
7. Dubrovin B.A., Novikov S.P., Fomenko A.T. Sovremennaya geometriya. Metody i prilozheniya [Modern geometry. Methods and applications] Second edition. Moscow, Nauka, 1986, 760 p.
8. Eskin D., Miller M.J. A model of non-Newtonian slurry flow in a fracture. Powder Technology. 2008. vol. 182. pp. 313-322.
interfaces // International Journal for Numerical Methods in Engineering. - 2007. - Vol. 69. - No. 4. - P. 772-793.
34. Dolbow J.E, Franca L.P. Residual-free bubbles for embedded Dirichlet problems // Computer Methods in Applied Mechanics and Engineering. - 2008. - Vol. 197. - P. 3751-3759.
35. Bechet E., Moes N., Wohlmuth B. A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method // Int. J. Numer. Meth. Engng. - 2009. -Vol. 78. - P. 931-954.
36. Hautefeuille M., Annavarapu C., Dolbow, J.E. Robust imposition of Dirichlet boundary conditions on embedded surfaces // Int. J. Numer. Meth. Engng. - 2012. - Vol. 90. - P. 40-64.
37. Nitsche J. Uber ein Variationsprinzip zur Losung von Dirichlet-Problemen bei Verwendung von Teilrümen, die keinen Randbedingungen unterworafen sind // Abh. Math. Sem. Univ. Hamburg. - 1971. - Vol. 36. - P. 9-15.
38. Dolbow J., Harari I. An efficient finite element method for embedded interface problems // International Journal for Numerical Methods in Engineering. - 2009. - Vol. 78. - No. 2. -P. 229-252.
39. Dziuk G., Elliott C.M. Surface finite elements for parabolic equations // J. Comput. Math. - 2007. - Vol. 25. -P. 385-407.
40. Variational problems and partial differential equations on implicit surfaces / M. Bertalmio, L.T. Cheng, S. Osher, G. Sapiro // J. Comput. Phys. - 2001. - Vol. 174. - P. 759-780.
41. Dziuk G., Elliott C.M. An Eulerian approach to transport and diffusion on evolving implicit surfaces // Comput Visual Sci. -2010. - Vol. 13. - P. 17-28.
42. Ratz A., Voigt A. PDEs on Surfaces: A Diffuse Interface Approach // Comm. Math.Sci. - 2006. - Vol. 4. - No. 3. - P. 575-590.
43. Olshanskii M.A., Reusken A., Xu X. A Volume Mesh Finite Element Method for PDEs on Surfaces // European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012), eds. J. Eberhardsteiner et.al. -Vienna, Austria, September 2012. - Р. 10-14.
44. Ruuth S.J., Merriman B. A simple embedding method for solving partial differential equations on surfaces // Journal of Computational Physics. - 2008. - Vol. 227. - P. 1943-1961.
9. Clifton R.J., Wang J.-J. Multiple fluids, proppant transport, and thermal effects in threedimensional simulation of hydraulic fracturing. SPE paper, 1988. 18198,
10. Clifton R.J., Brown U., Wang J.-J. Modeling of poroelastic effects in hydraulic fracturing. SPE Paper, 19912, 1871.
11. Cherepanov G.P. Mehanika hrupkogo razrushenija. [Mechanics of brittle fracture]. Moscow, Nauka, 1974, 600 p.
12. Kachanov L.M. Introduction to continuum damage mechanics. Dordrecht; Boston, Martinus Nijhoff 1986, 135 p.
13. Ramazanov M.M., Kritsky B.V., Savenkov E.B. Formulation of J-integral for Biot poroelastic medium. Journal of Engineering Physics and Thermophysics. 2018, vol. 91, no. 6. (accepted).
14. Linkov A.M. Speed Equation and Its Application for Solving Ill-Posed Problems of Hydraulic Fracturing. Doklady Physics. 2011, vol. 56, no. 8, pp. 436-438.
15. Detournay E., Peirce A. On the moving boundary conditions for a hydraulic fracture. International Journal of Engineering Science. 2014, vol. 84, pp. 147-155.
16. Garagash D.I., Detournay E. The tip region of a fluid-driven fracture in an elastic medium. ASME J. Appl. Mech. 2000, vol. 67, no. 1, pp. 183-192.
17. Garagash D.I. Propagation of a plane-strain fluid-driven fracture with a fluid lag: early-time solution. Int J Solids Struct. 2006, vol. 43, pp. 5811-5835.
18. Linkov A.M. Particle velocity, speed equation and universal asymptotics for efficient modelling of hydraulic fracturing. Journal of Applied Mathematics and Mechanics. 2015, vol. 79, no. 1, pp. 54-63.
19. Lewis R.W., Schrefler B.A. The Finite Element Method in the Static and Dynamic Deformation and Consolidation in Porous Media. J. Wiley, Chichester, 2nd ed., 1998.
20. Kim J., Tchelepi H.A., Juanes R. Stability, Accuracy and efficiency of sequential methods for coupled flow and geomechanics. SPE Paper, 2009, 119084.
21. Brezzi F., Fortin M. Mixed and Hybrid Finite Elements Methods. Springer, 1991.
22. Ern A., Meunier S. A-posteriori error analysis of Euler-Galerkin approximations to coupled elliptic-parabolic problems. ESAIM: M2AN. 2009, vol. 43, no. 2, pp. 353-375.
23. Murad M.A., Thomee V., Loula A.F.D. Asymptotic Behavior of Semidiscrete Finite-Element Approximations of Biot's Consolidation Problem. SIAM Journal on Numerical Analysis. 1996, vol. 33, no. 3, pp. 1065-1083.
24. Showalter R.E. Diffusion in Poro-Elastic Media. Journal of Mathematical Analysis and Applications. 2000, vol. 251, no. 1, pp. 310-340.
25. White J.A., Borja R.I. Stabilized low-order finite elements for coupled solid-deformation/fluid-diffusion and their application to fault zone transients. Comput. Methods Appl. Mech. Engrg. 2008, vol. 197, no. 49-50, pp. 4353-4366.
26. Commend S., Truty A., Zimmermann Th. Stabilized finite elements applied to elastoplasticity: I. Mixed displacement-pressure formulation.Comput. Methods Appl. Mech. Engrg, 2004, vol. 193, pp. 3559-3586.
27. Xia K., Masud A.A stabilized finite element formulation for finite deformation elastoplasticity in geomechanics. Computers and Geotechnics. 2009, vol. 36, pp. 396-405.
28. Chappele D., Bathe K.J. The inf-sup Test. Computers & Structures. 1993, vol. 47, no. 4-5, pp. 537-545.
29. Simo J.C., Rifai M.S. A class of mixed assumed strain methods and the method of incompatible modes. Int. J. Numer. Methods Engrg. 1990, vol. 29, pp. 1595-1638.
30. Oliver J., Cervera M., Manzoli O. Strong discontinuities and continuum plasticity models: the strong discontinuity approach. Int. J. Plasticity. 1999, vol. 15, pp. 319-351.
31. Peng D., Merriman B., Osher S., Zhao H., Kang M. A PDE-based fast local level set method. Journal of Computational Physics. 1999, vol. 155, pp. 410-438.
32. Chopp D.L. Computing minimal surfaces via level-set curvature flow. Journal of Computational Physics. 1993, vol. 106, pp. 77-91.
33. Mourad H.M., Dolbow J., Harari I. A bubble-stabilized finite element method for Dirichlet constraints on embedded interfaces. International Journal for Numerical Methods in Engineering. 2007, vol. 69, no. 4, pp. 772-793.
34. Dolbow J.E, Franca L.P. Residual-free bubbles for embedded Dirichlet problems. Computer Methods in Applied Mechanics and Engineering. 2008, vol. 197, pp. 3751-3759.
35. Bechet E., Moës N., Wohlmuth B. A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method. Int. J. Numer. Meth. Engng, 2009, vol. 78. pp. 931-954.
36. Hautefeuille M., Annavarapu C., Dolbow, J.E. Robust imposition of Dirichlet boundary conditions on embedded surfaces. Int. J. Numer. Meth. Engng, 2012, vol. 90. pp. 40-64.
37. Nitsche J. Uber ein Variationsprinzip zur Losung von Dirichlet-Problemen bei Verwendung von Teilrumen, die keinen Randbedingungen unterworafen sind. Abh. Math. Sem. Univ. Hamburg, 1971, vol. 36, pp. 9-15.
38. Dolbow J., Harari I. An efficient finite element method for embedded interface problems. International Journal for Numerical Methods in Engineering, 2009, vol. 78, no. 2, pp. 229-252.
39. Dziuk G., Elliott C.M. Surface finite elements for parabolic equations. J. Comput. Math, 2007, vol. 25, pp. 385-407.
40. Bertalmio M., Cheng L.T., Osher S., Sapiro G. Variational problems and partial differential equations on implicit surfaces. J. Comput. Phys, 2001, vol. 174, pp. 759-780.
41. Dziuk G., Elliott C.M. An Eulerian approach to transport and diffusion on evolving implicit surfaces. Comput Visual Sci, 2010, vol. 13, pp. 17-28.
42. Ratz A., Voigt A. PDEs on Surfaces: A Diffuse Interface Approach. Comm. Math.Sci, 2006, vol. 4, no. 3, pp. 575-590.
43. Olshanskii M.A., Reusken A., Xu X. A Volume Mesh Finite Element Method for PDEs on Surfaces. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012) J. Eberhardsteiner et.al. (eds.) Vienna, Austria, September 10-14, 2012.
44. Ruuth S.J., Merriman, B. A simple embedding method for solving partial differential equations on surfaces. Journal of Computational Physics, 2008, vol. 227, pp. 1943-1961.
