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

УДК 536.2

О. Ю. Чигирева

МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНОГО ПРОЦЕССА ТЕПЛОПРОВОДНОСТИ В СОСТАВНОМ ЦИЛИНДРЕ ПРИ ЛОКАЛЬНОМ ПЕРИОДИЧЕСКОМ ТЕПЛОВОМ ВОЗДЕЙСТВИИ

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

E-mail: [email protected]

Ключевые слова: математическая модель, нестационарный процесс теплопроводности, цилиндр, локальное периодическое тепловое воздействие, температурное поле, неидеальный тепловой контакт.

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

Постановка задачи и математическая модель процесса. Рассматривается нестационарный процесс теплопроводности в цилиндре радиусом г0, состоящем из двух слоев толщиной d и Н (рис. 1). Разогрев составного цилиндра осуществляется осесимметричным круговым (0 ^ г ^ г*) локальным периодическим источником теплоты, действующим на нижнем основании цилиндра, с плотностью теплового потока

nr

q (r,t) = qo ( 1 + cos — J ■ f] (t), t > 0, 0 < r < r*,

Рис. 1. Осевое сечение составного цилиндра

Рис. 2. График функции г/О

где г* — радиус пятна теплового контакта; п (¿) — ступенчатая периодическая функция с периодом по времени равным (рис. 2). Кроме того, на нижнем основании цилиндра вне зоны воздействия источника теплоты, происходит теплообмен излучением. Верхнее основание цилиндра охлаждается внешней средой по закону Ньютона с коэффициентом теплоотдачи а. Боковая поверхность цилиндра теплоизолирована. В начальный момент времени температура составного цилиндра постоянна и равна температуре внешней среды Тс. Предполагается, что тепловой контакт между слоями цилиндра является неидеальным [4], а теплофизические свойства материалов слоев зависят от температуры.

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

/ГТ1 dTi 1 d Л /ГТ1, дТЛ d Л /ГТ1, dTi

<T) "дГ = rsr (Al <T)r-w) + &(Ai <T) "z

t > 0, 0 < r < r0, 0 < z < d;

"T2 1 д / дТ2

р2С2 <T2) Ж = ГдГ A2 <T2)

д / дТ2

+ д* A2 <T2) "T

-Ai <Ti) -A2 <T2)

t> 0, 0 < r < r0, d < z < d + h; Ti <r, z, 0) = Tc, 0 < r < ro, 0 < z < d; T2 <r,z, 0) = Tc, 0 < r < ro, d < z < d + h; q <r,t) , t> 0, 0 < r < r,;

"Ti I

dz z=0 1

дТ2

dz z=d+h

dTi

dr r

dT2

dr r=r

ae <T4 - T4 <r, 0,t)), t > 0, r, < r < ro;

(1)

(2)

(3)

(4)

(5)

= a (T2 (r, d + h, t) - Ta), t > 0, 0 < r < ro; (6) = 0, t > 0, 0 < z < d; (7)

(8)

r=ro

= 0, t> 0, d < z < d + h;

- Л! (Ц)

dz

= -(Т (r,d,t) - T2 (r,d,t)) =

л —

z=d

- Л (T) dT2 -Л2 (T2) öz

, г > о, о < г < г0. (9)

Здесь приняты следующие обозначения: значения индекса 3 = 1 и 3 = 2 соответствуют слоям 1 и 2 составного цилиндра; Т (г, г, г) — искомые температурные поля; р, с и Л — плотность, удельная теплоемкость и коэффициент теплопроводности соответственно; а — постоянная Стефана-Больцмана; е — степень черноты излучающей поверхности; Я — термическое сопротивление на поверхности контакта.

Построение алгоритма приближенного решения. Для нахождения приближенного аналитического решения задачи (1)-(9) воспользуемся модификацией [5] метода, предложенного в работах [6,7]. Умножим обе части уравнений (1) и (2) на г и введем функции

Оз (Тз,г) = РзгС (Тз), Лз (Тз,г) = гЛз (Тз), 3 = 1, 2. Тогда уравнения (1) и (2) можно записать в виде дт1

Ог (Тг,г) = а1у (Л1 (Тг,г) gradТг), г > 0, 0 < г < г0, 0 < г < d;

г (10) -т2

О2 (т2,г) ^ = div (Л2 (Т2,г) gradТ2),

г > 0, 0 < г < г0, d < г < d + Н, (11)

где операции div и grad следует рассматривать как операции в прямоугольной декартовой системе координат (г, г).

Проведем дискретизацию временной переменной г системой точек = кт, к = 1, 2,... с достаточно малым шагом т > 0 и заменим в уравнениях (10), (11) производные по времени конечно-разностными отношениями

ÖTj T (r,z) - (r,z)

k = 1, 2,..., j = 1, 2,

-г т

где Т^ (г, г) — приближенные значения функций Тз (г, г, г) при г = Следует отметить, что согласно начальным условиям (3) и

(4) Т(0) (г, г) = Тс .

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

ОГ (г, г ) = Оз (т**-4 (г, г) ,г) , Л$й) (г, г) = Лз (т**-4 (г, г) ,г) , 3 = 1, 2;

д(к) (г) = q (г, ¿к).

Кроме того, на каждом временном слое £ = тепловые потоки в граничных условиях (5), (6) и условии сопряжения (9) вычислим, используя функции Т.(к 1 ) (г, г), найденные на временном слое £ = ¿к-1 :

rq(k) (r), 0 ^ r ^ r,

<r) =

aer ( T4 -

T(k-i) (r, 0) Q2k) <r) = ar (Vf-^ <r, d + h) - Tc) ; Q0k) (r) = R (zf-i) <r,d) - T2k-i) <r,d)) .

r, < r < r0;

Это позволяет записать дифференциально-разностный аналог начально-краевой задачи (1)-(9) в виде итерационной схемы (к =

= 1, 2,

решения двух краевых задач для линеиных эллиптических

уравнений с переменными коэффициентами Cj ) (r, z), Aj ) (r, z): - div (л<к) (r, z) gradT(k) (r, z)) + 1 Cf) (r, z) Tf0 (r, z) =

= 1 Cjk) (r, z) T(k-i) (r, z), 0 < r < r0, 0 < z < d; (12)

T

dTi

(k)

dr

= 0,

dTi

(k)

r=0

dr

= 0, 0 < z < d;

r=ro

-лik) (r, z)

dT(

(k)

dz

= Q ik) (r),

-Л ik) (r,z)

dT(

(k)

dz

z=0 (k)

(13)

(14)

= Q0k) (r), 0 < r < r0;

z=d

- div (Л2к) (r, z) gradT2k) (r, z)) + 1 C2k) (r, z) Tf (r, z) =

= - C2k) (r, z) Tf_i ;(r,z), 0 <r<r0, d < z < d + h; (15) т

- )

dT

(k)

dr

= 0,

dT2(k)

r=0

dr

= 0, d < z < d + h;

r=ro

^ „, "TW

-Л2к) (r,z)

dz

= Q0k) (r),

z=d

(k) _\ dT2

(k)

-Л2к) (r,z)

dz

(16)

(17)

= Q2k) (r), 0 < r < r0.

z=d+h

Отметим, что задачи (12)-(14) и (15)-(17) на каждом временном слое t = tk решаются независимо.

На k-м шаге итерации функции T^ (r, z) будем искать в форме разложения в двойные тригонометрические ряды Фурье

X те

j (r z) = Y1 Y1 Ymn^Amn (r, z) , j = 1 2,

m=0 n=0

Г 0,5, m = 0;

7mn = 7m7n, 7m = ^ . 0 по полным и ортогональным [5J в

I 1, ^^ > 0,

областях Qj, j = 1, 2, системам функций {Xmn (r,z)}X,n=0:

Xi,mn (r, z) = cos (pmr) cos (vi,nz), Qi = (0,ro) x (0,d);

X2,mn (r, z) = cos (^mf) cos (v2,n (z - d)), Q2 = (0,ro) x (d, d + h), mn nn nn

где ^m = -, v1,n = ~T, V2,n = ~T.

r0 d h

Для нахождения коэффициентов aj^ этих разложений умножим уравнения (12) и (15) на функции X1ps (r, z) и X2 ps (r, z) соответственно, а затем проинтегрируем полученные соотношения по областям Qi и Q2. С учетом граничных условий (13), (14) и (16), (17) приходим к бесконечным системам линейных алгебраических уравнений относительно искомых коэффициентов ajmn:

V^ V A(k) Y a(k) - b(k) p - n 1 2 s - n 1 2 1-12

/ j / j Aj,psmn fmnuj,mn — bjj,ps , p — n , 1 2 , • •• , s — n , 1 , 2 , • •• , J — ± , 2 ;

те те

1(й) ^ о^ = 6(й)

т=0 п=0

(18)

где

Л00 = т (и и + ^ ^ ) г +

+ т (итир — (^з()т-р|,п+в) — ^3, (т+р,|п_в|)) +

+ (*0 + (*0 + (*0 + (*0 + П3|(|т_р|,|п—в|) + Пз'|(т+р,|п_в|) + П3>(|т_ р|,п+в) + П3>(т+р,п+в);

те те

ь( *) = .( *о + ^ ^_г)х

ьз',рв = Л 3 + / у / у /тп°31тп л

т=0 п=0

/ ( *0 + ( *:) + ( *0 + (*0 V

\П3|(|т_р|,|п_в|) + П3|(т+р,|п_в|) + П3|(|т_р|,п+в) + П3>(т+р,п+в) I ;

=8т ((-1)^+г яро+)), /2 ^=8т ((-1)5+г ^+г).

Здесь <3(Р и г/3'(Р в) — коэффициенты Фурье разложений функций Л^) (г, г) и О((г, г) соответственно в двойные тригонометрические

ряды по системам функций {Xjips (r, z=0 ; ф,к), #Pk) и — коэффициенты Фурье разложений функций Qlk) (r), Q0k) (r) и Q2k) (r) соответственно в тригонометрические ряды по системе функций

{COS (^pr)}^=0-

Системы вида (18) можно преобразовать [5] к стандартному виду

Е

w=1

= , v =j =12, (19)

~(k) i(k)

где xj w и oj v — одномерные массивы, составленные из элементов

(k) (k) j,(fe) n(fe)

двумерных массивов xj ^n = Ymna},mn и соответственно, а D]^ — двумерный массив, составленный из элементов многомерного массива A(k)smn. Для решения бесконечных систем вида (19) применим метод редукции [8, 9], а решения конечных систем, полученных из (19) усечением, могут быть найдены методом квадратных корней [10].

В результате построен алгоритм нахождения приближенного аналитического решения задачи (1)-(9) в форме тригонометрических многочленов

N1 Ni-m

T (r, Z, tk) w ^ Ymn^l^mn COS (^mf) COS (v^Z) ;

m=0 n=0 N2 N2-m

T2 (r,Z,tk) ~ X] S

Ymna2,mn COS (^mr) COS (V2,n (Z - d)) ,

m=0 n=0

коэффициенты xjkmn = Ymnajkmn которых находят из конечных систем линейных алгебраических уравнений

Е

jС = С, v = 1,2....M,, j = 1,2.

•ш=1

Здесь М, = (Ж, + 1) (Ж, + 2) /2 , а значения N определяются согласно оценке Рунге [11].

Выбор шага по временной переменной осуществляется автоматически. Для обеспечения заданной точности решения используется правило двойного пересчета [10].

Результаты численных расчетов. Разработанный алгоритм использован для расчета температурного поля металлического цилиндра (слой 2) высотой к = 20 • 10-3 м, на нижнее основание которого нанесено поглощающее покрытие (слой 1) толщиной d = 0, 5 • 10-3 м [1]. Вычисления проводились при следующих значениях параметров задачи: г0 = 100 • 10-3 м, г, = 10 • 10-3 м, р1 = 3000 кг/м3, р2 = 7780 кг/м3, ТС = 300 К, д0 = 5 • 106 Вт/м2, а = 600 Вт/(м2• К), Л = 5 • 10-4 м2•К/Вт, а = 5,67 • 10-8 Вт/(м2•К4), е = 0,8, Д* = 2 с, * = 1,5 с.

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

Теплофизические свойства материалов

T, K 300 400 600 800 1000 1200 1400 1600 1800 2000

Аь Вт/(м-К) 32 28 20 16 7,5 6,4 5,6 5,4 5,2 5

c1, Дж/(кг-К) 860 875 943 1020 1086 1102 1140 1160 1170 1178

А2, Вт/(м-К) 48 47 41 37 32 23 21 20 — —

С2, Дж/(кг-К) 469 505 521 660 616 577 560 545 — —

На рис. 3 и 4 показана эволюция температуры в разных точках металлического слоя цилиндра. Периодическое внешнее тепловое воздействие приводит к возникновению колебаний температуры во всех точках цилиндра. Следует отметить, что, как и в линейных моделях [13, 14], период колебаний температуры совпадает с периодом воздействия источника теплоты. При удалении точек от источника в глубь материала наблюдается уменьшение амплитуды колебаний и сдвиг фазы колебаний. С момента времени £ = 90 с начинается процесс установления колебаний. Максимальный перепад температуры в точке по-

Т, К 1100

1000

900

800

700

600

500

400

300

1 : /л

А 2Л

\ 3

W i

10 12 14 16

t, С

Рис. 3. Эволюция температуры на начальном этапе разогрева (0 < t < 20 е) в

разных точках, лежащих на оси металлического цилиндра:

Н Н

1 — Т2(0, <1,1), 2 — т2(0,й +20,*), 3 — Т2(0 , 1 + —,1)

Рис. 4. Установившиеся колебания температуры в разных точках, лежащих на

оси металлического цилиндра:

Н Н

1 — Т2(0, <1,1), 2 — т2(0,1 +203 — Т2(0 , 1 + —

верхностного слоя металла для установившихся колебаний достигает величины 150 К.

СПИСОК ЛИТЕРАТУРЫ

1. Григорьянц А. Г., Шиганов И. Н., Мисюров А. И. Технологические процессы лазерной обработки. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2006.

2. Григорьянц А. Г. Основы лазерной обработки материалов. - М.: Машиностроение, 1989.

3. Углов А. А., С м у р о в И. Ю., Л а ш и н А. М., Гуськов А. Г. Моделирование теплофизических процессов импульсного лазерного воздействия на металлы. - М.: Наука, 1991.

4. К а р т а ш о в Э. М. Аналитические методы в теории теплопроводности твердых тел. - М.: Высш. шк., 2001.

5. Ч и г и р е в а О. Ю. Математическое моделирование процесса разогрева цилиндрической поверхности движущимся интенсивным источником тепла. // Инженерно-физический журнал. - 2006. - Т. 79, № 6. - С. 31-37.

6. М а л о в Ю. И., Мартинсон Л. К. Приближенные методы решения краевых задач. - М.: Изд-во МВТУ им. Н.Э. Баумана, 1989.

7. М а л о в Ю. И., М а р т и н с о н Л. К., Р о г о ж и н В. М. Математическое моделирование процессов тепломассопереноса при плазменном напылении // Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение. - 1994. - № 3. - С. 316.

8. К а н т о р о в и ч Л. В., Крылов В. И. Приближенные методы высшего анализа. -М.-Л.: Физматгиз, 1962.

9. Канторович Л. В., А к и л о в Г. П. Функциональный анализ. - М.: Наука, 1984.

10. А м о с о в А. А., Д у б и н с к и й Ю. А., К о п ч е н о в а Н. В. Вычислительные методы для инженеров. - М.: Высш. шк., 1994.

11. Чигирева О. Ю. Расчет оптимальной толщины слоя термоизоляции в многослойном цилиндрическом пакете // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2005. - № 1. - С. 94-101.

12. Ч и р к и н В. С. Теплофизические свойства материалов: Справочник. - М.: Физматгиз, 1959.

13. Лыков А. В. Теория теплопроводности. - М.: Высш. шк., 1967.

14. ТихоновА. Н., СамарскийА. А. Уравнения математической физики. - М.: Наука, 2004.

Статья поступила в редакцию 10.04.11

Ольга Юрьевна Чигирева родилась в 1979 г., окончила МГТУ им. Н.Э. Баумана в 2002 г. Канд. физ.-мат. наук, доцент кафедры "Математическое моделирование" МГТУ им. Н.Э. Баумана. Автор ряда научных работ в области математической физики и математического моделирования.

O.Yu. Chigiryova (b. 1979) graduated from the Bauman Moscow State Technical University in 2002. Ph. D. (Phys.-Math.), assoc. professor of "Mathematical Simulation" department of the Bauman Moscow State Technical University. Author of some publications in the field of mathematical physics and mathematical simulation.