Разработка термодинамически согласованных математических моделей и методов математического моделирования для анализа тепловых методов увеличения нефтеотдачи

Во введении диссертации описаны основные эффекты, сопровождающие
процесс термического воздействия на пласт, и особенности их моделирования.
Показана актуальность задачи комплексного моделирования пласта с учетом
деформации породы, фильтрации флюида, учета неизотермических эффектов
и разрушения породы. Описаны существующие подходы к моделированию, их
преимущества и недостатки. Сформулированы цели и задачи работы, а также
полученные результаты. Указаны сведения об апробации работы и публикациях
по теме исследования.
В первой главеописана физико-математическая модель разрушения в по-
роупругой среде, являющаяся обобщением пороупругой модели Био. Рассмат-
ривается элементарный объем пространства Ω, в котором содержатся две фазы
— пористый деформируемый скелет и однофазный слабосжимаемый флюид,
движущиеся со скоростями и соответственно, и одновременно занимаю-
щие один и тот же объем.
Обозначим величины, отнесенные к твердой и жидкой фазе, индексами « »
и « » соответственно. Базовая система уравнений модели состоит из фунда-
ментальных законов сохранения массы, импульса и энергии:

+ div( ) + div( ) = 0,
(︂)︂(︂)︂

div ( ) + −+ −= 0,

[︃ (︂(︃)︃]︃
2 2
)︂

++ ++(1)
22
[︃ (︃)︃ ]︃

[︂ (︂)︂ ]︂

+ div + + div + =
= div [ ( − ) + ] + ( + ) − div( + ),

где = , , — масса фазы в элементарном объеме среды, — плотность
флюида, — скорость фильтрации, — тензор полных напряжений, —
вектор внешних сил, — удельная внутренняя энергия фазы , — вектор
потока тепла фазы , оператор (·)/ = (·)/ + grad(·) — материальная
производная для фазы .
Для замыкания системы законов сохранения должны быть заданы опреде-
ляющие соотношения, описывающие ограничения на поведение материала под
влиянием внешнего воздействия. Данные определяющие соотношения для пол-
ной модели должны обеспечивать выполнение условий термодинамической со-
гласованности, то есть выполнение второго начала термодинамики для любой
истории состояний (процесса).
Для получения конкретного вида определяющих соотношений была исполь-
зована процедура Колмана-Нолла. Основная идея данного подхода заключается
в том, что при наличии функциональных связей между параметрами, описыва-
ющими термодинамический процесс, вид определяющих соотношений должен
быть таким, чтобы второй закон термодинамики выполнялся для любой после-
довательности состояний.
Считается, что в результате изменения напряженно-деформированного со-
стояния пласта в нем может возникнуть зона диффузного разрушения, под
которой понимается образование в твердой фазе зон микротрещиноватости. Су-
щественным является то, что характерный размер трещин существенно мень-
ше, чем размер представительного объема среды. По этой причине разрушение
описывается определенной в пространстве величиной (в простейшем случае —
скалярной), значения которой имеют смысл «степени поврежденности» среды.
Для построения определяющих соотношений рассматривается энтропийное
неравенство. В соответствии с классическими представлениями о механизмах
развития трещин, для образования единицы площади её поверхности необхо-
димо затратить определенную энергию. При этом, в связи с тем, что процесс
образования трещин является необратимым, изменение энтропии системы не
меньше, чем количество тепла, полученное системой, а также работа сил, отве-
чающих за разрушение среды.
Таким образом, второе начало термодинамики без учета внешних источни-
ков тепла постулируется в виде:

(︁ )︁ 1
+ + div− :⩾ 0,(2)

где — удельная энтропия фазы , — обобщенная сила, связанная с раз-
рушением материала, — тензор повреждаемости, «:» — оператор свертки
тензора.
Определим свободную энергию Гельмгольца для фазы как = −
. Воспользовавшись выражением для баланса энергии флюида =
− (1/ ) − , где — поровое давление, с учетом определения и нера-
венства (2) получим неравенство, в левой части которого находится сумма дис-
сипаций скелета , флюида и тепловой диссипации :

+ + ⩾ 0,(3)

где
(︃)︃

= − + + 2+

+ [ : grad ( − ) + : grad ( )] − :,

[︂(︂)︂]︂

= ( − ) −+ div( ) , = −grad( ).

Можно показать, что если течение флюида подчиняется закону Дарси, а
поток тепла описывается законом Фурье, то + ⩾ 0 и согласно процеду-
ре Колмана-Нолла можно предположить, что для выполнения неравенства (3)
достаточно, чтобы ⩾ 0.
Введем энергию Гиббса , зависящую от параметров = { , , , } та-
кую, что = − . Тогда выражение для диссипации скелета примет
вид:
(︂)︂(︂)︂

= − − + −

(︂)︂(︂)︂

− +− + ⩾ 0, (4)

где = / — пористость, — тензор малых деформаций.
Частным решением данного неравенства являются соотношения:

= ; = − ( , ) ; = −; = − ,(5)

обеспечивающие выполнение равенства = 0, что удовлетворяет принципу
термодинамической согласованности.
Предположим, что тензор упругих коэффициентов = 2 / 2 являет-
ся линейной функцией от параметра повреждаемости , а = − 2 / 2
— линейной функцией от тензора деформаций , причем должно выполняться
условие / = 3 / 2 = − / .
Будем считать, что энергия Гиббса может быть представлена в виде много-
члена относительно элементов вектора следующего вида:
⃒
1 ⃒⃒
∆ = · ∆ + ∆ · · ∆ +∆ · ∆ · ∆ ,(6)
22 2 ⃒ 0

где ∆ = − 0 для всех , = [ , , , ] , = / , = 2 / 2 .
Дифференцируя энергию Гиббса (6) в соответствии с выражениями (5) по-
лучим определяющие соотношения:

( 0 )
(︂)︂
∆ = ( ) +0
: ∆ : ∆ − ∆ − : ∆ − ( 0 ) : ∆ , (7a)

∆
∆ = : ∆ + − ∆ + : ∆ ,(7b)

∆ = : ∆ − ∆ + 0 ∆ + : ∆ ,(7c)
)︂
(︂
1 ( )
∆ = ( 0 ) +: ∆ : ∆ + ∆ + ∆ + : ∆ .(7d)
2

где — коэффициент Био, 1/ — модуль Био, — тензор термического
расширения, — коэффициент объемного термического расширения, —
теплоемкость скелета, , , , — тензоры коэффициентов, связанных с раз-
рушением, а также приняты следующие обозначения:
=+ , = + , = ,

где — модуль объемного сжатия флюида, — теплоемкость флюида,
— коэффициент теплового расширения флюида.
Внутренняя энергия скелета выражается через энергию Гиббса по формуле:

∆ = ∆ = ∆ + ∆( ) = ∆ + ∆ ( ) + ∆( ).

С учетом соотношений (6) и (7) выражение для внутренней энергии скелета
примет вид:
[︂]︂
1 0000
∆ = + : ∆ + + : ∆ +
[︂ (︂)︂(︂)︂]︂
+ − ∆ − − ∆ ∆ +
22
[︂(︂)︂(︂)︂]︂
1 1
+ − − ∆ + 0 − ∆ ∆ +
2 2
1 0
[︂(︂)︂]︂
− + +0
∆ ∆ − − + ∆ ∆ . (8)
2 2

Аналогично выражение для удельной внутренней энергии флюида:
{︂[︂(︂)︂(︂)︂]︂
1111
∆ = − ∆ − − ∆ ∆ −
22
[︂ (︂)︂(︂)︂]︂}︂
1 1
− − ∆ − − ∆ ∆ . (9)
2 02

Введем следующие дополнительные предположения:
ˆ плотность материала скелета пренебрежимо слабо зависит от давления;
ˆ внешние и инерционные силы пренебрежимо малы;
ˆ кинетическая энергия флюида пренебрежимо мала по сравнению с вели-
чиной внутренней энергии;
ˆ скорость движения скелета пренебрежимо мала по сравнению со скоро-
стью течения флюида;
ˆ параметр повреждаемости является скалярной величиной;
ˆ коэффициенты разрушения = 0, = 0, = 0 (таким образом, явной
зависимостью плотности флюида и энтропии системы от разрушения пре-
небрегаем);
ˆ тензор упругих коэффициентов зависит от повреждаемости по формуле:
˜ − ).
( ) = (1

С учетом данных допущений система уравнений (1) примет следующий вид:

+ div( ) = 0, div = 0,

( + )
+ div ( ) = div (− ) − div( ),(10)

= − grad( ), = − grad( ), = ( ),

которая замыкается определяющими соотношениями вида:

∆ = (1 − ) : ∆ − ∆ − : ∆ − : 0 : ∆ ,
∆
∆ = : ∆ + − ∆ ,

∆ = : ∆ − ∆ + 0 ∆ ,

(11)
(︂)︂
∆ = 0 + ∆ ∆ ,
∆ =∆ − ∆ ,

∆ = − ∆ + 0 ∆ .

Выражение для внутренней энергии скелета имеет вид:

∆ = ∆ + ∆ + ∆ + ∆ ,

где[︂]︂
∆ = 0 + : ∆ + 0 + 0 : ∆ ,
[︂ (︂)︂(︂)︂]︂
∆ = − ∆ − − ∆ ∆ ,
22
[︂(︂)︂(︂)︂]︂(12)
1 1
∆ = − − ∆ + 0 − ∆ ∆ ,
2 2
[︂(︂)︂]︂
∆ = − 0 + 0 + ∆ ∆ ∆ .
Выражение для внутренней энергии флюида имеет вид:
[︂(︂)︂(︂)︂]︂
11
∆ = − ∆ − − ∆ ∆ −
22
[︂(︂)︂(︂)︂]︂
1 1
− − ∆ − 0 − ∆ ∆ . (13)
2 2

В конце первой главы приведены различные варианты интерпретации па-
раметра повреждаемости, встречающиеся в литературе. В наиболее простом
случае, когда разрушение изотропно (то есть происходит одинаково во всех на-
правлениях), параметр повреждаемости является скалярной величиной и име-
ет смысл средней поверхностной плотности пересечений микротрещин с любой
плоскостью внутри тела.
Динамика изменения параметра повреждаемости может быть описана неко-
торым кинетическим уравнением вида / = ( , …), где зависит от и
таких переменных как тензор напряжений, температура, время и так далее. В
случае, если текущая реакция повреждающегося материала зависит от предыс-
тории состояний, то такую кинетику разрушения называют конечной. В ином
случае — мгновенной. В литературе описаны результаты различных экспери-
ментов по определению конкретного вида функции и значений входящих в
неё коэффициентов.
Во второй главеприводится описание численного алгоритма для реше-
ния системы уравнений (10)-(13). Данная система решалась методом конечных
элементов. Введем пробные функции , = { , , } в пространстве 0 ⊂ ,
где 0 = { ∈ : | Ω = 0}. После домножения уравнений (10) на пробную
функцию данная система уравнений примет вид:
∫︁∫︁[︂]︂

div Ω = 0, + div( ) Ω = 0,

ΩΩ
∫︁[︂(︂)︂]︂(14)
1
++ div + div ( ) + div( ) Ω = 0.

Ω

Введем дифференциальный оператор , который является матричным
представлением оператора деформации так, что = , где = [ , , ] .
Предположим, что начальный тензор напряжений 0 постоянен для всего
объема, а начальные деформации 0 равны нулю. Тогда, применяя формулу
Грина к уравнениям (10) и воспользовавшись определяющими соотношениями
для тензора напряжений и массы флюида, можно прийти к слабой постановке
задачи.
Для построения конечномерной задачи введем конечномерные пространства
( )( )
⊂ , 0, ⊂ 0 , причем = span( ), где — базисные функции. Тогда

∑︁
= ( )
, = , , ,
=1

Результирующая система уравнений для конечномерной задачи имеет вид:

− + = ,

+ + + = ,(15)

+ + + + + = ,

где используются следующие обозначения:
∫︁∫︁
)︀ )︀
= =
(︀(︀)︀(︀
= (1 − ) Ω, Ω,
ΩΩ
∫︁∫︁
)︀
(∇ ) / (∇ ) Ω,
(︀
= − (1 − ) : Ω, =
ΩΩ
∫︁∫︁
= −1 Ω, = − Ω,
ΩΩ
∫︁
0 + (1 − ) : + 0 + 0 : Ω,
(︀)︀ (︀ )︀
=
Ω
∫︁
−1 − Ω,
(︀)︀
=
Ω
∫︁
− + ( + ) / 0 Ω,
[︀]︀
=
Ω
∫︁
∇ [1/2 ( ) : ] Ω,
(︀ )︀
=
Ω
∫︁
∇ 0 / + [ / − ] ∆ − [ − / ] ∆ / ∇ Ω,
{︀}︀
=
Ω
∫︁∫︁∫︁
=∇ ∇ Ω, = / , = − / Ω,
Ω
∫︁∫︁ Ω∫︁ Ω

= − / − / − / .
Ω Ω Ω

Дискретизация по времени проводилась по неявной схеме для перемещений,
давления и температуры и явной для параметра повреждаемости по следую-
щему алгоритму:

1. Зачитывание данных модели, инициализация расчетной сетки.
2. Сборка матрицы системы с текущим значением параметра повреждаемо-
сти .
3. Решение системы уравнений относительно перемещений , давлений ,
температуры .
4. Явный расчет нового значения параметра повреждаемости .
5. Повторение пунктов 2-4 для следующего временного шага.

Решение нелинейной системы уравнений проводилось с использованием ме-
тода Ньютона. Для дискретизации уравнения по пространству использовалась
тетраэдральная сетка с квадратичными базисными функциями для переме-
щений и линейными для давления и температуры (элементы Тейлора-Худа).
Данный тип элементов обеспечивает выполнение inf − sup условий (условия
Ладыженской-Бабушки-Бреззи), которые необходимы для устойчивости реше-
ния уравнений пороупругости.
Для обеспечения устойчивости конечномерной задачи применялся ряд под-
ходов, в соответствии с которыми преобразовывалась матрица системы (15),
такие как метод диагонализации матриц масс, введение обезразмеривающих ко-
эффициентов и перестановка строк и столбцов по алгоритму Катхилла-Макки.
В третьей главеприводится описание разработанного программного ком-
плекса для расчета неизотермического течения в пороупругой среде с уче-
том разрушения породы в рамках разработанного численного алгоритма. Про-
граммный комплекс реализован на языке программирования С++ и состоит из
трех основных компонентов: препроцессор, вычислительное ядро и постпроцес-
сор.
Блок препроцессинга отвечает за построение расчетной сетки, зачитывание
и предобработки параметров модели, а также задание начальных и граничных
условий. В вычислительном ядре производится вычисление якобиана и правой
части, итерационное решение системы нелинейных уравнений, а также решение
системы линейных уравнений. В постпроцессоре производится выгрузка основ-
ных результатов (давление, температура, деформации, напряжения, параметр
повреждаемости, компоненты энергии и так далее) на каждый момент времени
для последующей визуализации и анализа.
В четвертой главеприведены результаты моделирования с использова-
нием разработанного программного комплекса. Валидация алгоритма проводи-
лась на ряде тестов (задача Терцаги, тест Манделя и тест на одномерное неизо-
термическое расширение) для которых известно аналитическое решение. Кроме
того, был проведен ряд расчетов, моделирующих воздействие на пласт добы-
вающих и нагнетательных скважин при различных условиях, с целью оценки
влияния геомеханических эффектов. В данных расчетах учитывались различ-
ные эффекты, характерные для процесса термического воздействия, такие как
изменение проницаемости при деформации породы, изменение вязкости флюи-
да, а также разрушение пласта. Для описания эволюции параметра повреждае-
мости была выбрана следующая зависимость, описанная в литературе, которая
основана на результатах лабораторных экспериментов по разрушению горных
пород:
⎧
⎪
⎪
⎪
⎪0,если ˜ < ˜ ,⎯ ⎪⎸ 3 ⎨ off ⎪ , если ˜ ⩽ ˜ ⩽ ˜off , ⎸∑︁ = ˜ − off ˜ = ⎷ ⟨ ⟩2 , (16) ⎪ ˜off − ˜ ˜off − ˜ ⎪ =1 ⎪ ⎪ ˜off если ˜ > ˜off ,
⎪
⎩ lim − ( lim − off )
⎪
˜

где ⟨ ⟩ = ( + | |) /2, — главные деформации.
В одном из проведенных расчетов моделируется реакция пласта при исполь-
зовании системы поддержания пластового давления, при котором нагнетатель-
ная скважина закачивает флюид при температуре 400 ∘ C с постоянной приеми-
стостью, а добывающая скважина работает с постоянным забойным давлением.
Пласт имеет размеры 100 × 100 × 1 м, скважины расположены на противопо-
ложных концах диагонали пласта. Полученные по результатам моделирования
распределения давления, температуры, параметра повреждаемости и проница-
емости через 6, 12 и 60 месяцев приведены на рисунке 1.
Также был проведен расчет данной модели без учета разрушения пласта.
Сравнение распределений компоненты тензора деформаций на момент вре-
мени 5 лет для случаев учета и в отсутствие учета повреждаемости приведено на
рисунке 2. Из данного рисунка видно, что учет разрушения породы существенно
влияет на расчет напряженно-деформированного состояния (для данного теста
различие в значениях тензора деформаций достигает 20%).
Дополнительно был проведен ряд расчетов по оценке зоне трещиноватости
в трехмерном пласте при закачке теплоносителя (рисунок 3) в условиях ани-
зотропии внешних напряжений.
В заключениисформулированы основные результаты работы.
Рис. 1. Распределение давления (1 ряд), температуры (2 ряд), параметра
повреждаемости (3 ряд), проницаемости (4 ряд) через 6 месяцев (слева),
12 месяцев (по центру) и 5 лет (справа).

Рис. 2. Компонента тензора деформаций на момент времени 5 лет при учете
(слева) и в отсутствие учета (справа) повреждаемости.
Рис. 3. Развитие зоны трещиноватости в трехмерном пласте с анизотропией
напряжений.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

Основными результатами выполненной работы являются:

1. Термодинамически согласованная физико-математическая модель разру-
шения пороупругой среды, учитывающая деформационные, фильтраци-
онные и неизотермические эффекты, пригодная для анализа эффектив-
ности современных и перспективных тепловых методов увеличения неф-
теотдачи.
2. Неявный численный алгоритм расчета эволюции термопороупругой среды
с учетом разрушения на основе метода конечных элементов с применением
неструктурированных тетраэдральных сеток.
3. Программный комплекс для моделирования термического воздействия на
пороупругую среду с учетом разрушения, а также результаты проведен-
ных валидационных и тестовых расчетов, подтверждающих корректность
разработанных алгоритмов.