Finite Element Analysis of Lightning Damage Factors Based on Carbon Fiber Reinforced Polymer

While carbon-fiber-reinforced polymers (CFRPs) are widely used in the aerospace industry, they are not able to disperse current from lightning strikes because their conductivity is relatively low compared to metallic materials. As such, the undispersed current can cause the vaporization or delamination of the composites, threatening aircraft safety. In this paper, finite element models of lightning damage to CFRPs were established using commercial finite element analysis software, Abaqus, with the user-defined subroutines USDFLD and HEAVEL. The influences of factors such as the structural geometry, laminate sequence, and intrinsic properties of CFRPs on the degree of damage to the composites are further discussed. The results showed that when a current from lightning is applied to the CFRP surface, it mainly disperses along the fiber direction in the outermost layer. As the length of the CFRP increases, the injected current has a longer residence time in the material due to the increased current exporting distance. Consequently, larger amounts of current accumulate on the surface, eventually leading to more severe damage to the CFRP. This damage can be alleviated by increasing the thickness of the CFRP, as the greater overall resistance makes the CFRP a better insulator against the imposed current. This study also found that the damaged area increased as the angle between the first two layers increased, whereas the depth of the damage decreased due to the current dispersion between the first two layers. The analysis of the electrical conductivity of the composite suggested that damage in the fiber direction will be markedly reduced if the conductivity in the vertical fiber direction increases approximately up to the conductivity of the fiber direction. Moreover, increasing the thermal conductivity along the fiber direction will accelerate the heat dissipation process after the lightning strike, but the influence of the improved thermal conductivity on the extent of the lightning damage is less significant than that of the electrical conductivity.


Introduction
In the aircraft industry, carbon-fiber-reinforced polymers (CFRPs) have gradually replaced traditional metallic materials in parts such as skins and fairings due to their excellent mechanical properties, including their high specific strength and corrosion resistance. More than 50% of the Boeing 787 Dreamliner and Airbus 350XWB aircrafts are made of composites [1,2]. However, it is inevitable that an aircraft will encounter lightning strikes during flights. According to statistics published by the Royal Canadian Air Force, a commercial plane is struck by lightning once a year on average, and the probability increases in thunderstorm-prone areas [3][4][5]. When a lightning strike occurs, up to 200 kA of electricity can be generated and released in microseconds [6][7][8]. Due to insulating resins that exist in CFRPs, the conductivity of CFRPs is much lower than that of metallic materials, and it transfer analysis. The thermal and electrical damage from the Joule heat generation process was calculated for each finite element.
The electric field in each grid unit was calculated from the Maxwell equation and shown as: where V is the unit volume, J is the unit current density, S is the unit cross-sectional area, r c is the unit current volume density, and n is the outward normal to S. Using this definition of the electrical field, the current density can be described with Ohm's law: where σ E is the electrical conductivity, E is the electrical field, and ϕ is the electrical potential. The Joule heat P ec can be calculated according to Joule's law as: The average Joule heat P ec over the entire time of the lightning strike can be expressed as: where E and σ E are the electrical field and electrical conductivity at time t, respectively, and the heat energy converted from the electricity can be written as: where the η v is the energy conversion factor. Therefore, the temperature field for the lightning damage model included lightning Joule heat generation, heat transfer, and dissipation. According to the energy balance, the heat conduction can be summarized as: where ρ, θ, C v , k, q, and r represent the density, temperature, specific heat, thermal conductivity, heat flux per unit area, and heat generation density, respectively.

Numerical Model
The lightning strike damage model for the CFRP was created using the commercial software ABAQUS, as shown in Figure 1. The model was based on a quasi-isotropic laminate with the sequence of [45/0/-45/90] 2s . The dimensions of the model were 150 × 150 × 4.6 mm 3 , and the thickness of each layer was set to 0.2875 mm. As the crosssection view shown in Figure 1, the mesh size in the thickness direction is the superposition of the thickness of each layer, instead of being divided according to the total thickness. The bottom 8 layers were simplified as an integrated layer to reduce the number of calculations, whereas the top half of the model was designed according to different laminate sequences, which was verified in simulations where the damage depth of composites was approximately 1.6 mm (5-6 layers) at 40 kA. The lightning current was injected into the model at the center node of the top surface using a node-loading form. The model was divided into a DC3D8E grid type, and a total of 23,868 solid elements were used to simulate the electrical-thermal propagation. The lightning current was injected into the model at the center node of the top surface, and the voltage was presumed to be zero at the subface and four side surfaces. a DC3D8E grid type, and a total of 2,3868 solid elements were used to simulate the electrical-thermal propagation. The lightning current was injected into the model at the center node of the top surface, and the voltage was presumed to be zero at the subface and four side surfaces. When the lightning current was injected, the heat was immediately transferred through the composite. Consequently, the thermal emissivity of the entire composite surface was set to 0.9, the convective heat transfer coefficient was 1, and the ambient temperature was 25 °C. The ablation behavior of the resin was considered a decay process, with a degree of pyrolysis C: In equation (7), W is the sample quantity; and i and f represent the initial and final states of the sample, respectively. Properties of each unit were modeled as a function of C, and the epoxy resin was assumed to be almost completely ablated at about 600 °C based on previous thermogravimetric analysis results evaluated from 25 °C to 1000 °C with a constant heating rate of 15 K/min. The thermal degradation was assumed to occur primarily in the carbon fibers. Therefore, the properties of the model CFRP after a degree of pyrolysis of 1 were assigned based on the measured material properties above 600 °C [20,21]. The specific thermo-electro-mechanical parameters are shown in Tables 1 and 2. Note: Cp is the specific heat capacity; is the coefficient of thermal expansion; k is the thermal conductivity; xx, yy, and zz represent the in-fiber direction, vertical fiber direction, and thickness direction, respectively. When the lightning current was injected, the heat was immediately transferred through the composite. Consequently, the thermal emissivity of the entire composite surface was set to 0.9, the convective heat transfer coefficient was 1, and the ambient temperature was 25 • C. The ablation behavior of the resin was considered a decay process, with a degree of pyrolysis C: In Equation (7), W is the sample quantity; and i and f represent the initial and final states of the sample, respectively. Properties of each unit were modeled as a function of C, and the epoxy resin was assumed to be almost completely ablated at about 600 • C based on previous thermogravimetric analysis results evaluated from 25 • C to 1000 • C with a constant heating rate of 15 K/min. The thermal degradation was assumed to occur primarily in the carbon fibers. Therefore, the properties of the model CFRP after a degree of pyrolysis of 1 were assigned based on the measured material properties above 600 • C [20,21]. The specific thermo-electro-mechanical parameters are shown in Tables 1 and 2. Note: Cp is the specific heat capacity; α T is the coefficient of thermal expansion; k is the thermal conductivity; xx, yy, and zz represent the in-fiber direction, vertical fiber direction, and thickness direction, respectively. However, the heating rate of the composites is extremely fast during lightning strikes, and resin matrix decomposition cannot be treated as an isothermal process [22], implying that the C value calculated from isothermal thermogravimetric analysis before 600 • C is inaccurate. Therefore, this study simulated the thermal degradation of the composite using kinetic equations, and the thermal degradation process was described as a n-level chemical kinetic rate equation as follows [23]: where t, n, and T are the pyrolysis time, reaction order of the pyrolysis process, thermodynamic temperature, respectively, and k(T) is the kinetic rate constant determined from the Arrhenius equation and related with activation energy E a : where A and R are the pre-exponential factor and gas constant, respectively. According to data published in the work by Ogasawaral et al., [14], the parameters used in the simulations were defined as n = 3.5, A = 3.0 × 10 15 s −1 , E a = 180 J/(mol·K), and R = 8.314 J/(mol·K). Figure 2 displays the simulated results for the resin ablation process at different heating rates using the chemical kinetics equation. It can be found that the shapes of the degree of pyrolysis curves at different heating rates were basically identical, but the onset temperature increased with increasing heating rates. Therefore, the onset temperature can be used to revise the non-isothermal chemical kinetics equation to better match the conditions of lightning strikes, as discussed below.
However, the heating rate of the composites is extremely fast during lightning strikes, and resin matrix decomposition cannot be treated as an isothermal process [22], implying that the C value calculated from isothermal thermogravimetric analysis before 600 °C is inaccurate. Therefore, this study simulated the thermal degradation of the composite using kinetic equations, and the thermal degradation process was described as a nlevel chemical kinetic rate equation as follows [23]: where t, n, and T are the pyrolysis time, reaction order of the pyrolysis process, thermodynamic temperature, respectively, and k(T) is the kinetic rate constant determined from the Arrhenius equation and related with activation energy Ea: where A and R are the pre-exponential factor and gas constant, respectively. According to data published in the work by Ogasawaral et al., [14], the parameters used in the simulations were defined as n = 3.5, A = 3.0 × 10 15 s −1 , Ea = 180 J/(mol·K), and R = 8.314 J/(mol·K). Figure 2 displays the simulated results for the resin ablation process at different heating rates using the chemical kinetics equation. It can be found that the shapes of the degree of pyrolysis curves at different heating rates were basically identical, but the onset temperature increased with increasing heating rates. Therefore, the onset temperature can be used to revise the non-isothermal chemical kinetics equation to better match the conditions of lightning strikes, as discussed below.
According to 30 kV/mm of the electrical breakdown strength of epoxy resin [24], the electrical breakdown of the composite with a 55% fiber content was set as 13.5 kV/mm. When the electrical breakdown intensity of grids was beyond 13.5 kV/mm, the electrical properties of the grids were converted to 0.1 S/mm. According to 30 kV/mm of the electrical breakdown strength of epoxy resin [24], the electrical breakdown of the composite with a 55% fiber content was set as 13.5 kV/mm. When the electrical breakdown intensity of grids was beyond 13.5 kV/mm, the electrical properties of the grids were converted to 0.1 S/mm.
In addition to resin ablation by the current, the sublimation temperature of the carbon fibers was set to 3316 • C by transforming the latent heat from 4.8 × 10 3 kJ/kg to 4.3 × 10 4 kJ/kg [16]. The lightning current waveform of component A in SAE ARP5412 with a modified peak of 40 kA was used in the experimental lightning tests on the CFRP specimens. According to this waveform, the corresponding theoretical waveform used in the present numerical model is plotted in Figure 3.
In addition to resin ablation by the current, the sublimation temperature of the carbon fibers was set to 3316 °C by transforming the latent heat from 4.8 × 10 3 kJ/kg to 4.3 × 10 4 kJ/kg [16]. The lightning current waveform of component A in SAE ARP5412 with a modified peak of 40 kA was used in the experimental lightning tests on the CFRP specimens. According to this waveform, the corresponding theoretical waveform used in the present numerical model is plotted in Figure 3.  Figure 4, the user-defined subroutines USDFLD and HEAVEL were employed to describe the thermal degradation field and electrical breakdown field, respectively. In each analytical step, the transient electrical and the thermal field analysis were performed sequentially. In the calculative process, the superposition calculations began with the degree of thermal degradation when the thermal field reached the initial ablation temperature value, and then the corresponding thermoelectric properties of the material for the given degree of thermal degradation were updated accordingly in each analysis step. Eventually, the obtained final model results for the degree of lightning damage were the thermal degradation field and the electrical breakdown field.  Figure 5 shows the simulated damage to the material with different onset temperatures of resin ablation. The in-plane damage area included the superposition of surface damage in each layer, which corresponded to the ultrasonic C-scan result; the depth damage area represented the damage area along the thickness direction, which corresponded to the ultrasonic B-scan result. The damage depth claimed the maximum damage depth in the thickness direction. In particular, the simulated damage results agreed well with  Figure 4, the user-defined subroutines USDFLD and HEAVEL were employed to describe the thermal degradation field and electrical breakdown field, respectively. In each analytical step, the transient electrical and the thermal field analysis were performed sequentially. In the calculative process, the superposition calculations began with the degree of thermal degradation when the thermal field reached the initial ablation temperature value, and then the corresponding thermoelectric properties of the material for the given degree of thermal degradation were updated accordingly in each analysis step. Eventually, the obtained final model results for the degree of lightning damage were the thermal degradation field and the electrical breakdown field. In addition to resin ablation by the current, the sublimation temperature of the carbon fibers was set to 3316 °C by transforming the latent heat from 4.8 × 10 3 kJ/kg to 4.3 × 10 4 kJ/kg [16]. The lightning current waveform of component A in SAE ARP5412 with a modified peak of 40 kA was used in the experimental lightning tests on the CFRP specimens. According to this waveform, the corresponding theoretical waveform used in the present numerical model is plotted in Figure 3.  Figure 4, the user-defined subroutines USDFLD and HEAVEL were employed to describe the thermal degradation field and electrical breakdown field, respectively. In each analytical step, the transient electrical and the thermal field analysis were performed sequentially. In the calculative process, the superposition calculations began with the degree of thermal degradation when the thermal field reached the initial ablation temperature value, and then the corresponding thermoelectric properties of the material for the given degree of thermal degradation were updated accordingly in each analysis step. Eventually, the obtained final model results for the degree of lightning damage were the thermal degradation field and the electrical breakdown field.  Figure 5 shows the simulated damage to the material with different onset temperatures of resin ablation. The in-plane damage area included the superposition of surface damage in each layer, which corresponded to the ultrasonic C-scan result; the depth damage area represented the damage area along the thickness direction, which corresponded to the ultrasonic B-scan result. The damage depth claimed the maximum damage depth in the thickness direction. In particular, the simulated damage results agreed well with  Figure 5 shows the simulated damage to the material with different onset temperatures of resin ablation. The in-plane damage area included the superposition of surface damage in each layer, which corresponded to the ultrasonic C-scan result; the depth damage area represented the damage area along the thickness direction, which corresponded to the ultrasonic B-scan result. The damage depth claimed the maximum damage depth in the thickness direction. In particular, the simulated damage results agreed well with the experimental results, when the simulated onset temperature was 320 • C. A detailed morphology comparison can be found in Figures 6 and 7 that the simulated in-plane and depth damage morphologies were directly compared with experimental results in reference [25]. In addition, the simulation setting condition corresponded to the experimental setup from reference [25]: the subface of the specimen was earth-grounded, and the distance between the tip of the discharge probe and the specimen surface was adjusted to about 2.0-3.0 mm. The lightning current waveform of t1/t2 of 4/20 µs with a peak of 40 kA was conducted in the experiment. The simulated results are in good agreement with the experimentally measured damage morphologies, indicating that Joule heating was the primary source of lightning damage to the composite. The red and blue regions of the diagram represent the fully and slightly decomposed regions of the composite, respectively, with C values between 0 and 1. However, other factors such as shockwaves, thermal expansion, pyrolysis gas, and electrical breakdown could result in further damage to the composite. Thus, the experimental damage is expected to be more severe than the simulated results [19]. Nevertheless, the non-isothermal pyrolysis method adopted in this study still showed superior simulated results compared to previous work based on isothermal pyrolysis results. Basically, the simulated result from the non-isothermal pyrolysis method predicted more severe damage. In particular, the simulated severely damaged region is closer to the experimentally observed morphology, because this non-isothermal pyrolysis method involves the overall consideration of the temperature-, time-, and space-dependent material properties as a result of a lightning strike.

As shown in
ence [25]. In addition, the simulation setting condition corresponded to the experimental setup from reference [25]: the subface of the specimen was earth-grounded, and the distance between the tip of the discharge probe and the specimen surface was adjusted to about 2.0-3.0 mm. The lightning current waveform of t1/t2 of 4/20 μs with a peak of 40 kA was conducted in the experiment. The simulated results are in good agreement with the experimentally measured damage morphologies, indicating that Joule heating was the primary source of lightning damage to the composite. The red and blue regions of the diagram represent the fully and slightly decomposed regions of the composite, respectively, with C values between 0 and 1. However, other factors such as shockwaves, thermal expansion, pyrolysis gas, and electrical breakdown could result in further damage to the composite. Thus, the experimental damage is expected to be more severe than the simulated results [19]. Nevertheless, the non-isothermal pyrolysis method adopted in this study still showed superior simulated results compared to previous work based on isothermal pyrolysis results. Basically, the simulated result from the non-isothermal pyrolysis method predicted more severe damage. In particular, the simulated severely damaged region is closer to the experimentally observed morphology, because this non-isothermal pyrolysis method involves the overall consideration of the temperature-, time-, and spacedependent material properties as a result of a lightning strike.

Influence of the Composite Structural Parameters on the Extent of Lightning Damage
The injected current from the lightning strike was primarily dissipated along the conductive fiber in the CRFP. In addition, changes in the structural geometry of the composite and laminate material sequences also determined the length of the path to the ground. In this section, the structural geometry parameters that can be tuned to improve the internal current export from the CRFP to reduce the associated lightning damage are discussed.

Influence of the Composite Structural Parameters on the Extent of Lightning Damage
The injected current from the lightning strike was primarily dissipated along the conductive fiber in the CRFP. In addition, changes in the structural geometry of the composite and laminate material sequences also determined the length of the path to the ground. In this section, the structural geometry parameters that can be tuned to improve the internal current export from the CRFP to reduce the associated lightning damage are discussed. With an increase in the composite thickness, the size of the in-plane damaged area decreased gradually, while the extent of damage in the thickness direction tended to increase (note: the damage in the sample with a 1 mm thickness penetrated the whole sample). In addition, in the CFRP specimen with a thickness of 7 mm, the breakdown damage area can barely be traced, because the specimen can withstand the voltage surge from the lightning current channel. The damaged area was circular in shape for small thicknesses (Figure 8c). With an increase in the composite thickness, the in-plane damage expanded along the fiber direction, while the damage vertical to fibers direction was gradually reduced.
Because of the insulating resin in the interply, the impedance in the thickness direction of the composite was assumed to decrease as the thickness of the resistive resin decreasing, which resulted in an incremental increase in the injected current in the thickness direction. When the CRFP thickness was small, it can be supposed that more of the injected current was conducted along the thickness direction, and thus, there were a less in-plane current and indistinguishable anisotropic in-plane damage. The simulation results for thinner composites showed evidence of a larger dielectric breakdown area, which supported that the current was concentrated in the attached region of the composite and conducted along the thickness direction. On the contrary, as the thickness increased, the associated greater impedance resulted in more of the current being dispersed along the in-plane fiber direction and more in-plane damage while less damage occurred in the depth direction [25].
for thinner composites showed evidence of a larger dielectric breakdown area, which supported that the current was concentrated in the attached region of the composite and conducted along the thickness direction. On the contrary, as the thickness increased, the associated greater impedance resulted in more of the current being dispersed along the inplane fiber direction and more in-plane damage while less damage occurred in the depth direction [25].

Side Length
In contrast to the trends seen with changes in the composite thickness, the in-plane damage decreased, as the side length increased (note: The composite with the side length of 37.5 mm exhibited the fully penetrated in-plane damage). Figure 9a,c shows that the degree of damage along the fiber direction first increased and then decreased as the side length increased. This non-monotonic trend was generated, because the current injected on the surface of the composite was not dissipated as the

Side Length
In contrast to the trends seen with changes in the composite thickness, the in-plane damage decreased, as the side length increased (note: The composite with the side length of 37.5 mm exhibited the fully penetrated in-plane damage). Figure 9a,c shows that the degree of damage along the fiber direction first increased and then decreased as the side length increased. This non-monotonic trend was generated, because the current injected on the surface of the composite was not dissipated as the conduction path increased with the increasing side length. As a result, the current accumulated in the area struck by lightning, inducing a strong electrical field and eventually causing a large dielectric breakdown, as shown in Figure 9b.
From the analysis above, it can be found that the structural geometry of the composite mainly affected the current export path. Therefore, it is essential to decrease the length of the export path to minimize the degree of damage to the composite. However, reducing the composite thickness did not fully alleviate the damage, because the low conductivity along the thickness direction led to more severe damage from the current. Due to the electrical conductivity of the carbon fibers, the conduction of the current along the fiber direction was promoted by decreasing the side length of the composite, thus promoting the dispersion of the current from a lightning strike. Therefore, shortening the dispersion path through a strategic structural design is beneficial for reducing lightning damage to composites, which can be applied to the fuel tanks fabricated by the CFRP. In general, subface ground is unfeasible and dangerous for fuel tanks, because it may introduce fuel ignition by lightning current injection. As a result, the distance of the grounded rivets installed in CFRPs should be decreased to shorten the dispersion path of the current. Moreover, thinner composites are more prone to dielectric breakdown. conduction path increased with the increasing side length. As a result, the current accumulated in the area struck by lightning, inducing a strong electrical field and eventually causing a large dielectric breakdown, as shown in Figure 9b. From the analysis above, it can be found that the structural geometry of the composite mainly affected the current export path. Therefore, it is essential to decrease the length of the export path to minimize the degree of damage to the composite. However, reducing the composite thickness did not fully alleviate the damage, because the low conductivity along the thickness direction led to more severe damage from the current. Due to the electrical conductivity of the carbon fibers, the conduction of the current along the fiber direction was promoted by decreasing the side length of the composite, thus promoting the dispersion of the current from a lightning strike. Therefore, shortening the dispersion path through a strategic structural design is beneficial for reducing lightning damage to composites, which can be applied to the fuel tanks fabricated by the CFRP. In general, subface ground is unfeasible and dangerous for fuel tanks, because it may introduce fuel ignition by lightning current injection. As a result, the distance of the grounded rivets installed in CFRPs should be decreased to shorten the dispersion path of the current. Moreover, thinner composites are more prone to dielectric breakdown.

Laminate Sequence
The electrical properties of different ply sequences can be derived using a transformation matrix R [26]:

Laminate Sequence
The electrical properties of different ply sequences can be derived using a transformation matrix R [26]: (12) [σ] = where σ ' is the conductivity of the origin ply, and σ is the conductivity of the transformed ply; the subscripts 11, 22, and 33 represent the fiber direction, the in-plane vertical fiber direction, and the thickness direction, respectively; m and n are parameters of the rotation matrix, namely m = cosα and n = sinα (α is the layup angle). This equation suggests that the laminate sequence can be designed to ensure multi-directional conductivity and create more current export paths to improve the internal current dispersion in composites.
To evaluate the influence of the laminate layer sequence on the current dissipation, the angles of the fibers in each layer were varied, namely [0] 16 4s . To simplify the labelling scheme, these models are noted using the 2nd layer sequence as 0, 30, 45, and 90, respectively.
The simulation results are shown in Figure 10, and models 0 and 30 had larger damaged areas than the other models. Because the current flow paths in the 1st layers for all models were identical, whereas the current paths in the adjacent 2nd layer distributed along another direction, which was dedicated to varying the current dispersion between adjacent layers. When the angle between adjacent layers increased, the current could disperse more widely, resulting in decreasing the extent of Joule heating. Therefore, the size of the damaged area reduced with an increase in the angle difference between the 1st and 2nd layers, but the depth of the damage in the transverse direction increased at the 1st layer and then dropped. From the perspective of the damage morphology, all the composite models with different layer sequences showed a distribution of morphologies parallel to the fiber direction in the 1st layer. Model 0 exhibited the maximal degree of damage, whereas model 90 showed the least damage. From the discussion above, the distance of the export path was a key factor in determining the Joule heating inside the model. Obviously, models 0 and 90 had the shortest exports paths. However, in model 0, the fibers are aligned in the same direction, resulting in the electrical current concentrating in the fiber direction [27]. From Figure 11, the nephogram of the electrical potential showed that the model 0 composite sequence had a higher concentration of the current compared to those of the other models, thus causing the accumulation of Joule heat in the 0 • direction and thereby leading to more severe damage. In contrast, the electrical potential area expanded with the angle difference between the 1st and 2nd layers, resulting in a wider distribution and a lower degree of electrical potential. Therefore, in addition to shortening the current export path, it is essential to disperse the electrical current inside the composite as much as possible to minimize the damage from a lightning strike. A sequence of laminates with orthogonal fibers is recommended to minimize damage from lightning compared to other ply sequences.

Influence of the Composite Properties on the Extent of Lightning Damage
Although many factors can influence the extent of the lightning damage to a composite, the thermal and electrical properties of materials undoubtedly play a dominant role. In this section, anisotropies in the thermal and electrical properties of the composites are discussed.

Electrical Conductivity
As shown in Figure 12, increasing the electrical conductivity along the fiber direction (σ11) effectively enhanced the current dissipation in the composite and reduced the extent of the in-plane damage, but the damage along the depth of the composite was unchanged.

Influence of the Composite Properties on the Extent of Lightning Damage
Although many factors can influence the extent of the lightning damage to a composite, the thermal and electrical properties of materials undoubtedly play a dominant role. In this section, anisotropies in the thermal and electrical properties of the composites are discussed.

Electrical Conductivity
As shown in Figure 12, increasing the electrical conductivity along the fiber direction (σ 11 ) effectively enhanced the current dissipation in the composite and reduced the extent of the in-plane damage, but the damage along the depth of the composite was unchanged. The damage was concentrated along the fiber direction, and the length of damage extended with increasing σ 11 ( Figure 13). Because the increased electrical conductivity facilitated the current propagation along the fiber direction, a lower current was conducted along other directions in the composite [18,28,29].

Materials 2021, 14, x FOR PEER REVIEW
The damage was concentrated along the fiber direction, and the length of dam tended with increasing σ11 (Figure 13). Because the increased electrical conductivi itated the current propagation along the fiber direction, a lower current was con along other directions in the composite [18,28,29].
Enhancing the electrical conductivity along the vertical fiber direction (σ22 a had little effect on the simulated lightning damage, until σ22 and σ33 were close to th of σ11. For instance, the extent of damage did not sharply improve, until σ22 increa to 1000 times (1.145 × 10 3 ·S/m). Because a slight improvement in the electrical condu along the transverse fiber direction did not affect the current conduction in the la the current was still mainly dissipated along the fiber direction. When the electric ductivity along the transverse direction was close to σ11, the current was transporte more paths. Consequently, the damage was more homogeneous, and the overall of damage was lower.     Enhancing the electrical conductivity along the vertical fiber direction (σ 22 and σ 33 ) had little effect on the simulated lightning damage, until σ 22 and σ 33 were close to the value of σ 11 . For instance, the extent of damage did not sharply improve, until σ 22 increased up to 1000 times (1.145 × 10 3 ·S/m). Because a slight improvement in the electrical conductivity along the transverse fiber direction did not affect the current conduction in the laminate, the current was still mainly dissipated along the fiber direction. When the electrical conductivity along the transverse direction was close to σ 11 , the current was transported along more paths. Consequently, the damage was more homogeneous, and the overall degree of damage was lower. Figure 14 reveals the influence of the thermal conductivity of the composite on the resulting lightning damage. It can be found that the thermal conductivity had a minimal effect on the degree of lightning damage. Because the lightning process is instantaneous, even if thermal conductivity increased by orders of magnitude, the high current conduction still resulted in significant amounts of Joule heating in microseconds, and heat dissipation in such a short time hardly affected the extent of the damage. However, increasing the thermal conductivity of the composite did promote thermal dissipation after the lightning strike, and as a result, the heat was less concentrated and cooling time was faster in the highly damaged areas near the lightning strike point (Figure 14b). resulting lightning damage. It can be found that the thermal conductivity had a minimal effect on the degree of lightning damage. Because the lightning process is instantaneous, even if thermal conductivity increased by orders of magnitude, the high current conduction still resulted in significant amounts of Joule heating in microseconds, and heat dissipation in such a short time hardly affected the extent of the damage. However, increasing the thermal conductivity of the composite did promote thermal dissipation after the lightning strike, and as a result, the heat was less concentrated and cooling time was faster in the highly damaged areas near the lightning strike point (Figure 14b). In Figure 15, because the high-temperature region induced by the Joule heating was distributed along the fiber direction, increasing the thermal conductivity along the fiber direction (K11) facilitated thermal conduction. Thus, high temperature was dissipated in shorter cooling durations, and thereby the size of the highly damaged area decreased with increasing K11. In contrast, increasing the thermal conductivity vertical to the fiber direction (K22) did not affect the heat transfer path. As a result, the size of the damaged area did not change with changes in K22. Increasing the conductivity along the thickness direction (K33) reduced the degree of the lightning damage. Heat transfer along the thickness direction is the shortest path compared to the other directions. Therefore, improving K33 promoted rapidly heat dissipation along the thickness direction, as seen in Figure 14a, where the depth of damage increased with increasing K33 as more heat was conducted along the thickness direction. In Figure 15, because the high-temperature region induced by the Joule heating was distributed along the fiber direction, increasing the thermal conductivity along the fiber direction (K 11 ) facilitated thermal conduction. Thus, high temperature was dissipated in shorter cooling durations, and thereby the size of the highly damaged area decreased with increasing K 11 . In contrast, increasing the thermal conductivity vertical to the fiber direction (K 22 ) did not affect the heat transfer path. As a result, the size of the damaged area did not change with changes in K 22 . Increasing the conductivity along the thickness direction (K 33 ) reduced the degree of the lightning damage. Heat transfer along the thickness direction is the shortest path compared to the other directions. Therefore, improving K 33 promoted rapidly heat dissipation along the thickness direction, as seen in Figure 14a, where the depth of damage increased with increasing K 33 as more heat was conducted along the thickness direction.

Specific Heat Capacity
As shown in Figure 16, the specific heat capacity of the composite determined the amount of heat energy required to raise the temperature of the material per unit of mass, and the composite was supposed to be more thermally stable as the specific heat capacity increased. Thus, the degree of damage was effectively reduced, as the specific heat capac-

Specific Heat Capacity
As shown in Figure 16, the specific heat capacity of the composite determined the amount of heat energy required to raise the temperature of the material per unit of mass, and the composite was supposed to be more thermally stable as the specific heat capacity increased. Thus, the degree of damage was effectively reduced, as the specific heat capacity of the composite increased. However, to be effective against a lightning strike, the heat capacity needed to be as high as 10 kJ/(kg·K), and only a few gases had such high specific heat capacities. Therefore, increasing the specific heat capacity of solid materials is very hard to realize. Figure 15. Morphologies of the lightning damage to composites with varying thermal conductivities.

Specific Heat Capacity
As shown in Figure 16, the specific heat capacity of the composite determined the amount of heat energy required to raise the temperature of the material per unit of mass, and the composite was supposed to be more thermally stable as the specific heat capacity increased. Thus, the degree of damage was effectively reduced, as the specific heat capacity of the composite increased. However, to be effective against a lightning strike, the heat capacity needed to be as high as 10 kJ/(kg·K), and only a few gases had such high specific heat capacities. Therefore, increasing the specific heat capacity of solid materials is very hard to realize.

Conclusions
The effects of structural geometry, laminate sequences, and intrinsic properties of composites on the degree of lightning damage were studied using numerical simulations, and the main conclusions were as follows: 1. The structural geometry of the composite determines the dissipation path of the current from a lightning strike. To improve the lightning strike protection, it is necessary to reduce the distance of the export path. In addition, the sequence of the laminates in the composites can affect the current dispersion and the export distance.
2. The current export path should be minimized in the design of composite structures for lightning protection. The sequence of [0/90] 4s is optimal for dispersing electrical potential and shortening the current export path.
3. The thermal and electrical conductivities are strongly influenced by the anisotropic structure of the composite. In terms of the fiber direction, improving the conductivity can facilitate thermal or electrical conduction to reduce the Joule heating. The conductivity along the transverse fiber direction should be increased to realize lightning current or energy path transformation.
4. Improving the electrical conductivity is a more effective means of enhancing the lightning protection of a composite than improving the thermal conductivity. Improving the electrical conductivity requires enhancing the current export path, whereas improving the thermal conductivity is beneficial for dissipating the heat after lightning strikes.
This work aimed to provide a reference for the design of composites for better lightning strike protection. In future work, experiments should be performed to determine if the predicted structural designs and property improvements better protect composite materials against lightning damage. Moreover, the presented work serves as a reference for the optimization of approaches to protect aircrafts from lightning damage.

Data Availability Statement:
The data presented in this study are available upon reasonable request from the corresponding author.