Применение метода Галеркина с разрывными базисными функциями к исследованию динамики изменения температуры и давления в пласте с нагнетательной скважиной и трещиной гидроразрыва

Полный текст

Введение

В настоящее время в связи с вводом в эксплуатацию месторождений с трудноизвлекаемыми запасами и значительной выработкой многих крупных месторождений развитие нефтегазовой промышленности России происходит на фоне внушительного падения запасов нефти и газа. Одним из важнейших факторов, оказывающих влияние на извлечение нефти из месторождения, является состояние призабойной зоны пласта (ПЗП). Важным источником информации о ПЗП являются гидродинамические исследования пластов и скважин (ГДИС).

В устоявшихся методах ГДИС анализируются кривые давления в бесконечном пласте при неустановившемся режиме радиальной фильтрации. Основные подходы (анализ данных по кривой падения давления и по кривой восстановления давления) базируются на решении уравнения пьезопроводности. Однако для более полного исследования скважин очень важно рассматривать, наряду с методами ГДИС, методы термометрии¹. Для скважин с гидравлическим разрывом пласта такие исследования особенно важны. Отсюда возникает необходимость в разработке математической модели для системы «скважина – трещина – пласт»². Использование методов термометрии скважин и пластов на сегодняшний день позволяет увеличить нефтеотдачу пластов за счет более эффективных мер по увеличению нефтедобычи.

Обзор литературы

Настоящая работа посвящена математическому моделированию процесса изменения температурного поля и поля давления в пласте с трещиной гидроразрыва под действием нагнетания охлаждающей жидкости в вертикальную скважину [1]. Для описания математической модели данного процесса используются уравнения конвекции-диффузии. В настоящее время существует множество подходов к решению этих уравнений. Одним из перспективных и активно развивающихся является метод Галеркина с разрывными базисными функциями [2–4], который прекрасно зарекомендовал себя для решения уравнений конвективного типа [5–8]. Также активно развиваются подходы к созданию лимитеров повышенного порядка точности, которые обеспечивают монотонность решения, полученного с помощью разрывного метода Галеркина [9–11]. Дальнейшее развитие метода Галеркина с разрывными базисными функциями привело к его модификации с использованием разнесенных сеток (Staggered Discontinuous Galerkin Method), которая объединяет хорошие качества этих способов [12–15]. К примеру, в ряде работ был построен оригинальный вычислительный алгоритм, в котором вспомогательные переменные, введенные для понижения порядка исходных уравнений переноса тепла, рассчитываются на двойственной сетке, представленной в виде медианных контрольных объемов вокруг узлов основной сетки [16–19]. Искомые величины аппроксимируются на основной неструктурированной треугольной сетке.

Материалы и методы

Для описания динамики изменения температуры, скорости и давления будем рассматривать следующие уравнения, которые подробно представлены в одной из наших работ³.

Процесс переноса тепла:

  $c\rho \frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(\lambda \frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(\lambda \frac{\partial T}{\partial y}\right)-c\rho \left(u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}\right),\left(x,y\right)\in D\mathrm{,0}<t\le T_{\text{max}},$ (1)

$T\left(x,y\mathrm{,0}\right)=T_0\mathrm{,0}<t\le T_{\text{max}},$

$-\lambda \frac{\partial T}{\partial n}=\mathrm{0,}\left(x,y\right)\in \partial D\mathrm{,0}<t\le T_{\text{max}},$

$\begin{array}{c}-\lambda \frac{\partial T}{\partial n}=\beta \left(T-T_{\Gamma }\right)n_x+\beta \left(T-T_{\Gamma }\right)n_y,\\ \left(x,y\right)\in \Gamma ,0<t\le T_{\text{max}}.\end{array}$

 

Процесс изменения давления [20]:

  $\begin{array}{c}{c}_r\frac{\partial p}{\partial t}-div\left(\frac{\kappa }{\mu }grad\left(p\right)\right)=\mathrm{0,}\\ \left(x,y\right)\in D,0<t\le T_{\text{max}},\end{array}$    (2)

$p\left(x,y\mathrm{,0}\right)=p_0\left(x,y\right)\mathrm{,0}<t\le T_{\text{max}},$

$-\frac{\kappa }{\mu }\frac{\partial p}{\partial n}=\mathrm{0,}\left(x,y\right)\in \partial D,0<t\le T_{\text{max}},$

$p=p_{\Gamma },\left(x,y\right)\in \Gamma ,0<t\le T_{\text{max}}.$

 

Для описания скорости течения жидкости используется закон Дарси⁴:

 $u=-\frac{\kappa }{\mu }grad\left(p\right),\left(x,y\right)\in D.$       (3)

Ранее авторами был представлен подробный вывод уравнений для решения систем (1)–(3) на неструктурированной двойственной сетке, здесь приведем лишь полученные выражения⁵.

Для решения системы (1) используются выражения:

$\begin{array}{r}c\rho {{\displaystyle \sum }}_{i=0}^2\frac{d{T}_{ij}}{dt}{{\displaystyle \int }}_{K_j}{\varphi }_i^j{\varphi }_k^{j}dS=-{\displaystyle {\oint }_{\partial K_j}{n}_x{\omega }_x^{\Gamma }{\varphi }_k^{j}dl}-\\ -{\displaystyle {\oint }_{\partial K_j}{n}_y{\omega }_y^{\Gamma }{\varphi }_k^{j}dl}+{{\displaystyle \int }}_{K_j}{\omega }_x\frac{\partial {\varphi }_k^j}{\partial x}dS+{{\displaystyle \int }}_{K_j}{\omega }_y\frac{\partial {\varphi }_k^j}{\partial y}dS-\\ -c\rho \left({\displaystyle {\oint }_{\partial K_j}{\left(uT\right)}^{\Gamma }{n}_x{\varphi }_k^{j}dl}+{\displaystyle {\oint }_{\partial K_j}{\left(vT\right)}^{\Gamma }{n}_y{\varphi }_k^{j}dl}\right)+\\ +c\rho \left({{\displaystyle \int }}_{K_j}{T}_j\frac{\partial \left(u{\varphi }_k^j\right)}{\partial x}dS+{{\displaystyle \int }}_{K_j}{T}_j\frac{\partial \left(v{\varphi }_k^j\right)}{\partial y}dS\right),\\ \forall {\varphi }_k^j\left(x,y\right),k=\mathrm{0...2,}\left(4\right)\end{array}$

 

   ${{\displaystyle \sum }}_{i=0}^2{\omega }_{xij}{{\displaystyle \int }}_{D_j}{\psi }_i^j{\psi }_k^{j}dS=-{\displaystyle {\oint }_{\partial D_j}{n}_x\lambda T^{\Gamma }{\psi }_k^{j}dl}+{{\displaystyle \int }}_{D_j}T\lambda \frac{\partial {\psi }_k^j}{\partial x}dS,$         (5)

${{\displaystyle \sum }}_{i=0}^2{\omega }_{yij}{{\displaystyle \int }}_{D_j}{\psi }_i^j{\psi }_k^{j}dS=-{\displaystyle {\oint }_{\partial D_j}{n}_y\lambda T^{\Gamma }{\psi }_k^{j}dl}+{{\displaystyle \int }}_{D_j}T\lambda \frac{\partial {\psi }_k^j}{\partial y}dS,$

 $\forall {\psi }_k^j\left(x,y\right),k=0\dots 2.$           (6)

Для решения систем (2), (3) используются выражения:

$\begin{array}{l}{c}_r{{\displaystyle \sum }}_{i=0}^2\frac{d{p}_{ij}}{dt}{{\displaystyle \int }}_{K_j}{\varphi }_i^j{\varphi }_k^{j}dS={\displaystyle {\oint }_{\partial K_j}{n}_x{u}^{\Gamma }{\varphi }_k^{j}dl}+{\displaystyle {\oint }_{\partial K_j}{n}_y{v}^{\Gamma }{\varphi }_k^{j}dl}-{{\displaystyle \int }}_{K_j}u\frac{\partial {\varphi }_k^j}{\partial x}dS-{{\displaystyle \int }}_{K_j}v\frac{\partial {\varphi }_k^j}{\partial y}dS,\\ \forall {\phi }_k^j\left(x,y\right),k=0\dots \mathrm{2,}\left(7\right)\end{array}$

 

${{\displaystyle \sum }}_{i=0}^2{u}_{ij}{{\displaystyle \int }}_{K_j}{\phi }_i^j{\phi }_k^{j}dS=-{\displaystyle {\oint }_{K_j}{n}_x}\frac{\kappa }{\mu }{p}^{\Gamma }{\phi }_k^{j}dl+{{\displaystyle \int }}_{K_j}p\frac{\kappa }{\mu }\frac{\partial {\phi }_k^j}{\partial x}dS,$    

$\forall {\phi }_k^j\left(x,y\right),k=0\dots 2$           (8)

${{\displaystyle \sum }}_{i=0}^2{v}_{ij}{{\displaystyle \int }}_{K_j}{\phi }_i^j{\phi }_k^{j}dS=-{\displaystyle {\oint }_{\partial K_j}{n}_y}\frac{\kappa }{\mu }{p}^{\Gamma }{\phi }_k^{j}dl+{{\displaystyle \int }}_{K_j}p\frac{\kappa }{\mu }\frac{\partial {\phi }_k^j}{\partial y}dS,$

$\forall {\phi }_k^j\left(x,y\right),k=0\dots \mathrm{2,}$              (9)

где $\left\{{\varphi }_i^j\left(x,y\right)\right\}$  , $\left\{{\psi }_i^j\left(x,y\right)\right\}$  – системы базисных функций, заданные на элементе K_j (элементы основной треугольной сетки) и D_j (элементы двойственной сетки) соответственно, в виде проекции на которые находятся температура, давление и компоненты вектора скорости⁶.

Для нахождения величин (uT)^Г и (vT)^Г на границах элементов в системе (4) используется потоковая функция Лакса – Фридрихса⁷. При вычислении потоковых величин $T^{\Gamma },{\omega }_x^{\Gamma },{\omega }_y^{\Gamma }$  на границе элементов в системах (4)–(6) применяется потоковая функция [21]. С учетом использования двойственных сеток потоковые величины представляются в виде:

Т^Г = Т,

${\omega }_x^{\Gamma }={\omega }_x-C_{11}\left(T^{+}-T^{-}\right)n_x$ ,

${\omega }_y^{\Gamma }={\omega }_y-C_{11}\left(T^{+}-T^{-}\right)n_y$ ,

где T⁺ – значение температуры из ячейки, для которой нормаль n = (n_x, n_y) является внешней, а T⁺ – значение температуры из ячейки, для которой нормаль n = (n_x, n_y) является внутренней; C₁₁ – стабилизирующая добавка.

На граничных ребрах, с учетом типа граничного условия, получаем:

$T^{Г}=T^{+}$ ,

${\omega }_x^{\Gamma }=\beta \left(T-T_{\Gamma }\right)n_x$ ,

${\omega }_y^{\Gamma }=\beta \left(T-T_{\Gamma }\right)n_y$ .

Для вычисления потоковых величин p^Г, u^Г, v^Г в системах (7)–(9) также используются стабилизирующие добавки, но в данном случае аппроксимация строится только на треугольной сетке. Вид потоковых функций представлен ниже:

$p^{\Gamma }=\frac{\left(p^{+}+p^{-}\right)}{2}$ ,

$u^{\Gamma }=\frac{\left(u^{+}+u^{-}\right)}{2}-C_{11}\left(p^{+}-p^{-}\right)n_x$ ,

$v^{\Gamma }=\frac{\left(v^{+}+v^{-}\right)}{2}-C_{11}\left(p^{+}-p^{-}\right)n_y$ ,

где верхний индекс «+» обозначает величины из ячейки, для которой нормаль n = (n_x, n_y) является внешней, а верхний индекс «–» – величины из ячейки, для которой нормаль n = (n_x, n_y) является внутренней; C₁₁ – стабилизирующая добавка.

На граничных ребрах получаем следующий вид:

$p^{\Gamma }=p^{+}$ ,

$u^{\Gamma }=q{n}_x$ ,

$v^{\Gamma }=q{n}_y$ .

В системах (4)–(9) необходимо с высокой точностью вычислять поверхностные и контурные интегралы. Для этого используются квадратурные формулы Гаусса [22]. Поверхностные интегралы вычисляются по трем точкам, контурные интегралы вычисляются с использованием двухточечного шаблона. Для подавления нефизических осцилляций используется лимитер TVD⁸. Для аппроксимации по времени используется явная схема Эйлера.

Результаты исследования

Описанный вычислительный алгоритм был реализован в виде программного пакета для расчета динамики изменения температуры и давления в нефтеносном пласте. Для сокращения времени расчетов была использована технология параллельных вычислений MPI.

Для анализа полученных результатов рассматривалась следующая постановка задачи: T₀ = 363 K, c = 2 000 Дж/кг·К, ρ = 950 кг/м³, λ = 2,5208, µ = 0,315 ∙ 10^–4 Па∙с, p₀ = 2,5 ∙ 10⁷ Па, p_Г = 2,5 ∙ 10⁷ Па, β = 150 Вт/м²·K. Для трещины были заданы следующие значения параметров: c_r = 4,18968 ∙ 10^–9 Па^–1, κ = 2,96 ∙ 10^–15 м². Для пласта были заданы следующие значения параметров: c_r = 0,2505 ∙ 10^–9 Па^–1, κ = 5 ∙ 10^–17 м².

Рассматривается область длиной 60 м и шириной 10 м. В центре области находится скважина с радиусом 0,025 м. Слева и справа к скважине симметрично примыкают трещины длиной 8 м и шириной 0,005 м каждая.

На рисунке 1 представлена расчетная сетка для описанной задачи. Вдоль трещины наблюдается заметное сгущение сетки. Расчетная сетка содержит 36 339 ячеек основной неструктурированной сетки.

 

 

 

Рис. 1. Сетка

Fig. 1. Mesh

 

Расчет производился с использованием параллельного комплекса программ на 12 процессорах [23]. На рисунке 2 представлена декомпозиция расчетной области по процессорам.

 

 

 

Рис. 2. Декомпозиция расчетной области

Fig. 2. Decomposition of the computational domain

 

На практике наибольший интерес представляет состояние призабойной зоны пласта. В связи с этим, а также из-за большого масштаба задачи дальнейшие рисунки представляют не всю расчетную область, а ее часть, приближенную к скважине. На рисунке 3 представлена расчетная сетка возле скважины, диаметр скважины 0,05 м.

 

 

 

Рис. 3. Сетка возле скважины

Fig. 3. Mesh near the well

 

На рисунках 4–6 представлено распределение поля давления в различные моменты времени. Из рисунков видно, что с течением времени вдоль трещины давление растет заметно быстрее, по сравнению с пластом, что согласуется с заданными параметрами задачи.

 

 

 

Рис. 4. Распределение давления, t = 1 с

Fig. 4. Distribution of the pressure, t = 1 s

 

 

 

Рис. 5. Распределение давления, t = 5 с

Fig. 5. Distribution of the pressure, t = 5 s

 

 

 

 

Рис. 6. Распределение давления, t = 10 с
Fig. 6. Distribution of the pressure, t = 10 s

 

На рисунках 7–9 представлены картины распределения поля температуры в различные моменты времени. Из рисунков видно, что холодная закачиваемая через вертикальную нагнетательную скважину жидкость охлаждает пласт. Можно отметить, что вдоль трещины охлаждение происходит немного интенсивнее, что согласуется с наблюдаемой картиной распределения давления. Существенное уменьшение температуры наблюдаются вблизи скважины и вдоль трещины, в частности на ее створках.

 

 

Рис. 7. Распределение температуры, t = 1 с
Fig. 7. Distribution of the temperature, t = 1 s

 

 

 

Рис. 8. Распределение температуры, t = 5 с
Fig. 8. Distribution of the temperature, t = 5 s

 

 

 

Рис. 9. Распределение температуры, t = 10 с
Fig. 9. Distribution of the temperature, t = 10 s

 

 

Обсуждение и заключение

В настоящей статье разработан и реализован вычислительный алгоритм для моделирования динамики изменения температуры и давления в нефтеносном пласте. Алгоритм построен на основе метода Галеркина с разрывными базисными функциями на разнесенных неструктурированных сетках с применением технологии параллельных вычислений MPI. С использованием разработанного программного кода была исследована задача закачивания в пласт охлаждающей жидкости через вертикальную нагнетательную скважину. Можно сделать вывод, что результаты моделирования показывают адекватные картины для температурного поля и поля давления в пласте, соответствующие заданным начально-краевым условиям. Для более точного моделирования рассматриваемого процесса в дальнейшем планируется решать данную задачу в трехмерной постановке на неструктурированных тетраэдральных сетках.  

 

 

¹           Гидродинамический разрыв пласта / Д. С. Кузнецов [и др.]. Томск, 2008. 114 с.

²           Чекалюк Э. Б. Термодинамика нефтяного пласта. М.: Недра, 1965. 238 с.; Руководство по исследованию и интерпретации. Термодинамические исследования при различных режимах работы скважин. Уфа, 2002. 248 с.

³           Жалнин Р. В., Масягин В. Ф., Пескова Е. Е. Применение разрывного метода Галеркина для математического моделирования динамики распространения температуры в пласте с нагнетательной скважиной // Аналитические и численные методы моделирования естественно-научных и социальных проблем: материалы XIV Междунар. науч.-техн. конф., 3–6 декабря 2019 г., Пенза / под ред. д-ра физ.-мат. наук, проф. И. В. Бойкова. Пенза: Изд-во ПГУ, 2019. С. 54–61. URL: https://dep_vipm.pnzgu.ru/files/dep_vipm.pnzgu.ru/konference/achm_2019.pdf (дата обращения: 12.02.2021).

⁴           Там же.

⁵           Там же.

⁶           Там же.

⁷           Shu C.-W. Numerical Methods for Hyperbolic Conservation Laws // Conference Proceedings (AM 257), 2007. 2007. 32 p. URL: https://mathema.tician.de/dl/academic/notes/257/257.pdf (дата обращения: 11.02.2021).

⁸           Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. 2-е изд., испр. и доп. М.: Физматлит, 2012. 656 с.

 

Об авторах

Елизавета Евгеньевна Пескова

ФГБОУ ВО «МГУ им. Н. П. Огарёва»

Email: e.e.peskova@mail.ru
ORCID iD: 0000-0003-2618-1674
ResearcherId: U-7971-2019

младший научный сотрудник, доцент кафедры прикладной математики, дифференциальных уравнений и теоретической механики факультета математики и информационных технологий, кандидат физико-математических наук

Россия, 430005, г. Саранск, ул. Большевистская, д. 68

Руслан Викторович Жалнин

ФГБОУ ВО «МГУ им. Н. П. Огарёва»

Email: zhrv@mrsu.ru
ORCID iD: 0000-0002-1103-3321
ResearcherId: Q-6945-2016

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

Россия, 430005, г. Саранск, ул. Большевистская, д. 68

Виктор Федорович Масягин

ФГБОУ ВО «МГУ им. Н. П. Огарёва»

Автор, ответственный за переписку.
Email: vmasyagin@mrsu.ru
ORCID iD: 0000-0001-6738-8183
ResearcherId: C-2439-2013

старший научный сотрудник, доцент кафедры прикладной математики, дифференциальных уравнений и теоретической механики факультета математики и информационных технологий, кандидат физико-математических наук

Россия, 430005, г. Саранск, ул. Большевистская, д. 68

Владимир Федорович Тишкин

ИПМ им. М. В. Келдыша РАН

Email: v.f.tishkin@mail.ru
ORCID iD: 0000-0001-7295-7002
ResearcherId: R-5820-2016

заведующий отделом, член-корреспондент РАН, доктор физико-математических наук, профессор,

Россия, 125047, г. Москва, Миусская пл., д. 4

Список литературы

Дополнительные файлы