CC BY
42
УДК 519.6
DOI: 10.53815/20726759_2021_13_3_90
И. М. Перепечкин
Московский физико-технический институт (национальный исследовательский университет)
Численное моделирование ГРП в постановке Planar3D
Гидроразрыв пласта является одним из основных мероприятий для интенсифики-цин добычи и увеличения коэффициента извлечения нефти. Для повышения предсказуемости операции ГРП и предупреждения осложнений требуется точный и быстрый инструмент. Рассматриваемый в данной статье программный комплекс, заключающий в себе быструю реализацию модели ГРП класса Planar 3D, предлагается в качестве подобного инструмента. Постановка Planar3D обладает необходимой физической достоверностью для точного предсказания геометрии плоской трещины, а применение нового подхода к построению численной схемы для реализации Planar3D модели позволяет добиться высокой скорости вычислений.
Построению численной схемы и численного алгоритма оценки критерия продвижения трещины посвящена данная статья.
Ключевые слова: ГРП, Planar3D, численное моделирование.
I. M. Perepechkin
Moscow institute of physics and technology
Numerical simulation of hydraulic fracturing in the Planar3D formulation
Hydraulic fracturing is one of the main measures for intensifying production and increasing the rate of oil recovery. Accurate and fast tools are required to increase the predictability of fracturing operations and prevent complications. The software package considered in this article, which includes the rapid implementation of the Planar3D hydraulic fracturing model, is proposed as such a tool. The Planar3D formulation possesses the necessary physical reliability to accurately predict the geometry of a plane fracture, and the use of a new approach to constructing a numerical scheme for implementing the Planar 3D model allows achieving high computational speed.
This article is devoted to the construction of a numerical scheme and a numerical algorithm for evaluating the criterion of crack propagation.
Key words: Hydraulic fracturing, Planar3D, numerical modeling.
1. Введение
На настоящий момент существует иерархия математических моделей. Данные модели описывают реакцию среды, течение жидкости в трещине и другие физические процессы, происходящие при ГРП с различной степенью физической достоверности. Простые аналитические модели способны предсказывать эволюцию трещины ГРП в частных случаях, самые подробные - ветвление трехмерных трещин. Увеличение подробности модели приводит к усложнению её реализации и осмысления. Начиная с некоторого уровня подробности модели, ГРП требуют привлечения методов численного моделирования.
© Перепечкин И. М., 2021
(с) Федеральное государственное автономное образовательное учреждение высшего образования
«Московский физико-технический институт (национальный исследовательский университет)», 2021
Всего модели, описывающие развитие трещины ГРП, можно разделить на четыре класса: аналитические, Pseudo3D, Planar3D и Full3D.
Двумерные аналитические модели используют решение, показанное в работе [4] для плоско деформированного состояния трещины в линейно-упругой среде. В работах [5, 6] предложена основопологающая аналитическая модель развития трещины ГРП в длину. Предположения модели следующие: трещина ограничена по высоте и подвергающейся горизонтальной плоской деформации. Предложенная модель известна как модель Христиановича-Гертсма-де Клерка (KGD). Предсказания модели KGD являются физически достоверными в случае, когда высота трещины много больше её длины, это выполняется в самом начале роста трещины от интервала перфорации.
Радиальная модель [7] описывает рост трещины ГРП в случае отсутствия контраста горных напряжений в плоскости трещины. На практике такое происходит при образовании трещины ГРП на малых глубинах, когда вертикальное составляющее тензора напряжений меньше горизонтальных составляющих. Трещина при этом развивается в горизонтальной плоскости, и горное давление, действующее на трещину, можно считать однородным. Радиальная модель также может использоваться для описания эволюции трещины в начальный момент времени при её инициализации из точечной перфорации.
Первой моделью, описывающей развитие трещины, длина которой много больше её высоты, является модель Перкинса-Керна-Нордгрена (PKN) [8, 9]. Модель описывает течение жидкости в трещине эллиптического сечения, рост трещины в длину и её раскрытие с учетом эффекта утечек. Трещина считается фиксированной по высоте.
Моделирование, позволяющее описывать вертикальных рост трещины в условиях высокой изменчивости минимальных сжимающих напряжений по пластам, может выполняться на моделях класса Pseudo3D и Planar3D. Модели Pseudo3D [10] можно разделить на два типа - параметрические (Lumped-P3D) [11,12] и кусочно заданные (Cell-based-P3D) [16]. В первой трещина описывается двумя сочлененными эллипсоидами, а во второй - набором связанных ячеек, каждая их которых характеризуется раскрытием. Все модели класса Pseudo3D отличает то, что они основаны на разбиении расчетной области на элементы, в пределах которых происходит решение аналитической модели, а между собой элементы связаны уравнением баланса массы.
Модели класса Pseudo3D обычно являются численными моделями, однако существует аналитическая модель класса Pseudo3D, основанная на получении приближенного численного решения в безразмерных параметрах [14].
Отличием Planar3D [15-18] модели от Pseudo3D является то, что она подразумевает решение трехмерной геомеханической задачи нахождения реакции среды на трещину и решения двумерного уравнения течения жидкости в трещине. При таком подходе отсутствует какое-либо выделенное направление, вдоль которого решаются упрощенные уравнения, как вертикальное направление в Pseudo3D модели. Planar3D модели считаются очень точными, но они требовательны к вычислительным ресурсам [19].
В полностью трехмерных (Full3D) моделях отсутствуют ограничения на геометрию трещины [20, 21], что делает их наиболее физически достоверными. Однако высокие вычислительные затраты ограничивают использование подобных моделей в инженерной практике, Full3D модели используются преимущественно в научных расчетах.
Классификация моделей на аналитические, Pseudo3D, Planar3D и Full3D главным образом определяет ограничение на геометрию трещины. Однако, кроме этого, модели можно разделить на то, с какой степенью физической достоверности они описывают различные эффекты, происходящие при гидроразрыве. С разной физической достоверностью могут быть описаны реакция среды, утечки жидкости в пласт, течение жидости в трещине и транспорт пропанта, разрушение породы и продвижение фронта трещины. Реакция среды обычно рассматривается в приближении однородной линейно упругой среды. Однако есть модели, рассматривающие эффекты реакции с учетом пластичности [22], пороупруго-сти [23] или термических эффектов [24].
С точки зрения баланса между физической достоверностью и скоростью вычислений, модель класса Planar3D является применимой в инженерной практике. Однако нерешенным вопросом остается численная реализация модели Planar 3D, обладающая скоростью расчета, близкой к продвинутым вариантам Pseudo3D модели, но обладающей куда большей точностью [25]. Разработке быстрой реализации Planar3D модели посвещана данная статья.
2. Модель Planar 3D
2.1. Положения и упрощения
В моделе Planar3D, рассматриваемой в данной статье, принимаются следующие положения и упрощения:
1) Считается, что гидроразрыв производится в упругой среде. Эффектами пороупруго-сти и пластической деформацией пренебрегаем. Трещина принимается за плоскую, и соответственно, она не ветвится и не изгибается. Плоскость, в которой расчет трещина выбирается такой, чтобы минимальные сжимающие напряжения в точке инициализации трещины были перпендикулярны к нормали плоскости.
2) Течение жидкости ГРП в трещине рассматривается в приближении смазочного слоя: предполагается, что раскрытие трещины много меньше чем линейные размеры трещины и не учитывается течение жидкости в направлениях, перпендикулярных плоскости стенок трещины.
3) Утечки в пласт моделируются, используя результат решения одномерной нестационарной линейной фильтрации в пористую среду.
4) Критерием разрушения породы на кончике трещины принимается силовой критерий
5) Транспорт пропанта моделируется в предположении теории суспензий. В описываемой модели учитывается бриджинг двух видов: бриджинг, связанный с превышением критической концентрации пропанта, и бриджинг, связанный с геометрическими размерами частиц пропанта.
2.2. Модель реакции среды
Принимается, что процессы, происходящие в упругой среде, на порядки быстрее процессов переноса жидкости и пропанат в трещине. Это допущение позволяет решать уравнение реакции среды в стационарной постановке. В связи с тем, что трещина плоская, то, можно считать, что она делит пространство на две половины, реакцию каждой из которых можно считать отдельно. Если считать, что упругие модули одинаковы во всем пространстве, то упругую реакцию стенки трещины допустимо оценить по следующей формуле:
где Р(х, у) - поле давления, Е - модуль Юнга, V - коэффициент Пуассона, О - область интегрирования, а(х,у) - поле раскрытия. Формула (1) позволят производить расчет реакции среды не только внутри трещины, но и за ее пределами. Это будет использовано при реализации критерия продвижения фронта трещины.
В уравнениях, приведенных далее, под Р будет пониматься суммарное давление от реакции среды из формулы (1) и минимальных сжимающих напряжений а:
Ирвина.
(1)
Р (х,у) = Ре(х,у) + а(х,у).
(2)
2.3. Модель течения жидкости и продвижения фронта трещины
Предполагается, что длина трещины много больше её ширины, и соответственно допустимо пренебречь течением жидкости в направлении, перпендикулярном стенкам трещины. В таком случае течение жидкости в трещине ГРП приводится в приближении смазочного слоя. Ниже показана зависимость потока от градиента давления в плоской трещине для неньютоновской жидкости:
ч К ™ ч 1 1 , ^
-(-УР + рд)п а п +2, (3)
2 п+!(I + 2)
где ц - поток, К - коэффициент густоты потока, п - степень жидкости, д - ускорение свободного падения. Вместе с уравнением переноса формулы (1) и (3) образуют замкнутую систему уравнений, описывающую течение жидкости и реакцию среды на трещину:
1
-(-УР + рд) п а п +2
2 п+1 1+2
п
(4)
р (х,у)=- щЬ?) //п ^ )2 ¿"'¿у.
Граничным условием уравнения (4) является то, что поток через границу трещины нулевой.
В качестве условия продвижения фронта трещины используется силовой критерий Ирвина [26]. Он формулируется следующим образом: трещина начинает распространяться тогда, когда коэффициент интенсивности напряжений достигает критической величины.
Кт = Кс.
Величину Кс ещё называют трещиностойкостью. Для определения коэффициента интенсивности напряжений К/ воспользуемся тем фактом, что асимптотически для плоской трещины напряжение за фронтом трещины выражается через коэффициент интенсивности напряжений следующей формулой [26]:
К1 ^ у2тгг
где г - расстояние от фронта трещины. Соответственно, для определения К1 за фронтом трещины необходимо произвести аппроксимацию найденного за фронтом трещины напряжения к формуле (5).
3. Численная реализация
Явные схемы требуют значительного измельчения шага по времени. В случае решения системы уравнения (4) явным методом длительность шага по времени может составлять сотые и тысячные доли секунды. Такие результаты были получены в численном эксперименте при попытке реализовать явную схему. Из малых шагов по времени следует пропорционально большее их количество. При этом в ходе каждого шага потребуется производить взятие интеграла (1), который, как будет показано далее, имеет вычислительную сложность порядка О (И4), гд е N - характерное количество ячеек в расчетной сетке вдоль одной из осей в плоскости трещины.
Еще одним вариантом численного решения является применение неявного метода, который позволяет повысить устойчивость решения и увеличить шаг по времени. Однако специфика построения неявной схемы приводит к необходимости решения системы нелинейных уравнений методом Ньютона, что в случае построения полностью неявной схемы ведет к решению СЛАУ с не разреженной матрицей. Для уравнения (4) матрица будет
полностью заполненной, так как давление в каждой точке трещины зависит от раскрытия по всей поверхности трещины. В результате, несмотря на сократившееся количество итераций, сложность самих итераций возрастет многократно.
В данной работе предлагается расщепить систему уравнений на каждой итерации по времени.
Система уравнений (4) состоит из уравнения реакции среды, уравнения для потока (3) и закона сохранения объёма.
Расчет каждого шага по времени предлагается организовать в следующей последовательности:
1) Производится расчет давления в трещине Р по формуле (2).
2) Решается система уравнений, состоящая из уравнения для потока (3) и закона сохранения объёма. Причем давление в уравнении потока учитывается как Р + йР, где йР -оценка изменения давления на следующем шаге по времени, выражение для которой будет приведено далее.
В разделении системы уравнений (4) на два выше обозначенных этапа и заключается расщепление по физическим процессам. Схематичное описание операций, производимых за один шаг по времени для явного, неявного методов и для варианта с расщеплением приведено на рис. 1, 2 и 3.
Нахождение реакции
Проверка критерия протрескивания и рост трещины
Решение уравнения гидродинамики вместе с полным уравнением
Проверка критерия протрескивания и рост трещины
Нахождение реакции среды
Решение уравнения идродинамики вместе упрощенной
Проверка критерия протрескивания и рос! трещины
Риг,. 1. Ятша.я схема.
Рис. 2. Неявная схема.
Рис. 3. Схема с расщеплением
3.1. Разностная схема
Для решения задачи в конечных разностях вводится двумерная регулярная сетка с квадратными ячейками. Шаг сетки по вертикальной и горизонтальной оси идентичен и обозначается как К. На сетке определяется дискретное поле раскрытия трещины а^. Пусть известно состояние (иоле раскрытия) трещины на шаге к расчета. Обозначим раскрытие в ячейке г,] расчетной сетки как а^. Тогда давление в каждой ячейке - результат численного интегрирования уравнения (1). Фактически полем раскрытия трещины определяется состояние трещины в каждый момент времени. При дискретизации решения уравнения (4) по времени за значение, которое необходимо найти на следующем шаге по времени, при-
Л+1
ячеек, а значения потоков для граней ячеек.
нимается поле раскрытия а^1. Далее значения раскрытия и давления даны для центров
а
к+1 пк пк+1 _ пк+1 к+1 + пк+1
'г, 1 - аг,3 + У-г+1/2,з %-1/2,з + %,^+1/2 + %,р-1/2 = ? к + ? к (
+ к + К ивоигсе,г,3 + и1еако//,г,3,
где а);, - раскрытие в ячейке с индексом 1,3 на шаге по времени к, 1 - раскрытие в ячейке с индексом г, з на шаге к + 1, т - шаг по времени, Ъ - шаг по пространству, Укяоигсе ^ ^
ь. ■ .
и р1еакоЦ г з ~ площадные источники и площадные утечки в ячеике с индексом г, 3 на шаге к соответственно. Обозначение с^Хх/^з означает поток между ячейками 1,3 и г + 1,] на шаге к + 1. Соответственно д^1 5 ^++4/2 и обозначают потоки между ячейками
1,3 иг - 1,_7, г,з +1, г,з — 1 на шаге к + 1. Значени е г/оигсе^^ отлично от нуля в ячейке инициализации трещины во время работы насоса и равно где Q - расход жидкости ГРП. Разностная схема для уравнения баланса жидкости из уравнения (4) приведена в формуле (6). Линеаризация потока на примере ^++1/2 ^ будет приведена далее. Для (¡1+\/2 р Я4++1/2 и
(¡1++-\/2 способ линеаризации аналогичен.
3.2. Численное интегрирование реакции среды
Предлагается метод численного взятия интеграла (1). Метод основан па том, что часть интеграла предлагается взять численно, а часть - аналитически. Похожее решение предлагалось в статье [27], по в пей рассматривался другой вид интеграла. При взятии интеграла (1) основной проблемой является особая точка х, у, подынтегральное выражение в которой стремится к бесконечности. Для решения этой проблемы предлагается отделить сингулярную составляющую от А а^представив интеграл в следующем виде:
Р,, —__
Р J 8тг(1 - 1/2)
А а ,з
г, 3
1
cell
у/ (хг' -x)2 + (Уз' -У ')2
dx dy ,
(7)
где сеIи, ^ обозначает область ячейки с координатами %, 3 на расчетной области, а А а*1 ^ -
,
тегральпое выражение:
1
cell
у/ (хг> -x )2 + (Уз' -У)2
dx dy .
(8)
Рис. 4. Обозначение углов ячейки
Интеграл в формуле (8) не зависит от раскрытия трещины и имеет аналитическую первообразную. Значение же интеграла зависит только от положения границ ячейки се1относительно координат ху,
Ii>,j' —
Уг2
1
d
У и
dx —
у/(хг! -Х )2 + (yf -у')2
— A(xi', yj', хг2 , У32 ) — A(xi', Uj', xii, У32 ) —
- А(хг,, у3' ,x%2, у л) + A(x,j, yr ,x%i, уп), (9)
X
i
где Хг1: Хг2, у^ и у^2 ~ координаты углов ячейки се^^- по осям х и у, как показано на рис. 4. Функция А(х, у, хс, ус) определяется так:
А(хс, ус, х, у) = -(у - ус)1 п(л/х2 - 2хсх + у2 - 2усу + х2 + у2 + хс - х)- (хс - х)1 п(л/х2 - 2хсх + у2 - 2усу + х2 + у2 - ус + у). (10)
Так как интеграл (9) зависит только от положения точки (х^ ,уу) относительно ячейки, то для случая регулярной прямоугольной сетки его можно рассчитать один раз, при хг' = 0 и = 0. В дальнейшем результат вычисления интеграла при х# = 0 и х# = 0 можно использовать для всех ячеек расчетной сетки. Для этого интегрирование необходимо представить в виде получения дискретного двумерного поля весов Ш, каждый элемент ,
= А(0,0,х^,уь) - А(0,0,х^1 ,ун) - А(0,0,х^,уу) + А(0,0,х^,у^). (11)
Дискретное поле лапласианов раскрытия, содержащую лапласианы раскрытия всей расчетной области, обозначим как А а. Дискретное поле А а вычисляется путем применения дискретного оператора Лапласа к дискретному полю раскрытия трещины:
&1-1,]-1 + 0,1-1^+1+0,1+1^-1 + 0.1+1,,]+1 4К2
аг-1,] + а1+1,] + аг,]-1 + а1,]+1 3аг,]
А а. . = ™г-1,]-1 | ^г-1,]+1 | ^г+1,]-1 I ^г+1,]+1 +
2К2 К2
(12)
Тогда двумерное поле реакции среды на трещину (перпендикулярная к трещине составляющая напряжения у поверхности трещины) вычисляется как дискретная свертка двух полей:
Е
р = - ¡ЙТ-гз) Ш * А°- (13)
Для расчета свертки в формуле (13) требуется 0(Ы4) операций, где N - характерный линейный размер расчетной сетки. Поэтому, данное уравнение ограниченно применимо для случая больших расчетных сеток, что важно при научных расчетах. Однако при инженерных расчетах такая асимптотика вычислительной точности расчета давления никак не оказывает решающего влияния - решение гидродинамической задачи занимает больше времени.
3.3. Линеаризация потока
Перепишем формулу (3) в новых обозначениях, где иод записью (—У(Р) + рд)« понимается (—V Р + рд)(| — УР + рд|)«-1. Для краткости обозначим далее (-УР + рд) как УР, а | — V Р + рд| как IVРВ новых обозначениях формула (3) обретает вид:
1
ТУ'--
^ " УР|УР| ъ-1а^+2. (14)
2 *+41 +2)
Пусть давление Р^ вычислено. Тогда, по формуле ниже вычислен численный градиент давления между ячейками V Р^+\/2] с поправкой на гравитацию:
р к _ р к Х7Рк = ^ + пп-
У ^+1/2,] = к +РП г,
где д1 - ускорение свободного падения вдоль индексов г расчетной сетки. Тогда аппроксимация потока из формулы (14) между двумя соседними ячейками г^ и г + Т, ] вычисляется
следующим образом:
„к = К " х/рк ({Х7р\к \ 1/п—1( + Ог+1,з , 1 +2 п ,
1/2,з = 1 +2( 1 , У^г+1/2,] 1 г+1/2,^ I 2 ) ,
2™ 1га
где иод IV Р!1+\/2з подразумевается численная оценка модуля градиента давления и гравитационного градиента между ячейками 1,3 и 1 + Т,]. Потоки между ячейками г,] и г — Т,], г, ] + Т Ц] — Т записываются по аналогии с (15). Был найден способ аппроксимации потока на следующем шаге по времени, который позволяет свести решение уравнения (4) неявным методом лишь к одному решению СЛАУ с разреженной матрицей, без необходимости итерационного решения нелинейных уравнений.
Для оценки потока на следующем шаге воспользуемся рядом Тейлора первого порядка для потока (3):
° ° ° ° да да
^ + М) = + = + дЦа + ■^pdVP. (16)
Частная производная потока по раскрытию ^ равна
д° = (Т + 2) 1 " (^Р + рд) 1а *+1.
да уп 2£ +2(га + 2)( Рд)
Принимается следующая численная аппроксимация частных производных ^ и ¿^тр-, которая для случая потока между ячейками г,] и г + Т,] равна:
(t W = (1 + 2)(|VP|Wj)i-(^^)i+ (17)
= - i KVP\+^,)i"'()i+2. dB)
9VPi+1/2,1 2i+2(i + 2)"" "+'''"" 4 2
Из формул (16), (17), (18) следует линеаризация потока (¿¡++1/2 у
К " wnfc ^Dlfc \i-1taXj +ai+l,j\i+2
q++1/2,j - 1 +2i i , о\ V Pi+1/2,j(1 V P1 i+l/2,j)n ( 2 ' '3) " +
2i l n + 2J 2
+(t+i/^^+i^+dvp++i/2*, (ig)
i+1/2,j
где da++i/2j = 0+++1^/2- — a++1/2 ^ а значение dPj++1/2- оценивается численно. Для этого дифференциал градиента давления представляется через дифференциал давления:
jwnfc dPi+i, j — dPi,j dV Pi+i/2, j = --" •
d P
трещины, это следует из формулы (1). Причем, вследствие линейности формулы (1) по раскрытию, зависимость будет аналогичной:
P(x,y, ti) — P(x,y, fc) = — j-^fl ^-y'M-^^vf dxdy, (20)
8^(1 — V2) JJq ^(x — X )2 + (y — y')2
где P (x, y, ti) — P (x, y, io) - изменение давления, a Aa(x ,y , ti) — Aa(x ,y , to) - изменение
P( x, )
ла (20) по всей области трещины. Однако для целей численной аппроксимации потока (19)
достаточно взять интеграл в некоторой области вокруг точки (х, у). Формула для численной оценки приведена ниже:
<Ргк, =--. Е 0, V АБау ?> [[ . 1 йхйу . (21)
— *2) 1,3 ^(хг — х' )2 + (у3—у ')2 У Х '
Формула (21) эквивалентна численному вычислению реакции упругой породы на трещину (7), за исключением того, что вместо раскрытия а берется изменение раскрытия Ба между соседними шагами по времени, а суммирование предлагается проводить в области, близкой к ячейке ( г, ]).
3.4. Программная реализация
Программный код модели реализован на языке программирования С++, для прото-типирования и отладки отдельных частей модели использовался язык программирования Python. Реализованные на С++ программные коды распараллелены с использованием ОрепМР. Программная реализация допускает возможность компиляции кода как на операционные системы Windows, так и Linux. API расчетного модуля реализовано таким образом, чтобы было возможно обращаться к нему из языка Python. Это облегчает подключение расчетного модуля к иным ПО и GUI.
4. Верификация реализованной модели
4.1. Верификация критерия продвижения фронта трещины
Критерий продвижения фронта трещины, реализованный по силовому критерию Ирвина, следует проверять вместе со всеми остальными уравнениями на таком кейсе, для которого точно известно аналитическое решение задачи гидроразрыва. Классические аналитические модели для этого не подходят, так как в них никак не учитывается влияние трещиностойкости. Однако была найдена постановка задачи, для которой можно найти простое аналитическое решение, учитывающее влияние трещиностойкости.
Постановка задачи следующая: пусть в однородный упругий пласт с нулевым коэффициентом утечек и с ненулевой трещиностойкостью происходит закачка текучей жидкости. В таком случае трещина окажется круглой. Какое-то время после окончания закачки она будет растекаться. Это будет происходить до тех пор, пока коэффициент интенсивности напряжений за фронтом трещины больше трещиностойкости. Когда коэффициент интенсивности напряжений станет равен трещиностойкости породы, рост трещины прекратится на определенном радиусе Rk, который можно оценить аналитически. При этом закачка жидкости должна происходить с достаточно высоким расходом, чтобы радиус трещины к моменту окончания закачки был хотя в несколько раз меньше Rk-
Опишем предложенную постановку задачи формально. Пусть в породу с модулем Юнга ^коэффициентом Пуассона v и трещиностойкостью Kj max была с высоким расходом закачена жидкость объема V. К тому моменту, как трещина в процессе растекания достигнет своего передельного радиуса, будем считать трещину эллиптической. Давление в такой трещине будет равно [28]:
рк:,,„. ^
ак^Е 8(1 — V2 )Кк
Раскрытие трещины ак вычисляется исходя из формулы для объема эллипсоида:
аК = . (23)
Интенсивность напряжений за фронтом эллиптической трещины радиуса Кк, давление в которой равно Рк, вычисляется по формуле [29]:
К = 2РкУКк. (24)
Воспользовавшись формулами (22), (23) и (25), нетрудно показать, что радиус трещины на момент окончания растекания, когда К\ = К\ тах, будет равен
Кк = (8 К^Е у2) 2 . (25)
Т а б .л и ц а 1
Основные параметры среды и закачки
Название параметра Обозначение Значение
Модуль Юнга среды Е 10 ГПа
Коэффициент Пуассона V 0
Объем закачки V 20 м3
Трещиностойкость К! тах 5е5, 1е6 и 2.5е6 Пау/м
Формула (25) позволяет проводить верификацию работы критерия продвижения фронта трещины путем сравнения Кк с результатом численного расчета.
Для сравнения использовались параметры, приведенные в табл. 1.
В численных расчетах закачка 20 м3 жидкости производилась за 500 секунд. На установление состояния трещины после окончания закачки отводилось 160 000 секунд. Вязкость закачиваемой жидкости равна 1 си. Размер расчетных ячеек 5 метров.
Рис. 5. К/ тах = 0,5 ГПал/м
Рис. 6. К/ тах = 1 ГПал/м
Рис. 7. К/ тах = 2,5 ГПа^м
На рис. 5, 6 и 7 приведено сравнение аналитического и численного решений. Как видно, наблюдается совпадение численного и аналитического радиуса трещины с точностью до размера ячейки, что можно интерпретировать как то, что удалось реализовать простой, но при этом достаточно точный критерий продвижения фронта трещины.
4.2. Сравнение с 1Ь8А
Для сравнения с 11.5А была выбрана статья за авторством Е. V. БоШяоу и А. Р. Рексе [27], в которой есть несколько кейсов, описывающих работу их Р1аиагЗО модели с П.5А критерием продвижения фронта трещины. Соответственно, сравнение результатов расчета реализованной в данной работе модели с кейсами из данной статьи позволит не только оценить адекватность реализации модели течения в трещине и расчета реакции среды, но и оценить обоснованность выбранного критерия продвижения фронта трещины.
Будут рассмотрены кейсы в трехслойной среде, схема которой с обозначениями приведена на рис. 8. Упругие модули среды и параметры закачки приведены в табл. 2.
Т а б л и ц а 2
Параметры среды и закачки
Название параметра Обозначение Значение
Модуль Юнга среды Е 9.5 ГПа
Коэффициент Пуассона V 0.2
Вязкость жидкости ГРП д 0.1 Па с
Расход жидкости Я 0.01 ^
Коэффициент утечек Сь 0.2605, 0.825 и 2.605 ^е-5 ' л/с
Время закачки жидкости Т 3600 с
Рис. 8. Трехслойная среда для сравнения с
Рис. 9. Сь = 2.605е-5 ^
Рис. 10. Сь = 0.825е-5 ^
Рис. 11. Сь = 0.2605е-5 ^
Для симметричного случая напряжения приняты такими: о\ = а3 = 7.75 МПа, ч2 = 7 МПа. Результаты расчетов на описанной в этой работе модели Р1апагЗБ при различных утечках на момент времени 3600 секунд от начала закачки приведены на рис. 9, 10 и 11. На рисунках показано поле раскрытия с размерностью в метрах. Геометрия трегцнн с точностью до размера одной-двух ячеек, составляющего 1.8 м, совпадает с результатами, приведенными в статье [27].
Для асимметричного случая с высоким контрастом между верхними слоями напряжения равны: о\ = 7.25 МПа, ч2 = 7 МПа, ч3 = 6.5 МПа. Результаты расчетов для этого случая приведены на на рис. 12, 13 и 14. Геометрия совпадает с той, что приведена на соответствующих рисунках в статье [27] с точностью до размера ячейки.
Совпадение геометрий трещины при моделировании на предлагаемой в этой работе модели Р1апагЗБ с результатами 1Ь8А Р1апагЗБ позволяет говорить о том, что в предложен-
ной модели корректно реализована гидродинамическая и гсомсханичсская часть, подтверждается корректность выбранного критерия продвижения трещины.
Также
Рис. 12. CL = 2.605е-5 ^
Рис. 13. CL = 0.825е-5 ^
Рис. 14. CL = 0.2605е-5 ^
5. Программная реализация
Программный код модели реализован на языке программирования С++ для прото-тинирования и отладки отдельных частей модели использовался язык программирования Python. Реализованные на С++ программные коды распараллелены с использованием ОренМР. Программная реализация допускает возможность компиляции кода как на операционные системы Windows, так и Linux. API расчетного модуля реализовано таким образом, чтобы было возможно обращаться к нему из языка Python. Это облегчает подключение расчетного модуля к иным ПО и GUI. В качестве бенчмарка скорости работы программной реализации и демонстрации эффективности численного алгоритма был выбран кейс, показанный на рис. 14. Кейс представляет собой закачку 36 м3 жидкости при размере расчетной ячейки 1.8 метров. Время расчета этого кейса на ПК с процессором Intel Core i5 7600k составляет 62 секунды в однопоточном режиме и 37 секунд при распараллеливании расчета на 4 потока.
6. Дальнейшее развитие модели
Модель, описанная в уравнении (4), может быть расширена для расчета транспорта пропанта и течения нескольких жидкостей в трещине.
6.1. Транспорт пропанта
Моделирование транспорта пропанта производится в предположениях теории суспен-
{—щ■ + div(qp) = ^ ^source prop, (26)
+ div(qf) = ^ Vsource fluid.
При этом пропант моделируется как отдельная фаза, и для него записывается отдельное уравнение баланса масс, как показано в уравнении (26), где С - объемная концентрация пропанта, qp - поток пропанта, qf поток жидкости.
К- й , __ ^ i i ,2, С rft
г, П г, 1 I
qs = qp + qf = 1 " ^ (-VP + pg) £ a * +2(1 - —) " . (27)
2 n+!( 1 + 2) ^max
В формуле (27) Стах - максимальная объемная концентрация пропанта, ш принимается равной 1,89. Потоки qp и д/ выражаются через средний поток смеси, выражение для которой приведено в формуле (27). Формула (27) аналогична формуле для потока жидкости с тем различием, что коэффициент густоты потока К8 для смеси пропант-жидкость предполагается равным:
/ С ч 5п
К = К (1 — —) 2 .
Стах
Для моделирования эффектов оседания пропанта вводится понятие отрыва скорости пропанта от жидкости:
Vp = У, + (1 - C)AV.,
где Vp = qp/а a Vs = qs/а. Для скорости отрыва пропанта предлагается следующее выра-
р
жение:
А V =
-gгаdP + Ppld2 ^ _
C
18^
Cm
d
Рис. 15. Моделирование оседания пропанта в трещине
На рис. 15 качественно приведено распределение пропанта для кейса, похожего на тот, что был приведен в статье за авторством Cohen Charles-Edouard et al. [25].
m
)
6.2. Течение нескольких жидкостей
В процессе проведения ГРП в трещину закачивается несколько типов жидкостей. Обычно это линейный гель, которым заполнена насосно-компрессорная труба (НКТ), и сшитый гель, которым производится накачка подушки и транспорт пропанта в трещину.
Течение нескольких жидкостей описывается системой уравнений баланса массы, состоящей из уравнений типа (29) для каждой жидкости.
даС^г ^^
^ = , (28)
где ^ таток г-й жидкости, а С/^ - её объемная концентрация. Основной проблемой при моделировании течения нескольких жидкостей в трещине ГРП является вопрос оценки потока смеси жидкости qf и потока отдельных жидкостей. Предложим следующую модель, построенную на следующих предположениях:
1) Жидкости не испытывают взаимного трения.
2) Каждая из жидкостей занимает все пространство между стенками трещины.
Тогда для каждой отдельной жидкости остается применимым приближение смазочного
Ч/г = Сг
__
Ки П'1 ^ -Х- +2 -(—V Р + рд) пПа ПП .
(29)
2 (^ + 2)
При этом поток смеси - сумма потов всех жидкостей:
(30)
Приведенная модель течения нескольких жидкостей позволяет с высокой точностью моделировать случай поршневого вытеснения, когда по ходу закачки вязкость увеличивается. При нестабильном режиме вытеснения, когда текучая жидкость вытесняет вязкую, вязкие пальцы не разрешаются, однако возникает эффект фильтрации текущей жидкости через вязкую - что качественно верно.
Рис. 16. Моделирование течения нескольких жидкостей в трещине
1
На рис. 16 показан пример расчета с течением нескольких жидкостей, закачиваемых в следующем порядке (в скобках указаны параметры степенной жидкости): (К = 0.05, п = 1.0), (К = 0.5, п = 0.9), (К = 5, п = 0.55). Так как вязкость жидкостей возрастает, то наблюдается поршневое вытеснение и возрастание забойного давления с каждой новой стадией.
7. Заключение
Разработан новый подход к численному моделированию ГРП в постановке Р1апагЗБ, заключающийся в расщеплении по физическим процессам. Реализованы вычислительные алгоритмы, решающие уравнения Р1апагЗБ согласно разработанному подходу путем решения уравнения явно-неявным методом. Разработан численный алгоритм, реализующий проверку критерия продвижения трещины по силовому критерию Ирвина. Верификация подтверждает корректность программной реализации и численной схемы.
Литература
1. Шевченко Д.В., Шевченко В.П. Выбор и оптимизация структуры построения автономных сейсмических средств обнаружения рубежного типа // Материалы VIII всероссийской научно-технической конференции «Современные охранные технологии и средства обеспечения комплексной безопасности объектов». 2010. С. 128-133.
2. Diallo M.S., Kulesh М., Holschneider М., Sherbaum F., Adler F. Characterization of polarization attributes of seismic waves using continuous wavelet transforms // Geophysics. 2006. V. 71, N 3. P. 67-77.
3. Лайонс P. Цифровая обработка сигналов. Москва : Бином, 2006. С. 361-369.
4. Sneddon I.N. The distribution of stress in the neighbourhood of a crack in an elastic solid // Proceedings of the Royal Society of London A: Mathematical Physical and Engineering Sciences. 1946. V. 187, N 1009. P. 229-260.
5. Khristianovich S.A., Zheltov Y.P. Formation of vertical fractures by means of highly viscous liquid // Proceedings of Fourth World Pet. Congress, Rome. 1955. V. 2, P. 579-586.
6. Geertsma J., Klerk F.A. Rapid method of predicting width and extent of hydraulic induced fractures // J. Pet Tech. 1969. V. 12. P. 1571-1581.
7. Zazovskii A.F. Growth of coin-shaped crack produced by hydraulic fracture in an impermeable rock 11 Mech. Solids. 1979. V. 2, N 14. P. 104-111.
8. Perkins Т.К., Kern L.R. Widths of Hydraulic Fractures //J. Pet. Technol. 1961. V. 13. P. 937-949.
9. Nordgren R.P. Propagation of a vertical hydraulic fracture // SPE 3009-PA. 1972.
10. Mack M.G., Warpinski N.R. Mechanics of hydraulic fracturing, Reservoir stimulation, 3rd ed. Chichester : Wiley, 2000. P. 856.
11. Cleary M.P. Analysis of mechanisms and procedures for producing favourable shapes of hydraulic fractures 11 SPE 9260-MS. 1980.
12. Meyer B.R. Design formulae for 2-D and 3-D vertical hydraulic fractures: model comparison and parametric studies // SPE 15240-MS. 1986.
13. Smith M.B., Klein H.A. Practical Applications of Coupling Fully Numerical 2-D Transport Calculation With a PC-Based Fracture Geometry Simulator // SPE 30505. 1995.
14. Шель E.B., Падерин Г.В. Модель для экспресс-оценок дизайна ГРП с использованием приближенного аналитического решения / / PRO НЕФТЬ. Профессионально о нефти. 2017. № 4(6). С. 40-43.
15. Barree R.D. A practical numerical simulator for three-dimensional fracture propagation in heterogeneous media 11 SPE 1227 3-MS. 1983.
16. Smith M.B., Klein H.A. Practical Applications of Coupling Fully Numerical 2-D Transport Calculation With a PC-Based Fracture Geometry Simulator // SPE 30505. 1995.
17. Adachi J., Siebrits E., Peirce A., Desroches J. Computer simulation of hydraulic fractures // International Journal of Rock Mechanics k, Mining Sciences. 2007. P. 739-757.
18. Аксаков А.В., Борщук О. С., Желтова И. С., Дедурин А.В., Калуджер 3., Пестри-ков А.В., Торопов К.В. Корпоративный симулятор гидроразрыва пласта: от математической модели к программной реализации // Нефтяное хозяйство. 2016. № 11. С. 35-40.
19. Nolte K.G., Economides M.J. Reservoir Stimulation Chichester, NY : John Wiley k, Sons, 2000.
20. Clifton R.J., Abou-Sayed A.S. On the computation of the three-dimensional geometry of hydraulic fractures // SPE 7943-MS. 1979.
21. Carter B.J., Desroches J., Ingraffea A.R., Wawrzynek P.A. Simulating fully 3D hydraulic fracturing // Modeling in geomechanics. 2000. V. 200. P. 525-557.
22. Zeng Q., Yao J., Shao, J. Effect of plastic deformation on hydraulic fracturing with extended element method // Acta Geotech. 2019. N 14. P. 2083-2101.
23. Baykin A.N., Golovin S. V. Modelling of hydraulic fracture propagation in inhomogeneous poroelastic medium 11 J. Phvs.: Conf. Ser. 2016. V. 722, N 012003.
24. Clifton, R.J., Wang J.J. Multiple Fluids, Proppant Transport, and Thermal Effects in Three-Dimensional Simulation of Hydraulic Fracturing // 63rd Annual Technical Conference and Exhibition in Houston, TX, 1988.
25. Cohen C.-E., Kresse O., Weng X. Stacked Height Model to Improve Fracture Height Growth Prediction, and Simulate Interactions With Multi-Layer DFNs and Ledges at Weak Zone Interfaces // SPE-184876-MS. 2017.
26. Матвиенко Ю.Г. Модели и критерии механики разрушения. Москва : Физматлит, 2006. С. 328.
27. Dontsov Е. V., Peirce А.P. A multiscale Implicit Level Set Algorithm (ILSA) to model hydraulic fracture propagation incorporating combined viscous, toughness, and leak-off asvmptotics // Computer Methods in Applied Mechanics and Engineering. 2017. V. 313. P." 53-84.
28. Martin P.A. On wrinkled penny-shaped cracks // Journal of the Mechanics and Physics of Solids. 2001. V. 49. P. 1481-1495.
29. Мураками Ю. Справочник по коэффициентам интенсивности напряжений: В 2-х томах. Т. 2. Москва : Мир, 1990. 1016 с.
References
1. Shevchenko D. V., Shevchenko V.P. Selection and optimization of constructing autonomous seismic detection struction a landmark type. Proceedings of the VIII Russian scientific conference «Modern security technology and means of complex security objectives». 2010. P. 128-133. (in Russian).
2. Diallo M.S., Kulesh M., Holschneider M., Sherbaum F., Adler F. Characterization of polarization attributes of seismic waves using continuous wavelet transforms. Geophysics. 2006. V. 71, N 3. P. 67-77.
3. Lyons R. Digital signal processing. Moscow : Binom, 2006. P. 361-369. (in Russian).
4. Sneddon I.N. The distribution of stress in the neighbourhood of a crack in an elastic solid. Proceedings of the Royal Society of London A: Mathematical Physical and Engineering Sciences. 1946. V. 187, N 1009. P. 229-260.
5. Khristianovich S.A., Zheltov Y.P. Formation of vertical fractures by means of highly viscous liquid. Proceedings of Fourth World Pet. Congress, Rome. 1955. V. 2, P. 579-586.
6. Geertsma J., Klerk F.A. Rapid method of predicting width and extent of hydraulic induced fractures. J. Pet Tech. 1969. V. 12. P. 1571-1581.
7. Zazovskii A.F. Growth of coin-shaped crack produced by hydraulic fracture in an impermeable rock. Mech. Solids. 1979. V. 2, N 14. P. 104-111.
8. Perkins Т.К., Kern L.R. Widths of Hydraulic Fractures. J. Pet. Technol. 1961. V. 13. P. 937-949.
9. Nordgren R.P. Propagation of a vertical hydraulic fracture. SPE 3009-PA. 1972.
10. Mack M.G., Warpinski N.R. Mechanics of hydraulic fracturing, Reservoir stimulation, 3rd ed. Chichester : Wiley, 2000. P. 856.
11. Cleary M.P. Analysis of mechanisms and procedures for producing favourable shapes of hydraulic fractures. SPE 9260-MS. 1980.
12. Meyer B.R. Design formulae for 2-D and 3-D vertical hydraulic fractures: model comparison and parametric studies. SPE 15240-MS. 1986.
13. Smith M.B., Klein H.A. Practical Applications of Coupling Fully Numerical 2-D Transport Calculation WTith a PC-Based Fracture Geometry Simulator. SPE 30505. 1995.
14. SheI E. V., Paderin G. V. The model for the express-estimation of the hvdrofracturing design with the approximate analytical solution. PRONEFT'. Professional'no o nefti. 2017. N 4(6). P. 40-43.
15. Barree R.D. A practical numerical simulator for three-dimensional fracture propagation in heterogeneous media. SPE 1227 3-MS. 1983.
16. Smith M.B., Klein H.A. Practical Applications of Coupling Fully Numerical 2-D Transport Calculation With a PC-Based Fracture Geometry Simulator. SPE 30505. 1995.
17. Adachi J., Siebrits E., Peirce A., Desroches J. Computer simulation of hydraulic fractures. International Journal of Rock Mechanics k Mining Sciences. 2007. P. 739-757.
18. Aksakov A. V., Borshuk O.S., Zheltova I.S., Dedurin A. V, Kaludzher Z., Pestrikov K. V., Toropov K. V. Corporate fracturing simulator: from a mathematical model to the software development. Neftvanoe khozvavstvo - Oil Industry. 2016. N 11. P. 35-40. (in Russian).
19. Nolte K.G., Economides M.J. Reservoir Stimulation Chichester, NY : John Wiley k Sons, 2000.
20. Clifton R.J., Abou-Sayed A.S. On the computation of the three-dimensional geometry of hydraulic fractures. SPE 7943-MS. 1979.
21. Carter B.J., Desroches J., Ingraffea A.R., Wawrzynek P.A. Simulating fully 3D hydraulic fracturing. Modeling in geomechanics. 2000. V. 200. P. 525-557.
22. Zeng Q., Yao J., Shao, J. Effect of plastic deformation on hydraulic fracturing with extended element method. Acta Geotech. 2019. N 14. P. 2083-2101.
23. Baykin A.N., Golovin S. V. Modelling of hydraulic fracture propagation in inhomogeneous poroelastic medium. J. Phvs.: Conf. Ser. 2016. V. 722, N 012003.
24. Clifton, R.J., Wang J.J. Multiple Fluids, Proppant Transport, and Thermal Effects in Three-Dimensional Simulation of Hydraulic Fracturing. 63rd Annual Technical Conference and Exhibition in Houston, TX, 1988.
25. Cohen C.-E., Kresse O., Weng X. Stacked Height Model to Improve Fracture Height Growth Prediction, and Simulate Interactions With Multi-Layer DFNs and Ledges at Weak Zone Interfaces. SPE-184876-MS. 2017.
26. Matvienko Y. .G. Models and criteria of fracture mechanics. Moscow : Fizmatlit, 2006. P. 328. (in Russian).
27. Dontsov E.V., Peirce A.P. A multiscale Implicit Level Set Algorithm (ILSA) to model hydraulic fracture propagation incorporating combined viscous, toughness, and leak-off asvmptotics. Computer Methods in Applied Mechanics and Engineering. 2017. V. 313. P.'53-84.
28. Martin P.A. On wrinkled penny-shaped cracks. Journal of the Mechanics and Physics of Solids. 2001. V. 49. P. 1481-1495.
29. Murakami Y. Handbook on stress intensity factors: In 2 volumes. V. 2. Moscow : Mir, 1990. 1016 p. (in Russian).
Поступим в редакцию 23.08.2021
