асп. ЧирковА.Е., д.т.н., доц. Салов А.Г.

Самарский государственный технический университет

Численное моделирование нестационарного температурного поля пластины от движущегося нормально распределенного источника тепла

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

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

Рассмотрим процесс пайки двух пластин лазерным источником энергии.

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

На рисунке 1 приведено распределение плотности мощности теплового источника.

Рис. 1. Распределение плотности мощности по сечению источника тепла

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

(1)

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

(2)

В начальный момент времени температура тела равна T₀

(3)

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

С целью определения температурного поля, возникающего в процессе лазернойобработкив пластине при действии нормально-распределённого источника тепла, воспользуемся уравнением теплопроводности в виде [1]:

,

(4)

где: T – температура тела, t – время, a – коэффициент температуропроводности тела, V – скорость движения источника тепла, c_γ – объемная теплоемкость тела, x, y, z – координаты.

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

Распределение плотности мощности нормально распределенного источника тепла описывается функцией:

,

(5)

где: q(r) – распределение плотности мощности по радиусу сечения источника тепла; q_m – максимальная плотность мощности в центре источника тепла; k – коэффициент кривизны нормального распределения; r–расстояние от центра источника тепла.

С помощью функции (5) можно описать нормально распределенный источник тепла любого радиуса. В первом приближении границей источника тепла будем считать радиус, на котором плотность мощности будет составлять 20% от максимального значения (то есть ограничим кривую нормального распределения по уровню 0.2). Тогда необходимый для достижения требуемого радиуса коэффициент кривизны можно рассчитать по формуле:

(6)

Обычно считается известной средняя мощность источника тепла по его сечению P. Используя соотношение объема гауссоиды V_г и объема цилиндра V_ц в уравнении (7) можно вычислить максимальное значение плотности мощности в середине источника тепла.

(7)

С учетом вышеизложенного, для источника тепла радиуса r = 0,1 мм и средней мощности P = 0,4 Вт, получаем необходимый коэффициент кривизны k = 160,944 и максимальную плотность мощности q_m = 0,789.

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

(8)

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

По принципу наложения температура в момент времени tпродолжающегося действия нормально-кругового подвижного источникаqравна сумме температур dTот всех элементарных источников dQ(t’), выделившихся за время от t’=0 до t’=t действия источника на всем пути его перемещения.

Таким образом, решением уравнения теплопроводности для случая нагрева полубесконечного тела движущимся нормально распределенным источником тепла будет функция следующего вида [2]:

(9)

где:T(x,y,z,t) – температура в точке (x;y;z) относительно начала координат в момент времени t; t – время воздействия источника тепла; с – теплоемкость материала; γ – плотность материала; a – температуропроводность материала;

P – мощность источника тепла; t₀ – постоянная времени, V – скорость движения источника тепла по оси X; T₀ – начальное состояние теплового поля.

Проблема определения температуры тела в конкретной точке по соотношению (9) заключается в решении численным методом интеграла, который аналитического решения не имеет.

Для численного решения интеграла преобразуем его, записав в виде:

(10)

Наиболее распространённым методом численного решения таких функций является метод Симпсона.

Суть метода заключается в описании подынтегральной функции на заданном участке интерполяционным многочленом второй степени. Для повышения точности численного интегрирования, лучше воспользоваться модификацией формулы Симпсона – формулой Котеса,которая имеет следующий вид:

(11)

где: – величина шага; x_k = a + k·h – узлы интегрирования, границы элементарных отрезков, на которых применяется формула Симпсона; N – число разбиений.

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

Фактически подынтегральная функция на интегрируемом отрезке аппроксимируется Nпараболами. Таким образом, на гладких функциях типа (9) требуемая точность достигается при наименьшем числе разбиений.

Подставляя формулу (10) в формулу (11), получаем итоговое уравнение для численного расчета теплового поля на ЭВМ

(12)

где:T(x,y,z,t) – температура в точке (x;y;z) относительно начала координат в момент времени t; t – время; с – теплоемкость материала; γ – плотность материала; a – температуропроводность материала; P – мощность источника тепла; t₀ – постоянная времени; V – скорость движения источника тепла по оси X; T₀(x,y,z) – начальное состояние теплового поля в точке (x,y,z);x, y, z – координаты; N – число разбиений.

Полученное соотношение (12) позволяет достаточно быстро рассчитать распределение температуры по оси X (вдоль движения источника тепла).

На Рис. 2 приведено рассчитанное распределение температуры по оси X при следующих параметрах расчета:

теплоемкость - 1052 Дж/(кг·К); плотность - 2201 кг/м³; коэффициент температуропроводности -5,96·10^-7 м²/с; мощность источника тепла - 0,4 Вт;

скорость движения источника тепла - 1,667·10^-4 м/с; радиус источника тепла -0,1·10^-3 м; расчетное время - 1 с.

Рис. 2. Распределение температуры по оси X

Анализ температурного поля показывает, что максимум температур следует за источником тепла, принимая максимальное значение на некотором расстоянии от него. (см. рис 2).

На Рис. 3 показано распределение температур по оси Y. Максимум температур располагается симметрично оси движения источника тепла.

Рис. 3. Распределение температуры по оси Y

Аналогичное распределение наблюдается и по оси Z.

Рис. 4. Распределение температуры по оси Z

На рис. 5 изображен термический цикл в точке, центр источника тепла достигнет спустя 0,6 с с момента начала процесса сварки. Именно в этой точке наблюдается максимум температуры. На стадии нагрева заметен крутой изгиб кривой температуры, это связано с нормальным распределением источника тепла, и с его движением.

Рис. 5. Термический цикл в точке (0,0001;0;0)

Выводы:

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

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

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

Литература:

[1] Лыков А.В. Теория теплопроводности. Москва, изд-во «Высшая школа», 1967 г.

[2] Рыкалин Н.Н. Расчеты тепловых процессов при сварке. Москва, государственное научно-техническое издательство машиностроительных технологий, 1951 г.