CC BY
41
УДК 532.546 - 622.276.6
ПРИМЕНЕНИЕ МЕТОДА АДАПТИВНЫХ РАЗНОСТНЫХ СЕТОК К ЗАДАЧЕ ДВУХФАЗНОЙ ФИЛЬТРАЦИИ ПРИ РАЗРАБОТКЕ НЕФТЯНЫХ МЕСТОРОЖДЕНИЙ
Дмитрий Сергеевич Евстигнеев
Институт горного дела им. Н. А. Чинакала СО РАН, 630091, Россия, г. Новосибирск, Красный проспект, 54, младший научный сотрудник, научно-инженерный центр горных машин и геотехнологий, тел. (983)127-88-52, e-mail: [email protected]
Реализован метод адаптивных разностных сеток, динамически связанных с решением уравнения Раппопорта-Лиса, учитывающей действия капиллярных и гравитационных сил.
Ключевые слова: двухфазная фильтрация, адаптивные сетки, нефть.
APPLICATION OF ADAPTIVE DIFFERENCE GRID TO TWO-PHASE FILTRATION PROBLEM IN OIL-FIELD DEVELOPMENT
Dmitry S. Evstigneev
Chinakal Institute of Mining SB RAS, 630091, Russia, Novosibirsk, 54 Krasny prospect, Junior Researcher of Mining Machinery and Geotechniques Research Center, tel. (983)127-88-52, e-mail: [email protected]
The researcher realizes the method for adaptive difference grids dynamically related to solution of Rappaport-Fox equation accounting for capillary and gravity force effects.
Key words: two-phase filtration, adaptive grids, oil.
В нефтяном пласте нефть движется по пласту-коллектору к добывающей скважине под действием перепада давления. Движение происходит при условии, что в пласте давление выше, чем на забое скважины. Впоследствии пластовое давление постепенно снижается и для его поддержания требуется проведение специальных мероприятий.
Восполнение энергии пласта и обеспечение продвижения нефти к добывающей скважине происходит методом искусственного поддержания внутрипластовой энергии. Основными источниками пластовой служат энергии:
- напора пластовой воды (краевой или подошвенной);
- расширения растворенного в нефти газа;
- упругости жидкости и породы;
- гравитационного напора нефти.
При водонапорном режиме эксплуатации залежи, в пласт закачивается вода под большим давлением и тем самым обеспечивается необходимый перепад давления для продвижения нефти к добывающей скважине. В зависимости от возникающего перепада давлений и фильтрационных характеристик пласта возникают области, в которых, при совместном движении несмешивающихся флюидов, капиллярными и гравитационными силами можно пренебречь. В этом случае Бакли и Левереттом была предложена модель фильтрации применяемая для прогноза добычи нефти на месторождениях [1]. Согласно предложенной модели, фронт насыщенности терпит разрыв на границе двух фаз: воды и нефти. Последующим усовершенствованием модели Ба-кли-Леверетта стала модель Раппопорта-Лиса. В ней действие капиллярных сил про-
является в основном вблизи фронта вытеснения, где градиенты насыщенности велики и изменяются непрерывно.
Для отыскания разрывных решений существуют два подхода: явный и неявный. Недостаток метода с явным выделением разрыва заключается в необходимости выделения положения разрыва на каждом временном слое [2,3]. Однако при больших деформациях фронта разрыва неравномерные расчетные сетки сильно искажаются, и возникает необходимость их периодической перестройки [2,4]. При неявном методе отслеживания фронта разрыва вводится некая пространственная функция-идентификатор [3,4], и положение разрыва совпадает с линией определенного уровня этой функции или находится в области ее максимального градиента. Существует группа методов, где в расчётную схему включен механизм размытия разрыва [5], а также схемы сквозного счёта типа С. К. Годунова [4,6], в которых потоки через боковые грани расчетной ячейки определяются из решения задачи о распаде произвольного разрыва.
В серии работ, посвященных построению адаптивных разностных сеток, динамически связанных с решением [7-11], предлагается метод, в котором с помощью взаимно-однозначного преобразования координат осуществляется переход от декартовой системы координат к криволинейной, где для численного решения задачи используются расчетные сетки с фиксированными узлами. В физическом пространстве им соответствуют сетки с тем же общим числом узлов, но с управляемым распределением, что позволяет концентрировать их в областях, содержащих большие градиенты и разрывные решения. Переход к нестационарной системе координат избавляет от проблем, связанных с подвижными границами, при этом решения не содержат разрывов, однако в этом случае искомыми являются не только сеточные функции, но и координаты узлов сетки. Распространим данный подход к описанию движения несмешивающихся флюидов в нефтяном пласте, описываемому уравнением Раппопорта-Лиса.
Систему уравнений двухфазной фильтрации составляют: закон Дарси, уравнения баланса масс, уравнения несжимаемости для общего потока и капиллярного давления. В случае одномерного движения она может быть сведена к одному уравнению для водонасыщенности параболического типа, получившему название уравнения Раппопорта-Лиса [1]:
где s - водонасыщенность, / (з ) = к2 (з)/(к2 (з) + //0-кх (з)) - потоковая функция Ле-веретта; ; Ау = у -у ; V(1) = щ + щ - суммарная скорость фильтрации флюи-
дов, не зависящая от пространственной переменной; т - коэффициент пористости; к - абсолютная проницаемость среды; рк - капиллярное давление, ^, ц2 - коэффициенты динамической вязкости воды и нефти; к (з), к2 (з) - относительные фазовые проницаемости по воде и нефти.
Капиллярное давление связано с водонасыщенностью соотношением:
где а - коэффициент поверхностного натяжения, 0 - краевой угол смачивания, J(з) -функция Леверетта.
(1)
Перепишем уравнение (1) в дивергентной форме с выделенными конвективным и диффузионным слагаемыми:
m* = + K(s), (2)
dt дх дх ^ dx ) k (s)
где введены функции: K(s) = k 1 , (p(s) = f (s)•\V(t)-K(s)•Arsina].
Mi
Положим, что в (2) V (t) = const (V (t) < 0), и введем новые независимые переменные:
х = х/L, t = V■ t/(Lm) , где L - характерная длина. Перепишем (2) в безразмерном виде:
ds dF(s) д — + —— = s — dt дх дх
cm §], (3)
где F(s) = f (s)
, kArsina , . . 1---—■ki(s)
, G(s) = k(s)f(s)p'k(s), s = acos6>Jmkj(m LV)>0
малый параметр.
Пусть (х, г) - исходные независимые переменные. Переход к произвольной нестационарной системе координат осуществляется с помощью замены переменных общего вида х = 0(д,т), г = т, имеющей обратное невырожденное преобразование
д = 0 х (х, г) г = т. Частные производные зависимых переменных выражаются следующим образом:
д д ^дд д д д х 1 д ^ д (4)
дг дт дг дд дт дт у/ дд дт у дд
д дд дд 1 д д2 1 д
'l д^
- - , -2 , (5)
дх дх дд у дд дх У дд у дд где у = дх / дд - метрический коэффициент или коэффициент трансформации, показывающий во сколько раз изменяется исходная область; дх / дт - скорость движения системы координат, подлежащая в дальнейшем определению. Связав движение системы координат с особенностями решения, задаваемых в виде некоторой функции Q, получим уравнение обратного преобразования:
Iх. (6)
дт
Функция Q фактически является параметром управления движения узлов. В (6) возьмем производную д / дд, в результате получим:
У = _дд_ ^
дт дд
Используя замены переменных х = 0(д, т) и г = т, запишем уравнение (3) в новых переменных д, т представляя дифференциальные уравнения в дивергентной форме:
дs дЖ д(д ■ s) д^
дт дд дд дд
Ж = е
G(s) дз
у дд' дщ _ дQ дт дд ' д х дд '
(9) (10) (11)
где д0 < д < дд при х < х < хл, ^ > 0. Предполагается, что на каждый момент времени существует невырожденное преобразование:
х=#(д,т), в(д0,т) = хо , о(дк,т) = х. (12)
Вид функции Q в общем случае произволен, однако условия (12) налагают на него следующие требования Q(g0,т) = Q(g^г ,т) = 0. В работе [9] Q(g,т) предлагатся задавать в виде:
Q(g,т) = -А Щ-Qо -д
дд дд
¥
дз
дд
(13)
Первый член в (13) ограничивает чрезмерное сближение расчетных узлов сетки, а второй обеспечивает их сгущение в области больших градиентов водонасыщенности з . ^ и Б0 - положительные произвольные константы.
Аппроксимацию системы уравнений (7)-(11), (13) осуществим по неявной разностной схеме интегро-интерполяционным методом (методом баланса) [12]. Для этого в расчетном пространстве О г введем расчетную сетку с шагом к по переменной
д и Ат по переменной т :
Г \
с =
(д> т),
д 1,т'
гн— V 2
к
, к = дм -дг, д 1 = дг +-, Ат = т+1 -т, г = о,1,...,к-1, у = о,1,...
. 1 " Ъ о
2 2
К целым узлам сетки (д,т) отнесем переменную х/ и сеточные функции Ж/, Q/, ^, а к дробным узлам (дг+1/2, Т) - сеточные функции
у+1 ы./+1 (V) 1+1 V) 1+1
У/+У2 ' з/+1/2
зу+! - зЛ М^1 (V)-1
,+2 г+1 жг+1 - (Q• -(Q• 4 ^+1 -
„У+1 (V1
Ат
к
-¥ 1
гн— гн—
2 2
Ат
к
{уун1 (v)j■н1
Qг+ 1 - Qг 2_
к
к
Т71 +1 _ Т7 ' +1
Хг+1 хг _ 1+1
1 ¥ ■ 1 '
к гн—
2
. з1, - в1 , . а1 г+ г-— Ж' -^, г = 1,...,К-2,
¥ г
к
¥] 1 - У1' 1 У1' 1 •' - з/ - ¥] 1 •з/ - '
г н— г— гн— 1 1 г— 1 1
$ =-О0—^-^ - Qо—2--2
к
к2
(14)
(15)
(16)
(17)
(18)
Значения функций и я} в целых узлах определим через значение этих функций в дробных узлах по интерполяционной формуле:
\ ГУ. 1+ ГУ. 1
1л— I-2 2
У = —^-■
У 1 + У 1
I — г'н— 2 2
Полученная система разностных уравнений (14)-(18) нелинейна и может решаться различными способами, например с помощью итерационного метода Ньютона [10] или с помощью матричной прогонки с итерациями (у) по нелинейности [8,12].
Опишем применение итерационного метода Ньютона к системе уравнений (14)-(18). Решение при переходе с временного слоя } на } +1 представляется в виде:
У У у (19)
[} = я} +ёэ}
Записывая разностные уравнения (14)-(18) относительно неизвестных приращений 8я} и 8\, получаем на каждой итерации линейную систему алгебраических уравнений, решаемую методом последовательных прогонок. На первой итерации у = 1 из уравнения \ с помощью прогонки определяем приращение 8\}. Затем решая разностное уравнение для водонасыщености я} определяем приращение 8я}. Заканчивается итерация определением величин \у+1 = \ + 8\, = яу + 8яу и осуществляется переход на следующую итерацию. Итерирование продолжается до выполнения условий:
[К
= \\ -\1\<а\1+£2
8я"\ = ку+1 - ¿у|<£,-я"
10-6 <(^1, ¿>2, )<10-
Численный эксперимент. Для тестирования построенной схемы примем следу-
ющие параметры [13]: г (я) =
^ 2+(1-* )2
ОД = 4т(1 - я),
е = 0.01 ■
Начальные условия г = 0: я(х, 0) = {1 - 3х, 0 < х < ±} ^ {0, ± < х < 1}. Граничные условия х = 0: я(0, г) = 1. Расчетные параметры: число точек на отрезке 0 < х < 1 п = 215, шаг по времени т = 0,001. Результат моделирования на момент безразмерного времени г = 0,2 приведен на рисунке.
Рис. Профиль водонасыщенности на момент безразмерного времени г = 0,2
я
Заключение. Применение метода адаптивных разностных сеток, динамически связанных с решением к задаче двухфазной фильтрации позволяет избежать появления разрывных решений в произвольной нестационарной системе координат и повысить точность расчёта в областях с большими градиентами и разрывными решениями.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований № 15-05-08824а.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Данаев Н.Т., Корсакова Н.К., Пеньковский В.И. Массоперенос в прискважинной зоне и электромагнитный каротаж пластов. - Алма-Ата, Казахский ун-т. 2005. 180 с.
2. Аганин А.А., Гусева Т.С. Численное моделирование контактного взаимодействия сжимаемых сред на эйлеровых сетках // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. 2012. Т. 154. Кн.4. С. 74-99.
3. Гусева Т.С. Численное решение задач взаимодействия жидкости и газа на эйлеровых сетках без явного выделения контактных границ // Вестник Казанского технологического университета. 2013. Т. 16. №15. С. 135-140.
4. Бураго Н.Г., Кукуджанов В.Н. Обзор контактных алгоритмов // Изв. РАН. МТТ. 2005. - №1. С. 45-87.
5. Boris Jay P., Book David L. Flax-Corrected Transport, I. SHASTA, a Fluid Transport Algorithm that Works // Journal of Computational Physics. -1973. -Vol.11. -PP. 38-69.
6. Годунов С.К. Уравнения математической физики. - М.: Наука, 1971. - 392 с.
7. Дарьин Н.А., Мажукин В.И. Об одном подходе к построению адаптивных разностных сеток. // ДАН СССР. 1988. Т. 289. - №1. С. 64-68.
8. Дарьин Н.А., Мажукин В.И. Математическое моделирование нестационарных двумерных краевых задач на сетках с динамической адаптацией // Матем. Моделирование. 1989. Т. 1. -№3. С. 29-43.
9. Дарьин Н.А., Мажукин В.И. Метод построения адаптивных сеток для одномерных краевых задач. - М.: ИПМ АН СССР. 1987. 26 с. (Препр. / ИПМ АН СССР; № 33).
10. Мажукин В.И., Самарский А.А., Кастельянос О., Шапранов А.В. Метод динамической адаптации для нестационарных задач с большими градиентами // Матем. Моделирование. 1993. Т. 5, -№ 4. С. 32-56.
11. Бреславский П.В., Мажукин В.И., Такоева Л.Ю. Математическое моделирование лазерного плавления и испарения однородных материалов. Пакет LASTEC-1. -М.: ИПМ АН СССР. 1991. 46 с. (Препр. / ИПМ АН СССР; № 22).
12. Самарский А.А. Теория разностных схем. - М.: Наука, 1983. - 616 с.
13. Liu Y., Shu Chi-Wang, Zhang M. High order finite difference WENO schemes for nonlinear degenerate parabolic equations // SIAM J. Scientific Computing, 33(2). 2011. pp. 939-965.
© Д. С. Евстигнеев, 2017
