Численное решение краевой задачи двумерного температурного поля с внутренними источниками теплоты Текст научной статьи по специальности «Физика»

УДК 621.1

ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ ДВУМЕРНОГО ТЕМПЕРАТУРНОГО ПОЛЯ С ВНУТРЕННИМИ ИСТОЧНИКАМИ ТЕПЛОТЫ

Докт. техн. наук, проф. ЕСЬМАН Р. И.

Белорусский национальный технический университет

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

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

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

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

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

Пусть размеры четверти формы а, х, Ь, а четверти плиты а0, х, Ь0, (рис. 1). Внутренняя поверхность формы покрыта масляной эмульсией

толщиной 5покр. Между плитой и формой возможно образование воздушного зазора, переменного вдоль линии контакта и изменяющегося во времени 5(У).

у у

Область I

I

О, Область ЛЛ

__Т

Область I

X

Т Ч

: >1'.1.-|.,.'1 1: х

1

Рис. 1. Эпюры температур и профильных температурных напряжений в сечениях металлической формы

Поле температур в плите и форме описывается дифференциальными уравнениями [1]:

с (Т )р (Т ) дТ1(Х,у,*) =-^^^ д1 дх

\ (Т)

дТД х, у, *) дх

_д_ ду

\ (Т1)

дТ1( х, у, *) ду

; (1)

С2(Т2)Р2(Т2 )

дТ2( х, у, *) д

д

дх

^2 (Т2 )

дТ2( х, у, *) дх

ду

^2 (Т2 )

дТ 2 (х, у, *) ду

(2)

где с1, р1, Х1, Т1 и с2, р2, X2, Т2 - теплофизические характеристики и температуры плиты и формы соответственно.

Уравнение (1) решается в прямоугольной области (0 < х < а0; 0 < у < Ь0), а уравнение (2) - в сложной области в виде угла, получаемого при вычитании из области (0 < х < а; 0 < у < Ь) области, занятой плитой.

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

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

Х дТ1 . дТ2

—Х1-= —Х2-

дх дх

при х = а0; 0 < у <

[ТД у, *) — Т2( у, *)]

х„

5( у, *)

-«л (у, *)

X

покр

покр

покр

Хв

5( у, *)

"ал (У, *)

(3)

Х дТ1 Х дТ2

—Х1-= —Х2-

ду ду

при у = Ь,; 0 < х < а0,

[Т1( х, *) — Т2( х, *)]

5( х, *)

-«л (х, *)

Х„

покр

Х„

Зпокр 8(х, *)

-«л ( х, *)

(4)

X

где Xпокр, Хв - теплопроводность покрытия и воздуха; 5(х, t) - зазор в контакте у = Ь0 в момент времени 5(у^) - то же х = а0 в момент времени V,

а л =в 1/2а(Т + Тп0кр )(Т? + Тп20кр), (5)

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

Т 2

Тпокр = Тр ^-1. (6)

Хпокр + 5 покр ^ + ал 1 5 л)

Хпокр 5 покр +Х. 5 - + ал

Значения Т1, Т2 и 5 берутся в соответствующей точке контактной поверхности, в которой определяются Тпокр и а л.

Напомним, что при записи (3) и (4) учитывали только процесс теплопроводности через покрытие и процессы теплопроводности и теплоизлучения через слой воздуха. Процессом конвекции в зазоре пренебрегаем.

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

дТ2

—1 = —2 = 0 при у = 0; дх дх

дТ1 дТ 2

—1 = —- = 0 при х = 0. ду ду

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

дТ.

-^^т2 = а(Т2 -Тш) при х = а; дх

дТ

= а(Т2 -Тш) при у = b,

дх

где Тш - температура внешней среды; а - коэффициент теплоотдачи с наружной поверхности формы.

Коэффициент а определяется способом охлаждения наружной поверхности формы.

При свободном охлаждении формы в безграничном пространстве коэффициент а характеризует собой теплоотдачу свободной конвекцией и излучением а = ак +ал. Тогда запишем:

ак = / (Ог,Рг);

ал =га(Т2 + Т2)(Т 2 + Тш).

При вынужденном охлаждении формы

а = / (Яе,Рг),

причем число Re вычисляется по толщине охладительной рубашки формы. Начальные условия для уравнений (1), (2) запишутся следующим образом:

ТДх, у 0) = Т1о; Т2(х, у 0) = Т2о,

где Т^ - температура заливки; Т2 - начальная температура равномерно

прогретой формы.

В качестве характерного размера выберем длину формы а, а в качестве характерной температуры - температуру окружающей среды Т0 = Тш. Перепишем задачу в безразмерных переменных:

ди

ди

д дх ^ 1 дх ) ду I 1 ду

Л ди К —

(7)

при 0 < х < а0; 0 < у < Ь0;

ди д

ди

а С2Р2^~ = ^г\ К2^Т 1 + ^Т

\ 2 I К 2-

д дх ^ дх) ду ^ ду

при 0 < х < 1; Ь0 < у < Ь; а0 < х < 1; 0 < у < Ь0;

(8)

- ди ди

-К1-= -К2-= —

дх дх

(и — ^ + ал

покр

Кпокр + 5покр 5

а

при х = а(); 0 < у < Ь();

(9)

г, ди ди

—К1-= —К2-= —

ду ду

к + ал ^а 5 л)б„

(и —и) — 1"покр

покр

при х = Ь0; 0 < х < 1;

Кпокр +

5 5

покр

а

ди ди

— = — = 0 при х = 0; дх дх

(10)

(11)

ди ди

= = 0 при у = 0;

ду ду

(12)

—К2 — = ааи при х = 1; дх

-Л2— = ааи при у = Ь; (14)

ду

и = ип при t = 0; ]

0 \ (15)

и = и0 при t = 0,

где и и и - безразмерные температуры в областях I и II соответственно. В уравнениях (9), (10): 5П01ф Ф 0; 5 Ф 0.

При 5 = 0; 5покр Ф 0 их следует заменить на условия соответственно:

дм, дх = -Л2 дм2 дх = (и -и) Лпокр 5покр

дм, ду = -Л2 ди2 ду = (и -и) Лпокр 5покр

при 5покр =0; 5 Ф 0 - на условия:

ди, дх . ди2 - -Л2- дх = (и ■ -и) ( Лв ч 5

ди, ду ди2 - -Л2- ду = (и -и)( X

(16) (17)

наконец, при 5покр =0 и 5 = 0 будем иметь:

-дм ди

Л, - = Л 2 - , М = О;

дх дх -дм ди

Л,-= Л 2-, М = и.

ду ду

Значения 5 в равенствах (9), (10), (16), (17) определяются исходя из рассмотрения упругих деформаций формы.

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

I

Л1Т8гай(Т1)| ^ -Л,ж8^(7; ^ = -ф—•

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

Р1 |>,+ *(т,- Тф )]дТ-=ддх (Л, ('«)

где

Р:, с1, =

[р1Г, с1Г, \т при Т1 ^ Тф; |р1ж, с1ж, пРи Т1 ^ тф.

Решение (18) производят путем сглаживания 5 -функции и теплофизи-ческих коэффициентов, осуществляя замену фронта фазового перехода на некоторую его область (Тф - А; Тф + А).

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

и у = —Ь0 лежат на узлах сетки. Пусть N и N - число узлов по горизонтали и вертикали соответственно, тогда шаг по горизонтали к = 1/ Ы1, а по вертикали - к2 = Ь /Ы2.

Предположим, что горизонтальная строка узлов на контактной поверхности имеет номер М1, а вертикальный столбец - М2. Будем решать задачу на фиктивной сетке с узлами:

1

х =\г+2)ъ; Уг =11 + 21-2

1

при г = —1,0,..., N1; j = —1,0,..., N2.

Неявные конечно-разностные уравнения, соответствующие выражениям (7)-(15), на узлах фиктивной сетки на шеститочечном шаблоне [1] в момент времени ^ = (I +1) т имеют вид

1+1 1

2 (1) (1) иг,]■ иг,]■ 1

ас- Г" = к

(

'+1 „ '+1 Л

1

+—

'+1 '+1 ... ^(1) г+1, ] г, 1 — ^(1) иг, ] иг— 1, ]

V 41' к 4; к )

'+1 _ '+1 Л

'+1 '+1 - ■ - _

^(1) иг, 1+1 иг, 1 — ^(1) иг, 1 Пг, 1—1

4 Ь2 к

V 2 2 2 1 )

(19)

где 1 = 0,1,2,..., М1 — 1; у = 0,1,2,..., М2-1; I = 0,1,2,...;

и1+1 —и1 1

2 (2) (2) г, 1 иг, 1 1 ас-Р-—=к

г

V 2

V1+1 +1 V1+1 — V1+1 ^

^(2) г+1,1 г, 1 — ^(2) иг, 1 иг—1, ]

г+1,1 к 1 к

+

V

и'+1 -и+1 и'+1 — и'+1 ^

^(2) г, 1+1 г, 1 — ^(2) иг, 1 г, 1—1 г,;+2 к2 г'4 к )

где

г = 0,1,2,..., N1 — 1; | г = Мх, Мх +1,..., N1 —1;| 1 = М2, М2 +1,..., N2,1 1 = 0,1,2,..., М2 — 1, )

' = 0, 1, 2;

1

и1+1 -и 1+1 и+1 -и+1

ИМ„ / ИМ„ / -1 = -^(2) ^Ч,/ Чц —1, / =

М1 —^ ь м—l2,/ И1

= к.

/'.г+1 , .л+1 ..I+1 , ..I+1 Л иЫ!, 1 + М^, 1 —1 Чц, 1 + Чц —1,1

где 1 =

—1,0,1,..., М 2 — 1;

и 1+1 — и 1+1 и+1 — и+1

^(1) г,М2 1,М 2 —1 = ^(2) 1,М 2 1,М 2 —1

г,М2 — 1 ^ г, М 2 — 2 Й2

Л.г+1 , +1 +1 , +1 Л

= к.

(21)

где I = —1,0,1,..., М1 — 1;

к, =

покр

V51

а

1 X X

покр + ^в

5 5

покр в

пРи 5покр *0 и 5в *0;

а

(22)

X

к 1 а при 5П0Кр *0, 5 = 0;

покр

К =7^ + ал, при 5покр =0, 5*0

Соответствующие формулы могут быть выписаны для к,. Величины ал и ал определяются по температуре в контакте на предыдущем временном слое

а л, = е1/2СТ0

^ и1,М 2 + и1,М2 —1

У

■( ипокрг + 1)2

(23)

^ и,М2 + иг,М, —1

-и + 2

покр,

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

68

покр!

1+1 л иг,М2 +иг,М2 -1 ^покр М 2 + и1,М 2 -1 ^в Л /

2 5 покр 2 V5!

^покр 5покр 5

(24)

Уравнения (23), (24) следует рассматривать как трансцендентные для определения ал , ипокр , которые решаются методом половинного деления.

При 5покр = 0 выражение (23) следует заменить на

ал( = 81/2сТ0

Г и\,Ы, + и',Ы-, -1

л2 v

У

г,М, г,М,-1

Ги'

(25)

Выражения для вычисления ал , ипокр составляются по аналогии с

(23)-(25). Вычисление 5г- и 5у в (22)-(25) производится исходя из рассмотрения напряженного состояния формы. Оставаясь в рамках одномерного расчета напряжений, предполагаем, что 5 изменяется по параболическому закону:

ботах 1 = 5

1 -

Г ^

V М 2 + N2 ,

ботах 1 = 52

1 -

Г 2г +1 Л2 М1 + N1

где 52ах и 52ах - деформации формы в точках 1 = 0; у = М2; 1 = М1; у = 0 соответственно. Получаем:

5тах =- 3 тОЩ^М^)2 .

2 2 [(N2 - М2) Н2 ]3 Е'

5тах =-3 тОт2 (2М2Н2)2

1 2 [(N1 - М1) Н ]3 Е'

где тОт1 и тОт2 - максимальные моменты изгиба в продольной и поперечной балках формы, определяются по формуле:

Т I ■ М 2 + N2

1 = ^ Е811 у —2—1

1 =М 2

где в1 =

2 пРи j=м2 или 1 = ; 1 - в остальных случаях;

N1

= * I

тот,,

г =М1

М1 + N в,1 г--1-к 1С,-;

в =

2 при г = М, или г = N; 1 - в остальных случаях,

где сг и с ■ - напряжения в сечениях балок по осям, определяемые так:

С =-вЕ

и,О +4,-1

+1

Т0 -Т2о

1 ^

— ¿вД.Е

Л Л гГг

N - М1 ,=М,

и,о + и,-1

+1

Т0 -Т2о

С, =-Р ,Е

^0, 1 +и-1, j

Л

м,

1 2

+—1— 1вР ,Е

м2-М. 11

2 1 =М 2

т - т

100 22„

и0,1 +и-1,1

т0 - Т2о

где Рг и р 1 - коэффициенты температурного расширения материала, вычисляются по температурам (иг 0 +иг -1)/2 и (и01- +и-11- )/2 соответственно. При вычислении интегралов использовалась формула трапеций.

Действительные величины и 5г считаются равными нулю, если соответствующие температуры плиты (игМг - игМг- )/2 и (им„; + ищ; )/2

больше или равны температуре затвердевания. В противном случае зазоры 61 и 5г определяются как остаточные деформации формы с использованием приведенных выше формул.

Аппроксимация остальных граничных условий и начальных условий запишется:

'+1 '+1 при / = -1,0,1,..., М2;

и0,1 = и-1, ]

и+1 = Ц.+1,. при 1 = М2, М2 +1, ..., м2; I

(26)

0, ]

-1,1

и = и при I = -1, 0, 1, ..., М1; и',+0 =и1+-11 при г = М„ М1 +1, ..., М1;1

и'+1 -и1+1 -^(2) иМ!,у uм1-l,j

N1 - 1 ] 12

= а 1а

и'+1 +и+1

N1,1 -1,1

где 1 = -1,0, ..., N2;

2 ?

и+1 _ и+1 +1 , +1

-Я(2) 1 ^2 Ц^-1 = аа-

г ,M2 -2 ^

2

где г = -1,0,..., И,;

где г = -1,0,...,М,; ] = -1,0,...,М2;

где г = М,, М, +1,..., И,; 7 = -1,0,..., М2 или /' = -1, 0,..., И,; 7 = М2,М2 +1,..., И2.

Выражения (19)-(21), (26), (27) с учетом (28) дают (И, + 2)(И2 + 2) + + 2 (М, + М2) алгебраических линейных уравнений для определения такого же количества неизвестных значений температур в узлах сетки. На каждом временном шаге (д + 1)т, I = 0, 1, ... их решение производится по методу продольно-поперечных направлений.

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

В Ы В О Д

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

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

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

Л И Т Е Р А Т У Р А

1. Е с ь м а н, Р. И. Расчеты процессов литья / Р. И. Есьман, Н. П. Жмакин, Л. И. Шуб. -Минск: Вышэйш. шк., 1977. - 264 с.

Представлена кафедрой промышленной теплоэнергетики

и теплотехники Поступила 3.03.2008