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

Вычислительные технологии

Том 7, № 5, 2002

ЧИСЛЕННЫЙ АНАЛИЗ НЕКОТОРЫХ МЕТОДОВ СОПРЯЖЕНИЯ ДВУХ МОДЕЛЕЙ ФИЛЬТРАЦИИ НЕСМЕШИВАЮЩИХСЯ ЖИДКОСТЕЙ

О.Б. БОЧАРОВ Институт водных и экологических проблем СО РАН Новосибирск, Россия e-mail: bob@ad-sbras.nsc.ru И. Г. ТЕЛЕГИН Горно-Алтайский государственный университет, Россия e-mail: dckt@bcenter.gorny.ru

By means of numerical methods the comparison of three approaches to the solution of the conjugate problem of two generally accepted models of filtration of immiscible fluids in porous media is fulfilled (Muskat — Leverett and Backley — Leverett models). The analysis of the solutions is performed. The a posteriori estimation of inaccuracy and convergence properties are provided.

В работе численно исследуются проблемы сопряжения двух наиболее употребляемых в компьютерных технологиях нефтедобычи моделей фильтрации двухфазной жидкости. В обводненной части нефтяного пласта, где подвижная нефть почти вытеснена, капиллярные силы оказывают слабое воздействие на процесс фильтрации двухфазной жидкости, и в этой области целесообразно использовать модель Баклея — Леверетта (БЛ-модель) с общим давлением для обеих фаз. В малообводненной части нефтяного пласта роль капиллярных сил при вытеснении нефти водой оказывается весьма существенной, что приводит к необходимости использования более сложной модели Маскета — Леверетта (МЛ-модель) с различными фазовыми давлениями. Однако, как показано А. Н. Коноваловым [1], за счет обращения в нуль функций относительных фазовых проницаемостей естественные граничные условия для этой модели являются плохо обусловленными (градиенты решения в окрестности границы становятся бесконечно большими). В работе О. Б. Бочарова [2] предложено в качестве граничного условия на эксплуатационной скважине рассмотреть уравнение модели Баклея — Леверетта, что решает эту проблему. В. Н. Монаховым в работе [3] БЛ-модель предложено применять в окрестности эксплуатационной скважины, тогда на скважине не требуется граничного условия. В итоге возникает задача сопряжения МЛ-и БЛ-моделей. В работе анализируются три подхода к решению этой задачи.

Уравнения моделей. Одномерная МЛ-модель фильтрации двухфазной жидкости в однородной изотропной пористой среде в отсутствие массовых сил при заданном суммарном расходе фаз имеет вид [4]

ms-t = (koao(s)ai(s)sx - Q(t)b(s))w, (1)

© О. Б. Бочаров, И. Г. Телегин, 2002.

где x £ [0,L] — пространственная координата, L — расстояние от нагнетательной скважины до эксплуатационной; t — время; s = (si — S0)/(1 — S0 — S0) — динамическая насыщенность смачивающей фазы, s1 — истинная насыщенность смачивающей фазы, (S0, S0) = const — остаточные водо- и нефтенасыщенности; m = m0(1 — S0 — S0), m0 — пористость коллектора; k0 = const — абсолютная проницаемость коллектора; a0(s) = a1(s) = —pc(s)/(^2( + ^ fc2)), pc(s) = (m0/k0)1/2Yj(s) — капиллярное давление, y — коэффициент поверхностного натяжения, j(s) — функция Леверетта; b(s) = fc1/( + ^ fc2) — коэффициент подвижности вытесняющей фазы, fc^ = fci(s) — относительные фазовые проницаемости, ^ = ^ — вязкости фаз (индекс i = 1 соответствует воде, i = 2 — нефти); Q(t) — расход воды на нагнетательной скважине. Здесь и далее для функций fa = df/da.

При этом функциональные параметры модели имеют следующие свойства [4]:

1) 0_< (fc1, fo) < 1 при s £ (0, l) ao(0) = ao(1) = M0) = fc2(1) = 0;

2) ( fc2s,pcs) £ C(G), где G — замкнутая область в пространстве переменных (x, s), причем M-1 < 0ub|u2, —pcs) < M, M > 0;

3) |ln(a1),b/fc1| < Mo(M).

К особенностям уравнения (1) следует отнести:

— вырождение типа при s = 0, s = 1,

— относительно малый коэффициент при s^ обусловливает появление внутренних пограничных слоев (областей больших градиентов насыщенности).

Полагая Q(t) = Q0 = const, введем безразмерные переменные x = x/L, t = tQ0/(mL), при этом уравнение (1) запишется в виде

st = (ea(s)s^ — b(s))x= — v*, (2)

где e = 7(m0fc0)1/2/(Q0L^2) — капиллярное число; a(s) = — a0js/(+ ^ fc2); v = v(x,t) — скорость фильтрации вытесняющей фазы. Значению параметра e = 0 соответствует модель Баклея — Леверетта. Обозначим через Г = {x,t|x = l,t > 0} линии, при l = 0 — это будет обозначение для нагнетательной скважины, l = 1 — эксплуатационной, l = (0,1) — линии в промежутке (0,1). Через Птга обозначим область {x, t|m < x < n, t > 0}, m £ [0, l), n £ (l, 1].

Сопряжение МЛ- и БЛ-моделей. Рассмотрим случай, когда вблизи эксплуатационной скважины капиллярные силы малы. Разобьем область П01 на две, в области П0г фильтрация двухфазной жидкости описывается МЛ-моделью, и соответственно этому в (2) e = 0, и в области 0,ц — БЛ-моделью, и в (2) e = 0. Предлагается для решения уравнения (2) в области П01 изучать следующую начально-краевую задачу:

s|x=0 = 1, s|t=0 = s0(x), x £ [0, 1], [s]x=1 = [v]x=i = 0, (3)

где [f]x= = {f(l — 0,t) — f(l + 0,t)} — скачок функции f(x,t) на линии Г^.

Из условий склейки на Г следует условие easx|x=i-0 = 0, и задача (2), (3) распадается на две подзадачи. В области

st = (easx — b)x, s|t=0 = s0(x), s|x=0 = 1, easx|x= = 0. (4)

В области фильтрация описывается БЛ-моделью и в соответствии с этим e = 0:

st = —bx, s|t=0 = s0(x), s|x=i = Sz(t), (5)

где функция Si(t) = s(x, t)|x=i определяется по решению s(x,t) задачи (4).

Задача (4), (5) состоит в сопряжении решений в = в(х, ¿), (х,£) £ П0г уравнения (4) и решений в(х,£), (х, ¿) £ Пц уравнения (5), эволюционных по переменной ¿. Задачу (4), (5) в дальнейшем будем называть в-склейкой.

В работах [3, 5] предложено сменить эволюционную переменную в области Пц и искомую насыщенность в на расход V. В этой ситуации легче обосновать непрерывность расхода. Тогда задача (5) принимает вид

Vx = VI *=о = ^(х), v|x=г = У (¿), (6)

где в = — функция, обратная зависимости V = Ь(в), зд(х) = Ь[в0(х)], Уг(¿) = Ь[$г(¿)].

Особенность задач (4), (6) заключается в склейке решений в = в(х,£), (х,£) £ Пог уравнения (4), эволюционного по переменной ¿, и функции в(х, ¿) = <^[и(х, ¿)], (х,¿) £ Пц, выражающейся через решения v(x,t) уравнения (6), эволюционного по переменной х.

Задачу (4), (6) условимся обозначать далее v-склейкой.

Вопросам существования обобщенного решения задачи сопряжения методом v-склейки посвящена работа [3], а первые предварительные расчеты проведены в [5]. Необходимо отметить также, что в [3, 5] изучалась параболически регуляризованная задача для уравнения (6).

Построение алгоритмов численного решения задачи (4), (6) затрудненно несколькими обстоятельствами:

1. Сменой эволюционной переменной, что приводит к необходимости усложнения алгоритма получения численного решения.

2. ^|а=о;1 = го, что затрудняет применение градиентных итерационных методов при линеаризации разностного уравнения и ухудшает сходимость численного метода.

3. Так как искомой функцией в (6) является v, для нахождения в необходимо обращение функции Ь в уравнении v = Ь(в), что в случае сложных функций /с1(в) и &2(в) — трудоемкая задача.

4. В настоящее время нахождение решения задачи (4), (6) для многомерного случая представляется авторам трудноразрешимой задачей.

В работе [6] проведен сравнительный анализ некоторых численных методов решения задачи v-склейки и выбрана оптимальная разностная схема для решения (6) в Пц:

И+1 — 1п+1 — „га Г.П+1 — пП

-Г--+ (1 - -Г- = --(1 - -, (7)

п п т т

где П — шаг по оси х; т = гП — шаг по ¿; — веса.

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

п2

Vx + = ^((2а1 - !)(^2V*)* + г(2^2 - 1)<Р«Vtt). (8)

Так как |^=0;1 = то, уравнение (8) является неравномерно параболическим. Также очевидно, что вблизи точек V = 0 и V = 1 разностное уравнение (7) не аппроксимирует уравнение (6) при П > 0, т > 0. С учетом всех этих обстоятельств в данной работе рассматривается в-склейка.

Перейдем к описанию алгоритма решения задачи сопряжения (4), (5).

Введем сетку в области П01 с распределенными узлами = {хг = гП, Ьп = пт, г = 0,Ж,п = 0,1, 2...}, П = 1/Ж — шаг по пространственной координате, т = КП2 —

шаг по временной переменной. Пусть при этом линии Г соответствует номер Мг = /Ж узла на сетке (предполагается, что М^ не является дробным числом).

Для нахождения численного решения использовались противопотоковые схемы [7] из соображений необходимости вычислять з(ж,£) на линиях склейки без привлечения дополнительных узлов сетки, а также из необходимости соблюдения балансовых соотношений. Шаг к брался равным 0.005 (Ж = 200), а т = 0.00025.

Уравнение МЛ-модели в Пог аппроксимировалось с помощью неявной разностной схемы первого порядка

„п+1 _ „га с г г _ ° I га „п+1 „га

к («п+ 2 „п+1 _ а- 2 0 _ 1, г = 1,м, _ 1, (9)

„И+1 _ „га 2с

М1 _ __пп „га+1 _ тга+1 га+1 _ п _ 1

т = к ММ -1 „х,мг %,Мг, „0 = „0 = 1 Для численного решения уравнения (5) использовалась следующая разностная схема:

„п+1 _ „п _

г т г = 1 _ (1 _ а)^, г = М + 1,М1; „М+1 = ЗГ+1, (10)

где 6П+1 линеаризовались как 6П+1 = ЬП + ЬПг(„П+1 _ „П). Схема (10) при а = 0.5 имеет погрешность аппроксимации 0(к + т2), при других значениях а £ [0,1] погрешность аппроксимации 0(к + т). Для согласования порядков аппроксимации (9) и (10) схема (10) использовалась с а = 1.

Для численного решения систем применялся метод правой прогонки. На каждом временном шаге вычисляли основные характеристики процесса вытеснения: положение жс(£) фронтовой водонасыщенности в БЛ-модели вс, которая определяется решением нелинейного уравнения

) = Ь(зс)/зс

с помощью метода деления пополам, п(£) — обводненность пласта

1

0

Интеграл в правой части уравнения вычислялся по формуле трапеций. Также контролировалась Ж/ (¿) — предельная точка распространения фронта водонасыщенности. В численных расчетах использовался следующий набор параметров:

кх = 52, к2 = (1 _ в)2, з = (1 _ в)/(& + в), & = 0.9,

с = 0.5, в0(ж) = 0, / = 0.9, ^ = 0.1.

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

Характерными особенностями решения задачи сопряжения являются подъем и выпола-живание графика в точке х = I, обусловливаемое наличием в этой точке краевого условия еазХ|Х=1 = 0 (рис. 1). Следующей особенностью решений задачи сопряжения является эффект недобегания, который выражается в том, что в формируется характерный для БЛ-модели скачкообразный профиль водонасыщенности, и поэтому приход воды в задаче сопряжения на эксплуатационную скважину происходит позднее, чем в задаче для МЛ-мо-дели во всей области. Следствием эффектов недобегания и выполаживания является тот факт, что после момента прихода воды на эксплуатационную скважину, рассчитанного по МЛ-модели, обводненность пласта п(£) в задаче сопряжения выше (максимальный дисбаланс по обводненности между решениями до прихода воды на эксплуатационную скважину составил 0.14%, а после прихода — 1.20%). На рис. 2 приведены графики решений з-склейки и решений по МЛ-модели при капиллярном запирании, в этом случае после прорыва воды дисбаланс обводненностей между моделями составил 5.78 %.

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

— при I = 0.99 дисбаланс по обводненности между решениями при £ = 1 составил 0.05 %, а максимальный дисбаланс был равен 0.14 %. Максимальная разница между решениями в норме С(0,1) составила 0.05, а при £ = 1 была порядка 0.002. Очевидно, что при таком I различия в гидродинамических характеристиках незначительны и проявляются только в модуле градиента водонасыщенности на эксплуатационной скважине;

— при I = 0.95 максимальный дисбаланс по обводненности между решениями до прихода воды по МЛ-модели (£ = 0.2855, время прорыва) составил 0.07%, а до £ = 1 максимальный дисбаланс был равен 0.65 %. Максимальная разница между решениями была равна 0.131, а при £ = 1 составила 0.008. Время безводной нефтеотдачи по сравнению с МЛ-моделью больше на 0.9 %;

— при I = 0.9 максимальный дисбаланс по обводненности между решениями до прихода воды составил 0.12%, а до £ = 1 был равен 1.20%. Максимальная разница между решениями составила 0.16, а при £ = 1 была равна 0.015. Время безводной нефтеотдачи по сравнению с МЛ-моделью больше на 4.9 %;

— при I = 0.8 максимальный дисбаланс по обводненности между решениями до прихода воды был равен 0.16%, а до £ = 1 составил 2.17%. Максимальная разница между

Рис. 1.

Рис. 2.

решениями была равна 0.191, при Ь = 1 она составила 0.025. Время безводной нефтеотдачи по сравнению с МЛ-моделью больше на 9.5 %. С уменьшением / < 0.95 эффекты недобегания и выполаживания в окрестности х = / усиливаются.

Анализ движения фронтовой водонасыщенности показывает, что график хс(Ь) для модели сопряжения состоит из трех участков (рис. 3):

— первый участок совпадает с графиком хс(Ь) для МЛ-модели в О01;

— второй участок характеризуется ускорением продвижения хс(Ь)(он догоняет Ж/ (Ь)). Это объясняется тем, что в течение некоторого времени формируется резкий профиль водонасыщенности, характерный для БЛ-модели. При этом происходит торможение х/(Ь), так как перестают работать капиллярные силы;

— третий участок — характерное для БЛ-модели равномерное движение хс (Ь).

На рис. 4 приведены графики движения фронта водонасыщенности х/(Ь). Видно, что график х/(Ь) в задаче сопряжения состоит из двух участков, которые соответствуют наличию капиллярных сил и их отсутствию.

0 0.2 0.4 0.6 0.8 г 0 0.2 0.4 0.6 0.8 I

Рис. 3. Рис. 4.

Б. Использование сквозного счета.

Для устранения выполаживания решения в окрестности х = / и уменьшения эффекта отставания фронта был применен другой подход к решению задачи сопряжения. Вместо схем (9), (10) соответственно в областях и использовалась единая разностная схема для уравнения (2) во всей области О01, но с разрывной функцией с(х) (метод сквозного счета)

„п+1 _ „п 1 , . _

г г сг+1 а", ! „Л1 _ £,•_ 1 аГ ! ^ _ Ьп+1, г = 1, N _ 1,

г г ( п + 1 ^^ п + 1

т = к V г+1 аг+1 ^ _ Сг_ 1 аг_ I

„п+1 _ „п 2с „ 1

% = _ N_I п „п+1 _ ьп+1 „п+1 = „п = 1

т = к _1 „х^ , „о = „о = 1

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

С е, х £ [0,/ _ Л/],

с(х) =\ е^ху, х £ (/_ Л/,/), (11)

[ 0, х £ [/, 1],

где А1 — ширина переходной зоны.

На рис. 5 приведены графики решений для случаев А и Б. Жирными линиями обозначены графики, соответствующие Б с а = 1, А1 = 0.06, а тонкими — графики для случая А. Видно, что при пересечении линии Г/ различия между решениями несущественны, но у графиков для случая Б нет выполаживания, хотя фронт водонасыщенности движется более медленно. При больших £ разница между решениями практически исчезает. Различия гидродинамических характеристик случаев А и Б невелики, только хс(£) на втором участке движется более плавно и тонкая линия на рис. 6 идет очень близко к жирной. Максимальный дисбаланс до прорыва был равен 0.13%, а при £ = 1 составил 0.14%. В варианте Б время безводной нефтеотдачи по сравнению со случаем А больше на 1.2%.

Были проведены многовариантные расчеты с различными значениями А1, а. Установлено, что уменьшение А1 приводит в пределе к случаю, близкому к А, а увеличение А1 растягивает область формирования характерного для БЛ-модели скачкообразного профиля водонасыщенности, при этом нефизичная полочка случая А исчезает. Была выбрана оптимальная ширина переходной зоны А1 £ (0.04, 0.06), что составило 8-12 узлов сетки.

Было также установлено, что а влияет на интенсивность формирования баклеевского фронта. Это выражается в том, что при а > 1 фронт формируется более интенсивно, а при а < 1 — медленнее и случай Б приобретает большее сходство с А.

0 0.2 0.4 0.6 0.8 х о 0.2 0.4 0.6 0.8 X

Рис. 5. Рис. 6.

Сравнение решений в-склейки в случае Б с А1 = 0.06 и решения по МЛ-модели во всей области показали, что дисбаланс обводненности до прихода воды на эксплуатационную скважину по МЛ-модели был равен ^10-8 %, а максимальный зарегистрированный дисбаланс составил 1.29%.

На рис. 6 приведены графики решений для случая Б (жирные линии) в сравнении с решениями для задачи ^-склейки (тонкие линии) с параметрами а1 = а2 = 1. Максимальный дисбаланс обводненности составил 1.25%. Так как у ^-склейки коэффициент при V = 0, V = 1 стремится к бесконечности, возникает искусственнное растяжение фронта и увеличивается дисбаланс. С учетом этого методы в-склейки предпочтительнее.

Если сравнивать случаи А и Б в-склейки, то очевидно, что случай Б более физичен, более прост в реализации, хотя и менее экономичен.

Сопряжение БЛ- и МЛ-моделей. Рассмотрим случай, когда вблизи нагнетательной

скважины капиллярные силы пренебрежимо малы. Предположение вполне оправданно, так как в окрестности нагнетательной скважины реализуется подобие задачи о точечном источнике, и следовательно, скорости здесь достаточно велики, а локальное капиллярное число очень мало, в результате этим слагаемым в уравнении (2) можно пренебречь.

Однако численные эксперименты убедительно показали, что постановка задачи сопряжения БЛ- и МЛ-моделей в том же виде, как и (4), (5), неудовлетворительна, так как дисбаланс обводненности на начальном участке мог превышать 50 %. Это объясняется тем, что в начальные моменты времени пропитка вносит существенный вклад в процесс вытеснения. Поэтому постановка задачи сопряжения БЛ- и МЛ-моделей отличается от постановки задачи склейки решений МЛ и БЛ тем, что в начальные моменты времени для нахождения решения з(х,Ь) применяется МЛ-модель во всей области П01 и только по прошествии определенного времени используется модель сопряжения БЛ, МЛ.

В соответствии с этим разобьем зону фильтрации линией Гг на две подобласти. Пусть в П0г процесс фильтрации описывается БЛ-моделью

„ = _ Ьх, з|4=о = Зо(х), з|х=о = 1, х £ (0,/(Ь)], (12)

а в О,ц вытеснение осуществляется в соответствии с МЛ-моделью

„ = (еазх _Ь)х, з|4=о = ^(х), з|ж=г = (Ь), еа5ж|ж=1 = 0, х £ (/, 1], (13)

вд = {0; '<*'

Таким образом, значение / будет зависеть от времени. В общем случае / может быть функцией от з|х=г или п.

Задача сопряжения решений „ = з(х,Ь) уравнений (12), (13) аналогична задаче сопряжения уравнений (4), (5).

В численных расчетах использовались схемы (9), (10) с теми же функциональными параметрами. Были проведены многовариантные расчеты с варьированием / и Ьг. Необходимо отметить, что при уменьшении Ьг и увеличении / дисбаланс обводненности увеличивается и различие решений по МЛ-модели и задаче сопряжения БЛ-, МЛ-моделей также увеличивается.

На рис. 7 приведены жирными линиями графики решения по задаче сопряжения (12), (13) с / = 0.2, Ьг = 0.1 и тонкими — решения по МЛ-модели на всей области П01(пунктир, как и ранее, линия сопряжения Гг). До прорыва воды на эксплуатационную скважину максимальный дисбаланс составил 0.27%, а при Ь =1 он был равен 1.45%. Дисбаланс уменьшается с увеличением d в функции Леверетта 3 („).

Максимальная разница между решениями в норме С(0,1) составила 0.08. Остальные гидродинамические характеристики процесса вытеснения по МЛ-модели и по моделям сопряжения БЛ, МЛ отличаются весьма незначительно.

Из рис. 7 видно, что структуры решений по модели сопряжения (12), (13) и МЛ-модели во всей области П01 в окрестности нагнетательной скважины различаются, это объясняется действием капиллярных сил. Однако если взять экспериментальные функциональные параметры моделей из [8]:

„3 т П ,3 ^ 0.16(1 _ „)

Т = „3, Т2 = (1_ „)3, 3 („)

(1.1 _ 0.8„)2(0.19 + 1.44„ _ 0.64„2):

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

Сравнение МЛ-модели и задачи сопряжения (12), (13) показывает, что можно заметно сократить время расчетов, при этом существенно не ухудшая баланса и точности.

0.8

0.6

0.4

0.2

л

i = 0 \

0.2

0.4

0.6

0.8

Рис. 7.

Рис. 8.

Апостериорная оценка порядков сходимостей. В силу нелинейности уравнений (4)-(6) вопрос о теоретическом исследовании порядка сходимости разностной задачи является весьма сложным. Использование при аппроксимации разностей против потока в конвективном члене приводит в общем случае к первому порядку. Исследуем вопрос о погрешности метода по апостериорному правилу Рунге: на сетке с шагом h погрешность R = max (|s(h) — s(h/2)|^ /(2p — 1), где p — порядок аппроксимации (равный 1 в нашем

случае). Задача (4)-(6) была просчитана на последовательности сеток с шагами h, h/2, h/4, h/8, h/16 (h = 1/N0, N0 = 10). Поведение погрешности в зависимости от величины N показано на рис. 8.

На основе решений Sh/2 (k = 0,1, 2, 3) численно получен порядок сходимости раз-

ГЛ D I (h/2k) (h/2k+1)|

ностной задачи по следующей методике. Определим величины Rk = max |si — s2i |

i

(k = 0,1, 2, 3). Очевидно,что Rk является максимумом разности на двух последовательных сетках. Предполагая, что Rk = C(h/2k)p, k > 0, C и p не зависят от k, вычислим величины pk = lg2(Rk/Rk+i), k = 0,1, 2. Это дает несколько оценок порядка классической сходимости для задач сопряжения. Чтобы объединить их в одну, возьмем среднее значение pi = I Pk I /3. Расчеты, проведенные по данной методике, дали для всех трех Vk=o J

рассмотренных методов на момент времени t = 0.35 значение p = 0.68 при h = 0.1 и т = 0.0005.

Список литературы

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

[2] Бочаров О. Б. О задаче с сосредоточенной емкостью для одномерных уравнений двухфазной фильтрации // Динамика сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ИГиЛ. 1985. Вып. 73. С. 149-155.

[3] ЖумАгулов Б. Т., Монахов В.Н. Гидродинамика нефтедобычи. Алматы: Казгос-ИНТИ, 2001. 336 с.

[4] АнтонцЕв С. Н., КАжихов А. В., МонАхов В. Н. Краевые задачи механики неоднородных жидкостей. Новосибирск: Наука, Сиб. отд-ние, 1983. 316 с.

[5] Телегин И. Г. Численная реализация сопряжения основных моделей фильтрации двухфазной жидкости // Динамика сплошной среды: Сб. науч. тр. / РАН. Сиб. отд-ние. ИГиЛ. 2000. Вып. 116. С. 107-111.

[6] Бочаров О. Б., Телегин И. Г. Экспериментальное исследование некоторых разностных схем при сопряжении различных моделей фильтрации двухфазной жидкости // Вычисл. технологии. 2002. Т. 7, №1. С. 34-40.

[7] Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971.

[8] КурБАнов А. К., КурАнов И. Ф. Влияние смачиваемости на процесс вытеснения нефти водой. М.: Гостоптехиздат, 1963. ВНИИ, НТС по добыче нефти, №24.

Поступила в редакцию 20 марта 2002 г.