К вопросу построения дискретной схемы для плоской задачи эволюции границы раздела различных жидкостей Текст научной статьи по специальности «Физика»

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

Том 13, № 4, 2008

К вопросу построения дискретной схемы для плоской задачи эволюции границы раздела различных

жидкостей

Д. Н. Никольский Орловский государственный университет, Россия e-mail: nikolskydn@mail. ru

In this work, a plane problem on an evolution of the interface between various liquids is numerically solved in frames of the Leibenzon—Muskat (piston replacement) formulation. Compared to the previous studies, an easily realizable discrete scheme is proposed.

1. Постановка задачи

Рассмотрим нестационарную фильтрацию в однородной изотропной пористой среде. Область течения D, не нарушая общности рассуждений, полагаем безграничной (границы области фильтрации могут быть учтены при помощи функции Грина). В некоторый момент времени t резкая подвижная граница rt делит область D на части D-(внутренняя) и D+ (внешняя), занятые жидкостями с вязкостями и соответственно. Границу rt полагаем гладкой. Нормаль n к ней направляем в область D+. Течение жидкости в области D возбуждается стоками (источниками).

Напорную фильтрацию жидкостей в каждой точке M = (x,y) £ D описывает потенциал р = p(M, t), удовлетворяющий уравнению Лапласа [3]:

Др = 0, р = -p/^,2 в D\rt. (1)

Здесь p — давление.

Потенциал р связан со скоростью фильтрации v = dr/dt = r следующим соотношением [3]:

v = Vp. (2)

Законы (1)-(2) представлены в безразмерных величинах. При этом характерные величины удовлетворяют соотношениям

Фо^о = KPo, V0L0 = Фо, aLo = VqTQ.

Здесь Ф0 — характерный потенциал, ^0 — характерная вязкость, K — коэффициент проницаемости пористой среды, P0 — характерное давление, V0 — характерная скорость, L0 — характерное расстояние, a — пористость, T0 — характерное время.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-96303), Президента РФ (грант № МК-491.2008.1) и Федерального агентства по образованию РФ.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

Полагаем, что вытеснение одной жидкости другой "поршневое" (модель Лейбензо-на—Маскета). В любой момент времени £ граница Г представима параметрически: г = г(£, 0). Ее перемещение описывается дифференциальным уравнением [1, 2]:

г= + V-) /2. (3)

Здесь и далее знаками "+" и " —" обозначены предельные значения соответствующих функций при приближении к границе Г из областей и . В начальный момент времени £ = 0 граница Г задана

г0 = г(О,0). (4)

На границе Г в любой момент времени £ выполняются условия непрерывности давлений и расхода жидкостей [3]:

р+ = р-, = V- на Г.

С учетом (1) и (2) эти условия примут вид

= = (0пр)- на Г4. (5)

Здесь дп — оператор (частного) дифференцирования. Потенциал <р(М, £) представим в виде

= ^о + , (6)

где — потенциал невозмущенного течения однородной жидкости, моделирующий стоки (источники), — потенциал возмущений, вызванных различием вязкостей жидкостей.

Условия (5) для потенциала примут вид

(1 — А)^+ — (1 + А)<£- = 2А^о, (дп^*)+ = (дп^*)- на Г*. (7)

Здесь параметр А = (^2 — ^1)/(^2 + ^1).

Для обеспечения единственности решения задачи на потенциал наложим условия регулярности [4]:

= 0(1/т^м), |У<р*| = 0(1/т2им), гИы ^ ж.

Дифференциальное уравнение движения границы Г (3), с учетом (6), примет вид

г = Vо + V*, Vо = У^о, V* = ((Ур*)+ + (У^*)-) /2. (8)

Итак, плоская задача об эволюции границы Г раздела различных жидкостей состоит в отыскании ее положения в любой момент времени £ > 0 из уравнений (1) и (8) для потенциала и радиуса-вектора г, при начальных и граничных условиях (4) и (7), с учетом условий регулярности.

2. Основная система интегрального и дифференциального уравнений

Потенциал р* ищем в виде потенциала двойного слоя, непрерывно распределенного по границе Г с плотностью д:

Р* = у g(N,t)ü(M,N)dlN, Ü(M,N) = dnNФ!(М,Ж), M / Г4, rt

где Ф1(M, N) — потенциал нормированного стока с расходом П = -1. Учтем известные свойства потенциала двойного слоя [5]:

= G*g ± g/2, (Vp*)± = VG*g ± Grad g/2 на rt.

Здесь в интегральном операторе G*

G*g = J g(N,t)Q(M, N)dlN

rt

(9)

(10)

(11)

интеграл понимается в смысле главного значения по Коши [5], Grad g — касательная составляющая градиента функции на кривой. Подставим (9) в (7) и (8), учтем (10), имеем

g — 2\G*g = 2Аро, r = vo + VG*g на Г4.

(12)

Таким образом, решение задачи о нахождении границы Г при £ > 0 сведено к решению системы интегрального и дифференциального уравнений (12) при заданном начальном условии (4).

3. Дискретная схема для основной системы

Несложно показать, что

0пмФ1 (М, N) = -01мФ2(М, N).

Здесь Ф2(М, N) — потенциал скорости нормированного вихря с циркуляцией Г2 = — 1, расположенного в точке N [3]. Подставляем последнее выражение в оператор (11), имеем

0*д = — У д^,1)31мФ2(М, N)Ш„. (13)

г4

Интегрируя по частям (13) [1], получим дифференциальное уравнение из системы (8) в виде

Г = vo + У 01м д^,г)У2(М^ )(Ц„, (14)

г4

где У2(М, N) = УмФ2(М, N) — скорость нормированного вихря.

В момент времени ¿^, 3 = 0,1,... , границу Г^. разобьем на п частей. Тогда ее положение будет определяться множеством точек ^(ж^Ут), т = 1,п, 3 = 0,1,... При этом начальное условие (4) примет вид

Гт,(хт,>Ут,), m 1,п-

Из (12) и (14) имеем систему линейных алгебраических уравнений (СЛАУ) и разностный аналог дифференциального уравнения движения границы

gm - 2А £ gk jAj = 2A^0m,

k=1 k = m

j = dxk ^m^m ,yk Kk+Byk ф^у- ,xk ^k Kk ,

n

Аг3 ___

= ^ + £ ))к лт, т = 1,п, з = 0,1,... (15)

к=1

В случае если на каждом шаге по времени граница Г4 переразбивается на равные по длине части, производные пХк, пУк и ) в (15) заменяются на разностные аналоги с

помощью центральных разностей [1, 2]. При этом для восстановления уравнения кривой Г в качестве интерполянта используются кубические и линейные сплайны. В случае если граница Г не переразбивается, для производных пХк, пУк и ) из (15) строятся

более точные разностные операторы [6].

Недостаток дискретной схемы (15) — необходимость аппроксимации производных пХк, пУк и ). От него можно избавиться. Для этого запишем дискретный оператор,

соответствующий (13), в виде

1к+1 П &

Gm ~ - gj 9in ф2(хт, ym >N )d/N

к=1

к-1

п

= Е9к (ф2т;к-1 - ф2т,к+1) ' 3 = 0,(16) к=1

С учетом (16) дискретная схема для основной системы (12) примет вид

п

9тт - 2^ 9к ^Ф2т,к— 1 - Ф2т,к+1) = 2А^0т, к=1

Ar-7

vo>

^f = Vom + £ gk (Vj^k-1 - V2m,k+ f) , ™ =1,П, j = 0, 1,... (17) k=1

Для решения (17) нет необходимости строить разностные аналоги производных.

Отметим, что если вычисление потенциала вихря представляется невозможным или неоправданно сложным (например в двумерном случае), то для решения задачи эволюции границы раздела жидкостей можно решать СЛАУ из (15) совместно с разностным аналогом дифференциального уравнения из (17).

n

n

4. Сопоставление с аналитическим решением

Границу Го моделируем окружностью радиуса R. Ее параметрическое представление

x = R cos в, y = R sin в, в = (2п, 0].

Потенциал невозмущенного течения р0(М, ¿) моделируем стоком, расположенным в точке М^х^ у1):

Я 1

Ро = — 1п ГМхМ, 2п

где я — дебит скважины (соответствует стоку при д < 0).

В этом случае на первом шаге по времени можно вычислить аналитически скорость перемещения границы раздела жидкостей. Сток расположим внутри окружности (х1 + у2 < Л2). Следуя [3], применяем метод изображений к внутренней (для стока в М\) и к внешней (для источника в бесконечно удаленной точке) областям окружности:

= 2П ((! + Л)1п гМхМ - А 1п ГМ) , Р2 = 2П (1п ГМ1М - А 1п тмгМ + А 1п , М1 ' (18)

Потенциалы (18) удовлетворяют условиям непрерывности давлений и расхода жидкостей (5) на Г0.

Используя (18), точно вычислим скорость смещения границы Г0:

V« = V = (1 + А) ^ - А + тМ^У (19)

V 2 ) V V тМ1М 2 \т" ткм)

При указанных условиях решение СЛАУ из систем (15) и (17) значительно упрощается. Так, после несложных эквивалентных преобразований СЛАУ из (15) имеем рекуррентные соотношения

n

сА 1 + - PEC*

у 1 -а ) 1 -a t2 , Ck - Ci

gi = -7-,-f-ñ-' gk = gi + ^-:

1 + a(n - 1) 1 - a

A

для СЛАУ из (17)

Ck = 2Ap0k, k =1,n, a =--

n

Ci (1 + a(n - 1)) - a£Ck

gi =--——, gk = gi + Ck - C1

1 + an

A

Ck = 2Ap0k, k =1,n, a =--

n

В таблице представлена зависимость ошибки ni и n2 численных расчетов скорости смещения Го в момент времени t = 0 по схемам (15) и (17) для параметра А = 0.5. Величины

nv = |1 - VvKI- 100%, V =1, 2.

Сопоставление результатов численного счета с аналитическим решением

n 50 100 200 400 800 1600

ni, % 3.9503 1.9245 0.9332 0.4586 0.2272 0.1131

П2, % 0.1502 3.943 ■ 10-2 9.866 ■ 10-3 2.467 ■ 10-3 6.168 ■ 10-4 1.542 ■ 10-4

Здесь vi,2 — модули скорости смещения границы Го, вычисленные по дискретным схемам (15) и (17) соответственно.

Из таблицы видим, что с ростом n ошибки ni,2 уменьшаются, что свидетельствует о

сходимости дискретной схемы. Также видим, что схема (17) более точная.

Список литературы

[1] Никольский Д.Н. Математическое моделирование работы системы скважин в однородных и неоднородных слоях с подвижной границей раздела жидкостей различной вязкости. Дис. ... канд. физ.-мат. наук. Орел, 2001. 191 с.

[2] Lifanoy I.K., Nikolsky D.N., Piyen' V.F. Mathematical modelling of the work of the system of wells in a layer with the exponential law of permeability variation and the mobile liquid interface // Russ. J. Numer. Anal. and Math. Model. 2002. Vol. 17, N 4. P. 381-391.

[3] ГолУБЕВА О.В. Курс механики сплошных сред. М.: Наука, 1971. 368 с.

[4] Пивень В.Ф. Единственность решения граничных задач сопряжения физических процессов в неоднородной среде // Тр. X Междунар. симп. "Методы дискретных особенностей в задачах математической физики". Херсон, 2001. С. 265-269.

[5] Довгий С.А., ЛиФАНОВ И.К. Методы решения интегральных уравнений. Киев: Наукова думка, 2002. 286 с.

[6] Никольский Д.Н. Математическое моделирование трехмерной задачи эволюции границы раздела жидкостей различной вязкости и плотности в однородном и безграничном грунте // Вычислительные методы и программирование. 2006. Т. 7, №. 2. C. 108-114. (http://num-meth.srcc.msu.su).

Поступила в редакцию 5 сентября 2007 г., в переработанном виде —16 ноября 2008 г.