Full Text

Введение В настоящей работе представлено описание математической модели и комплекса вычислительных алгоритмов для анализа динамики развития крупномасштабной флюидонаполненной трещины в пороупругой среде. Настоящая работа является продолжением работы авторов [1], которая посвящена детальному описанию математической модели. Вопросы описания флюидонаполнных трещин возникают в целом ряде приложений, см., например, [2, 3]. В частности, важной прикладной задачей является прогноз динамики трещины гидравлического разрыва пласта. Технология гидроразрыва пласта (ГРП) является одним из самых распространенных в настоящее время методов увеличения нефтеотдачи [4, 5]. Создание трещины гидроразрыва является экономически дорогим и технически сложным процессом. По этой причине важно обоснование применимости технологии в случае каждой отдельной скважины. Для этого необходимо оценить возможность создания трещины с заданными параметрами с учетом влияния свойств напряженно-деформированной среды на процесс развития трещины, ее конечные геометрические и фильтрационные характеристики, в том числе траекторию. Полноценное решение этих задач в настоящее время доступно только средствами математического моделирования. По этой причине анализу и математическому моделированию процесса развития флюидонаполненных трещин уделяется большое внимание, см., например, работы [6-14] и ссылки в них. Отметим, что используемые на практике модели развития трещин ГРП являются существенно упрощенными. Чаще всего они основаны на применении двумерных и псевдотрехмерных постановок, используются априорные допущения об ориентации и форме трещины, см. [4, 5]. В большинстве случаев они не учитывают сопровождающие развитие трещины фильтрационные процессы, для учета утечки флюида в трещине в пласт используются эмпирические модели. С алгоритмической точки зрения методы решения задач развития трещины в пороупругой среде (см., например, [6-10,12]) в подавляющем большинстве случаев либо используют достаточно сложные математические модели (однако соответствующие вычислительные алгоритмы не позволяют исследовать полный спектр актуальных постановок, поэтому часто рассматриваются упругие постановки с однородным или кусочно-однородным распределением свойств вмещающей трещину среды), либо основаны на применении нескольких независимых алгоритмов или пакетов программ, отвечающих за решение частных задач. В отличие от указанных подходов, используемая в настоящей работе математическая модель является трехмерной, полностью связанной и включает в себя группы уравнений пороупругости и течения в трещине, дополненные физически обоснованными критериями ее развития. В качестве модели вмещающей среды используется связанная модель Био для описания возмущений начального поля напряжений регионального масштаба, вызванных присутствием и развитием трещины и эффектами фильтрации. Для описания течения флюида в трещине использованы уравнения гидродинамики в приближении смазочного слоя. Срединная поверхность трещины является произвольной гладкой поверхностью с краем. Для описания процесса развития трещины использованы физически обоснованные критерии разрушения. Детальное описание модели можно найти в работе авторов [1], кратко она будет представлена в соответствующем разделе ниже. Вычислительный алгоритм основан на методе конечных элементов. Уравнения пороупругости в среде решаются с помощью оригинального варианта метода X-FEM [151-8], позволяющего учитывать наличие в расчетной области трещин, геометрия которых не согласована с геометрией расчетной сетки. Для представления срединной поверхности трещины при этом используется метод проекции ближайшей точки, с помощью которого также осуществляется расчет эволюции трещины [19] и интегрирование уравнений гидродинамики на ее поверхности [20-22]. Алгоритм решения полной задачи основывается на итерационном связывании групп уравнений, описывающих частные модели. Предложенный алгоритм является полностью эйлеровым и использует для решения задачи единую пространственную расчетную сетку. Описываемый подход был реализован в виде программного комплекса HFrac3D++ [23], разработанного в ИПМ им. М.В. Келдыша РАН и предназначенного для решения задач пороупругости, в том числе при наличии флюидонаполненных трещин. В заключительной части работы приводятся результаты демонстрационных расчетов для нескольких реалистичных постановок, показывающих применимость разработанной методики для рассматриваемого класса задач. 1. Математическая модель Описание динамики флюидонаполненной трещины в пороупругой среде сводится к решению сложной связанной задачи, включающей в себя: • систему уравнений пороупругости, описывающую эволюцию напряженно-деформированного состояния среды и полей давления флюида в ней в ходе развития трещины; • геометрическую модель представления и описания эволюции срединной поверхности трещины; • уравнение течения жидкости в трещине; • условия согласования между полями давления в трещине и в среде, потоками массы, импульса (и, в общем случае, энергии), а также кинематические условия, связывающие раскрытие трещины и перемещение точек пласта; • механические условия развития трещины, определяющие направление ее развития в каждой точке ее фронта. Далее кратко приведем уравнения модели для изотермического случая. Система уравнений пороупругости в трехмерной постановке [24-27], которая описывает состояние вмещающей трещину среды , имеет вид (1) где - плотность насыщенной флюидом вмещающей среды; - ускорение свободного падения; - тензор напряжений, который в общем случае определяется как ; - давление флюида, - тензор деформаций; - поле перемещений точек вмещающей среды; - симметричный тензор упругих коэффициентов; - тензор Био, - масса флюида в единице объема насыщенной среды (флюидосодержание); - плотность флюида; - пористость; - массовая плотность источников; для какой-либо величины , где индексом «0» обозначены опорные значения величины. Плотность потока массы определяется законом Дарси , - симметричный тензор коэффициентов проницаемости, - динамическая вязкость флюида. Геометрическая модель трещины предполагает, что срединная поверхность трещины является гладкой, односвязной поверхностью с краем («фронтом трещины»). Толщина (раскрытие) трещины является функцией точки , . Поверхность является срединной поверхностью трещины в том смысле, что боковые поверхности реальной, «физической» трещины находятся на расстоянии от срединной поверхности , . Здесь - нормаль к поверхности в точке , направленная в сторону . Естественными ограничениями на вид функции является ее неотрицательность и равенство нулю на крае . Второе допущение отражает тот факт, что в ходе развития трещины ее уже образованная поверхность не меняется, т.е. при имеем , где - время, - поверхность в момент времени . Помимо этого считается, что эволюция поверхности является гладкой. В этом случае семейство поверхностей может быть представлено как объединение начальной поверхности и следа «движения» ее края , Считается, что в каждый момент времени на линии задано векторное поле скорости описывающее ее эволюцию. Таким образом, движение точки границы описывается уравнением . Направление развития поверхности известно только в точках ее границы (фронта трещины) и определяется с помощью соответствующих критериев разрушения. Модель течения в трещине основана на приближении Рейнольдса и связывает раскрытие и давление в трещине при наличии источников массы [28]: (2) где - скорость течения флюида в трещине; - эффективная проницаемость трещины; обозначает двумерный поверхностный оператор градиента. В простейшем случае ньютоновской жидкости В общем случае неньютоновской жидкости Величина представляет собой мощность источников, описывающих поток флюида через скважину и утечку в пласт. Условия согласования, заданные на срединной поверхности трещины, состоят из двух групп соотношений. К первой группе относятся кинематические условия, связывающие перемещения боковых поверхностей трещины и ее раскрытие где , Ко второй группе уравнений относятся условия, выражающие непрерывность консервативных величин. Динамические условия согласования имеют вид условий непрерывности нормальной компоненты тензора напряжений на боковых поверхностях трещины: Для «гидродинамической» части задачи, включающей в себя уравнения фильтрации и уравнение течения в трещине, условия согласования связывают давление во вмещающей трещину среде и соответствующие потоки массы. Они будут иметь различный вид в зависимости от допущений, которым удовлетворяют условия течения. В простейшем случае идеального гидродинамического контакта между трещиной и вмещающей ее средой они имеют вид Условия непрерывности потока массы имеют вид Таким образом, величина потока из среды в трещину определяет правую часть (источник) закона сохранения массы флюида в трещине, а давление в трещине играет роль граничного условия для давления на боковых поверхностях трещины. Критерий разрушения, который определяет направление развития фронта трещины в его точках, формулируется с помощью векторного J-интеграла Черепанова - Райса [29, 30] для пороупругой среды. В частном случае плоской (двумерной) трещины выражение для векторного J-интеграла имеет вид где - двумерная область, окружающая кончик трещины; - вектор единичной нормали к ; - пористость; - свободная энергия скелета, определенная как В том случае, когда трещина развивается в направлении традиционный (скалярный) J-интеграл определяется выражением . Векторное выражение J-интеграла позволяет сформулировать следующий эмпирический критерий развития трещины: трещина развивается в направлении , таком, что , при этом для произвольного единичного вектора , где - критическое значение энергии, которое, в общем случае, может зависеть от направления. Для однородной изотропной среды и направление развития трещины совпадает с направлением вектора . Аналогично записываются соответствующие выражения для трехмерного случая. 2. Вычислительный алгоритм Предлагаемый в работе численный алгоритм решения задачи о динамике развития флюидонаполненой трещины в пороупругой среде основывается на итерационном связывании решений отдельных групп уравнений математической модели. Ниже сначала рассматриваются алгоритмы для решения частных задач, далее - алгоритм решения полной связанной задачи. В качестве базового алгоритма решения задачи во вмещающей трещину среде используется оригинальный вариант «расширенного» метода конечных элементов X-FEM/CP, позволяющий учитывать наличие в расчетной области трещин, геометрия которых не согласована с геометрией расчетной сетки. Для чисто упругой среды с трещиной детали алгоритма представлены в работе [15], в пороупругом случае аппроксимации строятся аналогичным образом. На первом этапе метода решения для системы (1)-(2) записывается слабая (вариационная) постановка задачи, условия на трещине в которой задаются с использованием метода штрафа в терминах множителей Лагранжа: (3) (4) (5) Здесь где - соответствующие пространства функций, которые могут быть охарактеризованы в терминах пространств Соболева нужной гладкости; - раскрытие трещины, определяет поток из вмещающей среды в трещину, и - скалярные произведения в и соответственно; и - параметры Био; - объемная деформация среды. Уравнения (3)-(4), описывающие напряженно-деформированное состояние вмещающей трещину среды, далее используются непосредственно. Уравнение для описания течения в трещине (5) применяются в несколько ином виде, основанном на применении метода проекции ближайшей точки для построения продолжения соответствующего уравнения из трещины во вмещающее ее пространство (однако общая структура задачи при этом сохраняется). На втором этапе метода X-FEM производится выбор дискретных пространств и аппроксимации задачи по пространству и времени, после чего соответствующая система уравнений для (3)-(4) в матричном виде может быть представлена как (6) где , , , , и - конечномерные аппроксимации соответствующих билинейных форм и правой части; и - зависящие от времени векторы узловых значений конечно-элементных аппроксимаций поля перемещений и давления; - шаг по времени, ^ обозначает величины на новом временном слое. Далее данная система решается подходящим методом, основанным на итерациях между группами уравнений теории упругости и фильтрации [31, 32]. Для решения возникающих на каждой итерации больших разреженных систем линейных алгебраических уравнений используется метод BiCGStab c предобуславливателем на основе неполного LU разложения [33]. Отметим, что, вообще говоря, исходная формулировка метода X-FEM не зависит от способа представления поверхности трещины, однако удобство представления является основой его эффективного применения. По этой причине единым способом представления срединной поверхности трещины был выбран метод проекции ближайшей точки. В рамках этого метода поверхность описывается неявно, путем задания проектора точки пространства на срединную поверхность трещины. Рассмотрим гладкую поверхность с краем, вложенную в трехмерную область . Пусть для произвольной точки , точка является ближайшей к ней точкой на поверхности . Точку будем называть проекцией точки на поверхность в смысле наименьшего расстояния, а соответствующий оператор проектирования будем обозначать , . В случае если поверхность и ее край являются гладкими, а область - «достаточно маленькая», то оператор однозначно определен [34] . Иными словами, каждая точка области однозначно проецируется в единственную точку на поверхности . Проектор однозначно описывает поверхность как множество своих неподвижных точек: . Проектор может быть использован как в континуальной, так и в дискретной постановке. В последнем случае он задается своими значениями в узлах расчетной сетки. Подробный алгоритм расчета эволюции зависящего от времени проектора в рамках геометрической модели трещины из раздела 1 можно найти в работе [20]. Метод на основе проекции ближайшей точки используется также при решении уравнения течения в трещине (5), см. [30]. Для этого уравнение, заданное на поверхности с краем, продолжается во вмещающее поверхность пространство. Для построения продолжения уравнения на поверхности в пространство используется оператор , который для произвольной заданной на поверхности функции определяется как . Тогда для определенного на поверхности уравнения [21, 22] имеем где - некоторый «цилиндр» вокруг . След решения этого уравнения на является решением исходного уравнения на поверхности, см. [34-37]. Отметим, что, в отличие от этих исследований, в настоящей работе используется предложенный в [20-22] конечно-элементный вариант метода проекции ближайшей точки. Таким образом, для решения как задач в объеме, так и на поверхности используются единая пространственная расчетная сетка, а также единое представление поверхности в рамках частных алгоритмов решения подзадач. Алгоритм решения полной связанной задачи в сокращенном виде может быть представлен следующим образом. Задача решается в области , в которой введена расчетная сетка , используемая для решения задачи пороупругости и течения в трещине. В каждый момент времени решение задачи описывается следующим набором переменных где - поле перемещений, - давление в области, - раскрытие трещины, - давление в трещине, - оператор проекции ближайшей точки, однозначно определяющий геометрию срединной поверхности трещины, - скорость движения фронта трещины. Рассмотрим задачу на одном слое по времени . Будем использовать обозначения . Расчет одного временного слоя соответствует переходу от набора переменных к набору переменных и состоит из следующих шагов: 1. Производится расчет дискретного оператора проекции ближайшей точки для новой геометрии трещины . Далее осуществляется определение типов узлов конечных элементов для алгоритма метода X-FEM, расчет параметров, необходимых для определения разрывных и сингулярных базисных функций метода. 2. Строится продолжение поля давления в трещине , полей и давления и перемещений в объеме с учетом новой геометрии трещины . 3. Производится инициализация внутренних итераций: . 4. Рассчитывается решение уравнения для поля перемещений при заданных полях давления в области и трещине 5. Рассчитывается раскрытие трещины по определенному полю перемещений . Далее производится решение уравнения для определения поля давления в трещине. 6. Производится решение уравнения для давления при заданном поле перемещений для определения поля давления . 7. Выполняется проверка критерия остановки итераций. В случае его невыполнения - переход на шаг 4, в противном случае - на шаг 8. 8. Определяются направления развития и скорости трещины в точках ее фронта на основе заданного критерия и скорости развития фронта трещины. 9. Окончание расчета временного слоя. При минимальных изменениях описанный выше алгоритм позволяет проводить расчеты для нескольких трещин в случае выполнения единственного условия: трещины расположены достаточно далеко друг от друга. С алгоритмической точки зрения это означает, что содержащие трещины сеточные области, в которых решается уравнение смазочного слоя и определены дополнительные степени свободы метода X-FEM, не пересекаются (т.е. отделены по крайне мере одним слоем конечных элементов). Формально это позволяет использовать описанные выше алгоритмы и структуры данных для случая одной трещины в каждой из таких сеточных областей. 3. Примеры расчетов В разделе приводятся результаты расчетов, которые демонстрируют возможности разработанных алгоритмов и программной реализации для решения связанных задач с учетом пороупругого поведения среды, течения в трещине и ее эволюции. Рассматриваются примеры со стационарной трещиной и эволюционирующей трещиной в среде с неоднородным распределением упругих и фильтрационно-емкостных свойств в естественном поле напряжений. Стационарная трещина. В этом расчете геометрия трещины (форма ее срединной поверхности) считается заданной и не меняется с течением времени. Вместе с тем задача решается в полностью связанной постановке - с учетом течения в трещине и эволюции полей давлений и напряжений в среде. Расчетная область представляет собой куб со стороной Для расчетов использовалась равномерная сетка из тетраэдров. Узлы сетки образуют кубическую решетку с . Таким образом, шаг сетки узлов м. Всего в сетке узлов и конечных элементов. Шаг по времени составлял c (15 мин), полное время расчета c (5 ч). Расчетная область разделена на две части вертикальной плоскостью, проходящей через ее центр перпендикулярно оси . Свойства среды в каждой из подобластей приведены в табл. 1, где - модуль Юнга, - коэффициент Пуассона, и - параметры Био, - коэффициент проницаемости. Вязкость и сжимаемость флюида задавались как и Стационарная трещина имеет форму плоского диска радиусом м, расположенного горизонтально в центре расчетной области. Распределение давления в трещине определяется решением полной системы уравнений, с учетом группы уравнений пороупругости и течения в трещине. В начальный момент времени давление внутри трещины принималось равным Па. В центральной части трещины давление считалось заданным в виде граничного условия Па. На границе области давление жидкости ставилось равным , нормальные компоненты поля перемещений на гранях задавались равными нулю. В начальный момент времени давление жидкости в области также равнялось . Таблица 1 Свойства модельной среды в расчете со стационарной трещиной Table 1 Properties of the model medium in the simulation with a stationary fracture На рис. 1 показано изменение давления в трещине и вмещающей ее среде с течением времени. Видно выполнение условий непрерывности давлений в среде и на срединной поверхности трещины. Распределение давления (равно как и раскрытие) в трещине непостоянно и несимметрично - это связано с отличием упругих и фильтрационных модулей среды в различных подобластях. На рис. 2 показано распределение компонент тензора напряжений на последнем шаге по времени, также хорошо заметно качественно отличающееся поведение полей напряжений, связанное с существенным отличием упругих свойств. Характер деформации боковых поверхностей трещины более подробно демонстрируют рис. 3 и 4. На этих рисунках трещина показана не как плоская срединная поверхность, а как трехмерный объект. Его толщина в различных точках срединной поверхности определяется раскрытием трещины в соответствующих точках (показано в измененном масштабе). На рис. 3 в различных ракурсах показано поле перемещений и раскрытие трещины . Наконец, на рис. 4 преставлена деформированная поверхность трещины в различных ракурсах с распределением давления в трещине. Рис. 1. Распределение давления в среде и трещине для Fig. 1. Pressure distribution in the medium and fracture at Рис. 2. Распределение давления в трещине и компонент тензора упругих напряжений при Fig. 2. Pressure distribution in fracture and elastic stress tensor components at Рис. 3. Деформированная геометрия трещины, поле перемещений и раскрытие трещины при Fig. 3. Deformed fracture geometry, displacement field and crack opening at Рис. 4. Деформированная геометрия трещины и распределение поля давления в ней при Fig. 4. Deformed fracture geometry and pressure field at Эволюционирующая трещина в естественном поле напряжений. Рассматривается задача об эволюции изначально плоской трещины в естественном поле напряжений с учетом неоднородности прочностных свойств среды. Это означает, что в различных подобластях вмещающей трещину среды коэффициенты трещиностойкости имеют разное значение, что приводит к тому, что скорость распространения фронта трещины в них также отличается. Развитие срединной поверхности трещины определяется как распределением поля напряжений, так и значением коэффициента трещиностойкости. Расчетная область представляет собой параллелепипед , где Для расчетов использовалась равномерная сетка из тетраэдров. Узлы сетки образуют кубическую решетку с и . Таким образом, шаг сетки узлов м. Всего в сетке узлов и конечных элементов. Шаг по времени составлял c (15 мин), полное время расчета c (6,5 ч). Расчетная область делилась на две части вертикальной плоскостью, проходящей через ее центр перпендикулярно оси . Свойства среды в каждой из подобластей приведены в табл. 2, где - модуль Юнга, - коэффициент Пуассона, и - параметры Био, - коэффициент проницаемости. Вязкость и сжимаемость флюида задавались как и Коэффициенты трещиностойкости, определяющие скорость развития трещины, принимают различные значения в подобластях, симметричных плоскости . Таблица 2 Свойства модельной среды в расчете с эволюционирующей трещиной Table 2 Properties of the model medium in the simulation with an evolving fracture В начальный момент времени зародыш трещины имеет форму плоского диска радиусом м, расположенного горизонтально в центре расчетной области. Распределение давления в трещине определяется решением полной системы уравнений, с учетом группы уравнений пороупругости и течения в трещине. В начальный момент времени давление внутри трещины принималось Па. В центральной части трещины давление считалось заданным в виде граничного условия Па. На границе области давление жидкости ставилось равным , нормальные компоненты поля перемещений на гранях задавались равными нулю. В начальный момент времени давление жидкости в области также равнялось . На верхней и нижней границах области прикладывались сдвиговые напряжения, имитирующие естественное поле напряжений в пластовых условиях. На рис. 5 показано изменение давления в трещине и вмещающей ее среде с течением времени. Результаты приведены для моментов времени . Видно выполнение условий непрерывности давлений в среде и на срединной поверхности трещины. Распределение давления (равно как и раскрытие) в трещине непостоянно и несимметрично - это связано с отличием упругих и фильтрационных модулей среды в различных подобластях. Более того, явно заметен несимметричный характер развития трещины - в одном из направлений она развивается быстрее, чем в другом. Более детально динамика развития срединной поверхности в трещины и распределения давления в ней показана на рис. 6. Хорошо заметны небольшие изломы на фронте трещины, которые образуются в местах разрыва прочностных свойств среды. Заключение Работа посвящена описанию методов математического моделирования динамики развития флюидонаполненной трещины в пороупругой среде в рамках полной связанной трехмерной постановки. Модель включает в себя группу уравнений пороупругости, описывающих связанные механические и фильтрационные процессы во вмещающей трещину среде, уравнение смазочного слоя, описывающее течение в трещине, группу условий согласования на границе «трещина - пласт» и граничные условия на фронте трещины и флюида в ней. С геометрической точки зрения моделью трещины является произвольная гладкая поверхность с краем. Вычислительный алгоритм основывается на итеративном связывании решения отдельных групп уравнений математической модели, каждая из которых представлена в работе. Уравнения пороупругости в среде решаются с помощью оригинального варианта расширенного метода конечных элементов X-FEM, позволяющего учитывать наличие в расчетной области трещин, геометрия которых не согласована с геометрией расчетной сетки. Для представления срединной поверхности трещины при этом используется метод проекции ближайшей точки, с помощью которого также осуществляется расчет эволюции трещины и интегрирование уравнений гидродинамики на ее поверхности. Предложенный алгоритм является полностью эйлеровым и использует для решения задачи единую пространственную расчетную сетку. Опыт программной реализации предложенного алгоритма показал его эффективность и надежность, прежде всего для решения задач с заранее неизвестной геометрией срединной поверхности трещины. Обобщение алгоритма на случай нескольких трещин может быть выполнено формально и не представляет существенных сложностей. Рис. 5. Распределение давления в среде и трещине для Fig. 5. Pressure distribution in the medium and fracture at Рис. 6. Эволюция срединной поверхности трещины и давления в ней для Fig. 6. Evolution of the fracture mid-surface and its pressure at Дальнейшее развитие алгоритма и его программной реализации будет связано с усложнением физической модели. В частности, для решения содержательных прикладных задач анализа развития трещин гидроразрыва пласта необходим учет неньютоновской реологии флюида внутри трещины. Помимо этого, существенный интерес представляют задачи развития трещине в среде с анизотропным распределением прочностных свойств.

About the authors

V. E Borisov

Keldysh Institute of Applied Mathematics RAS

A. V Ivanov

Keldysh Institute of Applied Mathematics RAS

B. V Kritsky

Keldysh Institute of Applied Mathematics RAS

E. B Savenkov

Keldysh Institute of Applied Mathematics RAS

References

Statistics

Views

Abstract - 256

PDF (Russian) - 156