CC BY
45
Вычислительные технологии
Том 13, № 1, 2008
Численное исследование разрывных решений задачи сопряжения двух неизотермических моделей фильтрации*
О.Б. Бочаров
Институт водных и экологических проблем СО РАН, Новосибирск, Россия
e-mail: obb-52@inbox.ru
И. Г. ТЕЛЕГИН
Научно-исследовательский и проектный институт “КогалымНИПИнефть”,
Тюмень, Россия e-mail: TeleginIG@tmn.lukoil.com
The problem of conjugation of two various mathematical models of a nonisothermal filtration in the homogeneous layer, without taking into account gravitational forces is considered. The one-dimensional case is studied when the aggregate rate of phase flux is given. With the help of numerical methods the discontinuous solutions are investigated which arise at conjugation of Muskat—Leverett and Bakley—Leverett temperature models.
Введение
В настоящее время задачи сопряжения различных математических моделей активно исследуются в связи с наличием в средах различных пограничных слоев, которые малы по размерам, но оказывают значительное влияние на общую структуру решения. Такими областями применительно к задачам нефтедобычи являются прискважинные зоны, т. е. участки вблизи нагнетательной и эксплутационной скважин. Здесь более физично использование моделей, в которых не учитываются капиллярные силы (например, модель Баклея—Леверетта (БЛ-модель)), в основной части пласта роль капиллярных сил оказывается существенной, что приводит к необходимости использования более сложной модели Маскета—Леверетта (МЛ-модель), учитывающей капиллярные эффекты.
Задача сопряжения моделей Маскета—Леверетта и Баклея—Леверетта достаточно подробно численно исследована в работах [1-3]. В данной работе с помощью численных методов исследуются разрывные решения, которые получаются при сопряжении температурной модели Маскета—Леверетта (МЛТ) и температурной модели Баклея— Леверетта (БЛТ).
*Работа выполнена при частичной финансовой поддержке СО РАН (интеграционный проект № 117).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
1. Описание моделей
В однородной изотропной пористой среде в отсутствие массовых сил одномерная модель неизотермической двухфазной фильтрации Маскета—Леверетта (МЛТ) имеет такой вид [4]:
mst = (koao(pCsSx + Рев0*) - Q(t)6)* = —v*, m
0t = (A0* - Q(t)0)*, (1)
где x G [0,L] — пространственная координата, L — расстояние от нагнетательной скважины до эксплуатационной; t — время; s = (si — S0)/(1 — Sj1 — SJ°) — динамическая насыщенность смачивающей фазы (воды), s1 — истинная насыщенность воды, (S0,S20) = const — остаточные водо- и нефтенасыщенности; m = m0(1 — Sj1 — S2°), m0 — пористость коллектора; 0 G [0min,0max] — температура; a0(s,0) = — h2b/^2; Pe(s,0) = (m0/k0)1/2Y (0)j (s) — капиллярное давление, y(0) — коэффициент поверхностного натяжения, j(s) —функция Леверетта; h0 —проницаемость; b(s,0) = k1 (s)/(k1(s) +
^(0)k2(s)) — коэффициент подвижности вытесняющей фазы, hi(s) — относительные
3
фазовые проницаемости, ^ = ^1/^2, ^i(0) — вязкости фаз; A(s,0) = ^2 aiAi/(picpi) —
i= 1
коэффициент температуропроводности смеси, а1 = m0s1, а2 = m0(1 — s1), а3 = 1 — m0 (индекс i =1 соответствует воде, i = 2 — нефти, i = 3 — коллектору), cpi = const — теплоемкость фазы при постоянном давлении, Ai = Ai(0) — коэффициенты теплопроводности, pi — плотности фаз; Q(t) — общий расход смеси; v — расход вытесняющей фазы.
Свойства функциональных параметров модели описаны в [4]. Отметим, что h1(0) = h2(1) = 0. Это приводит к вырождению типа уравнения для водонасыщенности (a0(O,0) = a0(1,0) = 0), что ведет к плохой обусловленности естественных граничных условий (градиенты решения в окрестности границы становятся бесконечно большими). Применительно к изотермическому случаю сложности, вызванные этим обстоятельством, описаны в [5].
Уравнение для s(x,t) системы (1), как следует из определения функциональных параметров, можно переписать в эквивалентном представлении:
mst = (ha Ре* — b)* = —v*. (2)
Положив Q(t) = Q0, введем безразмерные переменные: X = x/L, t = Q0t/(mL), 0 = (0 — 0min)/(0max — 0min), A = A/A0 (далее черта над безразмерными переменными опускается). В силу доказанного в [4] принципа максимума 0max и 0min достигаются на границах области при x = 0 и x = 1. С учетом представления (2) система уравнений
(1) в новых обозначениях запишется в следующем виде:
st = (£ap* — b)* = — V*^ (3)
0t = (£в A0* — m0)® = — 5*
где £ = 70(m0h0)1/2/(Q0L^0), a(s,0) = — h^/^(h1 + М2)), p(s,0)= j'7*, £e = mA0/(Q0L),
7* = 7/70, ^ = ^*/^2, ^* = ^1/^0, ^2 = ^2/^0, 70 = max (y(0)), ^0 = max (^(0)),
ве[0,1] ве[0,1]
A0 = A(0, 0). Звездочки у безразмерных величин в дальнейшем опускаются. При значении параметра £ = 0 будем иметь тепловую модель Баклея—Леверетта (БЛТ).
Обозначим через Гг = {x,t|x = l,t > 0} линию на отрезке x G [0,1], при x = 0 это будет координата для нагнетательной скважины, при x = 1 — координата добывающей скважины. Через Qzn обозначим область {x,t|z < x < n,t > 0}, z G [0,l), n G (l, 1].
2. Постановка задачи сопряжения
Нагнетание вытеснителя и добыча нефти производятся через скважины, которые моделируются точечными источниками и стоками. В результате этого в окрестности скважины могут развиваться достаточно большие скорости, которые сводят работу капиллярных сил к нулю. Рассмотрим случай, когда вблизи эксплутационной скважины капиллярные силы малы. Разобьем область П01 на две, в области П0г фильтрация двухфазной жидкости описывается МЛ-моделью и соответственно этому в (3) £=0, а в области 0,ц — БЛ-моделью и в (3) £ = 0. Предлагается для решения системы (3) в области П01 рассматривать следующую начально-краевую задачу:
s|£=0 1 0|*=0 01; [s]£=1 [v]£=1 [0]*=г [5]Ж=г 0; £в A0*|a;=1 0; (4)
s|t=0 = s0(x), 0|t=0 = 00 (x), xG(0,1],
где [/]*=г = {/(l — 0,t) — /(l + 0,t)} — скачок функции /(x,t) на линии Гг. Из условий склейки на Гг следует условие £арж|ж=г-0 = 0 и задача (3), (4) распадается на три
подзадачи. В области П0г уравнение для насыщенности имеет вид
st = (£ap* — b)* = —v*; s|*=0 = 1, £ар*|*=г = 0; s|t=0 = s0(x), x G (0,1]. (5)
В области Пг1 процесс фильтрации описывается БЛТ-моделью:
st = —b*; s|*=г = Si(t), s|t=0 = s0(x), x G (0,l], (6)
где Si^) определяется из решения (5).
Во всей области П01 решается уравнение для температуры:
0t = (£вA0* — m0)*; 0|ж=0 = 1 £вA0^|^=1 = 0; 0|t=0 = 0o(x), x G (0,1]- (7)
В отличие от изотермического случая в неизотермическом имеет место обратное влияние процессов в области Пг1 на область П0г посредством температуры, однако в силу малой величины коэффициента при 0** это влияние пренебрежимо мало.
В работе [1] с целью обоснования непрерывности расхода было предложено сменить эволюционную переменную в области Пг1 и искомую насыщенность s на расход v. Такой подход дает удовлетворительные результаты при ^ = const, в случае же, когда ^ = ^i(0), такая постановка становится слишком сложной. Это связано с тем, что в уравнении для насыщенности (7) появится член с 0t, что сильно осложняет нахождение решения.
3. Особенности численной реализации
Введем сетку с распределенными узлами = {xi = ih, tn = пт; i = 0, N, п = 0,1, 2,...}, h = 1/N — шаг по пространственной координате, т = Kh2 = rh — шаг по временной переменной. Пусть при этом линии Гг соответствует номер N = IN узла на сетке (предполагается, что N не дробное число). В численных расчетах шаг h брался равным 0.005 (N = 200), а т = 0.00025 (K =10, r = 0.05). В данной работе используются обозначения, принятые в [6].
Аппроксимируем уравнение для температуры (7) явной схемой второго порядка точности на гладких решениях:
0П+1т— 0П = (An+, 0n,i — An-10n,i) _ mAn+1/2 — A-1/2, i =1^—1; (8)
/in+1 nn n
0N 0N 2£e \n дп ™fln .
T = —~hAN -10S’N — m0*,N;
0,n+1 = 0n = 01; 00 = 00(xi), i = !7N,
где A'+1 = A((sn + Sn+1)/2), (0n + 0i+1)/2), разностный поток An+1/2 определяется следующим образом:
An+1/2 = °.5(0n + 0n+1 — гф”+1/2(0П+1 — 0n) — (1 — <Ai+1/2 )(0i+1 — 0n)), (9)
функция-ограничитель <Ai+1/2 была подобрана в работе [7]:
nn 0i 0i-1
(?Ai+1/2 = max[0, min(10ui+1/2,1), min(2, г/^^)], Mi+1/2 = jn-----------i-n. (10)
0i+1 0i
Отметим, что предпочтение явным схемам для уравнения температуры было вызвано, во-первых, простотой реализации схем высокого порядка точности, во-вторых, тем, что коэффициент вв мал и не играет роли серьезного ограничителя шага по времени. Уравнение (5) аппроксимировалось явно-неявной схемой второго порядка:
„га+1 _ „га е
-—1 = г (<п «1 р;+ 1 _ <_ 1 р;+1) + (1 _ ^1)(«?+1 р;,< _ <-1 р;,<))_
Т г ‘+ 2 ’ ‘ 2 ’ ‘+ 2 ’ ‘ 2 ’
—а^1 — (1 — о"1 )b*i, i = 1,Ni — 1; (11)
sn+1 — sn 2£
~ 1 = — (_Т (а1 °nNl-1 Pn+V1 + (1 — a1)aN;- 1Pn,N ) — ^ — (1 — a1)bn+V1;
sn+ = sn =1; s; = s0(xi), i = 1, N. ’i+1//2. (0,n+1 + 0n+11)/2), 6П+1 = bn + b"(sn
Здесь a”, I = a((s' + si+1)/2, (0n+1 + 0™+11)/2), b'+1 = bn + b™(sn+1 — s'), аналогично лине-
аризовывалось и р™+1, а1 — вес, при расчетах предполагалось а1 = 0.5. Для численного решения системы (11) применялся метод правой прогонки.
Уравнение для БЛТ-модели аппроксимировалось ТУБ-модификацией схемы Кран-ка—Николсона:
„га+1 _ „га 1
“ т ~ = _ г (а2(<^;+1/2 _ <^^г-1/2) + (1 _ а2)С+1/2 _ С+1/2^
г = ЖТ1Ж 1 = ^Г1, (12)
где а2 = 0.5 — вес, разностный поток <?’+1/2 определяется следующим образом:
3+1/2 = о.5(ь? + ь;+1 _ (1 _ ф;+1/2)(ь;+1 _ ь;)). (13)
Проводились численные эксперименты с различными функциями-ограничителями Фг+1/2. Была выбрана оптимальная функция (для схемы Кранка—Николсона) [8]:
Фг+1/2 = тах[0, ш1п(1,м^+1/2)]. (14)
Схема (12)—(14) неявная, линеаризация проводилась с использованием метода простой итерации. Уравнения решались итерационно, итерации прекращались при достижении точности 0.01т.
Для БЛТ-модели также использовалась явная ТУБ-модификация схемы Лакса— Вендроффа:
сп+1 сп пп — пп
3-----^ = —Уг+1/2^-1/2 , . = + 1, 5^+1 = 57+1, (15)
т К 1
в этом случае разностный поток <тП+1/2 определяется следующим образом:
С+1/2 = 0-5[6П + 6П+1 _ г^г+1/2(сП+1/2)2 _ |СП+1/2|(5П+1 _ ^ _ ^г+1/2)], (16)
"01+1/2 = <^г+1/2(5п+1 _
п = , (ЬП+1 _ ЬП)/(5П+1 _ О, (5П+1 _ 5П) = 0;
'т/2 I ^, (5П+1 _ О = 0.
В работе [9] были рассмотрены различные функции-ограничители ф для уравнения Баклея—Леверетта и была выбрана оптимальная функция для данной схемы:
1 _ |и*+1/21 + иг+1/2 _ ^ п „п \ к„п „п\ (17)
01+1/2 = ---:------Г”—;---, «г+1/2 = <А _ 5г-1)/(5г+1 _ )- (17)
|иг+1/2| + 1
На каждом временном шаге вычислялась п(£) — текущая обводненность пласта:
1
„(«) = / .(х,()<Ь ■ 100 %.
0
Интеграл в правой части вычислялся по формуле трапеций.
В численных расчетах использовался следующий набор модельных параметров:
к1 = 52, к2 = (1 _ з)2, ] = (1 _ з)/(0.9 + з); = 5°° = 0, з0(х) = 0;
^2(0) — ^2тах + (^2 тт ^2шах)^> Т(0) — Ттах + (Ттт Ттах)0, Ттах — 1,
£ = 0.5, ^1 = 0.1, т = 0.36, I = 0.8.
Выбор ^1 объясняется тем, что вязкость воды слабо зависит от температуры.
Остальные данные брались из [10]: А1 = 0.644 Вт/(м-К), А2 = 0.08 Вт/(м-К), Аз = 2.40 Вт/(м-К), р1 = 1000 кг/м3, р2 = 730 кг/м3, р1 = 1000 кг/м3, р2 = 730 кг/м3, р3 = 4216 кг/м3, ср1 = 4071 Дж/(кг -К), ср2 = 2100 Дж/(кг -К), ср3 = 920 Дж/(кг-К).
На рис. 1-8 толстыми линиями обозначены решения или характеристики, относящиеся к задаче сопряжения, тонкими — результаты контрольного счета по МЛТ-модели во всей области, тонкими линиями с кружками — профили температуры, пунктиром — линия сопряжения моделей.
4. Особенности разрывных решений задачи сопряжения (4)—(6)
Закачка горячей воды моделировалась с помощью задания = 1, 0о = 0. На рис. 1 даны графики решений по модели сопряжения и по МЛТ-модели, а также графики для температуры, при 7тіп = 0.5, ^2 = 1, є# = 10-4 (вариант 1). При прохождении линии сопряжения фронтом водонасыщенности в прискважинной зоне формируется баклеевский прямой скачок (рис. 1, і = 0.3). Отметим, что во всех просчитанных вариантах при прохождении фронтом насыщенности линии сопряжения формируется этот скачок, также во всех вариантах отмечаются подъем и выполаживание графика водонасыщенности в точке х = I — 0, обусловливаемые наличием в этой точке краевого условия єарх|х=г = 0, это показано в работах [2, 3]. Самые интересные изменения в решениях происходят при прохождении линии Гг температурным фронтом. На это и будем обращать внимание. Совпадение графиков для температуры объясняется малой величиной є#. В [10] отмечалось, что в случае, если от температуры зависит 7, перед температурным фронтом образуется локальный максимум (за счет источникового члена, пропорционального > 0, это видно из формулы (1)), а после температурного фронта — локальный минимум (за счет соответственно стокового члена, пропорционального
< 0). Наличие зоны действия БЛТ-модели в области Пгі изменяет структуру решения з(х,і), приближая друг к другу локальный минимум и локальный максимум, т. е. формируется обратный скачок насыщенности. Он возникает из-за немонотонности профиля 5г(і) и не связан непосредственно с температурным фронтом, поэтому при і = 2.4 отстает от температурного фронта и соответственно от структуры минимум-максимум в МЛТ-модели.
При БЛТ-модели в случае, если от температуры зависят только вязкости фаз, перед температурным фронтом на графике водонасыщенности образуется полочка (выполаживание), а после фронта — дополнительный вал вытеснения. МЛТ-модель сильно размазывает эти эффекты, и при больших є они слабо заметны. На рис. 2 приведены решения при следующих параметрах: 7 = 1, ^2тах =1, ^2тт = 0.2, є# = 10-6 (вариант 2).
0.0 0.2 0.4 0.6 0.8 X 0.0 0.2 0.4 0.6 0.8 X
Рис. 1. Закачка горячей воды: 7 = 7(0), ^2 = 1 Рис. 2. Закачка горячей воды: 7 = 1, ^2 =
М0)
Из рисунка видно, что отсутствие капиллярных сил, размазывающих области выпола-живания и дополнительного фронта в районе действия БЛТ-модели, приводит к появлению “стакана” на профиле водонасыщенности. Стенки “стакана” образованы прямым и обратным скачками з(х, і). Прямой скачок получается из дополнительного фронта вытеснения, его можно наблюдать только при малых є# и по схемам повышенного порядка
Рис. 3. Закачка горячей воды: 7 = 1, ^2 = ^2(0), решения по схемам первого порядка
Рис. 4. Графики в(1,і) к рис. 2
Рис. 5. Закачка горячей воды: ^2 = ^(0), 7 = 7(0)
Рис. 6. Графики в(1,і) к рис. 5
точности для БЛТ-модели и уравнения температуры. Например, на рис. 3 представлены результаты, полученные по неявным противопотоковым схемам первого порядка аппроксимации (при тех же параметрах). Видно, что дополнительный фронт вытеснения не носит разрывного характера, а “стакан” превращается в выемку. Обратный скачок на рис. 2 получается в результате интенсивного понижения водонасыщенности перед температурным фронтом. В БЛТ-модели нет капиллярных сил, и поэтому формируется решение с двумя разрывами. На рис. 4 нарисованы графики з(1,£). На этом рисунке хорошо виден “стакан” с двумя стенками.
Случай, когда от температуры зависят и капиллярное давление, и вязкости, показан на рис. 5 при следующих параметрах: 7тіп = 0.5, ^2тах = 1, ^2тіп = 0.2, є$ = 10-6 (вариант 3). При прохождении температурным фронтом линии сопряжения опять образуется “стакан” в профиле водонасыщенности. Но он имеет свои особенности, объединяя варианты 1 и 2. Прямой скачок формируется, как и в варианте 2, за счет влияния температуры на вязкости фаз и очень похож на соответствующий вариант. Обратный скачок формируется за счет понижения водонасыщенности перед температурным профилем и структуры минимум-максимум в профиле МЛТ-модели. Отметим, что использование явных схем повышенного порядка точности для уравнения Баклея—Леверетта (которые хорошо зарекомендовали себя в изотермическом случае) в вариантах 2 и 3 приводило к осцилляциям на разрывах и в районах понижения водонасыщенности, а при некоторых параметрах — и к разрушению решения. На рис. 5 пунктиром на графике водонасыщенности для задачи сопряжения указаны решения, полученые в области Пц по явной схеме (15)—(17). Видно, что решение осциллирует в районе понижения водона-сыщенности. Поэтому в данной работе было отдано предпочтение неявным схемам для уравнения Баклея—Леверетта.
На рис. 6 даны графики з(1,£), хорошо видно, что “стакан” шире, чем в варианте 2, а обратный разрыв в полтора раза больше. Интересно, что дно стакана не пологое, как
Рис. 7. Закачка горячей воды: ^2 = ^2(0), 7 = 7(0), £в = 0.001
Рис. 8. Закачка холодной воды: ^2 = ^2(0), 7 = 7(0), £в = 0.001
Дисбаланс и разница между решениями при t = 4
Вариант 0 1 2 3
Разница в норме С[0,1] 0.022 0.021 0.033 0.027
Максимальная разница в норме С[0,1] 0.191 0.192 0.220 0.388
Дисбаланс, % 1.00 0.90 1.97 1.23
Максимальный дисбаланс, % 2.23 2.56 6.05 9.36
в варианте 2, а со значительным скосом, это связано с тем, что минимум (сформированный в области П0г) сблизился с максимумом.
В варианте 3 при прохождении температурным фронтом линии сопряжения образуется “клювик” на профиле з(х,£) в точке х = I. Одна сторона “клювика” формируется условием Зх|Х=1 > 0, которое объясняется тем, что в краевом условии £арх1х=1 = 0 при прохождении температурным фронтом линии сопряжения производная 0Х|х= < 0. С другой стороны “клювика”, в области П^, перед температурным фронтом начинается интенсивное понижение водонасыщенности. В результате получаем структуру с острым “клювиком” в точке х = I. При £0 = 10-6 “клювик” имеет вид узкого всплеска, который быстро исчезает. При увеличении £0 (или при использовании схем первого порядка аппроксимации) “клювик” сохраняется дольше, но его амплитуда уменьшается. На рис. 7 приведен пример расчета варианта 3 с £0 = 0.001, на нем “клювик” показан при £ = 2.3.
На рис. 8 показан пример закачки холодной воды при в\ = 0, 00 = 1, 7т;п = 0.5, ^2тах = 4, ^2тт = 1, £0 = 0.001. В этом случае “клювик” образуется в точке х = I, но потом смещается вслед за температурным фронтом. В результате получается волновая структура в области Пц.
В таблице приведены данные по дисбалансу обводненности и разнице между решениями для рассмотренных вариантов (0 — изотермический вариант).
5. Апостериорный анализ сходимости
В общем случае точное решение задачи сопряжения построить очень сложно. Поэтому для определения порядков сходимости часто применяется правило Рунге, заключающееся в том, чтобы провести несколько расчетов задачи с разными шагами пространственной сетки:
hj = h/2j-1, Tj = rhj, j = 1, k, k = const.
Предполагая, что разностное решение Sj=shj при h^0 в некоторой точке Xi сходится с порядком p к искомому точному решению S, получим, что 5sj = Sj — s имеет порядок O(hp), т.е.
5sj = Sj(x,T) — s(x,T) « Chp(x’T),
где C — ограниченный функционал, не зависящий от h; T = пт — некоторый момент времени. Вычитая из 5sj представление для 5sj+1, получим 5*sj = sj(x, T) — sj+1(x, T) ^ C(hp — hp+1). Предполагая, что |C| > e > 0, выпишем формулу для экспериментального определения порядков сильной локальной сходимости разностного решения [11]:
|5Ч
|5*sj+i|
Pj~loS2^------Г, j = 1,k — 2. (18)
Рис. 9. Порядки сходимости в изотермическом случае При расчетах на нескольких сетках к > 3 можно найти осредненный порядок Р:
Р
к-2
р3
^к - 2 ■ 3 = 1
На рис. 9, а светлыми квадратиками показаны данные по сходимости для изотермического решения. На рис. 9, б приводится соответствующее изотермическое решение при Т = 0.15, к = 3, к\ = 0.01, з0(х) = х4. Видно, что в области действия МЛТ-модели решение сходится с порядком 2, в области разрыва водонасыщенности в БЛТ-модели порядки сильно осциллируют. Перед разрывом порядок равен 2.
Выводы
Модель сопряжения может сильно изменять структуру решения уравнения водона-сыщенности. При закачке горячей воды влияние зоны действия БЛТ-модели выражается в образовании разрывных решений з(х,£). Результат приобретает разнообразный характер, зависящий от различных параметров. Данная работа показывает, что разрывные решения можно получать только при использовании схем повышенного порядка аппроксимации.
На первый взгляд, полученные решения задачи сопряжения кажутся очень экзотическими, чтобы быть физичными. Однако это не так. Применительно к задачам сопряжения глинизированных нефтяных пластов в экспериментальной работе [12] были получены похожие результаты.
Список литературы
[1] МонАхов В.Н. Сопряжение основных математических моделей фильтрации двухфазной жидкости // Математическое моделирование. 2002. Т. 14, № 40. С. 109-115.
[2] Бочаров О.Б., Телегин И.Г. Численный анализ некоторых методов сопряжения двух моделей фильтрации несмешивающихся жидкостей // Вычислительные технологии. 2002. Т. 7, № 5. С. 11-20.
[3] Бочаров О.Б., Телегин И.Г. Численное исследование процесса вытеснения при сопряжении различных моделей фильтрации двухфазной жидкости // Наука, культура, образование: Труды ПАНИ. Горно-Алтайск: Изд-во ГАГУ, 2002. № 10-11. С. 118-125.
[4] Бочаров О.Б., МонАхов В.Н. Краевые задачи неизотермической двухфазной фильтрации в пористых средах // Динамика сплошной среды: Сб. науч. тр. Новосибирск: ИГиЛ СО АН СССР, 1988. Вып. 86. С. 47-59.
[5] Коновалов А.Н. Задачи фильтрации многофазной несжимаемой жидкости. Новосибирск: СО Наука, 1988. 166 с.
[6] Самарский А.А. Введение в теорию разностных схем. М.: Наука, 1971. 552 С.
[7] Бочаров О.Б., Осокин А.Е., Телегин И.Г. Сравнение разностных схем для одномерной модели переноса примесей в речных руслах // Сопряженные задачи механики, информатики и экологии: Матер. междунар. конф. Томск: Изд-во ТГУ, 2004. С. 39-40.
[8] Harten A. High resolution schemes for hyperbolic conservation laws // J. Comp. Phys. 1983. Vol. 49, N 3. P. 357-372.
[9] Бочаров О.Б., Телегин И.Г. Сравнительный анализ некоторых разностных схем для задач двухфазной фильтрации без учета капиллярных сил // Вычисл. технологии. 2003. Т. 8, № 4. C. 23-31.
[10] Бочаров О.Б., Телегин И.Г. О некоторых особенностях неизотермической фильтрации несмешивающихся жидкостей // Теплофизика и аэромеханика. 2002. Т. 9, № 3. С. 459-466.
[11] Остапенко В.В. Разностная схема повышенного порядка сходимости на нестационарной ударной волне // Сиб. журн. вычисл. математики. 1999. Т. 2, № 1. С. 49-54.
[12] ХАвкин А.Я., Чернышев Г.И. Исследование механизмов повышения нефтеотдачи с помощью метода радиоактивных индикаторов // Тр. VIII Междунар. выставки “Нефть, газ, нефтехимия — 2001”. Казань, 2001. Т. 2. С. 350-358.
Поступила в редакцию 17 августа 2007 г.
