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

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

applied, which has the second order of accuracy for all coordinates. Systems of linear algebraic equations, created by this finite difference scheme, are solved by modified Gauss method. Adaptive method is used for automatic design of variable difference mesh. It condenses nodes in high gradient zones and dilutes them in regions where temperature changes gradually. This gives the possibility for computation time shortening and receiving a result with predefined accuracy. Results of computer modeling show that taking into account dependence of material parameters on temperature produces more accurate results. However, it requires more steps and therefore it is more time-consuming. Keywords: 3D non-stationary quasi-linear heat equation, difference approximation method.

Введение

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

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

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

Исследуемый объект - металлическая пластина с геометрическими размерами Lx, Ly , Lz (рис. 1).

Лазер

Направление движения

L,

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

Будем считать, что пластина неподвижна. Все ее поверхности находятся в процессе теплообмена с окружающей средой, температура которой ис. Для деформации пластины на ее верхнюю грань воздействует луч лазера с плотностью мощности излучения д(х, у, /), который движется со скоростью V(/) параллельно оси ординат в течение времени Тк. Необходимо найти распределение температур в металлической пластине, на которую воздействует луч лазера.

Математические модели

Рассмотрим две математические модели данного процесса.

В первой модели не будем учитывать зависимость физических параметров металла от температуры. Процесс распределения температуры в пластине описывается нестационарным трехмерным уравнением теплопроводности, в котором теплоемкость, плотность и теплопроводность - константы [3]:

ди ( x, y, z, t ) „

GO--—--- = X

к дt

д 2U ( x, y, z, t ) д 2U ( x, y, z, t ) д 2U ( x, y, z, t )

дx2

дy 2

дz2

(1)

где с - теплоемкость; р - плотность; X - коэффициент теплопроводности материала, х е [0;Lx ],

.У е [0; ^ ] , г е[0; Lz ], t е[0;Т, ].

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

Х9и (хуХ') +а[ис-и (х, у, г, I)] = 0, (2)

дп

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

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

„ ди(х, у, 0, t) , ч л

д + ч( х, У, 0 = 0. (3)

дх

Начальная температура пластины равняется температуре окружающей среды: и(х, у, г,0) = ис.

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

д

c(U)p(U)dU%yZt) я dt dx

X(u )du ( x' y,z'0

o + —

Oy

X(U ) dBxyzH

dy

o + —

dz

xu ) dUxyzll

dz

.(4)

дх

Аналогично краевые условия принимают вид

Ш Л^Е^^И + а[ис-и (х, у, х, /)] = 0, (5)

дп

^(U)дU(XдyM + q(x,y,t) = 0. (6)

дг

Начальное условие остается без изменений.

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

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

1. дискретизация расчетной области;

2. замена дифференциального уравнения в частных производных системами алгебраических уравнений;

3. решение этих систем.

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

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

Пусть приближенное решение найдено на (к -1) -м слое на неравномерной сетке в узлах

(х1,, хт), где I = 0,...,п1, ] = 0,..., п2, т = 0,...,п3, п1 - количество узлов по оси абсцисс, п2 - количество узлов по оси ординат, п3 - количество узлов по оси аппликат. Обозначим величину шагов сетки в направлении оси абсцисс через // 1 = х1 -х-1, в направлении оси ординат - через / = —у -1, в направлении оси аппликат - через И3т = гт -хт_ 1. Обозначим также среднее арифметическое двух соседних шагов для узла (х1, у., гт ) в каждом из координатных направлений через /\с, /1с, /3 с; шаг по времени через х = хк = tk -tk-1 = (к+1 -(к; значение приближенного решения в точке (х1,,гт, 1к) - через и*т .

Рассмотрим метод разностной аппроксимации [8]. Заменив производные в уравнении (1) соответствующими разностными схемами, получим операторную форму:

к 1 к 1 к 1 к 1 к 1 к 1 к ^ к ^ и 3 -ик-1 и 3 + ик-1 и 3 -и 3 „и 3 + и 3 ик-и 3 ик + и 3

-= Л3-, -= Л9-, -= Л1-,

х 3 2 х 1 2 х 1 2

, 1 1 1 ,2,1 ,2,1 ,2 ,2

к+- , к+- , к+— к+- к+— к+- , ,, к+— , к+—

и 3 - ик и 3 + ик и 3 - и 3 А и 3 + и 3 и +1 - и 3 . и +1 + и 3

-= Л1-, -= Л 2-, -= Л3-,

т 2 т 2 т 2

где Л1, Л2 , Л3 - разностные операторы аппроксимации вторых производных на неравномерной сетке:

к _ Кг+1^-1, у, То - 2\Ат + \и+1,7, т А к _ К у+1и5-1, т - с^т + \Х;+1, т

Л К i,/-1-1 '-i,/, ru ini i,/ /-1-1, / ,/// » ir

Um _ а- i,/-— ' Л 2 u^ _ а

hAAi+1 h2,ch2,Jh2,

J' % J+1

. i hi.m+1Ui, j,m-1 2h3,cU»m + \mUi, j,m+1

AiU'Jm = Й-НПЛи-' (7)

cA mr i m+1

где a = —— коэффициент температуропроводности.

Ф

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

(7) коэффициент температуропроводности примет вид a (uk ) = x(uj)

a[Um) c( Ujm )p(Ujm ) .

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

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

Моделирование проводилось для тонкой пластины из стали углеродной 65Г [9] толщиной 0,5 мм, шириной 30 мм и длиной 50 мм. Мощность лазера - 0,2 кВт, скорость перемещения - 0,0167 м/сек, диаметр пятна - 1 мм. Температура окружающей среды - 27°С. Значения теплофизических параметров, используемых при моделировании, без учета их зависимости от температуры: a = 50Вт/(м2 • К),

c = 705 Дж/(кг • К), р = 7730 кг/м3 , X = 28 Вт/(м • К), что соответствует температуре 800°С. Расчет проводился на компьютере с процессором Intel Core Í7-3770 и тактовой частотой 3,4 ГГц под управлением ОС Windows Vista x64.

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

Обозначим температуру, измеренную верхней термопарой, через t1, а нижней - через 12 . В результате проведения эксперимента были получены значения t¡ = 440C и ¿2 = 800C .

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

3 - лазерный луч; 4 - фиксатор

Результаты для контрольной точки на верхней стороне пластины, полученные при компьютерном моделировании с учетом и без учета зависимости плотности, теплоемкости и коэффициента теплопроводности от температуры, отличаются примерно на 7% (таблица). Модель с переменными параметрами дает результаты на 5% ближе к результатам эксперимента по сравнению с моделью, в которой коэффициенты постоянны. Из таблицы видно, что температура в верхней контрольной точке для второго метода выше по сравнению с первым. Это объясняется прямо пропорциональной зависимостью теплоемкости от температуры и константным значением, соответствующим температуре 800°С.

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

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

Использованный метод Количество шагов t1, °С t2, °С t - н |'1 '1 1 t - te •100% t - te 2 '2 t - te 2 '2 •100%

te '1 te '2

Метод разностной аппроксимации с постоянными с, р , X 227 29,5 70 14,5 33% 10 12,5%

Метод разностной аппроксимации с переменными с, р , X 388 31,6 90,7 12,4 28% 10,7 13%

Таблица. Сравнение эмпирических и расчетных температур

а б

Рис. 3. Изотермы на нижней грани в разрезе xOy: с постоянными c, р , X (а); с переменными c, р , X (б)

Заключение

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

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

Литература

1. Price S. Laser forming // Journal of Manufacturing Science and Engineering. 2007. V. 129. P. 117-124.

2. Shi Y., Shen H., Yao Z., Hu J. Temperature gradient mechanism in laser forming of thin plates // Optics & Laser Technology. 2007. V. 39. N 4. P. 858-863.

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

4. Кутателадзе С.С. Основы теории теплообмена. М.: Атомиздат, 1979. 416 с.

А.А. Воронин, Г.Н. Лукьянов, Е.В. Фролов

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

6. Лук'яненко С.О. Адаптивш обчислювальш методи моделювання об'екпв з розподшеними параметрами. Ки1в: Полiтехнiка, 2004. 236 с.

7. Марчук Г.И. Методы вычислительной математики. М.: Наука, 1977. 456 с.

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

9. 65Г - сталь конструкционная рессорно-пружинная. Марочник стали и сплавов [Электронный ресурс]. Режим доступа: http://www.splav.kharkov.com/mat_start.php?name_id=265, свободный. Яз. рус. (дата обращения 07.11.2013).

10. Михайлова 1.Ю. Споаб функцюнального тестування точносп результапв, отриманих з використан-ням адаптивно! сггки // Математичне та комп'ютерне моделювання. Сер. Техтчт науки. 2012. В. 7. С. 124-132.

Лукьяненко Святослав Алексеевич -

Михайлова Ирина Юрьевна

Svyatoslav Luk'yanenko Irina Mykhailova

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

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

- D.Sc., Professor, National Technical University of Ukraine «Kyiv Polytechnic Institute», Kyiv, Ukraine, [email protected]

- assistant, National Technical University of Ukraine «Kyiv Polytechnic Institute», Kyiv, Ukraine, [email protected]

УДК 532.542.4:533.6.011.32:612.215.41

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТУРБУЛЕНТНОГО ПОТОКА ВОЗДУХА С ИСПОЛЬЗОВАНИЕМ МЕТОДА ОТСОЕДИНЕННЫХ ВИХРЕЙ А.А. Воронин3, Г.Н. Лукьянов3, Е.В. Фролов3

а Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, Санкт-Петербург, Россия, [email protected]

Приведен краткий обзор основных математических моделей описания потоков жидкостей и газов, от модели пути смешения Прандтля до одно- и двухпараметрических дифференциальных RANS-моделей, а также нестационарных моделей крупных LES- и отсоединенных DES-вихрей. Последний тип моделей, впервые предложенный П. Спалар-том в 1997 г., позволил объединить преимущества моделей LES и RANS, что впоследствии расширило область применения нестационарных математических моделей к описанию потоков жидкостей и газов. Авторы исследования приводят математическую формулировку моделей DES, подробно останавливаясь на принципиальных отличиях данных методов от моделей RANS.

Также приведен вид трехмерной геометрической модели носовой полости человека, полученной на основе данных компьютерной томографии с использованием программного пакета Mercury Amira. После сегментации данной модели с помощью программного пакета Altair Hypermesh была построена объемная нерегулярная сетка из 1,5*107 конечных элементов, на основе которой произведен нестационарный расчет параметров потока (программный пакет Ansys Fluent). Приведены полученные в результате расчета поля скоростей потока для вдоха и выдоха. Использование подробной расчетной сетки и нестационарной модели DES позволило выделить внутри потока отдельные мелкомасштабные вихревые структуры.

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

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

DETACHED-EDDY SIMULATION OF TURBULENT AIRFLOW A. Voronin3, G. Luk'yanov3, E. Frolova

a Saint Petersburg National Research University of Information Technologies, Mechanics and Optics, postgraduate, Saint Petersburg, Russia, [email protected]

A brief survey of the most significant mathematical models of the air and fluid flows from Prandtl's mixing length layer theory to RANS models with one and two differential equations, as well as unsteady LES and DES models is given. Detached-eddy simulation was first proposed by P.Spalart in 1997 and combined the main advantages of LES and RANS methods which gave the researchers the possibility to widen the sphere of such models application. The authors give the basic mathematical description of DES models paying special attention to the most important differences between DES and RANS models. The 3D geometrical model of human nasal cavities obtained from computer-aided tomography data using Mercury Amira program is also given. The 3D unstructured mesh with 1,5*107 finite elements was constructed after the segmentation using Altair Hypermesh software had been finished. The mesh was used to set up an unsteady simulation of airflow inside the obtained geometrical model. Application of DES method on the mesh of а good quality made it possible to distinguish the