CC BY
31
УДК 517:551.342
Численное решение задачи протаивания мерзлого грунта с учетом инфильтрации осадков
А.С. Алейников, А.Г. Петрова
Алтайский государственный университет (Барнаул, Россия)
Numerical Analysis of Frozen Soil Thawing under the Influence of Rainwater Infiltration
A.S. Aleinikov, A.G. Petrova Altai State University (Barnaul, Russia)
Работа посвящена исследованию влияния интенсивности осадков на скорость протаивания вечной мерзлоты. Изучение процессов протаивания мерзлых грунтов имеет большое значение как для северных территорий, так и для сельскохозяйственных районов Западной Сибири. В общем случае такие задачи требуют учета многих факторов. В данной работе рассматривается одномерная задача без учета силы тяжести, испарения, наличия вегетативного слоя и излучения. Доказывается, что область между дневной поверхностью и фронтом протаивания занята грунтом, который рассматривается как пористая среда с неподвижным скелетом и порами, заполненными воздухом и водой. В мерзлом грунте лед в порах и сам скелет неподвижны, а все поры заняты льдом.
Математическое описание процессов тепломас-сопереноса делается на основе методов механики сплошных сред и включает уравнения, являющиеся следствиями законов сохранения массы, импульса, энергии и замыкающих уравнений состояния. Модель строится в следующих предположениях: вода и лед несжимаемые; воздух — вязкий совершенный газ; температура и давление — общие для скелета и пор; поверхность грунта подвержена воздействию дождя, выпадающего с определенной скоростью и температурой. Скорость осадков, представляющая собой условие на градиент давления, задается на поверхности грунта.
Задача решается в полной постановке, без линеаризации на стационарном решении. Задача со свободной границей преобразуется в задачу с неподвижной границей, строится алгоритм численного решения. Анализируются результаты численных расчетов.
Ключевые слова: фазовый переход, инфильтрация, задача Стефана.
БОТ 10.14258Лгуа8и(2017)4-12
The paper is devoted to the study of rainwater infiltration influence on temperature, water saturation and permafrost. The investigation of thawing processes оf frozen soils is of great importance both for northern regions with permafrost and for agricultural areas in the West Siberia.
In general case, the problem is very complicated due to many factors to be taken into account. We consider here a one-dimensional problem without including gravity, evaporation, vegetative layer, and radiation. We suppose that 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 constant temperature and pressure and pores containing only ice. 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 velocity of rain on the ground surface is assumed to be given, and that leads to assuming the condition for the pressure.
In our paper, we solve nonlinear system of equation in full statement without linearization on stationary solution. Free boundary problem is transformed into the problem in fixed domain. The numerical simulation algorithm is elaborated, the analysis of numerical results is provided.
Key words: phase change, rainwater infiltration, Stefan problem.
1. Постановка задачи. Различные модели задач тепломассопереноса в грунте исследуются в работах [1-3]. Основная задача построения модели и численного расчета — оценить зависимость скорости протаивания вечной мерзлоты от скорости выпадающих на поверхность осадков. Решаем одномерную однофазную задачу. Ось 0Х направлена вниз. Область 0 < х < ф) — грунт, который рассматривается как пористая среда с неподвижными скелетом и порами, заполненными водой и воздухом. Область ф) < х < У — мерзлый грунт, в котором лед в порах и сам скелет неподвижны, а насыщенность равна единице, т.е. все поры заняты льдом.
Модель строится в следующих предположениях: вода и лед несжимаемые; воздух — вязкий совершенный газ; температура и давление общие для скелета и пор; поверхность грунта подвержена воздействию выпадающего с определенной скоростью и температурой дождя. Силой тяжести пренебрегаем. Определяющие уравнения в талом грунте представляют собой законы сохранения массы, энергии, дополненные обобщенным законом Дарси и уравнением состояния совершенного газа [4, 5].
В области инфильтрации осадков выполнены следующие уравнения относительно неизвестных функций —влагонасыщенности, плотности воздуха в порах ра, температуры Т и давления Р:
nPw +Pw = 0,n dt Pa (1 - Sw) + dX (pava) = 0,
L/i L/Л L/t С/Л
v = - ША (p -p g), P = pRT, j = a, w, U
(pC) — + — [P (v + v )] + (p C v + p C v )—T = —
yr 'm L 4 a 1 w'* yr w w w 1 ' a a a' о о
dt dt dx dx
л A T
m
OX
(1)
где
f (S ) = S , f (S ) = (1 - S ),
-'w4w/ w ^ J a x w ' 4 w' ^
лт = nSл + n(1 - )\a + (1 - n)Xs,
(pN)m = nStptCt + n(1 - St )paCa + (1 - n) psCs.
Здесь p, pw, pa — плотности льда, воды и воздуха; jUw, Ua — коэффициенты вязкости; Cw, Ca, Cs — удельные теплоемкости; k, q — проницаемость и удельная скрытая теплота фазового перехода; R—уни-версальная газовая постоянная; n — пористость.
Нижние индексы a, w, i, s относятся к характеристикам воздуха, воды, льда и скелета грунта соответственно, нижние индексы m,f указывают на усредненную (эффективную) характеристику соответственно талого и мерзлого грунта, которая здесь вычисляется как средневзвешенная [5].
На дневной поверхности заданы температура осадков, насыщенность
T(0, t) = T0, S(0, t) = S0;
и условие
kSwP0 dP
= N (t):
Л ^ Ф) дп
где — скорость дождя (объем воды, выпавший на единицу поверхности за единицу времени).
На границе фазового перехода х = ф) выполнены условия для температуры и давления
Т(фи) = Тр , Р(фи) = Рт .
Кроме того, выполнены следующие условия, являющиеся следствиями сохранения:
s -а
w
Pw
<(t) = - — fw (Sw )-дPm (ф), t) ;
npw dx
-fa (Sw ) d-Pm № ),t) dx
S -1)- i(t) = -—
д
nqpi(t) = -Лт — Tm (Ф), t) . dx
Последнее условие — это условие Стефана [6], в котором q—удельная скрытая теплота фазового перехода. Скорость фазовой границы в нашем случае равна ф).
Для того, чтобы замкнуть задачу, нам необходимо задать начальные распределения температуры, насыщенности и давления:
Т(х,0) = Т(х);
Б(х,0) = Б(х);
Р (х,0) = Р (х).
2. Преобразованная система. Будем искать реше- с криволинейной границей в другую задачу, опреде-ние, для которого все поры в мерзлом грунте заняты ленную в прямоугольной области путем замены пе-
льдом, а температура тождественно равна температуре фазового перехода (однофазная задача).
В целях численного исследования исходная система (1) будет приведена к безразмерному виду и рассмотрена в постановке без учета силы тяжести.
Применим метод выпрямления фронта [7]. Идея метода состоит в преобразовании задачи в области
ременной. Введем новую переменную
П =
Ф)
Введем безразмерные величины согласно [4]:
Т
С
Р =— ,Т = —, а: = ,С* = ,í =-,р* ,р =^ ; Р Т А ' С Г
аШ ±аЬ$ ь0
А р*
/ х* (рС) Р ■ С * *
(рС) = = пБ + п(1 -5 )-+ (1 -п)р С .
Ш „ п * * ЮТ „ п $
р*С* К1 ■ р*С*
Исходная система (1) примет следующий вид:
* * Р 1+ П -Т
Б' - а2ББ ■ Р"
I 2 п пп
а
Р ' - Рт' = а' РПП+ в'
1 >+П"Р, '= 71*Г+72Тп' + 7; + П**в *
(2)
где
1 а о (
щ =. АША* 10ТТ0
Qi(о2 ■ (ТТо - 7РРо)
Ж =-
X■ К ■Т
7а ¡оТР0
-Р'
в7 КТ
(ТТ0 - 7РР0) ^ О2 ■ (ТТ0 -7РР0) п (ТТ0 -7РР0)
Приведем систему (2) к виду
Б' - а,Б '= Б ■ Р"
2
ра
Р' = ер +ЩРР+е,
Т' = ж2т" + + щ
(3)
где
е,=-
а 7„
е =
£( О2 ■ (ТТ0 - 7РР0)
А А.гР
х ■ { р+ЁА - ЩР т"-М ■ ф)т>
0 р0 аа о щ а о п
Т1
(ТТ0 -7РР0)
3. Алгоритм численного решения. Для дальнейшего решения предлагается следующий алгоритм:
¡.Задаем начальные условия для неизвестных функций
Т(п>0) = VТ + (1 -п)'Т ;
На свободной границе:
т(1, г) = Тг, Р(1,г) = Рт.
3. Решаем задачу для насыщенности с условиями: На границах:
^ п = п' + (1 - п)' $0; Р(п,0) = п2' Р + п Рг + Р .
$(0, г) = $0; $(1, г) = .
2. Решаем задачу для температуры Т и давления Р с известной насыщенностью 5 и положением фронта 4 (1) при выполнении краевых условий: На дневной поверхности:
Т (0, г) = Т;
кБЕ дР
V.а 1К№) дп
= N (г)
Для существования решения необходимо взять начальную насыщенность в интервале (0.92..1)[8].
4. Находим новое положение свободной границы:
&) = дТ 1 г).
пчр1 дп
5. Возвращаемся к выполнению шага 2.
Для решения системы (3) воспользуемся неявной разностной схемой. Разностные методы подробно рассмотрены в [9,10].
Тп+1_тп тп+1_тп+1 тп+1_2ТП+1 + тп+1
т к к
РП+1 _РП РП+1 _РП+1 Пп+1 О Пп+! I пп+1
= WЛ
Рп+1 — о Рп+1 Рп+
^+е1 '+1 ' 2 + Р_
+е.
х=0
т
Будем искать решение в виде
Тп+ = а,..Тп++ + ви "21^1+1 + в21
111±1+1
Рп+Х = а. РР
Вычисляем Тп+1, Рп+1, используя значения 5п,
пп шп
ь 1,2,3.
В уравнении для 5 после перехода в фиксированную область мы должны проследить поведение харак-
с1г]
теристик • Другими словами, из уравнения перво-
го порядка не всегда можно найти значение насыщенности на фазовой границе. Дело в том, что даже при положительном наклоне, если наклон характеристик на линии раздела фаз меньше, чем наклон самой фазовой границы, то мы не сможем найти значение насыщенности, так как начальные условия задаются только на интервале (0, % (0)).
Рис. 1. Поведение характеристик (вертикальная — ось горизонтальная — ()
В зависимости от знака характеристик будем решать уравнение схемами «по» и «против» потока:
°г > 0
_n + 1 _n _n+1 _n+i ^
-- а
i. = s. P" ;
а2 <
- = oLs. p„ .
2 h n
Используя Tn+1, Pn+1, полученные выше, найдем значения Sn+1 и новое положение фазовой границы <(t).
0,855 0,85 0,845 0,84 0,835 0,83 0,825 0,82 0,815 0,81
• ^^
|06tia
54 мм/ч
12 14 £ (в часах) 11 мм/ч
Рис. 2. Результаты расчетов процесса движения фазовой границы в зависимости от времени для трех значений интенсивности дождя (0, 11 и 54 мм/ч)
4. Результаты расчетов. Результаты численных расчетов показали, что скорость выпадения осадков на поверхность оказывает значительное влияние на скорость протаивания грунта. Так, вечная мерзлота, расположенная на глубине 0,8142 м, в случае отсутствия осадков протаивает за двенадцать часов
на 2.244 см; небольшого дождя (11 мм/ч) — 2.517 см; сильного ливня (54мм/ч) — 3.562 см. Осадки, выпадающие на поверхность, имеют температуру 20°. Получено качественное совпадение движения фазовой границы с результатами [4].
Рис. 3. Результаты расчетов процесса движения фазовой границы в зависимости от времени для трех значений интенсивности дождя (0, 11 и 54 мм/ч)
При увеличении начальной глубины залегания до 6 м получены следующие результаты: в случае отсутствия осадков вечная мерзлота протаивает за двенадцать часов на 0.313 см; небольшого дождя (11 мм/ч) — 0.344 см; сильного ливня (54мм/ч) — 0.492 см.
Из полученных результатов видно, что большое влияние оказывают не только интенсивность осад-
ков, но и глубина залегания. В зависимости от глубины залегания скорость протаивания может изменяться в несколько раз.
Для расчетов выбраны значения исходных параметров, использованные в [4].
Библиографический список
1. Папин А.А. Краевые задачи двухфазной фильтрации : монография. — Барнаул, 2009.
2. Коробкин А.А., Папин А.А., Хабахпашева Т.И. Математические модели снежно-ледового покрова : монография. — Барнаул, 2013.
3. Ахмерова И.Г. Модельная задача фильтрации воды и воздуха в деформированном грунте // МАК-2015 : материалы Всероссийской конференции по математике. — Барнаул, 2015.
4. Петрова А.Г., Мошкин Н.П., Жирков А.Ф. Задача о возмущениях фазового фронта в ненасыщенном грунте под действием инфильтрации осадков // Известия Алтайского гос. ун-та. — 2015. — № 1/1 (85). DOI: 10.14258/ izvasu (2015)1.1-18.
5. Васильев В.И., Максимов А.М., Петров Е.Е., Цып-кин Г.Г. Тепломассоперенос в промерзающих и протаивающих грунтах. — М., 1996.
6. Мейрманов А.М. Задача Стефана. — Новосибирск, 1986.
7. Воеводин А.Ф., Гранкина Т.Б. Численное моделирование роста ледяного покрова в водоеме // Сиб. журн. ин-дустр. математики. — 2006. — Т.9, №1 (25).
8. Петрова А.Г., Алейников А.С., Бочкарева Ю.А., Ми-хина Д.Л. О точных решениях задачи протаивания грунта под действием инфильтрации осадков // МАК-2015 : материалы Всероссийской конференции по математике. — Барнаул, 2015.
9. Кузнецов Г.В., Шеремет М.А. Разностные методы решения задач теплопроводности. — Томск, 2007.
10. Крайнов А.Ю., Миньков Л.Л. Численные методы решения задач тепло- и массопереноса. — Томск, 2016.
