Моделирование температурного поля поверхности при электроискровом легировании металлов

Полный текст

Введение

Проблема улучшения эксплуатационных свойств (износо- и жаростойкости, а также коррозионной стойкости) рабочих поверхностей деталей машин, инструментов и технологической оснастки путем улучшения физико-химико-механических характеристик приобретает все большую актуальность. В настоящее время прогресс в машиностроении во многом связан с применением высокоэффективных методов модификации рабочих поверхностей деталей машин, инструментов и технологической оснастки, основанных на использовании потоков энергии с удельной мощностью в пятне нагрева более 10² Вт/мм². К их числу относится современный наукоемкий метод электроискрового легирования (ЭИЛ). Его достоинствами являются высокая прочность сцепления легированного слоя из любых токопроводящих материалов (в том числе тугоплавких металлов и сплавов) с обрабатываемым материалом, низкая энергоемкость процесса, простота выполнения технологической операции и др. Применение метода ЭИЛ для упрочнения рабочих поверхностей деталей, инструмента и оснастки обеспечивает повышение срока их службы в пять и более раз. В условиях современного машиностроительного производства метод ЭИЛ является востребованным, а его изучение с целью более эффективного использования при получении функциональных покрытий становится актуальным.

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

При каждом цикле ЭИЛ на катоде образуется лунка, заполненная материалом, полученным в результате взаимодействия катода, анода и межэлектродной среды [1; 2].

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

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

Обзор литературы

Многие российские и зарубежные исследования посвящены использованию метода ЭИЛ при создании упрочняющих покрытий на металлах и сплавах. Представлены результаты повышения физических, технологических и эксплуатационных характеристик покрытий, нанесенных методом ЭИЛ на титан и его сплавы [1; 2]; виден положительный эффект применения метода ЭИЛ на твердых сплавах [3; 4]. Кроме того, исследователи начинают проявлять интерес к ЭИЛ легких алюминиевых сплавов. Так, учеными [5] представлены характеристики микроструктуры и кавитационной эрозии покрытия из сплава Al–Si. Однако наиболее широкое распространение метод ЭИЛ получил при обработке поверхностей сталей [6–8]. В названных работах исследователи добились повышения термостойкости многослойного покрытия. Результаты зарубежных исследований свидетельствуют также об уменьшении коррозионной активности покрытия, полученного на нержавеющей стали методом электроискрового осаждения в расплавленном цинке [7].

В настоящее время проводятся исследования по разработке критериев эффективности [9; 10] и электрофизических моделей процесса ЭИЛ [11], которые позволяют найти зависимость критериев качества покрытия от технологических параметров процесса. Разработанные критерии метода ЭИЛ позволили выделить его в отдельный раздел материаловедения [12].

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

При выборе обрабатывающего материала для получения на обрабатываемом материале покрытий с заданными функциональными свойствами возникают сложности определения температурного поля в поверхностном слое при реализации метода¹ [13–15].

В работах зарубежных [16–19] и отечественных исследователей [20; 21] рассмотрены некоторые математические модели определения температурного поля в поверхностном слое катода в процессе ЭИЛ. Однако описанные модели не учитывают ряд факторов и сложны в реализации.

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

Материалы и методы

Рассмотрена задача нагрева катода, имеющего форму прямоугольного параллелепипеда (Q₀ = {(‒a,a)×(‒b,b)×(0,c)}), в декартовой системе координат при ЭИЛ (рис. 1). Здесь P₀ = $P_0=\left\{z=\mathrm{0,}x\in (-a,a),y\in (-b,b)\right\}$  – рабочая поверхность; q₀ = q₀{x,y,t} = q₀{x',t} – поверхностный тепловой поток; t – время; N – количество капель [20].

 

 

 

Рис. 1. Катод (прямоугольный параллелепипед Q₀) с каплями на его поверхности P₀;
a, b, c – размеры катода в декартовой системе координат XYZ

Fig. 1. The cathode (rectangular parallelepiped Q₀ ) with drops on its surface P₀;
a, b, c – the size of the cathode in the Cartesian coordinate system XYZ

 

 

Капли Q_i, равномерно заполняя некоторую часть грани {z = 0} (рис. 1), за один искровой разряд отдают посредством теплопроводности и излучения тепловой поток:

$q_{1i}={\alpha }_i(T_i-T_{sr})+\kappa \sigma T_i^4,i=\overline{\mathrm{1,}N},$

где α_i – коэффициент теплоотдачи капель; T_i – температура i-ой капли; T_sr – температура окружающей среды; k = λ/(cρ) – коэффициент температуропроводности; σ – универсальная постоянная Стефана-Больцмана [Там же].

В математическом плане исследование теплового процесса (нагрева параллелепипеда и остывания капель Q_i за некоторый фиксированный промежуток времени t^* от начала воздействия искрового разряда t₀) приводит к необходимости рассмотрения следующей начально-краевой задачи нелинейной зависимости:

${\lambda }_p\Delta T_0-c_p{\rho }_p\frac{\partial T_0}{\partial t}=0$ ,$\overrightarrow{x}\in Q_0$  ,$t_0<t<t*$  ,

${\lambda }_k\Delta T_i-c_k{\rho }_k\frac{\partial T_i}{\partial t}=0$ , $\overrightarrow{x}\in Q_i$ ,$t_0<t<t*$  ,

$T_0\left|{}_{t=0}=\right.{\Phi }_0(\overrightarrow{x})$ ,$\overrightarrow{x}\in Q_0$  ,$T_i\left|{}_{t=0}=\right.{\Phi }_i(\overrightarrow{x})$  ,

$\overrightarrow{x}\in Q_i$ ,$i=\overline{\mathrm{1,}N}$  ,              (1)

$\frac{\partial T_0}{\partial x}\left|{}_{x=\pm a}\right.=0$ , $\frac{\partial T_0}{\partial y}\left|{}_{y=\pm b}\right.=0$ ,$\frac{\partial T_0}{\partial z}\left|{}_{z=c}\right.=0$  ,

${\lambda }_p\frac{\partial T_0}{\partial z}+{\alpha }_p(T_0-T_{sr})=0$ ,$z=0$  ,

$(x,y)\notin {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\Omega }_i}$ ,

$T_0=T_i$ ,${\lambda }_p\frac{\partial T_0}{\partial z}={\lambda }_k\frac{\partial T_i}{\partial z}$  , $z=0$ ,

$(x,y)\in {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\overline{\Omega }}_i}$ ,

${\lambda }_k\frac{\partial T_i}{\partial \overrightarrow{n}}+{\alpha }_k(T_i-T_{ñð/sr})+\kappa \sigma T_i^4=q_0,$ ,$\overrightarrow{x}\in \partial Q_i$  ,$z<0$ 

где$\Delta =\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}$  – оператор Лапласа; T₀ – температура катода; T_sr – температура окружающей среды; Φ₀(x)  ${\Phi }_i(\overrightarrow{x})$ – известные заданные функции; λ_p, λ_k, c_p, c_k, ρ_p, ρ_k, α_p, α_k – коэффициенты теплопроводности, удельной теплоемкости, удельной плотности и теплоотдачи параллелепипеда (p) и капель (k) соответственно; $\overrightarrow{n}$ – внешняя нормаль к $\partial {\Omega }_i$  ; здесь и далее коэффициенты с индексом i = 0 относятся к катоду Q₀, с индексами $i\ne 0$ – к i-ой капле Q_i; ${\Omega }_i=Q_i{\displaystyle \cap \{z=0\}}$, ${\overline{\Omega }}_i={\Omega }_i{\displaystyle \cup \partial {\Omega }_i}$   [Там же].

Ввиду сложности задачи (1) ее решение требует перехода к новым переменным: U = T₀ ‒ T_sr; U_i = T_i ‒ T_sr;  $a_0^2=a_{0i}^2={\lambda }_i/(c_i{\rho }_i)$ . Тогда задача для параллелепипеда Q₀ получит следующий вид:

$a_0^2\Delta U-\frac{\partial U}{\partial t}=0$ ,$\overrightarrow{x}\in Q_i$  , $i=\overline{\mathrm{0,}N},$ 

$t_0<t<t^{*}$ ,

$U\left|{}_{t=t_0}=\right.\varphi (\overrightarrow{x})$ ,$\overrightarrow{x}\in Q_0{\displaystyle \cup \left({\displaystyle \underset{i=1}{\overset{N}{\cup }}{Q}_i}\right)}$  ,

$\frac{\partial U}{\partial x}\left|{}_{x=\pm a}\right.=0$ ,$\frac{\partial U}{\partial y}\left|{}_{y=\pm b}\right.=0$  ,$\frac{\partial U}{\partial z}\left|{}_{z=c}\right.=0$  , (2)

${\lambda }_p\frac{\partial U}{\partial z}+{\alpha }_{p}U=0$ ,$z=0$  ,$(x,y)\notin {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\Omega }_i}$  ,

${\lambda }_k\frac{\partial U_i}{\partial \overrightarrow{n}}+{\alpha }_k{U}_i+\kappa \sigma {(U_i+T_{sr})}^4=q_0,$

$z=0$ ,$(x,y)\in {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\overline{\Omega }}_i}$  ,

где  $\varphi (\overrightarrow{x})=\left\{\begin{array}{l}{\varphi }_0(\overrightarrow{x}),\overrightarrow{x}\in Q_0,\\ {\varphi }_i(\overrightarrow{x}),z=\mathrm{0,}(x,y)\in {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\Omega }_i,}\end{array}\right.$${\varphi }_0(\overrightarrow{x})={\Phi }_0(\overrightarrow{x})-T_{sr}$ , ${\varphi }_i(\overrightarrow{x})={\Phi }_i(\overrightarrow{x})-T_{sr}.$  [Там же].

Предположим, что размеры капель малы; тогда температура внутри объема капли будет почти постоянна, а существенно изменяться она станет только во времени t.

Будем считать, что при каждом t > 0 справедливы приближенные равенства, которые тем более верны, чем меньше размеры капли [Там же]:

$U(\overrightarrow{x},t)\approx const,$

$\partial U(\overrightarrow{x},t)/\partial t\approx const,$ $\overrightarrow{x}\in Q_i$  ,$i=\overline{\mathrm{1,}N}$ 

${\left.\partial U(\overrightarrow{x},t)/\partial z\right|}_{z=-0}\approx const,$ ,

$x^{\prime }\in P_i$ ,$i=\overline{\mathrm{1,}N}$  ,              (3)

где P_i – основание i-ой капли.

Это позволяет усреднить задачу по объемам Q_i, ограничившись только областью Q₀.

Так как для нашего исследования важна не конкретная геометрическая форма капли, а лишь ее размеры, для простоты вычислений предположим, что капля заключена в параллелепипед определенного размера (рис. 2) [Там же]:

$Q_i=\left\{(a_{1i},a_{2i})\times (b_{1i},b_{2i})\times (\mathrm{0,}d)\right\}$ ,

$\Delta a_i=a_{2i}-a_{1i}$ ,$\Delta b_i=b_{2i}-b_{1i}$  ;

тогда $\left|Q\right|=\Delta a_i\Delta b_{i}d$  ,$\left|P\right|=\Delta a_i\Delta b_i$ , $\left|S\right|=2\Delta a_{i}d+2\Delta b_{i}d+\Delta a{}_i\Delta b_i=2d(\Delta a_i+\Delta b_i)+\left|P\right|.$ 

 

 

 

Рис. 2. Расположение капли на катоде Q₀ в декартовой системе координат XYZ

 

Fig. 2. The location of a drop on the cathode Q₀ in the Cartesian coordinate system XYZ

 

 

Результатом математических преобразований становится следующая задача для параллелепипеда Q₀ [Там же]:

$a_0^2\Delta U-\frac{\partial U}{\partial t}=0$ ,$\overrightarrow{x}\in Q_0$  ,$t_0<t<t^{*}$  ,

$U\left|{}_{t=t_0}=\right.\varphi (\overrightarrow{x})$ ,$\overrightarrow{x}\in Q_i$  ,$i=\overline{\mathrm{0,}N},$ 

$\frac{\partial U}{\partial x}\left|{}_{x=\pm a}\right.=0$ ,$\frac{\partial U}{\partial y}\left|{}_{y=\pm b}\right.=0$  ,$\frac{\partial U}{\partial z}\left|{}_{z=c}\right.=0$  , (4)

$\frac{\partial U}{\partial z}+h_{0}U=0$ ,$z=0$  ,$(x,y)\notin {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\Omega }_i}$  ,

$\frac{\partial U}{\partial z}+h_{0}U=q_0-\left[c_0\frac{\partial U}{\partial t}+c_{1}U+c_{2}f(U)\right],$

$z=0$ ,$(x,y)\in {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\overline{\Omega }}_i}$  .

Данная задача определена только в области Q₀, так как граничное условие для температуры U заменяет все уравнения теплового баланса в каплях.

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

Для определения температуры катода в форме параллелепипеда Q₀ решение задачи (4) состоит из двух этапов [Там же].

Первый этап

Пусть U⁰ – решение задачи (5). С помощью найденного решения определим, что

$\psi (t)=q_0-\left[c_0\frac{\partial U^0}{\partial t}+c_1{U}^0+c_{2}f(U^0)\right]$ ;

$(x,y)\in {\Omega }_i$ .                  (5)

Второй этап

Сложный нелинейный теплообмен между каплями Q_i, $i=\overline{\mathrm{1,}N}$  и катодом  Q₀ в граничном условии задачи (4) заменяется линейным:

$\frac{\partial U}{\partial z}+h_{0}U=\psi (t)$ ,$(x,y)\in {\displaystyle \underset{i=1}{\overset{N}{\cup }}{\overline{\Omega }}_i},t>0.$ 

Для решения задачи применим вторую формулу Грина и функцию Грина $G(\overrightarrow{x},\overrightarrow{\xi },t-\tau ),$ которая является решением соответствующей задачи² [20], где $\overrightarrow{\xi }=(\xi ,\eta ,\zeta ),$ t, τ > 0. Тогда параллелепипед можно рассматривать как полупространство для капли, подставив в формулу (5) а = b = с = + ∞.

Наиболее простым для численной реализации является случай полупространства $\left\{z>\mathrm{0,}{x}^{\prime }\in R^2\right\}$ в качестве параллелепипеда Q₀, на границу {z = 0} которого нанесена единственная капля {N = 1}, и значения коэффициента теплоотдачи α_п = 0; следовательно, h₀ = 0, то есть учитывается только теплообмен между каплей и катодом.

Вторая формула Грина и данная функция Грина приводят задачу (4) к эквивалентному ей нелинейному интегральному уравнению следующего вида [20]:

$U(\overrightarrow{x},t)=U{}_{\wedge }(\overrightarrow{x},t,t^{\text{'}})-a_0^2\underset{t^{\text{'}}}{\overset{t}{\int }}d\tau \underset{0}{\overset{a^{\text{'}}}{\int }}\underset{0}{\overset{b^{\text{'}}}{\int }}G(\overrightarrow{x},\eta ,\zeta \mathrm{,0,}t-\tau )Q\left(U\right)d\eta d\zeta$ (6)

где $a^{\prime }=\Delta a/2,$ $b^{\prime }=\Delta b/2$  – половинные длины сторон капли;

$U_{\wedge }(\overrightarrow{x},t)=\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}G(\overrightarrow{x},\overrightarrow{\xi },t-t^{\text{'}})U(\overrightarrow{\xi },t^{\text{'}})d\overrightarrow{\xi }+a_0^2\underset{t^{\text{'}}}{\overset{t}{\int }}d\tau \underset{0}{\overset{a}{\int }}\underset{0}{\overset{b}{\int }}G(\overrightarrow{x},\xi ,\eta \mathrm{,0,}t-\tau )q_0(\xi ,\eta \mathrm{,0,}\tau )d\xi d\eta$

$Q\left(U\right)=c_0\frac{\partial U(\xi ,\eta \mathrm{,0,}\tau )}{\partial \tau }+c_{1}U(\xi ,\eta \mathrm{,0,}\tau )+c_{2}f\left(U\right(\xi ,\eta \mathrm{,0,}\tau \left)\right)$

 

Интегральные соотношения (6), как и задача (5), приближенно описывают процесс передачи тепла катоду каплей, помещенной на его границу.

В подынтегральном выражении уравнения (6) функция U = U(ξ,η,0) рассматривается в области, в которой она почти постоянна по ξ,η. Поэтому, подставив в уравнение (6) z = 0 и проинтегрировав его [Там же], получим нелинейное интегральное уравнение типа Вольтера:

$\Delta a\Delta bU\left(t\right)=U_{\wedge }\left(t\right)-a_0^2\underset{t^{\text{'}}}{\overset{t}{\int }}G(t-\tau )Q\left(U\right(\tau \left)\right)d\tau$         (7)

относительно усредненной функции

$U\left(t\right)=\frac{1}{\Delta a\Delta b}\underset{0}{\overset{a}{\int }}\underset{0}{\overset{b}{\int }}U(x,y\mathrm{,0,}t)dxdy$

Здесь

$G(t-\tau )=\underset{0}{\overset{a}{\int }}\underset{0}{\overset{b}{\int }}\underset{0}{\overset{a}{\int }}{\left.\underset{0}{\overset{b}{\int }}G\right|}_{\begin{array}{l}z=0\\ \zeta =0\end{array}}dxdyd\xi d\eta =+\frac{1}{\sqrt{\pi a_0(t-\tau })}\left[a\Phi \left(\frac{a}{\sqrt{a_0(t-\tau )}}\right)+\sqrt{\frac{t-\tau }{\pi }}\left(e^{-\frac{a^2}{a_0(t-\tau )}}-1\right)\right]\times$

$\times \left[b\Phi \left(\frac{b}{\sqrt{a_0(t-\tau )}}\right)+\sqrt{\frac{a_0(t-\tau )}{\pi }}\left(e^{-\frac{b^2}{a_0(t-\tau )}}-1\right)\right]$

 

$U_{\wedge }(x,y,z,t)=\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{a}{\int }}\underset{0}{\overset{b}{\int }}{\left.G\right|}_{\tau =0}d\xi d\eta dzdxdy$

 

где $\Phi (x)=\frac{2}{\pi }{\displaystyle \underset{0}{\overset{x}{\int }}{e}^{-t^2}}dt$  – интеграл вероятности³.

Для численного решения задачи (6) введем пространственную и временную сетку (x_i, y_j, z_k), $i=\overline{\mathrm{1,}{N}_1},$ $j=\overline{\mathrm{1,}{N}_2},$ $k=\overline{\mathrm{1,}{N}_3},$ $\{0=t_0<t_1<\mathrm{...}<t_{N_0}=t*\}$   .

Тогда $U_n^{ijk}=(x_i,y_j,z_k,t_n)$  – температура в узле (x_i, y_j, z_k) в момент времени t = t_n. Производную по времени заменяем разностным соотношением $\frac{\partial U}{\partial t}=\frac{U_n-U_{n-1}}{\Delta {\tau }_n},$  где  $U_n=U\left|{}_{t=t_n}\right.,$ $\Delta {\tau }_n=t_n-t_{n-1}.$ 

В уравнении (7) заменяем в подынтегральном выражении U(t) на $\frac{U_n+U_{n-1}}{2}$  и получаем:

 $\Delta a\Delta b{U}_n=U_{\wedge n}-G_n\left[c_0\frac{U_n-U_{n-1}}{\Delta \tau }+c_1\frac{U_n+U_{n-1}}{2}\right.\left.+c_{2}f\left(\frac{U_n+U_{n-1}}{2}\right)\right],$  (8)

где: $G_n=\underset{t_{n-1}}{\overset{t_n}{\int }}G(t_n-\tau )d\tau$ 

Решение нелинейного уравнения (5) описывает температурное поле катода под основанием капли в момент t = t_n. Поле внутри параллелепипеда определяется простой квадратурой (6), имеющей следующий вид [20]:

$U(x,y,z,t_n)=U_{\wedge }(x,y,z,t_n)-a_0^2\underset{t_{n-1}}{\overset{t_n}{\int }}G(x,y,z,t_n-\tau )Q_n\left(U_n\right)d\tau$   (9)

где

$U_{\wedge }(x,y,z,t_n)=\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}{\left.G\right|}_{\begin{array}{l}\\ \tau =t_{n-1}\\ t=t_n\end{array}}{U}_{n-1}(\xi ,\eta ,\zeta )d\xi d\eta d\zeta +\underset{t_{n-1}}{\overset{t_n}{\int }}G(x,y,z,t_n-\tau )q_0\left(\tau \right)d\tau$             (10)

$G(x,y,z,t-\tau )={\displaystyle \underset{0}{\overset{a}{\int }}{\displaystyle \underset{0}{\overset{b}{\int }}{\left.G\right|}_{\zeta =0}}}d\xi d\eta =\frac{\exp \left\{-\frac{z^2}{4a_0(t-\tau )}\right\}}{4\sqrt{\pi a_0(t-\tau })}\left[\Phi \left(\frac{x+a}{2\sqrt{a_0(t-\tau )}}\right)-\Phi \left(\frac{x-a}{2\sqrt{a_0(t-\tau )}}\right)\right]\times$

$\times \left[\Phi \left(\frac{y+b}{2\sqrt{a_0(t-\tau )}}\right)-\Phi \left(\frac{y-b}{2\sqrt{a_0(t-\tau )}}\right)\right]$

 

$Q_n\left(U_n\right)=c_0\frac{U_n-U_{n-1}}{\Delta \tau }+c_1\frac{U_n+U_{n-1}}{2}+c_{2}f\left(\frac{U_n+U_{n-1}}{2}\right)$

 

Для численных расчетов интеграл (10) заменим на тройную сумму следующего вида:

$\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}{\left.G\right|}_{\tau =t_{n-1}}{U}_{n-1}(\xi ,\eta ,\zeta )d\xi d\eta d\zeta =\sum _{i_1}^{N_1}\sum _{j_1}^{N_2}\sum _{k_1}^{N_3}{G}_{n{i}_1{j}_1{k}_1}^{ijk}{\overline{U}}_{n-1}^{i_1{j}_1{k}_1}$

где  $G_{n{i}_1{j}_1{k}_1}^{ijk}=\underset{x_{i_1-1}}{\overset{x_{i_1}}{\int }}\underset{y_{j_1-1}}{\overset{y_{j_1}}{\int }}\underset{z_{k_1-1}}{\overset{z_{k_1}}{\int }}G(x_i,y_j,z_k,\Delta t)d\xi d\eta d\zeta$ 

Среднеарифметическое значение функции относительно восьми узлов параллелепипеда  $(x_{i-1},x_i)\times (y_{j-1},y_j)\times (z_{k-1},z_k),$ $U_n^{ijk}=U(x_i,y_j,z_k,t_n),$ $i=\overline{\mathrm{1,}{N}_1},$ $j=\overline{\mathrm{1,}{N}_2},$ $k=\overline{\mathrm{1,}{N}_3},$ определяется по формуле:

$\begin{array}{c}{\overline{U}}_{n-1}^{ijk}=\frac{U_{n-1}^{i-\mathrm{1,}j-\mathrm{1,}k-1}+U_{n-1}^{i-\mathrm{1,}-\mathrm{1,}k}+U_{n-1}^{i-\mathrm{1,}j,k-1}+}{8}+\\ \frac{+U_{n-1}^{i-\mathrm{1,}j,k}+U_{n-1}^{i,j-\mathrm{1,}k-1}+U_{n-1}^{i,j-\mathrm{1,}k}+U_{n-1}^{i,j,k-1}+U_{n-1}^{i,j,k}}{8}\end{array}$  (11)

Так как решение $U_{n-1}(\xi ,\eta ,\zeta )$  является убывающим от нуля до бесконечности по переменным ξ,η,ζ, очевидно, что узловые точки ${\left\{x_i\right\}}_{i=1}^{N_1},$ ${\left\{y_j\right\}}_{j=1}^{N_2},$ ${\left\{z_k\right\}}_{k=1}^{N_3}$ следует уплотнить вблизи начала координат. С учетом вышеизложенного получаем окончательные выражения для определения температурного поля в узлах x_i, y_j, z_k параллелепипеда [Там же]:

$U_n^{ijk}=U_{\wedge n}^{ijk}-a_0^2{Q}_n\underset{t_{n-1}}{\overset{t_n}{\int }}G(x_i,y_j,z_k,t_n-\tau )d\tau$

 $U_{\wedge n}^{ijk}={\displaystyle \sum _{i_1=1}^{N_1}\sum _{j_1=1}^{N_2}\sum _{k_1=1}^{N_3}{G}_{n{i}_1{j}_1{k}_1}^{ijk}{U}_{n-1}^{i_1{j}_1{k}_1}+}{q}_0{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}{G}_n}(x_i,y_j,z_k,t_n-\tau )d\tau$    (12)

После этого получим все данные для перехода к следующему временному отрезку [t_n, t_{n + 1}].

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

$\Delta a_i\Delta b_i{U}_i\left(t\right)=U_{\wedge i}\left(t\right)-a_0^2\underset{t_{n-1}}{\overset{t_n}{\int }}\sum _{j=1}^N{G}_{ij}(t-\tau )Q\left(U_j\right(\tau \left)\right)d\tau$     (13)

где $G_{ij}(t-\tau )={\displaystyle \underset{a_{1i}}{\overset{a_{2i}}{\int }}{\displaystyle \underset{b_{1i}}{\overset{b_{2i}}{\int }}{\displaystyle \underset{a_{1i}}{\overset{a_{2i}}{\int }}{\left.{\displaystyle \underset{b_{1i}}{\overset{b_{2i}}{\int }}G}\right|}_{\begin{array}{l}z=0\\ \zeta =0\end{array}}}dx}}dyd\xi d\eta ,$  a_1i,a_2i,b_1i,b_2i – координаты i-ой капли; N = 4 (в данном случае).

Дискретизация по времени t преобразует систему (13) в систему нелинейных алгебраических уравнений:

 $\Delta a_i\Delta b_i{U}_i=U_{\wedge i}-a_0^2\sum _{j=1}^N{G}_{ij}Q\left(U_j\right)$    (14)

При этом поле внутри параллелепипеда определяется так:

$U(x,y,z,t)=U_{\wedge }(x,y,z,t)-a_0^2\sum _{i=1}^N\underset{t_{n-1}}{\overset{t_n}{\int }}\underset{a_{1i}}{\overset{a_{2i}}{\int }}\underset{b_{1i}}{\overset{b_{2i}}{\int }}{\left.G\right|}_{\begin{array}{l}\\ \zeta =0\end{array}}{\left.Q\left(U\right)\right|}_{\begin{array}{l}\\ \zeta =0\end{array}}d\xi d\eta d\tau$  (15)

Таким образом, определив поток Q(U) из одной капли в параллелепипед, можем начать рассмотрение упорядоченного (периодического по переменным х, у) множества капель [Там же].

Однако решение системы нелинейных уравнений вида (14) с последующим определением поля в области Q по формуле вида (15) становится трудоемким. Поэтому для решения задачи нелинейный поток из капель в параллелепипед заменяется линейным, что является более предпочтительным [Там же].

По представленному алгоритму создан пакет программ и проведены численные расчеты определения температурного поля параллелепипеда⁴.

С учетом функций в выражениях (7) и (9) для численного интегрирования применялись квадратурные формулы⁵с весовой функцией $1/\sqrt{t}$ . Пробные вычисления по разработанному алгоритму показали эффективность выбранного метода.

Расчеты проводились с учетом следующих данных [Там же].

1. Размеры капли: $\Delta a_i/2$ = 10^–2 см; $\Delta b_i/2$ = 2·10^‒2 см; d = 10^–2 см. Здесь размеры капли приняты неодинаковыми в ее основании, что соответствует более общей задаче – обработке посредством ЭИЛ не плоской (Δа = Δb), а криволинейной поверхности.
2. Теплофизические константы: T_cр = = 20 °C; T_o = 1 400 °С (начальная температура капли); q₀ = 0.
3. Материал капли – вольфрам (W):
   λ_k= 1,73 Вт/cм ∙ °C; ρ_к= 19,25 г/cм³; с_к = 0,15 Дж/г ∙ К; α_к = 0,001 Вт / (см² ∙ °С).
4. Материал катода – железо (Fe): λ_k= 0,733 Вт/cм ∙ °C; ρ_к= 7,87 г/cм³; с_к = 0,46 Дж/г ∙ К; α_к = 0,001 Вт/(см² ∙ °С).
5. Узлы пространственной сетки, см:

 $\left\{x_{i=1}^6\right\}$ = 0; 0,005; 0,01; 0,011; 0,013; 0,017;

 $\left\{y_{i=1}^6\right\}$ = 0; 0,01; 0,02; 0,022; 0,026; 0,034;

 $\left\{z_{i=1}^6\right\}$ = 0; 0,001; 0,003; 0,006; 0,01; 0,015.

Динамика температурного поля некоторых точек поверхности представлена в таблице. В круглых скобках даны координаты точек, обозначенных на рис. 3. Числа n и m в круглых скобках указывают на номер координаты вектора $\overrightarrow{x}=\left\{x_i\right\}$  и вектора $\overrightarrow{y}=\left\{y_i\right\}$  (например, точка (3, 1) таблицы соответствует точке (х₃, у₁) на плоскости x0y) [Там же].

 

Таблица Динамика температурного поля точек на поверхности в зависимости от времени

Table The dynamics of the temperature field of points on the surface depending on time

 

 

Время / Time	Координаты точек / The coordinates of the points
	(1,1)	(2,1)	(2,2)	(3,1)	(3,2)	(3,3)	(4,1)	(4,2)	(4,3)	(4,4)
t = 1 ∙ 10–5 c / t = 1 ∙ 10–5 s	1 078	1 043	1 042	541	540	270	294	293	147	65
t = 2 ∙ 10–5 с / t = 2 ∙ 10–5 s	904	835	790	467	439	228	323	299	156	96
t = 3 ∙ 10–5 c / t = 3 ∙ 10–5 s	775	712	660	419	386	209	310	281	155	105
t = 4 ∙ 10–5 c / t = 4 ∙ 10–5 s	670	617	568	375	342	193	287	258	148	106
t = 5 ∙ 10–5 c / t = 5 ∙10–5 s	585	538	495	334	305	177	260	234	139	103

 

 

 

Рис. 3. Схема расположения точек, в которых вычислялась температура

Fig. 3. The location of the points at which the temperature was calculated

 

 

На рис. 4 представлена динамика температурного поля и вычисляемого потока в соответствии с формулой (5). Кривые на рис. 4, а иллюстрируют убывание температурного поля начала координат (линия U) и потока (линия ψ). На рис. 4, b показана динамика температурного поля сечения y = 0, z = 0; на рис. 4, c ‒ динамика точек, расположенных вне основания капли [Там же].

 

 

 

 

Рис. 4. Динамика температурного поля и теплового потока: а) изменение температуры капли
(линия U) и теплового потока (линия ψ) во времени; b) зависимость температуры поверхности
октанта от координаты в сечении y = 0, z = 0; с) изменение температуры внутренних точек
октанта, расположенных вне основания капли, во времени

 

Fig. 4. The dynamics of the temperature field and heat flux: a) change in temperature of the droplet
(line U) and heat flux (line ψ); b) the dependence of the octant surface temperature on the coordinate
in the section y = 0, z = 0; c) the temporal variation of the temperature of the internal octant points
located outside the base of the drop

 

 

Обсуждение и заключение

Анализ полученных результатов показывает, что поле в точках, расположенных в основании капли, постоянно убывает (рис. 4, а); в точках, расположенных вне основания капли (как на поверхности, так и в глубине катода), сначала возрастает, а затем убывает вместе с полем капли (рис. 4, с). Это свидетельствует о том, что, пока капля достаточно горячая, тепла к данным точкам прибывает больше, чем убывает, что и приводит к увеличению температуры. Однако по мере охлаждения капли приток тепла становится меньшим, чем отток, что приводит к уменьшению температуры [Там же].

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

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

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

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

 

 

¹              ohnson R. N. Principles and applications of electro-spark deposition // Society of Vacuum Coaters 45^th Annual Technical Conference, Lake Buena Vista. 1987. URL: https://www.researchgate.net/publication/236399429_Principles_and_applications_of_electro-spark_deposition

²              Ладыженская О. А., Солонников В. А., Уральцева Н. Н. Линейные и квазилинейные уравнения параболического типа. М. : Наука, 1967. 736 c.; Friedman A. Partial differential equations of parabolic type. R. E. Krieger Pub. Co., 1983. 347 p.

³              Friedman A. Partial differential equations of parabolic type. R. E. Krieger Pub. Co., 1983. 347 p.; Санчес-Паленсия Э. Неоднородные среды и теория колебаний / Пер. с англ. М. : Мир, 1984. 472 с. URL: http://bookre.org/reader?file=441528

⁴              Программный комплекс для моделирования теплопереноса материала при электрофизическом воздействии : свидетельство о гос. регистрации прогр. для ЭВМ / Власенко В. Д., Колисова М. В. № 2017611479 ; заявл. 06.12.16 ; опубл. 03.02.17. URL: http://www1.fips.ru/fips_servl/fips_servlet?DB=EVM&rn=4256&DocNumber=2017611479&TypeFile=html

⁵              Крылов В. И., Шульгина Л. Т. Справочная книга по численному интегрированию. М. : Наука, 1966. 372 с

 

Об авторах

Виктор Дмитриевич Власенко

ФГБУН «Вычислительный центр Дальневосточного отделения Российской академии наук»

Email: vlasenko@as.khb.ru
ORCID iD: 0000-0001-7782-4532
ResearcherId: E-2432-2019

Ученый секретарь, кандидат физико-математических наук

Россия, 680000, г. Хабаровск, ул. Ким Ю Чена, д. 65

Валерий Игоревич Иванов

ФГБНУ «Федеральный научный агроинженерный центр ВИМ»

Email: tehnoinvest-vip@mail.ru
ORCID iD: 0000-0002-4568-8553
ResearcherId: H-4076-2018

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

Россия, 109428, г. Москва, 1-й Институтский проезд, д. 5

Вячеслав Федорович Аулов

ФГБНУ «Федеральный научный агроинженерный центр ВИМ»

Email: gosniti@mail.ru
ORCID iD: 0000-0001-6925-1260
ResearcherId: E-4179-2019

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

Россия, 109428, г. Москва, 1-й Институтский проезд, д. 5

Леонид Алексеевич Коневцов

ФГБУН «Институт материаловедения Хабаровского научного центра Дальневосточного отделения Российской академии наук»

Email: konevts@narod.ru
ORCID iD: 0000-0002-7212-3953
ResearcherId: H-4087- 2018

научный сотрудник, кандидат технических наук

Россия, 680042, г. Хабаровск, ул. Тихоокеанская, д. 153

Елена Геннадьевна Мартынова

ФГБОУ ВО «МГУ им. Н. П. Огарёва»

Автор, ответственный за переписку.
Email: el.mart2012@yandex.ru
ORCID iD: 0000-0002-6870-0498
ResearcherId: С-5023-2019

аспирант, кафедра технического сервиса машин

Россия, 430005, г. Саранск, ул. Большевистская, д. 68/1

Исмаил Халил Хасан

ФГБОУ ВО «МГУ им. Н. П. Огарёва»

Email: srorismael@gmail.com
ORCID iD: 0000-0002-4560-1016
ResearcherId: С-5025-2019

аспирант, кафедра физики твердого тела

Россия, 430005, г. Саранск, ул. Большевистская, д. 68/1

Список литературы

Дополнительные файлы