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

УДК 536.2

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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА РАЗОГРЕВА ДВУХСЛОЙНОГО ЦИЛИНДРА ДВИЖУЩИМСЯ КОЛЬЦЕВЫМ ИСТОЧНИКОМ ТЕПЛОТЫ

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

E-mail: k_fn12@bmstu.ru

Ключевые слова: двухслойный цилиндр, движущийся источник теплоты,

процесс разогрева, температурное поле, термическое сопротивление.

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

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

Постановка задачи и математическая модель процесса. Рассмотрим процесс разогрева металлического цилиндра радиусом r0 и высотой h, на боковую поверхность которого нанесен слой поглощающего покрытия [2] толщиной 8 (рис. 1). Разогрев цилиндра проводится кольцевым источником теплоты с плотностью теплового потока q0, совершающим импульсное движение по боковой поверхности цилиндра от нижнего основания к верхнему по закону

s (t) = (i - 1) (г - 1)Д t<t < iAt, г = 1, 2,..., I,

Рис. 1. Осевое сечение двухслойного цилиндра:

1 — металлический цилиндр; 2 — поглощающее покрытие

где £ — время; Д£ — продолжительность одного импульса; I — число импульсных шагов; I — ширина теплового контакта; в — безразмерный параметр. При 0 < в < 1 положения теплового контакта на г-м и (г + 1)-м импульсных шагах частично перекрываются. В момент времени = Д£ • I, когда тепловой контакт достигает верхнего основания цилиндра, процесс нагревания прекращается (д0 = 0 при £ > и начинается стадия остывания.

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

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

Цель исследований — определение глубины зоны локального прогрева металлического цилиндра до заданной предельной температуры Тп в зависимости от значений параметров (¡0 и в, определяющих интенсивность воздействия источника теплоты.

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

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

P1C1 (Tl) ИГ = ГдГ lAl (Tl) ^^

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

+ dZ(Ai (Ti) t

(1)

( dT2 1 д ( . дТ2 p2C2 (T2) -m = rör lA2 (T) rö7

t > 0, r0 <r <r0 + 8, 0 < z < h;

+ ÖZ(A2 (T2) £

(2)

Ti (r, z, 0) = Tc, 0 < r < ro, 0 < z < h; T2 (r, z, 0) = Tc, ro < r < ro + S, 0 < z < h;

Л2 (T2)

dTo

dr dTi'

qi (z,t), 0 <t < U;

r=ro+S ^2 (z,t) , t>t*;

dz

z=0, z=h

= 0, t> 0, 0 < r < ro;

dT2

dz

z=o,

z=h

= 0, t> 0, r0 < r < r0 + S;

(3)

(4)

(5)

(6) (7)

Ai (Ti)

dTi

dr

r=r 0

= ъ(Т2 (ro,z,t) - Ti (ro,z,t)) R

= A2 (T2)

dTo

dr

t> 0, 0 < z < h. (8)

r=r0

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

В граничном условии (5)

Г до, г е [8 (г) ,в (г) + I]; (г,г) = < а (Тс - Тщ (г,г)) +

( +ае (ТС - Т4 (г,г)), г е [0,к \ [в (г) ,8 (г) + I];

д2 (г, г) = а (Тс - Тщ (г, г)) + ае (ТС - Т^ (г, г)) , г е [0, к],

где Тщ (г, г) = Т2 (г0 + 5, г, г); а — постоянная Стефана-Больцмана; е — степень черноты излучающей поверхности.

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

Введем функции

С3 (Т,г) = р3гс3 (Т), Л3 (Т3,г) = гЛ3 (Т3), з = 1, 2, и запишем уравнения (1), (2) следующим образом:

dT

Ci (Ti, r) = div (Ai (Ti, r) gradTi),

t> 0, 0 < r <ro, 0 < z < h;

(9)

дт

С2 (Т2, г) = а1у (Л2 (Т2, г) gradТ2), (10)

£ > 0, г0 < г < г0 + 5, 0 < г < Л.

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

дТ3 j (r,z) - I(k-1) (r,z)

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

т

где Т(к) (г, г) — приближенные значения функций Т (г, г, ¿) при £ = £к.

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

С]к) (г,г) = С (т('-1) (г, г), г) , Л(к) (г, г) = Л^ (т(к-1) (г, г), г) , ; = 1, 2; 8(к) = 8 (£к);

#) ^ _ J q(k) (z) , 0 < tk < t*;

q2k) (z) , t > t*;

QW (z) =

Г (г0 + 5) ф), г е |>(к),8(к) + /] ; д(к) (г) = Г (г0 + 5) (а + 4теТ£) (тс - 1#-1) (г)) ,

[ г е [0,Л] \ |>(к),8(к) + /] ; ф(к) (г) = (г0 + 5) (а + 4<теТ£) (тс - Т^-1) (г)) , г е [0, Л];

д0к) (г) = Я (Т2(к-1) (г0,г) - Т(к-1) (г0,г)) .

Далее запишем дифференциально-разностный аналог начально-краевой задачи (1)-(8) в виде итерационной схемы (к = 1, 2,...) решения двух краевых задач для линейных эллиптических уравнений с переменными коэффициентами:

- Шу (Л1к) (г, г) grad Т(к) (г, г)) + 1 С((к) (г, г) Т(к) (г, г) =

= 1 С{к) (г, г) Т(к-1) (г, г), 0 < г < г0, 0 < г < Л; (11) т11

дТ{к)

0Г

= 0, л1к) (r,z) дТ(к)

r=0 0Г

= Q0k) (z), 0 < z < h; (12)

r=ro

öTi

(k)

dz

z=0, z=h

= 0, 0 < r < ro;

(13)

- div (A2k) (r, z) grad T2(k) (r, z)) + 1 C2k) (r, z) T2(k) (r, z) =

= 1 C2k) (r, z) T2(k-1) (r, z), ro < r < ro + 0 < z < h; (14)

Л2к) (r, z)

(k)

dr

(k)

(k) _\ dT2

Л2'') (r, z)

dT

dr

(k)

dz

z=0, z=h

= Qok) (z),

r=ro

= QWW) (z), 0 < z < h;

r=ro +5

= 0, r0 ^ r ^ r0 + 5.

(15)

(16)

Следует отметить, что на каждом временном слое £ = задачи (11)-(13) и (14)-(16) решаются независимо.

На первом шаге итерации, согласно начальным условиям (3) и (4), Т-0 (г, = Тс. На к-м шаге итерации функции Т^ (г, будем искать в форме разложения в двойные тригонометрические ряды Фурье

Tik) (r,z) = ^

ж ж

mnaj,jmnX'j,mn (r, z) , j 1 2

где Tmn TmTn, Tm

m=0 n=0

0,5, m = 0,

1, m > 0,

по полным и ортогональным

ж

m,n=0'

[4] в областях П, j = 1, 2, системам функций {Xj mn (r, z)}

Xi,mn (r, z) = cos (^I,mr) cos (vnz), П = (0, ro) x (0, h);

X2,mn (r, z) = COS (^2,m (r - Го)) COS (vnz) , П2 = (ro, Го + 0) X (0, h) . mn mn nn

Здесь = , ^2.m = f , = "T"•

r0 oh

Для нахождения коэффициентов a(kmn этих разложений умножим уравнения (11) и (14) на функции X1ps (r, z) и X2 ps (r, z) соответственно, а затем проинтегрируем полученные соотношения по областям П1 и П2. Используя граничные условия (12), (13) и (15), (16), приходим к бесконечным системам линейных алгебраических уравнений относительно искомых коэффициентов а

(k):

ж ж

VVA(k) Y a(k)

/ j / j j,psmn /mnU,ji,mn m=0 n=0

= b'

(k) ,

(17)

p = 0,1, 2,

s = 0,1, 2,..., j = 1, 2,

где

A(k) - т (и и + ) if(k) -Р(k)

Aj,psmn — 1 (u'j,mU'j,p + vnvs) ^Sj,(|m-p|,|n-s|) pj,(m+p,n+s)J + + т (uj,mUj,p — VnVü) ^Pj,()m-p|,n+s) — Pj, (m+p,|n—s|)) +

ъjXIm-pln+s) -j + (k) + (k) + (k) + (k) ; + П j,(lm-pl,ln-s\) + П j,(m+p,\n—s\) + П j,(lm—pl,n+s) + П j,(m+p,n+s);

<X <X /

b(k) — f (k) a(k-i) n(k) + bj,ps jj + 'mnaj,mn ylj^m-pHn-siy

m=0 n=0 ^

. (k) + (k) + (k) + n3,(m+p\n-s\) + Vj

f(k) — - (-I)p f(k) — 8-T ((-i)p4k) - ;

— ^1)p ^ f (k) = ^L

ro ( 1) ' f2 = 8 <p s) и nf) s) — коэффициенты Фурье разложений функций C(k (r, z) и Л<k (r, z) соответственно в двойные тригонометрические ряды по системам функций {XjiPS (r, z)}™s=0; ipS^ и — коэффициенты Фурье разложений функций Q0k) (z) и qWW (z) соответственно в тригонометрические ряды по системе функций {cos (vsz)}!;;=0.

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

Е D<&j = j), v = 1, 2,..., j = 1, 2, (18)

w=l ~(k) l(k)

где xjW и bjV — одномерные массивы, составленные из элементов

(k) (k) i(k) k) массивов x<< mn = YmnUjmn и b<< ps соответственно, а D<< vw — двумерный

(k)

массив, составленный из элементов многомерного массива AjJsmn.

Решения бесконечных систем вида (18) могут быть найдены методом редукции [7, 8].

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

правых частях уравнений (9), (10) являются симметрическими поло-

(k) (k)

жительно определенными, имеет место свойство Ajpsmn = Ajmnps, то

матрицы Djk), j = 1, 2, — симметрические положительно определенные. Это позволяет для решения конечных систем, полученных из (18) усечением, применить метод Холецкого [9].

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

N1 N-i-rn

Ti (r, z, tk) w E E Ymnd\lnn cos (vi>mr) cos (vnz);

r =0 n=0

N2 N2-m

T2 (r, Z, tfc) ~ ^ ^ Tmn^iL COS (r - Го)) COS (vnz) ,

m=0 n=0

коэффициенты x(kmn = T^j^n которых находят из конечных систем

Mj

Е ( xjkW = (J, V = 1, 2,...Mj, j = 1, 2.

w=1

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

При проведении численных расчетов выбор шага по временной переменной осуществлялся на основе сравнения в фиксированные моменты времени t распределений температуры на поверхности r = r0 при двух значениях шага (т и т/2):

\\Т- (ro,z,i) - Тт/2 (ro,z,i) || ^ ||Tr/2 (ro,Z,t)\\ ^

где е0 — заданное значение погрешности.

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

Расчеты проводились при следующих значениях параметров задачи: pi = 7780 кг/м3; р2 = 3000 кг/м3, Ta = 300K; Тп = 1500K; а = 50Вт/(м2-К); R = 10-4 м2-К/Вт; а = 5,6740-8 Вт/(м2^К4); е = 0,8; r0 = 24-10-3 м; S = 0,240-3 м; h = 12040-3 м; l = 1040-3 м; At = 0,5 c.

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

Таблица 1

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

Температура T, К

Свойства 300 400 600 800 1000 1200 1400 1600

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

еьДж/(кг-К) 469 505 521 660 616 577 560 545

Таблица 2

Теплофизические свойства материала поглощающего покрытия

Температура T, К

Свойства 300 400 600 800 1000 1200 1400 1800 2000

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

С2, ДЖ(кг-К) 860 875 943 1020 1086 1102 1140 1170 1178

Варьируя параметры задачи q0 и в в интервалах значений q0 = = (1... 5) • 107 Вт/м2 и 0 < в < 1, выбираем q0 = 1,5-107 Вт/м2 и

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

На рис. 2 и 3 представлены изотермы температурного поля (Т, К) в осевом сечении металлического цилиндра в различные моменты времени.

Как показали результаты расчетов, глубина зоны локального прогрева металлического цилиндра до температуры 1200 ... 1500 К при значениях параметров д0 = 1,5407 Вт/м2 и в = 0,5 составляет приблизительно 0,7мм, а ширина зоны прогрева не более 7мм.

Очевидно, что того же порядка значения глубины и ширины зоны локального прогрева могут достигаться и при иных значениях параметров, определяющих интенсивность воздействия теплового источника, например при д0 = 2,2407 Вт/м2 и в = 0,75.

60

40

20

310 ________

320 ... < \900

3JC к ( ™ i SQO^

400f Jllu \ 550 \ \ 60Q 650

О 5 10 15 20 г, мм

Рис.2. Изотермы температурного поля в осевом сечении металлического цилиндра в момент времени Ь = 5,5 c

Рис.3. Изотермы температурного поля в осевом сечении металлического цилиндра в момент времени Ь = 11,5 c

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

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

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

2. ГригорьянцА. Г., Сафонов А. Н. Основы лазерного термоупрочнения сплавов. - М.: Высш. шк., 1988.

3. Зарубин В. С. Инженерные методы решения задач теплопроводности. - М.: Энергоатомиздат, 1983.

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

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

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

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

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

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

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

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

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

Ольга Юрьевна Чигирева родилась в 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.