Проблема учета зависимости коэффициента объемной теплоемкости от температуры при моделировании лазерно-дуговой наплавки Текст научной статьи по специальности «Физика»

5 МАТЕМАТИЧЕСКОЕ И КОМПЬЮТЕРНОЕ

МОДЕЛИРОВАНИЕ MATHEMATICAL AND COMPUTER SIMULATION

УДК 519.633.6:517.958:536.2

ПРОБЛЕМА УЧЕТА ЗАВИСИМОСТИ КОЭФФИЦИЕНТА ОБЪЕМНОЙ ТЕПЛОЕМКОСТИ ОТ ТЕМПЕРАТУРЫ ПРИ МОДЕЛИРОВАНИИ ЛАЗЕРНО-ДУГОВОЙ НАПЛАВКИ С.А. Лукьяненко3, В.А. Третьяк3 а Национальный технический университет Украины «Киевский политехнический институт», Киев, Украина, [email protected]

На примере задачи моделирования динамики температурного поля при лазерно-дуговой наплавке композитных материалов исследуется проблема учета явлений плавления и кристаллизации сталей и сплавов. Математическая модель принадлежит классу задач Стефана и описывается нестационарным нелинейным трехмерным уравнением теплопроводности с граничными условиями. Задача решается конечно-разностным методом сквозного счета. Нелинейность задачи обусловлена зависимостью коэффициента объемной теплоемкости от температуры, которая имеет всплеск в интервале плавления-кристаллизации. Для учета этой зависимости предлагается метод линейно-экспоненциальной аппроксимации. За счет стремительности нагрева в некоторых узлах разностной сетки возникает проблема перехода температуры из диапазона твердой фазы в диапазон жидкой фазы за один шаг интегрирования по времени. Для решения этой проблемы предлагается алгоритмический метод учета всплеска зависимости коэффициента объемной теплоемкости от температуры в интервале плавления-кристаллизации, который основан на идее вложенных сеток. Предложенный метод увеличивает алгоритмическую сложность расчета, однако дает гарантию учета физики процесса плавления-кристаллизации. Приводится сравнение результатов расчетов, полученных при использовании кусочно-линейной и линейно-экспоненциальной аппроксимации зависимости коэффициента объемной теплоемкости от температуры. Показано, что в случае использования линейно-экспоненциальной аппроксимации процесс решения на адаптивной сетке ускоряется за счет уменьшения количества временных шагов интегрирования. Ключевые слова: уравнение теплопроводности, коэффициент объемной теплоемкости, задача Стефана, метод сквозного счета.

TEMPERATURE DEPENDENCE CONSIDERATION ISSUE FOR COEFFICIENT OF VOLUMETRIC HEAT CAPACITY IN SIMULATION OF LASER-ARC PAD

WELD PROCESS

S. Luk'yanenkob, V. Tret'yakb

b National Technical University of Ukraine «Kyiv Polytechnic Institute», Kyiv, Ukraine, [email protected] The paper deals with research of steels and alloys melting and crystallization on the example of temperature field dynamics simulation for the laser-arc pad weld of composite materials. Mathematical model is one of the Stefan problems for nonlinear transient three-dimensional Fourier equation of heat conduction with boundary conditions. The problem is solved by the finite-difference method with the pass-through calculation scheme. Nonlinearity of the problem is caused by the temperature dependency of the volumetric heat capacity coefficient that has a peak in the interval of phase transition (melting - crystallization). The linear-exponential approximation method is offered to consider this dependence. At some nodes of the discrete mesh the temperature changes from solid phase range to liquid phase at one time step on account of high heat velocity. The algorithm of temperature dependence consideration for volumetric heat capacity peak in the interval of melting - crystallization is suggested based on the nested mesh idea. Though the computational complexity is increased by suggested algorithm it provides more accurate solution and ensures account of phase transition phenomenon. Calculation results for linear and linear-exponential approximation of temperature dependence for volumetric heat capacity are compared. It is shown that in case of linear-exponential approximation the solution process on dynamic adaptive mesh is speeded up due to decreasing of time steps number. Keywords: heat capacity equation, volumetric heat capacity coefficient, Stefan problem, pass-through calculation scheme.

Введение

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

ПРОБЛЕМА УЧЕТА ЗАВИСИМОСТИ КОЭФФИЦИЕНТА ОБЪЕМНОЙ ТЕПЛОЕМКОСТИ...

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

Технологический процесс наплавки предполагает, что порошок наносится на поверхность потоком транспортирующего газа, при этом для увеличения эффективности он может подогреваться электрической дугой с распределением плотности мощности теплового потока qa, создавая плазму. Лазерный луч фокусируется на поверхности детали и завершает процесс плавления, пятно фокусировки луча может находиться перед областью действия дуги (на поверхности подложки) или после (на поверхности слоя наплавления), также они могут пересекаться. Более детально технологический процесс описан в [4]. Он включает множество сложных явлений [5, 6], все их невозможно учесть при моделировании.

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

Математическая постановка задачи

В работе зависимость теплофизических свойств подложки, а также коэффициента теплопроводности покрытия от температуры не учитываются. Считается, что слой напыления в виде идеально спрессованного композитного порошка моментально появляется под электрической дугой, т.е. область определения (ОО) задачи моделирования меняется во времени в соответствии с передвижением дуги (рис. 1). Положим размеры ОО Lx, Ly, Lz по соответствующим координатам, h - высота покрытия, Га - область нагрева дугой, Г, - пятно фокусировки лазера, Ц1 - неподвижный цилиндр, который образовался под дугой в первый момент времени, Ц2 - цилиндр, соответствующий текущему положению дуги, который перемещается со скоростью V вдоль оси Оу , оставляя за собой поверхностный слой, П - параллелепипедная часть покрытия. Таким образом, слой покрытия имеет переменную форму, образованную слиянием Ц1, П и Ц2.

х

Ча _У V

Рис. 1. Схематическое изображение переменной ОО

Распределение плотности мощности (РПМ) лазера по пятну фокусировки может быть различной конфигурации [7], чаще всего для моделирования РПМ лазера применяют функцию Гаусса [8]. В настоящей работе используются РПМ луча, описанные в [9]: гауссово, постоянное, кусочно-линейное и кусочно-постоянное. Хотя регулировать РПМ дуги так же, как и лазера, не представляется возможным, в предложенной модели считается, что конфигурация РПМ дуги может иметь такие же виды, как и луч. В местах пересечения луча и дуги РПМ комбинированного нагрева описывается предложенными в [10] формулами.

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

d

- (U • C(U )) = X dt

(d2U dU dU ^

dx2 dy2

dz2

Ч=0 = Uc-

X

dU ( x, y, z0, t )

dz

+ q(x, y, z, t) = 0, V(x, y) еГ = Га иГ,

X

dU dn

+ a[U - Uc ] = 0, V(x, y) ёГ = Г„ и Г,

(1)

где и - температура; С(Ц) - коэффициент объемной теплоемкости (КОТ); X - коэффициент теплопроводности; ис - температура окружающей среды; г0 - уровень верхней грани детали, соответствующий либо подложке (х0=И), либо покрытию (г0=0); ч(х, у, г, /) - функция РПМ нагревающих элементов; п - нормаль к поверхности детали; а - коэффициент теплоотдачи; О - граница области определения.

h

L

Z

z

+

ü

При этом X и С(П) зависят от пространственных координат. Если координаты соответствуют подложке, берутся значения коэффициентов материала подложки, которые считаются константами. Для покрытия берутся эффективные величины Хе, и Се, (Ц), соответствующие композитному материалу. Эффективный коэффициент теплопроводности, как и ранее [11], рассчитывается по формуле Максвелла:

1 + 2v/

f 1 Wf у Л 1 + 2 m

1

\ff = Xm -> V' If -Г"Y = COnSt^

У f

f J

1 _Vf

1 + 2

X /

где Хт, X, - коэффициенты теплопроводности матрицы и включений композита; Vf - объемная концентрация включений, не превышает 0,3. Определение значений КОТ рассмотрим детальнее.

Зависимость коэффициента объемной теплоемкости от температуры

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

= (ф)е# = ся а-V,)+,

где Ст и С, - эффективные значения КОТ матрицы и включений, которые рассчитываются согласно скрытой теплоте плавления:

Эф Эф f

Cm = Cm (U) _ Lm и Cf = Cf (U) _ Lf-/

m m\ / m QU J JV ' J QU

где Lm и Ь, - скрытая теплота плавления; фт и ф, - объемная концентрация жидкой фазы в материале матрицы покрытия и включений соответственно. Задача усложняется тем, что достоверной информации о значениях фт и ф, нет [3], фт и ф, зависят от состава материала, температуры и скорости нагрева или охлаждения. Например, в [12] предложена формула для мартенситных превращений. Однако удобнее выполнять аппроксимацию непосредственно зависимости КОТ от температуры. Ранее [11] зависимость КОТ от температуры определялась с помощью кусочно-линейной аппроксимации (КЛА). Очевидно [13], изменение КОТ во время плавления и кристаллизации имеет более плавный характер. В [10] предложено выполнять аппроксимацию квадратичным полиномом. В данной работе предлагается еще один способ аппроксимации КОТ - линейно-экспоненциальный (ЛЭА) (рис. 2), согласно которому КОТ вычисляется по формуле

с(Ц)=Ье-т -Ц™)2 + т+г, (2)

где Ь соответствует скрытой теплоте плавления материала; р определяет ширину интервала плавления-кристаллизации; Цпл - температура плавления, для сталей и сплавов считается медианой [Ц]; k -угловой коэффициент линейной составляющей функции; 2 - значение аппроксимирующей функции КОТ при температуре 0°С.

Стремительность возрастания С(Ц) при достижении Ц регулируется параметром у , который равен значению первого слагаемого формулы (2) при Ц = Ц :

Ье-р(Ц -Ц-)2 =у. Получаем формулу для расчета р :

' Ь

1п

Р=—^ •

(и S - Цпл)2

Проблема учета фазового перехода на уровне алгоритмического метода

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

ПРОБЛЕМА УЧЕТА ЗАВИСИМОСТИ КОЭФФИЦИЕНТА ОБЪЕМНОЙ ТЕПЛОЕМКОСТИ.

28 000 26 000 24 000 22 000 20 000 18 000 16 000 14 000 12 000 10 000 8 000 6 000 4 000

(U ; С(ик ))

200 400 600 800 1000 1200 1400 1600

U

sf

U

Lf

U,

s.

UL

Рис. 2. Линейно-экспоненциальная аппроксимация КОТ: штрихпунктир - материал матрицы покрытия (сталь 9Х18), штриховой пунктир - материал включений (бронза ОЦС 3-8-6), сплошная - график эффективного КОТ

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

В работе предлагается по окончании каждого шага интегрирования (перед формированием новой сетки) совершать поиск узлов, где имело место изменение фазы. Затем такие узлы объединяются в подобласти, для которых производится перерасчет с шагом в 2 раза меньше. На границах подобластей необходимо определить граничные условия. Поскольку значения в граничных узлах подобластей уже найдены на ^м и (k +1) -м шагах, можно использовать граничные условия первого рода. На дробном ^ +1/2) -м шаге значения в граничных узлах неизвестны, поэтому необходимо найти их при помощи интерполяции (используется метод Лагранжа 2-го порядка). Перерасчет предлагается выполнять рекурсивно, пока выделенные подобласти будут содержать узлы. Подобласти выделяются таким образом, что отстоящие друг от друга на 3 и менее пространственных шага узлы с фазовым переходом помещаются в одну подобласть. Более того, размерность подобласти обязательно должна быть не менее трех, что достигается путем добавления узлов по краю подобласти.

Исследование эффективности предложенных методов

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

Время, с Результат без пересчета U1, °С Результат с пересчетом u2, °С |u 1 - U2|,°C |u1 - u2 1 1 2,% U1

0,300 2300 2027 273 11,870

0,400 2410 2208 202 8,382

0,454 2330 2181 149 6,395

0,500 2311 2119 192 8,308

0,600 2407 2240 167 6,938

0,700 2325 2233 92 3,957

0,810 2430 2359 71 2,922

Средние значения 163,714 6,967

Таблица 1. Анализ метода учета фазового перехода путем пересчета по подобластям

Критерий сравнения КЛЛ (F1) ЛЭA( F 2) F1 - F 2 100 • (F1 - F2) I F1,%

Время расчета на равномерной сетке, с 153 187 -34 -22,2

Время расчета адаптивным методом, с 121 98 23 19

Количество временных шагов, шт. 80 62 18 22,5

Количество отказов, шт. 45 17 28 62

Таблица 2. Анализ влияния метода аппроксимации КОТ на характеристики расчета

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

Также был проведен сравнительный анализ результатов моделирования при условии использования КЛА (F1) и ЛЭА (F2) аппроксимации КОТ (табл. 2).

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

Заключение

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

Литература

1. Лазерш технологи та комп'ютерне моделювання / Шд ред. Л.Ф. Головка, С.О. Лук'яненка. Кшв: Вктка, 2009. 296 с.

2. LIA Handbook of Laser Materials Processing / Eds J. F. Ready, D. F. Farson. Orlando, USA: Laser Institute of America, Magnolia publishing, Inc., 2001. 715 p.

3. Князева А.Г., Крюкова О.Н., Букрина Н.В., Сорокова С.Н. Проблемы моделирования технологических процессов поверхностной обработки материалов и нанесения покрытий с использованием высокоэнергетических источников // Известия Томского политехнического университета. 2010. Т. 316. № 2. С. 93-101.

4. Третьяк В.А. Усовершенствование адаптивного метода построения сеток для задач теплопроводности с нестационарным источником энергии // Математичне та комп'ютерне моделювання. Сер. Техшчш науки. 2012. В. 7. С. 197-206.

5. Picasso M., Marsden C.F., Wagniere J.D., Frenk A., Rappaz M. A simple but realistic model for laser cladding // Metallurgical and Materials Transaction B. 1994. V. 25. N 2. P. 281-291.

6. Масленников А.В., Левин Ю.Ю., Ерофеев В.А. Моделирование нестационарного формирования шва при импульсной лазерной ND:YAG-сварке // Известия Тульского государственного университета. Технические науки. 2010. В. 4. Ч. 1. С. 255-265.

7. Гончарский А.А., Туницкий В.В. Об обратной задаче синтеза оптических элементов для лазерного излучения // Вычислительные методы и программирование. М.: НИВЦ МГУ, 2006. Т. 7. С. 138-162.

8. Qi H., Mazumder J., Ki H. Numerical simulation of heat transfer and fluid flow in coaxial laser cladding process for direct metal deposition // J. Appl. Phys. 2006. V. 100. N 2. P. 024903-1-024903-11 [Электронный ресурс]. Режим доступа: http://dx.doi.org/10.1063/L2209807, свободный. Яз. англ. (дата обращения 05.11.2013).

9. Головко Л.Ф., Лук'яненко С.О., Смаковський Д.С., Михайлова 1.Ю., Агеенко В.А. Моделювання температурного поля при змщненш матерiалiв лазерним випромiнюванням // Моделювання та шформацшш технологи. 2008. № 45. С. 28-35.

МЕТОДЫ МОДЕЛИРОВАНИЯ ТЕМПЕРАТУРНОГО ПОЛЯ.

10. Amara E.H., Hamadi F., Achab L., Boumia O. Numerical modelling of the laser cladding process using a dynamic mesh approach // Journal of Achievements in Materials and Manufacturing Engineering. 2006. V. 15. N 1-2. P. 100-106.

11. Головко Л.Ф., Лукьяненко С.А., Смаковский Д.С., Агеенко В.А., Михайлова И.Ю. Моделирование адаптивным сеточным методом температурного поля при лазерной наплавке порошковых материалов // Электронное моделирование. 2009. Т. 31. № 1. С. 21-32.

12. Mamat M., Tofany N., Kartono A. Numerical analysis of heat conduction and phase transformation in laser transformation hardening: influences of heating duration and laser beam intensity // Applied Mathematical Sciences. 2010. V. 4. N 61-64. P. 3019-3033.

13. Мышковец В.Н., Максименко А.В., Баевич Г.А., Грищенко В.В. Термические циклы в зоне рекристаллизации при импульсной лазерной наплавке среднелегированных высокопрочных сталей // Известия Гомельского государственного университета им. Ф. Скорины. 2011. № 6 (69). С. 105-109.

Лукьяненко Святослав Алексеевич - доктор технических наук, профессор, зав. кафедрой, Национальный технический университет Украины «Киевский политехнический институт», Киев, Украина, [email protected]

Третьяк Валерия Анатольевна - аспирант, Национальный технический университет Украины «Киевский

политехнический институт», Киев, Украина, [email protected]

SvyatoslavLuk'yanenko - D.Sc., Professor, Head of Department, National Technical University of Ukraine

«Kyiv Polytechnic Institute», Kyiv, Ukraine, [email protected]

Valery Tret'yak - postgraduate, National Technical University of Ukraine «Kyiv Polytechnic

Institute», Kyiv, Ukraine, [email protected]

УДК 519.63

МЕТОДЫ МОДЕЛИРОВАНИЯ ТЕМПЕРАТУРНОГО ПОЛЯ ПРИ БЕСКОНТАКТНОЙ ЛАЗЕРНОЙ ДЕФОРМАЦИИ ПЛАСТИНЫ С.А. Лукьяненко3, И.Ю. Михайлова3

а Национальный технический университет Украины «Киевский политехнический институт», Киев, Украина, [email protected]

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

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

METHODS OF TEMPERATURE FIELD MODELING FOR CONTACTLESS LASER

DEFORMATION OF A PLATE S. Luk'yanenkob, I. Mykhailovab

b National Technical University of Ukraine «Kyiv Polytechnic Institute», Kyiv, Ukraine, [email protected] Contactless laser deformation is a process of changing a geometrical form of an object by its heating and cooling without using mechanical forces. One of the factors, influencing deformation, is heating, which creates temperature field and overpatching speed of temperature field.

The article deals with a computer modeling result for temperature field of a plate irradiated by a moving laser beam. Comparison of results, obtained by using two mathematical models and experiments, is conducted. Such parameters as density, thermal capacity and thermal conductivity are constant values in the first model, and they are linearly dependent on temperature in the second one. Linear approximation of table values by least-squares method is used to define this dependence. Both models have the temperature field computed by 3D non-stationary heat equation: it is linear function in first model and quasilinear function in the second. To solve the equation a six- step implicit finite difference coordinate-wise splitting scheme is