Thermal Ablation Damage Analysis of CFRP Suffering from Lightning Based on Principles of Tomography

Coupled electrical–thermal finite element analysis (FEA) models are widely adopted to analyze the thermal ablation damage of carbon fiber reinforced polymer (CFRP) caused by lightning, but it is still difficult to analyze the ablation due to its complex space geometry. According to the principle of computerized tomography (CT), tomographic images of FEA models’ temperature fields with different thicknesses were obtained to calculate the mass loss and compare the damage morphology. The four areas including Area 0, Area I, Area II, and Area III; were separated from the temperature fields in terms of different vaporization and pyrolysis temperature ranges of carbon fiber (CF) and resin matrix. Ablation mass losses were calculated by pixel statistics and tomographic intervals, which were consistent with the experimental results. The maximum ablation area of unprotected CFRP was found on the tomography images of 50 μm rather than the surface by comparing tomographic images with different thickness due to the influence of the thermal radiation, but this effect was not found in CFRP protected by copper mesh. Some other phenomena, including continuous evolutions of ablation areas and the influence of the intersection angle on the direction of the ablation extension, were also discovered.


Introduction
The application research of carbon fiber reinforced polymer (CFRP) is in a stage of rapid development, and its application scope has achieved a leap from the secondary and small bearing component of aircraft to the main and large bearing component, and even the emergence of full composite aircraft [1][2][3]. Compared with metal materials, CFRP has better properties in terms of specific strength, specific modulus, corrosion resistance, fatigue resistance, and thermal stability. It is normal for an aircraft to be struck by lightning in flight, with the generally accepted probability of once every 3000 h. The traditional metal structures of the entire aircraft can overlap each other to form a good conductor, which can smoothly guide the lightning current, generate the "Faraday cage" effect, and protect the airborne equipment [4,5]. However, with the increasing application of the resin matrix composites in aircraft, its conductivity is gradually decreased, which makes these aircraft more vulnerable to lightning than before.
Carbon fiber and resin are the main components of CFRP, while some impurities can appear in CFRP such as curing agent and defoaming agent, etc. In fact, the amount of impurities is very According to the standard ASTM D7137 and principle of lightning ablation [25][26][27], two types of specimens were designed and fabricated, including unprotected CFRP laminates (Group A) and copper-mesh-protected CFRP laminates (Group B). The ply sequence of both types of laminates is [45/-45/0/90/90/-45/0/45/0/90/-45/45] with the single ply thickness of 0.15mm and the general size of 150mm×100mm×3.6mm. The copper mesh is vertically woven and its mesh number is 1100 with the copper wire diameter of 0.25mm. The appearance of the specimens is shown in Figure 1.  The boundary conditions and contact properties of FEA models were set according to the actual experimental conditions: The boundary heat flux of models was set to 0 W/m 2 , the surface thermal emissivity of CFRP laminates was 0.9 and that of the copper mesh was 0.78, the environmental temperature and the initial temperature of the specimens were both 25 • C, and the bottom and side potential in the models were 0V. There was often contact electrical and thermal resistance between different media. The electrical and thermal conductivity of the contact interface between copper mesh and CFRP laminates were reduced to 80% of the copper mesh at the same temperature [28]. The condition settings of the FEA models is shown in Figure 3.

Material Properties
The conductivity of CFRP is mainly determined by the carbon fiber (CF) inside, and the anisotropy of conductivity is determined by the different orientations of plies, but copper is isotropic material. The conductivity of CFRP and copper shows different trends with increasing temperature, which is because that the resin matrix in CFRP is an insulator which will melt when heated, while the copper is a conductor whose conductivity will decrease with increasing temperature. Property variations with temperature of CFRP and copper as shown in Tables 1 and 2.   Table 1. Property variations with temperature of CFRP [29]. T is the temperature, ρ is the density, C p is the specific heat, σ x , σ y , and σ z are, respectively, the longitudinal, the transverse, and the in-depth electrical conductivity; λ x , λ y , and λ z are, respectively, the longitudinal, the transverse, and the in-depth thermal conductivity.  Table 2. Property variations with temperature of copper [30]. T is the temperature, ρ is the density, C p is the specific heat, σ is the electrical conductivity, λ is the thermal conductivity.

The Waveform of Lightning Currents
On the basis of SAE ARP-5412B, the lightning current waveforms were classified into 4 types: Component A (the first return stock), Component B (the intermediate current), Component C (the continuing current), and Component D (the subsequent return stroke). Typical lightning current waveforms are shown in Figure 4 [30].  Table 2.
Component D was selected as the lightning current load to analyze the direct impacts of CFRP lightning strikes. Six combinations of currents and protection modes are shown in Table 3. Table 3. Combinations of currents and protection modes. Meanings of I p , t 1 , and t 2 are, respectively, shown in Figure 4b; E is the integral action; Q is the charge transfer.

Observation of the Ablation Damage
Ablation damages were first observed by temperature field of FEA models and CT, which are shown in Figure 5. The 3D reconstruction image of A20 was chosen, as other reconstructions were not sharp enough because of the warped layer and the complex structure of copper mesh both disrupting the propagation of the rays. The FEA temperature field of A20 showed clearly the inner complex shape of ablation rather than a taper, reflecting that a more precise observing method was needed. The CF inside the CFRP could be clearly found by 3D reconstruction based on CT, but only the first ply was reconstructed, as much more continuous CT images will be required for other unnecessary inner plies without damage. In the single CT image of Figure 5b, the absence of copper mesh with high light was found, but the CFRP was still intact, which was consistent with the temperature field. Therefore, combining the advantages of CT and FEA, a new analysis method of ablative damage observation will be described in the following.

Temperature Fields Zoning
Both CF and resin are organic polymer materials and belong to amorphous, whose vaporization and pyrolysis temperatures are not fixed. It is generally believed that 300 • C is the initial vaporize and pyrolysis temperature of the resin, and 573.5 • C is the completion temperature. The vaporize and pyrolysis of CF occur at 3042-3316 • C [31]. Four areas including Area I, Area II, and Area III were divided from the temperature fields calculated from the coupled electrical-thermal FEA models, and the resin content in CFRP is 40%. Properties of different temperature areas are shown in Table 4. The appearance of Area 0, Area I, Area II, and Area III also look different, 45 • unidirectional CFRP in different temperature areas are shown in Figure 6. Damage was not observed in Area 0, which still showed an intact appearance. Slight thermal ablation started to be established in Area I due to the fiber laying direction, which showed a part of resin was lost. Warpage and exposure of CF were found in Area II due to the full vaporization and pyrolysis of resin; the mechanical properties in this area had almost been lost, although the CF seemed to be intact. Breakages and losses of CF were found in Area III, which indicated that the temperature here had reached the vaporization pyrolysis point of CF causing the most serious damage.

Calculation of Ablative Mass Based on Tomographic Images
The lightning ablation region of CFRP extends along different directions in different plies due to the anisotropy of CFRP, which leads to the complex geometric space of ablation. Ultrasonic scanning is usually channeled into finding the internal damage, while some other researchers have simplified it into an inverted cone to calculate the damage volume [32]. The ablative space of CFRP is segmented along the thickness direction, as shown in Figure 7. Tomographic images at different thickness were extracted to discretize the continuous ablation space. The total ablative mass loss is expressed by Equation (1).
where m is the total ablative damage mass; Ι i S , ΙΙ i S , and i S ΙΙΙ are the areas of Area Ⅰ, Ⅱ, and Ⅲ in the image of Section i; ρ is the density of CFRP; hi is the distance between adjacent sections. If his are all equal between different adjacent sections, Equation (1) can be simplified as: Pixel statistics were used to calculate the area of different areas in temperature fields of FEA. The relationship between the number of the pixels and the area is expressed as: where i S Ι,ΙΙ,ΙΙΙ and t i S are the areas of Area Ⅰ, Ⅱ, or Ⅲ and the total area of Section I; i n Ι,ΙΙ,ΙΙΙ and t i n are the pixels' number of Area Ⅰ, Ⅱ, or Ⅲ and the total area of Section i. When Equation (3) is introduced into Equation (2), the ablation mass loss is expressed as:

The Ablation Mass Loss
The distance between adjacent sections was set as 50 μm in the calculation of temperature fields of Group A, which represents that three tomographic images will be obtained in every ply. Each tomographic image was numbered as "a-b" where a is the ply number between 1 and 24 and b is the Tomographic images at different thickness were extracted to discretize the continuous ablation space. The total ablative mass loss is expressed by Equation (1).
where m is the total ablative damage mass; S I i ,S II i , and S III i are the areas of Area I, II, and III in the image of Section i; ρ is the density of CFRP; h i is the distance between adjacent sections.
If h i s are all equal between different adjacent sections, Equation (1) can be simplified as: Pixel statistics were used to calculate the area of different areas in temperature fields of FEA. The relationship between the number of the pixels and the area is expressed as: where S I,II,III i and S t i are the areas of Area I, II, or III and the total area of Section I; n I,II,III i and n t i are the pixels' number of Area I, II, or III and the total area of Section i.
When Equation (3) is introduced into Equation (2), the ablation mass loss is expressed as:

The Ablation Mass Loss
The distance between adjacent sections was set as 50 µm in the calculation of temperature fields of Group A, which represents that three tomographic images will be obtained in every ply. Each tomographic image was numbered as "a-b" where a is the ply number between 1 and 24 and b is the sections' number in every ply between 1 and 3. The summary of tomographic images of Group A is shown in Figure 8.  The distance between adjacent sections was set as 5µm in the calculation of temperature fields of Group B due to the short damage depth of Group B and its significant variety along the thickness. A summary of the tomographic images of Group B is shown in Figure 9. Images of different areas were selected through color recognition based on Figures 8 and 9. The number of pixels was counted out and the area, volume, and mass of different areas were calculated according to Equations (2)-(4). The tomographic summary in different areas of A20 are shown in Figure 10. Results of pixel statistics and ablation mass losses obtained by Equations (2)-(4) are shown in Table 5. Table 5. Results of pixel statistics and ablation mass losses. S t is the total area; n t is the total number of pixels; n I , n II , and n III are, respectively, the total pixels' number of Area I, Area II, and Area III; ∆m CFRP and ∆m Cu are, respectively, the ablation mass losses of CFRP and copper mesh.

Number
S t (cm 2 ) n t n I n II n III ∆m CFRP (g) ∆m Cu (g) Mass losses of specimens were also obtained by weighing as shown in Table 6. Variations of the mass losses of the six specimens and their FEA models demonstrated the same trends with the variations of currents, which will increase with the increase in the peak current. In Group A, the FEA results were less than the experimental results, because there were still any other forms of energy dissipation including air exchange, thermal convection, etc., in addition to thermal radiation during experiments. In Group B, the FEA results were significantly less than the experimental results due to copper mesh falling off at a high temperature [33].

The Ablation Morphology
The maximum ablation area of CFRP in Group A was found on the tomographic images of 50 µm instead of the surface by comparing tomographic images in different thickness. This phenomenon was called "subsurface effect", whose main reason is the thermal radiation on surface; although the discharge time was just a few hundred microseconds, energy dissipation caused by the thermal radiation is considerable at a temperature difference of several thousand degrees Celsius. On the contrary, no thermal radiation to the external environment exists inside the CFRP, which counteracts the energy reduction caused by the weakening of the electric field, leading to the ablation area on the subsurface larger than on the surface. Surface temperature field images superimposed by tomographic images of 50 µm were found closer to the experimental results. The results of Group A are shown in Figure 11. Unlike Group A, the "subsurface effect" of Group B was not found in Figure 9, which means that the maximum ablation area of CFRP still appears on the surface by comparing the tomographic images. Furthermore, the ablation area extension in all directions of Group B was almost the same, while that of Group A was not. The reason why these two led to different results is that the conductivity of copper is much higher than that of CFRP, so the copper mesh conducted most of the currents during lightning. The isotropic conductivity of the copper mesh results in a nearly circular ablation area in Group B. The Joule effect is concentrated on the copper mesh and its contact interface, leading to almost no current passing through underground surface of CFRP. The results of Group B are shown in Figure 12.
Materials 2020, 13, x FOR PEER REVIEW 13 of 19 Figure 11. Results of Group A. Surface temperature field images of FEA, surface-subsurface superimposed temperature field images of FEA, specimen photos, and their ultrasonic C-scan images are shown from left to right.
Unlike Group A, the "subsurface effect" of Group B was not found in Figure 9, which means that the maximum ablation area of CFRP still appears on the surface by comparing the tomographic images. Furthermore, the ablation area extension in all directions of Group B was almost the same, while that of Group A was not. The reason why these two led to different results is that the conductivity of copper is much higher than that of CFRP, so the copper mesh conducted most of the currents during lightning. The isotropic conductivity of the copper mesh results in a nearly circular ablation area in Group B. The Joule effect is concentrated on the copper mesh and its contact interface, leading to almost no current passing through underground surface of CFRP. The results of Group B are shown in Figure 12. The ablation extension direction of copper-mesh-protected specimens is also influenced by the cross angle of copper wire. The conductivity of copper mesh will be anisotropic when the cross angle is not 90°, which will lead to different extension sizes of the ablation area in different directions. The ablation extension of rhombic braided copper mesh is exhibited in Figure 13 [30]. The ablation extension direction of copper-mesh-protected specimens is also influenced by the cross angle of copper wire. The conductivity of copper mesh will be anisotropic when the cross angle is not 90 • , which will lead to different extension sizes of the ablation area in different directions. The ablation extension of rhombic braided copper mesh is exhibited in Figure 13 [30]. The FEA result of B20 was in good agreement with its experimental result, followed by B40 and B60, with all ablation occurring on the copper mesh and the CFRP surface formed from the intact CF of 45° in ultrasonic C-scan images. The maximum temperature of B20 was 343.6 °C, which was much lower than the vaporization temperature of copper but higher than its oxidation temperature, being consistent with the specimen of B20 only showing oxidation discoloration [34]. The ablation area and The FEA result of B20 was in good agreement with its experimental result, followed by B40 and B60, with all ablation occurring on the copper mesh and the CFRP surface formed from the intact CF of 45 • in ultrasonic C-scan images. The maximum temperature of B20 was 343.6 • C, which was much lower than the vaporization temperature of copper but higher than its oxidation temperature, being consistent with the specimen of B20 only showing oxidation discoloration [34]. The ablation area and shape of B40 seemed the same in FEA and the experiments, while those of B60 were slightly different between the experiments and FEA. It can be inferred that although the maximum temperature of CFRP has not reached 3316 • C, the steam generated by the vaporization and pyrolysis of resin will cause the copper mesh to fall off at high temperatures and reduce its mechanical properties. Some other differences including the ablative pits of B20 and the oxidation discoloration and warping at the edge of B40 or B60 were also observed. The reason for the above phenomenon is that the actual current channel is a local energy concentration at the intersection of copper wires rather than a uniform plasma flow. The contact resistance between lapped copper wires or between specimens and fixtures also plays an important role, which will increase in ablation area of B40 and B60. The microscopic image of the central region of B20 specimen is shown in Figure 14.

Evolutions of Ablation Areas
In the past, the entire ply was usually scanned by ultrasound or sectioned by tomography images of FEA regardless of the ply thickness. Due to the limitation of the technical accuracy, discontinuity or mutation occurred in the evolution of ablation area, which was shown in Figure 15 [6,35]. Different from the above results, it was found that even in the same plies, the evolutions of ablation areas are continuous. The continuity of the evolution of ablation areas also conforms to the basic principle of "the macroscopic material change is continuous" in classical physics [36]. It was proved that the conductivity of different plies will affect its own electric field as well as the electric field of other plies, leading to the continuous change in the ablation area. Superimposed ablation area boundaries are shown in Figure 16. Different from the above results, it was found that even in the same plies, the evolutions of ablation areas are continuous. The continuity of the evolution of ablation areas also conforms to the basic principle of "the macroscopic material change is continuous" in classical physics [36]. It was proved that the conductivity of different plies will affect its own electric field as well as the electric field of other plies, leading to the continuous change in the ablation area. Superimposed ablation area boundaries are shown in Figure 16. Boundaries 1-1, 1-2, and 1-3 extended along the l1, becoming wider and shorter as the tomography deepens. The extended scale of Boundary 2-1 was located between Ply 1 and Ply 2, which looked like a square and symmetrical along l2, because it is located at the junction of the two plies. Boundary 2-2 and Boundary 2-3 were both symmetrical along l3, which become longer and narrower with the tomography deepening.

Fitting Analysis between Ablation Mass Losses and Currents
The integration of action and the charge transfer are the two main parameters of the currents, and the former can directly represent the current energy [37]. With the currents of same waveforms, the integration of action is proportional to the square of the amplitude, and the charge transfer is proportional to the amplitude. The quadratic curve fitting between ablation mass losses and peak currents can reflect the influence of the current action integral on the ablation damage. The quadratic curve fitting of mass losses between these two groups in experiments and FEA is shown in Figure 17. Boundaries 1-1, 1-2, and 1-3 extended along the l 1 , becoming wider and shorter as the tomography deepens. The extended scale of Boundary 2-1 was located between Ply 1 and Ply 2, which looked like a square and symmetrical along l 2 , because it is located at the junction of the two plies. Boundary 2-2 and Boundary 2-3 were both symmetrical along l 3 , which become longer and narrower with the tomography deepening.

Fitting Analysis between Ablation Mass Losses and Currents
The integration of action and the charge transfer are the two main parameters of the currents, and the former can directly represent the current energy [37]. With the currents of same waveforms, the integration of action is proportional to the square of the amplitude, and the charge transfer is proportional to the amplitude. The quadratic curve fitting between ablation mass losses and peak currents can reflect the influence of the current action integral on the ablation damage. The quadratic curve fitting of mass losses between these two groups in experiments and FEA is shown in Figure 17. The curve of "∆m = A·I peak 2 " is used for fitting in Figure 17. The parameter analysis is demonstrated in Table 7. The Coefficient A of the experiment was slightly less than that of the FEA by 14.75% in Group A, which indicates that the experimental results are in good agreement with the simulation results, while other forms of energy dissipation should be further considered in the FEA. However, the Coefficient A of the experiment is significantly greater than that of the FEA in Group B due to the surge of B60 results of the experiment, and the factors of copper mesh fracture and shedding should be further taken into consideration in the future.

Conclusions
(1) According to different degrees of vaporization and pyrolysis of carbon fiber and resin matrix at different temperatures, the temperature field of CFRP in lightning is divided into Area 0, Area I, Area II, and Area III with different density loss rates. The calculation accuracy will be improved if the temperature range is divided more finely without considering the simplicity of the calculation process. (2) The areas of Area I, Area II, and Area III in different tomographic images are calculated by pixel statistics. The total mass loss can be obtained by multiplying the area with the corresponding spacing, density, and mass loss rate. Ablation mass losses are well-matched with that of FEA. The calculation accuracy will increase by shortening the spacing without considering the simplicity of the calculation process. (3) The maximum ablative area of unprotected CFRP laminates appears underground rather than on the surface due to the surface thermal radiation, while this effect does not appear in copper mesh protected CFRP laminates, because the copper mesh bears most of the lightning current, leading to the ablation area being close to circular at the same time. (4) The ablative areas of CFRP gradually evolve with the increase in the thickness even in the case of the same ply, and the copper mesh can significantly reduce the anisotropy of ablation areas of CFRP. The anisotropy of ablation areas is closely associated with the cross angle of copper wire, where the ablation areas are close to circular if the angle is 90 • , and the ablation areas are close to ellipse if the angle is not 90 • . (5) The quadratic curves are used to fit ablation mass losses. Since the integral of current action is proportional to the lightning energy, the quadratic coefficient of the fitted curve can reflect the ablation energy.

Conflicts of Interest:
The authors declare no conflict of interest.