A Villa

R Malgesini
email: roberto.malgesini@erse-web.it

L Barbieri
email: luca.barbieri@rse-web.it

A multiscale technique for the validation of a numerical code for predicting the pressure field induced by a high-power spark ‡

A more precise knowledge of the pressure field induced by a high-power spark is essential to estimate the mechanical damages that a lightning strike can induce near the impact point. In this work we propose a multiscale approach to validate a numerical magneto-hydro-dynamic (MHD) model that can predict the pressure field when a very high-power discharge is considered. Two simplified models for the arc resistance are considered and their respective results are compared. A brief analysis regarding the numerical issues involved in the solution of a very high temperature gas is included. The numerical code has been validated against the experimental data of a short-arc discharge using a current waveform prescribed by the aeronautical standards.

Our study shows that a strong shock wave is generated in the first power peak and this travels away from the arc column maintaining a relatively high strength a few tens of centimeters away. The pressure in the arc region remains high for the whole discharge period.

Introduction

The numerical and experimental study of the evolution of electric arcs is a very active field of the applied physics. This research topic has found many applications such as the combustion ignition devices [START_REF] Reinmann | Temporal investigation of a fast spark discharge in chemically inert gases[END_REF][START_REF] Ekici | Thermal and flow fields modeling of fast spark discharges in air[END_REF], the plasma torch design [START_REF] Caillard | Experimental investigation of the overpressure generated by a low energy plasma igniter[END_REF] and the lightning protection [START_REF] Larsson | The lightning swept stroke along an aircraft in flight. part 1: thermodynamic and electric properties of the lightning arc channels[END_REF][START_REF] Larsson | The lightning swept stroke along an aircraft in flight. part 2: numerical simulations of the complete process[END_REF]. Regarding this latter field a lot of efforts have been devoted to develop some lightweight lightning-resistant materials for the aerospace industry. In the past years carbon-fiber was massively introduced in many structural components. The composite materials have many advantages over aluminum alloys from the structural viewpoint. However they are more fragile and therefore they are more vulnerable to transient loads such as the pressure shock waves induced by lightning strikes [START_REF] Ginestre | Assesment of lightning direct effects damages by modelling techniques[END_REF]. In particular we are interested in the highest-power phase of the discharge, which is also called phase-A, [START_REF] Chemartin | 3d simulation of electric arc coloumn for lighting aeroplane certification[END_REF]. The lightning discharge can be roughly divided in a sequence of different phases characterized by different peak currents and time evolutions. Aeronautical standards, see for instance [START_REF]Aircraft Lightning Test Method[END_REF][START_REF] Ihs | Aircraft lightning test methods[END_REF], prescribe a sequence of standardized waveforms that represent, for each phase, a real lightning discharge. These current waveforms represent some of the lightning characteristics: for instance they are designed to reproduce the lightning-induced damages. In phase-A a maximum current of two hundred thousand Ampere is reached in a few tens of microseconds. The Joule effect triggers a rapid temperature rise and the development of a pressure shock wave that can damage the external structure of the airplane. In other terms the arc inception and the subsequent shock wave generation trigger a pressure difference between the internal and external part of the aircraft structural panels. This results in a localized compression force that can puncture them. For these reasons in this work we are interested in mapping the pressure field induced by a typical phase-A laboratory test [START_REF] Chemartin | 3d simulation of electric arc coloumn for lighting aeroplane certification[END_REF]. The simulation of a real lightning strike is beyond the scope of this work, we concentrate on the simulation of a typical laboratory test. This kind of experiment can reproduce the pressure wave generated by the arc inception. In particular the aeronautical standards fix the current peak and the time to peak values. In the literature there are many numerical and experimental approaches to predict the arc-induced pressure field. An example are the studies of [START_REF] Graneau | The cause of thunder[END_REF] and [START_REF] Blanco | Modelling digital thunder[END_REF] regarding the thunder generation and propagation. However, it is very difficult to carry out precise measures on real lightning strikes since the impact point is often unknown. In some cases, provided that some favorable conditions exist, it is possibile to trigger the strike using sounding rockets [START_REF] Dwyer | Energetic radiation produced during rocket-triggering lighting[END_REF] or lasers [START_REF] Wang | A possible way to trigger lightning using a laser[END_REF]. Even in this case, virtually, no lightning is equal to another. As a consequence, the thunder measures are very scattered. Therefore, in many cases, the study of smaller-scale laboratory-made electric arcs is preferred. For instance in [START_REF] Uman | Four-meter sparks in air[END_REF] a four meter arc was produced and many electrical and spectral data were collected. A few years later a comprehensive theory of long sparks was introduced in [START_REF] Plooster | Numerical model of the return stroke of the lightning discharge[END_REF][START_REF] Plooster | Numerical simulation of spark discharges in air[END_REF]. The pressure field associated with a 20kA discharge was computed and a good matching with experimental data was shown. In [START_REF] Akram | The evolution of spark discharges in gases: 1. macroscopic models[END_REF][START_REF] Akram | The evolution of spark dischages in gases: 2. numerical solution of one-dimensional models[END_REF] the lightning model of [START_REF] Plooster | Numerical model of the return stroke of the lightning discharge[END_REF] was corrected and enhanced, and a small variation in the pressure profile was outlined. Unfortunately the experimental device used in this kind of researches can produce arcs with a length of a few meters and peak currents of a few tens of thousands of Amperes which is not enough to reproduce the phase-A. In a more recent work [START_REF] Aleksandrov | Effect of continuous current during pauses between successive strokes on the decay of the lightning channel[END_REF] a 100kA current is simulated. A fast rise, with a time derivative of approximately 10 7 A/s, is imposed and the effects of the turbulence in the arc channel are studied in relation to the channel extinction time. However, in this work there is no direct comparison against experimental data. Similar studies have also been performed in other fields such as the design of plasma torches. In [START_REF] Caillard | Experimental investigation of the overpressure generated by a low energy plasma igniter[END_REF] the overpressure generated by a plasma igniter has been measured. The plasma torches use relatively low currents. Also the fast-spark regime, see [START_REF] Bayle | Blast wave propagation in glow to spark transition in air[END_REF], has been deeply investigated. These kinds of discharges are used, for instance, in spark plugs and are characterized by a very short time to peak. In [START_REF] Reinmann | Temporal investigation of a fast spark discharge in chemically inert gases[END_REF][START_REF] Ekici | Thermal and flow fields modeling of fast spark discharges in air[END_REF] the complete fluidodynamic field has been characterized using numerical and optical techniques. Although fast-sparks can generate quite intense shock waves, their peak current is not even comparable with a lightning strike. The lack of pressure measures near the impact point of a high-power arc is due to the fact that a direct measure, using pressure transducers, is challenging since the pressure gauges can be easily damaged even by a small percentage of the current involved in the phase-A discharge. Moreover the induced currents are so high that the transducer signal could be heavily distorted. A possible work-around is to move it to a safe position far enough from the discharge and to connect it to the point where the pressure has to be measured with a tube. This technique is widespreadly used in the aeronautic field, see [START_REF] Tavoularis | Measurements in fluid mechanics[END_REF], and it is particulary suited for static measures. However, when a transient phenomenon is involved, the pressure at one end of the tube is different from the one at the other end. In other terms the fluid dynamics inside the tube makes the instantaneous measure of the pressure on the impulsed panel impossibile. Nevertheless these data can still be used to validate a numerical tool. The computer fluid dynamics (CFD) can be a valuable tool to predict the lightningrelated overpressure. For instance in [START_REF] Plooster | Numerical simulation of spark discharges in air[END_REF] a long-gap spark in air is simulated using a finite difference fluid dynamic solver. A flux-corrected transport method, see [START_REF] Boris | Flux-corrected transport 1. shasta, a fluid transport algorithm that works[END_REF], is used in [START_REF] Akram | The evolution of spark dischages in gases: 2. numerical solution of one-dimensional models[END_REF] to simulate fast discharges in an axis-symmetric geometry. A comprehensive review on the numerical methods used for this kind of problem is also included. A coupled fluid-electric MHD model is described in [START_REF] Ekici | Thermal and flow fields modeling of fast spark discharges in air[END_REF] and [START_REF] Jens | Ab inition modelling of thermal plasma gas discharges (electric arcs)[END_REF] where also the Lorentz force is taken into account. The magnetic effects are quite important when high currents are involved. The CFD finite-element analysis has also been applied in [START_REF] Trelles | Finite element modeling of flow instabilities in arc plasma torches[END_REF] to study the flow instabilities in a plasma torch. Although there is a very large number of works covering the numerical simulation of discharges, to the best of our knowledge, a few papers make an exhaustive comparison with experimental data when a very high current is considered.

In this work we model the fluid field with a three dimensional fluid dynamic solver coupled with a one-dimensional solver for the tubes. This multi-scale numerical scheme can generate an estimation of the pressure recorded by the gauges. With this new approach we can directly compare the experimental results with the numerical ones and we can validate our numerical scheme in a typical phase-A scenario. With our calibrated numerical model we can reliably compute the pressure field near the impact point even when a very high-power strike is considered. We stress that our numerical and experimental results are obtained in a laboratory setting using the damped oscillation current waveform prescribed by the aeronautical standards [START_REF]Aircraft Lightning Test Method[END_REF][START_REF] Ihs | Aircraft lightning test methods[END_REF]. The validation of a numerical tool in this kind of environment is quite important since this can cut the number of experimental tests that are needed to qualify new aircraft components. Let us now give a brief outline of this work. In section 2 we describe the experimental device, in 3 we review a suitable MHD model and we introduce a simplified one-dimensional model for the tubes. In section 4 we introduce a simple time-splitting technique for the solution of our physical model. Finally, in section 5 we compare the numerical estimates with the measures and we show the pressure loads on the panel surface.

The experimental device

The test discharge circuit is depicted in Figure 1(a). This testing device produces a damped oscillation which is one of the waveforms prescribed by the aeronautical standards [START_REF] Ihs | Aircraft lightning test methods[END_REF] to represent the first phase of the real lightning discharge. A capacitor (C) is charged to 73kV then a high power spinterometer (S) closes the circuit and the current starts to flow. From the electrode (E) and the plate (P) a 5cm arc is formed. The inductance (I) represents the equivalent inductance of the electric circuit. This can be modeled with the RLC scheme, where the capacitance is 26µF , the inductance is 2.9µH and the resistance is 24mΩ. The impulse waveform is a damped oscillation

I = I 0 exp(-αt) sin(βt)
where I is the current, I 0 = 2.18 • 10 5 A is the maximum current, α = R/(2L) is the damping factor, β = √ ω 2 -α 2 is the frequency, ω = 1/(LC) is the undamped frequency and t is time. The resistance is the sum of the capacitor internal resistance, the spinterometer arc resistance, the wire resistance and the electrode-to-panel arc resistance. This quantity is estimated along with the circuit inductance by fitting a current measure plot as is done in Figure 2. As we can see, the fitted waveform is a good approximation of the measured one. The discharge lasts about 0.8ms, although most of the power is released in the first oscillations.

The plate is 1cm thick and its surface is 50cm x 50cm. Three holes are drilled into it 5cm (A), 10cm (B) and 15cm (C) from the center of the panel where the current is injected by a 5mm-thick electrode. On the other plate side a plastic tube, with a length of 25cm and a diameter of 1cm, is installed, see Figure 1(b). At the end of the tube a Kistler 4043A5, 4043A10 pressure transducer is installed. The gauge has been protected by a metallic box that houses also the amplifying equipment and the optical fiber transmission device. This solution minimizes the measurement errors triggered by the induced currents.

The tube diameter we have chosen is a compromise between the aerodynamic interference between the tube and the panel and the necessity to have a low viscous pressure decay inside the tubes. In fact a thinner tube means a thinner hole on the panel and a weak disturb for the fluid flow on the panel. Nevertheless in thin tubes the viscous effects are very important and this curtails the frequency content of the pressure signal. We have found that a tube of 1cm in diameter is a good balance.

The physical model

Geometry

We model our discharge device with a two-domain scheme, namely the plate domain Ω p and the air domain Ω a . This latter is a 50cm x 50cm x 20cm box, see Figure 3, that has a 5mm hole in its center representing the electrode. This latter tip is l 0 = 5cm distant from the plate. The boundary of the air domain is divided in several regions: the plate-air common boundary Γ ap , the electrode part Γ a,E , the top of the box Γ a,T and its lateral contour Γ a,L . The plate boundary consists of the following parts: the plate-air common boundary Γ ap , the plate lateral contour Γ p,L and the plate bottom Γ p,B . On the various boundary parts we will apply different boundary conditions.

Three dimensional MHD model

Let us now introduce the MHD model:

                                           ∂ρ ∂t + ∇ • m = 0 in Ω a × (0, ϑ f ] ∂ m ∂t + ∇ • m m ρ + P Ī = J a ∧ B a in Ω a × (0, ϑ f ] ∂E t ∂t + ∇ • ( m (E t + P )) = m ρ • J a ∧ B a + Π a -Π I in Ω a × (0, ϑ f ] -∇ 2 A i = µ 0 J i in Ω i × (0, ϑ f ], i = a, p m • n a = 0, A a = A p , ∇ A a • n a + ∇ A p • n p = 0 on Γ ap × (0, ϑ f ] m • n a = 0, B a = B 0 on Γ a,E × (0, ϑ f ] P = P 0 , A a = A 0 on Γ a,L Γ a,T × (0, ϑ f ] ∇ A • n p = 0 on Γ p,L × (0, ϑ f ] ∇ A • n p = 0 on Γ p,B × (0, ϑ f ], (1)
where ρ is the density, m = ρ u is the momentum density, u is the speed, E t is the total energy per unit volume, P is the pressure, Ī is the identity tensor, J i , with i = a, p, is the current density vector of the air and plate respectively, B i = ∇ ∧ A i is the magnetic field, A i is the magnetic potential, Π a = η| J a | 2 is the power due to the Joule's effect, η is the channel resistivity and Π I is the power radiation term. Finally, n i , with i = a, p, is the outward pointing normal of the air and the plate domain respectively, A 0 and B 0 are suitable boundary conditions we will discuss below and ϑ f is the final time.

The power radiation term Π I has been computed using the net radiation model. In other terms only the power loss is modeled and no radiative heating is considered. The power radiation has been considered function of pressure and temperature and has been derived from [START_REF] Ekici | Thermal and flow fields modeling of fast spark discharges in air[END_REF][START_REF] Reynard | Modelisation tridimensionell de l'amorcage de l'arc electrique dans un disjoncteur basse-tension[END_REF][START_REF] Ph | One-dimensional hydro-kinetic modelling of the decaying arc in air-pa66-copper mixtures: 1. chemical kinetics, thermodynamics, transport and radiative properties[END_REF].

Regarding the boundary conditions for the magnetic problem we have imposed the continuity of the magnetic potential and of the magnetic potential gradient across the plate-air interface. For the remaining surfaces we have imposed either the magnetic field or its normal gradient. A 0 has been computed using the Biot-Savart solution considering the electrode a one-dimensional current carrying wire. We have imposed some non-penetration conditions on the solid walls Γ ap and Γ a,E and we have imposed the atmospheric pressure P 0 on Γ a,L and Γ a,T . Model ( 1) is a magneto fluid dynamic problem for a gas in chemical equilibrium with an imposed, gaussian-shaped, current density:

J a (t, r) = -I(t) πr 2 0 e -(r/r 0 ) 2 ẑ J p (t, r) =    I(t)r 2πr 2 0 τ r if r ≤ r 0 I(t) 2πrτ r if r > r 0 (2)
where r = | x -x c |, x is the position vector, x c = (0.25, 0.25, 0)cm is the center of the panel, r 0 is the arc radius, ẑ = (0, 0, 1) is a vertical unit vector, τ = 1cm is the plate thickness and r = (

x -x c )/| x -x c |.
The gaussian current profile (2) decays quickly with respect to r so that the greatest part of the current passes through the channel, i.e r ≤ r 0 . Moreover it can be easily shown that ∞ 0 J a (t, r) = I(t) . In [START_REF] Reinmann | Temporal investigation of a fast spark discharge in chemically inert gases[END_REF] we have also disregarded the viscous effects. In fact since the Reynolds number ρU L µ (where ρ ≈ 10 -2 is the characteristic density, U ≈ 10 2 is a representative speed modulus, L ≈ 10 -2 is a reference scale and µ ≈ 10 -3 -10 -4 is the air viscosity) is quite high 10 5 -10 6 and the time-scales of the phenomenon are small, the viscosity has effects only near the boundaries where a boundary layer is formed. In our case the boundary layer is completely neglected, however it is well known that the pressure is almost constant inside it. In other terms we expect, in this case, that the Euler and full Navier Stokes equations produce similar pressure results. In the same manner we have also neglected the thermal boundary layer. In fact the panel and the air have a very complicated thermal coupling dominated by many processes such as the surface melting and the thermal radiation. Fortunately, on these time scales, the heat conduction is quite limited, therefore the thermal boundary layer is very thin. Also the convective heat transfer is very limited since in the boundary layer the speed is small.

Resistivity field

Since we are interested only in reproducing the pressure effects, we use a couple of approximated models for the resistivity η. The resistivity coefficient η can be either estimated from the experimental data or computed as a function of the temperature and of the pressure, see [START_REF] Reynard | Modelisation tridimensionell de l'amorcage de l'arc electrique dans un disjoncteur basse-tension[END_REF]. In the first case short-circuiting the electrode-gap, we can derive the resistivity of the circuit alone. Comparing this datum with the resistivity estimated in the case when the electric arc is present, we can, by difference, compute the arc resistivity which nearly equals R a = 8mΩ. From this we can easily estimate the volume-averaged resistivity as

η = (πr 2 0 /l 0 )R a . (3)
We do not expect any major contribution to the resistance from the electrode sheaths since the ratio between the arc voltage drop and the sheath voltage drop is quite high. In fact, since R a = 8mΩ, we get that the sum of the two contributions is over one hundred volts for most of the discharge, in particular in the first peak of current a maximum of 1600V is reached. On the contrary, the sheath voltage drop is a few tens of volts as described in [START_REF] Djakov | A model describing the arc between closely spaced electrodes[END_REF][START_REF] Yokomizu | Electrode sheath voltages for helium arcs between non -thermionic electrodes of iron, copper and titanium[END_REF].

From Figure 2 we see that the current interpolation based on a RLC-behavior produces good results. This means that, from a circuital viewpoint and in the discharge timeframe, the arc behaves as a constant resistance. We expect that the estimate (3), which prescribes a constant resistivity in time, is good enough to get a consistent prediction of the pressure field. In fact the resistivity field may exhibit some variations during the discharge and a transitory phase on the first current rising front. However this transitory is particulary quick, in fact, in our variable resistivity simulations, described below, we have reached a high conductivity > 2500Sm -1 in less than one microsecond. The maximum conductivity, > 10000Sm -1 , is reached several microseconds later. However this relatively strong variability has not a major impact on pressure measurements as we will see in the result section 5. We have also checked the length of the transitory phase using the Toepler law:

R a = k t l 0 t 0 I(τ )dτ (4)
where k t = 5 • 10 -5 is a constant. The Toepler law can be used, see [START_REF] Kuffel | High voltage engineering: fundamentals[END_REF], to estimate the resistance of a few-centimeter electric arc in the first phases of the discharge. Also in this case we get that the resistance transitory phase lasts a few microseconds. On the other hand the resistivity-temperature curves found in the literature can provide a model more physically sound. In this case we expect to get some more refined temperature profiles, in particular in the arc region. However, this model is very sensitive to the temperature variation that in turn depends on the net radiation coefficient and on the thermodynamic model.

In the rest of the paper we refer to the fixed resistive model as the model that uses the approximation (3) where the calculated resistivity depends on the assumed value of r 0 but it does not depend on the thermodynamic variables. On the contrary, when we use the η = η(T, P ) relation, we refer to the variable resistivity model.

Here we have chosen the fixed resistivity model, however some extensive comparisons with the variable resistivity model are included in section 5.

The current shape (2) implies that the arc geometry is imposed a priori. This is a major simplification since the modeling of the geometrical evolution of the arc could be problematic. In particular, it would require the solution of a Laplace equation for the electric potential where the diffusion coefficient is represented by the conductivity coefficient 1/η. In our case the coefficient 1/η can vary by about eighteen orders of magnitude, see [START_REF] Ph | One-dimensional hydro-kinetic modelling of the decaying arc in air-pa66-copper mixtures: 1. chemical kinetics, thermodynamics, transport and radiative properties[END_REF][START_REF] Haidar | Non-equilibrium modelling of transferred arcs[END_REF]. The air resistivity is mainly dependent on temperature and its variation from the hot arc zone to the cold boundary is the cause of the resistivity changes inside the domain.

We have tried to solve the Laplace problem with the finite element method, see [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF]. This technique reduces the differential problem to a linear system like Ax = b where A is the stiffness matrix, x is the vector of the unknowns and b is a known vector.

The linear system has been solved using the iterative and direct solvers developed in the Trilinos package, see [START_REF] Heroux | An overview on the Trilinos project[END_REF]. In our numerical experiments we have experienced the numerical breakdown of the linear solvers when roughly more than 100000 degrees of freedom are used. This is somehow expected, in fact the conditioning number of the stiffness matrix A (i.e. a number that measures the stiffness of the problem) is bounded by

C(sup(η)/ inf(η)H -2 ) ( 5
)
where C is a constant [START_REF] Ern | Theory and practice of finite elements[END_REF]. This means that the matrix conditioning gets worse quickly as the mesh is refined and as the matrix size increases. This effect is greatly amplified by the variations of the conductivity coefficients and can make the solution of the linear system very difficult both using direct [START_REF] Davies | Direct methods for sparse linear siystems[END_REF] and iterative solvers [START_REF] Elman | Finite elements and fast iterative solvers[END_REF]. Actually the conditioning number affects the iterations needed to reach convergence when an iterative scheme is used, and also affects the precision of the direct solvers. In fact the solution of the linear system is affected by the discretized arithmetics used by the computers and it can be shown, see [START_REF] Quarteroni | Numerical mathematics[END_REF], that δx / x ≤ (2γκ(A))/(1 -γκ(A)) where δx is the variation between the solution of the linear system solved with exact and inexact arithmetics. Moreover γ is the roundoff and, using the double precision standard, is approximately 10 -16 , and κ(A) is the conditioning number of the matrix A. In other terms, if the conditioning number is big, i.e. above 10 16 , the solution of the linear system can be unreliable or even unstable.

In spite of this, many solutions have been obtained in axis-symmetric geometries, see for instance [START_REF] Zhu | A unified theory of free burning arcs, cathode sheaths and cathodes[END_REF][START_REF] Haidar | Non-equilibrium modelling of transferred arcs[END_REF][START_REF] Tanaka | Predictions of weld pool profiles using plasma physics[END_REF][START_REF] Lowke | Theory of free-burning arc coloumns including the influence of the cathode[END_REF][START_REF] Ekici | Thermal and flow fields modeling of fast spark discharges in air[END_REF][START_REF] Kovitya | Two-dimensional analysis of free burning arcs in argon[END_REF], or in three dimensional geometries, see [START_REF] Freton | Comparison between a two-and a three-dimensional ar plasma configuration[END_REF][START_REF] Speckhofer | Experimental and theoretical investigation of highpressure arcs-part 2: the magnetically deflected arc (three-dimensional modeling)[END_REF][START_REF] Yang | Numerical analysis of the influence of splitter-plate erosion on an air arc in the quenching chamber of a low-voltage circuit breaker[END_REF][START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part ii. the threedimensional modelinfluence of external forces on the arc column[END_REF][START_REF] Xu | Three-dimensional modeling of the plasma arc in arc welding[END_REF][START_REF] Blais | Three-dimensional numerical modelling of a megnetically deflected dc transferred arc in argon[END_REF]. Most of the algorithms are based on the SIMPLEC approach described in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] and then used in [START_REF] Zhu | A unified theory of free burning arcs, cathode sheaths and cathodes[END_REF][START_REF] Haidar | Non-equilibrium modelling of transferred arcs[END_REF]. In other cases the arc is simulated using the PHOENICS code, such as in [START_REF] Tanaka | Predictions of weld pool profiles using plasma physics[END_REF][START_REF] Goodarzi | The effect of the cathode tip angle on the gtaw arc and weld pool: I. mathematical model of the arc[END_REF], or FLUENT [START_REF] Reynard | Modelisation tridimensionell de l'amorcage de l'arc electrique dans un disjoncteur basse-tension[END_REF], while in other cases some dedicated codes are developed [START_REF] Trelles | Finite element modeling of flow instabilities in arc plasma torches[END_REF]. In almost all the cases a finite volume scheme on a structured grid is used however these schemes experience a similar behavior of the conditioning number with respect to the finite element method. In one dimension the finite element and finite volume methods are even the same [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF][START_REF] Quarteroni | Numerical mathematics[END_REF]. Moreover in almost all the three-dimensional cases the maximum number of degrees of freedom (dof) used is under 100000: in [START_REF] Speckhofer | Experimental and theoretical investigation of highpressure arcs-part 2: the magnetically deflected arc (three-dimensional modeling)[END_REF] 60025 dof are used, in [START_REF] Xu | Three-dimensional modeling of the plasma arc in arc welding[END_REF] 56000 dof are used and finally in [START_REF] Freton | Comparison between a two-and a three-dimensional ar plasma configuration[END_REF] 48216 dof are used. In [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part ii. the threedimensional modelinfluence of external forces on the arc column[END_REF] 162729 dof are used, which is higher than the 100000 threshold. We stress that the conditioning number depends on many factors: in (5) the constant depends on the discretization methods. Therefore different methods and even different grids may have different constants. For instance it is well known that, using un-structured grids (in our case tetrahedral elements are used), the mesh quality affects the conditioning number, see [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]. In other terms the constant in (5) depends on many factors that may contribute to keep lower values of the conditioning number. However we expect that any mesh refinement affects quadratically the conditioning of A Some dedicated techniques, mainly based on the domain decomposition theory, have been developed to tackle the conditioning problem, see [START_REF] Toselli | A feti domain decomposition method for edge element approximations in two dimensions with discontinuous coefficients[END_REF], however they can be applied only where the coefficient (in this case the electric conductivity coefficient) has isolated jumps aligned with the elements faces. Unfortunately this is not applicable to our case since the variations of η cannot be determined a priori and are, in general, not aligned with the grid.

The Laplace problem with a strong coefficient variability provides a tough mathematical problem. We have found that, if the grid is not very refined, some results can be obtained and there are many cases in the literature, see for instance [START_REF] Yang | Numerical analysis of the influence of splitter-plate erosion on an air arc in the quenching chamber of a low-voltage circuit breaker[END_REF], where the results have been validated against experimental data. However we have found that the algorithm can become unstable, particulary the arc inception phase, where the arc expands, requires a very small time interval. The inception phase is quite important in this case because the first and most important pressure shock wave is generated in the first tens of microseconds. Correlated difficulties have been experienced in [START_REF] Xu | Three-dimensional modeling of the plasma arc in arc welding[END_REF] where a small relaxation factor is used to reach the convergence of the SIMPLE algorithm or in [START_REF] Yang | Numerical analysis of the influence of splitter-plate erosion on an air arc in the quenching chamber of a low-voltage circuit breaker[END_REF] where an implicit solver is needed to get a relatively large time step of 2µs.

When the coefficient variability is so strong, the precision obtained when solving the linear system is better when H is large, on the contrary, the precision associated with the discretization of the partial differential equation (using finite elements of finite volume methods) gets better when H gets smaller. Due to this and since the simulation of the arc inception phase is sensitive to the time stepping, we preferred to use a more conservative approach and to impose a fixed arc radius. A dedicated study to this kind of mathematical problem involving a high coefficient variability has been postponed to allow further study. Moreover we think that our simplified model provides a quick and computationally sound manner to determine the effects of pressure. The resistivity physics is somehow different and more complicated however, as we will see in the result section 5, our model can represent the main factors that are needed to get a reliable pressure estimation.

Thermodynamic model

For the thermodynamic model we have chosen a reacting mixture in chemical equilibrium. We have used the 11-species model described in [START_REF] Mottura | An evaluation of roe's scheme generalizations for equilibrium real gas flows[END_REF] were O 2 , N 2 , O, N , NO, NO+, O 2 +, N 2 +, O+, N +, eare considered. The thermodynamical equations to be added to system [START_REF] Reinmann | Temporal investigation of a fast spark discharge in chemically inert gases[END_REF] are:

e = ns i=1 y i cv i (T )T + ns i=1 y i ∆f 0 i , P = ns i=1 ρ i R u w i T, (6)
where e = E t /ρ -(1/2)| u| 2 is the internal energy for unit mass, y i , with i = 1, . . . , n s , is the mass fraction of the i-th species, cv i (T ) is the heat specific coefficient at constant volume, n s = 11 is the number of the species, ∆f 0 i is the formation heat, R u is the universal gas constant, w i is the molar weight and T is the temperature. In our case ρ and e have been considered the independent variables subsequently, using the minimization of the Gibbs free energy, we have computed all the other thermodynamic variables. To be more precise, this problem is solved for several values of the couple ρ, e and the results are stored in a database. The fluid dynamic code only interpolates, from the database, the desired values reducing significantly the computational burden. We have also adopted the local thermal equilibrium (LTE) assumption since at the arc temperature, above 20000K, the deviation from the thermodynamical equilibrium is quite small, see [START_REF] Haidar | Non-equilibrium modelling of transferred arcs[END_REF], and a one-temperature field can be considered. In Figures 4,5 6 we have displayed some data provided by our database. As expected, see Figure 4, the temperature depends mainly on the internal energy nevertheless at very low densities the temperature becomes also density-dependent and, moreover, this dependence grows with the temperature. A similar behavior is outlined in Figures 5 and6 where the specific heat and the formation one are displayed. These latter variables are highly dependent on the dissociation reactions and show strong variations when a specific reaction takes place.

One dimensional model of the tubes

For the tubes we have used a simplified one-dimensional Euler equations employing a polytropic ideal gas:

                       ∂ρ ∂t + ∂m ∂x = 0 in [0, L] × (0, ϑ f ] ∂m ∂t + ∂ ∂x m 2 ρ + P = -sign(m) 1 2 λ D m 2 ρ in [0, L] × (0, ϑ f ] ∂E t ∂t + ∂ ∂x (m (E t + P )) = 1 2 λ D |m| 3 ρ 2 in [0, L] × (0, ϑ f ] m = 0 on x = 0 × (0, ϑ f ] P = P b , if u > 0 o n x = 0 × (0, ϑ f ] P = P b , ρ = ρ b , if u ≤ 0 o nx = 0 × (0, ϑ f ], (7)
where λ is the friction coefficient of the tube, D = 1cm is the tube diameter, L = 25cm is the tube length and P b , ρ g are the inflow boundary conditions:

P b = P g -1 2 κ m 2 ρ if m ≤ 0 P g if m ≥ 0 , ρ b = ρ g , (8)
where P g and ρ g are the pressure and density solutions of (1) on the panel at the conjunctions of the tubes and κ is the inlet pressure loss coefficient. System ( 7) is a one dimensional Euler problem with a source term that models the pressure loss by viscous friction in a simplified manner. The coefficient λ can be estimated using the Moody's table [START_REF] Lf | Friction factors for pipe flows[END_REF]. The Reynolds number in the tube ranges, during the first phases of the discharge, in the 1.2 • 10 7 -6 • 10 6 interval. Since the relative roughness is approximately 7 • 10 -4 we are in the complete turbulence region with an associated friction coefficient λ of approximately 0.018. This means that the viscous dissipation, due to the large diameter of the tube, is only a minor effect. Finally, the boundary conditions ( 8) model the pressure loss at the tube inlet and the coefficient κ can be estimated using some experimental data [START_REF] Batchelor | An introduction to fluid dynamics[END_REF].

The numerical scheme

Problem ( 1) can be divided in some more simple problems using the time splitting technique, see for instance [START_REF] Dehghan | Time-splitting procedures for the solution of the two-dimensional transport equation[END_REF]. We use a first order splitting method and we solve sequentially:

(i) a Laplace problem for the magnetic potential

-∇ 2 A i = µ 0 J i , i = a, p. (9)
The magnetic field B i = ∇ ∧ A i is then obtained through a post processing procedure;

(ii) a homogeneous Euler problem

       ∂ρ ∂t + ∇ • m = 0 ∂ m ∂t + ∇ • m m ρ + P Ī = 0 ∂E t ∂t + ∇ • ( m (E t + P )) = 0; (10)
(iii) an ordinary differential equation

∂E t ∂t = Π a -Π I ; ( 1 1 )
(iv) and a system of ordinary differential equations

   ∂ m ∂t = J a ∧ B a ∂E t ∂t = m ρ • J ∧ B a . (12)
This splitting technique separates the fluid dynamic part from the electric one and from the thermal forcing terms. This, as we will see below, has many advantages. The first step is a simple Laplace-type differential equation that can be solved with standard finite element techniques, see, for instance, [START_REF] Ern | Theory and practice of finite elements[END_REF]. The second step can be treated with finite volume techniques, in particular we have used the one described in [START_REF] Selmin | Unified construction of finite element and finite volume discretizatins for compressible flows[END_REF][START_REF] Guradone | Finite element/volume solution to axisymmetric conservation laws[END_REF] and we have used a Roe-type Riemann solver adapted to reacting gas mixtures, for details see [START_REF] Mottura | An evaluation of roe's scheme generalizations for equilibrium real gas flows[END_REF]. Problem [START_REF] Graneau | The cause of thunder[END_REF] can be formulated as a standard ordinary differential equation in fact, since in the second step ρ and m are constant, we get

ρ de(T ) dt = Π a (T ) -Π I (T )
and, using a few mathematical manipulations, we obtain an ordinary differential equation in the canonical form:

dT dt = ρ ∂e ∂T -1 (Π a (T ) -Π I (T )) .
This equation can be treated by standard ODE solvers such as the implicit Euler one. The choice of an implicit solver is mandatory since the variability of the net radiation coefficient and of the arc power term with temperature is very high.

For the fourth step we have used the following first order scheme:

         m n+1 = m n + ∆t n J n a ∧ B n a m n+1/2 = m n + ∆t n 2 J n a ∧ B n a E t,n+1 = E t,n + ∆t n m n+1/2 ρ J n a ∧ B n a . (13)
This splitting scheme has many advantages, first of all it guarantees the positivity of the internal energy e. In fact, from the internal energy definition and from the third equation of ( 13), we have

e n+1 = E t,n+1 ρ - 1 2 | u n+1 | 2 = E t,n + ∆t n m n+1/2 ρ J n a ∧ B n a - 1 2 | u n+1 | 2 .( 14
)
Then, using the first and the third of ( 13), we get

1 2 | u n+1 | 2 = 1 2 | u n | 2 + ∆t n ρ J n a ∧ B n a • m n + ∆t n 2 J n a ∧ B n a = 1 2 | u n | 2 + ∆t n m n+1/2 ρ J n a ∧ B n a
and, combining it with ( 14), we obtain

e n+1 = E t,n - 1 2 | u n | 2 = e n .
In other terms, the scheme ( 12) conserves the internal energy. The positiveness of the entire splitting scheme, i.e. the positiveness of the internal energy, depends on the positive treatment of all the splitting steps. In particular we have found that the most sensitive step is the solution of the Euler equation. It is known, see [START_REF] Einfeldt | On godunov-type methods near low densities[END_REF], that not all the approximated Riemann solvers guarantee the internal energy positivity, in particular the Roe approximation, described in [START_REF] Mottura | An evaluation of roe's scheme generalizations for equilibrium real gas flows[END_REF], may fail. A possible work around is to employ a positive approximated Riemann solver such as the HLLC described in [START_REF] Nishikawa | Very simple, carbuncle-free, boundary-layer-resolving, rotatedhybrid reimann solvers[END_REF] although this solver is very diffusive if compared to the Roe's scheme.

We have chosen the Roe's approximation since, for our particular geometry, it works well producing no un-physical solutions. We have also to stress that in more general cases, with more complicated geometries, the usage of positive solvers is mandatory. Another advantage of our splitting scheme is that the time-steps are constrained only by the cf l condition [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. In fact, the third step is treated with an implicit scheme and does not introduce a further time constraint. Finally, the fourth splitting part has a slow-varying time dependent forcing term and is stable also with long time steps. For the solution of the one dimensional model [START_REF]Aircraft Lightning Test Method[END_REF] we have used a similar first order technique that allows to separate the Euler equation from the viscous forcing terms.

The solution of the one dimensional Euler equations is obtained through a discontinuous Galerkin method, see [START_REF] Cockburn | Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws 2: general framework[END_REF], and the boundary conditions are enforced using the weak boundary imposition technique described in [START_REF] Villa | Convergence of weakly imposed boundary conditions: the one dimensional hyperbolic case[END_REF].

Results

Arc diameter estimation

Our model is dependent on the arc diameter. Some experimental and numerical results, see [START_REF] Plooster | Numerical model of the return stroke of the lightning discharge[END_REF][START_REF] Akram | The evolution of spark discharges in gases: 1. macroscopic models[END_REF], suggest that the arc is a few centimeter wide. To verify this, we have acquired some images at optical wavelengths using a streak camera, see Figure 7(a). There we have aligned the streak direction with the electrode direction: in this way we obtain the diameter of the arc versus time. Our data confirm the estimations found in the literature and so we have assumed an arc diameter of 4.0cm.

We have also performed a sensitivity analysis to study the effects of the radius diameter on the pressure profile. In Figure 7(b) we have depicted the pressure evolution at 5cm and 0cm from the center of the panel. The differences between the pressure evolutions are quickly reduced, in fact just at 5cm from the arc, the three curves corresponding to 1.5cm, 2cm and 2.5cm radii are almost undistinguishable. At the impact point the three curves show a more heightened difference. The spread in the results is quickly reduced after the first current peaks. The plot representing the evolution for the 1.5cm-radius case shows the greatest difference with respect to the 2cm and 2.5cm cases. However our streak images, see Figure 7(a), seem to confirm that in the first phases the arc radius equals or exceeds the 2cm threshold. In this case we can consider only the 2cm and 2.5cm cases and get a peak pressure variation from 7e 6 P a to 5e 6 P a which spans most of the previous estimates provided in the literature, see, for instance, [START_REF] Plooster | Shock waves from line sources. numerical solutions and experimental meausurements[END_REF][START_REF] Plooster | Numerical model of the return stroke of the lightning discharge[END_REF][START_REF] Plooster | Numerical simulation of spark discharges in air[END_REF].

Pressure comparison at gauges

We have simulated the discharge sequence using our three dimensional fixed resistance models and then we have used the computed data for imposing the boundary conditions of the one dimensional model of the tubes. We have obtained the expected pressure at the gauges and we have compared it with the experimental data shown in Figure 8 for three different positions. In figure 8(a) the data of the first hole, set at a distance of 5cm from the center of the panel, are displayed. Four pressure measures are depicted (meas1 through meas4) as in the nearest position some measure dispersion is present. Even in a laboratory environment, using the same charging potential for the capacitor, the electric arc effects are slightly different. However we can show that our numerical results can predict accurately the pressure peak.

In Figures 8(b) and 8(c) the same data are shown for the second and third hole and again we can show a good agreement between the computed and the measured data. The small oscillation of the measures just before the arrival of the shock wave is due to the discharge induced currents. However we stress that at 4ms, when the shock wave hits the pressure gauge placed in the hole A, the arc current is low.

As regards the tubes we have checked that the Mach number on the inlet remains subsonic. Otherwise a normal shock wave is formed and we have to use some relations more complicated than (8).

Variable resistivity model

The implementation of a variable resistivity model, even if the arc diameter is fixed, is more difficult than a mere substitution of the fixed resistive coefficient with a variable one that depends on the temperature and on the pressure. In fact, if we impose at atmospheric conditions even a relatively weak current, we get that an enormous amount of power is injected in air. As a matter of fact the air resistivity at ambient temperature is very high therefore if we impose the atmospheric conditions as the initial condition the Joule power would be very high. But this is un-physical, in fact the arc is anticipated by the streamer-leader phase in which a thin hot and conductive channel is created. A complete study of this very fast phase is beyond the scope of this work, more details can be found in [START_REF] Morrow | Streamer propagation in air[END_REF]. From our viewpoint it is only necessary to consider an initial condition in which the hot conductive channel has already been created. In our case we consider a pre-heating of 1000K in the arc region. As it is pointed out in [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part ii. the threedimensional modelinfluence of external forces on the arc column[END_REF] the initial condition has only a minor impact on the rest of the simulation. Some results for the pressure field are depicted in Figure 9. In Figure 9(a) we have depicted the pressure evolution at the gauge-A computed considering three arc radii.

The data are plotted together with the experimental data. As we can see, the variable resistivity model is much more dependent on the arc radius. For a comparison with the fixed resistivity model we recall that in that case the total power injected in the arc is independent from the arc radius. In fact the total power is (πr 2 0 /l 0 )R a | J| 2 (πr 2 0 l 0 ) where (πr 2 0 l 0 ) is the discharge volume. Moreover | J|πr 2 0 = I then the total power equals R a I 2 . On the contrary as we have already seen, in the variable resistivity case, the total power depends on the arc radius. In fact in that case the total power is and this depends on the arc radius.

A good agreement with the experimental data is found for an arc diameter ranging from 4cm to 5cm, however this value depends on many factors such as the radiated power, the thermodynamic model and the resistivity curve itself. In spite of all the uncertain parameters we got an arc radius which is at least comparable with the one used in the fixed arc case.

In Figure 9(b) we show the pressure evolution at the center of the arc channel. As the channel is pre-heated the initial pressure is far from the atmospheric one, in particular for r = 2.0cm and r = 2.5cm (that are the values that produce some pressure predictions close to the experimental values) we have an initial pressure close to 30e 5 P a. The variable resistivity model produces lower estimates for the pressure in the center of the arc channel yet also the fixed resistivity model produces a similar spread if different radii are considered. The lower pressure registered in our case is due to the fact that the pre-heated air begins to escape from the channel just from t = 0. In the fixed resistivity case the air is heated in a few tens of microseconds.

In Figure 10 we compare the temperature results for the fixed and variable resistive cases. The temperature field is the most sensible with respect to the model uncertain parameters. In fact in Figure 10 we outline a clear difference between the variable and the fixed resistance models. We have depicted the temperature evolution corresponding to the first two power peaks, i.e the first positive current peak and the subsequent first negative current peak. In both cases the temperature rises till the first negative current peak is reached, then the temperature remains at quite high values. This is due to the fact that the second power peak releases a power that is not much lower than the first one, however the density is much lower. After the second peak the temperature rises at values that are so high that almost all the power injected is radiated. This is a saturation mechanism that limits a further increase in temperature. Moreover a decrease in the temperature leads to a decrease in the radiation output and the temperature loss diminishes. In other words radiation effect tends to stabilize the temperature field between the power peaks. This is emphasized in the variable resistivity model where the power due to the Joule effect is also a function of the temperature: the higher the temperature the lower the Joule power. Since in our model there is no thermal coupling with the panel, in Figure 10 no thermal boundary layer is displayed.

Three dimensional field

We now give a few comments on the results obtained with our three dimensional MHD code. In particular we are interested in the mechanism of the shock wave formation and in the mechanical effects that the pressure field can induce on the panel. In Figure panel where the current is injected. In Figure 11(a) some positions very close to the arc are considered and the effects of the pulsed current are clearly visible. The pressure peak, at the center of the channel, is above 7 • 10 6 P a and this is in good agreement with the data reported in the literature, see for instance [START_REF] Plooster | Numerical model of the return stroke of the lightning discharge[END_REF][START_REF] Plooster | Numerical simulation of spark discharges in air[END_REF][START_REF] Akram | The evolution of spark discharges in gases: 1. macroscopic models[END_REF]. In Figure 11(b) the far pressure field is depicted: the oscillating behavior is quickly smoothed out. Also the pressure field quickly decays as the distance from the impact point increases as we can see in Figure 12, where we have produced a three dimensional pressure plot. Several iso-surfaces are shown for two different time steps: t = 1.4 • 10 -5 s in Figures 12(a) and 13(a) and t = 1.82 • 10 -4 s in Figures 12(b) and 13(b). As we have already stressed, the pressure shock wave forms almost immediately and propagates through the air but, just a few hundreds of milliseconds later, it is strongly weakened. In fact on the plate boundary a little overpressure is measured, as can also be seen in Figure 11. On the contrary, the pressure remains very high, exceeding 1 • 10 6 P a, in the center; this behavior is quite important for the lightning protection aspects as most of the mechanical damages are concentrated in a very small area corresponding to the impact point. Most of the pressure wave is created in correspondence of the first peak and no major contribution is given by the following peaks; this can be easily noted from the computed density fields. In Figure 14 channel at very high speed. After a few hundreds of milliseconds it reaches a cyclic equilibrium moving away and getting closer to the channel with a phase delay with respect to the current pulses due to the air inertia. This cycle is clearly depicted in Figure 16 (see also Table 1), in particular in Figure 16(a), even if the current reaches its local maximum, the air is still entering the arc channel. The current then decreases in Figure 16(b) and reaches a low value in 16(c). However the maximum speed is reached when the current is near its minimum value. Then the speed field reverses again in Figure 16(d).

Conclusions

In this work we have validated a MHD code in some extreme conditions. The multi-scale validation technique is innovative and has produced good results. The code has been used to investigate the pressure field near the impact point in order to give an estimate of the dynamic loads that a high power arc discharge can produce. Our results show that the pressure in the arc region remains high during the discharge phase while the shock wave strength quickly vanishes. This could produce some very localized damages in aluminium alloys but even more severe damages in more rigid materials such as the composite panels. However we stress that the interaction between the arc channel and a composite panel is much more complicated and a dedicated modeling is required in order to get some reliable data.

In our work the discharge mechanism is studied highlighting the importance of the time to peak. This aspect is important for the design of lightning test devices. Most of the strength of the shock wave is generated by the first peak of current, while the pressure field close to the arc is more influenced by the current waveform. In fact in [START_REF] Aleksandrov | Effect of continuous current during pauses between successive strokes on the decay of the lightning channel[END_REF], where a non-oscillating waveform is used, the pressure evolution in the center of the arc channel is quite different whereas a similar behavior is observed elsewhere.

In spite of the uncertain parameters our work gives a very clear picture of the pressure field at a distance of a few centimeters from the discharge point. The maximum spread of the computed results is in the center of the arc, however, even there we give some reliable upper and lower bounds that are generated by and extensive sensitive analysis. Unfortunately, it is almost impossible to get some reliable measures just inside the arc channel, therefore the validation procedure has to rely on indirect measures. For the arc radius we have also used some streak images to back our hypotheses.

In future developments some more sophisticated numerical techniques will be developed whose final goal will be to implement a simulation tool that can predict also the density current function.

The code we have developed is based on the open-source spark2 platform http://gitorious.org/spark2/spark2.

Figure 1 .

 1 Figure 1. (a) Electric circuit, (b) plate-tube junction and the position of the three holes.

Figure 2 .

 2 Figure 2. The measured current and a fitted version.

Figure 3 .

 3 Figure 3. The coupled plate-air geometry and the boundary subdivisions.

  ,

Figure 4 .

 4 Figure 4. The contour plots of the temperature plotted against the internal energy and the density.

Figure 5 .

 5 Figure 5. The contour plots of the specific heat plotted against the internal energy and the density.

Figure 6 .

 6 Figure 6. The contour plots of the formation heat plotted against the internal energy and the density.

Figure 7 .

 7 Figure 7. A streak picture of the arc (a), the first three current peaks are visible. The pressure evolution (b) at R = 5cm and R = 0cm from the center of the panel with different arc channel radii: r15 is the r 0 = 1.5cm case while r20 and r25 are the r 0 = 2.0cm and r 0 = 2.5cm cases respectively.

Figure 8 .

 8 Figure 8. Pressure at gauges measured at holes A (a), B (b) and C (c) (5cm, 10cm and 15cm distant from the center of the plate). The meas1, meas2, meas3, meas4 labels stands for the pressure plots obtained in four experimental tests while mun stand for the numerically-predicted pressure plot. In Figures (b) and (c) we have considered only two measurements.

Figure 9 .

 9 Figure 9. In Figure (a): the pressure evolution at the gauge A, comparison of the experimental data (four measurements: meas1 though meas4) with the data computed with the variable resistive model: the data obtained imposing r 0 = 20mm (num20) and r 0 = 25mm (num25) are displayed. In Figure (b): the pressure evolution in the center, comparison of the variable resistive model, with three arc diameters (r 0 = 20mm (var20), r 0 = 25mm (var25) and r 0 = 30mm (var30)), with the fixed resistive model ref .

Figure 10 .

 10 Figure 10. Temperature comparison for the fixed resistivity model (a-b) and the variable resistivity model (c-d). The first (a-c), 1.05e -5 s, and second, 4.84e -5 s, power peaks are depicted (b-d).

Figure 11 .

 11 Figure 11. The pressure evolution versus time for different positions near the impact point (a), and far from the impact point (b).

Figure 12 .Figure 13 .

 1213 Figure 12. Three dimensional pressure fields at t = 1.4•10 -5 s (a) and t = 1.82•10 -4 s (b). A cut of the iso-surfaces on the x-z plane at y = 25cm is displayed. The pressure is measured in P a.

Figure 14 .

 14 Figure 14. Three dimensional density fields at t = 4.2 • 10 -5 s (a) and t = 1.54 • 10 -4 s (b). A cut of the iso-surfaces on the x -z plane at y = 25cm is displayed. The density is measured in kg/m 3

Figure 15 . 3 Figure

 153 Figure 15. Contour plots of a slice on the x -z plane at y = 25cm of the density fields at t = 4.2 • 10 -5 s (a) and t = 1.54 • 10 -4 s (b). The density is measured in kg/m 3

Figure 16 .

 16 Figure 16. Four snapshots depicting the flow field cycle during a current half period. A cross section in the x -z plane at y = 25cm is depicted. The velocity is measured in m/s

Table 1 .

 1 Time and current data for the four snapshots in Figure16.

Acknowledgements

The authors wish to thank L. Barbareschi for her valuable contribution and suggestions.

‡ This work has been supported by the MIUR/FIRB grant: Tecnologie e processi innovativi per gondole motore di velivoli commerciali di nuova generazione, DM25837 Confidential: not for distribution.


S Sfikas

E Amanatides

D Mataras

D Rapakoulias

Simulation of Cylindrical Electron Cyclotron Wave Resonance Argon Discharges

A fluid model of a cylindrical electron cyclotron wave resonance (ECWR) Argon discharge is presented. The results for a 1 mTorr Argon discharge were checked against the analytical theory, simulation and experimental data. The basic plasma properties as power dissipation, magnetic field, electric potential, electron density and temperature were very well reproduced by using pre-defined boundary conditions for the magnetic potential. The results of this model were further used as inputs for the simulation of plasma expansion into a diffusion region, allowing thus fast and complete modeling of a typical ECWR plasma reactor.

Introduction

Electron Cyclotron Wave Resonance (ECWR) plasma sources have recently found several applications in materials processing [1]. Plasma Oxidation and Plasma Enhanced Chemical Vapor Deposition (PECVD) of a-C:H, c-BN, nc-Si:H thin films are a few examples of the processes where ECWR has already been applied successfully [2][3][4][5][6]. The main advantage of these types of sources is that radio-frequency coupling into the inductively coupled plasma is enhanced in this case by the ECWR effect, leading to higher plasma densities compared to capacitive or purely inductive discharges [1].

The ECWR effect in plasma reactors is induced by the superposition of a small static magnetic field perpendicular to the axis of the radio-frequency excited coil and in parallel to the inductively excited electron cyclotron waves [7]. ECWR utilizes the propagation branch for generator frequencies ( rf ) well below the cyclotron frequency ( ce ) in the dispersion diagram of the "whistler" or right-handed mode; thus it has to be well distinguished from the Electron Cyclotron Resonance (ECR) effect. The superposition of a small static magnetic field produces a slight anisotropy on the plasma electrical properties. Resonant plasma excitation by electron cyclotron waves is possible, when an odd number of half of the wavelength of the propagating cyclotron wave fits the plasma thickness in the direction of the wave vector (i.e. the formation of a standing wave pattern) [1,7,8]. The latter ensures a more efficient plasma heating and higher plasma densities compared to conventional capacitive and inductive sources.

In an effort to understand and explain the resonance effect, the rather complicated anisotropic electrodynamics of these systems were studied through various approaches. Thus, theoretical studies [9], analytical solutions [8,10] and simulations [11] together with experimental findings [12,13] have been used to establish the theory behind the ECWR effect. However, self-consistent modeling of such complex gas discharges is still a challenging task.

It is very important for further consideration of these complex plasma sources in industrial applications [14].

We extend our recently reported fluid model of an ECWR discharge in a 3d infinite slab [15] to more realistic cylindrical plasma reactor geometries including also the plasma expansion region, which is very important for the simulation of discharges used in materials processing. Model results for 1 mTorr Ar discharges are presented and compared to other theoretical and experimental findings.

Model Description

In figure 1 the reactor geometry is presented, including the bounding dielectric walls, together with the static (B st ) and high-frequency (hf) alternating magnetic field (H z ) directed along the x and z-axis, respectively. The geometry of the chamber, a cylindrical glass vessel of diameter d cyl = 0.095 m and height h=0.095 m filled with 1 mTorr Argon is similar to those in Ref's [10][11][12][13]. A single sector coil carrying a current of frequency rf =2 . . 27.12 MHz is positioned immediately outside the bounding walls. We introduced a Cartesian coordinate system with the x-axis in the direction of the stationary external field B st and the z-axis coinciding with the axis of the infinitely long cylinder.

In order to reproduce the electromagnetic field in this cylindrical configuration, we have adopted the estimation for the right-hand polarized (RHP) wave propagating along the field B st as in Ref. [9]:

) .y/d cos( x) cos(K B B cyl RHP RHP0 RHP π ⋅ ⋅ ⋅ = (1a)
where, RHP K and B RHP0 are the wavenumber and the peak value of the RHP wave. The value of B st in resonance can be calculated by the formula proposed in Ref. [10]. The peak value (B RHP0 ) of the RHP wave, was calculated to be nearly 64.16 % of the value of the alternating magnetic field B RF at the walls, using back-substitution of known values (K RHP and d cyl ) in eq. (9a) of ref [8]. In addition, the value of the alternating magnetic field B RF at the walls was set to 0.32 mTesla (3.2 Gauss), resulting from a of 24.2 Ampere (peak value) coil current.

) / ( ) ( ) ( ) ( ) ( ) ( 2 2 || * 2 2 * 2 || 3 * * 2 2 0 2 2 ω π λ β β λ ε π ⋅ ⋅ ⋅ ⋅ + - ⋅ - ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ = - e e i
For the magnetic vector potential, we have used a specific set of boundary conditions, under the following assumptions:

-The RHP component of the wave is responsible for electron heating.

-The z-component of the alternating magnetic field, induces a cross sectional electric field (i.e. field components in the xy-plane) which is of major importance compared to the axially induced field.

-Plasma conductivity in the direction of the stationary field ( xx σ ) is predominant:

Inductive power absorption occurs mainly in regions where the x-component of the induced electric field prevails.

-Modification of the effective collision frequency introduced in [15] can lead to the substitution of plasma tensor conductivity by a scalar.

The A x and A y components of the magnetic potential can be expressed as:

(y/r) .A | ) c r K sin(b | c -c A max 2 3 RHP x 2 1 x ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ = (2)
(x/r) .A | r) K sin(b | c c A max RHP y 2 1 y ⋅ ⋅ ⋅ ⋅ ⋅ = (3)
where the e j t time dependence has been eliminated and A max is the maximum value of the magnetic potential calculated to be 2.097 . The magnetic field can be also approximately reconstructed using first-order Bessel functions for the magnetic potential:

(y/r) A ) c r (K J c - A max 3 RHP 1 1 x ⋅ ⋅ ⋅ ⋅ ⋅ = (4) (x/r) A ) c / r (K J c A max 3 RHP 1 1 y ⋅ ⋅ ⋅ ⋅ = (5)
The spatial distribution of the high-frequency magnetic field B z , calculated from the difference y

A x A x y ∂ ∂ - ∂ ∂ in Eqs. (2, 3
) is presented in Fig. 2(a). The distribution has the form of the analytical expression given in Eq. (1a). The high-frequency magnetic field B z can be reproduced close to the x-axis using the analytical theory for B z , while in the direction perpendicular to B st it is close to the values predicted for the imaginary (B y ) component in ref. [8]. Moreover, close to the y-axis, the magnetic field resembles the simulation results of ref. [11] at the phase of vanishing coil current in the rf cycle. Nevertheless, Krimke et al [11] have introduced a discharge model with predefined plasma properties (N e , e-n , p ε t , p σ t ) , within a static magnetic field (B st ) of 2.55 mTesla, corresponding to res =2.63 which is quite different from the experimental value.

The modulus of the cross sectional induced electric field calculated as

2 2 2 2 y x RF y x RF A A E E E + ⋅ = + = ω
is shown in the map of Fig. 2(b). The coil surface charge irrotational field [11] will be shielded by the sheath and was not taken into account in the simulation since it is not contributing to electron heating. Two fields of dipole nature ( x E and y E components) are superposed in RF E so the quadrupolar character of the field that was revealed in [11] is still present-yet distorted.

The inductive power dissipation is calculated as:

2 ) Re( 2 1 RF p RF E P ⋅ ⋅ = σ (6)
Inductive power dissipation mainly occurs in regions of high plasma conductivity and field strength. Calculation of the components of conductivity tensor using the Appleton-Hartree theory [START_REF] Seshadri | Fundamentals of Plasma Physics[END_REF] [START_REF] Chen | [END_REF] to be =4 . 10 6 s -1 ).

The region of interest lies in the y-axis which is perpendicular to the static field B st . In order to incorporate this power absorption non-uniformity, we have inserted a function g(r) = (dE x /d r )/E x0 in eq. ( 4) of Ref. [15] for the effective collision frequency eff :

ν ν ω ⋅ ⋅ + = 2 ) ( 1 1 eRES e c eff N N r g (7)
where, g(r) is the normalized gradient of the x-component of the induced electric field E x , and it essentially represents the influence of the E x field in electron mobility: The x-component of the electric field tends to affect the electron motion in a direction parallel to the B st field; as a consequence, the electron mobility and plasma conductivity will be considerably increased.

The penetration distance for the E x field was set to =0.017 m, which corresponds to the normal skin depth value (non-locality parameter <<1 [18]). Modifications in eff are taking place only within this region of E x penetration. Moreover, the penetration depth for E x was assumed to vary in the interval [0;0.017] as =0.017 . |y/r|, because of the decreasing influence of the E x component as y 0.

For the case presented here, the domain was split into 1152 cells. The magnetic field was solved in the frequency domain; while for the electrostatic field and the charged species mass, momentum (for ions) and energy (for electrons) balance we have employed a time step of 1.8436 ns. Finally, for neutral species a time step of 184.36 ns was used. The system of equations for the particle, momentum and energy balances was obtained from moments of the Boltzmann transport equation and was coupled with Poisson (electrostatic field calculation) and Maxwell (hf electromagnetic field calculation) equations; finally, it was complemented by a set of boundary conditions concerning the species' densities as well as the electrostatic potential as in [15].

Results and discussion

In figure 3(a) we compare the results of our simulation for the electron density (spatial map with dark colored curves of constant density), with the experimental findings presented in fig.

7 of ref. [12] for N res =1 . 10 17 m -3 , =2.0, z=0: The variation of the ion current to the Langmuir probe between its mean and maximum value, is displayed by the superimposed light colored arcs for r =1.0, 2.0, 3.0, 4.0 cm, representing higher electron density regions. The alignment of the regions of maximum electron density in the simulation is in good agreement with experiment in the main volume of the discharge and the net result of the anisotropy in charged particle diffusion coefficient is closely reproduced. In ECWR sources, the electron density profile is a decreasing function of the radius r in the direction perpendicular to B st as resonance is approached, resulting in an increasing in-homogeneity of the extracted ion beam. This effect is demonstrated in figure 3(b): plasma density experimental data [12] collected in the center of the discharge along the diameter which is perpendicular to B st ( =0, 1.5, dotted lines) are compared to the simulation results ( res =2.13589, straight lines). The normalized electron density radial profile (N e (x), =2.13) along the static field is in fair agreement with the experimental data radial profile without the static field (N e (y), =0) as it should due to the unaffected electron motion along the lines of force. The radial profile of the simulated electron density (N e (y), =2.13) perpendicular to the static field is off the experimental values (N e (y), =1.5); nevertheless, due to the hindered electron motion perpendicular to the static magnetic field, the simulated profile for N e (y) is sharper than for N e (x), as it should.

The difference from the experimental findings can be attributed to the different plasma density symmetry properties at non-resonant values ( res > =1.5) where the alignment of maximum density regions is along the static magnetic field [12].

The simulation results for the inductive power dissipation along (P RF (x)) and perpendicular (P RF (y)) to the static field are presented in fig. 4(a): The radial profiles are entirely different due to the anisotropy of the conductivity in the plasma calculated as:

ω ν ω ε σ ⋅ + ⋅ = i eff p p 0 (8)
where p ω is the electron plasma frequency.

The map of P RF distribution is shown in fig. 4(b). The regions of maximum power dissipation are aligned along the y-axis; the curved skin layer resulted in a power absorption decreasing as r x ± → in a way comparable to the simulation results of Ref. [11]. A total power dissipation of 46 W was calculated from volume integration for the vessel of height h=0.095 m, assuming homogeneity in the cylinder's axial direction. This value is almost 60% of the 80 W of total dissipated power referred in experimental [12] and simulation [11] results; however in an experimental arrangement [13] a low level hf generator output power 50

≅ hf P
Watts was reported for the fundamental ( =2.0) resonance.

Furthermore, the model prediction of electron temperature T e is presented in fig. 5(a). The simulated profiles demonstrate a deep 4 eV minimum at the center of the discharge as a consequence of the hindered electron transport in the main volume. This value is close to experimental and simulation results reported in [14,19] for a (Ne ≈ 1 . 10 17 m -3 , p=1 mTorr, rf =2 . . 29 MHz, B st =20 Gauss, P RF =100 W) magnetized inductively coupled discharge configuration. The spatial anisotropy in energy transfer to electrons is reflected to the different maximum values and gradient of T e along (x) and transverse (y) to the static magnetic field B st . Although the power dissipation takes place in the direction perpendicular to the static magnetic field (fig. 4(b)), discharge heating seems to take place principally along the lines of force (T e (x) curve) where electron motion is undisturbed. The minimum in electron temperature at the center of the discharge can also be attributed to low energy electrons trapped in the electrostatic potential which exhibits spatial homogeneity and its' radial profiles are presented in fig. 5(b). It has a maximum value of max ≅ 36.5 Volts and it can be approximately reproduced by eq. ( 9) for the ambipolar potential to account for a nonuniform plasma density:

} 1} - d ]/ d ) {[(2 - 1 { ) ( 6 cyl cyl max + ⋅ ⋅ = x x , x ∈[-d cyl /2, d cyl /2] (9)
The coupling of the Poisson's relation to the fluid code was of great importance; otherwise the assumption of quasi-neutrality and ambipolarity led always to the underestimation of the electrostatic potential. The electrostatic potential distribution and range was calculated to be close to hybrid (fluid/kinetic) [20,21] simulation results (p=1 mTorr, Ne ≅ 10 17 m -3 ), of a nonmagnetized low pressure ICP discharge [START_REF] Kolobov | Simulation of Electron Kinetics in Gas Discharges[END_REF]. The spatial distribution and the range of the simulated electron temperature in fig. 5(a) are also comparable to the previous ref. [START_REF] Kolobov | Simulation of Electron Kinetics in Gas Discharges[END_REF]. Moreover, the model has been extended to include the plasma expansion into a diffusion region and towards the process chamber results for the 1 mTorr Ar ECWR plasma. The simulation geometry is similar with the one presented in ref. [START_REF] Vartolomei V Hannemann | Conf. on Phenomena in Ionized Gases[END_REF]. As an inlet of this diffusion regime, the bottom of the source (coil) was considered and the results of the previous analysis concerning the profiles along the x-axis of electron and ion density, electron temperature and potential were used as boundary conditions. Perfectly conducting wall at a distance of 25 cm from the bottom of the coil was considered in order to achieve a faster convergence. The results for the spatial distribution of electron density and potential from the center of the source towards the reactor wall are shown in figures 6 (a) and (b), respectively. In the same figures, a comparison of the data with the experimental results and the fitting of the spherical plasma expansion model [START_REF] Vartolomei V Hannemann | Conf. on Phenomena in Ionized Gases[END_REF][START_REF] Hippler | [END_REF] is included.

Even though the model nicely reproduces the rapid plasma degradation, electron density is overestimated especially for distances higher than 5 cm from the source (figure 6a) while a better agreement was found for the variation of potential (figure 6b) for distances higher than 10 cm. Consideration of the chamber and source walls as dielectrics seems to give slightly better results but it also leads to much higher computational time.

Further work is in progress with modified inlet conditions and more complicated gas phase and surface chemistry in order to validate and apply this fast converging code in more complicated chemical systems.

Conclusion

A fluid model of low pressure ECWR Argon discharges, based on a modified effective collision frequency and predefined boundary conditions for the magnetic potential was developed. The model predictions are in good agreement with other experimental measurements and simulations for Ar discharges, allowing for rapid estimation of basic plasma properties such as electron density, temperature and plasma potential. The simulation outputs can be used as inputs in a plasma diffusion model which are in reasonable agreement with experimental data and a spherical plasma expansion model [10].

The development of this fluid model is constantly under the perspective of its extension in order to simulate complex chemical gas phase and surface processes as deposition of amorphous and microcrystalline silicon from mixtures of silane in hydrogen.

FIGURE CAPTIONS

Figure 1: Geometrical model of the axially infinite plasma cylinder, with the axial component of the simulated alternating magnetic field H z and the superimposed static magnetic field B st .

Figure 2: Maps of (a) the simulated alternating magnetic field B z (normalized to the field value at the walls) and (b) the magnitude of the induced electric field, resulting from the specific boundary conditions used for the magnetic potential.

Figure 3: (a) Map of electron density N e as calculated and compared to the experimental curves [12] of higher density regions (light colored lines) at radius r = 1, 2, 3 and 4 cm and (b) Electron density profiles along (x) and transverse (y)direction to the static magnetic field B st as calculated in the present simulation and experimentally measured in [12]. [START_REF] Vartolomei V Hannemann | Conf. on Phenomena in Ionized Gases[END_REF][START_REF] Hippler | [END_REF].

  T e =80000 o K for the electron temperature. Setting the value of the desired electron density at resonance (N res ) to 1 . 10 17 m -3 , eq. (1b) results to res of 2.07 mTesla (20.7 Gauss) for the stationary magnetic field. Moreover, for K RHP= 51.88 m -1 , half a wavelength of the RHP wave fits into a slab of width d slab =

  10 -6 Weber/m[15]. The following values have been used for the dimensionless numerical constants b x , b y, c 1 , c 2 , c 3 in order to reproduce the magnetic field of Eq. (1a)

Figure 4 :

 4 Figure 4: (a) Inductive power dissipation PRF profile, along (x) and transverse (y) to the static magnetic field Bst.(b) Inductive power dissipation PRF spatial map.

Figure 5 :

 5 Figure 5: (a) Electron temperature T e and (b) electric potential profiles, along (x) and transverse (y) to the static magnetic field B st .

Figure 6 :

 6 Figure 6: (a) Electron density and (b) plasma potential as a function of the distance from the center of the coil towards as calculated in the present simulation and compared to the spherical expansion model and experimental measurements of[START_REF] Vartolomei V Hannemann | Conf. on Phenomena in Ionized Gases[END_REF][START_REF] Hippler | [END_REF].

Figure 1 :

 1 Figure 1: Geometrical model of the axially infinite plasma cylinder, with the axial component of the simulated alternating magnetic field H z and the superimposed static magnetic field B st .

Figure 2 ,Figure 2 :Figure 3 :

 223 Figure 2, S. Sfikas et al

Figure 4 :

 4 Figure 4, S. Sfikas et al

Figure 5 Figure 5 :

 55 Figure 5, S. Sfikas et al

Figure 6

 6 Figure 6, S. Sfikas et al

Figure 6 :

 6 Figure 6: (a) Electron density and (b) plasma potential as a function of the distance from the center of the coil towards as calculated in the present simulation and compared to the spherical expansion model and experimental measurements of[START_REF] Vartolomei V Hannemann | Conf. on Phenomena in Ionized Gases[END_REF][START_REF] Hippler | [END_REF].

  , for N res =1 . 10 17 m 3 , c = 363.95 MHz, rf =2 . . 27.12 MHz,

	leads to	σ	xx σ	yy	≅	2	.	235

p =17.8 GHz and the collision frequency as calculated in