Задача о возмущениях фазового фронта в ненасыщенном грунте под действием инфильтрации осадков Текст научной статьи по специальности «Физика»

УДК 517

Задача о возмущениях фазового фронта в ненасыщенном грунте под действием инфильтрации осадков

А.Г. Петрова1, Н.П. Мошкин2, А.Ф. Жирков3

1 Алтайский государственный университет (Барнаул, Россия)

2 Институт гидродинамики им. М. А. Лаврентьева Сибирского отделения Российской академии наук (Новосибирск, Россия)

3 Институт мерзлотоведения им. П. И. Мельникова Сибирского отделения Российской академии наук (Якутск, Россия)

Phase Front Perturbations in Unsaturated Soils under the Rainwater Infiltration

A.G. Petrova1, N.P. Moshkin2, A.F. Zhirkov3

1 Altai State University (Barnaul, Russia)

2 Lavrentyev Institute of Hydrodynamics, Siberian Branch of the Russian Academy of Sciences (Novosibirsk, Russia)

3 Melnikov Permafrost Institute, Siberian Branch of the Russian Academy of Sciences (Yakutsk, Russia)

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

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

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

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

The paper is devoted to the study of rainwater infiltration influence on temperature, water saturation, and permafrost thawing. In general case, the problem is very complicated due to a lot of factors to be taken into account. We consider here a one-dimensional problem without including evaporation, vegetative layer, and radiation.

The domain between a daylight surface and a thawing front is filled with soil that is considered as porous media with immovable skeleton and pores containing water and air. The domain ahead of phase change boundary is occupied with frozen soil with pores containing ice and air. The model is constructed under the following assumptions: water and ice are incompressible, air is a perfect viscous gas; temperature and pressure are common for skeleton and pores; the soil surface is under the action of falling rain with given temperature and rate.

The system of governing equations in frozen and thawing soil consists of mass and energy conservation laws completed with generalized Darcy's law and the Clapeyron equation. The conditions on phase change boundary that coincides with the line of discontinuity of water saturation are obtained in a standard way from the conservation laws under the assumption of temperature and pressure continuity. The velocity of rain on the ground surface is assumed to be given, and that leads to assuming the condition for the pressure. Nonlinear system of equation with the free boundary is considered to be linearized on the stationary solution of the problem with the absence of rain.

The numerical simulation algorithm for a linear problem is elaborated, and the numerical results are provided.

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

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

БОТ 10.14258/1гуа8и(2015)1.1-18

Key words: phase change, unsaturated soils, rainwater infiltration, Stefan problem.

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

Решаем одномерную задачу. Ось 0Х направлена вниз. Область 0<х<%(() занята грунтом, который рассматривается как пористая среда с неподвижным скелетом и порами, заполненными водой и воздухом. Область % (I)<х<У занята мерзлым грунтом, в порах находятся лед и воздух.

Модель строится в следующих предположениях: вода и лед несжимаемые, воздух — вязкий совершенный газ; температура и давление — общие для скелета и пор; поверхность грунта подвержена воздействию выпадающего с определенной скоростью и температурой дождя. Силой тяжести пренебрегаем.

Определяющие уравнения в мерзлом и талом грунте представляют собой законы сохранения массы, энергии, дополненные обобщенным законом Дарси и уравнением состояния совершенного газа [2].

В области 0<х<% выполнены следующие уравнения относительно неизвестных функций — влаго-насыщенности, плотности воздуха в порах ра, температуры Т и давления Р.

д о д / х

пр —& + р — (V ) = 0 ;

dt

dx

n—р (1 -S )+—(р v ) = 0 ; dt Л w) dx

(1.1)

v j = —

-Px, P = раRT , j = a,w ; (1.2)

(PC) f + 1 IP ( + "" )

д d

+(р C v + C pv )—T = —

1 V' w w w 1 vi a a J Qx dx

л A t

Ш о

dx

(1.3)

где

А„ = п&А + п(1-) +(1- п);

(рС) = пБкрС + п(1-Бы)раСу +(1-п)рС . (1.4)

В мерзлом грунте (область у (0 < х<У) лед в порах и сам скелет неподвижны, следовательно, выполнены следующие уравнения:

пр1 д^ = 0 , пдра (1- & ) + ^(раиа ) = 0 ; (1.5)

dt

dt

dx

k d ^ =~^a Й )dxP' P = PaRT >

) f + di(PVa ) +

(1.6)

+C d T d +C p v —T = —

1 vra a r\ r\

dx dx

Л -d t

dx

где

Л = nSt\ + n(l — S Л +(1- n)

(1.7)

(1.8)

(рС) = п8грсС + п(1- & )раС +(1- п)рС. (1.9)

Нижние индексы а, м>, ¡, s относятся к характеристикам воздуха, воды, льда и скелета грунта соответственно, нижние индексы т, / указывают на усредненную (эффективную) характеристику соответственно талого и мерзлого грунта, которая здесь вычисляется как средневзвешенная.

На верхней границе задана температура осадков, насыщенность

Т (0,0 = Г,, 8(0,0 = ; (1.10)

и условие

kfw (S, x) dP

Pw

dn

= N (t),

(1.11)

x=0

где N (г) — скорость дождя (объем воды, выпавшей на единицу поверхности x = 0 за единицу времени).

Условия на поверхности фазового перехода, которая считается совпадающей с линией разрыва водона-сыщенности, выводятся стандартным образом из законов сохранения в предположении непрерывности температуры и давления. Таким образом, на границе фазового перехода x = ) выполнены условия непрерывности температуры и давления

[Т ] = 0, [р ] = 0, условие фазового равновесия

т •= тщ (Р),

(1.12)

(1.13)

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

5 -

Б Р,

Ф) = -—/„ 5 )д-Рт (Ф и); д-

(1.14)

( -Б,)Ф) = ~/а Р -

-/ (Б )—Р

^ а * п ' т

д-

д-

(1.15)

Кроме того, имеет место условие Стефана

д

д

пБррЛф)=т т,г) - \т—тт (Ф),г). (1.16)

д-

д-

Ш1- р> 1р» )> /- О)». /и-.

(1.19)

Для линейных зависимостей /п (Б) = аБ, /(Б) = = Ь(1 - Б), а,Ь > 0 условие (1.19) выполнено при

а к. < 1-_р_

" р

Ь и

Поясним наличие трех условий (1.14)—(1.16) на фазовой границе для нашего случая, когда для одной из искомых функций (влагонасыщенность) имеет место уравнение первого порядка с положительным наклоном характеристик.

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

т к № ,) < а,)=,

/п (5) - Р 1 Рп

что равносильно неравенству

<

1

)

J п к п '

/п ) - 5Рг 1 Рп

(1.20)

Здесь q — удельная скрытая теплота фазового перехода. Скорость фазовой границы в нашем одномерном случае равна ^(г).

На неподвижной границе x = Y заданы температура и давление

Т(У,г) = т, Р(У,г) = Р/. (1.17)

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

Т(-,0) = Т., Р (-,0) = Р , Б(-,0) = Б... (1.18)

4 ' ' гпг ' 4 ' ' гпг ' 4 ' ' гпг

Относительно функций /п ($) и /а (Б), фигурирующих в законе Дарси, предполагается следующее: обе эти функции неотрицательны, дифференцируемы на (0,1) и ограничены при 0 < Б < 1; /п (Б) обращается в ноль при нулевом значении аргумента и монотонно возрастает; /а Б обращается в ноль при значении аргумента, равном единице, и монотонно убывает. Кроме того, потребуем выполнения следующего неравенства:

при выполнении которого для определения водонасы-щенности в области х>£ (0) необходимо знать значение этой функции на свободной границе. Если считать, что поры в мерзлой зоне заняты льдом, то условие (1.20) приобретает вид:

<

1

/'(Б )

/п ) - Рг 1 Рп и выполняется как для линейной, так и для степенной (со степенью, меньшей единицы) функций насыщенности.

На неподвижной границе х = Y задана температура (1.17).

Задача замыкается заданием начальных распределений температуры, давления и насыщенности (1.18).

2. Преобразованная система. Предполагая в дальнейшем, что воздух в порах мерзлого грунта отсутствует, и исключая из систем (1.1-1.4), (1.5-1.9) с использованием закона Дарси и уравнения состояния скорости и плотность воздуха, получим следующую систему нелинейных уравнений относительно насыщенности, давления и температуры в областях газ — вода:

к д/п д_ Р дБ =_^/пР- ; (21)

Условию (1.19) очевидно удовлетворяют функции

/п (Б) = Ба, / (Б) = (1 - Б)в ,0 < а, в < 1.

дБ

дг п»п дБ д- дх п»п

(1-8)дР-(1-8)РдТ - р дБ - кк/^РР =

• ' дг • 'Т дг дг п»а --

= — (/Р + Р

а х 1

пРа

дБ

- /аТтх)Рх; (2.2)

пБрС + п(1-Б)—С +(1-п)р,С,

дт

-к

/ /

^ ы + ^ а

Р, Ра

РРхх -(КТ ) =

= (п(Аы - Аа 8 + к

/ / С Р

^р С

гы ы 1 У» Т"1

Р, Ра К Т

т +

Ры

дБ

РБ

/ ^ + /

^ а 1 •> п>

Ра

Р .

(2.3)

В области мерзлого грунта останется одно уравнение

дТ

(прС +(1 - п)рС))-(А тх) =0.

(2.4)

пРа

дх

8 - А

т=-—/„ (Бт )д-Рттл (2.6)

прш дх

(1.10), (1.12), (1.13), (1.17), считая в (1.13) равновесную температуру равной нулю градусов Цельсия:

Лт =-Т/4,Вт = Т, Л; =Т} /(у

в, =-&, /(У

Положение фронта ^0 находится из условия Стефана (1.16), в котором левая часть равна нулю. В результате получим

КТУ

$0 =

А Т, - А Т

т , / /

На границе фазового перехода выполнены условия (1.14)—(1.16), записанные для случая отсутствия в порах мерзлой зоны воздуха

-к дР

ф)(-Бы ) = —/„ (Бы(Ф)Л; (2.5)

Заметим, что коэффициенты Ат, А/ в этой формуле зависят от влагонасыщенности. Постоянная влаго-насыщенность в решении (Р0, 50, у0, Т/т(х)) может выбираться любой. Давление можно взять атмосферное.

Вводя отклонения Т/ т = / - Т/т (х), Р = Р - Р0, Б = Б - Бд, линеаризуем задачу (2.1)-(2.7).

В талой зоне 0 < х < ф)

Т =—/ы ; (3.1)

пРы

(1-)Т-(1-50ТТ -Р0 'Р

Т п

/ /

J ы + J а

Ры Ра

пР а

= — / Б)—Т°Р ;

•> а 4 0' т 0 х х '

(3.2)

пряФ) = А —Т (Ф)> *) - К дЪТт ). (2.7)

д_ дх

Заметим, что из условий (2.5), (2.6) следует необходимое для разрешимости задачи в случае отсутствия воздуха в мерзлой зоне условие на насыщенность на фазовой границе, которое имеет вид:

(1 - )Ра = (- рг / р„ )Ры (28)

/ (Б ) / (Б ) '

Переписав (2.8) в виде

(1 - Б )Р / (Б )/ Р / (Б ) = (Б - р. / р ),

нетрудно проверить, что в условиях (1.19), которым удовлетворяют функции (Б), /а (Б), уравнение (2.8) имеет решение Б ы такое, что р{ / р ш < Б ы < 1.

3. Построение стационарного решения и линеаризация задачи вблизи него. Задача (2.1)-(2.7) в случае отсутствия дождя имеет стационарное решение с постоянным давлением Р0 , постоянной водона-сыщенностью Б0, постоянным положением фронта и линейным распределением температуры в обеих фазах Т/ т (х), удовлетворяющих постоянным условиям на известных границах и условию (2.7) с нулевой левой частью. Найдем это решение. Обозначая Т/т(х) = Лг тх + Вг т, найдем константы из условий

пБр С + п(1-С +(1 -п)рС

0гы п> 1 V 0 / РТТ V ^ 5

Т -

-к

/ы Б) , /а Б)

РпР -(^А, + п(1-)\ +

Р, Ра

+(- п)А )Тхх - Т(А - п\ +(1- п)А )Т0х =

(пА - Ха)Б + к

/ы Б)

р С +

I XV Я' 1

+

/а Б) СР0

Т 0

(3.3)

Ра КТ0 ,

В мерзлой зоне возмущенное давление равно нулю:

прС +(1-п)рС^)) - (А +(1-п)\)Т- = 0. (3.4)

(3.5)

На фронте х = % выполнены условия

Т. =Т = 0;

/ т '

к дРТ ^(1 - Б0) =--/а (Б,) (ф), t); (3.6)

пРа 0 дх

Б - А

Рп

к

т = -—/п Б)—Рт(£(г),г); (3.7)

д

пр, я£(г) = (Х +(1-п)\)—Т} (£(г), г) -

-(( + п(1- Б Ж +(- п)\

)т а(г ),г)

' д-

-Т(п\ - пХ. )

9Т0 (-)

д-

(4.4)

-( + п(1-)Х. +(1-п)Х ))Тт т,г) -

-т), г)(пХп - п\)—г ю- (3.8) а-

Очевидно, в решении, на котором линеаризуется задача, мы должны подчинить теперь условию (2.8):

(1 - Б0)»а = (Б0 - Рг 1 Рп )»п

/а Б)

/п Б)

Возмущенная насыщенность на фронте находится из условия (3.8).

Отметим, что линеаризованная задача (3.1)-(3.8) остается нелинейной в силу того, что условия (3.6)-(3.8) задаются на заранее неизвестной границе, таким образом, мы имеем дело с задачей со свободной границей. Даже если считать уравнения выполненными в фиксированной области, и граничные условия заданы на известной границе, рассматриваемая задача относится к классу начально-краевых задач общего вида для системы, параболической по Петровскому. Для установления ее разрешимости необходимо проверять условие Лопатинского аналогично тому, как это делалось в [3] для несколько более простой системы.

4. Алгоритм численного решения линеаризованной задачи. Для решения линеаризованной задачи предлагается следующий алгоритм:

1. Задаем начальные условия для возмущений неизвестных функций

Т(-,0) = а(-), 0 < - < £0; (4.1)

Тт(-,0) = От(-), Р(-,0) = ш(-). (4.2)

2. Решаем задачу для возмущения полей температуры и давления (3.2), (3.3) при выполнении краевых условий.

На дневной поверхности -= 0,

Т(0,г) = Тт,, - г) = щ(г). (4.3)

»п д-

На свободной границе, совпадающей с границей фазового перехода, - = £(г),

[Т = 0, [Т] = 0, Т(- = £(г)) = 0;

к г (? )дРт(Ф),г) --/п (Л)---пРгЯ =

п»п д-

8, - А

На нижней границе - = Ь

Т(1, г) = 0.

3. Решаем (3.1) с найденным возмущением давления в правой части. Краевых условий не используем. В результате найдем везде Б(-,г), 0 < - < £(г).

4. Используем (3.8) для нового положения границы фазового перехода - = £(г + Аг).

5. Продляем линейно Т(-,г) на 0 < - < £(г + Аг).

6. Возвращаемся к выполнению шага 2.

Уравнения (3.1)-(3.8) переписываются в безразмерном виде с использованием обезразмеренных переменных:

' = ^, г' = г/, р' = р/ ,

/Ь /го /р0

Т'= ТТ , Рг = '' /Р , С г =

Р.

(РС) т =

(Рс)

пБ + п(1 - Б)

РСУ

КТРпСп

РпСп

- +(1 - п)Р\С\

(РС)

(РС)/ = /Р С = пР+(1 -п)Р'^С\'

пп

Х' = Хт/ Х' = Х//

Х т /\ - Х / /к

Для решения используется метод спрямления фронта [4]. Перейдем к новым координатам

' - - £(г)

£(г) - - £(г)

-; 0 < - < £(г); £(г )< - < 1.

1- £(г)

Запишем уравнения (3.1)-(3.8) в новых коорди-

натах:

дТ д8 д2 Б ^

^7 + а1ТТ = а!2Т-1 + -1 < 7 < 0; дг оу ду

Ъ (1 + У)

; а,2 =0.0;

(Х +(1-п)Х,

) (£(г),г) д-

= < ) 1 дТ

а,3 = (Б0) »2^2.

(4.5)

дТ дР д2Р

* + ар1 9^ = ар2 дУ^ + ар3 (г,£г -1 < У < 0

ар1 =-

^ (1 + у) Ъ/а (Б0) Р0Т0 -

ар2 =

£ (1 - Ж

Ро(7а + Ъ/п Б)) ; £2(1 - Бе)

арз =

Р ¡дТ £ (1 + у) дТ

Т0 [дг £ ду

дТ дТ д2Т

+ а±1 = а±2 — + а±3(г,£,£,Т); -1 < у < 0; дг ду ду

^ (1 + у)

£ (рс....

(у -1)

(1 - £) (рС

, -1 < у <0;

, 0 < у < 1;

а±2 =

а- = ,в\ , -1 < у < 0;

Т2 £2 (Рс)т < у <

вХ /

(1 - £)2 (рс)

-, 0 < у < 1;

(рС)

/п Б) + /а Б) ^

Р0 д2 Р

£2 ду2

+трС-иен - \)1 дБ+-1)Х + )+

+щ \/ {8о)+/{8о)А1 дц^Ш.

1 'п^п4 0у 1 •> а 4 0У ^ т / \ о о

»а Т0(у) £ ду \ д-

Условия сопряжения на свободной границе у = 0 примут вид:

дР

Б - Рг)

ду

Х/ дТ/

(4.6)

п1Л /п (Б0) ч

-МК - \)

____Т

£ ду 1-£ ду дТ0(-У

д-

(4.8)

£г = —

1 ¡х/ дТ х дТ

I / / \п г

т т

Щ \ £ ду 1-£ ду

-ёп(\, - Ха)

дТ0(-)

д-

(4.9)

В приведенных уравнениях использованы следующие обозначения для безразмерных параметров:

^1

кр0г0

п»„Ь

2 * ^а

кР0г0

Хг

-00 в_ п 0

Р С Ь

I н' И'

2

К = -

крХ

Х = -

а-3 = 0.0. (4.7)

и р С Т,Ь Я р С Т,

При аппроксимации уравнений полученной системы строится равномерная по пространству сетка. Уравнения для давления и температуры решаются неявно с использованием метода прогонки. Шаг по времени брался постоянным. Предположим, что значения возмущений давления, температуры, влагонасыщенности и положение свободной границы известны в момент времени гп = гп-1 +Аг. Для перехода на следующий слой времени сначала решаем уравнение для температуры (4.7) с граничными условиями на границах у = -1, у = 0, у = 1. Используем

Результаты расчетов процесса движения фазовой границы в зависимости от времени для двух значений интенсивности дождя (10 и 50 мм/ч)

найденные значения при решении уравнения (4.6) при -1 < у < 0. На свободной границе выполняется разностный аналог уравнения (4.8). Находим распределение влагонасыщенности, используя только значения влагонасыщенности на предыдущем шаге и найденные значения давления. Привлекая уравнение (4.9), вычисляем смещение границы фазового перехода. На этом заканчивается расчет на одном шаге по времени.

5. Результаты расчетов. Для расчета были выбраны значения исходных параметров: р{ = 910.0; ры = 1000.0; ра =1.29; = 2000.0 [кг/ м3] — плотности льда, воды, воздуха и скелета соответственно; С = 2000.0; С = 4190.0; С = 1.009х103; С = 1.009х103:

ы а V

С5 = 920.0 [Дж/ кг К] — удельные теплоемкости льда, воды, воздуха и скелета соответственно; Аг = 2.23; А = 0.58; А = 2.370х10-2; А = 2.0 [вт/(м К)] — ко-

эффициенты теплопроводности льда, воды, воздуха и скелета; ры = 1.002х10-3; Ра = 0.0182х10-3 [кг/м с] — вязкости воды и воздуха; R = 287.0 [Дж/ кг] — универсальная газовая постоянная; q = 3.34х105 — удельная скрытая теплота фазового перехода; п = 0.2 — пористость среды; k = 10-12 — проницаемость. Относительные фазовые проницаемости были представлены следующим образом: /ы ($) = Б05, /а (Б) = 1 - Б. Температура плавления льда составляет 273.2 К. На границах талой и мерзлой областей брались значения температуры 290 и 270 К соответственно. Температура дождя равна 300 К. Нижняя граница мерзлой зоны находилась на глубине L = 10 м. Приведенным значениям параметров соответствует стационарное решение с положением фронта $(0) = 0.8142 и постоянной влагонасыщенностью Б0 = 0.9249.

Библиографический список

1. Петрова А. Г., Железняк М. Н., Янцен В.В. Автомодельные режимы протаивания насыщенного мерзлого грунта при выпадении дождя // Известия Алтайского гос. унта. — 2014. — № 1 (81).

2. Васильев В. И., Максимов А.М., Петров Е.Е., Цып-кин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. — М., 1996.

3. Gianni R., Petrova A. G. One-dimensional problem for heat and mass transport in oil-wax solution // Rend. Mat. Acc. Linceu. — (2005). — V. 9.

4. Воеводин А. Ф., Гранкина Т. Б. Численное моделирование роста ледяного покрова в водоеме // Сиб. журн. ин-дустр. математики. — 2006. — Т. 9, № 1 (25).