CC BY
19
УДК 625.711.84+625.31 А.С. Миляев
Санкт-Петербургский государственный лесотехнический университет
Миляев Александр Сергеевич родился в 1936 г., окончил в 1959 г. Ленинградское высшее военное инженерно-техническое училище ВМФ, доктор технических наук, профессор, засл. работник высшей школы РФ, зав. кафедрой теоретической и строительной механики С.-Петербургского государственного лесотехнического университета. Имеет более 150 научных работ в области механики деформируемого твердого тела, в том числе механики силового взаимодействия конструкций и сооружений с грунтом при статических и динамических нагрузках.
E-mail: icffi@home.ru
АЛГОРИТМ РАСЧЕТА ПРОМЕРЗАНИЯ СЛОИСТЫХ ОСНОВАНИЙ ЗИМНИХ ЛЕСОВОЗНЫХ ДОРОГ
Приводится алгоритм расчета глубины и скорости промерзания слоистых оснований зимних лесовозных дорог на базе текущих значений температуры воздуха и скорости ветра в данной местности. Алгоритм построен на численном решении нелинейной нестационарной задачи о теплопередаче в термодинамической системе «воздух -грунт» с учетом фазовых превращений поровой воды в лед, а также зависимости теп-лофизических характеристик грунтов от температуры.
Ключевые слова: зимние лесовозные дороги, теплопередача, промерзание оснований, фазовые превращения воды в лед, нелинейное нестационарное уравнение теплопроводности.
При создании блок-схемы алгоритма расчета промерзания оснований зимних лесовозных дорог рассмотрена замкнутая термодинамическая система (рис. 1), состоящая из слоев воздуха и грунта с известными физико-механическими и теплофизическими свойствами. Начальная температура и фазовое состояние системы предполагаются известными; температура воз-духа изменяется по заданному закону; поверхности y = 0, y = b, x = h4 тепло-изолированы. Требуется определить изменение температуры в системе в зависимости от координат, времени и изменения температуры воздуха и условий на ее границе.
Распространение теплоты в слоях грунта подчиняется известному уравнению теплопроводности [1]:
дИ _дТ_
дд ~ дх Iх дх
д ду
дТ д д Т
\ — +— X —
дх z ду _
qv,
(1)
где
Хэ
Н - энтальпия, Дж/м3; ^ - время, с;
Х2 - теплопроводность в направлении осей х, у, г, Вт/м2; Т - температура, °С;
ду - мощность внутренних источников энергии, Вт/м3.
© Миляев А.С., 2012
Рис. 1. Общая расчетная схема: 1-4 -слои грунтового основания; к1 - к4 -отметки подошв слоев; ТаГ - температура воздуха на поверхности основания,
Тагт Тагг(£)
Тепловые потоки через границы у = 0, у = Ъ, х = к4 отсутствуют:
О(х, у, 2, 0 = 0. Приращение энтальпии определяется выражением
Т
АН = | р С (Т)ёТ, (2)
Т
где ТЪ, Т - начальное и текущее значения температуры, °С; р - плотность грунта, кг/м3; С - удельная теплоемкость, Дж/(кг°С).
Внутренних источников энергии в системе нет = 0), однако в слоях грунта при замерзании поровой воды выделяется теплота кристаллизации (2ЯГ, Дж/м3), которая учитывается в расчете по СНиП 2.02.04-88 [5]:
= - жт)р», (3)
где О - теплота кристаллизации воды, О = 3,35 105 Дж/кг;
- суммарная влажность грунта, доли единицы;
Жт - влажность мерзлого грунта между включениями льда, доли единицы;
- плотность скелета грунта, кг/м3.
Принимаем, что выделяется в интервале температур от —1 до 0 °С; Wm = 0.
На границе х = 0 происходит конвективный теплообмен по следующему закону:
й = аЛ(Тшг - ТД (4)
где О — тепловой поток через границу х = 0;
а — коэффициент теплоотдачи (теплопритока); Л — площадь теплообмена (теплоотдачи); Т^ - температура грунта на его поверхности.
Коэффициент а на границе «воздух-снег» и «воздух-лед» вычисляется по формулам СНиП 2.06.04-82 [4].
Уравнение (1) решается численным методом конечных элементов (КЭ). Для обеспечения устойчивости численной процедуры необходимо соблюдать определенное соотношение между «шагами» по времени Дt и пространству Дк [5]:
ХМ 1
-< —
2
р С (АН )2
(5)
где X — теплопроводность в направлении распространения тепловой волны; Дк - приращение пространственной координаты в направлении распространения тепловой волны. Вариационное уравнение теплопроводности (1) с присоединенными граничными условиями для выделенного объема V имеет следующий вид:
5Т Г^ 1 + {ДТт8(5Т )([Л] {Ь}Т
= 1ЗТд * dS2 + | ЗТаЛ(Та1Г - Т)<^3, (6)
где
5Т — вариация температуры; Тгп8 — оператор транспонирования; Ч
* — тепловой поток через поверхность 52, Ч* = — {д}Тпк {п}; {Ч} = —[Л]{Ь}Т;
Ч — вектор теплового потока через поверхность 52; (п} — единичный вектор внешней нормали к поверхности 52; аА (Тсаг — Т) — конвективный тепловой поток через поверхность 53, оА(Та1Г — Т) = (п}Тппк [Л]{Ь}Т; [Л] — локальная матрица теплопроводности,
[Л] =
Х х (х, у, 2,Т) 0 0
0 Х у (х, у, г, Т) 0
0 0
Хг (х, у, 2,Т) {Ь} — оператор дифференцирования,
Тгпв
И = .
[дх ду дг ]
Уравнение теплопроводности для одного конечного элемента (КЭ) получается при интегрирования уравнения (6) по поверхности 52 и объему V,,:
| РС{Ф} {Ф}т™ {Ге} + { р[В]т™ [А][В]ОУв {Те} =
|{Ф}д*dS2 +|Та,1ГаЛ{Ф^5з -{аЛ{Ф}{Те^,
(7)
5
S
5
где р - плотность;
C — теплоемкость; {Ф} - интерполирующая функция элемента; Ve — объем конечного элемента;
е вектор производных по времени от узловых температур элемента,
дТ
{т;} = {Ф}г = {Ф}-;
[В] = ^}{Ф}Т™;
{Те} = {Ф}Т.
Уравнение (7) в матричном виде:
[Се ] {Г.} + ([ ^ ] + [к': ]) {Ге} = {ОТ } + {£ }, (8)
где [Се ^ = |рС {Ф}{Ф}TmsdV — матрица теплоемкости КЭ;
Ке = | р[ 5] Тгш[А][В] dV — матрица теплопроводности КЭ;
[Кес ^ = | аЛ {Ф} {Te}dS3 — матрица поверхностной конвекции КЭ;
{вС | = | {Ф} q*dS2 — вектор теплового потока;
{вС} = {ТаггаЛ {Ф}dS3 - вектор поверхностной конвекции КЭ.
^
Уравнение теплопроводности для всей расчетной области строится с помощью матрицы инцидентности узлов и сетки КЭ.
Рассмотрим определение матрицы инцидентности узлов КЭ для части расчетной области из пяти треугольных КЭ, представленной на рис. 2. Общий объем расчетной области — V. Номера КЭ (с 1 по 5) помещены в кружках; номера сетки узлов, образованных конечными элементами (с 1 по 6), - в квадратиках; номера узлов в каждом КЭ (показано только в первом КЭ) обозначены г, у, к. Отмеченный заливкой угол в КЭ имеет номер г. Обход узлов в КЭ — против хода часовой стрелки. Координаты сетки узлов в общей системе координат ОХУ обозначены X, Уг (г = 1, ..., 6). На части поверхности Sl расчетной области заданы температуры Т8, на части S2 - тепловой поток q*, на части S3 -
вектор поверхностной конвекции {вС }.
Температуры Тт (т = 1, 2, ..., 6) в узлах сетки КЭ принимают за искомые переменные, связь температур Т~е) в любой точке КЭ с температурами Тт в узлах сетки определяет с помощью интерполирующих функций фп (х, у), п = г,у, к:
Т(е) = ф/еТ + ф/е)Т; + фк(е)Тк, (9)
Рис. 2. Расчетная область из пяти треугольных КЭ
где
фп(е) — интерполирующие функции, фп(е) = ^ [аПе + ь<П) x + ^ у ] ;
а«(е), Ь„(е), Сп^ - коэффициенты,
а(e) = х(е) у,(е) — х£е) у(е) • ¿(е) = у(е) — Y(e', • с(е) = Х(е) — X(е) ■
г(е) = х(е) у(е) — X(е) У(е) • Ь(е) = У(е) — У(е) • с(е) = X(е) — X(е)
I к э у 1
и ,
е) к ■
а(е) = х (е) У (е) — X(е) у(е) • ¿¿е) = У(е) — У (е) • с(е) = X(е) — X (е) •
(10)
X®, уПе) — координаты узлов КЭ в общей системе координат ОХУ (т, п = 1, 2, ..., 6); Тп — температуры в узлах КЭ, п = 7,], к. Связь температур Т~е) в точках КЭ на рис. 2 с температурами Тт в узлах сетки определяется по формулам (9):
Т(1) = ф1(1)Т1 + ф2(1)Т2 + фб(1)Тб:
Т® = ф2(2)Т2 + фз(2)Тз + фб(2)Тб Т3) = фз(3)Тз + ф4(3)Т4 + фб(3)Тб
Т4) = ф4(4)Т4 + ф4(4)Т5 + фб(4)Тб
(11)
?5) = ф1(5)Т1 + ф5(5)Т5 + фб(5)Тб.
Здесь ф„ — интерполирующие функции, вычисляемые через координаты узлов КЭ в общей (глобальной) системе координат ОХУ:
Ф1(1) = ^ [Х276 - Х6У2 + (72 - ¥6)х + X - ХМ
Ф2(1) = ^ [Х7 -Х17б + (7б - 71)х + (XI -Хб)у];
Фз(1) = ^ [Х7 - Х7 + (71 - 72)х + (Х2 - Х1)у];
............................................................... (12)
Ф1(5) = Ъ [Хб75 - Х^7б + (7б - 7з)х + (Х5 - Хб)у];
Ф5(5) =^ [Х[7б - Хб71 + (71 - 7б)х + (Хб - ХОу];
Ф6(5) = Ъ [Х571 - Х175 + (75 - 71)х + (Х1 - Х)у].
В формулах (12) 2А равна удвоенной площади треугольного КЭ на рис. 2.
Поскольку интерполирующие функции фу обладают следующим свойством: ф;- (х, , уг) = 1, если I = у ; фу (хг- , уг) = 0, если I Фу, то связь температур Т'пе) в узлах КЭ на рис. 2 с температурами Тт в узлах сетки определяется с помощью матрицы инцидентности [Г]:
{Те)} = [Г]{Т}, {Т} = [Г]Тгпк(Т:е)}. (13)
{Те)} =_[Г]_т
Т1(1) 1 0 0 0 0 0
Т2(1) 0 1 0 0 0 0
Тз(1) 0 0 0 0 0 1
Т/2) 0 1 0 0 0 0
Т2(2) 0 0 1 0 0 0 Т,
Тз(2) 0 0 0 0 0 1 Т2
Т/3) 0 0 1 0 0 0 Т3
Т2(3) — 0 0 0 1 0 0 Т4
Тз(3) 0 0 0 0 0 1 Т5
Т1(4) 0 0 0 1 0 0 Тб
Т2(4) 0 0 0 0 1 0
Т3(4) 0 0 0 0 0 1
Т/5) 1 0 0 0 0 0
Т2(5) 0 0 0 0 1 0
Т3(5) 0 0 0 0 0 1
В компьютерных программах матрицы инцидентности строят, сравнивая числовые значения глобальных координат узлов сетки с глобальными координатами узлов КЭ.
Располагая матрицей инцидентности для всей расчетной области, можно составить матричное уравнение для определения вектора температур {Т} в узлах сетки КЭ (к — число КЭ во всей расчетной области):
[С']{Т } + + [£С]){Т} = {</} + {&}, (14)
где [С] = [Г]Тгп [[ С1 ], [С ],...,[ С£ ] |[Г] - матрица теплоемкости;
[Кь] = [Г]Тпш [[ К1Ь ], [ К2Ь ],...,[ К* ШГ] — матрица теплопроводности;
[ [ К{с ], [ К2С ],...,[ К{С ]][Г] — матрица поверхностной конвекции;
[КС] = [Г]Тпш 1
{0> = [Г]Тппк{{ }, { },.,{ }}ТпШ — вектор теплового потока; Шс} = [Г]Т^5{{ в1с}, { е2с },...,{ вс }}Тппк — вектор поверхностной конвекции. Ниже, на рис. 3, приведена укрупненная блок-схема алгоритма.
Ц Актуализация программы метода конечных элементов
2] Ввод исходных данных:
'.па: гшпД/; аугД*; тахД(; ^(Т1); 1о1(Д); (Т^-ГССТС, Г,.,,„); К.„,т, ХЛТ), Ст,,„(Т), СтЛТ), Нла(Т), ИЛТ)
Л Построение геометрической модели: координаты и номера узлов сетки КЗ, номера КЭ
Назначение физических свойств КЗ
I
_5] Задание начальных условий: Тл „ в узлах сетки при г = 0; Г111Г( 0 на дневной поверхности
"б] Задание граничных условий:
Т-лг -№\ (/-номер границы); йигг./-тепловые потоки через границы; - скорость ветра; б„гГ,„ - аА(Т^г -Т) — конвекция _на дневной поверхности_
71 Начало решения уравнения теплопроводности: г = 0
1азначение начального шага по времени: Д (= штД/
I
Вычисление текущего значения времени расчета: I = / + Д/
10|Начало итерационной процедуры /= 1
П] Вычисление текущего значения температуры * в узлах:
Ш =Л{ВД-.,Ж {С(ТМ)}, \Н(ТМ)}, ДА, д*)
I
12|Провера сходимости интерационной процедуры:
({Г,} = {Гн} + Ы{Г}) П ({Н,} = ({ДЧ} + 1о1{Л}) ?
13|Следующая интерация: 1 = 1 + 1
Рис. 3 Блок-схема алгоритма: /end - конечное значение интервала времени в расчете; iend - число итераций на каждом шаге по времени; hm - толщина слоя грунта основания; ттД/, аугД/, тахД/ - минимальное, среднее, максимальное значение шага по времени; (Tb, -1 °C, 0 °C, Tend) - таблица температур изменения теплофизических параметров грунта; Tb, Tend - начальная и конечная температура; tol(T), tol(^) - критерии сходимости итерационной процедуры по температуре и энтальпии
Представленный алгоритм расчета промерзания слоистых оснований зимних лесовозных дорог успешно использовался автором в работе [2]. Алгоритм позволяет на персональных компьютерах среднего класса производить теплотехнические расчеты слоистых оснований зимних лесовозных дорог по текущим значениям температуры воздуха и скорости ветра с учетом фазовых превращений воды в лед.
СПИСОК ЛИТЕРАТУРЫ
1. Исаченко В.П., Осипова В.А., Сукомел А.С. Теплопередача. М.: Энергоиздат, 1981. 416 с.
2. Миляев А.С. Влияние тепловыделения биомассы подстилающего слоя на промерзание оснований зимних лесовозных дорог//Лесн. журн. 2010. № 1. С. 53-58. (Изв. высш. учеб. заведений).
3. Рихтмайер Р., Мортон К. Разностные методы решения краевых задач. М.: Мир, 1972. 418 с.
4. СНиП 2.06.04-82. Нагрузки и воздействия на гидротехнические сооружения (волновые, ледовые и от судов)/Госстрой России. М.: ГУП ЦПП, 1982. 37 с.
5. СНиП 2.02.04-88. Основания и фундаменты на вечномерзлых грунтах /Госстрой России. М.: ГУП ЦПП, 2000. 52 с.
Поступила 17.10.11
A.S. Milyayev
St. Petersburg State Forestry Engineering University
Algorithm of Calculation of the Frost Penetration of the Winter Logging Roads Stratified Bases
The article presents an algorithm of calculation of frost penetration depth and rate of the winter logging roads based on current air temperature and wind speed in the given area. The algorithm is based on a computational solution of a nonlinear transient problem regarding the heat transfer in an "air-ground" thermodynamic system. Phase transformations "pore water - ice" as well as dependence of ground thermalphysic characteristics on temperature to be taken into consideration.
Key words: winter logging roads, heat transfer, base frost penetration, "water-ice" phase transformations, nonlinear non-stationary equation of heat conductivity.
