CC BY
89
ТЕХНИЧЕСКИЕ НАУКИ
УДК 519.63 А. Е. Колесов
ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ ПОРОУПРУГОСТИ ПРИ РАЗРАБОТКЕ НЕФТЯНЫХ МЕСТОРОЖДЕНИЙ
Рассматриваются проблемы численного решения задач пороупругости (фильтрационной консолидации), которые возникают при исследовании деформационных процессоров, происходящих при разработке месторождений нефти и газа. Система уравнений пороупругости включает уравнения Ламе для перемещений пористой среды и уравнения фильтрации для давления жидкости. Важнейшей особенностью математических моделей пороупругости является сильная связь между данными уравнениями. Уравнение упругости включает объемную силу, которая пропорциональна градиенту давления, а уравнение фильтрации содержит сжимаемость среды, пропорциональную дивергенции скорости перемещений. Вычислительный алгоритм основан на конечно-элементной дискретизации по пространству. Для аппроксимации по времени используются стандартные двухслойные связанные схемы с весами и трехслойные явно-неявные схемы расщепления по физическим процессам, когда переход на новый временной слой осуществляется решением отдельных задач для перемещений и давления. Для используемой схемы расщепления по физическим процессам дано условие его безусловной устойчивости. Представлены результаты численного решения трехмерной модельной задачи. Проведена оценка работоспособности схем расщепления по физическим процессам. Особое внимание уделено использованию высокопроизводительных вычислительных систем кластерной архитектуры. Вычислительная реализация осуществлена с использованием свободной библиотеки FEniCS для решения уравнений в частных производных методом конечных элементов.
Ключевые слова: пороупругость, фильтрационная консолидация, оседание земной поверхности, разработка нефтяных месторождений, метод конечных элементов, конечно-разностные схемы, схемы расщепления, регуляризация, высокопроизводительные системы, параллельные вычисления.
A. E. Kolesov
Numerical Solution of Poroelasticity Problems Associated with the Development of Oil Fields
We conisder the numerical solution of poroelasticity (filtration consolidation) problems, which arise in analysis of deformation processes during the development of oil and gas fields. The underlying system of equations consists of the Lame equations for the displacements and flow equation for the pressure. The most important feature of mathematical models for poroelasticity is that the system consists of coupled equations.
КОЛЕСОВ Александр Егорович - н. с. научно-исследовательской кафедры «Вычислительные технологии» ИМИ СВФУ им. М. К. Аммосова. E-mail: [email protected]
KOLESOVAlexandr Egorovich - research associate of research department of Computational technologies, Institue of mathematics and informatics, North-Eastern Federal University named after M. K. Ammosov. E-mail: [email protected]
The displacement equation in a poroelasticity problem includes a body force, which is proportional to the pressure gradient. The pressure equation includes a term describing the compressibility of the medium (the divergence of the displacement velocity). The numerical algorithm is based on a finite-element approximation in space. For time discretization we use the standard two-level coupled schemes with weights and splitting schemes with respect to physical processes, in which the transition to a new time level is associated with solving separate problems for the displacements and pressure. The results of numerical simulation of a 3D model problem are presented. The evaluation of performance of the splitting schemes with respect to physical processes are presented. Special attention is given to using high performance computing systems. Computational implementation is performed using the FEniCS library for solving partial differential equations.
Keywords: poroelasticity, flow consolidation, subsidence, oil field development, finite element methods, operator-difference schemes, splitting schemes, regularization, high performance systems, parallel computations.
Наличие движущейся жидкости в пористой горной породе оказывает существенное влияние на его механическое состояние [1]. В связи с этим недооценка данного процесса при разработке месторождений нефти и газа может привести к разрушающим явлениям со значительным экономическим и экологическим ущербом [2]. Деформации играют важную роль в процессах сжатия продуктивных пластов, оседания земной поверхности и морского дна, смятия обсадных колонн скважин и т. д. [3]. Также деформация горных пород может вызвать существенное изменение их фильтрационно-емкостных свойств [4].
Теория пороупругости изучает совместный механизм течения жидкости и деформирования пористых сред (фильтрационной консолидации) [5]. При этом базовые математические модели фильтрации состоят из уравнений упругости Ламе для перемещений среды и нестационарного уравнения фильтрации для давления жидкости. Важнейшей особенностью математических моделей пороупругости является сильная связь между данными уравнениями. Так, уравнение упругости включает объемную силу, которая пропорциональна градиенту давления, а уравнение фильтрации, в свою очередь, содержит сжимаемость среды, пропорциональную дивергенции скорости перемещений.
Вычислительные алгоритмы приближенного решения задач пороупругости базируются чаще всего на конечно-элементной аппроксимации по пространству [6-7]. А при аппроксимации по времени используются стандартные двухслойные схемы с весами [8]. Отдельного внимания заслуживают исследования схем расщепления по физическим процессам, когда переход на новый временной слой осуществляется последовательным решением отдельных задач расчета перемещений и расчета давления. Различные классы таких схем строятся на основе аддитивного расщепления оператора задачи [9-10].
В нашей работе численно исследуется схема расщепления по физическим процессам для задачи пороупругости [11-12], построенная на основе принципа регуляризации операторно-разностных схем А. А. Самарского [13-14]. В качестве модельной используется задача оседания нефтяных месторождений.
Постановка задачи
Математическая модель пороупругости состоит из уравнения фильтрации для давления р и уравнения Ламе для перемещений среды и:
Здесь а - тензор напряжения, который для изотропного тела задан через закон Гука:
а = + АсСгу и1,
где X - модули упругости, называемые коэффициентами Ламе, I - единичный тензор, е -тензор деформации:
Введение
div a(u) — agrad p = 0,
(1)
a
(2)
£ (и) = 1 (grad и + grad ит).
Также в уравнениях (1) и (2) а - коэффициент Био, характеризующий связь между смещением и давлением, в=1/М, М - модуль Био, k - тензор абсолютной проницаемости, п - коэффициент динамической вязкости жидкости, / - функция, описывающая интенсивность заданных источников (стоков) жидкости.
Учитывая специфику прикладной области, рассмотрим более подробно коэффициент Био а и модуль Био М [15]:
К 1 ф а-ф
а = 1--, — = -?— +-,
К М К, К
s .г s
где К, К , К{ - модули объемного сжатия пористой среды, твердой фазы (скелета) и жидкости, соответственно ф - коэффициент пористости, характеризующий относительное количество пор и твердой фазы.
Задача для системы уравнений (1), (2) рассматривается в ограниченной области О, на границе которой Г задаются следующие однородные условия первого и второго рода для перемещений:
и = 0, х е ГД, ип = 0, х е Г^, (3)
для давления:
р = 0, х е Г рв, - к др = 0, х еГр, (4)
П дп
где п — единичная нормаль к границе, Г = Тшв + Г^,= Грв + ГN, Г„ = ГN, Г^ = Грв.Кроме этого, задается начальное условие для давления:
Р (х,0) = Ро (х), х еО. (5)
Вычислительный алгоритм
Для численного решения задачи (1)-(5) получим дискретную задачу. Для аппроксимации по пространству используем метод конечных элементов. Перед проведением аппроксимации нам необходимо получить вариационную формулировку задачи. Для этого определим следующие пробные и тестовые функции для скалярных и векторных функций, соответственно:
Q = { е Я1 (О): q (х) = 0, X е грв },
V = { еН1 (О): V(х) = 0, X е г;}.
Умножив (1) и (2) на V е V и q е Q соответственно и проинтегрировав полученные уравнения по частям, получим
Ja(u)£(v)dx + ^а (gradр, ^)сСх = 0,
а а
Уа и qdx + ^fl — qdx + — gradp,grad q dx = ^fqdx,
а ^ а ^ п П п
где и е V, р е V. Для удобства введем следующие билинейные формы:
а (и, у) = ^а(ии)£(у)<1х, g (р, V) = а^gradpvdx, п п
с (р, q) = в^pqdx, d(и, q) = а ^div uqdx,
п п
■к
/к
— grad р qdx. а П
Тогда вариационная постановка читается: найти и е V, р е Q такие, что
a (u, v) + g (p, v) = 0, Vv € V,
d
du dt '
dp dt
+ b (p, q) = ( f, q), Vq G Q.
(6) (7)
Начальное условие (5) выражается как
(p(0), q) = (p0, q), q € Q. (8)
Далее проведем конечно-элементную аппроксимацию по пространству. Для этого сначала построим расчетную сетку Qh={m , т , ..., mN} в области Q. Здесь N - количество ячеек в Q h = maxhm, где hm - диаметр окружности, вписанной в ю. На построенной сетке определим [16] следующие конечно-размерные дискретные подпространства пробных и тестовых функций: Qh с Q, Vh С V. Чтобы получить дискретную задачу, ограничим вариационную задачу (6), (7) для данных подпространств [17]: найти uh е Vh, ph £ Qh такие, что
a U, v) + g (ph, v) = 0, Vv e Vh,
d
duh dt
dPh dt
+ b (Ph, q) = (f, q), Vq e Qh.
(9) (10)
Берем базисы {фj и {\у1 ^ для Ул и Qh соответственно, где N - размерности подпространств. Теперь пробные функции и и р можно аппроксимировать как:
: > Ph = £PiV,,
j=i
а тестовые функции берем и= ф., /=1,2, ..., Ы, д= £=1,2, ..., М. Здесь и. и р1 - векторы неизвестных, которые будут вычисляться. С учетом этого запишем (9), (10) как
N ^
^а () (V,)р, = 0, 1 =
У=1
.du, JX , dp,
Yd (ф Wk + Yc (Vi Wk )~dj- + Yb (Vt ,Vk )pt = (f ,Vk)
1=1
dt
k=1,2, ..., N .
Введя следующие линейные операторы (матрицы)
A = Ya(фф), G = Yg(Vlф) i = 1,2,...,N
(11) (12)
1=1
N N
D = ХУ() с = ^с(VI,¥к) в = (VI,¥к) к = 1,2,•••,N
j=l 1=1 1=1
перейдем к следующим дифференциально-операторной системе уравнений:
Au+Gp=0,
В — + С^- + Вр = f, dt dt
где вектор f = (fк ,k=1, 2, ..., N. Начальное условие возьмем в виде
" Р (0)=р0. (13)
Для дискретизации по времени будем использовать для простоты равномерную сетку
с шагом т>0. Пусть и"=ик(х,Г), р"=рА(х,Г), где У=т, и=0,1, ... Для приближенного решения
задачи (11)-(13) будем использовать двухслойную схему с весами:
Аии+^ри+1=0, (14)
где в - весовой параметр, 0<9<1, и
N
N
рв+ = в р'+ + (1-в)р', /в" + = /(х,).
Вычислительная реализация схемы (14)-(15) связана с определением м"+', р"+' на каждом временном слое и=0,1, ...как решения связанной системы уравнений, которая требует использования специальных вычислительных алгоритмов (см., например, [18]). В силу этого повышенное внимание уделяется построению схем расщепления [10]. В данной работе будем использовать трехслойную явно-неявную схему с весами [11]:
Au"+l+Gp"=0, (16)
{ п+1 п п п—1 '
Р_Р_+(1—е)Р-—Р—
БиП+1 — ип + с
-Врп+1 = Г+1. (17)
В этом случае сперва решается уравнение (16), в котором давление берется из предыдущего временного слоя, и находятся перемещения и"+'. Найденные перемещения используются при определении давления на новом временном слое р"+' (17).
Заметим, что в схеме (16)-(17) весовой параметр в отличается от веса в (14)-(15) и выбирается из условия устойчивости [11].
Утверждение 1. Схема расщепления (16), (17) устойчива при
2в>1+у (18)
В условии устойчивости (18) у определяется из условия как максимальное собственное значение спектральной задачи у=^ш!а:
GDu=MCu. (19)
Численные эксперименты
Рассмотрим численное решение трехмерной задачи пороупругости при добыче нефти из участка месторождения с протяженностью 1000,0 м по горизонтальным направлениям и 200,0 м по вертикальному направлению (рис. 1). В середине пласта пробурена добывающая скважина с радиусом г=0,1 м. Скважина несовершенная по степени вскрытия, т. е. она вскрывает пласт не на всю толщину, а на й=150,0 м от кровли пласта.
Пластовое давление в начальный момент времени равно р0=40,0 МПа. Добыча
Рис. 1. Расчетная область
производится при постоянном забойном давлении р0=30,0 МПа. Свойства продуктивного пласта и жидкости представлены в таблице 1.
Расчеты проводились на трех вычислительных тетраэдальных сетках с локальным сгущением вокруг скважины, которая вырезана из расчетной области (рис. 2). Сетки построены с помощью программы Gmsh. Количество вершин, ячеек и количество неизвестных в полях давления р и смещения и для каждой сетки представлены в таблице 2.
Таблица 1
Входные параметры
Параметр Символ Размерность Значение
Модуль сдвига М ГПа 0,284
Коэффициент Ламе Я ГПа 0,646
Модуль объемного сжатия пористой среды К ГПа 0,646
Модуль объемного сжатия твердой фазы К Б ГПа 40,0
Модуль объемного сжатия жидкости К ГПа 0,33
Коэффициент Био а - 0,98
Коэффициент пористости Ф - 0,2
Модуль Био М ГПа 1,598
Проницаемость к мД 10,0
Вязкость п Па-с 0,001
Таблицы 2
Параметры сеток
Сетка 1 Сетка 2 Сетка 3
Количество вершин 25107 84571 299713
Количество ячеек 108575 344867 1361790
Количество неизвестных р 25107 84571 299713
Количество неизвестных и 75282 253713 899139
Для решения линейной системы алгебраических уравнений используем обобщенный метод минимальных невязок (GMRES) с неполной Ши факторизацией в качестве предоб-уславливателя. Расчеты проводились на вычислительном кластере Северо-Восточного федерального университета имени М. К. Аммосова «Ариан Кузьмин».
В табл. 3 представлена зависимость общего времени расчетов в секундах от количества запущенных процессоров и сетки при использовании неявной схемы (14), (15) и схемы расщепления (16), (17). Здесь в общее время входят затраты на инициализацию расчетов, загрузку и разделения расчетной сетки, подготовку и сборку матрицы и вектора линейной системы, а также расчет 10 временных шагов при шаге по времени т=1 день.
Видно, что схема расщепления значительно сокращает время расчета (около 60 %). Также при согласовании размера используемой сетки и числа процессоров обеспечивается достаточно высокая степень эффективности параллелизации. При относительно малом количестве ячеек использование большого количества процессоров может замедлить процесс расчета, так как время расходуется на обмен информацией между процессорами.
Далее оценим величины параметра у и веса в (18), при которых схема расщепления (16), (17) безусловно устойчива согласно утверждению 1. Для этого решим спектральную задачу (19). Для рассматриваемых входных параметров (табл. 1) получим следующие результаты:
Рис. 2. Расчетные сетки. Сверху вниз: сетка 1, 2, 3
Таблица 3
Общее время расчетов
Кол-во проц-ов Общее время расчетов, сек
Неявная схема Расщепленная схема
Сетка 1 Сетка 2 Сетка 3 Сетка 1 Сетка 2 Сетка 3
1 58,21 219,98 1410,73 22,93 91,45 476,11
2 33,68 122,75 837,41 13,95 53,59 289,85
4 19,24 72,24 402,92 7,86 31,67 166,35
8 12,68 45,27 374,26 5,33 23,08 123,17
16 9,96 32,50 189,24 5,27 16,17 91,31
32 - 28,72 193,43 - 11,33 45,76
64 - - 152,32 - - 34,53
у=1,67, в>1,33.
Для исследования устойчивости схемы расщепления будем смотреть значение давления вблизи скважины в точке (450,0, 500,0, 200,0), чтобы уловить начальные колебания давления в случае неустойчивости. На рис. 3 представлены графики динамики давления при различных значениях веса в для шага по времени т=1 день. При увеличении веса до 1,4 схема становится устойчивой, что соответствует теоретическим результатам. Однако в начальный момент времени все еще наблюдаются колебания значения давления. Они полностью исчезают при в=1,8.
На рис. 4 и 5 показаны распределения давления и вертикальной компоненты вектора смещения в различные моменты времени (1 месяц, 1 год, 2 года, 4 года, 6 лет) соответственно. Данное решение получено с помощью схемы расщепления с весом в=1,8 на сетке 3 при т=1 день. Смещение представлено на деформированной области и умножено 100 для наглядности. Расчеты проводились на 64 процессорах. Время расчета составило 88 мин. За 6 лет разработки пластовое давление уменьшилось до 34,3 МПа и кровля пласта опустилась почти на 1 м в области скважины и на 0,93 м на границах.
40 (50
t, (Зауь
Рис. 3. Динамика давления для различных значений веса
Заключение
Проведено численное решение задачи пороупругости при разработке месторождений нефти и газа. Вычислительный алгоритм основан на аппроксимации по пространству методом конечных элементов. Проведен сравнительный анализ стандартной двухслойной связанной схемы с весами и схемы расщепления по физическим процессам. Установлено, что схема расщепления значительно сокращает время расчетов по сравнению с связанной схемой. Расчеты были проведены на вычислительном кластере Северо-Восточного федерального университета имени М. К. Аммосова «Ариан Кузьмин». При согласовании размера используемой сетки и числа процессоров обеспечивается достаточно высокая степень эффективности параллелизации.
Работа выполнена при финансовой поддержке гранта РФФИ № 13-01-00719.
30000000 40000000
Рис. 4. Распределение давления. Сверху вниз: через 1 месяц, 1 год, 2 года, 4 года, 6 лет
Л и т е р а т у р а
1. Wang H. F. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology / H. F. Wang - Princeton University Press, 2000.
2. Allen D. R. Environmental Aspects of Oil Producing Operations- Long Beach, California / D. R. Allen // J. Pet. Technol. - 1972. - Vol. 24, no. 2 - P. 125-131.
3. Detournay E. Poroelastic response of a borehole in a non-hydrostatic stress field / E. Detournay, A.-D. Cheng // Int. J. Rock Mech. Min. Sci. Geomech. Abstr. - 1988. - Vol. 25, no. 3 - P. 171-182.
4. Settari A. Geomechanics in Integrated Reservoir Modeling / A. Settari, V. Sen // Offshore Technol. Conf.
- 2008. - P. 5-8.
5. Biot M. A. General theory of three dimensional consolidation / M. A. Biot // J. Appl. Phys. - 1941.
- Vol. 12, no. 2 - P. 155-164.
6. Lewis R. W. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media / R. W. Lewis, B. A. Schrefler - Chichester: Wiley, 1998. - 508 p.
7. Haga J. B. On the causes of pressure oscillations in low-permeable and low-compressible porous media / J. B. Haga, H. Osnes, H. P. Langtangen // Int. J. Numer. Anal. Methods Geomech. - 2012. - Vol. 36, no. 12
- P. 1507-1522.
8. Gaspar F. J. A finite difference analysis of Biot's consolidation model / F. J. Gaspar, F. J. Lisbona, P. N. Vabishchevich // Appl. Numer. Math. - 2003. - Vol. 44, no. 4 - P. 487-506.
9. Vabishchevich P. N. Additive operator-difference schemes. Splitting schemes. / P. N. Vabishchevich - de Gruyter, 2014.
10. Kim J. Stability, accuracy, and efficiency of sequential methods for coupled flow and geomechanics / J. Kim, H. Tchelepi, R. Juanes // SPE J. - 2011. - P. 2-4.
11. Вабищевич П. Н. Схема пасщепления для задач пороупругости и термоупругости / П. Н. Вабищевич, М. В. Васильева, А. Е. Колесов // Журнал вычислительной математики и математической физики - 2014.
- Т. 54, № 8. - C. 1345-1355.
12. Kolesov A. E. Splitting schemes for poroelasticity and thermoelasticity problems / A. E. Kolesov, P. N. Vabishchevich, M. V. Vasilyeva // Computers and Mathematics with Applications. - 2014. - Vol. 67, no. 12. - P. 2185-2198.
13. Самарский А. A. Теория разностных схем / А. A. Самарский - М.: Наука, 1989.
14. Самарский А. A. О регуляризации разностных схем / А. A. Самарский // Журнал вычислительной математики и математической физики - 1967. - Т. 7, № 1. - C. 62-93.
15. Fjar E. Petroleum Related / E. Fjar, R. M. Holt, P. Horsrud, A. M. Raaen, R. Risnes - Elsevier, 2008. Ed. 2. - P. 492.
16. Hughes T. J. R. The finite element method: linear static and dynamic finite element analysis / T. J. R. Hughes - DoverPublications, 2012.
17. Logg A. Automated Solution of Differential Equations by the Finite Element Method / A. Logg, K.-A. Mardal, G. N. Wells. - Springer, 2012.
18. Gaspar F. J. An efficient multigrid solver for a reformulated version of the poroelasticity system / F. J. Gaspar, F. J. Lisbona // Comput. Methods Appl. Mech. Eng. - 2007. - Vol. 196, no. 8. - P. 1447-1457.
R e f e r e n c e s
1. Wang H. F. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology / H. F. Wang - Princeton University Press, 2000.
2. Allen D. R. Environmental Aspects of Oil Producing Operations- Long Beach, California / D. R. Allen // J. Pet. Technol. - 1972. - Vol. 24, no. 2 - P. 125-131.
3. Detournay E. Poroelastic response of a borehole in a non-hydrostatic stress field / E. Detournay, A.-D. Cheng // Int. J. Rock Mech. Min. Sci. Geomech. Abstr. - 1988. - Vol. 25, no. 3 - P. 171-182.
4. Settari A. Geomechanics in Integrated Reservoir Modeling / A. Settari, V. Sen // Offshore Technol. Conf. - 2008. - P. 5-8.
5. Biot M. A. General theory of three dimensional consolidation / M. A. Biot // J. Appl. Phys. - 1941. -Vol. 12, no. 2 - P. 155-164.
6. Lewis R. W. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media / R. W. Lewis, B. A. Schrefler - Chichester: Wiley, 1998. - 508 p.
7. Haga J. B. On the causes of pressure oscillations in low-permeable and low-compressible porous media / J. B. Haga, H. Osnes, H. P. Langtangen // Int. J. Numer. Anal. Methods Geomech. - 2012. - Vol. 36, no. 12
- P. 1507-1522.
8. Gaspar F. J. A finite difference analysis of Biot's consolidation model / F. J. Gaspar, F. J. Lisbona, P. N. Vabishchevich // Appl. Numer. Math. - 2003. - Vol. 44, no. 4 - P. 487-506.
9. Vabishchevich P. N. Additive operator-difference schemes. Splitting schemes. / P. N. Vabishchevich -de Gruyter, 2014.
10. Kim J. Stability, accuracy, and efficiency of sequential methods for coupled flow and geomechanics / J. Kim, H. Tchelepi, R. Juanes // SPE J. - 2011. - P. 2-4.
11. Vabishchevich P. N. Skhema passhchepleniia dlia zadach porouprugosti i termouprugosti / P. N. Vabishchevich, M. V. Vasil'eva, A. E. Kolesov // Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki
- 2014. - T. 54, № 8. - C. 1345-1355.
12. Kolesov A. E. Splitting schemes for poroelasticity and thermoelasticity problems / A. E. Kolesov, P. N. Vabishchevich, M. V. Vasilyeva // Computers and Mathematics with Applications. - 2014. - Vol. 67, no. 12. - P. 2185-2198.
13. Samarskii A. A. Teoriia raznostnykh skhem / A. A. Samarskii - M.: Nauka, 1989.
14. Samarskii A. A. O reguliarizatsii raznostnykh skhem / A. A. Samarskii // Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki - 1967. - T. 7, № 1. - C. 62-93.
15. Fjar E. Petroleum Related / E. Fjar, R. M. Holt, P. Horsrud, A. M. Raaen, R. Risnes - Elsevier, 2008. Ed. 2. - P. 492.
16. Hughes T. J. R. The finite element method: linear static and dynamic finite element analysis / T. J. R. Hughes - DoverPublications, 2012.
17. Logg A. Automated Solution of Differential Equations by the Finite Element Method / A. Logg, K.-A. Mardal, G. N. Wells. - Springer, 2012.
18. Gaspar F. J. An efficient multigrid solver for a reformulated version of the poroelasticity system / F. J. Gaspar, F. J. Lisbona // Comput. Methods Appl. Mech. Eng. - 2007. - Vol. 196, no. 8. - P. 1447-1457.
