Осредненное описание процессов разноплотностной фильтрации и массопереноса. 1. Уровень пор Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Том 14 7, кн. 3

Физико-математические пауки

2005

УДК 532.546

ОСРЕДНЕННОЕ ОПИСАНИЕ ПРОЦЕССОВ РАЗНОПЛОТНОСТНОЙ ФИЛЬТРАЦИИ И МАССОПЕРЕНОСА. 1. УРОВЕНЬ ПОР

Д.Е. Демидов, А.Г. Егоров

Аннотация

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

Введение

Безопасное захоронение химических и радиоактивных отходов является сегодня ключевой экологической проблемой. Среди прочих методов широко используется захоронение радиоактивных отходов в подземных каменно-соляных формациях. В подземных водах вблизи таких формаций наблюдаются высокие концентрации соли (вплоть до предела насыщения), соответствующие большим значениям плотности (до 1200 кг/м3). В случае утечки из такого хранилища наиболее вероятным механизмом доставки загрязнений к поверхности является перенос подземными водами. При этом возникают большие градиенты концентрации, и. как следствие, встает вопрос о применимости в рассматриваемой ситуации классических линейных реологических соотношений, лежащих в основе стандартных моделей фильтрационного переноса примеси: (I) закона Дарси, связывающего скорость фильтрации и градиент напора, и (и) зависимости Фика между дисперсионным потоком примеси и градиентом концентрации. Тщательные лабораторные эксперименты, проведенные в последние годы [1, 2], дали в целом отрицательный ответ на оба вопроса. Оказалось, что при росте разности плотностей воды и рассола наблюдается как снижение скорости размазывания концентрационных фронтов, так и ухудшение фильтрационных свойств среды. Считается, что данный эффект вызван стабилизирующим влиянием силы тяжести и может быть учтен зависимостью коэффициентов дисперсии и фильтрации от градиента концентрации. Вопрос о виде соответствующей зависимости до сих пор открыт. Более того, сама возможность этого далеко не очевидна.

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

масштаба лабораторных экспериментов (однородные пористые среды) до масштаба полевых испытаний (микронеоднородные среды).

Настоящая работа посвящена первому из указанных этапов. В этом случае базовые уравнения на микроуровне представляют собой связанную систему уравнений Стокса и уравнения конвективно-диффузионного массопереноса. Постановка возникающей при этом задачи на микроуровне приводится в разд. 1. Далее применяется метод гомогенизации [3. 4] для построения соответствующей осредненной системы уравнений (разд. 2). Коэффициенты этой системы (эффективный тензор дисперсии и коэффициент фильтрации) оказываются функциями градиента плотности и скорости фильтрации. Для вычисления этих функций на отдельном элементе пористой среды формулируются так называемые задачи на ячейке. В следующих двух разделах задачи на ячейке решаются для простейших одномерных моделей пористых сред, представляющих собой связку идентичных капилляров, и для более сложной и «реалистичной» двумерной решетчатой модели пористой среды. В разд. 5 проводится сравнение полученных результатов с известными данными лабораторных экспериментов [5]. В разд. 6 проводится обсуждение развиваемой теории и намечаются перспективы дальнейших исследований.

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

На микроуровне (уровень пор) сопряженные процессы фильтрации и массопереноса разноплотностной жидкости описываются в поровом пространстве П; следующей системой уравнений:

ди дт . .

др

и=£ (2)

А4Аи'=§г ~Р9-> др д . . В . . , „

т + д-х{ир) + Т^и'р)-с1-лр = 0- (4)

Здесь и, т - горизонтальная и вертикальная скорости движения жидкости, р -давление, р - плотность жидкости. Вертикальная координата г направлена вдоль силы тяжести; заданные константы д, ^ , ¿т определяют ускорение свободного падения, вязкость жидкости и молекулярный коэффициент диффузии. Уравнения (1) (3) сохранения массы и импульса записаны с использованием приближения Буссинеска [6] и в предположении о том, что инерционными силами можно пренебречь. На границе Г жидкой и твердой фазы ставятся условия прилипания и отсутствия потока примеси

др

Гя : и = ги = — = 0. (5)

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

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

являются давление р, плотность р и средние скор ости (и), (ад) движения жидкости

Асимптотический анализ микроскопической системы уравнений (1) (5) при стремлении к нулю малого параметра е = 1/Ь является типичной задачей теории гомогенизации [3. 4]. Результатом должна явиться замкнутая система так называемых осредненных уравнений для нахождения макроскопических характеристик в «эффективной» (неструктурированной) среде, учитывающая исходную информацию о микроструктуре посредством коэффициентных зависимостей. Вид получающейся системы существенно зависит от характерной величины уо скорости движения жидкости. Именно эта величина определяет два важнейших параметра, управляющих исследуемым процессом на микроуровне. Один из них. число Пекле

ре = —1

Я '

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

с= 9 12

ру о

определяет отношение гравитационных и вязких сил. Характерное изменение плотности 6р на ячейке периодичности может быть в силу макроскопичности р выражено через изменение Др плотности во всем пространстве как 6р = Ар (1/Ь), а значит.

с _ 9&Р13 Ь

В наиболее общей ситуации

Ре ж 1, С ж 1 (е ^ 0) (6)

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

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

дх

расположения менее плотной жидкости над более плотной.

Р-

Ь\

X

I, ................

а)

Ь)

Рис. 1. Макроскопическая схема рассматриваемых процессов: а) вертикальное движение концентрационного фронта: Ь) горизонтальный стационарный пограничный слой между смешивающимися жидкостями различной плотности

2. Процедура гомогенизации

Осреднение исходной системы (1) (4) микроскопических уравнений проводится однотипно для обоих типов процессов. Поэтому рассмотрим ее подробно лишь применительно к описанию горизонтального стационарного пограничного слоя. Перейдем в (1) (4) к безразмерным величинам

и

г>о'

т

«о'

Р - Р-

р

Р = —, ро

ро

А«'о

1е2

(7)

Р+ - Р-

и будем рассматривать их как функции двух типов пространственных координат

X =

I

Z =

I'

е =

Ь

п =

Х£

Т1

(8)

называемых, соответственно, быстрыми (X, Z) и медленными (£, ц). Каждую из

П

цией быстрых координат. Учитывая, что

дх I \дХ дп

д_

дг

перепишем (1) (4) в виде

ди д\¥ д\¥ ~дх + ~дг+£1й[ + £~

1 7

, ди

_д_ д_

дг+ед£

дг1

1 дР дР

1 дР 1 (дР

(9)

(Ю) (Н) (12) (13)

В этих уравнениях (X, Z) € П;

д2 д2 Л = А + 2е—— + 2е2

дедZ

+ £2

д2

д-цдХ де

дг!2

г

0

2

Д

пых. Граничные условия (5) преобразуются к

тт „т дС дС дС 2 дС

Гг: и = МГ = 0, пл—+пх—+£пл—+£-пх—= 0. (14)

В дальнейшем нам понадобится пара соотношений, получающихся из уравнений (10). (13) их осреднением по ячейке периодичности. Учитывая, что в силу периодичности и, Ш и условий прилипания (14)

/ ди дШ \ Г

Га+Гг

найдем из уравнения (10)

Щ<№+^(Ц)=0. (15)

Аналогично, хотя и несколько более громоздко, осреднение (13) с использованием (10). (14) приводит ко второму соотношению

Будем теперь искать решение задачи (10) (14) в виде асимптотического ряда

е

С = С0(С, п) + еС1 (X, 2; С, п) + е2С(X, 2; С, п) + ..., Р = Ро (С, п) + еР1 (С, п) + е2Р2 (X, 2; С, п) + • • •,

(17)

и = ио(Х, 2; С, п) + еи1(Х, 2; С, п) + • • •, Ш = Шо(Х, 2; С, п) + еШ1 (X, 2; С, п) + • • • Ре С е

стантамп.

Подстановка (17) в (15). (16) дает в главных членах откуда следует

(Шо) = 0. (18)

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

е

¿№ + ¿<^>=0, (19) Ре №) Со) + ^ (<^> Со)) + Щ = 0, (20)

где эффективный поток Q примеси в вертикальном направлении определяется как

Q = Ре (ШоС1)

/дел дСо

\ д2 / дС '

Подставим, наконец. (17) в (10) (13) и приравняем нулю получающиеся при различных степенях е члены. В главном члене (порядка е-1) уравнение (12) указывает на тот ожидаемый факт, что давление распределено по гидростатическому закону

^=СС0. (21)

ео

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

дЦ0 д\У0 = дХ дг

дС1 дС1 дС0

Гг : и0 = УГо = о, пй + пх —± + пй = 0.

дZ дХ де

Эта задача сформулирована в пространстве быстрых переменных и служит для определения в нем функций и<о, ЭД'о, Р2, С^. Величины дРо/дц, дР1/д£, дСо/д£ не зависят от быстрых переменных и, таким образом, являются параметрами задачи па ячейке. При фиксированных значениях дСо/де> 0, (Ц)), ^о) ш (С1) ее решение единственным образом определяет постоянные дРо/дц, дР1/д£ и функции ио, ^о, С^; функция Р2 находится с точностью до аддитивной константы.

Легко видеть, что одновременное изменение С1 и дР1/д£ на соответствующие константы А(£, п) и @А сохраняет задачу на ячейке. Поэтому при ее решении без

( С1 ) = 0

задача полностью характеризуется двумя параметрами (ио), дСо/д£ и после нормировки

ио = (ио), ^о = Vг (ио), Р2 = Р (ио), ^ = -тт. (г/0) , - (Сх) = (Ре (г/о»-1 ^с, дР1

_1-с(с1> = -7гг (г/0)

переписывается в виде

V- V =0, (22)

Д-и = VP - П - свг, (23)

Ре; (V • Vc) - Дс + Vг - (V,) = 0, (24)

Гг : г?=0, (Рег с + г) = 0, (25)

дп

(V,) = (с) = 0, («х) = 1. (26)

Здесь е2 = (0,1)т, а через Ре; и И.а; обозначены локальные числа Пекле и Релея

дС

Рег = Ре (11 о), Каг = РеС—

де

Решение задачи на ячейке (22) (26). с одной стороны, определяет нелинейный аналог закона Дарен

др

= -МРеьКаг)(£/о>, (27)

а, с другой стороны, конкретизирует эффективный поток примеси из (20) в виде

дС

д = -Р±(Рег,11аг)^ (28)

дС

с эффективным коэффициентом диффузии

Заметим, что произвольная константа в С1 те влияет та определение в силу (26). Как видно, замыкающие соотношения (27), (28) оказываются существенно нелинейными: и эффективный коэффициент диффузии , и коэффициент фильтрации = 1/пх зависят как от средней скорости (через параметр Ре;), так и от градиента концентрации (через параметр Иа,). Эти соотношения совместно с уравнениями (19) (21) представляют собой искомую макроскопическую модель для нахождения неизвестных величин Ро , (ио), (Ш.), Со, Q. Для наглядности выпишем эту модель в исходных, размерных переменных, сохраняя за компонентами средней скорости прежние обозначения и, т:

ди ди> дх дх

д'Р п

---рд = 0,

дх

д д д др

(29)

/2г &Р

и =--— .

р дх

Безразмерные коэффициенты эффективной диффузии и фильтрации являются, в отличие от классической погранслойной модели (см., например, [9]), функциями локальных чисел Пекле и Релея

и1 й/4 др

г = Т~' г =

дх

В классическом случае К^ - постоянная величина, К^ = К*!, а зависит лишь от локального числа Пекле. Соответствующий переход можно осуществить, положив в (23) Иа, = 0. При этом задача па ячейке расщепляется па две. Гидродинамическая задача (22)-(23) та зависит от Ре, и С. Решение этой задачи определяет константу пх = 1/К^ и поле скоростей У. Далее концентрационная задача (24) решается па фоне заданного поля V и дает искомую зависимость от Ре,. Обычно эту зависимость представляют в виде

Т>± (Ре,, 0)= С± + VI (Ре,), VI (0) = 0

и говорят о величинах = V± (0,0) и V0L = V± (Ре,, 0) — как о коэффициенте извилистости и коэффициенте дисперсии трассера соответственно. Далее будем придерживаться этой терминологии.

Интересно рассмотреть и другую предельную ситуацию, Ре; ^ 0. В этом случае задача на ячейке также расщепляется. Концентрационная задача для функции с = Ре; с

д

Д5 = 0, — (с + г)

дп

0

решается независимо от гидродинамической. Ее решение определяет эффективный коэффициент диффузии = 1 + (дс/дг) равным коэффициенту извилистости С.1. При задан ном с гидродинамическая задача

V •*?=(), Ы=УГ -С,сё, (С; = , (31)

= 0, (уг) = 0, (ух) = 1 (32)

линейна. Поэтому искомая величина пх является линейной комбинацией

,0

Пх = пх

1С; (33)

с константами , п1, определяемыми решениями, соответственно, задачи (30), (31) с С; = 0 и задачи

V-1 = 0, ДУ1 = УТ1 - П1 - О/ввг, (34)

-Чг. =0, (-1) = (-1) =0. (35)

В общем случае п1 = 0, и (33) является неклассическим обобщением закона Дарси вида

/от" дх 7Г° 1-1 дг

Однако наибольший интерес, как отражающие свойства реальных пористых сред, представляют «симметричные» пористые структуры, остающиеся инвариантными относительно преобразований вида 2 ^ и X ^ -X. Для таких структур замена

2 г с, V1, Т1 ^-Т1, п1 ^—п! (37)

П 1 = 0

а значит, для симметричных структур дополнительный фильтрационный поток (второе слагаемое в правой части (36)), индуцированный градиентом концентрации, невозможен. Такого типа среды в пределе Ре; ^ 0 независимо от значения Иа; описываются та макроуровне классическим законом Дарен с К± = = 1/п0 и стандартным уравнением массопереноса с постоянным и равным коэффициентом диффузии.

Аналогичные (37) замены показывают, что симметричные пористые структуры при любых Иа; инвариантны относительно изменения направления течения жидкости на противоположное:

(Ре;, Иа;) = (-Ре;, Иа;), (Ре;, Иа;) = (- Ре;, Иа;).

Поэтому в дальнейшем можно ограничиться рассмотрением случая Ре; > 0.

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

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

др

и = 0, ги = 1>0, ---рд = О,

* (38)

др др д др ~7Г7 + --ат—Т)п— = 0.

дг дг " дг

Эффективный коэффициент диффузии Х>ц, как и раньше, является функцией чисел Пекле и Релея. Необходимо лишь строить число Пекле не по горизонтальной, а по вертикальной скорости • Задача на ячейке (22)-(26) сохраняется, с заменой условий } = 0, («х) = 1 на } = 1, («х) = 0, а V] подсчитывается по формуле

Для симметричных пористых структур Х>ц, как и прежде, не зависит от знака Ре;, т. е. направления (вверх, вниз) движения фронта. Аналогично предыдущему определяются также коэффициент извилистости £ц и коэффициент дисперсии трассера р0(Ре;)- В общем случае они не сов падают с V0 (Ре;). Лишь в

дополнительном предположении об «изотропности» (инвариантности относительно преобразования (X, У) ^ (У, X)) микроструктуры £ц = . Но даже в этом случае V° (Ре;) = Р°(Ре;).

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

<«?-с> = ±(г?.Ус2) = ±<У(г>с)2)=0

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

V = (Ре; с + г)|2^ . (39)

Аналогичные выкладки после скалярного умножения (23) на V и интегрирования по ячейке периодичности дают

К-1 = (|^х|2 + |У«г12) - Ил; (с^} . (40)

3. Решение задачи на ячейке в одномерном случае

Конкретизация зависимостей , V]! от Ре;, Иа; предполагает задание микроструктуры пористой среды. Рассмотрим в качестве первого примера пористую среду в виде связки одинаковых, вертикально расположенных капилляров щелей раскрытия 21. В этом простейшем случае, очевидно, имеет смысл говорить лишь об определении продольного эффективного коэффициента диффузии V].

Решение задачи на ячейке здесь зависит только от поперечной координаты X, в силу чего Vх = 0, а функции (X), с(Х) удовлетворяют следующей системе обыкновенных дифференциальных уравнений

Vй = —п — Иа с, с'' = V - 1

и условиям

(v) = 1, c'(±l) = v(±1) = 0.

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

vIV +Ra v = Ra (41)

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

v(±1)= v'''(±1) = 0 (42)

для нахождения v(X). Условие (v) = 1 выполняется автоматически, в чем можно убедиться интегрированием (41) по интервалу [— 1,1] с использованием второго из условий (42). Продольный эффективный коэффициент диффузии определяется через решение задачи на ячейке формулой

D|| = 1 - Pe2 ((v - 1)c) .

Линейная задача (41), (42) допускает явное решение. Вид его, однако, громоздок, поэтому ограничимся асимптотическими представлениями. При Ra ^ 0 имеем

г. = -(1-Х2), С=--(Х4-2Х2 + -), 2?||=1 + —. 2 v ' ' 8 V 15 J 11 105

При Ra ^ ж имеем

v = 1 + e~s (sin s — cos s), с = _e-s (sin s + cos s),

Ra

1 —1/4

П, = 1 + -L Pe2 Ra-3/4, s = (1-Х2).

II 2 A/2 2 A/2 v '

v c Ra

рис. 2. В силу симметрии профили изображены для правой половины капилляра. Ra

щины порядка Ra-1/4, в котором резко изменяются скорость и концентрация. В

vc

стей потока в конечном итоге подавляет его дисперсивность. В пределе Ra ^ ж она отсутствует вовсе и эффективный коэффициент диффузии совпадает с молекулярным.

Обсудим более подробно свойства эффективного коэффициента диффузии, представив его предварительно в виде

D|| = Z|| + D0 (Pe) f (Ra). (43)

Коэффициент извилистости Z|| в рассматриваемом случае равен единице, продольный коэффициент дисперсии трассера

^(Ре)^Ре2

в точности определяет тейлоровский режим дисперсии, корректирующий множитель

f 105 К п \

/II = —— (('f - !)с)

2 1.5 1

0.5 0

-0.5 -1

Рис. 2. Профили V и с при различных Иа для связки капилляров щелевидпой формы

Рис. 3. Решетчатая модель пористой среды

учитывает влияние силы тяжести на дисперсивность. В силу нормировки /ц (0) = 1. С ростом И,а функция /ц монотонно убывает, выходя при И,а ^ то на асимптотический режим

/ (Еа) ^ Щ, Кд-з/4 _ 11 4 \/2

Интересно отметить, что множитель /ц и коэффициент Х>ц определены и при отрицательных значениях И,а вплоть до И.а = Касг < 0. Величина Касг = -31.285 есть наибольшее собственное число задачи (41), (42), подсчитываемое как максимальный корень уравнения

tg Асг = th Ас

А„г — Рас

При уменьшении И,а от нуля до Касг функция /ц (И,а) монотонно растет от единицы до бесконечности. Это означает, что даже при отрицательных значениях др/дг рассматриваемый процесс остается микроскопически устойчивым и может быть описан макроскопическим уравнением (38). Лишь при достижении др/дг критического значения, равного

др дг

-31.285^,

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

Аналогичные изложенным выше результаты могут быть получены для подобного типа сред с иной формой капилляров. Например, замена щелей круглыми капиллярами радиуса I (модель 2) не изменяет общей структуры (43). Коэффициент извилистости в этом случае по-прежнему равен единице, продольный коэффициент дисперсии трассера задается формулой

характер поведения /ц на бесконечности - формулой

/ц (11а) - 11а~3/4 (11а ->■ оо), 2

а критическое значение критерия Релея равно Касг = -104.36.

4. Решение задачи на ячейке в двумерном случае

Более сложной (и реалистичной) является еще одна рассматриваемая в данной работе решетчатая микромодель пористой среды. Она показана на рис. 3. Считается. что двумерная решетка одинаковых капилляров щелей повернута под углом 45° к вертикади. Раскрытие щелей равно 21, расстояние между центрами соседних щелей 2(а + 1)/. Параметр а определяет пористость

1 + 2а

т = -7г

(1 + о)

рассматриваемой системы. В представленных далее расчетах он выбирался равным 4, чтобы обеспечить типичное для пористых засыпок значение т = 0.36.

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

1. перейти от исходных (Х^) координат к координатам (у1,у2)

гу _ У2 - У1 V _У1 + У2

согласованным с осями решетки:

2. сформулировать гидростатическую часть задачи не для скорости и давления, как в (22)-(23), а в терминах функции тока (ф) и завихренности (ш)

дф дф дУ1 д'02

г'1 = тг—, г>2 = - —, --—-;

ду2 ду1 ду2 ду1

3. отыскивать вместо периодической функции с линейную комбин ацию Ре с + г.

Заметим, что именно функция Ре с + г, как можно видеть из (17), имеет простой физический смысл. За вычетом константы она совпадает с точностью до нормирующего множителя с главным членом в распределении плотности на ячейке периодичности. Подчеркивая это обстоятельство, будем обозначать далее в этом разделе статьи линейную комбинацию Ре с + г через р ■

В терминах р, ф, ш задачи на ячейке принимают вид

Дф = ш, (44)

дУ2 дУ1 дУ1 дУ2

где постоянная к равна нулю для ±-задачи и единице для || -задачи, = И.а;/Ре;. На твердой части границы Г8 = Га + Гв + Гс + Гд выполняются условия непроницаемости

г . ^Р = ^Ф = 0

дп дп

и задается значение функции тока

ф1га+Гс =0, ф|гв = -ф1гв = ф*, И: ^|Гв+гв=0, ф\Гл = -ф\Гс=ф*,

\/2 (1 + 2а)

ф*

1 + а

Рис. 4. Лилии тока (пунктирные лилии) и изосаты (сплошные лилии) при переносе высококонцентрированного рассола (Иа = 1000) и трассера (Иа = 0). Линии тока на (Ь) не показаны: они практически совпадают с представленными па (а)

Кроме того, требуется периодичность функций ш, р — ^ и функции ф — для ±-задачи либо функции ф + тХ для У-задачи. Сформулированные таким образом задачи решались конечно-разностными методами на равномерной по обоим направлениям сетке. Типичный шаг сетки составлял 1/64, что приводило к задачам, содержащим примерно 150000 узлов. Реалпзовывался итерационный метод с

р

и концентрационной задачи (46) при заданном ф. Каждая из задач решалась мно-госеточиым методом. Для подавления численной дисперсии при решении концентрационной задачи применялась специальная численная схема, обеспечивающая в конечном итого второй порядок аппроксимации коивективиых производных.

После решения задач на ячейке эффективные коэффициенты диффузии опре-

р

V = (|Ур|2) . (47)

Формула (40) для нахождения проницаемости преобразовывалась к виду

К!1 = (ш2) — Сг <(£ — р) Vг> .

Типичные результаты расчета полей плотности и скоростей движения жидкости для ||-задачи показаны та рис. 4. Для этой задачи прямая X = 0 (штрих-пунктирная линия на рис. 4) является осыо симметрии, поэтому без потери информации расчетные поля могут быть представлены в половине расчетной области. Важной особенностью рассматриваемого процесса является чисто канальная, без каких-либо вихревых зон, структура течения (рис. 4, а), что позволяет говорить о близости решетчатой модели к простейшим одномерным моделям. Близость подтверждается и структурой полей плотности. Как и в одномерном случае, с ростом И.а происходит их заметное выполаживание в ядре потока ортогонально силе тяжести. Гидродинамическое сходство объясняет практически полное совпадение коэффициента дисперсии трассера (рис. 5, а) для решетчатой и одномерных моделей. Различие в геометрии течения сказывается при вычислении коэффициента

10

а) Зависимость коэффициента дисперсии трассера от Ре

0.6

0.4

0.2

--- 10, роге

10, эМ

о 20, Ре = 10

□ 20, Ре = 20

о 20, Ре = 30

V 20, Ре = 100

200 400 600 800 ^

6) Зависимость /у (Яа)

Рис. 5. Сравнение характеристик для одномерных и решетчатой моделей

f

0

D

извилистости = 0.589 (в одномерном случае С = 1) и) чт0 более существенно, при нахождении поправочного множителя /ц. Для решетчатой модели, в отличие от одномерных, / является функцией обоих (Иа и Ре) параметров задачи (рис. 5, Ь). Однако в наиболее интересном с практической точки зрения случае

Ре /| Ре

и, в полном соответствии с результатами одномерного моделирования, поправочный коэффициент становится зависящим лишь от локального числа Релея. Более того, как можно видеть из рис. 5, Ь, зависимость /ц (Иа) является практически той же самой, что и в одномерном случае. Она может быть аппроксимирована формулой

1 + 3.4 • 10-6 Иа2

1\\ = -ТТ77- 48

11 1 + 4- Ю-3 11а + 10-711а11/4

Дадим пояснение, касающееся сравнения результатов, полученных для различных микромоделей пористой среды. Единственным параметром микромодели, входящим в управляющие комплексы Ре и Иа, является характерный размер / поры. Согласование моделей проводится за счет выбора / так, чтобы /ц (Иа) имел для всех моделей одно и то же асимптотическое представление при Иа ^ то. При выборе в качестве базового размера / радиуса поры в модели связки цилиндрических капилляров оказывается, что половина ширины щели в модели связки щелевидных капилляров и в решетчатой модели равны 1.2/ и / соответственно.

Обратимся к ±-задаче, определяющей поперечные коэффициенты диффузии и фильтрации. Ей можно придать наглядный физический смысл, разбив пористое пространство на совокупность горизонтально расположенных каналов, один из которых выделен на рис. 3. Течение жидкости осуществляется вдоль каналов: твердые стенки непроницаемы для примеси: на двух периодически чередующихся частях жидкой границы в силу симметрии задачи выполняются условия постоянства р. Можно принять, что при 2 = 0 (линия АС на рис. 3) р = 0, тогда на второй части жидкой границы р = л/2 (о + 1) в силу периодичности р — 2. Средняя плотность жидкости в таком канале равна, очевидно, р-£ = (о+ 1)/\/2. Линия 2 = 0 является линией симметрии, поэтому в нижнем соседнем канале средняя плотность есть —• Суммарный поток примеси между двумя каналами

па элементе периодичности определяется как

А

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

с

^ =/ Й ^ (49) А

на самом деле полностью эквивалентно (47), в чем можно убедиться, умножив (46) на р — рн и проинтегрировав результат по каналу периодичности.

В силу изотропности рассматриваемой модели ее поперечный коэффициент извилистости совпадает с продольным £ц; классический коэффициент фильтрации определяется равным 0.202. Коэффициент Х0 (Ре) трассера (рис. 5, а) монотонно растет с ростом Ре, выходя при Ре ^ то на предельное значение Х00 (то) « 16. Для объяснения такого поведения сравним между собой поля концентрации при И,а = 0 для умеренных и больших значений Ре. Они представлены на рис. 6, Ь, (1. По-прежнему изображается половина области течения вследствие симметрии ±-задачи относительно прямой Z = 0 (штрих-пунктирная линия на рис. 6).

При умеренных Ре (рис. 6, Ь) основное изменение р сосредоточено в узком диффузионном пограничном слое вблизи линии раздела двух потоков различной плотности (ось симметрии). Толщина диффузионного слоя пропорциональна Ре-1/2

Ре

Ре1/2. Это, однако, продолжается лишь до тех пор, пока соседние пограничные слои в узлах пересечения щелей не начинают взаимодействовать друг с другом посредством уходящих от них вниз по потоку концентрационных хвостов. Вследствие взаимодействия бывшие локальными изменения концентрации размазываются на весь поток (рис. 6, (I), рост Х0 (Ре) замедляется. В пределе при Ре ^ то толе р стабилизируется, а Х00 (Ре) выходит па постоянное значение.

Изучим изменение коэффициентов дисперсии и фильтрации с ростом безразмерного градиента плотности И.а. Из общих соображений [1] и качественного анализа лабораторных данных [10] следует ожидать уменьшения Х и К с увеличением И,а. Априори, однако, неясно, какие физические процессы лежат в основе этой закономерности, и лишь расчеты обнаруживают чрезвычайно интересный механизм ее реализации. Оказывается, что с ростом И.а происходит существенная перестройка структуры течения.

При небольших И,а течение — чисто канальное, а эффективные коэффициенты диффузии и фильтрации не отличаются от классических. Далее с ростом И.а возникает и укрупняется система возвратных течений. Динамика их развития показана на рис. 6, Ь / для умеренных (Ре = 10) и больших (Ре = 100) значений критерия Пекле. Указанная система отделяет друг от друга два основных потока существенно различной плотности (±рн). Чем шире разделяющая основные потоки зона возвратных течений, тем слабее диффузионный обмен между потоками, а

Рис. 6. Изосаты (а), (Ь) и линии тока (с)-(/) для ^-задачи. Рисунки (а), (с), (е) соответствуют умеренным. а (Ь), (г!), (/) большим значениям критерия Пекле

2.5

-2.5

в = 0 в = 10

в = 1000

-1

1 2

а) Ь)

Рис. 7. Профили р и |и| вдоль линии АВ (рис. 6, а); Ре = 10.

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

Динамику изменения с ростом И,а концентрации р и скорости движения жидкости V можно проследить и по рис. 7, па котором представлены профили изменения р и | VI вдоль о си X = 0 (линия АВ на рис. 6, а). Как видно из этих рисунков, при больших И,а основные потоки проходят в узкой зоне вблизи стенок канала. В пределе при И.а ^ то разность плотностей этих потоков равна 2р^. В остальной части сечения канала скорости течения жидкости малы, в силу чего основным механизмом переноса примеси является диффузия. Как следствие, при больших 11а профили плотности в этой части канала линейны, и градиент плотности равен др/дг = р^/\/2. Это позволяет, воспользовавшись (49), вычислить значение эффективного коэффициента диффузии при И.а ^ то:

Vо

1

2т

Обратим внимание, что в отличие от Х>ц, который обращался в нуль при И,а ^ то, поперечный коэффициент эффективной диффузии выходит па конечное значение . Поэтому поправочный та действие гравитации множитель /^ удобно определить как

- д - яг

Я

V0, - V0,

При таком определении с ростом И.а коэффициент /^ монотонно уменьшается от 1 до 0, а эффективный коэффициент дисперсии находится как средневзвешенное значение между предельными величинами

= а+я! /±+(1 - д).

Аналогичными свойствами обладает поправочный па действие гравитации множитель для коэффициента фильтрации /к = .

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

Р

0

□ Ре = 10

О Ре = 20 V Ре = 30 а Ре = 50

< Ре = 100

10-1 10° 101

102 О

0.8

0.6

0.4

0.2

0 10

□ Ре = 10

❖ Ре = 20

V Ре = 30

А Ре = 50

< Ре = 100

< А

а)

10 10

Ь)

102 О

Рис. 8. Зависимости fk(С;) (а) и f± (С¡) при фиксировашшх значениях Ре. Сплошные линии соответствуют аппрокснмациоппым формулам (50) и (51)

роль. Поэтому естественно ожидать, что в части гидродинамики основным управляющим параметром рассматриваемого процесса будет локальное гравитационное число , а те локальное число Релея И.а;, так в || -задаче. По этой причине удобно изображать поправочный коэффициент фильтрации Д как функцию от локального гравитационного числа С; при фиксированном значении Ре. Именно в таком виде он представлен на рис. 8. а. Как видно, в достаточно широком интервале изменения Ре, С; расчетные данные группируются вблизи общей кривой, которая может быть аппроксимирована выражением

1 + 0.08С;

fк =-(аО)

1 + 0.09С; + 0.0567 2

Для поправочного коэффициента дисперсии Д экспериментальные данные также оказывается возможным сгруппировать вблизи одной кривой, если предварительно нормировать локальное гравитационное число на эффективный коэффициент дисперсии трассера и изображать /0 как функцию от 6; = /Р0 . Соответствующие экспериментальные данные и результаты, полученные по аппрок-симационной формуле

1 + 0.1С, 1 + 0.25С?3/2

представлены на рис. 8, Ь. Заметим, что при больших значениях локального числа Пекле Р0 выходит та постоянную величину. Следовательно, при больших Ре параметр совпадавт с с точностью до постоянного множителя.

Формулы (50) и (51), совместно с формулой (48), и представляют, собственно, основной практический вывод проведенных исследований.

5. Сравнение с экспериментом

К сожалению, в настоящее время мы не располагаем достоверными данными по определению зависимости поперечных коэффициентов дисперсии и фильтра-

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

Тщательные эксперименты такого рода были проведены Н. Мобсг [2] и Б.Л. АпЛогбоп [11]. Их подробное описание, равно как и анализ результатов, можно найти в работах [5, 7]. Отсылая за подробностями к [5], скажем лишь, что экспериментальная установка представляла собой вертикально расположенную трубку, заполненную тщательно отсортированным грубозернистым (со средним диаметром зерна 0 = 0.73 мм) либо среднезернистым (0 = 0.27 мм) песком. В предварительно заполненную пресной водой систему через нижнее основание с постоянной скоростью нагнетался раствор МаС1 различной концентрации. Динамика размазывания концентрационного фронта наблюдалась в четырех портах, удаленных от входа на расстояние 10, 20, 30 и 40 см. Соответствующие динамические кривые являлись выходными данными каждого эксперимента. Нами были использованы все пригодные для анализа данные по крупнозернистому песку, представленные в [5]. В общей сложности это четыре эксперимента, проведенные при разности плотностей Др и средней скорости движения жидкости, равных, соответственно,

1. Др = 5 кг/м3, «о = 45 мм/мин.;

2. Др = 200 кг/м3, «о = 45 мм/мин.;

3. Др = 100 кг/м3, «о = 20 мм/мин.;

4. Др = 200 кг/м3, «о = 10 мм/мин.

Для численного моделирования эксперимента уравнение (38) дополнялось очевидными начальным и граничным условиями и решалось на равномерной сетке по явной схеме. Эффективный коэффициент диффузии принимался в виде (48): коэффициент извилистости £ц и коэффициент дисперсии трассера Х>о брались из

/| (Иа)

для согласования с 4 х 4 = 16 экспериментальными и теоретическими кривы/

Релея. Оказалось, что при выборе / = 0.18 мм согласие всех 16 экспериментальных н теоретических кривых становится вполне удовлетворительным. Качество совпадения иллюстрируется на рис. 9 для экспериментов 2 4. Данные по эксперименту 1 не приведены, так как этот случай близок к фильтрации трассера и совпадение результатов здесь гарантируется использованием экспериментальных значений Х>о •

/

усу поры. Можно ожидать, что отношение характерных радиусов зерна засыпки и поры примерно совпадает с отношением (1 — ш)/ш объемов, занимаемых твердой

/

0.64 мм, что разумно согласуется с экспериментальными данными (0 = 0.73 мм).

Аналогичным образом был обработан и единственный имеющийся в нашем распоряжении эксперимент для среднезернистого песка (Др = 100 к г/м3, «о = = 32.5 / /

шению размера зерна в эксперименте. И в этом случае достигалось удовлетворительное совпадение теоретических кривых с экспериментальными.

6. Обсуждение результатов

Макроскопические уравнения (29) и (38) были выведены методами теории гомогенизации для произвольных периодических пористых структур. Поэтому не

0

0

t - 'о (c)

4 - 'о (О

а)

Ь)

Рис. 9. Сравнение теоретических (сплошные лилии) и экспериментальных (пунктирные липли) кривых изменения концентрации со временем в датчиках, удаленных от точки входа на расстояние = 10 см (а) и г4 = 40 см (6). Сдвиг ¿0 = ^/г0 соответствует времени прохождения фронта через датчик

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

будут выполняться н для реальных пористых сред.

В лаборатории, как правило, реализуется случай больших Ре, поэтому величинами £ц, и обычно пренебрегают. 3начение К0 и зависимости Р°(Ре), Р0 (Ре) коэффициентов дисперсии трассера от скорости движения жидкости хорошо определяются в экспериментах. Основной практической проблемой является нахождение поправочных (на действие гравитации) множителей /ц, Д, /. Заманчиво предположить, что в факторизованных зависимостях (52) вся основная информация о структуре среды содержится в первых сомножителях, а поправочные коэффициенты /ц, Д, / универсальны для различных типов сред, что, в принципе, позволяло бы их вычислить на простых модельных средах. Именно это предположение было неявно сделано в предыдущем разделе статьи, когда при сравнении теоретических и экспериментальных результатов использовались экспериментальная зависимость Х>0(Ре) и теоретически обоснованная лишь для простейших моделей зависимость /ц(Ка) вида (48). Вполне удовлетворительное совпадение результатов вселяет определенную уверенность в истинности высказанного предположения, хотя конечно же не решает вопроса окончательно. Еще одним доводом в пользу данной гипотезы, казалось бы, является то, что зависимость /ц (И,а) одинакова (по крайней мере в некотором диапазоне изменения Ре) для одномерных и решетчатых пористых структур. Однако указанные структуры на самом деле однотипны. Это находит свое выражение, в частности, в одинаковом,

Р =С|| + Х>0(Ре)/||(аа, Ре), = с± + (Ре)/±(Иа, Ре) + (1 - Д (Ия, Ре)), К = К°Д(Ка, Ре)

(52)

тейлоровском, характере поведения коэффициента дисперсии V« ж

Ре2 реальных пористых засыпках наблюдается иной закон Х>о ж Реа (1 < а < 1-2, [13]).

Поэтому первоочередной задачей представляется проверка предположения об уни-

/|

тейлоровской зависимостью Х>о(Ре)-

Сказанное в полной мере относится и к поперечным коэффициентам , , /_!, /к. Особый интерес в данном случае представляет обнаруженный механизм самоорганизации внутрипоровых потоков под действием гравитационных сил. Есть основания полагать, что описанный в данной статье механизм понижения диспср-сивности и проницаемости за счет возникновения и роста вихревых структур является универсальным и должен проявляться, в той или иной мере, для любых типов пористых структур.

В заключение отметим, что учет гравитационных эффектов на уровне пор имеет смысл лишь в лабораторных экспериментах, где скорости фильтрации относительно велики, а ширины концентрационных фронтов малы. Для того чтобы понять это, достаточно оценить значение критерия Релея. При типичных значениях р = 10-3 кг/м с, йт = 10-9 м2/с для крупнозернистой среды с / = 10-4 м, при Др = 200 кг/м3 (предел насыщения для МаС1) и типичной для полевых условий ширине фронта Ь = 20 м [12] имеем Иа = 10-2, что практически соответствует случаю трассера (рис. 5, Ь). Другое дело, что в природных пористых средах появляется совершенно иной, на много порядков больший, чем пора, размер нсод-нородностсй: микронеоднородность поля проницаемости. Такие неоднородности, с одной стороны, генерируют на много порядков большую диспсрсивность, с другой стороны, за счет увеличения масштаба гравитационного воздействия обеспечивается ее более эффективное гравитационное подавление. Возникающая при этом задача макроскопического (уже не для пор, а для микронсоднородностсй) описания взаимосвязанных процессов разноплотностной фильтрации и массопсрсноса [13 15] близка к рассмотренной выше, однако ее обсуждение выходит за рамки данной статьи.

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

Работа выполнена при финансовой поддержке РФФИ (проект Л*1' 05-01-00516) и РФФИ-Х\УО (проект Л* 05-01-890001-Нв0).

Summary

D.E. Demidov, A.G. Egorov. Upscaling of density-driven flow and transport processes. 1. Porescale.

The macroscale equations describing coupled processes of density-driven flow and transport are thoroughly derived for the stable case (fresh water is above dense brine). Governing equations are Stokes equations and convection-diffusion mass transport equation, written at porescale. The main results are macroscale dependencies of dispersion and filtration coefficients on density gradient and filtration velocity.

Литература

1. Hassanizadeh S.M., Leijnse A, A non-linear theory of high-concentrat.ion-gradient. dispersion in porous media // Adv. Water Resour. 1995. V. 18, No 4. P. 203 215.

2. Moser Н. Einfluß der Salzkonzentration auf die hydrodynamische Dispersion im porösen Medium: Mitteilungen lieft. 128. Technische Universität Berlin. 1995.

3. Бахвалов H.C., Паиасеико Т.П. Осреднение процессов в периодических средах. Математические задачи механики композиционных материалов. М.: Наука, 1984. 352 с.

4. Беляев А.Ю. Усреднение в задачах теории фильтрации. М.: Наука, 2004. 200 с.

5. Watson S.J., Barry D.A., Schütting R.J., Hassanizadeh S.M. Validation of classical density-dependent, solute transport, theory for stable, high-concent.rat.ion-gradient. brine displacements in coarse and medium sands // Adv. Water Resour. 2002. V. 25, No 6. P. 611 635.

6. Bear J. Dynamics of fluids in porous media. N. Y.: Elsevier, 1972.

7. Schütting R.J., Moser H., Hassanizadeh S.M. High-concent.rat.ion-gradient. dispersion in porous media: experiments, analysis and approximations // Adv. Water Resour. 1999. V. 22, No 7. P. 665 680.

ß

schieden: Mitteilungen lieft. 60. Institut, für Wasserbau, 1985.

9. Thiele M. Gravity affected lateral dispersion and diffusion in a stationary horizontal porous medium shear flow // Transport, in Porous Media. 1997. V. 26. P. 183 204.

10. Leroy C., Hulin J.P., Lenormand R. Tracer dispersion in stratified porous media: Influence of transverse dispersion and gravity // J. Coiit.. Hydr. 1992. V. 11. P. 51 68.

11. Anderson S.J. An experimental investigation of high concentration displacements in saturated porous media: Pli.d. thesis. The university of Western Australia, 1985.

12. Sehelkes K., Vogel P. Paleoliydrological information as an important, tool for groundwater modelling of the Gorleben site // Paleoliydrological Methods and Their Applications for Radioactive Waste Disposal. Proc. OECD/NEA Workshop. Paris: OECD, 1992. P. 237 250.

13. Egorov A.G., Demidov D.E., Schütting R.J. On the interaction between gravity forces and dispersive brine fronts in micro-heterogeneous porous media // Adv. Water Resour. 2005. V. 28, No 1. P. 55 68.

14. Демидов Д.E., Егоров А.Г. Осреднение уравнений фильтрации рассолов в микропе-одпородпых пористых средах // Тр. Матем. центра им. Н.И. Лобачевского. Т. 16. Казань: Изд-во Казап. матем. о-ва, 2002. С. 163 174.

15. Демидов Д.Е., Егоров А.Г. Влияние гравитационных сил па дисперсионное размазывание фронтов при фильтрации рассола в микропеодпородпой пористой среде // Тр. Матем. центра им. Н.И. Лобачевского. Т. 27. Казань: Изд-во Казап. матем. о-ва, 2004. С. 107 114.

Поступила в редакцию 03.10.05

Демидов Денис Евгеньевич аспирант НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.

E-mail: Denis.DemiduvQksu.ru

Егоров Андрей Геннадьевич доктор физико-математических паук, старший научный сотрудник, заведующий кафедрой аэрогидромехапики Казанского государственного университета.

E-mail: Andrey.EgoruvOksu.ru