Содержание к диссертации
Введение
1 Определение гидравлических функций пористой среды 18
1.1 Модели невзаимодействующих капилляров 18
1.1.1 Модель связки цилиндрических капилляров 18
1.1.2 Модель Туллера-Ора 20
1.2 Модель взаимодействующих капилляров 25
1.2.1 Решётка 25
1.2.2 Капилляры 27
1.3 Гидравлические функции при дренаже 31
1.3.1 Капиллярная кривая 31
1.3.2 Относительная фазовая проницаемость 36
1.3.3 Верификация модели по экспериментальным данным для дренажа 40
1.4 Гидравлические функции при пропитке 45
1.4.1 Капиллярная кривая 45
1.4.2 Относительная фазовая проницаемость 52
1.5 Кривая первичной пропитки 55
2 Теоретический анализ релаксационной модели влагопереноса 61
2.1 Устойчивость уравнений влагопереноса: общие результаты . 62
2.1.1 Модель Ричардса 62
2.1.2 Анализ устойчивости модифицированной модели Ричардса 66
2.2 Релаксационная модификация модели Ричардса 68
2.2.1 Термодинамический анализ 68
2.2.2 Модели Р- и 5-релаксации 70
2.2.3 Учёт гистерезиса в релаксационной модели 71
2.3 Анализ релаксационных моделей на решениях типа бегущей волны 74
2.3.1 Модель Р-релаксации 74
2.3.2 Модель 5-релаксации 79
2.4 Анализ устойчивости модели Р-релаксации 81
2.5 Двумерное моделирование процесса пальцеобразования 86
3 Верификация релаксационной модели Ричардса 91
3.1 Эксперименты D. A. DiCarlo 92
3.2 Модель 97
3.3 Определяющие уравнения 99
3.4 Конкретизация гидравлических функций К и Р 100
3.5 Решение типа бегущей волны 102
3.6 Определение коэффициента релаксации 104
3.7 Результаты расчётов 106
Заключение ПО
Литература 112
- Модель связки цилиндрических капилляров
- Гидравлические функции при пропитке
- Анализ устойчивости модифицированной модели Ричардса
- Определяющие уравнения
Введение к работе
Угроза поверхностного заражения подземных вод токсичными веществами при авариях на промышленных предприятиях и складах химических продуктов является в наше время актуальной экологической проблемой. Попавшие на дневную поверхность водорастворимые загрязнения мигрируют с потоком влаги через зону аэрации к зеркалу грунтовых вод и могут в конечном итоге нанести непоправимый ущерб природным экосистемам обширного региона. Защитные свойства зоны аэрации, связанные с адсорбцией загрязнителя частицами почвы, во многом определяются режимом фильтрационных потоков.
Многочисленные эксперименты [55, 45, 83, 64] свидетельствуют о том, что в ненасыщенных грунтах движимый гравитацией однородный фронт пропитки, как правило, распадается на устойчиво развивающуюся систему потоков («пальцев»). Понимание причин, вызывающих неустойчивость таких фронтов, и учёт этого эффекта необходимы для верного предсказания интенсивности переноса влаги и водорастворимых загрязнений от дневной поверхности к зеркалу грунтовых вод. Действительно, наличие предпочтительных путей, порождаемых распадом фронта пропитки на отдельные пальцы, существенно уменьшает время миграции влаги в зоне аэрации, снижая тем самым её защитную роль.
Практическая важность и научная значимость феномена пальцеобразо-вания диктуют необходимость разработки адекватной математической модели этого явления. К сожалению, традиционные модели влагопереиоса (например модель Ричардса) оказываются в данной ситуации малопригодными, так как не принимают в расчёт динамических эффектов, играющих значительную роль при образовании и развитии пальцев. Многочисленные попытки релаксационной модификации модели Ричардса дали определённое понимание некоторых аспектов проблемы, однако построение целостной модели остаётся всё ещё делом будущего.
Одной из характерных особенностей пальцев, отмечаемой всеми экспериментаторами, является немонотонность профилей давления и насыщенности в их поперечном сечении. Эта немонотонность приводит к тому, что
Введение
насыщенность в теле пальца оказывается значительно выше, чем насыщенность вне его. В то же время поток влаги, определяемый градиентом давления, оказывается направлен к оси пальца. В этой ситуации высоконасыщенный палец не отдает влагу окружающей пористой среде, а, напротив, иссушает её.
На самом деле возможность движения влаги в пористой среде в направлении градиента влажности (т. е. из областей с малой влажностью — в области с большой) была впервые обнаружена не в экспериментах по паль-цеобразованию, хотя и приблизительно в то же время. По-видимому первые наблюдения этого явления были сделаны в работах Абрамовой [1,2] и Аллера [47, 48, 49]. Позднее опыты Абрамовой по сути были повторены Collis-George et al. [35]; аналогичные явления наблюдались в экспериментах Дмитриева [14] и Бондаренко [5], а также были воспроизведены Роде и Романовой [24]. Именно эти эксперименты впервые поставили под сомнение применимость модели Ричардса для описания переноса влаги в ненасыщенной пористой среде.
Опишем вкратце суть поставленных экспериментов. Во всех них исследователи изучали движение «подвешенной» (по терминологии А. Ф. Лебедева [19]) влаги в вертикальной колонке, заполненной образцом природного грунта. Грунт в колонке промачивали водой, после чего с её верхней поверхности было организовано испарение (как при помощи нагрева электрической лампой, так и естественным путём без нагревания при комнатной температуре). Давление снималось при помощи ряда установленных с определённым шагом тензиометров, влажность определялась по исследованию взятых с разной глубины образцов почвы.
Во всех экспериментах исследователи отмечали восходящее движение влаги из глубины почвы к поверхности испарения, которое, в особенности на первых стадиях процесса, совершалось при отсутствии градиента влажности и иногда даже в направлении возрастающей влажности.
Описанные эксперименты не могли быть объяснены в рамках традиционно используемой в почвоведении диффузионной теории. Согласно ей, перенос влаги в ненасыщенной пористой среде описывается уравнением Ричардса [77, 78]:
ds О ( ds\
Здесь s — влажность, t — время, х — пространственная координата, D(s) — коэффициент диффузивности, D(s) = Кдф/ds, где К — коэффициент вла-гопроводности, ip = i/)(s) — капиллярный потенциал влажности.
Кроме натурных экспериментов, о неприменимости диффузной теории по крайней мере для нестационарных процессов свидетельствуют и резуль-
Введение
таты численного моделирования. В качестве примера можно сослаться на работу Бондаренко и др. [6], в которой исследовалась степень отклонения численно рассчитанных по уравнению (1) значений влажности s от экспериментально наблюдаемых значений на различных стадиях сушки. Из сравнения этих величин видно, что расчётные значения даже качественно не совпадают с экспериментом на первых стадиях процесса.
Аллер был также первым, кто попытался теоретически обосновать описанный феномен [49, 50]. Для этого он модифицировал уравнение (1), заменив в нём капиллярный потенциал ф на так называемый эффективный потенциал влажности фе. Этот эффективный потенциал конструируется в виде суммы
фе = ф + А—,
где А — некий варьируемый коэффициент. После этой подстановки уравнение (1) принимает вид
Os д /_. . <9s л d2s \
Предложение Аллера вызывает два естественных вопроса: а) воспроизводит ли указанная модификация экспериментально наблюдаемые факты и б) насколько данная модификация является теоретически обоснованной.
Что касается первого вопроса, то ответ на него можно считать положительным. Численные эксперименты [22, 26], проведённые по модели Аллера для значений А = 0 и А = 100 показали, что при А = 100 происходит поток влаги вдоль градиента влажности, в то время как при А = 0 этот эффект не наблюдается. В другой работе [7] численные расчёты по модели Аллера сравнивались с экспериментальными профилями влажности на различных стадиях сушки. При А — const экспериментальные и рассчитанные кривые демонстрируют неплохое согласование на начальных стадиях процесса. При этом было отмечено, что величина А не остаётся постоянной в процессе сушки, хотя никаких предположений о виде зависимости А от времени не было сделано.
Несмотря на всё сказанное, против модели Аллера были выдвинуты многочисленные возражения. Критики указывали на то, что уравнение (2) и рассуждения, приводящие к нему, аналогичны таким, которые обычно используются для описания движения жидкости в трещиновато-пористых средах. Однако обычная почва имеет существенно более однородную структуру, и в ней едва ли можно выделить крупные транспортные каналы и мелкие питающие капилляры, как это сделал Аллер.
По мнению Роде [25], теоретическое обоснование, выдвигаемое Адлером, является чисто умозрительным. Роде считает, что Аллер попросту не
Введение
принял в расчёт гистерезис зависимости влажности от давления. Между тем отчётливо показано [66, 99, 15], что после того, как инфильтрация сменяется испарением, распределение давления вследствие влияния гистерезиса должно сразу же измениться. Это создаёт предпосылку к передвижению влаги к испаряющей поверхности не только без появления соответствующего градиента влажности, но даже в направлении этого градиента, хотя и против градиента давления. Имеются теоретические расчёты [17], показывающие, что обычная модель Ричардса с гистерезисом ip{s) даст движение воды в направлении градиента влажности.
Среди других подходов к модификации модели Ричардса нельзя не упомянуть работы А. В. Лыкова [20, 21]. Лыков исходит из совершенно других соображений, чем Аллер, хотя его уравнение также предназначено для более полного учёта нестационарности. Вывод Аллера базируется на рассмотрении капиллярной модели, тогда как Лыков получает своё уравнение на основе методов термодинамики необратимых процессов. Он рассматривает релаксационные процессы в бесконечно малом элементе среды и записывает уравнение движения жидкости в пористой среде в виде:
ds л d2s 0 ґп, ,ds\
Кулик [18] исследовал уравнения Аллера и Лыкова с точки зрения инвариантности относительно непрерывных групп преобразований. Его анализ показывает, что ни то, ни другое уравнение нельзя принять в качестве универсального уравнения движения влаги в почве. Модель Аллера, созданная специально для объяснения экспериментов по оттоку влаги при испарении, оказывется неспособна описать приток влаги в почву при пропитке. В то же время модель Лыкова, исходящая из экспериментов по инфильтрации жидкости в горизонтальную колонку, входит в противоречие с другой группой экспериментов. Кулик видит выход из этой ситуации в привнесении в модель гистерезисных свойств. По его мнению для того, чтобы быть равно пригодной как для пропитки, так и для дренажа, в модели Аллера коэффициент А должен испытывать очень сильный гистерезис. Коэффициент Лыкова AL должен либо быть очень малым, либо также быть подвержен сильному гистерезису. Кулик также предлагает своеобразный гибрид уравнений (2) и (3) в следующем виде:
Os я d2s д / ,ds\ л д2 fds\ д ,„ . где коэффициенты А и А] должны испытывать очень сильный гистерезис.
Введение
Возрождение интереса к модификации уравнения Ричардса было связано с появлением большого количества экспериментальных фактов, касающихся экспериментов по пальцеобразованию. Как и в случае с ранними экспериментами по движению влаги в природных грунтах, описать этот феномен в рамках традиционной модели не оказалось возможным. Интересно отметить, что и в этом случае при построении модели весьма продуктивным оказалось привлечение термодинамических соображений, а также учёт в явном виде гистерезиса гидравлических функций. Наибольший вклад в теорию нестационарных течений в ненасыщенных пористых средах в настоящее время принадлежит S. Majid Hassanizadeh. В серии работ [51, 52, 53, 54] он последовательно применяет методы неравновесной термодинамики к процессам перераспределения влаги в поровом пространстве. Записывая балансовые соотношения для свободной энергии Гельм-гольца на границе раздела фаз, Hassanizadeh получает необходимое ограничение, накладываемое на возможные релаксационные модификации модели Ричардса в виде энтропийного неравенства
Os SL(p-P(s))^0.
Им также была предложена одна из таких модификаций, известная как модель Р-релаксации:
ds
T—=p-P{s). (4)
Тщательный математический анализ этой модели был выполнен Cuesta et al. [36]. Предполагая, что зависимости P(s) и K(s) имеют степенной характер, ими было доказано существование и единственность решения типа бегущей волны уравнений (1), (4). Также был получен критерий немонотонности таких решений и исследовано их поведение при стремлении начальной водонасыщенности среды к нулю.
Привнесение гистерезисных свойств в предложенную Hassanizadeh модель было выполнено в работе Беляева и Hassanizadeh [30]. Авторы использовали простейшую модель гистерезиса (так называемую "play-type"), в которой все сканирующие кривые представляли собой вертикальные отрезки s = const.
Анализ устойчивости решения типа бегущей волны уравнения Ричардса, а также некоторых её модификаций (релаксационной и стефановского типа) был проведён Егоровым и Даутовым [13]. Они доказали, что модель Ричардса в своём исходном виде абсолютно устойчива, а её стефановская модификация — абсолютно неустойчива и, следовательно, оба этих случая не подходят для моделирования пальцев. В отличие от них, релаксационная модификация демонстрирует условную неустойчивость, т. е. устойчи-
Введение
вость на высокочастотных возмущениях и неустойчивость на низкочастотных. Это приводит к наличию минимума на дисперсионной кривой, положение которого определяет характерный масштаб системы пальцев, генерируемых изначально однородным фронтом пропитки. Эти же авторы предложили иную возможность релаксационной модификации модели вла-гопереноса — модель 5-релаксации. Вместо зависимости (4) в этом случае используется релаксационный закон вида
ds
Tdt = S(P) ~ S- (5)
Разумеется, всякая математическая модель должна быть тщательно протестирована и верифицирована по результатам натурных экспериментов. Только после этого можно делать окончательное заключение о её достоверности и практической применимости. С другой стороны, даже самый скрупулезный и аккуратный лабораторный эксперимент немного стоит без надёжного теоретического фундамента и непротиворечивых моделей. Накопленный массив экспериментальных данных с одной стороны и достижения теоретиков с другой позволяют в настоящее время значительно приблизиться к построению законченной модели влагопереноса в сухих средах.
Изложенное выше определяет основную цель данной работы, заключающуюся в построении и верификации математической модели влагопереноса в ненасыщенном грунте, способной адекватно описать процессы образования и развития пальцев.
Остановимся подробнее на содержании диссертации. Она состоит из введения, трёх глав, заключения, 61 рисунка, 7 таблиц. Список использованной литературы содержит 99 наименований.
Первая глава работы посвящена исследованию особенностей процессов перераспределения влаги в пористой среде на основе микромоделирования.
Первоначально исследуются процессы перераспределения влаги в масштабе пор. В параграфе 1.1 рассматривается модель пористой среды в виде связки цилиндрических капилляров различного размера, не взаимодействующих между собой. Подробно излагаются положения микромодели Туллера-Ора, в которой эти капилляры имеют полигональное поперечное сечение. Подход Туллера-Ора является более реалистичным по сравнению с моделью цилиндрических капилляров, поскольку позволяет учитывать плёночные и уголковые течения жидкости, доминирующие при малых насыщен-ностях пористой среды. Тем не менее игнорирование взаимодействия пор является его принципиальным недостатком.
Параграф 1.2 представляет модель, преодолевающую этот недостаток путём моделирования пористой среды в виде решётки со связями — цилиндрическими капиллярами различного поперечного сечения. Определя-
Введение
ющими параметрами для такой модели являются координационное число решётки Z, равное среднему числу выходящих из одного узла связей, и плотность распределения поперечных размеров связей (капилляров). Экспериментально установлено, что для типичных лабораторных песков, с которыми обычно имеют дело исследователи, это число в среднем равно 4.6. В качестве функции распределения пор по размерам использовалось гамма-распределение
fn{L) = [n{L I х, L0) = ,,.,., 77 exp(-L/L0),
І*+Г(х+1)
где Lq, x — параметры функции распределения. Выбор для нашей модели двухпараметрического закона распределения позволяет согласовать её с известными эмпирическими моделями Брукса-Кери и Ван-Генухтена, в которые также входят два определяющих параметра. Кроме того, дополнительным параметром модели, повышающим её адекватность реальным пористым средам является доля капилляров различного поперечного сечения. Для предлагаемой модели в качестве таких капилляров были выбраны криволинейный треугольник и круг.
Сформулированная выше микромодель была использована в параграфе 1.3 для вычисления гидравлических функций (капиллярного давления и гидравлической проницаемости) пористой среды. Эти функции определялись осреднением насыщенностеи отдельных капилляров по всему массиву пор с учётом взаимодействия между ними. Зависимость водонасыщенности среды от давления s(p) может быть представлена в виде суммы вкладов жидкости, сосредоточенной в заполненных капиллярах, уголковых менисках пор и в плёнках на твёрдой поверхности: s(p) = scap(p)4-smcn(p)-fSfiim(p). Асимптотический анализ поведения этих функций при р —> со показал, что влияние адсорбированных пленок начинает сказываться лишь при очень низких насыщенностях s ~ 10~4, и в типичных ситуациях им можно пренебречь. В то же время мениски в углах пор играют существенную роль при влажностях менее 20% и их влияние необходимо учитывать.
Другой важной гидравлической характеристикой пористой среды является зависимость относительной фазовой проницаемости К от насыщенности. Для её вычисления была применена известная электрическая аналогия, сводящая проблему отыскания фазовой проницаемости к задаче нахождения эффективной проводимости решётки по заданному распределению про-водимостей её связей. Считая, что проводимости отдельных пор распределены независимо друг от друга и используя хорошо известную аппроксимацию среднего поля, эффективная проводимость решётки К* (р) может быть
Введение
найдена как решение алгебраического уравнения
1- I k
Z~ \k + K,{Z/2- 1)
Здесь k — проводимость отдельной связи, а под средним понимается математическое ожидание относительно заданной плотности fn(k), которая восстанавливается по известной плотности /„(L) распределения пор по размерам. Полученная в результате решения зависимость относительной фазовой проницаемости К от давления р вместе с найденной ранее зависимостью водонасыщенности среды от давления s(p) даёт в итоге искомую функцию K(s).
Посі роенная микромодель верифицировалась по известным [82, 39] экспериментальным данным дренирования лабораторных засыпок 12/20,20/30, 30/40. Построенные гидравлические функции продемонстрировали вполне удовлетворительное согласование с экспериментальными данными как по капиллярной кривой, так и по относительной фазовой проницаемости для всех трёх песков. Дополнительным аргументом в пользу адекватности построенной модели можно назвать хорошую корреляцию между найденными в результате подбора параметров величинами среднего размера пор и известными из экспериментов средними размерами зёрен пористой засыпки. Кроме того, тестирование построенной микромодели на других (помимо гамма распределения) двухпараметрических плотностях распределения пор по размерам дало приблизительно те же значения параметров модели, что и гамма распределение.
Найденные в ходе верификации по данным дренажа параметры модели были использованы в параграфе 1.4 для определения гидравлических функций при пропитке. При построении модели пропитки дополнительно учитывался механизм кооперативного заполнения пор, а также наличие в среде защемлённого воздуха. Отличительной чертой предлагаемой модели является предсказание ею у капиллярной кривой остаточной по воздуху насыщенности sair < 1. Наличие такой предельной насыщенности — хорошо известный из опытов факт. Более того, вычисленные значения sair = 0.823,0.874,0.879 для сред 12/20, 20/30 и 30/40 соответственно, лежат в известной из экспериментов вилке 0.8 < sair < 0.9.
Наконец, в параграфе 1.5 строится модель первичной (быстрой) пропитки. Высокие темпы роста влажности изначально сухой пористой среды позволяют в этой ситуации пренебречь медленными (плёночные и уголковые течения) механизмами заполнения пор по сравнению с поршневым вытеснением в капиллярах и кооперативным заполнением. Особенностью этой модели является немонотонность кривой р — P(s) первичной пропитки.
Введение
Разумеется, в реальных натурных экспериментах подобная немонотонность не может быть реализована и переход из сухого во влажное состояние будет осуществляться скачком. Соответствующая этому скачку характерная полка давления на кривой первичной пропитки, соединяющая состояние нулевой насыщенности с близким к полному насыщению состоянием отмечается практически всеми исследователями, работающими с лабораторными песками.
Полученные результаты позволили сделать заключение о существенной зависимости общей картины протекания процесса пропитки от темпа его протекания и, следовательно, о необходимости модификации традиционных моделей влагопереноса учётом в них динамических (релаксационных) эффектов. В зависимости от темпа процесса пропитки элемент пористой среды в фазовой плоскости (s,p) должен следовать по траектории, промежуточной между первичной (бесконечно быстрый процесс) и главной (бесконечно медленный процесс) кривыми пропитки.
Вторая глава диссертационной работы посвящена анализу макроскопических уравнений влагопереноса и их возможных релаксационных модификаций.
В параграфе 2.1 проводится анализ устойчивости традиционной для механики грунтов модели Ричардса, которая при отсутствии источников или стоков жидкости имеет вид
- - V K(s)Vp - -gK(s) = 0, (6)
P = P(s). (7)
Здесь s — относительная водонасыщенность (s Є [0, 1]), р — давление жидкости, К — относительная фазовая проницаемость, z — вертикальная координата, противоположная по направлению силе тяжести. Функция P(s) определяет равновесное капиллярное давление, отвечающее заданному уровню водонасыщенности. Оказалось, что любое решение уравнения Ричардса является устойчивым, тем самым делая это уравнение непригодным для моделирования процессов пальцеобразования.
Естественно, поэтому, что следующим шагом исследования стал анализ модифицированной модели Ричардса. Обсуждаются такие модификации модели Ричардса, при которых изменению подвергается лишь соотношение между давлением и водонасыщенностью (7). Удалось доказать, что решения типа бегущей волны модифицированной модели могут быть неустойчивыми и получить критерий неустойчивости. Полученный критерий показывает, что неустойчивость фронтов пропитки для общего случая модификации модели Ричардса связана с немонотонностью профиля давления в базовом решении типа бегущей волны. Поэтому для того, чтобы моди-
Введение
фицированная модель Ричардса обладала свойством неустойчивости, необходимым для воспроизводства феномена пальцеобразования, нужно, чтобы оно продуцировало немонотонность в решениях типа бегущей волны. Это требование к гипотетической модели полностью согласуется с результатами экспериментальных работ.
Вид модифицированной модели конкретизируется в параграфе 2.2 на основе термодинамического подхода Hassanizadeh. Всякая релаксационная модификация модели Ричардса должна удовлетворять энтропийному неравенству, которое выводится на основе балансовых соотношений для свободной энергии Гельмгольца на границе раздела фаз и может быть записано в
виде
~(р - P(s)) ^ 0.
В качестве таких модификаций далее привлекаются модели Р- и S-релаксации. Первая из них заменяет зависимость (7) на релаксационный закон вида
ds
r^=p-P(s), (8)
а вторая — на соотношение
Ещё одним фактором, принятым в расчёт при построении модифицированной модели, стал учёт гистерезисного характера зависимости давления от насыщенности.
В параграфе 2.3 рассматриваются решения типа бегущих волн для описанных выше релаксационных моделей. Анализ фазовой плоскости модифицированных уравнений позволил получить критерий, связывающий скорость распространения волны пропитки и её форму (монотонная/немонотонная). Было показано, что при превышении скоростью волны некоторой критической величины, зависящей от начальных условий задачи, профили насыщенности и давления с необходимостью будут немонотонными. Это свойство решения, как было доказано выше, приводит к неустойчивости фронта пропитки при низкочастотных его возмущениях и к образованию пальцев. Здесь же демонстрируется важная роль гистерезиса в подавлении нефизичных (осциллирующих) решений.
Критерий устойчивости волны пропитки, полученный ранее, был выведен в предположении малости частоты возмущений ш- исходного решения типа бегущей волны. Параграф 2.4 обращается к определению устойчивости фронта пропитки во всём интервале изменения ал Для разрешения этого
Введение
вопроса приходится прибегать к численному решению спектральной задачи. Анализ решения показывает, что неустойчивость фронта пропитки при низких возмущающих частотах сменяется с ростом частоты устойчивостью. Имеется также характерное значение частоты uq, на которой возмущения растут быстрее всего. Естественно связывать эту величину с характерным расстоянием между пальцами, на которые распадается изначально однородный гравитационный фронт пропитки.
В параграфе 2.5 было проведено двумерное моделирование распространения одиночного пальца и распада фронта пропитки на систему пальцев на основе системы уравнений (6), (8). В качестве зависимости r(s,p) выбиралась простейшая функция вида т = то (—р). Численная реализация опиралась на схему, предложенную ранее Р. 3. Даутовым и А. Г. Егоровым. Численные расчёты подтвердили, что выявленная в предыдущих параграфах неустойчивость фронта пропитки приводит к его распаду на систему устойчиво развивающихся пальцев. Характерное расстояние между пальцами примерно соответствовало теоретически предсказанному. Размытие пальцев за счёт оттока влаги из них в окружающую сухую среду обеспечивается гистерезисом капиллярной кривой при дренаже и пропитке. Морфологические особенности пальца (область питания, тело пальца, головка пальца) качественно совпадают с выявленными в экспериментах, а профиль влажности вдоль оси пальца — с полученным в результате решения автомодельной задачи.
Заключительная третья глава диссертационной работы посвящена верификации построенной в предыдущей главе релаксационной модели влаго-переноса на основании экспериментов D. A. DiCarlo по одномерной гравитационной пропитке пористых сред. Результаты этих экспериментов в силу своей полноты и широты охвата предоставляют хорошую базу для тестирования теоретических моделей.
В параграфе 3.1 даётся описание экспериментов D. A. DiCarlo, процедур по обработке их результатов, указываются их характерные особенности. К таким особенностям, в частности, относятся:
зависимость вида профиля насыщенности от величины потока жидкости (монотонные профили для очень больших и очень маленьких потоков, немонотонные в промежуточном случае);
зависимость вида профиля насыщенности от величины начальной во-донасыщенности среды (немонотонные профили для сухой среды и монотонные для влажной);
почти постоянная скорость перемещения пальца.
Введение
В параграфе 3.2 обсуждается вид возможной релаксационной модификации модели Ричардса, в качестве которой была выбрана модель Р-релаксации
Далее, в параграфе 3.3 формулируется исходная постановка задачи одномерного влагопереноса в пористой среде, записываются начальные и граничные условия, приводятся соотношения для сканирующей дренажной кривой первого порядка.
Постановка задачи завершается в параграфе 3.4, в котором строятся аппроксимации главных кривых пропитки и дренажа, а также зависимость гидравлической проницаемости от насыщенности по экспериментальным данным D. A. DiCarlo [39] и Schroth et al. [82].
Параграф 3.5 описывает процедуру определения решения типа бегущей
волны модифицированной задачи. В связи с необходимостью учитывать ги
стерезис, задачу приходится решать в два этапа: первоначально находится
часть решения, соответствующая фазе быстрой пропитки; релаксация игра
ет здесь важнейшую роль. Дренирование, по сравнению с пропиткой, про
текает значительно медленнее и релаксационными эффектами на этом этапе
можно пренебречь, имея в виду, однако, что зависимость давления от насы
щенности будет определяться уже не главной, а сканирующей капиллярной
кривой. /
Параграф 3.6 посвящен определению вида коэффициента релаксации r(s,p). Для этого было использовано предположение о том, что он может быть представлен в виде произведения двух независимых функций насыщенности и давления: r(s,p) = rs(s)rp(p). Функции тр и rs определялись из условия согласования результатов вычисления зависимости насыщенности в голове пальца с соответствующими экспериментальными данными.
Адекватность построенной релаксационной модели проверялась в параграфе 3.7, где полученные с её помощью результаты сравнивались с результатами экспериментов D. A. DiCarlo. Сопоставление установившихся профилей насыщенности показало хорошее согласование рассчитанных и экспериментальных кривых во всём диапазоне изменения потоков. Столь же хорошее качество согласования продемонстрировали и зависимости насыщенности в хвосте и голове пальца от величины потока. Эти же зависимости обнаружили качественно сходное с экспериментами поведение при изменении начальной насыщенности среды, а именно исчезновение немонотонности профилей насыщенности по мере приближения этой величины к 0.02.
Основные результаты диссертационной работы:
1. Построена микромодель пористой среды, учитывающая взаимодействие пор и использующая приближение среднего поля для вычисления зависимости капиллярного давления и проницаемости от водона-
Введение
сыщенности. Применение этой модели к анализу процесса пропитки позволило сделать вывод о необходимости модификации традиционных моделей влагоперепоса учётом в них релаксационных эффектов.
Проведён анализ устойчивости модели Ричардса влагопереноса в пористой среде, а также её модификаций, учитывающих релаксационные свойства среды. Показано, что релаксационная модификация модели Ричардса способна воспроизвести экспериментально наблюдаемое явление пальцеобразования.
Построена и верифицирована по результатам экспериментов D. A. Di-Carlo по одномерной гравитационной пропитке модель влагопереноса в ненасыщенной пористой среде, учитывающая эффекты релаксации и гистерезисный характер зависимостей гидравлических функций пористой среды; получены аппроксимации коэффициента релаксации для различных лабораторных сред.
Основные результаты опубликованы в статьях [8, 16, 11, 12] и в тезисах [9, 46, 10]; они также докладывались и обсуждались:
на XVII сессии Международной школы по моделям механики сплошной среды (Казань, 2004);
на Международной конференции «Актуальные проблемы математики и механики» (Казань, 2004);
на Международном семинаре Recent Advances in Multi-phase Flow in Porous Media (Казань, 2004);
на Международной конференции 2004 American Geophysical Union Fall Meeting (San Francisco, 2004);
на Международном семинаре Upscaling Flow and Transport Process in Porous Media (Delft, 2005);
на Четвёртой молодёжной научной школе-конференции «Лобачевские чтения-2005» (Казань, 2005);
на Международном семинаре Summer School on Upscaling and Modelling of Coupled Transport Processes in the Subsurface (Utrecht, 2006);
на XVIII сессии Международной школы по моделям механики сплошной среды (Саратов, 2007);
Введение
на итоговых научных конференциях Казанского государственного университета (2003-2008);
на семинарах отделения механики НИИ математики и механики им. Н. Г. Чеботарёва Казанского государственного университета.
В совместных работах постановка задач, выбор направлений и методов исследований принадлежат соавторам, научному руководителю докт. физ.-мат. наук А. Г. Егорову и профессору Р. 3. Даутову.
Автор выражает глубокую благодарность своему научному руководителю докт. физ.-мат. наук А. Г. Егорову за предложенную тему исследований и поддержку при выполнении работы, профессору Р. 3. Даутову за плодотворные дискуссии и D. A. DiCarlo за предоставленные результаты лабораторных экспериментов.
Работа была выполнена при финансовой поддержке РФФИ (проект № 05-01-00516) и совместного проекта РФФИ и НВО (№ 05-01-890001-НВО).
Модель связки цилиндрических капилляров
Угроза поверхностного заражения подземных вод токсичными веществами при авариях на промышленных предприятиях и складах химических продуктов является в наше время актуальной экологической проблемой. Попавшие на дневную поверхность водорастворимые загрязнения мигрируют с потоком влаги через зону аэрации к зеркалу грунтовых вод и могут в конечном итоге нанести непоправимый ущерб природным экосистемам обширного региона. Защитные свойства зоны аэрации, связанные с адсорбцией загрязнителя частицами почвы, во многом определяются режимом фильтрационных потоков.
Многочисленные эксперименты [55, 45, 83, 64] свидетельствуют о том, что в ненасыщенных грунтах движимый гравитацией однородный фронт пропитки, как правило, распадается на устойчиво развивающуюся систему потоков («пальцев»). Понимание причин, вызывающих неустойчивость таких фронтов, и учёт этого эффекта необходимы для верного предсказания интенсивности переноса влаги и водорастворимых загрязнений от дневной поверхности к зеркалу грунтовых вод. Действительно, наличие предпочтительных путей, порождаемых распадом фронта пропитки на отдельные пальцы, существенно уменьшает время миграции влаги в зоне аэрации, снижая тем самым её защитную роль.
Практическая важность и научная значимость феномена пальцеобразо-вания диктуют необходимость разработки адекватной математической модели этого явления. К сожалению, традиционные модели влагопереиоса (например модель Ричардса) оказываются в данной ситуации малопригодными, так как не принимают в расчёт динамических эффектов, играющих значительную роль при образовании и развитии пальцев. Многочисленные попытки релаксационной модификации модели Ричардса дали определённое понимание некоторых аспектов проблемы, однако построение целостной модели остаётся всё ещё делом будущего.
Одной из характерных особенностей пальцев, отмечаемой всеми экспериментаторами, является немонотонность профилей давления и насыщенности в их поперечном сечении. Эта немонотонность приводит к тому, что насыщенность в теле пальца оказывается значительно выше, чем насыщенность вне его. В то же время поток влаги, определяемый градиентом давления, оказывается направлен к оси пальца. В этой ситуации высоконасыщенный палец не отдает влагу окружающей пористой среде, а, напротив, иссушает её.
На самом деле возможность движения влаги в пористой среде в направлении градиента влажности (т. е. из областей с малой влажностью — в области с большой) была впервые обнаружена не в экспериментах по паль-цеобразованию, хотя и приблизительно в то же время. По-видимому первые наблюдения этого явления были сделаны в работах Абрамовой [1,2] и Аллера [47, 48, 49]. Позднее опыты Абрамовой по сути были повторены Collis-George et al. [35]; аналогичные явления наблюдались в экспериментах Дмитриева [14] и Бондаренко [5], а также были воспроизведены Роде и Романовой [24]. Именно эти эксперименты впервые поставили под сомнение применимость модели Ричардса для описания переноса влаги в ненасыщенной пористой среде.
Опишем вкратце суть поставленных экспериментов. Во всех них исследователи изучали движение «подвешенной» (по терминологии А. Ф. Лебедева [19]) влаги в вертикальной колонке, заполненной образцом природного грунта. Грунт в колонке промачивали водой, после чего с её верхней поверхности было организовано испарение (как при помощи нагрева электрической лампой, так и естественным путём без нагревания при комнатной температуре). Давление снималось при помощи ряда установленных с определённым шагом тензиометров, влажность определялась по исследованию взятых с разной глубины образцов почвы.
Во всех экспериментах исследователи отмечали восходящее движение влаги из глубины почвы к поверхности испарения, которое, в особенности на первых стадиях процесса, совершалось при отсутствии градиента влажности и иногда даже в направлении возрастающей влажности.
Гидравлические функции при пропитке
Описанные эксперименты не могли быть объяснены в рамках традиционно используемой в почвоведении диффузионной теории. Согласно ей, перенос влаги в ненасыщенной пористой среде описывается уравнением Ричардса [77, 78]: ds О ( ds\ Здесь s — влажность, t — время, х — пространственная координата, D(s) — коэффициент диффузивности, D(s) = Кдф/ds, где К — коэффициент вла-гопроводности, ip = i/)(s) — капиллярный потенциал влажности. Кроме натурных экспериментов, о неприменимости диффузной теории по крайней мере для нестационарных процессов свидетельствуют и результаты численного моделирования. В качестве примера можно сослаться на работу Бондаренко и др. [6], в которой исследовалась степень отклонения численно рассчитанных по уравнению (1) значений влажности s от экспериментально наблюдаемых значений на различных стадиях сушки. Из сравнения этих величин видно, что расчётные значения даже качественно не совпадают с экспериментом на первых стадиях процесса. Аллер был также первым, кто попытался теоретически обосновать описанный феномен [49, 50]. Для этого он модифицировал уравнение (1), заменив в нём капиллярный потенциал ф на так называемый эффективный потенциал влажности фе. Этот эффективный потенциал конструируется в виде суммы фе = ф + А—, где А — некий варьируемый коэффициент. После этой подстановки уравнение (1) принимает вид Os д /_. . 9s л d2s \ Предложение Аллера вызывает два естественных вопроса: а) воспроизводит ли указанная модификация экспериментально наблюдаемые факты и б) насколько данная модификация является теоретически обоснованной.
Что касается первого вопроса, то ответ на него можно считать положительным. Численные эксперименты [22, 26], проведённые по модели Аллера для значений А = 0 и А = 100 показали, что при А = 100 происходит поток влаги вдоль градиента влажности, в то время как при А = 0 этот эффект не наблюдается. В другой работе [7] численные расчёты по модели Аллера сравнивались с экспериментальными профилями влажности на различных стадиях сушки. При А — const экспериментальные и рассчитанные кривые демонстрируют неплохое согласование на начальных стадиях процесса. При этом было отмечено, что величина А не остаётся постоянной в процессе сушки, хотя никаких предположений о виде зависимости А от времени не было сделано.
Несмотря на всё сказанное, против модели Аллера были выдвинуты многочисленные возражения. Критики указывали на то, что уравнение (2) и рассуждения, приводящие к нему, аналогичны таким, которые обычно используются для описания движения жидкости в трещиновато-пористых средах. Однако обычная почва имеет существенно более однородную структуру, и в ней едва ли можно выделить крупные транспортные каналы и мелкие питающие капилляры, как это сделал Аллер.
По мнению Роде [25], теоретическое обоснование, выдвигаемое Адлером, является чисто умозрительным. Роде считает, что Аллер попросту не принял в расчёт гистерезис зависимости влажности от давления. Между тем отчётливо показано [66, 99, 15], что после того, как инфильтрация сменяется испарением, распределение давления вследствие влияния гистерезиса должно сразу же измениться. Это создаёт предпосылку к передвижению влаги к испаряющей поверхности не только без появления соответствующего градиента влажности, но даже в направлении этого градиента, хотя и против градиента давления. Имеются теоретические расчёты [17], показывающие, что обычная модель Ричардса с гистерезисом ip{s) даст движение воды в направлении градиента влажности.
Среди других подходов к модификации модели Ричардса нельзя не упомянуть работы А. В. Лыкова [20, 21]. Лыков исходит из совершенно других соображений, чем Аллер, хотя его уравнение также предназначено для более полного учёта нестационарности. Вывод Аллера базируется на рассмотрении капиллярной модели, тогда как Лыков получает своё уравнение на основе методов термодинамики необратимых процессов.
Анализ устойчивости модифицированной модели Ричардса
Кулик [18] исследовал уравнения Аллера и Лыкова с точки зрения инвариантности относительно непрерывных групп преобразований. Его анализ показывает, что ни то, ни другое уравнение нельзя принять в качестве универсального уравнения движения влаги в почве. Модель Аллера, созданная специально для объяснения экспериментов по оттоку влаги при испарении, оказывется неспособна описать приток влаги в почву при пропитке. В то же время модель Лыкова, исходящая из экспериментов по инфильтрации жидкости в горизонтальную колонку, входит в противоречие с другой группой экспериментов. Кулик видит выход из этой ситуации в привнесении в модель гистерезисных свойств. По его мнению для того, чтобы быть равно пригодной как для пропитки, так и для дренажа, в модели Аллера коэффициент А должен испытывать очень сильный гистерезис. Коэффициент Лыкова AL должен либо быть очень малым, либо также быть подвержен сильному гистерезису. Кулик также предлагает своеобразный гибрид уравнений (2) и (3) в следующем виде:
Os я d2s д / ,ds\ л д2 fds\ д ,„ . где коэффициенты А и А] должны испытывать очень сильный гистерезис. Возрождение интереса к модификации уравнения Ричардса было связано с появлением большого количества экспериментальных фактов, касающихся экспериментов по пальцеобразованию. Как и в случае с ранними экспериментами по движению влаги в природных грунтах, описать этот феномен в рамках традиционной модели не оказалось возможным. Интересно отметить, что и в этом случае при построении модели весьма продуктивным оказалось привлечение термодинамических соображений, а также учёт в явном виде гистерезиса гидравлических функций. Наибольший вклад в теорию нестационарных течений в ненасыщенных пористых средах в настоящее время принадлежит S. Majid Hassanizadeh. В серии работ [51, 52, 53, 54] он последовательно применяет методы неравновесной термодинамики к процессам перераспределения влаги в поровом пространстве. Записывая балансовые соотношения для свободной энергии Гельм-гольца на границе раздела фаз, Hassanizadeh получает необходимое ограничение, накладываемое на возможные релаксационные модификации модели Ричардса в виде энтропийного неравенства Os SL(p-P(s)) 0. Им также была предложена одна из таких модификаций, известная как модель Р-релаксации: ds T—=p-P{s). (4) Тщательный математический анализ этой модели был выполнен Cuesta et al. [36]. Предполагая, что зависимости P(s) и K(s) имеют степенной характер, ими было доказано существование и единственность решения типа бегущей волны уравнений (1), (4). Также был получен критерий немонотонности таких решений и исследовано их поведение при стремлении начальной водонасыщенности среды к нулю.
Привнесение гистерезисных свойств в предложенную Hassanizadeh модель было выполнено в работе Беляева и Hassanizadeh [30]. Авторы использовали простейшую модель гистерезиса (так называемую "playype"), в которой все сканирующие кривые представляли собой вертикальные отрезки s = const.
Анализ устойчивости решения типа бегущей волны уравнения Ричардса, а также некоторых её модификаций (релаксационной и стефановского типа) был проведён Егоровым и Даутовым [13]. Они доказали, что модель Ричардса в своём исходном виде абсолютно устойчива, а её стефановская модификация — абсолютно неустойчива и, следовательно, оба этих случая не подходят для моделирования пальцев. В отличие от них, релаксационная модификация демонстрирует условную неустойчивость, т. е. устойчивость на высокочастотных возмущениях и неустойчивость на низкочастотных. Это приводит к наличию минимума на дисперсионной кривой, положение которого определяет характерный масштаб системы пальцев, генерируемых изначально однородным фронтом пропитки.
Определяющие уравнения
Разумеется, всякая математическая модель должна быть тщательно протестирована и верифицирована по результатам натурных экспериментов. Только после этого можно делать окончательное заключение о её достоверности и практической применимости. С другой стороны, даже самый скрупулезный и аккуратный лабораторный эксперимент немного стоит без надёжного теоретического фундамента и непротиворечивых моделей. Накопленный массив экспериментальных данных с одной стороны и достижения теоретиков с другой позволяют в настоящее время значительно приблизиться к построению законченной модели влагопереноса в сухих средах.
Изложенное выше определяет основную цель данной работы, заключающуюся в построении и верификации математической модели влагопереноса в ненасыщенном грунте, способной адекватно описать процессы образования и развития пальцев. Остановимся подробнее на содержании диссертации. Она состоит из введения, трёх глав, заключения, 61 рисунка, 7 таблиц. Список использованной литературы содержит 99 наименований. Первая глава работы посвящена исследованию особенностей процессов перераспределения влаги в пористой среде на основе микромоделирования. Первоначально исследуются процессы перераспределения влаги в масштабе пор. В параграфе 1.1 рассматривается модель пористой среды в виде связки цилиндрических капилляров различного размера, не взаимодействующих между собой. Подробно излагаются положения микромодели Туллера-Ора, в которой эти капилляры имеют полигональное поперечное сечение. Подход Туллера-Ора является более реалистичным по сравнению с моделью цилиндрических капилляров, поскольку позволяет учитывать плёночные и уголковые течения жидкости, доминирующие при малых насыщен-ностях пористой среды. Тем не менее игнорирование взаимодействия пор является его принципиальным недостатком.
Параграф 1.2 представляет модель, преодолевающую этот недостаток путём моделирования пористой среды в виде решётки со связями — цилиндрическими капиллярами различного поперечного сечения. Определя- ющими параметрами для такой модели являются координационное число решётки Z, равное среднему числу выходящих из одного узла связей, и плотность распределения поперечных размеров связей (капилляров). Экспериментально установлено, что для типичных лабораторных песков, с которыми обычно имеют дело исследователи, это число в среднем равно 4.6. В качестве функции распределения пор по размерам использовалось гамма-распределение и fn{L) = [n{L I х, L0) = ,,.,., 77 exp(-L/L0), І +Г(х+1) где LQ, x — параметры функции распределения. Выбор для нашей модели двухпараметрического закона распределения позволяет согласовать её с известными эмпирическими моделями Брукса-Кери и Ван-Генухтена, в которые также входят два определяющих параметра. Кроме того, дополнительным параметром модели, повышающим её адекватность реальным пористым средам является доля капилляров различного поперечного сечения. Для предлагаемой модели в качестве таких капилляров были выбраны криволинейный треугольник и круг.
Сформулированная выше микромодель была использована в параграфе 1.3 для вычисления гидравлических функций (капиллярного давления и гидравлической проницаемости) пористой среды. Эти функции определялись осреднением насыщенностеи отдельных капилляров по всему массиву пор с учётом взаимодействия между ними. Зависимость водонасыщенности среды от давления s(p) может быть представлена в виде суммы вкладов жидкости, сосредоточенной в заполненных капиллярах, уголковых менисках пор и в плёнках на твёрдой поверхности: s(p) = scap(p)4-smcn(p)-fSfiim(p). Асимптотический анализ поведения этих функций при р — со показал, что влияние адсорбированных пленок начинает сказываться лишь при очень низких насыщенностях s 10 4, и в типичных ситуациях им можно пренебречь. В то же время мениски в углах пор играют существенную роль при влажностях менее 20% и их влияние необходимо учитывать.
Другой важной гидравлической характеристикой пористой среды является зависимость относительной фазовой проницаемости К от насыщенности. Для её вычисления была применена известная электрическая аналогия, сводящая проблему отыскания фазовой проницаемости к задаче нахождения эффективной проводимости решётки по заданному распределению про-водимостей её связей. Считая, что проводимости отдельных пор распределены независимо друг от друга и используя хорошо известную аппроксимацию среднего поля, эффективная проводимость решётки К (р) может быть Здесь k — проводимость отдельной связи, а под средним понимается математическое ожидание относительно заданной плотности fn(k), которая восстанавливается по известной плотности /„(L) распределения пор по размерам. Полученная в результате решения зависимость относительной фазовой проницаемости К от давления р вместе с найденной ранее зависимостью водонасыщенности среды от давления s(p) даёт в итоге искомую функцию K(s).
Посі роенная микромодель верифицировалась по известным [82, 39] экспериментальным данным дренирования лабораторных засыпок 12/20,20/30, 30/40. Построенные гидравлические функции продемонстрировали вполне удовлетворительное согласование с экспериментальными данными как по капиллярной кривой, так и по относительной фазовой проницаемости для всех трёх песков. Дополнительным аргументом в пользу адекватности построенной модели можно назвать хорошую корреляцию между найденными в результате подбора параметров величинами среднего размера пор и известными из экспериментов средними размерами зёрен пористой засыпки. Кроме того, тестирование построенной микромодели на других (помимо гамма распределения) двухпараметрических плотностях распределения пор по размерам дало приблизительно те же значения параметров модели, что и гамма распределение.
Найденные в ходе верификации по данным дренажа параметры модели были использованы в параграфе 1.4 для определения гидравлических функций при пропитке. При построении модели пропитки дополнительно учитывался механизм кооперативного заполнения пор, а также наличие в среде защемлённого воздуха. Отличительной чертой предлагаемой модели является предсказание ею у капиллярной кривой остаточной по воздуху насыщенности sair 1. Наличие такой предельной насыщенности — хорошо известный из опытов факт. Более того, вычисленные значения sair = 0.823,0.874,0.879 для сред 12/20, 20/30 и 30/40 соответственно, лежат в известной из экспериментов вилке 0.8 sair 0.9.
Наконец, в параграфе 1.5 строится модель первичной (быстрой) пропитки. Высокие темпы роста влажности изначально сухой пористой среды позволяют в этой ситуации пренебречь медленными (плёночные и уголковые течения) механизмами заполнения пор по сравнению с поршневым вытеснением в капиллярах и кооперативным заполнением. Особенностью этой модели является немонотонность кривой р — P(s) первичной пропитки.
Разумеется, в реальных натурных экспериментах подобная немонотонность не может быть реализована и переход из сухого во влажное состояние будет осуществляться скачком. Соответствующая этому скачку характерная полка давления на кривой первичной пропитки, соединяющая состояние нулевой насыщенности с близким к полному насыщению состоянием отмечается практически всеми исследователями, работающими с лабораторными песками.
Полученные результаты позволили сделать заключение о существенной зависимости общей картины протекания процесса пропитки от темпа его протекания и, следовательно, о необходимости модификации традиционных моделей влагопереноса учётом в них динамических (релаксационных) эффектов. В зависимости от темпа процесса пропитки элемент пористой среды в фазовой плоскости (s,p) должен следовать по траектории, промежуточной между первичной (бесконечно быстрый процесс) и главной (бесконечно медленный процесс) кривыми пропитки.
