CC BY
102
ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
УДК 519.63
П. Н. Вабищевич, С. П. Варламов, В. И. Васильев, М. В. Васильева, С. П. Степанов
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕПЛОВОГО РЕЖИМА ЖЕЛЕЗНОДОРОЖНОГО ПОЛОТНА В УСЛОВИЯХ КРИОЛИТОЗОНЫ
Рассматривается численное моделирование температурного режима земляного полотна железной дороги в условиях криолитозоны. Численная реализация построена на основе метода конечных элементов, позволяющего производить численное моделирование в областях со сложной геометрией с учетом слоистости грунта и наличия теплоизоляции. Проведено численное сравнение влияния сезонных колебаний температуры окружающей среды, снежного и напочвенного покровов на температурный режим грунтов основания железной дороги. Представлены результаты численного расчета для различных геометрических форм насыпей земляного полотна с учетом теплоизоляции пеноплексом.
Ключевые слова: железнодорожное полотно, задача Стефана, фазовый переход, метод конечных элементов, численное моделирование, многослойный грунт, снежный покров, теплоизоляция, криолитозона.
P N. Vabishchevich, S. P Varlamov, V I. Vasiliev, M. V Vasilieva, S. P Stepanov
Mathematical Modeling of the Thermal Regime of a Railway Line in Conditions of Cryolithozone
The numerical modeling of the thermal regime of a roadbed in conditions of cryolithozone is observed. The numerical implementation was made on the base of finite elements approach, with a help of which the numerical modeling in the field of complex geometry taking into account the ground banding and presence of thermal covering can be produced. The numerical comparison of the influence of seasonal fluctuations of temperature of the environment, snow and soil cover on the thermal regime of railway subsoil is held. The results of numerical calculation for various geometrical shapes of soil cover banket taking into account thermal covering with penoplex are represented.
Key words: roadway, Stefan problem, phase change, finite elements approach, numerical modeling, multilayer ground, snow banket, thermal covering, cryolithozone.
ВАБИЩЕВИЧ Петр Николаевич - д. ф.-м. н., профессор, зав. лабораторией Института проблем безопасного развития атомной энергетики РАН, ведущий ученый, научный руководитель ЦВТ ИМИ СВФУ им. М.К. Аммосова.
E-mail: vabishchevich@gmail.com
ВАРЛАМОВ Степан Прокопьевич - к. г. н., и. о. зав. лабораторией геотеплофизики и прогноза Института мерзлотоведения им. П. И. Мельникова Сибирского отделения Российской академии наук (ИМЗ СО РАН).
E-mail: vsp@mpi.ysn.ru
ВАСИЛЬЕВ Василий Иванович - д. ф.-м. н., профессор ИМИ, проректор по научной работе СВФУ им. М.К. Аммосова.
E-mail: vasvasil@mail.ru
ВАСИЛЬЕВА Мария Васильевна - к. ф.-м. н., начальник отдела параллельных технологий Центра вычислительных технологий Института математики и информатики СВФУ им. М.К. Аммосова.
E-mail: vasilyevadotmdotv@gmail.com
СТЕПАНОВ Сергей Павлович - аспирант Института математики и информатики СВФУ им. М.К. Аммосова.
E-mail: cepe2a@inbox.ru
Исследование изменений температурного режима грунтов под насыпью железнодорожного полотна является необходимым элементом инженерногеокриологического обоснования строительства железнодорожного пути в районах вечномерзлых грунтов. При сезонном оттаивании мерзлых грунтов изменяются его физико-механические свойства, что приводит к осадке оснований, нарушению устойчивости земляного полотна и увеличению затрат на ремонтные работы.
Для численного моделирования теплового режима железнодорожного полотна в условиях криолитозоны используется уравнение теплопроводности с учетом фазовых переходов поровой влаги вода - лед [1, 2]. Для моделирования в сложной геометрической области насыпей земляного полотна необходимо строить неструктурированные сетки, а для численного моделирования использовать конечно-элементные аппроксимации для возникающей начально-краевой задачи. При построении математической модели также необходимо учитывать сезонные колебания температуры воздуха, высоту снежного покрова и радиационный баланс подстилающей поверхности, влияющие на глубину протаивания грунта под насыпью железной дороги.
В работе рассмотрена математическая модель промерзания-протаивания грунта. Вычислительный алгоритм базируется на конечно-элементной аппроксимации температурного поля по пространству. Для аппроксимации по времени строится стандартная чисто неявная разностная схема. Возможности вычислительного алгоритма иллюстрщэуются расчетами модельной задачи для различных участков трассы железной дороги. Рассмотрено влияние теплоизоляции пеноплексом при различном его расположении.
Математическая модель. Рассмотрим математическую модель, описывающую распределение температуры при наличии фазовых переходов твердая фаза - жидкая фаза при некоторой заданной темпера-^ре фазового перехода Т* в области П = П~ U П+ (рис. 1) [1-6].
°д 1сь r + (t)- облаеть, занятая жидкой! фазой, где темперотура превышает темперааур= фазового перофода
/2+(t) = {х|х 6 /2,Н(хД) > Н*} и C_1t) е обббсть,занятн1 твердой фазой,
р_Пре ы {р|р £ рз, аПр, ре < а*}.
Фаооваш герехоа щюисходит на границе раздела
Г
п
г3
Рис. 1. Расчетная область
фаз S = S(tO-
Для моделирования процессов теплопереноса с фаоовыми перехоогми используется классическая моделt Стефана, описываюнкп тепловые процессы, сопровождающиеся фазовыми превращениями среды, поглощеси ем и выделением скр ытой теплоты
^ , зн
(а(ф) + р+сф к — - div(Я(фК grad НК = 0, (1)
где L - удельная теплота фазового перехода.
Для коэффициентов уравнения с использованием фун ющ и Хэвисайда ф ифеем следующие соотношения
афф) = р~с~ + ф(р+с+ — р~с~),
жф) = я + <-а+ — я)
, __ го, т < т*,
Ф ~ (1, т > т*,
где р+, с+ и р~, с~ -плотность и удельная теплоемкость талой и мерзлой зоны соответственно.
Поркоську рассматривается процесс распространения тепла в пористой среде, то для коэффициентов имеем:
рсс~ = К & — т)Я)С + тЯ+,
р+с+ = (1 — т)Я)С + тЯ,, где т - пористость. Индексы sc, w, i обозначают соответственно каркас портстой среды, воду и лед. Для коэффициентов теплопроводности в талой и мерзлой зоне имеем аналогичные соотношения
Я = (1 — т)Я)С + тЯ+,
&+ = (1 — т)Я)С + тЯ,.
На практике фазовые превращения не происходят мгновенно и могут' происходить I! малом интервале температуры [Т* — Д,Т* + Д]. В качестве функции ф можно взять ф2 [1, 6]
О, Т <Т* -А,
)Т — Т* + А Фь = {----—------, Т*-Д<7'<7'* + А,
> (а(фд) + РіЬф\) v dx + > (Л(фд) grad Г, grad v) dx Q(l- A) + I — а(Т-Гаіг)
I.
aRH+ l
- vds = 0, (5)
2A
1, T > T* + A,
О, Т <T* — А,
= '2А' Т* ~ А< Т <Т* + Л’
0, Т >Т* -1-А.
Тогда получим следующее уравнение для температуры 15 области 0::
■Те
Ма 6Я1(1).
Здесь Н1(П')- пространство Соболева, состоящее из функций а таких, что у2 и |Ка|2 имеем конечный интеграл.
Определим равномерную (для простоты) сетку по времени
"т = {и = П-Т,П = 0,1,..., Л тЛо = £тах} и проведем аппроксимацию по времени с использованием стандартной чисто неявной схемы [12]. Для линеаризации уравнения воспользуемся простейшей ланенризадией, ногда конффициенты зависят от значения функции с предыдущего временного слоя:
—И г
[аефа) + РіЬф а] — - diл,[Яeфa) grad И] н ф. (2- ф (а(фр + р1ф д)
+г!1 +г
дґ
Полученное уравнение (2) является стандартным нелинейным параболігаес ким уравнением.
Урнвнееие (Р) дополняется начальнымтслояилм
% dx ■
+ f (Д(фД) grad Tn+1, grad %) Cx
д 0 Q(1-A)+ I- a(dn+1 - dair) Jr aff + 1
(6)
lids = О,
T(x, вв = Г0, x Е Д
(3)
и граничными условиями
dT Q(1-a)+ I г a(TrTair)
-kd+=---------------атп-, иЄ Гі'
dT
-k — = С, и Є Гг, dn
(4)
Т = Т0, & е Г3,
где Р - суммарная коротковолновая радиация, А -альбедо, а - коэффициент конвективного теплообмена,
I - длинноволновое излучение, Т.г - еемпература наружного воздуха, Я - термическое сопротивление наземного покрова (зимой - снега) [ 7-9].
Конечно-элементная дискретизация. Для численного решения задачи проведем аппроксимацию уравнения (2) с учетом граничных условий с использованием метода конечных элементов [10, 11]. Умножим уравнение для температуры на функцию V и проинтегрируем с использованием формулы Грина
где Тп = Т(£п) 6 Я^О).
Таким образом, мы придем к следующей злассичес-кой вариационтой цистгшовке задачи: определить
Ы 6 Я"(СИ, (Г(х, 7) — Ыс(х, ЗИ)|г; еН"(СИ такую, что
- 9 (а(ФЫР + Рі^Ф-)^1 % сРх
т
Г х, Г аТп!"
+ $ 7’тг+1, ёг-а- г?) ф а3 + 1
= -[ (а(ф%) + р1Ьф'Зф)Тп % (іх
I
vds
(7)
+
Q(l-/)+ l + aTa
aR + 1
■vds,Vv й 71 (8).
Д/ численного решениямы должны перейти от непрерывной вариационной задачи (7) к дискретной [13]. Введем койечномейные пространства
и определим в них следующую задачу: определить Г? бК? , такую, что
7 [ >(0д) + Рг^0'д)а?г+1 %? т -'л
+ [ ()(Ф2) ёга-ТЦ+1,*га<1 %? )&" +п
# і
аТ?1+1
аЯ + 1
= -[ (а"ф1) + т ■)(
+ I
(8)
.(1-/)+ I + аТа
аЯ + 1
-%(І8, <&< Є V?.
Заметим, что выбор пространства @ непосредственно вытекает из типа применяемых конечных элементов.
Вычислительный эксперимент. Проведем численное моделирование температурного поля под железнодорожным полотном Амуро-Якутской магистрали «Беркакит - Томмот - Якутск». Расчетная область представлена из нескольких слоев грунта, состоящих из песка, суглинка и щебенистого грунта (рис. 2) [14, 15]. Грунт в слое годовых тепло-
оборотов имеет начальную температуру -3 оС. Температура на поверхности задавалась с учетом амплитуды колебания температуры воздуха, характерной для участка Томмот. Теплофизические параметры грунтов и тело насыпи приведены в таблицах 1 и 2. Будем предполагать, что фазовый переход при оттаивании-промерзании грунта происходит при
Рис. 2. Расчетная область и сетка
температуре 0 оС. Расчеты проводились на 3 года с шагом 1 день.
Результаты численных расчетов представлены на рис. 3. Иллюстрируется динамика температурного грунта под железнодорожным грунтом для различных сезонов (весна, осень) на третий год моделирования.
При моделировании учитывался не только конвективный теплообмен с окружающей средой, но также излучение, радиация, альбедо, высота и теплопроводность снежного покрова. Суммарная коротковолновая радиация Q, альбедо А, длинноволновое излучение I бралось в зависимости от сезонов (табл. 3). Высота снежного покрова изменялась по времени и достигала своего максимального значения в середине марта (0,35 метра).
Приведем результаты моделирования для различных условий на верхней границе (рис. 4). Сравнение проводилось при задании в качестве граничного условия только конвективного теплообмена с
Теплофизические параметры грунта основания насыпи
Таблица 1
Объемная теплоемкость [Дж/(м3 Град)*106] Теплопроводность [Вт/м Град] Уд.тепл. таяния [кДж/кг]
Грунт Талый Мерзлый Талый Мерзлый
Сугл. Заторф. 3,9 2,7 0,93 1,4 15000
Суглинок 3,15 2,35 1,51 1,7 71957
Супесь 1,64 1,47 0,64 0,7 50000
Песок 2,01 1,67 1,51 1,86 60437
Теплофизические параметры отсыпки насыпи и теплоизоляции
Таблица 2
Материал Теплопроводность [Вт/м Град] Объемная теплоемкость [Дж/(м3 Град)*106]
Тело насыпи Талый Мерзлый Талый Мерзлый
2,67 3,37 2,13 2,09
Пеноплекс 0,03 0,06
ВЕСНА ОСЕНЬ
Рис. 3. Распределение температуры по сезонам на третий год
Таблица 3
Параметры граничных условий
Периоды года Суммарная коротковолновая радиация ^) [Дж/Град м2] Длинноволновое излучение (I) [Дж/Град м2] Число Альбедо (А)
КО (Теплый период) (У-1Х месяцы) 198,93 63,44 0,3
КО (Холодный период) (Х-1У месяцы) 61,43 20,48 0,67
-4 -2 0 2 4
1
і
5 5
Рис. 4. Сравнение температуры осенью третьего года (1 - конвективный теплообмен, 2 - полный тепловой поток)
окружающей средой и для полного учета теплового потока для открытой местности.
Для увеличения охлаждения грунта на телах бермы и насыпи в практике используются теплоизоляция пеноплексом или сезонно-охлаждающие устройства (СОУ) (рис. 5). Результаты моделирования при различном расположении теплоизоляции и СОУ представлены на рис. 6. Охлаждение грунтов основания насыпи происходит при заложении пеноплекса под бермами. Заметим, что заложение пеноплекса выше основания оказывает более высокое охлаждающее
влияние на грунты основания. Из рис 6 d видно, что СОУ имеет такое же охлаждающее свойство, как пеноплекс, но установка СОУ имеет более высокие затраты.
Приведем также численное моделирование для различных геометрических форм полотна, представленных на рис. 7 и 8.
Представленные результаты численных расчетов иллюстрируют процесс сезонного оттаивания и промерзания грунтов для различных участков трассы железной дороги с использованием теплоизоляции
Рис. 5. Расчетные области при различном расположении пеноплекса Ь, c и СОУ)
Рис. 6. Распределение температуры осенью третьего года при различном расположении пеноплекса
Ь, c и СОУ)
Рис. 7. Распределение температуры по сезонам на третий год
пеноплексом и сезонно-охлаждающих устройств.
Работа выполнена при финансовой поддержке грантов РФФИ (проекты №12-01-98514 и №13-01-00719А).
Л и т е р а т у р а
1. Вабищевич П. Н., Самарский А. А. Вычислительная теплопередача. - М.: Едиториалл УРСС, 2003.
2. Васильев В. И., Максимов А. М., Петров Е. Е, Цыпкин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. - М.: Наука, 1996.
3. Andersland O. B., Ladanyi B. An introduction to frozen ground engineering. Chapman & Hall, 1994.
4. Harris J. S. Ground Freezing in Practice. Thomas Telford Limited, 1995.
5. Alexiades V., Solomon A. D. Mathematical modeling of melting and freezing processes. Hemisphere Publishing Corporation, 1993.
6. Samarskii A. A., Vabishchevich P. N., Iliev O. P., Churbanov A. G. Numerical simulation of convection/diffusion phase change problems - a review // Int. Journal of Heat and Mass Transfer. 1993. Vol. 36. No. 17. P. 4095-4106
7. Материалы Международной научно-практической конференции по инженерному мерзлотоведению, посвященной 20-летию создания ООО НПО «Фундаментстройаркос». - Тюмень: Сити-Пресс, 2011.
8. Павлов А. В., Перльштейн Г. З., Типенко Г. С. Актуальные аспекты моделирования и прогноза термического состояния криолитозоны в условиях меняющегося климата // Криосфера Земли. Академ. изд-во “Гео” - 2010. -Т. 14. - № 1. - C. 3-12.
9. Павлов А. В. Теплофизика ландшафтов. -
Новосибирск: Наука, 1979, 284 с.
10. Васильева М. В., Павлова Н. В. Конечно-элементная реализация задачи замораживания фильтрующих грунтов // Математические заметки ЯГУ. - 2013. - Т. 20. - № 1. -С. 202-212.
11. Зенкевич О., Морган К. Конечные элементы и
аппроксимация: Пер. с англ. - М.: Мир, 1986.
12. Самарский А. А. Теория разностных схем. - Наука. Глав. ред. физико-математической лит-ры, 1989.
13. Anders Logg, Kent-Andre Mardal, Garth N. Wells
Automated Solution of Differential Equations by the Finite
Element Method. The FEniCS Book. 2011.
14. Software package Gmsh. http://geuz.org/gmsh/
15. Software package Paraview. http://www.paraview.org/
ll
