CC BY
18
УДК 532.546
БОТ: 10.15587/2312-8372.2018.148386
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ЗАДАЧИ НЕЛИНЕЙНОЙ ТРЕХФАЗНОЙ ФИЛЬТРАЦИИ
Гасымов С. Ю., Мамедов Р. С., Керимова С. Р.
1. Введение
Как известно [1], с проблемой многофазной фильтрации связаны важные практические задачи нефтегазодобычи, гидротехники, химической технологии и т. д. Наиболее актуальной является эта проблема для теории и практики добычи нефти и газа.
Задачи многофазной фильтрации обладают определенными специфическими особенностями, которые не дают возможность всегда применять при их численном решении ставший классическим метод конечных разностей. Поэтому возникает необходимость в разработке разностных схем на адаптивных сетках с учетом особенностей решения [2, 3].
Адаптивные сетки уменьшают искусственную вязкость и осцилляцию численного решения. Если даже вычислительная сетка имеет минимальные точки, они опять же дают удовлетворительные результаты по количеству и качеству для всей области (даже при включении зон, обладающих особенностями решения).
Поэтому актуальным является исследование процессов совместной фильтрации нефти, газа и воды в пористой среде с учетом:
- реальных геолого-промысловых данных (упругие свойства и неоднородность пласта);
- свойств фильтрующихся жидкостей и газа (сжимаемость, зависимость вязкостей от давления);
- относительных фазовых проницаемостей и капиллярных сил.
2. Объект исследования и его технологический аудит
Объектом исследования является численное моделирование процесса фильтрации нефти, газа и воды на адаптивных сетках с учетом некоторых свойств жидкостей во время их совместного течения.
Задачи многофазной фильтрации обладают специфическими особенностями, которые не дают возможность всегда применять при их численном решении ставший классическим метод конечных разностей.
Одним из наиболее проблемных мест в теории многофазной фильтрации является то,что шаги по пространственной переменной должны измельчаться в областях резкого изменения не только градиента водонасыщенности, но и градиента газонасыщенности. Это объясняется тем, что из-за очень малой вязкости свободный газ под действием градиента давления обгоняет остальные компоненты смеси - воду и нефть. Поэтому возникает необходимость в разработке разностных схем на адаптивных сетках, который дает возможность учитывать особенности решения.
3. Цель и задачи исследования
Цель исследования - разработка эффективных и достаточно универсальных численных методов решения задачи одномерной трехфазной фильтрации.
Для достижения поставленной цели необходимо выполнить такие задачи:
1. Построить математическую модель процесса одномерной трехфазной фильтрации в пористой среде.
2. Построить адаптивную сетку, учитывающую измельчение шагов в областях резкого изменения не только градиента водонасыщенности, но и градиента газонасыщенности.
3. Разработать алгоритм разностно-итерационного метода в подвижных сетках.
4. На основании численных экспериментов обсудить вопросы вычислительной реализации предложенной методики и дать практические рекомендации по ее применению.
4. Исследование существующих решений проблемы
При проведении гидродинамических расчетов, связанных с проведением и эксплуатацией нефтегазовых месторождений при водонапорном режиме, в основном используется двухфазная модель Рапопорта-Лиса. Однако, из-за того, что во многих месторождениях происходит трехфазное течение, возникает необходимость рассмотреть модель трехфазного течения и оценить влияние газовой фазы на процесс вытеснения. А это, в свою очередь, приводит к необходимости разработать численный метод решения задач трехфазной фильтрации, описываемых системой нелинейных дифференциальных уравнений в частных производных. Для исследования этого процесса в основном используются численные методы [4, 5].
В работах [6, 7] используется модель Баклея-Леверетта, а в работе [4] учитывается также влияния капиллярных сил. И для численного решения полученных задач применяется идея метода [8]. Однако этот метод не учитывает особенности многофазной фильтрации, и в ряде случаев дает качественно и количественно неверные результаты [4].
Среди основных направлений решения этой проблемы, выявленных в ресурсах мировой научной периодики, могут быть выделены работы [9-11]. Однако в работе [9] вытеснение предполагается поршневым, и получены приближенно-аналитические решения. В работе [10] предлагается численный метод, идейно отличающийся от метода конечных разностей и конечных элементов. Необходимо отметить, что модификации этих работ широко используются многими авторами [1 1].
В настояшей работе предлагается эффективный разностно-итерационный метод в подвижных сетках, обладаюший свойством адаптируемости к особенностям задач и отличающийся высокой точностью [2].
5. Методы исследования
5.1. Математическая постановка задачи
Допустим, рассматривается круговой горизонтальный пласт радиуса Я и мощностью Н. Предполагается, что пласт не пропускает жидкость снизу и
сверху, в центре пласта находится гидродинамическая совершенная эксплуатационная скважина радиусом г=гс, на контуре г=Я расположена батарея скважин, закачивающих в пласт сжимаемую жидкость.
Предположим, что в рассматриваемом пласте имеется трехфазная сжимаемая жидкость, фазы несмешивающиеся и вначальный момент времени 1=0 находятся в состоянии капиллярного равновесия. На основании модели Рапопорта-Лиса для жидкостей с приведенными условиями изотермическая плоско-радиальная нелинейная задача фильтрации описывается системой дифференциальных уравнений в частных производных в безразмерном виде:
и д г — д (г,з) Р1 дг 01
и д г — д с)т) ГР^ПР^Б 2,5:1) - дг _
-1
а
дг
(г^)еПТ={г<гс<Я , 0<t<T}.
Здесь использовались следующие обозначения: - коэффициент подвижности фаз:
\ Ра)
к(г) - проницаемость; г - полярная координата; ? - время. Функция:
(1)
- ' \%га(1Р\>С О, ^гиЛР\<С„
учитывает вязкопластичное свойство нефти [1, 12]. Здесь О} - начальный градиент давления.
К системе уравнений (1) необходимо добавить уравнения состояния жидкостей и газа:
Ра=РАР«)> « = 1,2,3,
(2)
<
и соотношения:
^+52+53=1.
(3)
Отметим, что вообще в задачах многофазной фильтрации давления Ра(а=1,2,з) не равны друг другу и различаются на величину соответствующего капиллярного давления, т. е.:
Р\ РГ1 = Рк 1 ( ' ' )'
где ры - капиллярное давление между нефтью и водой, а рк2 - капиллярное давление между газом и нефтью.
Зная рк1 и рк2, можно при помощи них выразить капиллярное давление между газом и водой:
Рз-Р2=Ри+Р/<2-
Необходимо отметить, что в задачах трехфазной фильтрации функции /«(51'52'5з)> (а = 1>2,3) обладают следующими свойствами [6, 7]:
У ) ^ = 0, ^ ,
/2 > ' ) = — '
Уд > ' ) = 53 ~ 53'
где значения
называются остаточными насыщенностями.
а = 1,3
V У
Предположим, что неизвестные в системе (1) функции:
/ п
Ра(Г>Щ а = 1>3
Определим их из начальных и краевых условий.
Предположим, что в момент времени t=0, т. е. до начала момента эксплуатации, значения определяемых функций известны, т. е.:
ра(г,0) = Р°Аг),г 1г]К,а = 1,2Л (5)
Для рассматриваемой задачи краевые условия могут быть в следующем виде.
Для эксплуатационной скважины (т. е. при г=гс) предполагается, что из пласта добывается только нефть, и в этот момент скачок капиллярного давления не учитывается:
др1 _ др2
дг дг
дРз = Ф2 дг дг
г = г,0^<Т.
На внешнем контуре (т. е. при г=Я) предполагается, что в пласт закачивается только сжимаемая жидкость (вода) и соответственно течения первой и третьей фаз нет:
дг
Р2(г,г) = д>2(г)
Фз
дг
= 0, г = Я, О<Ь<Т.
(7)
<
Таким образом, согласно условиям (5)-(7), задачу можно поставить следующим образом. Необходимо найти такие функции:
ра(г,£)(а = 1,2,3),
чтобы они удовлетворяли бы системе уравнений (1) и начальным краевым условиям (5)- (7).
Предполагая существование и единственности решения задачи (1)-(7), дадим приближенный метод ее решения.
5.2. Численный метод решения задачи
Прежде всего, отметим одну специфическую особенность задач многофазной фильтрации. Установлено, что в процессе всего исследования, обычно вокруг скважин происходит резкое изменение фазовых давлений. Поэтому при моделировании процессов многофазной нелинейной фильтрации краевые условия необходимо аппроксимировать с высокой точностью, т. е. вокруг скважин шаг сетки необходимо умельчать [2, 3]. Таким образом, покроем область:
Q = {(r,t):rc<r<R,0<t<T}
неравномерной сеткой соьг =с0кхс0т следующим образом. Узлы неравномерной сетки со/г по г определяются так:
ъ = п_± + к , I = 1, М - 1; г0 = тс , гм = Я ,
^¿+1 = ,1 — 1,2,..., ¿о 2,..., ¿0 + п ; /11 = г0 ,
- _ | 2Л*0_1 ' £г > ,
" 10,5(1! - г0 - Е^о £г < 2/ч0_1
4+1 =
/¿¿0 , £г < ,
^¿0+2 — п» ' ^¿о+п+1 + 1 — ^¿о+п—¿+2 » ^ — 1'2, ..., ¿0 + п + 1
Здесь
'1 ,1 = 1,2,...,п, ¿о + 2,..., ¿о + п,
п' 1 2 ,1 = п + 1, п + 2,.., ¿о — 2 ¿о = шах{1 : г0 + ЕЙ кк + <
£г = Ь± — Г0 — £¿0 кк — 2/1;о_2 ,
где Ж - количество разбиений отрезка [г0, ц ];
п п п
п- = го - количество постоянных шагов; п - количество разбиений отрезка [¿р^].
Для более лучшей сходимости и точности результатов при решении рассматриваемой задачи итерационно-разностным методом вычисления целесообразно прово -пь с ^озра.хающим по времени шагом. Таким образом, вычисления начинаются с очень малым шагом по ? и постепенно
п
увеличиваются до заданного Т [2]. В соответствии с этим узлы сетки (От определяются следующим образом:
т ={к- тп-1 Л<п< nfc_1 , П Ьтах > Пк <П< N .
Здесь
а А = о,
и п
к = к_
п к
щ = шж[п '.егптп < тпшх},
где \п\ - целая часть; кс - целое положительное число.
Таким образом, увеличение шага по ? происходит через каждый кс
о
временной слой. При вычислениях брали кс=100 и гшч=4-10" .
Используя неявную консервативную схему, дифференциальную задачу (1)-
(7) аппроксимируем на неравномерной сетке о)ь следующей дискретной задачей:
1
А
(гр^щ)
Рим -Ри
,+-> ¡1 2 "7+1
гт
г дР2 , Г
Ун У12 "'"У 13
дЬ ш о1
\ГР№\\.
Фз
Ри ~Ри-\ к
1+1
1
А
Р 2./+1 Р'1,1
1——-1
2 П;
■{гр2Ъ)г
Ри -Р2,1-1
*7+1
/— 2
/г,-
/+1
= Г-171
г др^.г др^.г Фз
У21 У22 У 23
¿ж д1
А
)/+1-Е--(г АА)/_!-7-
/г;
/+1
к
/+1
= Г;Ш
", , { Фз
Уз1 """У 32 ~|~Узз от д/ д!:
Ра М) = Л ,п = 0 ,0 <1< М ,а = 1,2,3 ,
(8) (9)
С1- =
V (^11 - ) = ¿Г1« - й1)' * = о,о < «Ж
(10)
Ь~Члп+1 -2/г+1 "1 = 0
=ср?\ 1 = М ,Ъ<п<Ы, (10)
/тЧ»'г+1 - п"+] "1 = 0
пм м гш-1) Здесь
А и = [ А - 52 - 53 ) - А52 + А53 ],> /21 = А52> /13 = - А52>
/21 = Р2^2 >./21 = Р2^2 ~ Р2^2 1/21 = ^ ' /31 = _А53 >./з2 = ^ 'У32 = А53 ~~ Р353 '
Ра=^ ,ОС = 1,2,3 ;52 = 52(ри) ,53 = 53(р/г2), ф«
<> = -— Ч
Ф/,1 Ф/,2 Фи Фи
др
Аппроксимируя производные —- в праюй части системы (8) следующим образом:
др рп+1-рп-= а = 123
Ы т
и проведя некоторые элементарные преобразования, приведем (8) к следующему виду:
Г1/*:1 (гМ^П/Ри-! - + (грЛпХ.\
1+1 /+1
V 1 + 1
+
+гтт
У",у; + -г/у/г ;/;,Х:/
= -?)тт~] (/и.рЧ. -/¡зХ/)>
+
¿"V - (г РоА Г*+(ГР'А> )•
)/+1
ч/+1 /-1/
+
+г1тт~%2. ~]р£ + +Г%+\ (гр,Л2 р!2% -
Г V(гр.л,р13% -Г1 (г/?:,Я,+ к +г,тт >£ + (гр,Я3-
-г^г-Уз^КУЙ1 = (/змК/ +/ззХ/)-
/+1
+
(12)
Таким образом, заменив р"^ = У7«./ ,// ( = К, ( , введем нижеследующие матрицы и векторы:
В0 = "-1 0 0" 1 -1 0 II '0 0 0" 1 -1 0 п 0 = -(р! 0 1
0 1 -1 0 0
'А 0 0 " 'К ' 0 0"
А = 0 Л 0 и = п ¥г = ,С = 0 с2 0 ?
0 0 ^3 Рз 0 0 С3_
В =
А 0 0" '10 0" '"
0 К 0 А -= о 0 0 >Вм = 0 1 0
0 0 Ь3 0 0 1 0 0 1
$ О
Система (12) нелинейная, определяемые функции входят в коэффициенты системы. Для того чтобы линеаризовать нелинейные коэффициенты, используем метод простой итерации, т. е. в каждый новый момент времени эти коэффициенты считаются известными от предыдущего момента времени [13]:
1//-1
Т1+\г,р)*Т1(г,р).
Принимая во внимание вышесказанное, разностную задачу (12) можно записать в следующем трехдиагональном каноническом виде:
„(')$('+!) ,,(')$('+!) п С)
-Во У о +Со У1 =-,Ро ,/ = 0,
„(')$('+!) ,?(0$(/+1) „(/)$(/+!) „О
Аг Ум -5, У,- + С/ У/+1 = -,Р/ ,1ЩМ-1,
(13)
Лм Уми -Вы Ум
-^М, г =М.
Элементы матриц Л, в и С, а также вектора имеют следующий вид:
Аг =
л..=
3,г
Р
1,г
/7(0 _
2,/ ~~
/7(0 _
3,/ ~~
а
и
-Чп^{1иУРи +/12X1-/13X1);
-^"1(/21,Х:+/22Х,0;
п
-1
Кх {гРАп 1XIV + К1 (гРЛп1)% + г!тт~1/п,1
у2
-1/ С)
ки =ггтти/г
Ь^; =
'13,/
¿от , = ^
22,/
^ (гаЛ +К1 (гр-А +г<пп~Аш
ь2и=^К,=-г,тт~ Ъ и-
^ , =
32,/
ь..
33,/
г1 0.
(г аЛ +(г АЛ )!!!/+'умА/,,.
<
Таким образом, задавая начальные приближения в виде (9), чтобы найти из системы (13) ¡+1-ое приближение, применяем матричный вариант метода
прогонки (алгоритма Томаса) [13, 14]. В каждый момент времени итерационный процесс продолжаем до такого значения 10, чтобы условия:
max
О <i<M
yv о 1,/
max
0<i<M
max
0 <i<M
+1) у{1о) 1li 1,
У (¿0+1) у(к) 12,( 12,(
3,/ 3,/
<s.
3 ,
удовлетворялись бы одновременно. Здес s1, s2, s3 - предварительно заданные точности для нахождения /-ой компоненты приближенного решения.
6. Результаты исследования
Предложенный выше итерационно-разностный метод был применен для численного моделирования одного процесса трехфазной фильтрации. Для проведения расчетов были использованы следующие начальные данные:
Я = 100 м, Я = 10 м, гс = 0.1 м; ш = 0.2; k = 10~12 м2;
цл = 0.3 пуаз; /у2 = 0.01 пуаз; ¿и3 = 0.0013 пуаз;
срх (t) = 139.5 атм; ср2 {£) = 140 атм; р\ (г) = 140 атм;
S\ (г) = 0.25; S" (г) = 0.45; 53° (г) = 0.3; G7=100 дин/см3;
Rq = 100 м; k0 = 10"12 м2; р{) = 0.01 пуаз; Т0 = 108 сек;
р0 = 140 атм; р0 = 1 г/см
Здесь R0, k0, Т0,р0, р0 - характеристические величины измерения, используемые для приведения первоначальной задачи к безразмерному виду. Уравнения состояния фаз даются в следующем виде [15]:
A (ft) = 0.000853а + 0.82592, р2 ( р2) = 0.01035а + 0.99989, р3(р3) = 0.063р3.
До начала процесса насыщенности фаз составляли соответственно £1=0.45, 52=0.25, £?=0.30. Зависимости относительных проницаемостей и капиллярных сил от насыщенностей брались, следуя [15].
Проведенные численные эксперименты показывают, что в начале процесса воздействия в пласте начинает формироваться увеличенная зона нефтенасыщенности - нефтяной вал, размеры которого увеличиваются и который, вытесняя газ, движется по направлению к эксплуатационной скважине. Вытеснение нефти же происходит практически после нефтяного
вала. Под «шириной» и «приростом» профиля «нефтяного вала» понимаются следующие величины:
/., =те5^г :ге[гс,Щ ,51(г)>510},
кг = тах 51 (г) - ,
ге[гс,В].
На рис. 1 показаны результаты расчетов при (р2 -ср{ =0.5 кгс./см2 в различные моменты времени для различных значений вязкости нефти с использованием адаптивной сетки. Здесь сплошной линией показаны графики, отражающие состояние нефтяного вала в различные моменты времени, а штриховой линией обозначены графики, показывающие водонасыщенность.
Необходимо отметить следующую особенность рассматриваемого процесса. При приближении нефтяного вала к эксплуатационной скважине из пласта выходит только газ. Из рис. 1 видно, что по мере уменьшения вязкости нефти уменьшается время приближения нефтяного вала к эксплуатационной скважине.
Рис. 1. Состояние нефтяного вала в различные моменты времени
1а хис. 2 показано влияние вязкости нефти на размеры и «прирост» нефтяного вала. Здесь сплошной линией показан график нефтенасыщенности и
штриховой линией график газонасыщенности в момент времени 1=3 месяца (ф—
2
ф1=0.7 кгс/см ).
Б *
0.55
0.5
0.45
0.4
0.35
0.3
0.25 0.2
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Г
Рис. 2.Влияние вязкости нефти на размеры и «прирост» нефтяного вала
На рис. 2 начальное распределение нефте- и газонасыщенности показано точками. Кривые, отмеченные цифрами 1, 2, 3, 4, 5, 6 являются результатами проведенных вычислений, соответственно при /^ = 0,01; 0,02; 0,03; 0,1; 0,2; 0,4 пуаз.
Как видно из рис. 2, при увеличении вязкости нефти длина ¡у и «прирост» Ну нефтяного вала уменьшаются, причем уменьшение ¡у по сравнению с Ну происходит с большей скоростью. Например, из табл. 1 видно, что при увеличении значения ¡л1 в 2, 10 и 40 раз ¡у уменьшается примерно на 41 %, 73 %, 86 %, а Ну, соответственно, уменьшается на 5 %, 36 %, 65 %.
Таблица 1
Влияние вязкости нефти на размеры нефтяного вала
9 Я 9 t=5 месяцев; ф2-ф1=0.7кгс/см , к=210' см , =0,2
пуаз 0.01 0.02 0.03 0.1 0.2 0.4
¡у(м) 60.01 3з.02 25.29 14.49 11.21 9.40
Ну 0.0079 0.0739 0.0701 0.0501 0.0427 0.0269
Таким образом, исходя из сказанного выше, при наличии фазы свободного газа геометрические размеры нефтяного вала характеризуют полноту вытеснения нефти, так что увеличение размеров вала является результатом вытеснения нефти водой.
Теперь рассмотрим влияние скорости вытеснения на процесс трехфазного вытеснения (табл. 2).
Таблица 2
Влияние скорости вытеснения на процесс фильтрации
Пф 0.5 кгс/см 2 1 кгс/см 2
r Sj S3 Sj S2 Г S "1
0.1 0.4500 0.2500 0.3000 0.4511 0.2500 0.2980
0.2 0.4500 0.2500 0.3000 0.4750 0.2503 0.2747
0.3 0.4500 0.2500 0.3000 0.5446 0.2506 0.2048
0.4 0.4525 0.2502 0.2974 0.5446 0.2506 0.2048
0.5 0.5199 0.2505 0.2296 0.5320 0.2632 0.2048
0.6 0.5370 0.2560 0.2070 0.3936 0.4016 0.2048
0.7 0.4329 0.3624 0.2048 0.3197 0.4755 0.2048
0.8 0.3617 0.4336 0.2048 0.2579 0.5373 0.2048
0.9 0.3139 0.4813 0.2048 0.2150 0.5802 0.2048
1.0 0.2726 0.5226 0.2048 0.1804 0.6148 0.2048
Результаты вычислений показывают, что при увеличении скорости фильтрации геометрические размеры и «прирост» нефтяного вала резко увеличиваются. Как видно из табл. 2, при увеличении разности давлений до 0.4 кгс/см длина вала увеличивается на 35.4 %, а прирост на 9.6 %.
7. SWOT-анализ результатов исследований
Strengths. Предложенный алгоритм можно использовать для разработки определения оптимального режима и времени процессов теплопроводности в неоднородных средах.
Weaknesses. К слабым сторонам исследования относится трудоемкость вычисления.
Opportunities. Используя полученные исследования, можно повысить нефтеотдачу месторождения.
Threats. Для внедрения полученных исследований необходимо приобрести дополнительное оборудование, что означает денежные затраты.
8. Выводы
1. Построена математическая модель процесса одномерной трехфазной фильтрации в пористой среде, которая более адекватно описывает рассматриваемый процесс. Особенностью данной модели является то, что она учитывает многие из вышеупомянутых свойств фильтрующихся жидкостей и газа, например, сжимаемость, вязкопластичность, зависимость вязкостей от давления, относительные фазовые проницаемости и капиллярные силы (модель Рапопорта-Лиса).
2. Дан алгоритм построения адаптивной сетки, учитывающий измельчение шагов в областях резкого изменения решения, и который дает возможность более точно аппроксимировать производные искомых функций.
3. Разработан алгоритм разностно-итерационного метода в подвижных сетках, который учитывает особенности решения и производит автоматическое сгущение в областях резкого изменения решения.
4. На основании численных экспериментов установлено, что в начале процесса воздействия в пласте начинает формироваться нефтяной вал, размеры которого увеличиваются и который, вытесняя газ, движется по направлению к эксплуатационной скважине. А вытеснение нефти происходит практически после нефтяного вала.
Литература
1. Мирзаджанзаде А. Х., Ковалев А. Г., Зайцев Ю. В. Особенности эксплуатации месторождений аномальных нефтей. М.: Недра, 1972. 200 с.
2. Пирмамедов В. Г. Об одном разностно-итерационном методе в подвижных сетках решения некоторых нелинейных задач теории фильтрации и теплопроводности. Деп в ВИНИТИ, 11.07.1975, № 2027-75.
3. Касумов С. Ю., Сулейманов С. Г., Нифтиев Я. М. О применении адаптивной сетки к одной радиальной задачи трехфазной фильтрации // Изв. АН Аз.Р, сер. физ. тех и мат.наук. 1996. Т. 139, № 1-2.
4. Peery J. H., Herron E. H. Three-Phase Reservoir Simulation // Journal of Petroleum Technology. 1969. Vol. 21, Issue 2. P. 211-220. doi: http://doi.org/10.2118/2015-pa
5. Sonier F., Ombret O. A Numerical Model of Multiphase Flow Around a Well // Society of Petroleum Engineers Journal. 1973. Vol. 13, Issue 6. P. 311-320. doi: http://doi.org/10.2118/3627-pa
6. Клевченя А. А., Таранчук В. Б. Численное моделирование процесса неустойчивого вытеснения неньютоновской нефти // Динамика многофазных сред. 1981. С. 193-198.
7. Шалимов Б. В. Численное моделирование одномерной трехфазной фильтрации // Изв. АН СССР, МЖГ. 1975. № 26. С. 59-66.
8. Douglas J. Jr., Peaceman D. W., Rachford H. H. A method for calculating multidimensional immiscible displacement // Trans. SPE of AIME. 1959. Vol. 216. P. 297-308.
9. Pascal H. Dynamics of moving interface in porous media for power law fluids with yield stress // International Journal of Engineering Science. 1984. Vol. 22, Issue 5. P. 577-590. doi: http://doi.org/10.1016/0020-7225(84)90059-4
10. Elnaggar H., Karadi G., Krizek R. J. Effect of non-darcian behavior on the characteristics of transient flow // Journal of Hydrology. 1971. Vol. 13. P. 127-138. doi: http://doi.org/10.1016/0022-1694(71 )90210-1
11. Shutler N. D. Numerical, Three-Phase Simulation of the Linear Steamflood Process // Society of Petroleum Engineers Journal. 1969. Vol. 9, Issue 2. P. 232-246. doi: http://doi.org/10.2118/2233-pa
12. Лапин В. Об исследовании некоторых нелинейных задач теории фильтрации // ЖВМ и МФ. 1979. Т. 19, № 3. С. 689-700.
13. Abou-Kassem J. H., Farouq-Ali S. M. Petroleum Reservoir Simulations. Gulf Publishing Company, 2006. 445 p.
14. Самарский А. А. Теория разностных схем. М.: Наука, 1984. 656 с.
15. Killough J. E. Reservoir Simulation With History-Dependent Saturation Functions // Society of Petroleum Engineers Journal. 1976. Vol. 16, Issue 1. P. 3748. doi: http://doi.org/10.2118/5106-pa
