Account of near-cathode sheath in numerical models of high-pressure arc discharges

Impressive advances achieved in the modelling of high-pressure arc discharges have been summarized in recent review papers [1, 2]. Significant improvements have been acknowledged, attained over the last 15 years in the modelling of cathode and anode phenomena in order to couple the electrodes with the bulk plasma in a general model. On the other hand, a clear need for further improvement has been pointed out, in the first place, as far as the cathode phenomena are concerned. It is worth noting that the near-cathode space-charge sheath plays an important role in the energy balance of arc discharges: substantial electric power is deposited by the arc power supply into the sheath. Therefore, an adequate account of the near-cathode sheath in the modelling of high-pressure arc discharges is indispensable for the modelling to be physically justified.

A number of different models to account for near-cathode sheaths in numerical simulations of high-pressure arc discharges have been described in the literature. The aim of this work is to systematize these models and to compare the results wherever possible.

The outline of the paper is as follows. Different approaches to accounting for near-cathode sheaths in numerical models of high-pressure arc discharges are discussed in section 2. Characteristics of the cathodic part of arc discharges predicted by the different approaches are compared in sections 3 and 4. Conclusions are summarized in section 5. Boundary conditions describing positive and negative space-charge sheaths in non-equilibrium models of arc discharges are derived in the appendix.

The physics of current transfer to electrodes of high-pressure arc discharges is rather complex and involves processes characterized by substantially different length scales, such as the Debye length, the mean free paths of the ions, electrons and neutral atoms, the lengths of Maxwellization of emitted electrons and the relaxation of their energy, the ionization length etc. There are a large number of publications dedicated to different aspects of near-electrode physics in arc discharges, e.g. the reviews [3–12]. Without trying to provide another review, we will rather focus on a particular aspect: methods of quantitative description of charge separation in the near-cathode region.

2.1. Unified modelling

2.2. Separate treatment of the plasma bulk and space-charge sheaths

The rest of this work is concerned with comparing characteristics of the cathodic part of arc discharges computed by means of three approaches: the unified modelling, the NLTE-sheath approach and the model of nonlinear surface heating. Since the unified modelling has been performed until now only in 1D cases [13–15], while the NLTE-sheath approach has been realized only for the 2D (axially symmetric) case [28], we cannot compare all three approaches at once. In this section, we compare the unified modelling and the model of nonlinear surface heating. The NLTE-sheath approach and the model of nonlinear surface heating are compared in the next section.

Let us consider an arc cathode having the shape of a right cylinder with the base being externally cooled; the electrical conductivity of the cathode material is high enough so that Joule heat production in the cathode body is negligible; the lateral surface is thermally and electrically insulated, so the energy flux and the electric current from the plasma enter the cathode through the front surface and leave it through the base. The heat conduction equation in the cathode body admits a 1D solution in this geometry, with the temperature distribution varying only in the axial direction and all parameters being constant along the cathode surface. This makes such geometry a convenient test case for the unified modelling.

More precisely, we compare characteristics of the cathodic part of the discharge computed by means of two approaches. In the framework of the first approach, the near-cathode region, separating the cathode from the bulk plasma where LTE holds, is simulated by means of the 1D unified modelling code [13]. The cathode bulk is described by the 1D solution of the heat conduction equation, with the energy flux from the plasma to the cathode surface being matched to the heat flux from the surface into the cathode bulk. The other approach is the model of nonlinear surface heating, based on the model of the near-cathode layer [29–31] (a summary of equations of this model can be found in [43]) and realized by means of a tool available on the internet [44]. While the tool is capable of computing both the diffuse and spot modes of current transfer, the results given here refer to the diffuse mode, where the distribution of temperature in the cathode is 1D.

All the results of this work refer to the cathode material being tungsten and the plasma-producing gas being argon at 1 bar, except for figure 6 which refers to mercury vapor at 200 bar. Simulations reported in this section refer to a cathode height equal to 1 cm with the bottom of the cathode maintained at room temperature.

The choice of calculation results to compare is not quite trivial since the unified modelling [13] and the model of the near-cathode layer [29–31] involve different sets of variables. The most natural choice seems to be as follows. The results will be compared for the same value of the current density, which is equivalent to the same value of the arc current in the 1D case being considered. The parameters to be compared are the cathode surface temperature ${{T}^{\left(\text{s}\right)}}$; the near-cathode voltage U; the sheath voltage U_sh; the average electron temperature in the ionization layer, T_eil; the densities of the ion current to the cathode surface, j_iw, the electron emission current, j_em, and the current of fast plasma electrons counter-diffusing to the cathode surface against the sheath electric field, j_CD; the density of the net energy flux to the cathode surface, $q_{\text{w}}^{\left(\text{s}\right)}$; densities of energy fluxes transported to the cathode surface by the heavy particles and the fast plasma electrons, q_hw, q_CD, and from the cathode surface by the emitted electrons, q_em. All these quantities appear in a natural way in the model of the near-cathode layer [29–31], but most of them have not been expressly specified in the unified modelling [13] and for completeness need to be related to the variables employed there. The quantities j_CD, q_em and q_CD are expressed as

$$\begin{eqnarray}&&{{j}_{\text{CD}}}=\frac{e{{n}_{\text{ew}}}}{4}{{\left(\frac{8k{{T}_{\text{ew}}}}{\pi {{m}_{\text{e}}}}\right)}^{1/2}},\;\;\;\;{{q}_{\text{em}}}=\frac{{{j}_{\text{em}}}}{e}\left({{A}_{\text{f}}}+2k{{T}^{\left(\text{s}\right)}}\right),\\ &&{{q}_{\text{CD}}}=\frac{{{j}_{\text{CD}}}}{e}\left({{A}_{\text{f}}}+2k{{T}_{\text{ew}}}\right),\end{eqnarray} \tag{ 1 }$$

where A_f is the work function of the cathode material and the subscript 'w' in the unified model means that the corresponding quantity (in this case, the electron number density n_e or temperature T_e) is evaluated at the cathode surface. The sheath was identified in the results of unified modelling as a near-cathode region beyond which the charge separation is below $10\%$. The ionization layer was identified as a region adjacent to the sheath beyond which the deviation of the charged particle density from the value n_b obtained from the Saha equation, evaluated in terms of local values of T_e and T_h, is below $10\%$. (More precisely, ${{n}_{\text{b}}}$ was obtained from the assumption of a local balance of ionization and recombination with account of radiation escape. It deviates from the Saha value if the electron number density is of the order of 10²¹m⁻³ or smaller; see figure 5(a) of [13].) U_sh and U have been evaluated as potential differences between the edge of the sheath or the ionization layer, respectively, and the cathode. T_eil was evaluated as the average of the values of T_e at the edge of the space-charge sheath (which represents the inner boundary of the ionization layer) and at the edge of the ionization layer.

Densities of energy fluxes transported by each particle species in the direction of the cathode in the framework of the unified model are defined as (see equation (31) of [13])

$$\begin{eqnarray}&&{{q}_{\text{a}}}=\frac{5}{2}k{{T}_{\text{h}}}\,{{J}_{\text{a}}}+{{\kappa}_{\text{a}}}\frac{\partial {{T}_{\text{h}}}}{\partial x},\quad {{q}_{\text{i}}}=\left(\frac{5}{2}k{{T}_{\text{h}}}+{{A}_{\text{i}}}-{{A}_{\text{f}}}\right)\,{{J}_{\text{i}}}+{{\kappa}_{\text{i}}}\frac{\partial {{T}_{\text{h}}}}{\partial x},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{q}_{\text{e}}}=\left(\frac{5}{2}k{{T}_{\text{e}}}+{{A}_{\text{f}}}\right)\,{{J}_{\text{e}}}+{{\kappa}_{\text{e}}}\frac{\partial {{T}_{\text{e}}}}{\partial x}-{{A}_{\text{itd}}},\end{eqnarray} \tag{ 3 }$$

where the x-axis is directed from the cathode into the plasma, ${{J}_{\alpha}}$ and ${{\kappa}_{\alpha}}$ are the density of flux in the direction of the cathode and the thermal conductivity of species α ($\alpha =\text{a},\text{i},\text{e}$), A_i is the ionization energy and A_itd designates the term accounting for the effect inverse to the thermal diffusion (the second term on the rhs of equation (15) of [13]). The density of energy flux transported in the direction of the cathode by the heavy particles is ${{q}_{\text{h}}}={{q}_{\text{a}}}+{{q}_{\text{i}}}$.

A comparison of the above-listed quantities as functions of the current density is shown in figure 1. The agreement between the two approaches is good, except for a discrepancy in the current density and energy flux transported by the fast plasma electrons ( j_CD, q_CD) occurring in the range ${{j}_{\text{w}}}\lesssim 4\times {{10}^{6}}$ A m⁻². This discrepancy, however, has little relevance since contributions of the plasma electrons in this range are small anyway: ${{j}_{\text{CD}}}\ll {{j}_{\text{em}}}$, ${{q}_{\text{CD}}}\ll {{q}_{\text{hw}}}$.

Zoom In Zoom Out Reset image size

The net electron energy flux to the cathode surface is ${{q}_{\text{ew}}}={{q}_{\text{CD}}}-{{q}_{\text{em}}}$. (In the framework of the unified model, this expression conforms to equation (3) evaluated at x  =  0, which can be verified with the use of the boundary condition (25) of [13].) It is seen from figure 1(c) that ${{q}_{\text{em}}}>{{q}_{\text{CD}}}$ for both models, i.e. cooling of the cathode by thermionic emission prevails over the heating by plasma electrons, so q_ew is negative and the net electron energy flux is directed from the cathode surface into the plasma.

The net energy flux to the cathode surface in both models is evaluated as $q_{\text{w}}^{\left(\text{s}\right)}={{q}_{\text{hw}}}+{{q}_{\text{ew}}}$. However, the energy flux transported by the heavy particles is evaluated differently: account is taken of only the ion contribution (${{q}_{\text{hw}}}={{q}_{\text{iw}}}$) in the model of the near-cathode layer [29–31] and of both ion and neutral atom contributions (${{q}_{\text{hw}}}={{q}_{\text{iw}}}+{{q}_{\text{aw}}}$) in the unified modelling. In fact, q_aw exceeds q_iw in the framework of the unified model as seen from figure 2. In spite of this difference, values of the net heavy-particle energy flux q_hw predicted by the two models agree quite well, as seen from figure 1(c).

Zoom In Zoom Out Reset image size

It is seen from figure 2 that T_eil, the average temperature of electrons in the ionization layer, is close to T_ew, the electron temperature at the cathode surface in the unified modelling, for ${{j}_{\text{w}}}\gtrsim 4\times {{10}^{6}}$ A m⁻². (Note that the line T_eil in figure 2 represents the same data as the similar line in figure 1(b).) This suggests that the electron temperature does not vary much across the near-cathode layer. It is unsurprising that the integral-balance evaluation of T_eil, employed in the model of the near-cathode layer [29–31], is accurate enough under these conditions. On the other hand, one can see from figure 1(b) that values of T_eil predicted by the two models agree quite well also for ${{j}_{\text{w}}}\lesssim 4\times {{10}^{6}}$ A m⁻², where T_eil is significantly lower than T_ew and the electron temperature distribution across the near-cathode layer is therefore strongly non-uniform.

In view of the above, the reasons for the agreement seen in figure 1 are not obvious. In order to figure out these reasons, it is useful to look at examples of distributions of parameters in the near-cathode layer obtained from the unified modelling and shown in figures 3 and 4. Here n_i is the ion number density and E is the x-projection of the electric field. The vertical dotted lines here and elsewhere represent the edges of the space-charge sheath and the ionization layer, identified as described above.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The quantity $e{{J}_{\text{i}}}/{{j}_{\text{w}}}$ in figures 3(a) and 4(a) has the meaning of fraction of electric current transported by the ions. In the bulk plasma, this fraction is low, as should have been expected. On the cathode surface, this fraction attains a value of about 0.86 in the case ${{j}_{\text{w}}}={{10}^{5}}$ A m⁻² and 0.22 for ${{j}_{\text{w}}}={{10}^{7}}$ A m⁻². The ion current is generated primarily in the ionization layer, as expected.

In the case ${{j}_{\text{w}}}={{10}^{7}}$ A m⁻² (figure 4(b)), the electron temperature in the near-cathode layer does not change much. In contrast, the variation of T_e is very large in the case ${{j}_{\text{w}}}={{10}^{5}}$ A m⁻² (figure 3(b)). However, the variation of T_e in the ionization layer remains relatively small. This explains why the integral-balance evaluation of the average electron temperature in the ionization layer, employed in the model [29–31], is accurate over the whole range of j_w shown in figure 1.

Neglecting variation of T_e in the case ${{j}_{\text{w}}}={{10}^{7}}$ A m⁻², one can assume that the density of plasma electrons in the sheath follows the Boltzmann distribution evaluated at the temperature T_eil. Then one can write ${{n}_{\text{ew}}}=n_{\text{i}}^{\left(\text{pl}\right)}\exp \left(-e{{U}_{\text{sh}}}/k{{T}_{\text{eil}}}\right)$, where $n_{\text{i}}^{\left(\text{pl}\right)}$ is the density of the ions (or the electrons) at the sheath edge, and the expression for j_CD in equation (1) will coincide with the corresponding expression in the model [29–31]; see equation (A.7) of the appendix or equation (5) of [43]. This explains the good agreement in current density and energy flux transported by the fast plasma electrons ( j_CD, q_CD) seen in figures 1(b) and (c) for ${{j}_{\text{w}}}\gtrsim 4\times {{10}^{6}}$ A m⁻². For lower ${{j}_{\text{w}}}$, T_e in the sheath significantly exceeds T_eil, and the electron density in the sheath significantly exceeds values given by the Boltzmann distribution evaluated at the temperature T_eil. This explains why, for ${{j}_{\text{w}}}\lesssim 4\times {{10}^{6}}$ A m⁻², j_CD and q_CD given by the unified modelling significantly exceed the values given by the model of the near-cathode layer [29–31], as seen in figures 1(b) and (c).

One can see in figures 3(b) and 4(b) that in the bulk plasma the energy fluxes of the heavy particles, q_a and q_i, are small compared to the electron energy flux, as should have been expected. The net energy transport is directed not from the bulk plasma into the near-cathode layer, but in the opposite direction, which is also expected [34]. q_a and q_i grow in the near-cathode layer.

Equation (11) of [13], which describes conservation of energy of heavy particles in the unified model, may be transformed to

$$\begin{eqnarray}&&-\frac{\text{d}{{q}_{\text{h}}}}{\text{d}x}=e{{J}_{\text{i}}}\frac{\text{d}\varphi}{\text{d}x}-\left({{A}_{\text{i}}}-{{A}_{\text{f}}}\right)\frac{\text{d}{{J}_{\text{i}}}}{\text{d}x}+\frac{3{{n}_{\text{e}}}{{k}^{2}}{{T}_{\text{e}}}}{{{m}_{\text{i}}}n}\left({{T}_{\text{e}}}-{{T}_{\text{h}}}\right)\left(\frac{{{n}_{\text{a}}}}{{{D}_{\text{ea}}}}+\frac{{{n}_{\text{i}}}}{{{D}_{\text{ei}}}}\right),\end{eqnarray} \tag{ 4 }$$

where φ is electrostatic potential, D_ea and D_ei are binary diffusion coefficients, and n is the total number density of the plasma. The terms on the rhs of this equation account for, respectively, work of the electric field over the ions, energy gained by the ion species due to ionization, and the energy received by the atoms and ions in elastic collisions with the electrons. Integrating equation (4) from x to the outer boundary L of the calculation domain, one finds the integral balance of energy of heavy particles in the layer $[x,L]$:

$$\begin{eqnarray}&&{{q}_{\text{h}}}-{{q}_{\text{h}L}}={{W}_{\text{i}J}}+{{W}_{\text{ion}}}+{{W}_{\text{ea}}}+{{W}_{\text{ei}}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{W}_{\text{i}\,J}}={\int}_{x}^{L}e{{J}_{\text{i}}}\frac{\text{d}\varphi}{\text{d}x}\text{d}x,\;\;\;\;{{W}_{\text{ion}}}=\left({{A}_{\text{i}}}-{{A}_{\text{f}}}\right)\left({{J}_{\text{i}}}-{{J}_{\text{i}L}}\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{W}_{\text{ea}}}={\int}_{x}^{L}\frac{3{{n}_{\text{e}}}{{k}^{2}}{{T}_{\text{e}}}}{{{m}_{\text{i}}}n}\left({{T}_{\text{e}}}-{{T}_{\text{h}}}\right)\frac{{{n}_{\text{a}}}}{{{D}_{\text{ea}}}}\text{d}x,\\ &&{{W}_{\text{ei}}}={\int}_{x}^{L}\frac{3{{n}_{\text{e}}}{{k}^{2}}{{T}_{\text{e}}}}{{{m}_{\text{i}}}n}\left({{T}_{\text{e}}}-{{T}_{\text{h}}}\right)\frac{{{n}_{\text{i}}}}{{{D}_{\text{ei}}}}\,\text{d}x.\end{eqnarray} \tag{ 7 }$$

An example of the distribution of the terms on the rhs of equation (5) in the near-cathode layer is shown in figure 5. The energy flux at the outer boundary of the calculation domain is relatively small, ${{q}_{\text{h}L}}=3\times {{10}^{3}}$ W m⁻², and cannot be seen on the graph. As x decreases, q_h increases, first due to energy transfer from electrons to the ions in ionizing and elastic collisions (W_ion, W_ei) and then due to work of the electric field over the ions (${{W}_{\text{i}J}}$). The energy W_ea transferred from electrons to the atoms is small. Hence, the increase in q_h in the near-cathode layer occurs due to the energy supplied to the ions. Part of this energy is subsequently transferred to the atoms through ion–atom collisions. The latter process affects the distribution of the energy between the ions and the atoms, but does not affect the total energy of the heavy particles. This explains the very good agreement between values of q_hw given by the two models seen in figure 1(c).

Zoom In Zoom Out Reset image size

Note that the use of a single heavy-particle temperature in the unified modelling implies that thermal equilibrium between the ions and the atoms is established instantaneously. In reality, the (translational) energy is transferred from the ions to the atoms over distances of the order of the mean free path for collisions between ions and atoms, ${{\lambda}_{\text{ia}}}$. Let us consider as an example a high-current arc in atmospheric-pressure $\text{Ar}$ and assume that ${{T}_{\text{e}}}=15\,000$ K and ${{T}_{\text{h}}}=3500$ K in the near-cathode plasma. An estimate for these conditions gives ${{\lambda}_{\text{ia}}}\approx 2.5$ μm. It follows from figure 5 that most of the energy that the ions receive from the electric field in the near-cathode layer cannot be transferred to the atoms. Thus, values of the ion energy flux density q_iw, predicted by the unified modelling and the model of nonlinear surface heating, should be considered as, respectively, lower and upper estimates.

The model of nonlinear surface heating, as it is used in this work, comprises the model of the near-cathode layer [29–31], based on the assumption that the ions cross the near-cathode space-charge sheath without collisions. If the pressure of plasma-producing gas is much higher than atmospheric, the latter assumption becomes unrealistic. However, the model of nonlinear surface heating comprising the model of the near-cathode layer [29–31] may still remain useful. As an example, a comparison similar to that of figure 1 is shown in figure 6 for an arc burning in mercury vapor at 200 bar, which is a representative plasma-producing gas in projection lamps. (The dependences $q_{\text{w}}^{\left(\text{s}\right)}$ computed by means of the two models coincide to the accuracy of the graph, so the corresponding curves in figure 6(c) overlap.) There is a significant discrepancy in the current density and energy flux transported by the fast plasma electrons (j_CD, q_CD), but all the other quantities are close to each other.

Zoom In Zoom Out Reset image size

Let us consider a free-burning arc in atmospheric-pressure argon. The cathode is a tungsten rod of radius 1 mm and height 12 mm, with a hemispherical tip. The anode is planar and made of copper. A numerical investigation of this arc by means of the 2T-ionization layer-sheath approach was performed in [34, 35] and by means of the NLTE-sheath approach in [28]. In particular, the computed arc voltage was in good agreement with the experimental data [45]. In this work, we do not consider characteristics of the arc as a whole and will rather focus on comparison of characteristics of the cathodic part of the arc given by the two approaches.

The NLTE-sheath approach is used as described in [28]. The model of nonlinear surface heating, to which the 2T-ionization layer-sheath approach is reduced, is used as described in [34, 35]; the only difference from the previous section is that the nonlinear boundary-value heat conduction problem in the cathode body is solved by means of the steady-state solver of COMSOL Multiphysics instead of a home-made FORTRAN code implemented in the internet tool [44].

Distributions along the cathode surface of the surface temperature and the density of electric current to the surface, computed by means of the NLTE-sheath approach for the arc current I  =  100 A and three different values of the interelectrode distance (7, 10 and 13 mm), are plotted in figure 7 as solid lines. (The distance l is measured from the centre of the front surface of the cathode along the generatrix of the cathode surface.) The lines corresponding to the different interelectrode distances coincide. This is a manifestation of plasma–cathode interaction being to the first approximation independent of processes in the arc bulk, which constitutes the basis of the model of nonlinear surface heating as discussed in section 2.2.1. The results given by the model of nonlinear surface heating (dotted lines in figure 7) are in excellent agreement with those given by the NLTE-sheath approach.

Zoom In Zoom Out Reset image size

A comparison of distributions along the cathode surface of characteristics of the cathodic part of the discharge, computed with the use of the two models for three values of the arc current I  =  60, 100 and 200 A, is shown in figure 8. The distributions computed by means of the NLTE-sheath model for different interelectrode distances (7, 10 and 13 mm) are very close to each other; we note for definiteness that the data shown in figure 8 refer to the interelectrode distance of 10 mm.

Zoom In Zoom Out Reset image size

In most cases, the agreement is good, especially at the front surface of the cathode tip ($l\lesssim 1.6$ mm), which is the main part of the arc attachment. Distributions of the surface temperature and density of electric current to the cathode surface (figure 8(a)) agree quite well, similarly to those in figure 7. Distributions of the energy flux to the cathode surface (figure 8(b)) are in good agreement for I  =  60 A, but deviate for 100 A and especially for 200 A. Note that one could think that this deviation is a consequence of the Joule effect inside the cathode, which was taken into account in the NLTE-sheath approach but not in the realization of the model of nonlinear surface heating being employed; however, this has been disproved by additional simulations performed by means of the model of nonlinear heating with account of the Joule effect as in [46].

Values of the density of ion current to the cathode surface (figure 8(b)) given by the NLTE-sheath approach are close to those given by the model of nonlinear surface heating at the centre of the front surface of the cathode, but decrease faster with increasing l. Values of the sheath voltage (figure 8(c)) given by the two models are close to each other at the centre of the front surface of the cathode. As l increases, the values given by the NLTE-sheath model increase weakly, while those given by the model of nonlinear surface heating decrease weakly; however, the difference between the two is not very significant.

In the framework of the NLTE-sheath model, T_e in figure 8(c) represents the electron temperature on the plasma side of the interface describing the sheath. In the nonlinear surface model, T_e has the meaning of the average electron temperature in the ionization layer as in section 1. The two distributions are not very different, except for the maxima occurring in the range $4~\text{mm}\lesssim l\lesssim 6~\text{mm}$ in the NLTE-sheath model, which, however, have little relevance since the current density at these l values is already quite small.

Three approaches to describing separation of charges in near-cathode regions of high-pressure arc discharges are considered. The most straightforward one is the unified modelling, which does not rely on a priori division of the interelectrode gap into quasi-neutral plasma and space-charge sheaths.

The second approach, which may be called the NLTE-sheath approach, is based on separate descriptions of the bulk quasi-neutral plasma and space-charge sheaths formed near solid surfaces (electrodes and insulators). The description of the bulk plasma is fully non-equilibrium, i.e. it does not rely on assumptions of thermal (${{T}_{\text{e}}}={{T}_{\text{h}}}$) or ionization equilibrium. The appropriate description of space-charge sheaths is developed in the appendix.

The third approach employs separate descriptions of the bulk plasma, where deviations between T_e and T_h are taken into account but deviations from ionization equilibrium are not, of the ionization layer and of the near-cathode space-charge sheath. This approach may be called the 2T-ionization layer-sheath approach. Since processes in the bulk have little effect over the cathodic part of the arc, calculation of the cathodic part (the cathode, the sheath and the ionization layer) is decoupled from calculation of the bulk plasma in the framework of this approach. The reduced version, in which only the cathodic part is simulated, is sometimes referred to as the nonlinear surface heating model.

Since the unified modelling has been performed until now only in 1D cases while the NLTE-sheath approach has been realized only for the axially symmetric case, we cannot compare all three approaches at once. Results given by the unified modelling are compared with those given by the model of nonlinear surface heating for the simple 1D test case of a rod cathode with thermally and electrically insulated lateral surface. Results given by the NLTE-sheath approach are compared with those given by the model of nonlinear surface heating for the (axially symmetric) test case of a free-burning atmospheric-pressure argon arc with a rod cathode with a hemispherical tip. It is found that the results given by different models are generally in good agreement, and in some cases the agreement is even surprisingly good.

On one hand, this finding is consistent with estimates given in the appendix (see also section 1 of [28]), which indicate that for typical high-pressure arcs the Debye length is much smaller than the ionization length, which is in turn much smaller than characteristic arc dimensions; hence the necessary conditions of applicability of both the NLTE-sheath and 2T-ionization layer-sheath approaches are fulfilled. On the other hand, the above finding shows that the predicted integral characteristics of the arc-cathode interaction are not strongly affected by details of the model (such as the ion motion in the sheath being collision-free or collision-dominated, variability of T_e and T_h in the near-cathode layer, distribution of the heavy-particle energy flux from the plasma to the cathode surface between the atoms and the ions) provided that the basic physics is right.

The unified modelling approach is at present prohibitively intensive computationally in situations of practical interest that require multidimensional simulations. If the main objective is to simulate the cathodic part rather than the arc as a whole, then it seems natural to employ the model of nonlinear surface heating, which is the simplest one and is ready for use for a wide range of plasma-producing gases and cathode materials (e.g. [44]). This model is a natural first step also in simulations of the arc as a whole, which can be performed by means of either the NLTE-sheath approach or the 2T-ionization layer-sheath approach. The former is the method of choice in cases where deviations from ionization equilibrium occurring in the vicinity of the anode and in the arc fringes are of interest. The current density distribution along the cathode surface, given by the model of nonlinear surface heating, may be used as an initial approximation in this approach, which will contribute to the solution's stability [28].

The 2T-ionization layer-sheath approach may be used in cases where deviations from ionization equilibrium occurring in the vicinity of anode and in the arc fringes are not of primary interest. The solution given by the model of nonlinear surface heating provides boundary conditions at the cathode surface for the bulk plasma equations in this approach.

Up to now, both the NLTE-sheath and the 2T-ionization layer-sheath approaches have been applied only to cases where current transfer to the cathode occurs in the diffuse mode. It will be very interesting to extend the modelling to spot modes. Since the model of nonlinear surface heating allows a routine calculation of spot modes, calculation of corresponding states of the arc bulk is straightforward in the 2T-ionization layer-sheath approach and seems to be feasible in the NLTE-sheath approach. Another effect to be treated in the future is motion of the molten metal in the cathode, which is especially important at high currents.

The work at Universidade da Madeira was supported by FCT—Fundação para a Ciência e a Tecnologia—of Portugal through the projects PTDC/FIS-PLA/2708/2012 and Pest-OE/UID/FIS/50010/2013. The work at INP Greifswald e.V. was supported in part by the DFG (German Science Foundation) under grant UH106/11-1. The collaboration between INP Greifswald e.V. and Universidade da Madeira has been supported in part by funding from the European Union Seventh Framework Programme under grant No. 316216.

The calculation domain in fully non-equilibrium models of arc discharges includes the quasi-neutral plasma and eventually the adjacent solid bodies (electrodes and insulators), separated by interfaces describing the space-charge sheaths. The boundary conditions at these interfaces are formulated in what follows.

The electrostatic potential is discontinuous at the interface:

$$\begin{eqnarray}&&{{\varphi}^{\left(\text{pl}\right)}}-{{\varphi}^{\left(\text{s}\right)}}={{U}_{\text{sh}}}.\end{eqnarray} \tag{ A.1 }$$

Here U_sh has the meaning of the sheath voltage, which varies along the solid surface and has to be found as part of the solution, the upper indices 'pl' and 's' denote values on the quasi-neutral plasma side of the interface and on the side of the solid, respectively. The electric current is continuous at the interface:

$$\begin{eqnarray}&&j_{\text{w}}^{\left(\text{pl}\right)}=j_{\text{w}}^{\left(\text{s}\right)},\end{eqnarray} \tag{ A.2 }$$

where $j_{\text{w}}^{\,\left(\text{pl}\right)}$ and $j_{\text{w}}^{\,\left(\text{s}\right)}$ are the normal components of the densities of electric current from the quasi-neutral plasma to the interface and from the interface into the solid. (If the solid represents a dielectric, then $j_{\text{w}}^{\,\left(\text{s}\right)}=0$ and this boundary condition implies that temporal variations of surface charge on the dielectric are neglected.) In the following, let us replace both designations $j_{\text{w}}^{\,\left(\text{pl}\right)}$ and $j_{\text{w}}^{\,\left(\text{s}\right)}$ with j_w. Again, j_w varies along the surface and has to be found as part of the solution. One can write

$$\begin{eqnarray}&&{{j}_{\text{w}}}=e\left(J_{\text{iw}}^{\left(\text{pl}\right)}-J_{\text{ew}}^{\left(\text{pl}\right)}\right).\end{eqnarray} \tag{ A.3 }$$

Here $J_{\text{iw}}^{\left(\text{pl}\right)}$ and $J_{\text{ew}}^{\left(\text{pl}\right)}$ designate normal components of the number densities of, respectively, the ion and electron fluxes coming to the interface from the quasi-neutral plasma, which are governed by the hydrodynamic (diffusion) equations describing the plasma.

Let us consider as an example a high-current arc in atmospheric-pressure $\text{Ar}$ and assume that ${{T}_{\text{e}}}=15\,000$ K and ${{T}_{\text{h}}}=3500$ K in the near-cathode plasma. An estimate for these conditions gives ${{\lambda}_{\text{D}}}\approx 1.5\times {{10}^{-2}}$ μm for the Debye length and $d\approx 26$ μm for the ionization length. Since ${{\lambda}_{\text{D}}}\ll d$, ionization and recombination are negligible inside the sheath. Hence, the densities of ion and electron fluxes across the sheath are constant and equal to $J_{\text{iw}}^{\left(\text{pl}\right)}$ and $J_{\text{ew}}^{\left(\text{pl}\right)}$, respectively.

If the solid represents a dielectric, ${{j}_{\text{w}}}=0$, or an electrode operating as a cathode, ${{j}_{\text{w}}}>0$, then ${{U}_{\text{sh}}}>0$. If the solid represents an electrode operating as an anode, ${{j}_{\text{w}}}<0$, then U_sh may be either negative or positive. Note that in terms of the so-called near-anode voltage (see, e.g. definition in [14]) ${{U}_{\text{sh}}}<0$ corresponds to positive near-anode voltage and vice versa.

It is natural to neglect variations of the heavy-particle temperature inside the sheath, which gives a boundary condition

$$\begin{eqnarray}&&T_{\text{h}}^{\left(\text{pl}\right)}={{T}^{\left(\text{s}\right)}}.\end{eqnarray} \tag{ A.4 }$$

The energy flux is discontinuous at the interface and the corresponding boundary condition is represented by the equation of integral balance of energy:

$$\begin{eqnarray}&&q_{\text{w}}^{\left(\text{s}\right)}-q_{\text{w}}^{\left(\text{pl}\right)}={{j}_{\text{w}}}{{U}_{\text{sh}}}-{{j}_{\text{w}}}\frac{{{A}_{\text{f}}}}{e}+J_{\text{iw}}^{\left(\text{pl}\right)}{{A}_{\text{i}}}-{{q}^{\left(\text{rad}\right)}},\end{eqnarray} \tag{ A.5 }$$

where $q_{\text{w}}^{\left(\text{pl}\right)}$ and $q_{\text{w}}^{\left(\text{s}\right)}$ are normal components of the energy fluxes coming, respectively, from the quasi-neutral plasma to the interface and from the interface into the solid, A_f is the work function of the electrode material, A_i is the ionization energy and ${{q}^{\left(\text{rad}\right)}}$ is the net density of energy flux removed from the solid surface by radiation. Note that $q_{\text{w}}^{\left(\text{pl}\right)}$ comprises the heat fluxes caused by heat conduction and by the effect inverse to thermal diffusion, the enthalpy transport by the diffusion fluxes, but not transport of production enthalpies and radiation. The term ${{j}_{\text{w}}}{{U}_{\text{sh}}}$ on the rhs of equation (A.5) represents power deposited into the sheath by the arc power supply, ${{j}_{\text{w}}}A/e$ is cooling of the electrode surface due to passage of the electrons from metal into the plasma, and $J_{\text{iw}}^{\left(\text{pl}\right)}{{A}_{\text{i}}}$ is the power deposited at the solid surface due to neutralization of the ions coming from the plasma.

It may be convenient for numerical implementation to introduce a new variable, which equals the temperature in the solids and the heavy-particle temperature in the bulk plasma, and to treat the equation of energy of heavy particles in the quasi-neutral plasma jointly with the equations of heat conduction inside the solids. In this connection, it may be convenient to transform the boundary condition (A.5) with the use of the relation $q_{\text{w}}^{\left(\text{pl}\right)}=q_{\text{hw}}^{\left(\text{pl}\right)}+q_{\text{ew}}^{\left(\text{pl}\right)}$, where $q_{\text{hw}}^{\left(\text{pl}\right)}$ and $q_{\text{ew}}^{\left(\text{pl}\right)}$ are normal components of energy fluxes transported from the quasi-neutral plasma to the interface by heavy particles and electrons, respectively, and to rewrite this boundary condition as

$$\begin{eqnarray}&&q_{\text{w}}^{\left(\text{s}\right)}-q_{\text{hw}}^{\left(\text{pl}\right)}=q_{\text{ew}}^{\left(\text{pl}\right)}+{{j}_{\text{w}}}{{U}_{\text{sh}}}-{{j}_{\text{w}}}\frac{{{A}_{\text{f}}}}{e}+J_{\text{iw}}^{\left(\text{pl}\right)}{{A}_{\text{i}}}-{{q}^{\left(\text{rad}\right)}}.\end{eqnarray} \tag{ A.6 }$$

The form of the rest of the boundary conditions depends on the sign of the sheath voltage. Let us consider first the case of a positive sheath, ${{U}_{\text{sh}}}>0$. In this case, the electrons from the quasi-neutral plasma are repelled by the sheath electric field and follow in the sheath the Boltzmann distribution in space with a temperature equal to $T_{\text{e}}^{\left(\text{pl}\right)}$, the electron temperature in the quasi-neutral plasma estimated at the interface. The density of the flux of plasma electrons counterdiffusing to the surface against the retarding electric field is given by the usual expression

$$\begin{eqnarray}&&{{J}_{\text{CD}}}=\frac{n_{\text{i}}^{\left(\text{pl}\right)}}{4}{{\left(\frac{8kT_{\text{e}}^{(\,pl)}}{\pi {{m}_{\text{e}}}}\right)}^{1/2}}\exp \left(-\frac{e{{U}_{\text{sh}}}}{kT_{\text{e}}^{\left(\text{pl}\right)}}\right).\end{eqnarray} \tag{ A.7 }$$

If the solid represents an electrode and its surface is hot enough, there is a second group of electrons in the sheath: those emitted by the electrode surface. Therefore, in the general case one can write

$$\begin{eqnarray}&&J_{\text{ew}}^{\left(\text{pl}\right)}={{J}_{\text{CD}}}-{{j}_{\text{em}}}/e,\end{eqnarray} \tag{ A.8 }$$

where j_em is the density of electron emission current.

Integral balance of electron energy in the sheath yields

$$\begin{eqnarray}&&\frac{{{j}_{\text{em}}}}{e}2k{{T}^{\left(\text{s}\right)}}-{{J}_{\text{CD}}}2kT_{\text{e}}^{\left(\text{pl}\right)}+q_{\text{ew}}^{\left(\text{pl}\right)}-J_{\text{ew}}^{\left(\text{pl}\right)}e{{U}_{\text{sh}}}=0.\end{eqnarray} \tag{ A.9 }$$

The first term on the lhs represents the density of flux of energy removed from the solid surface by the emitted electrons, evaluated under the conventional approximation of the velocity distribution of the emitted electrons being the half-Maxwellian function with the surface temperature. ${{J}_{\text{CD}}}2kT_{\text{e}}^{\left(\text{pl}\right)}$ is the flux of energy delivered to the surface by the plasma electrons, so the first and second terms jointly represent the net electron energy flux from the solid to the sheath. The third and fourth terms represent, respectively, the electron energy flux from the quasi-neutral plasma to the sheath and the work of the electric field over the electrons inside the sheath.

Lowering of the work function due to the presence of an electric field (the Schottky effect) is appreciable for arc electrodes if ${{U}_{\text{sh}}}>0$. It is logical to assume in such conditions that the energy gained by an emitted electron on crossing the space-charge sheath is equal to $e{{U}_{\text{sh}}}-\Delta A$, where $\Delta A$ is the Schottky correction. The energy gained (lost) by a singly charged ion (plasma electron) moving to the electrode is the same. Therefore, in the following we will assume that the factor U_sh in the last term on the lhs of equation (A.9) is replaced by ${{U}_{\text{sh}}}-\Delta A/e$. Note that there is no need for such replacement in equation (A.5) or (A.6), where U_sh appears in the combination $e{{U}_{\text{sh}}}-{{A}_{\text{f}}}$.

Relaxation of the two-group electron distribution function in the sheath to the Maxwellian distribution in the bulk plasma occurs on the length scale ${{\lambda}_{\text{ee}}}$, the mean free path for electron–electron collisions. Under the above-mentioned conditions of a high-current arc in atmospheric-pressure $\text{Ar}$, ${{\lambda}_{\text{ee}}}=0.43$ μm. Hence, another zone distinguished by specific physics is relevant: the layer of electron relaxation, which is much thicker than the sheath but much thinner than the ionization length. The fact that equations (A.8) and (A.9) involve electron particle and energy fluxes to the interface, $J_{\text{ew}}^{\left(\text{pl}\right)}$ and $q_{\text{ew}}^{\left(\text{pl}\right)}$, which are governed by the hydrodynamic (diffusion) equations, implies that the layer of electron relaxation is implicitly included in the interface jointly with the space-charge sheath.

The number density of charged particles in the sheath is much smaller than that in the quasi-neutral plasma, so the boundary condition for the distribution of the charged particles in the quasi-neutral plasma to be applied at the interface could be written simply as $n_{\text{i}}^{\left(\text{pl}\right)}=0$. However, we consider the case where the sheath is collision-free for ions, and in this case the use of this boundary condition would require one to expressly treat the Knudsen layer, i.e. a region of quasi-neutral plasma that is adjacent to the sheath and has a thickness of the order of the mean free path for ion–neutral collisions, where the ions are accelerated up to the Bohm speed ${{u}_{\text{B}}}=\sqrt{k\left(T_{\text{e}}^{\left(\text{pl}\right)}+T_{\text{h}}^{\left(\text{pl}\right)}\right)/{{m}_{\text{i}}}}$ (in other words, the Knudsen layer plays the role of a presheath). In order to avoid such treatment, one can assume that the acceleration of the ions up to the Bohm speed may be described to a reasonable accuracy by the hydrodynamic (diffusion) equations and can use at the interface the boundary condition (e.g. [29])

$$\begin{eqnarray}&&n_{\text{i}}^{\left(\text{pl}\right)}=J_{\text{iw}}^{\left(\text{pl}\right)}/{{u}_{\text{B}}}.\end{eqnarray} \tag{ A.10 }$$

The electric field at the electrode surface needed for evaluation of the Schottky correction may be estimated with the use of equation (10) of [29] (which was derived for the case where the sheath is collision-free for ions). Then relations (A.1) and (A.2); (A.4) and (A.5); (A.9); and (A.10) represent a complete set of boundary conditions at the interface for the differential equations of, respectively, current continuity in the quasi-neutral plasma and the electrode; heavy-particle energy in the plasma and heat conduction in the solid; electron energy in the plasma; and ambipolar diffusion in the plasma. Since the sheath voltage U_sh remains indeterminate in this reasoning, an additional relationship is needed: equation (A.8).

Let us proceed to the case of a negative sheath, ${{U}_{\text{sh}}}<0$. In this case, the ions in the sheath follow the Boltzmann distribution with a temperature $T_{\text{h}}^{\left(\text{pl}\right)}$ and one can write, similarly to equation (A.7),

$$\begin{eqnarray}&&J_{\text{iw}}^{\left(\text{pl}\right)}=\frac{n_{\text{i}}^{\left(\text{pl}\right)}}{4}{{\left(\frac{8kT_{\text{h}}^{\left(\text{pl}\right)}}{\pi {{m}_{\text{i}}}}\right)}^{1/2}}\exp \left(\frac{e{{U}_{\text{sh}}}}{kT_{\text{h}}^{\left(\text{pl}\right)}}\right).\end{eqnarray} \tag{ A.11 }$$

Most of the electrons eventually emitted by the solid cannot overcome the potential barrier in the sheath and immediately return to the surface. Therefore, electric current transported by the emitted electrons can be neglected. In particular, one should drop the term involving j_em and set ${{J}_{\text{CD}}}=J_{\text{ew}}^{\left(\text{pl}\right)}$ in equation (A.9) (of course, there is no need to replace U_sh with ${{U}_{\text{sh}}}-\Delta A/e$ in this equation in the case ${{U}_{\text{sh}}}<0$).

The electrons from the quasi-neutral plasma enter the sheath with the electron Bohm speed ${{u}_{\text{eB}}}=\sqrt{k\left(T_{\text{e}}^{\left(\text{pl}\right)}+T_{\text{h}}^{\left(\text{pl}\right)}\right)/{{m}_{\text{e}}}}$; see discussion in [47]. Therefore, equation (A.10) is replaced with

$$\begin{eqnarray}&&n_{\text{i}}^{\left(\text{pl}\right)}=J_{\text{ew}}^{\left(\text{pl}\right)}/{{u}_{\text{eB}}}.\end{eqnarray} \tag{ A.12 }$$

Equation (A.9) with the above-described changes and equation (A.12) represent boundary conditions at the interface for the differential equations of electron energy and ambipolar diffusion in the plasma in the case ${{U}_{\text{sh}}}<0$. The additional relationship is supplied by equation (A.11).

Note that in finite-element codes such as COMSOL Multiphysics the Dirichlet boundary condition for the ambipolar diffusion equation may be less favorable from the point of view of numerical stability than the Neumann condition. In this connection, it may be appropriate to use as a boundary condition for the ambipolar diffusion equation in the plasma not equation (A.12) as it is, but rather a combination of this equation with (A.11):

$$\begin{eqnarray}&&J_{\text{iw}}^{\left(\text{pl}\right)}=\frac{J_{\text{ew}}^{\left(\text{pl}\right)}}{4{{u}_{\text{eB}}}}{{\left(\frac{8kT_{\text{h}}^{\left(\text{pl}\right)}}{\pi {{m}_{\text{i}}}}\right)}^{1/2}}\exp \left(\frac{e{{U}_{\text{sh}}}}{kT_{\text{h}}^{\left(\text{pl}\right)}}\right).\end{eqnarray} \tag{ A.13 }$$