The Solution Structure of the Elenbaas–Heller Equation for Inductively Coupled Plasmas and Wall Stabilized Arcs

1. Introduction

Radio frequency (rf) inductively coupled plasmas (ICP) since their first practical implementation by Reed [1] in 1961 are continuously improving and are encountered in many industrial and scientific applications [2,3]. A salient feature of ICP is the generation of a plasma torch in an electrodeless manner, free of particle contamination due to electrode erosion, via electro-magnetic induction [4,5,6]. Therefore, ICP sources are widely used in a variety of applications such as synthesis of nanopowders and carbon nanotubes [7], semiconductor technology [8,9], plasma spray processes [3,10,11], ground testing facilities for simulating atmospheric entry vehicles thermal loads (hypersonics) [12,13], ion thrusters for space propulsion [14,15], environmental [16,17] and analytical chemistry [18,19,20] which is the largest area of application for this technology.

Experimental methods for measuring the ICP parameters are time consuming and require specialized and expensive instrumentation so numerical modeling and simulation is an indispensable tool in understanding key features of the plasma flow and optimizing existing devices and developing new. Detailed numerical simulation of ICP requires the simultaneous solution of Navier-Stokes, energy and Maxwell’s equations coupled through nonlinear transport and kinetic coefficients. Over the years several numerical models have been developed for the description of the chemical processes and the magneto-hydrodynamic phenomena inside the torch; from one-dimensional [4,21,22,23,24,25] to multi-dimensional configurations [13,26,27,28,29]. Most of the simulations assume local thermodynamic equilibrium (LTE) which is a very good approximation for ICP torches, since at atmospheric pressure the collisional frequency between the gas particles is sufficiently high to maintain local equilibrium. Higher accuracy may be obtained by non-LTE models with state-to-state resolution of the internal energy states at the atomic level (translational, rotational, electronic, vibrational modes) but are computationally very expensive [13,30,31,32,33].

The governing equations for the plasma torch are complex and highly nonlinear. As a result bifurcations [34], instabilities [35], multiple solutions [36,37], and pattern formation [38] is a characteristic of the numerical solutions. For example the Elenbaas–Heller equation [39,40] describing the temperature profile in a wall stabilized arc with cylindrical symmetry [41,42,43,44,45,46,47] admits multiple solutions. Phillips [48] essentially calculated multiple temperature profiles when switching from an initial dc to an ac excitation of the arc column. Lowke [41] obtained a multivalued current-voltage characteristic for an air plasma and explained the multiple solutions on the non-linear and non-monotonic behavior of the ratio between the radiation losses and the electric conductivity. Recently Zhovtyansky et al. [47] calculated multiple solutions using the Elenbaas-Heller model for an air-Cu plasma. The multiplicity was explained on the basis of the non-monotonic thermal conductivity–temperature curve. For an inductively coupled argon plasma Mensing and Boedeker [4] calculated two solutions one stable and one unstable. Similar nonlinear phenomena are observed in electrodeless radiofrequency discharges which exhibit bistability i.e. two modes of operation: a higher density mode which is an inductive discharge known as the H-mode and a low-density known as the E-mode in which the power is capacitively coupled to the plasma. As a result of the multiple steady states due to the multipoint equilibrium between absorbed and dissipated power the transition between these modes is hysteretic, i.e. the E to H transition occurs at a different current than the reverse H to E transition [49,50,51,52,53].

The aim of the present work is to present a systematic analysis of the solution structure for ICP torches and wall stabilized arcs and it is organized as follows: In Section 2 the electrothermal model for the ICP and the wall stabilized arc are developed in dimensionless form. Multiplicity analysis is carried out in Section 3, first for the simpler case of wall stabilized arcs and then for the ICP. In Section 3.1.1 argon plasma is considered as the working gas and in Section 3.1.2 hydrogen. The solution structure for argon and hydrogen ICP is presented in Section 3.2.

2. Analysis

The operating principle of an rf-ICP is schematically shown in Figure 1. A rf current passes through an externally wounded coil around a dielectric tube. The flowing gas, argon or hydrogen in our case, absorbs the electromagnetic power supplied from the coil and it is subsequently heated and ionized (torch generation). In the following analysis it is assumed that the length of the dielectric tube is sufficient for the flow to be fully developed, while the plasma is assumed to be optically thin and in LTE condition.

2.1. Electromagnetic Field

Maxwell’s equations for a conductor with negligible displacement current and unit permeability are given by [54]

$$\nabla \times \mathbf{E}+{\mu }_0{\displaystyle \frac{\partial \mathbf{H}}{\partial t}}=0,$$

(1)

and

$$\nabla \times \mathbf{H}-\sigma \left(T\right)\mathbf{E}=0,$$

(2)

where $\sigma \left(T\right)$ is the temperature dependent electric conductivity. Expanding and considering the cylindrical symmetry for the electric field $\mathbf{E}(\rho ,\theta ,z)=\left(0,E_{\theta },0\right)$ with $E_{\theta }=E_{\theta }\left(\rho \right)$ and the magnetic field $\mathbf{H}(\rho ,\theta ,z)=\left(\mathrm{0,0},H_z\right)$ with $H_z=H_z\left(\rho \right)$ yields

$${\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial \left({\rho E}_{\theta }\right)}{\partial \rho }}+{\mu }_0{\displaystyle \frac{\partial H_z}{\partial t}}=0$$

(3)

and

$${\displaystyle \frac{\partial H_z}{\partial \rho }}-\sigma \left(T\right)E_{\theta }=0$$

(4)

A time harmonic electro-magnetic field is assumed in the form [4]

$$\begin{array}{cc}{E}_{\theta }=E_p{e}^{i(\omega t-{\varphi }_E)},& {\varphi }_E={\varphi }_E\left(\rho \right)\end{array}$$

(5)

$$\begin{array}{cc}{H}_z=H_p{e}^{i(\omega t-{\varphi }_H)},& {\varphi }_H={\varphi }_H\left(\rho \right)\end{array}$$

(6)

where $E_p$ and $H_p$ are the amplitudes of the electric and magnetic field respectively and ${\varphi }_E$ and ${\varphi }_H$ are the corresponding phase angles. Substituting Eqs.(5,6) into Eqs.(3,4) and separating real and imaginary parts, the following system of equations is obtained:

$$\begin{array}{c}{\displaystyle \frac{d{H}_p}{d\rho }}=\sigma \left(T\right)E_p\mathrm{c}\mathrm{o}\mathrm{s}({\varphi }_H-{\varphi }_E)\\ -H_p{\displaystyle \frac{d{\varphi }_H}{d\rho }}=\sigma \left(T\right)E_p\mathrm{s}\mathrm{i}\mathrm{n}({\varphi }_H-{\varphi }_E)\\ \begin{array}{c}{\displaystyle \frac{d{E}_p}{d\rho }}+{\displaystyle \frac{E_p}{\rho }}=-\omega {\mu }_0{H}_p\mathrm{s}\mathrm{i}\mathrm{n}({\varphi }_H-{\varphi }_E)\\ E_p{\displaystyle \frac{d{\varphi }_E}{d\rho }}=\omega {\mu }_0{H}_p\cos \left({\varphi }_H-{\varphi }_E\right).\end{array}\end{array}$$

(7)

The boundary conditions are:

$$\begin{array}{ccc}{H}_p\left(0\right)=H_c,& E_p\left(0\right)=0,& \begin{array}{cc}{\varphi }_H=0,& {\varphi }_E=\pi /2\end{array}\end{array}.$$

(8)

For the wall stabilized arc the rf coil is absent since a dc electric field is established between the anode and cathode electrodes. Expanding Eq.(3)

$$\nabla \times \mathbf{E}=\left({\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial E_z}{\partial \theta }}-{\displaystyle \frac{\partial E_{\theta }}{\partial z}}\right){\mathbf{a}}_{\rho }+\left({\displaystyle \frac{\partial E_{\rho }}{\partial z}}-{\displaystyle \frac{\partial E_z}{\partial \rho }}\right){\mathbf{a}}_{\theta }+\left({\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial \left({\rho E}_{\theta }\right)}{\partial \rho }}-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial E_{\rho }}{\partial \theta }}\right){\mathbf{a}}_z=0,$$

(9)

and taking into consideration the symmetry we are left with $\partial E_z/\partial \rho =0$ so $E=E_z=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$ The current–voltage relationship is determined from the temperature field and the below integral equation:

$$I=E{\int }_0^{R}2\pi r\sigma \left(T\right)dr.$$

(10)

2.2. Temperature Field

For the fully developed flow the energy balance between heat generation (Joule heating) and heat dissipation by conduction and radiation is given by

$${\displaystyle \frac{1}{\rho }}{\displaystyle \frac{d}{d\rho }}\left[\rho k\left(T\right){\displaystyle \frac{dT}{d\rho }}\right]=Q\left(T\right)-\sigma \left(T\right)E^2,$$

(11)

where $Q=4\pi {ϵ}_N$ is the radiative heat flux density, ${ϵ}_N$ is the net emission coefficient and $k$ is the thermal conductivity. Equation (11) is also known as the Elenbaas–Heller equation [39,40]. For the ICP $E=E_{\theta }$ and for the wall stabilized arc $E=E_z$. The boundary conditions at the axis of symmetry and the wall are:

$$\begin{array}{cc}{\left.{\displaystyle \frac{\partial T}{\partial \rho }}\right|}_{\rho =0}=0,& T\left(R\right)=T_w\end{array}.$$

(12a,b)

2.3. The Electro-Thermal Problem in Dimensionless Form

Introducing the below dimensionless variables

$$\begin{array}{ccc}r={\displaystyle \frac{\rho }{R}},& \mathsf{\Theta }={\displaystyle \frac{T}{T_{\mathrm{r}\mathrm{e}\mathrm{f}}}},& \begin{array}{ccc}\lambda ={\displaystyle \frac{k}{k_{\mathrm{r}\mathrm{e}\mathrm{f}}}},& s={\displaystyle \frac{\sigma }{{\sigma }_{\mathrm{r}\mathrm{e}\mathrm{f}}}},& \begin{array}{cc}q={\displaystyle \frac{Q}{Q_{\mathrm{r}\mathrm{e}\mathrm{f}}}},& \begin{array}{cc}e={\displaystyle \frac{E_p}{E_{\mathrm{r}\mathrm{e}\mathrm{f}}}}& h={\displaystyle \frac{H_p}{H_{\mathrm{r}\mathrm{e}\mathrm{f}}}}\end{array}\end{array}\end{array}\end{array},$$

(13)

the electromagnetic equations take the form:

$$\begin{array}{c}{\displaystyle \frac{dh}{dr}}=ues\left(\mathsf{\Theta }\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\delta \varphi \right)\\ -h{\displaystyle \frac{d{\varphi }_H}{dr}}=ues\left(\mathsf{\Theta }\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\delta \varphi \right)\\ \begin{array}{c}{\displaystyle \frac{de}{dr}}+{\displaystyle \frac{e}{r}}=-uh\mathsf{\Omega }\mathrm{s}\mathrm{i}\mathrm{n}\left(\delta \varphi \right)\\ e{\displaystyle \frac{d{\varphi }_E}{dr}}=uh\mathsf{\Omega }\cos \left(\delta \varphi \right),\end{array}\end{array}$$

(14)

where $\delta \varphi ={\varphi }_H-{\varphi }_E$ and $\mathsf{\Omega }=\omega /{\omega }_{\mathrm{r}\mathrm{e}\mathrm{f}}$ with ${\omega }_{\mathrm{r}\mathrm{e}\mathrm{f}}=Q_{\mathrm{r}\mathrm{e}\mathrm{f}}/{\mu }_0{\left(\sigma kT\right)}_{\mathrm{r}\mathrm{e}\mathrm{f}}$. The energy equation becomes:

$${\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}\left[r\lambda \left(\mathsf{\Theta }\right){\displaystyle \frac{d\mathsf{\Theta }}{dr}}\right]=u^2\left[q\left(\mathsf{\Theta }\right)-s\left(\mathsf{\Theta }\right)e^2\right],$$

(15)

The dimensionless boundary conditions for the electromagnetic field are now:

$$\begin{array}{ccc}h\left(0\right)=h_c,& e\left(0\right)=0,& \begin{array}{cc}{\varphi }_H=0,& {\varphi }_E=\pi /2\end{array}\end{array},$$

(16)

and for the temperature field

$$\begin{array}{cc}{\left.{\displaystyle \frac{d\mathsf{\Theta }}{dr}}\right|}_{r=0}=0,& \mathsf{\Theta }\left(1\right)={\mathsf{\Theta }}_w\end{array}.$$

(17a,b)

An important parameter that appears in Eqs.(14,15) is the conduction-radiation parameter (CRP) which is defined as

$$\begin{array}{c}u=R{\left({\displaystyle \frac{Q}{kT}}\right)}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{1/2}=R{\left({\displaystyle \frac{\sigma E}{H}}\right)}_{\mathrm{r}\mathrm{e}\mathrm{f}},\end{array}$$

(18)

where $H_{\mathrm{r}\mathrm{e}\mathrm{f}}={\left(\sigma kT\right)}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{1/2}$ and $E_{\mathrm{r}\mathrm{e}\mathrm{f}}={\left(Q/\sigma \right)}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{1/2}$. CRP is a measure of the heat transferred by radiation to the heat transferred through conduction in the radial direction and it is encountered in electrothermal models of metallic conductors [56] and superconductors [55].

3. Results and Discussion

The two-point boundary value problem described by Eqs.(14,15) has been solved numerically with a multi-shooting algorithm, Ascher et al. [57]. Continuation along the various branches and the calculation of the singular points has been carried out along the lines suggested by Seydel [58]. For argon the thermal conductivity was fitted from tabulated data from Boulos et al. [2], the electric conductivity from Cressault et al. [59] and the net emission coefficient from Cressault and Gleizes [60]. For hydrogen the thermal and the electric conductivity was adopted from Boulos et al. [2] and the radiation correlation from Gueye et al. [46]. For the scaling of the variables the following reference values were used: $T_{\mathrm{r}\mathrm{e}\mathrm{f}}=1000$ [K], $k_{\mathrm{r}\mathrm{e}\mathrm{f}}=5.87\times {10}^{-2}$ [W/mK], ${\sigma }_{\mathrm{r}\mathrm{e}\mathrm{f}}=1$ [S/m] and $Q_{\mathrm{r}\mathrm{e}\mathrm{f}}=1$ [MW/m³]. To get a better understanding of the solution structure the simpler case of the wall stabilized arc, where the electric field is uniform along the radius, will be considered first for argon and hydrogen plasmas respectively.

3.1. Wall Stabilized Arcs

3.1.1. Argon

Selecting $e$ as the bifurcation parameter, the projection of the solution structure on the $\left(e,{\mathsf{\Theta }}_c\right)$ plane for argon when the radiation is negligible, is shown in Figure 2. For a certain range of the field intensity $e$ up to three solutions exist denoted as A, B, and C. For instance, for $e=0.46$, A is the cold branch, B is the intermediate, and C is the hot one. Solution A is in fact a trivial one since the flow is isothermal, determined by the wall temperature and the gas in this state is rather an insulator (no plasma). Thus, a very high field is required to increase the electric conductivity for the current to start flowing, just after passing the first singular point $L_1$. This can be clearly seen on the current-voltage characteristic curves in Figure 3, where the field intensity is rapidly decreasing while the current is increasing as soon as the ionization starts after $L_1$. As expected, the lower the wall temperature the higher the field required to break through the singular point $L_1$. Although the cold solution is of no practical interest, the continuation algorithm is robust and able to capture it. After the second singular point $L_2$ the wall temperature has no influence on the center temperature and the current-voltage characteristics as all curves are merging into one.

The multiplicity structure when radiation is taken into consideration is shown in Figure 4 on the $\left(e,{\mathsf{\Theta }}_c\right)$ plane and in Figure 5 on the $\left(i,e\right)$ plane respectively. The radiation does not affect the bifurcation pattern as it does not introduce additional solutions, and the solution structure remains essentially the same. Thus, the existence of multiple solutions is the result of the non-linear and non-monotonic dependence of the thermal conductivity on the temperature as shown on the right of Figure 2, where the singular points correlate to the local inflexion points on the $k-T$ diagram. The comparison with the experimental data shown in Figure 6 for $R=1$cm [61,62] and $R=2$cm [61,62] and in Figure 7 for $R=0.5$cm [61,62,63] and $R=0.25$cm [64] suggest that radiation is significant after the second singular point $L_2$ where higher temperatures are encountered. It is worth pointing out that despite the simplicity of the theoretical model the overall agreement with the measurements is very good, inspiring confidence in the methods and in the results of the bifurcation analysis. Figure 8 shows the three solutions denoted as A, B and C and the corresponding temperature profiles for the non-radiating and the radiating cases. As expected, radiation affects the hot branch C, where a significant reduction in the temperature and the current is observed because of the heat losses towards the wall. While the temperature profile in the intermediate branch B may be approximated by a parabolic curve the hot branch profile is highly nonlinear with an almost constant temperature in the conducting core which abruptly falls to the wall temperature within a thin non-conducting boundary layer.

3.1.2. Hydrogen

Let us now consider the more complicated case of a hydrogen arc column as depicted in Figure 9 without radiation with ${\mathsf{\Theta }}_w$ as a parameter. Comparing with Figure 2 for argon, on the $\left(e,{\mathsf{\Theta }}_c\right)$ diagram there are now four singular points $L_1$ to $L_4$. For a certain range of the field intensity $e$ up to four solutions exist, the one with lower temperature corresponding to the trivial (no plasma) solution. The multiplicity is clearly the result of the non-monotonic thermal conductivity as the singular points appear in the neighborhood where strong gradients in the $k-T$ curve exist. The function $k=k\left(T\right)$ is also multivalued. For instance, the equation $k\left(T\right)=5$ [W/mK] is satisfied for four different temperatures simultaneously. As the heat dissipation due to conduction is proportional to the thermal conductivity the balance between Joule heating (heat generation) and heat dissipation is satisfied at multiple points. The complicated structure of current-voltage characteristic curve is shown in Figure 10. The existence of multiple solutions for the Elenbaas-Heller equation has been reported by Phillips [48] who essentially calculated multiple temperature profiles when switching from an initial dc to an ac excitation of the arc column. Lowke [41] obtained a similar current-voltage characteristic as the one in Figure 10 for an air plasma and explained the multiple solutions on the non-linear and non-monotonic behavior of the ratio between the radiation losses and the electric conductivity $4\pi {ϵ}_N/\sigma$ (Figure 13 in [41]). However, the present study clearly showed that multiple solutions is a characteristic feature of both radiating and non-radiating arcs. Recently Zhovtyansky et al. [47] calculated multiple solutions using the Elenbaas-Heller model for an air-Cu plasma. The multiplicity was explained on the basis of the non-monotonic $k-T$ curve which exhibits several peaks as the one for hydrogen in Figure 9 (Figure 3 in [47]). On the other hand, this does not necessarily mean that radiation does not play a role in the multiplicity. For instance, the net emission coefficient for hydrogen is also a non-monotonic function of the temperature as can be seen from Figure 4 in [46], where after a peak value at high temperatures ${ϵ}_N$ decreases. The effect is shown in the high current regime of Figure 11 in the present study, where the current-voltage curve becomes non-monotonic and multivalued as another singular point appears at the “tail” of the curve. It is worth comparing with Lowke’s [41] Figure 12 where a monotonic current-voltage curve has been calculated in the high temperature and high current regimes.

3.2. Inductively Coupled Plasmas

For this case the complete set of Eqs. (14,15) is solved with $h_c$ as the bifurcation parameter. The projection of the solution with radiation effects on the $(h_c,{\mathsf{\Theta }}_c)$ plane for argon is shown in Figure 12. There exist three solutions separated by the two singular points $L_1$ and $L_2$. The corresponding temperature and magnetic field profiles are shown in Figure 13 in agreement with the measurements of Eckert and Pridmore–Brown [24]. The one with the lowest temperature (A) is the trivial one (no plasma, uniform magnetic field). For argon the existence of multiple solutions has been recognized by Mensing and Boedeker [4], who reported profiles (B) and particularly (C) with the characteristic off-center peak in the temperature, in detail. The multiplicity structure becomes much more complicated compared to argon, when hydrogen is used as the working gas for the plasma torch as can be seen from Figure 14 where the solution is projected on the $(h_c,{\mathsf{\Theta }}_c)$ plane with the wall temperature as a parameter. The temperature and the magnetic field profiles are presented in Figure 15 for selected values of the center line magnetic field $h_c$ (points A to G). To the right side of Figure 14 at $h_c={10}^6$ the curve ${\mathsf{\Theta }}_w=3$ admits four solutions, points A to D, with point A being the trivial one. As the wall temperature increases $L_1$ remains but the two other singular points $L_2$ and $L_3$ disappear. This is because the temperature field is restricted above the first major peak of the thermal conductivity at about 3.8kK and consequently the number of solutions are reduced. Interestingly another limit point $L_4$ appears to the left at about $h_c=700$ where all curves seem to converge as shown in the inset in Figure 14. The detail also suggests that there exist another solution branch, practically attached to G–F section but it is very difficult to numerically distinguish. Along the path B–G a two–zone structure is observed. A high temperature core surrounded by a low temperature section with an abrupt temperature drop. The temperature increases as the core expands outwards up to point G where it becomes dominant with an off–center temperature maximum (Figure 15).

Figure 13. Temperature profiles for inductively coupled argon plasma. The three solutions correspond to points A, B and C shown in Figure 12.

Figure 13. Temperature profiles for inductively coupled argon plasma. The three solutions correspond to points A, B and C shown in Figure 12.

Figure 14. Bifurcation diagram on the $\left(h_c,{\mathsf{\Theta }}_c\right)$ plane for inductively coupled hydrogen plasma with radiation effects.

Figure 14. Bifurcation diagram on the $\left(h_c,{\mathsf{\Theta }}_c\right)$ plane for inductively coupled hydrogen plasma with radiation effects.

Figure 15. Temperature profiles for inductively coupled hydrogen plasma. Center line magnetic field values $h_c$ A to G, according to Figure 14.

Figure 15. Temperature profiles for inductively coupled hydrogen plasma. Center line magnetic field values $h_c$ A to G, according to Figure 14.

4. Summary and Conclusions

A detailed multiplicity analysis has been carried out for inductively coupled plasmas and wall stabilized arcs. Two working gases have been considered, argon and hydrogen. For the bifurcation analysis a one–dimensional model for the electromagnetic and the temperature field has been adopted. The agreement with the experimental data is good. An important finding is that the solution structure depends on the working gas and it is primarily affected by the particular shape (non–monotonicity) of the thermal conductivity versus temperature curve. As a result for a certain range of the parameters, three solutions (one being trivial, no plasma) have been identified for argon. For hydrogen a more complicated solution structure emerges, up to four solutions, because of the highly non–linear and non–monotonic thermal conductivity curve.