1. J. Wendelstorf: Two-temperature, two-dimensional modeling of cathode-plasma interaction
   in electric arcs. In P. Pisarczyk, T. Pisarczyk and J. Wolowski (Ed.):
   Proceedings of the XXIV International Conference on Phenomena in Ionized Gases, Warschau, Polen, 11-16 July 1999 (ISBN 83-902319-5-6), Bd.II, S.227-228
   (best computational poster award)
   
2. J. Wendelstorf, H. Wohlfahrt und G. Simon: Two-temperature, two-dimensional modeling of cathode hot spot formation including sheath effects - application to high pressure xenon arc lamps. In: C. Deeney (Ed.): 1999 IEEE Int. Conf. on Plasma Science, Monterey, California, USA, 20-24.6.1999. IEEE Cat. No. 99CH36297 (ISBN 0-7803-5224-6), S.240 (5P03)
   
3. J. Wendelstorf, I. Decker, H. Wohlfahrt und G. Simon:
   Modeling of arc-electrode systems - preliminary results for high pressure xenon arc lamps.
   8th Int. Symposium on the Science & Technology of Light Sources (LS-8), Greifswald,
   30.8.-3.9.1998 (ISBN 3-00-003105-7), S.382-383
   
4. J. Wendelstorf, I. Decker, H. Wohlfahrt und G. Simon: Investigation of cathode spot behavior of atmospheric argon arcs by mathematical modeling.
   In: G. Babucke (ed.), XII Int. Conf. Gas Discharges and their Applications, Greifswald,
   8.-12.9.1997 (ISBN 3-00-001760-7), Bd.I (1997), S.62-65
   
5. I. Decker, O. Voß, J. Wendelstorf und H. Wohlfahrt: Nutzen numerischer Simulationen für den Einsatz realer Schweißprozesse. DVS-Berichte, Bd.186 (1997), S.34-38
   
6. J. Wendelstorf, I. Decker, H. Wohlfahrt und G. Simon: TIG and plasma arc modelling.
   In: H. Cerjak (ed.), Mathematical Modelling of Weld Phenomena 3 (ISBN 1-86125-010-X),
   The Institute of Materials, London, 1997, S.848-897
   
7. B. Rethfeld, J. Wendelstorf, T. Klein und G. Simon: A self-consistent model for the cathode fall region of an electric arc. J. Phys. D: Appl. Phys., Bd.29 (1996), S.121-128
   
8. J. Wendelstorf, I. Decker, H. Wohlfahrt und G. Simon:
   Berechnung anwendungsrelevanter Zielgrößen von thermischen Lichtbögen:
   Kathodenfallmodelle und deren Integration in die Lichtbogenmodellierung.
   8te Bundesdeutsche Fachtagung Plasmatechnologie, Dresden, 14-17.7.1997, S.17
   
9. J. Wendelstorf, G. Simon, I. Decker und H. Wohlfahrt: Untersuchungen zur Berechenbarkeit technologischer Parameter thermischer Lichtbögen am Beispiel des atmosphärischen Argonbogens. Verhandl. DPG (VI), Bd.32 (1997), S.7
   
10. J. Wendelstorf, I. Decker, H. Wohlfahrt und G. Simon:
    Modellierung thermischer Lichtbögen zur Optimierung von technologischen Prozeßparametern.
    7te Bundesdeutsche Fachtagung Plasmatechnologie, Bochum 13-14.3.1996, S.118
    
11. J. Wendelstorf, I. Decker, H. Wohlfahrt und G. Simon:
    Modellierung des Kathodenfalls zur Simulation thermischer Lichtbögen.
    7te Bundesdeutsche Fachtagung Plasmatechnologie, Bochum 13-14.3.1996, S.37
    
12. B. Rethfeld, T. Klein, C. Tix, J. Wendelstorf und G. Simon:
    Modellierung der katodennahen Plasmarandschichten eines WIG-Lichtbogens.
    Verhandl. DPG (VI) Bd.30 (1995), S.243
    

Thermische Plasmen wurden in Form von Lichtbogenentladungen Anfang des neunzehnten Jahrhunderts entdeckt. Die fundamentalen physikalischen Mechanismen dieser Gasentladungen wurden bereits im Jahre 1903 von J. Stark beschrieben. Während die grundlegende Physik dieses Phänomens somit weitgehend bekannt ist, war es bis heute praktisch kaum möglich, wesentliche Entladungsparameter quantitativ vorauszuberechnen.

Ein fundamentales Problem der Modellierung lag in der Beschreibung der Nichtgleichgewichtsgebiete (Plasmaschichten) vor den Elektroden. Das andere war die Tatsache, daß es sich bei diesen elektrischen Entladungen um dissipative, selbstorganisierende Systeme handelt, deren Eigenschaften sich erst aus der komplexen Wechselwirkung ihrer Teilsysteme ergeben (Emergenz). Diese Probleme werden im Rahmen dieser Arbeit diskutiert und weitgehend gelöst.

Basierend auf der Motivation des Themas und der Geschichte seiner Untersuchung wird zunächst ein Konzept zur vollständigen quantitativen Vorausberechnung dieses Typs von Gasentladungen entwickelt. Es folgt eine physikalische Beschreibung der Lichtbogenteilsysteme Elektrodenfestkörper, Elektrodenoberfläche, Raumladungs- und Ionisationsschichten und der Plasmasäule. Die zur Lösung notwendigen Plasmaparameter und Transportkoeffizienten werden für den Fall des partiellen lokalen thermodynamischen Gleichgewichts (pLTG) in Abhängigkeit von Elektronen- und Schwerteilchentemperatur und Entladungsdruck berechnet. Mit diesem ab initio Gesamtmodell werden dann zunächst die wesentlichen Eigenschaften der Entladung im Detail berechnet. Eine Variation der Entladungsparameter Gas, Druck, Strom und Katodendurchmesser weist im folgenden die quantitative Berechenbarkeit der fundamentalen Verhaltensweisen von Lichtbögen im Entladungsdruckbereich von 0.1 bis 8 MPa und Strömen oberhalb von 1A für Argon, Xenon und Quecksilber als Plasmamedium nach. Durch das rechnerische Ausschalten einzelner physikalischer Effekte wird deren konkreter Einfluß auf das Verhalten der Entladung untersucht. Diese Analyse gestattet zusammen mit einer im Bezug auf die benötigten Stoßquerschnitte durchgeführen Sensitivitätsanalyse eine Bewertung des Vergleichs mit anderen Modellen und den wenigen quantitativen Messwerten, die für derartige Hochdruckbogenentladungen bisher publiziert wurden.

Für sehr unterschiedliche Entladungsparameter wird ein Vergleich mit dem Experiment vorgenommen, welcher insbesondere die hohe Genauigkeit des entwickelten Kathodenfallmodells nachweist und erste Validierungsaussagen liefert.

Die Zusammenfassung stellt die wesentlichen Neuerungen des Berechnungskonzeptes und die grundlegenden physikalischen Prozesse kurz dar und beschreibt die große Zahl der möglichen Anwendungsbereiche des Modells. Es werden zudem Detailverbesserungen angeregt und Hinweise zur Durchführung von weiteren Validierungsexperimenten gegeben.

Bei dem vorliegenden Berechnungsverfahren wurde weitgehend auf ungerechtfertigte Vereinfachungen verzichtet und durch Kopplung der modellmäßigen Erfassung der Teilsysteme erstmalig die weitgehend vollständige quantitative Vorausberechenbarkeit dieses Typs von Gasentladungen nachgewiesen.

Thermal plasma gas discharges (electric arcs) were discovered at the beginning nineteenth century. The fundamental physical mechanisms of such gas discharges were described first by J. Stark in 1903. While the basic physical laws governing this phenomenon are mostly well known, up to now, a quantitative prediction of the fundamental discharge parameters was practically impossible.

One fundamental modelling problem was the description of the non equilibrium layers (plasma sheaths) in front of the electrodes. Additionally, such electric discharges are dissipative self organizing systems. Their properties emerge from the complex interaction of their parts. These problem will be addressed and solved by this work.

Based on a motivation of the objective and the history of its investigation, a concept of a complete self consistent quantitative prediction of such electric arcs was be developed. A physical description of the physical partial systems electrode solid, electrode surface, space charge and ionization sheaths and plasma column is provided. For the case of partial local thermodynamic equilibrium (pLTE), the plasma parameters and transport coefficients as a function of electron- and heavy particles temperature are calculated. Using this ab initio model, the properties of arc discharges are computed. A variation of the discharge parameters gas, pressure, current and cathode diameter establishes the quantitative predictability of the fundamental arc behaviour for a discharge pressure range of 0.1 to 8 MPa and arc currents above 1 A for argon, xenon and mercury fillings. By computational disableing individual physical effects, their specific influence on the discharge behaviour will be investigated. Together with the sensitivity analysis performed with respect to the required cross section data, this analysis allows for an assessment of the comparison with other models and the few available experimental results, actually published for such high pressure arc discharges.

For a number of very different discharge parameters an experimental validation is provided. Especially the high precision of the cathode layer model becomes evident leading to a first positive model assessment.

Finally, the main innovations of the computational concept and the basic physical processes together with the broad range of application areas of the model are summarized. Some enhancenments are proposed and the realization of validation experiments will be discussed.

The present computational scheme cuts out most unjustified simplifications and, by a numerical coupling of the partial systems models, a complete quantitave predictability of this type of electric gas discharges is established.

Contents

Chapter 1
Introduction

From the viewpoint of a scientist, atmospheric pressure gas discharges are very easy to obtain and thus understanding and modelling is a matter of scientific interest since the first detection of the arcing phenomenon itself in the nineteenth century. During the last decades, modern computing hardware and advanced software provides the opportunity to compute the overall discharge behaviour in advance and in a quantitative way. This work is dedicated to contribute to the understanding of the overall discharge properties by using these new modelling tools. This objective is reached by the development and application of accurate and validated modelling software.

The investigations are applicable to the so-called thermal plasmas, which are distinguished from other types of plasmas by the following characteristics:

• All species have Maxwellian velocity distribution. If they are all according to the same temperature the plasma is named to be in Local Thermodynamic Equilibrium (LTE), otherwise the equilibrium is stated to be partial (pLTE).
• The excited energy states are populated to a Boltzmann term.
• Composition can be derived from Local Chemical Equilibrium (LCE) or demixing effects have to be calculated from an additional set of diffusion equations.
• Temperatures are less than a few eV (1eV=11604K).
• The plasma is not confined by an external magnetic field.
• The discharge is stabilized by its own magnetic field, by external gas flow, the electrodes or the surrounding solid or liquid walls.

In the American and European literature these plasmas are often called hot plasmas (the heavy particles are hot compared to low pressure gas discharges). In the Russian literature and in Germany they are often called low temperature plasmas, because their electron temperature is small compared to nuclear fusion plasmas.

During the last decades it became increasingly clear that the existence of LTE in a plasma is rather an exception than a rule. The excited states population often deviates from Boltzmann equilibrium. This is of major importance for spectroscopical plasma diagnostics, but for most practical applications at least the partial Local Thermodynamic Equilibrium (pLTE) concept can be used for the following reasons:

• The contribution of the excited states to the total plasma enthalpy is very small.
• The local energy balance of the discharge is dominated by Joule heating, radiative,
  convective and conductive source terms.
• The plasma composition and transport coefficients mainly depends on the excitation (electron) temperature and the fundamental gas properties.
• There is no direct impact of specific excited states populations on global arc discharge properties.
• Large errors in the calculation of the excited states density highly affect experimental methods, were these quantities are required for data interpretation, but in the self consistent system calculation, the effect is limited to very small changes in the calculated self-consistent plasma temperature distribution.

Typical discharge parameters are listed in table 1.1. The values are reflecting the typical application areas of these arc discharges: lighting and material processing, synthesis and destruction.

In the following section, a short historical and application overview is provided. This is followed by the definition of ab initio modelling, the objective of the following chapters, as summarized in section 1.4.

1.1  A short history of thermal plasma gas discharges

High intensity gas discharges are also called electric arcs, because the first discharge observed in detail was a carbon arc burning in horizontal position in the ambient atmosphere and thus bent by natural convection.

In table 1.2, a small collection of the milestones in the 200 years history of the arcing phenomenon is sketched. The arc discharge not only forms a basis for efficient discharge lamps, it is also an initiator for the development of plasma physics itself in the 1920's.

From discovery of a certain phenomenon (e.g. the use of metal halides for lighting purposes) to it's industrial use it can take more than 50 years and up to now, most practically important discharges cannot be computed in advance and thus may not be regarded as fully understood.

Figure 1.1: Development of HID science and technology.

Figure 1.2: Applications of high intensity gas discharges.

Concluding from this long lasting history of the subject, every fifty years, there was a qualitative advance in the overall development paradigma (see figure 1.1). This work is dedicated to the actual change from trial & error supported by a first quantitative understanding to a computation of specific discharge features in advance.

1.2  Fields of model applications

As demonstrated in figure 1.2, discharge powers range from a few W to several MW and discharge pressures are atmospheric in most cases. For lighting and (underwater) materials processing purposes, hyperbaric applications became more and more important in the last decades.

This work is dedicated to the computation and quantitative understanding of arc discharges. As a first step, the models will be applicable to stationary discharges showing cylindrical symmetry:

• Arc lamps, at low frequency alternating current (AC) or direct current (DC), especially for the optimization of
  • electrode geometry and erosion,
  • global energy balances,
  • radiative properties.

• Welding arcs (see [127])
  • optimized arc configurations will permit higher processing quality and performance,
  • the impact of new gas mixtures and electrode materials can be determined in advance.

• Arc furnaces (steel and high melting point materials synthesis)
  • optimized high power torches will allow to replace conventional carbon electrodes.

• Plasma torches for enhanced or even new applications:
  • waste destruction,
  • gas heating for reentry simulation and space propulsion,
  • synthesis of nanostructure materials (fullerenes),
  • plasma chemistry and materials deposition (plasma spraying) or removal (plasma etching).

1.3  A practical definition of ab initio modelling

Figure 1.3: Information flow and subtasks of ab initio modelling.

The common meaning of ab initio modelling is computation based on first principles and fundamental constants. The practical definition used for this work is sketched in figure 1.3:

An ab initio model is based on fundamental physics and any material, gas- or plasma property determined by a discharge configuration independent method. The model should not contain any parameter or input value which has to be fitted to experimental discharge data for every individual configuration under investigation. A practical ab initio model should allow for the computation of important discharge characteristics with an accuracy better than experimental error or sufficient for automatic optimization of discharge configurations.

The underlying motivation for model development with such high accuracy is the extrapolation capability these models should have. Additionally the strong nonlinear character of the arc discharge implies a need for relative high computation accuracy. The model validation by quantitative comparison with experimental data has to be supplemented by a sensitivity analysis: Some experimental quanitities are rather insensitive to changes in the arc configuration (e.g. doubling the current often implies only a few percent change in plasma temperature), while other are easy to determine (e.g. voltage-current characteristics) and also sensitive to important nonlinearity effects like mode transitions.

Sensitivity analysis has to be undertaken within all fields important for the development of the models and their validation:

• Modelling accuracy can be restricted in principle, by the numerical techniques applied and by the available computing resources.
• Validateability of the physical effects included will require detailed knowledge of the transport coefficients involved or the need to compute critical effects like mode transitions. The validation of physical details may be completely hidden by configuration- or measurement-inaccuracies.
• Analysis of the numerical schemes and their parameters allows to increase computation speed as well as to obtain important information on the principle accuracy of the concept. As a result, some of the involved numerical schemes may be less accurate (10^-3), while other have to solve stiff problems at high accuracy ( < 10^-8).

One of the most important experience from modelling such systems can be summarized in one sentence:

The validity of the physical description of the system can be justified a posteriori (afterwards), but not a priori (in advance).

Another is the importance of a sensitivity analysis in order to obtain information about the role of a specific input parameter (e.g. electrode work function) on the overall or local discharge behaviour.

The modelling approach is therefore relying on the comparison with selected reproduceable experimental data of known accuracy. The answers of this validation and sensitivity analysis are not always positive:

• The specific experimental data may be not accurate enough to validate a specific modelling approach, e.g. a validation of a cathode layer model by plasma temperature measurements may require unaffordable spatial resolutions and accuracies - cathode layer models are better justified by an electrical determination of the cathode fall voltage.
• The physical detail may be hidden by a material property not known with sufficient accuracy (e.g.  cathode material work function).
• The physical or application detail may be superimposed by the dependency of the results on numerical parameters (e.g. grid size).
• The computing resources may not allow for a sensitivity analysis at all.
• The behaviour of the discharge system itself can be out of the range of the specific modelling approach (e.g. time dependent and 3-D, chaotic or stochastic).

The model will enhance the understanding of the arc phenomenon itself, but the overall objective is quantitative prediction of specific discharge configurations. Some sort of understanding of such complex systems can be obtained by many different plausible explanations or models with internal fit parameters. But these will fail if new discharge configurations are to be computed in advance.

1.4  Outline of this work

Mathematical or physical modelling conventionally consists of a derivation of the equations describing the system and then solving them. This requires making all assumptions and approximations in advance. As the systems under investigation become more complex, this procedure becomes inapplicable in terms of effort and due to the nonlinearity of nature. Initially made assumptions need to be controlled afterwards and the modelling work becomes an iterative process itself.

Most of all important arc phenomena occur as an emerging result of the interaction of different physical processes and discharge regions. The accurate modelling of the individual parts of the discharge is a prerequisite, but the iterative linkage of these models is of the same importance.

As sketched in table 1.3, the overall modelling approach is discussed first in chapter 2. After sectioning the discharge into regions governed by different physical processes, the global iteration algorithms needed for the complete discharge description are presented. Within chapter 3, the physical models of the plasma column, the electrodes and the boundary layers are presented. The computation of the basic plasma properties and transport coefficients is summarized in chapter 4.

In chapter 5, the model proves its ability to reproduce all major discharge features. Initially, an easy computable arc configuration is selected, and the detailed results available from the modelling data are provided. Additionally, the impact and origin of the flow pattern formation are discussed. For this model lamp, the general arc behaviour is computed for a large range of external parameters. The rest of the chapter is dedicated to the validation of the model by comparison with available experimental data.

The last chapter (6) gives a summary of the achievements, it applications and the possible future enhancements of the model.

Chapter 2
Conceptual framework for modelling electric arc discharges

The objective of this chapter is to develop the concept for an ab initio model of the overall thermal plasma gas discharge behaviour. As sketched in figure 2.1, there are at least three physically different regions:

• The (solid or liquid) electrodes and their surfaces.
• The thermal plasma column and its interaction with the environment and the layers.
• The nonequilibrium boundary layers providing the transition from a solid surface at 1¼4 kK to an (partial) equilibrium plasma at 4¼100 kK.

Looking at the details of these regions, one realizes even finer structures, e.g. the boundary layer split into a space charge sheath and ionisation nonequilibrium presheath.

As a conclusion, the overall electric discharge system can not be described by a unified mathematical model, or such a model will be useless and untreatable. One has to realize a hierarchy of physical processes occuring in different regions, at different scales and with different relevance to specific features of the discharge. The electric arc is an arrangement of parts, so intricate as to be hard to understand or deal with - a complex system. Its behaviour emerges from a self organization of its parts.

For this reason, we have to look at the general behaviour and modelling tasks for general complex systems first, as provided in the next section. Within this framework, in section 2.2 the arc discharge is identified to show the behaviour of a Complex Adaptive System.

The specific modelling approaches for the individual arc regions are discussed in section 2.3. Here, we also provide the algorithms required to link physically and numerically different models without loosing their individual modelling accuracy. This and alternative approaches are discussed in section 2.3, followed by a summary of the whole chapter.

The overall target is not to develop a model of all possible discharge situations, but to develop a framework which may be regarded as sufficient for some simple arc configurations and open for further improvements. We treat the arc discharge as a dissipative complex system and we will get a comprehensive and simple picture of some discharges together with the identification of new complex phemomena relevant to others ¼ (just another iteration towards the impossible complete model).

2.1  Introduction to complex adaptive systems

This section is an adaption of the arc discharge relevant material of an introduction provided by Badii and Politi [9]:

Research on complex systems means searching for general patterns in a number of very distinctive systems. Characterizing complexity in a quantitative way is also a vast and rapidly developing field. In general, complexity is not just another phenomenon, it may be regarded as another approach, contradictionary to the reductionist approach currently used in most areas of science (figure 2.2):

Whenever substantial disagreement is found between theory and experiment, the system has been observed with an increased resolution in the search for its elementary constituents. Matter has been split into molecules, atoms, nucleons, quarks, thus reducing reality to the assembly of a huge number of bricks, mediated by only three fundamental forces: nuclear, electro-weak and gravitational interactions.

The discovery that everything can be traced back to such a small number of different types of particles and dynamical laws is certainly gratifying. Can one thereby say, however, that one understands the origin of the arc discharge? Well, in principle, yes. One has just to fix the appropriate initial conditions for each of the elementary particles and insert them into the dynamical equations to determine the solution. Without the need of giving realistic numbers, this undertaking evidently appears utterly vain, at least because of the immense size of the problem. An even more fundamental objection to this attitude is that a real understanding implies the achievement of a synthesis from the observed data, with the elimination of information about variables that are irrelevant for the sufficient description of the phenomenon. For example, the equilibrium state of a gas is accurately specified by the values of only three macroscopic observables (pressure, volume and temperature), linked by a closed equation. The gas is viewed as a collection of essentially independent subregions, where the internal degrees of freedom can be safely neglected. The change of descriptive level, from the microscopic to the macroscopic, allows recognition of the inherent simplicity of this system.

These examples introduce two fundamental problems concerning physical modelling: the practical feasibility of predictions, given the dynamical rules, and the relevance of a minute compilation of the system’s features.

In fact, as the study of nonlinear systems has revealed, arbitrarily small uncertainties about the initial conditions are exponentially amplified in time in the presence of deterministic chaos (as in the case of a fluid). This phenomenon may already occur in a system specified by three variables only.

The resulting limitation on the power of predictions is not to be attributed to the inability of the observer but arises from an intrinsic property of the system. In section 2.2, as an example, cathode phenomena in low pressure arc dicharges are identified to show such a behaviour.

Nature provides plenty of patterns in which coherent macroscopic structures develop at various scales and do not exhibit elementary interconnections. They immediately suggest seeking a compact description of the spatio-temporal dynamics based on the relationships among macroscopic elements rather than lingering on their inner structure. In a word, it is useful and possible to condense information.

Similar structures evidently arise in different contexts, which indicates that universal rules are possibly hidden behind the evolution of the diverse systems that one tries to comprehend. Many systems can be characterized by a hierarchy of structures over a range of scales. The most strong evidence of this phenomenon comes from the ubiquity of fractals (Mandelbrot, [86]), objects exhibiting nearly scale-invariant geometrical features which may be nowhere differentiable.

Hierarchical structures appear to be a general characteristic of nature. The difficulty of obtaining a concise description may arise from fuzziness of the subsystems, which prevents a univocal separation of scales, or from substantial differences in the interactions at different levels of modelling.

The concept of complexity is closely related to that of understanding, in so far as the latter is based upon the accuracy of model descriptions of the system obtained using a condensed information about it. Hence, a theory of complexity could be viewed as a theory of modelling, encompassing various reduction schemes (elimination or aggregation of variables, separation of weak from strong couplings, averaging over subsystems), evaluating their efficiency and, possibly, suggesting novel representations of natural phenomena. It must provide, at the same time, a definition of complexity and a set of tools for analysing it: that is, a system is not complex by some abstract criteria but because it is intrinsically hard to model, no matter which mathematical means are used. When defining complexity, three fundamental points ought to be considered:

1. Understanding implies the presence of a subject having the task of describing the object, usually by means of model predictions. Hence, complexity is a function of both the subject and the object.
2. The object, or a suitable representation of it, must be conveniently divided into parts which, in turn, may be further split into subelements, thus yielding a hierarchy. Notice that the hierarchy need not be manifest in the object but may arise in the construction of a model. Hence, the presence of an actual hierarchical structure is not an infallible indicator of complexity.
3. Having individualized a hierarchical encoding of the object, the subject is faced with the problem of studying the interactions among the subsystems and of incorporating them into a model. Consideration of the interactions at different levels of resolution brings in the concept of scaling. Does the increased resolution eventually lead to a stable picture of the interactions or do they escape any recognizable plan? And if so, can a different model reveal a simpler underlying scheme?

2.1.1  Complexity measures

After observing a new phenomenon, research initially attempts to find a quantitative measure for it. For complexity this is historically strongly related to discrete mathematics, coding and compression theory. Details may be found in the literature [9,110]. For our discharge modelling context, a specific meaning of the complexity measures is proposed as follows:

• Algorithmic complexities:
  • Algorithmic Information Content, Lempel-Ziv complexity or randomness:
    The amount of modelling software (typically 10^4¼6 lines of code) required to get the results on one hand, and the compressabilty (the modelling data compressability by standard methods is larger for more simple discharge configurations) of the results on the other hand.
  • Logical depth: The computing time (typically more than 10¹¹ floating point operations).
  • Sophistication: The modelling and numerical concepts and the management of their realization.
  • Grammatical, regular language and set complexities:
    The programming concepts and languages used.

• Hierarchical scaling complexities:
  The physical submodels and their interaction.

In the long term, the historical development of discharge models is clearly attracted towards descriptions with balanced complexity. In the short term, the development follows the path of least resistance, where increasing logical depth is simple and increasing sophistication is most difficult.

Before affordable high performance computers emerge in the 1980'ties, the development was focused on the mathematical description of specific physical processes relevant for specific arc configurations or regions. Obtaining numerical solutions was possible for highly simplified models only. Nowadays, special care has to be taken, when using measured quantities (like plasma conductivity) obtained using oversimplified discharge models.

2.1.2  Emergence: self generated complexity

Figure 2.2: The emerging system properties replace the inability of the observer to understand the system in the classical (reductionism) sense.

We speak of self-generated complexity whenever the (infinite) iteration of a few basic rules causes the emergence of structures having features which are not shared by the rules themselves. Simple examples of this scenario are represented by various types of symmetry breaking (superconductors, heat convection) and long-ranged correlations (phase transitions, cellular automata). The relevance of these phenomena and their universal properties discovered by statistical mechanical methods indicate self-generation as the most promising and meaningful paradigm for the study of complexity. It must be remarked, however, that the concept of self-generation is both too ample and too loosely defined to be transformed into a practical study or classification tool.

In any case, two extreme situations may occur: the model consists either of a large automaton, specified by many parameters, which achieves the goal in a short time, or of a small one which must operate for a long time. Consequently, the complexity of the pattern is often visualized either through the size of the former automaton or through the computing time needed by the latter.

Another aspect of the problem emerges, namely, that complexity is associated with the disagreement between model and object rather than with the size of the former. This, in turn, calls for the role of a subject (the observer) in determining the complexity itself through its ability to infer an appropriate model. Of course, a system look: complex as long as no accurate description has been found for it.

These considerations and the limited domain of applicability of all existing complexity measures strongly suggest that there cannot be a unique indicator of complexity, in the same way as entropy characterizes disorder, but that one needs a set of tools from various disciplines (e.g. plasma physics and computer science) and - of course - a significant amount of manpower and computer capacity. As a result, complexity is seen through an open-ended sequence of models and may be expressed by numbers or, possibly, by functions. Indeed, it would be contradictory if the complexity function, which must be able to appraise so many diverse objects, were not itself complex!

As sketched in figure 2.2, the global structure of a system is not unidirectional determined by the local interaction of its parts. The global structure may change local interactions. The objective of this work is not to change the viewpoint from reductionism (escape into detail) to holism (ignore the detail), but to apply both approaches as they are two sides of the same medal.

2.1.3  Is the arc discharge a complex system?

The electric arc discharge is driven by an external supply of electrical energy and permanently exporting heat and radiation to its environment. The plasma electrode system is far from thermodynamic equilibrium and exporting entropy. The discharge process is irreversible and can become chaotic. The distribution of the hot spots across the electrode surface or the discharge itself can show a fractal nature [45]. As shown in the next chapter, the effecive number of degrees of freedom can be reduced by exploiting its hierarchic structure. Such far-from-equilibrium systems manifest self-organization (synergetics). There are mathematical concepts for analysing time series data or nonequilibrium phase transitions by analytical methods [60]. In this work, there will be no attempt to extend these concepts to an only numerically treatable hierarchical system like thermal plasma gas discharges. An analytical treatment of the emerging discharge properties and thus the application of existing synergetic concepts is assumed to be possible by changing the type of the model from ab initio to heuristic.

A complex adaptive or self organizing system is not a mathematically well defined object. The arc discharge simply gets this characterization by another complex adaptive system - the human researcher investigating it. By accepting the complex details and their emerging properties, an increased understanding will finally allow the reader to characterize the arc as a less complex system.

Figure 2.3: The different levels of complexity emerging in electric arc discharges.

2.2  The overall behaviour of arc discharges

Now, we are prepared to realize the different levels of complexity appearing within the parameter range of high intensity gas discharges. This will enable us to define the arc systems predictable by actual modelling concepts (section 2.3) and possible extensions for the future.

In figure 2.3, the different levels of complexity emerging in high intensity gas discharges are sketched. The physical details and models for them are discussed in chapter 3, here we discuss the three major spatial and temporal appearences of the overall discharge:

1. Stationary discharge in cylindrical symmetry:
   This kind of discharge appearance can be found in Tungsten Inert Gas welding arcs [95,65,118] and high pressure lamps [63]. Such configurations are under investigation within current modelling attempts (see chapter 5). The cathode material should be refractory (i.e. not evaporating) and with homogenous emission properties. The electrode geometry must be properly selected (e.g. tipped and/or rounded). Stabilization to cylindrical geometry is provided by magnetic compression, electrode geometry or outer walls. There will be a single hot spot at the electrodes, often called diffuse. Variation of total arc current may lead to a sudden change of the total spot area, a mode transition. The smaller hot spot area mode is called spot mode.

2. [Quasi]Stationary non symmetric discharge:
   Using geometrically problematic electrode configurations (e.g. an unrounded stick electrode) will force the discharge to break cylindrical symmetry. Lower currents support for that by decreasing magnetic compression. Higher currents or hollow electrodes may lead to rotation of the hot spot, as found in arc gas heaters [113].

3. Transient non symmetric discharge:
   Decreasing discharge pressure, the arc starves for ionizable material. Decreasing current provides less heating and thus less electron emission in the normal hot spot mode. The system reacts by forming a variety of smaller and smaller cathode spots on different time and length scales. Additional support for such microspot formation comes from spatial of temporal inhomogenities of the electrode surface, or such inhomogenities result from burning the discharge for a significant time. Over the years, increased experimental resolution support the assumption of some fractal regimes of cathode spot formation (both in time and space).
   Finally we arrive at the vacuum arc with its characteristic transient and non-symmetrical behaviour [71,72,103,69].

• The arc appearence (1-3) depends on all major parameters stated in figure 2.3.
• There are some plausible rules (see above), but prediction of arc appearance for a specific configuration will require models suitable for all possible symptoms.
• It will be a positive quality factor, if e.g. a 2-D stationary model will fail to deliver results for a configuration showing only 3-D instationary arc appearances.
• Deviations of the experimentally realized configuration from the anticipated idealized state may totally hide any other parameter dependencies one tries to observe (e.g. the cathode surface exhibits an inhomogenous work function distribution with dramatic impact on the discharge appearance).
• Discharges may be similar, but not because of the aggreement of some external parameters, but because of similar physics and symmetry.

As a rule of thumb, the arc reacts on harder burning conditions with increased spatial and temporal complexity. The prediction of optimum burning conditions allowing for stationary cylindrical symmetric discharges is one of the objectives of this work.

2.3  Numerical modelling concepts

A physical self consistent decription of the overall discharges (as described in chapter 3) needs to be based on a numerical solution concept, which is realistic in terms of development and computing time (algorithmic complexity and logical depth).

Modern methods of scientific computing have to be evaluated with respect to availability, costs and numerical effort. For electric arc modelling, the plasma is described as a conducting fluid. Mainly the concepts of Computational Fluid Dynamics (CFD) have to be applied. Because of its tremendous complexity and the amount of practical experiences needed for an understanding of such concepts, the following section can provide a summary only. Understanding will require at least the knowlegde of the basic concepts as described in standard textbooks [98,46,5,133].

Figure 2.4: Finite volume (this work, Braunschweig: BS) and finite element (Lichttechnisches Institut, Karlsruhe: LTI) grids used for complete arc discharge modelling (only one half of the FEM grid is shown).

Regarding the numerical solution of the mathematical plasma electrode model, two different approaches are currently in use:

1. The discharge geometry is divided into finite volumes using a non equidistant grid spacing and volume boundaries perpendicular to the spatial coordinates {r,z} (BS, left part of figure 2.4).
2. The computational region is divided into a solution adapted finite element grid based on triangles. (LTI, right part of figure 2.4) [128,50].

In their current implementation (BS: this work, LTI: [128,50]), both concepts have their specific drawbacks.

The orthogonal structured grid of the BS concept wastes a lot of grid volumes because the refinement depends only on the r and z coordinate. A solution adaptive refinement is implemented for the finite volume method by commercial CFD-codes (unstructured grids). The present approach was selected in order to limit code development time to a few years and because of the straightforward fluid flow implematation as described by Patankar [98]. Similar numerical concepts were implemented at CNRS/ENSCP in France [37], CSIRO in Australia [132,82], at UMN/ME in the USA [64,67] and several other groups.

The finite element arc modelling approach (LTI) is described in detail by Wiesman [128]. There is no inclusion of fluid flow phenomena in this implementation. The major drawback is not a numerical one, it is a result of the physical model used by the LTI approach. The neglection of one dimensional sheath effects and the description of the non equilibrium electrode layers by a 2D model implies the requirement for a very fine numerical mesh near the electrodes (see figure 2.4). The mesh elements scale down to a nanometer scale leading to a very large number of finite elements.

Both implementations (LTI and BS) may be regarded as less professional than commercial CFD/FEM software packages, but they allow for the prediction of electric arc behaviour, while commercial packages actually do not support this application.

In the future, easy to use and efficient arc modelling should be based on such commercial tools. The problem of these large codes is the unavailablity of the sources and their algorithmic complexity (more than 10⁵¼10⁶ lines of code compared to the several 10⁴ lines developed at LTI and BS). For the development of the basic technology, the test and evaluation of new physical concepts and in order to have development times of only several man years, only concepts like those realized by this work and the groups cited above, are realistic. The following section will provide a summary of the general numerical concept for electrode plasma linkage (transfer function) developed within this work.

The orthogonal FVM discretization used for this work was selected for a maximum efficiency in code development time. This work is focussing on the investigation of physically self consistent and accurate cathode and anode layer models within the framework of a self consistent calculation of the overall arc discharge including convection effects.

2.4  The transfer function concept

The basic algorithm for solving coupled fluid flow and heat transport problems by the finite volume method is perfectly described by Patankar [98]. The early implementations of this approach for an overall arc electrode model by a so-called conjugate heat transport method can be found in the literature [82,132]. The solution of the 2D plasma equations (electric potential, flow- and temperature-fields) as well as the heat (and current) transport inside the electrodes can directly follow the Patankar approach. His method solves a system of equations of the type

¶ ¶t (rF)+ ® Ñ    · æ è r ® v   F ö ø = ® Ñ    ·GF ® Ñ    F+ SF

(2.1)

to which the set of arc plasma equations can be reorganized, e.g. (see section 3.5):

Figure 2.5: Integration of local 1-D sublayer models by the effective conductivity approach.

Figure 2.6: Electrical conductivity within the thermal plasma cathode system.

Figure 2.7: Transfer function algorithm (see text).

The major problem of the numerical concept is the development of a physically realistic and treatable boundary layer algorithm. In front of the electrodes, sheath effects occur (see chapter 3) and these are very difficult to implement into a 2D model. Because the size of the nonequilibrium layer is only about 10 mm (or even a few nanometers if the space charge layer is regarded as the boundary layer only), it can be treated as a 1D-layer forming a skin in front of the electrodes.

The finite volumes in front of the electrodes (resulting from the 2D meshing) require special attention. A numerically stable approach is to use a 1D layer model for the calculation of the effective electric and heat conductivity (figure 2.5). The result is a non smooth electric conductivity at the plasma electrode interface (figure 2.6). Negative anode sheath voltages (as resulting from the model described in the next chapter), will give even negative electric conductivities and can imply some numerical problems.

The effective conductivity approach works fine for some high current atmospheric argon arcs (see section 5.3.6), but it is not very realistic. Therefore the problem was generalized and found to be a special case of a general iterative scheme for the calculation of related boundary conditions for multidimensional regions linked by surfaces with additional physical processes to be modelled. This so called transfer function approach (figure 2.7) is very general and may also be used for other application areas or numerical schemes.

It is straightforward to get solutions for partial differential equations defined on multidimensional regions with well defined boundary conditions. Thus the electrode plasma system has to be divided into the 2D arc plasma and 2D electrode solid regions. The boundary conditions used for the 2D solutions are heat flux and voltage or current density. The layer between these regions is regarded as a black box transferring actual plasma and solid surface parameters to new boundary condition values for the next iteration step. As a result, the multidimensional regions are linked by a general scheme allowing for an implementation independent of a specific physical model of the processes inside the layer. It is a local concept allowing for a variation of the layer parameters and boundary conditions along the electrode surfaces. Regarding the plasma and solid heat balances, it simply transfers current values of the local plasma and surface temperatures to new heat flux boundary conditions.

2.5  Summary of the concept

As a matter of fact, electric arc discharges are not homogenous physical systems allowing for a description by a single set of equations. The internal boundaries between regions of different physical provenience and the long range interaction of the discharge parts form a self-organizing complex system. In the next chapter, the physical processes governing the individual regions will be described by mathematical modelling. The numerical linkage of the regions is sketched in this chapter. Finally the software implementation of this concept will deliver a computer model allowing to predict most of the characteristic features of the electric arc discharge. Through the emergence of more and more complex structures resulting from specific operating conditions and depending on discharge pressure, electrode geometry and discharge medium, the real arc may be still regarded as unpredictable for a randomly choosen configuration. But the remainder of this work will show the predicablitiy of a large class of discharge parameters: the stationary, radial-symmetric electric arc.

Chapter 3
Physics of the electric arc discharge

3.1  Introduction

The overall aim of this chapter is to develop a numerically treatable ab initio model for a (large) class of thermal plasma gas discharges. Such electric discharges are also called electric arcs. Within the variety of gas discharges, the electric arc is defined by the following typical behaviour:

• The discharge is self-sustaining with a relatively low voltage caused by a cathode fall voltage of only several volts and a high conductivity of the thermal plasma in the column.
• Special discharge ignition techniques like high voltage pulses, short circuit or the increase of pressure at constant current are generally needed.
• The cathode is at thermionic emission temperatures provided by the energy flow from the near cathode plasma or external heating.
• The current density in the cathode hot spot is 100-1000 times that of the corresponding glow discharge.
• The total discharge current is always above the corresponding glow discharge current
  (figure 3.2). There is a non-smooth transition from glow to arc with increasing current.

An introductionary discussion of the different types of arc discharges can be found in the literature [102,47]. Not all of them are treatable by the model described in this chapter or even any ab initio model to be developed in the near future. The model does not attempt to describe non-stationary multi-cathode-spot discharges (like vacuum arcs). Nevertheless, it can deliver the physical basis for future attempts to predict alternating current (AC) arcs used for lighting or welding.

3.1.1  Basic discharge features

As sketched in figure 3.1, a typical arc discharge consists of a (highly luminous) plasma column linked to the electrodes by small (skin like) non equilibrium layers. The basic arc behaviour observed experimentally and discussed in the following sections will be predictable by the quantitative ab initio model developed in this chapter.

3.1.1.1  Arc stabilization

In order to obtain numerical solutions for the model developed within this chapter, the arc discharge has to be stationary and cylindrical symmetric. Most practical applications require such stable burning conditions which are realized by the following means:

• The arc column is confined, i.e. by a quarz or alumina tube with a radius not larger than the natural arc radius, e.g. most HID-lamps, where the electrode separation is much larger than the electrode diameter.
• The arc is stabilized by the electrodes, i.e.  the high current densities of the cathode and anode hot spots generate high magnetic compression forces. The resulting expansion jets dominate the arc flow field and temperature map stabilizing the overall discharge against external pertubations (e.g. in high pressure short arc lamps).
• The arc is stabilized by natural convection, e.g. the horizontally burning arc in air.

3.1.1.2  Current voltage characteristics

As shown in figure 3.2, the arc discharge is distinguished from the glow discharge by much lower voltages at higher current. The transition from glow to arc is in general non-smooth.

The current voltage characteristics of the arc is generally falling. Above a certain current (about 10 times the glow to arc transition current), the variation of the cathode fall is rather small. The anode fall voltage is nearly constant with current. The following model will allow for the calculation of all these voltage variations.

3.1.1.3  Arc column properties

As shown in figure 3.3, the temperature profile in the column of (long) arc discharges is almost parabolic. Deviations of electron and heavy particles temperatures can be located (at least) in the arc fringes.

The electric field in the column of arc discharges is much lower than in glow discharges. This is due to the different ionization mechanisms. In the kinetic regime of the glow discharge plasma, the electric field must deliver electron energy up to the region of the ionization energy (E_ion), while in the equilibrium arc discharge, the ionisation is due to the high energy tail of the electron energy distribution function, thus the electric field must only deliver electron energy of k_bT_e << E_ion.

As shown in figure 3.4, the electric field in the center of the arc discharge increases with total discharge pressure (the mean free path to gain energy from the electric field decreases with increasing pressure) and the total current (heat loss by conduction, radiation and convection).

3.1.1.4  Cathode hot spot and fall voltage

Glow discharges need a high cathode voltage drop, because the electrode is not at thermionic emission temperatures. The ions are accelerated in the cathode fall in order to provide an electron current sustained by secondary emission from the cathode surface.

The cathode voltage drop of the arc discharge is much smaller and its task is to heat the cathode surface to thermionic emission temperatures and to accellerate the emitted electrons for gaining energy needed for ionization within the near cathode plasma layer. Below a certain saturation current, the cathode fall voltage strongly depends on total arc current and decreases with increasing pressure.

As an example, the strong current dependencies of the arc voltage, cathode- and anode fall voltages for a low current non LTE arc are plotted in figure 3.5. It can be seen, how the cathode heating required to reach thermionic emission temperatures is provided by the increasing fall voltage (located in the space charge layer) for total arc currents below 10A.

3.1.1.5  Anode hot spot and fall voltage

Analogous to the situation at the cathode, the anode attachment can be diffuse, i.e. the current is spread over a relatively large area at a density of about 10^5¼7 A/m². The energy flux densities are not very high compared to the spot mode and the material erosion is small as long as the evaporation temperature of the anode material is not reached.

If the increasing current is forced to occupy the edges or the anode surface structure is non-homogenous, the current density can increase locally by several orders of magnitude forming a single or numerous small hot spots.

Again, directly in front of the anode surface, a small space charge layer develops. If there is no need to supply additional energy for sustaining the electron current, i.e. the current density is below the electron saturation current density of the undisturbed column plasma, the fall voltage becomes negative.

In the case the anode surface area is smaller than the cross section of the column, there is also a plasma constriction near the anode resulting in an additional geometric fall voltage.

In most cases, the anode fall voltage does not vary appreciably with total discharge current (see figure 3.5).

3.1.2  Modelling tasks

In order to reproduce the basic features summarized above and to get accurate modelling software, at least the following physical phenomena have to be included in the discharge model:

1. Heat conduction within the (solid) electrodes.
2. Electron emission from the cathode surface.
3. Electrical and thermal transition from the electrode surface to the
   (equilibrium) thermal plasma.
4. Current and heat transport within the arc plasma column
   (plasma description as a conducting fluid).
5. The symmetry of arc discharges is at least cylindric
   - an overall one-dimensional description is not possible.

As elaborated in the following sections, there are different models for different regions of the arc plasma. Very thin (skin like) non equilibrium layers in front of the electrodes are of special importance for the existence and formation of the arc itself. At least one of these layers, the space charge sheath, is small enough to be treated one-dimensional. Summarizing the results, the overall arc discharge model will consist of three sub-models linked by an iterative procedure:

• 2-D description of the heat conduction within the electrodes: section 3.2.
• 1-D description of the cathode sheath and presheath: section 3.3.
• 2-D description of the arc plasma (current transport, heat and fluid flow): section 3.5.
• 1-D description of the anode sheath and presheath: section 3.4.

Within the iteration algorithm, the data of the plasma and electrode surface is transferred to new boundary conditions for the next iteration cycle. The concept is thus called transfer function approach (see section 2.4). This approach can be implemented into the numerical simulation software independently of the specific details of the layer (1-D) or plasma/electrode (2-D) modelling and independent of the numerical procedures to be used for the 2-D regions (like finite volume or finite element methods).

3.2  Physical processes inside the electrodes and at the electrode surface

The electrodes are heated within a small current-carrying area, the so-called hot spot. This heat load is dissipated by conduction and thermal radiation of the surface. These losses and sources are boundary conditions for the heat conduction equation to be solved within the electrodes:

rcp ¶T ¶t = ® Ñ    · é ë l(T) ® Ñ    T ù û + ® j 2 s(T)  .

The Joule heating term was included into the calculations, but was found to be negligible for most high pressure discharge situations, while it becomes important for very high current densities.

A typical welding cathode tip is shown in figure 3.6. The active (electron emitting area) is easy to identify, because the electrode surface is influenced by the high heat loads during arcing. Because electron emission workfunction is one of the most important electrode material parameters, real electrodes are often doped with low workfunction materials. As a consequence, surface and bulk diffusion of the emitter material [112] becomes important for the overall discharge behaviour and cathode lifetimes. Such non homogenous surface states imply spacially non uniform thermal radiation emissivities and work functions. In principle, the local radiation emissivity and electron emission work function therefore depends on local temperature and electrode composition. For a detailed investigation of the cathode surface structures and processes see [111].

The heat and electric conductivity of typical electrode materials are plotted in figure 3.7. They are several orders of magnitude larger than the plasma conductivities. As a conclusion, current transport within the electrodes is not important for the overall discharge behaviour, while electrode geometry and heat conductivity significantly influences the hot spot formation at the cathode and thus the cathode fall voltage.

3.3  The cathode-plasma transition layers

As discussed above and sketched in figure 3.8, the transition from the electron emitting cathode surface to a plasma in (partial) LTE can be further divided into a space charge layer (sheath) and an ionization layer (presheath). The interface between the thermal plasma and the presheath is associated with the subscript P, physical quantities at the sheath-presheath interface will be identified by the subscript SE, surface quantities by the subscript S. The sheath and presheath sublayers and their linkage are described in the following sections.

3.3.1  Electron emission of the cathode surface

The local electron emission of the cathode surface depends on the local surface temperature T_S and the material property work function F. It is given by the Richardson emission formula [106]:

jRS = AR TS2 exp æ ç è - Feff kbTS ö ÷ ø

with the effective work function

Feff = F+ dFRS

and the Richardson constant A_R. Its theoretical value is 4 pk_B² m_ee/ h³ » 1.2·10⁶A/(K² m²).

For practical applications, the influence of the surface electric field E_S on electron emission is given by the Schottky correction formula (in eV) [4]:

dFRS = -   æ  ú Ö eES 4 peo

The work functions of some non refractive cathode materials used for high intensity arc discharges are summarized in table 3.1. The F and A_R values in this table where determined experimentally (best fit to the Richardson equation). There seems to be a noticeable variation of A_R. The values provided for doped tungsten materials are estimated and may also depend on the surface state of the electrode and the local surface concentration of the dopand.

Secondary electron emission by the impinging ion current j_i is given by j_SEE = g_i ·j_i where g_i = 0.01¼0.3. It can be neglected for most arc discharge configurations. For low pressure discharges, g_i becomes important. For argon with tungsten cathodes a value of 0.07 can be estimated [100].

In the case of high surface electric fields and cathode hot spot current densities, the emission current has to be calculated by the full thermo-field emission theory (see figure 3.9 and [99]).

The strong temperature and work function dependence of the emission current density determines the cathode hot spot temperature and may result in a bifurkation of the cathode hot spot system: Choosing optimum cathode geometry and a reasonable large discharge current, the arc contricts to a diffuse hot spot with current densities of 10⁶ ¼10⁸ A/m² and several hundred microns in diameter (diffuse mode).

If the cathode is cooled by external means or by geometrical design, the hot spot may constrict below a critical current density and field electron emission becomes increasingly important (spot mode).

3.3.2  Space charge layer (sheath)

For the discharge the space charge layer (sheath) has the following functions:

• limiting the back diffusion of the plasma electrons by a retarding electrical field which also isolates thermally the electron fluid from the surface - i.e. the diffusion of the electrons from the plasma is inhibited by a potental barrier,
• accelerating the ions heating the surface to thermionic emission temperatures,
• accelerating the emitted electrons in order to gain energy for ionization in the presheath,
• increasing electron emission by the high electric field at the surface.

The detailed potential and electric field distribution within the sheath is not relevant to the hot spot formation or overall discharge behaviour. The most important quantity, the voltage drop accross the layer U_S, is calculated by the energy balance at the sheath edge in section 3.3.3.

For the calculation of the electric field at the surface a simplified model is used. In the case of a collision-free space charge sheath (l_ia >> l_Debye), the Poisson equation for the sheath can be solved analytically for the electric field strength at the cathode surface (see [84,125]):

ESMK =   æ  ú  ú Ö 4 eo æ ç ç ç è jion ·   æ  ú Ö miUS 2 Zeff e   - jRS ·   æ  ú Ö meUS 2 e   ö ÷ ÷ ÷ ø    .

Within an accuracy of Ö{m_e/ m_i} this simplifies to the so called McKeown formula

ESMK = æ ç è 8  Zeff  mi jion2  US e eo2 ö ÷ ø [ 1/ 4]    .

The influence of multiple charged ions was taken into account by using the effective ion charge Z_eff. As a conclusion, the electric field at the cathode surface depends on the voltage drop accross the sheath U_S and the ion current density at the sheath edge j_ion. For a collision-dominated sheath (ion atom mean free path l_ia << l_Debye), the electric field can be also be calculated by integrating a simplified Poisson equation for the sheath [123,80]

ESW = æ ç ç ç è 5  jion  US 3  eo x   æ  ú Ö mi e lia   ö ÷ ÷ ÷ ø [ 2/ 5]

with x = 1.143.

For the overall current balance, the electron back diffusion current is calculated by

je,SEback = 1/4 e ne ve,SEth ·exp æ ç è - e US kbTe,SE ö ÷ ø

where T_e,SE is the electron temperature at the sheath edge and n_e the electron density at this position. The mean thermal velocity is given by

vth =   æ  ú Ö 8 kbT pm    .

3.3.3  Energy balance and current continuity at the sheath edge

As known since the beginning of the century [117], one of the main purpose of the space charge layer is to accelerate the emission electrons in order to sustain the ionization within the presheath. Approximately, the energy needed is given by the ion current density at the sheath edge times the ionization energy j_i,SE ·E_ion / e . Within this estimate, the sheath voltage drop is given by

US = ji  Eion e jRS

Benilov and Marotta [15] formulated a more accurate electron energy balance at the sheath edge which is used for this work:

There are two sources, the flux of energy brought into the layer by the emitted electrons accelerated in the space charge sheath j_RS/e·(2 k_bT_S + eU_S) and the work of the electric field over the electrons inside the layer calculated from the mean presheath electron current to [`(j_e,PS)] ·U_PS, where U_PS is the voltage drop accross the presheath. In the following the electron temperature T_e means the electron temperature at the sheath edge.

The energy losses at the sheath edge are coming from the electron back diffusion current leaving the presheath for the sheath j_e,back / e·(2 k_bT_e+ eU_S) and the electron current leaving to the bulk plasma j_tot ·(⁵/₂ + D_e^T) k_bT_e/ e and the losses due to inelastic collisions required to sustain the ion current j_i / e·E_ion. Together with the current continuity condition, the following equation system has to be solved for the current density j_tot and the space charge sheath voltage drop U_S:

jtot = jion - je,back + jRS (3.1) jRS/e·(2 kbTS + eUS) + je,PS    UPS = je,back / e·(2 kbTe+ eUS) (3.2) + jtot / e·(5/2 + DeT) kbTe + ji / e·Eion

This equation system is solved for all positions along the cathode hot spot surface.

For single charged ions, using D_e^T » 0.7 and [`(j_e,PS)] = ¹/₂ (j_tot + j_RS - j_e,back), the space charge layer voltage drop can be calculated for a given ion current densitiy j_ion, emission current density j_RS and electron back diffusion current density j_e,back, sheath edge electron temperature T_e and presheath voltage drop U_PS

US = jion ·( Eion + 3.2 kbTe- 0.5 eUPS ) e·(jRS-je,back) - jRS ·(2 kbTS + eUPS - 3.2 kbTe) e·(jRS-je,back) - je,back·(1.2 kbTe- eUPS ) e·(jRS-je,back)  .

3.3.4  Ionization layer (pre-sheath)

A detailed one-dimensional model of the cathode layer was developed by Rethfeld and Wendelstorf [105] based on the work of Hoffert and Lien [62] also continued by Hsu and Pfender [66,64]. The presheath was identified as a stiff boundary value problem. The spatial variation of the plasma parameters can be adapted to almost every set of plasma and sheath boundary conditions.

For modelling the overall cathode hot spot self consistently, only two parameters are important: the presheath potential drop and the sheath edge ion current density. First, the actual ion current j_i at the sheath edge can be calculated from the local current density in the plasma j_P, the electron emission current density j_RS and the electron back diffusion current j_e,b

ji = jP - jRS + je,b  .

The velocity of the ions at the sheath edge is given by the Bohm criterion [19,107]:

vi, Bohm =   æ  ú Ö kb( Ti+ Zeff Te) mi    .

The effective ion charge Z_eff is again used as a crude approximation for the influence of multiple charged ions. Using this ion current, the electron and ion density at the sheath edge can be calculated:

ni,SE = ne,SE = jion e vi,Bohm  .

Assuming the electrons are in Boltzmann equilibrium with the presheath electric field, the presheath potential drop can be calculated:

eUPS kbTe = log ne,P ne,SE  .

(3.3)

3.3.5  A simple presheath control model

The cathode layer model sketched above allows for the calculation of all discharge relevant data without prior knowledge of the detailed spatial variation of the physical quantities in the presheath. For controlling the physical validity of the results, the presheath voltage drop can be computed by the model of Benilov and Marotta [15]. The presheath voltage drop is again calculated from equation 3.3, but the sheath edge electron density is calculated by a presheath model to

ne,SEBM = ne,P · 0.8 2+a

with

a =   æ  ú Ö kbTh miDia arek ne,P2    .

The heavy particles temperature T_h and ion-atom diffusity D_ia are calculated for an average presheath heavy particles temperature T_h. The recombination coefficient a_rek is calculated in section 4.10. Using this sheath edge electron density, the sheath edge ion current density becomes

jSEBM = e ne,SEBM ·vBohm  .

The ion current density calculated from this equation agrees perfectly with the ion current densities calculated self consistently using the scheme from the previous subsection (see also figure 5.9).

3.3.6  Sheath-surface boundary conditions

The local net heating or cooling of the cathode surface is given by a summation over all possible contributions

qS = qS,i + qS,em + qS,rad + qS,e,back + qS,a + qS,a,back

(3.4)

with

qS,i = ji / e·( 2 kbTh,SE + 1/2 Zeff kbTe,SE + Zeff eUS + Eion - Zeff Feff ) (ion heating) qS,em = - jRS  Feff (emission cooling) qS,rad = eSB CSB TS4 (thermal radiation cooling) qS,e,back = je,back / e· ( 2 kbTe,SE + Feff ) (electron back diffusion heating) qS,a = - la ® Ñ    Th (neutrals heat conduction) qS,a,back = - ji/ e·(2 kbTS) (back flow of recombined ions) (3.5)

The Stefan-Boltzmann constant is given by C_SB » 5.67·10^-8 W/m² and the temperature dependence of the (tungsten) surface emissivity is given by [131]:

eSB = -0.0266 + 1.8197·10-4 ·TS - 2.1946·10-8 ·TS2

where the local surface temperature T_S is given in K.

Within the cathode hot spot, q_S is mainly given by ion heating balanced by emission cooling. The typical behaviour of this net energy flux density for different values of the sheath voltage drop is plotted in figure 3.10. From these plots, one can estimate the upper limit for the cathode hot spot temperature that is possible for a specific work function of the cathode material: 3630¼4610K for thoriated tungsten (F = 2.7eV) and 4340¼4905K for pure tungsten (F = 4.5eV) [43].

3.3.7  Presheath-plasma boundary conditions

While the boundary conditions for the cathode surface are determined by equation 3.4, the electron energy balance of the plasma in front of the cathode surface is dominated by the Joule heating term and thus the major boundary condition is coming from the calculation of the local layer current density (equation 3.1).

Additionally, the following source term is included into the energy balance within the cathode layer:

S = jRS (2 kbTC / e+ US) + jtot UPS -ji Eion / e -je,b (2 kbTe,SE / e+ US)

Within the anode layer, the source term

S = -je,b (5/2 kbTe,SE / e+ US)

is used.

3.4  The anode-plasma transition layers

Analogous to the cathode layer, the anode layer is divided into a sheath and a presheath structure. The distribution of the electric potential within this system is a bit different to the one at the cathode. As sketched in figure 3.11, the voltage drop of the space charge layer is negative. Assuming the anode is not emitting electrons, the current transport is mainly due to the plasma electrons diffusing to the anode surface against the (negative) space charge voltage drop (the nomenclature is the same as for the cathode):

je,SEback = 1/4 e ne ve,SEth  exp æ ç è - e US kbTe,SE ö ÷ ø

The voltage drop U_S is about several volts (3-8V) because the electron saturation current (U_S = 0) is above the usual current densities in the plasma.

As can be seen from figure 3.11, there must be a location in the layer, where the electric field is zero. At this position, the electric current is only diffusive (see also section 3.5.3):

je = s· kbTe e · ® Ñ    ne ne  .

Assuming a linearly descending n_e(x), the sheath edge electron density becomes

ne,SE = ne,P · æ ç è 1+ e dlayer  jtot kb s Te ö ÷ ø -1    .

The presheath voltage drop U_PS is then calculated by equation 3.3 and the sheath voltage drop U_S reads

US = kbTe e log æ ç è e  ne,SE  veth 4 jtot ö ÷ ø  .

Finally, the local net heating of the anode surface is analogous to the cathode given by

qS = qS,rad + qS,e,back + qS,a

(3.6)

with

qS,rad = eSB CSB TS4 (thermal radiation cooling) qS,e,back = je,back / e· ( 2 kbTe,SE + Feff ) (electron back diffusion heating) qS,a = - la ® Ñ    Th (neutrals heat conduction) (3.7)

The surface energy balance within the hot spot region of the anode is dominated by condensation of the current carrying electrons at the surface (q_S,e,back).

Compared to the cathode layer model in section 3.3, this anode layer model is rather crude. For enhancement, there may be the need to perform a detailed local 1-D calculation instead of the pure analytical model used here.

Additionally, the role of energetic photons (UV resonance radiation, see also [35]) in the cathode and anode layers is currently not included into any self-consistent model, but can change the ionization balance in the presheath. The models used in this section does not require a detailed calculation of the presheath properties.

3.5  The thermal arc plasma

Electric arcs are rather complex entities because of their spatial inhomogenities and the amount of plasma physics necessary for a description of the discharge. While low pressure discharges require a more or less kinetic treatment, because at least the electron distribution function is non-Maxwellian, the thermal plasma of the arc can be treated by a multi-fluid approach. Additionally, the full multi-fluid theory [24] can be simplified to a one-fluid treatment regarding flow phenomena and a coupled two-fluid model regarding energy transport.

While the layer physics presented in the preceeding sections, provide the boundary conditions for the thermal arc plasma, the multidimensional description of this region (often called the arc column, but also valid for the near electrode plasma regions) is provided in the following subsections.

The thermal arc plasma is treated as a quasineutral (n_e = n_i) conducting fluid. Except for the energy balance, the electrons and heavy particles are treated as a single fluid.

3.5.1  Modelling the flow of the conducting plasma fluid

3.5.1.1  Mass continuity

3.5.1.2  Momentum balance

From Newtons second law, the momentum equation is derived ([97] and [18]):

The total time derivative is defined as usual¹. The strain-rate tensor \sf p is described in [97]. For cylindrical coordinates, \sf p can be found in [18]. The dynamic viscosity h is the proportionality factor between momentum current and velocity gradient and is calculated and discussed in section 4.6. Natural convection forces become important for high pressure arc lamps or gas shielding effects occuring in welding applications.

Regarding momentum transport, the type of the flow is determined by the viscosity h, velocity u and characteristic dimension l of the problem. From these quantitites, a characteristic dimensionless number can be derived, the Reynolds number:

Re = r u  l h

Flows with equal Reynolds numbers are similar. In thermal plasma gas discharges, Re is not small ( > 1), but well below the critical values for turbulence ( » 1000).

3.5.2  Modelling the energy transport within the plasma

3.5.2.1  Equilibrium (LTE) heat balance

Heat is generated in the plasma by Joule heating and carried away by radiation (net radiation loss coefficient S_R, section 4.9), conduction (heat conductivity l, section 4.7) and convection. The compression work [(v)\vec]  · [(Ñ)\vec] p and viscous dissipation [(v)\vec]:\sf p can be neglected, because the Mach-Number is about 0.1 and Prandtl-Numbers² Pr = c_p h/l are small (0.03¼0.7). Because the total pressure in the discharge region is mainly constant, the energy balance is formulated in terms of the plasma enthalpy h and reads

(For the Joule heating term one can often find the equivalent forms j²/s or sE²). The heat flux is given by conduction, enthalpy transport by the (current carrying) electrons and the diffusion-thermo (Dufuor) effect from the (de)mixing of different plasma components [44]:

® q   = l cp   ® Ñ    h + 5kb 2 e   ® j cp  h + ® Ñ    · é ê ë å i  æ ç è rD - l cp ö ÷ ø ( hi ) ® Ñ    Ci ù ú û .

(3.10)

For simple one gas component plasmas, the Dufuor effect vanishes and for numerical convenience, the temperature formulation is introduced ([(Ñ)\vec] h = c_p [(Ñ)\vec] T):

rcp D Dt T = ® Ñ    · æ ç è l ® Ñ    T + 5 kb 2 e   ® j   ·T ö ÷ ø + ® j   · ® E   - SR

Even for modelling discharge configurations showing negligible deviations from LTE, the electron enthalpy flow term 5 k_b/ (2 e) [(j)\vec] ·T_e becomes undefined near the cathode layer, where a temperature split (T_e ¹ T_h) takes place in every discharge situation and the LTE model has to fullfill boundary conditions at the electrode coming from the heavy particle / electrode surface interaction³, while the electrons are thermally isolated by the space charge layer voltage drop.

3.5.2.2  Two temperature (pLTE) heat balances

Typical for thermal plasmas is the phenomenon of temperature split at pressures below 0.1 MPa, at temperatures below 9000K and near the electrodes. The electrons become decoupled from the heavy particles. For this reason, accurate modelling requires the electrons and the heavy particles to be treaten as different fluids (two-fluid model). The energy balance of the electrons becomes

re D Dt he = ® Ñ    · ® qe   + ® j   · ® E   - SR - . Eeh

(3.11)

with the electron heat flux

® qe   = le cp ® Ñ    he + 5 kb 2 e ® j cp  he

(3.12)

and the energy exchange term [(E)\dot]_eh as discussed in section 4.4.

Analogously, neglecting transport of ionization energy [75], the enthalpy balance of the heavy particles becomes

r D Dt hh = ® Ñ    · ® qh   + . E   eh

(3.13)

with the heavy particles heat flux

® qh   = lh cp ® Ñ    hh + ® Ñ    · é ê ë å i  æ ç è rD - l ch ö ÷ ø ( hi ) ® Ñ    Ci ù ú û .

The heat conductivity is discussed in section 4.7. Incidentally, one has to distinguish between the general (3D) formulation given above, the (published) 2-D equations and the equations really treated by the individual computer code.

For this work (chapter 5), the full set of the equations is solved (without the Dufuor term, because no complex gas mixtures are investigated).

3.5.3  Modelling the electric current transport

Regardless the importance of hydrodynamic transport phenomena, the arc is mainly an electrical phenomenon. Taking into account that the plasma column is electrical neutral (quasineutrality, r_el » 0), the generalized Ohm's law becomes [47,122]

® j   = s· ì í î ® E   + ® v   × ® B   - 1 nee ® j   × ® B   + 1 2 ene ® Ñ    pe + 0.41· kb 2 e ® Ñ    Te ü ý þ ,

(3.14)

The second and third terms (induction and Hall current) can be safely neglected. The fourth term, the diffusion current can influence the current transport within the electric arc [47]:

• In the outer regions of the arc, there is a diffusive current towards the axis and thus a space charge compensating electric field.
• On the axis, the diffusion current contributes about 10% to the overall electric current density.
• Near the anode, the diffusive current density can become as large as the normal one.

Because of the numerical difficulties of using this generalized Ohm's law and because diffusion current contribute less than 10% to the total current density (except near the anode), most authors use the conventional Ohm law

® j   = s ® E

where s denotes the electric conductivity depending on the local electron temperature. Recently, 2-D modelling results without this simplification where published [2]. The results give additional evidence of the diffusive current transport for the plasma anode transition region (anode attachment modes), while it can be neglected in the main plasma column and the cathode attachment region. Thus, current continuity law

Ñ· ® j   = 0

(3.15)

gives a Laplace-type equation

® Ñ    · æ è s ® Ñ    F ö ø = 0.

(3.16)

for the electric potential F. This equation is solved in the main arc plasma throughout this work.

3.5.3.1  Magnetic field calculation

Starting with Ampere's law

® Ñ   ×  ® H   = ® j   + . ® D

and using

Ñ· ® B   = 0   ,    ® B   = m0mr ® H      ,    ® B   = ® Ñ   ×  ® A

the well known equation for the vector potential [(A)\vec] is found in the form

\triangle ® A   = - m0 ® j

with the solution

® A   ( ® r   ) = m0 4 p ó õ ® j   ( ® r¢   ) | ® r   - ® r¢   | d3 r¢  .

The displacement current may be neglected as |[([(D)\vec])\dot]| << |[(j)\vec]| and the Coulomb gauge

® Ñ   ×  ® A   = 0

is used. Application of curl yields Biot-Savart's law

® B   ( ® r   ) = m0 4 p ó õ ® j   ( ® r¢   ) ×( ® r   - ® r¢   ) | ® r   - ® r¢   |3 d3 r¢  .

For this work, the plasma geometry is radial-symmetric and the axial current density component j_z is much larger than the radial one j_r. The magnetic field is simply given by [(B)\vec] = {0,0,B_q }. Neglecting the displacement current and with m_r » 1, Ampere's law reduces to

1 r ¶ ¶r (r Bq) = m0jz .

This equation can be integrated from r = 0 (where B_q = 0) to r = r. Finally, the azimutal magnetic field B_q(r,z) is obtained by a simple integration over the axial electric current density:

Bq(r,z) = m0 r ó õ r 0  jz(r,z)  r dr .

(3.17)

3.5.4  Diffusion and other transport phenomena

3.5.4.1  Ambipolar diffusion

In high pressure discharges, diffusion of charged particles is not free. There is a strong coupling between the electrons and ions. This phenomenon is called ambipolar diffusion. For the electric arc, the corresponding diffusity D_amb is almost identical to the ion-atom diffusity D_ia (see section 4.5). The additional equation to be solved within the (2 or 3D) plasma region reads:

D Dt ne+ ® Ñ    · é ê ë e neDia kbTh æ ç è ® E   - 1 e ne ® Ñ    pi ö ÷ ø ù ú û = . n   e   .

The heavy particles (ions and atoms) are described by a single local temperature T_h = T_i = T_a. p_i = n_ik_bT_h is the partial pressure of the ions (n_i = n_e, [41]) and [(n)\dot]_e the net ionization rate (see section 4.10). The major application of this additional diffusion equation is detailed modelling of the anode attachment [2].

An alternative approach more applicable for arcs producing large plasma jets or influenced by external gas flow, is to model the neutral atom balance equation

D Dt n0 = - . n   e

together with Dalton's law and charge neutrality condition (see section 4.2 and [59]).

3.5.4.2  Turbulent plasmas

Some applications like plasma gas cutting result in turbulent plasma flows. The inclusion of turbulence into the arc column model is reviewed in [127] and not used or discussed here.

3.6  Summary of the arc discharge physics

Applying the physical models of the discharge regions sketched above, means making a number of basic assumptions, which are well justified and accepted:

• The fluid flow in the arc column is laminar (Re < 400).
• The arc column is in partial local thermodynamic equilibrium (pLTE), i.e. locally determined by the variables velocity [(v)\vec], electron temperature T_e, heavy particles temperature T_h, pressure p and the electromagnetic quantities current density [(j)\vec], electric potential F and magnetic field [(B)\vec]. This is called two-fluid (T_e([(r)\vec]) ¹ T_i,a º T_h([(r)\vec])) hydrodynamic treatment.
• The arc column plasma is quasi-neutral, i.e. free of space charges.
• The non equilibrium layers (sheaths) in front of the electrodes are small against the overall electrode arc attachment region - they are skin like. This is always valid for the space charge layer (sheath), but sometimes questionable for the ionization layer (presheath).

These assumptions are valid for most discharge situations not too far away from LTE. This physical model requires additional computation of the plasma composition and transport coefficients depending mainly on electron temperature, as described in the next chapter.

Finally, the model will allow predictions for the overall discharge behaviour and the detailed spatial distributions of the physical parameters like plasma and electrode temperature. Results for a number of typical applications are provided in chapter 5.

It is evident from the different physical processes governing the main arc plasma and the sheaths, that the discharge cannot be described by a single set of equations. For physical self consistence, a proper definition of the internal boundaries and an iterative linkage of the submodels is cruical - as provided by this and the preceeding chapter.

Chapter 4
Plasma properties and transport coefficients

4.1  Introduction

The composition, thermodynamic properties and transport coefficients of thermal plasmas strongly vary with temperature. Their calculation in local thermodynamic equilibrium (LTE, T_e = T_h) is based on the work of Hirschfelder, Curtis and Bird [61] using tabulated partition functions [42]. Significant efforts where dedicated to the transport theory foundations since the early work of Devoto [38]. Starting 1977 [70], non equilibrium effects within the partial LTE concept (pLTE, T_e ¹ T_h) were introduced [64,20].

Focussing the efforts on modelling the overall arc discharge, the results of these references were selected and applied on the basis of the current knowledge from major working groups. The approach is to focus on LTE dicharges, but to use the two-temperature model in order to identify parameter ranges, where identical electron and heavy particles temperature in the main arc plasma (column) is not a valid assumption.

An ab initio decription of arc discharges more far away from LTE requires a closer look to the excitation equilibria [121] and is regarded not to be within the range of current modelling resources.

As summarized below, the plasma composition and transport coefficients show a strong variation with electron temperature, while the effects of temperature split (Q = T_e/T_h) > 1 are relatively weak. This does not justify the use of LTE models for discharges not in LTE, because deviations of electron temperature from the heavy particles temperature will imply important variations of plasma parameters and transport coefficients through the T_e dependencies. As an example, the LTE electric conductivity of an atmospheric argon plasma at 5000K is almost zero, but allowing for an electron temperature of 10000K will give a conductivity of 2000 S/m. Such temperature split effects are important for all electric arc discharges within a small layer in front of the electrodes, where LTE is almost impossible (except for some sodium vapour lamps) because the heavy particles equilibrate with the surface but the electrons are thermally isolated through a potential barrier (space charge sheath). Using a two temperature model will permit a self-consistent computation of this thermal boundary layer and enables the inclusion of electron enthalpy flow into the overall energy balance, because this term requires the electron temperature gradient, especially near the electrodes.

In the following sections, the plasma composition and thermodynamic properties are calculated (section 4.2) and the basic collision integrals are provided (section 4.3). These are needed to calculate the transport coefficients used within the macroscopic models discussed in chapter 3:

• The energy exchange rate coefficient linking the electron and heavy particles energy balances (4.4).
• The diffusities required to calculate demixing, ambipolar diffusion and used as a basis for the other transport coefficients (4.5).
• The dynamic viscosity required to calculate the hydrodynamic plasma flow (4.6).
• The electron- and heavy particles heat conductivity (4.7).
• The variation of electric conductivity with electron temperature (4.8).
• The net plasma radiation emission coefficient is briefly discussed (4.9).
• The net recombination coefficient needed to calculate deviations
  from ionization equilibrium (4.10).

4.2  Plasma composition and thermodynamic properties

The foundation for the calculation of all plasma transport coefficients is the determination of the plasma composition as a function of electron temperature T_e and heavy particles temperature T_h. The particle densities of the individual components are determined by the condition of quasineutrality, Daltons law and the Saha-Eggert equations.

4.2.1  The condition of quasineutrality

Deviations from electrical neutrality (space charges) will produce high electric fields efficiently restoring electrical neutrality. On macroscopic scales, the plasma is thus electrical neutral (quasineutral). The electron density is equal to the sum of all ion densities times their ionization state:

ne = å (i,j)  Zi,j ·ni,j

where i is the sum over all plasma components and j the sum over all ionization states taken into account.

4.2.2  Daltons law

Real gas effects (virial corrections) can be neglected within the temperature and pressure range of thermal plasmas [23]. The total discharge pressure is described by the ideal gas law. The summation of the contributions from the individual plasma components gives:

p = å j  nj kbTh+ nekbTe  .

Due to the effectiveness of charge-exchange collisions and the low ionisation degree of most high pressure discharges, ion temperature is assummed to be equal to the gas temperature.

At total discharge pressure levels above 10 MPa, the plasma becomes non-ideal. For a weak non-ideality, the Debye pressure correction can be introduced [55]

pDebye Correction = - e2ne 12pe0lDebye

The total discharge pressure is given by external conditions and thus one of the input parameters. The hydrodynamic pressure variations are usually less than 1% and therefore neglected for the calculation of the plasma composition and transport coefficients.

4.2.3  Saha-Eggert equation

The ionization equilibrium of the plasma is described by a system of Saha-Eggert equations [121,120]

ne æ ç è nz nz-1 ö ÷ ø = Ze Zz Zz-1 exp æ ç è - Ez - dEz kbTe ö ÷ ø

describing the ratio of the densities of the different ionization states by a Boltzmann term with the excitation temperature equal to the electron temperature. This approximation is still under discussion [32,10], but well accepted by most groups. Some minor modifications of the excitation temperature can be introduced, but are not cruical for the overall discharge behaviour [53].

The partition function Z is given by its translational, rotational, vibrational and internal contributions:

Z = Ztr ·Zrot ·Zvib ·Zint

where rotational and vibrational contributions are unity for the monatomic gases used for this work. The translational partition function is

Ztr = æ ç è me kbTe 2 p(h/2p)2 ö ÷ ø 3/2

and the internal partition function is determined by the Energy E_n and degeneration d_n of the excited levels of the individual atom or ion:

Zint = En £ Ez - dEz å n  dn ·exp æ ç è - En kbTe ö ÷ ø .

The exact energy levels of the ions and atoms within the plasma cannot be computed ab initio, thus the undisturbed energy levels of the free atoms and ions are used. The divergence of the partition function is avoided by a so called cut-off criterion (figure ) [56], allowing to calculate the lowering of the ionization energy:

dEz = Z e2 4 peolDebye

with the Debye length of the plasma

lDebye = æ ç è e2 eokb å a  Za2 na Ta ö ÷ ø -1/2   .

(4.1)

The calculated particle densities are in good agreement with the values published in the literature (figure 4.2). The cut-off criterion has significant effects on viscosity and reactive heat conductivity only [29,90].

4.2.4  Mass density

As shown in figure 4.3, the mass density is sensitive to the non equilibrium parameter Q = T_e/T_h. This is due to Daltons law. As the total particle number increases with the onset of ionization, the ideal gas scaling r ~ 1/T is not valid above the onset of ionization (e.g. above T_e = 9000K for 0.1 MPa argon).

4.2.5  Enthalpy

The high temperature enthalpy (figure 4.4) is dominated by chemical and translational enthalpy, while excitational enthalpy can be neglected.

4.2.6  Specific heat

The temperature dependence of the specific heat reflects the different ionization levels within the plasma (figure 4.5).

4.3  Particle collision integrals

The basic particle interactions (collisions) within the plasma are modeled by the calculation of so-called collision integrals. The interaction between the different particles within the plasma (atoms, ions and electrons) can be described by an interaction potential depending on the relative (kinetic) energy of the collision partners. The calculations assume binary instantanous collisions. Details are described by Hirschfelder, Curtis and Bird [61]. The collision integral of the order (l,s) is defined by [7]:

W   (l,s)   = Q (l,s) p = 2(l+1) p(s+1)! [2l+1-(-1)l] × ó õ ¥ 0  exp(-X) Xs+1 Ql (E) dX

where

X º E/ kbT  , E º relative energy of the collision partners, Ql (E) º elastic collision cross section.

For a number of simple interaction potentials, the collision integrals are tabulated relative to the hard sphere values or given by analytic approximations. Otherwise, the collision integrals are computed by numerical integration. In the following subsections, the different collision categories are discussed.

4.3.1  Neutral particle interaction

Collisions between neutral particles are modeled by a Lenard-Jones interaction potential

V(r) = 4 kbtLJ [ (rLJ/r)12 - (rLJ/r)6 ]

using tabulated values for the parameters t_LJ and r_LJ [61]. The advantage of this crude approximation is the availability of the interaction potential parameters for a large number of different plasma gases.

Compared to other estimates, the collison integrals for temperatures above 1600K (for Argon) deviate by less than 20% (see table 4.2 and [76]). This will give slightly wrong values for the low temperature viscosity but has no remarkable effects regarding electrical and heat conductivity. For this interaction potential, the collision integrals can be related to the corresponding hard sphere values

Q   (l,s)   = prLJ2 · W(l,s)* æ ç è T tLJ ö ÷ ø

using the W^(l,s)* function tabulated in [61].

4.3.2  Charged particle interaction

Collisions between charged particles are modeled using a screened Coulomb potential. The collision integrals are tabulated [87,40]. Allowing for a deviation of about 15%, an analytical approximation can be used [78]:

Q   (1,s) ei  = 4 pb02 s (s+1) æ ç è log 2 lDebye b0 -2 g   + s-1 å i = 1  1 i - 1/2  ö ÷ ø Q   (2,s) ei  = 12 pb02 s (s+1) æ ç è log 2 lDebye b0 -2 g   + s-1 å i = 1  1 i - 1 ö ÷ ø

with [`(g)] = 0.577¼, b₀ = [(Z_i Z_j e²)/(8 pe_ok_bT_e)] and l_Debye from equation 4.1.

4.3.3  Ion atom interaction

Charge exchange collisions between ions and atoms can be modeled by the heuristic interaction potential

Qin(E[eV]) = Q0,CT · æ ç è 1 + aCT ·log E0 E ö ÷ ø 2

where Q_0,CT and a_CT are tabulated [4] (E₀ = 1eV). While the sources of these data are not known to the author, the same principal dependencies investigated by [104] and compiled by [130] are used for this work.

These collisions dominate the collision integrals [`(Q)]^(l,*) for even l. For odd l, Devoto has computed some values for argon. A comparison with the even l charge-exchange collision integrals gives a scaling factor of ¹/₄ , avoiding the extensive calculations done by Devoto and without the requirement of detailed atomic data. Alternatively, a simple polarisation model (dipol interaction) can be used [130].

4.3.4  Electron atom interaction

Electron atom collisions are highly affected by quantum mechanical matter wave inflection (Ramsauer effect). The cross sections can be measured [68] or calculated [3]. The collision integrals are calculated numerically using cross section data from the literature (see figure 4.6).

4.3.5  Summary of particle collision modelling

Using the simplifications discussed above and the fundamental data provided in table 4.1, reasonable agreement with the more detailed computations of Devoto where found (table 4.2). This gives a sufficient accuracy for a large number of plasma gases where the basic data needed is available, while more sophisticated calculations like those of Devoto will require unaffordable efforts or will be even impossible due to missing fundamental data.

4.4  Energy exchange term for the pLTE model

4.4.1  Elastic energy exchange rate

The elastic energy exchange between electrons and heavy particles is computed by a mean collision frequency approach [17,31]:

. E   elastic eh  = 3/2 kb( Te- Th) ·ne· _ v   th e  · å i  æ ç è 2 me mi · ni Q   ei  ö ÷ ø  .

Due to their small mass, the mean thermal velocity of the electrons is much higher and the heavy particles can be treated at rest. The sum is taken over all atom and ion plasma components and [`(v)]_e^th = Ö{8k_bT_e/ (pm_e)} is the electron mean thermal velocity.

The energy exchange (see figure 4.7) between the electron fluid and the heavy particles drastically increases with pressure.

4.4.2  Inelastic energy exchange rate

With increasing electron temperature, inelastic energy exchange between electrons and heavy particles should be taken into account. The detailed computation of this rate coefficient is extremely difficult. The principal dependency is

. E   inelastic eh  ~ neTe-1/2  Ek

where E_k is the atomic or ionic energy level to be excited [17]. Because these processes are important only for highly ionized plasmas, the total energy exchange rate is assumed to be increasing with temperature. The elastic decrease beyond a certain temperature (typically 15000K, see figure 4.7) is thus corrected.

4.5  Diffusities

The principal parameter dependencies of the diffusities are [130]:

D ~ 1 p ·Q   æ  ú Ö T3 m   » 10-5 m2/s

The binary diffusion coefficient of two plasma species is given by [61]:

Dijbin = 3 8 · Ö p p   Q   (1,1) ij  ·   æ  ú Ö (kbT)3 2m

with the reduced mass m = (m_i·m_j)/(m_i + m_j) and the total pressure p. From the binary diffusities, the total diffusity can be computed [61]:

Dij = N å k = 1  xk  mk mj · |Kji| - |Kii| | K |

using

Kij = ì ï ï í ï ï î xi Dijbin + mj mi å k ¹ i  xk Dikbin for i ¹ j 0 for i = j

|K_ij| is the sub determinant of the K matrix constructed by canceling the ith column and the jth row and multiplying with (-1)^j+i.

Ambipolar diffusity

Because macroscopic space charge formation is not possible (quasineutrality condition), the flux of electrons and ions out of any region must be equal (neglecting the net total current). This coupling avoids independent diffusion of charged particles. The coupled diffusity is called ambipolar diffusity and given by [79]:

DA = Di me + De mi me + mi  .

Using the Einstein relation for the mobility and diffusity ratio m_a/D_a = e/(k_bT) and the fact of much higher electron mobility m_e >> m_i, we get

DA = Di (1 + Te/Ti)

For practical applications in high pressure discharges, the ambipolar diffusity is almost identical to the ion-atom diffusity D_ia.

4.6  Viscosity

The proportionality factor between momentum current and velocity gradient is called dynamic viscosity. The principal dependence on plasma parameters is as follows [130]:

h ~   __ Öm T Q » 10-5 kg/(m s)

(viscosity is releated to the self diffusity by D/h = r). The ratio of the dynamic viscosity and the mass density n = h/ r is called kinematic viscosity. Here, we use the first approximation given by Hirschfelder et.al. [61]:

h = N å i = 1  xi2 xi2 hi + 1.385 N å [(k = 1) || ((k ¹ i))]  xi xk kbTh p  mi  Dik

(4.3)

where p is the total pressure, x_i the molar fraction of the component i and h_i the viscosity of the pure gas component i.

hi = 5 16 Ö   pkbThmi   · 1 Q (2,2) ii

The electronic contribution is less than 2 % and thus neglected [23]. Due to this small influence of the electrons, viscosity decreases with increasing non-equilibrium parameter Q (figure 4.8). The difference of the calculated viscosity to that provided by Boulos [23] is due to different treatment of the atom-atom interaction.

4.7  Heat conductivity

Heat transport in thermal plasma gas discharges can be split into an electron contribution l_e dominating the high temperature regime and a heavy particles contribution l_h dominating the low temperature regime:

l = le + lh

Using only the dominating translational heat conductivity, the principal dependence on plasma parameters is as follows:

l ~ 1 Q ·   æ  ú Ö T m       » 1 W/(m K)

4.7.1  Electron heat conductivity