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

•перепад температур в поперечном сечении стержня или по толщине пластины отсутствует;

•теплопроводность материала стержня или пластины не зависит от температуры;

•стержень или пластина имеют постоянное сечение.

Рассмотрим температурное поле стержня.

Выделим элементарную область длиной ∆х по направлению движения теплового потока. Запишем составляющие теплового баланса для элементарной области.

Q1 – тепловая мощность, подводимая теплопроводностью к элементарной области; Q2 – тепловая мощность, отводимая теплопроводностью из элементарной области.

Тогда тепловая мощность, отдаваемая теплоносителю на элементарном участке, равна

Q1 − Q2 = α f (t (х)−tt ),	(7.1)
где α – коэффициент теплоотдачи от поверхности стержня к теплоносителю; f =П∆x	– омываемая пло-

щадь поверхности элементарной области; П – периметр элементарной области; t(х) – текущая температура стержня; tt – температура теплоносителя; х – координата, направленная по длине стержня.

С другой стороны

Q1 −Q2 = F (q1 −q2 ),

(7.2)

где F – площадь поперечного сечения стержня; q1, q2 – плотности тепловых потоков, соответственно подводимых к элементарной области и отводимых от нее теплопроводностью.

Устремляя ∆х к нулю, имеем:

dq (x)

d

dt (x)

d 2t (x)

(7.3)

q1 −q2

= −

∆x = −

−λ

dx = λ

dx.

dx

dx

dx2

dx

	Тогда	λ	d 2t (x)		F dx = αП(t (х)−tt )dx.				(7.4)
			dx2
	Введя обозначения T (x)= t (x)−tt и k 2 =						αП	,	(7.5)
							λF
	окончательно получим:				d 2T (x)	−k 2T (x)= 0.			(7.6)
					dx2

Общее решение этого уравнения имеет вид:

T (x)= C1 сh (k x)+C2 sh (k x).

(7.7)

При использовании краевых условий вида:

−λ

dT (0)

= q0 ;

(7.8)

dx

−λ

dt (L)

= αt (t (L)−tt ),

(7.9)

dx

где q0 – входящий тепловой поток; L – длина стержня; αt – торцевой коэффициент теплоотдачи,

αt tt +

q0

(αt sh (kL)+λk сh (kL))

имеем:

C1

=

λk

;

(7.10)

αt сh (kL)+λk sh (kL)

C2	= −	q0	.	(7.11)
		λk

Теперь рассмотрим температурное поле пластины. Обозначим толщину пластины h.

При подводе тепловых мощностей к соседним торцам пластины получаем двумерное температурное поле.

Выделим элементарную область размерами ∆х × ∆у × h.

Запишем составляющие теплового баланса для выделенной элементарной области.

Qх1 – тепловая мощность, подводимая теплопроводностью к элементарной области в направлении координаты х;

Qх2 – тепловая мощность, отводимая теплопроводностью из элементарной области в направлении координаты х;

Qу1 – тепловая мощность, подводимая теплопроводностью к элементарной области в направлении координаты у;

Qу2 – тепловая мощность, отводимая теплопроводностью из элементарной области в направлении координаты у.

Тогда тепловая мощность, отдаваемая теплоносителю на элементарном участке, равна

Qx1 −Qx2 +Qy1 −Qy2 = (α1 +α2 )f (t (х, y)−tt ),

(7.12)

где α1 и α2 – коэффициенты теплоотдачи от наружных поверхностей пластины к теплоносителю; f = ∆x ∆y – омываемая площадь поверхности элементарной области на каждой из плоскостей пластины;

t(х, у) – текущая температура пластины; tt – температура теплоносителя; х, у – координаты. С другой стороны,

Qx1 −Qx2 = Fx (qx1 −qx2 ), Qy1 −Qy2 = Fy (qy1 −qy2 ), (7.13)

где Fx =h ∆y , Fy =h ∆x – площади поперечного сечения элементарной области в направлениях х и у соот-

ветственно; qх1, qу1, qх2, qу2 – плотности тепловых потоков, соответственно подводимых к элементарной области и отводимых от нее теплопроводностью в соответствующих направлениях.

Устремляя одновременно ∆х и ∆у к нулю, имеем:

∂q (x, y)

∂

∂t (x, y)

∂2t (x, y)

qx1 −qx2

=−

∆x =−

−λ

dx =λ

dx;

∂x

∂x

∂x2

∂x

(7.14)

∂q (x, y)

∂

∂t (x, y)

∂

2

t (x, y)

qy1 −qy2

=−

∆y

=−

−λ

dy.

2

∂y

∂y

dy =λ

∂y

∂y

Тогда

∂2t (x, y)

∂2t (x, y)

λ

+

h dx dy =(α

1

+α

2

)(t (х, y)−t

)dx dy. (7.15)

2

2

∂x

∂y

t

Введя обозначения

T (x, y)= t (x, y)−tt и

k 2 =

(α1 +α2 )

,

(7.16)

λh

окончательно получим:

∂2T (x, y)

+

∂2T (x, y)

−k 2T (x, y)= 0.

(7.17)

∂x2

∂y2

При граничных условиях вида:

		−λ	∂T (0, y) = qx0 ;			(7.18)
			∂x
		−λ	∂T (x ,0) = qy0 ;			(7.19)
			∂y
	−λ	∂T (Lx , y)			= αxt T (Lx , y);	(7.20)
		∂x
	−λ	∂T (x, Ly )			= αyt T (x , Ly ),	(7.21)
		∂y

где qх0, qy0 – входящие тепловые потоки по направлениям х и у; Lх, Ly – длины пластины в направлениях х и у; αхt, αуt – торцевые коэффициенты теплоотдачи, получаем задачу, которая решается методом конечных интегральных преобразований.

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

d 2 P(y)+µ2 P(y)= 0 dy2

с однородными граничными условиями:

dPdy(0) = 0 ;

−λ dPdy(Ly )= αyt P (Ly ).

С точностью до постоянного множителя

P (y)= cos (µ y),

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

µλsin (µLy )=αyt сos (µLy ).

Переход к изображениям выполняется по формуле

Ly

U (x)= ∫T (x, y)ρ(y)P (y)dy ,

0

где ρ(y)= 1.

Обратный переход выполняется по стандартной формуле

T (x, y)= ∑∞ U (x)P (y),

n=1 N

где суммирование ведется по значениям µn и

	Ly	Ly		1	(µLy +sin (µLy )cos (µLy )).
	N = ∫ρ(y)P2	(y)dy = ∫cos2	(µ y)dy =
				2µ
	0	0

(7.22)

(7.23)

(7.24)

(7.25)

(7.26)

(7.27)

(7.28)

Переходим к изображениям задачи (7.17) – (7.21):

d 2U (x)		−ν2U (x)+	qy0	= 0,	ν2 = µ2 +k 2 ;
dx2			λ

−λ dUdx(0) = Qx0 ;

							−λ	dU (Lx )		= αxt U (Lx ),
								dx
		Ly			µ
		∫
где	Qx0 =		qx0	P (y)dy =	qx0	sin (µLy ).	(7.33)
		0

Решение задачи имеет вид:

U (x)= Ash (νx)+B сh (νx)+ λνqy02 .

А и В определяются из краевых условий:

	A = −	Qx0	;
		νλ
	A(λνch (νLx )+αxt sh (νLx ))+			αxt qy0
B = −					λν2		.
	αxt ch (νLx )+λνsh (νLx )

Таким образом, все компоненты решения (7.28) определены.

8 ДИФФЕРЕНЦИАЛЬНОЕ УРАВНЕНИЕ ПЕРЕНОСА ТЕПЛА ЖИДКОСТЬЮ, ДВИЖУЩЕЙСЯ В РЕЖИМЕ ИДЕАЛЬНОГО ВЫТЕСНЕНИЯ ПО КАНАЛУ

(7.29)

(7.30)

(7.31)

(7.32)

(7.34)

(7.35)

(7.36)

Рассмотрим случай, когда движущаяся жидкость омывает стенки канала с двумя разными температурами. Эта ситуация соответствует случаям движения теплоносителя по кольцевому каналу рубашки емкостного аппарата или по межтрубному пространству кожухотрубчатого аппарата. Случай, когда движущаяся жидкость омывает стенки канала с одной температурой, рассмотрен в [7].

При выводе уравнения используются следующие допущения:

•температура жидкости по сечению канала не меняется;

•жидкость несжимаемая;

•теплофизические характеристики жидкости не зависят от температуры;

•канал имеет постоянное сечение.

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

Рассмотрим нестационарный температурный режим.

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

Примем следующие обозначения:

х – пространственная координата по направлению движения потока; τ – время;

t(x, τ) – текущая температура жидкости; G – массовый расход жидкости;

с – теплоемкость жидкости;

qi – плотность теплового потока через стенку канала;

∆Fi = Пi∆х – площадь i-й стенки канала элементарной области; Пi – омываемый периметр i-й стенки канала;

αi – коэффициент конвективной теплоотдачи от i-й стенки канала к жидкости; tFi(x, τ) – температура i-й стенки канала;

i = 1, 2 – индекс поверхности канала.

Выделим элементарную область длиной ∆х по направлению движения потока. Запишем составляющие теплового баланса для элементарной области.

Тепло, привносимое потоком жидкости за время ∆τ

Q0 = G ct(x,τ)∆τ ;		(8.1)
тепло, уносимое потоком жидкости за время ∆τ
Q3 = −G c t(x + ∆х, τ)∆τ;		(8.2)
тепло, отдаваемое элементарной области первой поверхностью
Q1 = q1∆F1	= α1 (tF1 (x, τ)−t (x, τ))∆F1 ∆τ ;	(8.3)
тепло, отдаваемое элементарной области второй поверхностью
Q2 = q2 ∆F2	= α2 (tF (x, τ)−t (x, τ))∆F2 ∆τ .	(8.4)
	2

Составляющие элементарного теплового баланса можно записать через приращения независимых переменных:

Q1	= α1	(tF (x, τ)−t (x, τ))П1∆х∆τ;	(8.5)
		1
Q2	= α2 (tF2 (x, τ)−t(x, τ))П2 ∆х∆τ.		(8.6)

Изменение теплосодержания жидкости в элементарной области за время ∆τ:

Q = S ∆xρc(t (x,τ+∆τ)−t (x,τ)),

(8.7)

где S – площадь сечения канала; ρ – плотность жидкости.

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

или	Q = Q0 +Q1 +Q2 +Q3	(8.8)
	S ∆xρc (t (x, τ+∆τ)−t (x, τ))=
= G ct(x, τ)∆τ+α1(tF1 (x, τ)−t (x, τ))П1∆х∆τ+		(8.9)
+α2	(tF2 (x, τ)−t (x, τ))П2∆х∆τ−G ct(x + ∆х, τ)∆τ.

Разделив уравнение почленно на произведение ∆x∆τ и выполнив предельные переходы, получим:

Sρc	∂t (x, τ)	= −Gc	∂t (x, τ)	+			(8.10)
	∂τ		∂x			(tF
		+α1 (tF (x, τ)−t(x, τ))П1			+α2		(x, τ)−t(x, τ))П2 .
			1			2

Далее, поделив уравнение почленно на произведение Sρc и перегруппировав слагаемые, получим

окончательный вид уравнения:
			∂t (x, τ)	+W		∂t (x, τ)	+ α1П1 +α2П2	t	(x, τ)=
			∂τ			∂x	Sρc				(8.11)
			=		α1П1tF (x, τ)+α2П2tF (x, τ)				,
						1	2
							Sρc
где	W =	G	– скорость движения жидкости.
		Sρ
	Запишем задачу в виде:
			∂t (x, τ)+W ∂t (x, τ)+K t (x, τ)= F(x, τ);								(8.12)
			∂τ			∂ x
где			t(0, τ)= t0 (τ);				t(x, 0)= f (x),				(8.13)
						K = α1П1 +α2П2 ;
							Sρc				(8.14)
			F(x, τ)=			α1П1tF	(x, τ)+α2П2tF (x, τ)
						1	2			.
							Sρc

Полное решение уравнения (8.12) приводится в [1].

При использовании решения для моделирования температурных полей элементарных областей возможно применение более простого решения задачи (8.12) – (8.13), полученного при использовании ко-

нечно-разностного аналога частной производной ∂t (∂xτ, τ):

∂t (x, τ)	≈	t (x, τ)−t (x, τ−dτ)	.	(8.15)
∂τ
		dτ

При этом для фиксированного значения времени dτ внутри каждого временного интервала температурное поле потока является функцией только координаты х и описывается уравнением:

dt (x)

+Pt (x)= V (x),

(8.16)

dx

где

K dτ+1

V (x)=

1

f (x)

(8.17)

P =

;

;

W dτ

F (x, dτ)+

d τ

W

f (х) – температурное поле теплоносителя в элементарной области в начальный момент. При начальном условии вида t (0)= t0 имеем решение уравнения (8.16):

						x			(8.18)
			t (x)= еxp (−Px) t0			+∫V	(x)еxp (P x)dx .
						0
	Средняя температура жидкости на участке длиной ∆х равна
		1	∆x	1	∆x		x
	t =		∫t(x)dx =		∫еxp (−P x) t0		+∫V (ξ)еxp (P ξ)dξ	dx.	(8.19)
		∆x		∆x
			0		0		0

Если канал образован одной стенкой замкнутого периметра с температурой tF (x, τ), то

αП

F(x, τ)=

α

ПtF (x, τ)

K = Sρc ;

Sρc

,

где П – периметр сечения канала; α – коэффициент теплоотдачи.

В стационарном случае задача упрощается:

dt (x)

+K t (x)= S (x),

где

dx

K = α1П1 +α2П2 ;

Gc

S(x)=

α1П1tF (x)+α

2П2tF

(x)

.

1

2

Gc

При начальном условии вида t(0)= t0 решение уравнения (8.21) имеет вид

x

t (x)= exp (−Kx) t0

+∫S (x)exp (Kx)dx .

0

Если

температуры стенок

постоянны, т.е. tF

(x)= const = tF

и

tF

(x)= сonst = tF

,

t (x)= V +(t0 −V )еxp (−K x),

1

1

2

2

(8.24)

где

V =

α1П1tF (x)+α2П2tF

(x)

(8.25)

1

2

.

α1П1 +α2П2

(8.20)

(8.21)

(8.22)

(8.23)

то

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

		1	∆x	(t0 −V )
	t =		∫t(x)dx = V +		(1−exp (−Kx)).	(8.26)
		∆x
			0	K∆x

Если канал образован одной стенкой замкнутого периметра с температурой tF (x), то

					x
				t (x)= exp (−K1x) t0	+K1 ∫tF (x)exp (K1x)dx	,
					0
где	K1 =	αП	.	(8.28)
		Gc

Если температура стенки постоянна, т.е. tF (x)=сonst =tF , то

t (x)= tF +(t0 −tF )exp (−K1x),

а средняя температура жидкости на участке длиной ∆х равна

t =	1	∆x t(x)dx = tF +	(t0 −tF )		(1−exp (−K1x)).
	∆x ∫		K	∆x
	0		1

(8.27)

(8.29)

(8.30)

Таким образом, решения задач (4.1) – (4.5), (5.1) – (5.6) и (8.12) – (8.13) составляют минимальный набор аналитических решений задач теплопроводности, необходимый для математического моделиро-

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

9 РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ

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

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

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

Классическая инженерная методика расчета коэффициентов конвективной теплоотдачи, базирующаяся на теории подобия, основана на использовании критериальных уравнений алгебраического типа, обобщающих экспериментальные данные по различным веществам, выступающим в роли теплоносителей, для каждого набора условий протекания теплообменного процесса. Поэтому использование критериальных уравнений, являющихся по сути результатом многомерной аппроксимации, приводит в каждом конкретном расчете к погрешностям, не поддающимся оценке. Как правило, погрешность расчета коэффициентов конвективной теплоотдачи по критериальным уравнениям составляет 30 – 50 %. По замечанию А.А. Гухмана, высказанному в частной беседе, эта погрешность может достигать 900 %.

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

Существует ряд методик расчета коэффициентов теплоотдачи по экспериментальным данным, это:

1)расчет непосредственно по определению коэффициента теплоотдачи как удельного количества тепла, приходящегося на единицу площади поверхности теплообмена в единицу времени и отнесенного

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

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

3)прямой расчет коэффициента теплоотдачи по результатам решения обратной задачи теплопроводности для исследуемого процесса.

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

Второй способ является «кустарным», т.е. не имеет строгого математического обоснования, хотя он относительно прост и нагляден. Этот способ целесообразно использовать в оценочных расчетах. Третий способ математически строг, но сложен и специфичен. Решению обратных задач теплопро-

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

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

В[7] приводится методика, основанная на использовании интегральных преобразований Лапласа для решения обратной задачи теплопроводности в цилиндрических координатах.

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

Пусть измеряется температура во времени в какой-то точке пластины. Температурное поле такой пластины описывается решением следующей задачи теплопроводности:

∂t (x, τ)		= a		∂2t (x, τ)				, 0 ≤x ≤R,	τ≥0;	(9.1)
∂τ					∂x2
с начальным условием
	t (x, 0)= t0						= сonst			(9.2)
и граничными условиями
λ	∂t (R, τ)				+α(t (R, τ)−tc )=0 ;					(9.3)
	∂x
			∂t (0, τ)				=0 ,			(9.4)
					∂x

где х – координата, направленная по нормали к поверхности пластины; τ – время; t(x, τ) – температурное поле пластины; a = cλρ – коэффициент температуропроводности материала пластины; λ, с, ρ – теп-

лопроводность, теплоемкость и плотность материала пластины соответственно; t0 – начальная температура пластины; tс – температура окружающей среды; R – полутолщина пластины; α – коэффициент теплоотдачи.

Путем замены переменной T (x, τ)= t(x, τ)−t0 можно получить задачу с нулевым начальным распреде-

лением:	∂T (x, τ)=a	∂2T (x, τ)
			;	(9.5)
	∂τ	∂x2
	T (x, 0)= 0 ;			(9.6)
λ	∂T (R, τ)+α(T (R, τ)−tc +t0 )= 0 ;			(9.7)
	∂x
	∂T (0, τ)	= 0 .		(9.8)
	∂x
Выполним преобразование Лапласа данной задачи по формуле
	∞
	S (x, p)=∫T (x, τ)exp(− p τ)dτ.			(9.9)
	0

Получаем задачу в изображениях:

	a	d 2S (x, p)		− p S (x, p)= 0 ;							(9.10)
		dx2
	∂S	(R, p)				tc +t	0			= 0 ;	(9.11)
λ
		∂x	+α S (R, p)−			p
		∂S (0, p)			= 0 .						(9.12)
			∂ x

Решение этой задачи в изображениях имеет вид

S (x, p)=		p			p		,	(9.13)
	Ash		x	+B сh		x
		a			a

где sh(z) и ch(z) – гиперболические функции.

A и B определяются из граничных условий (9.3), (9.4):

λ

p

p

R

+ B sh

p

R

+

a

Ach

a

a

(9.14)

p

p

t

−t

+B

−α

c

0

= 0;

+ α Ash

a

R

сh

a

R

p

λ

p

p

0

p

0

= 0 ,

(9.15)

a

Ach

a

+B sh

a

А = 0,

откуда

(9.16)

B =

α(tc −t0 )

.

(9.17)

p

p

p

p λ

sh

R +αсh

R

a

a

a

α(t

)сh

p

−t

x

a

c

0

Тогда

S (x, p)=

.

(9.18)

p

p

p

p

λ

sh

R

+αсh

R

a

a

a

Если известно изображение температурного поля S(x, p), значение коэффициента теплоотдачи может быть найдено аналитически:

S (x, p)pλ

p

p

sh

R

a

a

α =

.

(9.19)

p

p

(t

c

−t

0

)сh

x

− p S (x, p)сh

R

a

a

Изображение температурного поля S(x, p) определяется по формуле