CC BY
28
Вычислительные технологии
Том 15, № 1, 2010
Регуляризация дискретной схемы для плоской задачи об эволюции границы раздела различных жидкостей*
Д. Н. Никольский Орловский государственный университет, Россия e-mail: nikolskydn@mail.ru
Для численного решения задачи об эволюции границы раздела различных жидкостей построена регуляризованная дискретная схема. Регуляризация проводится методом сглаживания ядра сингулярного интеграла, входящего в дифференциальное уравнение подвижной границы. Схема апробирована на конкретных задачах.
Ключевые слова: регуляризация сингулярного интеграла, эволюция границы раздела жидкостей, совместная фильтрация различных жидкостей.
1. Основная система интегрального и дифференциального уравнений
В работе [1] показано, что решение задачи об эволюции границы rt раздела жидкостей различных вязкоетей вне и внутри) в однородной безграничной пористой среде состоит в совместном решении уравнений
Д^ = 0 в R2\ rt, (1)
dr/dt = v0 + (v+ + v-) /2 на rt (2)
при граничных условиях
= , = на rt, (3)
— 0 при r —> oo (4)
и начальном условии
r = ro(0) при t = 0, в — некоторый параметр. (5)
Здесь <р — потенциал скорости фил ьтрации v; r — радиус-век тор; t — время; знаками "+" и "—" обозначены предельные значения соответствующих функций при подходе к границе rt против и вдоль нормали n; — потенциал скорости возмущения, вызванного внесением границы раздела жидкостей rt в область протекания процесса. Выражения (1)-(5) записаны в безразмерных величинах.
Используя потенциал двойного слоя, решение задачи (1) и (2) при условиях (3)-(5) сводится к системе интегрального и дифференциального уравнений:
g — 2AG*g = 2А^0, dr/dt = v0 + VG*g на rt, (6)
*Работа выполнена при финансовой поддержке Президента РФ (проект МК-491.2008.1). © ИВТ СО РАН, 2010.
где С*д = J д(^ £)П(М, N— потенциал двойного слоя плотности д; П = (Е, п^); г4
Е = VNФ(М, N); Ф — потенциал нормированного стока, мощность которого равна — 1; _ потенциал скорости невозмущенного течения, параметр А = (^2 — ^1)/(^1 + ^2)-В каждый момент времени tj, j = 0, ,1 разобьем границу Г^. на п частей и получим множество точек {г3т\т = 0,п — 1}, Элементы этого множества при £ = 0 вычисляются по формуле (5), Методом дискретных вихревых пар выполним дискретизацию основной системы (6) [1], в результате получим итерационную последовательность для системы линейных алгебраических уравнений (СЛАУ) и выражение для вычисления элементарных смещений частиц границы раздела жидкостей:
?т+1/2 = Адт+1/2 + /т+1/2> /т+1/2 = 2Ат^0 т+
gm+l/2 — Agm+l/2 + fm+l/2' fm+1/2 — m+1/2'
n— 1
Agm+ 1 /2 — ^ gfc+1/2 (Fm+1/2,fc+1/2' nfc+1/^ A/k+1/2 + (1 — T)gl+1/2>
k=0 k=m
n—
(7)
Arm/Atj — V0m + +1/2 (Vjmfc - Vjm,fc+0 ' — ro(^m),
k=0 k=m
j — 0, J, m — 0,n — 1, i — 1,1,
где пк+1/2 = Т 1к+1/2/А1к+1/2; Т = о^ _ матрица поворота на 90° против часовой стрелки (единичные векторы касательной и нормали образуют правую двойку);
A/j A/fc+1/2
fc+1/2
lk+1/2 — rk+1 — rl rk+1/2 — (ri + 4+0А причем rn = r0; i
число итераций для решения СЛАУ, которое определяется из условия || g1 — g1 —< (1 — ||А||)/||А||б, e — заданная точность; V£mfc — ©£Vmfc, ©£(r) — в1(г/г£) — сглаживающая функция [2], удовлетворяющая следующим условиям: в1(0) — ©1(0) — ©1(0) — 0, ©'/'(r) < Cr2 для V r £ [0, то), r£ — эффективный радиус вихря, C — константа, ©1'(r) — 1 при r > 1; Arm — rm+1 — rJm] Atj — tj+1 — tj V — скорость нормирован-
-1T
обеспечивающий выполнение необходимого и достаточного условия сходимости итерационной последовательности gl,i = 1,1 из (7): р(А) < 1 [3], здесь р(А) — спектральный радиус матрицы A
Схема (7) отличается от ранее исследованных дискретных схем наличием сглажи-©
скорости вихря каждый отрезок rm,m+1 смещается за его концы, а не за центр. Ниже рассмотрены практические задачи, в которых применение данной дискретной схемы позволяет произвести расчет. Для численных расчетов использовался один из возможных вариантов сглаживающей функции ©£ — (63(r/r£)5 — 90(r/re)7 + 35(r/re)9) /8 при r < rf.
2. Эволюция в поступательном потоке
Рассмотрим эволюцию границы Г раздела различных жидкостей, представляющей собой в начальный момент времени £ = 0 окружность радиуса Я с центром в точке (хс, ус)
выбранной декартовой системы координат xOy, Движение границы раздела происходит под действием поступательного потока, скорость которого па бесконечности равна uex, Так как область протекания процесса без гран ична, функции <р0, v0, F и V из дискретной схемы (7) примут вид |4|
^0 = ux, vo = uex, F = (2n)-1r nm/t2Nm , V = (2n)-1r nm/t2Nm ■
В начальный момент времени t = 0 можно точно вычислить скорость смещения частиц границы Г0, Действительно, применяя теорему об окружности [4], получим потенциалы, описывающие течение вне и внутри окружности:
^ = (1 + XR2/(x2 + y2))ux, = (1 - A)ux. (8)
Непосредствеппой подстановкой убеждаемся, что потенциалы (8) удовлетворяют уравнению Лапласа (1) и граничным условиям (3).
Поело дифференцирования (8) получим точное значение скорости смещения частиц границы Г0:
va = (1 — Xx2/R2) uex — Xuxy/R2ey на Г0. (9)
Множество точек, моделирующее окружность Г0, получим, используя ее параметрические уравнения:
Х<пг = хс + R cos tm, Ут = Ус + R sin tm, 7П = 0, П ~ 1,
t = {тг{')п/с)р\ m = 0, с — 1} и {-7г(т/с)р| m = с- 1,0} , с = п/2
Здесь параметр p характеризует неравномерность разбиения окружности. При p = 1 и p = 1 имеем систему точек, разбивающую окружность Г0 соответственно на равные и неравные но длине части.
Введем погрешность nn = |1 — v/va| 100% оде v — модуль скорости смещения частиц границы Г0, полученный численно из системы (7), va — модуль вектора скорости, вычисленный по точной формуле (9), Для определенности выберем R = 1, u = 1, X = 0.5, xc = 0 и yc = 0,
Рис. 1. Зависимость погрешности п200 (я) и п400 (б) от радиуса вихря r£; p = 1 (i), 2 (2), 3 (3), 4 (4), 5 (5)
(10)
На рис, 1 представлены зависимости погрешности вычислений цп от радиуса вихря ге для п = 200 и п = 400, Анализируя рисунки, видим, что при разбиении границы раздела жидкостей Г4 на неравные по длине части результаты численных расчетов имеют большую погрешность. Однако эта погрешность может быть устранена путем использования сглаживающей функции ве из (7) и подбора оптимального значения радиуса вихря ге,
В табл. 1 представлена зависимость погрешности пп от степени неравномерности разбиения окружности р при оптимальном значении радиуса вихря ге для числа точек разбиения п = 200, 400, 600, 800, Видно, что с увеличением п при оптимальном значении эффективного радиуса вихря ге погрешность численного решения задачи уменьшается при р = 1 и 2 для всех п, а при р =3, 4 и 5 для п = 200 и 400,
На рис, 2 показаны положения границы раздела жидкостей Г при £ =1.8 для случаев ге = 0 (без сглаживавия) и ге = 0.064 (со сглаживанием), В начальный момент времени множество точек границы раздела жидкостей Го было получено из (10) при р = 1, т.е. при £ = 0 точки разбивали границу раздела жидкостей на равные по длине дуги части. Затем, с течением времени, в силу касательных смещений частиц границы это множество стало образовывать неравные по длине отрезки, что при ге = 0 привело к неустойчивым численным расчетам и не физическому разрыву границы раздела жидкостей (рис. 2, о).
Таблица 1. Погрешность пп при оптимальном значении ради уса вихря те
V 1 2 3 4 5
Ге 200 0 0.0137 0.0021 0.0069 0.0077
Г]200,% 0.7682 1.6253 2.5511 3.5296 4.5476
Ге 400 0 0.0102 0.017 0.0227 0.0028
Щоо,% 0.3796 1.1885 1.4668 1.8086 2.1244
Ге 600 0 0.0082 0.0204 0.0291 0.0330
Г]600,% 0.2520 0.9601 1.5486 1.9858 2.2762
Ге 800 0 0.0071 0.0234 0.0306 0.0357
Г]800,% 0.1886 0.8179 1.6049 2.0945 2.3675
0.5
-0.5
0 0.5 1 1.5 X
а
0.5
-0.5
0 0.5 1 1.5 X
б
Рис. 2. Граница раздела жидкостей при те = 0 (а) и те = 0.064 {б)
0
0
V 12
10
8
6
4
2
0
О 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 t
Рис. 3. Изменение числа обусловленности v с течением времени t\r£ = 0 (i), 0.064 (2)
Для наглядности выбрано n = 200, Оптимальное значение параметра т = 0.8, На каждом временном шаге для решения СЛАУ из системы (7) с точностью е = 10-6 выполнялось 13-14 итераций. Расчет для случая т£ = 0 был выполнен с приближенным условием остановки итерационного процесса из (7): ||д1 — д1-1\\ < е.
На рис, 3 показано изменение числа обусловленности v матрицы СЛАУ из (7) с течением времени t. Линия 1 на рисунке соответствует расчету без сглаживания при т£ = 0 линия 2 — расчету со сглаживанием при т£ = 0.064, Видно, что при достижении некоторой степени неравномерности, в случае, когда сглаживание пе выполняется, число обусловленности v резко возрастает. Последнее не позволяет считать результаты расчетов достоверными и объясняет отсутствие физического смысла у результатов расчета, представленных па рис, 2, а. Применение регуляризации позволило избежать
v
логичные выводы при исследовании задачи Лейбепзопа получены в |5|,
Отметим, что согласно (8) скорость движения частиц жидкости, находящихся внутри границы раздела жидкостей, постоянна и равна v2 = (1 — A)uex, В исследованной задаче модуль средней скорости движения частицы жидкости, имеющей в начальный момент времени координату (0,1), равен 0.513 и превышает v2 = 0.5 на несколько процентов, что согласуется с данными |6|,
3. Эволюция к стоку
Центр первоначальной окружности Го радиуса R поместим в точке (а, 0), Невозмущенное течение моделируем стоком, расположенным в начале координат. Его потенциал и скорость v0 равны
^о = q(2n)-1 ln тм, vo = q(2n)-1rM/т2м,
где q — дебит скважины (соответствует стоку при q < 0),
При параметре А = 1 задача об эволюции стока к окружности имеет точное решение
(а, 0)
t
ЖС ) = АС/(1 — atZ ) + YtZ,
(11)
причем коэффициенты at, в и It определяются системой уравнений
(.R2 - qt/ir)2w3 - (2а2Д2 + 2a2qt/n + a4) w + 2а4 = 0, at = ±Vw,
et = (1 - at2) (a/at + at/a (R2 - qt/n)) /2, (12)
Yt = (a/at - at/a (R2 - qt/n)) /2.
Полагаем R = 1, a = 0.5, q = -n, t = 0.1. Численно решая систему (12), получим параметры ao.i = 0.56991 во.1 = 0.71947, Yo.i = -0.18824. При решении кубического уравнения был выбран корень, имеющий физический смысл. Подставляя найденные параметры в (11), находим координаты частицы, смещающейся вдоль направления прорыва воды в скважину в момент времени ta = 0.1 — (-0.27004, 0.0).
Определим время достижения границей раздела жидкостей rt точки (-0.27004, 0.0) численно путем решения системы (7) при Л =1. Первоначальную границу раздела жидкостей Г0 моделируем системой точек, вычисленией из (10) при xc = 0.5 и yc = 0.0. В табл. 2 приведена зависимость величины цп = |1 - t/ta|100% от неравномерности разбиения первоначальной окружности p. При подборе г£ выполнялся визуальный контроль за границей раздела жидкостей. При этом значения г£, соответствующие деформации границы раздела жидкостей, не имеющие физического смысла и не совпадающие с аналитическим решением, были отброшены.
Анализируя табл. 2, видим, что при разбиении начальной границы раздела жид-
p = 1
требуетея и численное решение сходится к точному с ростом числа и, при разбиении Таблица 2. Влияние неравномерности разбиения на погрешность Пп при dt = 0.002 и 0.001
V 1 2 3
dt = 0.002
Ге 200 0 0.038 0.036
Г]200,% 3.7223 5.1999 5.1855
Ге 400 0 0.030 0.031
Щоо,% 2.9215 4.4799 4.5013
Ге 600 0 0.028 0.028
Г]600,% 2.2861 4.4199 4.4509
Ге 800 0 0.026 0.028
Г]800,% 1.6953 4.3746 4.4307
dt = 0.001
Ге 200 0 0.037 0.036
Г]200,% 3.58345 5.0499 5.0468
Ге400 0 0.030 0.031
Щоо,% 2.23943 4.3413 4.3626
ГебОО 0 0.028 0.028
Г]600,% 1.87199 4.0107 4.3047
Ге 800 0 0.026 0.028
Г]800,% 1.55051 3.9655 4.2922
на неравные по длине части регуляризация необходима и позволяет добиться практической сходимости численного решения к точному.
Из проведенных вычислительных экспериментов следует, что при численном решении задач об эволюции границы раздела различных жидкостей методом дискретных особенностей в случае разбиения подвижной границы на неравные по длине части в силу отсутствия в дискретной схеме механизмов саморегуляризации возникает и быстро нарастает численная неустойчивость, приводящая к результатам, не имеющим физического смысла. Эта численная неустойчивость устраняется методом сглаживания ядра сингулярного интеграла в дифференциальном уравнении движения границы раздела жидкостей. Необходимость сглаживания, а также ее параметры (в нашем случае вид функции в£ и радиус вихря г£) предлагается устанавливать в ходе численного эксперимента путем визуального наблюдения за подвижной границей и контроля характеристик дискретной схемы (7), например, числа обусловленности СЛАУ,
Список литературы
[1] Никольский Д.Н. К вопросу построения дискретной схемы для плоской задачи эволюции границы раздела различных жидкостей // Вычисл. технологии. 2008. Т. 14, № 4. С. 89-94.
[2] Кирякин В.Ю., Сетуха A.B. О сходимости вихревого численного метода решения трехмерного уравнения Эйлера в лагранжевых координатах // Дифференц. уравнения. 2007. Т. 43, № 9. С. 1263-1276.
[3] Формален В.Ф., Ревизников Д.Л. Численные методы. М.: Физматлит, 2006. 400 с.
[4] Голувева О.В. Курс механики сплошных сред. М.: Наука, 1971. 368 с.
[5] Никольский Д.Н., Дорофеева В.И. Исследование дискретных схем для задачи об эволюции границы раздела жидкостей в постановке Лейбензона // Сб. Междунар. школы-семинара "МДОЗМФ-2008". Вып. 6. Орел: Картуш, 2008. С. 73-77.
[6] Воинов В.В. О точных решениях задачи движения границы раздела несмешивающихся жидкостей в пористой среде // Прикл. механика и техн. физика. 1991. № 1. С. 68-71.
[7] Варченко А.Н., Этингоф П.И. Почему граница круглой капли превращается в инверсный образ эллипса. М.: Наука, 1995. 80 с.
Поступила в редакцию 20 апреля 2009 г., в переработанном виде — 18 августа 2009 г.
