CC BY
632
НАУКА и ОБРАЗОВАНИЕ
Эя №ФС 77 - 3056-9. Государствен над регистрация №0421100025.155Н 1994-0405_
Математическое моделирование температурных полей с учетом фазовых переходов в криолитозоне 77-30569/354740
# 04, апрель 2012 Крылов Д. А.
УДК. 536.42
МГТУ им. Н.Э. Баумана [email protected]
Будущее топливно-энергетического комплекса России связано с освоением нефтегазовых ресурсов арктических шельфов. Задачи, связанные с изучением, прогнозом и управлением температурными полями в грунтовых средах, в настоящее время являются одними из наиболее важных и востребованных в сфере инженерной геологии и строительства. В современных условиях, когда на территории Российской Федерации происходит разработка новых месторождений нефти и газа, все чаще поднимается проблема строительства в районах вечной мерзлоты, ведь они занимают 64 % от площади РФ и 25 % площади всего мира. В их строении выделяются ярусы охлажденных ниже 0°С и многолетнемерзлых пород. Составление прогноза изменений температурного режима грунтов является необходимым элементом инженерно-геологического обоснования строительства (реконструкции, расширения) объектов нефтегазовой промышленности в районах распространения вечномерзлых грунтов. Целью данной работы является построение математической модели, позволяющей составлять подобные температурные прогнозы. Эта модель должна учитывать фазовые переходы в трехмерной области.
Прогноз изменений температурного режима грунтов может выполняться приближенными аналитическими методами, с помощью моделирования, а также численными методами. Задача определения температурного поля в процессе термостатирования многофазных веществ является нелинейной и относится к классу задач со свободной границей, в которых часть границы не задана и должна быть определена вместе с решением системы дифференциальных уравнений теплопроводности.
Научная новизна работы заключается в следующем. С помощью интегро-интерполяционного метода (метода контрольных объемов) была получена дискретная
постановка задачи Стефана, сводящаяся к задаче нахождения решения системы линейных алгебраических уравнений с трехдиагональной матрицей. Решение такой системы было получено с помощью метода прогонки (для одномерного случая) и метода переменных направлений (для двух- и трехмерного случаев). Задача была приведена к безразмерному виду.
Температурный режим (совокупность последовательных температурных полей в грунтовом массиве, соответствующих любым заданным моментам времени от начала расчета) рассчитывается как результат задаваемых на весь период расчета прогноза тепловых воздействий на верхней, боковых и нижней границах грунтового массива.
Расчетный алгоритм составлен и может быть применен для расчетных областей, ограниченных прямоугольной верхней границей, в контурах которой фиксированы тепловыделяющие объекты. В трехмерном случае область исследования параллелепипед, в двухмерном — прямоугольник и в одномерном — прямая.
В разработанном алгоритме учитываются условия на границе слоев, условия на границе фазового перехода воды в лед, физические свойства грунтов: суглинок бурокоричневый, супесь буровато-коричневая с прослоями суглинка и песка, песок пылеватый. Теплофизические характеристики корректируются по ходу вычислений при изменении влажности и температуры. Величина солнечной радиации определяется в зависимости от месяца и географической широты.
Численная реализация выполнена на языке Fortran. Программа способна рассчитать температурное поле в заданной области исследования при наличии определенного набора входных данных: физические характеристики области
исследования, геометрическое разбиение области исследования, начальное распределение температуры, граничные условия.
Следует отметить, что большая часть работ, связанных с температурными полями в таких средах носит экспериментальный характер, поэтому численное моделирование прогноза физических характеристик является особенно актуальным. Предшествующие работы посвящены решению подобных задач в двумерных и одномерных областях, либо построению моделей трехмерных моделей без учета фазовых переходов. Автором была рассмотрена подобная задача в энтальпийной постановке, с которой можно ознакомиться в [1]. Разработанная модель может найти применение в таких отраслях как строительство в зоне вечной мерзлоты, геология, космонавтика и др.
1 Теплопроводность твердых тел
Процесс теплообмена представляет собой перенос энергии, происходящий между телами (средами), имеющими различную температуру. Тепловое поле на данный момент времени t определяется распределением температуры по телу, т.е. функцией и(x,y, z,t), где x, y, z — декартовы координаты. В простейшем случае тепловой поток направлен из области с более высокой температурой в область с более низкой температурой.
Основное положение теории теплопроводности, известное как закон Фурье, состоит в предположении пропорциональности теплового потока градиенту температуры в однородной среде: q = -Л grad и, где Л - коэффициент теплопроводности. Для получения уравнения переноса тепла запишем закон сохранения энергии в виде
сРд~ = -divq + s, (11)
где s определяет мощность внутренних источников теплоты, с - удельная теплоемкость, р - плотность среды.
Подстановка выражения для потока в (1.1) позволяет записать основное дифференциальное уравнение теплопроводности
ди
ср— = div (Л grad и) + s . (1.2)
dt
Коэффициенты и правая часть уравнения могут зависеть как от точки
пространства (при неоднородной среде), так и от температуры.
Если теплофизические свойства среды постоянные (однородная среда), то
уравнение теплопроводности (1.2) упрощается и принимает вид
ди . s
— = а Ди +----,
dt ср
Л
где а =-----коэффициент температуропроводности, Д = div grad - оператор Лапласа.
сР
Другой частный случай уравнения теплопроводности (1.2) соответствует описанию установившихся тепловых полей. Стационарное уравнение теплопроводности имеет вид div (Л grad и) = -s и является эллиптическим уравнением второго порядка.
При моделировании теплопереноса в движущейся среде (движущейся системе координат) уравнение теплопроводности (1.2) соответствующим образом трансформируется. Переход основан на замене частной производной по времени d/dt на
полную производную по времени (субстанциональную производную) d = — + v grad,
dt dt
где v - локальная скорость среды. Поэтому уравнение теплопроводности в движущейся среде записывается в виде
(ди Л
ср\-----+ v ■ grad и I = div(A, grad и) + s. (13)
J
Здесь член с v ■ grad и определяет изменение температуры за счет конвективного переноса.
Распределение температуры в точках среды в различные моменты времени определяется из уравнения теплопроводности. Для однозначного определения температуры поля помимо уравнения необходимо сформулировать дополнительные (замыкающие) соотношения, так как решения уравнений в частных производных определяются с точностью до некоторых произвольных функций.
Для этого дифференциальное уравнение в общем виде дополняется: уравнениями состояния, уравнениями неразрывности, условиями в начальный момент времени, условиями на границах тела, данными о геометрии (в частности, условиями симметрии), о теплофизических свойствах материала, а иногда и другими замыкающими задачу сведениями.
Отдельно рассмотрим условия на границе контакта двух сред с различными теплофизическими характеристиками. Прежде всего, остановимся на вопросе, какие условия сопряжения естественны для уравнения теплопроводности, и, следовательно, их можно явно не формулировать. В этом случае условия контакта обуславливаются только разрывами теплофизических характеристик при переходе границы сред.
Будем считать, что контакт двух однородных сред проходит по плоскости x = 0 . Теплофизические характеристики, отнесенные к среде, занимающей полупространство x > 0 , будем помечать плюсом, а для второй среды (когда x < 0) будем использовать обозначения с минусом. Рассмотрим уравнение теплопроводности в виде
ди д ( ди Л д ( ди Л д ( ди Л ( ,
ср— = — \Л—I + — Л— + — \Л—I + s, (1.7)
dt dx ^ dx J dy ^ dy J dz ^ dz J
причем коэффициенты разрывны при x = 0 :
|с+,р+ ,Л+,s +, x > 0, с, р,Л, s = \ (1.8) [с , р , Л , s , x < 0.
С учетом такого разрыва коэффициентов уравнения (1.7) естественно считать, что само решение уравнения (температура и) непрерывно, а вот его первые производные уже
разрывны. Поэтому напишем условие непрерывности температуры при переходе из одной среды в другую в форме
[и ] = 0, X = 0, (1.9)
где через [ ] обозначен скачок при переходе границы контакта. В рассматриваемом случае [и] = и(х + 0, у, г) - и(х - 0, у, г). Осталось сформулировать условия сопряжения для первых производных температуры.
Для того, чтобы получить естественные условия сопряжения для уравнения теплопроводности (1.7) с разрывными коэффициентами (1.8), можно поступить следующим образом. Выделим на границе контакта некоторую ограниченную область 38 и проинтегрируем исходное уравнение теплопроводности (1.7) по области шириной 2е
={(х,у,г)\ -е < х <е,(у,г) &38}.
С учетом непрерывности температуры (условие (1.9)), конечности разрывов при е ^ 0 получим
ди ди
[| a+ — -A~—\ dydz = 0. Dx Dx >
SS4
В силу произвольности выбора элемента 38 получим условие сопряжения в виде
= 0, х = 0 . (1.10)
A^U.
6x
Условие сопряжения (1.10) отражает непрерывность теплового потока.
Естественные условия сопряжения записываются аналогично (1.9), (1.10) на произвольной криволинейной границе сопряжения двух сред £ внутренней границе области Q. Они отражают непрерывность температуры и теплового потока и записываются в виде
[и] = 0, (х,у,z) є £, (111)
A^U
6x
= 0, (x, y, z) є S. (1.12)
Уравнение (1.11) характеризует условие непрерывности температурного поля, а уравнение (1.12) - закон сохранения энергии на поверхности соприкасающихся тел при условии идеального теплового контакта. Подчеркнем еще раз, что условия (1.11) и (1.12) являются естественными для уравнения (1.7) с разрывными коэффициентами (1.8), они включены в само уравнение и нет необходимости каждый раз формулировать их отдельно.
Если на границе двух сред не выполнены условия идеального контакта (непрерывности температуры и теплового потока), то необходимо специально формулировать условия сопряжения, которые дополняют уравнение теплопроводности, записанное в каждой отдельной среде. Эти условия сопряжения отражают особенности тепловых процессов на границе контакта, особенности поведения решения задачи при переходе границы контакта и не всегда могут включаться в само уравнение теплопроводности.
Пусть на границе контакта имеется поверхностный тепловой источник мощности д5. Тогда справедливо предположение о непрерывности температуры, т. е. выполнено
условие (1.11), а вот тепловой поток претерпевает разрыв. Вместо однородного условия сопряжения (1.12) необходимо записать неоднородное условие сопряжения вида
хди_
дх
(1.13)
Условия сопряжения могут включаться в само уравнение теплопроводности, которое записывается для всей расчетной области без выделения границы контакта 8. Поверхностный источник тепла учитывается дополнительным слагаемым в уравнении (1.7)
ди д (. ди Л д (. ди Л д (. ди Л
ср— = — \Л—I +---Л —
дї дх ^ дх) ду ^ ду )
(1.14)
где 58 есть поверхностная 3 -функция, определяемая так, что для любой функции Г(х, у, г) справедливо соотношение
|35 Г(х,у, г) ёхёу ёг = |Г(х,у, г) ds .
В вычислительной теплофизике выделяют и другие типы условий сопряжения. В качестве примера отметим условия типа сосредоточенной теплоемкости, когда также имеет место непрерывность температуры, а разрыв теплового потока определяется соотношением
ди
Лд~и
дх
= с8р^т. (хy,z) є £, дї
(1.15)
где с8 — сосредоточенная теплоемкость контакта. Условия сопряжения (1.11), (1.15) с
использованием поверхностной 3 -функции включаются в уравнение теплопроводности подобно уравнению (1.14).
В прикладных исследованиях большого внимания заслуживают условия неидеального контакта. Такой случай реализуется, например, при недостаточно плотном соприкосновении шероховатых твердых тел. При неидеальном контакте тепловой поток непрерывен, что является отражением закона сохранения энергии. Таким образом, одно условие сопряжения мы имеем — условие (1.12). Температура при переходе границы неидеального контакта терпит разрыв, пропорциональный тепловому потоку
[м] = — А—, (х,у, z) е £. а дп
Здесь коэффициент контактного теплообмена а связан с условиями контакта.
2 Математическая постановка задачи
В классической задаче Стефана рассматриваются фазовые превращения с участием твердой фазы — твердое тело—жидкость. Примером таких превращений являются процессы затвердевания и плавления в металлургии. Соответствующие математические модели характеризуются наличием подвижных заранее неизвестных границ фазового перехода — задач со свободными (неизвестными) границами.
Основная предпосылка при моделировании фазовых превращений твердое тело— жидкость состоит в том, что фазовый переход происходит при заданной постоянной температуре фазового перехода и*. Пусть фазовый переход происходит на границе раздела фаз, которую мы обозначим S, причем S=S(t). Эта граница разделяет расчетную область Q на две подобласти. Область Q+ (t), занятая жидкой фазой, где температура превышает температуру фазового перехода, есть Q+ (t) = {(x, y, z) e Q, u(x, y, z, t) > и *}. Соответственно, область Q- (t), занятая твердой фазой,
Q- (t) = {( x , y, z) e Q, u(x, y, z, t) < и *}. Аналогичные обозначения используем и для теплофизических величин в каждой отдельной фазе.
Выпишем соответствующие уравнения теплопроводности. В твердой фазе имеем
с Р ~~ = div(z~grad и”)+ s_, (x, y, z, t) e 0-, (2.1)
где 0-= {(x, y, z, t)|(x, y, z) e Q- ,0 < t < tmax }. Учитывая конвективный перенос в жидкой фазе, получим
ди+ ^ \
+ ~+--------+ v ■ grad и+ = div(X~grad и+ )+ s+, (x, y, z, t) e 0+, (2.2)
V dt J
c p
где 0+ ={(^ У, ^ t)|(^ ^ z) eQ+ ,0 < t < tmax }.
Нас интересуют условия на границе фазового перехода 8. Прежде всего на этой границе контакта двух сред справедливы предположения о непрерывности температуры (1.11), т. е.
[и] = 0, (х,у,і) є 8. (2.3)
Фазовый переход сопровождается выделением/поглощением определенного количества тепла. Поэтому тепловой поток на границе фазового перехода разрывен и определяется величиной
дх
где Q - теплота фазового перехода, а ¥п - скорость движения границы фазового перехода по нормали.
Как мы уже выше отмечали, предполагается, что фазовый переход происходит при постоянной температуре и *. Поэтому сама граница фазового перехода 8 определяется на каждый момент времени следующим образом: 8 = 5 (г) = {( х, у, і) є О, и(х, у, і, г) = и *} или, в другой форме, на границе фазового перехода выполнены условия первого рода
и(х, у, і, г) = и *, (х, у, і) є 8(г). (2.5)
Условия (2.3)-(2.5) есть условия Стефана, а соответствующая задача для уравнений (2.1), (2.2) называют задачей Стефана. Рассматриваемая задача
характеризуется тем, что исследуются процессы в обеих фазах, поэтому в этом случае говорят о двухфазной задаче Стефана. Предельная ситуация характеризуется тем, что тепловое поле в одной из фаз известно (температура равна температуре фазового перехода). Поэтому рассматривается тепловое поле только для одной из фаз — однофазная задача Стефана [6]. В этом случае неизвестная граница фазового перехода 8 является не внутренней, а внешней.
Например, будем считать, что область О занята твердой фазой при температуре и *. Тогда для нахождения температуры в жидкой фазе используется уравнение (2.2) в переменной области О+ (г) со следующими условиями на 8 :
и + (х, у, і, г) = и *, (х, у, і) є 8(г). (2 6)
ди+
Л+^— = Ч2К, (х,у,і)єБ(г). (2.7)
дх
Условия типа (2.6), (2.7) характеризуют однофазную задачу Стефана.
С точки зрения построения эффективных вычислительных алгоритмов чрезвычайно важное значение имеет тот факт, что задача Стефана допускает обобщенную формулировку, при которой условия (2.3)-(2.5) включаются в само уравнение теплопроводности. Включение неоднородных условий сопряжения на заданной границе раздела в само уравнение обсуждалось выше. В задаче Стефана ситуация осложняется тем, что граница раздела фаз S сама неизвестна, и ее нужно найти.
Поясним переход от уравнений (2.1), (2.2) с условиями (2.3)-(2.5) к одному уравнению теплопроводности. В соответствии с (2.4) и (1.13)-(1.14) уравнения (2.1), (2.2) можно записать в виде одного уравнения (ди ^
cp\-----+ v■ grad и I = div(A grad и)-8SQV + s , (x,y,z,t) e 0 . (2.8)
\dt J
Вблизи границы фазового фронта введем локальную ортогональную координатную систему (x',y',z'), метрические коэффициенты которой равны единице, В
этих новых координатах поверхностная 8 -функция 8S есть 8S =8(x' - x0), где уравнение
x ' = x0 определяет границу S . Аналогично для скорости движения свободной границы
имеем Vn = dx'/dt. Условие Стефана (2.6) соответствует тому, что в новых координатах
и = u(x\y', z', t), и (x0, y', z') = и *. С учетом этого получим
8sVn =8( x - x0) ^ = 8(и - и*) . (2.9)
dt dt
Подстановка (2.9) в (2.8) дает уравнение
ди
(cp + Q8(u - и*))— = div(A grad и) + s, (x, y, z, t) e0 , (2.10)
dt
Уравнение (2.9) примечательно тем, что в него не входит явно сама неизвестная граница фазового перехода. Фронт фазового перехода в такой постановке задачи находится как изотермическая поверхность и = и* = const, положение которой в пространстве, а в общем случае и форма, изменяются с течением времени. Учет теплоты фазового перехода эквивалентен заданию эффективной теплоемкости: ceff = c + plQ8 (и - и*).
Квазилинейное уравнение теплопроводности (2.10) лежит в основе эффективных процедур приближенного решения задач Стефана. Нелинейности в этой задаче изменяют не только количественный характер тепловых процессов, но и качественную картину их протекания. Они значительно усложняют математические модели тепловых процессов,
причем во многом эти трудности связаны с невозможностью применения для нелинейных задач принципа суперпозиции решений.
В общем виде (2.10) такая задача не имеет полного аналитического решения. Однако найдены точные аналитические решения для некоторых частных случаев (например, для постоянных коэффициентов и для коэффициентов, зависящих только от температуры). Подробнее с ними можно ознакомиться в [3-4, 7-11].
3 Расчетная часть
Требуется найти распределение температуры в среде с фазовыми переходами твердое тело - жидкость. Тепловое состояние такой среды с учетом теплоты фазового перехода описывается уравнением теплопроводности вида
[с(х, у, г, и)р(х, у, г,и) + QS(u - и*)]дИ = ^^(х, У,2,и) +
д ( ди ^ д ( ди ^
+----Я(х, У,г, и)----- +-|Я(х, У, г, и)-------- 1 + ^(х, У, г, і),
ду ^ ду) дг \ дг )
(3.1)
где с - удельная теплоемкость; р - плотность; Я - коэффициент теплопроводности; и(х, у, г, і) - температура среды; и * - температура фазового перехода; Q - теплота фазового перехода; ^ - мощность внутренних источников тепла; 5 (и - и*) - дельтафункция.
Решение и требуется найти в ограниченной области (рис. 1) - прямоугольном параллелепипеде Б = {0 < х <Ьх, 0 <у <Ьу, 0 < г <Ьг}, удовлетворяющее начальному
условию и (х, у, г,0) = ф( х, у, г). На границе 2 = 0 с температурой происходит конвективный теплообмен со средой, имеющей температуру 0(і). Плотность теплового потока на этой границе задается через коэффициент теплоотдачи к в виде
ди
Я — = к(в(і) - и(х, у,0, і)).
дг
Коэффициент теплоотдачи h записывается в виде h = (а 1 + R) 1, где а -коэффициент конвективного теплообмена и R - коэффициент термического сопротивления.
На границе z = Lz поддерживается постоянный тепловой поток (поток из недр Земли согласно [12])
Л— = а„ = const.
dz E
Боковые границы области D теплоизолированы:
дич _ ди ,т . Л
— (О, y,z,t) = 0, —(Lx, y,z,t) = 0, dx dx
£ (x,0, z, t) = ^ ^ (x, Ly, z, t) = °.
dy dy
В качестве примера расчета рассматривается задача о распределении температурных полей в грунте в районах с сезонными промерзаниями и оттаиваниями почвы (аналогично РСН 67-87 ([29]) и [30]). Рассмотрим здание со свайным фундаментом и термостабилизаторами. Такой тип фундамента часто встречается при строительстве в зоне вечномерзлых грунтов. Предположим, что сверху на грунте построено здание. Необходимо найти распределение температурных полей под зданием, чтобы узнать глубину оттаивания грунтов под ним. Если под зданием через некоторый промежуток
времени будут скапливаться оттаявшие грунты, здание может осесть в котлован, образовывающийся под тяжестью самого здания, а затем и вовсе разрушиться [31].
Термостабилизатор - это парожидкостное устройство для охлаждения грунтов, представляющее собой металлическую герметично запаянную, заправленную хладагентом трубку диаметром от 36 до 57 мм, длиной от 6 до 10 м и более, состоящее из конденсатора с оребрением (надземной части длиной в пределах 1—2,5 м) и испарителя (подземной части длиной от 5 до 9 ми более). Работа осуществляется без внешних источников питания, только за счет законов физики — переноса тепла вследствие испарения в испарителе хладагента и его поднятия в конденсаторную часть, где пар конденсируется, отдавая тепло, и стекает по внутренним стенкам трубы вниз. В качестве хладагента применяется экологически безвредный, не воспламеняющийся и взрывобезопасный фреон R22. В расчетах будем использовать гладкостенные термостабилизаторы СОУ СГВ-100-40/9 с глубиной погружения в грунт г0 = 8 м., общей длиной 9 ми мощностью
ккал
50 = 28.3932----- [32, 35]. Термостабилизатор работает с 1 октября по 31 марта и не
м • ч
работает с 1 апреля по 30 сентября каждого каждого года. Термостабилизатор рассматривается как бесконечно тонкий стержень:
5( х ^ Z, г) = ^(г) • 3( х - х0) •5( у - У0) • /[0,г0] (г) .
В задаче мы будем рассматривать самый простой тип термостабилизаторов: прямой, без изменения формы для теплоотвода. У здания нет вентилируемого подполья, поэтому термостабилизаторы находятся вне границы здания.
Пусть область исследования Д имеет вид, представленный на рис. 2. Плоскость Оху соответствует поверхности земли. Ось Ог направлена вглубь грунта. Закрашенный прямоугольник на верхней границе области соответствует четверти здания,
расположенному над этой областью. Размер этого прямоугольника (горизонтального сечения здания) составляет 10мх4м (рис. 3.1). Характеристики области Д: Ьх =21 м,
Ь =15 м, Ь =17 м. Область Д состоит из 3 литологических слоев: Д, Д и Д
у 5 г 1 ’ 2 3
(соответственно, области 1, 2 и 3 на рисунке 1). Эти множества задаются следующим образом:
Д1 = (УМ(х, у, г) е Д : 0 < г < 2};
Д2 = (УМ(х, у, г) е Д : 6 < х < 21, 2 < г < 6} ;
Д3 = Д \ (Д1 и Д2).
Область исследования Д содержит 6 свай и 5 термостабилизаторов, расположенных как показано на рис. 2.
Рис. 2. Верхняя граница расчетной области
В зоне 2 (на рис. 2 - это заполненный текстурой «кирпичная стена» прямоугольник) располагается производственное здание, физические характеристики в этой зоне остаются постоянными во времени. Условия в зонах 1 и 3 изменяются по интервалам времени - месяцам года (см. таблицу 1). За начало отсчета принимается 01 января текущего года. Через год эти условия повторяются, т.е. условия на верхней границе являются периодическими функциями времени с периодом Т = 8760 ч (1 год). Боковые границы области теплоизолированы, и условия на них не зависят от времени. Поток на боковых границах в плоскостях х = 0 и у = 0 считаются равным нулю в силу симметрии (см. п. 3.1.1), а на двух других плоскостях равен нулю, т.к. эти грани удалены достаточно далеко и не оказывают влияния на распределение температурных полей под зданием. Начальное распределение температуры по глубине представлено в таблице 3.3.
Физические характеристики на верхней границе расчетной области
Зоны Характерис- тика Месяц
I II III IV V VI VII VII I IX X XI XII
1 Температура в ,°С -24,6 -23,8 -19,2 -9,5 3,7 13,1 18,0 12,4 4,6 -4,6 -16,4 -22,4
Термическое сопротивление Я, м2-ч-°С/ккал 2,6 2,7 2,8 3,0 1,7 - - - - 1,3 2,1 2,3
Коэффициент конвективного теплообмена а, ккал/(м2-ч-°С) 13,8 12,3 12,3 12,0 12,1 12,6 13,4 13,0 14,3 13,4 13,5 12,7
2 Температура в, °С в =20°С= соті
Термическое сопротивление я, м2-ч-°С/ккал Я=1,8= соті
Коэффициент конвективного теплообмена а, ккал/(м2-ч-°С) а =20.0= соті
3 Температура в, °С -24,6 -23,2 -19,2 -9,5 3,7 13,1 18,0 12,4 4,6 -4,6 -16,4 -22,4
Термическое сопротивление я, м2-ч-°С/ккал 4,4 4,4 4,7 5,0 2,82 - - - - 2,14 3,57 3,59
Коэффициент конвективного теплообмена а, ккал/(м2-ч-°С) 13,8 12,3 12,3 12,0 12,1 12,6 13,4 13,0 14,3 13,4 13,5 12,7
Согласно РСН 67-87, СНиП 2.02.04-88, и ГОСТ 25100-82 ([29], [32] и [33] соответственно) можно написать для слагаемого cp в левой части уравнения (3.1) выражение
сР=
рй ^ Сй + cice (Wtot - Ww ) + CwWw + k
Ра (са + ),
дW
ди
w 1 и < и*;
(3.2)
и > и *.
г, л ккал л ккал
Здесь cd = 0,22----------удельная теплоемкость сухого грунта, cice = 0,49---------
кг • °С кг• °С
ккал
удельная теплоемкость льда, cw = 1,006--------- - удельная теплоемкость воды, pd -
кг •°С
плотность сухого грунта, Wtot - суммарная влажность грунта в долях к массе сухого грунта, Ww = Ww (и) - доля незамерзшей воды по отношению к массе сухого грунта,
к = 79,4 ккал - удельная теплота фазового перехода льда, и* = 0 - температура фазового
кг
перехода.
Формула для объемной теплоемкости (3.2) учитывает фазовые переходы в области отрицательных температур. Функцией Ww = Ww (и) описывается массовая доля (по
отношению к массе сухого грунта) незамершей воды при температуре и [29, 34]. Подробнее нахождение этой функции описано ниже.
Коэффициент теплопроводности Я в уравнении (3.1) определяется выражением
[Лы, и < и*,
ЯЧ M (3.3)
[Ят, и > и*,
где Ям - коэффициент теплопроводности мерзлого грунта, Ят - коэффициент теплопроводности талого грунта.
Коэффициент теплоотдачи h записывается в виде h = (а-1 + R) 1, где значения коэффициента конвективного теплообмена а и термического сопротивления R находятся по таблице 1 (см. [29]).
Теплота фазового перехода Q принимается равной количеству теплоты, необходимой для таяния льда (замерзания воды) в единице объема грунта и вычисляется по формуле Q = kPd (Wtot - Ww (и*)) в соответствии со СНиП 2.02.04-88 [32]. Внутренние источники тепла в рассматриваемой задаче отсутствуют, поэтому s(x, y, z, t) = 0 .
Физические характеристики литологических слоев, необходимые для расчета, приведены в таблице 2 (см. РСН 67-87, ГОСТ 25100-82, СП 32-101-95 [29, 33, 35]).
На нижней границе области z = Lz (см. рис. 2) зададим поток из недр Земли [12]
„ „ ккал qE = 0,043-^.
м ч
Физические характеристики литологических слоев
Номер слоя Характеристика
Р, кг/м3 сй. ккал/(кг-°С) W^ot, доли аи. ккал/(м-ч-°С) Ат. ккал/(м-ч-°С) Ww (и), доли
° О Ww (и)
-0,3 0,14
1 1390 0,22 0,25 1,3 1,15 -1,0 0,12
-10,0 0,08
-0,3 0,12
2 1520 0,24 0,22 1,55 1,55 -1,0 0,08
-10,0 0,05
0 0
3 1500 0,23 0,27 2,35 2,15 0 0
0 0
Верхняя граница области г = 0, представленная на рисунке 3.4, разбита на 3 зоны:
1, 2 и 3. Будем рассматривать режим теплообмена на прилегающей территории, наименее благоприятный для эксплуатации здания. Предположим, что снежные надувы образовывались симметрично со всех сторон здания. В результате определялся температурный режим под четвертью частью здания, образованной осями симметрии. В процессе эксплуатации здания с одной из его сторон на части прилегающей территории в зимнее время будут образовываться снежные надувы (увеличение максимальной мощности снега на 0,4 м по сравнению с естественным). Зона снежных надувов обозначена на рис. 2 цифрой 3.
Таблица 3
Начальное распределение температуры по глубине
Г лубина, м 1 3 5 7 9,5 12,5 15,5
°С -0,1 -0,4 -1,0 -1,2 -1,4 -1,4 -1,4
Дополнительно к разбиению области расчета на блоки, в каждом блоке по всем трем направлениям х, у иг вводилась равномерная сетка: длины блоков по каждому направлению делились на 5 частей. Таким образом, в результате число контрольных объемов по х, у и г составило, соответственно, 40, 30 и 35, т.е. общее число элементов в
расчетной области (это же число равняется числу внутренних расчетных точек в области Б ) равняется 40 х 30 х 35 = 42000. В ходе методических расчетов проводилось сравнение вариантов разбиения каждого блока на 5, 7, 9 и 11 частей. Установлено, что разбиение на 5 частей позволяет с точностью до 10"5 находить распределение температур. Остальные варианты не дают ощутимый выигрыш в точности, в то время как эти варианты существенно снижают время расчетов. Минимальный безразмерный шаг по пространству
Актт =^ 2^^ = 0,019, АЛтах =^3^^ = 0,029. В качестве рассматриваемого интервала времени выбран интервал 5 лет = 5 • 8760ч = 43800 ч . Число итераций по времени выбрано равным 6000. Тогда безразмерный шаг по времени Аt = ^ 43800 ) Т = 0,00083.
2
Безразмерная полуширина сглаживания 8 -функции А =-----------------= 0,08. В качестве
^тах
показателя даты выбран стандартный формат «дд.мм.гг», где годы будет иметь вид 01, 02,
03, 04 или 05.
На вертикальных линиях, проходящих параллельно оси Ог через точки наблюдения М1(1.5,1), М2(9,5) (см. рис. 2), выведены на печать в определенных точках по глубине зависимости температуры и в этих точках от времени t (рисунки 4-5). По графикам наглядно видно, что в точках М1 и М2 термостабилизаторы постепенно снижают температуру.
Видно, что к концу рассматриваемого периода (5 лет) температура под зданием становится отрицательной уже на глубине порядка 2 м. Это означает, что большая часть грунта находится в твердом состоянии, и зданию не угрожает разрушение. Для свай глубиной залегания 8 м. это вполне допустимая глубина промерзания. Три четверти длины сваи находится в твердом грунте, что позволит зданию устоять без просадки. Хотя следует отметить, что для большей безопасности необходимо использовать термостабилизаторы строго под зданием рядом со сваями. Таким образом, можно подморозить грунт в непосредственной близости от свай и почти полностью сваи будут находиться в твердом грунте. Это особенно важно при строительстве крупных зданий. Также есть смысл строительства зданий с вентилируемым подпольем, за счет чего значительно уменьшится эффект накопления тепла под зданием. В данном случае граничные условия не будут принимать постоянное значение. Однако построение такой модели требует учесть движение воздушных масс и их прогревание, а также их
специфический теплообмен с дисперсной средой. По данному направлению ведутся исследования, достигнутые результаты изложены в [37].
Рис. 3. Зависимость температуры от времени в точке М1 на глубинах: 1 — 1 м., 2 - 3 м., 3 —
5 м., 4 — 7 м., 5 — 15,5 м.
и,°с | | | | | Г
12______I_____I_____I_____I______I_____I_____I_____I_____I_____
О 0.5 1 1.5 2 2.5 3 3.5 4 4.5 I. год
Рис. 4. Зависимость температуры от времени в точке М2 на глубинах: 1 — 1 м., 2 - 3 м., 3 —
5 м., 4 — 7 м., 5 — 15,5 м.
Выводы
В результате выполнения данного проекта была разработана эффективная математическая модель для определения температурных полей в низкотемпературных двухфазных средах. В ходе работы были получены результаты прогноза распределения температурных полей для зданий с разным типом фундамента: свайный, ленточный и без фундамента; рассмотрен вариант охлаждения грунта под зданием с помощью термостабилизаторов.
Кроме того была создана программа, которая реализует численное решение класса трехмерных обобщенных задач Стефана в безразмерном виде. Программа позволяет прогнозировать распределение температурных полей в низкотемпературных двухфазных средах во времени. Разработанная программа может найти применение в таких отраслях промышленности как космонавтика, геология, строительство в зоне вечной мерзлоты и др.
Список литературы
1 Крылов Д.А., Мельникова Ю.С. Математическое моделирование распределения температурных полей в криолитозоне // Студенческий научный вестник. Сборник статей
четвертой научно-технической выставки «Политехника». М.: Изд-во МГТУ им.
Н.Э. Баумана, 2009. С.94-97.
2 Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. М.: Книжный дом «ЛИБРОКОМ», 2009. 784 с.
3 Карслоу Г., Егер Д. Теплопроводность твердых тел. М.: Наука, 1964. 488 с.
4 Лыков А.В. Теория теплопроводности. М.: Высшая школа, 1967. 599 с.
5 Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах. Т.1. - М.: Мир, 1991. - 504 с.
6 Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах. Т.2. М.: Мир, 1991. 552 с.
7 Тихонов А.Н., Самарский А.А. Уравнения математической физики: учеб. пособие. М.: Изд-во МГУ, 1999. 799 с.
8 Мартинсон Л.К., Малов Ю.И. Дифференциальные уравнения математической физики: учебник для студентов вузов. М.: Изд-во МГТУ им. Н.Э. Баумана, 1996.
9 Фридман А. Вариационные принципы и задачи со свободными границами. М.: Наука, 1990. 536 с.
10 Власова Е.А., Зарубин В.С., Кувыркин Г.Н. Приближенные методы математической физики: учеб. для вузов. М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. 704 с.
11 Синай Я.Г. Теория фазовых переходов. Строгие результаты. М.: Наука, 1980. 208 с.
12 Короновский Н.В. Общая геология: учебник. М.: Изд-во «Книжный дом университет», 2006. 528 с.
13 Седов Л.И. Методы подобия и размерности в механике. М.: Наука, 1977. 440 с.
14 Зарубин В.С. Математическое моделирование в технике: Учеб. для вузов. М.: Изд-во МГТУ им. Н.Э. Баумана, 2001. 496 с.
15 Ferziger Joel H., Petric M. Computational methods for fluid dynamics. 3 rev. ed. - New York: Springer, 2002. - 423 p.
16 Самарский А.А., Гулин А.В. Численные методы: Учеб. пособие для вузов. - М.: Наука, 1989. - 432 с.
17 Галанин М.П., Савенков Е.Б. Методы вычислений. Уравнения в частных производных. Интегральные уравнения: Учеб. пособие. М.: Изд-во ИПМ им. М.В. Келдыша РАН, 2006. 72 с.
18 Аттетков А.В., Галкин С.В., Зарубин В.С. Методы оптимизации: Учеб. для вузов. М.: Изд-во МГТУ им. Н. Э. Баумана, 2003. 440 с.
19 Корн Г., Корн Т. Справочник по математике (для научных работников и инженеров). Определения, теоремы, формулы. Спб.: Изд-во «Лань», 2003. 832 с.
20 Галанин М.П., Савенков Е.Б. Методы вычислений. Обыкновенные дифференциальные уравнения. Разностные схемы: Учеб. пособие. М.: Изд-во ИПМ им. М.В. Келдыша РАН, 2006. 90 с.
21 Отчет по теме «Фундаментальные основы управления нестационарным температурным режимом в криолитозоне» / МГТУ. Руководитель темы Н.И. Сидняев. ГР № 01201151171, Инв. № 02201157901. М., 2010, 306 с.
22 Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.
23 Патанкар С.В. Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах. М.: Изд-во МЭИ, 2003. 312 с.
24 Веденеева Е.А. Применение метода контрольного объема к решению задач динамики вязких теплопроводных жидкостей. М.: Изд-во МГУ, 2009. 87 с.
25 Сборник задач по методам вычислений: Учеб. пособие для вузов / Под ред. П.И. Монастырского. М.: Физматлит, 1994. 320 с.
26 Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.: Высшая школа, 1994. 544 с.
27 Галанин М.П., Савенков Е.Б. Методы вычислений. Алгебра. Интерполяция: Учеб. пособие. М.: Изд-во ИПМ им. М.В. Келдыша РАН, 2006. 80 с.
28 РСН 67-87. Инженерные изыскания для строительства. Составление прогноза изменений температурного режима вечномерзлых грунтов численными методами. М., 1987.
30 Крылов Д.А., Сидняев Н.И., Федотов А.А. Математическое моделирование распределения температурных полей в вечномерзлых грунтах // Материалы Международной научно-практической конференции по инженерному мерзлотоведению, посвященной ХХ-летию создания ООО НПО «Фундаментстройаркос»» - Тюмень. Изд-во «Сити-пресс», 2011. С. 126-132.
31 Цытович Н.А. Механика мерзлых грунтов. М.: Высшая школа, 1973. 448 а
32 СНиП 2.02.04-88. Основания и фундаменты на вечномерзлых грунтах. М., 1998.
33 ГОСТ 25100-82. Грунты. Классификация. М., 1982.
34 Роман Л.Т. Механика мерзлых грунтов. М.: МИАК «Наука/Интерпериодика», 2002. 426 с.
35 СП 32-101-95. Проектирование и устройство фундаментов опор мостов в районах распространения вечномерзлых грунтов. М., 1996.
36 Попов А.П., Милованов В.И. и др. К вопросу о типовых технических решениях по основаниям и фундаментам для криолитозоны // Инженерная геология. 2008. Сентябрь. -С. 22-38.
37. Крылов Д.А., Сидняев Н.И. Метод расчета массовой кристаллизации многофазных реологических сред // Материалы Четвертой конференции геокриологов России. МГУ им. М.В. Ломоносова, 7-9 июня 2011 г. Т. 1. Часть 1. Физико-химия, теплофизика и механика мерзлых пород. М.: Университетская книга, 2011. С. 129-136.
electronic scientific and technical periodical
SCIENCE and EDUCATION
__________El. .Vs KS 77 - 30569. -V«042l 100025. ISSN 1994-0408_
Mathematical simulation of temperature fields considering phase transitions in cryolithic zone 77-30569/354740
# 04, April 2012 Krylov D.A.
Bauman Moscow State Technical University [email protected]
One of the most dangerous processes that reduce stability of constructions in permafrost is seasonal soil thawing and heat accumulation beneath buildings. This paper describes a threedimensional mathematical model for analyzing and forecasting temperature distribution under buildings subject to seasonal temperature fluctuations. An example of calculation of a temperature field under a building is included in the article. The developed model and software package can be used for geotechnical investigations at any stage of construction projects in the permafrost zone.
Publications with keywords: mathematical model, ground, the forecast, temperature field, thermostabilizer Publications with words: mathematical model, ground, the forecast, temperature field, thermostabilizer
References
1. Krylov D.A., Mel'nikova Iu.S. Matematicheskoe modelirovanie raspredeleniia temperaturnykh polei v kriolitozone [Mathematical modeling of distribution of temperature fields in the permafrost zone]. Studencheskii nauchnyi vestnik. Sbornik statei chetvertoi nauchno-tekhnicheskoi vystavki «Politekhnika» [Student scientific Journal. Collection of articles of the 4th Science and Technology Exhibition "Politehnica"]. Moscow, Bauman MSTU Publ., 2009, pp. 94-97.
2. Samarskii A.A., Vabishchevich P.N. Vychislitel'naia teploperedacha [Computational heat transfer]. Moscow, Knizhnyi dom «LIBROKOM», 2009. 784 p.
3. Karslou G., Eger D. Teploprovodnost' tverdykh tel [Thermal conductivity of solids]. Moscow, Nauka, 1964. 488 p.
4. Lykov A.V. Teoriia teploprovodnosti [Theory of heat conduction]. Moscow, Vysshaia shkola, 1967. 599 p.
5. Fletcher K. Vychislitel'nye metody v dinamike zhidkostei: V2-kh tomakh. T.1 [Computational approach in Fluid Dynamics. In 2 vols. Vol. 1]. Moscow, Mir, 1991. 504 p.
6. Fletcher K. Vychislitel'nye metody v dinamike zhidkostei: V2-kh tomakh. T.2. [Computational approach in Fluid Dynamics. In 2 vols. Vol. 1]. Moscow, Mir, 1991. 552 s.
7. Tikhonov A.N., Samarskii A.A. Uravneniia matematicheskoifiziki [The equations of mathematical physics]. Moscow, MSU Publ., 1999. 799 p.
8. Martinson L.K., Malov Iu.I. Differentsial'nye uravneniia matematicheskoi fiziki [Differential equations of mathematical physics]. Moscow, Bauman MSTU Publ., 2006. 368 p.
9. Fridman A. Variatsionnyeprintsipy i zadachi so svobodnymi granitsami [Variational principles and free boundary problems]. Moscow, Nauka, 1990. 536 p.
10. Vlasova E.A., Zarubin V.S., Kuvyrkin G.N. Priblizhennye metody atematicheskoi fiziki [Approximate methods of mathematical physics]. Moscow, Bauman MSTU Publ., 2004. 704 p.
11. Sinai Ia.G. Teoriia fazovykh perekhodov. Strogie rezul'taty [Theory of phase transitions. Rigorous results]. Moscow, Nauka, 1980. 208 p.
12. Koronovskii N.V. Obshchaiageologiia [General Geology]. Moscow, «Knizhnyi dom universitet», 2006. 528 p.
13. Sedov L.I. Metody podobiia i razmernosti v mekhanike [Methods of similarity and dimensionality in mechanics]. Moscow, Nauka, 1977. 440 p.
14. Zarubin V.S. Matematicheskoe modelirovanie v tekhnike [Mathematical modeling in engineering]. Moscow, Bauman MSTU Publ., 2001. 496 p.
15. Ferziger Joel H., Petric M. Computational methods for fluid dynamics. New York,
Springer, 2002. 423 p.
16. Samarskii A.A., Gulin A.V. Chislennye metody [Numerical methods]. Moscow, Nauka, 1989. 432 p.
17. Galanin M.P., Savenkov E.B. Metody vychislenii. Uravneniia v chastnykhproizvodnykh. Integral'nye uravneniia [Methods of computation. Equations in partial derivatives. Integral equations]. Moscow, IPM im. M.V. Keldysha RAN Publ, 2006. 72 p.
18. Attetkov A.V., Galkin S.V., Zarubin V.S. Metody optimizatsii [Optimization procedure]. Moscow, Bauman MSTU Publ., 2003. 440 p.
19. Korn G., Korn T. Mathematical Handbook for Scientists and Engineers: Definitions,
Theorems, and Formulas for Reference and Review. New York, Dover Publications, 2000. 943 p. (Russ. ed.: Korn G., Korn T. Spravochnikpo matematike (dlia nauchnykh rabotnikov i inzhenerov). Opredeleniia, teoremy, formuly. Spb., Lan', 2003. 832 p.).
20. Galanin M.P., Savenkov E.B. Metody vychislenii. Obyknovennye differentsial'nye uravneniia. Raznostnye skhemy [Methods of computation. Ordinary differential equations. Difference schemes]. Moscow, IPM im. M.V. Keldysha RAN Publ., 2006. 90 p.
21. Sidniaev N.I. Otchetpo teme «Fundamental'nye osnovy upravleniia nestatsionarnym temperaturnym rezhimom v kriolitozone» [Report on the theme "Fundamental principles
for management of the non-stationary temperature regime in the permafrost zone"]. Moscow, MGTU im. N.E. Baumana [Bauman MSTU], 2010. 306 p. Not published.
22. Patankar S. Chislennye metody resheniia zadach teploobmena i dinamiki zhidkosti [Numerical methods of solving tasks of heat transfer and fluid dynamics]. Moscow, Energoatomizdat, 1984. 152 p.
23. Patankar S.V. Chislennoe reshenie zadach teploprovodnosti i konvektivnogo teploobmena pri techenii v kanalakh [Numerical solution of heat conduction and convection heat transfer for flow in channels]. Moscow, MEI Publ., 2003. 312 p.
24. Vedeneeva E.A. Primenenie metoda kontrol'nogo ob"ema k resheniiu zadach dinamiki viazkikh teploprovodnykh zhidkostei [Application of the method of control volume to the decision of problems of dynamics of viscous heat conducting fluids]. Moscow, MSU Publ., 2009. 87 p.
25. Monastyrskii P.I., ed. Sbornik zadach po metodam vychislenii [Collection of problems on methods of calculation]. Moscow, Fizmatlit, 1994. 320 p.
26. Amosov A.A., Dubinskii Iu.A., Kopchenova N.V. Vychislitel'nye metody dlia inzhenerov [Computational approach for engineers]. Moscow, Vysshaia shkola, 1994. 544 p.
27. Galanin M.P., Savenkov E.B. Metody vychislenii. Algebra. Interpoliatsiia [Methods of computation. Algebra. Interpolation]. Moscow, IPM im. M.V. Keldysha RAN Publ., 2006. 80 p.
28. Republ. Build. Code 67-87. Moscow, Gosstroi RSFSR, 1987.
29. Krylov D.A., Sidniaev N.I., Fedotov A.A. Matematicheskoe modelirovanie raspredeleniia temperaturnykh polei v vechnomerzlykh gruntakh [Mathematical modeling of the distribution of temperature fields in the permafrost soils]. Materialy Mezhdunarodnoi nauchno-prakticheskoi konferentsii po inzhenernomu merzlotovedeniiu, posviashchennoi KhKh-letiiu sozdaniia OOO NPO «Fundamentstroiarkos» [Proc. of the International Scientific and Practical Conference on Permafrost Engineeringdevoted to 20th anniversary of the establishment of a Limited Liability Company Scientific Production Association "Fundamentstroyarkos"]. Tiumen', Siti-press, 2011, pp. 126-132.
30. Tsytovich N.A. Mekhanika merzlykh gruntov [The mechanics of frozen soils]. Moscow, Vysshaia shkola, 1973. 448 p.
31. Construction Norms and Regulations 2.02.04-88. Moscow, 1998.
32. The State Standard 25100-82. Moscow, 1982.
33. Roman L.T. Mekhanika merzlykh gruntov [The mechanics of frozen soils]. Moscow,
MIAK «Nauka/Interperiodika», 2002. 426 p.
34. Sets of Rules to construction 32-101-95. Moscow, 1996.
35. Popov A.P., Milovanov V.I., et al. K voprosu o tipovykh tekhnicheskikh resheniiakh po osnovaniiam i fundamentam dlia kriolitozony [On the question of typical technical solutions on the basis and foundation for cryolithozone]. Inzhenernaia geologiia, 2008, September, pp. 22-38.
36. Krylov D.A., Sidniaev N.I. Metod rascheta massovoi kristallizatsii mnogofaznykh reologicheskikh sred [The method of calculating of mass crystallization of multiphase rheological environments]. Materialy Chetvertoi konferentsii geokriologovRossii. MGU im. M.V. Lomonosova, 7-9 iiunia 2011 g. T. 1. Chast' 1. Fiziko-khimiia, teplofizika i mekhanika merzlykh porod [Proc. of the Fourth Conference of geocryologists Russia. Lomonosov Moscow State University, 7-9 June 2011. Volume 1. Part 1. Physico-chemical, thermal physics and mechanics of frozen soils]. Moscow, Universitetskaia kniga, 2011, vol.
1, pt. 1, pp. 129-136.
