Численное решение одной обратной задачи фильтрации Текст научной статьи по специальности «Математика»

____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 157, кн. 4 Физико-математические науки

2015

УДК 519.6

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОЙ ОБРАТНОЙ ЗАДАЧИ ФИЛЬТРАЦИИ

П.Н. Вабищевич, В.И. Васильев, М.В. Васильева, Д.Я. Никифоров

Аннотация

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

Ключевые слова: фильтрация, обратная задача, дебит, забойное давление, метод конечных элементов, неструктурированные сетки.

Введение

В теории и практике решения уравнений с частными производными большое внимание уделяется обратным задачам [1]. Из них выделим и рассмотрим обратную задачу по восстановлению зависимости правой части параболического уравнения от времени при известном распределении по пространству. Такая линейная обратная задача относится к классу корректных в классическом смысле задач математической физики. Дополнительная информация наиболее часто дается как значение функции в некоторых внутренних точках или как среднее значение, получаемое в результате интегрирования по области. Для численного решения таких обратных задач предложен и развивается вычислительный алгоритм, основанный на специальной декомпозиции задачи на каждом временном слое [2-5].

Разработка реальных месторождений углеводородного сырья приводит к необходимости решения задач в сложных расчетных областях. В настоящей работе задачи фильтрации изучаются в сложных многомерных областях с использованием достаточно подробных расчетных сеток с наличием сгущения в местах особенностей задачи. Поэтому для дискретизации по пространственным переменным целесообразно использовать метод конечных элементов. Так, в работе [6] предложен вычислительный алгоритм трехмерных задач фильтрации жидкости в напорнобезнапорном коллекторе со сгущающимися участками сетки в прискважинных подобластях.

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

79

80

П.Н. ВАБИЩЕВИЧ И ДР.

1. Постановка задачи

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

дтр

~дГ

+ div pu = f (x, t)

(1)

и законом Дарси

k

u = (grad p + pg). (2)

P

Здесь p - давление, u - скорость фильтрации флюида в пористой среде, к - проницаемость пористой среды, р - вязкость флюида, g - вектор ускорения свободного падения и f - плотность внутренних источников/стоков.

Подставляя закон Дарси (2) в уравнение неразрывности (1), учитывая слабую сжимаемость флюида, упругую деформируемость коллектора и пренебрегая действием гравитационных сил ввиду незначительной мощности нефтяных месторождений [7-9] (толщина пласта на 2-3 порядка меньше характерных размеров пласта в горизонтальной плоскости), получаем следующее уравнение для давления:

др [к

- div ( p gradp ) = f (x,t),

(3)

где в - коэффициент совместной упругоемкости флюида и коллектора.

В вычислительной практике решения задач фильтрации правую часть f (x,t) можно задать в виде суммы мощностей источников/стоков

др (к

в ----div ( — grad р

dt \p

N

J2qi(t)^i(x), x e Q, t e (0,T],

i=l

(4)

где qi для двумерного случая, когда задаем точечный источник/сток, - дебит i -й скважины (см. рис. 1, а), а в случае горизонтальной скважины (см. рис. 1, б и в) -приток флюида, приходящийся на единицу поверхности ствола i -й скважины, ^i (x) - неотрицательные весовые функции, Nq - количество скважин, T > 0 и Q e Ra, а = 2, 3.

Уравнение (4) дополняется соответствующими граничными и начальным условиями

-рй = '° x e Г- t e (0T (5)

p(x, 0) = po(x), x e Q, (6)

где Г - граница Q, n - внешняя нормаль к Г.

Таким образом, требуется найти функцию p(x,t), t e (0, T], T > 0, удовлетворяющую параболическому уравнению, граничным и начальным условиям (5), (6) при заданных входных данных к, р, c, po(x), qi(t), ^i(x), i = 1, 2,..., Nq. Начально-краевая задача (4)-(6) относится к классу прямых задач.

Рассмотрим теперь обратную задачу, где входящие в правую часть коэффициенты qi, i = 1, 2,..., Nq, уравнения (4) подлежат определению. Дополнительную информацию зададим следующим образом:

p(x,t)^i(x)= <fii(t), x e Si, i = 1, 2,..., Nq, t e (0,T]. (7)

Обратная задача состоит в определении функций p(x,t), qi(t), i = 1, 2,..., Nq, удовлетворяющих начально-краевой задаче (4)-(6) при дополнительном условии, что само решение p(x,t) на забоях скважин предполагается известным (7). Отметим, что поставленная задача является корректной, и в настоящей работе рассматривается вопрос о ее эффективной вычислительной реализации.

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОЙ ОБРАТНОЙ ЗАДАЧИ ФИЛЬТРАЦИИ 81

а)

б)

Рис. 1. Расчетные области О с различными расположениями скважин

2. Аппроксимация по времени и пространственным переменным

Введем равномерную сетку по времени

и т = ит [J {T} = {tn = пт, п = 0,1,...,NT, tNt = T}

и обозначим pn = p(tn). При использовании чисто неявной схемы по времени для уравнения (4) приближенное решение на каждом временном слое определяется как решение дифференциально-разностного уравнения [10]

pn+i _ pn / — \ Nq

в------------div [ — grad pn+M = q

т Vp J 7—1 '

,n+l

i

фг (X), X G Q.

(8)

Для численного решения поставленной задачи будем использовать дискретизацию по пространственным переменным методом конечных элементов. Запишем вариационную постановку прямой задачи (4)-(6): найти p(x,t) G V такую, что

\pn+l —рп , +

в--------v dx +

/(

— grad pn+1, grad v P

dx =

Г Nq

pq’n

о г—1

qn+1 фг (x)vdx

(9)

V v G H1 (Q),

где H1 (Q) - пространство Соболева, состоящее из функций v таких, что v2 и \Vv|2 имеют конечный интеграл в Q.

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

N

p(x) = ^2 pi Фг,

i—0

где фг - стандартные линейные базисные функции, определенные в Uh; Uh - разбиение области Q на конечные элементы - треугольники (рис. 2 и 3), тетраэдры (рис. 4), N - количество узлов сетки.

Тогда уравнение (9) можно представить в следующей матричной форме:

pn+1

Nq

M-

+ Apn+1 =£ qn+1 фг (X),

n

p

(10)

где

M = [mij] = J сфгфj dx, A

—

P

УфгVфj dx.

Таким образом, на каждом временном шаге мы будем решать систему линейных алгебраических уравнений вида (M + TA)pn+1 = тf n+1 — Mpn с квадратной

82

П.Н. ВАБИЩЕВИЧ И ДР.

матрицей размерностью N х N. Матрица M - диагональная с положительными элементами, A - разреженная, положительно определенная, симметричная матрица, обладающая свойством диагонального преобладания.

3. Вычислительный алгоритм решения обратной задачи

Рассмотрим приближенное решения обратной задачи (4)-(6), где дополнительное условие (7) в точке наблюдения приобретает вид

p(x*,t)= (fii(t), i = 1, 2,...,Nq, t e (0,T ], (11)

или в случае линии, моделирующей криволинейную (горизонтальную) скважину, задается в виде

<Pi(t) = J p(x,t) dx J dx, i = 1, 2,...,Nq, t e (0,T ], (12)

Si Si

где Si - линия, моделирующая i-ю скважину. Таким образом, <^i(t) - среднее значение давления на стволе горизонтальной скважины на каждом временном слое.

Для нахождения решения на новом временном слое pn+1 используем декомпозицию, предложенную в работе [5]:

N

pn+1(x,t) = yn+1(x,t) + ^ qn+ (t) wn+1(x). (13)

i=1

Подставим линейную форму (13) в уравнение (4) и получим уравнения для определения yn+1(x,t) и wrn+l(x).

Таким образом, для определения функции yn+1(x,t) получим следующее уравнение:

• yn+1 _ pn в-----------v dx +

к \

— gradyn+1, grad v ) dx = 0 Vv e H 1(0). M J

(14)

Для вычисления сеточных функций wi(x), i = 1, 2,..., Nq, решаем следующие уравнения:

в л + — wiv dx +

— grad wi, grad v M

dx

!Mx}vdx Vv e H,(n).

Q

(15)

Отметим, что функции не зависят от времени w'n+1(x) = wi(x). Для определения сеточных функций yn+1(x,t), t e шт, wi(x), i = 1, 2,. .. ,Nq, x e шh, приме-

ним метод обобщенных минимальных невязок (GMRES) с предобуславливателем ILU [11] (методы крыловских подпространств для систем линейных алгебраических уравнений с разреженными матрицами).

Из дискретного аналога дополнительного условия (11) определяем qrn+1:

n+1 = рП+1 - yn+1 (*l)

i wi(x*i)

1, 2,

Nq

i

(16)

В случае горизонтальных скважин (условие (12)) сначала найдем средние значения функций yn+1(x,t) и wi(x) на забоях скважин

пП+1 = yn+1 dx / dx, 6i = wi dx / dx,

(17)

Si

Si

Si Si

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОЙ ОБРАТНОЙ ЗАДАЧИ ФИЛЬТРАЦИИ 83

Рис. 2. Сетка расчетной двумерной области, содержащая 7602 вершин и 15042 треугольных элементов со сгущением в окрестности скважин

затем из дискретного аналога находим q,П+1 - объемный приток флюида, приходящийся на единицу поверхности ствола i-й скважины:

п+1

qn+1 —

qi

n+1

■Vi

i — 1, 2,

,Nq.

i

(18)

Для определения Qi(t) - дебита i-й криволинейной (горизонтальной) скважины - нужно вычислить следующий интеграл:

Qn +1 — J qn +1 dx, i — 1, 2,...,Nq. (19)

Si

Следует отметить, что для обеспечения корректности данного вычислительного алгоритма необходимо, чтобы 0i — 0 или wi(x*) — 0. Эти условия выполняются в силу неотрицательности финитных функций ^i (x).

4. Численные результаты

Рассмотрим численное решение обратной задачи с использованием предложенного алгоритма идентификации правой части параболического уравнения. Нами будут рассмотрены три обезразмеренные модельные задачи фильтрации слабосжимаемого флюида в упруго-деформируемой пористой среде:

84

П.Н. ВАБИЩЕВИЧ И ДР.

Рис. 3. Сетка расчетной двумерной области, содержащая 44189 вершин и 87872 треугольных элементов со сгущением в окрестности скважин

Рис. 4. Сетка расчетной трехмерной области, содержащая 170363 вершин и 981878 тетраэдральных элементов со сгущением в окрестности скважин

1) двумерная задача с одиночной нагнетательной скважиной в центре расчетной области и четырьмя добывающими скважинами, расположенными по бокам от добыващей скважины. Расчетная область Q = (0,1) х (0,1) (рис. 2);

2) двумерная задача с двумя криволинейными добывающими и одной вертикальной нагнетательной в центре области скважинами. Расчетная область Q = = (0,1.5) х (0,1) (рис. 3);

3) трехмерная задача с одной нагнетательной и четырьмя вертикальными добывающими скважинами. Расчетная область Q = (0, 3) х (0, 3) х (0,0.5) (рис. 4).

Расчеты проведены до момента времени T = 0.1 c временным шагом т = 0.001, с начальным условием po(x) = 1.0 и со следующими гидродинамическими характеристиками: в = 1.0, к = 1.0 и р = 1.0. В качестве граничных условий использовались нулевые однородные граничные условия второго рода (Неймана).

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОЙ ОБРАТНОЙ ЗАДАЧИ ФИЛЬТРАЦИИ 85

Рис. 5. Функции qi(t) и qi(t), i = 2, 3, . . . , при Nq = 5, справа при Nq = 3

Рис. 6. Дополнительные условия (забойные давления скважин) ^i(t) и ^i(t), i = = 2,..., 5. Тестовая задача 1

q, при разных количествах скважин: слева

Рис. 7. Абсолютные погрешности идентификации дебитов скважин. Тестовая задача 1

При проведении квазиреального эксперимента коэффициент qi (t) - дебит скважины (объемный приток флюида, приходящийся на единицу поверхности ствола i-й скважины) - задавался в виде следующей функции, зависящей от времени:

qi(t) =

20t

(1 + е7(t-0.5T) )’

qi (t) = -

qi(t)

Nq - 1 ’

i = 2, 3,..., Nq, t e (0,T]. (20)

Здесь при больших значениях параметра у функция qi (t), i = 1, 2,..., Nq, стремится к разрывной функции c точкой разрыва t = 0.5T. В настоящей работе у = = 2000.

Для демонстрации работоспособности предложенного алгоритма сначала решим прямую задачу, заключающуюся в определении зависимости забойных давлений уi (t) от времени (рис. 6, 8 и 10) по заданным функциям qi (t) (рис. 5). Полученные таким образом значения забойных давлений в дальнейшем используем для проверки точности предлагаемого алгоритма для численной реализации поставленной неклассической задачи. Во второй и третьей тестовых задачах значения забойных давлений находились посредством вычисления значений интегралов в уравнении (12) на каждом временном слое.

На рис. 7, 9 и 11 представлены абсолютные погрешности вычисления дебитов скважин (приток флюида) в каждый момент времени. Заметим, что они лежат в диапазоне 10-7^10-4. Приведенные значения абсолютной погрешности свидетельствуют о достаточно точном восстановлении искомых функций qi, i = 1, 2,... ,Nq.

Дебиты горизонтальных скважин находим из уравнения (19). На рис. 12 мы видим, что отношения между дебитами добывающих и нагнетательных скважин различаются не в два или четыре раза, как это следовало ожидать. Причинами тому являются длины добывающих горизонтальных скважин - они по площади

86

П.Н. ВАБИЩЕВИЧ И ДР.

Рис. 8. Дополнительные условия (забойные давления скважин) ф\ (t) и фi (t), i = = 2, 3. Тестовая задача 2

Рис. 9. Абсолютные погрешности идентификации дебитов скважин. Тестовая задача 2

Рис. 10. Дополнительные условия (забойные давления скважин) ф\ (t) и фi (t), i = = 2,, 5. Тестовая задача 3

Рис. 11. Абсолютные погрешности идентификации дебитов скважин. Тестовая задача 3

150

100

50

-50

-100

-150

1 1 1 1 ^1 1 1 1 1

V 4i

v q2 ...........

Чз

( .

1 1 1 1 J 1 1 1 1 -

0.02

0.04

0.06

0.08

а) Тестовая задача 2

20

-20

- 1 1 1 1 1 1 1 1 г qi 42 -

- Чз -

44

45

- 1 1 1 1 1 1 1 1 1 -

0.02 0.04 0.06 0.08 1

б) Тестовая задача 3

Рис. 12. Дебиты Qi горизонтальных скважин

поверхности превышают площади поверхности нагнетательных скважин, что влияет на вычисление значения интеграла в уравнении (19) по области скважины.

Далее, поставляя найденные коэффициенты qi (t), i = 1, 2,..., Nq , и сеточные функции y(x,t) и w(x) в каждый момент времени в уравнение (13), находим распределение давления p(x, t) по всей области. Распределения давления в разные моменты времени для первой модельной задачи представлены на рис. 13. На рис. 14 и 15 приведены графики с заливкой между линиями уровней в те же моменты времени для численных решений второй и третьей модельных задач соответственно.

Заключение

Предложен вычислительный алгоритм решения обратной задачи определения дебитов скважин по заданным забойным давлениям. Для дискретизации по пространственным переменным использован метод конечных элементов, позволяющий

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОЙ ОБРАТНОЙ ЗАДАЧИ ФИЛЬТРАЦИИ 87

t = 0.2T t = 0.49T

Рис. 13. Решение p(x, t) обратной задачи в разные моменты времени. Тестовая задача 1

p(x,t) p(x,t)

■6.79 - -6.79

t = 0.2T t = 0.49T

Рис. 14. Решение p(x, t) обратной задачи в разные моменты времени. Тестовая задача 2

t = 0.2T

t = 0.49T

Рис. 15. Решение p(x,t) обратной задачи в разные моменты времени. Тестовая задача 3

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

Работа выполнена при частичной финансовой поддержке Российского научного фонда (проект № 15-11-10024) и Российского фонда фундаментальных исследований (проект № 14-01-00785).

88

П.Н. ВАБИЩЕВИЧ И ДР.

Summary

P.N. Vabishchevich, V.I. Vasil’ev, M.V. Vasil’eva, D.Ya. Nikiforov. Numerical Solution of an Inverse Problem of Filtration.

A numerical method is suggested to solve the multidimensional inverse problem for the fluid flow in a porous media in order to determine the flow rates by the given bottomhole pressure. Approximation by space is based on the finite element method that allows simulation in unstructured grids with refinement near the location of wells. Discretization by time is performed with the use of implicit difference approximation. The numerical results for two-and three-dimensional test problems are presented.

Keywords: filtration, inverse problem, debit, bottomhole pressure, finite element method, unstructured grids.

Литература

1. Самарский А.Л., Вабищевич П.Н. Численные методы решения обратных задач математической физики. - М.: Изд-во ЛКИ, 2009. - 480 с.

2. Borukhov V.T., Vabishchevich P.N. Numerical solution of the inverse problem of reconstructing a distributed right-hand side of a parabolic equation // Comput. Phys. Commun. - 2000. - V. 126, No 1-2. - P. 32-36.

3. Вабищевич П.Н., Васильев В.И. Вычислительная идентификация младшего коэффициента параболического уравнения // Докл. РАН. - 2014. - Т. 455, № 3. -C. 258-260.

4. Vabishchevich P.N., Vasil’ev V.I. Computational algorithms for solving the coefficient inverse problem for parabolic equations // Inverse. Probl. Sci. Eng. - 2014. - doi: 10.1080/17415977.2014.993984.

5. Вабищевич П.Н., Васильева М.В., Васильев В.И. Вычислительная идентификация правой части параболического уравнения // Журн. вычисл. матем. и матем. физики. - 2015. - Т. 55, № 6. - С. 1020-1027.

6. Морозов П.Е., Хайруллин М.Х., Шамсиев М.Н. Численное решение прямой и обратной задачи при фильтрации флюида к горизонтальной скважине // Вычисл. методы и программирование. - 2005. - Т. 6. - С. 262-268.

7. Щелкачев В.Н. Упругий режим пластовых водонапорных систем. - М.: Гостоптех-издат, 1948. - 144 c.

8. Булыгин В.Я. Гидромеханика нефтяного пласта. - М.: Недра, 1974. - 232 с.

9. Васильев В.И., Попов В.В., Тимофеева Т.С. Вычислительные методы в разработке месторождений нефти и газа. - Новосибирск: Изд-во СО РАН, 2000. - 126 c.

10. Самарский Л.Л. Теория разностных схем. - М.: Наука, 1977. - 656 с.

11. Саад Ю. Итерационные методы для разреженных линейных систем. Т. 1. - М.: Изд-во Моск. ун-та. 2013. - 344 с.

Поступила в редакцию 08.09.15

Вабищевич Петр Николаевич - доктор физико-математических наук, профессор, заведующий лабораторией, Институт проблем безопасного развития атомной энергетики РАН, г. Москва, Россия.

E-mail: [email protected]

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОЙ ОБРАТНОЙ ЗАДАЧИ ФИЛЬТРАЦИИ 89

Васильев Василий Иванович - доктор физико-математических наук, профессор, заведующий кафедрой, Северо-Восточный федеральный университет им. М.К. Аммосова, г. Якутск, Россия.

E-mail: [email protected]

Васильева Мария Васильевна - кандидат физико-математических наук, старший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН, г. Новосибирск, Россия.

E-mail: [email protected]

Никифоров Дьулустан Яковлевич - магистрант, Северо-Восточный федеральный университет им. М.К. Аммосова, г. Якутск, Россия.

E-mail: [email protected]