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

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

CC BY
181
37
i Надоели баннеры? Вы всегда можете отключить рекламу.
Область наук

Аннотация научной статьи по физике, автор научной работы — Бочаров О. Б., Телегин И. Г.

Рассматривается классическая задача вытеснения нефти водой без учета капиллярных и гравитационных сил в одномерном случае при заданном суммарном расходе фаз(задача Баклея-Леверетта). На примере этой начально-краевой задачи, используя методы апостериорного анализа, сравниваются шесть явных разностных схем. Приводятся также две новых модификации известных ранее схем.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Comparative analysis of some difference schemes for two-phase filtration problems without capillary forces

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.

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

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

Том 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), где

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

*)={!!? >- ¿Г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 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.