Математическая модель промерзания ненасыщенного грунта при наличии капиллярных сил

Математические заметки СВФУ Апрель—июнь, 2017. Том 24, № 2

УДК 532.546

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПРОМЕРЗАНИЯ НЕНАСЫЩЕННОГО ГРУНТА ПРИ НАЛИЧИИ КАПИЛЛЯРНЫХ СИЛ Г. Г. Цыпкин

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

Б01 10.25587/8УРи.2017.2.9249

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

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

д'д _ д дЬ дх

Здесь •& — влажность грунта или масса воды, отнесенная к единице массы сухого образца, О — коэффициент диффузии влажности, определяемый из эксперимента. Подход, основанный на эффективных уравнениях, широко распространен и использовался в большом числе работ (например, [4-6]).

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

Работа выполнена при финансовой поддержке гранта РФФИ № 16—01—00363.

© 2017 Цыпкин Г. Г.

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

В [9] сформулирована задача о плавлении льда в мерзлом грунте при наличии капиллярного давления в области вода — воздух. Выведена полная система граничных условий на фронте фазового перехода. В [10] рассмотрена задача о плавлении льда в грунте, изначально насыщенном льдом и воздухом. В нелинейном случае проведено численное исследование зависимости водонасыщенно-сти от вида функции Леверетта, величины капиллярного давления и внешних градиентов давления и температуры.

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

1. Основные уравнения. В положительном диапазоне температур будем моделировать грунт недеформируемой пористой средой, насыщенной гетерогенной смесью воды и воздуха. Предположим равенство температур компонент в областях их совместного существования. С учетом сделанных предположений уравнения неразрывности для газа и воды, обобщенный закон Дарси для каждой из компонент и уравнение состояния газа имеют вид [10]

= 0,

д

=

= АР3, з=9,™,

Рд = Рд ЯГ.

Здесь Р — давление, Т — температура, К — газовая постоянная, V — вектор скорости фильтрации, р — плотность, Б — насыщенность, к — проницаемость, ф — пористость, р — вязкость, / — относительная фазовая проницаемость. Индексы: д — газ (воздух), ад — вода.

Рассмотрим влияние капиллярного давления на движение смеси газ — вода в пористой среде. Капиллярное давление по определению равно разности давлений между газом и водой:

РС = Рщ - Рд .

В случае смачиваемой среды Рс < 0.

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

, ВТ дРд

+

Исключая плотность газа и скорости фильтрации воды и газа, в области воздух — вода получаем систему четырех уравнений для неизвестных Бщ, Т, Рщ и Рд:

дБ к

сИУ[/ш

д£ фр.

' Р,

д^ (л я \ 1 дРа ^ (л я \ 1 дТ ± к Т А-

^"аГ +(1 - ^тЖ + лгдтд а1У

Т

= 0,

дТ дБ д Р

к

г, /ш (Бщ) , , СР Рд /д (Бщ) , „

о-шС-и,-grad Рги + -—-ф —-grad Рд

Рш К Т Рд

grad Т

к2 рш

Рш Рд + Рс ( Бш ) 1

Ат = фБщАщ + ф(1 - Бш)Ад + (1 - ф)А«, (РС )т = фБщ Рш Сщ + ф(1 - Бд )Рд Ср + (1 - ф)рз С«.

В пористой среде, насыщенной гетерогенной смесью воды и газа, капиллярное давление зависит от водонасыщенности и определяется соотношением Леверетта [11]:

Рс{Бш)=х\1|</

Параметр х зависит от смачиваемости пород, 3(Бщ) — функция Леверетта.

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

Рс = —2— cos в.

r

Здесь a — коэффициент поверхностного натяжения, в — контактный угол, r — радиус кривизны капиллярного мениска.

Величина среднего капиллярного давления в пористой среде следует из оценки характерного размера (диаметра) поры: d = \Jф/k. Радиус капиллярного мениска совпадает с радиусом поры и равен половине диаметра r = d/2.

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

Рс = -Act cos eJj-J

V k

Поскольку измерение контактного угла представляет трудности, зачастую при расчетах в качестве параметра задается величина капиллярного давления или значение параметра х, который, как следует из (1.10), х = -4a cos в. Для смачиваемой среды х < 0 [12,13].

Положим, что капиллярное давление по величине не превосходит атмосферного давления и давления в жидкости и газе являются величинами одного порядка: Pw ~ Pg Pc |.

2. Формулировка задачи. Упростим задачу, используя существенную разницу физических свойств воды и воздуха. Поскольку вязкость газа на два порядка меньше вязкости воды, то скорость фильтрации газа значительно выше, если различие в фазовых проницаемостях невелико. Тогда можно положить, что давление в газе устанавливается мгновенно по отношению к характерному времени перераспределения функции водонасыщенности и в первом приближении считать его постоянным: Pg = const. В этом случае из соотношения (1.5) следует grad Pw = grad Pc и скорость перераспределения давления в жидкости совпадает со скоростью перераспределения капиллярного давления, т. е. функции водонасыщенности. Тогда функция водонасыщенности удовлетворяет нелинейному уравнению диффузии вида

dSw k п.

div

dPc

fw(Sw)-^- gradSV, dSw

0,

dt фj

а уравнение энергии принимает вид

dT dS

= -pwCw^^^gra,dSw-gra,dT-—fw

fj-w dSw fw \dSU

Уравнение (2.1) можно записать в виде

и.т

дБ«,

шу

дЬ

\ ф йБш

Уравнение (2.3) по форме совпадает с уравнением (1), если положить, что эффективный коэффициент влагопереноса определяется выражением

/ к /ш

к = Х\Пк--Тя-•

V ф йош

Поскольку производная функции Леверетта всегда отрицательна: " < 0, то для смачиваемой пористой среды (х < 0) коэффициент диффузии влаги к всегда положителен. В случае несмачиваемой среды, когда х > 0, коэффициент диффузии отрицателен и уравнение переноса влаги сводится к некорректному уравнению диффузии.

Уравнение диффузии влажности (2.3) для движения двухкомпонентной смеси вода — воздух, выведенное из закона сохранения массы, закона Дарси и выражения для функции капиллярного давления, аналогично эффективному уравнению диффузии влаги. Уравнение (2.3) получено в предположении, что градиент давления в газовой компоненте мал и относительные фазовые проницаемости различаются незначительно. Если учесть движение воздуха, то в этом случае задача сводится к системе двух нелинейных уравнений, и уже одно уравнение диффузии не будет адекватно описывать процесс.

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

§ = а1АТ< а1 =

А/ = фБгАг + ф(1 - Бг)\д + (1 - ф)А*, (рС)/ = фБгргсг + ф(1 - Бг)рдСр + (1 - ф)р,С8,

В талой области перед фронтом кристаллизации перераспределение воды описывается уравнением диффузии влаги (2.3). Рассмотрим задачу в линейном приближении, когда отклонение функции водонасыщенности от невозмущенного состояния невелико. Невозмущенным является состояние, когда вода неподвижна и начальная водонасыщенность равна Бшо, а для невозмущенного состояния за фронтом насыщенность льда определяется выражением р¿о = Бшорш/р%.

Для простоты аппроксимируем функцию относительной фазовой проницаемости и функцию Леверетта линейными функциями: /ш — Бш и 3(Бш) — 1 Бш. Тогда в области смеси вода — воздух уравнение диффузии влаги принимает вид

дБш кБшо <1РС / к 5Шо г \

-— = кАБ-и,, к=~--тег = 4а сов 6* --.

дЬ фрш <БШ у ф

Если капиллярное давление меньше атмосферного: |РС| < Ра — 105 Па, то второе слагаемое в левой части уравнения (2.2) много меньше первого слагаемого и в линейном приближении в области вода — воздух уравнение энергии сводится к обычному уравнению теплопроводности

^ = а2АТ, а2= ,

дЬ (рС)т

Ат — фБшоХш + ф(1 - Бшо)\д + (1 - Ф)Аа, (рС)т — фБшОРшСш + ф(1 - Бшо)РдСр + (1 - ф) р3С3.

Граничные условия на подвижной поверхности образования льда состоят из условий сохранения энергии и массы Н20 при известной температуре кристаллизации Т* — 273.15 К в предположении, что капиллярные силы не изменяют температуру фазового перехода.

Уравнение баланса энергии имеет вид

фБlplqVn — Аf

Здесь теплота фазового перехода лед — вода q — Ьш - Ь,^ а V — Х(Ь) — скорость фронта плавления. Индексы плюс и минус — значения величин соответственно перед и за фронтом кристаллизации, п — нормаль.

Из баланса массы Н20 следует соотношение

--) Уп = --Ри))п+

Рш ) ФРш

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

= (2.8) Рш ) Бш0

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

Сформулированная задача содержит три неизвестных параметра: скорость подвижной границы фазового перехода V, насыщенность льда и значение водонасыщенности Бш+ на поверхности раздела. Для их определения имеется два соотношения на границе кристаллизации — (2.7), (2.8), и задача является недоопределенной. Для замыкания системы соотношений на подвижной границе требуется выставить дополнительное условие, связывающее искомые параметры. Естественно предположить, что весь растущий кристалл льда должен

омываться водой. В той части, где кристалл не контактирует с водой, а контактирует с воздухом, лед формироваться не будет. Из этих соображений следует дополнительное условие равенства насыщенностей воды и льда Б^- = Бш+ = , которое замыкает систему соотношений на поверхности кристаллизации.

3. Автомодельное решение. Рассмотрим задачу в одномерном приближении. Пусть в начальном состоянии талый грунт занимает полупространство х > 0 с начальной водонасыщенностью Бшо. Если температура на неподвижной стенке х = 0 опускается ниже температуры фазового перехода вода — лед, то вправо распространяется фронт кристаллизации, двигающийся по закону х = X(£). Если начальная температура и водонасыщенность, а также температура охлаждающей стенки являются постоянными величинами, то сформулированная задача имеет автомодельное решение вида

Т = Т

С = х/уД.

В области перед фронтом решения уравнения теплопроводности и диффузии водонасыщенности имеют вид

> X

Бw — Бш0 +

ет{с(0/2^/01)'

ег&(/3/2Л/к)'

В области за фронтом распределение температуры и водонасыщенности определяется выражениями

0 < х < X

ей

Подставляя решения (3.1), (3.2) в соотношения (2.7) и (2.8), получаем трансцендентные уравнения для определения безразмерной автомодельной скорости фронта образования льда 7 = ¡3/2^/а{ (£ = £/2^/07) и значения насыщенности Б :

^/ттфБ^РгЧах /оГ АТ1* 1

-г-7- Л/\/— тр - 1

Т0 V а-2 \Т0

1

ехр(-72а!/а2) ег£(7)ч/а1/а2 ' Т*

— Ат

Т0 Т0

Т0^ ехр

Б* ~ Б° ехр

= 0.

= 0,

ров:

Система (3.4), (3.5) решалась численно при следующих значениях парамет-

ф = 0.25, ра = 2 • 103 кг/м3, рг = 0.91 • 103 кг/м3 рш = 103 кг/м3, Аз = 2Вт/(м • К), Аг = 2.23 Вт/(м • К), Аш = 0.58Вт/(м • К), Сз = 103 Дж/(кг • К), С = 2 • 103 Дж/(кг • К), Сш =4.2 • 103 Дж/(кг • К), д = 3.34 • 105 Дж/кг, = 2 • 10-3 Па • с.

х

7

0.54 р 1

0.49 -1-1-1-1

0 0.5 1 1.5 £ 2

Рис. 1. Характерные распределения водонасыщенности воды и льда в зависимости от величины капиллярного давления. То = 280 К, Т0 = 260 К, £—о = 0.5.

Кривые 1 и 2 соответствуют коэффициентам диффузии к = 10-7 и к = 10-6

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

Поскольку при кристаллизации объем Н2 С возрастает, то значение водо-нащенности на фронте должно быть больше, чем начальное значение водонасыщенности. Отсюда следует, что согласно диффузионному уравнению (2.5) жидкая фаза должна оттекать от фронта и наибольшая водонасыщенность не должна превосходить невозмущенного значения насыщенности льда Б го, которое определяется из условия неподвижности воды. В противном случае должен реализовываться подток влаги к фронту, что не обеспечивается уравнением (2.5) при условии максимального значения на фронте кристаллизации.

Численные эксперименты показали, что представленная модель адекватна и удовлетворяет этим требованиям, а водонасыщенность , равная насыщенности льда Б г, изменяется в пределах от начального значения Бы0 до значения Бго = Бы0/рг. При больших значениях коэффициента к, соответствующего высокой проницаемости грунта, характерное время перераспределения воды перед фронтом существенно меньше характерного времени перераспределения температуры и водонасыщенность на фронте лишь немного превышает начальную водонасыщенность. Так, например, при к = 105 м2/с будет Б* = 0.509 при Бы0 = 0.5. При уменьшении коэффициента диффузии влажности отток от фронта замедляется и значение водонасыщенности на фронте возрастает (рис. 1), приближаясь к невозмущенному значению объемной доли льда.

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

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

0.55 |-

S

I

0.54 -

0.53 -

0.52 -

0.51 -

0.5 -

220 230 240 250 260 270 ^ 280

Рис. 2. Влияние температуры охлаждающей стенки на количество образовавшегося льда: To = 280 К, Swо = 0.5. Кривые 1 и 2 соответствуют коэффициентам диффузии к = 10-7 и к = 10-6

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

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

ЛИТЕРАТУРА

1. Цытович Н. А. Механика мерзлых грунтов. М.: Высшая школа, 1973.

2. Жесткова Т. Н. Формирование криогенного строения грунтов. М.: Наука, 1982.

3. Hoekstra P. Moisture movement in soils under temperature gradient with the cold-side temperature below freezing // Water Resources Res. 1966. V. 2. N 2. P. 241-250.

4. Zhang Y., Carey S. K., Quinton W. L. Evaluation of the algorithms and parametrizations for ground thawing and freezing simulation in permafrost regions //J. Geophys. Res. 2008. V. 113. D17116.

5. Li Q., Sun S., Dai Q. The numerical scheme development of a simplified frozen soil model // Adv. Apmosph. Sci. 2009. V. 26. N 5. P. 940-950.

6. Akbari G., Tabrizi H. B., Damangir E. Numerical and experimental investigation of variable phase transformation number effect in porous media during freezing process // Heat Mass Transfer. 2009. V. 45. N 4. P. 407-416.

7. Menot J. M. Equations of frost propagation in unsaturated porous media // Eng. Geology. 1979. V. 13. N 1-4. P. 101-109.

8. Цыпкин Г. Г. Линейная задача о фазовых переходах вода — лед в ненасыщенных грунтах // Изв. АН СССР. Механика жидкости и газа. 1990. № 3. С. 68-73.

9. Цыпкин Г. Г. Влияние капиллярных сил на распределение влагонасыщенности при про-таивании мерзлого грунта // Изв. РАН. Механика жидкости и газа. 2010. № 6. С. 130-140

10. Васильев В. И., Попов В. В., Цыпкин Г. Г. Нелинейная задача протаивания ненасыщенного мерзлого грунта при наличии капиллярных сил // Изв. РАН. Механика жидкости и газа. 2012. № 1. С. 132-141.

11. Leverett M. C. Capillary behavior in porous solids // Trans. AIME. 1941. V. 142. N 1. P. 152-169.

12. Bear J. Dynamics of Fluids in Porous Media. New York: Acad. Press, 1972.

13. Pruess K., O'Sullivan M. Effect of capillarity and vapor adsorbtion in the depletion of vapor-dominated geothermal reservoirs // Proc., 17th Workshop on Geothermal Reservoir Engineering. Stanford, CA. USA: 1992. P. 165-174.

Статья поступила 12 марта 2017 г. Цыпкин Георгий Геннадьевич

Институт проблем механики им. А. Ю. Ишлинского РАН, пр. Вернадского, 101-1, Москва 119526 tsypkin@ipmnet.ru

Математические заметки СВФУ Апрель—июнь, 2017. Том 24, № 2

ypK 532.546

A MATHEMATICAL MODEL OF FREEZING OF UNSATURATED SOILS IN THE PRESENCE OF CAPILLARY PRESSURE G. G. Tsypkin

Abstract. A mathematical model of freezing of soils saturated with a heterogeneous mixture of water and air in the presence of capillary pressure is proposed. The derivation of the diffusion equation for the redistribution of moisture from the laws of conservation of mass and momentum is given. This allows us to define the diffusion coefficient through the parameters of the porous medium and fluids. Balance relations through the water crystallization front is derived. A self-similar solution of the problem in the linear approximation is obtained. It is shown that the growth of capillary forces reduces the amount of ice formed, and a more intensive freezing regime leads to an increase in the saturation with ice.

DOI 10.25587/SVFU.2017.2.9249 Keywords: unsaturated soil, capillary pressure, filtration, freezing, ice.

REFERENCES

1. Tsytovich N. A., The Mechanics of Frozen Ground, McGraw-Hill, New York. 1975.

2. Zhestkova T. N., Formation of the Cryogenic Structure of Soils [in Russian], Nauka, Moscow (1982).

3. Hoekstra P., "Moisture movement in soils under temperature gradient with the cold-side temperature below freezing," Water Resources Res., 2, No. 2, 241—250 (1966).

4. Zhang Y., Carey S. K., and Quinton W. L., 11Evaluation of the algorithms and parametrizations for ground thawing and freezing simulation in permafrost regions," J. Geophys. Research, 113, D17116 (2008).

5. Li Q., Sun S., and Dai Q., "The numerical scheme development of a simplified frozen soil model," Adv. Apmosph. Sci., 26, No. 5, 940-950 (2009).

6. Akbari G., Tabrizi H. B., and Damangir E., "Numerical and experimental investigation of variable phase transformation number effect in porous media during freezing process," Heat Mass Transfer, 45, No. 4, 407-416 (2009).

7. Menot J. M., "Equations of frost propagation in unsaturated porous media," Eng. Geology, 13, No. 1-4, 101-109 (1979).

8. Tsypkin G. G., "Linear problem of water-ice phase transitions in unsaturated soils," Fluid dynamics, 25, No. 3, 384-389 (1990).

9. Tsypkin G. G., "Effect of the capillary forces on the moisture saturation distribution during the thawing of a frozen soil," Fluid dynamics, 45, No. 6, 942-951 (2010).

10. Vasil'ev V. I., Popov V. V., Tsypkin G. G., "Nonlinear problem of unsaturated frozen soil thawing in the presence of capillary forces," Fluid dynamics, 47, No. 1, 106-113 (2012).

The work was supported by the grant from the Russian Foundation for Basic Research (project code 16-01-00363).

© 2017 G. G. Tsypkin Г. Г.

11. Leverett M. C., "Capillary behavior in porous solids," Trans. AIME, 142, No. 1, 152-169 (1941).

12. Bear J., Dynamics of Fluids in Porous Media, Acad. Press, New York (1972).

13. Pruess K., O'Sullivan M., "Effect of capillarity and vapor adsorbtion in the depletion of vapor-dominated geothermal reservoirs," Proc., 17th Workshop on Geothermal Reservoir Engineering. Stanford, CA. USA, 165-174 (1992).

Submitted March 12, 2017. Georgii G. Tsypkin

Ishlinsky Institute for Problems in Mechanics of the Russian Academy of Sciences, pr. Vernadskogo 101-1, Moscow 119526, Russia tsypkin@ipmnet.ru