CC BY
61
Вычислительные технологии
Том 19, № 4, 2014
Математическое моделирование
>к
искусственного замораживания грунтов*
П.Н. ВлвищЕвич1, М.В. Васильева2, В.Ф. Горнов3, Н.В. Павлова2
1 Институт проблем безопасного развития атомной энергетики РАН,
Москва, Россия
2 Северо-Восточный федеральный университет им. М.К. Аммосова, Якутск, Россия 3ОАО "ИНСОЛАР-ИНВЕСТ", Москва, Россия e-mail: vabishchevich@gmail.com, vasilyevadotmdotv@gmail.com, insolar-invest@yandex.ru, npav@rambler.ru
Вабищевич П.Н., Васильева М.В., Горнов В.Ф., Павлова Н.В. Математическое моделирование искусственного замораживания грунтов // Вычислительные технологии. 2014. Т. 19, № 4. С. 19-31.
Рассмотрен процесс термостабилизации грунтов вблизи фундаментов зданий и сооружений с использованием вертикальных замораживающих колонок. Математическая модель процесса описывается нестационарным уравнением теплопроводности с учётом фазовых переходов поровой влаги. Численная реализация базируется на методе конечных элементов с использованием библиотеки научных вычислений FEniCS. На модельной двумерной задаче приводится сравнение изотерм температуры и динамики погрешности на различных расчётных сетках в зависимости от выбора шага по времени и ширины интервала фазового перехода. Приводятся результаты численного моделирования трёхмерных температурных полей при наличии 39 замораживающих колонок с учётом температуры атмосферного воздуха на вычислительных системах кластерной архитектуры.
Ключевые слова: численное моделирование, искусственное замораживание, теп-лоперенос, фазовый переход, задача Стефана, метод конечных элементов, вычислительный кластер.
Vabishchevich P.N., Vasilyeva M.V., Gornov V.F., Pavlova N.V. Mathematical modeling of the artificial freezing of soils // Computational Technologies. 2014. Vol. 19, No. 4. P. 19-31.
Research of changes of soil temperature regime is necessary for the geotechnical studies of the construction objects in regions where permafrost soils are presented. During seasonal thawing of frozen soil, their physico-mechanical characteristics are changed, which leads to the violation of their bearing capacity. Artificial freezing of soils using cooling devices is applied near the piles for sustainability foundations of buildings and structures. It was done by creating a large mass or piles of frozen ground around boulders, which prevents the soil from thawing during the summer period.
For numerical simulation of the process, we design the mathematical model of heat and mass transfer that takes into account the freezing temperature of the liquid in the seasonal cooling devices and typical outside temperature for Yakutia (from -55 °C to +35 °C). Mathematical model of the process is described by nonstationary heat
* Работа выполнена при финансовой поддержке РФФИ (гранты № 12-01-98514 и 13-01-00719А) и Министерства науки и образования России (контракт RFMEFI57914X0026).
equation with phase transitions of pore water. The numerical implementation of the problem is based on the finite element method using FEniCS library of scientific computing. Comparison of the temperature isotherms and error dynamics at different computational grids is given for the two-dimensional problem depending on the time step and the phase transition width.
Results of numerical simulation of the three-dimensional temperature fields after five years (spring, summer, autumn, winter) are given in the presence of 39 freezing columns. The calculations were performed on a computational cluster of NEFU. The graphs illustrate the effectiveness of the use of freezing columns for sustainability the foundations of buildings and structures.
Keywords: numerical modeling, artificial freezing, heat transfer, phase transition, Stefan problem, finite element method, the computing cluster.
Введение
Исследование изменений температурного режима грунтов является необходимым элементом инженерно-геологического обоснования строительства объектов в районах залегания многолетнемёрзлых грунтов. При сезонном оттаивании изменяются физико-механические свойства мерзлых грунтов, что приводит к нарушению их несущей способности [1-4]. Для стабилизации или укрепления зданий и сооружений применяется метод замораживания грунтов. Искусственное замораживание грунтов с использованием охлаждающих устройств применяется вблизи свай для обеспечения устойчивости фундаментов зданий и сооружений за счёт создания вокруг сваи глыбы замерзшего грунта большой массы, которая предохраняет грунт от размораживания в течение летнего периода [5-8].
Замораживание грунтов может проводиться с помощью специальных холодильных установок или сезонных охлаждающих устройств, не требующих затрат электрической энергии в условиях резко континентального климата Крайнего Севера. Использование сезонных охлаждающих устройств позволяет также производить охлаждение грунтов в районах, где электричество не доступно.
При проектировании инженерных сооружений на вечномерзлых грунтах проводится научное обоснование технических решений. Для этого применяются методы математического моделирования технологических процессов взаимодействия оснований (фундаментов) зданий и сооружений с многолетнемерзлыми грунтами [9, 10]. Для численного моделирования рассматриваемого процесса строятся математические модели тепломас-сопереноса с учётом температуры замораживающей жидкости в сезонных охлаждающих устройствах и температуры атмосферного воздуха, задаваемого с учётом её типичного для Якутии сезонного колебания (от -55 до +35 °C) [11-13].
Подобные прикладные задачи связаны с необходимостью решения краевых задач в реальной геометрии. Расчётная область является трёхмерной и достаточно сложной. Для учёта геометрических факторов должны использоваться достаточно большие нерегулярные расчётные сетки, что учитывается, в частности, путём применения конечно-элементных аппроксимаций по пространству [14-16].
Особенностью моделирования термостабилизации фильтрующих грунтов является ярко выраженная геометрическая разномасштабность моделируемых объектов. При прикладном моделировании, даже в случае использования существенно неравномерных расчётных сеток, приходится ориентироваться на сетки достаточно больших размеров:
типичная из них содержит десятки миллионов ячеек. Численное решение таких задач в настоящее время невозможно без применения вычислительных систем параллельной архитектуры [17-19].
1. Математическая модель
Основная предпосылка при моделировании фазовых превращений твёрдая фаза — жидкая фаза состоит в том, что фазовый переход происходит при заданной постоянной температуре фазового перехода Т* [20 - 22]. Пусть фазовый переход происходит на границе раздела фаз S = S(£). Эта граница разделяет расчётную область П на две подобласти (рис. 1). Температура в подобласти П+(£) (талая зона) выше температуры фазового перехода Т*:
П+(г) = {х | х 6 П, Т(х,г) > Т*}, в подобласти П-(£) (мерзлая зона) — ниже температуры фазового перехода:
П-(£) = {х | х 6 П, Т(х, ¿) < Т*}.
Тепловое состояние грунта может быть охарактеризовано двумя выражающимися друг через друга функциями: энтальпией Н(х, ¿) и температурой Т(х, ¿). Энтальпия является монотонно возрастающей функцией температуры, имеющей скачок в точке замерзания грунтовой воды Т = Т*:
Н(Т) =
Н = р с Т,
Т < Т*,
Н + = р+Ь + р-с- Т* + р+с+(Т - Т*), Т > Т*,
где Ь — удельная теплота фазового перехода, р+, с+ и р -, с - — плотность и удельная теплоёмкость соответственно талой и мерзлой зоны. Поскольку рассматривается процесс распространения тепла в пористой среде, то для этих величин имеем следующие зависимости:
С-р- = (1 - m)Cscрsc + тС;р;, с+р+ = (1 - m)cscрsc + тстр-,
Рис. 1. Расчётная область
где т — пористость. Индексы зс, означают соответственно каркас пористой среды, воду и лед. Для коэффициентов теплопроводности в талой и мерзлой зонах используем
аналогичные выражения
Введём переменную ф как
и запишем энтальпию в виде
Л = (1 — т)Л5С + тЛ,, А+ = (1 — т)Л5С + тЛад.
ф
0, х е П-,
1, х е П+
н = н + р+ьф,
н = I р с Т, х е П ,
н1 = ^ р-с- Т* + р+с+(Т — Т*), х е П+.
Для Н1 верны следующие соотношения:
дН1
р с -
р+с+
дТ
дТ 3?
х е П -, х е П+,
дН1 дж,-
_ _дТ
р с -, г = 1, 2, 3, х е П
дж,
^дТ
р+с+ —, г = 1, 2, 3, х е П+ дж,-
Представим уравнение сохранения энергии в энтальпийной формулировке во всей области П = П - и П+:
дн
+ и grad Н — div(Л grad Т) = 0, (1)
где и — вектор скорости течения в пористой среде. Уравнение (1) можно записать через температуру в виде
а(Ф) (дТ + и grad Т^ + р+Ь ^дф + и grad ф^ — div(Л(ф) grad Т) = 0, (2)
где
а(ф) = р - с - + ф(р+с+ — р - с -), Л(ф) = Л- + ф(Л+ — Л-),
Г 0 Т < Т *
ф = ф(Т) = Т>Т*'
В вычислительной практике используется модель, основанная на том, что фазовый переход происходит в интервале температур [Т* — А, Т* + А]. В этом случае функция ф заменяется на ступенчатую функцию фд:
ф
д
0, Т < Т* — А,
Т — Т * + А
-, Т* — А <Т<Т* + А,
2А
1,
Т > Т* + А.
Далее имеем
где
дф _ дфдТ _ Ф дТ Ж _ дТЖ _ ФаЖ'
0, Т < Т* - А'
ФА _ { , Т* - А < Т < Т* + А, (4)
0, Т > Т* + А.
Подставляя выражения (3)-(4) в (2), получим следующее уравнение для температуры:
(«(Фа) + Р+^ФА) (5Т + и §Гаа т) - ^(Л(ФА) grad Т) _ 0, (5)
которое является стандартным нелинейным параболическим уравнением.
При отсутствии течения в пористой среде и _ 0 и уравнение теплопроводности с учётом фазовых переходов примет вид
5Т
(«(ФА) + Р+¿ФА) ^ - МЛ(ФА) gгadТ) _ 0. (6)
Уравнение (6) дополняется начальным условием
Т(х, 0) _ То, х € П,
и граничными условиями
Т _ Тс, х € Гд, дТ
-к— _ 0, х € Г/Гд. дп
Здесь Гд — место контакта с замораживающей скважиной.
2. Конечно-элементная аппроксимация
Для численного решения задачи проведём аппроксимацию уравнения (6) с использованием метода конечных элементов [15, 16]. Умножим уравнение для температуры на функцию V и проинтегрируем с использованием формулы Грина. Это дает
[ дТ
(«(ФА) + Рг^ФА) V
+ У (Л(фА)gгadТ,gгadV) ^г _ 0, Vv € Я,1(П). (7)
п
Здесь Я1 (П) — пространство Соболева, состоящее из функций V таких, что V2 и |Vv|2 имеют конечный интеграл в П и Я1 (П) _ {V € Я 1(П) : v|ГD _ 0}. Определим равномерную для простоты сетку по времени
Шт _ {Г _ П • Т, П _ 0, 1, ..., N0, Т^о _ ¿шах}
и проведем аппроксимацию по времени с использованием стандартной чисто неявной схемы. Для линеаризации уравнения воспользуемся простейшей линеаризацией, когда коэффициенты берутся с предыдущего временного слоя:
/Тп+1 _ Тп Г
(а(фД) + рг^фд)-ь^ж + (А(ф^гааТп+1, gradь) ¿ж = 0, (8)
п п
где Тп = Т(Г) е Н1 (П).
Таким образом, придем к следующей классической вариационной постановке задачи: определить функцию Т е Н 1(П), (Т(х, ¿) — Тс(х, ¿)) |Гд е Н^П) такую, что
1/ (а(Фд)+ Рг^фД^) Тп+1 ь^ж + J (А(фД)gгadТп+1, gгadь) ¿ж = пп
= ^ (а(фд) + Рг^Ф'Д) Тп ь^ж, Уь е Н^П). (9)
п
Для ее численного решения необходимо перейти от непрерывной вариационной задачи (9) к дискретной. Введем конечномерные пространства УН е Н 1(П), УН е Н^(П) и поставим в них следующую задачу: необходимо найти функцию ТН е УН такую, что
1/ (а(фд)+ РгЬфд) Тп+1 ьн ¿ж ^ (А(фД) gгad Т^1, gгad ьЛ) ¿ж = пп
= 1 / (а(фд) + Рг¿фд) ьн ¿ж, Уьн е Ун. (10)
п
Заметим, что выбор пространства УН непосредственно зависит от типа применяемых конечных элементов.
3. Численное исследование параметров задачи
При разработке вычислительного алгоритма основное внимание уделялось расчету температуры с учетом фазового перехода, включающему подбор ширины интервала фазового перехода А, шага по времени и определение их влияния на точность приближенного решения при использовании различных расчетных сеток. Такой методический прием для методов сквозного счета является стандартным (см., например, [21]).
Для построения двумерной геометрической модели и генерации неструктурированной сетки применяется программа ОМБИ [23]. Численная реализация производится с использованием свободно-распространяемого вычислительного пакета РЕшСБ, который реализует метод конечных элементов [16], являющийся одним из наиболее эффективных численных методов решения математических задач и позволяющий производить расчет в сложных геометрических областях. В качестве линейного решателя использован итерационный метод минимальных невязок (GMR.ES) с ЬИ-предобуславливателем вычислительного пакета ЕЕшСБ.
Рассмотрим численное решение задачи искусственного замораживания грунтов в двумерной расчётной области протяжённостью 2 м по горизонтальному и 2 м по вертикальному направлениям. В левом углу располагается замораживающая колонка радиусом 0.1 м, заглубленная в грунт на 1 м (см. рис. 1). В таблице представлены входные параметры задачи.
Приведём результаты численных расчётов в зависимости от вычислительных параметров задачи: величины шага по времени, расчётной сетки и параметра сглаживания А. Расчётные сетки представлены на рис. 2 и содержат 1014, 4144, 9206 треугольных ячеек. Расчёт проводился на 30 дней работы замораживающей колонки.
В качестве эталонного решения Те для вычисления погрешности было использовано решение на достаточно подробной сетке, содержащей около 30 тысяч треугольных ячеек при А _ 1, Т _ 0. 1 сут. Для сравнения ошибки вычислялась относительная погрешность решения
^ _ ^ТТ-р, ЦТе - Т||2 _/№ - Т)2^Х,
п
где Те — эталонное решение, Т — вычисляемая температура.
Параметры задачи
Обозначение Значение Единица Описание
измерения
Тс -30.0 град Температура замораживающей колонки
То 2.0 град Начальная температура
Т* 0.0 град Температура фазового перехода
ь 33.5 кДж/кг Удельная теплота фазового перехода
Свс 1.0 кДж/(кг • град) Удельная теплоёмкость грунта
С/ 4.2 кДж/(кг • град) Удельная теплоёмкость воды
Сг 2.1 кДж/(кг • град) Удельная теплоёмкость льда
Рзс 1500 кг/м3 Плотность грунта
Р/ 1000 кг/м3 Плотность воды
Рг 910 кг/м3 Плотность льда
Авс 1.1 Вт/(м •град) Коэффициент теплопроводности грунта
А/ 0.56 Вт/(м •град) Коэффициент теплопроводности воды
Аг 2.26 Вт/(м •град) Коэффициент теплопроводности льда
Рис. 2. Расчётные сетки. Слева направо: 1014, 4144, 9206 треугольных ячеек
Рис. 3. Динамика погрешности решения £ при т = 0.1, 0.2 и 0.5 сут
Рис. 4. Динамика погрешности решения £ на различных расчётных сетках: 1014, 4144 и 9206 ячеек
Рис. 5. Динамика погрешности решения £ при А = 0.5, 1 и 2
Ниже приводится сравнение динамики погрешности численного решения при разном т (рис. 3). Расчет проводился при А = 1 на расчетной сетке, содержащей 9206 треугольных ячеек. Видно, что шаг по времени достаточно сильно влияет на погрешность решения и расчет с минимальным шагом по времени т = 0.1 дает наиболее точное решение.
На рис. 4 представлено численное сравнение влияния размерности сетки на погрешность получаемого решения. Расчет проводился при временном шаге т = 0.1 сут, А = 1. Графики показывают сходимость решения при увеличении размера расчетной сетки.
Влияние ширины интервала фазового перехода на погрешность решения иллюстрирует рис. 5. Сравнение проводилось с решением задачи на сетке, содержащей 9206 ячеек при т = 0.1 сут. Из графиков видно, что большой интервал сглаживания, как и слишком малый, приводит к потере точности решения.
4. Моделирование термостабилизации оснований зданий и сооружений
Приведём результаты численного моделирования замораживания грунтов оснований зданий и сооружений. Расчётная область имеет протяженность 40 м по одной горизонтали, 50 м — по другой и 15 м — по вертикали (рис. 6). Вокруг здания расположено 39 замораживающих колонок радиусом 0.1 м, заглублённых в грунт на 8 м. Расчётные параметры задачи аналогичны представленным в таблице.
Особенностью моделируемой области является ярко выраженная разномасштабность моделируемых объектов: большая протяжённость расчётной области и малые размеры замораживающих колонок. При моделировании даже при использовании существенно неравномерных расчётных сеток приходится ориентироваться на сетки достаточно больших размеров. Для численного решения поставленной задачи была построена сетка, содержащая приблизительно 19 млн тетраэдральных ячеек [23-25].
Трёхмерная расчётная область строится с использованием графической программы NETGEN [24]. Процесс численного решения задачи представим в виде следующей схемы. До тех пор, пока t < tmax:
1) определяем текущее время t: t = t + т;
2) сохраняем значение с предыдущего временного слоя tprev = t;
3) вычисляем температуру наружного воздуха tair
tair = 41 sin((2n(t + 250))/365) - 10.2;
4) если температура в грунте меньше tair, то сезонно охлаждающие устройства перестают работать, иначе — температура в них будет равна tair;
5) вычисляем распределение температуры с использованием итерационного решателя;
6) сохраняем распределение температуры для текущего временного слоя в файл формата vtk.
Для визуализации полученных численных результатов использовалась программа PARAVIEW [25]. Вычислительная реализация основана на вычислительном пакете FEniCS, который использует технологию MPI для обмена данными при проведении параллельных вычислений, а параллельная реализация итерационных методов основана на библиотеке PETSc. Параллельная реализация программы запускалась на вычис-
Рис. 6. Расчётная область и сетка
лительном кластере Ариан Кузьмин Северо-Восточного федерального университета. Кластер состоит из 160 узлов (1920 ядер) и имеет пиковую производительность 23.5 ТПОРЯ. Это один из мощнейших гибридных суперкомпьютеров в РФ и самый мощный на Дальнем Востоке.
Рис. 7. Распределение температуры (0С) на пятый год при постоянной температуре на замораживающих колонках
Расчёты проводились на вычислительном кластере на пять лет работы замораживающих колонок с шагом по времени т _ 1 сут. Для численного решения задачи использовались 64 параллельных потока, время счёта составило около 25 мин.
Рис. 8. Распределение температуры (0С) на пятый год при наличии сезонных охлаждающих устройств
Численное моделирование проводилось с учётом сезонного колебания температуры воздуха на верхней границе расчётной области и работы замораживающих колонок, в которых с использованием холодильных установок поддерживается постоянная отрицательная температура (T = — ЗО 0C). На рис. Т показано распределение температуры в вертикальном срезе на пятый год для различных сезонов: весна, лето, осень и зима, на рис. 8 продемонстрировано моделирование температурного поля с использованием сезонных замораживающих устройств.
Представленные результаты численных расчётов иллюстрируют процесс искусственного замораживания грунтов с использованием замораживающих колонок для обеспечения устойчивости фундаментов зданий и сооружений за счёт создания вокруг оснований (фундамента) зданий ледопородной завесы, предохраняющей грунт от размораживания в течение летнего периода.
Список литературы
[1] Цытович Н.А., Сумгин М.И. Основания механики мерзлых грунтов. M.: Изд-во АН СССР, 1937.
[2] Докучаев В.В. Основания и фундаменты на вечномерзлых грунтах. M.: Гос. изд-во лит-ры по строительству, архитектуре и строительным материалам, 1963.
[3] Велли Ю.Я. Здания и сооружения на Крайнем Севере: Справ. пособие. M.: Гос. изд-во лит-ры по строительству, архитектуре и строительным материалам, 1963.
[4] Войтковский К.Ф. Фундаменты сооружений на мерзлых грунтах в Якутии. M.: Наука, 196S.
[5] Хакимов Х.Р. Вопросы теории и практики искусственного замораживания грунтов. M.: Изд-во АН СССР, 1957.
[6] Хрусталёв Л.Н., Никифоров В.В., Войтковский К.Ф. Стабилизация вечномерзлых грунтов в основании зданий. Новосибирск: Наука. Сиб. отд-ние, 1990.
[7] Andersland O.B., Ladanyi B. An Introduction to Frozen Ground Engineering. Chapman к Hall, 1994.
[S] Harris J.S. Ground Freezing in Practice. Thomas Telford Limited, 1995.
[9] Alexiades V., Solomon A.D. Mathematical Modeling of Melting and Freezing Processes. Hemisphere Publ. Corporation, 1993.
[10] Павлов А.В., Перльштейн Г.З., Типенко Г.С. Актуальные аспекты моделирования и прогноза термического состояния криолитозоны в условиях меняющегося климата // Криосфера Земли. 2010. Т. 14, № 1. C. 3-12.
[11] Бондарев Э.А., Васильев В.И. Искусственное замораживание фильтрующих грунтов // Численные методы решения задач фильтрации многофазной несжимаемой жидкости. Новосибирск: ИТПM СО АН СССР, 19S7. С. 3S-47.
[12] Васильев В.И., Максимов А.М., Петров Е.Е., Цыпкин Г.Г. Тепломассоперенос в промерзающих и протаивающих грунтах. M.: Наука, 1996.
[13] Васильева М.В., Павлова Н.В. Численное решение нестационарной задачи искусственного замораживания фильтрующих грунтов // Mатематические заметки ЯГУ. 2010. Т. 17, № 2. С. 14S-15S.
[14] Васильева М.В., Павлова Н.В. Конечно-элементная реализация задачи замораживания фильтрующих грунтов // Математические заметки ЯГУ. 2013. Т. 20, № 1. С. 195-205.
[15] Зенкевич О., Морган К. Конечные элементы и аппроксимация: Пер. с англ. М.: Мир, 1986.
[16] Anders Logg, Kent-Andre Mardal, Garth N. Wells. Automated Solution of Differential Equations by the Finite Element Method. The FEniCS Book, 2011.
[17] Vabishchevich P.N., Vasilyeva M.V., Pavlova N.V. Mathematical Modeling of Thermal Stabilization of Vertical Wells on High Performance Computing Systems. Arxiv Preprint, arXiv:1304.1625, 2013.
[18] Эндрюс Г.Р. Основы многопоточного, параллельного и распределенного программирования. М.: Вильямс, 2003.
[19] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб.: БХВ-Петербург, 2002.
[20] Samarskii A.A., Vabishchevich P.N., Iliev O.P., Churbanov A.G. Numerical simulation of convection/diffusion phase change problems — a review // Intern. J. of Heat and Mass Transfer. 1993. Vol. 36, No. 17. P. 4095-4106.
[21] Вабищевич П.Н., Самарский А.А. Вычислительная теплопередача. М.: Едиториал УРСС, 2003.
[22] Цыпкин Г.Г. Течения с фазовыми переходами в пористых средах. М.: Наука, 2009.
[23] Software package GMSH. (http://geuz.org/gmsh/)
[24] Software package NETGEN. (http://sourceforge.net/apps/mediawiki/netgen-mesher/)
[25] Software package PARAVIEW. (http://www.paraview.org/)
Поступила в 'редакцию 6 ноября 2013 г., с доработки — 4 апреля 2014 г.
