CC BY
28
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Том 14 7, ки. 3
Физико-математические пауки
2005
УДК 532.546
ВЫЧИСЛЕНИЕ ГИДРАВЛИЧЕСКИХ ФУНКЦИЙ НА РЕШЁТЧАТОЙ МОДЕЛИ ПОРИСТОЙ СРЕДЫ В ПРИБЛИЖЕНИИ СРЕДНЕГО ПОЛЯ
Б.Э. Гоголашвили, А.Г. Егоров
Аннотация
В приближении среднего поля изучены процессы влагоперепоса па решётчатой микромодели пористой среды. В рамках микромодели учтены капиллярная, пленочная и уголковая компоненты влагосодержапия и основные механизмы изменения состояния отдельного капилляра (поршневое вытеснение, "snap-off", кооперативное заполнение). На основе предложенного подхода определены кривые дрепажа/пропитки (главные и первичные) и кривые фазовых пропицаемостей микромодели. За счет выбора параметров микромодели удалось добиться хорошего согласования построенных кривых с экспериментальными данными для отсортированных лабораторных песков.
Введение
Явление пальцеобразования при пропитке сухих пористых сред известно достаточно давно. Повышенный интерес к этому феномену в последнее время вызван экологическими проблемами и связан с потенциальной возможностью быстрого переноса загрязнений с поверхности почвы к зеркалу грунтовых вод по системе вертикальных пальцев [1 3].
Теоретическое изучение процесса пальцеобразования во многом лимитируется тем обстоятельством, что. с одной стороны, он критическим образом зависит от вида функций капиллярного давления Р(в) и относительной фазовой проницаемости К (в) в области малых влажностей в [4, 5], а с другой стороны, определение этих функций при малых в встречает известные экспериментальные трудности. Возможным выходом из этой ситуации может стать микромоделирование процесса перераспределения влаги в пористой среде с тем. чтобы получить гидравлические кривые Р (в) и К (в) во всём диапазоне изменения насыщенности [6].
Стандартным средством микромоделирования в последние годы стали решётчатые модели пористой среды [7 14]. Они позволяют учесть как основные составляющие влагосодержапия (капиллярную, уголковую, плёночную), так и основные механизмы изменения состояния пор (узлов решётки) и соединяющих их капилляров (связей решётки). Анализ процесса перераспределения влаги на решётке проводится, как правило, с помощью прямого вычислительного эксперимента. Возможной альтернативой этому может служить рассматриваемый в данной работе аналитический подход, опирающийся на известное приближение среднего поля [15].
1. Микромодель
1.1. Решётка. Пористая среда моделируется решёткой со связями цилиндрическими капиллярами различного поперечного сечения. Узлы решётки являются простым пересечением капилляров. Объём узлов игнорируется, перераспределяясь между соответствующими связями, а размер узла считается совпадающим с максимальным размером выходящих из узла капилляров.
Рис. 1. Типичпые сочопия капилляров
Двумя наиболее важными характеристиками решётки являются координационное число 2, определяющее среднее число выходящих из одного узла связей, и плотность распределения поперечных размеров связей (капилляров). Основной интерес для данной работы представляют используемые в лабораторных экспериментах пористые засыпки из неконсолидированных сферических частиц. Ряд исследователей [7. 16] определяли среднее координационное число для таких сред и нашли его равным 2 = 4.6. Поэтому в дальнейших расчетах величина 2 будет варьироваться вблизи указанного значения 2 € [4, 6].
Эмпирические модели для определения гидравлических характеристик пористых сред, наиболее известными из которых являются модели Брукса Кори [17] и Ван-Генухтена [18]. включают в себя два параметра. Естественно поэтому использовать двухпараметрические законы распределения пор по размерам. В данной работе в качестве такого закона используется гамма-распределение
/п(Ь) = X), /лг(у,х) = , -.4 ехр(-у).
Г(х + 1)
Здесь Ь0, х _ параметры функции распределения, fn(L) - частичная плотность вероятностей:
число пор с размерами из интервала [Ь, Ь + 3,Ъ]
/п(-ь) аЬ =---.
общее число пор
Соответствующие данным плотностям / интегральные функции распределения будем обозначать через
1.2. Капилляры. Ограничимся рассмотрением капилляров двух различных типов: криволинейный треугольник и окружность (см. рис. 1). Их выбор в качестве базовых, на наш взгляд, в большей степени отвечает внутренней геометрии сферических упаковок, чем правильные многоугольники, используемые при моделировании консолидированных пористых сред [19, 20]. Размер капилляра Ь определялся как радиус вписанной в него окружности. Доли капилляров первого и второго типов (а и 1 — а соответственно) фиксировались одними и теми же независимо от Ь.
Степень Б заполнения водой поры данного размера зависит от уровня давления р < 0 в системе. При низком давлении величина Б определяется толщиной к(р) адсорбированных на твердой поверхности водных пленок и радиусом кривизны Гтеп(р) жидких менисков в углах поры (рис. 1). Толщина плёнок вычисляется через константу Гамакера Ан < 0 [21], а радиус кривизны - через поверхностное натяжение а на границе вода - воздух по формулам
Ч'Р) = (/тг^ гтсп(р) = --. у 6пр р
Вклад 5fiim жидких плёнок в S = Sfiim + Smen подсчитывался как
с / \
¿film (Р) = Ol——,
где константа ai равна, соответственно, 3.014 и 2 для пор первого и второго типов. Вклад Smen менисков для треугольной поры вычисляется из простых геометри-ческтх соображений. Соответствующая формула может быть аппроксимирована выражением
Sme п(е) = a2i£1'5 + a22£2 + a23£2'5, £ = rmen/L = —a/pL, a21 = 1.067, a22 = -0.699, a23 = 0.173,
с относительной погрешностью не более 1.4% во всём интервале [0,1] изменения £.
С ростом давления вплоть до некоторого критического значения p+(L) величина S непрерывно растёт. При p = p+ (L) граница раздела между водой и воздухом в поре становится неустойчивой, и пора «мгновенно» заполняется водой (так называемый механизм "snap-off" заполнения). Согласно [22]
а
Р+ =
а константа a+ равна единице для обоих типов пор.
Аналогичным образом происходит и осушение поры при снижении давления. Отличие состоит в том, что переход от режима полного заполнения поры к плёночно-уголковому режиму идет по иному сценарию, за счет движения мениска вдоль капилляра (режим поршневого вытеснения), и при ином значении критического давления
а
р~ = ~а-Т
Величина p_ есть не что иное, как давление, при котором внутри поры существует равновесный мениск, разделяющий в поперечном направлении воздух и воду. Сумма двух главных кривизн такого мениска, очевидно, должна совпадать с кривизной 1/r_ уголковых менисков. Следствием этого является простое уравнение
A(r_ )= г_П(г_) (2)
для нахождения критического радиуса r_ уголковых менисков (а, значит, и критического давления p_ = — a/r_). Здесь через A и П обозначены площадь и периметр заполненной воздухом центральной части капилляра. Для треугольной поры решение (2) даёт a_ = 1.75, а для круглой поры - a_ =2.
S
её проницаемость k. При полном заполнении поры проницаемость пропорциональна квадрату её размера: k(L) = a3L2. При частичном заполнении k складывается из проницаемости жидких плёнок kfijm = a5h3/L и проницаемости менисков в углах поры (для треугольной поры). Последняя вычислялась на основе решения методом конечных элементов задачи стоксовского течения внутри жидкого мениска. На твёрдых границах задавались условия прилипания, на внешней границе нормальная производная скорости принималась равной нулю, что отвечает контакту с невязкой средой (воздухом). Результаты расчёта оказалось возможным с относительной погрешностью менее 1% аппроксимировать приближённой формулой
kmen = L2A4(p) = L2(a4i£3'5 + a42£4 + a43£4'5 + a44£5). (3)
Коэффициенты в представленных выше зависимостях равны a3 = 1/8, a5 = = 4/3 (для круглой поры); a3 = 0.134, a5 = 2.01, a41 = 0.1524, a42 = —0.29 93, a43 = 0.25 34, a44 = —0.0857 (для треугольной поры).
2. Гидравлические функции при дренаже
2.1. Капиллярная кривая. Процесс осушения каждого капилляра зависит не только от общего давления в системе, но и от состояния соседних с ним капилляров. Если оба узла, на которые опирается данный капилляр, заполнены водой, то капилляр тоже останется заполненным даже при снижении давления ниже критического. Осушение в этом случае произойдёт лишь тогда, когда воздух проникнет в один из соседних узлов, иными словами, когда хотя бы один из 2(Z — 1) капилляров, являющихся непосредственными соседями данного капилляра, осушится.
Обозначим далее через Z(р) и n(p, L) числовую долю заполненных капилляров при данном уровне давления среди всех капилляров образца и среди капилляров размера L соответственно. Тогда вероятность того, что все 2(Z—1) соседей данного капилляра заполнены, равна Z2(Z-1), в силу чего
fl, p > P-(L),
n(P, L) = < (4)
(z2(Z-1), p < P-(L).
Z
POO
Z(p)=/ n(p,L)fn(L) dL. (5)
Jo
Подставляя (5) в (4) и обозначая через L-(p) обратную к p-(L) функцию L-(p) = = — a-a/p, получим для нахождения Z следующее алгебраическое уравнение
Z = b +(1 — b)Z2(Z-1), b(p) = F„(L-(p)). (6)
Это уравнение имеет при b < bz = (2Z — 3)/(2Z — 2) два корня. Один из них равен единице: физический интерес представляет второй расположенный между нулём и единицей. С превышением b критического значения bz единственным вещественным корнем уравнения (6) остаётся Z = 1 • Соответственно, зависимости Z(b) терпят излом в точке bz, продолжаясь правее неё прямой Z = 1 • Для различных значений Z эти графики представлены та рис. 2. Отметим, что кривые Z(b)
Z(b) = b
ют с пей при малых b. Последнее есть следствие того, что при b ^ 0 решенне Z
b
малости Z = b + O(b2(Z-1)), b ^ 0.
После определения величины Z функция n(p, L) иодсчитывается согласно (4), а полная водопасыщеппость s находится осреднением по всем порам:
s = s(p) = j™ (n + (1 — n)(es(1)(p, L) + (1 — e)S(2)(p, L))) fv(L) dL. (7)
Здесь через fv(L) обозначена объёмная плотность распределения пор по размерам, которая с точностью до нормирующего множителя совпадает с L2fn(L). Соответствующую данной плотности fv функцию распределения будем обозначать через Fv ■ Верхний индекс 1, 2 указывает на треугольные и круглые поры соответственно, объёмная доля в треугольных пор выражается через их числовую долю а как
(1)
о _ аа0_ m _ о 7оо (2) _ ^
Р — —щ-;-—(Г), ао -6.738, а0 -тт.
аад + (1 — а)а0
Коэффициент а0г) определяет площадь aO^L2 поперечного метения капилляра г-го типа.
Рис. 2. Зависимость доли заполненных пор Z °т параметра b при различных значениях Z для процесса первичного дренажа
Обозначая для простоты записи 1 — Z2(Z-1) через ^ и проводя интегрирование в (7). получаем окончательную расчетную формулу для определения капиллярной кривой в виде
s(p) = Scap(p) + «те n(p) + Sfilm (p), (8)
scap = 1 — M(z) + M(z)FV (a-z), (9)
Sülm = (a^ß + af^l - /3))At(z) v^G^a-z), (10)
3
sme„ = a21fc)z1+fc/2Gi+fc/2(a-z), (И)
fc=i
S= ' 2, Gk{z)= y~kfv{y)dy,
rmenL0 Jz
где fV (y) = L0fv (L0y). Здесь и далее через z обозначается обратное безразмерное давление z = — <r/pL0, коэффициенты a^ определены согласно (1).
Структурно правая часть (8) представляет собой сумму вкладов от полностью заполненных капилляров, уголковых менисков и плёнок. Типичные графики суммарной зависимости s(p) и её структурных составляющих scap, smen, sfiim представлены на рис. 3.
Наличие предельного pe значения давления, при превышении которого s = = 1, т. н. "air entry pressure", очевидным образом связано с критическим значением bz в уравнении (6): оно может быть найден о из условия b(pe) = bz- В терминах
z
27 3
FN(o,-ze) = ^--. (12)
Как видно из рис. 3, влияние адсорбированных пленок начинает сказываться лишь при очень низких насыщенностях s ~ 10-4, и в типичных ситуациях его можно игнорировать. В то же время учёт наличия менисков в углах пор необходим, особенно при рассмотрении малых (менее 20%) влажностей. Ещё одной характерной особенностью представленных графиков является резкое изменение вида зависимости s(p) при значениях насыщенности порядка 10%, связанное с переходом от режима полного заполнения пор к уголковому режиму.
а) Во всем диапазоне изменения насыщенно- 6) В области умеренных влажностеи
Рис. 3. Зависимость пасыщеппости от давления (в см водного столба) для грунта 20/30
2.2. Относительная фазовая проницаемость. Для вычисления относительной фазовой проницаемости К решётки капилляров можно воспользоваться известной электрической аналогией, сводящей проблему отыскания фазовой проницаемости к задаче нахождения эффективной проводимости решётки по заданному распределению проводимостей её связей. В таком контексте проводимость к каждой связи есть коэффициент пропорциональности между объёмным потоком влаги через соответствующий капилляр и перепадом давления в узлах, соединяемых данной связью. Следовательно, с точностью до несущественного множителя, проводимость определяется через проницаемость капилляра и площадь его поперечного сечения аоЬ2 так к = аоЬ2&.
В свою очередь, проницаемость капилляра зависит от его формы, размера и состояния (заполнен он полностью нлн нет):
а3 'Ь2, капилляр заполнен,
Ь2А4 ) (р) + а5г)^3(р)/Ь, в противном случае.
Здесь верхний индекс (г) указывает та тип капилляра, коэффициент а5 и функция
определяют, соответственно, плёночную и уголковую составляющие проницаемости.
Считая далее проводимости пор распределёнными независимо друг от друга и пользуясь хорошо известной аппроксимацией среднего поля [15], найдём эффективную проводимость решётки к*(р) как решение следующего алгебраического уравнения
2 = / к я = \ к + к„да -1)
Под средним здесь понимается математическое ожидание относительно заданной плотности /и(к), которая полностью восстанавливается по известным плотности /П(Ь) распределения пор по размерам, долям а, 1 — а пор различных типов и доли п(р, Ь) заполненных пор данного размера. С точностью до множителя, не зависящего от давления (который не существенен и определяется на завершающей стадии вычислений по условию нормировки К(1) = 1), найденная величина к*(р) совпадает с искомой зависимостью К (р).
(13)
Рис. 4. Зависимость относительной фазовой проницаемости от насыщенности
Переходя далее в правой части (13) к интегрированию по Ь, получим после элементарных преобразований следующую безразмерную форму алгебраического уравнения для нахождения К(р)
с Ж
7?=/ Ро1м{у)с1у-ф) I Ро/м{у)с1у + ф) I Р1/М{у)ау, (14)
^ «/0 «/ а— г а—г
(1) (1) 4 (2) (2) 4
*Ь = *Ь(у; К) = а +(!-«) У
а,
а01)а31)у4 + к ' 4 а02)а32)у4 + к' 01)у44(1)(./у) + Л1)^).™, а(2)а(2)
у4А4)(^/у) + а0 )а5 ^гу а0 )а^ )егу
Л = ЗД; К) = а , :» ■» 4 - ' » 5 ^ + (1 -
а)
а01)у4А41) (г/у) + а01)а51)егу + К а02)а52)егу + К'
Численное решение этого алгебраического уравнения не вызывает трудностей. Типичная зависимость относительной фазовой проницаемости К (в) представлена на рис. 4 сплошной линией. На этом же рисунке штриховой линией показана та же самая зависимость при игнорировании потока по адсорбированным жидким плёнкам (а551) = а52) = 0). Обе кривые практически совпадают друг с другом вплоть до достижения в некоторого порогового значения. Будем обозначать это пороговое значение через вге8 и говорить о нём как об остаточной насыщенности. Наличие порогового значения вГе8 обусловлено тем, что по его достижению система хорошо проводящих капилляров (полностью заполненные капилляры и капилляры с менисками) теряет связность. Поэтому в пренебрежении проницаемостью плохо проводящих капилляров потоки влаги в среде оказываются вообще невозможными (штриховая линия на рис. 4), а с её учётом они определяются плёночной составляющей проницаемости и чрезвычайно малы (сплошная линия на рис. 4).
Величина , отвечающая остаточной насыщенности в-8 при дренаже, может быть найдена из уравнения (14), если положить в нём а5 = 0, К = 0:
2
= 1 - (1 - (1 - FN{a-z-s)) .
Учитывая, что в силу (6) 1 — ^^(а_г) = (1 — £получим имеющее очевидный физический смысл соотношение
= !-(!_ а)(1-С(г-8)). (15)
а) Зависимость давления от насыщенности; сплошная линия дренаж, штрих-пунктирная линия главная кривая пропитки, пунктирная линия кривая первичной пропитки, маркеры экспериментальные данные для дренажа [23], sa[r — остаточная по воздуху насыщенность
6) Зависимость фазовой проницаемости от насыщенности; сплошная линия дренаж, штриховая линия пропитка, маркеры экспериментальные данные для дренажа [24]
Рис. 5. Сравнение теоретических гидравлических функций с экспериментальными данными для грунта 20/30
Его левая часть указывает па перколяционный порог решётки, вычисленный на основе приближения среднего поля. Правая часть подсчитывает этот порог как долю хорошо проводящих пор (величина (1 — а)(1 — С) есть доля плохо проводящих пор с исключительно плёночным режимом течения жидкости). Заметим далее, что найденное из (15) значение С мало, и, следовательно, хорошим приближением для С ) (см. рис. 2) является величина Ь = ^^ (а_я_8). Поэтому (15) можно упростить, придав ему аналогичную (12) форму
, _ 2 —
Рм{а-г1С5) = --—. (16)
(1 — а)^
3. Верификация модели по экспериментальным данным для дренажа
Построенная микромодель верифицировалась по известным [23, 24] экспериментальным данным дренирования лабораторных засыпок 12/20, 20/30, 30/40 со средним размером зерна 1.105, 0.713 и 0.532 мм соответственно. Код засыпки отвечает размерам сит, через которые просеивается песок при подготовке образца.
При заданном координационном числе решётки микромодель полностью определяется тремя параметрами: параметры Ь0 и х функции распределения и доля а треугольных капилляров. Значения Ь0, х> а подбирались так, чтобы согласовать наилучшим образом теоретические и экспериментальные кривые фазовой проницаемости и капиллярного давления. Результаты такого согласования для засыпки 20/30 при Z = 4 представлены на рис. 5 сплошными линиями. Как видно, совпадение теоретических кривых с экспериментальными данными более чем удовлетворительно. Оно обеспечилось выбором Ьо = 0.138 мм, х = 33, а = 0.36. При варьировании координационного числа Z в диапазоне [4, 5.5] параметры Ь0, х> а микромодели изменялись незначительно.
Столь же хорошее согласование между теоретическими и экспериментальными данными было получено и для двух других видов засыпок. При этом оптимальные значения параметров микромодели составили Ьо = 0.198 мм, х = 9, а = 0.205 для
засыпки 12/20 и L0 = 0.096 мм, х = 56, а = 0.33 для засыпки 30/40. Отметим, что найденные в результате подбора параметров оптимальные значения среднего размера пор L0 хорошо коррелируют со средним размером зёрен. Первые из них соотносятся между собой в пропорции 1:0.7:0.48, в то время как вторые в пропорции 1:0.65:0.48. Кроме того, относительное среднеквадратичное отклонение (х + 1)-1/2 размеров пор от среднего значения растёт при переходе от песков 30/40 к 20/30 и далее к 12/20, что также согласуется с увеличением разброса диаметров частиц, составляющих указанные среды.
Для проверки чувствительности микромодели к выбору функции распределения пор по размерам помимо гамма-распределения тестировалась и другая (лог-нормальная) двухпараметрпческая функция распределения. Средний размер зёрен, среднеквадратичное уклонение и доля треугольных капилляров при этом изменялись несущественно. Представленные факты дополнительно свидетельствуют в пользу адекватности предложенной мнкромоделн.
4. Гидравлические функции при пропитке
4.1. Капиллярная кривая. Пропитка является более сложным процессом, чем дренаж. Это связано с тем, что условия заполнения водой узлов при пропитке отличаются от условия проникновения в узел воздуха при дренаже. Различие связано с тем, что размер узла не меньше размера выходящих из него капилляров. Поэтому при дренаже воздух автоматически проникнет в узел, как только он окажется в одном из выходящих из него капилляров, а при пропитке узел может остаться незаполненным, даже если несколько выходящих из него капилляров (меньшего размера) заполнены водой. Это обстоятельство ясно осознано, начиная, по-видимому, с работ [8 10]. В это же время введено понятие кооперативного механизма заполнения узлов. Считается, что чем больше выходящих из узла капилляров заполнено водой, тем мениск меньшего радиуса способен проникнуть в узел, а, значит, тем при меньшем уровне давления в системе произойдёт заполнение узла. Формализуя этот факт, вводят ранжированную последовательность Pi < p2 < ... < Pz-1 пороговых давлений, при которых происходит заполнение узла при наличии Z — i заполненных капилляров, выходящих из пего (события I1, /2,..., Iz-1). Величина p1 обычно определяется [11] условием существования равновесного мениска в узле, p1 = — а_ст/Ь8не. В силу того, что размер LSjte узла в предлагаемой модели совпадает с максимальным размером выходящих из него капилляров, имеем
P1 = max{p-,.. .p-}, (17)
где через p- обозначены критические значения p- давления для выходящих из узла капилляров. Информация о последующих пороговых значениях давления носит менее точный характер. Предлагаемые в различных работах [7, 11, 13] формулы для их вычисления вряд ли можно назвать строго обоснованными.
p2 p1 поэтому во многих ситуациях события /2,/3,... разумно попросту игнорировать. Применяя, например, соотношение из [11] к случаю однородных (с порами одного
p2 p1 = p-
в 2.74 раза. Это значение заметно превышает критическое давление p+ для пор любой формы. Ясно, что в этом случае событие 12 попросту не может наступить, так как значительно раньше с ростом давления сработает "snap-off" механизм заполнения пор. Учитывая изложенное, ограничимся в дальнейшем рассмотрением лишь механизма 11 кооперативного заполнения узла.
Ещё одним обстоятельством, которое будет приближённо учтено в развиваемой нами микромодели, является наличие в среде защемлённого воздуха. В процессе пропитки пора, или некоторая система связанных пор (кластер), может оказаться окруженной со всех сторон водой. Заполнение водой такой поры (кластера) в дальнейшем оказывается невозможным. Полный учет эффекта защемления воздуха является сложной задачей, которую можно решить, видимо, лишь на основе непосредственного численного моделирования процесса пропитки на решётке капилляров. Соответствующие численные эксперименты [13. 14] показывают, однако, что для пористых сред с малой относительной дисперсией распределения пор по размерам воздушные кластеры невелики, а большинство из них состоит лишь из одной поры. Поэтому в первом приближении разумно учитывать лишь эту возможность защемления воздуха.
Дополнительный учет механизма /1 кооперативного заполнения и наличия защемлённого воздуха приводит к следующему выражению для зависимости п(р, Ь) доли полностью заполненных водой пор данного размера Ь от давления р
(1 — с22-2(р+(ь)), р > р+(Ь), п(р, Ь) = < 2£2_1(р)(1 — С^(р)), р- <р<р+, (18)
[о, р<р-(Ь).
Здесь С(р) по-прежнему обозначает общую долю заполненных пор при данном р
р<
< р - ( Ь )
ностыо заполнен. Средняя строка (с учётом (17)) указывает на то. что необходимым и достаточным условием заполнения выделенной поры в диапазоне давлений р € [р_,р+] является требование того, чтобы либо все 2 — 1 её «правых» соседей, либо все её 2 — 1 «левых» соседей (по ни те и ни другие одновременно) оказались заполненными водой. Под «правыми» и «левыми» соседями здесь понимаются капилляры, входящие в одни из двух узлов, на которые опирается выделенный капилляр. Указанное условие при этом есть ни что иное, как условие срабатывания /1
поры. Наконец, верхняя строка в (18) указывает на то, что
• с переходом давления через критическое значение р+(Ь) заполненным будут
Ь
непосредственными соседями:
• при р > р+ (Ь) пет механизмов, изменяющих состояние поры, в силу чего величина п будет оставаться неизменной с дальнейшим ростом давления.
Для нахождения £(р) по-прежнему можно воспользоваться соотношением (5), получив в результате интегрирования (18) уравнение следующего вида
Ф)= Г(1 — С2^_2Ы)/№(у) ф + 2С^_1(г)(1 — С^М) /2 /м(у) ф. (19) J0 иг
Через г по-прежнему обозначается безразмерное обратное давление. При численном решении удобно дифференцированием по г свести уравнение (19) к обыкновенному дифференциальному уравнению первого порядка. Результаты такого решения для 2 = 4 и значений параметра плотности распределения х = 9) 33, 56, соответствующих лабораторным пескам 12/20, 20/30, 30/40, представлены на рис. 6.
3|
Рис. 6. Зависимость доли полностью заполненных пор С °т безразмерного давления при различных значениях параметра х плотности распределения пор по размерам
Как видно из этого рисунка, рассматриваемая микромодель предсказывает наличие предельного значения < 1, а значит, и остаточной по воздуху насыщенности ва;г < 1. При достаточно больших значениях х (для достаточно однородных сред) значение Са1г зависит лишь от координационного числа 2 решётки. Это можно наблюдать на рис. 7 (для 2 = 4 с изменением х от 9 до 56 остаточное значение изменяется слабо, оставаясь равным приблизительно 0.89).
Обратимся к вычислению главной кривой пропитки в = Б(р). Будем игнорировать при этом для простоты плёночную составляющую насыщенности, что вполне допустимо, если только не работать в области сверхмалых влажпостей (в < 10-4). Разобьём при данном давлении р всё множество капилляров па три класса: заполненные водой, неблокированные капилляры с уголковым режимом заполнения
Б(р)
вкладов
всар = П(р (Ь) в = в Пте„(р,Ь)Б(1)(р,Ь)/„(Ь) ¿Ь. (21)
Здесь через Птеп обозначена числовая доля неблокированных пор с уголковым режимом заполнения. Она за вычетом доли блокированных пор пгар совпадает со значением 1 — ц\
в — Б (р) — всар + втеп + strap,
подсчитаем их независимо друг от друга.
Первые два вклада вычисляются аналогично случаю дренажа:
(20)
С2^-2(р+(ь)), р >р+(ь), С22-2(р), р < р+(Ь).
Вычисляя интегралы (21), получим
вСарМ= Г (1 — С22-2Ы)/у(у) ¿у + 2С^-1(^)(1 — С^ЖРу(а_г) — Ру(г)), (22) Jo
Рис. 7. Распределение по размерам остаточной доли пор с защемлённым воздухом (заштрихована) для двух различных пористых сред
3
smen(z) = Р ^ ^ z + / х k=1
х [(1 - CZ-1(z))2(Gi+k/2(z) - Gi+k/2(a-z)) + (1 - C2Z-2(z))Gi+k/2(a_z)]. (23)
Определённая сложность подсчета strap вызвана тем, что входящие в число %rap(p, L) капилляры были блокированы при различных уровнях давления p' < p,
а значит, содержат различное количество влаги S (1)(p' , L). Поэтому Strap зависит
p
Strap = /3^00(|P gW^^^'^'^d^/^L)^. (24)
Переставляя здесь пределы интегрирования и подсчитывая интеграл по L явно, получим
Strap(z) = /г z1+k/2Gi+k/2(z)dC2Z-2(z). (25)
k=i Jo
Результаты проведённых по формулам (20), (22), (23), (25) расчётов для грунта 20/30 представлены на рис. 5, а штрих-пунктирной линией. Параметры микромодели Lo, х а были определены ранее по данным дренажных испытаний.
В отличие от главной дренажной кривой (сплошная линия), кривая пропитки имеет характерной особенностью наличие остаточной по воздуху насыщенности sa¡r < 1. В рассматриваемом случае sa¡r = 0.85; для других пористых засыпок 12/20 и 30/40 значение sa¡r составляет 0.82 и 0.87 соответственно. Во всех случаях sa¡r несколько меньше, чем соответствующее значение Zair • Объёмная доля sa¡r пор с защемлённым воздухом оказывается меньше числовой доли в силу того, что воздух блокируется преимущественно в больших порах. Этот факт иллюстрируется рис. 7, заштрихованная часть которого отвечает блокированным порам на момент
s
вестной (см., например, [25, 26]) экспериментальной вилке 0.8 < sa¡r < 0.9.
Как видно из рис. 5, а, в области больших и умеренных насыщенностей имеет место примерно двукратное отлично капиллярных кривых пропитки и дренажа, сохраняющееся и для других типов засыпок. Оно объясняется примерно двукратным
отличием в значениях р+, р_ критических давлений и находит своё подтверждение в экспериментах [25].
Обратим внимание, что ярко выраженный при умеренных в гистерезис капиллярных кривых становится пренебрежимо малым в области малых насыщенно-стей, где обе кривые оказывается возможным описать единой зависимостью. Этот теоретический результат приобретает особое значение в связи с отсутствием достоверных экспериментальных данных по измерению капиллярных кривых при малых в [27-29].
4.2. Относительная фазовая проницаемость при пропитке. После подсчёта доли п(р, Ь) заполненных пор относительная фазовая проницаемость при пропитке определялась подобно тому, как это делалось для случая дренажа. В конечном итоге приходилось решать аналогичное (14) алгебраическое уравнение.
Результаты решения для среды 20/30 представлены штриховой линией на рис. 5, Ь. Как видно, характер зависимости К (в) при пропитке аналогичен наблюдаемому при дренаже (сплошная линия). Отличие проявляется в области малых насыщенностей. Оно связано главным образом с иным (большим) значением остаточной насыщенности, > в-. В случае пропитки значение оказывается возможным подсчитать по формуле
аналогичной формуле (16). Из сравнения (26) и (16) следует, что = Простые оценки показывают, что капиллярная составляющая всар предельной насыщенности практически одна и та же в обоих случаях. Количество же втеп жидкости, содержащейся при г = в уголковых менисках, при пропитке заметно больше, чем при дренаже. Причина - в сильной зависимости втеп от г и существенной (в а_ « 2 раза) разнице между и 8.
Как правило, в физике почв игнорируют различие между кривыми фазовой проницаемости для дренажа и пропитки. Представленные на рис. 5, Ь результаты в целом оправдывают это упрощение вне области малых насыщенностей. Принципиально важным для ряда теорий (см., например, [30]) является, однако, то, что гистерезис всё-таки есть, и кривая фазовой проницаемости при пропитке лежит ниже соответствующей дренажной кривой. Именно такое поведение фазовых пронпцаемостей наблюдается и в тщательно проведённых лабораторных экспериментах [25, 31].
Первичная пропитка характеризуется высоким темпом роста влажности изначально сухой пористой среды. В то же время увеличение плёночной и уголковой составляющих насыщенности, а также механизм "snap-off" заполнения пор являются медленными процессами. Скорость их протекания лимитируется большим гидравлическим сопротивлением плёнок и менисков по сравнению с гидравлическим сопротивлением заполненных пор. Естественно поэтому при описании первичной пропитки пренебрегать указанными медленными механизмами повышения влажности, считая, что они не успевают проявиться в течение изучаемого процесса в силу его высокого темпа. Единственными учитываемыми механизмами заполнения капилляров н узлов решётки при этом остаются поршневое вытеснение в капиллярах и кооперативное заполнение узлов.
В отличие от предыдущего, описание процесса при этом удобнее проводить, рассматривая в качестве первичного элемента узел, а не капилляр. Функция Fs(L)
(26)
5. Кривая первичной пропитки
распределения узлов по размерам строится, исходя из аналогичной функции .Р„(Ь) для капилляров:
^(Ь) = ^ (Ь). (27)
Подобно тому, как это было сделано ранее для капилляров, введём в рассмотрение функции С(р) и п(р, Ь) как численные доли заполненных узлов при данном уровне давления среди всех узлов образца и среди узлов размера Ь соответственно. Функции п(р, Ь) и п(р, Ь) можно связать друг с другом следующими простыми
Ь
Ь ^ Ь и Ь2 ^ Ь. При этом
|о, Ь' <Ь,
вероятность {Ьк < Ь'} = < к = 1,2. (28)
(/„^(Ь'), Ь' > Ь,
В силу того, что заполнение узла эквивалентно заполнению всех поросокаю-
Ь
тогда, когда хотя бы один из узлов Ьх, Ь2 заполнится. Отсюда, с учётом (28) следует, что
с ж
П(Р,Ь) = 1 - (1 - х)2, х = ^-1(Ь)П(Р,Ь) + ^ П(Р,Ь') ^-1(Ь'). (29)
В свою очередь, функция п(р, Ь) определяет, как и прежде (см. (7)), насыщенность образца при данном уровне давления
8 = 5(р)= / п(р,Ь)/«(Ь) ¿Ь. (30)
ио
В формуле (30) по сравнению с (7) иренебрежено дополнительно плёночной и уголковой составляющими насыщенности.
Таким образом, для вычисления кривой первичной пропитки по (29), (30) остаётся определить п(р, Ь). С этой целью обозначим через
л(С) = с|(1- скг=
k!(Z - k)!
вероятность того, что 2 — к узлов, соседствующих с данным, заполнены (событие /к). Представ им п(р, Ь) в виде
L) = <
/i + /2 + ... + /z-i, P>PZ-i,
/i + /2 + ... + /k, Pk <p<pk+i,
ib Pi <P <P2,
0, P<Pi,
(31)
обозначив pk = —«k^/L. Воспользовавшись далее тем, что
POO
C(p)=/ n(P,L)dFf (L), (32)
/0
Ь
£(р). В прежней безразмерной переменной г = —ст/рЬ0 оно записывается как
z-i
С=Е /k (OFz К z). (33)
k=i
Рис. 8. Зависимость доли заполненных узлов £ от безразмерного давления z
Можно показать, что при любом значении Z и произвольном выборе ai > > > ... > az-i > 0 решения этого уравнения устроены одинаково. Помимо тривиального решения £ = 0, имеющегося при любом значении z, уравнение (33) обладает ещё одним решением при z > z* и двумя решениями па интервале z ££ (z**, z*).
Величину z* можно найти, устремив в (33) £ к нулю. Учитывая, что fz-i =
= Z£ + O(C2) и fk = O( С2) ^и k < Z — 1, получим следующее условие для z*
FZ (az-iz*) = 1/Z.
С* z
чипы F^(akz) в правой части (33) обращаются в единицу, сама правая часть после суммирования оказывается равной 1 — £Z — (1 — c)z. В результате С* находится как решение уравнения
С* = 1 — CZ — (1 — C*)z.
С* z*
ak
Типичные результаты расчета функции £(z) представлены па рис. 8. Координа-Z ak
ai = 1.7, a2 = 1.15, аз = 0.7.
Насыщенность s выражается через решение £(z) уравнения (33) по формулам (29) (31).
С
интересным фактом. Она приводит к неоднозначной зависимости насыщенности от давления, или, что то же самое, к немонотонности кривой p = P(s) первичной пропитки. Соответствующая кривая, построенная для среды 20/30, представлена на рис. 5, а штриховой линией. Очевидно, что возрастающая часть этой кривой неустойчива и не может быть физически реализована. Переход из сухого во влажное состояние может осуществиться лишь скачком, переводящим изиачальио сухую среду в состояние, близкое к полному насыщению. Насыщенность при этом вычисляется явно через предельное значение С* как s* = 1 — (1 — £*)2. При Z = 4 имеем s* = 0.923, а с изменением координационного числа величина s* меняется несущественно; так при Z = 5 и 6 значение s* равно 0.94 и 0.95 соответственно. Во
s* s
главной кривой пропитки.
Характерная полка давления на кривой первичной пропитки, соединяющая состояние нулевой насыщенности с близким к полному насыщению состоянием отмечается практически всеми исследователями, работающими с лабораторными песками. Развитая здесь теория даёт этому факту физическое объяснение.
Авторы благодарны профессору Р. ДиКарло за представленный экспериментальный материал и профессору Дж.Л. Нибору за плодотворные обсуждения результатов работы.
Работа выполнена при финансовой поддержке РФФИ (проекты Х- 04-01-00484. 05-01-00516).
Summary
B.E. Gogolashvili, A.G. Egorov. Application of effective medium theory for obtaining hydraulic functions on pore networks.
Effective medium theory on pore networks was employed for modeling hydraulic characteristics of porous media. Different flow regimes (flow in completely filled pores, corner and film flow in partially filled pores) and different mechanisms of liquid displacement (piston-type displacement, "snap-off", cooperative wetting) was taken into account in modeling of single pore behavior at microlevel. Proposed micromodel was then upscaled to obtain macroscopic hydraulic characteristics for typical laboratory silica sands. All calculated hydraulic functions (drainage curve, hydraulic conductivity, main wetting curve, and primary wetting curve) showed good agreement with experimental data.
Литература
1. Dautov R.Z., Egorov A.G., Nieber J.L., Sheshukov A.Y. Simulation of two-dimensional gravity-driven unstable flow // Proc. 14t.li Int. Conf. on Сотр. Motli. in Wat. Res. Delft. The Netherlands, 2002. V. 1. P. 9 16.
2. Egorov A.G., Dautov R.Z., Nieber J.L., Sheshukov A.Y. Stability analysis of traveling wave solution for gravity-driven flow // Proc. 14t.li Int. Conf. on Сотр. Motli. in Wat. Res. Delft, The Netherlands, 2002. V. 1. P. 120 127.
3. Егоров А.Г., Даутоо Р.З. Моделирование процессов пальцеобразовапия при влагопе-репосе в ненасыщенных пористых средах // Тр. Матом, центра им. Н. И. Лобачевского. Т. 15. Модели механики сплошной среды. Обзорные доклады и лекции XVI сессии Международной школы по моделям механики сплошной среды, Казань, 27 июня
3 июля 2002 г. Казань: Изд-во Казан, матем. о-ва, 2002. С. 28 60.
4. Nieber J.L., Dautov R., Egorov A., Shehsukov A. Dynamic capillary pressure mechanism for gravity-driven flows: review and extension to dry conditions // Transport Porous Media. 2005. V. 58. P.147 172.
5. Cuesta C., van Duijn C.J., Hulshof J. Infiltration in porous media with dynamic capillary pressure: travelling waves // Europ. J. Appl. Math. 2000. V. 11. P. 381 397.
6. Гоголашаили Б.Э. Определение равновесных гидравлических функций па основе микромоделировапия пористой среды // Тр. Матем. центра им. Н. И. Лобачевского. Т. 27. Модели механики сплошной среды. Материалы XVII сессии Межд. школы по моделям механики сплошной среды, Казань, 4 10 июля 2004 г. Казань: Изд-во Казап. матем. о-ва, 2004. С. 88 95.
7. Вакке S., 0геп Р. 3-D pore scale modelling of heterogeneous sandstone reservoir rocks and quantitative analysis of the architecture, geometry and spatial continuity of the pore network // Soc. Pet. Eng. 1996. SPE paper 35479.
8. Lenormantl R., Zareone С., Sarr A. Mechanisms of the displacement, of one fluid by another in a network of capillary ducts // J. Fluid Mecli. 1983. V. 135. P. 337 353.
9. Lenormantl R., Zareone C. Role of roughness and edges during imbibition in square capillaries // Proc. of 59t.li SPE Annual Technical Conference and Exhibition. Houston. TX. 1984. SPE paper 13264.
10. Lenormantl R. Liquids in porous media // J. Pliys.: Condens. Matter. 1990. V. 2. P. 79 88.
11. Patzek T. W. Verification of a com pi et e pore network model of drainage and imbibition // Soc. Pet. Eng. 2000. SPE paper 59312.
12. Hughes R.G., Blunt M.J. Pore scale modeling of rate effects in imbibition // Transport Porous Media. 2000. V. 40(3). P. 295 322.
13. Blunt M.J., Seher H. Pore-level modeling of wetting // Pliys. Rev. E. 1995. V. 52(6). P. 6387 6403.
14. Панфилов М.Б., Панфилова И.В. Осредпёппые модели фильтрационных процессов с неоднородной внутренней структурой. М.: Наука, 1996. 383 с.
15. Kirkpatriek S. Percolation and conduction // Rev. Mod. Pliys. 1973. V. 45. P. 574 588.
16. Yanuka M., Dullien F., Eniek D. Percolation processes and porous media // J. Colloid Interface Sci. 1986. V. 112. P. 24 41.
17. Brooks R.H., Corey A.T. Hydraulic properties of porous media // Hydrol. Pap. 3. Colorado State Univ., Fort Collins. 1964.
18. van Genuehten M.T. A closed form equation for predicting the hydraulic conductivity of unsaturated soils // Soil Sci. Soc. Am. J. 1980. V. 44. P. 892 898.
19. Tuller M., Or D. Unsaturated hydraulic conductivity of structured porous media: A review of liquid configuration-based models // Vadose Zone J. 2002. P. 14 37.
20. Mason G., Morrow N.R. Capillary behavior of a perfectly wetting liquid in irregular triangular tubes // J. Colloid Interface Sci. 1991. V. 141. P. 262 274.
21. Iwamatsu M., Horii K. Capillary condensation and adhesion of two wetter surfaces // J. Colloid Interface Sci. 1996. V. 182. P. 400 406.
22. Alt W., Luekhaus S., Visintin A. On nonstationary flow through porous media // Ann. Math. Рига Appl. J. 1984. V. 136. P. 303 316.
23. Sehroth M.H., Ahearn S.J., Selker J.S., Istok J.D. Characterization of Miller-similar silica sands for laboratory hydraulic studies // Soil Sci. Soc. Am. J. 1996. V. 60. P. 1331-1339.
24. DiCarlo R. Experimental measurements of saturation overshoot on infiltration // Water Resour. Res. 2004. V. 40(4).
25. Ustohal P., Stauffer F., Draeos T. Measurement and modeling of hydraulic characteristics of unsaturated porous media with mixed wettability // J. Cont. Hydrol. 1998. V. 33. P. 5 37.
26. Stauffer F., Kinzelbaeh W. Cyclic hysteretic flow in porous medium column: model, experiment, and simulations // J. Hydrology. 2001. V. 240. P. 264 275.
27. Ross P.J., Williams J., Bristow K.L. Equation for extending water-retention curves to dryness // Soil. Sci. Soc. Am. J. 1991. V. 55. P. 923 927.
28. Morel-Seytoux H.J., Nimmo J.R. Soil water retention and maximum capillary drive from saturation to oven dryness // Water Resour Res. 1999. V. 35(7). P. 2031 2041.
29. Toledo P.G., Novy R.A., Davis Н.Т., Seriven L.E. On the transport properties of porous media at low water content., in Indirect methods for estimating the hydraulic properties of unsaturated soils. U.S. Salinity Laboratory, Riverside, CA, 1992. P. 97 114.
30. Kaeimov A.R.,Yakimov N.D. Nonmonotonic moisture profile as a solution of Richards' equation for soils with conductivity hysteresis // Adv. Water Resour. 1998. V. 21(8). P. 691 696.
31. Wong P-Z. Flow in porous media: permeability and displacement patterns // MRS Bulletin. 1994. V. 19(5). P. 32 38.
Поступила в редакцию 12.09.05
Гоголашвили Вулат Эдуардович младший научный сотрудник НИИ математики и механики им. И.Г. Чеботарева Казанского государственного университета.
E-mail: Bulat. GogolashviliQksu.ru
Егоров Андрей Геннадьевич доктор физико-математических паук, старший научный сотрудник, заведующий кафедрой аэрогидромехапики Казанского государственного университета.
E-mail: Andrey.EgorovQksu.ru
