An improved equation of state for air plasma simulations
AAn improved equation of state for air plasma simulations
An improved equation of state for air plasma simulations
F. Träuble, a) S. T. Millmore, b) and N. Nikiforakis Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom (Dated: 12 January 2021)
This work is concerned with the development of a novel, accurate equation of state for describing partially ionised airplasma in local thermodynamic equilibrium. One key application for this new equation of state is the simulation oflightning strike on aircraft. Due to the complexities of species ionisation and interaction, although phenomenologicalcurve fitting of thermodynamic properties is possible, these curves are intractable for practical numerical simulation.The large difference in size of the parameters (many orders of magnitude) and complexity of the equations meansthey are not straightforward to invert for conversion between thermodynamic variables. The approach of this paperis to take an accurate 19-species phenomenological model, and use this to generate a tabulated data set. Coupledwith a suitable interpolation procedure this offers an accurate and computationally efficient technique for simulatingpartially ionised air plasma. The equation of state is implemented within a multiphysics methodology which can solvefor two-way coupling between a plasma arc and an elastoplastic material substrate. The implementation is validatedagainst experimental results, both for a single material plasma, and an arc coupled to a substrate. It is demonstrated thataccurate, oscillation-free thermodynamic profiles can be obtained, with good results even close to material surfaces.
I. INTRODUCTION
Numerical simulation of air plasma interacting with solidsubstrates provides insight into the physical mechanismswhich lead to structural damage due to lightning strike, whichis a key safety concern for aircraft . The aluminium skin tra-ditionally used for aircraft rapidly dissipates the energy inputfrom the lightning strike, due to its high electrical and ther-mal conductivity. With the current trends towards lightweightcomposite skins, with lower conductivity, there is greater en-ergy input local to the lightning attachment point. To accu-rately understand the physical effects at this attachment, anaccurate description of the ionisation processes within theplasma arc is essential. This requires an equation of state(EoS) which can account for the complex thermodynamic andelectromagnetic effects arising from partially ionised mate-rial. The multiple atomic species comprising an air plasma areeach governed by different characteristic ionisation energies,governing their initial emergence and further ionisation. Thisleads to a non-linear relationship between intensive thermo-dynamic variables within the EoS for an air plasma. Compu-tation of these quantities, which on a continuum scale modelof a plasma are local physical properties, and are independentof system size, must still take into account transport proper-ties arising from molecular interactions between the differentspecies present.It is a key challenge in developing an EoS for air plasmato balance accuracy with computational efficiency. One ofthe simplest techniques is to use a surrogate ideal gas modelwith an adiabatic index appropriate for much of the ionisedregime, for example Ekici et al. use a an effective adiabaticexponent of γ = .
16. In a lightning arc simulation, for whichthere is the transition of unionised air to a plasma, an idealgas model has limited validity. An improved model was de-veloped by Plooster and has been applied in several plasma a) Also at Max Planck Institute for Intelligent Systems, Tübingen, Germany b) Electronic mail: [email protected] discharge models. It is a simplified equation of state for asurrogate diatomic gas, which takes into account dissociationeffects. However, due to the different dissociation energiesfor N and O , key components of an air plasma, and the factthat nitric oxide ( NO ) is not considered, this approach is lim-ited to lower temperature mixtures, below 9 , . Villa etal. developed an EoS which models 11 of the most commonspecies within an air plasma, allowing for single ionisation of N and O , as well as the dissociated atoms. Each species istreated as an ideal gas which are then combined as a gener-alised ideal gas mixture. They start from the work of Mot-tura et al. who assume a simplified heat capacity of the form C v ( T ) with T being the temperature, neglecting effects froma variation in pressure. The equation of state for each speciesassumes a linear relation between the specific heat contribu-tion and species’ dependent formation heats. To fully describea given set of thermodynamic data, a system of coupled lin-ear equations must be solved. In order to maintain computa-tional efficiency, Villa et al. tabulated the thermodynamic datafor given pairs of density and specific internal energy, and in-terpolate values between these points. This tabulated formallows for a complex description of the plasma composition,though, in common with other equations of state available inthe open literature, it does not consider plasma-specific ef-fects arising from electromagnetic interaction on microscopicscales.In this work, an improved equation of state is developedbased on an accurate air plasma model covering a wide rangeof variables from the theoretical model of D’Angola et al. .In section 2 the thermodynamic theory of air plasma, uponwhich the equation of state is based, is detailed. Addition-ally, a mathematical formulation for describing the governingequations of motion for an air plasma, which will be used tovalidate the improved EoS, are given, using the model devel-oped by Millmore and Nikiforakis . The generation of thenew EoS is detailed in section 3 and section 4 presents val-idation of the improved EoS within a state-of-the-art multi-physics multimaterial lightning code. It is demonstrated thatan accurate EoS is required to model the arc attachment point a r X i v : . [ phy s i c s . c o m p - ph ] J a n n improved equation of state for air plasma simulations 2to a substrate, and thus is essential in order to gain a deeperinsight into this highly complex and nonlinear phenomenon.Conclusions and further work are presented in section 5. II. THEORY AND MATHEMATICAL FORMULATIONA. Air plasma model
Atmospheric air is a mixture of many different particlespecies which, over the temperature and pressure range ex-perienced in a lightning strike, will interact in a complexmanner . The composition of dry air per volume is primar-ily made up of 78 . . . . .An appropriately chosen theoretical model of air plasmahas to cover the physics between the limits of unionised low-temperature air and highly ionised air plasma. The generationof the EoS presented here is based on an accurate theoret-ical model for a 19-species air plasma model developed byCapitelli, Colonna, D’Angola and others in . Thisair plasma model considers the following 19 air species: N , N + , N , N + , N ++ , N +++ , N ++++ O , O + , O − , O , O − , O + , O ++ , O +++ , O ++++ NO , NO + , e − Molecular species with more than two atoms are not consid-ered, as they only occur in low quantities, so their effect onthe application presented here is negligible, and the theory forexcited energy levels of these species is complex. The chem-ical processes among the different species are assumed to bein detailed balance at all times with the only exception beingradiation processes, which involve the emission and absorp-tion of photons. This assumption is valid as long as molecu-lar timescales are significantly smaller than fluid and gas flowtimescales; this is satisfied under the thermodynamic condi-tions of interest. In this case, a unique temperature and pres-sure can be assigned to a given composition. As such inten-sive variables are inhomogeneous in time and space for theapplications of interest in this paper, the concept of local ther-modynamic equilibrium (LTE) needs to be defined. LTE re-quires that the particles’ energies in a local neighbourhoodobey Maxwell-Boltzmann distributions, so that a local tem-perature and pressure can be defined. This is not always per-fectly justified, but work from Haidar indicates that it is suf-ficiently accurate to predict the expected behaviour in a light-ning attachment.
1. Thermodynamics
In dynamical equilibrium, every reaction is governed by alaw of mass action, taking the form K Pr ( T ) = N total ∏ s = P ν r , s s r = , , ..., N R (1)where K Pr ( T ) is a temperature dependent equilibrium constantof the r th reaction and P s is the partial pressure of the s th species. N total is the number of considered species and ν r , s the stoichiometric coefficients of the s th species in the r th re-action. The total pressure of a gas mixture is then given byDalton’s law P = N total ∑ s = P s + P DH (2)where the correction term P DH , derived from the Debye-Hückel theory, considers additional effects for ionised gasmixtures, with further details below. Every species is par-tially approximated by P s = n s k B T with n s the particle numberdensity of the s th species, k B the Boltzmann constant and T the temperature. The particles are modelled as dimensionlesshard spheres and virial corrections are neglected. This is jus-tified as they are usually only important for T < P ≈ , i.e. unrealistic conditions for thearc channel plasma and surrounding air.For given pressure and temperature, the system of nonlinearequations arising from (1) and (2) determines the exact com-position of all the species present in the air model. A furthercondition is given by the conservation of the predeterminedtotal mass of each atomic species, adding two additional con-straints to the 17 equations to close the system. It is challeng-ing to solve such a nonlinear system numerically, since equi-librium constants differ by many orders of magnitude. Thecomposition used for the numerical fits in is based on anaccurate hierarchical algorithm for solving the equilibrium ofa given reactive mixture .Once the composition { n s } s = ,..., N total is known, the ther-modynamic potentials are calculated using the respective par-tition function for each species . Due to the formation ofa plasma potential in ionised gases, and plasma in general,the thermodynamic properties are modified when comparedto an ideal gas. The theory and framework to consider sucheffects were developed in the Debye-Hückel theory . Thisintroduces a cut-off energy, making the sum concerning theinternal partition functions finite and further correction termsfor the pressure and potential . Finally, the heat capacitiesare given by C P = (cid:32) ∂ H ∂ T (cid:33) P (3) C V = (cid:32) ∂ U ∂ T (cid:33) V (4)with H being the enthalpy and U the internal energy.n improved equation of state for air plasma simulations 3
2. Transport properties
In the presence of gradients in the spatial distribution ofthermodynamic quantities, transport effects emerge which actto reduce these gradients. Each thermodynamic behaviour hasa different responsible transport mechanism. Ordinary diffu-sion is associated with a concentration difference and thermaldiffusion connected to a gradient in temperature. Viscosityis the mechanism describing the transfer of momentum due tovelocity gradients, while thermal conductivity contains effectsthat lead to transport in thermal energy such as thermal gradi-ents, chemical reactions or redistributions of internal degreesof freedom. Electrical conductivity is related to the transportof charged particles due to differences in the electric poten-tial. At energies relevant to air plasmas within this work, it isjustified to assume that the transport coefficients of electronsand heavy atoms or ions decouple from each other . Nu-merical fits regarding the transport coefficients in have beencalculated using third-order approximated coefficients derivedfrom the Chapman-Enskog theory . The Chapman-Enskogtheory provides a set of equations for calculating the trans-port properties of a multi-species gas mixture under the as-sumption that thermodynamic states can be described in LTE.Its starting point is the Boltzmann equation formulation andit assumes some justified microscopic models for the binarycollision term . The main approach of the theory representsa Chapman-Enskog expansion where the probability densityfunction is expanded perturbatively in a small parameter ε .This will result in a general set of Navier-Stokes equationsfor this system, including terms for transport coefficients. Adetailed derivation can be found in Capitelli et al. , and thisderivation provides expressions for characteristic transport co-efficients to an arbitrary degree of order ξ . In this work,the transport coefficients of interest are electrical and thermalconductivity . Corresponding collision integrals, as calcu-lated by D’Angola et al., take into account a variety of species-depended interaction approximations. B. Mathematical formulation of a plasma discharge
To model the dynamics of a plasma discharge for lightningstudies, a magnetohydrodynamic formulation is used. Thisassumes the plasma is modelled as a single fluid with a cor-responding equation of state that considers thermodynamicsof interacting species; this is appropriate for an air mixture inLTE. The governing equations are ∂∂ t ρ + ∇ · ( ρ u ) = ∂∂ t ( ρ u ) + ∇ · ( ρ u ⊗ u + P I ) = J × B (6) ∂∂ t E + ∇ · ( u ( E + P )) = u · ( J × B ) + ρ res | J | − S r (7) ∇ A = − µ J (8) with density, ρ , velocity u , total energy E = ρ e + ρ u , mag-netic potential, A , magnetic field B = ∇ × A , current den-sity, J , electrical resistivity ρ res = / σ , electrical conduc-tivity σ , permeability of free space µ , and a source termgoverning radiative losses, S r . As the system is underdeter-mined (9 unknowns in 8 equations), one further equation isneeded for closure. This is provided by an equation of statefor air plasma, providing the mathematical relation betweenspecific internal energy ( e ) , pressure ( P ) and mass density, inthe form e = e ( ρ , P ) . Within this framework, a cylindricallyaxisymmetric geometry is considered, with the arc connec-tion at the centre of the domain, which is sufficient to capturethe bulk behaviour of the arc-substrate interaction . The sys-tem, which is representative of a typical set up in lightning-protection laboratory experiments, is schematically shown inFigure 1. Current profiles representative of those used in labo- FIG. 1. Schematic showing the axisymmetric plasma dischargesetup. A plasma arc is generated at the electrode at the top of theimage, and connects to the substrate beneath it. A cut-through ofthe air/plasma region is shown, illustrating the arc; the plasma doesoccupy the full domain, to the edge of the substrate. Beneath thesubstrate is a subsequent air region (not shown), an approximate two-dimensional representation of this setup is shown in Figure 2. ratory lightning testing are included within the model throughthe current density. Two different models are studied: a sim-plified 1D model and a more realistic axisymmetric model.Within the 1D model, radial arc profiles are reproducedby assuming translational invariance along the z-coordinate.Consequently, the current density has a component only inthe z -direction and is approximated by a Gaussian profile assuggested in J ( r , t ) = − I ( t ) π r e − ( r / r ) e z (9)where r is the characteristic length scale of the plasma arcchannel radius and the total current I ( t ) is a given input func-tion. By assuming such a predetermined profile, the resultingmagnetic field is not coupled to the fluid evolution itself.Within the axisymmetric model, the geometry of an elec-trode is explicitly included and the current density is com-puted dynamically. Hence, this approach is more accurate as itn improved equation of state for air plasma simulations 4calculates the magnetic field in a plasma arc self-consistentlysince it is coupled with the evolution of the air plasma. It alsoallows for interaction with the boundaries to affect the currentprofile. The current density is connected to the electric field E or its corresponding potential φ through J = σ E = − σ ∇ φ (10)and due to the conservation of charge, it holds ∇ · J =
0. Asa consequence, elliptic PDEs for φ and A have to be solved,given by ∇ · (cid:0) σ ∇ φ (cid:1) = ∇ A = µ σ ∇ φ (12)which represent some form of generalised Poisson’s equa-tions on the domain. Figure 2 shows the domain and its ge- air plasma substrate (aircraft skin) air electrode z r
1 2 3 4 5 7 6 φ A r A z ∂φ∂ n = − σ I ( t ) π r A r = ∂ A z ∂ n = ∂φ∂ n = A r = ∂ A z ∂ n = ∂φ∂ n = ∂ A r ∂ n = A z = φ = ∂ A r ∂ n = A z = ∂φ∂ n = ∂ A r ∂ n = A z = ∂φ∂ n = A r = ∂ A z ∂ n = ∂φ∂ n = ∂ A r ∂ n = ∂ A z ∂ n = ometry for the 2D axisymmetric model. Applied boundaryconditions for the electromagnetic potentials are summarisedin the table beneath the diagram. If an aluminium panel isused as the aeronautical substrate, an electrical conductivity of σ = . × Sm − is used, i.e. higher than typical air plasmaconductivity within the arc channel. Alternatively, this sub-strate will be modelled as a carbon fibre composite material.A complete implementation of an anisotropic equation of state for carbon fibre composites is beyond the scope of this work.Instead, an approximation to CFRP is used, following , whichdefines the material as a isotropic material, which can be em-ployed to look at behaviour either perpendicular or parallelto the CFRP weave. This material differs from aluminium inelectrical conductivity ( σ = . × Sm − ) and mass den-sity ( ρ = / m ) .The model described above is implemented within a moreextensive simulation code. This includes explicit thermome-chanical coupling with various aircraft layer configurations ina full multimaterial model . This uses the Riemann GhostFluid Method (GFM) to accurately model the multimaterialboundary conditions, i.e. propagation and reflection of shockwaves across a multimaterial interface. The location of the in-terface is tracked by a level set function. The original GFMwas developed by Fedkiw et al. in , and a detailed introduc-tion to Riemann GFM is given in the work of Sambasivanet al. . Different materials are evolved independently fromeach other on their respective domains, and a unique equa-tion of state can be assigned to each material. In this code,the elastoplastic substrate and electrode are described usingthe Eulerian framework as presented by Schoch et al. andMichael et al. , based on the formulation of Godunov andRomenskii . Plasticity effects are incorporated following thework of Miller and Collela . Details are discussed in Mill-more and Nikiforakis , which validate the equations from theprevious section. III. METHODOLOGYA. Generating data for an improved equation of state
D’Angola et al. calculated a wide range of thermody-namic and transport quantities for given temperature and pres-sure (mixture composition, molar mass, enthalpy, heat ca-pacities, electrical and thermal conductivity), based on theair plasma model described in section II A. Therewith an im-proved equation of state for a 19-species air plasma model wasdeveloped. The equation of state is inserted into the plasmadischarge model from section II B by the help of a tabulardatabase, denoted by EOS19. In principle, EOS19 can be usedby any simulation code that attempts to model air plasma ap-plications in the future. Every data point in the database gen-erated contains the following quantities:• mass density ρ i • pressure P i • specific internal energy e i • temperature T i • speed of sound c i • adiabatic index γ i • molar fractions of each species χ i • electrical conductivity σ i n improved equation of state for air plasma simulations 5• thermal conductivity κ i The actual relation for the equation of state is thengiven through discrete tuples of the form ( ρ i , P i , e i ) with e ( ρ i , P i ) = e i . To generate this data the remaining quantitiesnot provided as analytical fits in need to be calculated. Themass density of air plasma ρ is calculated with the help of themean molar mass ¯ M via the formula ρ = PRT ¯ M (13)with R being the gas constant. Specific internal energy is de-rived from the specific enthalpy with the help of a Legendretransformation according to e = h − P ρ (14)Although negligible in this model, thermal conductivity fromanalytical fits in is included in the database, such that it canbe used to considering associated effects in the future. Finally,the adiabatic index and corresponding speed of sound, a , forair plasma as suggested in are calculated as γ = C P C V ρ P (cid:16) ∂ P ∂ ρ (cid:17) T (15) a = (cid:115) γ P ρ (16)with the heat capacities defined in equations (3) and (4).For typical choices of conserved variables, it is desirableto generate the thermodynamic EoS quantities for given pairsof density and pressure, instead of temperature and pressure.Therefore, the temperature has to be obtained inversely forgiven density and pressure in the form T ( ρ , P ) with the helpof the known relation in (13) ρ ( P i , T i ) = ρ i (17)where the tuple ( T i , P i , ρ i ) represents a thermodynamical state.Finding the solution for T i in (17) reduces to a root-findingproblem, and a simple bisection method is applied. It isnot important to have an efficient algorithm here since thetabulated data is generated only once. However, even witha more advanced method, this root-finding process is time-consuming, which is the reason for creating a tabulateddatabase instead of calculating parameters dynamically dur-ing simulation. From this data, the equation of state isinterpolated. The layout of EOS19 follows the SESAMEDatabase standard for equations of state as introduced atLos Alamos National Lab . It contains data points rangingfrom 0 . to 10kg/m in mass density and 0 . B. Numerical approach for solving the plasma discharge
For each computational time-step ∆ t , the system of equa-tions (5) - (8) is evolved in four steps, following the approachin .Step 1 The Poisson’s equation in (8) for the magnetic vectorpotential A is solved using the Thomas-Algorithm in1D and the finite element solver framework FEniCS in axisymmetry. Both algorithms are an appropriatechoice regarding accuracy and efficiency, however anysuitable solver could be used. In axisymmetry, themagnetic vector potential has the non-zero components A z ( r , z ) and A r ( r , z ) . The corresponding magnetic field,pointing in the poloidal direction, is then given by B φ = (cid:0) ∇ × A (cid:1) φ = ∂ A r ∂ z − ∂ A z ∂ r (18)with A r = ∂∂ t ρ + ∇ · ( ρ u ) = ∂∂ t ( ρ u ) + ∇ · ( ρ u ⊗ u + P I ) = ∂∂ t E + ∇ · ( u ( E + P )) = .Step 3 To incorporate Joule heating and radiation source termsin the energy balance equation (7), an ordinary differ-ential equation of the form dEdt = ρ res | J | − S r (22)has to be solved, with details regarding the implemen-tation of the radiation model given in .Step 4 To include Lorentz force effects, the ODEs given by ddt ( ρ u ) = J × B (23) dEdt = u · ( J × B ) (24)are solved by using a standard second-order ODE solverwith two half-time step updates that conserves the en-ergy and momentum . IV. RESULTS
In this section, data from the improved, newly generatedequation of state is shown together with its implementation inthe 1D and axisymmetric plasma models.n improved equation of state for air plasma simulations 6
EOS19 data
The specific internal energy generated by EOS19 for agiven mass density and pressure is shown in Fig. 3. These
FIG. 3. Heat map of the specific internal energy for EOS11 (top) andEOS19 (bottom) values are compared to the publicly available data of the equa-tion of state developed by Villa, denoted as EOS11 . Asthere are no classical phase transitions observed in air plas-mas, a physically accurate description requires smooth con-tour lines in any physical variable. Qualitatively, both equa-tions of state show similar behaviour with smooth contoursin EOS19. However, there are accuracy issues in EOS11,which generate unphysical unsteady features at low densitiesand high pressure as will be shown shortly. The missing re-gion in EOS11 is due to a lack of sufficient data points forreliable interpolation.
1D Model
Simulations in 1D were calculated on a radial domain with r ∈ [ , . ] , with initial data is chosen to represent a real-istic air state with mass density ρ = . , velocity u = P ( r ) = P + P ( r ) with P = , P = × Pa · exp ( − ( r / r ) ) with r = T peak ≈ , V current are shown in Figure 4 for FIG. 4. Radial profiles ( N = ) for a component V current forEOS11 (left) and the new EOS19 (right). In both cases, the resultsqualitatively agree, however, the oscillations visiible at about 2,000 Kin the EOS11 results do not appear when using EOS19. EOS19 and EOS11 in comparison. The mathematical formof a component V current profile, and others applied in thiswork, are detailed in the Table I. Both models show similarbehaviour, within the first 10 µ s, the plasma arc channel heatsup quickly, and the pressure increases rapidly reaching peakvalues of P peak = P peak = , , Component
D V (Villa) I ( t ) I (cid:0) e − α t − e − β t (cid:1)(cid:0) − e − γ t (cid:1) I e − α tsin ( β t ) I ,
405 218 , α [ s −
1] 22 ,
708 4 , . β [ s −
1] 1 , ,
530 114 , γ [ s −
1] 10 , ,
100 -time to peak [ µ s ] 3 .
18 13 . µ s ] 34 . . ( nd peak ) TABLE I. Current components applied in simulations lightning striketesting. Component D corresponds toa a single high-amplitude pulsewith decreasing peak current as defined in . Additionally, a compo-nent V current profile has been used in and reflects a high-amplitudeoscillatory exponentially decaying current profile. high. EOS11 displays oscillatory spikes in the pressure pro-file behind the leading shock wave; similar oscillations canbe seen in the corresponding temperature profiles. These fea-tures occur in a temperature range between T < , T < , , where pressure loading dueto the expansion of a cylindrical plasma arc is recorded.Specifically, in their experimental set up, they drilled threeholes into the substrate at radial distances 5cm, 10cm and15cm from the arc attachment point. These holes wereconnected by a 25cm tube to pressure transducers whichrecorded the effects of the passing shock wave and sub-sequent post-shock behaviour. Villa et al. compared theseexperimental results with their numerical simulations usinga one-dimensional approximation for the tube; here this testis used to validate data for the equation of state and plasmamodel. Results for the pressure at the three locations usingthe model presented here, both with EOS11 and EOS19,are shown in Figure 5. The two models show qualitativelysimilar behaviour, although EOS19 predicts a higher pressurepeak than EOS11 upon arrival of the initial wave. However,the subsequent pressure decrease after this peak is then moreconsistent with experimental measurements from Villa forEOS19. The finite width of the tubes is not modelled inthe one-dimensional simulation of the behaviour, hence itis likely that the experimental pressures near the peak areslightly flattened due to geometric effects not captured by thenumerical models. Multimaterial lightning code
Having validated EOS19 in isolation, it is now usedwithin the multimaterial model introduced in Millmore andNikiforakis . Using this, it is possible to investigate the effect FIG. 5. Pressures calculated at the end of transducer tubes of length L = . of the EoS on simulations of arc attachment, both qualitativelyand quantitatively.The experiments of Martins provide imaging of the arcevolution attaching to various substrates, as well as measure-ments of the expansion of both the shock wave generated bythe attachment, and of the arc itself. These tests use a com-ponent D current, as defined in table I, and a simulation inline with a typical laboratory set up is considered, a conicalelectrode geometry with a flattened tip of radius 2 . δ = .For all tests, initial data in the air plasma, outside of thearc, is ρ = . m , u = m/s and P = , , r = P = × , σ ≈ − . Thedomain beneath the substrate is air, modelled as an ideal gaswith γ = . FIG. 6. Results for the plasma arc interacting with an aluminiumpanel (left half) and CFC panel (right half) at t = µ s. The top plotshows temperature and the bottom plot shows pressure. files for the arc attachment to both substrates at t = µ s usingEOS19.In order to understand the difference between the two equa-tions of state for arc attachment, the initial data describedabove is simulated with both EOS11 and EOS19. Resultsare shown for arc attachment to an aluminium substrate attwo time instances, t = . µ s and t = µ s, for pressureand temperature in Figures 7 and 8 respectively. The firstof the two times corresponds to the peak current input timeand there are similar pressure profiles for both EOS11 andEOS19. It is noted that EOS11 results in a lower pressuredirectly above the substrate than for EOS19, whilst the pres-sure within the shock wave appears to be higher, indicatedby the contour plotted at P = × Pa. This higher pres-sure in the peak is likely a result of the oscillatory behaviourwithin the shock wave shown in Figure 4. However, thereare clear differences in the temperature profiles between thetwo equations of state, even at this early time. EOS11 dis-plays unphysical hot spot-like domains, with temperatures ofup to 150 , , FIG. 7. Pressure distribution and stresses for a component D plasmaarc interacting with an aluminium panel. The top image is plottedafter 3.18 µ s, and the bottom plot after 25 µ s. of Martins, where high radiation associated with high temper-ature would be expected to be visible. Additionally, EOS19predicts a more uniformly cylindrical shape of the arc chan-nel, including at the attachment point, as is the expected forarc attachment to aluminium .At t = µ s the current flowing through the arc has de-creased substantially. A characteristic pressure shock wavehas emerged, which moves radially outwards, as can be seenin Figure 7. At this time, there are clear differences in the theshock front comparing EOS11 and EOS19. For EOS 11, theshock wave has not travelled as far (approximately r ≈ r ≈ . P = × Paclearly shows the unphysical oscillatory behaviour seen in theplasma discharge results in Figure 4. As expected, EOS19gives smooth profiles without oscillatory regions throughoutthe shock wave.There are similar differences in the temperature profile inFigure 8. For EOS11, the hot spots observed at earlier timesare still present, and have expanded in size. These hot spotssignificantly alter the behaviour of the attachment point tothe aluminium substrate, effectively widening the attachmentregion. For both equations of state, these results also showthe aluminium substrate mechanically deforming due to thepressure loading on the panel.n improved equation of state for air plasma simulations 9
FIG. 8. Temperature distribution and stresses for a component D plasma arc interacting with an aluminium panel. Temperatures athot spots in EOS11 are substantially higher than 80 , µ s, and the bottom plot after 25 µ s. Figure 9 compares the temperature distribution for thetwo different substrate materials; aluminium and the low-conductivity material. At t = . µ s, the peak temperaturewithin the arc is similar for both material, around 70 , t = . µ s corresponds similarly to a lower pressure. Awayfrom the attachment point, the arc profiles are similar for thetwo substrate materials.Figures 9 and 10 also show the pressure within the sub-strate; pressure rise is caused by a direct loading from the arc,but also Joule heating within the substrate. It is clear thatthe magnitude of the pressure rise is significantly higher forthe low conductivity substrate, indicating the vulnerability togreater damage displayed by such materials.After t = µ s, the pressure shown in Figure 10 shows thatthe radially expanding shock wave is comparable in structurefor both substrates. However, as with the earlier time, at thecentre of the arc, pressure is lower directly above the sub-strate. The temperature plot in Figure 9 shows a significantdifference in the arc root structure between the two substrates. FIG. 9. Temperature distribution and stresses for a component D plasma arc interacting with a carbon composite using EOS19. Thetop image is plotted after 3.4 µ s, and the bottom plot after 46 µ s. This shows that for the low conductivity substrate, the arc rootis significantly wider at the attachment point. This qualita-tively distinct behaviour of the plasma arc root channel be-tween the two substrates is observed experimentally . Ofparticular interest for the attachment to the low conductivitysubstrate is the feature which arises in arc root, close to thepanel, particularly visible in the temperature profile shown inFigure 9; an emerging filament, which separates from the pri-mary arc channel. The physical mechanism behind this is thatthe plasma arc provides a favourable path with lower electricalresistance opposed to the direct path to the substrate. Such afilament is a characteristic phenomenon for carbon compositepanels, which has been observed in experiments by Martins ,where the evolution of the visible arc root radius was mea-sured over time for both substrate materials.The large pressure rise in the low conductivity materialis again visible at after t = µ s, and due to the wider arcin this case, the extent of the high pressure region is widerin comparison to the aluminium substrate. Additionally,throughout the evolution, the lower density of the low con-ductivity material can be seen to result in greater deformationof the panel.In addition to the qualitative improvement shown by EOS19for arc attachment simulations, a comparison is made to themeasurements of Martins . The expansion of both the arcand the shock wave was measured through optical detectionof the attachment process. Figure 11 compares these quan-n improved equation of state for air plasma simulations 10 FIG. 10. Pressure distribution and stresses for a component D plasmaarc interacting with a carbon composite using EOS19. The top imageis plotted after 3.4 µ s, and the bottom plot after 46 µ s. tities for an aluminium substrate; it is clear that both equa-tions of state follow similar patterns, with EOS19 demonstrat-ing an improvement in accuracy. In order to compare the arcand shock expansion to a low conductivity substrate, follow-ing Millmore and Nikiforakis , two measurements are possi-ble; arc growth parallel and perpendicular to the compositeweave. The symmetry in these cases makes the cylindricallysymmetric model presented here valid. The electrical conduc-tivity used for this comparison is 3 × S/m parallel to theweave and 100S/m perpendicular to it. Results for these testsare shown in Figure 12, where it is noted that in the case oflow ‘against weave’ electrical conductivity in the substrate,EOS11 failed to produce a stable simulation. For the parallel‘with weave’ case, results are similar to those in Figure 11,with EOS19 demonstrating an improved fit to the experimen-tal results over EOS11. It is noted that the underprediction inshock radius at early times in the ‘against weave’ case is dueto both the arc itself and the shock being close together, henceit is difficult to measure an accurate radius. At later times,when the two features can be more clearly distinguished, thesimulation data matches the experiment well.
V. CONCLUSIONS
The primary objective of this work was to implement animproved equation of state for air plasma, which can be usedfor various applications involving the partial ionisation of air.
FIG. 11. Comparison of the arc evolution for attachment on alu-minium for EOS11 and EOS19. It is clear that EOS19 captures theevolution behaviour well, with the faster expansion of the shock wavedemonstrated in Figure 7 matching experimental values. The top plotshows the expansion of the arc radius, whilst the bottom plot showsthe shock radius.
In this paper, emphasis was put on the computational multi-physics modelling of solid-plasma interactions. Such numer-ical tools are of high interest for aeronautics to study light-ning strike on new aircraft materials. The equation of state(EOS19) based on the theoretical results of D’Angola et al. ,who considered a 19-species air plasma. Their model is basedon the Debye-Hückel theory for (partially) ionised gas mix-tures, and transport coefficients have been calculated usingthe Chapman-Enskog theory. Comparison to existing open-source EoSs, e.g. Villa et al , show EOS19 can accuratelycapture the evolution of a plasma arc in a non-oscillatory man-ner. Profiles of various air plasma properties are shown to beaccurate over a broad temperature and pressure range. Ad-ditionally, the transport coefficients for thermal conductivityhave been included in the database for usage in other modelsin the future.To validate EOS19, an axisymmetric plasma discharge codewas used, simulating a plasma arc channel under laboratoryconditions. Results are presented in 1D, based on a modelwith fixed current density from Villa et al. and a more real-istic model in 2D axisymmetry, which generates the currentand magnetic field self-consistently based on the plasma evo-lution. The results have been validated against experimentaland numerical results of others . It was addition-ally identified that small errors in the levels of ionisation inan equation of state can lead to oscillatory behaviour; the ac-curacy of EOS19 is demonstrated through the smooth profileswhich arise across the temperature and pressure range.To further test EOS19, it was also applied within a state-of-the-art multimaterial lightning simulation code, which ad-n improved equation of state for air plasma simulations 11 FIG. 12. Comparison of the arc evolution for attachment on the lowconductivity substrate for EOS11 and EOS19. Results are shownfor two cases, with electrical conductivity corresponding to paralleland perpendicular alignment with a composite weave. In the parallel‘with weave’ case, the improvement shown by EOS19 is clear. Inthe ‘against weave’ case, EOS11 was unable to simulate this case,whilst EOS19 again shows good agreement. The top plot shows theexpansion of the arc radius, whilst the bottom plot shows the shockradius. ditionally computes the thermomechanical coupling betweenthe plasma arc root and solid aeronautical skin layers. Resultsare computed for an aluminium and quasi-isotropic CFC panelunder extreme lightning testing conditions by inducing a stan-dardised arc current that is used in lightning-protection testingin laboratories. Results for both materials indicate good agree-ment with experimental measurements by Martins , and nounphysical high-temperature hot spots close to interfaces areobserved when applying EOS19. Additionally, a quantitativeinvestigation of the expansion of both the shock wave causedby attachment, and the arc itself, show that EOS19 achievesimproved accuracy on EOS11, matching the measured be-haviour of these features well.Future work will focus on improved radiative source termwithin the arc in a computationally efficient manner. Ad-ditionally, one may include possible effects associated withthermal conductivity. Finally, the implementation of a fully3D anisotropic equation of state for carbon composite materi-als will be of high interest for the future development of verycomplex materials in aeronautics.To the best of the authors’ knowledge, this state-of-the-artlightning code in combination with EOS19 is the only mul-tiphysics code which considers aeronautical materials andtheir coupled thermomechanical interaction with a plasma arcchannel within a single simulation. Hence, the multimateriallightning code with EOS19 functions as a reliable computa-tional tool to support the design and validation of aeronautical materials concerning lightning protection. ACKNOWLEDGMENTS
The authors would also like to thank Micah Goldade ofBR&T, and Carmen Guerra-Garcia of BR&T and the Mas-sachusetts Institute of Technology, for their input in improv-ing the plasma model. D. Morgan, C. Hardwick, S. Haigh, and A. Meakins, “The interaction oflightning with aircraft and the challenges of lightning testing,” AerospaceLab Journal , 1–10 (2012). O. Ekici, O. Ezekoye, M. Hall, and R. Matthews, “Thermal and flow fieldsmodeling of fast spark discharges in air,” Journal of Fluids Engineering ,55–65 (2007). M. N. Plooster, “Shock waves from line sources. numerical solutions andexperimental measurements,” The physics of fluids , 2665–2675 (1970). M. Akram and E. Lundgren, “The evolution of spark discharges in gases:I. macroscopic models,” Journal of Physics D: Applied Physics , 2129(1996). A. Villa, R. Malgesini, and L. Barbieri, “A multiscale technique for thevalidation of a numerical code for predicting the pressure field induced bya high-power spark,” Journal of Physics D: Applied Physics , 165201(2011). L. Mottura, L. Vigevano, and M. Zaccanti, “An evaluation of roe’s schemegeneralizations for equilibrium real gas flows,” Journal of ComputationalPhysics , 354–399 (1997). A. D’angola, G. Colonna, C. Gorse, and M. Capitelli, “Thermodynamicand transport properties in equilibrium air plasmas in a wide pressure andtemperature range,” The European Physical Journal D , 129–150 (2008). A. D’Angola, G. Colonna, C. Gorse, and M. Capitelli, “Thermodynamicproperties of high temperature air in local thermodynamic equilibrium: Iiaccurate analytical expression for electron molar fractions,” The EuropeanPhysical Journal D , 453–457 (2011). S. Millmore and N. Nikiforakis, “Multi-physics simulations of lightningstrike on elastoplastic substrates,” Journal of Computational Physics ,109142 (2020). M. Capitelli, D. Bruno, and A. Laricchiuta,
Fundamental aspects of plasmachemical physics: transport , Vol. 74 (Springer Science & Business Media,2013). P. Brimblecombe,
Air composition and chemistry (Cambridge UniversityPress, 1996). W. J. O’Sullivan et al. , “Standard atmosphere - tables and data for altitudesto 65,800 feet, report 1235,” Tech. Rep. (International Civil Aviation Orga-nization and Langley Aeronautical Laboratory, 1956). M. Capitelli, G. Colonna, and A. D’Angola,
Fundamental aspects ofplasma chemical physics: thermodynamics , Vol. 66 (Springer Science &Business Media, 2011). G. Colonna and A. D’angola, “A hierarchical approach for fast and accurateequilibrium calculation,” Computer physics communications , 177–190(2004). G. Colonna, “Improvements of hierarchical algorithm for equilibrium cal-culation,” Computer physics communications , 493–499 (2007). J. Haidar, “Non-equilibrium modelling of transferred arcs,” Journal ofPhysics D: Applied Physics , 263 (1999). J. Hilsenrath and M. Klein, “Tables of thermodynamic properties of air inchemical equilibrium including second virial corrections from 1500K TO15,000 K,” Tech. Rep. (National Bureau of standards Gaithersburg MD,1965). M. Capitelli and E. F. Varracchio, “Thermodynamic properties of Ar-H2plasmas,” Rev. Int. Htes Temp. Et. Refract , 195–200 (1977). R. Sevast’yanov and R. Chernyavskaya, “Virial coefficients of nitrogen,oxygen, and air at temperatures from 75 to 2500 K,” Journal of engineeringphysics , 851–854 (1986). L. D. Landau and E. M. Lifshitz, “Statistical physics,” American Journal ofPhysics , 371–372 (1959). M. Capitelli, G. Colonna, C. Gorse, and A. d’Angola, “Transport proper-ties of high temperature air in local thermodynamic equilibrium,” The Eu- n improved equation of state for air plasma simulations 12 ropean Physical Journal D-Atomic, Molecular, Optical and Plasma Physics , 279–289 (2000). R. Devoto, “Transport properties of ionized monatomic gases,” The Physicsof Fluids , 1230–1240 (1966). S. Chapman and T. Cowling, ““the mathematical theory of non-uniformgases,” third edition,” Journal Vol (1970). R. S. Martins,
Étude expérimentale et théorique d’un arc de foudre et soninteraction avec un matériau aéronautique , Ph.D. thesis, Université Paris-Saclay (2016). F. Tholin, L. Chemartin, and P. Lalande, “Numerical investigation of the in-teraction of a lightning and an aeronautic skin during the pulsed arc phase,”International Conference on Lightning and Static Electricity (Toulouse,France) (2015). R. P. Fedkiw, T. Aslam, B. Merriman, and S. Osher, “A non-oscillatoryeulerian approach to interfaces in multimaterial flows (the ghost fluidmethod),” Journal of computational physics , 457–492 (1999). S. K. Sambasivan and H. Udaykumar, “Ghost fluid method for strong shockinteractions part 1: Fluid-fluid interfaces,” AIAA Journal , 2907–2922(2009). S. Schoch, K. Nordin-Bates, and N. Nikiforakis, “An eulerian algorithmfor coupled simulations of elastoplastic-solids and condensed-phase explo-sives,” Journal of Computational Physics , 163–194 (2013). L. Michael and N. Nikiforakis, “A multi-physics methodology for the sim-ulation of reactive flow and elastoplastic structural response,” Journal ofComputational Physics , 1–27 (2018). S. Godunov and E. Romenskii, “Nonstationary equations of nonlinear elas- ticity theory in eulerian coordinates,” Journal of Applied Mechanics andTechnical Physics , 868–884 (1972). G. Miller and P. Colella, “A high-order eulerian godunov method forelastic–plastic flow in solids,” Journal of computational physics , 131–176 (2001). G. Colonna, A. D’Angola, and M. Capitelli, “Electronic excitation andisentropic coefficients of high temperature planetary atmosphere plasmas,”Physics of Plasmas , 072115 (2012). M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, “The FEniCSproject version 1.5,” Archive of Numerical Software , 9–23 (2015). E. F. Toro,
Riemann solvers and numerical methods for fluid dynamics: apractical introduction (Springer Science & Business Media, 2013). A. Villa, R. Malgesini, and L. Barbieri, “Source Code correspond-ing to paper: A multiscale technique for the validation of a numericalcode for predicting the pressure field induced by a high-power spark,”http://gitorious.org/spark2/spark2 [Accessed: 2018, August 14]. “Aircraft lightning test methods, ARP5416, SAE,” Tech. Rep. (SAE, 2013). M. N. Plooster, “Numerical model of the return stroke of the lightning dis-charge,” The Physics of Fluids , 2124–2133 (1971). M. N. Plooster, “Numerical simulation of spark discharges in air,” ThePhysics of Fluids14