Туголуков Е.Н. - Матем. моделирование технологического оборудования [2004]

	∞
	S (x, p)=∫	(ti (x, τi )−t0 )exp(− p τi )dτi =
	0			(9.20)
	τ		tk (x, τk )−t0
	= ∫k	(ti (x, τi )−t0 )exp(− p τi )dτi +		exp(− p τk ),
			p
	0

где ti(x, τi), i = 1, 2, …, k – экспериментальные значения температур, измеренных в фиксированной точке в различные моменты времени.

Можно показать, что расчетный результат не зависит от численного значения параметра интегрального преобразования Лапласа р, которое может быть как действительным, так и комплексным.

Однако практически метод может быть использован только для регулярного температурного режима из-за его ярко выраженной неустойчивости по отношению к погрешности измерений температур при нерегулярных температурных режимах.

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

Проиллюстрируем этот метод на том же примере для пластины, температурное поле которой описывается решением задачи (9.1) – (9.4).

Введем другую замену переменной

	T (x, τ)= t(x, τ)−tс ,				(9.21)
позволяющую перейти к задаче с однородными граничными условиями:
	∂T (x, τ)	= a	∂2T (x, τ)	;	(9.22)
	∂τ		∂x2
	T (x, 0)= Т0		= t0 −tc ;		(9.23)
λ	∂T (R, τ)+αT (R, τ)= 0 ;				(9.24)
	∂x
	∂T (0, τ)		= 0 .		(9.25)
	∂x
Применим к этой задаче конечное интегральное преобразование
	R
W (τ)= ∫T		(x, τ)ρV (x)dx ,			(9.26)
	0

где функция V(х), являющаяся ядром интегрального преобразования, определяется как решение вспомогательной задачи Штурма – Лиувилля:

	d 2 S (x)				+µ2 S (x)= 0 ;			(9.27)
		dx2
λ		dS (R)			+αS (R)= 0 ;			(9.28)
		dx
			dS (0)				= 0 .	(9.29)
			dx

Решение этой задачи с точностью до постоянного множителя имеет вид

S (x)= sin (µx+ϕ).

(9.30)

Тогда

dS (x)

= µcos(µx +ϕ).

(9.31)

dx

µсos(ϕ)= 0 ,

Из граничного условия (9.29)

(9.32)

откуда (минимальный положительный корень) ϕ =

π

,

(9.33)

2

следовательно

S (x)= сos(µx ).

(9.34)

Из граничного условия (9.28) имеем

−λµsin (µR )+αcos(µR )= 0 ,

(9.35)

откуда в конечном итоге и определяется искомый коэффициент теплоотдачи

α = λµtg (µR).

(9.36)

Весовая функция ρ определяется, как решение уравнения

′

(9.37)

a ρ = 0,

откуда ρ = сonst, в частности, ρ = 1.

(9.38)

Выполняем переход к изображениям:

dW (τ)

+a µ2 W (τ)

= 0 ;

(9.39)

dτ

изображение начального условия (9.23):

W (0)=

R T cos (µx)dx =

T0

sin (µR).

(9.40)

µ

∫

0

0

Решение задачи (9.22) – (9.23) в изображениях:

W (τ)= W (0)exp (−µ2 a τ)=

T0

sin (µR)exp (−µ2 a τ).

(9.41)

µ

Изображение температурного профиля находим по формуле:

R

W (τ)= ∫Ti

cos(µx)dx ,

(9.42)

0

где Ti = T (xi ,τ) – массив экспериментальных значений температур;

i = 1, 2, …, N – номер точки изме-

рения температуры.

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

– подынтегральная функция; xi – координата i-й точки изме-

При наличии не менее четырех точек измерения температуры по толщине пластины с равномерными интервалами удобна схема Эйлера третьего порядка точности [2]

R	N
I = ∫Ti cos(µx)dx ≈ h(G2 −7	(G1 +GN )+GN −1)/12+h∑Gi , (9.43)
0	i =1

где h – величина интервала; Gi = Ti сos (µxi ) рения температуры.

Значение µ определяется как любой положительный корень уравнения

	I =	T0	sin (µR)exp (−µ2 a τ).	(9.44)
		µ

Лучше выбирать первый корень во избежание накопления погрешности компьютерного счета. Теперь для момента времени τ из уравнения (9.36) можно найти искомое значение коэффициента

теплоотдачи α.

Устойчивость метода обусловлена сглаживанием значений экспериментальных температур при интегрировании по толщине пластины.

10 МЕТОДИКА РАСЧЕТА ТЕПЛООБМЕННОГО ОБОРУДОВАНИЯ

Рассмотрим элементарную область одноходового кожухотрубчатого теплообменника, работающего в стационарном температурном режиме.

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

Принятое допущение позволяет выделить элементарную область длиной ∆х вдоль оси аппарата, охватывающую все его поперечное сечение, включая трубный пучок, корпус и теплоизоляционное покры-

тие (рис. 10.1).

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

∆x

РИС. 10.1 ЭЛЕМЕНТАРНАЯ ОБЛАСТЬ КОЖУХОТРУБЧАТОГО ТЕПЛООБМЕННИКА

Примем допущение о постоянстве теплофизических характеристик теплоносителей внутри элементарной области. Значения теплофизических характеристик определяются сначала температурами теплоносителей на входе в элементарную область, а затем, в ходе расчета, средними температурами теплоносителей в элементарной области.

Температуры стенок трубок и корпуса принимаются постоянными по длине элементарной области. По толщине их температуры меняются.

Пусть горячий теплоноситель находится в трубном пространстве, и его характеристики имеют индекс 1.

Пусть G1, G2 – массовые расходы теплоносителей; t1, t2 – температуры теплоносителей на входе в элементарную область; α1, α2, αк, αoc – коэффициенты теплоотдачи соответственно от внутренней и внешней поверхностей трубки к теплоносителям, от внутренней поверхности корпуса, от наружной поверхности теплоизоляции; с1, с2, ρ1, ρ2, λ1, λ2 – соответственно теплоемкости, плотности и теплопроводности теплоносителей; ст, ск, си, ρт, ρк, ρи, λт, λк, λи – соответственно теплоемкости, плотности и теплопроводности материала трубок, корпуса и теплоизоляционного покрытия; dк, dт – соответственно внутренние диаметры кожуха аппарата и трубки трубного пучка; δт, δк, δи – соответственно толщины стенок трубок, стенки корпуса и теплоизоляционного покрытия; tос – соответственно температура окружающей среды.

При известных температурных полях стенки трубки, корпуса и слоя теплоизоляции расчет тепловых потоков и температур теплоносителей внутри элементарной области и на выходе из нее не представляет трудностей, поэтому рассмотрим возможности расчета этих полей.

Структура элементарной области кожухотрубчатого теплообменника, работающего в стационарном температурном режиме, представлена на рис. 10.2.

Здесь стенка трубки трубного пучка рассматривается как полый неограниченный цилиндр, в частном случае однослойный.

G1	t1			t1
		t1с		α1	rт
	δ1	λк	U(r)
G2	t2	t2с			α2
				t2
					rк
		tк		αк
	δк	λк		S1(r)
	δи λи			tгр S2(r)
		tос		tи
				αос

∆x

Рис. 10.2 Структура элементарной области кожухотрубчатого теплообменника, работающего в стационарном температурном режиме

Теплоизолированная стенка корпуса аппарата также рассматривается как полый неограниченный цилиндр, в частном случае двухслойный.

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

Это обосновано следующими соображениями:

• продольные тепловые потоки в тонких стенках и стенках с низкой теплопроводностью при небольших перепадах температур пренебрежимо малы; так, в кожухотрубчатом теплообменнике тепловая мощность продольного теплового потока в стенках трубок на 2 – 3 порядка меньше тепловой мощности, передаваемой от горячего к холодному теплоносителю;

•для элементарной области входящие и выходящие продольные тепловые потоки близки, в результате чего в элементарном тепловом балансе они практически взаимно компенсируются.

Таким образом, температурное поле элементарной области описывается следующими функциями:

•t1(х) – температурное поле потока в трубном пространстве;

•U(r) – температурное поле стенки трубки трубного пучка;

•t2(х) – температурное поле потока в межтрубном пространстве;

•S1(r1) – температурное поле стенки корпуса аппарата;

•S2(r2) – температурное поле слоя теплоизоляционного покрытия.

Эти функции являются решениями соответствующих задач теплопроводности, совокупность ко-

торых составляет математическую модель температурного поля элементарной области.

dt1 (x)

+K1t1 (x)= V1 (x),

0 ≤x ≤∆ x ;

(10.1)

dx

t1 (0)= t10 ;

(10.2)

d 2U (r)

+

1

dU(r)

= 0,

r ≤r ≤r +δ

т

;

(10.3)

dr2

r

dr

т

т

λ

dU (rт )

+α (U (r )−t ) = 0 ;

(10.4)

т

dr

1

т

1

λт

dU (rт +δт )

+α2 (U (rт +δт )−t2 )= 0 ;

(10.5)

dr

dt2 (x)

+K2 t2 (x)= V2 (x),

0 ≤x ≤∆ x ;

(10.6)

dx

t2 (0)= t20 ;

(10.7)

d

2S

(r )

1 dS

i

(r )

Ri −1 ≤ ri ≤ Ri ;

(10.8)

i

i

+

i

= 0,

i =1, 2,

dr

2

r

dr

i

i

i

λ

dS1(R0 )

+α

(S (R )−t

)

= 0 ;

(10.9)

к

dr

к

1

0

2

1

λ

и

dS2 (R2 )

+α

oc

(S

2

(R

2

)−t

oc

)= 0 ;

(10.10)

dr

2

S (R )= S

(R

);

λ

dS1(R1 )

= λ

dS2 (R1 )

,

(10.11)

1

1

2

2

к

dr

и

dr

1

2

где

K

1

= α1П1

,

V

(x)= K t

(x);

(10.12)

G1c1

1

1 F

1

K2 =

α

к

П

к

+α

2

П

2 ,

V2 (x)

=

αкПкtF

(x)+α

2

П2tF

(x)

;

(10.13)

к

2

G2c2

G2c2

= 2π(rт +δт )n,

П1 = 2πrтn,

П2

Пк

= 2πrк ;

(10.14)

R0

= rк,

R1 = rк +δк,

R2

= rк +δк +δи.

(10.15)

Если температуры поверхностей стенок трубок и корпуса в

элементарной области посто-

янны, т.е. tF

(x)= сonst = tF ,

tF

(x)= сonst = tF

и

tF (x)= сonst = tF

, то искомые функции с учетом граничных

1

1

2

2

к

к

условий имеют вид:

t1(x)= tF +(t10 −tF )exp (−K1x);

(10.16)

1

1

U (r)= AU +BU ln (r);

(10.17)

t2 (x)= V2 +(t20 −V2 )exp (−K2 x);

(10.18)

Si (ri )= Ai

+Bi

ln (ri ),

i = 1, 2,

(10.19)

где

V2

=

αкПкtF +α2П2tF

;

(10.20)

к

2

αкПк +α2П2

BU =

t2 −t1

;

(10.21)

δ

1

1

т

ln

1+ r

+λт α

2

(r +δ

т

)−

α r

т

т

1 т

λ

т

;

(10.22)

AU = t1 −BU ln (rт )+

α r

1 т

B1

=

t2 −toc

; (10.23)

λк

λк

λи

1

1

ln(R0 )+ R α

к

−λ

и

ln(R2 )+ R α

oc

+λк ln(R1)

λ

к

−λ

0

2

и

			λ		к			;
A1	= t2
		−B1 ln(R0 )+	α	к	R
						0

B2 = λк B1 ;

λи

λ

A2

= A1

к

+B1 ln(R1 ) 1−

λ

.

и

Средняя температура жидкости в элементарной области:

(t10 −tF

)

1

∆x

(x)dx = tF1 +

(1−exp

(−K1∆x));

t1

=

∫t1

1

∆x

K

∆x

0

1

1

∆x

(t20 −V2 )

t

=

t2 (x)dx = V2 +

(1−exp (−K2 x)).

2

∆x ∫

K

2

∆x

0

Сопряжение температурных полей происходит на поверхностях трубок и корпуса:

tF = U

(rт );

tF = U (rт +δт ); tF = S1 (rк ).

1

2

к

(10.24)

(10.25)

(10.26)

(10.27)

(10.28)

(10.29)

Численное значение коэффициентов теплоотдачи α1 и αк отрицательно для учета направления тепловых потоков и сохранения общности подхода к решению задач теплопроводности, принятому в данной работе.

В общем случае:

α1 = α1 (t1,U (rт )), α2 = α2 (t2 ,U (rт +δт ));

αк = αк (t2 ,V1 (rк )), αoc = αoc (toc ,V2 (rк +δк +δи )).

Явный вид этих зависимостей определяется соответствующими критериальными уравнениями. При наличии в межтрубном пространстве перегородок или экранов, вводятся соответствующие по-

правки.

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

1 Используя для расчета коэффициентов теплоотдачи температуры поверхностей стенок трубки и корпуса из предыдущей элементарной области, находим средние температуры теплоносителей в текущей элементарной области по формулам (10.27), (10.28).

2По полученным значениям средних температур и коэффициентов теплоотдачи рассчитываем температуры поверхностей стенок трубки и корпуса по формулам (10.17) и (10.19).

3Для найденных температур поверхностей стенок трубки и корпуса уточняем коэффициенты теплоотдачи и средние температуры теплоносителей в текущей элементарной области.

4Повторяем пункты 2 и 3.

5Вычисляем конечные температуры теплоносителей по формулам (10.16), (10.18).

Тепловой баланс элементарной области может быть использован для независимой проверки и оценки качества расчета (как и тепловой баланс по всему аппарату в целом).

∆Q1 = ∆Q2 + ∆Qп,

(10.30)

где ∆Q1 – тепловая мощность, отдаваемая горячим теплоносителем; ∆Q2 – тепловая мощность, получаемая холодным теплоносителем; ∆Qп – тепловая мощность потерь в окружающее пространство.

Здесь	∆Q1 = G1c1 (t1k −t10 );	(10.31)
		∆Q21 = G2c2 (t2k −t20 );	(10.32)
		∆Qп = αoc 2πR2∆x (V2 (R2 )−toc ).	(10.33)

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

∆G =	∆Q ,	(10.34)
	rf

где rf – удельная теплота фазового перехода.

Выполняется пересчет теплофизических характеристик теплоносителя, частично совершившего фазовый переход.

Проверяется, весь ли теплоноситель совершил фазовый переход.

После того, как весь теплоноситель совершит фазовый переход, продолжается дальнейший пересчет его температур.

Таким образом, расчет элементарной области включает следующие действия:

•расчет текущих значений теплофизических характеристик теплоносителей в зависимости от их текущих температур и агрегатного состояния;

•расчет скоростей теплоносителей и коэффициентов теплоотдачи;

•расчет средних значений температур теплоносителей в элементарной области;

•итеративный расчет поверхностных температур стенки, разделяющей теплоносители;

•итеративный расчет поверхностных температур теплоизолированного корпуса аппарата;

•расчет выходных температур теплоносителей и тепловых потерь для элементарной области. Результирующими данными расчета элементарной области являются:

•t1k, t2k – конечные температуры горячего и холодного теплоносителей;

•Gф1, Gф2 – массовые расходы части теплоносителей, совершившей фазовый переход;

•∆Qп – тепловая мощность, теряемая в окружающую среду.

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

1Теплообмен при отсутствии фазовых переходов в теплоносителях.

2Теплообмен при конденсации одного из теплоносителей.

3Теплообмен при испарении одного из теплоносителей.

4Теплообмен при одновременных фазовых переходах в обоих теплоносителях.

В этом случае расчет представляет собой последовательное использование перечисленных процедур в заданной последовательности.

Этим исключается необходимость задания избыточных исходных данных для расчета частных случаев теплообменных процессов.

Теперь рассмотрим элементарную область кожухотрубчатого теплообменника при нестационарном температурном режиме работы.

Пусть горячий теплоноситель также находится в трубном пространстве и его характеристики имеют индекс 1. Также сохраняются все остальные обозначения, принятые при моделировании стационарного режима.

При нестационарном температурном режиме расходы теплоносителей и их начальные температуры могут оставаться постоянными или меняться во времени.

В этом случае элементарная область представляет собой элементарный внутренний объем аппарата, определяемый так же, как в стационарном случае, но рассматриваемый в течение элементарного интервала времени ∆τ.

Допущение о постоянстве теплофизических характеристик теплоносителей внутри элементарной области, принятое для элементарной области при стационарном режиме, остается в силе.

Тепловой баланс элементарной области в нестационарном температурном режиме:

∆Q1 = ∆Q2 + ∆Qт + ∆Qк + ∆Qи + ∆Qп ,

(10.35)

где ∆Q1 – тепловая мощность, отдаваемая горячим теплоносителем; ∆Q2 – тепловая мощность, получаемая холодным теплоносителем; ∆Qт – тепловая мощность, затрачиваемая на нагрев стенок трубок; ∆Qк – тепловая мощность, затрачиваемая на нагрев стенки корпуса; ∆Qи – тепловая мощность, затрачиваемая на нагрев слоя теплоизоляции; ∆Qп – тепловая мощность потерь в окружающее пространство.

Температурное поле элементарной области при нестационарном температурном режиме также моделируется решением сопряженной задачи теплопроводности, включающей задачи нестационарной теплопроводности для полого многослойного неограниченного цилиндра с произвольным начальным распределением и несимметричными неоднородными граничными условиями 3-го рода на внутренней и внешней поверхностях (температурные поля стенок трубок и теплоизолированного корпуса), и уравнениями движения жидкости, движущейся в режиме идеального вытеснения по каналу (температурные поля теплоносителей в трубном и межтрубном пространствах).

Необходимость использования постановки задач теплопроводности с произвольным начальным распределением диктуется методикой использования температурного поля элементарной области для моделирования температурного поля аппарата, при которой для очередной элементарной области прихо-

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

При отсутствии загрязнений стенка трубки рассматривается как полый однослойный неограниченный цилиндр, при наличии – соответственно как многослойный. Теплоизолированный корпус аппарата в зависимости от наличия загрязнений, покрытий, экранов, также может рассматриваться как двухили многослойный полый неограниченный цилиндр.

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

∂t1 (x, τ)

+W1

∂t1 (x, τ)

+ K1t1

(x, τ)= F1 (x, τ),

0 ≤ x ≤∆ x;

(10.36)

∂τ

∂x

t1(0, τ)= t10 (τ);

t1(x, 0)= f1(x);

(10.37)

∂U (r, τ)

= a

2

∂2U (r, τ)

+

1

∂U (r, τ)

,

r

≤r ≤r +δ

т

,

τ >0;

(10.38)

∂τ

t

∂r

2

r

∂r

т

т

U (r,0)= ft (r);

(10.39)

λт

∂U (rт, τ)

+α1 (U (rт, τ)−t1 )= 0;

α1 <0;

(10.40)

∂r

λ

т

∂U (rт +δт , τ)

+α

2

(U (r +δ

т

, τ)−t

2

)= 0;

(10.41)

∂r

т

∂t2 (x, τ)

+W2

∂t2 (x, τ)

+K2t2 (x, τ)= F2 (x, τ),

0 ≤ x ≤∆ x ;

(10.42)

∂τ

∂x

t2 (0, τ)= t20 (τ);

t2 (x,0)= f2 (x);

(10.43)

∂S

(r ,

τ)

2

∂2 S

i

(r , τ)

1

∂S

(r , τ)

i

i

i

i

i

Ri −1 ≤ri ≤Ri ;

(10.44)

= ai

2

+

,

i =1, 2,

∂τ

∂r

r

∂r

i

i

i

Si (ri

, 0)= ϑi (ri );

(10.45)

λ

dS1(R0 , τ)

+α

(S (R , τ)−t

)= 0 ;

(10.46)

к

к

2

dr1

1

0

λи

dS2 (R2 , τ)

+αoc (S2 (R2 , τ)−toc )= 0 ;

(10.47)

dr2

S

(R

, τ)

= S

2

(R

2

, τ);

λ

к

dS1(R1 , τ)

= λ

и

dS2 (R1 , τ)

,

(10.48)

1

1

dr1

dr2

где

K

1

= α1П1

,

F

(x, τ)= K t

F

(x, τ);

(10.49)

G1c1

1

1

1

K2

=

α

к

П

к

+α

2

П

2 ,

V2 (x, τ)=

αкПкtF

(x, τ)+α

2П2tF

(x, τ)

;

(10.50)

к

2

G2c2

G2c2

П1 = 2πrтn,	П2	= 2π(rт +δт )n,			Пк = 2πrт ;				(10.51)
R0 = rк ,	R1 = rк +δк,			R2 = rк +δк +δи .					(10.52)
Здесь t1, t2 – средние температуры теплоносителей в элементарной области, tF					(x, τ), tF	(x, τ),		tF	(x, τ) –
				1	2			к
температуры внутренней и наружной поверхностей трубок, внутренней поверхности корпуса.
Решения уравнений (10.36) – (10.48) имеют следующий вид.
			x				,		(10.53)
t (x, dτ)= exp(−P x) t (dτ)+			∫	θ (x, dτ)exp(P x)dx
1	1	10		1	1
			0

где

K

dτ+1

1

f

1

(x)

P =

1

;

θ (x, dτ)=

F (x, dτ)+

.

W1dτ

dτ

1

1

W1 1

U (r,τ)

∞

ξ(νn ,τ)ζ(r,νn ),

= A+Bln(r)+∑

n =1

Zn

где

B =

t1 −t2

;

(10.56)

1

1

α

− (r

+δ

)α

ln(rт )−ln(rт +δт )+λт r

1

т

2

т

т

A = t1 −B ln(rт )+ αλ1 сrт ;

ξ(νn , τ)= ξ(νn ,0)exp(−ν2n τ);

ζ(r, νn )= J0 νanсr + Dn Y0 νanсr ;

rт +δт

ξ(νn , 0)

=

∫r(ft (r)− A −B ln (r))ζ(r, νn )dr ;

rт

λ

т

ν

n

ν

n

ν

n

−α1

at

J1

at

rт

J0

at

rт

Dn

=

;

ν

n

λ

т

ν

n

ν

n

−

α1 Y0

at

rт

at

Y1

at

rт

νn – последовательные положительные корни уравнения

ν(rт +δт )

νλт

ν(rт +δт )

J

J1

0

at

−

at α2

at

+

ν(rт +δт )

νλт

ν(rт +δт )

−

Y1

= 0;

+D Y0

at

at

at α2

(10.54)

(10.55)

(10.57)

(10.58)

(10.59)

(10.60)

(10.61)

(10.62)