CC BY
26
Вычислительные технологии
Том 7, № 1, 2002
ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ НЕКОТОРЫХ РАЗНОСТНЫХ СХЕМ ПРИ СОПРЯЖЕНИИ РАЗЛИЧНЫХ МОДЕЛЕЙ ФИЛЬТРАЦИИ ДВУХФАЗНОЙ ЖИДКОСТИ *
О.Б. БОЧАРОВ Институт водных и экологических проблем СО РАН Новосибирск, Россия e-mail: bob@ad-sbras.nsc.ru И. Г. ТЕЛЕГИН Горно-Алтайский государственный университет, Россия e-mail: dckt@bcentr.gorny.ru
Numerical analysis of difference schemes for solving the problem on conjugation of different two-phase fluid filtration models is made. As a result of numerical experiments on different grids the optimal scheme is chosen.
При изучении вытеснения нефти водой из пористых коллекторов с учетом капиллярных сил наиболее распространенной является модель Маскета — Леверетта (МЛ-модель). Однако, как показано А. Н. Коноваловым в [1], вследствие обращения в нуль функций, описывающих относительные фазовые проницаемости, естественные граничные условия для этой модели являются плохо обусловленными(градиенты решения в окрестности границы становятся бесконечно большими). В работе О. Б. Бочарова [2] предложено в качестве граничного условия на эксплуатационной скважине рассмотреть уравнение модели Ба-клея — Леверетта (БЛ-модель). В работах В. Н. Монахова, И. Г. Телегина [3, 4] эту модель предложено применять в окрестности эксплуатационной скважины. В итоге возникает задача сопряжения обеих моделей. Целью данной работы является сравнительный анализ некоторых численных методов решения данной задачи сопряжения.
Уравнения моделей. Одномерная МЛ-модель фильтрации двухфазной жидкости в однородной изотропной пористой среде в отсутствие массовых сил при заданном суммарном расходе фаз имеет вид [5]
msi = (k°a°(s)ai(s)sx — Q(t)b(s))w, (1)
где x £ [0,L] — пространственная координата; L — расстояние от нагнетательной скважины до эксплуатационной; t — время; s = (si — S°)/(1 — S° — S°) — динамическая насыщенность смачивающей фазы; si — истинная насыщенность смачивающей фазы;
* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-00622.
© О. Б. Бочаров, И. Г. Телегин, 2002.
) = const — остаточные водо- и нефтенасыщенности; m = m0(1 — S0 — S0), m0 — пористость коллектора; k0 = const — абсолютная проницаемость коллектора; a0(s) = fci k2, ai(s) = —pcs/(^2(k1 + ^k2)), pc(s) = (m0/fc0)1/27j(s) — капиллярное давление; 7 — коэффициент поверхностного натяжения; j(s) — функция Леверетта; b(s) = k1/(k1 + — коэффициент подвижности вытесняющей фазы; k = kj(s) — относительные фазовые проницаемости; ^ — вязкости фаз (индекс i =1 соответствует воде, i = 2 — нефти). Здесь и далее для функций fa = df/da.
При этом функциональные параметры модели имеют следующие свойства [5]:
1) 0 < (k1,k2) < 1 при sG(0,1), ao(0) = a0(l) = k1(0) = k2(1) = 0;
2) (fc1s,fc2s,pcs)gC(G), где G — замкнутая область в пространстве переменных (x, s), причем M—'Pcs)<M, M > 0;
3) | ln(a1),b/fc1|<Mo(M).
Полагая Q(t) = Q0 = const, введем безразмерные переменные: x = x/L, t = tQ0/(mL), при этом уравнение (1) запишется в виде
st = (ea(s)sx — b(s))x= — (2)
где e = y(m0fc0)1/2/(Q0L^2) — капиллярное число; a(s) = — a0js/(k1 + ^fc2); v = v(x,t) — скорость фильтрации вытесняющей фазы.
Значению параметра e = 0 соответствует модель Баклея — Леверетта. Обозначим через Г0 = {x,t|x = 0,t>0} нагнетательную скважину, через Гг = {x,t|x = М>0} — линию сопряжения моделей, через Г1 = {x,t|x = 1,t>0} — эксплуатационную скважину, а через Птга — область {x, t|m < x < n, t > 0}, m = 0; /, n = /; 1.
Сопряжение решений МЛ и БЛ-моделей с использованием смены эволюционной переменной. Рассмотрим случай, когда вблизи эксплуатационной скважины капиллярные силы малы. Разобьем область П01 на две, в области П0г фильтрация двухфазной жидкости описывается МЛ-моделью и соответственно этому в (2) e=0, а в области 0,ц — БЛ-моделью и в (2) e = 0. Предлагается для решения уравнения (2) в области П01 изучать следующую начально-краевую задачу:
s|x=0 = 1; s|t=0 = S0(x), x G [0,1]; [s]x= = [v]x= = 0, (3)
где [f]x= = {f(/ — 0,t) — f(/ + 0,t)} — скачок функции f(x,t) на линии Гг.
Из условий склейки на Гг следует условие easx|х=г-0 = 0, и задача (2), (3) распадается на две подзадачи. В области П0г
st = (easx — 6)л; s|t=0 = S0(x), s|x=0 = 1, easx|x= = 0. (4)
В области 0,ц фильтрация описывается БЛ-моделью, в соответствии с этим e = 0:
St = —bx; s|t=0 = S0(x), s|x=i = Sz(t). (5)
В работах [3, 4] предложено сменить эволюционную переменную в области 0,ц и искомую насыщенность s на расход v. В этом случае легче обосновать непрерывность расхода. Тогда задача (5) принимает вид
Vx = —^t(v); v|t=0 = v0(x), v|x=i = vi(t), (6)
где зависимость s = <^(v) является обратной по отношению к зависимости v = b(s), v0 = b[s0(x)], v(t) = b[S^(t)], функция S^(t) = s(x,t)|x= определяется по решению s(x,t) задачи (4).
Особенность задач (4), (6) заключается в склейке решений в = в(ж,£), (ж,£) € П0г уравнения (4), эволюционного по переменной ¿, и функции в(ж, ¿) = ^[^(х, ¿)], (ж, ¿) € Пц, выражающейся через решения ^(ж,^) уравнения (6), эволюционного по переменной ж.
Вопросам существования обобщенного решения задачи сопряжения МЛ- и БЛ-моделей (4), (6) посвящена работа [3], а первые предварительные расчеты проведены в [4]. Необходимо отметить, что в [3, 4] изучалась параболически регуляризованная задача для уравнения (6).
Построение алгоритмов численного решения задачи (4), (6) затруднено несколькими обстоятельствами:
1) заменой эволюционной переменной;
2) условием ^|^=0;1 = го, которое осложняет применение градиентных итерационных методов при линеаризации разностного уравнения.
Перейдем к описанию алгоритма решения задачи сопряжения (4), (6). Вначале произведем регуляризацию используемой в дальнейшем функции ^
( ^(в*), в€[0, в*],
^(в) = \ ^(в) [ ^(1 - в*), в€[1 - в*, 1].
Значение в* подбиралось по результатам численных экспериментов, черта над ^ (в) в дальнейшем изложении опускается.
Введем основную сетку в области П01 с распределенными узлами ш^т = (ж = гЛ, ¿п = пт; г = 0, У, п = 0,1, 2...} (Л = 1/У — шаг по пространственной координате, т = гЛ — шаг по временной переменной). Пусть при этом линии Гг соответствует номер Мг = /У узла на сетке (предполагается, что Мг не является дробным числом).
Также введем вспомогательную сетку в области Пг1 ш^ Тг, — (жг — гЛ^, tn = пт^; г = Мг, М1, п = 0,1, 2...} (Л^ = Г1Т^ — шаг по координате ж, т,и — шаг по времени), М1 соответствует линии Г1 или эксплуатационной скважине.
При этом г и г1 выбирались так, чтобы на линии Гг узлы сетки ш^^ совпали с узлами сетки ш^т.
Для нахождения численного решения использовались противопотоковые схемы из соображений необходимости вычислять в(ж,£) на линии склейки без привлечения дополнительных узлов сетки. В дальнейшем используются обозначения, принятые в [6].
Уравнение МЛ-модели в П0г аппроксимировалось с помощью неявной разностной схемы первого порядка
„п+1 _ „га р 'Ч___пп еп+1 „»
Л
г+ 1 г-1
^ _ ЬП++1, г =1,Мг _ 1, (7)
„п+1 _ „п 9 в м, в м 2е
Мг °Мг = _ ¿Ь „п+1 _ ьп+1 „п+1 = п =1
т = Л "мг -1 вх,мг ЬХ,Мг, в0 = в0 = 1-Для численного решения уравнения (5) использовалась разностная схема
„п+1 „п
I — в _
г L = 1 _ (1 _ ^п,г, г = Мг + 1,М1, вМ+1 = 57+1, (8)
т
где Ьп+1 линеаризовывались как Ьп+1 = Ьп + Ьпг(„п+1 _ „п). Схема (8) при а = 0.5 имеет погрешность аппроксимации О (Л + т2), при других значениях а € [0,1] погрешность аппроксимации О(Л + т).
Для численного решения систем (7), (8) применялся метод правой прогонки. Для аппроксимации уравнения (6) использовались два семейства разностных схем
_ „.п+1 ..га _ ..га ?.п+1 _ „.га _ „га а1-7--+ (1 _ -7- = _а2^г--(1 _ а2^г -' (9)
/¿V т,и
г = М ,М1 _ 1; <1 = 6(^Г+1);
1п+1 _ 1п+1 31п+1 _ + 7,п-1 31п+1 _ + 7,«-1
^г+1 1 _ ^ ,га+131г+1 41г+1 + ^г+1 п + 1 ПП^ -7- = _ а3^г -2--(1 _ -2-' ( )
» ¿V 2 / V 2 / V
г = М ,М1_ 1; 1М+1 = 6(^П+1)-
Схема (9) при а1 = а2 = 0.5 имеет погрешность аппроксимации 0(72 + т2), при а1 = 0.5 и 0"2 = 0.5 погрешность аппроксимации 0(7^ + т2), при а1 = 0.5 и а2 = 0.5 — 0(7^ + ), в остальных случаях погрешность аппроксимации 0(7^ + ).
Схема (10) при а3 = 0.5 аппроксимирует исходное уравнение со вторым порядком по обеим переменным, при а3 = 0.5 погрешность аппроксимации 0(7^ + т2). Схема (10) является трехслойной, поэтому на первом шаге, т. е. когда вода подходит к линии склейки, применялась двухслойная схема (9) с а1 = 1 и а2 = а3.
В численных расчетах использовался следующий набор параметров:
= в2; = (1 _з)2; 3 = (1 _в)/(0.9 + в); I = 0, 8; е = 0.5; ^ = 0.1; 5* = 0.05.
Выбор функциональных параметров был обусловлен простотой вычисления обратной функции ^(г>) = 6(-1)(^) для уравнения (6).
Ниже на рисунках жирными кривыми обозначены решения или характеристики, относящиеся к задачам сопряжения, а светлыми — результаты контрольного счета для БЛ-модели с а = 0.5 во всей области П01, пунктиром — линия Г.
Результаты расчетов. Первоначально тестировали разностные схемы на сопряжение смены эволюционной переменной, т. е. во всей области решали задачу для БЛ-модели без линии сопряжения и с ней.
Рассмотрим случай применения единой сетки ш^т на всей области. В численных расчетах использовались шаги сетки 7 = 7^ = 0.01 (Ж = 100), т = т^ = 0.001, М^ = 80, М1 = 100. Были проведены многовариантные расчеты при различных а, а1, а2, а3, з0(ж).
Так как на нагнетательной скважине задан расход, а жидкости принимались несжимаемыми, то можно вести контроль решения по водному балансу.
Сравнение численного решения задачи (4) - (6) с решением по БЛ-модели во всей области показывает, что после прихода воды на линию Г при малых з0(ж) имеет место дисбаланс порядка 2.3-2.4%, дисбаланс уменьшается с увеличением з0(ж) и /. Так, при 50(х) = 0.05 дисбаланс составил 0.6-0.7%. При з0(ж) = 0 фронт водонасыщенности по модели сопряжения движется быстрее. Для схемы (10) с а3 = 1 и I = 0.8 это приводит к уменьшению безводной нефтеотдачи на 6%. Еще одной особенностью решений задачи сопряжения является то, что при уменьшении а1, а2, а3 в схемах (9), (10) возникают осцилляции перед фронтом водонасыщенности. Они имеют место уже при а1 = 0.5, а2 = а3 = 0.72. С уменьшением з0(ж) осцилляции появляются даже при больших весах. Наилучшие результаты получились при решении задачи по схеме (10) при а3 = 0.74 и 50(х) = 0.05 (рис. 1).
Опыт использования разностных схем для задачи (4) - (6) на основной сетке показал, что оптимальной является более простая схема (9), так как, варьируя а1 и а2, можно
получить решение, близкое к тестовому. В то же время эта схема несколько хуже сохраняет баланс, нежели схема (10). Однако этот дисбаланс ложится в рамки допустимых пределов. График функции обводненности
1
„№ = ш>% /
0
для худшего варианта, когда дисбаланс достигает 2.45%, приведен на рис. 2. Интеграл в правой части уравнения вычислялся по формуле трапеций.
Теперь рассмотрим случай использования в области Пц вспомогательной сетки Тг,. В численных расчетах использовались шаги сетки Л = = 0.01, т = = 0.001 (Ж = 100), Мг = 80, Мх = 280.
Сравнение численного решения задачи (4) - (6) с решением по БЛ-модели во всей области показывает, что дисбаланс несколько меньше, чем при расчете задачи (4) - (6) на основной сетке, при в0(ж) =0 он составляет 1.8-2.6%, а при в0(ж) = 0.05 — порядка 0.29-1.0%.
0 0.2 0.4 0.6 0.8 х
Рис. 1.
Рис. 2.
О 0.2 0.4 0.6 0.8 X
Рис. 3.
Рис. 4.
При использовании схемы (9) с а < 0.57, а2 = 0.5 и в0(ж) = 0.05 имели место осцилляции после фронта водонасыщенности, при этом ухудшался баланс. С уменьшением в0(ж) осцилляции появляются даже при больших весах. Опять имело место опережающее движение фронта при в0(ж) = 0, как и на основной сетке. Расчет по схеме (10) при а3 = 1 и в0(ж) = 0.05 приводил к решению без осцилляций, но с нефизичной полочкой в окрестности фронтовой насыщенности. Наилучший результат получился при расчете по схеме (9)
при а1 = 0.58 и а2 = 0.5, з0(ж) = 0.05 (рис. 3).
Опыт использования вышеописанных разностных схем для задачи (4) - (6) на вспомогательной сетке показал, что оптимальной является схема (9).
На рис. 4 приведены результаты расчетов задачи (4)-(6) с е = 0.5, сравнение проводилось с решением по МЛ-модели во всей области П01. Расчет в области Пц проводился по схеме (9) при а1 = 0.58 и а2 = 0.5 на вспомогательной сетке. Вследствие условия еазх|х=г = 0 слева от Г образуется полочка. Капиллярные силы в МЛ-модели растягивают фронт водонасыщенности, поэтому время безводной нефтеотдачи по модели сопряжения больше, чем по МЛ-модели. Дисбаланс до прорыва воды составлял 0.6 %, а наибольший дисбаланс возникает после прорыва воды по модели сопряжения 1.75% при £ ^ 0.27. Такое расхождение вполне приемлемо, тем более, что с уменьшением е отличия в решениях уменьшаются, а в пластовых условиях обычно е < 0.1.
Список литературы
[1] Коновалов А. Н. Задачи фильтрации многофазной несжимаемой жидкости. Новосибирск: Наука. Сиб. отд-ние, 1988. 166 с.
[2] Бочаров О. Б. О задаче со сосредоточенной емкостью для одномерных уравнений двухфазной фильтраци // Динамика сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Ин-т гидродинамики. 1985. Вып. 73. С. 149-155.
[3] МонАхов В. Н. Сопряжение основных математических моделей фильтрации двухфазной жидкости // Мат. моделирование. 2001 (в печати).
[4] Телегин И. Г. Численная реализация сопряжения основных моделей фильтрации двухфазной жидкости // Динамика сплошной среды: Сб. науч. тр. / РАН. Сиб. отд-ние. ИГиЛ. 2000. Вып. 116. С. 107-111.
[5] АнтонцЕв С. Н., КАжихов А. В., МонАхов В. Н. Краевые задачи механики неоднородных жидкостей. Новосибирск: Наука. Сиб. отд-ние, 1983. 316 с.
[6] Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971.
Поступила в редакцию 24 октября 2001 г.
