Конечно-элементное моделирование соударения капли расплава с подложкой при плазменном напылении Текст научной статьи по специальности «Физика»

Конечно-элементное моделирование соударения капли расплава с подложкой при плазменном напылении

О.П. Солоненко, Э.П. Шурина, А.А. Головин

Институт теоретической и прикладной механики СО РАН, Новосибирск, 630090, Россия

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

Условные обозначения

и = (иг, иг) — вектор скорости; Т — температура;

Ек — кинетическая энергия;

Ер — потенциальная энергия поверхностного натяжения; Еф — энергия вязкой диссипации; а — коэффициент поверхностного натяжения; р — удельная плотность; ц — коэффициент динамической вязкости;

V =ц/р — коэффициент кинематической вязкости;

С — удельная теплоемкость;

X — коэффициент теплопроводности; а = Х/рС — температуропроводность;

Ь — теплота плавления; Ор — диаметр капли;

Я — радиус эквивалентного цилиндра; h — высота затвердевшего слоя;

I — время; г, г — радиальная и продольная координаты;

№е = Ррт Л,ир20 / а® — число Вебера;

Ке = °рироР Р» / Цр!, — число Рейнольдса;

Рг = Ц?» С®,/ Хрт — число Прандтля;

Ре = Re Рг — число Пекле; Fo = ар» — число Фурье;

С(1) Т

^ = рт рт------число Стефана;

Ьр»

Ки£! = Ь /[ер» Тр»] — критерий Стефана-Кутателадзе;

К^,р) = (Х'^»/ХР»^ ар»/ а)® — критерий тепловой активности материала основы к расплаву частицы;

5р,ь = Ррт Ьрш/ Рь» Ь» — критерий относительной интенсивности плавления (затвердевания); п — внешняя нормаль;

2 г _ 2

G(г) =^= I е а За — интеграл ошибок.

^ о

Безразмерные переменные:

t = tUpо/Ор — время; Ф = ТТр» — температура;

; = ^°р, г = г/°р, ; = ^°р.

Верхние индексы:

s, 1 — твердое, жидкое состояние;

~ — безразмерная величина.

Нижние индексы:

р, Ь — частица, основа;

0 — начальное значение параметра; т — значение параметра при температуре плавления.

1. Введение

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

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

© Солоненко О.П., Шурина Э.П., Головин А.А., 2001

В работе [1] была предложена модель, в которой капля в начальный момент времени заменяется эквивалентным цилиндром, сохраняющим ее массу, а также кинетическую и потенциальную энергию. Высота деформирующегося слоя расплава зависит только от времени. Для определения фронта затвердевания используется решение задачи Стефана, не учитывающей отвод тепла в подложку и вклад конвекции расплава в теплообмен с основой. Эволюция жидкого цилиндрического слоя определяется из уравнения энергии, при этом используется простейшее поле скоростей, удовлетворяющее уравнению неразрывности. Автор отмечает, что в такой постановке невозможно точно определить параметры эквивалентного цилиндра, а полученное приближенное решение может вносить ошибку порядка 10 % в начальное значение потенциальной энергии.

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

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

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

ячейка полностью занята жидкостью, и нулю, когда ячейка свободна. Ячейки со значением 0< F < 1 содержат свободную поверхность. Подобным образом определенная функция удовлетворяет уравнению переноса

З F|dt + (V -V) F = 0.

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

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

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

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

Цель настоящей работы — исследование взаимодействия капли расплава с подложкой с учетом протекающих при этом фазовых превращений посредством решения двумерной задачи теплопроводности и привлечения модели эквивалентного цилиндра, предложенной в [1, 2] для учета гидродинамических процессов; разработка и тестирование алгоритма триангуляции, ориентированного на построение сеток в деформирующихся областях; исследование возможности использования аналитических приближений при задании начальных условий в системе капля расплава - подложка.

2. Физическая модель

Будем считать, что капля несжимаемой жидкости с начальной скоростью ир0, диаметром Ор и темпера-

Таблица 1

Возможные соотношения между контактной температурой и температурами плавления материалов частицы и основы

Номер и описание сценария Соотношение температур

1. Деформация и одновременное затвердевание капли на твердой подложке Т 3 > < Ь-а 3

2. Деформация и одновременное затвердевание капли расплава с подплавлением подложки в пятне ее контакта с частицей Т 3 > > Ь-а 3

3. Деформация капли без затвердевания и без подплавления основы Трт < Тсо < ТЬт

4. Деформация капли без затвердевания с одновременным подплавлением основы т Л V Є

турой Тр0 соударяется с подложкой, начальная температура которой ТЬ0. В данной работе рассматриваются все четыре сценария возможного взаимодействия “капля расплава - основа” [10]. Эти сценарии определяются соотношением между температурой Тс0, устанавливающейся в контакте, и температурами Трт и ТЬт плавления материалов частицы и основы.

Описание сценариев представлено в таблице 1, где

То =

То./АгРьСу + Тьот/АЇРС7

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

Как и в [2], будем предполагать, что при ? = 0 капля мгновенно принимает форму эквивалентного цилиндра радиусом Яо и высотой ^0), который имеет ту же массу, потенциальную и кинетическую энергию, что и капля перед соударением. Далее жидкость начинает растекаться, при этом возможны процессы фазовых превращений как в капле, так и в подложке.

Примем следующие допущения:

- линия раздела капля - подложка не деформируется;

- высота слоя жидкости является функцией времени, но она постоянна в пространстве в каждый момент времени;

- свойства материалов не зависят от температуры, но различны для разных агрегатных состояний;

- радиационный теплообмен капли с окружающей средой пренебрежимо мал;

- на линии соприкосновения капли с подложкой устанавливается идеальный тепловой контакт;

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

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

Распределение температур в капле и подложке подчиняется уравнению переноса тепла

7\Т

Р с — + р С и grad(T) = Егаё(Т)}, (1)

дt

где А, р, С — свойства, зависящие от материала.

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

- на линии фазового перехода в капле г =

ьр8) > 0

дТ ■рш дп

-о дп+ о

dЬS) “рт^рт' dt

где п — внешняя нормаль к границе затвердевания;

- на линии фазового перехода в подложке г = Ь®, Ь® < 0

дТ ,Л(1) дТ = г _(1) Т = Т . (3)

АЬт д + АЬт д іЬтРьт , , 1 ТЬт’ (3)

дп-о дп+о dt

- условие идеального контакта на границе контакта капли с подложкой г = 0,0<г<Я

А» дТ +а(8) дТ = о Т = Т АЬт дп-о +Арт дп+о ’ Тп-о Тп+о’

(4)

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

- на остальных границах задается однородное условие второго рода

дТ

дп

= о.

(5)

Начальные условия для уравнений переноса тепла: Т = Тро, г = 0, г > 0,

Т = Тм, г = 0, г < 0. (6)

Для различных сценариев взаимодействия капли и подложки существуют аналитические оценки, позволяющие оценить текущую толщину затвердевшего слоя в капле йр8) (Ро), подплавленного слоя в подложке йь(1) (Ро) и температуру на границе контакта капли с подложкой Фс. Согласно [10-13] эти оценки имеют вид:

~р(8) (Ро) = л/Р0, ~ь(1) (Ро) = с^ТРО. (7)

Значения коэффициентов с^ и с^, определяющих динамику затвердевания расплава и возможного под-плавления подложки в пятне контакта ее с частицей, а также безразмерная температура в контакте вс = = Тс/Трт, отвечающие конкретному сценарию формирования сплэта, представлены в табл. 2.

Величина с?, характеризующая скорость роста фронта затвердевания, позволяет получить оценки конечной толщины и радиуса сплэта [10]:

Ь?* = с,

•р с^ЫFo = (1 - PeFo ),

6{&)* = 2Д/бЬр(8)*

(8)

(9)

где

Fo =

с?| ф + 4Ре/ с^ -1 | /2Ре

Для решения задачи нужно также определить динамику деформации жидкого диска. В нашей постановке это означает определение закона изменения величин

й®, R, dR/dt, полностью определяющих параметры эквивалентного цилиндра. Необходимые соотношения могут быть получены из уравнения сохранения энергии [1, 2]:

где

— ^ ^ ч —ЕФ

— (Ек + Ер) = - —Ф dt р —

Я 1

Ек = 2п| - рр“(и2 + 42 )—zr—r, о2

Ер =а«(пЛ2 + 2пЛЬ(1)),

(10)

—Е,

—

-=/ м-рт Ф—и

Ф = 21

диг

дг

+ 21

диг

дг

(диг ди2

+ | —- + —2 I дг дг

+ 2

Таблица 2

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

Номер

сценария

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

Р =

. = РІ^/Ї+^Р7-1]/2, с^= о

п Ар’^Ки® + 2(1 + ас)К?,р)(Оро -1)

л/П К <Ь,р>Ки®

Q=-

2А(р,р) (1 -вьо)

Кир1»

(1 + ас)(«ро -1)

(1 -«Ьо)К £

(Ь,р)

«с = ^П+А(р;;)К^Ь,р)«ьоссУ(^+А®;)К[Ь,р)сс)

сс= р - 1]/2,

с^= 5р,ь [ + 2(—2 - —,)/Ки®]

Р = 2[—2 (25р,ь +АЬр>) - —.5 р,ь]

Кир1)(5р,ь + <>)

= 2 [2—2(—2 - —1)5р,ь -А(Ь:рКир1>(1 -«Ьт)]

К® )2(5Р,Ь +А(Ь:;>)

—1 = ^(«Ьт -«Ьо)Д/П , —2 = (1 + Яс)(«ро - 1)/т/П

«с = 1 -

а(Ь;;>(1 -«Ьт)С

Ьт)сС

2(—2 - —1)

Кир

1+

А(М) 'ї А Ь,р

5 р,ь

С£ _ о, С\~ о, «с _ Тсо

= о, с = Р|У 1 -4^Р2 -1^2,

_л/ПаЬр_ + А(Ь,р) (1 -«ь^«ьт)

1 + ас 4ПК <Ь-р)КиЬ8)

АМ(1 -«Ы>/ «Ьт)

(1 + ас)К<Ь-р)КиЬ8>

1 - (1 + ас)(«ро/«Ьт -1) . КРЬ'р>(1 -«Ь^«Ьт)

«с = ^«Ьт + (1 + ас)Ар,11)«роС5 ] / + (1 + ас)Ар,1і)С

]

Здесь ^с = 0.26 — коэффициент, определяющий относительный вклад конвекции расплава в окрестности точки его торможения в сравнении с кондуктивным теплообменом между растекающейся частицей и основой

3. Вычислительная схема

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

-1

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

1. Задание входных данных (начальная скорость, диаметр и температура капли, начальная температура подложки, материалы частицы и основы, их термодинамические свойства и параметры переноса).

2. Определение начальных параметров эквивалентного цилиндра.

3. Пока не выполнено условие завершения счета /гр(1) < eps:

3.1. определение величины текущего шага по времени;

3.2. пока не найдено точное положение линий фазовых переходов для текущего шага по времени:

3.2.1. определение параметров эквивалентного цилиндра;

3.2.2. построение сетки;

3.2.3. расчет поля температур;

3.2.4. уточнение положения линий фазовых переходов.

Рассмотрим некоторые этапы более подробно.

3.1. Выбор шага сетки и шага по времени

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

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

Еще одним моментом, тесно связанным с выбором шага по времени, является наличие эффекта отката (изменение направления скорости растекания эквивалентного цилиндра). Это требует точного определения момента остановки, для которого dR/dt = 0. Для решения этой задачи мы использовали следующий алгоритм определения ^ого шага по времени:

1) определяется величина Д^, обусловленная процессами деформации и фазовых превращений;

2) экстраполированием оценивается

t =t■_1 +Д‘.

„ч dRi _ dR

3) если —и имеют одинаковый знак, то dt dt

Дti =Д^ и на этом определение текущего шага по времени заканчивается;

dRi_l dRi

4) если —и —~ противоположно направлены, dt dt

то экстраполированием определяется ДtiR, такое, что

= 0, и тогда Дti = ДtiR.

t = и. +Дti,

3.2. Триангуляция области решения

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

3.3. Решение уравнения теплопереноса

При решении уравнения теплопереноса (1) использовалась неявная схема дискретизации по времени:

- Тп+1 _ Тп

рС

+ рШ Ега^т; +1) = +1)), О1)

_ + рС-

+ рШ grad(7;+l) _рС-

; +1 ;

= 0.

(12)

рс

Введем обозначения и = Тп+,, F = Тп, у = -

л оч 1п+1 _1п

и тогда (12) примет вид

-drv(Agrad('o)) + уи + рСUgrad('o) _ YF = 0. (13)

Условия на границах фазовых переходов (2) и (3) примут вид

Л (8) , -1(1) = Г р(1) ^ П +1 ^ П

_ т л г™ _ -^ршгрш

дп.

-0

дп.

+0

.(8) ди , .(1) ди АЬш ^ ГАЬш 3

дп-0 дп

(1)

’ ^ЬшрЬ]

(1) hb ' п+1 _ hb ' п

(1)

(14)

(15)

+0

Помимо (14), (15) граничными условиями для (13) будут условия (4), (5), записанные относительно и.

Для решения сформулированной эллиптической краевой задачи с несамосопряженным оператором использовался метод конечных элементов в слабой форме (метод Галеркина) [15], с квадратичными базисными функ-

‘п+1 _ tn

+

‘'п+1 _ t

циями [16]. Дискретизация области решения выполнялась треугольными конечными элементами (рис. 2).

3.4. Определение точного положения границ фазового перехода

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

дТ Л(|) дТ (1) dhPs)

'рш дП_

1) Арт

-+

А®.

= L о® рт Эп+0 ртРрт dt

п_0 ^п+0

2) Т = Трш.

Записанные условия относятся к границе затвердевания в капле. Аналогичные свойства имеют место и на границе подплавления в подложке.

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

Для обеспечения требуемой точности мы использовали итерационную процедуру определения координат 2 каждого узла:

гкп +1 = гкп - f (гкп )

к к-1 _______________

-----———----------- п = 1 N

4-І кч /V к-К ’ п 1 ’

Ґ(2п ) - Ґ(2п )

где Ґ(г’к) = Т(гп, г*к) - Трт. В качестве начального приближения г0 использовалась координата фронта с предыдущего временного слоя, а в качестве г на первом шаге по времени использовались аналитические оценки (7), а на последующих шагах — оценка с помощью скорости фронта на предыдущем шаге по времени

И 00 - - И00 . ,

Ир п,1 Ир п,1 -1 , Л

■ +------------------------VI+1 -11 Л

^1-11 -1)

где п — номер точки в линии затвердевания; г — номер шага по времени.

3.5. Обработка особых точек

Определение положения фронта затвердевания при растекании капли на ее периферии имеет особенность. Она касается определения скорости затвердевания в точке А (рис. 4).

Дело в том, что при натекании расплава на холодную подложку на линии соприкосновения расплава с подложкой мгновенно устанавливается контактная температура (на рис. 4 областью соприкосновения расплава с подложкой является отрезок (%, +1)). В зависимости

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

Рис. 4. Особая точка при определении скорости роста фронта

ление его скорости как

dh

(®)

dt

г>) _ ГОО

р 1+1--— может быть

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

йр(8)^о) = С;^ ~р(8)(~) = С;

1

а°) t

рт

= Сг

D„

а®7

“рт1

р0

йЬ~_

dt

р0

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

dhp(s) = 1 Ч'Г 1 о а р (1 )

dt ч1+1 _ ч1 ч1 ро

V ~+1 _ ~ ар(1т) ^рт 1 а® ^рт

ч+1 _ч ^ ^>ир0 д/~+1 _ 1 DpUро

3.6. Тестирование вычислительной схемы

Для тестирования правильности построения вычислительных схем использовалась задача в следующей постановке: неподвижный цилиндрический столб жидкости высотой Н и радиусом R в начальный момент времени имеет температуру Тр0, при t> 0 поверхность подложки ^ = 0) поддерживается при нулевой температуре (рис. 5).

Распределение температуры подчиняется уравнению теплопроводности:

рС дТ = гігу(^гагі(Г)). дЧ

(16)

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

дТ

Т = Т-0- Щ = °.

(17)

При выполнении этих предположений границу z = = Н можно считать бесконечно удаленной. Кроме этого, в области решения должны выполняться следующие краевые условия:

Рис. 5. Схематическое изображение области решения тестовой задачи

- на линии фазового перехода г = hPs), hPs) > 0 Х(8) 0) _Э^ = ь л) т = т

рт рт ^ ^ртКрт ^ , 1 Трт’ (18)

_0

рт дп+0 рт"рт dt

- на линии z = 0

Т = 0, (19)

- на остальных границах задается однородное условие второго рода

дт

дп

- = 0.

(20)

Начальным условием задачи является Т = Тро, t = 0.

В такой постановке для задачи существует аналитическое решение (решение Неймана) [17]. Для высоты фронта затвердевания

№(1) = 2р(др^Г )1/2.

р V / к\~рт

Для температуры

Тр

Т = -

(21)

О(в)

(арт t)

1/2

при г < hPs)(t),

Т = Тро _

Т _Т

Тр0 Трт

G *[р(ар)/ арт)12]

2(арт t)12

при г > h.(s) (t),

(22)

(23)

где параметр в определяется из следующего уравнения

е_р2 С арт12(Тро _ Трт)ехр[_ арт в1/а

(1)]

рт

О(в) Є арт12Трт е*[в(арт/арт)12]

= Р^ртПІ

С (^т

СртТрт

где О (х) = 1 _ О(х).

чі+і _ ч

Рис. 6. Высота фронта затвердевания для задачи Неймана на всем протяжении тестирования (а), в начальные моменты времени (б)

Для определения точности решения были выбраны следующие критерии: ЕгИ = (Ир* - И®)/И^*, где

hPs») — точное значение текущей высоты затвердевшего слоя; ЕгТ = ||Т» - Т||2/||Т»||2 , где Т» — точное значение температуры.

Кроме того, использовался критерий, позволяющий контролировать консервативность используемого численного метода, Ег^) = Qe(t)/Qe(0), где

Se(t) = }}(^р С■^дГ - div( grad(T)) йО йт. (25)

Очевидно, что для точного решения Qe (т) = 0, а для приближенного величина отклонения от нуля будет характеризовать степень соблюдения законов сохранения, выраженных уравнением теплопроводности (16).

При проведении тестовых расчетов в качестве материала использовался алюминий с начальной температурой Тр0 = 1033 К.

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

Видно, что погрешность счета в начальные моменты времени максимальна и достаточно быстро убывает. Обоснование этому можно найти на рис. 7. Здесь ошибки в расчете высоты фронта затвердевания и температур имеют сходный характер, а величина ErQ, характеризующая сохранение тепловой энергии, на всем протяжении времени остается практически постоянной, близкой к нулю величиной. Очевидно, что в момент установления на нижней границе нулевой температуры там появляется тонкий слой расплава (или уже затвердевший слой), где градиент температуры очень велик. Ясно, что чем ближе мы будем подходить к моменту времени t = = 0, тем тоньше будет этот слой и тем больше градиент. Следует также отметить, что скорость изменения поля температур (а следовательно, и скорость роста фронта затвердевания) максимальна в начальный момент и с

течением времени убывает. При вычислениях это предъявляет определенные требования к шагу по времени и по пространству. С одной стороны, большая скорость изменения поля температур требует малого шага по времени; с другой стороны, чем меньше начальный шаг по времени, тем больше градиент температур и тем мельче должен быть шаг по пространству. Так как невозможно измельчать сетку до бесконечности, то для решения этой проблемы необходим разумный компромисс. Следует учесть возможности вычислительной техники и определить момент времени, с которого нужно иметь достаточно точное решение. Из рис. 7 видно, что полный баланс тепла в системе сохраняется хорошо. Он является стабилизирующим фактором и обеспечивает быстрое уменьшение ошибки в поле температур и высоте фронта, и начальная погрешность не оказывает существенного влияния на результаты последующих вычислений.

4. Результаты вычислительных экспериментов

Численные эксперименты проводились для различных сочетаний материалов капли и подложки, их основные теплофизические свойства приведены в табл. 3.

Рис. 7. Динамика изменения величин, характеризующих точность решения задачи Неймана

Л» •

л;.

М*’

6

6 -

4 -

6

Рис. 8. Сравнение вычисленного диаметра сплэта с аналитическими оценками (8), (9) (а) и аналитическими оценками (26) (б)

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

В качестве базовых материалов исследования быши вы-

браны никель для частицы и медь для подложки. При численном моделировании варьировались скорость частицы перед ударом (30, 75, 100 м/с), диаметр капли (0.025, 0.05, 0.1 мм), температура капли (1726, 2000,

2250, 2500, 2750, 3000, 3250, 3500 К). Температура подложки перед ударом оставалась постоянной и равнялась 293 К.

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

На рис. 8 представлены сравнение результатов вычислений диаметра сплэта DC с аналитическими оценками (8), (9) DS, и оценки диаметра и высоты сплэта согласно [18] DJ:

А =4 4/(3Яе14).

(26)

Рис. 9. Критериальное обобщение результатов выгаислений диаметра сплэта

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

Таблица 3

Теплофизические свойства материалов

Материал р®,кг/м3 р«кг/м3 с®, Дж/кг • К с?, Дж/кг • К Л®, Вт/м• К Л1», Вт/м• К Lm, Дж/кг а®, Н/м ц®,Па • с

А1 2 370 2 500 1 095 890 88 220 396 500 933 0.914 0.00125

Си 8 000 8 500 495 490 166 244 206150 1 356 1.28 0.0045

№ 7 790 8 902 640 430 72 67 303 000 1 726 1.62 0.00472

РЬ 10 686 11 058 142 146.6 31.6 15.5 23 000 600 0.468 0.00256

8п 6 980 7184 248 262 30 60.3 58 977 505 0.544 0.00185

8Ш 8 060 7 284 502 639 30 31 276 785 1 683 — —

Zn 6 575 6 878 481 452 54.5 90.7 111 000 692.6 0.782 0.00375

0.25 0.5 0.75 1 1.25 1.5 1.75 2 г

Рис. 10. Этапы формирования сплэта (Гр0 = 2500 К, ир0 = 30 м/с, Пр0 = 1 мм)

лее значительный разброс в определении диаметра сплэта дает сравнение с аналитическими оценками (8), (9).

На рис. 9 представлено критериальное обобщение результатов вычислений и аналитических оценок (8), (9) для диаметра сплэта.

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

О 0.25 0.5 0.75 1 1.25 1.5 Г

г ' 0.25 -

О 0.25 0.5 0.75 1 1.25 1.5 7

О 0.25 0.5 0.75 1 1.25 1.5 1.75 Г

Рис. 11. Этапы формирования сплэта (Гр0 = 3500 К, ир0 = 30 м/с, ,Ор0 = 0.5 мм)

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

При исследовании взаимодействия капли расплава с подложкой важны не только конечные параметры сис-

темы, но и ее эволюция. На рис. 10, 11 представлены этапы формирования сплэта при различных параметрах взаимодействия.

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

Таблица 4

Сравнение результатов вычислительного эксперимента при различных сочетаниях материалов капли и подложки и параметрах взаимодействия с результатами физического эксперимента

Частица Основа Тро,К Яо,мм и ро,м/С PeFo* Dъ Dc Ds

Sn Си 508 1.023 2.86 0.600 2.39 3.532 2.46 0.107

Sn Си 510 0.896 4.27 0.653 2.50 3.733 2.63 0.094

Sn Си 535 0.897 7.80 0.813 2.94 4.223 3.12 0.066

Sn Си 542 0.62 6.00 0.665 2.77 3.708 2.70 0.089

Sn Sus 544 0.49 8.56 0.856 3.24 4.221 3.39 0.056

Sn Sus 544 0.458 9.04 0.851 3.15 4.215 3.38 0.057

Sn Sus 544 0.462 10.80 0.899 3.49 4.351 3.53 0.052

Sn Sus 657 1.034 4.28 1.047 3.00 4.538 3.26 0.061

РЬ Си 625 0.566 6.60 0.594 2.69 3.531 2.40 0.112

РЬ Си 645 0.399 13.10 0.680 2.79 3.783 2.64 0.092

РЬ Си 657 0.315 17.10 0.695 2.84 3.834 2.69 0.089

РЬ Си 655 0.381 12.60 0.670 2.75 3.754 2.64 0.092

РЬ Си 710 0.744 6.11 0.703 2.93 3.826 2.72 0.087

РЬ SUS 625 0.54 10.20 0.744 3.01 3.970 2.86 0.079

РЬ SUS 630 0.599 6,22 0.659 2.80 3.706 2.61 0.094

РЬ SUS 645 0.553 6.60 0.661 2.70 3.717 2.64 0.093

РЬ SUS 665 0.401 12.70 0.765 3.11 4.016 2.96 0.073

РЬ SUS 710 0.547 4.93 0.679 2.90 3.704 2.75 0.085

Zn SUS 690 1.29 2.41 0.687 2.26 3.606 2.27 0.123

Zn SUS 700 0.985 6.17 0.854 3.39 4.209 3.26 0.060

Zn SUS 783 3.47 3.69 1.165 4.44 4.985 4.64 0.030

Zn SUS 784 3.59 3.68 1.217 4.67 5.099 4.69 0.029

Zn Sn 783 3.14 3.68 1.199 3.61 5.125 5.46 0.021

4 -

2 -

У*

4 -

2 -

#

Г

Рис. 12. Сравнение результатов непосредственных измерений диаметра сплэта с аналитическими оценками (8), (9) (а) и с результатами численного эксперимента (б)

D

5 -4 -

3 -

I а

О % ЛА, А

Ъ

AiAi

4%

i

*

л с + j

♦ Е

О S

D

5 Н 4

3 -

2 -

♦

О

А ^ АА *

А

д с

+ j

♦ Е

❖ S

500

1000

1500

2000

We

10000

20000

Re

Рис. 13. Зависимости диаметра сплэта от чисел Вебера (а) и Рейнольдса (б) для результатов физического эксперимента и различных методов расчета

раметрах взаимодействия, для которых имеются данные прямого физического эксперимента. Сравнение расчетов с экспериментальными данными является основным критерием оценки адекватности используемой физической модели и вычислительной схемы. Данные расчетов и непосредственных измерений результатов физического эксперимента 3Е [19] приведены в табл. 4.

На рис. 12 представлено сравнение результатов непосредственных измерений диаметра сплэта с аналитическими оценками (8), (9) и результатами численного эксперимента.

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

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

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

5. Заключение

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

существенное влияние на результат. Упрощенная модель эквивалентного цилиндра в своих разновидностях не дает удовлетворительных оценок параметров системы сплэт - основа.

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 98-02-17810), а также Сибирского отделения РАН (Междисциплинарная интеграционная программа на 2000-2002 гг., проект № 45).

Литература

1. Madejski J. Solidification of droplets on cold surface // Int. J. of Heat Mass Transfer. - 1976. - V. 19. - P. 1009-1013.

2. Rangel R.H., Bian X. Metal-droplet deposition model including liquid deformation and substrate remelting // Int. J. Heat Mass Transfer. -1997. - V. 40. - No. 11. - P. 2549-2564.

3. Zhao Z., Poulikakos D. Heat transfer and fluid dynamics during the collision of a liquid drop on a substrate // Int. J. Heat Mass Transfer. -1996. - V. 39. - No. 13. - P. 2771-2789.

4. Fukai J., Zhao Z., Poulikakos D., Megaridis C., Mityatake O. Modeling of the deformation of a liquid droplet impinging upon a flat surface // Physics Fluids. - 1993. - V. A5. - P. 2588-2599.

5. Pasandideh-Fard M., Mostaghimi J. On the spreading and solidification of molten particles in a plasma spray process: effect of thermal contact resistance // Plasma Chemistry and Plasma Processing. -1996. - V. 16 (Supplement). - P. 83S-97S.

6. Liu H., Lavernia E.J., Rangel R.H. Numerical simulation of impingement of molten Ti, Ni, and W droplets on a flat substrate // J. Thermal Spray Technology. - 1993. - V. 2. - No. 4. - P. 369-378.

7. Pasandideh-Fard M., Bhola R., Chandra S., Mostaghimi J. Deposition of tin droplets on a steel plate: simulations and experiments // Int. J. Heat Mass Transfer. - 1998. - V. 41. - P. 2929-2945.

8. Mostaghimi J. Modeling droplet impact in plasma spray processes // Pure and Appl. Chem. - 1998. - V. 70. - No. 6. - P. 1209-1215.

9. Bertagnolli M., Marchese M., Jacucci G., Doltsinis I. St., Noelting S. Thermomechanical simulation of the splashing of ceramic droplets on a rigid substrate // J. Computational physics. - 1997. - V. 133. -P. 205-221.

10. Солоненко О.П., Смирнов А.В., Клименов В.А., Бутов В.Г., Иванов Ю.Ф. Роль границ раздела при формировании сплэтов и структуры покрытий // Физ. мезомех. - 1999. - № 1-2. - C. 123-140.

11. Solonenko O.P., Smirnov A.V Generalized map of the plasma sprayed splats formation // Proc. of 3rd Asia-Pacific Conference on Plasma Science and Technology, Tokyo, Japan 15-17 July, 1996. - 1996. -P. 247-252.

12. Solonenko O.P., Smirnov A.V Comparative analysis and testing of different theories characterizing a diameter and thickness of plasma sprayed splats // Proc. of 12th Intern. Symp. on Plasma Chemistry. -Mineapolis, USA, 1995. - P. 921-926.

13. Жуков М.Ф., Солоненко О.П., Федорченко А.И. Равновесная кристаллизация расплавленных частиц на поверхности при плаз-

менном напылении // Доклады Академии наук СССР. - 1990. -Т. 314. - № 2. - C. 369-374.

14. Markworth A.J., Saunders J.H. An improved velocity field for the Madejski splat - quench solidification model // Int. J. Heat Mass Transfer. - 1992. - V. 35. - No. 7. - P. 1836-1837.

15. Съярле Ф. Метод конечных элементов для эллиптических задач. -М.: Мир, 1980. - 512 с.

16. СегерлиндД. Применение метода конечных элементов. - М.: Мир, 1979. - 392 с.

17. Карслоу Г., Егер Д. Теплопроводность твердых тел / Пер. с англ. -М.: Наука, 1964. - 487 с.

18. Bennet T., Poulikakos D. Splat - quench solidification: estimating the maximum spreading of a droplet impacting a solid surface // J. Materials Science. - 1993. - V. 28. - P. 963-970.

19. Смирнов А.В. Экспериментальное исследование взаимодействия капель металлических расплавов с основой / Дис. ... канд. физ.-мат. наук. - Новосибирск: ИТПМ СО РАН, 2000. - 203 с.