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

MS С 76М25

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ СОВМЕСТНОГО ДВИЖЕНИЯ ДВУХ НЕСМЕШИВАЮЩИХСЯ НЕСЖИМАЕМЫХ ЖИДКОСТЕЙ В ПОРИСТОЙ СРЕДЕ НА МИКРОСКОПИЧЕСКОМ УРОВНЕ

O.A. Гальцева, О.В. Гальцев

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

Аннотация. Работа посвящена численному исследованию моделей фильтрации двух нееме-шивающихся несжимаемых жидкостей различной плотности, разделенных свободной границей в iiopoyiipyi'OM пространстве. Приводятся численные результаты по разработанным а.;п'орит-мам в случае структуры поровох'о пространства в виде изолированных капилляров и несвязных элементов пористой среды.

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

1. Введение. Проблемы моделирования физических процессов различной природы в различных средах возникают в механике жидкости и газа, в механике твердого те.иа, электродинамике и многих других областях. При этом общей проблемой является соотношение микро-и макроскопических подходов их описания. Часто требуется построить модель среды, локальные свойства которой резко меняются, поэтому удобнее перейти от микроскопического ее описания к макроскопическому, то есть рассматривать усредненные характеристики такой среды. Во многих случаях рассматриваемые физические процессы в сильно неоднородных средах описываются уравнениями с частными производными, причем сильная неоднородность этих сред приводит к дифференциальным уравнениям с резко изменяющимися коэффициентами. Такие задачи возникают в теории упругости и гидродинамике, в теории гетерогенных сред и композитных материалов, теории фильтрации и других задачах физики и механики. Непосредственное численное решение таких задач, как правило, затруднительно даже на современных ЭВМ. Поэтому возникает вопрос о построении моделей для сильно неоднородных сред, приводящих к более простым дифференциальным уравнениям, которые называются усредненными. Часто такие дифференциальные уравнения имеют постоянные коэффициенты. Усредненные уравнения позволяют определить с большой точностью эффективные характеристики первоначальной среды. Это условие обеспечивается основным требованием, которому должны удовлетворять усредненные уравнения - близость решений соответствующих краевых задач для исходных и усредненных уравнений. Математическое описание сильно неоднородных сред часто основано на предположении о

наличии у таких сред какой-либо унорядочешюй микроструктуры (например, периодической, квазинериодической, случайной однородной и др.).

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

Уравнения нороунругости, полученные К. фон Терцаги |1| и М, Био |2| долгое время являлись общепринятыми и служили основой дня решения практических задач нороунругости. Эти уравнения учитывают перемещение не только жидкости в норах, но и твердого сколота. Предлагаемые К. фон Терцаги и М. Био модели называют феноменологическими: в них постулируются свойства смеси твердой и жидкой компонент. Позже, ряд авторов (Р. Барридж и Дж. Келлер |3|, Э. Сапчес-Палопсия |4|, Т. Лови) предложили вывод уравнений нороунругости на основе основных законов механики сплошных сред и методов усреднения. Это было вполне естественно, сначала описать совместное движение упругого сколота и жидкости в норах на микроскопическом уровне, используя классические законы механики сплошных сред, а затем найти соответствующие аппроксимирующие модели с помощью теории усреднения (усредненные уравнения).

Так, совместное движение в области П описывалось ими математической моделью [5]

+ У • (рьОь-хЬ + ^-х)^) =рР, (1)

^ + V■(pv) = 0, (2)

где V • и - дивергенция и, матрица а ® Ь определяется как

(а ® Ь) • с = а(Ь • с),

для векторов а, Ь, и с, х есть характеристическая функция порового пространства П/, Р/ и Р5 - тензоры напряжения жидкой и твердой компонент соответственно, V -скорость среды, р - плотность среды и^ - заданный вектор распределенных массовых

Уравнения (1) и (2) понимаются в смысле теории распределения (как соответствующие интегральные тождества) и содержат динамические уравнения дня жидкой компоненты

dv „ ,_, ¿р „

в П/ при Ъ > 0, динамические уравнения для твердой компоненты

¿V „ ,_, ¿р „

р— = У-Р, + р^, -£+рУ-*; = 0 аъ аъ

в П8 при Ъ > 0, и условие непрерывности нормальных напряжений

(Рв - Р/) • п = 0

на общей границе «поровое пространство - твердый скелет» Г(£), оде п есть единичная нормаль к Г(£).

Приведенная задача сильно нелинейная и содержит еще одну неизвестную величину

- границу раздела норового пространства и твердого сколота. Главный постулат здесь

- твердая и жидкая компоненты не смешиваются. Таким образом неизвестная (свободная) граница Г(£) является поверхностью контактного разрыва [5], которая определяется из задачи Коши

для характеристической функции X в области П при £ > 0.

Дня описания совместного движения двух неоднородных жидкостей в упругом сколото динамическая система уравнений дополняется уравнением переноса для плотности р£(х,£) смеси жидкой и твердой компонент:

В частности, эта система описывает движение двух несмошивающихся жидкостей различной (и постоянной) плотности р++ и р- в несжимаемом твердом скелете постоянной плотности р8, если уравнение (4) дополнить начальным условием

где есть область занятая жидкостью плотности р±± в начальный момент времени.

Очевидно, что даже такое упрощение по делает задачу намного легче. Более того, если задачу (1), (3)-(5) удастся решить, то эта математическая модель будет бесполезна для практического применения, так как функция X меняет свои значения от 0 до 1 в масштабе нескольких микрон в то время как задача в целом рассматривается в области в несколько десятков (сотен) метров. Поэтому наиболее подходящим здесь будет усреднение рассматриваемой модели. Но в этом случае задача (1)-(3) становится совершенно неразрешимой. Чтобы получить что-то подходящее для практических нужд и все еще разумное с теоретической точки зрения, воспользуемся схемой, предложенной в |3,4| и линеаризуем основную динамическую систему.

А именно, аппроксимируем характеристическую функцию X жидкой части П/ ее значением в начальный момент времени

~Т7 ^¡7 • V \ • т 0. \(х.0) \„(х)

(3)

р(х, 0) = р8, х е П, р(х, 0) = р±, х е

(5)

где р/ и р, есть плотности жидкости в порах и твердого скелета соответственно, и

Р/ = 2р V) - р I, (6)

Р, = 2А Р(ж, т) - р I. (7)

Здесь 0(я, V) есть симметричная часть У^ I - единичный тензор, т вектор перемещения среды, р - динамическая вязкость, V - объемная вязкость, и А - постоянная упругости Ламе.

Предположение 1. Пусть х(у) есть 1-периодическая по переменной у функция, такая, что х(у) = У ^^ 7 С У, х(у) 0. у € У8 = Y\Yf, У - единичный квадрат Я2.

1) Множество У/ открытое и 7 = ЗУ/ П ЗУ8 есть Липшицева поверхность.

2) Пусть У/ есть периодическое повторение в Я2 элементарной ячейки £ У/. Тогда У/ есть связное множество с лппшпцевой границей 0У£, которая является периодическим повторением границы £ 7.

3) П С Я2 - ограниченная область с Липшицев ой границей Я = ОП и П/ = П П У£ есть норовое пространство, П^ = П \ Щ - твердый скелет, I '0 I'; дЩ П д0,£3 - общая граница « твердый скелет - норовое пространство».

Пусть ((х) характеристическая функция области П. Тогда х£(х) = С(х)х(х/£) будет характеристической функцией жидкой области П£, В безразмерных переменных

х т Ь ^ Г

ж ->• — , го — , £->•-, .Р ->• — , Ь Ь т д

где Ь есть характерный размер физической области, т есть характерное время физи-д

Тогда динамическая система примет вид

02т

атд£-^- = Ч- Р+£^Р, (8)

От

Р = Xе Щх, + (1 - \ )'1 ,\ Щх, го)-р1, (9)

У ■ т = 0 . (10)

В (8) - (10) £ = 1/Ь есть безразмерный размер пор, I - средний размер пор,

Ь 2р 2А

о 1 ^и Г п '

дт2 ' м тЬд р0 ' Ьдр0 '

о£ = в/ х£ + о, (1 - х£) >

в/ и О« _ безразмерные плотности жидкости в порах и твердого скелета, соотнесенные р0

Различные частные случаи линеаризации (1) - (3) изучались многими авторами: Бьюкеиен - Гилберт-Лии |6,7|, Бакингем |8|, Барридж-Келлер |3|, Клопиу-Ферри-Гил-берт-Микелич-Паоли |9-11|, Лови |12|, Hrvcrceiir |13|, Саичес-Хыоберт |14|, Сапчес-Палеисия |4|.

Наиболее ио.нно задача усреднения общей линеаризованной системы для сжимаемых сред была исследована в работах A.M. Мейрмаиова |15,16|,

В частности, в |15,16| была предложена классификация физических процессов и физических сред в зависимости от значений безразмерных критериев

lim ат(е) = т0 , lim aß(e) = ß0 , lim ад(е) = A0 .

e\0 e\0 e\0

Дня очень медленных процессов, таких как фильтрация жидкости, скорость среды составляет всего 3-6 метров в год. Поэтому характерное время процесса очень велико и ат ~ 0. Для быстропротекающих процессов таких, как акустика или гидравлический удар, ат ~ 1, ил и ат ~ то,

В этом случае инерционными слагаемыми в (6) можно пренебречь и ограничиться уравнением

V- P + q£F = 0 . (11)

При описании совместного движения двух неоднородных жидкостей в упругом сколото система уравнений (10) - (11) дополняется уравнением переноса

для плотности р£ смеси жидкой и твердой компонент и начальным условием

р£(х, 0) = ps, x е Qs, Ps(x, 0) = р±, x e Q± . (13)

Самым простым случаем системы (10) - (12) является случай, когда твердый сколот предполагается абсолютно твердым толом. Он характеризуется равенством

А0 = то.

Соответствующая система уравнений состоит из уравнений Стокса

V- v = 0 , (14)

V- {aß D(x, v) - p I) + Qf F = 0 (15)

для скорости v и давления p жидкости в области Qf при t > 0 и равенства

v=0

в твердом скелете Qs (более точно, система (14)-(16) получается из системы (10)-(12) предельным переходом при а\ ^ то).

Отмстим, что нетривиальное усреднение системы (14) - (16) возможно только при условии

ам = £2р0, р0 = const > 0 . (17)

(см. |15,16|),

При описании совместного движения двух неоднородных жидкостей в абсолютно твердом сколото система уравнений (14)-(17) дополняется уравнением переноса

^ + v ■ Vpf = 0, ж G О/ (18)

для плотности р£ смеси жидкой и твердой компонент и начальным условием

Pf (x, 0) = р±, x е . (19)

В исходной постановке дня точных уравнений динамики (1), (3), (4), задачи (12)-(3) и (18)-(19) эквивалентны и дня упругого сколота, когда равенство (15) может и не выполняться. Но при линеаризации динамических уравнений мы потеряли важное свойство границы r(t) «поровое пространство - твердый скелет» быть поверхностью контактного разрыва. Грубо говоря, жидкая и твердая компоненты среды теперь могут «протекать» через ее аппроксимацию - границу Г0, Поэтому на участках, где жидкая компонента «втекает» в поровое пространство, необходимо задавать значение плотности, что не вызвано существом дола.

Вообще говоря, такая постановка задачи, когда решением является искомая функция pf (а не комбинация р£, в которой плотность твердого скелета ps известна) более естественна и привычна дня исследователя. Поэтому, чтобы «подправить» последнюю постановку задачи о транспорте массы, воспользуемся диффузионным приближением

^L + vVpf = D0Apf (20)

с достаточно малым коэффициентом диффузии D0, который впоследствии можно устремить к пуню.

Данное уравнение допускает два вида краевых условий на границе Г0, не требующих знания искомой плотности pf

(DoV pf - pf v) • n = 0 , (21)

DoV pf • n = 0 , (22)

где n - вектор внешней нормали к границе Г0,

Условие (21) является точным дня линеаризованных динамических уравнений с фиксированной границей раздела Г0, Здесь под «точным условием» мы понимаем условие, запрещающее перетоки жидкости из норового пространства в твердый сколот или перетоки из твердого сколота в норовое пространство. С другой стороны, условие (22) является точным дня точных понипоаризоваппых динамических уравнений с

неизвестной границей раздела Г(£). В этом случае условие на сильном разрыве для уравнения диффузии (22) примет вид |5|

Pf (Уп - V • п) + Д,У Pf • п = 0 , (23)

где Уп есть скорость движения поверхности Г(£) в направлении нормали п. Поскольку данная поверхность Г(£) есть поверхность контактного разрыва, то

Уп = V • п

и мы приходим к (22).

2. Разработка алгоритмов численного решения математических моделей

1 и 2. Задача численного усреднения математической модели 1 совместного движения жидкости и упругого сколота решалась для двух различных геометрий элементарной ячейки норового пространства: несвязные капилляры (Рис. 1 - геометрия «а») и несвязный упругий сколот (Рис. 1 - геометрия «б»).

Уз

а) б)

Рис. 1. Структура поровшх) пространства: а) в виде изолированных капилляров, б) в виде несвязных элементов твсрдох'о скелета.

Предположение 2.

1) Пусть «жидкая компонента» У^ квадрата У = (0,1)2 С Я2 является открытым множеством, а К = У/Уу - «твердая компонента» в У. 7 = <91/ П дУ8 -липшицева граница.

У = Yf и У8 и 7 и дУ , У = (0,1) х (0,1), еУ = (0,е) х (0,е).

2) Пусть У^ есть периодическое повторение в Я2 элементарной ячейки еУf, с лип-шпцевой границей дУу, которая является периодическим повторением границы е^ У/

- периодическое повторение в Я2 элементарной ячейки еУ8.

3) Пусть П € Я2 - ограниченная область с Липшицев ой границей Б, П = Щ и Г£ и Q£s, где Щ = П П Уу - норовое пространство, = П \ Щ - твердый скелет, п I ; дЩ П дП£3

- общая граница «твердое тело - норовое пространство».

Модель 1. Совместное движение вязкой жидкости в норовом пространстве и упругого сколота на микроскопическом уровне описывается системой

д*щ£

V • (х^оЮ>0г, —) + (1 - АсЩх, ги£) - р!) + р^Р = 0 , (24)

V ■ 'ш£ = 0 , (25)

др£ д,ш£

+ — • V Р/ = о, Р£ = хеР/ + (1 - \ ;•/',, (26)

которая состоит из стационарной системы уравнений Стокса

V • (р0 Щх, -Ур£-р}Г = 0,У-го£ = 0, (27)

для скорости д,ш£/дЬ и давления р£ неоднородной несжимаемой жидкости в области Щ£ при Ь > 0 и стационарной системы уравнений Ламэ

V ■ (АоО(ж, ю£)) -Vp£ - р£ ^ = 0, V ■ ю£ = 0 , (28)

для перемещения и;£ и давления р£ упругого несжимаемого скелета в области Щ£ при Ь > 0, транспортного уравнения (26). На общей границе «твердый скелет - поровое пространство» Г£ выполнено условие непрерывности перемещений и нормальных напряжений.

Система (26)-(28) дополняется начальным и граничным условиями

,ш£(х, 0) = 0, х е Щ , (29)

ю£ = 0, х е дЩ, г> 0, (зо)

р(х, 0) = ро(х), х е Щ£ (31)

и условием нормировки

/ р(х,г) йх = 0 , г> 0 . (32)

ип

Процесс численного усреднения (е \ 0) моделируется увеличением количества капилляров дня геометрии «а», или количества элементов упругого сколота для геометрии «б». Таким образом, можно считать, что при достаточно малом е система (26)-(31) описывает задачу вязкоупругой фильтрации.

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

величина шага по времени, и - число шагов, известно поле скоростей V = д'ш/дЬ, дав-рр

можпо представить в следующем виде:

Этап 1. Решается система уравнений Ламэ (28) с граничными и начальными условиями для перемещения ш (30), (31). На первом шаге по времени находим промежуточное значение решая уравнение

Ао Д' - Ур = - раГ , (33)

дополненное граничным условием

W Isur^ = 0

В данном случае условие несжимаемости (25) не удовлетворяется. Поэтому будем искать поправки к нолю перемещения и давления в следующем виде

w = W + W, p = p + ps, (34)

где p|t=0 = 0 .

Подставив (34) в (28), получим уравнения для коррекции давления и перемещения упругого сколота

Ао AWs - Vps = 0 , (35)

V ■ Ws = —V ■ Wos, (36)

Ws 1 Sur£ =0 •

Выразим градиент Vps из системы (35) и применим дивергенцию к обеим частям, получим уравнение для поправки давления

Aps = Ао V ■ (AWs) • (37)

Учитывая, что V ■ (AW) = A(V ■ W) и равенство (36), уравнение (37) примет вид

Aps = —Ао A(V ■ Wos), (38)

которое дополним с.:юдующими граничными условиями

Ур5 • = 0 ,

Ур8 • п|Г£ = -АД' • п|Г£.

Найдя р3, легко решим (35) и, соответственно, найдем окончательные значения перемещения и давления в (34).

Этап 2. Используя найденные значения и р3, находим нормальные напряжения на границе Г

(Ао»(ж, ш) - р1)п = А , (39)

где п есть единичный вектор нормали к границе Г£.

Этап 3. Так как на общей границе между упругим сколотом и жидкостью выполняется условие непрерывности перемещений и нормальных напряжений, то решаем систему уравнений Стокса (27) в Щ£, учитывая граничное условие на об щей границе Г£

(роО(х, V) - р1)п = А , (40)

где правая часть равенства известна с предыдущего этана.

В виду того, что р уже известно из (40), уравнение для V = д'ш/дЬ на первом шаге будет иметь следующий вид

ро ЛV - Vp = р£Е , V|S = 0 . (41)

Затем аналогичные действия повторяются, а именно, находим поправку дня ноля скоростей и давления

V = V + V)£, р = р + р£ . (42)

Подставляя (42) в (27) и в уравнение неразрывности, получим уравнения дня коррекции скорости и давления

ро ЛV£ -VP>£ = 0 , (43)

V- V)£ = -V- V , (44)

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

V|5 = 0, (ро О (х, V) - р>£ I) п = 0 на г£ .

Выразим градиент Vp£ и применим дивергенцию к обеим частям в (43). Получим уравнение дня поправки давления

Ар£ = -ро Д^- Vn) , (45)

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

^р>£ ■ п)18 = 0, Vp£ ■ п = -р0ДV ■ п на Г£.

Найдяр^ решим (43), Учитывая известные V V и р£ легко найдем V, р£. Нормальное напряжение на общей границе Г£ находится, как

В = (ро О (х, V) - р)п .

Этап 4. Используя уже известное значение скорости жидкости, находим значение

£

плотности на следующем временном слое (рП+1), решая численно уравнение

Рп+1 - Рп

^-^ = -V ■ Чр'} . (46)

Этап 5. На этом текущий цикл закапчивается. На последующих временных слоях алгоритм нахождения неизвестных в упругом сколото решаем уравнение

АоД'П - Ур= -р.^ (47)

с граничным условием

(А0О(ж,шп) - Ур) • п = В ,

найденным на Этапе 3, и условием на Б

ш X = 0.

1 for каждой ячейки сотки do

2 if ячейка принадлежит область упругого сколота then

3 находим w из уравнения (33).

4 При p|t=0 = 0 ищем Ws, решая систему уравнений (35)-(36), Получаем ps из (37).

5 if Ap's = —Л0 A(V ■ Ws) then

6 решаем (34).

7 end

8 Находим нормальное напряжение на Г из (39).

9 end

10 else

11 Ищем v£ с учетом найденного значения нормального напряжения на общей границе «упругий сколот - норовое пространство», соблюдая аналогичный порядок действий, что и дня расчета перемещения в упругом сколото.

12 end

13 Находим р™+1 (значение плотности на следующем временном слое).

14 end

Алгоритм 1. Решение системы уравнений (26)-(28). Поправки дня давления и перемещения будем искать как

wns = Wn + Wn, p = p + РП. (48)

Уравнения для коррекции Ws и ps будут

ЛoAWS — VPs = 0 , (49)

V ■ WS = — V ■ Wns (50)

со следующими граничными условиями

гЪП|^ = 0, (Ло О(х, гг8) - Vs I) п = 0 на Г£. Поправку для давления находим из уравнения

Дрп = -ЛоД^ ■ гп) (51)

Vps ■ V^ = 0, VVs ■ п|Г£ = -ЛоДгП ■ п на Г£.

Как и ранее, находим ггП, затем г'П и р и повторяем этапы 2-4. Процедура повторяется, пока не получим сходящееся решение.

Модель 2. Математическая модель 2 движения жидкости в абсолютно твердом скелете решалась численно в области П, которая имеет структуру элементарной ячейки как геометрии «а», так и геометрии «б». Движение жидкости в одном капилляре описывается системой уравнений Стокса

V-(ро е1 О(х, v£) - р£1)+ р£Е = 0 , X е П£ ,г е (0,Т), (52)

V- v£ = 0, х е П£ ,г е (0,т), (53)

которая дополняется уравнением переноса

^ + = рЦх,0)=р{р(х),хеЩ,1е(0,Т) (54)

р££

I р£(х, t)dx = 0 , (55)

где ро = а^/е1.

Выполняется следующее граничное условие

V" = 0 , х е 5 и Г£, t> 0 , (56)

где Б = дП,

Процесс численного усреднения (е \ 0) моделируется увеличением количества капилляров дня геометрии «а», или количества элементов твердого сколота для геомете

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

Этап 1. На первом временном слое находим промежуточное значение скорости V в Щ из уравнения

^Аг, = -р1Г, ъ\8игв = 0,р1 = р°г (57)

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

V1 = V + V, р1 = р + р. (58)

Находим поправки скорости из уравнений

2

АчЗ — Ур = 0 , (59)

У • V = -У • V , (60)

^Г* = 0 .

Ур

дня поправки давления

2

У-(А^). (61)

Учитывая У • (Д V) = Д(У • V) и равенство (60), найдем р

АР = ^ V • (Д «), (62)

зная которое, решим (59).

Используя уже известные значения V и V, найдем V1, р1 и вычислим (р2) на следующем шаге

¿^А = . ур1. (63)

т

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

2|

АЬп - Vр(га"1} = -рпГ , Ьп\3иГ£ = 0 . (64)

Находим поправки дня скорости и давления, как на предыдущем шаге

Vй = Vп + Vй , рп = р(п-1) + рп . (65)

Подставив (58) в исходное уравнение Стокса и уравнение неразрывности, получим

уравнения дня коррекции скорости и давления

2

Айп - Vрп = 0 , (66)

У • ¿п = -У • V"- . (67)

1

Арп = V • (АЬп) , (68)

Затем решим (65).

И, наконец, находим значение (рп+1), численно решая уравнение

„п+1 - 0п

-= -Vй • Урга . (69)

1 for каждой ячейки сетки do

2 for первого каждого временного слоя n =1 v из (57). do

3 Находим v из (59).

4 Находим p из (62).

5 Определим v1 и p1,

6 Ищем р2 из (63) для следующего шага по времени.

7 end

8 for второго и последующего временного слоя находим vn из (64). do

9 Находим поправки скорости и давления,

10 Решаем (65),

11 Ищем pn+1 из (69),

12 end

13 end

Алгоритм 2. Решение системы уравнений (52)-(54).

3. Численное решение модели 1. В качестве численного метода решения моделей 1 и 2 был выбран метод VOF (volume of fluid).

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

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

которая касается перемещения элемента жидкости в новое положение. В Лагранжевом случае сетка просто перемещается с рассчитанными скоростями элемента потока, тогда как в Эйлеровом случае необходимо вычислить поток жидкости через сетку. Этот поток требует усреднения параметров жидкости дня всех элементов жидкости, которые находятся в данной ячейки сетки.

Именно этот процесс усреднения присущ аппроксимации конвективного потока. Это самый большой недостаток метода Эйлера. Конвективное усреднение приводит к сглаживанию всех вариаций в потоке и, в частности, к размытию поверхности разрыва (свободной поверхности). Единственный способ преодолеть эту проблему на границе -это ввести некоторую специальную обработку, которая позволила бы избежать разрыва. Сравнения относительных преимуществ и недостатков уже существующих методов преодоления этой проблемы привели к новой технике, которая является простой, но эффективной. Это метод объема жидкости в ячейке (УСЖ),

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

В основе УОЕ-метода .нежит схема дробного объема жидкости дня отслеживания свободной поверхностью. В этой технике, функция р(х 1,х2^) равна единице если ячейка занята жидкостью полностью и нулю, если ячейка пуста, когда усредненное но клеткам вычислительной сетки, среднее значение р в ячейке равно дробному объему ячейки, занятой жидкостью. В частности, единичное значение р соответствует ячейке, заполненной жидкостью, а нулевое значение указывает, что ячейка не содержит жидкости.

р

Как отмечалось выше, процесс численного усреднения моделируется увеличением количества капилляров дня геометрии «а», или количества элементов упругого сколота для геометрии «б» (е \ 0) .

А

/ N V

( \

V )

К \ /

>ь

х% V

РЩ РШ

Рич Ш/+1

Ру Ри

Ри-1 РФ

Рш. I Р'-и

■ ]+Ш

' ]+1/2

П/2

' ]-3/2

¡-3/2 Ьг г-1/2

1+1/2 1+3/2

Рис. 2. Шаблон расчетной сетки. В силу выбранного нами метода решения задачи исследуемая область течения но-

крывается равномерной по х1 и х2 сеткой ячеек (смотрите рисунок 2, «а» и «б»)

где Л,1; - размер сетки, М2 - количество ячеек сетки, соответственно, в направлении х1 и х2 (точка с координатами (г,Э) совпадает с центром ячейки).

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

Здесь, как и в исходном методе УСЖ, будем использовать «шахматную» сетку. Это дает возможность четко интерпретировать каждую ячейку, как элемент объема, который характеризуется рассчитываемыми давлением и плотностью в его центре.

Дискретные значения переменных расположены в ячейке, как показано на рисунке 2 «б». Здесь свободная граница (г,Э) определяется как ячейка, содержащая ненулевое значение р и имеющая по крайней мере одну соседнюю ячейку (г ± 1, э) или (г,Э ± 1),

плотность в которой равна 0. Ячейки с нулевым р содержат вторую жидкость. Ячейки

р

первой жидкостью.

Кратко, основную процедуру решения модели 1 за один шаг но времени можно представить в виде трехэтаиной схемы:

(1) Используя явные приближения уравнений (27),(28) и, следуя алгоритму, изложенному выше, находим значения скорости и перемещения.

(2) Используя разработанный алгоритм, удовлетворяем уравнению неразрывности (25) в каждой ячейке, так как из-за изменения давления в одной ячейке, нарушается банане в четырех соседних клетках. В стандартном методе УСЖ предполагается, что на свободной границе раздела жидкостей выполняется условие сжимаемости среды. Так как нам необходимо учитывать поведение еще и упругого тола, то необходимо удовлетворять условие (25) в каждой ячейке всей области. Это небольшая модификация метода позволяет избежать непредусмотренных возмущений при переходе на новый временной слой.

р

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

Конечно-разностные аналоги пространственных переменных дня соответствующих производных, входящих в исходную систему уравнений, центрируются в соответствии с выбранным шаблоном. Например, слагаемые с градиентом давления вычисляются с помощью односторонних разностей но формулам вида

гН1) Н1 > 0; г = 0,1,..., М1; ]к2, Н2 > 0; э = 0,1,..., N2;

др _ р+ы - р^ др _ рл - р.

(70)

дх1 Н1 ' дх2 к2

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

еги

д'2и _ Щ+3/2,7 - 2 "¿+1/2,7 + »'¿-1/2,7 ^

5x1

д2и _ Щ+1/2,7+1 - 2%1/2,3- + Щ+1/2,7-1

5x2 Л2

(72)

Более по.нпые выражения конечно-разностных аналогов соответствующих слагаемых двухмерных и трехмерных уравнений движения приведены, например, у О.М. Бе-.иодерковского |17|, Используемая конечно-разностная схема аппроксимирует рассматриваемые уравнения с первым порядком точности но времени и со вторым порядком точности по пространственным переменным 0(Дт, Л2) и можно показать, что она устойчива. Подставим конечно-разностные формулы в исходную систему уравнений движения (24)-(26). Тогда поело простых преобразований получим их дискретные аналоги для х1 и х2 направлений соответственно.

Полученные разностные алгебраические уравнения, разрешенные относительно компонент перемещения и^ ¿+^,7, и2 ¿7+^2 и дополненные уравнением неразрывности, преобразуются к следующему конечно - разностному виду

Л (¿+3/2,7 - 2 ¿+1/2,7 + ^"¿-/2,7 .

Ло ( К2 +

1

+

¿+1/2,7+1 2 ¿+1/2,7 + ¿+1/2,7-1

Л2

-^--г Ра ¿+1/2,7 ^ '

л (и2 ¿+1,7+1/2 2 и2 ¿,7+1/2 + и2 ¿-1,7+1/2 Ло I ц +

и'2 ¿,7+3/2 2 и,2 ¿,7+1/2 + и,2 ¿,7-1/2 N

ц }

и1 ¿+1/2,7 и1 ¿-1/2,7 и2 ¿,7+1/2 и2 ¿,7-1/2

Л1 Л

2

^¿,7+1 ^,7 , га т-г /7л

-^- + Ра ¿,7+1/2 -^2 , (' I)

= 0 , (75)

1 +1 1 р1 р1

Ра ¿,7 Ра ¿,7 _ „ ¿+1/2,7 ¿-1/2,7

; --«'1 ¿+1/2,7 ^

р1 р1

ра ¿,7+1/2 ра ¿,7-1/2 , , - «'2 ¿,7+1/2-^- • ( < 6)

где п - номер шага по времени, и+ и2 - компоненты вектора перемещений, р3 - плотность упругого скелета.

Величины с дробными индексами относятся к границам ячеек, например

Wl ^ + Wl <Ш2 г^ + ^2 г,з + 1

™И+1/2,з — -2-' ^¿,¿+1/2 - -2- •

В случаях, когда нам необходимо определить значения функции в точках сетки, необходимо воспользоваться средним арифметическим, например

Р1+1/2,3 = + р1,.р\,з +1/2 = ^ '

Легко видеть, что эти конечно-разностные схемы (так и нижеследующие) аппроксимируют систему уравнений (26)-(28) с погрешностью порядка О(т, К2), где К = шах(К, К2). Таким образом, приведенная разностная схема имеет второй порядок точности по пространству.

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

шах

V П+1 - Vп

,п+1

< а.

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

На рис. 3 представлены результаты численного усреднения задачи (2б)-(31) по разработанному алгоритму для случая элементарной ячейки «б».

ЦЦ

1 = 860

1 = 1200

1 = 2600

1 = 4800

Рис. 3. Совместное движение жидкости и упругого скелета в виде изолированных капилляров.

4. Численное решение модели 2. Процесс аппроксимации системы уравнений (52)-(56) разностным представлением ничем не отличается от описанного выше, так как модель 2 является подмоделью модели 1, где отсутствуют перемещения скелета и, соответственно, на границе жидкости и твердого скелета выполняется условие непротекания.

Полученные разностные алгебраические уравнения, разрешенные относительно компонентов перемещения м»+1/2,7, ^+1/2 и дополненные уравнением неразрывности, преобразуются к следующему конечно-разностному виду

РФ2 ( '"'¿+3/2,7 ~ 2 Ц"+1/2,7 + иг-/2,з

^2

_ 2 I

¿+1/2,7+1 - "¿+1/2,7 ^ "¿+1/2,7-1

Л2

2

РП+1,7 _ Р"

¿7

+ Р?+1/2,7 , (77)

2 ^I™ _ О,,« I ,

^'¿+1,7+1/2 ^ 1м+1/2 +

Л?

+

''"7+З/2 2 г'"7+1/2 + г м-1/2 \

Р™ _р™

ы

+ Р^+1/2 ^2 , (78)

л«" _ 7/™ 7 _ 7

"¿+1/2,7 "г-1/2,7 , ^»,7+1/2 %-1/2

Лт /ь

0, (79)

Рг,7 Рг,7 , /г Рг+1/2,7 Рг-1/2,7 ,, Рг,7+1/2 Рг,7-1/2 , , --- = -«¿+1/2,7-^--¿>7+1/2-^- • (Ь0)

При

Зс1М6Н6 дифференциальной задачи на конечно-разностное иредстав.:1ение необходимо обратить внимание на приближение граничных условий, так как их конкретное приближение прямо влияет на корректность метода и устойчивость схемы, а также скорость сходимости. Сформулируем граничные условия, введя серию фиктивных ячеек (так, чтобы каждая вычисляемая точка становилась внутренней и описанный алгоритм сохранялся дня всех ячеек). Дня первого порядка аппроксимации достаточно одного с.ноя, дня второго порядка аппроксимации - два слоя. Так в случае, когда боковые стенки есть твердая поверхность, то условие прилипания и ненротекапия представляется в виде

V-1/2,7+1/2 = 0 (условие прилипания) , " г-1/2,7 = 0 (условие непротек ания). Из условия прилипания можно получить

0,7+1/2 - —+ —+ т—¡^Ъ + .

Дня расчета ноля плотности в «приграничных» ячейках рассматривается следующее представление (смотрите рисунок 4):

п+1,к+1 _ п

р03 = р03 — т

(Ри)П/2,3 , ^0,3+1/2 - 1/2

+ 0(к2). (81)

приведенные краевые условия дня расчета скорости и плотности в «приграничных» ячейках имеют порядок точности но пространственным переменным не ниже второго.

-1/2 1

I'

/ /

1

1/2

3/2

V0,j+í/l.

¡£р,;+1

^Яу+1/2

Ра,; Е • 1 Ро] и-1/2,] » Ри и1/2,]

?>„.,• 1 Ра,]-1

7 + 1

1-1

Рис. 4. Расчетная сетка у границы «жидкость упругое тело» («абсолютно твердое тело»).

Численное решение этой задачи в одном капилляре в абсолютно твердом сколото показало совпадение с результатами |18|, На рис. 6 можно увидеть гладкую свободную границу в капилляре в различные моменты времени.

Процесс численного усреднения (е \ 0) моделируется увеличением количества капилляров дня геометрии «а», или количества элементов твердого сколота для геометрии «б». Таким образом, можно считать, что при достаточно малом е система (52), (56) описывает классическую задачу Маскета.

Легко заметить, что дня одного и того же времени процесса в задаче Маскета возникает переходная фаза в случае абсолютно твердого сколота (см. |16|), в то время, как при вязкоунругой фильтрации свободная граница сохраняется.

Литература

1. Terzaghi К. Die Bereehnung dcr Durehlassigkeitsziffer des Tones aus dcm Verlauf der hyd-rodynamisehen Spaimungsereheinungen /7 Sitzung beriehte. Akademie der Wissensehaften, Wien Mathematieseh-Naturwissensehaftliehe Klasse. 1923. 132, №IIa. P.104-124.

2. Biot M. Theory of elasticity and consolidation for a porous anisotropic solid /7 .Journal of Applied Physics. 1955. 26. P.182-185.

3. Burridge R., Keller .J.B. Poroelastieitv eciuations derived from microstructure // .J. Acoust. Soe. Am. 1981. 70, №4. Р.1140-П46.

4. Sanchez-Palencia E. Non-Homogeneous Media and Vibration Theory /7 Lecture Notes in Physics. Berlin: Springer, 1980.

5. Овсянников Л.В. Введение в механику сплошных сред. Часть 2 / Новосибирск: НГУ, 1977. 69 с.

6. Buchanan .J.L., Gilbert R.P. Transition loss in the farfield for an ocean with a Biot sediment over an elastic substrate// ZAMM. 1997. №77. P.121-135.

7. Gilbert R.P., Lin .J.Z. Acoustic waves in shallow inhomogeneous oceans with a poro-elastic seabed // ZAMM. 1997. №4. P.l-12.

8. Buckingham M..J. Seismic wave propagation in rocks and marine sediments: a new theoretical approach /7 Underwater Acoustics. 1998. 1. P.299-300.

9. Clopeau Th., Ferrin .J.L., Gilbert R.P., Mikclic A. Homogenizing the acoustic properties of the seabed. Part II /7 Mathematical and Computer Modelling. 2001. 33. P.821-841.

10. Ferrin .J.L., Mikclic A. Homogenizing the acoustic properties of a porous matrix containing an incompressible inviseid fluids /7 Math. Meth. Appl. Sei. 2003. 26. P.831-859.

11. Mikclic A., Paoli L. Homogenization of the inviseid incompressible fluid flow trough a 2D porous medium /7 Proceedings of the AMS. 1999. 17. P.2019-2028.

12. Lew T. Acoustic phenomena in elastic porous media // Mech. Res. Comm. 1977. №4. P.253-257.

13. Nguetseng G. Asvmptotie analvsis for a stiff variational problem arising in mechanics // SIAM .J. Math. Anal. 1990. 21. P.1394-1414.

14. Sanchez-Hubert .J. Asvmptotie studv of the macroscopic behavior of a solid-liquid mixture // Math. Methods Appl.'Sci. 1980. '2. P.l-18.

15. Менрманов A.M. Метод двухмаештабнон сходимости Нгустсснга в задачах фильтрации и сепсмо-акустики в упругих пористых средах // Сибирский Математический Журнал. 2007. №3." С.645-667.

16. Мейрманов А. М. Вывод уравнений еейемоакуетики и уравнений фильтрации в упругих пористых средах через усреднение периодических структур // Труды семинара имени И.Г. Петровского. М.: Наука, 2008. С.178-238.

17. Белоцерковекий О.М. Численное моделирование в механике сплошных сред / М.: Физ-матлит, 1994. 448 с.

18. Dalv В..J. Numerical studv of two fluid Ravlaigh-Tavlor instability // Phvs. Fluids. 1967. №2.' P.297-307.

NUMERICAL SOLUTION OF TWO IMMISCIBLE INCOMPRESSIBLE FLUIDS SIMULTANEOUS FLOW IN A POROUS MEDIUM AT MICROSCOPIC LEVEL PROBLEM

O.A. Galtseva, O.V. Galtsev

Belgorod State University, Pobedy St., 85, Belgorod, 308007, Russia, e-mail: [email protected]

Abstract. Two models of filtration of two immiscible incompressible fluids of different density separated by free boundary in a porous elastic medium is studied numerically. Numerical results are given according to developed algorithms in the case when the structure of pore space has the form of isolated capillaries and disconnected elements of porous medium.

Key words: free boundary problem, liquid filtration, numerical averaging.