Two-way interaction between switching arc and solid surfaces: distribution of ablated contact and nozzle materials

This paper is devoted to two-way plasma-surface interactions by investigating how the plasma arc ablates the nozzle and contacts and how the distribution of ablated materials changes the plasma parameters. For this purpose, a two-dimensional time-dependent model, in axial symmetric coordinates, for an arc at atmospheric pressure burning within a polytetrafluoroethylene nozzle is created. A computational fluid dynamics equations system is solved for plasma velocity, pressure, temperature, magnetic vector potential, and electrical potential. Radiation is modeled based on net emission coefficient and contacts, and nozzle ablation is also considered to better describe the arc formation, contact cooling, and arc temperatures, more precisely. The sublimated materials from contacts and nozzle will be used to calculate the distribution of plasma composition (i.e. ablated mixture ratio). The calculated ratio is used to change the plasma parameters, and data processing techniques are utilized to derive particle distribution and temperature profiles of the arc to investigate its thermo-electrical behavior. The simulation results show good agreement with the measurements obtained in an experimental setup already designed and published. This study provides support to the experimental work and contrariwise. The presence of ablated points on nozzle and contacts, which further modify plasma parameters and therefore the arc voltage are discussed.


Introduction
Although many investigations have been performed on MHD arc simulations for low voltage miniature circuit breakers (MCB) as well as for medium and high voltage SF 6 circuit breakers, there are almost no investigations on MHD arc simulation considering the simultaneous interaction of material on plasma and vice versa.
Most thermal plasma applications involve interactions between the plasma and solid or liquid surfaces. Computational models of thermal plasma applications typically consider the influence of the plasma on the surface, but often neglect the effect of the surface on the plasma. In many cases, however, it is not possible to model the plasma process accurately without taking the two-way interactions into account [1].
There are few studies about these interactions, but the research trend in this field is increasing recently [2][3][4][5][6][7], but most are focused on SF 6 [4,7,8] and CO 2 , and some of them have investigated both contact and nozzle interactions with plasma [9][10][11]; fewer of them have studied the distribution of particles inside the plasma [12].
There is still more work required to obtain precise plasma parameters, including the admixtures from contacts and nozzle material. Moreover, further research is necessary to evaluate the distribution of these admixtures in the plasma and to understand its impact on radiation and conductance of the plasma as well as on the interruption capability. For the dielectric recovery, precise breakdown models and validations are needed. The work for identification of alternatives to SF 6 is ongoing [13], where one of the environmentally-friendly gases in nature is air. So, we focus in this paper on investigating the numerical studies of switching arc characteristics in air.
To overcome the difficulties of simulating seven coupled physics simultaneously with moving mesh, the plasma is modeled in the local thermodynamic equilibrium in a layered and very thin state to simulate the arc in two-dimensions (2D); it is a trade-off between the time, accuracy and complexity. 2D simulation results are accurate enough when there are no strong vortices in flow [14]. A three-dimensional (3D) plasma can be assumed as 2D, if it is in local thermodynamic equilibrium (LTE), and it is in a very thin state so that the effects of assumption of symmetry on an asymmetric turbulent state can be ignored and it can be modeled as a laminar state.
There is a border between the thermal plasma column inside the arc and the turbulent gas flow around it [15]; so no large vortices field has been found inside the arc column. If the mesh is fine enough to resolve the size of the smallest vortices in the flow, the laminar model and time-dependent Naiver-Stokes (TDNS) will be suitable to reduce the simulation time.
The approach and outcomes of the present paper are shown in figure 1. The numerical model is based on the magnetohydrodynamics and the radiative transfer equation through a net emission coefficient (NEC), together with a detailed treatment of heat transfer in contacts and polytetrafluoroethylene (PTFE) nozzle considering the contact and nozzle ablation. The material sublimation is simulated using a FEM model and imported as an initial condition to the next study utilizing fluid and particle interactions to calculate the distribution of molecules and ions inside the arc chamber. The results of this study are used to modify the radiation and the electric conductivity of the nonhomogenous mixture with the time-variable mixing ratio of ablated materials to calculate the voltage and temperature profile of the electric arc in the next run of the MHD model.
Measured arc current is as input to the FEM model, and the simulated voltage will be one of the outputs that will be compared with the measured one.

Geometry, initial, and boundary conditions
The simulated geometry, including initial and boundary conditions, is shown in figure 2(a). A precise triangle mesh, including 37 000 elements with dense configuration near the contacts and nozzle, is applied as shown in figure 2 To get accurate results, it is important to have the minimum and average mesh element quality in general bigger than 0.1, and the ideal quality is one [16]. For the physics and solver used in this study, the minimum and average quality of 0.323 and 0.856 is used.

Governing equations
Here, a short description of the physics and relevant equations are presented; more details can be found elsewhere [17]. Frequently used notations are presented in table 1. The rest will be defined in the text.
2.2.1. Current conservation equations. The first group of the modeled equations in plasma is related to the current conservation equations: The measured current ( I measured ) divided by the cathode area is applied to the cathode as a current source J n (t) and the anode is grounded.

Magnetic field equations.
Magnetic field equations in plasma include the Ampere laws and Lorentz force, as shown in the following equations: − → J e in equations (3) and (7) is defined in equation (9) as the static current density resulted from the electric field. A realistic simulation of the arc splitting process needs special treatment of the arc roots [18]. One of these requirements is the modeling of electron emission from electrodes. A comparison between different electron-emission equations for different pressures and surface temperatures [19] shows that equations (10)- (12) for the cathode are still valid for the pressure and surface temperature range (1 bar and 1400 K) in this work. Although the influence of thermionic electron emission is not significant on the current density as at 1356 K will be 4.16 × 10 −5 A m −2 while the current density at the cathode is in the order of 9.3 × 10 8 A m −2 . Increasing T does not ensure observable electron emission since copper will melt and then evaporate or decompose when heated but it affects the contact temperature which is essential for the ablation model.
The required parameters of copper contacts can be found elsewhere [20,21]. These are φ eff = 4.94 V, A R = 120.2 A (K · cm) −2 and V ion = 15.7 V.

Heat transfer equations.
Heat transfer equations in plasma are shown in equations (13) and (14).
According to equation (15) for the cathode, a change in the cathode temperature is occurred by absorption of energy and electron emission from the cathode, i.e. overcoming the space charge energy level, and overcoming the potential of plasma ionization layer of positive ions near the cathode region [22].
According to equation (16) at the anode surface, the change in the anode temperature is a result of absorption of electrons at the anode, i.e. overcoming the energy level due to the accumulation of electrons near the anode.
The outer boundaries of the simulated region are always considered at ambient temperature and pressure.

Flow equations in plasma.
All gases and most nonviscous liquids are Newtonian materials [23]. In Newtonian fluids, if the Mach number is more than 0.3 or the temperature inside the fluid changes significantly and rapidly, the fluid is considered compressible. In our study, although the temperature rises sharply, the changes are not in a vast range, and calculation shows that the Mach number is not more than 0.3. Also, the Reynolds number is relatively small, so the fluid flow in each layer is laminar. The plasma core temperature is in the range of 10-18 kK, and by reducing the nozzle diameter reaches 26 kK. Densities of ions and electrons at this temperature are high enough to consider the air plasma as fully ionized. Therefore, Navier-Stokes equations for an incompressible fluid, including mass and momentum conservation in addition to Lorentz force are applicable.
It is important to mention that the eliminated time derivative of the density in equation (17), is only valid if the flow is incompressible. In the case studied in this paper, this assumption is, however, justified, because of low velocity of the flow within the nozzle compared to the thermal velocity of the ablated particles; the assumption of the incompressibility reduces the calculation costs significantly while keeping good accuracy.
2.2.5. Ablation model. Ablation helps cooling and affects arc voltage characteristics. The PTFE ablation occurs in different temperatures with different products. The CF 4 and C 2 F 4 , dominant decomposition products, equalize at a temperature higher than about 2000 K, and the full equilibrium state for all species is reached around 3500 K, being a typical value for the vapor temperature of PTFE ablation [24,25]. The details of polymer decomposition are too complicated for describing with a single distribution function as an ablation model. Dominant features of ablating nozzle system are shown to be based on 11 molecular reactions [8] and ten ionizations defined for PTFE [9], but a widely used simplified ablation equation for the nozzle can also be applied [26]: If ablation is not too strong and the Mach number (χ) is much smaller than one, then χ 2 can be neglected. The   vaporization enthalpy h v that is necessary to produce gaseous C 2 F 4 at 1000 K is set to 2.6 MJ kg −1 [25]. In addition, the specific heat dh PTFE to heat up the PTFE to the effective ablation temperature of 3500 K is taken into account [26]. The cooling effect of ablation is modeled through the boundary heat source Q surface as described in equation (21); it is coupled with heat equations through T.
T ablate is the ablation temperature which is 1356 K for Cu. δ is the Dirac delta function. Moving mesh is coupled to calculate the sublimated material in this physics.
2.2.6. Particle tracing model. Almost none of the other researches mentioned above has simulated the distribution of admixtures created from PTFE and metal vapor in detail. In one of the previous studies [26], the authors address this topic, but no detail is shown. There are only some practical measurements of one-way [27] or two-way measurements [1] of plasma effects on nozzle materials through high-speed imaging and considering just drag and gravity forces or experimental injection of particles to the plasma. In this work, a detailed distribution of species movement based on equations (22)-(30) is presented. It is assumed that the impact of the particles on the flow field is negligible. This assumption is valid when the mass of particles is less than 1% of fluid and particle-particle reactions, including Coulomb, linear elastic, Lennard-Jones, are ignored, as particle distribution is sparse [28]. Motion process and heating process of particles and simple charging process of Cu atoms, by arc plasma, is considered in the model. transient behaviour of ablation considering mass reduction of macroparticles during first 20 µs of ablation process is discussed in a recently published literature [28] but it is mentioned in another publication [27] that PTFE does not generate macroparticles and we are focusing on plasma effects on arc voltage profile; so, decomposition of macroparticle of copper contact in first 20 µs of ablation is not important. Then we assume the PTFE and Cu evaporate in molecules and atoms. First, the flow field is computed, and then, the motion of the particles is calculated as an analysis step. The total force, − → F t , in the second law of Newton (equation (22)) is the sum of all forces acting on the particles consisting of drag force, The latter is dependent to temperature. m p is the particle mass.
The drag force represents the force that a fluid exerts on a particle due to a difference in velocity between the fluid and the particle. It includes the viscous drag, the added mass, and the Basset history term [29]. u is the fluid velocity from flow equations and v is the particle speed.
τ p is the particle velocity response time for the stokes model whereas applicable to the low Reynolds numbers ( R e 1) and spherical particles. µ is the fluid viscosity.
Critical Reynolds number for each configuration must be determined through experiments or numerical simulations. The studied case is, however, similar to flow inside a cylinder, whereas Re ⩽ 0.03 was calculated and a steady, predictable laminar flow is observed [30].
M is wall correction factor [31], L is the distance to the nearest wall, and P(n) is the projection onto the vector normal to the wall. When the particle contacts the walls, it takes the wall temperature and is scattered diffusely.
The other forces are shown in equation (27). ρ p and r p are the particle density and radius, respectively. The particle mass and size are calculated based on the molar weight and density.
− → g is the gravity vector k B is Boltzmann constant and ζ is a Gaussian random number with zero mean and unit variance. The random direction of the Brownian force was accounted for by evaluating both the r and z components of − → F b at each time step using independent values of ζ in the two directions [32] and Δt is the time step defined in equation (31). It is assumed that particles will separate from the wall and are distributed with a random thermal velocity having a Maxwell-Boltzmann distribution [28].
Knudsen's cosine law shown in equation (29) applies to the particle speeds. Γ is a random number.
Particles temperature changes due to equation (30). If the particle temperature reaches the ionization temperature, it will be considered as a positive ion in the field.
h is the heat transfer coefficient, T p is the particle temperature and A p the particle surface area. The Biot number (Bi), a simple index of the ratio of the heat transfer resistances inside of a body and at the surface of a body, is much smaller than 1. So, it is assumed that there is no resistance to heat transfer inside the particles, i.e. particle temperature is uniform [28]. The smaller the particle, the smaller time steps, and an increase in the number of particles is equivalent to a higher degree of freedom (DOF) of the equations. Equation (31) shows the relationship between maximum time step, Δt, and the particle time constant.
The time steps of at least one picosecond at the start of the solution and a maximum of 0.2 µs were selected; this is vital for precise modeling of the arc voltage transients. Jacobian matrix was updated at each time step and considered absolute tolerance of simulation was 5 × 10 −4 .
To have an estimate of the particle density distribution in time, ∂N d ∂t , an accumulator is defined in equation (32) based on the total number of species in each mesh, R i , divided by mesh volume, V mesh , that is equal to the mesh area in 2D axial symmetry.
The obtained variable is utilized to modify the temperature-dependent properties of the plasma mixture.

Radiation model.
An exact formulation of the energy transport by radiation is a very complicated procedure, but three methods are mostly used named P1, ordinary discrete method (ODM) and NEC where the latter shows the difference between the radiated and absorbed power. A compariso n between these methods is discussed elsewhere [33]. In the range of temperatures and thicknesses of the modeled arc, NEC has acceptable accuracy and much less complexity [34].
In this paper, the model will be applied to an arc burning inside nozzle with a radius of 1 and 2 mm (this is the same geometry used in experimental studies reported earlier [35]). This implies that the maximum radius of the core varies between 1-2 mm. So, the radiation absorption radius, R p , in calculating the NEC is assumed 2 mm. The arc in the experiments were initiated by an ignition wire between two identical copper spherical head electrodes of 2 mm diameter, which is a usual arc initiating method in practical experiments [36][37][38]. The insulating cylindrical tubes with a length of 20 mm and different inner diameters are used. We have modeled air-Cu mixture and absorption effect for NEC. If absorption is ignored (in thin layered arcs), then most of the metal vapors increase the NEC. An increase in arc radius decreases NEC [39]. The ablated nozzle does not participate in radiation but exploding wire of size 25 µm generates 1% Cu vapor-99% air mixture [40] in a pressure tank of 15.7 l at the start of the simulation. The air data is modified based on results reported elsewhere [41]. Thermal radiation of PTFE is neglected since T p is not so high [8], and its effects on the arc through affecting transport properties of the mixture are considered as described in literature [10]. Surface to surface radiation between electrodes and nozzle and from the inner walls of the nozzle were modeled and is ignorable compared with the arc radiation.

Material properties
2.3.1. Hot air plasma. Thermodynamic properties of the pure hot air in constant pressure were presented by Yos in 1963 [42] then revised in 1967 [43] and 1990 [44]. Some other  calculations were published by Murphy [45,46], Capitelli [47][48][49] and Boulos [50]. Comparisons between the mentioned reports were presented in [48,51,52]. Effects of pressure on thermodynamic properties were introduced in [53,54]. Simulation results show that the pressure change is less than 0.15% in our model, so the impact of pressure on thermodynamic properties can be ignored here. Thermal plasma properties in mixtures of air and metallic vapors have been discussed in an earlier study [41]. Based on the results presented on thermal plasma properties for mixtures of air and copper [53], the metal vapor has a considerable influence on the electric conductivity and on the radiation but a negligible effect on the thermal conductivity and heat capacity [55]. The presence of metal vapors increases the electric conductivity of the plasma for temperatures lower than 10 000 K. So, its effect on electrical conductivity is preponderant in low-current arcs whereas at high currents, the impact on radiation outweighs [56]. To obtain the thermodynamic properties of hot air, the practical results in previous studies [48,53] have been used. Material properties of pure dried air including density, heat capacity, dynamic viscosity, thermal conductivity, electrical conductivity, and NEC are presented in figure 3 for the temperature range up to 25 000 K at the constant atmospheric pressure. The data is modified through results from particle tracing to estimate the properties of the plasma mixture and to achieve the maximum accuracy for arc voltage simulation.

Solid parts.
Frequently used parameters for modeling the solid parts, including copper contacts and PTFE are presented in table 2. The rarely used parameters will be defined in the text.

Arc Ignition
To initiate an arc in the air through a narrow tube and getting convergence on equations, the minimum conductivity of air is increased to 1 S m −1 , then a voltage ramp of 500 V µs −1 to 10 µs, 100 V µs −1 to 20 µs, 40 V µs −1 to 22 µs, and 1 V µs −1 to 50 µs then fixed voltage to 66 µs is applied.
As shown in figure 4, it results in an area with temperature higher than 2830 K and up to 7500 K around the contacts (green borders). Next, 300 ADC current is applied up to 100 µs. This opens a conductive channel with conductivity higher than 2000 S m −1 and the temperature higher than 8700 K. Finally, 150 A alternating current is applied from 100 µs.

Results
Sublimated surface and mass fraction curves, velocity contours in plasma, core temperature profiles and particles distribution and densities are presented. Simulated arc voltage based on these results is compared with the measured.

Sublimation
The sublimated surface of the cathode at 1356 K at six different times, temperature contours for T = 300, 600, 900, and 1200 K and heat flux surfaces are shown in figure 5(a). The rate of mass ablation at edges of the nozzle tube is around 30 times higher than the inner surface which is in line with previous measurements [57].
The nearest corner to the cathode is ablated 1.5 times more compared with the other corner close to the anode in first half cycle from the cold state but in the next half cycle the other side sublimates 10 times faster due to its preheated state. The sublimated surface of the nozzle edges are shown in figures 6(a) and (b). The cooling effect of ablation is clear from figure 6. The temperature contours of T = 300, 700, 1100 and 1500 K are shown at nozzle surface in cyan. Just 150 µm of the inner thickness of the nozzle has a temperature higher than 300 K. The sublimated surface inner side of the nozzle is very narrow, and it is shown in figure 6(c). On each cycle, the electrodes interchange their roles as anode or cathode, when the current reverses [58]. Due to the direction of current injection in a simulated half cycle, the top contact acts as anode and the lower contact is the cathode. Figure 7 shows the effects of ablation and nozzle size on the maximum temperature of the cathode. Without the ablation cooling, the cathode temperature can reach 2300 K while it is melted at 1356 K. The cathode starts ablating at around 700 µs (zero crossing) and almost at the same time in both cases with different nozzles. The latent heat of fusion for Cu is 205 (kJ kg −1 ), and the heat of vaporization is 300.3 kJ mol −1 ; so, the heat of sublimation for Cu at 3081 K is 313.3 kJ mol −1 . Figure 8 shows the effects of nozzle ablation and size on the maximum temperature of the nozzle surface. Without the ablation cooling, the nozzle surface for D Nz = 2 mm can reach temperatures as high as 4000 K while it is ablated at 2800-3300 K. The nozzle surface starts ablating at around 1300 µs for D Nz = 2 mm, and in case of the 4 mm nozzle there almost is not ablation as the nozzle temperature does not reach more  than 1200 K. The experiments reported earlier confirm this results [40].
The simulation shows that ablation only occurs on the cathode in each half cycle and as it is shown in figure 9, the cathode sublimation is 135 µg, which is equal to 316 µg C −1 for a 350 Hz and 150 A peak current in half cycle. Ablation mass for copper electrodes in low currents is reported around 100-300 µg C −1 [59], which is in line with our findings.
The ablated mass of 2 mm PTFE nozzle in the first half cycle from the cold state is 680 µg while in the 200 µs of the next half cycle due to its hot state; it will be 800 µg. Figure 10 shows the velocity surface in the plasma inside the electrical conductivity contour with sigma >500 S m −1 . It is revealed that the flow velocity in the arc core (σ > 500 S m −1 ) inside the nozzle is in range of 3.6-4.4 m s −1 at the arc center and it decreases to less than 1 m s −1 at 250 µm from the walls. So, the flow will not move PTFE particles more than 8.8 mm in 2 ms from their positions. Out of the nozzle and near to its edges, velocity is in range of 8-40 m s −1 except near the cathode where it reaches locally 132 m s −1 due to its high temperature. The flow velocity in the entrance of the nozzle at both ends is higher than 40 m s −1 and its maximum happens at the nozzle's throat near to the anode, which is 149 m s −1 .

Plasma velocity
At the arc core center and inside the inner gray contour in figure 11, although the flow speed is in range of 40-132 m s −1 , it is clear from the cyan streamlines that the flow direction reverses in short distances from the cathode. So, the effect of the flow speed in total is ignorable.  Our simulation with the diffusion and gravity as only forces applied shows a maximum displacement of ablated contact material in the order of 100 µm which was reported elsewhere [4]. The effect of the flow speed is ignorable, and therefore, the drag force is not dominant, but the thermal velocity is the primary reason for spreading the particles in the plasma region. So, the assumption of equation (17) is valid Figure 12 shows the particle distribution in the plasma region in three different times. Particles from the cathode, nozzle edges and the inner surface of the tube are displayed with different thermal speeds they reach. Particles have a Maxwellian speed distribution in the range of 0-500 m s −1 and mostly around 150-250 m s −1 . So, the small gas flow inside the nozzle and around the edges has no important effect on the particle distribution.

Particle trajectories and densities
The calculated mixture densities of the copper and PTFE are shown in figure 13. Other studies showed when the current is lower than one kA; metal vapor resides only in two small regions (1-3 mm) in front of the two electrodes because of limited electrode vaporization [56,60]; our simulations confirm this reported fact. The ablated metals remain mostly near to the surfaces; then, arc voltage will not be so sensitive to the presence of metal vapour [60] from ablation. Ablated cathode particles that are ionized due to  temperature increase, are exposed to the electric field and cannot move away.

Temperature profile
Profiles of temperature and labelled electrical potential contours along the core are shown in figures 14(a) and (b). The highest temperature and the maximum conductivity have a local peak near the anode and cathode, but at the current zero, these points will be cooler than the middle of the arc column. This temperature drop increases the resistance and voltage drop at both ends. The voltage profile inside the nozzle is linear. The cooling resulted from ablation affects the temper ature profile at the nozzle surface and affects the nozzle temperature. There are local temperature extremums between contacts and nozzle.
The study shows that the maximum temperature happens further from the nozzle with reducing the nozzle diameter, but the minimum occurs closer to the nozzle; it means the temper ature gradient increases by decreasing the diameter. The reason is an increase in flow speed and Lorentz force and moving the hot plasma toward the inside of the tube. It means that the pinch effect of the narrow nozzle has increased the arc column temperature to points higher than the cathode region.    This observed change in the simulated temperature with considering the ablation effects is in line with other studies [61]. Figure 17 provides a comparison between the measured and simulated arc voltages. It can be seen that both simulated and  measured voltages jump from zero to 350-380 V, then step down to about 100 V and 150 V for 4 mm and 2 mm nozzles, respectively. The arc voltage shape complies with the measurements except at current zeros where the temperature falls to below 9 kK, and therefore the plasma cannot be considered as fully ionized at these temperatures; so, the arc behavior suddenly changes.

Arc voltage versus time
The maximum current density was simulated near the cathode equal to 9.3 × 10 8 A m −2 . Current density at the cathode is independent of current magnitude and mostly depends on the cathode material [62].
The arc current is similar for both cases with 2 mm and 4 mm nozzles, but the arc voltage is almost 50% higher with D Nz = 2 mm and temperature near the cathode is lower due to ablation explained in 3.4. So, the maximum temperature inside the core increases from 18 kK for D Nz = 4 mm to 26 kK for D Nz = 2 mm to keep the arc energy balance. Figure 18 shows the core temperature at the midpoint of the nozzle for 2 mm and 4 mm nozzle diameters. Numbers in the figure legend are different radial positions.

Arc temperature and core radius
It is evident that the arc temperature at the peak of the current is 26 kK and 16 kK, respectively. The core radius is 0.97 mm and 1.85 mm by considering 6000 K as arc wall temper ature. The calculated arc temperature and core radius for 4 mm nozzle diameter based on simplified models [63,64]   are 14 kK and 1.4 mm; this arc radius corresponds to the arc wall temperature of 11 000 K.

Conclusion
A detailed 2D FEM model of air arc in PTFE nozzle was presented. The presented numerical model was built basis on magneto-hydrodynamics and the radiation-absorption through NEC in addition to ablation and particle tracing and was validated theoretically and experimentally. Sublimation of contact and nozzle materials was revealed, and the effects of ablated particles on arc parameters were investigated for a low current air arc burning inside a PTFE nozzle. All assumptions and simplifications were explained in detail and were validated through comparison with other researches and physical modeling. Thermodynamics and electrical behaviors of arcs as well as sublimated particle's distributions were extracted through simulation and were validated through physical measurements. It is observed that simulated values of temperature, voltage, and arc behavior are all consistent and confirm measured and theoretical predictions for similar setup [40].