МНОГОМАСШТАБНОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ОБРАБОТКИ ПОРОШКОВЫХ МАТЕРИАЛОВ ДЛЯ АДДИТИВНОГО ПРОИЗВОДСТВА В ИНДУКТИВНО-СВЯЗАННОЙ ПЛАЗМЕ Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

2023, Т. 165, кн. 1 С. 82-100

ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)

ОРИГИНАЛЬНАЯ СТАТЬЯ

УДК 519.63, 532.61, 66.088 10.26907/2541-7746.2023.1.82-100

МНОГОМАСШТАБНОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ОБРАБОТКИ ПОРОШКОВЫХ МАТЕРИАЛОВ ДЛЯ АДДИТИВНОГО ПРОИЗВОДСТВА В ИНДУКТИВНО-СВЯЗАННОЙ ПЛАЗМЕ

И. В. Цивильский, А. С. Мельников, А. Х. Гильмутдинов

Казанский национальный исследовательский технический университет им. А.Н. Туполева - КАИ, г. Казань, 420111, Россия

Аннотация

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

Ключевые слова: математическое моделирование, порошковые материалы, индуктивно-связанная плазма, фазовые переходы, сфероидизация, частицы «ядро-оболочка»

Введение

С развитием аддитивного производства к исходному порошкообразному сырью стали предъявляться всё более высокие требования для обеспечения изготовления продукции требуемого качества. Способ синтеза и обработки порошковых частиц в индуктивно-связанной плазме (ИСП) открывает большие перспективы и возможности в создании порошковых материалов с соответствующими свойствами. Технология плазменной обработки металлических и керамических порошков позволяет получать частицы идеальной сферической формы и необходимого размера, с заданными гранулометрией, химическим составом, морфологией. Кроме того, большой интерес представляет возможность создания композитных частиц «ядро-оболочка» из различных материалов.

ИСП достигает температуры плавления и испарения порошкового сырья, а также является химически чистым (безэлектродным), что обеспечивает идеальные условия для газофазной обработки частиц без образования оксидов и вредных примесей. Схематично плазменный реактор изображён на рис. 1. Основным элементом является кварцевая горелка 2, где возникает ИСП 4 при воздействии высокочастотного переменного магнитного поля от индукционной катушки 3. Через трубку 1

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

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

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

1. Обзор исследований и моделирования процессов плазменной

обработки

Первые математические модели процессов плазменной обработки порошков были разработаны исследовательской группой Вои1ои8 а1. [1,2]. Они изучали нена-груженную аргоновую плазму без учета движения частиц и подачи порошкового материала в реактор. В работе Ргои1х а1. [3] рассчитан характер течений газопорошковой смеси в плазме с учетом теплопереноса, но без учета испарения частиц. Группой ВегпаМ1 а1. [4] была разработана трехмерная модель нестационарных процессов, сопровождающих горение ИСП, и впервые продемонстрировано влияние асимметрии индукционной катушки и турбулентности потока на характер течений в конденсационной камере. В указанных работах не рассматривалась возможность плазменной обработки микроразмерных композитных частиц и не моделировалась их внутренняя структура.

В работе Aghaei et al. [5] подробно изучена индуктивно-связанная плазма, используемая для спектрометрического анализа материалов. Также авторами разработана численная модель динамики наночастиц меди в ИСП. В работе Benson et al. [6] разработаны двумерные модели для расчёта динамики микрочастиц в ИСП с учётом их взаимных столкновений и уменьшения размеров за счёт испарения.

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

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

2. Основные решаемые уравнения

2.1. Моделирование динамики частиц порошка на макромасштабе

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

1) Перенос частиц

Движение отдельных частиц рассчитывается по уравнению переноса

где ир = (ир- скорость частицы, ир и гир - её осевая и радиальная составляющие, й = (и, г') - скорость течения газа (плазмы), и и V - осевая и радиальная составляющие скорости газа (плазмы), |?7гег| = ^(ир — и)2 + (г>р — г>)2 - относительная скорость, р -- локальная плотность газа (плазмы), рр - плотность материала частицы, ¿р - диаметр частицы, д - ускорение свободного падения, / -дополнительные силы, т - масса частицы порошка.

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

24 Re '

Re < 0.2,

^(1 + JLRe), 0.2 < Re < 2,

с = I Ке у 16 / ' ^ — ' (2)

^ М(1+олтеа81), 2 < Ре <21, Ы

Ц (1+0.189Ре°'62) , 21 < Ре < 200,

где Ие = ^-\иге1\ - число Рейнольдса, I] - динамическая вязкость газа (плазмы).

Под дополнительными силами / в правой части уравнения (1) подразумеваются броуновское движение Fbr, сила термофореза Ftf и силы взаимодействия заряженных частиц порошка друг с другом и с электромагнитным полем (Fq, Fe, FL):

/ = FЬr + Ftf + Fq + Fe + FL.

Броуновское движение позволяет учесть хаотический характер движения частиц порошка:

о 216икв ¿о

где С _ случайное число с нулевым средним, <Ьп = 7Т2рс1Р(£р-)2Сс > 17 ~ кинематическая

вязкость, кв - постоянная Больцмана, Тд - локальная температура газа (плазмы), Сс - поправочный коэффициент Каннингема, ДЬ - шаг по времени.

Термофоретическая сила, придающая дополнительное ускорение частицам, движущимся в переменном температурном поле, определяется уравнением

67Г г/2 С,д (А' + С(Кп) 1

р{ 1 + ЗСтКп)(1 + 2К + 2С(Кп)

р =__"""Г ^'ЬЧ-^ Т у-'^уц;__—Х7Т (4)

„(Л I IV \ /1 I О V I оп Т.',-Л ,™Т ' ^ '

где Сд = 1.17, Сг = 2.18, Ст = 1.14, К = , к = -^¿«Д, кр - теплопроводность частицы, Кп = — число Кнудсена, А - длина свободного пробега молекул газа.

Я р

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

= кф (5)

где ке = 4^-, £о) ~~ электрическая постоянная, и цо - заряды частиц порошка.

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

Ёе = ЧЁ, (6)

где Е - напряженность электрического поля.

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

Ёь = ч[Ир,Б], (7)

где Б - напряженность магнитного поля. 2) Нагрев частиц

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

д = не (Тд - Тр) + Базв £(Тр4 - Т^), (8)

где /? = - коэффициент конвективного теплообмена, N11 - число Нуссельта,

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

Нагрев частиц в твердом и жидком состоянии в определенном интервале температур рассчитывается с учетом (8) как

тСр^~ = д' (9)

где ср - удельная теплоемкость материала частицы.

3) Плавление частиц

Когда температура частицы достигает температуры плавления Тт материала, начинается процесс перехода частицы в жидкое состояние. С учётом (8) необходимо решить уравнение

тНт^=Я, (10)

где Нт - удельная теплота плавления материала частицы, у - безразмерный параметр, соответствующий объёмной доле материала частицы в жидком состоянии; его значение (у = 0) соответствует полностью твёрдому состоянию частицы, (у = 1) -полностью расплавленному. Температура частицы при плавлении остаётся постоянной до полного её перехода в жидкое состояние (у > 1).

4) Испарение частиц

После полного расплавления частицы ( у = 1 ) начинается вынос паров с её поверхности при температуре выше температуры испарения Ту . Молярный поток паров определяется уравнением

I = кС(в - Сд ), (11)

где кс - коэффициент массопереноса, с = р"- концентрация паров над поверхностью расплавленной частицы, сд = Х^1^" - концентрация паров в окружающем газе, Хг - мольная доля компонента, рат - атмосферное давление, Е - универсальная газовая постоянная [9]. Давление паров руар над искривлённой поверхностью частицы определяется уравнением Кельвина

. 4аИ . . .

Риар =р8а(ехр(—--—), (12)

ЕТ р рр &р

где psat - давление насыщенных паров над плоской поверхностью, выражаемое из уравнения Антуана

в

]Е(рааг) = А-—, (13)

Тр

где А и В - коэффициенты для определённых материалов.

Снижение температуры при испарении описывается с учётом (11) выражением

Щ = 1БМНи ¿1 тср '

где М - молярная масса материала, Н - удельная теплота парообразования.

5) Кипение частиц

По достижении температуры кипения Ть начинается процесс парообразования по всему объёму частицы, при этом температура частицы остаётся постоянной до полного её испарения. Диаметр частицы изменяется с учётом (8) согласно выражению

(15)

6) Динамика паров материала

Облако паров испаряющихся частиц распространяется согласно уравнению переноса плотности пара

+ {йу • У)ру = БрАри + Б „ар, (16)

где иу - скорость паров, Бр - коэффициент диффузии паров, Буар(гр; Тр) -источниковое слагаемое, отвечающее за образование паров (при испарении частиц) и его исчезание (при конденсации паров на поверхностях).

2.2. Моделирование динамики частицы на мезомасштабе

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

[Ж + (« ' V}« = + »Аи + Г,

• и =0,

где и = (и,у) - поле скоростей смеси, р - плотность, р - динамическое давление, V - кинематическая вязкость, р - внешнее силовое поле.

Гравитация учитывается в соответствии с приближением Буссинеска: для расчета относительного ускорения дг вводится относительная плотность смеси, чтобы не вносить неустойчивость в течения [10]

-.Р ~ Р0 (1о\

9г = д-, (18)

р

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

Рьн = РоаТ,

где ро - стандартное давление, а - коэффициент теплового расширения материала, Т - температура смеси. В уравнении (17) термофорез учитывается в виде ускорения

VРгк У(роаТ) Роа^ пп,

ам =--=--=--VI. (19)

Р Р Р

Ускорение от конвективных течений вносит вклад, когда локальная температура

бинарной смеси Т отличается от комнатной температуры То :

аса™ = -ад(Т - То). (20)

Поверхностное натяжение вызывает в смеси ускорение

а= акп, (21)

где а - коэффициент поверхностного натяжения, к - кривизна поверхности раздела компонентов смеси, Ёп - нормаль к этой поверхности. Таким образом, слагаемое / в уравнении (17) представляет собой сумму ускорений (18)-(21):

/ = дг + аЬН + асаиу + авг-

В мезомасштабной модели такие свойства материалов частиц и газа (плазмы), как плотность (р), теплопроводность (к), удельная теплоемкость (ср) и кинематическая вязкость (V) задаются соответственно линейными комбинациями, зависящими от безразмерного фазового параметра <:

X = Х1(1 - <) + Х2<,

где Х1 и Х2 - любое из свойств (р, к, ср, V) соответствующего компонента бинарной смеси, представленное кусочно-заданными функциями температуры.

Для решения уравнения (21) необходимо определить границу раздела фаз и её кривизну. Для двумерной модели кривизна по Менгеру может быть представлена величиной, обратной радиусу окружности, проходящей через смежные точки границы раздела. Предлагается следующий алгоритм расчёта силы поверхностного натяжения: для каждого узла расчётной сетки определяются восемь ближайших точек (рис. 2). Радиус-вектор силы усредняется по этим смежным точкам, причём в материале частицы (^ = 1), а в газе (плазме) ^ = 0:

ег=1 ^ ег=1

Рис. 2. Иллюстрация к алгоритму определения силы поверхностного натяжения. Синяя линия соответствует границе раздела фаз

Затем вычисляем силу поверхностного натяжения, направление которой совпадает с направлением вектора г:

а вызванное ею ускорение

г =

ая+.

е;

рУ

где V - объём ячейки расчётной сетки.

Когда материал частицы достигает температуры плавления Тт, всё переданное ему тепло расходуется на разрушение кристаллической решётки, а энергия этого процесса равна удельной теплоте плавления. Теплопроводность к и удельная теплоёмкость ср существенно изменяются при фазовом переходе, а температура остаётся неизменной. Влияние удельной теплоты плавления на распределение температуры учитывается путём введения эквивалентной теплоёмкости, которая значительно возрастает вблизи точки плавления [11]. Увеличение теплоёмкости реализуется путём добавления кривой Гаусса к интерполированным табличным значениям теплоёмкости в некотором интервале ДТ вблизи температуры плавления. Аналогичный скачок теплоёмкости добавляется вблизи температуры кипения Ть. Таким образом, выражение для эффективной теплоёмкости можно записать в следующем виде

= т ЦТ) + #аехр(-^-^)2 + я6

ДТ

ДТ

1

1

где int(T) - функция, содержащая интерполированные табличные значения температуры, Ha, Hb - высоты графиков функции Гаусса (рис. 3), при которых интеграл (22) равен удельной теплоте фазового перехода:

/Tm+AT/2 T — T

На ехр(--fdT,

Hb =

rTb+AT/2

>Tb

(22)

Рис. 3. Геометрическое описание теплоёмкости материала частиц вблизи точек фазового перехода (синяя линия); интерполированные табличные значения теплоёмкости в зависимости от температуры показаны красной линией

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

3. Аналитические оценки

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

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

Использовав данные из [8] и предположив, что Ие = 0.5, С = 48, р =

100 кг/м3, рр = 850 кг/м3, ¿р = 1 • 10-5 м, и = 1 м/с, ир = 0, [7ге/ = 1 м/с,

д = 30000 • 1.6 • 10-19 Кл, ке = 8.99 • 109 Н • м2 /с2, В = 5 • 10-3 Тл, Е = 1 • 104 В/м,

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

-л Cd{u - up)\Urel\(-?-)m = h8(

100

Pp dp

г«):

850n

Н, = 1.88 • 10-7 Н,

4 850 • 10-5' 6(10-5)3

кулоновскую силу взаимодействия частиц между собой:

я2 9 (30000 • 1.6 • 10-19)2 21

Д, = = 8.99 • 109 • ---- Н, = 3.69 • 10~21 Н,

9 ("о" )2

а также силы Кулона и Лоренца, действующие со стороны электромагнитного поля:

Де = яЕ = 30000 • 1.6 • 10-19 • 1 • 104 Н, = 4.81 • 10-11 Н,

Дь = яПрБ = 30000 • 1.6 • 10-19 • 1 • 5 • 10-3 Н, = -2.4 • 10-17 Н.

Как можно заметить, влияние электростатических сил для достаточно крупных порошковых частиц (100 мкм) будет пренебрежимо мало по сравнению с газодинамическими силами (газодинамическая сила превышает силы электромагнитного воздействия более, чем в 10 000 раз), поэтому эти явления в данной работе также не учитываются. Однако в случае мелких частиц (менее 10 мкм) электростатическая сила может оказывать решающее влияние на их траектории, что отмечается другими авторами [8].

4. Программная реализация и результаты численного моделирования

4.1. Макромасштаб на ИЛТЬЛБ Уравнения (1), (8)-(10), (15)-(16) решим методом конечных разностей по неявной схеме.

По уравнению (1) введём обозначения

А=-С,\иге1\(^-), В = 3(1--)-4 ррар рр

Тогда

и в конечных разностях

¿п.

^ = А(и - г/,р) + Д,

пр" - иРп-1) = (Ап - Апр" + Б)^.

Пусть С = Ап + Б. Итерационным методом Гаусса-Зейделя будем решать следующее выражение для скорости

(п) = 4"-1) + СШ "р 1 +АЛ '

По уравнениям (8)-(9) введём обозначения

I / ьСр

Тогда

(Щ

1м

или в конечных разностях

£ = А(Тд - Тр) + Д,

Тр(п) - Тр(п-1) = (АТд - АТр(п) + Б)Л.

Пусть С = АТр + В. Итерационным методом Гаусса-Зейделя будем вычислять следующее выражение для нагрева частиц

T (n)

Tp(n-1) + Cdt

p 1 + Adt

По уравнению (10) введём обозначения

Тогда

А = 0, В =--—, C = A<f + B.

mHm

V) + B = B = C,

или в конечных разностях

^(п) _ ^(п-1) =

Итерационным методом Гаусса-Зейделя вычислим следующее выражение для плавления частиц

(в) = ^+СсЫ * 1 + АсЫ

По уравнению (15) введём обозначения

Я

A = 0, B = C -

2 РрНу Тогда

ddp _

или в конечных разностях

dpn) - dpn-1) = Cdt.

Итерационным методом Гаусса-Зейделя вычислим следующее выражение для кипения частиц

(n) = + ^

'р 1 + Adt ' По уравнению (16) введём обозначения

D v (Рж+1'2) - pvx-1>z)) v (P*,^1) _ p(x.z-1))

л , Dp 7 . vvx(pv pv ) , vvz (pv pv )

A = 4_ advect = ---+---

rri (x+1,z) (x-1,z) (x,z+1) (x,z-1)

B = ^-advect, С = А—---±»L

dtlg

Тогда

(п) _ (п-1)

21-(!1— = с- Ар(»)

dt Ру

Итерационным методом Гаусса-Зейделя найдем следующее выражение для динамики паров

{п) = + ски Ри 1 + АМ '

Таким образом, основные уравнения сводятся к выражению

а(п-1) + С ¿г

a(n)

1 + Adt

для которого может быть реализована универсальная функция-решатель на языке МЛТЬЛБ.

С помощью разработанной математической модели была выполнена серия вычислительных экспериментов для оценки траекторий и динамики движения частиц порошка, их нагрева, плавления, испарения и конденсации в конденсационной камере при известных полях распределения температуры, давления, скоростей газовых потоков. Рассматривались частицы порошка никеля диаметром 40 мкм и оксида алюминия диаметром 70 мкм. Начальная скорость частиц (скорость инжекции) 60 м/с. Рассматриваемая конфигурация технологической ИСП характеризуется параметрами [12], указанными в табл. 1. Давление в реакторе поддерживается атмосферное.

Табл. 1

Расходы газа в технологической ИСП

Газовый поток Расход, л/мин

Транспортный (несущий) газ Плазмообразующнй (вспомогательный) газ Защитный (охлаждающий) газ 3 20 30

На рис. 4 представлены результаты моделирования динамики частиц с учётом уравнений (1)-(16): для наглядности показано пять частиц оксида алюминия слева и пять частиц никеля справа. Заметно отклонение частиц от начальной траектории транспортного газового потока, связанное с движением частиц в условиях резкого градиента температуры; при этом по мере испарения частицы теряют массу и отклоняются сильнее из-за меньшей инертности. Также приведено поле распределения паров материала испарившихся частиц. Из-за характерных течений в конденсационной камере облако паров приобретает сердцевидную форму.

Тд [К] ру [кг/м3] х10-,„ с1 [мкм] Тр [К] Ф с1 [мкм]

Рис. 4. Поле температуры газа (плазмы) Тд, поле плотности паров ру материала частиц, траектории частиц в конденсационной камере. Цвет графиков соответствует диаметру частиц ^, их температуре Тр вдоль траекторий, агрегатному состоянию ^ (0 - твёрдое, 1 - жидкое, 2 - газообразное)

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

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

4.2. Макромасштаб в ANSYS

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

Использована модель дискретной фазы (Discrete Phase), а испаряемый материал частиц был представлен в виде смеси газа и пара (Species Model). Вязкость жидкостей и газа рассмотрена в рамках стандартной k - е модели турбулентности. Для учёта давления насыщенных паров при испарении решается расширенное уравнение Антуана (13) с дополнительными коэффициентами. Электромагнитные эффекты решалось с использованием ранее разработанной модели ИСП [12], основанной на пользовательских функциях ( User-Defined Functions).

Рассмотрена динамика температуры ненагруженной ИСП и ИСП с введёнными в неё частицами порошка из никеля. На рис. 5 приведены распределения полей массовой доли никеля и температуры. На рис. 5а-б наблюдается влияние потока частиц и их паров на температуру плазмы. Также заметно вращение плазмы и паров вокруг оси реактора. На рис. 5 в представлено горение плазмы без введения частиц. В этом случае быстрого охлаждения плазмы не происходит.

Рассмотрение задачи в стационарном режиме (рис. 5 г-д) подразумевает, что частицы подаются непрерывным потоком с расходом в 3 • 10-5 кг/с. При движении смеси газа и частиц учитывается взаимодействие частиц с плазмой: горячий поток нагревает и испаряет частицы порошка, а они, в свою очередь, локально охлаждают плазму вдоль своих траекторий, особенно в области интенсивного испарения частиц. В случае с введёнными частицами в усреднённом температурном поле по центру плазменного потока наблюдается холодный след, обусловленный локальным охлаждением плазмы облаком паров от частиц.

Динамика параметров частиц (диаметр, температура, модуль скорости и время пребывания) вдоль их траектории (рис. 6д-ж) демонстрирует сходство с результатами, полученными из расчётов в MATLAB (раздел 4.1).

4.3. Мезомасштаб на Python

На мезомасштабе выполнены расчёты процессов плавления и сфероидизации частиц. Рассматриваются отдельные частицы никеля неправильной (несферической) формы (рис. 7). Граничные условия следующие: на боковых и верхних стенках расчётной области задана постоянная температура 3000 К; на нижней стенке -граничное условие Неймана, имитирующее локальное охлаждение газового потока из-за присутствия частиц. Под действием силы поверхностного натяжения частицы становятся более округлыми; при этом они быстрее оплавляются в верхней части за счёт градиента температуры. Таким образом имитируется режим обтекания нисходящим потоком горячего газа (плазмы).

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

1 = 0.15 [с] 1 = 0.27 [с] 1 = 0.42 [с] 1 = 0.54 [с] 1 = 0.66 [с] 1= 1.005 [с]

Рис. 5. Рассчитанная массовая доля никеля а) и температура плазмы б) при испарении частиц; в) температура плазмы без испарения частиц; усреднённое по времени поле температуры ненагруженной плазмы г) и плазмы с введёнными частицами порошка никеля д)

паров микрокапли материала (зародыши) осаждаются на более холодной поверхности частицы-ядра, образуя тем самым однородное плёночное покрытие.

Для достижения гетерогенной конденсации (рис. 8 а) были заданы соответствующие граничные условия: концентрация паров материала частиц на верхней стенке выражена через значение фазы (р = 0.8). Таким образом может быть получено достаточно однородное плёночное покрытие. Толщина плёнки сверху и снизу частицы оказалась разной из-за разницы в скоростях течений газа (плазмы).

Обнаружен и негативный эффект в условиях повышенной концентрации насыщенных паров в окружающем частицу газе (плазме): конденсация паров происходит не на поверхности, а в объёме газа (плазмы) в виде мелких капель. За счёт тер-мофоретической силы эти капли притягиваются к относительно холодной частице, образуя так называемые «сателлиты» неправильной формы. Подобные частицы представляют собой дефект производства и негативно влияют на технологические свойства порошкового материала (рис. 8 б).

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

Рис. 6. Рассчитанные поля в реакторе: а) электрическое; б) магнитное; в) температура; г) модуль скорости. Параметры частиц вдоль их траекторий: д) диаметр; е) температура; ж) массовая доля никеля

5. Обсуждение результатов

В результате исследования была разработана математическая модель на двух масштабных уровнях. Моделирование на макромасштабе выполнено на языке программирования MATLAB. Реализован конечно-разностный решатель дифференциальных уравнений, описывающих эволюцию дискретных частиц в конденсационной камере плазменного реактора. Учтено влияние броуновской силы и термофо-реза. Процесс парообразования рассмотрен с точки зрения классического подхода (кипение происходит во всём объёме частицы при критической температуре) и согласно кинетической теории (испарение паров материала частицы с поверхности полностью расплавленной частицы при любой температуре).

Выполнено сопоставление с моделью в ANSYS: сравнивалось горение ненагру-женной ИСП и ИСП с введёнными частицами порошка никеля. Наблюдалось локальное охлаждение газа (плазмы) за счёт испарения частиц. Эти результаты качественно подтверждают результаты модели на MATLAB.

Модель на мезомасштабе реализована на языке программирования Python. Частица в окружающем газе рассматривалась как непрерывная бинарная смесь. Рассмотрены явления сфероидизации, коагуляции нескольких частиц, фазовые переходы, а также образование частиц-сателлитов и частиц «ядро-оболочка» с плёночным покрытием.

Заключение

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

Рис. 7. Результаты моделирования коагуляции и сфероидизации двух частиц в разные моменты времени. Представлены поля фазы ^, агрегатного состояния и температуры Т

Рис. 8. а) гетерогенное осаждение паров на частицу-ядро диаметром 200 мкм. Представлены поля фазы ^ и агрегатного состояния; б) формирование частицы-сателлита при гомогенной конденсации паров

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

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

Благодарности. Научные исследования проведены при финансовой поддержке Министерства науки и высшего образования Российской Федерации (Рег. номер НИОКТР АААА-А20-120122490071-1).

Литература

1. Boulos M.I. The inductively coupled R.F. (radio frequency) plasma // Pure Appl. Chem. 1985. V. 57, No 9. P. 1321-1352. doi: 10.1351/pac198557091321.

2. Boulos M.I., Gagne R., Barnes R.M. Effect of swirl and confinement on the flow and temperature fields in an inductively coupled r.f. plasma // Can. J. Chem. Eng. 1980. V. 58, No 3. P. 367-376. doi: 10.1002/cjce.5450580313.

3. Proulx P., Mostaghimi J., Boulos M.I. Plasma-particle interaction effects in induction plasma modeling under dense loading conditions // Int. J. Heat Mass Transfer. 1985. V. 28, No 7. P. 1327-1336. doi: 10.1016/0017-9310(85)90163-2.

4. Bernardi D., Colombo V., Ghedini E., Mentrelli A. Three-dimensional modeling of inductively coupled plasma torches // Pure Appl. Chem. 2005. V. 77, No 2. P. 359372. doi: 10.1351/pac200577020359.

5. Aghaei M., Bogaerts A. Particle transport through an inductively coupled plasma torch: Elemental droplet evaporation // J. Anal. At. Spectrom. 2015. V. 31, No 3. P. 631-641. doi: 10.1039/C5JA00162E.

6. Benson С.M., Gimelshein S.F., Levin D.A., Montaser A. Simulation of droplet heating and desolvation in an inductively coupled plasma - Part I // Spectrochim. Acta, Part B. 2001. V. 56, No 7. P. 1097-1112. doi: 10.1016/S0584-8547(01)00233-6.

7. Shemakhin A.Yu., Zheltukhin V.S. Mathematical modelling of RF plasma flow at low pressures with 3d electromagnetic field // Adv. Mater. Sci. Eng. 2019. V. 2019. art. 7120217. doi: 10.1155/2019/7120217.

8. Файрушин И.И. Аналитический расчет состава термической пылевой плазмы с металлическими частицами // Химия высоких энергий. 2020. Т. 54, № 6. с. 497-500. doi: 10.31857/S0023119320060042.

9. Sect. 16.4: Mixture Model Theory // ANSYS FLUENT 12.0 Theory Guide. Ch. 16: Multiphase Flows. Canonsburg, PA: ANSYS, Inc., 2009. URL: www.afs.enea.it/project/neptunius/docs/fluent/html/th/node308.htm.

10. Gray D.D., Giorgini A. The validity of the Boussinesq approximation for liquids and gases // Int. J. Heat Mass Transfer. 1976. V. 19, No 5. P. 545-551. doi: 10.1016/0017-9310(76)90168-X.

11. Luo Z., Zhao Y. A survey of finite element analysis of temperature and thermal stress fields in powder bed fusion Additive Manufacturing // Addit. Manuf. 2018. V. 21. P. 318-332. doi: 10.1016/j.addma.2018.03.022.

12. Tsivilskiy I.V., Gilmutdinov A.Kh., Nikiforov S.A., Rublya R.S., Khamidullin B.A., Melnikov A.S., Nagulin K.Yu. An experimentally verified three-dimensional non-stationary fluid model of unloaded atmospheric pressure inductively coupled plasmas // J. Phys. D: Appl. Phys. 2020. V. 53, No 45. art. 455203. doi: 10.1088/1361-6463/aba45f.

Поступила в редакцию 16.11.2022 Принята к публикации 24.04.2023

Цивильский Илья Владимирович, кандидат технических наук, доцент кафедры лазерных и аддитивных технологий

Казанский национальный исследовательский технический университет им. А.Н. Туполева - КАИ

ул. К. Маркса, д. 10, г. Казань, 420111, Россия E-mail: icivv@mail.ru

Мельников Антон Сергеевич, аспирант, инженер кафедры лазерных и аддитивных технологий

Казанский национальный исследовательский технический университет им. А.Н. Туполева - КАИ

ул. К. Маркса, д. 10, г. Казань, 420111, Россия E-mail: ASMelnikov@kai.ru Гильмутдинов Альберт Харисович, доктор физико-математических наук, профессор, заведующий кафедрой лазерных и аддитивных технологий

Казанский национальный исследовательский технический университет им. А.Н. Туполева - КАИ

ул. К. Маркса, д. 10, г. Казань, 420111, Россия E-mail: Albert.Gilmutdinov@kai.ru

ISSN 2541-7746 (Print) ISSN 2500-2198 (Online) UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)

2023, vol. 165, no. 1, pp. 82-100

ORIGINAL ARTICLE

doi: 10.26907/2541-7746.2023.1.82-100

Multiscale Modeling of Powder Materials Processing for Additive Manufacturing

in Inductively Coupled Plasma

I.V. Tsivilsky*, A.S. Melnikov**, A.Kh. Gilmutdinov***

Kazan National Research Technical University named after A. N. Tupolev - KAI,

Kazan, 420111 Russia E-mail: *icivv@mail.ru, ** ASMelnikov@kai.ru, ***Albert.Gilmutdinov@kai.ru,

Received November 16, 2022; Accepted April 24, 2023 Abstract

This article analyzes the approaches to macro- and mesoscale computational modeling of the dynamics of metal powder particles motion in the condensation chamber of a plasma reactor, spheroidization, coagulation, and phase transitions in particles. The features of different regimes of vaporization and condensation were described. The influence of phenomena such as Brownian motion and thermophoresis on the process was explored. The parameters of the process at which the formation of core-shell particles occurs were determined. The model can be used to optimize and select the effective regimes for the processing and synthesis of powder materials in inductively coupled plasma.

Keywords: computational modeling, powder materials, inductively coupled plasma, phase transitions, spheroidization, core-shell particles

Acknowledgments. This study was supported by the Ministry of Science and Higher Education of the Russian Federation (R&D project no. AAAA-A20-120122490071-1).

Figure Captions

Fig. 1. Left part: plasma reactor diagram: 1 - gas-powder mixture supply tube, 2 - quartz plasma torch burner, 3 - inductive coil, 4 - inductively coupled plasma area, 5 - condensation chamber; right part: diagram of gas supply.

Fig. 2. Illustration for the algorithm to determine the surface tension force. The blue line corresponds to the phase boundary.

Fig. 3. Geometric description of the heat capacity of the particle material near the phase transition points (blue line); the red line shows interpolated reference values of the heat capacity depending on temperature.

Fig. 4. Gas (plasma) temperature field Tg , vapor density field pv of the particle material, particle trajectories in the condensation chamber. The color of the graphs corresponds to the diameter of particles d, their temperature Tp along the trajectories, and the state of aggregation p (0 - solid, 1 - liquid, 2 - gaseous).

Fig. 5. Calculated mass fraction of nickel a) and plasma temperature b) during particle evaporation; c) plasma temperature without particle evaporation; time-averaged temperature field of unloaded plasma d) and plasma with introduced nickel powder particles e).

Fig. 6. Calculated fields in the reactor: a) electric; b) magnetic; c) temperature; d) speed module. Parameters of particles along their trajectories: e) diameter; f) temperature; g) mass fraction of nickel.

Fig. 7. Simulation results of the coagulation and spheroidization of two particles at different times. The fields of phase p , state of aggregation, and temperature T are shown.

Fig. 8. a) heterogeneous vapor deposition onto the core particle with a diameter of 200 pm. The fields of phase p and state of aggregation are shown; b) formation of a satellite particle during homogeneous condensation of vapors.

References

1. Boulos M.I. The inductively coupled R.F. (radio frequency) plasma. Pure Appl. Chem., 1985, vol. 57, no. 9, pp. 1321-1352. doi: 10.1351/pac198557091321.

2. Boulos M.I., Gagne R., Barnes R.M. Effect of swirl and confinement on the flow and temperature fields in an inductively coupled r.f. plasma. Can. J. Chem. Eng., 1980, vol. 58, no. 3, pp. 367-376. doi: 10.1002/cjce.5450580313.

3. Proulx P., Mostaghimi J., Boulos M.I. Plasma-particle interaction effects in induction plasma modeling under dense loading conditions. Int. J. Heat Mass Transfer., 1985, vol. 28, no. 7, pp. 1327-1336. doi: 10.1016/0017-9310(85)90163-2.

4. Bernardi D., Colombo V., Ghedini E., Mentrelli A. Three-dimensional modeling of inductively coupled plasma torches. Pure Appl. Chem., 2005, vol. 77, no. 2, pp. 359372. doi: 10.1351/pac200577020359.

5. Aghaei M., Bogaerts A. Particle transport through an inductively coupled plasma torch: Elemental droplet evaporation. J. Anal. At. Spectrom., 2015, vol. 31, no. 3, pp. 631-641. doi: 10.1039/C5JA00162E.

6. Benson C.M., Gimelshein S.F., Levin D.A., Montaser A. Simulation of droplet heating and desolvation in an inductively coupled plasma - Part I. Spectrochim. Acta, Part B, 2001, vol. 56, no. 7, pp. 1097-1112. doi: 10.1016/S0584-8547(01)00233-6.

7. Shemakhin A.Yu., Zheltukhin V.S. Mathematical modelling of RF plasma flow at low pressures with 3d electromagnetic field. Adv. Mater. Sci. Eng., 2019, vol. 2019, art. 7120217. doi: 10.1155/2019/7120217.

8. Fairushin I.I. Analytical calculation of the composition of thermal dusty plasma with metal particles. High Energy Chem., 2020, vol. 54, no. 6, pp. 477-479. doi: 10.1134/S0018143920060041.

9. Sect. 16.4: Mixture Model Theory. ANSYS FLUENT 12.0 Theory Guide. Ch. 16: Multiphase Flows. Canonsburg, PA, ANSYS, Inc., 2009. URL: www.afs.enea.it/project/neptunius/docs/fluent/html/th/node308.htm.

10. Gray D.D., Giorgini A. The validity of the Boussinesq approximation for liquids and gases. Int. J. Heat Mass Transfer, 1976, vol. 19, no. 5, pp. 545-551. doi: 10.1016/0017-9310(76)90168-X.

11. Luo Z., Zhao Y. A survey of finite element analysis of temperature and thermal stress fields in powder bed fusion Additive Manufacturing. Addit. Manuf., 2018, vol. 21, pp. 318-332. doi: 10.1016/j.addma.2018.03.022.

12. Tsivilskiy I.V., Gilmutdinov A.Kh., Nikiforov S.A., Rublya R.S., Khamidullin B.A., Melnikov A.S., Nagulin K.Yu. An experimentally verified three-dimensional non-stationary fluid model of unloaded atmospheric pressure inductively coupled plasmas. J. Phys. D: Appl. Phys., 2020, vol. 53, no. 45, art. 455203. doi: 10.1088/1361-6463/aba45f.

Для цитирования: Цивильский И.В., Мельников А.С., Гильмутдинов А.Х. Мно-I гомасштабное моделирование процессов обработки порошковых материалов для ( аддитивного производства в индуктивно-связанной плазме // Учен. зап. Казан. \ ун-та. Сер. Физ.-матем. науки. 2023. Т. 165, кн. 1. С. 82-100. doi: 10.26907/25417746.2023.1.82-100.

For citation: Tsivilsky I.V., Melnikov A.S., Gilmutdinov A.Kh. Multiscale modeling of / powder materials processing for additive manufacturing in inductively coupled plasma. \ Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2023, vol. 165, no. 1, pp. 82-100. doi: 10.26907/2541-7746.2023.1.82-100. (In Russian)