CC BY
128
ТЕХНИЧЕСКИЕ НАУКИ
УДК 519.63
В. И. Васильев, М. В. Васильева, Д. Я. Никифоров
РЕШЕНИЕ ЗАДАЧ ОДНОФАЗНОЙ ФИЛЬТРАЦИИ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ НА ВЫЧИСЛИТЕЛЬНОМ КЛАСТЕРЕ
Рассматривается численное решение задачи однофазной линейной фильтрации в пористой среде в многомерной постановке с учетом действия гравитационных сил. Для аппроксимации по пространственным переменным используется метод конечных элементов, который позволяет использовать неструктурированные сетки со сгущением в местах особенностей задачи. Аппроксимация по времени строится с использованием чисто неявной разностной схемы. Проведено численное исследование итерационных методов с предобуславливанием. Приводятся результаты моделирования задачи в трехмерной постановке на вычислительном кластере СВФУ «Ариан Кузьмин».
Ключевые слова: однофазная фильтрация, пористая среда, дебит, скважина, аппроксимация, метод конечных элементов, неструктурированные сетки, подпространство Крылова, вычислительный кластер, параллельные вычисления.
V. I. Vasil'ev, M. V. Vasil'eva, D. Ia. Nikiforov
Solving One Phase Filtration Problems Using Finite Element Method on Computing Cluster
In this paper, numerical solutions of one phase linear filtration problems in porous media in multidimensional formulation are considered taking into account the gravity forces effect. For the approximation over spatial variables finite element method that allows using unstructured grids with condensation in particular place
ВАСИЛЬЕВ Василий Иванович - д. ф.-м. н., проф., зав. каф. вычислительных технологий ИМИ СВФУ им. М. К. Аммосова.
E-mail: vasvasil@mail.ru
VASIL'EV Vasilij Ivanovich - Doctor of Physical-Mathematical Sciences, Professor, Head of the Department of Computational Technologies of the Institute of Mathematics and Informatics of the M. K. Ammosov NorthEastern Federal University.
ВАСИЛЬЕВА Мария Васильевна - к. ф.-м. н., с. н. с. Института вычислительной математики и математической геофизики СО РАН, г Новосибирск.
E-mail: vasilyevadotmdotv@gmail.com
VASIL'EVA Marija Vasil 'evna - PhD Candidate, Senior Researcher, Institute of Computational Mathematics and Mathematical Geophysics, Siberian Branch of RAS in Novosibirsk.
НИКИФОРОВ Дьулустан Яковлевич - аспирант кафедры вычислительных технологий ИМИ СВФУ им. М. К. Аммосова.
E-mail: dju92@mail.ru
NIKIFOROVDjulustan Jakovlevich - PhD Student of Department of Computational Technologies of the Institute of Mathematics and Informatics, M. K. Ammosov North-Eastern Federal University.
of the problem is used. The approximation over time is constructed using fully implicit difference scheme. The numerical research of iterative method with preconditioning has been done. The results of modeling of problem in three-dimensional formulation are given on the computing cluster of NEFU named "Arian Kuzmin".
Keywords: one phase filtration, porous media, debit, well, approximation, finite element method, unstructured grids, Krylov subspace, computing cluster, parallel computing.
Введение
Математическое моделирование течения жидкостей и газов в природных коллекторах связано с необходимостью расчета в сложных геометрических областях, как следствие, метод конечных элементов является основным инструментом в их численных исследованиях [1, 2].
Теоретические основы математического описания движения жидкостей и газов в пористых средах заложены в трудах [3-12]. Базовые математические модели фильтрации флюида включают уравнение неразрывности и закон Дарси. После дискретизации задачи методом конечных элементов наступает этап решения системы линейных уравнений. Для каждой задачи итерационный метод и предобуславливатель для вычислений систем уравнений выбираются индивидуально [13].
Особую сложность представляет решение трехмерной задачи. Как правило, при решении таких нестационарных задач нужно вычислять в каждый момент времени большую систему линейных уравнений, что практически вынуждает решать ее параллельно на нескольких процессах.
В данной работе рассматривается численное моделирование процессов однофазной фильтрации в пористой среде методом конечных элементов. Приводятся основные уравнения и ставятся задачи в многомерных случаях. В двумерной постановке проводятся численные исследования погрешности решения на разных неструктурированных пространственных сетках и на разных временных шагах. Также проводится численное сравнение решателей систем линейных уравнений, получающихся при дискретизации рассматриваемой задачи с учетом гетерогенности нефтяного пласта. Каждое численное исследование дополняется обсуждением и анализом полученных результатов. В трехмерной постановке строятся подробные неструктурированные сетки, чтобы достичь повышенной точности решения, что приводит к использованию для численных расчетов вычислительного кластера «Ариан Кузьмин». В конце работы приводятся основные результаты.
Постановка задачи
Течение однофазной жидкости в пористой среде описывается законом сохранения массы (уравнением неразрывности) и законом Дарси
д(фр)
dt
- + div (ри ) = 0, х ей, I е( 0, Т ], (1)
и = -kgrad(р + pgz) = 0, х еП, (2)
И
где ф - пористость среды, р - плотность, k - коэффициент проницаемости пористой среды, ^ - коэффициент вязкости фильтрующейся жидкости, р - давление, g - ускорение свободно -го падения и г - вектор, соответствующий вертикальной координате.
В данной работе мы рассматриваем процессы фильтрации несжимаемой жидкости в слабосжимаемом коллекторе, т. е. р=const и
дф дф др др
— = —— = сг—, (3)
дt др дt г дt у '
где сг - коэффициент сжимаемости коллектора.
Таким образом, для давления запишем следующее параболическое уравнение:
др (к *) сг—— divI — grad(р + р^) I = 0, х ей, г е (0,Т].
Уравнение (4) дополним граничными условиями второго рода (Неймана)
к д(р + pgz)
(4)
^ дп к д(р + pgz)
= 0, х еГл, (5)
= qi (), х еГ,, , = 1,2.....Nq (6)
р дп
и начальным условием
Р (х,0) = Ро (х), х еП, (7)
где п - внешняя нормаль к границе = Г0 и Г] и Г2 и.. .и ГК , Гд - внешняя граница области, исключающая все скважины, Г. - граница .-й скважины, к - количество скважин. Функция д() в граничном условии (6) - приток флюида, приходящийся на единицу поверхности ствола .-й скважины.
Аппроксимация по времени и пространству
Для аппроксимации по времени сначала построим равномерную сетку
ат = ат и{т} = {п = пт,п = 0,1,...,Nt,тN = т}
и используем обозначение р=р(Г).
При переходе с одного временного слоя на другой используем чисто неявную разностную схему по времени для уравнения (4)
п+1 _ п / ^ \
С р--_р— I _gШd (п+1 +pgz )1 = 0 (8)
с начальным условием р0=р0(х).
Для дискретизации по пространственным переменным используем метод конечных элементов. Для записи вариационной постановки задачи умножим уравнение (8) на тестовую функцию V и проинтегрируем по области &
п+1 _ п / — \
\ сг Р-уСх _ \ div I -grad (п+1 + pgz) I vdx = 0, Уу е V. (9)
Далее, используя формулу интегрирования по частям для оператора Лапласа с учетом граничных условий (5)-(6), получим
рп+1 _ (k , . , \ ^ - д(рп+ +pgz)
I агР-^ vdx +}1 - у(рп+1 +pgz), Уу I dx -XI - —-= 0, Уу е V. (10)
п т И у ' ) 1=1Г1М дп
В силу условий (6) на границах скважин получаем
п+1 _ п / — \ Ыу
| сг Р-vdx +}1 -у(рп+1 +pgz), Уу I dx = £/ Уу е К, (11)
п т И у ' ) ;=1гг
где F={v6Я1(&)}, Н\Ц) - пространство Соболева, состоящее из функций V таких, что V2 и VI2 имеют конечный интеграл в
В методе конечных элементов решение представленной задачи ищем в виде линейной комбинации заданных базисных функций
м
Р (х ) = £ Р№,
i=1
где ф. - стандартные линейные базисные функции, определенные в х (разбиение области О на конечные элементы) и N - количество узлов сетки.
Тогда уравнение (11) представим в следующей матричной форме:
M
pn+1 - pn
тг П+1 rn+1
KP =f ,
(12)
где
М = \ау ] = 1 К = \Ъу ] = I -Чф^.
п п И
Таким образом, на каждом временном шаге решаем систему линейных уравнений вида
(M + тК )pn+1 = tf"+1 + Mp
(13)
с матрицей размерностью
Численное моделирование в двумерной постановке
В данной главе рассмотрим численное решение задачи (4)-(7) в двумерной многосвязной области (рис. 1). Область является квадратной со сторонами по 4 км и радиусами скважин, равными 0,1 м.
Входные данные возьмем следующими: с=10~9 Па-1, &=10"13 м2, ,м=240-9 Па-с, «0=36 МПа, 7=1 год, ql=10 "4 м3/с, q=q/(N -1), .=2, ..., N, N=5.
т
or. 1 о
Гд
о о
Рис. 1. Двумерная область Q
Для исследования влияния шага по пространству на решение задачи решим задачу на разных пространственных сетках (табл. 1). Для этого построим четыре разные сетки и сетку с названием mesh 0 примем за эталонную, т. е. решение на этой сетке сравним с решениями на остальных сетках в каждый момент времени. Норму относительной погрешности вычислим в пространствах L2(Q) по формулам
:(t ) =
1(Р (x t)- Pe (x, t)) dx
J pe (x, t) dx
•100%, x eQ, t e (0, T] и H(Q),
:(t ) =
|(Vp (x, t)-Vpe (x, t)) dx
\Vpe (x, t) dx
•100%, x eQ, t e (0,T],
где p (x,t) - решение на эталонной сетке.
Q
Q
Таблица 1
Сетки с разными количествами узлов и элементов
Название сетки mesh 0 mesh 1 mesh 2 mesh 3
Количество узлов 186046 31827 7931 1087
Количество элементов 368540 62242 15150 2062
Рис. 2. Сетки из табл. 1
Из рис. 3 видно, что погрешность более четко выделяется в пространстве H(Q) и то, что сетка mesh 2 является оптимальной, т. к. дальнейшее сгущение приводит к незначительным уменьшениям погрешности с учетом времени, которое тратится на вычисления на ЭВМ. В дальнейших исследованиях все задачи в двумерной постановке решим на сетке mesh 2.
Для исследования шага по времени решим задачу со следующими разными временными шагами:
0.006 0.005 0.004
^ о.ооз 0.002 0.001 о
tau0=1 час, tau1=1 день, tau2=10 дней, tau3=1 месяц,
35
esh 1
- ; ____________mesh 2
/
/ - - -
- —^-
30 25 ^20 15 10 5
v— mesh 1
mesh 2
___________j___________{___________|____________i____________
4 6 8 t, месяцы
4 6 8 t, месяцы
Рис. 3. Относительная погрешность в пространстве Ь2(О) (слева) и Н(О) (справа) в зависимости от сетки
эталонный шаг возьмем tau0, чтобы сравнить ее с остальными. На рис. 4 видно, что шаг tau2 является оптимальным, как и в случае с пространственным сетками, поэтому в дальнейших исследованиях используем ее в качестве шага по времени.
Исследование итерационных методов
Для вычисления результирующей матрицы (13) и для того, чтобы выявить более эффективный решатель, применим разные итерационные методы подпространств Крылова с предобусловливанием и без него: метод обобщенных невязок (GMRES), метод сопряженных градиентов (CG) и стабилизированный метод бисопряженных градиентов (BICGSTAB). В качестве предобуславливателей используем неполную LU факторизацию (ILU) и метод релаксации (SOR).
Как следует из табл. 2, для данной задачи метод BICGSTAB с предобусловливателем ILU является наиболее эффективным. Из той же таблицы видно, что чем больше значение отношения параметров проницаемости пласта kjk (рис. 5), тем сильнее увеличивается количество итераций для всех решателей.
О 2 4 6 8 10 12 0 2 4 6 8 10 12
t месяцы 1, месяцы
Рис. 4. Относительная погрешность в пространстве ¿2(0) (слева) и Я(Щ (справа) в зависимости от шага по времени
_I №
■*>
Рис. 5. Метки подобластей
Таблица 2
Количество итераций в зависимости от решателей
kA GMRES CG BICGSTAB
ILU 1 369 126 90
10 459 177 117
100 500 186 127
1000 1012 207 149
SOR 1 381 132 95
10 487 192 121
100 537 194 138
1000 737 222 141
NONE 1 2637 289 224
10 7100 730 583
100 18211 1480 1090
1000 168207 3558 4385
к /к =100
к /к =1000
Рис. 6. Распределение поля давления и(х,Г): слева, когда =1 месяц, и справа, когда =1 год - 37
На рис. 6 представлено распределение поля давления в гомогенной и гетерогенных средах. Заметим, что скорость фильтрации (2) больше в подобластях с большей проницаемостью.
Численное моделирование в трехмерной постановке
Рассмотрим теперь задачу (3)-(7) в трехмерной области, где, как и в двумерной задаче, имеем пять скважин: одна вертикальная нагнетательная и четыре кривых горизонтальных добывающих скважин (рис. 7).
Входные данные возьмем из двумерной задачи с учетом гравитационных сил, где р=863 кг/м3, g=9,8 м/с2. Параметры области следующие: ширина и длина по 2571 м, высота области 250 м, длина центральной скважины 200 м, длина наклонных скважин 314 м, расстояние между центральной и наклонными скважинами 750 м, радиус скважин 0,1 м.
Область для достижения приемлемой точности разбивается на достаточно мелкие тетраэдры. Для исследований сгенерируем следующие две сетки: сетку mesh 1, имеющую 710 099 узлов и 3701704 элементов (рис. 8), и сетку mesh 2, имеющую 1407467 узлов и 7874635 элементов (рис. 9). Для параллельного вычисления варьируем количество процессоров, чтобы найти оптимальное их количество по времени вычисления для данной задачи и данных сеток.
Рис. 8. Трехмерная сетка mesh 1
Рис. 9. Трехмерная сетка mesh 2
В табл. 3 приведена зависимость времени решения от количества процессоров, где видим, что 64 процессоров являются оптимальными для обеих сеток, т. к. дальнейшее его увеличение приводит к увеличению времени счета за счет того, что тратится дополнительное время на обмен данными между процессорами. На рис. 10 изображено распределение давления для трехмерного случая на сетке mesh 2.
Таблица 3
Зависимость времени решения от количества процессоров
Количество процессоров 1 2 4 8 16 32 64 72 128
Время в минутах mesh 1 20,2 12 6,8 4,5 3,5 2,3 1,7 2,1 2,3
mesh 2 177 44,2 53 33 26 13,2 6 6,2 6,7
3.4е+07 3.71 е+07
Рис. 10. Распределение поля давления u(x,t) на сетке mesh 2: слева, когда t=1 месяц, и справа, когда t=1 год
Заключение
В данной работе рассмотрено численное решение задач однофазной нестационарной фильтрации в пористой среде в двумерной и трехмерной постановках. В двумерной задаче были исследованы итерационные методы, где метод BICGSTAB с предобуславливателем ILU показал наибольшую эффективность. Трехмерная задача смоделирована с учетом гравитационных сил на вычислительном кластере «Ариан Кузьмин». Приводятся результаты некоторых исследований и их анализ.
Л и т е р а т у р а
1. Logg A. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book / A. Logg, K. A. Mardal, G. Wells. - Springer, 2012.
2. Вабищевич П. Н., Васильев В. И., Васильева М. В., Никифоров Д. Я. Численное решение одной обратной задачи фильтрации // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2015. - Т. 157, кн. 4. - С. 79-89.
3. Азиз Х., Сеттари Э. Математическое моделирование пластовых систем. - М.: Недра, 1982. - 408 с.
4. Афанасьева Н. М., Васильева М. В., Колесов А. Е. Математическое моделирование фильтрации: учебное пособие. - Якутск: Издательско-полиграфический комплекс СВФУ, 2011. - 86 с.
5. Барренблатт Г. И., Ентов В. М., Рыжик В. М. Теория нестационарной фильтрации жидкости и газа. - М.: Недра, 1972. - 288 с.
6. Басниев К. С., Дмитриев Н. М., Розенберг Г. Д. Нефтегазовая гидромеханика. - М.-Ижевск: Институт компьютерных исследований, 2005.
7. Васильев В. И., Попов В. В., Тимофеева Т. С. Вычислительные методы в разработке месторождений нефти и газа. - Новосибирск: Изд-во СО РАН, 2000. - 126 с.
8. Лейбензон Л. С. Нефтепромысловая механика. Ч. 2. Подземная гидравлика воды, нефти и газа.
- М.-Грозный-Ленинград-Новосибирск: Горгеолнефтеиздат, 1934. - 352 с.
9. Леонтьев Н. Е. Основы теории фильтрации: Учебное пособие. - М.: Изд-во Центра прикладных исследований при механико-математическом факультете МГУ, 2009. - 88 с.
10. Нигматулин Р. И. Динамика многофазных сред. Ч. 2. - М.: Наука, 1987.
11. Чарный И. А. Подземная гидрогазодинамика. - М.: Государственное научно-техническое издательство нефтяной и горно-топливной литературы, 1963.
12. Cheng Z., Huan G., Ma Yu. Computational methods for multiphase flows in porous media // SIAM.
- Philadelphia, 2006.
13. Саад Ю. Итерационные методы для разреженных линейных систем. Том 1. - М.: Издательство Московского университета, 2013. - 344 с.
R e f e r e n c e s
1. Logg A. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book / A. Logg, K. A. Mardal, G. Wells. - Springer, 2012.
2. Vabishchevich P. N., Vasil'ev V. I., Vasil'eva M. V., Nikiforov D. Ia. Chislennoe reshenie odnoi obrat-noi zadachi fil'tratsii // Uchen. zap. Kazan. un-ta. Ser. Fiz.-matem. nauki. - 2015. - T. 157, kn. 4. - S. 79-89.
3. Aziz Kh., Settari E. Matematicheskoe modelirovanie plastovykh sistem. - M.: Nedra, 1982. - 408 s.
4. Afanas'eva N. M., Vasil'eva M. V., Kolesov A. E. Matematicheskoe modelirovanie fil'tratsii: uchebnoe posobie. - Iakutsk: Izdatel'sko-poligraficheskii kompleks SVFU, 2011. - 86 s.
5. Barrenblatt G. I., Entov V. M., Ryzhik V. M. Teoriia nestatsionarnoi fil'tratsii zhidkosti i gaza. - M.: Nedra, 1972. - 288 s.
6. Basniev K. S., Dmitriev N. M., Rozenberg G. D. Neftegazovaia gidromekhanika. - M.-Izhevsk: Institut komp'iuternykh issledovanii, 2005.
7. Vasil'ev V. I., Popov V. V., Timofeeva T. S. Vychislitel'nye metody v razrabotke mestorozhdenii nefti i gaza. - Novosibirsk: Izd-vo SO RAN, 2000. - 126 s.
8. Leibenzon L. S. Neftepromyslovaia mekhanika. Ch. 2. Podzemnaia gidravlika vody, nefti i gaza.
- M.-Groznyi-Leningrad-Novosibirsk: Gorgeolnefteizdat, 1934. - 352 s.
9. Leont'ev N. E. Osnovy teorii fil'tratsii: Uchebnoe posobie. - M.: Izd-vo Tsentra prikladnykh issledovanii pri mekhaniko-matematicheskom fakul'tete MGU, 2009. - 88 s.
10. Nigmatulin R. I. Dinamika mnogofaznykh sred. Ch. 2. - M.: Nauka, 1987.
11. Charnyi I. A. Podzemnaia gidrogazodinamika. - M.: Gosudarstvennoe nauchno-tekhnicheskoe izdatel'stvo neftianoi i gorno-toplivnoi literatury, 1963.
12. Cheng Z., Huan G., Ma Yu. Computational methods for multiphase flows in porous media // SIAM.
- Philadelphia, 2006.
13. Saad Iu. Iteratsionnye metody dlia razrezhennykh lineinykh sistem. Tom 1. - M.: Izdatel'stvo Moskovskogo universiteta, 2013. - 344 s.
Учебно-научная лаборатория «Современный якутский язык: перево-доведение и лингвостилистика» предлагает услуги профессионального перевода с русского на якутский, с якутского на русский язык и редактирования текста в публицистическом, официально-деловом, научном стилях речи. Цены устанавливаются в зависимости от сложности текста. Обращаться в 219 каб. ГУК СВФУ Телефон: +7 (4112) 496-750. E-mail: ksip219@mail.ru.
