CC BY
36УДК 622.276:531.71.001.572
МОДЕЛИРОВАНИЕ МНОГОФАЗНОЙ МНОГОКОМПОНЕНТНОЙ ФИЛЬТРАЦИИ ПРИ РАЗРАБОТКЕ МЕСТОРОЖДЕНИЙ УГЛЕВОДОРОДОВ
И.В. Афанаскин, А.В. Королев, С.Г. Вольпин, П.В. Ялов
Федеральное государственное учреждение «Федеральный научный центр Научно-исследовательский институт системных исследований Российской академии наук», ivan@afanaskin.ru
Разработана трехфазная многокомпонетная (композиционная) модель фильтрации для моделирования разработки месторождений жидких и газообразных углеводородов. Приведены некоторые подходы к расчету фазового равновесия углеводородных систем в пласте. Дан пример расчетов по описанной модели.
Ключевые слова: разработка нефтяных месторождений, гидродинамическое моделирование, уравнение состояния, композиционная модель, многокомпонентная модель.
SIMULATION OF MULTICOMPONENT MULTIPHASE FILTRATION FOR HC DEPOSIT DEVELOPMENT
I.V. Afanaskin, A.V. Korolev, S.G. Volpin, P.V. Yalov
System Research Institute, Russian Academy of Sciences, ivan@afanaskin.ru
A three-phase multicomponent (composed) flow simulator for oil and gas fields development has been developed. Some approaches to phase equilibria analysis in hydrocarbon reservoir systems have been shown. An example of analysis with the proposed simulator is presented.
Keywords: oil field development, hydrodynamic simulation, equations of state, composed simulator, multicomponent model.
Введение
В течение последних 50 лет мы наблюдаем значительный прогресс в области математического и компьютерного моделирования процессов разработки нефтяных и нефтегазовых месторождений, базирующихся на научных методах подземной гидродинамики.
Первые серьезные успехи в этой области датируются началом 70-х гг. прошлого века [1, 2], которые сделали возможным создание универсальных 2D и 3D численных моделей, описывающих большинство процессов разработки месторождений углеводородов, в том числе и современных методов повышения углеводородоотдачи пластов.
Важность создания единой углеводородной модели, описывающей процессы подземной гидродинамики, одним из первых осознал известный американский специалист Coats K.H., который сделал первый значительный шаг в этом направлении (1979-1980 гг.) [3, 4]. Отметим также многоцелевой симулятор MARS (автор Kendall и др.) [5].
С математической точки зрения модель фильтрации (подземной гидродинамики) представляет собой систему дифференциальных уравнений в частных производных и алгебраических уравнений (ограничений) с соответствующими начальными и граничными условиями. В изотермической постановке уравнения в частных производных описывают условия материального баланса для сложных многофазных многокомпонентных течений, базирующихся на обобщенном законе Дарси. Могут учитываться также диффузионные и адсорбционные процессы, капиллярные силы, гравитация и т.д.
Алгебраические уравнения описывают термодинамические условия, в частности фазовые равновесия многофазных многокомпонентных систем, в том числе и взаимодействие пластовых флюидов со скелетом пористой среды.
Основное отличие между традиционным и универсальным симуляторами заключается в количестве компонентов, которые учитываются для адекватного описания процессов повышения нефтеотдачи пластов, а также в способе учета условий
термодинамического равновесия. Обычный симулятор базируется на модели «черной нефти» (Black Oil), широко применяемой до сих пор. Термодинамические условия в данном случае определяются только давлением. Это, как известно, позволяет эффективно использовать последовательную процедуру IMPES (неявный учет давления и явный -насыщенности при конечно-разностной аппроксимации дифференциальных уравнений). Также с большой эффективностью применяется и полностью неявный метод (FULLY IMPLICIT). Однако в случае большого числа компонентов и сеточных блоков полностью неявный метод оказывается громоздким и требует огромных затрат машинного времени и оперативной памяти даже на современных ЭВМ, что делает симулятор малоэффективным в таких расчетах.
Поэтому при большом количестве блоков и компонентов применяют последовательные методы IMPECS (неявное давление, явные концентрации компонентов и насыщенности).
В современных широко распространенных в мире симуляторах ECLIPSE и STARS предусмотрена возможность применения целого набора численных методов, в т.ч. IMPES, FULLY IMLICIT и их комбинаций.
1. Система уравнений трехфазной многокомпонентной фильтрации нефти, газа и воды
Композиционные модели фильтрации используются при подробном моделировании залежей, содержащих легкие углеводороды (конденсат и газ) в том случае, когда необходимо тщательно описывать массообмен между фазами, либо когда пластовые флюиды содержат ценные неуглеводородные компоненты. Также такие модели являются основой для изучения множества методов увеличения нефтеотдачи.
Одна из первых изотермических композиционных моделей трехмерной трехфазной фильтрации (1980 г.) была предложена в работе Coats K.H. [3], численное решение системы нелинейных уравнений было осуществлено с помощью полностью неявной (SS) вычислительной схемы (Simultaneous Solution, современное название - Fully Implicit).
Система уравнений представляет собой комбинацию уравнений материального баланса и некоторых ограничений:
j^koSoX +^gSgyt )]+у[%0хгЖ0+ j + q = 0, i = 1,2,..., N, (1)
ñ Í ^Л
] + v[íwWo\ + qw = 0, (2)
ñt'
fO - fg = 0, i = 1,2,..., N, (3)
Wa=-kkra£ (VPa+PagVD), « = o, g, w, (4)
= 1, (5)
i=1 N
I y = 1, (6)
i=1
б0++ ^ = 1, (7)
Р - Р = Р , (8)
о м смо ' \ /
Р - Ро = Рсо, , (9)
Р - Рм = Рсмо + РСО, , (10) здесь N - общее число компонентов в нефтяной и газовой фазах, к(X, Р) и ф(Х, Р) -
соответственно тензор абсолютной проницаемости и пористость пласта, X = (х, у, I) -
точка в трехмерном пространстве, t - время, Ра - фазовое давление, - фазовая
насыщенность, x и y - мольная концентрация i -ого компонента соответственно в
нефтяной и газовой фазе, £а - мольная фазовая плотность, /ла - динамическая вязкость фазы а, ра - массовая плотность фазы а, (Бк) - капиллярное давление в системе нефть-вода, р ) - капиллярное давление в системе нефть-газ, кга - относительная фазовая проницаемость, ч - компонентная плотность источников и стоков, ^ -плотность источников и стоков по воде, 0(х) - превышение точки пласта над некоторой горизонтальной плоскостью, g - ускорение силы тяжести, и - летучести компонентов нефтяной и газовой фазы, V - оператор Гамильтона.
В системе (1)-(10) количество уравнений составляет + 9, что соответствует
^ ^ ^
количеству неизвестных хх,х2,-, хМс; У\,у^-, Уыс; Ро,Р,Р; 8о,^,^; К,Wg• Количество уравнений и неизвестных можно сократить до 2 N + 4, подставляя (4), а затем (8)-(10) в систему (1), (2). За опорное давление принято давление газа Р = Р . Таким образом, в качестве неизвестных функций приняты:
Х1, Х2хшс; У1' у2' Ушс; р; ^, ^, ^ • (11)
Исходная система (1.1)-(1.4) преобразуется к виду:
УМ^Мом: ^р +PogvD)+ + ^УкКХ (vР + PggVD)]+ дг = 0, , = 1,2,..., Мс,
| [ФРАккм(VР ^о ^Р^о дК = 0. (13)
Система уравнений и ограничений (12), (13), (5)-(7) содержит N = 2ЖС + 4 неизвестных функций (11) и описывает математическую модель многокомпонентной трехфазной и, в общем случае, трехмерной фильтрации.
Система уравнений (12), (13) и ограничений (3) нелинейна относительно искомых функций. Она решается с помощью полностью неявных вычислительных схем методом Ньютона. Нелинейные члены уравнений и ограничений линеаризуются на итерациях
N
иМ =1 Рз)+ и , (14)
1 =1 \Ф у
где и нелинейный член, I - номер итерации, р - искомая функция р = ^ Р2 = Х2 ,..., РN = Sw •
Линеаризованная система уравнений с помощью Гауссова исключения редуцируется к системе, содержащей N +1 искомых функций [у3,у4,..., ум ,Р,Б0,^}•
Все вышесказанное описывает общий случай, когда в каждой точке пласта существуют обе углеводородные фазы: нефть и газ. Если в какой-то области пласта ^ = 0, вместо 2N + 4 уравнений надо рассматривать N + 3 уравнений с N + 3
неизвестными функциями {хг}, Р, , ^. Отпадает необходимость в N уравнениях для летучести компонентов и ограничения (6).
В ряде работ [6, 8] уравнения материального баланса для углеводородных фаз представляются в виде:
^[ф 2, (§аБо + )\^[^0х1кктм-1 ^Р-VРCgo )+
+ V + Ч, = 0, , = 1,2,..., Nc,
£ V £ 8 Nc
2 = хЬ + уУ, Ь =-^-, V =-^-, у 2 = 1 (15)
где zi - мольная концентрация г -го компонента в углеводородной системе.
Указанный подход удобен в случае применения вычислительной схемы IMPECS (неявное давление, явные концентрации компонентов и насыщенности). Уравнение для давления получают следующим образом. Суммируют по г уравнения (15)
+ 1 (Vp + р^™)]+ до + д8 - 0.
Далее уравнение (16) суммируется с уравнением (13). Иногда уравнение (13) берется с весом в [6, 8]:
)]-v[aw (vp - vpcgo +
+ ея**гх£ ^Р + (vP Pcgo-Р^о (17)
+ до + д& +вд„ - о.
Уравнение (17) нелинейно относительно давления Р. Поскольку применяемая разностная схема неявная, используются итерационный метод Ньютона. Затем по явной схеме определяются zi. Наконец с помощью методов расчета фазовых равновесий углеводородная смесь при известных Р, ^ расщепляется на нефтяную и газовую фазы и определяются х и у . Следует отметить, что при линеаризации уравнений для давления Р производными по мольным концентрациям пренебрегают [6, 8].
Для моделирования трехфазной многокомпонентной фильтрации была создана программа NCOMP. Использована конечно-разностная схема явная по концентрациям и неявная по давлению. Дискретизация производных по пространству и аппроксимация межблочных проводимостей реализована следующим образом: слабые нелинейности -среднее арифметическое, сильные нелинейности - среднее гармоническое. Относительные фазовые проницаемости отнесены вверх по потоку. Аппроксимация проводимостей по времени реализована явным способом. Аппроксимация слагаемых, учитывающих источники и стоки, реализована явно по концентрациям и неявно по давлению. Для линеаризации использован метод Ньютона. Разностная схема консервативная, условно устойчивая.
2. Уравнение состояния в рамках модели многокомпонентной фильтрации
Для расчетов фазовых равновесий можно использовать различные уравнения состояния. Рассмотрим уравнение состояния Пенга-Робинсона, которое часто используется для моделирования углеводородных систем и хорошо себя зарекомендовало. Оно имеет вид:
г3 + Е2Х2 + ЕхХ + Е0 - 0, (18)
где Е2, Е и Е0 - коэффициенты, X - фактор сверхсжимаемости фазы.
Максимальный положительный корень уравнения (18) г равен коэффициенту сверхсжимаемости газа, а минимальный положительный корень 2о - коэффициенту сверхсжимаемости нефти.
Для расчета коэффициентов уравнения (18) для каждого компонента необходимо задать пластовую температуру Т и следующие константы:
1. Т - критическая температура компонента г;
2. р - критическое давление компонента г;
3. щ - ацентрический фактор компонента г;
4. с. - коэффициенты попарного взаимодействия компонентов г и у .
как:
где
Приведенные давление и температура для каждого компонента рассчитываются
(19)
Р Т
Р = —, Т. = —.
Г1 -р " Г1 г^г
В параметр Оа вычисляется как:
Па(Г,,) = 0,457235529[1 + (0,37464 +1,5422®,. - 0,26992^г2)(1 -Т-1/22)]2. Теперь определим следующие коэффициенты:
Р Р , ч _
4 = ^ (—, ¡)—, В = 0,07796074, = (1 - у4д.
гг г,
Теперь рассчитаем параметры уравнения состояния:
N N «
4=ЕЕ , 5=Ё ,
]=1 к=1 ]=1
Е = ^ - 1,
Е = А - 2Б2 - 2В,
Е0 =-АВ + В2 (В +1). (25)
Коэффициенты летучести определяются с помощью следующих выражений:
(20)
(21)
(22)
(23)
(24)
1п
Рх
V Гхг
= - 1п(10 - В)
+
_А_Г 25+В, щ
2,828427125В V А В
(10 - 0,414213562В Л V ^ + 2,414213562В у
+
+ В Со - ^
(26)
1п
л
V РУг У
1n(Zg - В)-
А
-е - В)+ ,
g 7 2,828427125В V А
25 + В 11п
В
(- 0,414213562ВЛ
7 + 2,414213562В
V ^ ' у
+
(27)
+ В (7 -1)
в\ g ь
«с
5 = £ . (28)
]=1
Напомним, чтобы термодинамическая система пребывала в равновесии, необходимо, чтобы летучести компонентов в жидкой и газовой углеводородных фазах были равны для каждого компонента:
у:0 -= о, г = 1,2,..., .
3. Пример моделирования многофазной многокомпонентной фильтрации
В качестве примера использования рассмотренного подхода приведем моделирование водогазового воздействия на отложения ачимовской свиты одного из месторождений Западной Сибири. Для этого была использована изотермическая композиционная секторная гидродинамическая модель в симуляторе КСОМР. В качестве прототипа принят пласт Ачз месторождения ХХХ. На этом объекте проводятся работы по водогазовому воздействию (закачка в нефтяной пласт воды и углеводородного газа) в районе скважины № 480. В секторной модели представлен участок обращенной девятиточечной системы разработки, включающий скважину № 480, рис. 1. В модель вошли скважины №№ 480, 482, 484, 485, 494, 496, 762, 763, 764. Отметим, что номера скважин изменены для сохранения коммерческой тайны недропользователя.
Рис. 1. Участок работ на карте эффективных нефтенасыщенных толщин
При моделировании рассмотрено влияние вертикальной анизотропии проницаемости на показатели разработки.
Проницаемость и пористость приняты постоянными по площади и переменными по толщине. Геолого-статистический разрез (ГСР) приведен на рис. 2. Он соответствует распределению фильтрационно-емкостных свойств по стволу скважины №480. ГСР характеризуется наличием высокопроницаемого пропластка (240 мД против 20 мД выше по разрезу) в нижней части интервала перфорации.
Рис. 2. Геолого-статистический разрез (ГСР), интервал перфорации 2879-2889 м.
Рассмотрены значения коэффициента анизотропии проницаемости (ак=к/кГ) 0,1; 0,05 и 0,01 д.ед.
Перфорирована верхняя, нефтенасыщенная половина толщины пласта. Нижняя, водонасыщенная половина не вскрыта.
Начальная пластовая температура 92,2°С, начальное пластовое давление 255 ат. Вязкость нефти 0,62 сПз, объемный коэффициент нефти 1,16 ед., плотность нефти 847
кг/м3. Объемный коэффициент воды 1,02 ед. Плотность газа 1,17 кг/м3. Давление насыщения нефти газом 98 ат., растворимость газа в нефти 65 м3/м3. Для моделирования РУТ-равновесия использована модель Пенга-Робинсона. При расчете используются следующие компоненты:
1. N2 - молекулярная масса 28,013,
2. СО2 - молекулярная масса 44,01,
3. С3 - молекулярная масса 44,097,
4. С4 - молекулярная масса 58,124,
5. Н2О
и псевдокомпоненты:
1. С1+С2 - молекулярная масса 17,693,
2. С5+С6+С7 - молекулярная масса 81,739,
3. С8+ - молекулярная масса 209,512.
Нормализованные относительные фазовые проницаемости и капиллярное давление в системе нефть-вода представлена на рис. 3, нормализованные относительные фазовые проницаемости в системе нефть-газ представлены на рис. 4. Капиллярным давлением в системе нефть-газ пренебрегается. Критические точки и значения относительных фазовых проницаемостей в критических точках задавались в виде констант и представлены в табл. 1. Для моделирования относительной фазовой проницаемости по нефти в трехфазном случае использована первая модель Стоуна.
Начальное распределение давления и насыщенности соответствует состоянию на момент начала закачки газа в скважину № 480.
Начальные геологические запасы нефти 1,459 млн. м3.
Все внешние границы модели считаются непроницаемыми.
1,0 0,9 0,8
0,7
Ч 0,6 ш
* 0,5 с
О 0,4
0,3 0,2 0,1 0,0
3.0
2.7
2.4
2.1
1.8
1.5 1,2 0,9 0,6 0,3 0,0
ш
О) С т га Ч
О)
о
X
о. к с с
2
0,0 0,1
0,2 0,3 0,4 0,5 0,6 0,7 Водонасыщенность, д.ед.
0,8 0,9
1,0
ОФП воды
ОФП нефти
Капиллярное давление в системе нефть-вода
Рис. 3. Нормализованные относительные фазовые проницаемости (ОФП) и капиллярное давление в системе нефть-вода
1,0 0,9 0,8
0,7
Ч 0,6 а>
* 0,5 с
О 0,4
0,3 0,2 0,1 0,0
0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0
Газонасыщенность, д.ед.
ОФ П газа ОФ П нефти
Рис. 4. Нормализованные относительные фазовые проницаемости (ОФП)
в системе нефть-газ
Добывающие скважины управляются дебитом жидкости. Нагнетательные скважины управляются закачкой. Дебиты и закачка соответствуют фактическим данным. Скважины в углу модели имеют дебит жидкости 1/4 от фактического. Скважины на ребре модели имеют дебит 1/2 от фактического. Компенсация отбора закачкой на начальный момент времени 100%. Водогазовое отношение 0,005 м3/м3. Закачивается растворенный в нефти газ после сепарации. Состав закачиваемого газа приведен в табл. 2. Начальные мольные концентрации компонентов в углеводородной системе приведены в табл. 3.
Табл. 1
Критические точки и значения относительных фазовых проницаемостей в критических точках_
SWU 1,0 Максимальная водонасыщенность
8т 0,332 Минимальная водонасыщенность
SWCR 0,335 Критическая водонасыщенность
SOWCR 0,236 Критическая нефтенасыщенность в системе с водой
SGU 0,668 Максимальная газонасыщенность
SGL 0,05 Минимальная газонасыщенность
SGCR 0,05 Связанная газонасыщенность
SOGCR 0,220 Критическая нефтенасыщенность в системе с газом
SGWCR 0,02 Критическая газонасыщенность в системе с водой
KRO 0,5 Максимальная ОФП для нефти
KRORW 0,5 ОФП для нефти при критической насыщенности водой
KRORG 0,5 ОФП для нефти при критической насыщенности газом
KRW 0,3 Максимальная ОФП для воды
KRWR 0,07 ОФП для воды при остаточной водонасыщенности
KRG 0,8 Максимальная ОФП для газа
KRGR 0,2157 ОФП для газа при остаточной газонасыщенности
Табл. 2
Состав закачиваемого газа (мольные концентрации)
N2 С02 С1+С2 С3 С4 С5+С6+С7 С8+
Попутный газ 0,0191 0,0 0,9412 0,0122 0,0169 0,0106 0,0
Табл. 3
Начальные мольные концентрации компонентов в углеводородной системе_
N2 СО2 С1+С2 С3 С4 С5+С6+С7 С8+
Углеводородная система 0,0039 0,0053 0,2615 0,0603 0,0679 0,1476 0,4535
Закачивается смесь попутного газа и воды. Расход газа 20000 м3/сут, расход воды 100 м3/сут.
Ограничения на работу добывающих скважин следующие:
1. Минимальный дебит нефти 1 м3/сут.
2. Максимальная обводненность 98 %.
3. Максимальный газовый фактор 1000 м3/м3.
Срок разработки - не более 30 лет.
Всего было рассмотрено 3 варианта в зависимости от коэффициента анизотропии проницаемости, табл. 4.
Вариант 1. Коэффициент анизотропии проницаемости 0,1 д.ед. В результате получен КИН (коэффициент извлечения нефти) 0,310 д.ед. и срок разработки 30 лет. Причина окончания разработки - истечение предельного срока 30 лет. На конец разработки газовый фактор по группе скважин составляет 402 м3/м3, обводненность 62 %, суммарный дебит нефти 14,8 м3/сут, накопленная добыча нефти 452138 м3. Суммарный дебит нефти, накопленная добыча нефти, газовый фактор и обводненность во времени по группе скважин показаны на рис. 5-8.
Вариант 2. Коэффициент анизотропии проницаемости 0,05 д.ед. В результате получен КИН 0,320 д.ед. и срок разработки 30 лет. Причина окончания разработки -истечение предельного срока 30 лет. На конец разработки газовый фактор по группе скважин составляет 320 м3/м3, обводненность 54 %, суммарный дебит нефти 31,1 м3/сут, накопленная добыча нефти 466407 м3. Суммарный дебит нефти, накопленная добыча нефти, газовый фактор и обводненность во времени по группе скважин показаны на рис. 5-8.
Вариант 3. Коэффициент анизотропии проницаемости 0,01 д.ед. В результате получен КИН 0,377 д.ед. и срок разработки 30 лет. Причина окончания разработки -истечение предельного срока 30 лет. На конец разработки газовый фактор по группе скважин составляет 366 м3/м3, обводненность 76 %, суммарный дебит нефти 23,6 м3/сут, накопленная добыча нефти 508356 м3. Суммарный дебит нефти, накопленная добыча нефти, газовый фактор и обводненность во времени по группе скважин показаны на рис. 5-8.
Табл. 4
Основные результаты моделирования
Вар. кг!кг, д.ед. Срок разработки, годы КИН, д.ед. Газовый фактор*, м3/м3 Обводненность*, % Суммарный дебит нефти*, м3/сут
1 0,1 30 0,310 402 62 14,8
2 0,05 30 0,320 320 54 31,1
3 0,01 30 0,377 366 76 23,6
* - на конец разработки.
Резкие изменения параметров на рис. 5 -8 объясняются послойным прорывом закачиваемого газа и воды в скважины, а также отключением некоторых добывающих скважин из-за достижения экономических ограничений. При уменьшении анизотропии проницаемости КИН увеличивается, что косвенно говорит о корректности модели.
1_ 140 -Т------
0 Н------
0 5 10 15 20 25 30
Время, годы
Рис. 5. Суммарный дебит нефти по группе скважин
3 550000 500000
ш Я
£ 450000
§ 400000
№
Я
§5 350000
ш
О 300000
я
250000 200000 150000 100000 50000 0
0 5 10 15 20 25 30
Время, годы
Рис. 6. Накопленная добыча нефти по группе скважин
а
о н X ез
■е
«
ва о
РО
ез
и-
1000 900 800 700 600 500 400 300 200 100 0
10
15
20
25 30
Время, годы
Рис. 7. Газовый фактор по группе скважин
.а н и о Я Я о Я
п о ва »о О
100 90 80 70 60 50 40 30 20 10
-2 -3
--- '—
\ \
V/л
10
15
20
25 30
Время, годы
Рис. 8. Обводненность по группе скважин
Заключение
Разработана трехфазная многокомпонетная (композиционная) модель фильтрации для моделирования разработки месторождений жидких и газообразных углеводородов. Приведены некоторые подходы к расчету фазового равновесия углеводородных систем в пласте. Описана численная схема. Создана альфа-версия программы для моделирования под названием NCOMP. В качестве примера использования рассмотренного подхода
0
5
0
0
5
приведены результаты моделирования водогазового воздействия на отложения ачимовской свиты одного из месторождений Западной Сибири. Получено косвенное подтверждение корректности модели
Работа выполнена при поддержке гранта РФФИ 16-29-15105 офи_м.
Литература
1. Peaceman D.W. Fundamentals of Numerical Reservoir Simulation. Amsterdam - Oxford -New York: Elsevier Scientific Publishing Co., 1977. 176 p.
2. Aziz K, Settari A. Petroleum Reservoir Simulation. London: Applied Science Publishers Ltd. 1979.
3. Coats K.H. An Equation of State Composition Model // Soc. Pet. Eng. J. (Oct. 1980) P. 363376.
4. Coats K.H. Reservoir Simulation: A General Model Formulation and Associated Physical // Numerical Sources of Instability, Boundary and Interior Layers. Computational and Asymptotic Methods / ed. J.J.H. Miller. Dublin: Boole Press (198o) P. 62-76.
5. Kendal R.P. Development of a Multiple Application Reservoir Simulator for Use on a Vector Computer // SPE Middle East Oil Technical Conference and Exhibition. Bahrain: Manama, March 14-17. 1983.
6. Ngiem L.X., Fond D.K., Aziz K. Compositional Modeling with an Equation of State // Soc. Pet. Eng. J. (Dec. 1981) P. 687-698.
7. Coats K.H., Thomas L.K., Pierson R.G. Compositional and Black Oil Reservoir Simulation // SPE Reservoir Evaluation & Engineering (August 1998) P. 372-379.
8. Kazemi H., Vestal C.R., Shank G.D. An Efficient Multicomponent Numerical Simulator // Soc. Pet. Eng. J. (Oct. 1978) P. 355-368.
