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

У статті розглянуто результат комп’ютерного моделювання температурного поля пластини під дією лазерного променя з урахуванням залежності густини, питомої теплоємності та теплопровідності матеріалу від температури. В основу моделі покладено тривимірне нестаціонарне квазілінійне рівняння теплопровідності, що розв’язується за допомогою методу покоординатного розщеплення з використанням адаптивної сітки

Ключові слова: метод покоординатного розщеплення, адаптивний метод побудови сітки, метод балансу

□--------------------------------------□

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

Ключевые слова: метод покоординатного расщепления, адаптивный метод построения сетки, метод баланса -------------------□ □-------------------------

УДК 519.63

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

С. А. Лукьяненко

Доктор технічних наук, професор, завідуючий кафедри* E-mail: lukian@aprodos.kpi.ua И. Ю. Михайлова

Аспірант* E-mail: imikh@aprodos.kpi.ua *Кафедра автоматизації проектування енергетичних процесів і систем Національний технічний університет України «Київський політехнічний інститут» вул. Політехнічна, 6, м. Київ, Україна, 03056

1. Введение

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

В работе рассматриваются модели распространения тепла в результате локального нагрева металли-

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

2. Постановка задачи

Исследуемый объект - металлическая пластина с геометрическими размерами Lx , Ly , Lz (рис. 1). Будем считать, что пластина неподвижна. Все ее поверхности находятся в процессе теплообмена с окружающей средой, температура которой ис. На ее верхнюю грань воздействует луч лазера с плотностью мощности излучения q(x,y,t) , который движется со скоростью V(t) параллельно оси ординат в течение времени Тк.

Рис. 1. Схематическое изображение расположения луча лазера относительно пластины

© С. А. Лукьяненко, И. Ю. Михайлова, 2013

Необходимо найти распределение температур в металлической пластине, на которую воздействует луч лазера.

3. Математическая модель

Рассмотрим математическую модель процесса нагрева детали лазерным излучением. В ней будем учитывать зависимость физических параметров металла от температуры. Эти зависимости примем линейными, полученными в результате среднеквадратичной аппроксимации экспериментальных табличных данных. Процесс распределения температуры будем описывать с помощью нестационарного квазилинейного трехмерного уравнения теплопроводности [4, 5]:

с(и)р(и)

Э

Эу

ОД

Эи(х, у, z,t) дї “

дЩх,у^,ї)

_Э_

Эх

ад

дЩх,у^,ї)

Эу

Э

ад

Эх

дЩх,у^,ї)

дz

(1)

где с(и) - удельная теплоемкость, р(и) - плотность, Х(и) - коэффициент теплопроводности материала,

х е[0;Е ], у ф'Е ], z е[0Е ], * Ф'Л ].

Краевые условия вне зоны действия лазера моделируют теплообмен с окружающей средой по закону Ньютона (с учетом закона Фурье) [6]:

х(и)

Эи(х, у, z,t) Эп

+ а[ис -Щх,у^,ї)] = 0,

(2)

. .... дЩх,у,0,ї) . _

Х(и)—У 7 + q(x,y,t) = 0. Эz

(3)

Начальная температура пластины равняется температуре окружающей среды:

Щх,у^,0) = ис.

(4)

ной матрицей. Для формирования систем уравнений на каждом шаге метода покоординатного расщепления будем использовать метод баланса [8].

Пусть приближенное решение найдено на (к -1) -ом слое на неравномерной сетке в узлах (х^у^т), где і = 0,...,^, ] = 0,...,п2, т = 0,...,п3, п1 - количество узлов по оси абсцисс, п2 - количество узлов по оси ординат, п3 - количество узлов по оси аппликат.

Обозначим величину шагов сетки в направлении оси абсцисс через Ь1 і = хі - хі-1, в направлении оси ординат - через = уі - уі_1, в направлении оси аппликат -

через Ь3т = zm -zm_1. Обозначим также среднее арифметическое двух соседних шагов для узла (х^у^т) в каждом из координатных направлений через Ь1с, Ь2с, — с; шаг по времени через т = тк = їк - їк-1 = їк+1 - їк ; значение приближенного решения в точке (xi,Уj,Zm,tk) -через икт ; точного решения - через и.т .

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

Используя неявную линейную схему метода баланса, получим уравнения первых двух этапов метода расщепления:

к-1

Чт

с (ик-1 )р(икт1) =-І- х

Т — 3,с

к-2 к-2

и.. 3, - и.. 3

і]гп+1_____і|т _____

где п - нормаль к поверхности, а - коэффициент теплоотдачи.

Краевое условие в зоне воздействия лазерного луча

[4]:

к-1 к-1

и.. ,+ и..

і]т+1 і]т

-X

( к-1 к-1 Л и.. + и.. „ і]т і] т-1 к-2 к-2 Л и.. 3 - и.. 3 і]т і]т-1

2 V ) — 3,т

к-2^ 3

V

і]т ( /

/

, \ к-1 к-2 2Л„ 3. "

V

/

к-2 к-2

и.. 3 + и.. 3

і]+1т і.т

к-1 к-1

и.. 3 -и.. 3

і]+1т________і.т

2]+1

4. Методы решения

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

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

----- /

к-2 к-2

и.. 3 + и.. 3

і]т і]-1т

к-1 к-1

и„ 3 -и„ 3

і]т і]-1т

Аналогично записываются уравнения для остальных четырех этапов.

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

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

т

2,с

Е

u

v

ijm / /

/

V

U-- 3 — uk 1

ijm ijm

/

— X

k—3 k—3

u.. 3 + u.. 3

ijm+1 ijm

k—3 k—з

u.. 3 + u.. 3

ijm ijm—1

3,c

k—3 k—3

u.. 31 — u.. 3

ijm+1________ijm_______

3,m+1

k—3 k—3

u.. 3 — u.. 3

ijm ijm—1

/

V

k—1 k—3

3 — її 3

/

.—1 .—

u.. 3 + u.. 3

ij+1m ijm

k—1 k—1

u., 3 — u., 3

ij+1m________ijm

3,j+1

/

k—1 k—1

u.. 3 + u.. 3

ijm ij—1m

k—1 k—1

u.. 3 — u.. 3

ijm ij—1m

5) .

Таким образом, алгоритм перехода с (к-1) -го слоя на (к +1) -й может быть описан следующим образом: Вычисляем с помощью разностной схемы с шагами Ь1,Ь2,Ь3,х значения иЬ — для каждого узла (к + 1)-го слоя.

Дважды применив разностную схему с шагами Ь Ь2 Ь3 т

2^~2^~2^ вычисляем значения иЬ т для каждого

узла (к +1) -го слоя. ь

Вычисляем с шагами —,Ь2,Ь3,т значения иЬ для 2 Т—

каждого узла (к + 1)-го слоя^

Вычисляем с шагами Ь^—,Ь3,т значения иЬ для

2 Т’т

каждого узла (к + 1)-го слоя. ь

Вычисляем с шагами Ь^Ь,,—,т значения иЬ для

2 Т’т

каждого узла (к + 1)-го слоя.

Из системы уравнений

U — uh т = С0т3 + C1h12x + C2h3x + C3h3x,

Аналогично записываются уравнения для остальных четырех этапов.

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

К полученным системам алгебраических уравнений добавляются уравнения, соответствующие краевым условиям.

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

(н+| Ml

можно представить в виде:

||h3||) [8], то погрешность ehT

п = Uijm — uijm = Т (С0т + C1h3 + C3h3 + C3h3 ),

где С!, 1 = 0...3 и и^т - неизвестные. Эта величина используется для контроля точности и выбора сетки.

Для оценки погрешности и построения новой сетки переход на очередной временной слой следует выполнить несколько раз, например, с такими величинами шагов [9]:

1) Ь1,Ь2,Ь3,т;

2) ^2 ^3 — ;

2 ’ 2 ’ 2 ’2

3) ь1,Ь2,Ь3,т;

4) Ь1,Ь2,Ь3,т;

U — u,, =

С0т3 + C^+ C3h3T + C3h3т

3 4 4 4 ,

U — uh =C0t3 + C1h т + C3h3x + C3h3x,

^f,T 4

U — % = C0T 3 + C1h3T + ^4-^ +

3,T 4

U — uh3 =C0t3 + C1h13x + C3h3x + С^зТ 1TT 4

находим неизвестные:

C = 3

C0 = T

3uh,T — uh, — uh3 — % + uh.

u!

3" 3’3'

C1 31i13t

C = 4

3h3x

3h3x

V 3’

V 3

V 3

U = 3u, + u, ----------------------

h t h,T q

3,3 3

V 3 3

3

ОтсюДа eh,T= U — uJ =

3uh T —

h t q

3,3 3

V 3

(5)

3

Если полученная погрешность меньше допустимой, то результаты шага по времени принимаются и по формуле (5) уточняем температуру.

1. Если полученная погрешность больше допустимой, то результаты аннулируются, шаг по

х

т

3,c

3

3.

7.

времени уменьшается и формируется новая функция на уплотненной сетке.

Для построения новой сетки по формуле

1

а ¡jm=т

4Cn

рассчитываем коэффициент из-

менения временного шага для каждого узла и находим минимальный из них а .

По формуле тк+1 = а-тк-1 находим новый шаг по времени.

4. По формулам e1 =

;- C0 (ат) (ат)

р.. = — ~Єч

Pljm h^3CJ

5. По

рассчитываем коэффициенты изменения шага в направлении оси Ох для каждого узла х;, по которым строим неравномерную сетку в этом направлении [9].

е - С0 (ат)2 - С (РЬ1 )2

формулам

e, =-

Y lim =

1

M2|C2|

(ат)

рассчитываем коэффициенты

6.

изменения шага в направлении оси Оу для каждого узла у^ и формируем новую неравномерную сетку в этом направлении.

По формулам

\2 ^ \2 ^ )2

Єо =

e - Co (ат)2 - Cl (phi )2 - C2 (Ph2

(ат)

1 Ге

= — і—3 рассчитываем коэффициенты h3 V C3

изменения шага в направлении оси Ог для каждого узла zm и формируем новую неравномерную сетку в направлении Ог. Если текущее время не превышает Тк, то проводится интер-

нительно, измерение проводилось на некотором расстоянии. Для этого к образцу припаивались две термопары (рис. 2): одна на лицевой стороне образца на расстоянии

7 мм от края обрабатываемой зоны, вторая - на обратной стороне непосредственно на центральной оси прохода.

Рис. 2. Схема измерения температуры образца во время обработки: 1 — пластина; 2 — термопары; 3 — лазерный луч; 4 - фиксатор

Обозначим температуру, измеренную верхней термопарой через ^ , а нижней - через ^ . В результате проведения эксперимента были получены значения ^ = 44°С и t2 = 80°С .

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

Таблица 1

Сравнение эмпирических и расчетных температур

поляция полученных значении на новую сетку и алгоритм повторяется с п. 1.

Для интерполяции полученных значений используется аппроксимация сплайнами Коханэка-Бартелса [10], что позволяет избежать осцилляций искомой функции по сравнению с использованием многочленов Лагранжа третьего порядка [5].

Метод балан- са Метод интер- поля- ции Кол- во шагов ti, °С it — tei t1 - te 1 1 11 100% t2 °С Í2 - t2 It, -12 1 2 2i 100%

te 4 te 2

линей- ный мно- гочлен Ла- гранжа 355 29,4 14,6 33,2 70,3 9,7 12,1

сплайн 348 32,0 12,0 27,3 70,5 9,5 11,9

не- линей- ный мно- гочлен Ла- гранжа 364 29,4 14,6 33,2 71,5 8,5 10,6

сплайн 359 31,9 12,1 27,5 72,2 7,8 9,7

е

2

5. Результаты моделирования

Моделирование проводилось для тонкой пластины из стали углеродной 65Г толщиной 0,5 мм, шириной 30 мм и длиной 50 мм. Мощность лазера - 0,2 кВт, скорость перемещения - 1 м/мин, диаметр пятна -1 мм. Температура окружающей среды - 27 °С. Значения теплофизических параметров материала даны в [11]. Допустимая погрешность - 10 %. Расчет проводился на компьютере с процессором Intel Core i7-3770,

8 Гб ОЗУ и тактовой частотой 3,4 ГГц под управлением ОС Windows Vista x64.

Поскольку экспериментальное определение температуры непосредственно в зоне обработки крайне затруд-

Из табл. 1 видно, что данные, полученные с использованием различных методов аппроксимации функции при переходе на новую сетку, отличаются на 2-8 %. Это можно объяснить тем, что сплайн-функции, кроме того что проходят через заданные точки, также являются устойчивыми к выбросам, т.е. не осциллируют в окрестностях точки, существенно отличающейся своим значением от соседних, в отличие от многочлена Лагранжа. Это свойство сплайнов также позволяет решить проблему получения в результате моделирования температур, ниже начальной.

Из результатов для обеих контрольных точек (табл. 1), полученных при компьютерном моделировании, видно, что использование сплайнов Коханэка-

Е

Бартелса дает на 2-6 % более близкий к экспериментальному результат в сравнении с использованием многочлена Лагранжа. Кроме того, из табл. 1 видно, что расхождение между 14е и ^ для обоих методов больше чем между 12 и 12 .

Если сравнивать результаты, полученные при использовании линейной и нелинейной схемы метода баланса, то видно, что нелинейная схема дает на 2 % результат ближе к экспериментальному. Это можно объяснить тем, что в нелинейной схеме используются значения физических параметров материала на новом слое, а не на предыдущем, как в линейной.

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

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

7. Выводы

Результаты компьютерного моделирования, проведенного на основе описанных методов, показали, что нелинейная схема метода баланса с использованием сплайнов Коханэка-Бартелса для интерполирования значений температуры при переходе на новую разностную адаптивную сетку дает более точный результат. Отметим также, что нелинейная схема метода баланса несколько более трудоемка по сравнению с линейной, так как требует выполнения итераций метода Ньютона, что увеличивает время расчета примерно на 7 %.

Література

1. Price, S. Laser forming [Текст] / S. Price // Journal of Manufacturing Science and Engineering. - 2007. - Vol. 129. - P. 117-124.

2. Shi, Y. Temperature gradient mechanism in laser forming of thin plates [Текст] / Y. Shi, H. Shen, Z. Yao, J. Hu // Optics & Laser Technology. - 2007. - Vol. 39(4). - P. 858-863.

3. Головко, Л. Ф. Моделювання температурного поля при зміцненні матеріалів лазерним випромінюванням [Текст] / Л. Ф. Головко, С. О. Лук’яненко, Д. С. Смаковський, І. Ю. Михайлова, В. А. Агеєнко // Моделювання та інформаційні технології: Збірник наукових праць Інституту проблем моделювання в енергетиці НАНУ- К.: ИПМЭ НАНУ - 2008. -С. 28-35.

4. Михайлова, І. Моделювання температурного поля з урахуванням залежності фізичних характеристик від температури [Текст] / І. Михайлова // Технологічний аудит та резерви виробництва. - 2013. - T. 5, N 4(13). - С. 12-15.

5. Калиткин, Н. Н. Численные методы [Текст] / Н. Н. Калиткин - М.: Наука - 1978. - 512 с.

6. Кутателадзе, С. С. Основы теории теплообмена [Текст] / С. С. Кутателадзе - М.:Атомиздат - 1979. - 416 с.

7. Марчук, Г. И. Методы расщепления [Текст] / Г. И. Марчук - М.: Наука. Гл. ред. физ.-мат. лит. - 1988. - 264 с.

8. Тихонов, А. Н. Уравнения математической физики [Текст] / А. Н. Тихонов, А. А. Самарский — 5-е изд. — M.: Наука - 1977. — 735 с.

9. Лук’яненко, С. О. Адаптивні обчислювальні методи моделювання об’єктів з розподіленими параметрами [Текст] / С. О. Лук’яненко — К.: ІВЦ «Видавництво «Політехніка»» - 2004. — 236 с.

10. Kochanek, D. H. U. Interpolating splines with local tension, continuity and bias control [Текст] / D. H. U. Kochanek, R. H. Bartels // ACM SIGGRAPH. - 1984. - Vol. 18. - No. 3. - P. 33-41.

11. 65Г - сталь конструкционная рессорно-пружинная - Марочник стали и сплавов [Электронный ресурс]. - Режим доступа: \ www/ URL:http://www.splav.kharkov.com/mat_start.php?name_id=265/ 07.11.2013 г. - Загл. с экрана.

3