Усовершенствованный попеременно-треугольный метод численного решения пространственно-трехмерных задач фильтрации двухфазной несжимаемой жидкости Текст научной статьи по специальности «Математика»

МЕХАНИКА MECHANICS

УДК 519.6 DOI 10.12737/ 18131

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

А. И. Сухинов1, Л. А. Григорян2, А. А. Сухинов3

1Донской государственный технический университет, г. Ростов-на-Дону, Российская Федерация 2Северо-Кавказский федеральный университет, г. Ставрополь, Российская Федерация 3Южный федеральный университет, г. Ростов-на-Дону, Российская Федерация

Advanced alternate triangular method of numerical solution of spatial three-dimensional problems of two-phase filtration of incompressible fluid ***

A. I. Sukhinov1, L. А. Grigoryan2, A. A. Sukhinov3

'Don State Technical University, Rostov-on-Don, Russian Federation 2North Caucasus Federal University, Stavropol, Russian Federation 3Southern Federal University, Rostov-on-Don, Russian Federation

Целью данной работы является рассмотрение пространственно-трехмерной модели Баклея — Леверетта фильтрации двухфазной несжимаемой жидкости «нефть — вода» в приближении Буссинеска и ее численная реализация усовершенствованным вариантом модифицированного попеременно-треугольного метода (МПТМ), учитывающим специфику функций точечных источников, имеющих дельтаобразный характер. Для предложенного усовершенствованного варианта МПТМ выполнено построение набора итерационных параметров и получены необходимые оценки спектральных параметров метода. Оценка скорости сходимости итерационного процесса асимптотически эквивалентна оценке для базового итерационного МПТМ, предложенного ранее А. А. Самарским и Е. С. Николаевым. Построенный вариант итерационного метода был применен для решения модельной задачи фильтрации — вытеснения нефти методом заводнения в существенно неоднородном пласте, для которого значения коэффициента проницаемости менялись более чем в 104 раза. Что касается вычислительной эффективности такого подхода, то число итераций сократилось в 10-70 раз по сравнению с традиционно используемыми итерационными методами — Зейделя и верхней релаксации с шахматной упорядоченностью узлов.

The work objective is to consider a space 3D Buckley -Leverett filtration model of the two-phase incompressible liquid "oil-water" to Boussinesq approximation, and its numerical solution by the modified alternate triangular method (MATM) with account for the specificity of the point source functions with deltalike character. For the proposed advanced MATM version, a set of iteration parameters is built, and the needed estimates of the method spectral parameters are obtained. The estimation of the convergence rate of the iteration process is asymptotically equivalent to the assessment for the base iterative MATM proposed earlier by A.A. Samarsky and E.S. Nikolaev. The constructed iteration method modification is applied to solve a model task of filtering-displacement of oil by the flooding method in the essentially heterogeneous stratum for which permeability coefficient values have been changed by more than 104 times. As for the computational efficiency of this approach, the number of iterations is reduced by 10-70 times as compared with the widely used techniques - Seidel method, and the overrelaxation with checkered nodes ordering.

Ключевые слова: модели фильтрации Баклея — Леверетта, двухфазная несжимаемая жидкость, разностные схемы, итерационный попеременно-треугольный метод, итерационные параметры.

Keywords: Buckley - Leverett filtration models, two-phase incompressible liquid, difference schemes, iteration alternate triangular method, iteration parameters.

Работа выполнена при финансовой поддержке РФФИ (грант 15-01-08619-а), а также в рамках Программы Президиума РАН «Фундаментальные проблемы математического моделирования».

** E-mail: sukhinov@gmail.com, honey.lusine@mail.ru, andreysukhinov@gmail.com

*** The research is supported in part by RFFI (grant 15-01-08619-а) and is done within the frame of the RAS Presidium Program "Fundamental research challenges of mathematical modeling".

eö «

S X

eö

*

(U

Введение. В современных условиях для решения значимых задач проектирования разработок нефтяных месторождений (РНМ) и совершенствования технологии добычи углеводородного сырья используется методология математического моделирования и комплексы программ, реализующих современные вычислительные методы [1, 2]. Важный класс задач, возникающих при проектировании РНМ, — задачи фильтрации многофазных жидкостей в пористых средах, сформулированные при тех или иных упрощающих предположениях. Данный класс задач после дискретизации приводит к сеточным уравнениям, которые могут иметь размерности 10 5-109 и должны решаться многократно (многие десятки, сотни тысяч раз) для временных промежутков, составляющих десятки месяцев. В практике вычислительных методов решения задач данного класса широко применяется IMPES-метод (IMplicit on Pressure and Explicit on Saturation) — неявный по давлению и явный по насыщенности [3]. Поэтому основной объем вычислительной работы приходится на решение уравнений эллиптического типа в относительно небольшом числе узлов сетки с функциями источников, отличными от нуля и имеющими характер дельта-функций (функций точечных источников) для определения давления. Следует учитывать данную специфику сеточного оператора задачи. В противном случае, а также при недостаточной обусловленности системы разностных уравнений, вызванной высокой размерностью и существенной неоднородностью характеристик пластовой системы, снижается качество вычислений. В первую очередь следует отметить, что замедляются итерационные процессы, а при использовании техники чебышевского ускорения наблюдается потеря сходимости метода, обусловленная ошибкой в определении спектральных характеристик операторов. В данной работе указанная проблема снимается за счет построения оператора-предобусловливателя, содержащего член, обусловленный наличием источников и корректных спектральных двухсторонних оценок для оператора-предобусловливателя.

1. Непрерывная пространственно-трехмерная модель двухфазной фильтрации. Нефтеводоносный пласт может существенным образом менять глубину залегания и (или) иметь существенную толщину и неоднородность по вертикали. В таких случаях требуется использование пространственно-трехмерной модели, учитывающей влияние силы тяжести, толщины пласта и неоднородности его параметров по вертикали. Будем использовать Декартову систему координат — ось Ox направлена на север, ось Oy — на восток, Oz — вертикально вниз. Процесс фильтрации рассматривается в пространственно-трехмерной области G — цилиндрической области с боковой поверхностью £(x, y), (x, y)еГ, Г = D \ D, ограниченной снизу поверхностью подошвы пласта, задаваемой функцией поверхности

z = Hb (x, y), (x, y) e D, D e R2, а сверху — поверхностью кровли пласта z = Hr (x, y), (x, y) e D.

В отношении граничных поверхностей области G предполагается их достаточная гладкость, необходимая для корректной постановки задачи. Заметим, что в этом случае граничная поверхность Г области G состоит из объединения достаточно гладких поверхностей Г = Hb UHr U G = GU Г.

Для проницаемой кровли пласта уравнения, описывающие процесс фильтрации в переменных «давление — водонасыщенность», в случае несжимаемости сред имеют вид:

д dx

irf1(s)kh t f2(s)kh )др \ d ff fi(s)kh t fi( s)kh W

Mi M2 )dx) dy 111 M1 Mi )dy

+

+1 Вт d+-)+M (§+-11+«1=(i) md^= д_ffi(s)kh др\ д ff2(s)kh dp

dt dx ^ m dx) dy ^ m dy

+ 3( /^fdP+^lU. (2)

дк ^ /л2 ^ дк

и В уравнениях (1), (2): я = я(х, у, г, /) — водонасыщенность; /¡(я), /2(я) — относительные фазовые проницае-

й

-§ мости для нефти и воды соответственно; т — пористость пласта; /и1, /и2 — динамические вязкости нефти и воды со-

■М

•д ответственно; р1, р2 — плотности нефти и воды соответственно, а рь р2 — их фазовые давления соответственно;

1л

«и кн, ку — коэффициенты абсолютной проницаемости соответственно в горизонтальном и вертикальном направлениях.

Здесь мы отметим следующее: как правило, кк значительно больше ку и ку существенно изменяется в зависимости £ от координаты г.

Будем считать, что на скважинах задаются либо дебиты, либо забойные давления. В качестве функциональных зависимостей для задания / (я), I = 1,2 будем использовать полиномы второго порядка относительно переменной я.

Граничные условия будем рассматривать в связи с тем, является ли граница Г области G проницаемой. На

dp

непроницаемой границе поток по нормали должен быть равен нулю, что приводит к условию — = 0 . На проницаемой

dn

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

wn = Wj + w2,

f (s)

где w1, w2 — потоки нефти и воды соответственно, удовлетворяющие условиям: wi = -к —'— grad p.

Mi

Граничные условия 2-го рода реализуются, когда задан поток, закачиваемый (отбираемый) на части цилиндрической границы Sj:

wn(x,y,z, t) = W(x,y,z, t), (x,y, z) e Sj, 0 < t < t0, s(x,y,z,t) = S(x,y,z,t), (x,y,z) e Sj, 0 < t < t0, а также для отбора, заданного на части цилиндрической границы S2:

wn (x, y, t) = W (x, y, t), (x, y, z) eS2, 0 < t < t0,

или давления

wn (x, y, t) = W (x, y, t), (x, y, z) eS2, 0 < t < t0, p(x,y,t) = P(x,y,t), (x,y,z) eS2, 0 < t < t0. Разумеется, на части цилиндрической поверхности S3 также могут быть заданы условия непроницаемости.

2. Дискретизация пространственно-трехмерной модели фильтрации на основе метода неявного по давлению и явного по насыщенности. Ранее МПТМ для двумерных задач фильтрации с огрубленными спектральными оценками рассматривался в работе [4], а для уравнений параболического типа — в монографии [5]. МПТМ для случая несамосопряженного оператора сеточной задачи при ограниченном значении числа Пекле, предназначенный для решения задач морской гидродинамики, был построен и исследован в работах [6, 7]. Далее рассматривается как наиболее распространенный алгоритм дискретизации задачи, базирующийся на IMPES — методе, неявном по давлению и явном по насыщенности. Также будем считать, что G( x, y, z) — область, в которой требуется найти решение задачи (1)-(2), — параллелепипед. В области G используем равномерную пространственную сетку ah :

ah ={ x = ihx, yj = jhy, zk = khz ,0 < i < N ,0 < j < Ny ,0 < к < Nz} и неравномерную временную сетку al.

Аналогично предыдущему определим пространственно-временную сетку a)hT = ahxa>T.

Разностные аппроксимации уравнений (1) и (2) во внутренних узлах сетки ah имеют вид для уравнения (2):

2 (Pi- 1,j,k Pi, j,k ) B., ,2 (Pi, j,k Pi+1,j,k ) B4., +

i—,j,k h ' ' ' i+-,j,k 2 x 2

+ 77 (Pi, j-1, к Pi, j, к ) B ( ) l ~2 (Pi, j+1, к Pi, j, к ) B ( ) l + h;V 7 ' ' i, J--, к J ' i, j+ -, к

+ ТГ(Pi,j,k-1 -Pi,j,k)B(3\ \-~T(Pi,j,k-Pi,j,k+1)B У 1 +

... L 1 Z \ ^Ji^ 1 * / . . , 1

i, j,k— h i, j ,k+-

2 z 2

+ (0,5 ( f1i, j,k-1 + f1i, j,k ) <vi, j ,k-0,5 -0,5 ( f1i, j,k+1 + f1i, j,k ) <vi, j,k+0,5 ) +

M1hz

+ g^2 (0,5 ( f2 i, j ,к-1 + f2 i, j,k ) <vi, j ,k-0,5 -0,5 ( f2 i, j,k+1 + f2 i, j,k ) <vi, j,k+0,5 ) = М2 hz

0 - если узел вне скважины,

-B . ., IP . ., - P. . , - если узел находится на нагнетательной скважине; Bgi,j,k[Po i,j,k Pi,j,kk)-если узел находится на эксплуатационной скважине. ( )

ев

и —

X ев

X

К

Следует отметить, что значения функций относительной фазовой проницаемости в уравнении (3) и в уравнении (4) в зависимости от значений функции водонасыщенности вычисляются на предыдущем временном слое (метод явный по функции водонасыщенности):

тН

+1 _ 1

л1,},к л1,] ,к 1

"(31 /2 ¡_1,],ккМ_1,],к + 31 /2 ¡, ],ккМ, ],к )(р _1, ],к р, ],к ) +

+ 2 (32 /2и],ккЫ,},к + 32 /21+1,j,kkйi+1,},к )'(р,},к р+1,],к )

М-2 "х

72 (З3 /2 i,]_1,ккМ,]_1,к + З3+/2 i,],ккЫ,],к )(р,]_1,к _р,],к ) +

М-2 "у

+ 72 (34 /2 i,],ккЫ,],к + 34 /2 i,] +1,ккЫ,]+1,к )(р,],к _ р,]+1,к )_

М-2 "у

72 (35 /2 i, ],к_1куг', ] ,к _1 + 35+/2 и ],кКи ],к )(р, ],к _1 _ р, ],к ) +

М-2 "г

+ 72 (Зб /2 i, ] ,kkvi, ],к + Зб /2 i, ],к+1куг', ],к+1 )(р, ],к _р, ],к+1 )_

М-2""г

(35 /2 i, ] ,k_1kvi, ],к _1 + 35+ /2 i, ],кkvг, ],к ) + (Зб /2 и j,k+1kvi, ],к+1 + Зб /2 i, ],кkvi, ],к ) :

^2 "г М-2 "г

0 _если узел вне скважины;

В . . , I Р . . , _Р. . , - если узел находится на нагнетательной скважине; wi,J,k i,J,k) '

В . . ,(р . . ,_Р. )с

(4)

gi,J,k[ о i,J,к i,J,k) i,J,k _ если узел находится на эксплуатационной скважине.

Забойные давления на нагнетательной и эксплуатационной скважинах соответственно:

Р Р '

w /,],к' о /,],к'

Bgг,],к -Щг,],к"г \ -+-

. 2п

, Л-~к-■

i,j,k 1п--1,51п2 _ с

Г

(5)

Символы

31+ , 32 , З3+ ,3+,35+ , Зб 31 , 32 , З3 ,34 ,35 , Зб

<5+ =(1-<5/"Х (>. =1,2,3,4,5,6

3 -

3 -

3 -

1- если р_1] к _ р ],к > 0 - в противном случае;

1- если р ]_1,к _ р,],к > ° 0 - в противном случае;

1- если р ] к_1 _ р ]к > 0, 0 - в противном случае;

определяются из условий:

з -

з -

з -

1- если р ],к _ р+1,],к >0,

0 - в противном случае;

I- если рр, ],к _ р, ] +1,к > 0 0 - в противном случае;

1- если р],к _ р,] к+1 > 0

0 - в противном случае.

(б)

с о тз

"¡3

-м

ел

Л

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

В

(1)

' ± ],к 2

^ 1

сХ

В

(2)

¿, ] ±—,к 2

"х Х к" (x, У] )(^1/2 ( s(x, У] ^ )) + М2 /l(x, У] ,2к ))

1 у_М1М2 сУ_

"у У1 К (Х¿, У, гк ) (М'1/2 (Х¿, У, гк ) + И2 М^ , У, гк ) )

В

■(3)

],к ±-

1 2к±1 1 1

К ч К (х1, у}, г) ( ь<х, , У],г) ) + ^2 /1(хг , Уу ,г) )

(7)

(8) (9)

3. Усовершенствование модифицированного попеременно-треугольного метода решения пространственно-трехмерной разностной задачи двухфазной фильтрации. Будем для простоты рассматривать граничные условия 1-го рода. Систему разностных уравнений (6) запишем в стандартном виде [8]:

-2 (°а Уха )_ + Ч (х)У (х) = ^(х) , х 6®,

а=1 ха

у (х) = ц(х), х 6 у, у 6 Г , 0 < с1 < аа (х) < с2 ч (х) ^ 0, Ух 6 ® . Здесь Г — граница прямоугольного параллелепипеда О = {0 < ха < !а, а = 1,2,3}, ю = ю + у — равномерная сетка, покрывающая область О, у — множество граничных узлов сетки.

Связь между коэффициентами уравнений (8) и (3) задается очевидным образом

а1 = а1 (хИ + К1, х2 ] , х3 к )- В ( \ . , а1 = а1 (хИ , х2 ] , х3 к )- В () ,

г +—, ],к 2

г—, ],к 2

= а2 (х1г, х2] + К2 , х3к ) - В ^ 1 . а2 = а2 (х1г, х2] , х3 к ) - В( .) 1

г, ] +— ,к 2

г,1Ч— ,к 2

= а2 (х1г, х2 ] , х3 к + К3 ) - В ( ) 1 . а3 = а3 (х1г, х2 ] , х3 к ) - В ^ ^

г, ],к +— г, 1,к +—

2 2

х1 6®, х2 6®2, х3 6®3,

®а={ха= гаКа , 1 < 'а < ^ - О, а= 1,2,3, ® = ®1 х ®2 х®.

Предполагается, что имеет место связь между искомыми сеточными функциями уравнений (10) и (11) и сеточными функциями правых частей этих же уравнений:

У (х1г, х2] , х3к ) - Рг,],к , Ч (х1г, х2] , х3 к ) - Bg г,],к ,

Я>( хг, х2 ], хз к )-- Ц" (0,5 (/\]Л-1 + Л] ) К

%г,],к-0,5 '

"°,5 (/ 1г,],к+1 + / 1г,],к )Кг,],к+0,5 ) ^^ (0,5 (/ 2 г,],к-1 + / 2 г,],к )к1 '(/ 2 г, ],к+1 + / 2 г, ],к ) Кг, ],к+0,5 )--^г, ], кр

Уг, ],к-0,5

-0,5 (

ТТ И и и

Здесь значения относительных фазовых проницаемостей: / ц у к_1, / ц у к, / 2 г ] к-1' / 2 г ] к входят в соответствии с равенствами (10) в функции правых частей сеточного уравнения (8) и рассчитываются на основе значений сеточной функции водонасыщенности на предыдущем п-м временном слое. При этом

Р . ., - если узел находится на нагнетательной скважине; ™ г,],к 7

, Ро г у к~ если узел находится на эксплуатационной скважине;

Р

gi,],k

0 - если узел не совпадает ни с одной из скважин. Будем использовать смещенные сетки

®+ ={х = г К, 1 < \ < }, ={х2 = К., 1 < г2 < N2 }, ={х3 = г3К3,1< г3 < N3}

и скалярные произведения

(и,у) = 2 м (х)V(х)К, (и,V) + = 2 м (х)V(х)К ,

(м,у)® = 2 м (х)у (х)К2, (м,у)®+ = 2 и (х)у (х)К2,

х26®2 z х26®2+

(м,у)® = 2 и (х)у (х)^3, (м,у)®+ = 2 и (х)у (х)К

х26®3 3 х2 6®3

И К

в

к

(и,V) = Е Е Е и(х)х)Н1И2И3,

х1е®+ х2е®2 х3е®3

(и,= Е Е Е и(х)у(х)Н1И2ИЪ,(и,у)3 = Е Е Е и(х)у(х)Н1Н2И3.

х1е® х2е®2 х3е®3 х1е® х2е®2 х3е®3

Представим систему разностных уравнений (8), (9) в виде операторного уравнения с однородными граничными условиями, изменив соответствующим образом правую часть:

Ay = u,

Ay = -T\aa yx | + qy, xe®,

(10)

7 (х) = 0, х еу. (11)

Представим схему итерационного двухслойного модифицированного попеременно-треугольного метода в операторном виде:

п+1 - п

(D + аЯг) D- (D + юК2 )^-?— + Ayn = f,

vn+1

где операторы, входящие в уравнение (12), запишутся

(

R y = Z

а=1

Л

аа aaxa 1

yx +—- y+-qy

h - 2 h 6

V - ^"a u у

( „+1

R2 y = -Z

a=1

a- a-xa 1

v ha a 2 ha 6 у

(12)

(13)

(14)

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

к = - к;.

Очевидно, что

А = К + К2.

Получим необходимые оценки констант, входящих в операторные неравенства [8]:

Ю < К + К2, К1 Б^Щ <^(к1 + К2 ) .

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

RD~% y, y| =

f , ( „+1

o o } 13

-rZ

d a=

V a

f

a- o aax- o 1 o

yx +—^^ qy К a 2h„ 6

3

Z

a=1

(„+1

w

a- o aaxa o 1 yx +—-y— qy

h a 2h 6

V a u у у

+1

3

Z

V-=1 V

a - o a—x- o 1 yx +—- y— qy

h„ a 2h„ 6

2 Л

у у

Будем использовать лемму: «Для любых чиселpa, ua, va, ra имеет место следующее неравенство:

п 2 п п IГ I

Z[ P-ua+ vara ] ^Z Z S-ß(\ P- +Xa\ ra\) |p- u- + v- ~

a=1 ß=1

X

- у

с о тз

"¡3

-м

ХЛ (U

£ л

где £aß > 0,£aß = —, £aa = 1, > 0».

Доказательство этой леммы основывается на е-неравенстве и здесь не приводится [7]. Положим в неравенстве (17)

P- =-

1 0

u-= yx-, v-= - =

tl,v

^Ч _М 2 6 '

"-ß

( x ) =

aßß hß + Xß aßxß 2hß 6

a- + X- a-x- q

h- 2h- 6

eß(x-, xs,) xß, xs2)'

где ß = 1,2,3, - = 1,2,3,

(15)

(16)

(17)

0 (х2 , Х3 ) , (х2 , Х3

)е®2 х®з' если Р =

вр(^ха, х^) = -в2 (х1; х3), (х1; х3) е юх хю3, если р = 2;

в3 (х1, х2 ), (х1, х2 )е®1 х®2, если ¡3 = 3;

01 (х2 , х3 ) , (х2 , х3 )е®2 х«3, если а = 1; ва (хр, х8^ = ■ в2 (х1, х3), (х1, х3) е а>1 х«3, если а = 2;

в3 (х1, х2 ), (х1, х2

)е«1 х®2, если а = 3. Функцию ^(х), соответствующую оператору Б, определим так:

й (х) = 2

а=1

а+а+ьа Ха

аха Ч 6

2 к

а

Из соотношений (16) и (17) получаем неравенство

Б-% у, у |< 2

а=1

о 2 ^

- у^ ,1

+ 2

а=1

Заметим, что справедливо неравенство

^+1 о 2 ^ (

— ух ,1 в

V а У V

1 ааха

ва Ха 2 ка

>

,1

У а

х ею.

о 2 ^

у ,1

У

Из этого и соотношения (18) следует:

Я,Б-1К2 у, у |< 2

а=1

(а о2 ^ ^ ,1

У а

1 ааха

в Х аЛ-а 2Ьа

о 2 1 у ,1

Будем считать, что уа (х) — решение краевой задачи

(а \ Ч а

ад -3^

"аха Я

2К 6

V = 0, ха = ° ха = 1а.

а = 1,2,3, х„ е ю

Положим

са(хр,х3)= тах Vх (х), х = (х1,х2,х3)ею,

где а = 1,2,3.

са (хр, хз) =

с1 (х2, х3 ), (х2, х3 ) е ю2 хю3, если а = 1; с2 (х1, х3 ), (х1, х3 ) е ю1 хю3, если а = 2; с3 (х1, х2 ), (х1, х2 )ею1 хю2, если а = 3.

(18)

(19)

(20)

(21)

(22)

(23)

Используя лемму 14, приведенную в п. 4, § 2, гл. V [8], получаем:

ааха Ч

2И„ 6

о 2 А

у ,1

< с„

^ о 2 о 2А ^

о а о аа у~ха+ 3 у

1

У У

Умножая обе части последнего равенства на произведение сеточных функций ваха и суммируя скалярно по Декартову произведению сеток юрхю5,

где

а = 1,2,3;

(Р, 8) =

(2,3), если а = 1; (1,3), если а= 2; (1,2), если а = 3,

получаем

ваХа

аха Ч

2И„ 6

о 2 ^ (

у ,1

ваХаса

аа ух.. + Т у

2

1

У У а

(24)

Рассмотрим теперь другую совокупность трехточечных разностных задач:

и к

X а

X к

11

Пусть

(а \ Ч а аа а г\ п 1 л л о

аа™ха) ~ м = -¡Г' М = ° ха = ° ха = -, а = 1,2,3;

а 'х„ 3 /И

'х- 3 ка

ха еа-' х = ( х1' х2' х3 )ею

Ьа(хр,х8^) = тах м- (х), х = (х1,х2,х3) е ю, а = 1,2,3,

(25)

(26)

где

Ьа (хр, х$) =

Ь1 ( х2, х3 ), (х2, х3 )ею2 хю3, если а = 1; Ь2 (х1, хз ), (х1,хз)е®1 х®з, если а = 2; Ь3 (х1, х2 ), (х1, х2 )ею1 хю2, если а= 3.

Базируясь на уже упоминавшейся лемме [8], из соотношений (25) и (26) получаем оценку вида

<аЛ '2 ^

а у•'

V а

< Ь

(( о 2 ^ о „о

аа У-а + 3 у

а = 1,2,3.

у

Умножим обе части последнего равенства на сеточную функцию 9а, суммируя скалярно по Декартову произведению сеток юр х щ8, где

а = 1,2,3;

(Р, 8) =

(2,3), если а = 1; (1,3), если а = 2; (1,2), если а = 3,

получим неравенство вида

( а+1 о 2 Л ( ( У ,1 <

а к2 V а

2 о Л Л

о а о аа у- + 3 у V V у У

ОаЬа

(27)

Сложим (24) и (27) и с учетом формулы (19) получим

( (

( „2 Л

йу ,1

< 2

а=1

(Ьа+Уаса)ва

\ \

2 а °

аа Уха + 3 У

(28)

й о чз

"¡3

и

и

О, С Л

Построим функцию

Здесь

а = 1,2,3, х еюр, х8 е ю8;

(хр>хз)= , 4 ' Ь + У с

иа ^ Л—а

(Р, 8) =

(2,3), если а = 1; (1,3), если а = 2; (1,2), если а = 3.

Тогда неравенство (28) запишется в виде

(о2 Л Г о о Л

йу , 1 <1 Ау,у I.

V У V У

Следовательно, 8 = 1. Рассуждая подобным способом, придем к оценке следующего вида

1 ааха

®ауа 2к-

о 2 ^

у ,1

^ауа

аа у* у

Л Л 1

у у а

(29)

(3°)

Неравенство (21) с учетом оценки (27) запишем

R D-1К2 У, "у 1< Z

а=1

Га o 2 ^ ^ ,1

( (л

/а

аа Ух, +Т У

^

1

J У

(31)

Предположим, что

q = max

X,,Х-.,х,

[q (X1, х2, хз)}

Исходя из соотношения (32), равенство (31) можно усилить и записать так:

RD-lR2 у, У 1< Z

а=1

Га o2 ^ (

агУж. ,1

У а

( „2 _„2Л Л

\°аХа \

аа Ух„ + Т У

1

J J

Применяя лемму 17, п. 5, § 2 гл. V [8], получаем неравенство:

( „2 Л /„2 Л

аа уХ„ ,1

\ Ут,

< Та

У ,1

\ Jm„

где

Га = Уа(хр>Xs);а = 1,2Д хр х,

(А 8) =

8 ш8>

(2,3), если а = 1; (1,3), если а = 2; (1,2), если а = 3;

У1 = —max— (0,х2,х3),а1 (lx,х2,х3), max

4 I

У2 =~7max | а2 (X1,0, х2 ), а1 (^ ^ х3 )

0<х <1 -h

max

0< х2 <l2-h2

4 -

у3 = — max < а3 (х1, х2, 0), а3 (х1, х2, l3), max

h, j 0<х3 <l3-h3 Ч

Опираясь на неравенство (31), имеем цепочку неравенств:

( „2 Л

а1 (х1, х2, х3

) + а1 (х1 + h1, х2, х3 )

0*2 (х1, х2, х^ ) + а2 (х1, х2 + h2, х^ )

> (х1, х2, х3 ) + а3 (х1, х2, х3 + h )

ааУх„ ,1 \ J ш

< Га

(„2 \ А

= Га ~

q

у ,1

\ J ш

fv ,1

а = 1,2,3.

ш

Будем задавать постоянные величины в соответствии с равенствами:

k = Уа На - 1 '

У а + 3 q

а = 1,2,3.

Базируясь на соотношениях (35), (36), приходим к цепочке неравенств:

Г „2 Л ( „2 Л ( „2 Л

аа Уха,1

= к1а

аа Ух,1

а s ха

\ уш

+ (1 - к1а)

аа Ух,1

< к1а

= К

( o 2 ^

аа ^ха ,1

( o 2 ^

аа ^ха ,1

+ (1 - к1а)-Уа

q

а ха

\ уша

2

qУ ,1 3

+ к1а

-о 2 ^

q Л

—у ,1

3

У ш„

а = 1,2,3.

У ш

Рассуждая аналогичным образом и учитывая неравенство (34), приходим к соотношению

fa о2 ^ (

ат.Ух. ,1

( „2 _„2Л Л

%1а

аа Ух„ + ТУ

1

J J

а = 1,2,3.

(32)

(33)

(34)

(35)

(36)

(37)

ей И К X а X к

(38)

Воспользуемся для преобразования неравенства (38) соотношением (33):

КБ-1Я2 "у,у|< £

а=1

( 2 \ ( ( 2

п г о 77 о

аа

Уа \в<хХа V

аа ух,. + ~у

3

< £

а=1

(1. ( 1а

Ь О2 о2

ааух_ у

-„2Л А ( ( Ч

V " V

,1 + са

@аХа

у

2 2 о 2 а о

аауа 3 у

У У 1

У У

Используя соотношение (29), из последнего неравенства получаем оценку вида:

( ' „ V о 2 „о 2 Л Л

о о \ 3

КБ-1Я2 у, у |< £

а=1

(Ьа+Хаса)

Ь1а + '

V Ха У

о а о

аа у 7 + Т у

1

У У а

Чтобы получить в неравенстве (39) наиболее точную оценку для скалярного произведения

К Б-1Я2 у, у

необходимо найти минимум выражения

(Ьа+Хаса)

На + V Ха У

который будет достигнут, если

Ха{ хр> хз) =

Ьа (хр, хз)

Ь1а (хр, х3 )

; а = 1,2,3; х е®р; х5 еа>5;

(2,3), если а = 1; (1,3), если а = 2; (1,2), если а = 3.

Неравенство (39) примет вид

к Бу, у|< £

а=1

(

( „2 _„2Л Л

(Ьа+Хаса)

Ь1а + "

V Ха У

о а о

аа у* у

1

У У

Соответствующее выражение для определения функции й(х) запишется следующим образом:

й (х)= £

а=1

„+1 Г А

а V Ь1а У

а а

2к„ 6

Ьа + са

1

^Ь. Л2

V Ь1а У

Из соотношения (37) и неравенства (36) получаем оценку для параметра А:

Д = тах

а=1,2,3

тах

xрe®р' х5чш5

(Ьа (xР' хз) Ь1а (xР'х^)) 2 + са (xР' хв)

а = 1,2,3, х е ар, х8 е а8,

(Р, ^) =

(2,3), если а = 1 (1,3), если а = 2 (1,2), если а = 3.

й о чЗ

"¡3

и

и

О, С -Й

(39)

(4°)

(41)

(42)

Ранее было показано, что имеет место оценка 5 = 1. Таким образом, определить параметр ю° можно с помощью следующего равенства:

2

т - =-

(43)

2

ТД.

Поскольку сеточный оператор А, входящий в операторную постановку задачи, задаваемую равенствами (8), (9), является самосопряженным, то в итерационной формуле (12) итерационные параметры тп + 1, п = 0Д,...и°(е) - 1 14 следует для ускорения сходимости выбирать совпадающими с чебышевским набором. Тогда для числа итераций п°(е),

требуемого для достижения заданной точности е, справедлива оценка [8]:

4! 1п р ^

I \ V«

(е) = ~йТ

(44)

Если = тах {Л^, И2, N },

где максимум берется по числу узлов сетки по всем трем координатным направлениям, то нетрудно показать, что оценка (44) принимает вид:

' 2'

>(«) = о

N0 1п| -

(45)

4. Результаты численных экспериментов. Для того чтобы построить модельную задачу, будем ориентироваться на среднюю проницаемость нефтеносных пластов, порядок которой составляет 10-12 м2. Это значение взято в качестве базовой проницаемости. Далее базовая проницаемость умножается на неоднородный коэффициент проницаемости, показанный на рис. 1 (сгенерирован случайно из диапазона 0,01-100 для области 150*150 м). Таким образом, разброс значений проницаемости составляет 4 порядка, что соответствует наиболее сложным с вычислительной точки зрения постановкам задач [9]-[12], на обусловленность которых влияют как количество используемых узлов сетки, так и разброс коэффициентов сеточных уравнений.

150

125

100

100.00 51.79 26.83 13.89 7.20 3.73 1.93 1.00 0.52 0.27 0.14 0.07 0.04 0.02 0.01

Рис. 1. Распределение проницаемостей пласта

Будем считать, что:

— толщина пласта — 10 м;

— горизонтальные размеры области — 150*150 м;

— граничные условия — непроницаемые стенки.

Зададим суточные дебиты скважин, в одну из которых (нагнетательную) вода закачивается с дебитом 300 м3/сутки, а из другой (эксплуатационной) — выкачивается нефтеводоносная смесь с отрицательным дебитом 300 м3/сутки. Проницаемость пласта в зависимости от вертикальной координаты не меняется. Исходные данные для наглядности представлены в табл. 1.

Таблица 1

Исходные данные к модельной задаче

Наименование параметра Значение Единицы измерения

Длина области по направлению ОХ 150 м

Длина области по направлению ОУ 150 м

Мощность пласта 10 м

Пористость 0,2

Проницаемость В соответствии с рис. 1

Начальная водонасыщенность 0,1

Дебит на нагнетательной скважине 300 м3/сут

Дебит на эксплуатационной скважине -300 м3/сут

ев И К X ев X

К

Использовалась пространственно-трехмерная сетка с общим количеством шагов: 300*300*20, шаг по времени составлял 0,003-0,005 сут. Сеточное уравнение для функции давления решалось двумя методами: Зейделя (МЗ) и методом верхней релаксации с шахматной упорядоченностью узлов (МВРШУ) [1], а также построенным в данной работе усовершенствованным МПТМ.

На рис. 2 представлено пространственное распределение функции водонасыщенности в нефтеводоносном пласте в зависимости от горизонтальных координат по истечении 20 и 50 суток, а на рис. 3 — распределение давлений. Следует отметить, что полученные результаты удовлетворяют закону материального баланса на дискретном уровне и качественно согласуются с результатами моделирования подобной задачи другими авторами [9].

■ 0.75

Ш 0.70

— 0.65

— 0.60

— 0.55

— 0.50

— 0.45

— 0.40

— 0.35

— 0.30

а 0.25 0.20 0.15

■ 0.75

Ш 0.70

— 0.65

— 0.60

— 0.55

— 0.50

— 0.45

— 0.40

— 0.35

— 0.30

00.25 0.20 0.15

а) б)

Рис. 2. Распределение функции водонасыщенности: через 20 суток (а); через 50 суток (б)

й О чЗ

М

'й

и

и

О, £ -Й

150

125

100

150

Рис. 3. Распределение функции давления через 50 суток (изменяется в соответствии с палитрой от 100 атмосфер до 250 атмосфер)

Результаты численных экспериментов демонстрируют значительное преимущество построенного варианта МПТМ по сравнению с широко используемыми методами — число итераций сокращается в 10-70 раз (табл. 2).

Таблица 2

Количество итераций, требуемое для перехода на новый слой

Номер временного слоя Метод решения сеточных уравнений для функции давления

МЗ МВРШУ МПТМ

1 37535 7232 752

2 23347 3243 314

3 14138 2172 196

4 13190 1918 188

5 12932 1863 184

Заключение. Для 3Б-модели фильтрации двухфазной жидкости выполнена ее дискретизация на основе метода, явного по насыщенности и неявного по давлению. Для полученной системы сеточных 3Б-уравнений построен вариант модифицированного попеременно-треугольного метода (МПТМ). Он учитывает наличие 5-образных источников (скважин) как в структуре оператора-предобусловливателя, так и в оценке энергетической постоянной, определяющей минимально необходимое число итераций при использовании чебышевского набора итерационных параметров. Данный вариант МПТМ применен для решения модельной пространственно-трехмерной задачи в пласте с существенно (на 4 порядка) изменяющейся проницаемостью на сетках, включающих до 200 000 узлов. В сравнении с другими методами данный вариант МПТМ требует меньшего числа итераций (в 10-70 раз). Таким образом, его можно рассматривать как базовый итерационный метод в алгоритме типа IMPES — неявном по давлению и явном по насыщенности.

Библиографический список

1. Азиз, Х. Математическое моделирование пластовых систем / Х. Азиз, Э. Сеттари. — Москва ; Ижевск : Институт компьютерных исследований, 2004. — 416 с.

2. Коновалов, А. Н. Задачи фильтрации многофазной несжимаемой жидкости / А. Н. Коновалов. — Новосибирск : Наука, 1988. — 166 с.

3. Вабищевич, П. Н. Явно-неявные вычислительные алгоритмы для задач многофазной фильтрации / П. Н. Вабищевич // Математическое моделирование. — 2010. — Т. 22, № 4. — 118-128.

4. Сухинов, А. И. Модифицированный попеременно -треугольный метод для задач теплопроводности и фильтрации / А. И. Сухинов // Вычислительные системы и алгоритмы. — Ростов-на-Дону : Ростовский государственный университет, 1984. — С. 52-59.

5. Сухинов, А. И. Двумерные схемы расщепления и некоторые их приложения / А. И. Сухинов // Москва : МАКС Пресс. — 2005. — 408 с.

6. Сухинов, А. И. Модификация метода минимальных поправок для решения сеточных уравнений с несамосопряженным оператором // А. И. Сухинов, А. Е. Чистяков, А. В. Шишеня // Известия Юж. федер. ун-та. Техн. науки. —2013. — № 4. — С. 194-202.

7. Сухинов, А. И. Адаптивный попеременно-треугольный метод для решения сеточных уравнений с несамосопряженным оператором / А. И. Сухинов, А. Е. Чистяков // Математическое моделирование. — 2012. — Т. 24, № 1. — С. 3-20.

8. Самарский, А. А. Методы решения сеточных уравнений // А. А. Самарский, Е. С. Николаев. — Москва : Наука, 1978. — 592 с.

9. Two-Phase Porous Media Flow Simulation on Hybrid Cluster / M. Trapeznikova [et al.] // Lecture Notes in Computer Science. — 2012. — № 7116. — P. 644-651.

10. Бервено, Е. В. Фильтрация двухфазной жидкости в неоднородной среде на компьютерах с распределенной памятью / Е. В. Бервено, А. А. Калинкин, Ю. М. Лаевский // Вестник Томск. гос. ун-та. Математика и механика. — 2014. — № 4 (30). — С. 57-62.

11. Лаевский, Ю. М. Об одном вычислительном алгоритме решения уравнений Баклея — Леверетта /

Ю. М. Лаевский, С. А. Литвиненко // Сибирский журн. индустр. матем. — 2013. — Т. 16. — № 3. — С. 106-115. й

12. Mathematical modeling of flows in porous media / V. R. Dushin [et al.] // WSEAS Transactions on Fluid Mechan- Щ ics. — 2014. — Vol. 9. — P. 166-130. S

¡5

References

1. Aziz, K., Settary, E. Matematicheskoe modelirovanie plastovykh system. [Mathematical modeling of tabular systems.] Moscow; Izhevsk: Institute of Computer Science, 2004, 416 p. (in Russian).

2. Konovalov, А.К Zadachi fil'tratsii mnogofaznoy neszhimaemoy zhidkosti. [Filtering problems of multiple incom-

pressible fluids.] Novosibirsk: Nauka, 1988, 166 p. (in Russian).

3. Vabishchevich, P.N. Yavno-neyavnye vychislitel'nye algoritmy dlya zadach mnogofaznoy fil'tratsii. [Explicit-implicit numerical algorithms for porous media multiphase flow problems.] Mathematical Models and Computer Simulations, 2010, vol. 22, no. 4, pp. 118-128 (in Russian).

4. Sukhinov, AI. Modifitsirovannyy poperemenno-treugol'nyy metod dlya zadach teploprovodnosti i fil'tratsii. [Modified alternating triangular method for heat conduction problems and filtering.] Vychislitel'nye sistemy i algoritmy. [Computation systems and algorithms.] Rostov-on-Don: Rostov State University, 1984, pp. 52-59 (in Russian).

5. Sukhinov, А.I. Dvumernye skhemy rasshchepleniya i nekotorye ikh prilozheniya. [Two-dimensional splitting schemes and some of their applications.] Moscow: MAKS Press, 2005, 408 p. (in Russian).

6. Sukhinov, А.I., Chistyakov, A.E., Shishenya, A.V. Modifikatsiya metoda minimal'nykh popravok dlya resheniya setochnykh uravneniy s nesamosopryazhennym operatorom. [Modification of minimal residuals iterative method for solving grid equations with nonselfadjoint operators.] Izvestiya SFedU. Engineering Sciences, 2013, no. 4,pp. 194-202 (in Russian).

7. Sukhinov, А.I., Chistyakov, A.E. Adaptivnyy poperemenno-treugol'nyy metod dlya resheniya setochnykh uravneniy s nesamosopryazhennym operatorom. [Adaptive analog-SSOR iterative method for solving grid equations with nonselfadjoint operators.] Mathematical Models and Computer Simulations, 2012, vol. 24, no. 1, pp. 3-20 (in Russian).

8. Samarsky, А.А., Nikolaev, E.S. Metody resheniya setochnykh uravneniy. [Solution methods of finite-difference equations.] Moscow: Nauka, 1978, 592 p. (in Russian).

9. Trapeznikova, M., et al. Two-Phase Porous Media Flow Simulation on Hybrid Cluster. Lecture Notes in Computer Science, 2012, no. 7116, pp. 644-651.

10. Berveno, Е.V., Kalinkin, A.A., Laevskii, Y.M. Fil'tratsiya dvukhfaznoy zhidkosti v neodnorodnoy srede na komp'yuterakh s raspredelennoy pamyat'yu. [Two-phase fluid filtration in nonuniform media on clusters.] Tomsk State University Journal of Mathematics and Mechanics, 2014, no. 4 (30), pp. 57-62 (in Russian).

11. Laevskii, Y.M., Litvinenko, S.A. Ob odnom vychislitel'nom algoritme resheniya uravneniy Bakleya — Leveretta. [On a numerical algorithm for solving the Backley-Laverett equations.] Journal of Applied and Industrial Mathematics, 2013, vol. 16, no. 3, pp. 106-115 (in Russian).

12. Dushin, V.R., et al. Mathematical modeling of flows in porous media. WSEAS Transactions on Fluid Mechanics, 2014, vol. 9, pp. 166-130.

Поступила в редакцию 02.12.2015 Сдана в редакцию 02.12.2015 Запланирована в номер 22.01.2016

й о чЗ

"¡3

и

<и

>

Л £ -Й