РЕШЕНИЕ ЗАДАЧИ СТЕФАНА ПРИ ПРОМЕРЗАНИИ ТРУБОПРОВОДОВ
Парфентьева H.A., Самарин О.Д.
Одним из практических приложений задачи Стефана является исследование процессов, происходящих при промерзании трубопроводов водопроводных и тепловых сетей в аварийных режимах, например, при отключении сетевых насосов или выходе из строя источника теплоты в холодный период года.
Рассматриваемая задача решалась рядом авторов, в том числе в последние годы с использованием численного моделирования на современных ПЭВМ [1], [2]. Однако получаемые результаты обладают рядом недостатков из-за неполного учета факторов, влияющих на протекание процесса.
Поэтому представляется целесообразным проведение исследования, позволяющего оценить влияние дополнительных факторов на время промерзания и предложить достаточно простые соотношения, пригодные для ориентировочных расчетов. Математическая постановка данной задачи в простейшем случае может быть описана дифференциальным уравнением теплопроводности в цилиндрических координатах для слоя льда на стенках трубопровода [3]:
а д ( дг ^ _ дг
— Т ; (1)
Эх
от начала промерзания на расстоянии r, м, от
3
г дг^ дг
где 1, К - температура льда в момент времени т оси трубопровода; а = Х/ср, м2/с - коэффициент температуропроводности льда, причем р, кг/м, с Дж/(кг-К) и X, Вт/(мК) - соответственно плотность и удельная теплоемкость льда и коэффициент его теплопроводности, принимаемые по данным [4]. Тогда при 1 = 0оС а = 2.25/(2260-920) = 1.08-10"6 м2/с. Уравнение (1) решается при начальных и граничных условиях:
т = 0: t = 0, rF = rD; г = r0: t = tH; (2)
4vrF
r = rF: pr„
drF _ (^dt_
dx l dr
(условие Стефана) [3].
Здесь го, м - радиус трубопровода; гр, м - радиус фронта промерзания; гпл - удельная теплота плавления воды, равная 3.34105 Дж/К [4]; 1„, К - температура окружающей среды, в первом приближении принимаемая равной температуре грунта на уровне заложения трубопровода при подземной прокладке. В этом случае уравнение (1) решается при внешнем граничном условии 1-го рода, что упрощает исследование. Величина qv, Вт/м3, представляет собой объемную плотность потока теплоты от трения при движении воды в трубопроводе. Несложно показать, что данный параметр может быть выражен как
ч, =-
4r
-(2a),
где рж, кг/м и V, м/с - соответственно плотность воды и скорость ее движения, а \ - безразмерный коэффициент гидравлического трения. Поскольку шероховатостью слоя льда можно пренебречь, ^ вычислялся как для гидравлически гладких труб по формуле Блазиуса [4]:
Е, = 0.3164 -Ке"0-25 =0.3164 ■
4 (2b). 2rw
Здесь Ке = - безразмерное число Рейнольдса, V - кинематическая вязкость воды. При тем-
пературе около 0оС ее значение составляет порядка 1.78-10"6 м2/с [4].
Для получения наиболее достоверных результатов уравнение (1) с условиями (2) решалось численными методами с применением ПЭВМ. При этом необходимо использовать конечно-разностную аппроксимацию дифференциального уравнения теплопроводности и условия Стефана. Поскольку решается задача с подвижной границей, шаг по времени должен быть очень малым, а в этом случае для аппроксимации можно использовать явную схему [4]. Помимо существенного упрощения вычислений, такая схема позволяет за счет надлежащего подбора соотношения шагов по времени Ат и по простран-
ственной координате Ь повысить порядок аппроксимации до уровня Лт2 + Ь4. Как показано в [3], для этого достаточно выполнения равенства аЛт/Ь2 = 1/6.
Тогда температура 1-го слоя льда в j+1-й момент времени определится через температуры 1-го, 1-1-го и 1+1-го слоев для j-го момента по выражению [3], [4]:
1 6
Т =■
11, .+1
-1.1
+ 4ТГ. + Т+1. 2 '.■'
. (3)
2/ - 2 '^ 2/ - 2_ Нумерация слоев принята в данном случае от оси трубопровода.
Перемещение фронта промерзания при этом можно вычислить по конечно-разностной аппроксимации условия Стефана:
дуг Ах
АгР =-
( ХАг Л ЙРГ,„
/+1..
2р>;,
(3а)
Следует, однако, заметить, что величина qv является переменной, поэтому в качестве характерного значения при теоретическом анализе следует принимать параметр qmax, вычисленный при г = г0 и начальной величине V, а в расчетах на ПЭВМ необходим пересчет qv на каждом временном шаге.
Приближенную аналитическую зависимость для величины гР (без учета тепловыделений от трения) можно получить путем решения системы (1)-(2), если пренебречь влиянием теплоемкости ледяного слоя и принять распределение температуры в нем по выражению, соответствующему мгновенно-стационарному режиму. Для этого в (1) необходимо положить дШт = 0, откуда нетрудно найти соотношение для линейной плотности теплового потока ql, Вт/м, через поверхность фронта промерзания [4]:
<7/ = -—Г.Т- (4)
1п (гР )'
Здесь г р = гр/г0 - безразмерный радиус фронта промерзания. Из условия Стефана (2), имея в виду, что удельная поверхность границы раздела фаз на 1 пог.м трубопровода равна 2ягр, и если пренебречь значением qv, та же самая величина ql запишется так:
<7/ = /;,/>
2 ' йгр г г -—
■щ'о'г ,
а т
(4а).
= 2прг
1 Г }П
ах
Система из двух уравнений (4) и (4а) имеет два неизвестных: ql и г поэтому она замкнута и имеет единственное решение. Отсюда дифференциальное уравнение для основного параметра г р после разделения переменных приобретает следующий вид :
-гр 1п [гЕ = йГо (4Ь).
Параметр Ро'' представляет собой модифицированное число Фурье (безразмерное время), определяемое в данном случае по выражению (4с):
(4с)
Го =-
рГ
г2
т о
Интегрируя (4Ь) в пределах соответственно от 1 (поскольку при т = 0 гР = го) до г Р и от 0 до т, окончательно находим:
2Го =1 + Г 2 Г
1п (ГГ)" 2
(4ё)
Подобная зависимость была получена авторами в работе [5] для граничных условий 3-го рода. Совершая в ней предельный переход по условию стремления к бесконечности коэффициента теплообмена между трубопроводом и окружающей средой, также приходим к формуле (4ф.
Если теперь положить гР = 0, для времени полного промерзания получаем Ро'' = Учитывая влияние теплоемкости ледяного слоя, искажающей температурное поле промерзшей зоны и замедляющей фазовые превращения, можно найти зависимость (5):
Го = (1+V4; (5)
где безразмерный параметр К, определяющий соотношение теплоаккумуляции за счет теплоемкости промерзшего слоя и за счет фазовых переходов, определяется по выражению (5а):
К = с^Л/гт . (5а)
Очевидно, что такая же поправка (1 + К/2) может быть введена и в правую часть зависимости (4) для текущего радиуса фронта промерзания.
В качестве примера можно привести результаты расчетов промерзания трубопровода радиусом г0 = 0.1 м при температуре окружающей среды = -20оС и скорости движения воды w = 6.5 м/с. Графики соответствующих зависимостей приведены на Рис.1. Таким образом, неучет теплоемкости приводит к занижению времени промерзания на 12.4%, что является уже достаточно заметной величиной.
0,70
0,60
0,50
0.40
сч
0,30
0,20
0,10
0,00
4 --■•в! 3
2
1
0.00
0,10
0,20
0.Х
0,40
0.50
0,60
0,70
0,80
0,90
1,С
Рис.1. Связь гГ и Ро" для трубопровода: по формуле (4) - кривая 1; то же с поправкой на теплоемкость (5) - кривая 2; численный расчет - кривая 3; то же с учетом внутреннего трения - кривая 4.
Для учета влияния внутренних тепловыделений от трения была проведена серия расчетов на ПЭВМ по схеме (3)-(3а) при w от 0.5 до 12.5 м/с. Полученные значения продолжительности промерзания сопоставлялись с решением для w = 0, и рассчитывалось значение поправки Т к значению Fo , вычисляемому по формуле (5). При этом за определяющий параметр был принят безразмерный критерий А, представляющий собой отношение плотностей теплового потока на границе раздела фаз, создаваемого соответственно тепловыделениями от трения и за счет теплопроводности в промерзшем слое (см. формулу (6)).
Аналитически определить величину Т можно, если подставить в условие Стефана (2) зависимость для qv по выражениям ^^(2^, приняв связь между текущей скоростью воды и значением гр, исходя из постоянства перепада давлений на концах трубопровода. Тогда в режиме гидравлической гладкости можно получить, что w ~ г/'7, и окончательное выражение для Т с учетом (4)-(4d) будет иметь следующий вид:
- 1П (тр У'рйгР
т=41-
1 + А ■ г
1п (гР)
- А =
9 г2 ± max о
2Ц,
0.3164 р
У2Ч07Ч75
2'' */./
(6)
При выводе (6) учтено, что А, в отличие от 1„, неотрицательно. Интеграл в формуле (6) в элементарных функциях не выражается, но его легко можно найти численными методами. Его аппроксимация при А < 4 может быть представлена формулой (6а):
Г= 1 + А + _А1; 9 12 1200
(6а).
На Рис.2 приведены зависимости текущих значений Fo'' от с учетом влияния трения, которые можно получить из (6) заменой верхнего предела интегрирования с 1 на и исключением множителя, равного 4.
Существует некоторое предельное значение w, при котором промерзания вообще не происходит. Ее можно определить, если в условии Стефана (2) приравнять производную ^рМт нулю при = 1. При этом плотность теплового потока через поверхность трубопровода, заложенного в грунт, вычисляется по формуле Форхгеймера [6]. Тогда для критической скорости получаем:
0.58Лг//н
Чп
(7).
0
ш =
2.75
Г
о
1,0
0,8
0.7
0,5
о.з
0,2
У
/
/ /
1 / / /
> 1 в
J
# у Г У>
Л--;-;. _ - - ---- -----
— —'—
-А а 0 -А = 2 А = 4 -А = 6 -А = 7 ■А = 7.378
0,00
0.10
0,20
0,30
0,40
0,50
Г-ТГ
0,60
0,70
0.80
0,90
1.00
Рис.2. Связь между критерием Ро" и величиной 1 - гТ при различных значениях параметра А.
Здесь - коэффициент теплопроводности грунта, Вт/(мК); Ь - глубина заложения трубопровода до его оси, м. Формула (7) справедлива при Ь/г0 > 2. Принимая для влажного грунта = 0.658 Вт/(мК) [4], = -20оС как для предыдущего примера, Ь = 0.7 ми г0 = 0.25 м, находим = 2.50 м/с. Но эта скорость существенно выше характерной для городских тепловых сетей (1 - 1.5 м/с), хотя и остается еще допустимой [6]. Однако защита от замерзания теплообменных аппаратов в системах вентиляции осуществляется именно созданием скорости теплоносителя выше критической (здесь это технически возможно) [7].
При воздушной прокладке условие незамерзания становится более жестким, поскольку здесь остается только термическое сопротивление тонкого слоя льда на стенках трубопровода. Тогда из выражения (6) имеем А > 19е/7 = 7.378 (см. также Рис.2), или в нашем примере = 11.24 м/с.
Тем не менее, полученные результаты дают возможность располагать заметным резервом времени при ремонте аварийных участков тепловой сети в холодный период года. Соответствующие зависимости имеют достаточно простой вид и удобны в инженерной практике.
БИБЛИОГРФИЧЕСКИЙ СПИСОК
[1] - Рымаров А.Г. Стационарная и нестационарная теплопередача теплоизолированного трубопровода. (Сб. докл. конф. НИИСФ, 2001, с. 125 - 130).
[2] - Миханькова Ю.О. Моделирование и идентификация тепловых режимов трубопроводов систем теплоснабжения. Автореферат дисс. ... к.т.н. - Челябинск, 2002, 19 с.
[3] - Беляев Н.М., Рядно А.А. Методы теории теплопроводности. Кн. 2. - М., Высшая школа, 1982, 304 с.
[4] - Теория тепломассообмена. / Под ред. А.И.Леонтьева. - М., МВТУ, 1997, 684 с.
[5] - Парфентьева Н.А., Самарин О.Д. К расчету регенеративных теплоутилизаторов с фазовыми превращениями. // Энергосбережение и водоподготовка., 1999, №4, с. 43 - 46.
[6] - Ионин А.А. и др. Теплоснабжение. - М., Стройиздат, 1982, 336 с.
[7] - Автоматика и автоматизация систем теплогазоснабжения и вентиляции. / под ред. В.Н.Богословского. - М., Стройиздат, 1986, 479 с.