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

ИНФОРМАТИКА, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И УПРАВЛЕНИЕ INFORMATION TECHNOLOGY, COMPUTER SCIENCE, AND MANAGEMENT

УДК 519.6

DOI 10.12737/20944

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

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

1 Донской государственный технический университет, г. Ростов-на-Дону, Российская Федерация

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

3 Южный федеральный университет, г. Ростов-на-Дону, Российская Федерация

Mathematical modeling of two-phase compressible fluid filtration based on modified adaptive method of minimum

amendments***

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

1 Don State Technical University, Rostov-on-Don, Russian Federation

2 North Caucasus Federal University, Stavropol, Russian Federation

3 Southern Federal University, Rostov-on-Don, Russian Federation

Ö О T3

M

'S

M

(U

Ü £ Л

Целью работы является построение и исследование адаптивного метода минимальных поправок (МАММП), предназначенного для численного моделирования процессов фильтрации двухфазной сжимаемой жидкости в пористых средах. Описанный подход позволяет преодолеть известные ограничения использования других методов решения сеточных уравнений, определяемые такими условиями, как: значительные перепады давления, действующего на нефтеводоносный пласт; сжимаемость среды при значительном содержании газа в нефтяной фазе.В качестве основы исследования выбран метод аппроксимации — явный для определения функции водонасыщенности и неявный для расчета функции давления.При постановке начально-краевой задачи и ее дискретизации рассмотрен процесс фильтрации сжимаемой двухфазной жидкости в пространственно-трехмерной области с боковой поверхностью, ограниченной снизу поверхностью подошвы пласта, а сверху — поверхностью его кровли. Построен двухслойный итерационный метод вариационного типа — модифицированный метод минимальных поправок, адаптированный для решения сеточных уравнений двухфазной сжимаемой жидкости с несамосопряженным оператором при самых общих предположениях относительно свойств оператора сеточной задачи. Показано, что МАММП обладает асимптотической скоростью сходимости, характерной для «классического» ПТМ, не использующего технику чебышевско-го ускорения и применяемого для задач с самосопряженным оператором. Численные эксперименты подтвердили высокую эффективность МАММП. Установлено, что для достижения заданной точности число итераций при МАММП сокращается в 3-20 раз по сравнению с методом Зейделя и методом верхней релаксации.

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

The work objective is to build and investigate the modified adaptive method of minimum amendments (MAMMA) which is destined for the numerical simulation of the two-phase compressible fluid filtration in porous media. This approach allows overcoming the known use limitations of other methods of the finite-difference equations solution, such as: crucial differential pressures acting on the oil-and-water bearing formation; and the compressibility of the medium at the considerable gas content in the oil phase. An approximation method - an explicit one for defining the function of water saturation, and an implicit one for the pressure function computation - is selected as the research basis. When setting the initial boundary value problem and its sampling, the process of the two-phase compressible fluid filtration in the space-dimensional domain with the lateral area bounded below by the subface of stratum, and above - by the bed top, is considered. A two-layer iterative method of the variational type - a modified method of minimal amendments adapted for solving finite-difference equations of the two-phase compressible fluid with a non-selfadjoint operator under the most general assumptions on the properties of the grid-problem operator is built. It is shown that a MAMMA has the asymptotic convergence rate characteristic of the "classical" alternate triangular method that does not use the Chebyshev acceleration technique and can be applied to the problems with a self-adjoint operator. Numerical experiments have confirmed the high efficiency of MAMMA. It is established that to achieve the specified accuracy, the number of iterations at the MAMMA reduces to 3-20 times as compared to the method of Seidel and the overrelaxation method.

Keywords: two-phase compressible fluid filtration problems, finite-difference equations, non-selfadjoint operator, adaptive method of minimum amendments.

* Работа выполнена при финансовой поддержке грантов РФФИ: № 15-01-08619, № 16-07-00100, № 15-07-08408, № 15-07-08626.

** E-mail: [email protected], [email protected], [email protected]

*** The research is supported by RFFI grants (nos. 15-01-08619, 16-07-00100, 15-07-08408, 15-07-08626).

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

В данной работе, как и в [1], авторы ориентируются на метод аппроксимации — явный для определения функции водонасыщенности и неявный для расчета функции давления (/МРЕХ-метод) [2]. В этом случае в построении и численной реализации дискретной модели для определения функции давления возникают дополнительные трудности, связанные с изменением плотности более сжимаемой нефтяной фазы. Это, в свою очередь, приводит к изменению свойств оператора уравнения — он становится несамосопряженным [3], [4], и в общем случае невозможно гарантировать ограниченность сверху величины сеточного числа Пекле, например, константой 2. Поэтому ни один из вариантов ПТМ, построенного в работах [1], [5], [6], не пригоден для численной реализации модели. Сравнение построенного метода с употребительными подходами к решению сеточных уравнений для восстановления функции давления (метод Зейделя, метод верхней релаксации) показывает существенное преимущество модифицированного адаптивного метода минимальных поправок (МАММП). Выигрыш в числе итераций составляет 3-40 раз в зависимости от числа ячеек используемых сеток, что позволяет использовать МАММП в качестве базового алгоритма решения сеточных задач проектирования разработок нефтяных месторождений.

1. Постановка начально-краевой задачи и ее дискретизация. Как и в [1], рассмотрим процесс фильтрации сжимаемой двухфазной жидкости в пространственно-трехмерной области О с боковой поверхностью £( х, у), (х, у )ег, г = Б \ б, ограниченной снизу поверхностью подошвы пласта

7 = Нь (х,у), (х,у) е Б, Б е Я2,

<и К

а сверху — поверхностью кровли пласта щ

<и

ч ю ей

Запишем уравнения неразрывности для каждой из двух фаз в отдельности в дивергентной форме. ^

5(тРЛ ) + ^ Г и 1=_р.д,. , I = ^ (1)

z = Hr ( х, y ), ( х, y )e d.

Я

I ' 11 | ' I J-1 ' ' ' v ' ^

Ot { ) И

к

где индексы i = 1,2 обозначают соответственно нефтяную и водную фазы; рг-, si — плотность и насыщенность соот- ^

<U

- к

ветствующих фаз; ui — скорость их фильтрации; qi — мощность объемных источников (стоков) каждой из фаз. <§

К

Будем ориентироваться на закон Дарси, который запишем в виде уравнений для каждой из фаз в отдельности: ¡2

н

U = - k • — | gradp - Pig I, i = 1,2, ч

(2) !

3 ю

g и

где k — коэффициент абсолютной проницаемости; f = f (s2) — относительная проницаемость соответствующей g

|

фазы; g — вектор ускорения свободного падения; - • gradp = f kh —, kh —, kv — | , - •g =(0,0, kvg)T . &

I oX oX Oz J

X

Уравнения (1) и (2) дополняются уравнениями состояния для каждой из фаз: К

= -k •f f gradP - P i g ), i = 1,2,

g T

k = (kh, kh, kv )

Pi = Pоi [ 1 + Pi (P - Po )], i = 1,2, (3)

где р0/, р0 — соответственно характерные значения плотности и давления, р1 — коэффициенты сжимаемости сред.

Предположим, для определения функций / (¿2) и ^(¿г) используются многочлены третьего порядка относительно функции водонасыщенности ¿2 [2]. Чтобы получить уравнение для определения давления, подставим в каждое из двух уравнений (1) выражение для скорости фильтрации соответствующей фазы из уравнений (2). Таким образом, придем к соотношениям

д(рт ) = ^ ( р1/1 (¿2)

д/

-к-I gradp-Р1 g

л

д(р2т'Ъ ) = л f Р2Л2 (¿2 )

дл

-кgradp-Р2 g

т

Ж

т

(4)

(5)

Разделим уравнение (4) на произведение функций тр1, а уравнение (5) — на тр2 с учетом уравнений (2). Складывая преобразованные уравнения, получим:

Ср^ Ф2. \дР = —с1„

р1 Ср р2 Ср ) д/ тр1

1

Р1Л (¿2 )

М>1

к -I gradp- р- g

1

тР2

С/у

Р2 Л2 (¿2 )

Д 2

к -I gradp-р2 g

(6)

р-т

Ч2 р 2 т

При естественном характере сжимаемости имеет место неравенство

Срг

Ср

> 0, / = 1,2.

(7)

При этом условии соотношение (7) является уравнением параболического типа, содержащим в общем случае младшие производные (первого порядка по пространственным переменным) относительно плотности. Учитывая (3), уравнение содержит первые производные относительно давления, а также переменные множители вида 1/тр/, i = 1,2 при операторах

С/у

р,/1 (¿2 )

к -1 gradp-рг g

/ = 1,2.

Итак, дифференциальный оператор в правой части уравнения (6) является несамосопряженным [3]. Получена система выражений, включающая уравнение (5) или вместо него — (4), уравнение (6), уравнение состояния (3). К данной системе необходимо добавить соответствующие начальные и граничные условия, сформулированные в [1], а также уравнение (4). Для численной реализации данной модели целесообразно провести линеаризацию задачи (по временной переменной), выделить в операторе задачи на верхнем временном слое симметричную и кососимметрич-ную части, а затем выполнить дискретизацию. Построим неравномерную временную сетку:

шт = |^ =Д тК 1 < п < N

Для водонасыщенности s2 опустим нижний индекс и построим явное разностное уравнение (по времени). Будем считать пористость т не зависящей от давления.

Из уравнения состояния (3) получено равенство

+ ¿2. (1 - ¿)¿Р2

р1 СР р 2 СР 1 + Р1 (Р - Р0 ) 1 + Р 2 (Р - Р0 )

(8)

Преобразуем первое слагаемое в правой части уравнения (7), выполняя линеаризацию на временной сетке сСТ. При этом отнесем линеаризованные коэффициенты на нижний временной слой (/ = 1п ), а члены, определяющие давление, — на верхний временной слой. Для определенности / = /п+1 в следующих ниже выражениях эти члены отмечены символом «а» над функциями).

й о тз

'й -м

и

<и

а £ .й

1

С/у

тр1

р1Л1 (5)

=-р1 С/у

тр1

ш

к -I gradp - р1 g

к-1 gradp- р01 [1 + Р1 (р - Р0 )] g

Ро 101 к' 8гаФ- ~ Ро 1 [! + Р1 (Р ~ Ро )] 8 ■

V

1

т

Р018 д

т дг

ЛМ №

к • gradp

Р1Р018 д (, Л (5)

т дг

к V Д

-(р - Ро )

Р018 д Л (5) ] + к Р1Р01 Л (5) др др 1 ^ Р1Р01 Л (5) др др

Д1

к ' + кк ' ■ тр1 ц дх дх тр1 ц ду Йу

^ р^ т др .др - к, р^АМ [1+в1 (р - Р0 )]др=

тр1 ц дг дг тр1 ^ дг

д_

дх

V V

/1 (5) др И дх

д + —

ду

/1 (5) др И ду

д + —

дг

/1 (5) др И дг

Р1Р01 Л (5) др др , ,, Р1Р01 Л (5) др др+ к Р1Р01 Л (5)

+к^^01^^^ .^1 + кк

тр1 и дх дх тр1 и ду ду

т и

1 др Р1 дг

др д

Р1Р018 д( Л (5)

т дг

к

V И

р - К

Р012Р128 Л (5) др

тР1

И дг

(

р +

+к Р012 (Р1 р0 -1)8 Л (5) др+ (Р1 р0 - 1)Р018 д [ь Л (5)

тР1

И дг

дг

к

V И

(9)

Аналогично предыдущему приходим к равенству:

1

тР2

к / 8~аф - Р 2 я 1

Д2 V ))

Р2Л2 (5)

( я (

дх.

V V

Л (5) др

д + —

ц 2 ^) ду V Д 2 ду

Л (5) др

д + —

дг

+к Р2Р02 Л2 (5) др. др+ к Р2Р02 Л (5) др. др+ к Р2Р02 Л^С51 др

к тР2 ц2 йх дх к тР2 ц2 ду ду

./2 (5) др Д2 дг 1

Р2 дг

Д 2

др д

Р2Р028 д (к )

V У Д 2

т дг

\ к Р022р228 У2(5)дрр +

тР2

Д 2 дг" (

+к Р022 (Р2рр -1)8 Л (5) др + (Р2р0 -1)Р028 д /2 (5)

к

V Д2

(10)

тР2 ц2 дг т дг

Подставляя (8), (9) и (10) в соотношение (6), получаем линеаризованное дифференциально-разностное уравнение относительно давления с младшими производными:

г (1 -

*Р2

1 + в (р - р0 ) ' 1 + Р2 (р - р0 ).

- р

2

= 2

г=1

дх.

V V

Л (5) др

дх

д + —

ду

(

Л (5) др

Дг ду

д + —

дг

(

Л (5) др

+к Р^ЛМ др. др+ к Р^ЛМ др. Ф+ к Рг Р0г Л (5)

к тРг дх дх к тРг ду ду У т

РгР0г8 д(, Л (5)

г г-г

( Г1Л\

Дг дг,

' А. др

V Рг дг

др ~дг

(11)

т дг

р - к

Р0г'Рг28 Лг (5) др

тРг дг

р +

+кУ

Р0г2 (Ргр -1)8 Л (5) др , (Рг-р0 - 1)Р0.8 д (. Л (5)

тР.

Дг дг

дг

Дг

Чг

Ргт

Для функции водонасыщенности построим на временной сетке явное относительно 5 уравнение:

т (Р25 - Р25^ ( Р2Л2 (5)

^+1

к ( к

Д 2

к А 8radp - Р 2 8

42 т

которое можно переписать в развернутом виде:

щ

5 I

(и

4

И ей

а

с

^

5 ей И 5

1

X Щ

Н «

ей X Л

4 (и н

5

4 О

5

(Г

2 СП

ей И 5

<3

а о

-е х X

Тп+1 т

к

дх

V V

/ (^) др

М2 дх

ду

/ (^) др

М- 2 дУ,

+ — дх

/2(^) др

М 2 дх

+к р^Л^)др.+к Ё2Р01 /2(5)др.+к Р2Р02 Л2(^)( 1 др ^

к тр2 м2 дх дх к тр2 м2 ду ду у т

М 2

Р2 5х

др) ~дх

вм ^(к Л2 (^)

V У М 2

т дх

\ к Р022Р22 g /2 (¿)дрр +

р ку ^Гр +

тР2

М 2 дх"

(12)

+к Р022 (Р2р0 ~ !) g Л2 (^) др | (Р2р0 - 1) Р02g д Л2 (^)

т

дх

\

к,

V М 2

Ч2

Р2т

(13)

тР2 М 2 дх

Для вычисления плотностей водной среды на нижнем и верхнем временных слоях используем уравнения состояния вида (3):

р2 = Р02 [1 + Р2 (р - р2 )] > Р2 = Р0 / [1 + в2 (р-р2 )] .

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

' (1 - ¿Ж , ^2

1 + Р1 (р-р0) 1 + Р2 (р-р0)

2

+2

/=1

Р/Р0/^ д

дх V V

т дх

/ (')'

/ 2п 2

п+1

р+ку № др р -

тр/ М/ дх

(14)

к'

■Ш др

М/ дх

+— ду

к

V

/ (*) др р

М дУ

+— дх

/ (^) др М/ дх

Ч/

Р/т

являются самосопряженными. Данное свойство сохраняется при использовании дискретизаций этих выражений на пространственных сетках — например, так, как это описано в [7]. Аппроксимация выражения (14) на сетке и одновременно определение самосопряженной части оператора сеточной задачи на сетке имеет вид:

(-*Р и ■ Ь ВЩ, (^ - р) вх,

х 2 2

12В3.,(Р+и,к - РМ) + к12В'(2-- 1,к (Р1,к - Р}-1,к )-

"х 2

--В(2). (Р,1+1,к - Р, 1,к ) +

К

/,],к р,],к+1 Iв 1

' ',] ,к+- к 2

+—(Р 1к - Р к

В(3)

1т (

х

( (1 - 1 )Р1 т_

ч 1 + Р1 (Р,М - р0 ) + 1 + Р2 (Р,М - р0 )

1 (ДМ-1 - Ам) в(3) 1 +

],к Р2

т,

/, 1 ,к

п+1

Й О ТЗ

М

"¡3

ел

и

О, £ -Й

2

+ 2

а=1

РаРоа g ( кУ', 1,к+1Ла /, 1,к+1 куг, ¡,к-\Ла/, ] ,к-1 2кх Ма

2к,

РИк +

Р,1,к +

+к

Р0/ Ра g /а1,к

( Р

ра 1,к Ма

/, 1,к+1 /, 1,к-1

Р

2к

/, 1,к

0,

В& 1 кР11 1 к, если узел находится на нагнетательной скважине; В&, 1,кРг,] к, если узел находится на эксплуатационной скважине.

При произвольных значениях шага по времени тп+1, а также при неизвестных заранее знаках производных

вида

Л (5)'

др

и — нельзя в общем случае гарантировать положительность сеточного оператора (15), аппрокси-

дг

мирующего соответствующий дифференциальный оператор. Достаточное условие положительности сеточного оператора, аппроксимирующего дифференциальный оператор в соотношении (14), принимает вид:

' (1" *г,М )в , " ^

5г', ],к в2

т п+1 ^

min —

(хг у,гк )6ю £

а=1

1 + в1 (рг,],к - р0 ) 1 + Р2 (рг,],к - р0 )

РаРоа8

куг,],кЛа г,],к

^ р0аРарг, г,к ^ 1 +--—

Ра ¡, г ,к

],к У

=о кг.

(16)

кгтг, ] ,к Да

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

Обозначим

к Рг^ М*) др. др+ к РгР01 ЛМ др. др.

к тр. д. дх дх к тр1 д. ду ду

+к

Рг Р0г Л (5)

т Дг

' 1 др.

V Рг дг

др, г = 1,2.

дг

(17)

Ьа,1 (х,у,г,^) = кк ^АМ дХ, (х,у,г)б О, Ра Да дг

Ьа,1 (х + 0,5кх, У, г, 1п ) = Ьа,1 (х. + 0, 5кх, у] , гк , 1п ), (хг, у] , гк ) 6 юк

'др (х. + 0,5кх, у], гк, ^ )1 Р+1,],к - Рг,],к

юк = ю1 х ю2 х ю3,

дх

Ьа,1 (х - 0,5кх, У, г, гп ) = Ь1 (хг - 0,5кх, у}, гк, гп ), (хг, у}, гк )

( -л / « .

: Юг,

ю к+ = ю1 х ю 2 х ю 3,

ф (х - 0,5кх, у], гк, гп) _ Р,],к - Рг-1,],к

дх

К

а = 1,2,

Ьа,2 (х,у, г,Гп ) = кк ^^др, (х,у, г)е О,

Ра Да 5у

Ьа,2 (х, у + 0, 5ку, г, гп ) = Ьа,2 (хг, у} + 0, 5ку, гк, гп ), (хг, у} , гк ) (

ю

юк = ю1 х ю2 х ю4,

др(х.,у] + 0,5ку,гк, 1п)^ Рг ]+1,к -РР,],к ду

Ьа,2 (х, у - 0,5ку, г, гп )= Ь1 (хг, у] - 0,5ку , гк, гп ), (хг, у], гк )6 Юк

Рг,],к - Рг,]-1,к

ю к = ю1 х ю 2 х ю 3,

Гдр(х,у} -0,5ку,гк,)Л

ду

а = 1,2.

(18)

(19)

и К X <и ч и

ей Л

С

^

К ей И К

X

*

н ч

ей X л ч и н К

ч о К

Е и

ей

и

К ей

л о

X

К

. / ,\,а fa (5) f 1 Ф 1 / ч Г

ba,3 (X У, z, tn ) = kvßaP0a-— I -Т--g , (X z )6 G

Да I Ра dz )

ba,3 (X У, z + 0,5hz , 1п ) = ba,3 (X , У] , zk + 0, 5hz , tn ), (Xi, У ] , zk ) 6

wh = fflj x ю2 x ю3,

dp (xt, У], + 0,5hz, tn)

5z

P, ],k+1 Pi, j,k

(20)

ba,3 (X У, z " 0, 5hz , tn ) = b1 (Xi, У] , zk " 0, 5hz , tn ) , (Xi, У] , zk )

6 Wi,

w = w1 x w 2 x w3

dp (X, У] - 0,5hz , zk , tn ) N = P,],k - P,],k-1 dz "

a = 1,2.

С учетом введенных обозначений (18)-(20) сеточный аналог непрерывного выражения (17), умноженный на функцию т,-, у, к, принимает вид:

(AP)

2

i, ],k a=1

ßaP0a fa (5) dp dp ßgpQa fa (5) dp dp

• _ -

+kvßaP0<

fa ( 5 )f 2. dp

Pa dz

dp) dz

Sr dx Pa ^a dV dV

1 [(ba,1 (X + 0,5ÄX, У, z, tn ) Px +

2

= z

a=1

(21)

(х - 0,5кх, У, z, 1п ) Р-х ) + (^ (x, у + 0,5ку, х, 1п ) Ру + Ьа 2 (^ у - 0,5ку, х, 1п ) Ру ) + + (Ьа,з (х,у, х + 0,5кх, /п )Р2 + Ьа,з (х,у, х - 0,5кх,/п )Рх )]], (х/, у} , хк )6 ®1 Х Ю2 Х Ю3 .

Таким образом, в операторном виде сеточную задачу, аппроксимирующую уравнение (11), можно записать в

виде:

(22)

AP = A P + A1P = f (P,s),

A = A0 + Aj, A0 = A0*, A = -A*.

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

2. Построение адаптивного метода минимальных поправок. Далее построен двухслойный итерационный метод вариационного типа — модифицированный метод минимальных поправок, адаптированный для решения сеточных уравнений двухфазной сжимаемой жидкости с несамосопряженным оператором [7] при самых общих предположениях относительно свойств оператора сеточной задачи. Требуется лишь, чтобы в явном виде были выделены симметричная и кососимметричная часть оператора сеточной задачи. Сформулируем задачу в общем виде. В конечномерном Гильбертовом пространстве H рассматривается задача на обнаружение решения операторного уравнения:

AP = f, A : H ^ H , (23)

где A — линейный, положительно определенный оператор, определенный выше соотношениями (15), (21), (22),(A>0).

Для нахождения решения задачи (23) используем неявный итерационный процесс:

й о тз

"¡3

и (U

Ü £ -Й

B

P m +1 - P п

. + APm = f, B : H ^ H .

(24)

Здесь т — номер итерации; ат+1 > 0 — итерационный параметр; В — легко обратимый оператор. Последнее означает, что обращение оператора В в (24) должно быть существенно проще, чем обращение исходного оператора А. В нестационарном итерационном методе вариационного типа

B

P m+1 _ P п

. + APm = f

(25)

итерационные параметры ат+1 выбираются из соображений минимизации скалярного произведения для погрешности решения в энергетическом пространстве ИП [7]:

а

а

гт = Рт — Р

т+1 т+1

гт+1 = Р"+1 — Р

( Бг

) — min, Б = Б > 0.

Учитывая условие (26), а также представление гт+1 = гт — ат+1В 1 Azт , приходим к задаче

которая равносильна задаче

(Бгт - ат+11 Агт, гт - ат+1 1Лгт ) — тт,

(Бгт, гт ) — ти+1(РВТ 1Azт, zт ) — Т"+¿Р*", В А" ) + + т"+1 (РВ ^1А2т, В 1А2 т ) — тт.

(26)

(27)

(28)

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

ат+1 ="

(БВ 1 Агт, гт ) + (Бгт, В Аг" ) 2( БВ 1 Агт, В ^Аг" )

(29)

Дальнейшее зависит от выбора энергетического пространства НБ . В частности, итерационные методы вариационного типа (27), (29), для которых Б = Лл"В 1А приводят к следующим равенствам:

ат+1 ='

(А В-1 АВ—Агт, гт ) + (А*В-1 Аг", В-1 Агт ) 2(А* В 1АВ1 Агт, В Агт )

_ (АВ1 Агт, В-1Агт ) + (А*ВАгт, В^Аг" ) _ = 2( В 1АВ1 Агт, АВ1 Агт ) =

(АВ—Агт, В 1 Агт ) + (АВ—Агт, В 1 Агт ) (АВ Агт, В ^ Агт )

2( В 1 АВ"1 Агт, АВ1 Агт )

(В-1 АВ-1 Агт, АВ- Агт )

Следовательно, итерационные параметры определяются из соотношений:

(Лwm, wm )

=

(В- А^т, А^т)'

В^т = Ахт - Л т = 0,1,...

(30)

3. Исследование сходимости метода минимальных поправок. Если, как это принято, w = В- Агт — вектор поправки, а векторы погрешностей определены равенствами (26), то из уравнений (23), (24) получаем:

Вгт+1 = Вгт - т^, Azт .

(31)

Представим последнее соотношение в виде:

wm+1 = wm - ат+В1 Awm .

Выполним замену у" = в12wm, С = В 12АВ 12 и из последнего равенства получим:

что равносильно равенству

В-\пут+1 = в-1/2 Ут —ат+В 1 АВ-1/2 V"

у"+1 = у" -а

(32)

!СУ" =(Е — ат+1С)У" . Используем оценку ||ут+^|:

||у"+1 = ||( е — ат+1с) у"| = ||((0 т+1е — ат+1с0 ) + ((1 —0 т+1) е — ат+1с1)) ут\ |.

Выполняя очевидные преобразования, из последнего равенства получим С = С0 + С1, где С0 = С0*, С1 = —С1*. Имеем очевидную оценку:

\у"+11 <

_ т+1

е —

С

т+1

'0

'т+1 У

+ ||((1 —0т+1) е — ат+1с1) у"\ i.

(33)

Проведем оценку первого слагаемого в правой части неравенства (33):

и

к

X <и ч ю ей Л

С

^

к

ей И

к

X

*

н ч

ей X л ч и н к ч о

к Е

ю ей

и к

ей

а о

X

К

((

\\

( а т+1

Е -

V 0 т+1

В1'2 - а т+1

0 т+1

Е — а"+1 с

V т+1 У

Л (

В1'2 w"

w

\ \

В1/2 В -1/2 А0

А 0

"т+1

w

/ У

= ( В1/2 wm, В1/2 wm )- 2 -а т+1

0 т+1

(а^т,wm)+ а-] (лоwm,b-1лоwm).

(34)

Примем во внимание, что минимум дроби ат+1/0 1 достигается при

(А0 wm, wm )

0т+1 (В 1А Wm , А0Wm )

Тогда оценку (34) запишем:

( а ^ Е -^"+1 С,

А 0

V т+1 У

'(Б

/ о1/2 т т-А/2 т\ -ч

= (В w ,В w )-2

(л т т\

|А W , W I ( В- А0 wm, А wm )

(л т т \

А^ , W I -

(А0т

)

(В~1А0 ^^,, wm) = (у", у" )

(Л0wm ,т )2

(А wm, В-1 А0 wm ) = ( В1/2 wm, В1/2 wm '

(В-1А0 ^^,, wm )

(сУ" , У" С у" , С0 у" )

)2

1 -Л

С0 у",У")

(С0 у", С0 У" )(У",У")

(35)

Таким образом, из (35) получаем равенство — оценку для первого слагаемого, стоящего в правой части неравенства (33):

( а Л

Е-с0 0 0

V т+1 у

1

(С0У",Ут )2 , (с0 у", с0у")(у",у")'

(36)

Из соотношения (33) для второго слагаемого с учетом равенства (су" ,у" ) = 0 получаем: ||((1 - 0 т + 1) Е - а"+С ) У"|2 =((1 - 0 "+1) у" , (1 - 0 т +1) у" )-

- (ат+1С1У" ,(1-0" + 1) У")+(ат + 1С1У" ,ат+Су" ) = = (1 - 0т+1 )2 (у", У" ) + 0"

32 I ат+1 п т ат+1 п т

^т+1 IТ С1У С1У .

(37)

т+1 т+1

С учетом полученных неравенств (36), (37) из соотношения (33) имеем:

\\ут+1 < (

( а Л

Е-^т+1 с

0 0

V т+1 У

+ ((1 - 0„+1) Е - а „+1С1) Ут\\ =

й о

ТЗ

'й

и

и

О, £ -Й

1

(с0у",у")

2 Л

(с0 у", С0 у")(у",у")

у +

+ 1(1 - 0т+1 )2 (У" , У" ) + 0

2

т+1

(

т+1 г<ат+1 /-<,,т

Л

Су" _m±L Су

1 0 1

V ^ т+1 °т +1 У

а

}2

т+1,

, (

2

Сгут, ут )

(С0ут, С0ут)(ут,ут)

V" +

(1 - 0т+1 )2 + 0

2 , п2 т+1

ат+1 Сут ат+1 Сут

т+1

(ут,ут)

Введем обозначения:

т+1

1 <

С0ут,ут)

2 Л

(с0 ут, С0 ут)(ут,ут)

1 т+1

ат±1_ су , 5«±1 С1ут ® т+1 ® т+1

(ут, ут)

Определим итерационный параметр следующим образом:

(Ао (йт ,(йт) 0

т+ (В- А0 шт, А0 шт)' Оценка для ||У™+1|| с учетом соотношений (38) и (39) примет вид:

||т|| = (е т+1 ¿т +1 +71 - 20т +1 + 0 2+1 (1 + У„ +1 ))| .

1 - п

Положим, 0т+1 =-т+1. Тогда из последнего равенства получим оценку:

1 + У т+1

V =

Таким образом, последнее неравенство запишем в виде

т+1 ^ V <

1 ^т+1 „ , ¡Ут+1 + 'Пт+1

' »-> 111^1 ~Г А

(38)

(39)

(40)

1 , т-п -I/ ,

1 + Ут+1 \ 1 + Ут+1

\ у

Чтобы выбрать оптимальное значение параметра цт+1, вычислим производную от выражения, входящего в

правую часть неравенства (40), и приравняем ее нулю:

1 - Пт+1

1 + У т+1*

¡Ут+1 + Пт+1

I 1 + Ут+1

"У

Лт+1

Пт+1

Пт+1

1 + 1т+1 ^ + У т+1)( У т+1 + П1+1)

= 0.

(41)

Заметим, что

1 У т+1 „ , У т+1 + Пт+1 , , , т+1 , ,

1 + кт+1 V 1 + Ут+1

> 0.

Пт+1 ' Пт+1

Тогда точка минимума выражения, входящего в неравенство (40), достигается при оптимальном значении параметра Пт+1 :

Пт+1 =

¿т+1Ут+1

Ц1 + У т+1 - ¿т+1)

(42)

и

к

X <и ч и

ей Л

С

^

<й И

к

X

*

н «

ей X Л

ч и н к ч о

к

Е

и

<й и к

<3

а о

X

К

С учетом полученного оптимального значения цт+1, подставляя (42) в (40), получим оценку для У

,т+1

\\ут+1 <

1 -

5т+1Ут+1 1 + Ут+1 - 5т+1

1 + У т+1

+

Ут+1 +

5т+1Ут+1 1 + Ут+1 - 5т+1

1 + Ут+1

У =

5т+1 5т+1

Ут+1

1 + Ут+1 5т+1

т +1

1 + Ут+1

1 + Ут+1 А'т +1

+«

Ут+1 (1 + Ут+1 5т+1)

1 + Ут+1

Константа р , определяющая скорость сходимости итерационного метода, оценивается сверху следующим

образом:

Р<

5т+1

Ут+1 (1 + У т+1 5т+1)

1 + т+1

Скорость сходимости метода минимальных поправок зависит от:

(43)

5т+1 =

1-^

С0 у",у")

2 Л

(с0 у", С0у" )(у",у")

1

(А wm, wm )2

(А0 В- 1А0 wm, wm )(Bwm, wm )

Пусть Ат1п (С0), Л,тах (С0) — соответственно минимальное и максимальное собственные числа оператора С0 а v = Xтах (с0 )/X(с0) — число обусловленности оператора С0.

Применим неравенство ху < (ах + у / а)2 /4, справедливое для а ф 0 . Тогда

5 т+1 < 1-

4 (А0 wm, wm )2

( а (А В-1А0 wm, wm ) + (Bwm, wm ) / а )

< 1 —

\2 '

(а Лпах (С0 ) + 1/ «Лтт (С0 ))

Положим, а = ^Хт1П (С0) / Хтах (С0) . Тогда выражение, стоящее в правой части соотношения (43), меньше единицы при 5т+1 < 1. Для 5т+1 и ут+1 имеют место оценки:

~4 V-1

й о тз

"¡3

и

и

О, С -Й

5т+1 <

1

1

(v1/2 +v-1/2 )2 V + 1

(44)

У т+1

(в-1 , ЛlWm ) (В-1 ЛlWm, ЛlWm )

-4 (1-52+1 )<У '

(В-1А wm, А wm )'

(-1 т л т\

В Л0W , Aw )

Введем обозначение

(Б1Л1wm, )

кт+1=(Бzcлоwm:^wm), у-+1=кт+1 (1-5"+1).

ч

/

4

Выполним минимизацию оценки сверху постоянной р. Вычислим производную по 'т+1 от выражения в правой части (43) и, учитывая (45), получим:

( . Л 2 \ Г,-П-ТТ V

'т+1 +(1 'т+1 )\/кш+1 (кт+1 + 1)

1 + к

т+1

(1 — +1 )

(1 + кт+1 ) 2 'т+1 д/кш+1 (кт+1 + 1) + кш+1'

2

т+1

(1 + кт+1 (I- 'Ш+1 ))2

> 0.

(46)

Из неравенств (43) и (46) получаем оценку:

' тах + ^

Ут+1 (^ + Ут+1 'тах )

Р " ,

1 + Ут+1

Выполняя тождественные преобразования — умножая на ^ 'тах у т+1 + т+1 (1 + У т+1 - ¿тж) ^(1 + У т +1) литель и знаменатель выражения в правой части неравенства (47), получим:

(47)

чис-

Р ^

|^'тах + ^ У т+1 ( 1 + Ут+1 'тах ) 1^'тах Ут+1 + ^ У т+1 ( 1 + У т+1 'тах ) |)/(1 + У т+1 )

'тах У т+1 + \1 У т+1 1 I1 + У т+1 — '2 ) тах

[ У т+1 (1 + У т+1) + 'тах (У т+1 + 1) ^ Ут+1 ( 1 + У т+1 'тах ) |]/(1 + Гт+1 )

'тах У т+1 + М У т+1 (1 + Ут+1 'тах )

В результате из неравенства

Р ^ '„

У т+1/ 'тах ^У т +1 (1 + У т+1 'тах )

'тах У т+1 + ^ У т+1 ( 1 + У т+1 — 'тах )

в силу (44) получим:

^Ут+1 (1 + У т+1 5тах ) + Ут+1 ^У т+1 (1 + Ут+1 — 'тах ) — Ут+1

Ут+1 (1 + Ут+1 ¿тах ) + Ут+1 ^Ут+1 (1 + Ут+1 — ¿тах ) — Ут +1

+1

Введем обозначение

V = V

У т+1 (1 + У т+1 'тах ) + У т+1

^У т+1 (1 + У т+1 'тах ) У т+1

Тогда с учетом оценки (48) и введенного обозначения получаем

* 1

V -1

Р<-

*

V +1

Равенство (49) с учетом (45) примет вид:

„*_„м

кт+1 (1 'т+1 + кт+1 (^ 'т+1) 'тах )+ кт+1 (^ 'т+1)

У кт+1 (^ 'т+1 + кт+1 (^ 'т+1) 'тах ) кт+1 (^ 'т+1)

\]кт+1 0 + кт+1 ) + кт+1 (^кт+1 (1 + кт+1 ) + кт+1 ) , г—-

' / =-= ^-к--+ кт+1 +

Л1кт+1 (! + кт+1 ) — кт+1 кт+1

■\1кт+1 )

(48)

(49)

(50)

и К X <и ч и

ей Л

С

^

ей И К

X

*

н ч

ей X Л

ч и н К

ч о К

Е и

ей И К

Й

а о

X

К

k _Ш1 <

Km+1~ ||2 -

Таким образом, в постоянную р, определяющую скорость сходимости метода минимальных поправок в соответствии с неравенством (50), входит величина

v* =v(VlT

km+1 +1 ) , (51)

где v — число обусловленности матрицы C0, а km+1 характеризует отношение нормы кососимметрической части оператора A к симметрической:

II н2

' m .«М"

(АwmfB_1' \\C0Л,

На этом построение и оптимизация МАММП завершены.

Приведем результаты численных экспериментов, используя в качестве тестовой модельную задачу, приведенную в [1]. В качестве характерного значения коэффициента сжимаемости возьмем ß1 = 10_3 / МПа . Остальные данные — такие же, как для тестовой задачи в статье [1]. В таблице приведены результаты численных экспериментов — количества итераций, необходимые для достижения заданной точности е = 10-6 тремя методами: Зейделя (МЗ), верхней релаксации с шахматной упорядоченностью узлов (МВРШУ) и модифицированного адаптивного метода минимальных поправок (МАММП).

Таблица 1

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

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

МЗ МВРШУ МАММП

1 32535 7342 2752

2 20347 3453 1014

3 13138 2387 696

4 12190 2123 688

5 12132 2135 692

Заключение. В статье описано построение и исследование модифицированного адаптивного метода минимальных поправок — МАММП, предназначенного для численного моделирования пространственно-трехмерных процессов фильтрации двухфазной сжимаемой жидкости в пористых средах. Такие процессы являются типичными в задачах моделирования нефтяных месторождений — в случае, когда нефтеводоносный пласт находится под воздействием давления, изменяющегося значительным образом, и необходимо учитывать сжимаемость среды. Для данного класса сеточных задач с несамосопряженным оператором построен МАММП. Проведена его оптимизация и доказана сходимость метода, который обладает асимптотической скоростью сходимости, характерной для «классического» ПТМ, не использующего технику чебышевского ускорения и применяемого для задач с самосопряженным оператором. Результаты численных экспериментов продемонстрировали эффективность метода. Число итераций уменьшилось в 3-20 раз по сравнению с другими часто применяемыми подходами к решению сеточных задач фильтрации (методом Зейделя, методом верхней релаксации и другими).

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

3 1. Сухинов, А. И. Усовершенствованный попеременно-треугольный метод численного решения простран-

§ ственно-трехмерных задач фильтрации двухфазной несжимаемой жидкости / А. И. Сухинов, Л. А. Григорян,

А. А. Сухинов // Вестник Дон. гос. техн. ун-та. — 2016. — Т. 16, № 1. — С. 5-18. 3 2. Коновалов, А. Н. Задачи фильтрации многофазной несжимаемой жидкости / А. Н. Коновалов. — Новоси-

<3 бирск : Наука, 1988. — 166 с.

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

4. Mathematical modeling of flows in porous media / V. R. Dushin [et al] //WSEAS Transactions on fluid mechanics. — 2014. — Vol. 9. — P. 130-166

43

a £ л

5. Sukhinov, A. I. Adaptive modified alternating triangular iterative method for solving grid equations with a non-self-adjoint operator / A. I. Sukhinov, A. E. Chistyakov // Mathematical Models and Computer Simulations. — 2012. — Vol. 4, iss. 4. — P. 398-409.

6. Sukhinov, A. I. Increasing efficiency of alternating triangular method based on improved spectral estimates / A. I. Sukhinov, A. V. Shishenya // Mathematical Models and Computer Simulations. — 2013. — Vol. 5, iss. 3. — P. 257-265.

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

2013. — № 4. — С. 194-202.

References

1. Sukhinov, А!., Grigoryan, L.A., Sukhinov, A.A. Usovershenstvovannyy poperemenno-treugol'nyy metod chislennogo resheniya prostranstvenno-trekhmernykh zadach fil'tratsii dvukhfaznoy neszhimaemoy zhidkosti. [Advanced alternate triangular method of numerical solution of spatial three-dimensional problems of two-phase filtration of incompressible fluid.] Vestnik of DSTU, 2016, vol. 16, no. 1, pp. 5-18 (in Russian).

2. Konovalov, А.К Zadachi fil'tratsii mnogofaznoy neszhimaemoy zhidkosti. [Problems of multiphase filtration of incompressible fluid.] 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 multi-phase filtration problems.] Mathematical Models and Computer Simulations, 2010, vol. 22, no. 4, pp. 118-128 (in Russian).

4. Dushin, V.R., et al. Mathematical modeling of flows in porous media. WSEAS Transactions on fluid mechanics,

2014, vol. 9, pp. 130-166.

5. Sukhinov, A.I., Chistyakov, A.E. Adaptive modified alternating triangular iterative method for solving grid equations with a non-self-adjoint operator. Mathematical Models and Computer Simulations, 2012, vol. 4, iss. 4, pp. 398-409.

6. Sukhinov, A.I., Shishenya, A.V. Increasing efficiency of alternating triangular method based on improved spectral estimates. Mathematical Models and Computer Simulations, 2013, vol. 5, iss. 3, pp. 257-265.

7. Sukhinov, A.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 non-selfadjoint operators.] Izvestiya SFedU. Engineering Sciences, 2013, no. 4, pp. 194-202 (in Russian).

<u

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

<о

Сдана в редакцию 20.06.2016 ч

Запланирована в номер 07.07.2016 ¡^

^

К

ей И

к

к

*

<а н ч

ей К Л

ч <и н к

4 о

к

Е 3 и

ей И К

fe

5

Л

о

-е к К