О численном моделировании задачи Маскета со свободной границей Текст научной статьи по специальности «Физика»

УДК 517.958:531.72, 517.958:539.3(4)

О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ЗАДАЧИ МАСКЕТА СО СВОБОДНОЙ ГРАНИЦЕЙ ^ О.В. Гальцев, О.А. Гальцева

Белгородский государственный университет,

ул. Победы, 85, Белгород, 308015, Россия, e-mail: [email protected]

Аннотация. Работа посвящена задаче фильтрации двух несмешивающихся несжимаемых жидкостей различной плотности, разделенных свободной границей. Среди математических моделей совместного движения двух несмешивающихся жидкостей хорошо известна задача Мас-кета. В работе приведены результаты численной аппроксимации известной математической модели Маскета, используемой для описания фильтрации указанного типа, точными микроскопическими моделями со свободной границей для структуры порового пространства в виде изолированных капилляров.

Ключевые слова: задача Маскета, задача со свободной границей, фильтрация жидкости.

1. Введение. В работе рассматривается численная аппроксимация модели Маскета [9], описывающей совместное движение двух несмешивающихся жидкостей в пористой среде, например, вытеснение нефти водой. Она является одной из большого числа моделей, предназначенных для решения указанной задачи описания и удовлетворяющих основным ограничениям, вытекающим из её физической постановки. Модель Маскета описывает фильтрацию несмешивающихся несжимаемых жидкостей различной вязкости и различной плотности, разделенных движущейся границей (свободной границей).

Движение одной из жидкостей, имеющей постоянную вязкость р+ и постоянную плотность р+, происходит в области О+(1). Оно описывается системой уравнений фильтрации Дарси

к

v+= —-Ур| + р+^, У-г;+ = 0, хеП+(*), (1)

р+ J J

для макроскопической скорости течения этой жидкости и её давления .

Соответственно, движение другой жидкости в области О-{Ь), имеющей постоянную вязкость р- и постоянную плотность р-, описывается аналогичной системой уравнений

к

=-----V • V- = 0, хеП"(*), (2)

р

для её скорости течения V- и давления р-.

1 Работа выполнена в рамках ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009-2013 годы (госконтракт к 02.740.11.0613).

Будем полагать, что на свободной границе Г^) = дО+ Р| дО (t), разделяющей обе жидкости, давления и нормальные скорости непрерывны:

р+ = р- х е ГС0, (з)

v+ ■ п = V- ■ п = УП, х е Г(t), (4)

где п есть единичный вектор нормали к границе Г^) в точке х е Г(t) и Уп есть скорость границы Г(t) в направлении нормали п в точке х е Г(t).

Условие (4) означает, что граница Г(t) представляет собой поверхность, которая в

течение всего движения состоит из одних и тех же материальных точек. Это условие

позволяет сформулировать понятие слабого решения задачи Маскета. А именно, определим давление р неоднородной жидкости

р/ = р+, если х е О^^, и р/ = р-, если х е О-^) ,

и скорости V

V = v+, если х е О^^, и V = V-, если х е О-^).

Тогда неизвестные функции V, р/ и р/ будут решениями системы уравнений фильтрации Дарси

к

v =-----Vpf + pfF, V • V = 0, ж е П , £>0, (5)

р

и уравнение переноса

^ _ ур/ = о, жбО, £>0. (6)

^ 6 д 6

Первое уравнение в (5) (закон Дарси) понимается в обычном смысле почти всюду в От = О х (0,Т), а второе уравнение (уравнение неразрывности) понимается в смысле теории распределений, как интегральное тождество

/ V -V ф йхдьЪ = 0

о Пт

для произвольной гладкой функции ф. Аналогично, понимается уравнение переноса (6). А именно, используя равенство

р v ■ V р/ = V ■ (v р/ р) - р/ v ■ V р - р/ р V ■ v , и уравнение неразрывности получим, после умножения уравнения (6) на произвольную гладкую финитную в области От функцию р и интегрирования по частям, интегральное тождество

[ йр

Рг -г- а хж = 0 ,

Упт ^

которое выполняется для всех таких функций р.

Постановка задачи дополняется однородным граничным условием

V ■ п = 0, х е 5 = дО , t > 0, (7)

где п есть нормальный вектор к границе 5, и начальным условием

р/(х, 0) = р0(х), х е О , (8)

где

0 (х

и

р0 (ж) = р+ = const > 0 , если ж G П+ = П+(0) ,

р0 (ж) = р- = const > 0 , если ж G П- = П-(0).

Развитие методов математического анализа в последние десятилетия позволило решить вопрос об адекватности некоторых моделей подземной гидродинамики. На главный вопрос о выполнимости закона Дарси был дан положительный ответ Л.Тартаром (Ь.ТаЛаг, 1980г.). Он изучил движение вязкой жидкости, описываемое на микроскопическом уровне уравнениями Стокса

aмДv — Vр + р^ = 0, V- V = 0 (9)

(движение в периодических порах размера е) и, используя методы теории усреднения, строго вывел закон Дарси при е \ 0 для случая абсолютно твёрдого скелета грунта. Таким образом, математическая модель фильтрации жидкости, описываемая системой уравнений Дарси, нашла своё обоснование на микроскопическом уровне.

Естественно ожидать, что аналогичным образом может быть обоснована модель Маскета. С этой целью, задачу со свободной границей о совместном движении в периодических порах размера е двух несмешивающихся жидкостей необходимо рассмотреть на микроскопическом уровне и затем перейти к усреднённой задаче с последующим предельным переходом при е \ 0. В нашей работе мы попытались реализовать эту схему численно на примере задачи о неустойчивости Рэлея-Тейлора, рассмотрев движение двух жидкостей под действием силы тяжести.

В классической гидродинамике идеальных жидкостей, если более тяжелая жидкость находится сверху, то наблюдается так называемая неустойчивость Рэлея-Тейлора. При её изучении прослеживаются следующие ярко выраженные стадии: линейная, промежуточная, регулярная, асимптотическая и турбулентная [2],[3]. Неустойчивость Рэлея-Тейлора наиболее исследована для случая плоской поверхности раздела и стремящегося к бесконечности отношения плотностей тяжелой и легкой жидкостей. Линейная стадия подробно изучена в классических работах Рэлея, Тейлора и Льюиса [4],[5],[6], регулярная асимптотическая - в работах Биркгофа [7], в [8] развита феноменологическая теория турбулентной стадии, а в [3] высказаны некоторые соображения о механизме её образования.

Однако, аналитического математического аппарата для анализа в целом рэлей-тей-лоровской неустойчивости в идеальной жидкости недостаточно. Экспериментальные

же исследования весьма трудоемки. В такой ситуации наиболее полная информация об этой неустойчивости может быть получена на основе численных расчётов.

В случае медленных движений вязких жидкостей, описываемых системой уравнений Стокса, задача о поведении границы раздела двух жидкостей изучена аналитически [1]). Вне зависимости от того, какая жидкость находится сверху, граница раздела остаётся гладкой и не теряет свою регулярность.

Для макроскопической модели Маскета никаких результатов аналитического или численного анализа неустойчивости Рэлея - Тейлора автору не известно (есть несколько работ о существовании классического решения в малом отрезке времени или глобального решения по времени, находящегося вблизи некоторого точного решения [4],

[12], [14]).

В настоящей статье приведены результаты, полученные при решении системы уравнений

а^ Ди — Vp — ре2 = 0 , V • и = 0 ,

в отдельно взятом капилляре. Это решение есть аппроксимация задачи Маскета со свободной границей для структуры порового пространства в виде изолированных капилляров.

Расчёты двумерного течения при наличии рэлей-тейлоровской неустойчивости показали хорошие совпадения с результатами [13]. На рис. 1 приведены хорошо различимые формы контактной поверхности в различные моменты времени

*=0. 003 *=0. 255 *=0. 525 *=1. 035

Рис. 1. Форма контактной поверхности в различные моменты времени при двумерной рэлей-тейлоровской неустойчивости

При увеличении количества капилляров (е \ 0) чёткой границы раздела двух жидкостей не наблюдается. При этом становится различимой переходная фаза (рис.2) р- < р < р+.

*=0.623

*=3.562

*=8.560

Рис. 2. Периодические изолированные капилляры

Заметим, что при задании начального возмущения в расчетах со сжимаемыми средами важно везде соблюдать условие V ■ и = 0 (за исключением поверхности раздела), иначе возникнут возмущения, которые могут исказить картину развития неустойчивости Рэлея-Тейлора.

2. Вычислительный алгоритм

Для численного моделирования задачи Маскета был выбран случай абсолютно твёрдого скелета грунта. А именно, на микроскопическом уровне была рассмотрена задача со свободной границей совместного движения в периодически расположенных порах размера е двух несмешивающихся жидкостей. Поры были выбраны в виде изолированных капилляров.

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

и уравнение переноса

при условии, что

ам Ди — Vp — ре2 = 0 и=0

(10)

(11)

(12)

р(х, 0) =

р+, Х2 > #0^1) р- , Х2 < #0(Х1)

Система уравнений (2.1)-(2.4) дополняется однородными условиями:

р+ = р- , х е Г^) , v+ ■ п = V- ■ п = УП , х е Г^)

V ■ п = 0, х е 5 = дО,

t > 0.

.к?

П+(4)

П-(*)

Рис. 3. Область 0±(£)

При моделировании задачи использовался модифицированный метод расщепления. Выбор данного метода обусловливается тем, что его разностная схема позволяет рассчитывать течения вязкой несжимаемой жидкости без использования граничного условия для вихря на твёрдой поверхности и при всех прочих равных условиях обладает большей эффективностью.

Расчёты для одного капилляра размера е проводились в области масштаба 1 х 200. При этом очень важно, что для фиксированного е > 0 (е - характерный размер поры) доказано [1], что существует единственное решение задачи, в которой граница раздела двух жидкостей есть липшицева поверхность.

Основная трудность при численном решении системы уравнений (2.1),(2.2) связана с расчётом поля давления. Первый значительный успех в преодолении отмеченной трудности был достигнут благодаря использованию идеи искусственной сжимаемости [12]. Существо этой идеи состоит во введении в уравнение неразрывности дополнительного д

члена вида — (р+и /2). В результате получается модифицированная система уравнений вида * ^

др

777 + V

— Ур — ре2 = 0

ді

и

0.

В нашем случае уравнение неразрывности будет иметь вид

0.

др

— + ІУсііу и ді

(14)

В отличие от стандартного метода расщепления, мы не будем рассчитывать отдельно промежуточное поле скоростей и подправлять его с учётом градиента давления, а, заменив на начальном шаге по времени величину р в уравнении (2.1) выражением ^т£0, где £0 - малая функция, получим легко решаемое уравнение

а^Ди + т^Уо — ре2 — 0 ,

(15)

где N - большой параметр, ам = 2р/тЬдр0, Ь - характерный размер рассматриваемой области, д - сила тяжести, р есть вязкость жидкости.

Приведенные итерации повторяются до тех пор, пока во всей области не будет полностью вычислено поле скоростей на текущем временном шаге. Затем на следующем временном шаге полагается ^+1 = Нг+1. Если | ^+1 | > 0 вплоть до указанного порядка, то необходимо увеличить N и повторить цикл. При N ^ то выполняется переход от слабосжимаемой жидкости к несжимаемой. Подставив полученные значения скорости в уравнение переноса, находим плотность жидкости в следующий момент времени.

После расчёта поля скоростей мы переходим к решению уравнения переноса с соответствующими значениями компонент вектора скорости на текущем временном шаге.

При использовании в численном моделировании прямоугольной системы координат исследуемая область покрывается равномерной по х и у сеткой ячеек

^ = /х^+1 = (г + 1)к, к> 0 , г = 0,1,...,Ь; (Ь + 1)к = Хт£

7 \У?+1 = 3 + 1)^, к> 0 , 3 = 0, 1,...,М; (М + 1)к = Ут;

где к - размер шагов сетки, Ь и М - числа ячеек сетки соответственно в направлениях осей х и у.

Здесь используется последовательное вычисление значений скоростей и плотностей жидкостей на отдельных накладывающихся сетках. Это означает, что каждая компонента вектора скорости вычисляется на разнесённой сетке отдельно, после чего вычисляется изменение плотности на следующем временном шаге в соответствие с изменением значения скорости.

В случае декартовой системы координат и равномерной сетки двумерная разностная схема имеет следующий вид:

_ аиг^+1 + 1 + Щ+1,^+1 — 1,г-1^+1 — Уг+1^-1 + Щ-1,^ — 1 + &Мг+1^ + ЪЩ-\^

2Ь~+2а ’ [Ь)

_ (Щ+1,;/ + + Ьгц^-\

Щ'3 2Ь~+2а ' 10

Так как начальные значения плотностей для жидкостей известны, то для уравнения переноса была выбрана явная схема:

1

2 р

Р%4 — 7Г \ти(Рг+1^ ~ Рг-I,]) + Tv(piJ+1 — Рг^-1)] , (18)

где а = 4кам, Ь = 4к(ам — 1), и и V - узлы ячеек с координатами (г,3), т - шаг по времени.

При замене дифференциальной задачи конечно-разностным представлением особое внимание следует уделять аппроксимации граничных условий, так как конкретная аппроксимация последних влияет на точность метода, устойчивость схемы, а также на скорость сходимости.

В случае, когда боковые стенки - твёрдая поверхность, то условие непротекания представляется в виде

V*, ]—1 = 0 , (19)

а условие прилипания - в виде

иг±1,7-1 = 0 . (20)

Из последнего условия следует, что в случае прямоугольной декартовой системы коор-

динат скорость жидкости у твёрдой границы будет

(п) (п)

ИЙ.1] = + Л3 • (2!)

Учитывая поведение жидкости на границе с твёрдым скелетом и, решая систему

уравнений (2.1)—(2.4), получим картину изменения положения двух жидкостей относительно друг друга в отдельном капилляре.

Для получения приближенного решения задачи Маскета необходимо рассмотреть периодическую структуру, состоящую из изолированных капилляров и тем самым усреднить задачу при е \ 0.

Литература

1. Antontsev S., Meirmanov A., Yurinsky B.V., A Free Boundary Problem for Stokes Equations: Classical Solutions jj Interfaces and Free Boundaries. - 2000. - 2. - P.413-424.

2. Биркгоф Г. Неустойчивость Гельмгольца и Тейлора j В кн.: Гидродинамическая неустойчивость. - М.: Мир, 1964. - C.68-94.

3. Inogamov N.A. Turbulent phase of the Rayleigh-Taylor instability j j N.A. Inogamov. - 1978. -Черноголовка: Институт теоретич. физ. им. Л.Д.Ландау АН СССР, 1978.

4. Lord Rayleigh. Theory of sound V.2 j Lord Rayleigh. - N.Y.:Dover Publications Inc., 1984.

Б. Taylor G. The instability of liquid surface when accelerated in a direction perpendicular to their plans. I jj Proc. Roy. Soc. - 19Б0. - Ser. A. - 201,1065. - P.192-196.

6. Lewis D.G. The instability of liquid surface when accelerated in a direction perpendicular to their plans. II. jj Proc. Roy. Soc. - 19Б0. - Ser.A. - 202,1068. - P.81-96.

7. Birkhoff G. Los Alamos Scientific Lab., Rept. №LA-1982. - Los Alamos, 19ББ.

8. Беленьких С.З., Фрадкин Е.С. Теория турбулентного перемешивания jj Тр. ФИАН СССР. - 196Б. - 29. - C.207-238.

9. Fahuai Yi, Global classical solution of Muskat free boundary problem // J. Math. Anal. Appl. - 2003. - 288. - P.442-461.

10. Radkevich E., On the spectrum of the pencil in the Verigin-Muskat problem // Sbornik: Mathematics. - 199Б. - 80,1. - P.33-74.

11. Siegel M., Caflish R.E, Howison S. Global existence, singular solutions, and ill-posedness for the Muskat problem jj Comm. on Pure and Appl. Math. - 2004. - LVII. - P.1-38.

12. Владимирова Н.Н., Кузнецова Б.Г., Яненко Н.Н. Численные рассчеты симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости //В кн.: Некоторые вопросы вычислительной и прикладной математики. - Новосибирск: Наука, 1996. -C.186-192.

13. Daly B.J., Numerical study of two fluid Raylaigh-Taylor instability // Phys. Fluids. - 1967. -10,2. - P.297-307.

ABOUT THE NUMERICAL SIMULATION OF MUSKAT PROBLEM О.V.Galtsev, O.A.Galtseva

Belgorod State University,

Pobedy St., 85, Belgorod, 308015, Russia, e-mail: [email protected]

Abstract. The filtration of two immiscible incompressible liquids with various viscosity and density divided by free boundary is studied. The Muskat model being most plausible and physically correct is considered. The problem connected with the filtration modelling in the porous space structure is investigated.

Key words: Muskat problem, free boundary problem, liquid filtration.