Вычислительные технологии
Том 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 г.