Numerical investigation of transient, low-power metal vapour discharges occurring in near limit ignitions of flammable gas
Rajiv Shekhar Sergey Gortschakow Holger Grosshans Udo Gerlach Dirk Uhrlandt
physics.app-ph











]Numerical investigation of transient, low-power metal vapour discharges occurring in near limit ignitions of flammable gas

^1Physikalisch-Technische Bundesanstalt (PTB), Bundesallee 100, 38116 Braunschweig, Germany
^2School of ITEE, The University of Queensland, Qld, Australia
^3Leibniz Institute for Plasma Science and Technology (INP), 17489 Greifswald, Germany
rajiv.shekhar@uqconnect.edu.au




  August 2017



This article presents an investigation of a transient (30s – 5ms) electrical discharge in metal vapour with low voltage (≤ 30 V) and current (≤ 200 mA), drawn between two separating electrodes. Discharges of this type are rarely studied, but are important in electrical explosion safety, as they can ignite flammable gasses. An empirical model is developed based on transient recordings of discharge voltages and currents and high speed broadband image data. The model is used for predicting the electrical waveforms and spatial power distribution of the discharge. The predicted electrical waveforms show good accuracy under various scenarios. To further investigate the underlying physics, the model is then incorporated into a simplified 3-D gas dynamics simulation including molecular diffusion, heat transfer and evaporation of metal from the electrode surface. The local thermodynamic equilibrium (LTE) assumption is next used to calculate electrical conductivity from the simulated temperature fields, which in turn is integrated to produce electrical resistance over time. This resistance is then compared to that implied by the voltage and current waveforms predicted by the empirical model. The comparison shows a significant discrepancy, yielding the important insight that the studied discharge very likely deviates strongly from LTE.



52.80, 64.70,  47.11
















§ INTRODUCTION


Electrical discharges have been studied extensively as a means of igniting flammable gasses in the field of combustion science. Previous works have typically focussed on high voltage spark discharges, possessing substantially different characteristics to the type of discharge considered here. Although scientific literature on the subject is limited, these low power discharges, referred to variously as break arcs, break sparks, or contact arcs, are well known in electrical explosion safety. Industry standards prescribe methods for managing the explosion risk they pose <cit.>, but require experiment based methods for establishing safety limits which have been proven to be unreliable <cit.>. There is thus an increasing demand for a reliable alternative based on sound science to be developed, which requires a fundamental understanding of the relevant phenomena.

The specific scenario being investigated is described in figure <ref>. A tungsten wire anode is positioned with its tip in contact with the surface of a cadmium block cathode. The wire moves along the surface until it reaches the edge of the block, where contact is broken. As the two electrodes separate, the discharge is drawn between them, reaching a typical maximum length of around 300m. The electrode materials and geometry are based on a reference scenario for standardised explosion safety testing <cit.>. Additionally, the electrodes are energised by an electronic circuit which adjusts output voltage so that a programmed constant current flow is maintained.



A previous investigation applied a combination of optical diagnostic experiments and reactive computational fluid dynamics simulation to the problem, focussing on the combustion resulting from the discharge <cit.>. Here, the discharge was represented by an energy source term derived from an empirical model of the discharge. The empirical model was based on a quasi-static correlation between voltage, current and length of the discharge. This work presents an expanded derivation of this empirical model, including several improvements and the incorporation of new experimental data. A partial validation is also presented utilising electrical measurement data. A key challenge of the empirical approach is quantifying how power is spatially distributed within the discharge. This is currently done using high speed broadband image data, which is difficult to analyse. Additionally, this aspect of the model is not easy to validate.

The direct imposition of a spatial power distribution could be avoided if electrical conductivity within the discharge region were known. This calculation would be significantly simplified by the assumption of local thermodynamic equilibrium (LTE), i.e.: that all species present at any point in the plasma have kinetic energies described by a Maxwell-Boltzmann distribution for a given single temperature. This has the important implication that the parameters of the plasma such as ionisation fraction and electrical conductivity are defined uniquely by this temperature. Previous studies of high voltage sparks conducted simulations under the LTE assumption, where tabulated data for plasma parameters as functions of temperature were used to couple the equations of gas dynamics to those of electromagnetics <cit.>.

Applying a similar approach to modelling the low power discharge in question would require some indication of whether the LTE assumption is valid. In models of discharges in switching electrodes at high currents, the bulk of the plasma is often considered to be in LTE <cit.>. Quasi-static empirical models are also used for equilibrium arcs, and have shown some success in this application <cit.>. The radiation emitted by the discharge was also confirmed to comprise predominantly of atomic lines of cadmium <cit.>.

Characterising this discharge as an equilibrium arc would, however, be questionable, considering the extremely small energy and time scales involved. This work therefore additionally aims to provide insight into this matter. Given the substantial difficulties involved in modelling a non-equilibrium plasma, it is desirable to know beforehand whether doing so is necessary. For this purpose, the developed empirical model is used to describe the energy input of the discharge into a volume of gas. This model reflects measurements of the discharge and does not make assumptions with regard to the underlying physics. Next, assuming LTE, fields for temperature and electrical conductivity are obtained from the simulation. The latter is then integrated spatially to calculate electrical resistance, which can be compared to that predicted directly by the empirical source model. The consistency (or otherwise) between these two calculations provides an indication of the LTE assumption's validity.

The article is presented as follows: section <ref> describes the development of the empirical discharge model, section <ref> describes the simulated solid, gas phase and LTE plasma physics, and section <ref> presents simulation results and comparisons with the empirical model. Conclusions from the work and avenues for future development are presented in sections <ref> and <ref>.



§ EMPIRICAL DISCHARGE MODEL


The empirical part of the model defines the total power input into the discharge over its duration, as well as how that power is distributed over the discharge region. This section outlines the derivation of these two relationships. A similar approach to Ref. <cit.> is taken, applying it to the recent dataset of Ref. <cit.> and adding a more detailed consideration of varying arc geometry. The dataset consists of voltage and current waveforms, together with corresponding high speed image recordings from which length variation of the discharge over time is measured. The data is thus a series of voltage/current/length triples – one for each recorded time instant of each discharge.

An example of this data for one discharge event is shown in figure <ref>. As explained in section <ref>, the current remains roughly constant over the duration of the discharge as regulated by the electronic source circuit, while the voltage and length of the arc increase continuously, the former up to a maximum of 30V. Data for four different current values (70, 100, 150 and 250mA) are used in the proceeding analyses. Images are recorded every 10 , of which five are shown. A voltage and current value is assigned to each image as shown.





 §.§ Voltage and current

The quasi-static relationship between voltage (v), current () and length (l) in m of the discharges takes the form of (<ref>) <cit.>.

    v(i,l) =  +α l (1+β/^n);   < t < 
    l(t)= t

The length is modelled as increasing linearly over the duration of the discharge ( to , as shown in figure <ref> with velocity ). The quantity 
is known as the fall voltage, and is the minimum voltage for which the discharge can exist. Fall voltage is a property of the electrode materials, and is given a constant value of 10 V based on experimental observation. As this value compares favourably to cathode fall voltage values for cadmium given in the literature <cit.>, it is assumed that  consists entirely of a cathode fall and does not arise from anode effects. The remaining parameters α, β and n are fitted from the data.

Discharge lengths are determined using simple image analysis techniques. Here, and for the subsequent analyses of section <ref>, it is assumed for simplification that the light intensity of the discharge recorded over the spectral sensitivity of the high-speed camera (mainly the optical range) scales in good approximation with the power density of the discharge, and therefore indicates the active discharge area. Details like the change of atomic and ionic radiation intensities with species temperatures in the plasma are not considered.
In the first step of the analysis, the greyscale image of the discharge is converted to monochrome, setting to white all parts of the image with intensity greater than 20% of the maximum value. Next, small white areas of the image are removed using a non-linear filtering operation known as morphological opening <cit.>. The largest contiguous white region of the image is then selected, and presumed to be the discharge area. The centroid and major and minor axes of this region are then calculated, with the major axis taken as the estimate of discharge length. An example of this analysis procedure is shown in figure <ref>. Here, the ellipse corresponding to the identified major/minor axes is shown, superimposed on the original greyscale image.



The empirically determined parameters of (<ref>) and a systematic length error (l_err) comprise a vector 

    =[α  β   n    l_err]^T

whose value can then be estimated by non-linear least squares minimisation of the function

    f()=V̂(θ_1,θ_2,θ_3,I,L+θ_4)-V

where V, I and, L are vectors of all measured voltage data, and their corresponding currents and lengths. V̂ is the voltage predicted by the model (<ref>), and is therefore a function of the fit parameters θ_1–θ_4 as well as measured current and length data. Here, the systematic length error (θ_4=l_err) is added to the measured length. The measured data together with the fitted voltage/length relationship is shown in figure <ref> for current values of 100 and 250 mA.

By linearising f about its least squares minimiser , asymptotic estimates of the standard error for each parameter θ_i can be calculated as <cit.>

    _θ_i = f() /√(N-P)√(( J_f()   J_f()^T )_i)

where N and P are the sizes of V and  respectively. The fitted parameter values with asymptotic estimates of the associated standard error are given in table <ref>. The corresponding data and fitted length/voltage lines for two of the current values are shown in figure <ref>. The relatively large errors are a reflection of the difficulties involved in the experiment and the image based evaluation methods, as discussed in <cit.>.







 §.§ Spatial power distribution


The spatial power distribution is accounted for by another empirical relationship, namely

    q_src(r,x) = (v-)/lσ^22π(- r^2/2σ^2)   < t < ,    0 < x < l

This describes the discharge as cylindrically symmetric, for axial coordinate x and radial coordinate r. The distribution of power is thus Gaussian in the radial direction and a uniform in the axial direction. This Gaussian approximation has been used in classical models of electric arcs <cit.>. Although it is known that power density in the axial direction is likely to be non-uniform due to electrode effects, the non-uniformity cannot yet be reliably quantified from the data.



The parameter σ of (<ref>) is derived from line profiles of intensity taken from the images (figure <ref>). The Gaussian function is fitted to this profile data, providing an estimate of σ. This analysis is carried out across the image dataset to determine the relationship between σ and the length and current of the discharge. In the absence of a physical model, the dependence of σ on both quantities is presumed linear, with an expression of the form

    σ=σ_0+a_σ l + b_σ

for length coefficient a_σ and current coefficient b_σ. This expression is fitted to the data using a conventional linear least squares method, and standard errors of the parameters calculated in the usual manner <cit.>. The data and fit are show in figure <ref>, and parameters with corresponding standard errors of regression are shown in table <ref>.







 §.§ Model Validation

The validity of the electrical part of the empirical model can be assessed by examining voltage and current predictions under various circumstances. In section <ref>, it was noted that the electrodes were connected to an electronic source circuit designed to maintain constant current. By replacing this special circuit with a simple resistive or inductive source, a discharge with time varying current is produced, and its electrical waveforms recorded. These circuits, together with a component described by (<ref>), can then be implemented in a common commercial circuit simulator to produce a prediction for comparison to the data. Note that only the parameter u_arc must be adjusted for these simulations, with the others taken from table <ref>. A comparison of predicted and actual discharge waveforms for the resistive and inductive source circuits are shown in figures <ref> and <ref> respectively.





The ability of the model to provide a good quantitative prediction under conditions outside of the original dataset used for parameter fitting is a good indication of its validity for describing the electrical characteristics of the discharge. The validity of the spatial power distribution is a more difficult question, as spatially resolved measurements of power, or even temperature, cannot be easily made. Hence, the proceeding sections attempt to shed some light on this issue by consideration of the underlying physics.



§ PHYSICS


The empirical source term of section <ref> is simulated in a gas dynamics model together with a simplified consideration of heat conduction and evaporation in the solid cadmium electrode via boundary conditions. A gas mixture of 21% hydrogen in air is simulated, this being a standard mixture used to investigate ignitions by the discharges <cit.>.

The geometry and dimensions of the computational domain are described in figure <ref>. This is intended to represent the scenario of figure <ref>, excluding the moving tungsten wire electrode. The relevant conservation equations and boundary conditions are discussed in this section, together with details of the calculation of electrical conductivity based on an equilibrium plasma composition.






 §.§ Conservation equations


The conservation equations of the model are given in (<ref>)–(<ref>). Conservation of mass is stated by,

    ρt + ÷(ρu) = 0

for density ρ and velocity u. Additionally, conservation of momentum is given by the Navier-Stokes equations

    ρut + ÷(ρuu) + p = ÷τ + ρg

for pressure p and gravity vector  g = [0    -9.8    0]^T. The viscous stress tensor τ is calculated from Newton's law

    τ=-μ[u+(u^T)-2/3(÷u)I]

for dynamic viscosity μ and identity matrix I. Conservation of enthalpy (h) is given by

    ρ(h+e_k)t + ÷(ρu(h+e_k) +j_q) - pt = q_src

for heat flux j_q, specific kinetic energy e_k=1/2uu and q_src as per (<ref>). Finally, conservation of the ith chemical species is given by

    ρ Y_it+÷(ρuY_i + j_i)  =  0

for species mass fraction Y_i and diffusive flux j_i. The chemical species considered are _2, _2, _2 (from the gas mixture), and Cd (evaporated from the electrode). No source terms are present, as chemical reactions are not considered.

In the above equations, pressure is related to density by the ideal gas equation p=ρ R T for specific gas constant R and temperature T. Temperature is derived from enthalpy by numerically inverting the equation h=∫_0^Tc_p T, where c_p is the temperature dependent mixture specific heat, obtained from a thermodynamic database <cit.>. The diffusive fluxes are given by

    j_i = -ρ D_i Y_i - D_T_i/TT

and

    j_q = -kT + ∑ h_i j_i

for mixture averaged diffusion coefficients D_i, thermodiffusion coefficients D_T_i and thermal conductivity k. These transport parameters, together with viscosity μ are calculated as a function of temperature using the first Chapman-Enskog approximation <cit.> for the pure species, and then combined using mixing rules <cit.>.



 §.§ Boundary conditions

The open boundaries of the computational domain are simulated using a Dirichlet boundary condition for pressure p=1 atm and applying Neumann boundary conditions to all other variables.

Heat conduction at the electrode surface is approximated by solving the 1-D diffusion equation in the perpendicular direction for each point on the surface

    ρ_s c_p_sTt = qx_s
    
    q=k_sTx_s     x_s > 0

where x_s is the depth.  Here ρ_s, c_p_s, and k_s and the density, specific heat and heat conductivity of cadmium, assumed constant for the simulation. The equation is solved up to a depth of 100 μm, corresponding to a 2 ms characteristic time of heat conduction. As suggested by <cit.>,  heat flux and constant temperature boundary conditions are used for the solid (<ref>) and gas (<ref>) respectively, with the value for each derived from the other's solution.

The heat flux at the surface q_x_s=0 comprises of three components, namely

    q_x_s=0  = q_cond + q_fall - q_vap

The first is the heat conducted from the gas domain.

    q_cond = - kT_gasx_s

The second represents the heating produced in the cathode fall region,  which is assumed to possess the same radially symmetric spatial distribution as the power within the discharge – as described in (<ref>)

    q_fall = v_fall/σ^2√(2π)(- r_s^2/2σ^2)

where r_s represents distance along the surface from the origin of the discharge. It is also assumed here that this energy is completely transferred towards the cathode surface and not to the plasma. The last component  represents heat lost due to evaporation

    q_vap  =	(L_v+c_p_s T_x_s=0)J_m

where L_v, is the specific heat of vaporisation. J_m is the surface normal mass flux  due to evaporation. This is calculated with

    J_m = -ρ DCdYCdx_s + α_ep_v/√(2 π RCdT_x_s=0)

where the first term represents diffusion flux and the second bulk evaporation according to the Langmuir relation <cit.>. Here, DCd, YCd and pCd represent the diffusion coefficient, mass fraction and partial pressure of cadmium at the surface. RCd is the specific gas constant for cadmium. An accommodation coefficient α_e is calculated according to the procedure of <cit.>. The equilibrium vapour pressure p_v is estimated as a function of temperature from the August equation coefficients given in <cit.>.

To account for the evaporation process, Dirichlet boundary conditions on the surface are also imposed for velocity in the surface normal direction (|u| = J_m/ρ) and cadmium mass fraction (YCd = p_v MCd/(pM)), where MCd and M are the molar mass of cadmium and and the mixture average molar mass respectively.



 §.§ Plasma


For the calculation of electrical conductivity, a pure cadmium plasma at a pressure equal to the partial pressure of cadmium in the gas mixture has been assumed. This assumption is based on the previous results showing the discharge radiation to comprise solely of cadmium lines, as mentioned in section <ref>. Determination of the transport properties of this cadmium plasma has been performed using the procedure presented in <cit.>. First, the plasma composition has to be calculated. Under the assumption of local thermodynamic equilibrium (LTE), the temperature dependent densities of plasma species, namely, electrons, atoms and ions, can be predicted using the method of the minimisation of generalised Gibbs free energy <cit.>. The following components have been taken into account: electrons e, atoms Cd, ions Cd^+, Cd^2+, Cd^3+, Cd^4+ and Cd^5+. The densities have been calculated for various vapour pressures. When the densities are known, transport parameters for the corresponding pressures and temperatures can be calculated.

Similarly to the calculation of gas transport parameters (section <ref>), the Chapman-Enskog theory is utilised for obtaining transport properties of the thermal plasma in equilibrium state <cit.>. In the specific case of a monatomic gas, transport parameters can be obtained using the approach developed by Devoto <cit.>. The starting point for the determination of transport coefficient is the calculation of kinetic integrals which requires the knowledge of the scattering cross sections as input data. Kinetic integrals for Cd-Cd collisions have been calculated using the Lennard–Jones model potential with parameters from <cit.>. For the determination of the data for the collisions between electrons and Cd atoms scattering cross sections from <cit.> have been adopted. In order to describe the ion-atom interaction the common assumption that the interaction cross section for the ion-atom pair is generally defined by polarisation ion capture by the atom has been used.
The polarisability of cadmium  has been taken from <cit.>. Finally, interaction between charged particles has been characterised by the Coulomb cross-section <cit.>. Figure <ref> shows the resulting temperature dependence of electrical conductivity for different vapour pressures.





§ SIMULATION RESULTS


The equations of section <ref> are implemented in the OpenFOAM finite volume method solver <cit.>. The pressure implicit with splitting of operator (PISO) method is used to solve for pressure and velocity <cit.>. Transport parameters in the conservation equations are calculated using the Cantera software <cit.>.

Simulations are conducted for the discharge parameter values given in table <ref>, these being representative of the data collected in <cit.> and of the application in general. Here, the simulated discharges are specified in terms of current () and duration (), from which velocity  is calculated using (<ref>). For each simulation, the time varying fields of mass fraction for cadmium are used to calculate the corresponding partial pressure. The partial pressure and temperature fields are then used to calculate an electrical conductivity field using the data of section <ref>. Example fields for temperature, cadmium partial pressure and electrical conductivity are shown in figure <ref>.





A notable characteristic of figure <ref> is the elongated distribution of temperature and cadmium partial pressure. Note that the discharge length at this point is around 170 m, however the high temperature region extends well beyond this length. This is due to a jet of evaporated material produced by the cathode fall region heating, which can be seen in the velocity distribution of figure <ref>.



The electrical conductivity field is then sampled to produce an averaged distribution in 2-D cylindrical coordinates. An example of the axial and average radial profile of electrical conductivity is shown in figure <ref>.



The sampled axisymmetric distribution of electrical conductivity can be numerically integrated to estimate overall resistance Z.

    Z = ∫_L_0^L_1dl/∫_0^R_1 2π r σ_e dr

The radial limit of integration R_1 can simply be made large enough to include the conductive region (200 m is selected here). The axial limits L_1 and L_0 require consideration. As the axial profile of figure <ref> shows, conductivity is negligible near the ends of the discharge, particularly on the cathode (left) side. Since electrode layer effects are not modelled, including the areas near the electrodes in the integration would not be meaningful. The analysis is thus conducted by integrating only over the region where conductance is no less than 20% of its maximum value. Under these conditions, the evaluation of (<ref>) will provide a lower bound on possible resistance.


The resistance at each time step of the simulated cases was calculated both (<ref>) and the empirical relationships of (<ref>), as Z=v(i,l)/. A comparison of the results from both calculation methods can be seen in figures <ref> and <ref>, for the simulated discharges of 200 and 1000 s respectively. Here, the expected values are those predicted by (<ref>) as described above. The comparison provides an indication of consistency of the model under the LTE assumption. The results show an inconsistency of at least an order of magnitude, which is particularly stark at the beginning of the discharge period.







§ DISCUSSION


The results show that an LTE temperature–electrical conductivity relationship  cannot adequately describe the discharge in question. Even under the favourable conditions assumed here, (i.e.: no radiation losses, near electrode regions assumed perfectly conducting), the discrepancy between electrical resistance predicted this relationship, and that predicted by an empirical correlation is not less than one order of magnitude. In addition to the difference in magnitude of the electrical resistance, the way in which this resistance varies over the duration of the discharge is also noteworthy. Under the LTE assumption, the resistance would be expected to decrease over time. This is because the average temperature increases over time, and the electrical resistance decreases near-exponentially with temperature (and increases only linearly with length). The data and empirical correlations of section <ref>, however, show that voltage is linearly proportional to length. This implies a constant electrical resistance per unit length over the duration of the discharge.

From these observations, it can be reasonably concluded that the properties of the plasma are not (or at best, very weakly) dependent on the temperature of the gas. This is a defining characteristic of a non-equilibrium (or non-thermal) plasma. Although further investigation is required to determine the exact nature of the plasma, it could possibly be described as a type of glow discharge. An obvious question arising from this hypothesis is how the empirical model of section <ref> is able to describe the electrical characteristics of the discharge with reasonable accuracy. This is surprising, given that the hyperbolic relationship between voltage and current is based on LTE, specifically, the idea that higher currents produce a broader and hotter (therefore more conductive) discharge column. It is possible that the apparent applicability of the model is coincidental. One theory of stable glow discharges suggests that a constant current density at the cathode is maintained, with increases in current only increasing the size of the electron emitting area (or cathode spot) <cit.>. It is conceivable that this would, in turn, also increase the width of the discharge column, leading to a qualitatively similar voltage/current relationship as in the case of a stable arc. It should be stressed that these ideas are speculative, and require confirmation through experiment and/or more detailed simulations.



§ CONCLUSION


The results have significant implications for the simulation of combustion initiated by the electrical discharges in question. Thus far, the empirical modelling approach has shown some success in predicting electrical properties of the discharge. In the prediction of flammable gas ignition, however, the gas temperature, or more accurately, the heavy particle translational temperature, is also an important parameter. The current empirical approach detailed in sections <ref> and <ref> assumes that all of the electrically measured power is thermalised, and distributed spatially in a similar manner to the discharge's observed radiation intensity. This is problematic, firstly because the power loss due to this radiation is disregarded, and secondly, because the distribution is difficult to accurately measure by experiment. Although previous research <cit.> found that qualitative features of the flammable gas ignition were adequately predicted by this empirical approach, a more physically realistic approach would be more satisfactory, and could potentially approach quantitative accuracy without arbitrary model tuning. Conventional simulation methods which use a temperature–electrical conductivity curve to couple hydrodynamics to electromagnetic equations cannot be used to solve this problem, as they rely on the LTE assumption.

An improvement on the experimental side would require more sophisticated measurement apparatus, such as absolutely calibrated spectroscopy or spectrally resolved imaging. On the numerical side, improving the physical correctness of the model would only be possible by accounting for the non-equilibrium thermodynamics and associated radiative processes. A more rigorous treatment of the electrodes, especially electron emission, is also required. These extensions to the model pose significant challenges due to the uncommon material (cadmium) involved, for which limited published data exists.



§ ACKNOWLEDGEMENTS

The authors wish to thank Carsten Uber and Steffen Pohl for the provision of the experimental dataset, Mohammed Abdul Moiz for the creation of a program to automate the image based measurements, and Thomas Uehlken for the development of the experimental apparatus.




§ REFERENCES


iopart-num

