CC BY
91
УДК 517.958:531.72, 517.958:539.3(4)
НЕУСТОЙЧИВОСТЬ РЭЛЕЯ-ТЕЙЛОРА В ЗАДАЧЕ МАСКЕТА СО СВОБОДНОЙ ГРАНИЦЕЙ 8) О.В.Гальцев
Белгородский государственный национальный исследовательский университет, ул. Студенческая, 14, Белгород, 308007, Россия, e-mail: [email protected]
Аннотация. Настоящая работа посвящена моделированию совместного движения двух несмешивающихся несжимаемых жидкостей в пористых средах. Жидкости имеют различную плотность и изначально разделены свободной границей. Рассматриваются результаты численного моделирования задач со свободной границей на микроскопическом уровне для абсолютно твердого и упругого скелета различной геометрии. Эти задачи имеют естественный малый параметр е, который представляет собой отношение среднего размера поры к размеру рассматриваемой области. Если уменьшать параметр е, то решения микроскопических уравнений моделируют усредненную картину движения жидкостей. Для движения жидкостей в упругом грунте присутствует поверхность контактного разрыва, в то время, как при движении жидкостей в абсолютно твердом скелете вместо свободной границы имеется переходная фаза, занятая смесью двух жидкостей.
Ключевые слова: задача Маскета, задача со свободной границей, фильтрация жидкости.
1. Проблема Маскета
Неустойчивость Рэлея-Тейлора (см. [1], [2]) представляет собой неустойчивость границы раздела двух жидкостей различной плотности, которая возникает, когда более легкая жидкость толкает более тяжелую. Эквивалентная ситуация возникает, когда сила тяжести действует на две жидкости разной плотности - более плотная жидкость находится над жидкостью меньшей плотности. В настоящей работе мы будем рассматривать рэлей-тейлоровскую неустойчивость неоднородной жидкости, заполняющей поры в твердом скелете. В дополнение к несомненному теоретическому интересу, эта задача имеет важное практическое значение. Например, описание вытеснения одной жидкости другой в пористой среде. Решение этой проблемы связано с выбором правильной математической модели.
Существует много различных типов математических моделей, но мы будем иметь дело только с классом моделей, которые мы называем физически корректными. Феноменологическую модель будем считать физически корректной, если она является одной из базовых моделей механики сплошных сред (как, например, уравнение Стокса для медленного движения вязкой жидкости, или уравнение Ламэ для перемещений упругого скелета), или асимптотически близкой к какой-либо физически корректной феноменологической модели.
8 Работа выполнена в рамках ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы (госконтракт №02.740.11.0613).
Рис. 1.1- область П+(£), 2 - область О (£), 3 - свободная граница Г(£).
Среди математических моделей, описывающих совместное движение двух несмеши-вающихся жидкостей, наиболее верной (физически корректной) является задача Мас-кета, предложенная М. Маскетом [3]. Эта модель описывает процесс фильтрации двух несмешивающихся несжимаемых жидкостей различной вязкости и плотности, разделенных некоторой подвижной (свободной) границей. Движение первой жидкости под воздействием силы тяжести в области П+(£) с постоянной вязкостью Ц и постоянной плотностью р+ описывается системой уравнений Дарси
V
+
к.
+ р}
+
е,
V ■ v+ = 0, х € П+(£),
для микроскопической скорости v+ и давления р+ жидкости. Единичный вектор е совпадает с направлением силы тяжести.
Соответственно, движение второй жидкости под воздействием силы тяжести в области П_(£) с постоянной вязкостью I и постоянной плотностью р- описывается системой уравнений Дарси
к
V- = Vр- + р-е, V ■ V- = 0, х € П-(£), (2)
для микроскопической скорости V- и давления р-. На общей свободной границе
г(*) = ш+(£)р| д п-(г)
давления жидкостей и их нормальные скорости непрерывны:
р+ = рп х € г(t),
V
+
п
V
п
Уп, х € Г^),
(4)
где п есть единичный вектор нормали к границе Г^) в точке х € Г^) и Уп - скорость в нормальном направлении к границе Г^) в точке х € Г^).
Условие (4) означает, что граница Г^) есть материальная поверхность. Иными словами, во время движения, она состоит из постоянного количества материальных точек.
Этот факт позволяет сформулировать слабую постановку задачи Маскета. Определим давление р/ неоднородной жидкости
п+ если
плотность р/ и скорость V
р/ = р+ если х Е П+(Ь), р/ = р- если х Е П (Ь),
р/ = р+ если х Е П+(Ь), р/ = р- если х Е П-(Ь),
V = v+ если х Е П+(Ь), V = V- если х Е П-(Ь).
Таким образом неизвестные функции V, р/, и р/ будут решениями системы уравнений фильтрации Дарси
к
V = —Vр/ + р/ е, V- V = 0, х Е П, Ь > 0 , (5)
ц
и уравнения переноса
йр/ _ др/ йЬ дЬ
+ V -V р/ = 0, х Е П, Ь > 0 , (6)
Первое уравнение в (5) (закон Дарси) следует понимать в обычном смысле почти всюду в Пт = П х (0,Т), а второе уравнение (уравнение неразрывности) понимается в смысле теории распределения. Уравнение переноса понимается также в смысле теории распределения, если использовать равенства
V - Vр = V - ^р), V -Vр/ = V - ^р/).
Проблема дополняется однородным граничным условием
V - п = 0, х Е Б = дП, Ь> 0 , (7)
где п есть нормальный вектор к границе Б, и начальным условием
р/(х, 0) = р0(х), х Е П , (8)
с разрывными начальными значениями:
р0 (х) = р+, если х Е П+ и р0 (х) = р-, если х Е П-.
Считается, что задачу Маскета можно рассматривать в двух случаях. В обоих случая проблему легко сформулировать, но почти невозможно решить. По этой причине известно очень мало работ о ее классическом или слабом решении. Существует несколько работ о классической разрешимости локальной или глобальной во времени задачи, но до сих пор нет каких-либо результатов по слабой разрешимости (см. [4], [5], [6]).
Следуя схеме, предложенной Р. Барриджем и Дж. Келлером [7], попытаемся найти более общую физически корректную математическую модель, описывающую тот же
процесс. Чтобы объяснить идею, во-первых, мы рассмотрим систему фильтрации Дарси, которая отвечает за динамику в задаче Маскета. Хорошо известно, что эта система есть асимптотический предел уравнения Стокса для несжимаемой вязкой жидкости, когда размер пор стремится к нулю (см. [7], [8]). Р.Барридж и Дж.Келлер предложили рассматривать более общую систему
д2ю ( ( д ю) )
ат р£~д№ = ^ (х^0! Х’дг) +(1 - х£)ах°(х’ ю) - Р Ч + р£ е , (9)
р + арЧ - ю = 0 , (10)
для перемещения ю и давления р. Микроскопические системы (9), (10) описывают совместное движение вязкой жидкости в поровом пространстве и упругого скелета и понимаются в смысле теории распределения. Эта система содержит уравнение Стокса для вязкой жидкости в поровом пространстве П/, уравнение Ламэ для твердого скелета в П и граничное условие (непрерывность нормальных напряжений) на общей границе Г = дП/ П дП3. В (9) 0(х, ю) есть симметричная часть Vю, х£ - характеристическая функция порового порового пространства П/, е = 1/Ь - безразмерный размер пор,
Ь 2р 2А
ат 2 , а№ т , а^
дт2 тЬдро Ьдр0
р£ = р/Х£ + р8(1 — Х£), ар = С/Х>£ + С2(1 — х£) ,
I - средний размер пор, Ь - характерный размер рассматриваемой области, т есть характерное время процесса, р/ и р3 - безразмерные плотности жидкости в порах и твердого скелета, отнесенные к средней плотности воды р0, с/ и с3 есть, соответственно, безразмерные скорости звука в жидкости и твердом скелете, д есть значение ускорения силы тяжести, р - вязкость жидкости, и А - постоянная упругости Ламэ. В дальнейшем будем предполагать, что структуры областей П/ и П периодические с периодом е. То есть Х£(х) = х(х/е) с 1 - периодической функцией х(у).
Теоретически система (9), (10) с соответствующими граничными и начальными условиями наиболее точно описывает данный физический процесс, но они все еще не имеют никакого практического значения. Оно появляется только после усреднения. Все предельные режимы системы (9), (10) зависят от безразмерного параметра ат, а№, а\ [9]. Например, для филтрации ат ~ 0, для абсолютно твердого скелета а\ ~ то, и для несжимаемой среды ар ~ то. Таким образом, мы можем рассматривать эти параметры, как некоторые переменные значения в процессе перехода к пределу, и вместо стандартного предельного режима (когда все параметры фиксированы) получим несколько осредненных моделей, асимптотически близких к основной. Все эти модели описывают тот же физический процесс, но с разной мерой приближения.
Для фильтрации несжимаемой жидкости в несжимаемом твердом скелете стандартная усредненная модель для фиксированных а№ = р0, а\ = Ао, и ат = т0 имеет вид
,8 2ю ^
То = V- Р + £ е , £ = тр/ + (1 - т)р3,
дь2
^ ( дю\ Ґ
Р = р0Ш1 : О ( х, + Ао^2 : О(х, ю) + J Ж.3(і — т) : О(х, ю(х, т))йт .
Система (9), (10) на микроскопическом уровне имеет другие асимптотические пределы, чем система фильтрации Дарси (ат = 0, ар = то, а\ = то, а^ = р1 є2):
д ю 1 * . „ . „ д ю
Ж = •(—Ур + р<е)' V• аТ = 0'
или система Био-Терцаги для пороупругости (ат = 0, ар = то, а\ = Л0, а^ = р1є2)
(см. [9], [10], [11]):
дю ди 1 ^ ,
~Ш = т Ж + 7і ’ •(—Ур + р‘е) ■
дю ди
+(1 — т) а) = 0 ■
V • (Л0^0 : О(х, и)) — Vр + д е = 0 , или система вязкоупругой фильтрации
V • Р — Vp + £ е = 0 , V • ю = 0 ,
Р = р0Ш4 : О ^х, ~дЮ^ + А0^5 : 0(х, ю) + J Ж.е(Ь — т) : 0(х, ю(х, т))йт .
Эту схему мы можем применить к задаче со свободной границей, которая на микроскопическом уровне для фильтрации в несжимаемой среде будет иметь вид
V- ^х£ам0^ х,8т^ +(1 — х£)ал0(х, ю£) — р£+ р£ е = 0 , (11)
V • ю£ = 0 , (12)
К
Ї ■ Ж + = „. „
£
Задача наиболее полно изучена в [12]. Основные результаты говорят о том, что для любой е > 0 существует слабое решение {ю£, р£, р£} задачи со свободной границей (11)
- (13) и для а№ = р0, ал = А0 соответствующие последовательности сходятся при е \ 0 к решению {ю, р, р} (усредненной) задачи Маскета для вязкоупругой фильтрации
V- Р — Vp + р е = 0 , V- ю = 0 , = 0 . (14)
аЬ
При доказательстве этого результата было использовано понятие двух - масштабной сходимости [13].
м її
*=0.003
4=0.255
*=0.525
*=1.035
Рис. 2. Поверхность контактного разрыва в различные моменты времени.
С другой стороны формальный предел є \ 0 в (11) - (13) для а\ ^ то и а^ = р1є2 ведет к задаче Маскета (5), (6). Понятно, что тот же формальный предел можно получить, если в качестве базовой модели на микроскопическом уровне будем рассматривать задачу со свободной границей для системы Стокса
д?п£
V • —) — рЧ) + р£
е
0
х Є Щ , і > 0,
V • ю£ = 0 , ^ = 0
Iи
х
Є Щ*, і > 0 , ю£(х) = 0 , х Є дЩ£* , і > 0,
:15)
16)
только в поровом пространстве П/. Последняя проблема была изучена в [14], где авторы доказали, что для любой е > 0 задача (15) имеет единственное классическое решение {ю£, р£, р£}.
Целью настоящей работы является компьютерное моделирование задачи (11) - (13) с а№ = р0, ал = А0, и задачи (15), (16) для двух различных структур порового пространства: несвязные капилляры и несвязный твердый скелет в единичном квадрате К2.
Численное решение этой задачи в одном капилляре в абсолютно твердом скелете показало совпадение с результатами [15]. На рис. 2 можно увидеть гладкую свободную границу в капилляре в различные моменты времени.
Мы показали, что движение жидкостей в (15), (16) зависит от трех факторов: отношения 8 = р+ /р- плотностей р+ и р- верхней и нижней жидкостей, вязкости р жидкостей и размера пор е. Изменяя эти параметры, можно получить различные сценарии развития неустойчивости Рэлея-Тейлора. Для задачи (11) - (13) процесс зависит от тех же параметров 8, р, е, как и ранее, и дополнительно параметра коэффициента упругости Ламэ А.
*=50 *=860 *=2631 *=3012 *=4873
Рис. 3. Несвязные капилляры: численное моделирование для абсолютно твердого скелета (сверху) и для упругого скелета (снизу).
Процесс численного усреднения (є \ 0) моделируется увеличением количества капилляров для первой геометрии, или количества элементов твердого скелета для второй геометрии. Таким образом, можно считать, что при достаточно малом є система (15), (16) описывает классическую задачу Маскета, а система (11) - (13) описывает задачу Маскета вязкоупругой фильтрации.
На рисунке рис. 3 показаны результаты численного моделирования для первой геометрии задачи (15), (16) движения в абсолютно твердом скелете (сверху), и задачи (11)
- (13) движения в упругом скелете (снизу).
На рис. 4 можно видеть случай упругого скелета в большем масштабе, чтобы лучше рассмотреть изменение свободной поверхности раздела фаз.
Кроме того, проводились численные расчеты для различных значений Л и 8:
• для є = 2 * 10“5, Л \ 0 и 8 = 1.25 наблюдается изменение поверхности раздела жидкостей;
• для є = 2 * 10“5, Л = 0.5 и 8 \ то наблюдается изменение поверхности раздела жидкостей;
• для є = 2 * 10“5, Л = 0.5 и 8 \ 1 изменения поверхности раздела жидкостей не наблюдается.
Такой же вывод справедлив для случая второй геометрии (несвязный твердый скелет). На рис. 5 показаны сравнительные результаты для аналогичных значений 8, р, д, р3, Л, є (см. Таблица 1) и тех же начальных значений.
*=50
*=860
*=2631 *=3012
Рис. 4. Случай упругого скелета в большем масштабе.
*=50 *=860 *=2631 *=3012 *=4873
Рис. 5. Несвязный скелет: численное моделирование для абсолютно твердого скелета (сверху) и для упругого скелета (снизу).
Таблица 1
Значения параметров для случаев абсолютно твердого скелета и упругого скелета
Параметр Абсолюно твердый скелет Упругий скелет
Р+ 998.2 998.2
. 1 — 800 800
Ps — 2000
f!+ 10-2 10-2
Ц~ 9 * 10-1 9 * 10-1
A — 0.5
£ 2 * 10-5 2 * 1 О 1 Сл
L 100 100
Легко заметить, что для одного и того же времени процесса в задаче Маскета возникает переходная фаза (см. [16]), в то время как при вязкоупругой фильтрации, свободная граница сохраняется.
Как и в случае геометрии первого типа, были проведены численные расчеты для различных значений A and 8:
• для 8 = 1.25, £ = 2 * 10-5 и A \ 0, наблюдается изменение поверхности раздела жидкостей;
• для 8 \ то, £ = 2 * 10-5 и A = 0.5, наблюдается изменение поверхности раздела жидкостей и процесс фильтрации протекает быстрее с увеличением 8;
• для 8 =1.25, £ = 2 * 10-5 и A \ то, наблюдается изменение поверхности раздела жидкостей и процесс фильтрации протекает быстрее с увеличением A.
На основе численного моделирования можно сделать вывод, что:
1) при движения жидкости в абсолютно твердом скелете вместо свободной границы возникает переходная фаза, в то время, как при движении жидкостей в упругом скелете сохраняется свободная граница;
2) за одно и то же время развития процесса, жидкости в абсолютно твердом скелете полностью поменялись местами, в то время, как жидкости в упругом скелете все еще сохраняют свое положение.
Таким образом, решение модели вязкоупругой фильтрации - классическое и обладает гладкой и устойчивой свободной границей, в то время, как решение задачи Маскета в лучшем случае будет только обобщенным с переходной фазой вместо свободной границей.
2. Постановка задачи
2.1. Абсолютно твердый скелет. Предположим, что П единичный квадрат (0 < Х\ < 1) х (0 < х2 < 1), П = Щ и Б£ и П*, где П поровое пространство (область, занятая жидкостью), П£ - твердый скелет и Б£ есть поверхность «твердый скелет -поровое пространство».
Для геометрии задачи первого типа поровое пространство представляет собой набор изолированных капилляров
п— 1
є
Qf = U {єк < x1 < є(к + 1)m) x (0 < x2 < 1).
k=0
Здесь е = 1/п, т - пористость, 0 < т < 1.
Для геометрии второго типа твердый скелет есть объединение несвязных квадратов
п1
Щє = U (є(к+1)в < x1 < є(к+1)(1—в)) x (є(к+1)в < x2 < є(к+1)(1—в)) , 2в = 1 —л/m
k=0
Как мы уже отмечали выше, две несмешивающиеся жидкости моделируются неоднородной жидкостью, где плотность р может принимать только два значения р+ или р-. Скорость v, давление р, и плотность pf неоднородной жидкости в поровом пространстве Щ описываются системой уравнений Стокса
V • (a^D(x, v) — pf I) — pf e2 = 0 , V - v = 0 , (17)
и уравнением переноса
df + v •V p, = 0. (18)
Дифференциальные уравнения (2.1) и (2.2) дополняются условием нормализации
/ р,(x, t)dx = 0 , (19)
Jnf
и однородным граничным условием
v = 0 , x e S U S£, t> 0 , (20)
где S = дЩ, и начальным условием
Pf (x, 0) = p0(x), x e Щ , (21)
pf (x) = p+ = const > 0, если x e = Щ£ Л ,
и
р^(x) = p— = const > 0, если , x Є Щ— = Qf П Щ .
2.2. Упругий скелет. Совместное движение вязкой жидкости в поровом пространстве и упругого скелета на микроскопическом уровне описывается системой (10) - (13), которая состоит из стационарного уравнения Стокса
V- ^-^Р/ - Р/ 62 = 0 , V- Wf = 0 , (22)
для перемещения W/ и давления р/ в неоднородной жидкости, стационарного уравнения Ламэ
V - (a\B(ws)) - Vps - Рз 62 = 0 , V - Ws = 0 , (23)
для перемещения ws и давления рз в упругом скелете, и двух граничных условий на
общей границе «твердый скелет - поровое пространство» Б£:
w/ = Ws, (24)
(а^О() - р/1) - п = (ахВ^з) - Рз1) - п , (25)
где п есть единичный вектор нормали к границе Б£.
Система (14) - (17) дополняется начальным и граничным условиями
W/(х, 0) = 0 , х Е П£ , (26)
w = 0 , х Е Б, г > 0 , (27)
где
w(x,t) = W/(х,г) , если х Е П/ ,г> 0 ,
и
w(x, г) = ws(x, г), если х е П, г > 0, и условием нормализации
/ р(х,Ь)йх = 0 , (28)
./П£
где
р(х,г) = р/ (х,г), если х е П/ ,г> 0,
и
р(х,г) = ps(x,t), если х е Пs, г > 0,
Таким образом имеем задачу Коши для уравнения переноса (10). На микроскопическом уровне уравнение переноса для плотности жидкости определяется только в поро-вом пространстве П/. Для существующей микроскопической модели (см. [9]) характеристическая функция х£ порового пространства есть неизвестная функция и определяется как решение задачи Коши
1к + Ж- V^£ = 0 ’ Х£(х’ 0) = х£(х) •
Это означает, что поверхность между твердым скелетом и жидкостью Бє в этой модели есть материальная поверхность и нам не нужно граничное условие для плотности жидкости на Бє. Но наша основная динамическая система здесь линейная, где хє = Xе известная функция. Таким образом поверхность Бє не больше материальной поверхности и нам необходимо граничное условие для плотности жидкости на части Бє, где жидкость «проникает» в поровое пространство. Чтобы избежать этого, распространим задачу Коши (10) на всю область П. Во первых, предположим, что функция (х) определена во всей области П и
р®(х) = р+ в П+, р0(х) = р- в П-.
Перепишем задачу Коши для плотности жидкости как задачу Коши для плотности смешивания
Р = ХєРї + (1 - Хє)Рз
в виде
дР + '^Р = 0 ’ (х>*) Є Пт ; р(х, 0) = р0(х) , х Є П , (29)
где
Р(0)(х) = ХєР°(х) + (1 - хє)р
-/в •
3. Вычислительный алгоритм
3.1. Абсолютно твердый скелет. В качестве численного метода решения системы уравнений (9) - (10), дополненной условиями (11) - (13), был выбран метод крупных частиц [17]. Этот метод основан на расщеплении исходных дифференциальных уравнений по физическим процессам. Метод крупных частиц - это усовершенствованный метод «частицы в ячейке» разработанный Харлоу. Этот метод используется для решения определенного класса уравнений в частных производных, в том числе неустойчивости Рэлея-Тейлора для задачи Маскета со свободной границей.
Процесс решения системы (9) - (10) очень трудный из-за наличия стратификации. Неоднородность жидкости требует дополнительного расчета поля плотностей. Чтобы найти неизвестные функции, процесс вычисления можно представить в виде трехэтапной схемы.
Этап 1. Вычисляется промежуточное значение скорости V из уравнений
V - (амО(£)) - ре = 0 , (30)
и промежуточное значение плотности р из уравнения
+ V - р = 0 , (31)
Этап 2. Основная трудность при численном решении уравнений (2.1), (2.2) связана с вычислением поля давления. Успех в преодолении этих трудностей был достигнут благодаря идее искусственной сжимаемости [18]. Этот метод очень важен для вычисления неизвестного давления р из уравнения
др _ дг + ^
V
0 •
(32)
Здесь необходимо соблюдать условие соленоидальности (V - V = 0), что достижимо при увеличении ср каждый шаг по времени, пока равенство не будет выполняться с заданной точностью.
Этап 3. Вычисляется итоговое значение скорости V из уравнения
V - (амО^)) - Vp - ре = 0
и давления Р из уравнения
др
дг + v - р
0,
(33)
(34)
Наконец, выпишем основные разностные схемы. Область интегрирования покрывается стационарной (Эйлеровой) разностной сеткой произвольной формы (прямоугольная сетка в нашем случае, см. рис.6):
П
/
х1г+1/2) = гк1, Л1 > 0; г
0,1,-
х2?+1/2) = уЬ2, Л2 > 0; ] = 0,1, ••
N1;
,N2;
;;
где Л1, Л2 - размер сетки, N1, N - количество ячеек сетки, соответственно, в направлении х1 и х2 (точка с координатами (%,]) совпадает с центром ячейки). Здесь, как и в исходном методе расщепления, будем использовать «шахматную» сетку. Это дает возможность четко интерпретировать каждую ячейку, как элемент объема, который характеризуется рассчитываемыми давлением и плотностью в его центре.
Рис. 6. Шаблон расчетной сетки.
Изменяются только величины, относящиеся к ячейке в целом. По этой причине пренебрегаем конвективными члены, относящимися к поступательным эффектам. В остальных уравнениях р выносится за знак дифференцирования, и уравнения (2.1)—(2.2) решаются для и, V по времени.
Простейшая конечно-разностная аппроксимация уравнения (3.1) будет иметь вид
~.п _ 2 ~п I ~.п ~.п _ 2 ~.п I ~п
и 1+3/2,] 2иг+1/2,] + аг-/2,] _ иг+1/2,]+1 2 ui+1/2,j + иг+1/2,]-1
а н? =-а1 ’
п <,п _ О п чп | /1 чп п чп _ О /1 <,п | /1 чп
Vi+1,j+1/2 2 Vi,j+1/2 + Vi-1,j+1/2 Vi,j+3/2 2 Vi,j+1/2 + Vi,j—1/2
----------------- = -а11-----------------Щ-----------------р^+1/2 е2
где п - шаг по времени, ~, V - промежуточные значения скорости течения.
Величины с дробными индексами относятся к границам ячеек, например:
Vi,j + Vi+1,j Ч+1/2,] = -----2------ ’
В методе крупных частиц необходимо вычислить плотность потока (3.2) через границы ячейки. Это говорит о том, что плотность большой частицы находится в движении только за счет нормального компонента скорости к границе. Эти значения параметров течения на следующем временном слое вычисляются в соответствие с уравнениями (направление потока слева направо и снизу вверх):
ру1 - = <р~‘)?++11/2,] - <р~‘)Г-+11/2,] (р~)п+|11/2 - {р~)п+-1/2 , .
Р = н Н2 • (35)
В случаях, когда необходимо определить значения функции в точках сетки (см. рис. 6), достаточно воспользоваться средним арифметическим, например:
рi+1/2,j = 2 (рг+1,^ + р^) •
На завершающем шаге (при использовании дискретной модели среды) проводится дополнительный расчет плотности, что сглаживает колебания и повышает точность вычислений. Комбинируя различные этапы, получаем ряд разностных схем, которые и составляют метод крупных частиц. Выбранный метод может быть интерпретирован с различных точек зрения: метод расщепления, смешанный метод Эйлера-Лагранжа, вычисление в локальных лагранжевых координатах с масштабированием на предыдущей сетке, различные обозначения закона сохранения для жидкого элемента (большой частицы), эйлерова разностная схема.
При замене дифференциальной задачи на конечно-разностное представление необходимо обратить внимание на приближение граничных условий, так как их конкретное приближение прямо влияет на корректность метода и устойчивость схемы, а также скорость сходимости. Сформулируем граничные условия, введя серию фиктивных ячеек (так, чтобы каждая вычисляемая точка становилась внутренней и описанный алгоритм сохранялся для всех ячеек). Для первого порядка аппроксимации достаточно одного слоя, для второго порядка аппроксимации - два слоя. Так в случае, когда боковые стенки есть твердая поверхность, то условие непротекания представляется в виде
и-1/2^ = 0
(36)
а условие прилипания в виде
Vi-1/2,j+1/2 = 0 • (37)
В случае задачи на плоскости, геометрические характеристики дробных ячеек могут быть определены непосредственными измерениями. В случае аксиально-симметричной задачи необходимо провести дополнительные вычисления для нахождения расстояния от каждой дробной ячейки до оси симметрии. Разностные уравнения для дробных ячеек получаются путем небольшой модификации разностного уравнения для целой ячейки.
3.2. Упругий скелет. Компьютерное моделирование совместного движения вязкой жидкости в поровом пространстве и упругого скелета (на микроскопическом уровне) есть численное решение системы, которая состоит из стационарной системы уравнений Стокса (14) для перемещения Wf и давления гр$ неоднородной жидкости и стационарного уравнения Ламэ (15) для перемещения ws и давления р3 упругого скелета с соответствующими условиями на общей границе.
В качестве численного метода для моделирования неустойчивости Рэлея-Тейлора с учетом упругой составляющей скелета был выбран аналогичный метод, описанный выше (метод крупных частиц). Рассматриваемая область заменяется системой жидких частиц, которые совпадают с ячейкой эйлеровой сетки (см. рис. 6). Предлагаемый алгоритм состоит из четырех главных этапов:
Этап 1. Решается система уравнений Ламэ (15) с заданными граничными и начальными условиями для перемещения ws и р3.
ws = Wf, х е Б£, ws = 0, х е Б^ (38)
Решение этой системы ничем не отличается от решения системы Стокса (п. 3.1), то есть вводится искусственная сжимаемость и вычисляется неизвестное ps из уравнения
+ CpsV • Ws = 0 • (39)
dps dt
Ввиду того, что жидкость несжимаема, необходимо также увеличивать cp,s в каждый момент времени до тех пор, пока не будет выполнено условие соленоидальности (V • Ws = 0).
Следует отметить, что процедура составления двумерной разностной схемы уравнения Ламэ идентична, описанной ранее для уравнения Стокса. Простейшая конечноразностная аппроксимация уравнения (2.7) будет иметь следующий вид:
wli+3/2,j — 2 wli+1 /2,j + w1i—/2,j wn+1 /2,j+1 — 2 wli+ 1 /2,j + Wli+ l/2,j-1
ax h = -«* h ’
i1 2 ян1 —I— ян1
2 i+1,j+1/2 2 w2 i,j+1/2 + w2 i— 1,j+1/2 _
h =
w2i,j+3/2 — 2 w2i,j+1/2 + w2i,j — 1/2
a\-
= -а\ Н2 Рsi,j+1/2 е2
где п - шаг по времени, т1, и>2 значения перемещения, рs плотность скелета.
Этап 2. Используя найденные значения ws и ps, находим нормальное напряжение на границе Б£
где п единичный вектор нормали к границе Б£.
Этап 3. Далее решается система уравнений (14) в П£, как и для случая абсолютно твердого скелета, повторяя этапы 1 - 3 из предыдущего пункта (3.1), с условием на общей границе 5£
в котором, на основе результатов расчета на предыдущем этапе, нам известна правая часть равенства.
Это граничное условие предполагает, что вектор смещения и давления удовлетворяет непрерывности нормальных напряжений на общей границе между жидкостью и упругим скелетом.
Найденное значение скорости жидкости подставляем в уравнение переноса (29) в ^, что позволяет вычислить значение плотности pf на следующем шаге.
Этап 4. Зная скорость дWf /дЬ в жидкой части, находим ws на следующем шаге по времени из условия непрерывности нормальных напряжений (16).
Таким образом, учитывая поведение жидкостей на границе «жидкость-упругий скелет» и, решая систему уравнений (14), (15), (29) с соответствующими начальными и граничными условиями (16) - (28), получаем численную аппроксимацию совместного движения жидкости и упругого скелета.
В настоящей работе показано, как необходимо производить численное моделирование процессов совместного движения несмешивающихся жидкостей. Мы начали с задачи движения двух жидкостей со свободной границей в рамках микроскопической модели. Теоретически, эта математическая модель наиболее подходящая для описания рассматриваемого физического процесса. Однако, она не имеет никакого практического значения, так как приходится решать задачу в физической области в несколько сотен метров, в то время, как коэффициенты колеблются от физического размера в несколько микрон. Практическая значимость модели появляется только после усреднения. В свою очередь, усреднение имеет по крайней мере два уровня аппроксимации, которые зависят от безразмерного критерия физической задачи. Первый уровень приближения - хорошо известная задача Маскета. Второй уровень приближения - задача со свободной границей на микроскопическом уровне (задача Маскета для вязкоупругой фильтрации). На основе численных расчетов для периодической структуры было смоделировано усреднение за счет увеличения числа элементарных ячеек в единице
(40)
4. Заключение
объема. Результаты показывают, что решение проблемы Маскета неустойчиво, свободная граница резко меняет структуру и превращается в некоторую области (переходную фазу), занятую смесью из двух жидкостей. Между тем, в случае вязкоупругой фильтрации решение Маскета остается классическим, то есть обладает гладкой свободной границей. Поэтому задача Маскета для вязкоупругой фильтрации является естественным обобщением классической задачи Маскета, которая все еще остается нерешенной, как математическая задача, и очень трудно реализуема численно.
Литература
1. Rayleigh, Lord (John William Strutt). Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density // Proceedings of the London Mathematical Society. - 1883. - 14. - P.170-177.
2. Taylor, Sir Geoffrey Ingr.m The instability of liquid surfaces when accelerated in a direction perpendicular to their planes // Proceedings of the Royal Society of London. - Series A, Mathematical and Physical Sciences. - 1950. - 201(1065). - P.192-196.
3. Muskat M. Two-fluid system in porous media. The encroachment of water into an oil sand // Physics. - 1934. - 5. - P.250-264.
4. Fahuai Yu. Global classical solution of Muskat free boundary problem // J. Math. Anal. Appl. - 2003. - 288. - P.442-461.
5. Radkevich E. On the spectrum of the pencil in the Verigin-Muskat problem // Sbornik: Mathematics. - 1995. - 80;1. - P.33-74.
6. Siegel M., Caflish R.E., Howison S. Global existence, singular solutions, and ill-posedness for the Muskat problem // Comm. on Pure and Appl. Math. - 2004. - LVII;1. - P.38.
7. Burridge R., Keller J.B. Poroelasticity equations derived from microstructure // Journal of Acoustic Society of America. - 1981. - 70;4. - P.1140-1146.
8. Sanchez-Palencia E. Non-Homogeneous Media and Vibration Theory // Lecture Notes in Physics: Springer-Verlag,1980. - №129.
9. Meirmanov A. Nguetseng’s two-scale convergence method for filtration and seismic acoustic problems in elastic porous media // Siberian Mathematical Journal. - 2007. - 48. - P.519-538.
10. Terzaghi K. Die Berechnung der Durchlassigkeitsziffer des Tones aus dem Verlauf der Hyd-rodynamischen Spannungsercheinungen // Sitzung berichte. Akademie der Wissenschaften, Wien Mathematiesch-Naturwissenschaftliche Klasse. - 1923. - 132;IIa. - P.104-124.
11. Biot M. General theory of three dimensional consolidation // Journal of Applied Physics. -1941. - 12. - P.155-164.
12. Meirmanov A. The Muskat problem for a viscoelastic filtration // Accepted for publication to Interfaces and Free Boundaries.
13. Nguetseng G. A general convergence result for a functional related to the theory of homogenization // SIAM J. Math. Anal. - 1989. - 20. - P.608-623.
14. Antontsev S., Meirmanov A., Yurinsky B. A Free Boundary Problem for Stokes Equations: Classical Solution // Interfaces and Free Boundaries 2000. - 2. - P.413-424.
15. Daly B.J. Numerical study of two fluid Raylaigh-Taylor instability // Phys. Fluids. - 1967. -
2. - P.297-307.
16. Meirmanov A. The Stefan problem / 1992. - Walter de Gruyter, Berlin-New York. - 244 p.
17. Белоцерковский O.M. Численное моделирование в механике сплошных сред / М.: Наука, 1984. - C.510-520.
18. Владимирова Н.Н., Кузнецова Б.Г., Яненко Н. Н. (1996) Численные рассчеты симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости // Некоторые вопросы вычислительной и прикладной математики. - Новосибирск, 1996. - C.186-192.
THE RAYLEIGH-TAYLOR INSTABILITY FOR THE MUSKET FREE BOUNDARY PROBLEM
O.V. Galtsev
Belgord State University,
Belgorod, 308015, Russia, e-mail: [email protected]
Abstract. Joint motion of two immiscible incompressible liquids in porous media is studied. Liquids have different densities and initially separated by a surface of strong discontinuity (free boundary). Results of numerical simulations for exact free boundary problems on the microscopic level for the absolutely rigid solid skeleton and for the elastic solid skeleton of different geometries are discussed. Problems under consideration have the natural small parameter which is the ratio of average pore size to the size of the physical domain. If the parameter e formally tends to zero, solutions of microscopic equations describes the averaged motion of liquids both in the case of the Muskat problem with the absolutely rigid solid skeleton and in the case of the viscoelastic Muskat problem with the elastic solid skeleton. The last model preserves a free boundary during the motion while in the first model, instead of the free boundary, the mushy region appears occupied by a mixture of two fluids.
Key words: Musket’s problem, the Rayleigh-Taylor instability, Stokes and Lame’s equations.
