Mathematical modeling of non-stationary heat transfer in selective laser melting based on machine learning

Full Text

Введение

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

Тепловой расчёт процесса печати в трёхмерной постановке ввиду нестационарного характера и нелинейных эффектов в материале сопровождается значительными вычислительными затратами. В частности, опыт использования программного комплекса Ansys Additive¹ показывает, что время расчёта реальных изделий может измеряться сутками. Поэтому актуальной является проблема ускорения численного моделирования процесса 3D-печати.

Задача упрощается для расчётов однотипных изделий. В этом случае целесообразно обучить нейронную сеть (НС) и в последующем использовать её для быстрого определения температурного поля при различных комбинациях исходных параметров. НС в настоящее время используются в различных предметных областях: от генерации текста и мультимедиа контента [4, 5] до аппроксимации сложных математических моделей (ММ) физических систем [6].

НС и машинное обучение в области аддитивных технологий применяются следующим образом [7]: при проектировании изделий для 3D-печати – это синтез материалов с заданными свойствами и проведение топологической оптимизации; при выборе параметров технологического процесса – оптимизация и мониторинг дефектов печати в реальном времени; при планировании аддитивного производства – для контроля качества готовых изделий.

Для оценки принципиальной возможности использования НС для приближённого описания процессов теплопередачи при 3D-печати в данной работе рассматривается одномерная идеализация теплового процесса, когда учитывается только распространение тепла вдоль оси печати (оси Z). Использование такой модели представляется вполне оправданным в случаях, когда высота изделия значительно превышает два других габаритных размера. Одномерная идеализация при моделировании теплопереноса в аддитивном процессе использована в работе [8]. Решение нестационарной тепловой задачи оптимального управления технологическим процессом индукционного нагрева заготовки рассмотрено в работах [9, 10].

1 Математическая модель теплопередачи в процессе 3D-печати

В данной работе рассматривается ММ теплопередачи в процессе послойного синтеза в одномерной постановке, когда изменение температуры учитывается только в направлении оси печати Z. При этом геометрическая модель изготавливаемой конструкции может быть описана в виде стержня с изменяющимся вдоль оси Z поперечным сечением $A=A\left(z\right)$.

Количество теплоты Q, необходимое для нагрева массы m с некоторой температуры $T_{*}$ до температуры T, вычисляется по известной формуле: $Q=cm\left(T-T_{*}\right)$, где с – удельная теплоёмкость материала. Для вычисления количества теплоты в единице объёма стержня необходимо продифференцировать обе части этого уравнения по объёму V:

$\frac{dQ}{dV}=c\rho \left(T-T_{*}\right)$,

где $\rho =dm/dV$ – плотность материала. Тепловая мощность в единице объёма находится путём дифференцирования полученного соотношения по времени t, что с учётом независимости $T_{*}$ от t приводит к формуле:

$\frac{d\dot{Q}}{dV}=c\rho \frac{\partial T}{\partial t}$,

где точкой над Q обозначена производная по времени. Тепловой поток q в стержне, по определению, представляет собой тепловую мощность, приходящуюся на единицу поперечного сечения: $q=\dot{Q}/A$, откуда $\dot{Q}=q\cdot A$. Учёт того факта, что элемент объёма dV выражается через площадь сечения как $dV=A\cdot dz$ даёт

$\frac{1}{A}\cdot \frac{\partial }{\partial z}\left(q\cdot A\right)=c\rho \frac{\partial T}{\partial t}$.

Связь градиента температуры и теплового потока определяется законом Фурье:

$q=k\frac{\partial T}{\partial z}$,

где $k$ – коэффициент теплопроводности. Подстановка этого выражения в предыдущую формулу в предположении постоянства k по координате z приводит к дифференциальному уравнению теплопроводности для стержня переменного сечения следующего вида:

$\frac{\partial }{\partial z}\left(A\frac{\partial T}{\partial z}\right)=\frac{\rho cA}{k}\frac{\partial T}{\partial t}$. (1)

Следует заметить, что при постоянной площади сечения величина A выносится из-под оператора дифференцирования по z и, после очевидного сокращения, формула (1) трансформируется в классическое одномерное уравнение нестационарной теплопроводности. Граничное условие (ГУ) на «подвижном» (наплавляемом) конце стержня соответствует поддержанию заданной максимальной температуры:

${\left.T\left(z,t\right)\right|}_{z=L\left(t_i\right)}=T_{\max }.$ (2)

ГУ на «неподвижном» конце стержня задаётся в виде уравнения конвекции:

$h\left({\left.T\left(z,t\right)\right|}_{z=0}-T_{\min }\right)=k\frac{\partial T}{\partial z}.$ (3)

Здесь T_min – температура окружающей среды, $h$ – коэффициент конвективной теплоотдачи. При $h\to \infty$ данное ГУ трансформируется в условие поддержания постоянной T_min на неподвижном конце стержня. Принимая для $h$ конечные значения, можно обеспечить физически адекватное значение температуры на неподвижном конце стержня.

Уравнение теплопроводности (1) решается итерационно в цикле по временным шагам. На каждом шаге $\Delta t$ происходит увеличение длины стержня на величину $\Delta L$. Длина стержня $L$ увеличивается с течением времени, что соответствует пошаговому добавлению материала при 3D-печати. Начальное распределение температур задаётся на основе поля температур на предыдущей итерации:

${\left.T\left(z,t\right)\right|}_{t=t_i}={\left.T\left(z,t\right)\right|}_{t=t_{i-1}}$.

Для решения уравнений в частных производных в Matlab используется функция «pdepe».

Правила использования данной функции предполагают, что уравнение должно быть представлено в виде [11]:

$c\left(z,t,T,\frac{\partial T}{\partial z}\right)\frac{\partial T}{\partial t}=z^{{}^{-n}}\frac{\partial }{\partial z}\left(z^{{}^n}f\left(z,t,T,\frac{\partial T}{\partial z}\right)\right)+s\left(z,t,T,\frac{\partial T}{\partial z}\right)$. (4)

где функции $c$, $f$ и $s$, а также коэффициент n задаются таким образом, чтобы уравнение (4) в точности соответствовало решаемой задаче. В данном случае для величин с, m, f, s получаются следующие значения:

$с\left(z,t,u,\frac{\partial T}{\partial z}\right)=\frac{\rho c_{{}_p}A}{k},f\left(z,t,T,\frac{\partial T}{\partial z}\right)=A\frac{\partial T}{\partial z},s=\mathrm{0,}n=0.$. (5)

Подстановка этих значений в уравнение (4) приводит к исходному уравнению (1). Общий вид выражения для ГУ в Matlab следующий:

null. (6)

Здесь функции $p_i$ и $q_i$ должны быть заданы для обоих концов стержня: условно левого (индекс l) и условно правого (индекс r) – неподвижного и подвижного концов стержня соответственно. Для реализации ГУ по обеспечению заданной температуры необходимо на подвижном конце стержня задать

$\begin{array}{l}{p}_r\left(z,t,T\right)=T-T_{\max },\\ q_r\left(z,t\right)=0.\end{array}$

где T – текущая температура стержня, T_max – температура на подвижном конце стержня.

ГУ конвекции на неподвижном конце стержня задаётся как

$\begin{array}{l}{p}_l\left(z,t,T\right)=A_{\text{plate}}h\left(T_{\min }-T_0\right),\\ q_l\left(z,t\right)=-k_{\text{plate}}.\end{array}$

Здесь $A_{\text{plate}}$ – площадь поперечного сечения плиты, $k_{\text{plate}}$ – коэффициент её теплопроводности, T₀ – температура окружающего плиту пространства. Подстановка этих значений в (6), с учётом выражения для функции f в формуле (5), приводит к условиям (2) и (3) соответственно.

Зависимость площади сечения от координаты z задаётся функцией

$A\left(z\right)=\left\{\begin{array}{l}{A}_{\text{plate}},z\le L_{\text{plate}},\\ A_{\text{part}},z>L_{\text{plate}}.\end{array}\right.$

Здесь $A_{\text{part}}$ – площадь поперечного сечения детали, $L_{\text{plate}}$ – толщина опорной плиты.

В качестве тестовой задачи использован образец, представленный на рисунке 1, на котором показан процесс печати в текущей стадии, где итоговые очертания детали показаны пунктиром. Размеры на рисунке указаны в миллиметрах. Свойства материала соответствуют алюминиевому сплаву: ρ = 2700 кг/м³, с_p = 900 Дж/(кг·К), k = 122 Вт/(м·К). Коэффициент теплоотдачи h принят равным 500 Вт/м².

 

Рисунок 1 – Геометрия образца для тестовой задачи

 

Максимальная температура, поддерживаемая на верхнем конце детали, соответствует температуре плавления типового алюминиевого сплава и составляет 660 °C.

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

$T\left(z\right)=\left\{\begin{array}{l}{T}_{{}_{\min }}+\frac{T_{{}_{\text{mid}}}-T_{{}_{\min }}}{z_{\text{plate}}}\cdot z,0\le z\le z_{{}_{\text{plate}}}\\ T_{{}_{{}_{\text{mid}}}}+\frac{T_{{}_{\max }}-T_{{}_{\text{mid}}}}{z_{{}_{\text{part}}}-z_{{}_{\text{plate}}}}\cdot \left(z-z_{{}_{\text{plate}}}\right),z_{{}_{\text{plate}}}\le z\le z_{{}_{\text{part}}}\end{array}\right.$ (7)

где $z_{\text{plate}}=L_{\text{plate}}$ – толщина опорной плиты, $L_{\text{part}}$ – высота детали, $z_{\text{part}}=L_{\text{plate}}+L_{\text{part}}$.

 

Рисунок 2 – Распределение температуры в конечный момент времени

 

Для вычисления температуры на границе раздела «плита-деталь» $T_{\text{mid}}$ и температуры в основании опорной плиты $T_{\text{min}}$ можно воспользоваться законом сохранения энергии:

$\left\{\begin{array}{l}{q}_{\text{plate}}{A}_{\text{plate}}=q_{\text{part}}{A}_{\text{part}},\\ q_{\text{conv}}{A}_{\text{plate}}=q_{\text{plate}}{A}_{\text{plate}}.\end{array}\right.$

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

$\left\{\begin{array}{l}-\frac{A_{{}_{\text{plate}}}}{L_{{}_{\text{plate}}}}{T}_{{}_{\min }}+\left(\frac{A_{{}_{\text{plate}}}}{L_{{}_{\text{plate}}}}+\frac{A_{{}_{\text{part}}}}{L_{{}_{\text{part}}}}\right)T_{{}_{\text{mid}}}=\frac{A_{{}_{\text{part}}}}{L_{{}_{\text{part}}}}{T}_{{}_{\max }},\\ \left(h+\frac{k_{{}_{\text{plate}}}}{L_{{}_{\text{plate}}}}\right)T_{{}_{\min }}-\frac{k_{{}_{\text{plate}}}}{L_{{}_{\text{plate}}}}{T}_{{}_{\text{mid}}}=h{T}_{{}_0}.\end{array}\right.$ (8)

В результате получается система двух линейных уравнений относительно $T_{\text{mid}}$ и $T_{\text{min}}$. Ввиду громоздкости общего решения системы (8), в программном коде она решается численно методом Гаусса.

Таким образом, системы уравнений (7) и (8) определяют ММ теплового процесса, которая используется для обучения НС.

2 Модель теплопереноса в 3D-печати при помощи НС

Рассмотрена задача построения приближённого представления ММ, описанной в разделе 1, при помощи искусственной НС. В качестве параметров для обучения НС и получения распределения температуры $T\left(z\right)$ исследованы четыре геометрических размера: площадь сечения плиты $S_{\text{plate}}$ и образца $S_{\text{part}}$, их длины $L_{\text{plate}}$ и $L_{\text{part}}$. Свойства материала и ГУ остаются неизменными. Таким образом, температура является функцией пяти аргументов: $T=T\left(z,S_{\text{plate}},L_{\text{plate}},S_{\text{part}},L_{\text{part}}\right)$. Заданы следующие диапазоны изменения этих параметров: $p_{i\min }=0.5p_{i0},p_{i\max }=1.5p_{i0},$ где для краткости обозначено $p_i=S_{\text{plate}},S_{\text{part}},L_{\text{plate}},L_{\text{part}}$, а $p_{i0}$ – некоторое начальное значение параметра. Диапазон изменения координаты определяется максимальными значениями длин плиты и детали: $z_{\text{m}in}=\mathrm{0,}{z}_{\max }=L_{\text{plate max}}+L_{\text{part max}}.$

Для обучения НС необходимо в диапазоне изменения параметров сгенерировать обучающую выборку. Это набор комбинаций параметров и соответствующей им температуры, вычисленной при помощи ММ. Здесь можно воспользоваться планом эксперимента, например, таким как центральный композиционный план или план Бокса-Бенкена [13]. В данной работе использован подход, в котором количество точек в пространстве параметров ограничено заданным числом n_train = 1000. Обучающая выборка строится следующим образом:

$\left\{p_i\right\}=\text{rand}\left(p_{i\min },p_{i\max },n_{train}\right)$, $\left\{z\right\}=\text{rand}\left(\mathrm{0,}{z}_{\text{m}ax},n_{train}\right)$,

где «rand» – функция, формирующая массив случайных чисел размерности n_train в заданном диапазоне. Построение обучающей выборки завершается вычислением температуры по формулам (7) и (8) для каждой комбинации параметров p_i и координаты z.

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

Архитектура НС показана на рисунке 3.

 

Рисунок 3 – Архитектура используемой нейронной сети

 

Программная реализация НС осуществлена на языке Python с использованием библиотеки TensorFlow³. Архитектура НС определена на основе примеров из документации к данной библиотеке. В качестве функции активации для всех слоёв использован tanh(x). Для обучения НС использован алгоритм оптимизации Adam [15]. Максимальное количество «эпох» – просмотров обучающей выборки – равно 100. Тестовые расчёты показывают, что данный параметр в совокупности с типом функции активации и методом оптимизации оказывают наиболее существенное влияние на точность нейросетевой аппроксимации. В качестве функции потерь принято среднеквадратичное отклонение. Для ускорения обращения к обученной НС реализовано сохранение архитектуры сети и её весов в иерархическом формате данных (HDF5) с последующим импортом перед использованием.

3 Использование НС для определения температур в образце

С использованием обученной НС построено распределение температур в конечный момент времени в изготовленном образце. Набор входных параметров определён следующим образом. Задано количество тестовых точек (комбинаций значений параметров) n_test =50. Определены векторы входных параметров {p_i} размерности n_test и постоянные значения p_i0, например, $L_{plate0}=\mathrm{0,04}$ м, $L_{part0}=\mathrm{0,1}$ м, $S_{plate0}=\mathrm{0,0144}$ м², $S_{part0}=\mathrm{0,0016}$ м². Вектор координат {z} заполняется арифметической прогрессией от 0 до $L_{plate0}+L_{part0}$. Полученный набор значений формирует тестовую выборку, которая нормируется на отрезке [-1,1] и подаётся на вход НС. Вычисляются температуры по заданным входным параметрам. Полученный на выходе вектор температур приводится к исходному масштабу, поскольку на этапе обучения он подвергался нормированию.

Результаты расчёта температур в образце с использованием НС приведены на рисунке 4.

 

Рисунок 4 – Распределение температур по математической модели и нейросети

 

Выводы

Сопоставляя распределения температур по высоте плиты и детали, полученные с использованием ММ и НС (рисунок 4), можно сделать вывод, что НС обеспечивает достаточно высокую точность. Этот результат показывает принципиальную возможность нейросетевой аппроксимации ММ физических систем, поведение которых описывается дифференциальными уравнениями в частных производных. Деталь, выращиваемая на 3D-принтере, в совокупности с процессами теплопереноса, относится к таким системам.

 

¹ Модуль для моделирования процессов 3D-печати Ansys Additive. https://www.ansys.com/products/additive.

² Библиотека численных методов NumPy для языка Python. https://numpy.org/doc/stable/reference/generated/numpy.interp.html.

³ Библиотека Tensor-Flow (https://www.tensorflow.org/) для языка Python.

About the authors

Evgeniy A. Kishov

Samara University (Samara National Research University named after academician S.P. Korolev)

Author for correspondence.
Email: evgeniy.kishov@ssau.ru
ORCID iD: 0000-0002-6893-1894
Scopus Author ID: 57039050400

PhD, Associate Professor of the Department of Aircraft Design and Construction

Russian Federation, Samara

References

Supplementary files