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

УДК 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