Нечаев Моделирование процессов ядерной технологии 2007

ля и экранов не меняются во времени, т.е. потоки тепла между всеми объектами равны между собой:

Qн,1э = Q1э,2э = ... = Q(N–1)э,Nэ = QNэ,c = Qэ.

Поэтому в рассматриваемой задаче число неизвестных величин принимаем равным N+1, а именно: N температур экранов плюс одно значение величины стационарного потока.

Для каждой пары соседних поверхностей с номерами i, i+1 (в том числе поверхности нагревателя (i=1) и стенки камеры (i=N+2)) запишем через их свойства N+1 уравнений для расчета теплового

потока Qэ в системе:

Qэ = Сi, i+1. [(0.01. Ti)4 – (0.01. Тi+1)4] (i=1,2,...N+1) [Вт], (4.2.9)

где

Сi, i+1 =

5.67 πDiH

.

1

Di

(

−1) +(

1

−1)

+1

ε

(T )

ε

i+1

(T

)

D

i

i

i+1

i+1

4.3. Алгоритмчисленногорешения

Общая методика решения. Численные методы. Проблемы численного решения. Начальное приближение. Выбор рекуррентной записи. Окончание счета.

Если бы степень черноты не зависела от температуры, то система (4.2.9) весьма просто решалась бы аналитически. Поделив левые и правые части уравнений (4.2.9) на свои постоянные в этом случае коэффициенты Сi, i+1 и просуммировав их, для расчета величину потока тепла Qэ в системе с экранами получаем формулы, в которые входят только температуры нагревателя и стенки:

.

N +1

1

.

4

.

4

Qэ

∑

= [(0.01

Tн)

– (0.01

Тс)

],

(4.3.1)

C

+1

i=1

i,i

.

4

.

4

N +1

1

или

Qэ = [(0.01

Tн)

– (0.01

Тс)

] / ∑

.

(4.3.2)

C

+1

i=1

i,i

Tемпературы экранов Ti+1 (i = 1,2,...N), необходимые для коррекции выбора материала для их изготовления, последовательно

рассчитываются по N уравнениям системы (4.2.9):

Тi+1 = 100. [Сi,i+1. (0.01. Ti)4 – Qэ]0.25. (4.3.3)

Но поскольку степень черноты имеет температурную зависимость: εi = ai + bi. Ti , то при расчете тепловых потерь излучени-

51

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

Одним из простейших численных методов, а именно методом

последовательных приближений, называемом также методом про-

стых итераций. Этим методом можно пользоваться, если уравнение вместо обычного в классической математике вида F(X) = 0 записать в рекуррентной формe: X = f ( X ).

Из (4.2.9), (4.3.1), (4.3.2) следует, что в системе с экранами можно представить в виде рекуррентных функций двa типа величин:

1)величину потока тепла: Qэ = Q{ Ci { ε [Ti(Qэ) ] } } ;

2)температуры экранов : Ti = Т{ Q { C[ ε (Τi ) ] } } .

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

за нулевое приближение можно принять любые числа. Например, положить все Тi(0) равными (Т1+ТN+2)/2.

Задавшись некоторыми начальными значениями температур экранов Тi(0) (нулевое приближение), вычисляем последовательно :

εi(0) = ai + bi. Тi(0) ;

(4.3.4)

(0)

5.67 πDi H

Ci,i=1 =

;

(4.3.5)

1

1

Di

(

−1) +(

−1)

+1

εi(0)

εi(+01)

Di+1

Q(0)

= [(0.01. Tн)4 – (0.01. Тс)4] / N∑+1

1

,

(4.3.6)

э

i=1 C

(0)

i,i+1

и, наконец, получаем следующее (первое) приближение для значений температур:

T (1)	= 100. [С(0)i,i+1. (0.01. Ti(0))4 – Qэ(0)]0.25 .	(4.3.7)
i+1

Повторяя расчет по формулам (4.3.4)–(4.3.7) с новыми значениями температур, получаем все более точные значение для иско-

мого потока тепла Qэ(1) … Qэ(m) , Qэ(m+1) . В качестве условия нахож-

дения решения принимается условие достижения разностями значений расчетных величин на двух последовательных (m и m+1) итерациях величин, меньших некоторой наперед заданной точно-

52

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

mod( Qэ(m) – Qэ(m+1) ) < δ;

или относительной:

mod[( Qэ(m) – Qэ(m+1) )/ Qэ(m) ] < δ.

На практике ограничиваются относительной точностью в 1%.

Контрольныевопросы

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

2.Чем отличаются законы излучения для абсолютно черных и серых тел?

3.Сформулируйте и докажите закон Кирхгофа.

4.В чем заключается закон Ламберта?

5.В чем природа угловых коэффициентов? Сформулируйте для них правила: замыкаемости, взаимности, невогнутости.

6.Почему тонкий экран эффективно уменьшает теплопотери?

7.Выведите формулы теплообмена излучением в системе коаксиальных цилиндров.

8.Что такое «приведенные коэффициенты»?

9.В чем заключается общая методика решения системы уравнений для теплообмена излучением в системе коаксиальных цилиндров?

10.Каковы основные проблемы численного решения?

53

5. НЕСТАЦИОНАРНЫЙ ТЕПЛООБМЕН

5.1. Математическаямодельнестационарноготеплообмена

Второй закон Фурье (Фика). Коэффициент температуропроводности. Его физический смысл и связь с теплопроводностью. Условие теплоотдачи (массоотдачи) по Ньютону.

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

Для простоты рассуждений будем исследовать одномерный вариант. Рассмотрим некоторый элементарный объем около точки xо, который, имея размеры ∆x·∆S, содержал количество тепла :

Т·Cp·ρ ·∆x·∆S, (5.1.1)

где Cp – теплоемкость, [кДж/кГ/К], ρ – плотность, [кГ/куб.м]. Из-за разности градиента температуры у левой и правой стенки

элемента за время ∆τ в этот элемент вольется тепла:

	∂T	∆S·∆τ,	(5.1.2)
Q+ = –λ·
	∂x x=x
	0

а выйдет из него другое количество:

	∂T	·∆S·∆τ .	(5.1.3)
Q – = –λ·
	∂x x=x +∆x
	0

Изменение теплосодержания ∆Q = Q+ – Q– дает приращение температуры ∆Т, вычисляемое из уравнения:

∆Q = Q+ – Q –= ∆Т·Cp·ρ·∆x·∆S. (5.1.4)

Выражения (5.1.2) и (5.1.3) отличаются только тем, что относятся к разным точкам пространства. Тогда, используя прием разложения функций в ряд Тейлора, выразим значение Т(xо+∆x) через Т(xо), ограничившись первым членом разложения:

54

∂T

·∆x +

... ,

(5.1.5)

Т(xо+∆x) = Т(xо) +

∂x

x=x

0

откуда для первой производной имеем:

∂T (x

0

+ ∆x)

=

∂T (x

0

)

+

∂2T (x

)

·∆x

+ ...

(5.1.6)

∂x

0

∂x

∂x2

Подставляя (5.1.6) в (5.1.3), а (5.1.2) и (5.1.3) в (5.1.4), имеем:

∆Q = –λ·

∂T (x0 )

·∆S·∆τ

· ∂T (x0 )

·∆S·дτ –λ·

∂2T (x0 )

·∆x·∆S·дτ) =

∂x

–(–λ

∂x

∂x2

= λ·

∂

2T (x

)

·∆x·∆S·∆τ = ∆Т·Cp·ρ·∆x·∆S .

(5.1.7)

0

∂x2

Разделив во второй строке равенства (5.1.7) левую часть на правую и проведя перегруппировку членов, получим уравнение в частных производных, которое описывает изменение температуры во времени в каждой точке пространства:

∆T	=	∂T	=	λ	∂2T .	(5.1.8)
∆τ		∂τ		Cpρ ∂x2

Из аналогичных рассуждений для концентраций можно полу-

чить II закон Фика – уравнение нестационарной диффузии:

	∂C	= D·	∂2C	,	(5.1.9)
	∂τ		∂x2

где D – коэффициент диффузии, С – концентрация. Для полной аналогии уравнений диффузии и теплопроводности перепишем

(5.1.8) в виде:

∂T	= а	∂2T	,	(5.1.10)
∂τ		∂x2

где через а обозначен коэффициент температуропроводности:

	a =	λ	,	(5.1.11)
		Cpρ

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

Уравнения (5.1.9) или (5.1.10) описывают связь функции – температуры или концентрации – с аргументами – временем и координатой – в самом общем виде без учета конкретных условий нагрева

55

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

Если такая связь линейна, то такой тип граничных условий называют условием Ньютона. Для тепловых задач граничное условие по Ньютону записывается как уравнение:

	–λ·	∂T (x = S,τ)	= α·[Тсреды – Т(x=S,τ)]·k ,	(5.1.12)
		∂x

где α – коэффициент теплоотдачи, показывающий количество теплоты, получаемой (или отдаваемой) единицей поверхности в единицу времени при разности температур поверхности и окружающей среды в один градус. Через k в (5.1.12) обозначен некоторый коэффициент пропорциональности, учитывающий разницу величин:

a) потока тепла, поступающего на поверхность тела извне:

qпов(τ) = α·[Tсреды – Т(x=S,τ);

б) потока тепла, уходящего от поверхности тела в его глубину:

qвнутр(τ) = –λ· ∂T (x∂= S,τ) . x

Именно различие величин этих потоков и делает задачу нестационарной – некоторая часть тепла накапливается (или уходит) из слоя толщиной dx, изменяя его температуру во времени.

Для диффузионных задач граничное условие по Ньютону записывается аналогично:

jпов = β·(Cвнешн – Cпов),

(5.1.13)

где β – коэффициент массоотдачи.

5.2. Моделированиенагреватвердыхтел

Тонкие имассивные тела. Вывод критерияБио. Его физическийсмысл. Лимитирующая стадия. Вывод уравнения нагрева тонкого тела. ЗадачаЭйлера. Качественнаякартинанагревамассивноготела.

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

56

	∂T
q = λ·			= α·(Тпечи – Тпов).	(5.2.1)
	∂x пов

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

половину толщины тела L:
λ·(Тпов – Тцентр)/L = α·(Тпечи – Тпов).	(5.2.2)

Еще более загрубим (5.2.2), заменив разность температур на максимальную температуру в системе среда – тело (разность двух

величин всегда меньше максимального из них):
λ·Тпечи/L =α·Тпечи .	(5.2.3)

Равенство (5.2.3) выглядит еще более справедливым, если учесть, что в конце процесса нагрева наступит состояние, в котором:

Тпов = Тцентр = Тпечи .

(5.2.4)

Поделив правую часть (5.2.3) на его левую, получим безразмер-

ную величину, называемую безразмерным критерием Био Bi :

	Bi =	α L	,	(5.2.5)
		λ

смысл которого легко понять, если (5.2.5) переписать в виде:

Bi =	L / λ	=	Internal Heat Resist .
	1/ α		External Heat Resist

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

Если критерий Био очень большой (Bi > 10), то можно решать уравнение теплопроводности с граничным условием 1-го рода, приравняв Тпов и Тпечи с самого начала процесса нагрева.

Если критерий Био мал (Bi < 0.25), то внутреннее тепловое сопротивление мало, градиенты температуры внутри него малы, тело

57

хорошо прогревается, и перепада температур между центром и поверхностью почти не образуется. Такое тело называется тонким в тепловом отношении. Т.е., оно воспринимает весь поток, который на него приходит – это соответствует граничному условию 2-го рода. Для изделий, имеющих большой коэффициент теплопроводности, малое значение критерия Био Bi получается даже при их большой физической толщине. В то же время, если поместить "тонкое" в одних условиях тело в более мощный источник тепла, то критерий Био увеличится и тело надо будет рассматривать как

массивное.

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

Выведем формулу для расчета времени нагрева тела от начальной температуры То до заданной температуры ТК тонкого тела массой М, теплоемкостью Ср и с тепловоспринимающей поверхностью F. В произвольный момент времени тело, имеющее температуру Т за элемент времени dτ получит из окружающей среды элементарное количество тепла δQ:

δQ = a·(Tпечи – T(τ)) ·F dτ ,

которое нагреет тело на dT градусов:

dT·(M·Cp)= a·(Tпечи – T) ·F dτ . (5.2.6)

Перепишем (5.2.6) в форме уравнения для dτ :

	M C p
dτ =			dτ.	(5.2.7)
	a (T	−T )
	печи

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

58

	τ	M Cp	TK	dτ
	∫dτ =		∫		,
		α F		Tпечи −T
	0		T0

откуда для времени нагрева получаем следующую формулу:

	M Cp		T	−T
τ =		ln	печи	0	.	(5.2.8)
	α F		T
				−T
			печи	K

Как видим из (5.2.8), время нагрева тонких тел во многом зависит от разности температур источника (печи) и конечной (требуемой) температуры. Для нагрева тела до нужной температуры за конечное время необходимо, чтобы температура печи была выше, чем требуемая.

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

Впервые уравнение в частных производных, иначе – уравнения математической физики, составил и решил в Петербурге Леонард Эйлер, рассмотрев задачу о промерзании грунта при годичных колебаниях температуры. Как известно, начиная с глубины 30 метров температура Земли практически постоянна и не зависит ни от времени года, ни от географического положения. Однако сезонные колебания температуры на поверхности Земли проникают на некоторую глубину с течением времени. Эйлер к уравнению нестационарной теплопроводности (5.1.10) записал следующие граничные условия.

На поверхности Земли (х = 0) температура изменяется периодически, по синусоидальному закону: T(0,τ) = sin(τ).

Вглубине Земли (условно х= ∞ ) температура равна некоторой константе: T( ∞ ,τ) = const .

Вначальный момент времени примем, например, что температура всюду постоянна и также равна const:

T(x,0) = const.

(5.2.9)

Обращаем внимание, что при больших временах решение вообще не зависит от начальных условий!

Эйлер предложил решать уравнения в частных производных

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

59

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

T(x,τ) = ϕ(x)·f(τ). (5.2.10)

Подставляя (5.2.10) в уравнение теплопроводности (5.1.10), получаем:

	ϕ(x)	∂f (τ)	= a· f(τ)	∂2ϕ(x)	.	(5.2.11)
		∂τ		∂x2

Разделяем переменные, поделив правую и левую части (5.2.11) на

(5.2.10):

∂f (τ) 1		/f(τ) = a	∂2ϕ(x) 1		= –ξ	2	,
∂τ	f (τ)		∂x2 ϕ(x)

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

∂f (τ)		– ξ·f(τ) = 0		с условием	f(0) = const;
	∂τ
a·	∂2ϕ(x)		– ξ·ϕ(x) = 0	с условиями:	ϕ(0) = sin(τ);
	∂x2

ϕ( ∞ ) = const,

методы решения которых достаточно разработаны.

5.3. Безразмерныекритерии

Метод Гухмана. Безразмерные время и координата. Критерий Фурье. Безразмерная температура. Особенности нагрева цилиндрических тел. Практический пример.

Заметим, что условие (5.2.10) очень редко выполняется в элементарных функциях (как, например, в задаче Эйлера с синусоидальной функцией). Решение обычно находится в различных классах специальных (неэлементарных) функций, представляемых бесконечными рядами разложения по ортогональным функциям или полиномам. Для плоских задач это прежде всего ряды Фурье по тригонометрическим функциям. Для цилиндрических задач более удобными, дающими более быстро сходящиеся числовые ряды, оказываются функции Бесселя. Последние часто так и называются цилиндрическими функциями. Для задач со сферической симметри-

ей следует пользоваться сферическими полиномами Лежандра.

60