Thermal damage analysis of aircraft composite laminate suffered from lightning swept stroke and arc propagation

Xingteng MA, Fusheng WANG,*, Hn CHEN, Donghong WANG, Bin XU

a School of Mechanics, Civil Engineering and Architecture, Northwestern Polytechnical University, Xi’an 710129, China

b Shanxi Key Laboratory of Electromagnetic Protection Material and Technology, 33th Institute of China Electronics Technology Group Corporation, Taiyuan 030032, China

KEYWORDS Ablation damage;Composite laminates;Coupling data transfer;Element deletion;Lightning swept stroke

Abstract Lightning strike is a complicated process involving multi-field coupling. In order to investigate the thermal damage mechanism of carbon fiber reinforced composites subject to lightning swept stroke, a complete numerical method is presented. Numerical model of lightning discharge is established based on Magneto Hydro Dynamics (MHD) and calculated by FLUENT secondary development technology. Considering aerodynamic flow effect, channel formation and evolution process during lightning discharge is analyzed for lightning current waveform A.Thermal-electric Coupling model is presented according to Radial Basis Function (RBF) interpolation theory, which is implemented by compiling program to make lightning current and heat energy inject into composite laminate. Consequently, damage mechanism of composite laminate under lightning swept stroke is studied based on the coupled numerical model and element deletion method. Ablation damage morphology of composite laminate is analyzed to understand plasma expansion and reattachment in arc root. The results show that aerodynamic flow makes the lightning channel move fast and composite laminate is deteriorated due to thermal damage.

1. Introduction

Lightning discharge will release electrical energy,generate high temperature and cause strong electromagnetic radiation,which usually occurs within cloud to cloud discharge,intra cloud discharge or could to ground discharge. It also happens between thunderstorm cloud and aircraft when an aircraft intercepts with a lightning channel branch or is triggered by itself.1As electric charge accumulates in cloud, air will be broken down and a leader channel will be generated when the electric field reaches 25-30 kV/m in thunderstorm. With electric field increasing gradually, a stepped leader gradually propagates downward driven by electronic avalanche. Downward leader and the induced upward leader are connected to form a discharge channel. Lightning flash possesses strong destructive capacity and high potential hazard, which will seriously affect flight safe of vehicles due to lightning direct effects such as thermal ablation damage and indirect effects such as electromagnetic interference on electronic equipment.

As external electromagnetic environment of aircraft, lighting discharge research is very important and primarily concerned about the electrification mechanism, channel formation,heat conduction and electromagnetic field distribution. Although many theories can be used to explain thunderstorm cloud electrification, the actual electrification mechanism is very complicate and involves various electrification mechanisms simultaneously. Lightning discharge direction is affected due to electric charge distribution, humidity and pressure. The main discharge channel is usually shown as irregular bending shape with numerous small branches around it for the reason of partial discharge. Since discharge time is extremely short, current rise rate is generally over 10,000 A/ms that results in electromagnetic field changing rapidly and the temperature ranges from 8000 K to 30000 K owing to Joule heating effect. In this case, air in lightning channel is transformed into thermal plasma because of ionization caused by the interaction of high temperature and electric field. To some extent, lightning channel can be regarded as a thermal plasma channel.

Lightning observation technology is usually utilized to investigate the characteristics of lightning plasma channel.The temperature, electrical conductivity and electron thermal diffusivity can be calculated through lightning spectral information.2Lightning parameters are interesting in engineering application such as current rise time, current duration and peak current, all of them can be available from direct current measurements.3In order to further study lightning discharge mechanism, long-gap discharge test in laboratory and rockettriggered lightning experiment are carried out to duplicate natural lightning discharge due to their similarity. Bazelyan4and Rakov5et al. summarize all aspects of lightning discharge in detail, including lightning physics, hazards, interaction with different objects and protection techniques based on lightning discharge observation.

Many numerical models of lightning discharge have been proposed to study space evolution characteristics, electric distribution and heat transfer. Fractal properties of lightning branched discharge can be described by Dielectric Breakdown Model (DBM)6,7and discharge propagation direction is random depending on electric field distribution. These DBMs are more focused on simulating discharge path. But lightning discharge involves multi-physics response such as heat, electricity and magnetism that are particularly concerned in lightning strike. The physical characteristics of plasma or conductive fluid can be described by Magneto Hydro Dynamics(MHD)method,which is originally applied in studying arc plasma8and gradually introduced to simulate lightning channel. A numerical model based on MHD is established by Gonzalez et al.9to study free burning arc and analyze the influence of cross flow or external magnetic field on electric arc. Chemartin et al.10first adopts MHD method to simulate partial lightning channel to understand discharge physics.Numerical results are verified through lightning discharge experiment conducted by Tanaka et al.11and the distortion of arc column is characterized by normalized length and expansion radius.Overall,MHD method is a new and feasible approach to study the characteristic, mechanism and multiphysics field interaction in lightning discharge.

The interaction mechanism between lightning channel and aircraft is complex, especially when composite material is widely used in aircraft design. Chemartin et al.12investigate direct effect of lightning strike on aircraft structure and explain damage characteristic of composite materials through thermal,electrical and mechanical coupling mechanism. High intensity current is conducted into composite through lightning attachment point.However,high energy is so difficult to be dispersed rapidly that results in thermal damage due to the anisotropic electric conductivity of Carbon Fiber Reinforced Plastics(CFRP).There are at least two attachment points such as entry and exit points when aircraft becomes a part of the current path during lightning strike process.1Lightning channel will move along aircraft surface due to aerodynamic flow effect.Lightning swept stroke makes the damage region enlarge owing to multiple reattachment in dwell time.13Hirano et al.14summarize three kinds of damage morphology such as fiber fracture, resin deterioration and delamination according to lightning strike test of composite laminate. Li et al.15consider that carbon woven fabric laminate has less damage compared with prepreg taped material under lightning strike. Thermal damage behaviors of composite laminate can be predicted by thermal-electrical coupled analysis,which is helpful for understanding the thermal damage process.16Besides, Wang17considers the electric-thermal-mechanical-chemical coupling effects to study the lightning damage of carbon/glass fiber reinforced composite.Nevertheless,lightning current is applied on a point or some points to represent arc attachment region in present works.18-20It is means that only the current is injected into composite.But heat flux in lightning channel is neglected.Besides,lightning attachment region will change due to plasma expansion especially under swept stroke. Therefore, lightning channel should be treated as plasma column to study dynamic thermal damage of composite laminate under lightning strike.

In reality, aerodynamic flow can enlarge thermal damage region when aircraft suffers from lightning strike. Research work on lightning swept stroke and thermal damage on composite structures is seldom reported especially for aircraft composite laminates.21,22This paper is intended to investigate the characteristics of lightning channel evolution in air considering aerodynamic flow effect based MHD method. A complete numerical method for composite laminate under lightning strike is proposed, which include MHD model of lightning channel, interface data transfer based on Radial Basis Function (RBF) method and multi-fields numerical model of composite laminates. It is also study the thermal damage characteristics of composite laminate exposed to lightning swept stroke during dwell time.

2. MHD model of lightning arc plasma

2.1. Assumptions

Lightning arc belongs to magnetic fluid category and MHD theory has been established to study the discharge evolution mechanism. The influence of atmospheric pressure and external electromagnetic is neglected in this model. Local Thermodynamic Equilibrium (LTE) state23is formed in lightning channel owing to high energy transfer promoted by electrons collision. In reality, air plasma departs from LTE slightly and transport coefficients are hard to obtain in this state considering the temperature of electrons is higher than heavy particles.It is commonly assumed that air plasma is in LTE state24for calculation convenience. Besides, the following assumptions are proposed when MHD model is established.

(1) Discharge channel is considered as single phase flow.

(2) Arc plasma is regarded as laminar and Newtonian flow.

(3) Arc plasma is treated as continuous and incompressible.

(4) Gravity effect is neglected.

2.2. Governing equations of lightning arc plasma

MHD equations consist of Navier-Stokes equations and Maxwell equations, which can be used to calculate the electromagnetic field, flow field and temperature field distribution in discharge channel. The continuity equation is written as Eq.(1).Compared with general fluid,Lorenz force J×B and Joule heat J·E are considered in plasma due to electromagnetic effect. Both of them are added in MHD equations. The momentum conservation equation is written as Eq. (2) and energy conservation equation is expressed as Eq. (3),respectively.

In which ρ, t, ν, μ, P, J, B, H, κ, Cpand SMare density,time, velocity vector, dynamical viscosity, pressure, current density,inductive magnetic field,enthalpy,thermal conductivity, specific heat and momentum source, respectively. The energy loss caused by radiation should not be ignored during lightning discharge process. Heat radiation loss term SRof thermal plasma can be expressed by 4π·εN.Net Emission Coefficient(NEC)εNof air thermal plasma is a function of temperature.25The transport properties of air depend on temperature at atmospheric pressure such as dynamical viscosity μ,thermal conductivity κ, specific heat Cpand electrical conductivity σ,which are given in Table 126.

Lorenz force J×B and Joule heat J·E are relevant to current density, electric field and magnetic field, which can be obtained from Maxwell equations. Electric field E mainly depends on electric potential gradient and the relationship is expressed in Eq. (4). Ohm’s law and current conservation law are written as Eq. (5) and Eq. (6), respectively. Consequently, electric potential φ can be calculated by Eq. (7).

According to Maxwell-Ampere law, magnetic field B is related to vacuum permittivity μ0and current density J expressed as Eq. (8). On the other hand, the curl of magnetic vector potential A is equal to magnetic field B based on vector analysis theory, which is written in Eq. (9). Therefore, Kulum Condition is presented by Eq. (10) and magnetic potential A can be obtained through Poisson Eq. (11), in which vacuum permittivity μ0is 4π×10-7N/A2.

The evolution and dynamic characteristic of discharge plasma can be described by Eqs. (1)-(11). There are many methods adopted to solve MHD equations such as Finite Elements Method (FEM),27Element Free Method (EFM)28and Finite Volume Method (FVM).29In this paper, numerical investigation of lightning channel evolution is calculated by FVM in FLUENT software. It is necessary to expand thefunctions in FLUENT software achieved through further development in order to add electromagnetic source terms to Navier-Stokes equations.30

Table 1 Air thermodynamic properties varying with temperature26.

Maxwell equations are defined as additional scalar transport equations in this model by means of User Defined Scalars(UDS). The transport coefficients of thermal plasma varying with temperature are defined by macro function. Momentum source and energy source are calculated by User Defined Functions (UDF). The Governing equations are calculated by Pressure-Based Solver to realize spatial discretization because plasma is treated as incompressible laminar flow. The second-order upwind discretization scheme is selected to calculate electric potential φ and magnetic vector potential A.

3. Evolution of lightning discharge channel

3.1. Lightning discharge model

Based on the lightning strike experiment,20three-dimensional model of lightning discharge is established. The calculation region is composed of cathode, anode plate and air zone shown in Fig. 1. Cathode tip is flattened and its angle is 60°.According to SAE standard,31discharge gap between cathode tip and anode surface is not less than 50 mm.Size of air zone is 500 mm×250 mm×65 mm and the thickness of anode plate is 2 mm. This section focuses on studying lightning channel evolution in air in order to understanding discharge leader development and arc attachment.

The whole computational domain is divided by hexahedral elements and performed on ANSYS/ICEM software. Structured mesh scheme is adopted to ensure mesh quality and calculation accuracy,which is illustrated in Fig.2.Since discharge channel bending mainly occurs in the center region and tip discharge effect occurs at cathode tip, the meshes in these places need to be refined for the purpose of improving calculation accuracy. There are 312,700 hexahedra meshes in the whole computational domain and the minimum cubic mesh is 1.4×10-3mm3.

3.2. Boundary conditions of lightning arc plasma

Fig. 2 Mesh generation for lightning arc plasma calculation.

surfaces are set as open boundary type and electrode surfaces are defined as walls. Temperature of cathode tip is set to 3500 K and the Neumann boundary conditions (∂ /∂n) are applied on other surfaces.Considering lightning channel propagating in open air, boundary conditions of velocity are defined as Dirichlet boundary conditions and the operating pressure is set as standard atmospheric pressure. The components of velocity vector ν in x, y, z direction are denoted as,u, v and w, respectively.

Magnetic vector potential A is also defined as Neumann boundary condition in this model. The components of magnetic vector potential A in x, y, z direction are denoted as Ax, Ayand Az, respectively. Lightning current is imposed on the cathode tip.The simulated current waveform A is adopted to study the direct effects during first return strike. It can be defined as double exponential expression shown in Fig. 3.Where t1is the front rise time and t2is the time to 50%of current peak Ip. The relationship between current density J and potential φ can be expressed as Eq. (12). Current potential is set to 0 V at anode and current density is imposed on cathode tip. The current density J(t) is expressed by Eq. (13).

Here,lightning current I(t)changes along with time and the radius r of cathode tip is 0.15 mm.

3.3. Lightning plasma evolution in air

Boundary conditions of lighting discharge are mainly referred to those of argon arc model proposed by Hsu et al.32The detailed boundary conditions are given in Table 2. Air zone

Fig. 1 Three-dimensional model of lightning discharge.

Table 2 Boundary conditions of three-dimensional lightning arc model.

Fig. 3 Lightning current component A waveform.

The transient evolution characteristics and temperature distribution during lightning discharge are shown in Fig. 4.A downward leader is formed rapidly and develops in space with lightning current increasing.Lightning strike process consists of corona inception, downward propagating, initial attachment and return strike. The surrounding air is initially broken down because of high intensely current on cathode tip. Temperature in plasma increases due to Joule heating effect and the conductivity of air plasma also increase in hot region, which makes the discharge process accelerate. A plasma channel is generated in the vicinity of cathode tip.Temperature of arc channel decreases from the center to periphery and the hottest region appears in the center. Lightning channel evolves from cathode tip toward anode plate under the interaction of Joule heat and Lorenz force.

The discharge channel is similar to a column with a slight twisting and bending before initial attachment,which is shown in Fig.4(a)and(b).Due to plasma impact and stagnation,the arc column expands in attachment region shown in Fig. 4 (c).An upward leader propagate in opposite direction after attachment, which is positive discharge process. The interaction between two discharge leaders with opposite direction can lead to arc channel bending and twisting.Temperature in the vicinity of cathode tip and close to anode plate is relatively higher than that of other regions because of electrode effect. Consequently, attachment area gradually expends because of heat transfer and current conduction on anode surface.

It takes about 25 μs when the arc firstly attaches to anode surface within 50 mm discharge gap,which means propagation velocity of plasma channel is up to almost 2 km/s. When plasma impact on anode plate with supersonic speed, the air is compressed and results in overpressure.Pressure distribution on interface between air and anode plate is shown in Fig. 5.The high pressure on anode surface occurs in arc attachment area and high temperature region.

Current density along the channel is not uniform when plasma arc attach the anode surface. Current density distribution along x axis on interface between air and anode plate is similar to Gaussian distribution shown in Fig.6.As discharge duration increases, arc attachment area increases and current density decreases. Lightning current is delivered into anode plate and Joule heat is related to current density,which means thermal damage mainly depends on current density profile in arc attachment region.

Fig. 4 Evolution characteristics and temperature profile when discharge channel propagates in space.

Fig. 5 Pressure profile in lightning attachment region.

Fig. 6 Current density along radial direction in lightning attachment region.

3.4. Plasma channel propagation during lightning swept discharge

Evolution process of lightning channel has been simulated and studied above. But it ignores aerodynamic flow effect. In reality, there is a relative motion between flying aircraft and arc.Lightning arc attaching to aircraft surface usually lasts for a few milliseconds. Therefore, the arc root can be displaced on aircraft surface and forms a swept channel when aircraft moves at a typical flight speed. Lightning swept discharge model is established to further investigate the evolution characteristics of arc plasma channel and study the damage mechanism of anode plate during arc root moving. It is assumed that aircraft flies at a speed of 100 m/s. Anode plate is considered immobile in calculation domain.The incoming flow has a velocity equal to 100 m/s in the y direction and is applied with a Blasius profile. In this case, boundary layer thinness δ is set to 10 mm.Because the boundary layer around a flying aircraft is about 10-20 mm.33

Evolution and temperature distribution during discharge is shown in Fig. 7, which is effected by aerodynamic flow in y direction. Discharge channel is bending at arc root because of boundary layer effect and the arc reattaches to plane surface. For this reason, the whole arc column gradually moves along the y direction and is affected by cathode tip at the same time. Discharge channel is extended by swept stroke and plasma expands near cathode tip at t=0.892 ms. As energy loss is caused by aerodynamic flow, temperature profile gradually decreases in dwell time.

4. Damage of composite laminate under lightning swept discharge

4.1. Interface data transfer model

Strong current and high temperature will be transferred into composite laminate through reattachment points during swept stroke process. Here, numerical method is proposed to study the damage characteristics of composite laminate subject to lightning swept stroke. Generally, the meshes generated in fluid domain are more refined than solid domain,which results in mismatch at interface. It is difficult to ensure the load data accurately transfer at interface between lightning arc plasma and composite laminate. There are many interpolation methods34-37for fluid-structure coupling calculation and Radial Basis Function (RBF) as a global interpolation approach has the higher accuracy.37Here, transient data such as pressure,heat flux and current density are transferred to composite laminate to investigate the dynamic response in lightning stroke based on RBF interpolation method.

The fluid meshes and structure meshes on interface are shown in Fig. 8. There are a total of Nfnodes on the fluid interface and Nsnodes on structure interface, which can be expressed as Eqs. (14) and (15) in three-dimensional.

Some physical quantities such as current density, pressure and heat flux of the nodes on fluid interface and structure interface are defined as gs（xsi） and gf（xf）, respectively. The relationship between them can be expressed as Eq. (16). Interface data transfer is realized by means of transfer matrix Hsf.According to RBF method, transfer matrix Hsfcan be written as Eq. (17), in which coefficient matrix Asf, M, P and Mpare detailed described in Ref. 37. Interpolation calculation is realized by compiling the program on MATLAB.

4.2. Finite element model and boundary conditions

Carbon fiber/epoxy resin composite laminate is made of 16 layers having a set of ply orientations [45/-45/02/45/90/-45/0]sand each layer is 0.125 mm. Its size is 500 mm×250 mm×2 mm. Composite material properties of IM600/133 are listed in Table 3 including density,specific heat,thermal conductivity(κx,κy,κz)and electrical restively(RSVx,RSVy, RSVz), which are related to temperature.

Fig. 7 Evolution and temperature distribution of lightning swept channel under aerodynamic flow.

Fig. 8 Finite element mesh on fluid-structure coupling interface.

Because epoxy resin matrix is insulate at low temperature,electrical conductivity of composite material is poor along transverse and in-depth direction.Decomposition temperature of epoxy resin matrix is about 400°C and carbon fiber is seriously carbonized and fractured when temperature exceeds 3000°C. Epoxy resin matrix will be gasified rapidly in lightning arc attachment region under high temperature and lightning current mainly conducts along carbon fiber before vaporization. Transverse and in-depth conductivity of composite laminate is increased when temperature exceeds 400°C. The ablation temperature of Composite laminate defined as 3316°C, which means carbon fiber are ablated or vaporized. It can be assumed that composite laminate has excellent electrical conductivity due to phase transition when the temperature exceeds 3316°C.19

Finite element model of composite laminate and boundary conditions are shown in Fig.9,which is established in ANSYS software. Hexahedral mesh is generated in each layer of composite laminate and SOLID5 element is selected for multi-field coupling calculation.There are total 28800 elements and 32147 nodes. Arc initial attachment occurs in the central region and moves along y direction. Electrical potential is set as 0 V at bottom face and side faces of composite laminate to represent ground connection. Thermal radiation is 0.9 at upper surface and side faces.The environment temperature is 25°C and bottom face is defined as adiabatic boundary condition. Load data of nodes are obtained from RBF interpolation method and imposed on upper layer surface. Here, nonlinear solution scheme is selected in view of the material properties depending on temperature.

4.3. Damage characteristics of composite laminate

Current and heat conducts rapidly along the fiber direction after lightning arc attachment. The main energy source within composite laminate originates from plasma heat flux in lighting arc channel and Joule heat caused by conductivity anisotropy. It is commonly suggested that composite laminate can be defined as pyrolysis damage when temperature between 400-600°C, which means only epoxy resin is pyrolyzed while carbon fiber can still conduct current in this state. In this coupled analysis model, elements will be considered invalid or deleted when both carbon fiber and matrix are pyrolyzed.Because its average temperature exceeds ablation temperature 3316°C.

Table 3 Thermal and electrical properties of carbon fiber/epoxy composite laminate.16,19,22

Fig. 9 Finite element model and boundary conditions of composite laminate.

Ablation area of composite laminate propagates with arc moving on surface and forms an ablation path shown in Fig. 10 at typical time. Initial attachment occurs in the center of upper surface at 0.0241 ms and current mainly conducts along carbon fiber with 45° direction. Because dwell time of arc plasma is extremely short, the heat energy generated by Joule heating effect and heat flux is not enough to make epoxy resin matrix decompose. Arc plasma channel expands in attachment region with lightning current injecting continuously. At this time, the high temperature region on composite laminate surface increases. Ablation damage area appears at 0.166 ms because of carbon fiber and resin are vanished in high temperature region,in which elements are deleted.In addition,pyrolysis vanished is accelerated through reattachment in lightning arc root due to aerodynamic flow effect. The radius of arc root and reattachment region will decrease as a result of the thermal energy in arc channel will be weakened during swept process. Finally, an ablated path on CFRP surface is generated by arc swept stoke at 1.078 ms.

In order to understand the internal decomposed degree,temperature profile and ablative damage area of each layer are shown in Fig. 11 at different typical time. Joule heat is mainly generated along the fiber direction on each layer.Although damage area increases, the expansion rate is lower due to energy dissipation in arc plasma.There is no doubt that most serious damage occurs on up layer, which suffers from heat flux and Joule heat at the same time and the maximum ablation damage area is up to 1589.97 mm2.Ablation damage area of each layer decreases from top to bottom due to the poor electrical conductivity and thermal conductivity in depth direction.

Damage behaviors such as delamination in composite laminate can be study by analyzing the damage morphology in cross-section. A cross section perpendicular to x-axis is selected to investigate internal pyrolysis damage. For convenient observation, arc attachment region in composite laminate is partially enlarged for detailed study ablation process along depth direction. A series of sectional calculation results are shown in Fig.12 at different typical times.It is found that the ablation path is extended when arc root is displacing and the ablation length is up to 63.1 mm within 1.078 ms. Meanwhile, the pyrolysis of carbon fiber and matrix expands to fourth layers in depth direction.The maximum ablation thickness is up to 0.375 mm during swept process of arc channel.

Lightning current mainly conducts on the upper surface along carbon fibers and is hard to conduct in depth direction due to the poor conductivity of epoxy resin when lightning arc initially attaches. Internal temperature of composite laminate increases rapidly as dwell time prolonging.Lightning current can conducts to the second layer owing to the temperature high enough to decompose epoxy resin. From the calculation results, it is shown that ablation damage occurs in initial attachment region at 0.166 ms.Nevertheless,injection position of lightning current changes when arc root is displaced during swept stroke process. Lightning current rapidly conducts in depth direction subject to high temperature. Ablation damage decreases with the lightning current intensity weakening.

5. Conclusions

Fig. 10 Ablation path propagation of composite laminate under lightning swept stroke.

Fig. 11 Temperature profile and ablation damage area of each layer.

Fig. 12 Ablation propagation in depth direction of composite laminate during sweeping process.

In this study, an arc plasma channel is obtained by numerical method to understand lightning discharge process. Dynamic and evolution characteristics of arc plasma channel are investigated when the speed of airflow is 100 m/s. BRF interpolation method is used to ensure the data information transfer accurately in fluid-structure coupling interface. It finally realizes the thermoelectric coupling numerical calculation when composite laminate exposed to arc plasma channel during swept stroke process.

(1) When the heat energy accumulates,lightning plasma leader develops rapidly from cathode tip to anode plate with the air electrical conductivity increasing Thermal plasma impingement and stagnation on anode plate surface makes arc root expand after initial attachment. At the same time,overpressure occurs in attachment region.

(2) Lightning discharge channel obviously distorts and expands in sweeping process. Driving by airflow, reattachment in arc root causes the whole arc plasma channel displaced from initial attachment position.Temperature in arc plasma channel decreases over time.But it is not less than 10000 K.

(3) An ablation path is generated on upper layer due to high energy injection from thermal plasma channel consistently. Thermal ablation damage in-plane and in-depth direction is analyzed in order to understand damage mechanism and estimate the decomposition degree of composite laminate under lightning swept stroke.

Acknowledgements

This study is supported by the National Natural Science Foundation of China(No.51875463)and the Natural Science Basic Research Plain in Shaanxi Province of China (No.2018JM1001)and the Research Funds of Shanxi Key Laboratory of Electromagnetic Protection Material and Technology(No. 33ZD1807KF001C).