Содержание к диссертации
Введение
1 Макроскопическая модель фильтрации и массопереноса в однородных пористых средах 19
1.1 Постановка задачи 20
1.2 Процедура гомогенизации 23
1.3 Одномерные модели 32
1.4 Решетчатая модель 36
1.5 Модель идеального перемешивания 48
1.6 Сравнение с экспериментом 55
1.7 Обсуждение результатов 63
2 Макроскопическая модель фильтрации и массопереноса в микронеоднородных пористых средах 66
2.1 Постановка задачи и процедура гомогенизации 67
2.2 Многомерные слабонеоднородные среды 72
2.3 Сравнение с вычислительным экспериментом 82
2.4 Обсуждение результатов ., 89
3 Макроскопическая модель массопереноса в слоистых пористых средах 93
3.1 Постановка задачи и вывод макроскопических уравнений 94
3.2 Решение задачи на ячейке 97
3.3 Сравнение с вычислительным экспериментом 102
3.4 Фазовый переход в процессах массоперепоса при фильтрации разноплотностной жидкости 111
3.5 Обсуждение результатов 119
Список иллюстраций 124
Список таблиц
- Процедура гомогенизации
- Модель идеального перемешивания
- Многомерные слабонеоднородные среды
- Сравнение с вычислительным экспериментом
Введение к работе
Безопасное захоронение химических и радиоактивных отходов является сегодня ключевой экологической проблемой. По общему мнению наиболее приемлемым методом конечного захоронения высокорадиоактивных, генерирующих тепло отходов является их размещение в подземных геологических формациях. Такой способ имеет ряд преимуществ перед наземным складированием. Геологические образования, такие как соляные формации, существовали на протяжении многих миллионов лет. Они не терпели существенных изменений на протяжении длительных периодов времени. Таким образом, они позволяют изолировать высокорадиоактивные отходы (с крайне продолжительными периодами полураспада) на время, превышающее 100 000 лет. Более того, благодаря глубине залегания таких образований и малым скоростям фильтрации в них, радиоактивные изотопы при возможной утечке потеряют свою активность к тому времени, когда они достигнут биосферы.
Подземные каменно-соляные формации имеют ряд преимуществ в качестве возможных хранилищ радиоактивных отходов [65]. Во-первых, каменная соль практически непроницаема в силу своих пластических свойств и отсутствия трещиноватой пористости. Вследствие этого внутри каменно-соляных формаций практически отстутствует циркуляция подземных вод. Во-вторых, высокая теплопроводность позволяет рассеивать тепло, порождаемое высокорадиоактивными отходами.
В случае утечки из такого хранилища наиболее вероятным механиз-
мом доставки загрязнений к поверхности является перенос подземными водами [59,01]. Следовательно, для каждого потенциального места захоронения должна быть проведена оценка последствий возможного выброса изотопов [74]. Такая оценка должна учитывать тот факт, что в подземных водах вблизи соляных формаций наблюдаются высокие концентрации соли (вплоть до предела насыщения NaCl), соответствующие большим значениям плотности (до 1200 кг/м3) [34]. При этом возникают большие градиенты концентрации. Так, например, вблизи Горлебенского соляного купола в Германиии наблюдалось изменение плотности в 150 кг/м3 на расстояниях менее 30 м [21]. Наличие больших градиентов концентрации вызывает вопрос о применимости в рассматриваемой ситуации классических линейных реологических соотношений, лежащих в основе стандартных моделей фильтрационного переноса примеси: (г) закона Дарси, связывающего скорость фильтрации и градиент давления, и (гг) зависимости Фика между дисперсионным потоком примеси и градиентом концентрации.
В более общем виде проблема заключается в построении математических моделей фильтрации и массопереноса разноплотностных жидкостей, учитывающих наличие больших градиентов концентрации. Ее решение принципиальным образом должно зависеть от того, рассматривается ли неустойчивая (градиент плотности направлен против силы тяжести) либо устойчивая (градиент плотности сонаправлен силе тяжести) ситуация. В первом случае развитие неустойчивости приводит за счет процессов пальцеобразования к активному макроскопическому перемешиванию жидкостей; при этом концентрационные фронты размываются со скоростью, на много порядков превышающей предсказанную линейной теорией. Возможные подходы к моделированию таких процессов предложены в рабо-
тах [12,13,56,60]. Вторая ситуация реализуется вблизи соляных формаций и представляет основной интерес с точки зрения оценки экологических рисков создания подземных хранилищ радиоактивных отходов. Именно она и будет рассматриваться в данной работе. Ранние лабораторные эксперименты [46,66,68] качественно подтвердили тот физически ожидаемый факт, что увеличение градиентов концентрации при устойчивом расположении высококонцентрированного рассола и пресной воды приводит как к существенному снижению скорости размазывания концентрационных фронтов, так и к некоторому ухудшению фильтрационных свойств среды. В целом это дает отрицательный ответ на вопрос о возможности использования классических линейных реологических соотношений. Последующие тщательные лабораторные эксперименты [18,35,39,43,51] подтвердили это заключение. Авторы пришли к выводу, что наблюдаемые эффекты вызваны стабилизирующим влиянием силы тяжести и могут быть учтены зависимостью коэффициентов дисперсии и фильтрации от градиента концентрации. Вопрос о виде соответствующей зависимости, однако, до сих пор остается открытым. Наиболее популярная полуэмпирическая формула Hassanizadeh-Leijnse [28]
(l+/?|J|).7=-BVp,
связывающая между собой величину дисперсионного потока примеси J, градиент плотности V/? и заданный (классический) тензор дисперсии В, требует экспериментального определения в каждом конкретном случае параметра модели (3. В свою очередь, этот параметр оказывается [64,72] зависящим от скорости фильтрации и других параметров среды и процесса. Подчеркнем, что указанные лабораторные эксперименты имели дело с однородными пористыми средами. Подавление дисперсивности здесь
происходит уже на уровне пор. Аналогичные, но более масштабные эффекты имеют место и в реальных неоднородных пористых средах. Дисперсив-ность потока, вызванная неоднородностью среды, и здесь должна стабилизироваться гравитационными силами. Результаты, полученные в [76,77] на основе методов стохастической гидрогеологии, подтвердили это предположение. Для частного случая слабонеоднородных сред с экспоненциальной корреляционной функцией поля проницаемости были получены зависимости продольного [76] и поперечного [77] коэффициентов дисперсии от градиента концентрации. Полученные результаты, однако, нуждаются как в экспериментальной проверке, так и в обобщении на другие тины неодно-родностей.
Изложенное выше определяет цель диссертационной работы. Она заключается в выводе макроскопических уравнений, адекватно описывающих процессы фильтрационного переноса примеси при наличии больших градиентов плотности в устойчивой ситуации. Вообще говоря, выделяются два этапа такого исследования. На первом этапе необходимо провести осреднение базовых уравнений с масштаба нор до масштаба лабораторного эксперимента. Второй этап исследования предполагает осреднение уравнений фильтрации и массопереноса с масштаба лабораторных экспериментов (однородные пористые среды) до масштаба нолевых испытаний (микронеоднородные среды). На обоих этапах средством построения макроскопических уравнений служит метод осреднения дифференциальных уравнений с быстроосциллирующими коэффициентами (метод гомогенизации). Он был разработан в 70-х годах прошлого века в работах Н. С. Бахвалова, С. М. Козлова, В. В. Жикова, О. А. Олейник, J. L. Lions [2,11,17,20], и с тех
пор неоднократно применялся для анализа процессов, происходящих в пористых средах [4,15,38]. Процессы разноплотностпой фильтрации, однако, с этих позиций до сих пор не рассматривались.
Остановимся подробнее на содержании диссертации. Она состоит из введения, трех глав и списка использованных источников, содержит 137 страниц сквозной нумерации, в том числе 30 рисунков, 2 таблицы. Список использованной литературы насчитывает 77 наименований.
Первая глава работы посвящена первому из выделенных выше этапов исследования, а именно, осреднению с масштаба пор до масштаба лабораторного эксперимента базовых уравнений, описывающих фильтрационное движение и массоперенос на уровне пор.
Базовые уравнения в этом случае представляют собой связанную систему уравнений Стокса и уравнения конвективно-диффузионного массо-переноса. Постановка соответствующей задачи, анализ размерностей и схематизация исследуемых процессов приведены в параграфе 1.1. В результате анализа размерностей выделены управляющие параметры — число Релея и число Пекле — изучаемого процесса
\i&L а
Здесь д — ускорение свободного падения, Др — характерное изменение плотности, \i — вязкость жидкости, d — коэффициент молекулярной диффузии, г;о ~~ характерная скорость движения жидкости, величины I и L определяют соответственно характерные масштабы пор и процесса.
В параграфе 1.2 проводится процедура гомогенизации базовых (микроскопических) уравнений, состоящая в их асимптотическом анализе при
стремлении к нулю малого параметра є = l/L. Результатом этой процедуры является взаимосвязанная система нелинейных макроскопических уравнений фильтрации и массоиереноса. В отличие от классических моделей, коэффициенты этой системы — тензоры фильтрации и дисперсии оказываются функциями локальных чисел Пекле и Релея
a [id oz
Здесь через z обозначена вертикальная координата, а угловые скобки означают операцию осреднения. Показано, что в предельных случаях Ре, —> 0 и Ra, —» 0 полученная макроскопическая модель переходит в известные модели. Для вычисления зависимости тензоров фильтрации и дисперсии от Ре^, Raj в результате процедуры гомогенизации оказываются сформулированными так называемые задачи на ячейке. В следующих трех параграфах эти задачи конкретизируются и решаются для различных типов пористых сред.
В параграфе 1.3 рассмотрены простейшие одномерные модели пористой среды в виде связки одинаковых вертикально расположенных капилляров щелевидной или круглой формы. Эти модели позволяют рассчитать только продольные коэффициенты дисперсии Рц, но зато допускают простые аналитические решения. Результатом соответствующего анализа оказывается простое соотношение
2>ц (Ре,, Ra,) = d (С + V\ (Реї) /и (Ra,J J , Ra, = Ra, "^ . (0.1)
Здесь С — постоянный коэффициент извилистости, V?,(Pei) — коэффициент дисперсии трассера, вычисленный в пренебрежении плотностными эффектами, /||(Ra/) представляет собой поправочный на действие гравитации
множитель, /|j(0) = 1. С ростом Raj поправочный множитель монотонно стремится к нулю.
В параграфе 1.4 рассмотрена более реалистичная решетчатая модель пористой среды. Задачи на ячейке в этом случае двумерны; они решались численно с помощью специально разработанной вычислительной процедуры. Реализовывался итерационный метод расчета с последовательным решением гидродинамической (записанной в терминах функция тока — завихренность) подзадачи на заданном поле плотности и концентрационной подзадачи на заданном поле скоростей. Эта процедура позволяет рассчитать не только продольные, но и поперечные коэффициенты дисперсии Т>± и фильтрации К.±. Что касается поперечных коэффициентов, то из общих соображений и качественного анализа лабораторных данных следовало ожидать их уменьшения с ростом числа Релея. Априори, однако, было неясно, какие физические процессы лежат в основе этой закономерности и лишь расчеты обнаружили весьма специфичный механизм ее реализации. Оказалось, что с ростом Raj происходит существенная перестройка внут-рипорового течения с образованием и развитием зон возвратных течений. Именно они, с одной стороны, сужают живое сечение поровых каналов, а с другой стороны, препятствуют эффективному массообмену в жидкости. Результаты расчетов 2}j_ (Ре/, Ra/), JC± (Ре/, Ra/) в наиболее интересном диапазоне больших значений Ре/ оказалось возможным апроксимировать простыми зависимостями вида (0.1).
Аналогичные апроксимации продольного коэффициента дисперсии привели к результату, полностью аналогичному (0.1). Более того, как коэффициент дисперсии трассера V?, (Ре/), так и поправочный на действие гра-
витации множитель /||(Ra/) оказались практически теми же самыми, что и в одномерном случае. Это объясняется однотипностью этих моделей в части продольного массоперепоса. Однотипность находит свое отражение, в частности, в одинаковом, тейлоровском, поведении коэффициента дисперсии трассера при больших значениях критерия Пекле (^(Ре) ос Ре2). Известно, что в реальных пористых засыпках наблюдается иной закон роста коэффициента дисперсии трассера V2 ос Реа, 1 < а < 1.2. При этом возникает вопрос о степени универсальности найденной зависимость /ц(Иа/).
Для того, чтобы ответить на этот вопрос, в параграфе 1.5 рассматривается принципиально иная модель пористой среды в виде «цепочки бусинок». Эта модель отличается наличием зон идеального перемешивания в точках контакта бусинок и обладает «реальным» (Pj? ос Ре In Ре), а не тейлоровским видом зависимости коэффициента дисперсии трассера от средней скорости фильтрации. Численные расчеты, проведенные для этой модели, показали, что и в этом случае продольный коэффициент дисперсии может быть представлен в виде (0.1); поправочный множитель /ц при этом начинает существенно отличаться от вычисленного ранее лишь при /ц < 0.3. Последнее является следствием различного характера затухания /||(Ra) с ростом Ra^: /ц ос Ra; для одномерных и решетчатой моделей, и /ц ос Ra-2'3 для модели бусинок.
Для проверки адекватности построенных макроскопических моделей в параграфе 1.6 проводится сравнение полученных с их помощью результатов с результатами лабораторных экспериментов по одномерному замещению пресной воды концентрированным рассолом [72]. При этом для всех типов моделей коэффициент извилистости и коэффициент дисперсии
трассера брались из экспериментов. Единственным параметром для согласования 16 экспериментальных кривых с теоретическими был характерный размер I пор, участвующий в определении безразмерного градиента концентрации Ra/. Оказалось, что оба класса моделей, как «тейлоровские», так и модель идеального перемешивания, вполне удовлетворительно согласуются с экспериментальными данными. Это объясняется тем, что эксперименты были проведены в области малых и средних градиентов плотности, там, где обе модели дают практически одинаковые результаты. Окончательный выбор, очевидно, требует дополнительных экспериментов в области больших Ra/. В то же время отмечается, что дополнительным преимуществом модели идеального перемешивания является то, что она хорошо согласуется с экспериментами также и в части зависимости эффективного коэффициента дисперсии трассера от критерия Пекле.
Наконец, в параграфе 1.7 проводится обсуждение развиваемой теории и намечаются перспективы дальнейших исследований.
Процедура гомогенизации
Однако наибольший интерес, как отражающие свойства реальных пористых сред, представляют «симметричные» пористые структуры, остающиеся инвариантными относительно преобразований вида Z — —Z и X — —X. Для таких структур замена Z -Z, с--с, vl-vl Vl- -V\ тг- -т (1.24) оставляет систему уравнений (1.17), (1.21), (1.22) неизменной. Как следствие, 7Г = 0, а значит, для симметричных структур дополнительный фильтрационный поток (второе слагаемое в правой части (1.23)), индуцированный градиентом концентрации, невозможен. В пределе Ре/ — 0 такого типа среды независимо от значения Ra/ описываются на макроуровне классическим законом Дарси с JC± = K?L = l/7r и стандартным уравнением массоиереноса с постоянным и равным (ц. коэффициентом диффузии.
Аналогичные (1.24) замены показывают, что симметричные пористые структуры при любых Ra/ инвариантны относительно изменения направления течения жидкости на противоположное: К± (Ре/, Ra/) = /CJL (-Ре/, Ra/), V± (Ре/, Ra/) - V± (-Ре/, Ra/). Поэтому в дальнейшем можно ограничиться рассмотрением случая Ре/ 0.
Обратимся теперь к рассмотрению процесса одномерного вертикального движения с заданной скоростью VQ горизонтального фронта между пресной водой и высококонцентрированным рассолом. Процедура гомогенизации в этом случае проводится в подвижной, движущейся со скоростью Vo, системе медленных координат и приводит к следующим макроскопическим уравнениям: др и = 0, w = vo, - pg = 0, (1.25а)
Эффективный коэффициент диффузии 2?ц, как и раньше, является функцией чисел Пекле и Релея. Необходимо лишь строить число Пекле не по горизонтальной, а по вертикальной скорости VQ. Задача на ячейке (1.13) сохраняется, с заменой условий (vz) = 0, (vx) = 1 на (vz) = 1, (vx) — 0, а Х подсчитывается по формуле І):=1 + РС/( )"РЄ?(( "1)С Для симметричных пористых структур ц, как и прежде, не зависит от знака Ре/, т. е. направления (вверх, вниз) движения фронта. Аналогично предыдущему определяются также коэффициент извилистости ц и коэффициент дисперсии трассера ТЯ(Ре/). В общем случае они не совпадают с Сі и Dj_(Pei). Лишь в дополнительном предположении об «изотропности» (инвариантности относительно преобразования (X, У) — (У, X)) микроструктуры ц = (j_. Но даже в этом случае V±(Pci) DJ?(Pej).
В заключение данного параграфа укажем полезные в дальнейшем альтернативные формы записи эффективных коэффициентов диффузии и фильтрации. Для этого умножим (1.13с) нас и проинтегрируем результат по ячейке периодичности. С учетом того, что (cv-c) = \(v-V,?)=l-{V(vcf) = 0 в силу несжимаемости поля скоростей, условий прилипания и периодичности, получим, как в продольном, так и в поперечном случаях V = (\V(VelC + Z)\2Y (1.26)
Аналогичные выкладки после скалярного умножения (1.13b) на v и интегрирования по ячейке периодичности дают /CI1 = (V 2 + \Vvz\2) - 11 (cvz). (1.27)
Конкретизация зависимостей JC±, Т ±, Т \\ от Ре/, Raj предполагает задание микроструктуры пористой среды. Рассмотрим в качестве первого примера пористую среду в виде связки одинаковых, вертикально расположенных капилляров — щелей раскрытия 21. В этом простейшем случае, очевидно, имеет смысл говорить лишь об определении продольного эффективного коэффициента диффузии Т)
Модель идеального перемешивания
Рассмотренные выше модели пористой среды, на самом деле, однотипны (действительно, решетчатую модель в силу ее симметричности можно представить в виде системы изогнутых поровых каналов). Это находит свое выражение, в частности, в одинаковом, тейлоровском, характере поведения (Т 9, ос Ре2) коэффициента дисперсии при больших значениях Ре. В реальных пористых засыпках, для интересных с практической точки зрения значений Ре = 101 -flO2, наблюдается иной закон роста коэффициента дисперсии [72]: ТІЇ ос Реа (1 а 1.2). Различие в поведении jj(Pe), однако, еще не означает, что вся информация, полученная на основе анализа одномерных и решетчатой модели, представляет лишь теоретический интерес. В первую очередь это относится к выявленному выше характеру зависимости поправочных коэффициентов /ц, /j_, fk от гравитационных сил. Возникает вопрос о степени универсальности зависимостей, справедливых для «тейлоровских» модельных сред. Естественно адресовать этот вопрос в первую очередь таким моделям, которые демонстрируют «реальный», а не тейлоровский закон изменения коэффициента дисперсии в зависимости от критерия Пекле.
В данной главе мы ограничим соответствующий анализ продольным поправочным коэффициентом /ц, вычислив его для модели «цепочки бусинок», представленной на рис. 1.9 и имеющей близкую к наблюдаемой в экспериментах зависимость ){?(Ре). Особенностью этой модели является наличие зон идеального перемешивания жидкости в точках контакта бусинок — кругов заданного радиуса I. Элементом периодичности здесь
Задача (1.40)-(1.42) требует для своего решения применения численных методов. Алгоритм численного решения практически совпадает с использованным ранее при решении задач на ячейке для решетчатой модели и состоит в последовательном решении гидродинамической и концентрационной подзадач. Единственное отличие заключается в выборе сетки. При ее построении расчетная область (круг) отображалась на прямолинейный канал с помощью конформного отображения [14]
а в качестве узлов сетки выбирались точки пересечения равномерно распределенных линий уровня вещественной и мнимой частей отображения (рис. 1.9). Преимуществом такой сетки является то, что линии тока всюду близки к линиям = Const. Данное обстоятельство позволило упростить апроксимационную схему для конвективных членов в (1.40с). В поперечном () направлении они апроксимировались центральными разностями. Точность апроксимации в продольном направлении (rj) не являлась критичной из-за отсутствия в этом направлении пограничного слоя. Поэтому продольный конвективный член в (1.40с) аироксимировался направленными разностями.
Модель идеального перемешивания дает в части эффективного коэффициента дисперсии трассера результаты, близкие к экспериментальным. Коэффициент дисперсии V? трассера монотонно возрастает с ростом Ре, выходя при Ре — со на асимптотический режим ТУ» ос Peln(Pe). Во всем диапазоне изменения критерия ПеклеТ 9, может быть апроксимирован формулой
Многомерные слабонеоднородные среды
При заданных средних скоростях (U ) и заданном макроскопическом градиенте дС /дуі 0 плотности она имеет единственное решение, од нозначно определяющее функции Л)(У), С (У), P 2\Y) И константу П через параметры задачи на ячейке. Далее, определив можно по формулам (2.8) найти Q и П как функции дС /дуі, замкнув тем самым систему осредненных уравнений (2.7).
В общем случае нелинейная система уравнений (2.9) сложна для сколь нибудь исчерпывающего анализа. Существует, однако, достаточно широкий класс слабонеоднородных сред, для которых задача на ячейке линеаризуется и поэтому допускает достаточно полное исследование.
Под слабонеоднородными понимаются среды, проницаемость К(У) которых может быть представлена в виде K{Y) = 1 + 51C, ()=0, (JC2) = 1.
Малый параметр 5 характеризует среднеквадратичное отклонение проницаемости от своего среднего значения. Введенное ранее характерное значение &о проницаемости в этом случае есть ни что иное как средняя про ницаемость. Малость 5 позволяет отыскивать решение задачи на ячейке в виде 0W = V + 6u, PW = 5V, С = 6Рес -, (2.10) ду где для простоты записи через V обозначено среднее значение скорости (СЛЛ. Будем считать V = 1.
Этого всегда можно добиться, задав характерное значение qo как модуль средней скорости фильтрации. Конкретизируя еще одно характерное значение d = dlL(qo), запишем (с точностью 0(5)) локальный тензор дисперсии в виде 4Ы - dl±(q0) Dij = 6ij + dViVj, d = . , d±{Qo) и, удерживая в (2.9) главные по 5 члены, получим следующую линейную задачу на ячейке для определения й, V, с: V-u = 0, (2.11а) и = KV - VV + Raceb (2.11b) PeV Vc - Ac - dViVj— - + «i = 0. (2.11c)
Дополнительные условия (и) = (с) = 0 здесь выполняются автоматически, в чем можно убедиться, взяв среднее от (2.lib), (2.11с)Будем считать ее в дальнейшем зависящей лишь от \Y — Y"\, что отвечает рассмотрению макроскопически изотропных сред. Фигурирующее в (2.12) среднее (щс) может быть выражено через 1Z в виде (ию) = — JJn (У - Y"\) Фс (Г) ФХ (Y") dY dY" = = -Л72 /// (И) e Y -Y 4c (Г) Фс (Y") dY dY"du = = - 1п(\и\)Фс(иЩ(и;) du. (2.14) Здесь звездочка означает комплексное сопряжение, а через it обозначена трансформанта Фурье корреляционной функции, так называемая спектральная функция. Аналогичные выкладки для среднего {К,й) дают {Кит) = —i [K{\u \)&m(u)dLj, т = 1,2,3. (2.15) (27Г) J
Таким образом, задача нахождения эффективных характеристик среды свелась к задаче вычисления интегралов (2.14), (2.15) в пространстве изображений. Далее ограничимся тем, что представим результаты вычислений для V\ = ±1, V2 = V3 = 0, что отвечает вертикальному макроскопически однородному фильтрационному потоку. В этом случае интеграл (2.14) принимает вид (ию) = -— Jii(\Lj\) pc(u) du, (2.16) ..лл_ ( + 4)( -offa Ve2u\u + [о;2 (w2 + du\) + Ra (a;2 - a;2)]2 Важно отметить, что практический интерес имеют лишь большие значения критерия Ре. В этом случае для исследования интеграла (2.16) используется стандартный подход. Он опирается на то, что основной вклад в (2.16) вносит малая окрестность точки ш\ = 0. Используя этот факт, оказывается удобным перейти в (2.16) от интегрирования по переменным о; к интегрированию по переменным 5 = РеО 1, Г = л/U)\ + 6 3 Считая далее s, г величинами порядка единицы, и преобразуя с точностью 0(Ре-2) подынтегральные функции к виду .
Сравнение с вычислительным экспериментом
Для проверки адекватности полученной макроскопической модели проводились специальные вычислительные эксперименты, основанные на прямом решении системы уравнений (3.2). Расчетная область, с учетом симметричности поля проницаемости, представляла собой полосу 0 X 7Г, -со Z со, на боковых сторонах которой задавались условия симметрии дС Х = 0,тг: Ux = 0, - 7 = 0.
Начальные условия отвечали острому фронту между пресной водой и концентрированным раствором при устойчивом их расположении. При проведении экспериментов помимо гармонической (JC(X) = cosX) исследовались и другие, в частности, кусочно-постоянная, типы неод пород ностей. Все расчеты проводились при одном и том же значении Ре = 200/7г; значения d, д, G варьировались в широких пределах.
В типичных расчетах использовалась неоднородная прямоугольная сетка, состоящая из 129 х 2049 узлов. Высота Z сеточной области \Z\ Zoo увеличивалась по мере размазывания концентрационного фронта с тем, чтобы выдержать с заданной погрешностью е = 10_6 условия C\z=Zoo = 0, C\z=_z \ (3.20) на верхней и нижней границах расчетной области. На каждом временном шаге расчет состоял в последовательном решении двух подзадач. Вначале на заданном с предыдущего временного слоя поле плотности решалась гидродинамическая задача (3.2а), (3.2b) нахождения поля скоростей. Использовался известный метод [16,54] быстрого преобразования Фурье вдоль вертикальной координаты с последующим обращением методом прогонки серии трехдиагональных матриц. Далее, найденное поле скоростей использовалось для пересчета на текущем временном слое поля плотности. Пересчет основывался на решении уравнения (3.2с) конвективно-диффузионного мас-сопереноса. Использовался метод расщепления по физическим процессам.
Конвективная часть задачи решалась методом ENO третьего порядка [36]. Диффузионную часть задачи, в силу слабости диффузионного механизма (малость Ре- ), оказывалось возможным решать по явной схеме. Для сравнения с ММ результаты каждого эксперимента осреднялись вдоль горизонтальной оси X и представлялись в виде 7Г (С) {Z,T) = - fc(X,Z,r) dX. о Типичные результаты расчетов, включающие микроскопические поля скоростей фильтрации и плотности, а также осредненный профиль плотности (С), представлены на рис. 3.2.
Макроскопическая нелинейная система уравнений (3.17), (3.18), дополненная очевидными начальным и граничным условиями, решалась по чисто неявной схеме методом Ньютона. Шаг по времени в начале расчетов принимался равным Ат = 0.1 (AZ) , затем, по мере размазывания фронта, динамически изменялся так, чтобы метод Ньютона сходился за 2 -т- 8 итераций. Длина расчетной области, так же, как и в вычислительном эксперименте, увеличивалась по мере размазывания концентрационного фронта с тем, чтобы с заданной погрешностью є выдержать условие, аналогичное (3.20).
Перейдем к сравнению результатов вычислительных экспериментов и макроскопической модели. Прежде всего, отметим, что задача на ячейке (3.7) инвариантна относительно направления Z. Следовательно, во-первых, ММ не зависит от направления (вдоль или против силы тяжести) средней скорости движения жидкости, и, во-вторых, изначально симметричные профили плотности остаются таковыми и далее.
Следовательно, профили Ra и J для обратных направлений потока являются зеркальным отображением друг друга по отношению к Z = 0. Сами профили при этом асимметричны (чем больше число гравитации, тем сильнее эта асимметрия). Если бы эта асимметрия была значительной, ММ (3.17), (3.18) не могла бы дать адекватного описания рассматриваемой задачи и для достижения желаемой точности необходимо было бы решать вместо (3.7) задачу на ячейке второго (или более высокого) порядка по S. Однако, при реальных значениях гравитационного числа (G 5) асимметрия оказывается незначительной, а профили осредненной плотности (С), безразмерного градиента плотности Ra и среднего потока примеси J, полученные для противоположных направлений скорости фильтрации, близкими друг к другу. Соответствующие профили ММ лежат между экспериментальными, полученными для противоположных направлений движения фронта. Типичная ситуация представлена на рис. 3.3(a) Дополнительно на рис. 3.3(b) представлена на тот же момент времени фазовая плоскость. Сплошная и пунктирная линии на этом рисунке соответствуют множеству пар {(Ra(Z), J{Z)) — оо Z оо} для экспериментальных и теоретических профилей плотности.
