Решение задачи Стефана при промерзании трубопроводов

РЕШЕНИЕ ЗАДАЧИ СТЕФАНА ПРИ ПРОМЕРЗАНИИ ТРУБОПРОВОДОВ

Парфентьева H.A., Самарин О.Д.

Одним из практических приложений задачи Стефана является исследование процессов, происходящих при промерзании трубопроводов водопроводных и тепловых сетей в аварийных режимах, например, при отключении сетевых насосов или выходе из строя источника теплоты в холодный период года.

Рассматриваемая задача решалась рядом авторов, в том числе в последние годы с использованием численного моделирования на современных ПЭВМ [1], [2]. Однако получаемые результаты обладают рядом недостатков из-за неполного учета факторов, влияющих на протекание процесса.

Поэтому представляется целесообразным проведение исследования, позволяющего оценить влияние дополнительных факторов на время промерзания и предложить достаточно простые соотношения, пригодные для ориентировочных расчетов. Математическая постановка данной задачи в простейшем случае может быть описана дифференциальным уравнением теплопроводности в цилиндрических координатах для слоя льда на стенках трубопровода [3]:

а д ( дг ^ _ дг

— Т ;

Эх

от начала промерзания на расстоянии 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;

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 '.'

.

2/ - 2 '^ 2/ - 2_ Нумерация слоев принята в данном случае от оси трубопровода.

Перемещение фронта промерзания при этом можно вычислить по конечно-разностной аппроксимации условия Стефана:

дуг Ах

АгР =-

( ХАг Л ЙРГ,„

/+1..

2р>;,

Следует, однако, заметить, что величина qv является переменной, поэтому в качестве характерного значения при теоретическом анализе следует принимать параметр qmax, вычисленный при г = г0 и начальной величине V, а в расчетах на ПЭВМ необходим пересчет qv на каждом временном шаге.

Приближенную аналитическую зависимость для величины гР (без учета тепловыделений от трения) можно получить путем решения системы (1)-(2), если пренебречь влиянием теплоемкости ледяного слоя и принять распределение температуры в нем по выражению, соответствующему мгновенно-стационарному режиму. Для этого в (1) необходимо положить дШт = 0, откуда нетрудно найти соотношение для линейной плотности теплового потока ql, Вт/м, через поверхность фронта промерзания [4]:

<7/ = -—Г.Т-

1п (гР )'

Здесь г р = гр/г0 - безразмерный радиус фронта промерзания. Из условия Стефана (2), имея в виду, что удельная поверхность границы раздела фаз на 1 пог.м трубопровода равна 2ягр, и если пренебречь значением qv, та же самая величина ql запишется так:

<7/ = /;,/>

2 ' йгр г г -—

щ'о'г ,

а т

(4а).

= 2прг

1 Г }П

ах

Система из двух уравнений (4) и (4а) имеет два неизвестных: ql и г поэтому она замкнута и имеет единственное решение. Отсюда дифференциальное уравнение для основного параметра г р после разделения переменных приобретает следующий вид :

-гр 1п [гЕ = йГо (4Ь).

Параметр Ро'' представляет собой модифицированное число Фурье (безразмерное время), определяемое в данном случае по выражению (4с):

Го =-

рГ

г2

т о

Интегрируя (4Ь) в пределах соответственно от 1 (поскольку при т = 0 гР = го) до г Р и от 0 до т, окончательно находим:

2Го =1 + Г 2 Г

1п (ГГ)" 2

Подобная зависимость была получена авторами в работе [5] для граничных условий 3-го рода. Совершая в ней предельный переход по условию стремления к бесконечности коэффициента теплообмена между трубопроводом и окружающей средой, также приходим к формуле (4ф.

Если теперь положить гР = 0, для времени полного промерзания получаем Ро'' = Учитывая влияние теплоемкости ледяного слоя, искажающей температурное поле промерзшей зоны и замедляющей фазовые превращения, можно найти зависимость (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) учтено, что А, в отличие от 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 с.