Сравнительный анализ некоторых разностных схем для задач двухфазной фильтрации без учета капиллярных сил Текст научной статьи по специальности «Физика»

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

Том 8, № 4, 2003

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

О.Б. БОЧАРОВ Институт водных и экологических проблем СО РАН,

Новосибирск, Россия И. Г. Телегин Горно-Алтайский государственный университет, Россия e-mail: bob@ad-sbras.nsc.ru

The classical oil-water displacement problem with given full discharge and without capillary and gravitational forces in one-dimensional case is considered (The Buckley — Leverett problem). Six explicit difference schemes are compared using this initial-boundary value problem and methods of a posteriori analysis. Two new modifications of the well-known schemes are suggested.

В настоящее время широкое распространение получили разностные схемы повышенной точности для сквозного расчета разрывных решений гиперболических уравнений. Разрабатываются новые принципы построения разностных схем и методов их анализа [1, 2]. Развиваются эти подходы, как правило, для уравнений газовой динамики, "мелкой воды" и их модельных одномерных аналогов вида

ди + df (и) = 0 dt дх

где f (и) = 0.5u2 — вогнутая. Развитие численных методов подземной гидродинамики традиционно шло вслед за вышеуказанными разделами вычислительной механики c учетом той особенности, что в задачах двухфазной фильтрации f (u) не является ни выпуклой, ни вогнутой. В настоящей работе сравниваются результаты применения некоторых классических разностных схем и их модификаций к решению задач фильтрации.

Одномерная модель фильтрации двухфазной жидкости в однородной изотропной пористой среде в отсутствие капиллярных и массовых сил при заданном суммарном расходе фаз Q(t) (модель Баклея — Леверетта (БЛ)) имеет вид [3]

mst + Q(t)bx = 0, (1)

где х £ [0,L] — пространственная координата; L — расстояние от нагнетательной скважины до эксплуатационной; t — время; s = (si — S0)/(1 — S0 — S0) — динамическая насыщенность порового пространства смачивающей фазой, s1 — истинная насыщенность смачивающей фазы, (S0, S0) = const — остаточные водо- и нефтенасыщенности; m = m0(1 — S0 — S0),

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

т0 — пористость коллектора; Ъ(в) = ^1/(^1 + — коэффициент подвижности вытесняющей фазы (функция Баклея — Леверетта), кг = кг(в) — относительные фазовые проницаемости, ^ = ^г — вязкости фаз (индекс г = 1 соответствует воде, г = 2 — нефти).

Функциональные параметры модели описаны в [3] и определяются по экспериментальным данным. Наиболее типичными свойствами функции Ъ(в) являются следующие: 6(0) = 0, 6(1) = 1, Ъ'(в) > 0, Ъ'(0) = Ъ'(1) = 0. Существует в* € (0,1) такая, что Ъ'' (в)(в* - в) > 0.

_ г

Введем безразмерные переменные: ж = ж/Х, I = (тХ)-1 / (далее черта над

о

безразмерными переменными опущена), при этом уравнение (1) запишется в виде

в4 + Ьх = 0. (2)

Для уравнения (2) будем рассматривать следующую начально-краевую задачу:

в|х=о = 1, в|г=о = во(ж). (3)

Введем сетку с распределенными узлами UhT = {xi = ih, tn = пт, i = 0, N, n = 0,1, 2...}, h = 1/N — шаг по пространственной координате; т = rh — шаг по временной переменной, r = т/h. В численных расчетах использовались шаги сетки h = 0.01 и т = 0.001. В дальнейшем в разностных схемах используются обозначения, принятые в [4].

На каждом временном шаге вычисляются такие характеристики решения, как положение xc(t) — фронтовой водонасыщенности в БЛ-модели sc, которая определяется решением нелинейного уравнения

( \ -h i \ b(sc)

dS(sc) = b*(sc) = ~

1

с помощью метода деления пополам, n(t) — обводненность пласта n(t) = 100% -J s(x, t)dx.

0

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

k1 = s2, k2 = (1 - s)2, s0(x) = 0, j = 0.1.

При этом параметров значение фронтовой водонасыщенности sc равно 0.30145.

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

hj = h/2j-1, Tj = rhj, j = 1, k, k = const.

Предполагая, что разностное решение Sj = s^ при h ^ 0 в некоторой точке xi сходится с порядком p к искомому точному решению s, получим, что Ssj = Sj — s имеет порядок O(hp), т.е.

Ssj = sj(x,T) — s(x,T) к Chp(x'T),

где C — ограниченный функционал, не зависящий от h; T = пт — некоторый момент времени. Вычитая из Ssj представление для Ssj+1, получим S*sj = sj(x,T) — sj+1(x,T) к

С(Щ — Предполагая, что |С| > е > 0, выпишем формулу для экспериментального

определения порядков сильной локальной сходимости разностного решения [2]:

1

Рз ™10§2 ^з! 3 = 1,к — 2- (4)

При проведении расчетов на нескольких сетках (к > 3) можно найти осредненный порядок

к-2

Р = Е Рз

к — 2 "=1

Аналогично вводится формула для нахождения порядков слабой локальной сходимости

10§2 ^-^3 = 1,к — 2 (5)

где а € [0,1], ^"(^Т) = / 8"(£,Т— интеграл от кусочно-линейной по х функции 8":

8" (х, Т) = 8™ + (х — Х"г)(з™+1 — 8™г)/й, х € [Х"г, Х"г+1], I = 0, N — 1. Осредненный порядок слабой сходимости определяется по формуле [1]

к-2

* = £

гз

к—2 з=1

На рисунках данные по слабой сходимости показаны темными кружками, для сильной локальной сходимости — светлыми квадратами. В нижней части рисунка приводится решение, для которого рассчитывались порядки при Т = 0.3, к = 5.

Отметим, что для задачи (2), (3) при 8о(х) = 0 легко построить точное разрывное решение, так как на момент времени £ координата точки разрыва равна хс = Ъ3(8С)£, а для координаты х3 насыщенности 8 € [8С, 1] имеем формулу х8(£) = Ъ8(8)£. Это решение на рисунках, где приведено семейство графиков решений, отмечено тонкой линией.

Двухслойная схема с центральной разностью по пространству (ВВЦП) и ее модификация

Запишем разностную схему с аппроксимацией конвективного члена центральной разностью в следующем виде:

г , I т , г, т

^ + Ъ0,, = / ,г = 1, N — 1, (6)

где /П имеет вид (С™+1/2ЪП,г — 1/2Ъ™,г)/й, еП+1/2 = с((8? + функЦия е(8) выбрана

следующей

' (8 — (^2 — 8)а2 ,8 <

Ф)" 1 0,8 >&.

Первое дифференциальное приближение для схемы (6) имеет вид

84 + ъж = й ((ее — Ъл) .

Разностная схема (6) аппроксимирует исходное уравнение с погрешностью 0(еИ + И2 + т).

Проведены расчеты с разными параметрами, хорошие результаты получены при 51 =0, S2 = 0.4, «1 =0.1, а2 = 1. При этом выбирались так, чтобы £ [51, S2].

На рис. 1 жирными линиями показано решение с использованием разностной схемы (6) при е = 0, пунктиром — решение при е(И) = 10И (в этом случае схема имеет второй порядок аппроксимации по пространству) и тонкими линиями — точное решение.

Как видно, разностная схема (6) в отсутствие искусственной вязкости дает нефизичное численное решение с выполаживанием (образованием полки) вблизи й = й* = 0.5 > , также наблюдается отставание фронта Xf (£). Квадратичная вязкость с е(И) улучшает вид решения, но незначительно. Поэтому е выбиралось в интервале [2, 3], такая регуляризация схемы (6) позволяет избавиться от выполаживания и получить решение, близкое к точному (рис. 2 — пунктир). В численных экспериментах с разными значениями 52 установлено, что с его увеличением точность определения решения падает. Так, на рис. 2 жирными линиями показаны решения при 52 = 0.8, ориентировочно лучше брать 52 не больше, чем + 0.1.

Отметим, что явная двухслойная схема с центральной разностью (е = 0) для уравнений газовой динамики является абсолютно неустойчивой и не используется в численных расчетах. Относительная устойчивость этой схемы в данной задаче объясняется специфическим видом функции Ь(й) (наличием зоны, где Ь < 0, при й > й*), а также начальным приближением й0(х) = 0. На гладких начальных данных типа й0(х) = (1 — х)5 неустойчивость схемы ярко проявляется и в задачах фильтрации мощными осцилляциями в окрестности фронта й = .

На рис. 3 приведены порядки локальных сходимостей для схемы (6) при е = 2. Как видно, в окрестности фронта порядки осциллируют, что является известным фактом [2]. Немонотонный характер порядков сходимостей наблюдается и при й £ [52,0.6), т.е. там, где сопрягаются гиперболическое (st + Ьх = 0, с(й) = 0, й £ [52,1)) и параболическое (st + Ьх = еИ(с(й)Ьх)х), й £ [0, 52] уравнения. Особенно это заметно на рис. 4, где приведены

Рис. 3.

Рис. 4.

графики локальных сходимостей при 80(х) = (1 — х)5, Т = 0.2. Видно, что в среднем схема имеет первый порядок аппроксимации. Интересно отметить, что для гладких начальных данных осцилляции в порядках слабой сходимости на рис. 4 выделяют кроме разрыва еще и точку выключения регуляризатора, график порядка сильной сходимости выделяет также точку возникновения градиентной катастрофы х « 0.26. Последнее для компактных схем отмечено в [1].

Схема Лакса и ее модификация

Схема Лакса получается заменой 8™ в схеме ВВЦП на среднеарифметическое значение:

sn+1

'L + bn =0,i = 1,N - 1, (7)

T x,i v 7

где zn = (sn+1 + sn-i)/2.

Для схемы (7) первое дифференциальное приближение имеет вид

h2 T

St + bx = Sx - -bsbx)x.

2T 2

Эта схема аппроксимирует исходное уравнение с погрешностью O(h2 + h2/T + т) и условно устойчива при T < h/ max bs.

Аппроксимация уравнения БЛ-модели разностной схемой (7) имеет условный характер. Так, при h/T = const аппроксимируется исходное уравнение, а при h2/T = const производится аппроксимация некоторого параболического уравнения st + bx = £Sxx.

Решение, полученное по схеме Лакса (жирная линия при r = 0.1 на рис. 5), сильно размазывает скачок водонасыщенности и завышает решение, что приводит к большому дисбалансу по водосодержанию (площадь подграфика). Пунктиром отмечено решение при r = 0.2. Видно, что при h/T = 1/r ^ 0 дисбаланс уменьшается и решение приближается к точному. То есть с уменьшением временного шага T увеличиваются дисбаланс и расхождение с точным решением, что характерно для схем с условной аппроксимацией и является главным их недостатком: схема хороша при больших r, а устойчива при r < b-1.

Отметим одну из возможных модификаций схемы Лакса. Для того чтобы уменьшить размазывание скачка и подъем графика вблизи х = 0, выберем специальное усреднение вида

= (в15-1 + в25? + вз5+1 )/в4, в1 + в2 + вз = в4. (8)

Первое дифференциальное приближение для схемы (7) с усреднением (8) имеет вид

+ Ь Т (ЬЬ ) ■ (вз - в1) Н + Н2 (в1 + вз) + 0(Н2 + 2)

5 + Ьх = --(6А)х + ----5Ж + -----5ЖЖ + 0(Н + Т ).

2 в4 т 2Т в4

Мы можем сохранить второй порядок по Н и убрать один из членов, дающих условную аппроксимацию, положив в1 = в3, тогда получим полное дифференциальное приближение вида

н 2 в1

* + Ч- 2Ь*)х = *

Для уменьшения эффекта условной аппроксимации можно выбрать /= СтЬ2, что приводит нас к схеме типа Лакса — Вендроффа

* =( - 2) bs= Ch2( t1 - 2<У ЬК)х '

недостатки которой обсуждаются ниже. Как говорилось ранее, схема Лакса хороша при больших r и плоха при малых, а отвечает за это коэффициент С1 = —в1. Если взять

Г в4

в1 = 1, в2 = 1/r, то получим С1 = h/(2r + 1), который имеет порядок не ниже первого равномерно по г. На рис. 6 решение по этой схеме показано жирной линией. Дальнейшие экперименты показали, что близость численного решения к точному в зоне волны разряжения можно улучшить, если взять в1 = d(s), где

*)={!!? >- ¿Г1""£ *•

выбиралось так, чтобы > 8С.

Тогда вк имеют вид

в1 = вз = ф),в2 = 1/г. (9)

При нашем выборе параметров хорошие результаты получены при = 0.4, а = 0.1.

На рис. 6 пунктиром показаны решения, полученные по схеме Лакса со специальным усреднением (8), (9). Видно, что усреднение (9) позволило уменьшить размазывание разрыва, улучшить баланс и достичь хорошей близости к точному решению. Соответствующие графики локальных сходимостей представлены на рис. 7. Схема имеет ярко выраженный первый порядок сходимости. Но следует отметить, что локальные порядки сильной и слабой сходимости ведут себя более плавно, с меньшими осцилляциями, нежели для схемы (6).

ТУО-модификация схемы Лакса — Вендроффа

Если в схему ВВЦП ввести искусственную вязкость, гасящую ее схемную вязкость первого порядка, то получится одношаговый вариант схемы Лакса — Вендроффа [5]:

Т

т

_и Ьп — —(Ьп Ьп —Ьп Ьп \ + жг = 2й (Ъ®г+1/2Ъх,г Ъвг—1/2Ъх,г)'

(10)

Схема (10) имеет погрешность аппроксимации 0(й2 + т2). К недостаткам схемы Лакса — Вендроффа следует отнести сильное выполаживание решения вблизи 8* > 8С, хотя баланс обводненности соблюдается неплохо, при этом имеет место отставание фронта Хf (£). Для того чтобы улучшить схему (10), используют ТУБ-модификацию схемы Лакса — Вендроффа [6]. Представим схему в виде

8п+1

+5+1/2

е

г—1/2

т

= 0,

(11)

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

5п+1/2 = 2 (Ъп + Ъп+1 — йфг+1/2Ъп+1/2Ъп,г + (1 — фг+1/2)(Ъп+1 — Ъп)

8п

ьп I ьп

гп _ ^вг 1 ^¿+1

°вг+1/2 — ^ ' ( )

В литературе предлагаются различные функции-ограничители фг+1/2 [6, 7]:

Фг+1/2 — тах[0, тт(1, иг+1/2)], (13)

Фг+1/2 — тах[0, т1п(2иг+1/2, 1), тт(2, иг+1/2)], (14)

, К+1/21 + иг+1/2 ( Л

фг+1/2 — -¡—:-, (15)

|иг+1/2| + 1

иг+1/2 — (5П - 5П-1 )/(^П+1 - 3П).

Проведены расчеты с разными функциями фг+1/2 и установлено, что наилучшая точность определения решения в окрестности разрыва была при использовании формулы (15). Решения, полученные с применением (13), оказались ближе к точному решению в гладкой части, но на разрыве хуже, и, наконец, формула (14) была наихудшей.

На рис. 8 жирными линиями показаны решения по схеме Лакса — Вендроффа, пунктиром обозначено решение для ТУБ-модификации (11), (12), (15) и тонкими линиями с кружками — решение для ТУБ-модификации (11), (12), (14). Как видно, использование ТУБ-модификаций позволило избавиться от выполаживания решения. Контроль баланса при этом также несколько улучшился по сравнению со схемой Лакса — Вендроффа (см. таблицу). Соответствующие порядки локальных сходимостей для схемы (11), (12), (15) приведены на рис. 9. Поскольку начальный профиль имеет разрыв при х — 0, схема имеет первый порядок аппроксимации. При начальном приближении, более гладком, 50(х) — (1-х)5, Т — 0.2, локальные порядки оказались приблизительно равными 2 (рис. 10).

В задачах фильтрации важен контроль за массовым балансом и фронтом вытеснения. Данные, полученные по этим показателям, для рассмотренных схем приведены в таблице. Для контроля за фронтом приведено отклонение расчетного положения Xf от теоретического хс в единицах шагов сетки К.

Дисбалансы водосодержания на два момента времени, % и тах К(¿)

Название схемы t = 0.3 t = 0.6 K

ВВЦП 0.37 10.15 13

Модификация схемы ВВЦП 0.01 0.35 4

Схема Лакса 23.41 12.67 29

Схема Лакса с усреднением (8), (9) 4.09 0.65 6

Схема Лакса — Вендроффа 1.62 3.47 5

Схема (11), (12), (15) 1.03 0.59 3

Таким образом, стандартные схемы для расчета гиперболических уравнений в задачах фильтрации оказываются малопригодными. Это связано с появлением таких эффектов, как осцилляции, выход на полку, поднятие решения, опережающее движение фронта ж/(¿). Все это приводит к необходимости модифицировать стандартные схемы для успешного их применения. Метод создания ТУО-схем оказался эффективным для построения монотонных разностных схем, но в то же время и относительно сложным. Данное исследование показывает, что из рассмотренных схем конкурировать между собой могут лишь ВВЦП с регуляризатором и ТУБ-модификация схемы Лакса — Вендроффа с корректором потоков (15) [7]. Первая схема проста в реализации и дает меньший дисбаланс. Вторая меньше размазывает скачок, но дает больший дисбаланс и сложнее в реализации.

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

[1] ВоЕводин А.Ф., Остапенко В.В. О расчете прерывных волн в открытых руслах // Сиб. журн. вычисл. математики. 2000. Т. 3, № 4, С. 305-321.

[2] Остапенко В.В. Разностная схема повышенного порядка сходимости на нестационарной ударной волне // Сиб. журн. вычисл. математики. 1999. Т. 2, № 1. С. 49-54.

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

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

[5] Lax P.D., Wendroff B. Systems of conservation laws // Comm. Pure Appl. Math. 1960. Vol. 13. P. 217-237.

[6] Harten A. High resolution schemes for heperbolic conservation laws //J. Comp. Phys. 1983. Vol. 49, № 3. P. 357-372.

[7] Van Leer B. Toward the ultimate conservative difference sheme. V. A second-order sequel to Godunov's method //J. Comp. Phys. 1979. Vol. 32, № 1. P. 101-136.

Поступила в редакцию 10 февраля 2003 г.