Sage Journals: Discover world-class research

Abstract

In this work, three-dimensional (3D) finite element simulations were undertaken to study the effects of lightning strikes on the microscale behaviour of continuous fibre-reinforced composite materials and to predict and understand complex lightning damage mechanisms. This approach is different from the conventional mesoscale or macroscale level of analysis, that predicts the overall lightning damage in composite laminates, thus providing better understanding of lightning-induced thermo-mechanical damage at a fundamental level. Micromechanical representative volume element (RVE) models of a UD composite laminate were created with circular carbon fibres randomly distributed in an epoxy matrix. The effects of various grounding conditions (one-, two-, and four-side grounding), fibre volume fractions (V_f = 55, 60%, and 65%), and peak current amplitudes (10, 20, and 40 kA) on microscopic damage in RVE models were characterised to understand fibre, epoxy matrix, and fibre-matrix interfacial damage associated with a lightning strike. Thermal damage, estimated based on epoxy matrix decomposition temperature, was largely constant up to 20 kA before increasing significantly at 40 kA. Severe thermal damage steadily decreased with increasing V_f. Thermo-mechanical damage was predicted using a ductile plasticity model with Drucker–Prager yield criterion for epoxy matrix failure, and cohesive surfaces for fibre-matrix interface debonding. Thermal strain had the largest contribution to thermo-mechanical damage, while dynamic pressure loading was negligible in all RVE models. The RVE model proposed in this work, to the best of the authors’ knowledge, is the first model predicting lightning-induced fibre-matrix interfacial damage.

Keywords

Polymer-matrix composites thermomechanical micro-mechanics finite element analysis lightning damage prediction

Introduction

Lightning strikes present a potential risk to the structural integrity of composite materials, particularly those used in primary aircraft structures. Lightning strikes produce both direct effect damage (physical damage near the lightning attachment) and in-direct effects (appreciable, non-structural damage associated with coupling with on-board electronics) on aircraft composites. Lightning direct damage is caused by a combination of various thermal and mechanical loading, i.e., Joule heating emanating from lightning current flow, direct plasma heat fluxes impinging the lightning attachment location, and dynamic mechanical loading.^1,2 Therefore, research into lightning strike damage and protection systems, has been increasing for many years and has followed parallel experimental and simulation pathways.^3–12

Experimental research has been conducted with the primary focus on idealised standard lightning waveforms with ≥100 kA peak current amplitudes, named Waveforms A or D^13–15 since these waveforms produce ample electrical current, which is instantaneously dissipated as heat (i.e., thermal damage). Direct lightning damage in typical carbon fibre-reinforced polymer (CFRP) composites includes a number of significant damage modes including fibre damage (i.e., fray, fracture, or blow-out), matrix damage (i.e., cracking or thermal decomposition), delamination, and through-thickness damage penetration.¹⁶

Modelling of lightning strikes has typically fallen into two categories focussing on either (i) the external plasma conditions at the lightning attachment point, due to rapid ionisation of the high temperature and high pressure air in the arc channel^17–20 or (ii) the composite material internal damage mechanisms divided into three broad areas; thermal-electric,^5,7,16,19,21 mechanical pressure loading^10,22–24 and thermal expansion^25,26 analyses. The majority of this work, presented in open literature, has been conducted at the mesoscale level, where each composite ply is modelled as a homogenous medium (thus, no damage characteristics of each constitute) and inter-laminar interactions are considered between plies using gap conductance or cohesive zone models. These models have almost exclusively used a simulated specimen measuring 150 mm × 100 mm^{5,16,19,21,27} and a specimen layup of [45/0/-45/90]_4s,^{16,21,27–29} as inspired by the pioneering work of Hirano et al.¹³

Sun et al.³⁰ recently studied the lightning‐induced conduction characteristics of CFRP composites. The authors suggested that at lower electric field strength (the intensity of an electric field at a particular location) electrons cannot overcome the interfacial barrier in the “fibre‐resin-fibre” unit and are therefore not involved in the conduction process. Rather, the low‐amplitude direct current (DC) or alternating current (AC) conduction paths are controlled only by the contact of neighbouring fibres. However, for higher intensity current loading (50–1500 A), contact between neighbouring fibres was not sufficient and therefore mixed conduction paths “fibre–fibre” and “fibre‐resin‐fibre” both contribute to current flow. The authors used their findings and assumptions to inform mesoscale modelling of resulting lightning strike damage. While this approach could inform mesoscale modelling, representative volume element (RVE)-based microscale modelling can characterise unique damage mechanisms of each constitute and better illuminate the current flow through and along constitutes under realistic lightning strike loading.

Representative volume element modelling

The RVE method, using a small representative volume to represent the whole material, is widely applied to characterize the macroscopic response of composites and primarily used for effective property calculation. RVE modelling has been widely used in the study of two-dimensional (2D) and three-dimensional (3D) woven composites^31–35 and also in multi-scale modelling of unidirectional (UD) composites.^36–41 Totry et al.⁴² used a 2D RVE model to predict the failure of a composite laminate under transverse compression and out‐of‐plane shear loading, while Sun et al.⁴³ used a 3D RVE model to predict the failure of a composite laminate under combined transverse compression and out‐of‐plane/in‐plane shear bi-axial loading. Several studies have typically used RVE modelling to obtain effective material properties from a representative unit cell. Chen et al.⁴⁴ proposed a 3D RVE micromechanical model to capture the coupled effects of fibres, matrix, and fibre-matrix interface on failure under different multiaxial loading conditions. The authors used an RVE model (fibre diameter d_f = 7 μm and fibre volume fraction V_f = 60%) measuring 50 μm × 50 μm × 10 μm. The interface thickness was identified using an artificial neural network (ANN) and interface properties were taken from literature. González and Llorca⁴⁵ modelled the mechanical behaviour of UD carbon or glass fibre polymer-matrix composites. A square RVE model (d_f = 10 μm and V_f = 50%) was used with randomly distributed parallel fibres created using a modified random sequential adsorption (RSA) algorithm.

RVE modelling has been limited within the field of lightning strike research. Guo et al.⁴⁶ and Karch et al.²³ both proposed RVE-based micromechanics models to estimate the effective properties of expanded metal foils, the most common lightning strike protection (LSP) layers. In the work done by Guo et al.,⁴⁶ the RVE models consisted of a 6 × 6 unit-cell, where one unit was a single grid in the long width of the diamond (LWD) and short width of the diamond (SWD) directions. The effective electrical conductivities were calculated to evaluate the lightning shielding performance of expanded foils subjected to 100 kA Waveform D. Karch et al.²³ used similar RVE models, but with periodic boundary conditions (PBCs), to conduct a microstructural homogenisation to determine the effective mechanical properties of an expanded foil LSP layer. More recently, Wang et al.⁴⁷ developed RVE models to calculate material properties of NCF/CIPCF film using a RSA algorithm. Chippendale and Golosnoy⁴⁸ used an RVE model to determine the effect of percolation phenomena and the percolation threshold on the bulk electrical conductivity of CFRP.

Background summary and paper aim

In summary, all existing lightning strike simulation research has largely focussed on a singular modelling strategy - homogenised mesoscale modelling to predict the overall lightning damage in composite laminates in the presence or absence of an outer LSP layer. However, this modelling approach is not capable of capturing unique damage mechanisms of each constitute (i.e., matrix and carbon fibre). Instead, RVE-based micromechanics models can be developed to understand lightning damage mechanisms at a smaller length scale. As previously mentioned, a few RVE studies were performed within the field of lightning strike research but were limited to the calculation of effective properties for composite plies and LSP layers. Understanding lightning-induced damage in composites at a microscale level is still lacking. It is still not known how each constitute is damaged due to lighting or what the dominant damage mechanisms are. In addition, it is not known how lightning waveform properties (i.e., peak current) influence these damage mechanisms or where the weakest (critical) location for damage initiation is. Therefore, the motivation of the current study is to address these unknowns.

This paper proposes, for the first time, micromechanics-based 3D RVE models to characterize the progressive failure of CFRP composites under thermal, mechanical, and coupled thermo-mechanical loading, consistent with actual lightning direct effects. Microscopic fibre damage, matrix damage, and fibre-matrix interface failure were predicted from the RVE models, consisting of randomly distributed carbon fibres in an epoxy matrix. Various electric potential boundary conditions were imposed to understand the effects of grounding on microscopic thermal damage in RVE models. The effects of fibre volume fraction and peak current amplitude on microscopic damage evolution were also thoroughly studied. In this work, the matrix failure was predicted using a ductile plasticity model with Drucker–Prager (DP) yield criterion, and the fibre-matrix interface debonding was characterized using cohesive surfaces with a bilinear traction-separation law. This work has attempted to propose an efficient and robust modelling framework to understand localised damage in CFRP composites caused by lightning, at the microscale level. This work is fundamentally different from the typical mesoscale modelling used to predict the overall lightning damage in composite laminates. Rather, this work focusses on predicting the microscopic lightning damage to understand complex lightning damage mechanisms near the lightning attachment.

Methodology

Representative volume element modelling

The 3D RVE model of a UD composite laminate included parallel, circular cross-sectioned carbon fibres randomly distributed in an epoxy matrix. The RVE models with three fibre volume fractions (V_f = 55%, 60%, and 65%) were developed with a total volume of 50 μm × 50 μm × 50 μm. Fibres had diameters of 5.3–7 μm and were randomly generated and positioned such that the distance between neighbouring fibres was ≥0.035d_f, in line with Ref. 45, to ensure an adequate discretisation of this region. A periodic microstructure was maintained by splitting those fibres which intersected the RVE boundary and copying the overlapping portion of the fibre to the opposite side or corner of the RVE model, in line with Refs. 44,45. Note that a periodic microstructure differs from PBC connecting two opposite RVE edges by constraints. The resulting RVE models are shown in Figure 1. No designated interface region was modelled in the RVE models. Rather, cohesive surfaces were applied at fibre-matrix interfaces to understand interfacial debonding.

Figure 1.

Representative RVE models with V_f = 55%, 60%, and 65%.

Modelling of lightning strike effects on the RVE model was completed in two stages. The first stage was a thermal-electric model predicting the temperature distribution in the RVE model primarily due to Joule heating. The second step predicted mechanical damage resulting from temperature, thermal expansion and/or pressure loading. This two-step sequential modeling procedure enables the thermal and mechanical damage to be decoupled from the total (combined) lightning damage in the RVE model. However, it also allows temperature, thermal expansion and pressure loading to be applied together.

To study both the thermal and mechanical damage within the RVE models, there was a need to assemble an appropriate set of simulation methods based on existing literature to predict damage with different boundary and loading conditions. Therefore, an analysis workflow was assembled for the work herein, as outlined in Figure 2, which was adapted from Refs. 26,49,50. The time-dependant temperature profile predicted in a thermal-electrical model was transferred to a thermo-mechanical model using python scripts and the method developed by Millen et al.,^26,49,50 in which a temperature profile is recorded and transferred for each node from one model to the next. This approach allowed the electro-thermal coupling effect to be replicated within the thermo-mechanical RVE simulation through the transfer of the temperature profile. Thus, the temperature prediction of each thermal-electric simulation was passed to an equivalent thermo-mechanical simulation. Therefore, the time dependent variation in temperature is represented throughout the duration of the thermo-mechanical RVE simulation when predicting specimen damage. Subsequently, two sequential simulations could be completed for each RVE model with time periods of 30 μs for both analyses. Both model types will be now discussed in more detail.

Figure 2.

Simulation flowchart.

Thermal-electric representative volume element model

In the thermal-electric model, electric current loading was applied to the top surface of the RVE shown in Figure 3(a). The arc radius was assumed to be constant given the very small size of the RVE model. Therefore, electrical loading (i.e., surface current) was applied uniformly to the entire top surface. Zero electrical potential boundary conditions, representing the grounding terminal, were imposed on the bottom surface of the baseline RVE model (Figure 3(a)). PBC were not considered in this work. Note that PBC are typically used to calculate the effective properties of a RVE model and avoid potential issues associated with boundary effects. This work intends for microscopic characterisation of lightning-induced damage in a RVE model, thus PBC are not required. While the recommended Waveform A is 200 kA (T₁/T₂ = 6.4/69 µs),⁵¹ the typical laboratory and simulation lightning test Waveform A (T₁/T₂ = 4/20 µs)^13,51 with three peak current amplitudes of 10 kA, 20 kA and 40 kA and a time period of 30 µs was applied to each RVE model, as shown in Figure 3(b).

Figure 3.

Representative volume element boundary conditions and current waveforms used in this work from Refs. 13,52.

Electrical current flow strongly depends on both the optimal conduction path and grounding conditions, i.e., current preferentially flows through the least resistance path in a material to the grounding terminal. Thus, careful consideration of electrical potential was made in this work. In addition to the baseline RVE model (Figure 3(a)), additional parametric studies were performed to study the effects of grounding conditions on microscopic lightning thermal damage in RVE models (V_f = 55%) subjected to 40 kA peak current load. Three different sets of ground boundary conditions were considered herein and are shown in Figure 4, including (i) bottom surface only, (ii) bottom surface and two side surfaces (where two side surfaces corresponds to the surfaces parallel to the fibre direction), and (iii) bottom surface and four side surfaces.

Figure 4.

Grounding conditions used in thermal-electric representative volume element models.

In the thermal-electric models, the epoxy matrix was assumed to be isotropic, while the carbon fibres were assumed to be transversely isotropic. Thermal and electrical material properties, at room temperature, for the fibres and matrix are shown in Table 1. The epoxy matrix was assumed to be a very poor conductor at room temperature. The value of room temperature electrical conductivity of the epoxy matrix used herein was derived from the value used for through-thickness conductivity in mesoscale models (in “resin-rich” regions) in preceding research.^2,4,19 As the temperature of the epoxy matrix increases, so does the conductivity before reaching a char state with a conductivity of 1 × 10⁵ 1/Ω.mm.² Thermal and electrical conduction between the fibres and epoxy matrix was included using the surface-to-surface contact functionality with *Gap Conductance and *Gap Electrical Conductance in ABAQUS.

Table 1.

Thermal-electric material properties.

Material	Density (kg/mm³)	Specific heat (J/kg^oC)	Thermal conductivity (W/mm.K)			Electrical conductivity (1/Ω.mm)
Material	Density (kg/mm³)	Specific heat (J/kg^oC)	k₁₁	k₂₂	k₃₃	σ₁₁	σ₂₂	σ₃₃
Fibre	1.52 × 10⁻⁶	1065	0.008	0.00067	0.00067	35.97	0.001145	0.001145
Matrix	1.33 × 10⁻⁶	1065	0.00067			1.79 × 10⁻⁶

In Table 1 the 1-direction is the fibre axial direction, while the 2- and 3-directions correspond the radial direction.

After an initial mesh sensitivity study, the RVE models were meshed using 3D linear coupled thermal-electrical solid (DC3D8E) elements with an approximate element size of 1.5 μm. A total of between 60,000 and 70,000 elements were used to balance the result accuracy and computational cost. A full breakdown of the elements used in the thermal-electric RVE models are given in Table 2.

Table 2.

Mesh details for representative volume element models.

V_f (%)	Total elements	Matrix elements	Fibre elements
55	62,370	25,212	37,158
60	65,373	25,410	39,963
65	74,844	27,984	46,860

Thermo-mechanical model

Thermo-mechanical analysis was completed with a dynamic, temperature-displacement, explicit analysis step in ABAQUS/Explicit using 3D trilinear displacement and temperature, reduced integration elements, with hourglass control (C3D8RT).

Pressure loading was applied to the same specimen surface as the current load previously. To determine the equivalent mechanical (shockwave) loading from the peak current waveform a conversion procedure was required. In practice, the magnitude of mechanical loading associated with high current/voltage impulse waveforms (consistent with natural lightning), can be calculated either from a classical arc channel expansion theory or by solving arc channel energy balance equations.⁵³ Each method requires several thermo-physical properties of ionised air and arc channel/spark energy, thus mechanical loading calculation is not straightforward. Lee et al.^4,54 proposed a simple conversion procedure from a peak current amplitude to equivalent mechanical loading, valid for lower (<50 kA) peak current amplitude waveforms, and provided a conversion table. Figure 5(a) shows equivalent lightning mechanical loading as a function of the peak amplitude of a current impulse waveform. In this work, the magnitudes of mechanical loading corresponding to peak current amplitudes of 10 kA, 20 kA and 40 kA, were calculated as 2.8 MPa, 4.1 MPa, and 6.7 MPa, respectively. Figure 5(b) shows the dynamic lightning mechanical loading which was assumed to follow identical temporal characteristics of the current waveforms (Figure 3(b)).

Figure 5.

Graphs of conversion from peak current amplitude to lightning mechanical pressure.

Constitutive material models were developed to predict fibre and matrix damage at the microscale level. Table 3 provides the temperature-dependent elastic properties of carbon fibre and epoxy matrix constituents used in the thermo-mechanical models. The fibres were assumed to be linearly elastic. Pressure loading calculated up to 40 kA peak current and associated temperature increases were not significant, thus no damage criteria was applied to the fibres.

Table 3.

Elastic properties of carbon fibre and epoxy matrix.

Material	Temperature (^oC)	Elastic moduli (GPa)		Poisson’s ratio		Shear moduli (GPa)
Material	Temperature (^oC)	$E_{11}$	$E_{22} = E_{33}$	$ν_{12}$	$ν_{23}$	$G_{12}$	$G_{23}$
Fibre	25	230	186	0.255	0.405	20.5	7
	190	230	186	0.255	0.405	20.5	7
	300	230	186	0.255	0.405	20.5	7
	500	230	186	0.255	0.405	20.5	7
	3316	230	186	0.255	0.405	20.5	7
	>3316	0.23	0.046	0.025	0.405	5.125	2.2
Matrix	25	3080		0.35
	190	1105		0.35
	300	1105		0.35
	500	1105		0.35

The subscript 1, 2, and 3 each refers to the longitudinal (or fibre), transverse, and thickness direction.

The epoxy matrix was modelled as isotropic and elastoplastic with a Drucker-Prager (DP) plasticity model coupled with progressive damage. A DP model was originally developed for geopolymer materials (i.e., soil, concrete) based on the Morh-Coulomb theory, but also widely used for modelling the failure of a polymer matrix with yield (hardening) behaviour. In this work, the ABAQUS built-in, modified DP model⁵⁵ was used to capture matrix failure. The DP yield function F_s is defined as

F_{s} = t - p \tan β - d = 0

(1)

with

t = \frac{1}{2} σ_{e q} [1 + \frac{1}{K} - (1 - \frac{1}{K}) {(\frac{J_{3}}{q})}^{3}]

(2)

where

σ_{e q}

is the von-Mises equivalent stress;

K

is the flow stress ratio;

J_{3}

is the third invariant of the deviatoric stress tensor. Three material parameters, which can be determined from uniaxial tests, were required to capture the yield and hardening criteria and the flow rule: (i) angle of friction

β

, (ii) flow stress ratio

K

, and (iii) dilation angle

ϕ

β

indicates the slope of the linear yield surface in the meridional (p–t) plane representing a sensitivity of material’s yield stress to hydrostatic pressure;

K

is a ratio of the flow stress in triaxial tension to compression;

ϕ

is a suitable parameter describing inelastic volume changes, i.e., a positive angle indicates plastic volume expansion. In this work,

β

= 0.137,

ϕ

= 0 assuming incompressible flow, and

K

= 0.954 (for typical epoxy matrices) were used.^56,57 Note that

β

and

K

govern material yielding and hardening, while

ϕ

determines plastic flow rule.

The evolution of DP yield surface was assumed to be controlled by the yield stress under uniaxial compression. The summary of the compressive yield stress and absolute plastic strain used in this work are tabulated in the following: [(30.095, 0), (61.93, 0.03), (61.93, 0.25)]. Hooputra et al.⁵⁸ proposed a continuum ductile damage model predicting the ductile nature of noticeable plastic deformation of an epoxy matrix prior to damage initiation. Two material parameters are required in a ductile damage model to predict the onset of damage ${\bar{ε}}_{D}^{p l}$ : stress triaxiality, $η$ and equivalent plastic strain, ${\bar{ε}}_{0}^{p l}$ . $η$ is a ratio of hydrostatic stress $σ_{m}$ to von-Mises equivalent stress $σ_{e q}$ and found to be 0.333 under uniaxial loading.⁵⁹ At the onset of damage, ${\bar{ε}}_{0}^{p l}$ reaches a threshold value, often known as fracture strain. Assuming strain-rate independent epoxy matrix, the ductile damage initiation was implemented with $η$ = 0.333 and ${\bar{ε}}_{0}^{p l}$ = 0.0875. A linear damage evolution law with a fracture energy of $G_{f}$ = 10 J/m²⁶⁰ was also used. More details on damage initiation and evolution criteria for ductile materials are available in reference 55.

The interfaces between fibres and matrix were modelled using cohesive surfaces incorporating the bi-linear traction-separation law, as in previous work.²⁶ Penalty interface stiffness values

K

, which controlled the linear elastic behaviour of the cohesive interfaces, were left default.⁵⁵ In this work, the onset of fibre/matrix interfacial damage was governed by the quadratic stress criterion and damage occurred when the sum of nominal stress ratios reached 1, as follows:

{(\frac{〈 t_{n} 〉}{t_{n}^{0}})}^{2} + {(\frac{t_{s}}{t_{s}^{0}})}^{2} + {(\frac{t_{t}}{t_{t}^{0}})}^{2} = 1

(3)

where

t_{n}^{0}

t_{s}^{0}

and

t_{t}^{0}

are the strengths (or peak contact stress) in each orientation, normal

(n)

and two shear (

s

and

t

) directions. The symbol

〈 〉

represents the Macaulay bracket. Dissipated fracture energy is then used to determine the separation displacement at failure. The Benzeggagh and Kenane (B-K) criterion^22,61 was adopted for fracture energy dissipation during damage propagation. The criterion reads

G^{C} = G_{n}^{C} + (G_{s}^{C} - G_{n}^{C}) {(\frac{G_{s} + G_{t}}{G_{n} G_{s} + G_{t}})}^{η_{B K}^{C}}

(4)

where

η_{B K}^{C}

is the mixed-mode interaction and

G_{n}^{C}

and

G_{s}^{C}

are the critical fracture energies required to cause failure in the normal and shear directions, respectively. A relationship for the temperature dependence of the interface has recently been proposed⁶² and applied within mesoscale models.^26,49 This same relationship was used herein. Room temperature properties were adapted from literature based on established assumptions for fibre-matrix interface failure in both normal and shear modes.^44,63–65 Therefore, in this work, fracture energies of 2 J/m² and 100 J/m² were used for normal and shear failure modes, respectively, with further assumed temperature dependent properties shown in Table 4.

Table 4.

Interface mechanical properties.^44,62–65

Temperature (^oC)	$t_{n}^{0}$ (MPa)	$t_{s}^{0}$ = $t_{t}^{0}$ (MPa)	$G_{n}^{C}$ (J/m²)	$G_{s}^{C}$ = $G_{t}^{C}$ (J/m²)
25	58.0	92.0	2	100
300	1 × 10⁻⁴	1 × 10⁻⁵	1 × 10⁻⁶	1 × 10⁻⁷
3000	1 × 10⁻⁴	1 × 10⁻⁵	1 × 10⁻⁶	1 × 10⁻⁷

The thermo-mechanical models also incorporated thermal expansion of each constituent in a RVE model. The coefficients of thermal expansion (CTE) of the carbon fibres and epoxy matrix, taken from literature,¹⁹ are shown in Table 5.

Table 5.

Thermal expansion material properties from Ref. 19

Temperature (^oC)	Carbon fibre		Epoxy matrix
Temperature (^oC)	α₁₁ (×10⁻⁸)	α₂₂ = α₃₃ (×10⁻⁵)	α (×10⁻⁵)
25	1.80	2.16	2.16
200	5.40	3.78	3.78
260	5.40	3.78	3.78
600	5.40	3.78	3.78
3316	5.40	3.78	3.78
>3316	5.40	3.78	3.78

Model summary

This modelling approach and its derivation is presented in further detail in Figure 6. The typical mesoscale model of a lightning strike on CFRP considers a specimen measuring 150 mm × 100 mm with a typical ply thickness ranging from 0.1-0.15 mm.^2,16 The typical lightning arc is modelled with a diameter of 10 mm,¹⁶ shown in Figure 6(a). The typical size of an RVE modelled in literature is around 50 μm × 50 μm × 50 μm with fibre diameters of 5.3–7 μm.^44,45 Therefore, the models herein use this size of RVE and assume it is located at the top surface of the top ply in a typical laminate, shown in Figure 6(b). As such, the typical 10 mm lightning arc used in mesoscale lightning strike models will completely cover the top surface of the RVE, shown in Figure 6(c). Therefore, rather than applying the entire lightning waveform to the 0.0025 mm² area at the top of the RVE, a scaling factor is used to account for the difference in area between the typical 10 mm arc and the size of the RVE. As a result, only the fraction of the current that would flow into the RVE is applied in each case. If this approach was not used, the applied current density and associated Joule heating would be significantly exaggerated and therefore exaggerate the predicted damage mechanisms. In fact, the actual gradients in the current density and initial arc channel radius are not known since they are strongly influenced by various factors, including the peak current amplitude, temporal characteristics, surrounding air conditions, and arc channel expansion velocity.^1,54

Figure 6.

Summary of representative volume element modelling approach and derivation of boundary conditions and thermo-mechanical loading from typical meso-scale models.

In addition, the boundary conditions studied herein have been chosen in order to capture the actual behaviour within the ply of a mesoscale level analysis and translate this to the RVE. When the current enters the ply (and the RVE) it can flow in three possible directions, along the fibres, transversely or through-thickness. This depends on the potential difference between the attachment (current injection location) and the grounding terminal, and electrical conductivity in each orthogonal direction. Each situation is therefore considered before choosing the most appropriate boundary conditions for further study.

While the typical mesoscale models of a lightning strike consider the specimens with multiple plies and the potential delamination damage between the neighbouring plies, this is beyond the scope of the models herein. Rather, the proposed RVE models can predict the fibre-matrix interfacial debonding in a small unit of a ply.

Lightning produces thermo-mechanical damage in composite specimens. Mechanical damage is highly localised and confined inside a lightning arc attachment. In addition, the mechanical damage area is fully surrounded by the thermal damage area, thus it is not straightforward to decouple pure mechanical damage from overall lightning damage since thermal, electrical, and mechanical damage occur simultaneously within the arc attachment area.⁴ Therefore, in order to determine the contribution of different physics, pressure loading and thermal expansion, these were applied separately and then together. Therefore, three thermo-mechanical RVE models were completed with pressure loading only, temperature loading (associated with thermal expansion) only, and combined pressure-temperature loading. All simulations were completed on the Kelvin-2 High Performance Computing (HPC) cluster at Queen’s University Belfast.

Thermal damage results

In this section, coupled thermal-electric model results are discussed by focusing on the temperature distribution predicted within the RVE model. Typical epoxy matrix degradation occurs within the temperature range of 300–800°C, while carbon fibre damage occurs at temperatures greater than 3000°C. Therefore, these temperature ranges were used to characterise lightning-induced thermal damage within the RVE constituents. In general, since the carbon fibre damage region is totally encompassed by the widespread matrix damage region, lightning thermal damage in composites often corresponds to the extent of the matrix damage region. Along with this definition, lightning thermal damage calculated in an RVE model was estimated from the volume of the epoxy matrix region with temperatures ≥300°C or ≥500°C.

Effects of grounding condition

Three different electric potential boundary conditions were compared using a RVE model (V_f = 55%) subjected to a 40 kA peak current load. Figure 7 shows predicted temperature distributions, indicative of lightning thermal damage in the epoxy matrix, within the RVE models when zero electric potential ground boundary conditions were applied to different RVE surfaces. Note that zero electric potential indicates a grounding terminal. As mentioned earlier, thermal decomposition occurs in typical epoxy matrices in the temperature range of 300 – 800°C. In this thermal decomposition temperature range, the epoxy matrix is progressively and irreversibly damaged due to pyrolytic matrix decomposition.⁶⁶ It can be seen that the addition of more grounding terminals reduced the predicted epoxy matrix damage significantly, limiting the peak temperature within the RVE model. This is reasonable since more grounding terminals allow electric current to exit, thus guaranteeing better current flow distribution (less Joule heating) over the RVE model.

Figure 7.

The effect of grounding conditions on representative volume element thermal damage prediction at t = 30 µs.

However, the most representative and interesting case is the RVE model with only the bottom surface grounded (Figure 7(a)). This RVE boundary condition can be assumed to represent the centre of the top ply on a large specimen under the direct influence of the lightning arc. Since even the smallest lightning arc channel (radius R ≥5 mm) is much larger than the RVE length (50 μm), the current likely tends to travel through-the-thickness of the RVE model rather than in-plane directions due to large potential difference, equivalent to voltage. As a result, only the bottom surface of each RVE will be grounded in the following simulations.

Initial assessment of the temperature profile in the carbon fibres of the RVE model with bottom surface grounding is based on Figure 7(a). Predicted temperature increases in carbon fibres are well below their sublimation temperature (∼3000°C), indicating no thermal damage. Note that the predicted maximum fibre temperature was 788°C in Figure 8. This can be explained by their comparatively large electrical conductivity compared to the epoxy matrix (Table 1). Therefore, further discussion of lightning thermal damage is primarily focused on the epoxy matrix regions of the RVE model.

Figure 8.

Fibre temperature distribution under 40 kA peak current amplitude at t = 30 µs.

Effects of peak current amplitude

Figure 9 shows the predicted temperature distribution resulting from 10 kA, 20 kA and 40 kA peak current waveforms on the RVE models (V_f = 55%) and corresponding epoxy matrix regions. In all cases considered, the most severe damage occurs at the top surface of the RVE model, where a uniform electric current is applied. In general, the carbon fibres are largely unaffected by the current load, however, the fibres which intersect with the loading surface do see a modest temperature increase, relative to their sublimation temperature. As the applied peak current increases, the extent of thermal damage onset (specifically damage to the matrix through Joule heating where the range 200–400°C represents moderate damage and 600°C and above represents matrix ablation), taken as the volume of the epoxy matrix with temperatures greater than 300°C, was nearly constant under 10 kA and 20 kA peak amplitudes, i.e., 1057 μm³ (10 kA) and 976 μm³ (20 kA), both around 2% of the total epoxy matrix volume. However, thermal damage onset under 40 kA peak amplitude increased significantly to 3920 μm³, corresponding to around 7% of the initial epoxy matrix volume. The extent of more severe thermal damage, taken as the volume of the epoxy matrix region with temperatures greater than 500°C increased steadily with increasing peak current from 545 μm³ at 10 kA to 2211 μm³ at 40 kA.

Figure 9.

The effect of a peak current amplitude on representative volume element thermal damage prediction at t = 30 µs.

Figure 10 shows the primary current flow (or conduction) paths through each RVE model with varying peak current amplitudes (10–40 kA), established by inferring current flow through temperature contours, due to the Joule heating effect. It can be seen that the nearly identical current flow paths are produced in all cases. However, the amount of current flowing in each conduction path varies with peak current amplitude. This current was calculated by the sum of all nodes within each conduction path on the bottom surface of the RVE, shown on Figure 10(c). For the 55%-10 kA case (Figure 10(a)), 50% of the current leaves through the central flow path, Path 1, while this increases to 68% for the 55%-20 kA case (Figure 10(b)) and 80% for the 55%-20 kA case (Figure 10(c)). This suggests that a primary conduction path is mainly dictated by the fibre distribution, but its magnitude is strongly controlled by the incident peak current amplitude. This hypothesis is further explored when varying V_f in the following section.

Figure 10.

Primary conduction paths through representative volume element models with varying peak current amplitudes at t = 15 µs.

Figure 11 shows the development of the primary (red) and secondary (black) conduction paths, i.e. the larger and smaller of the main current flow paths, through the RVE model (V_f = 55%, 40 kA) at three time points, t = 6 μs, 12 μs and 30 μs, respectively. These times were chosen as they correspond to 20% (just after $T_{1}$ of the waveform), 40% and 100% of the total waveform time period. These images show that even at t = 6 μs, the primary conduction paths have been determined and the magnitude of current through these paths steadily increases through the simulation. The secondary condition path is formed at t = 6 μs due to increasing loading across the width of the top surface where current is applied.

Figure 11.

Primary and secondary conduction paths in the representative volume element model (V_f = 55%, 40 kA).

Effects of fibre volume fraction

Figure 12 shows the temperature distribution predicted from 40 kA peak current on the RVE models with V_f = 55%, 60%, and 65%, respectively. As observed previously, the carbon fibres are largely unaffected, even at higher V_f. As V_f increased, the extent of thermal damage onset in the RVE models increased from 3920 μm³ (55%-40 kA case, Figure 12(a)) to 5121 μm³ (60%-40 kA case, Figure 12(b)) to 5918 μm³ (65%-40 kA case, Figure 12(c)). These thermal damage volumes correspond to 7%, 10%, and 14% of the pristine (undamaged) epoxy matrix volumes, respectively. However, the volume of more severe damage (taken as >500°C) decreased with increasing V_f from 2211 μm³ (55%-40 kA) to 285 μm³ (65%-40 kA). The predicted peak temperature within each RVE model decreased with increasing V_f from 788°C (Figure 12(a)) to 590°C (Figure 12(c)). This makes sense since more carbon fibres spread electric current more rapidly, causing less Joule heating within the epoxy matrix and, therefore less temperature increase in each RVE model. When the threshold temperature was 300°C, predicted thermal damage increased steadily with increasing V_f. In contrast, with a threshold temperature of 500°C, predicted thermal damage decreased with V_f. Therefore, this somewhat counterintuitive behaviour may be a result of the chosen threshold temperatures.

Figure 12.

The effect of V_f on representative volume element thermal damage prediction at t = 30 µs.

Figure 13 shows the predicted current conduction paths, from the RVE models with V_f = 55%, 60%, and 65%. While the nearly identical current flow path was predicted at a fixed V_f (Figure 10), modifying V_f produced significant changes in the flow of current. For example, when V_f = 55%, the flow of current was divided over two paths (Figure 13(a)). The rightmost flow path took around 37% of the current, while the central path took around 63% of the current. Increasing V_f to 60% produced a more uniform current spread over three broad columns of fibres (Figure 13(b)). Finally, increasing V_f to 65% spread the current over the majority of the RVE model with one localised region of high current density down the approximate centre of the RVE model (Figure 13(c)), accounting for around 24% of the current flow. The nature of this primary conduction path formation is highly influenced by the distance between each fibre, associated with fibre diameter characteristics (i.e., mean and standard deviation).

Figure 13.

Primary conduction paths generated in the representative volume element models with varying V_f at t = 15 µs.

Thermo-mechanical damage results

This work proposed a sequentially coupled thermo-mechanical modelling framework to predict lightning thermo-mechanical damage in a RVE model. Heat transfer primarily takes place via conduction within a RVE model and requires some time (likely a few seconds) to achieve thermal equilibrium. Thermal damage, resulting from Joule heating due to the application of current, may be accelerated due to thermal strain - the product of the coefficient of thermal-expansion (CTE) and a temperature change. In addition, lightning damage may become worse with dynamic mechanical loading (i.e., lightning shockwave overpressure). The following sections highlight thermo-mechanical damage prediction in the RVE model each obtained from dynamic pressure loading, temperature loading (due to thermal strain), and its combined effects.

Pressure loading only

Figure 14 shows the stress generated in the matrix of the 65%-40 kA RVE model under pressure loading only. Note that a 40 kA peak current amplitude is equivalent to 6.7 MPa lightning shockwave overpressure (Figure 5(a)). It can be seen in Figure 14 that the magnitude of stress reached a maximum of 2.7 MPa, well below that which would cause both carbon fibre and epoxy matrix damage. This indicates that dynamic pressure loading equivalent to a 40 kA peak current amplitude does not induce severe mechanical stress within the RVE. This is consistent with work done by Refs 4,24,26,49,54, who despite using much larger specimens, e.g. 150 mm × 100 mm, also found negligible damage from lightning mechanical loading alone. This is an important finding as lightning mechanical damage does not influence the overall lightning damage in a RVE model, so there is no need to decouple it from lightning damage assessment.

Figure 14.

Stress predicted in the matrix region of the representative volume element model (V_f = 65%) under only pressure loading equivalent to 40 kA peak current at t = 30 µs.

Temperature loading only

The temperature loading defined herein indicates the temperature profile transferred from the preceding thermal-electric analysis. The effects of thermal strains (associated with thermal expansion during thermal equilibrium) on thermo-mechanical damage can be evaluated in an RVE model. Results for temperature loading are more interesting and severe than pressure loading only. Figure 15 shows the predicted equivalent plastic strains (PEEQs) contours of the 55%-40 kA and 65%-40 kA RVE models. The RVE deformation and corresponding PEEQ strains are higher for the 55%-40 kA RVE model than the 65%-40 kA RVE model due to a large volume fraction of relatively “weak” matrix (i.e., the epoxy matrix is mechanically weaker than the carbon fibres). In addition, the CTE mismatch between carbon fibres and epoxy matrix contributes to significant RVE deformation, i.e., 5.40 × 10⁻⁸ (carbon fibre) versus 3.78 × 10⁻⁵ (epoxy matrix) at 260°C from Table 5. The 55%-40 kA RVE model exhibits a somewhat continuous strain path developed within the narrow regions between fibres (highlighted in green in Figure 15(a)).

Figure 15.

Predicted PEEQ distributions on the 55%-40 kA and 65%-40 kA representative volume element models at t = 30 µs.

Figure 16 shows the predicted ductile matrix damage for both the 55%-40 kA and 65%-40 kA RVE models. Comparing Figure 16 and Figure 12, it can be seen that the ductile matrix damage closely matches the thermal damage profile, resulting from Joule heating, predicted in the corresponding thermal-electric models, particularly for the 55%-40 kA model. As with PEEQ, there was significantly more ductile damage in the 55%-40 kA RVE model due to a large volume of (relatively weak) matrix-rich areas. An amount of ductile matrix damage may be associated with the mean and standard deviation (STD) of fibre-to-fibre distances. This is greatly dependent on the effect of fibre positioning and variability in the fibre-to-fibre distance associated with V_f. For example, looking at the widest region of ductile damage on the 55%-40 kA RVE model, bottom right of Figure 16(a), the average fibre-to-fibre distance was 2.1 μm while in the same region of the 65%-40 kA RVE model the average fibre-to-fibre distance was 0.8 μm.

Figure 16.

Predicted ductile matrix damage in (a) 55%-40 kA and (b) 65%-40 kA representative volume element models at t = 30 µs.

Comparing Figure 13 with Figure 16, the primary conduction paths generally match the severe damaged matrix regions (in red) with relatively large fibre-to-fibre distances. More investigation is required to further understand the effects of fibre distribution on subsequent ductile matrix damage development.

As a quantitative measure of overall thermo-mechanical damage in the RVE models, the total damage volume was calculated based on the value of the ductile damage initiation criteria (DUCTCRT), i.e., ≥0.99. The resulting matrix damage for the 55%-40 kA RVE model was 17,523 μm³ (14% of the full RVE model and 31% of the matrix region) and 4514 μm³ (4% of the full RVE model and 8% of the matrix region) for the 65%-40 kA RVE model.

Figure 17 shows the stress distributions predicted within the matrix and fibre regions for both 55%-40 kA and 65%-40 kA RVE models. In both cases, the stresses within the matrix were relatively low, only with mild stress concentration. In the 55%-40 kA RVE model, the peak fibre stress occurred in the central fibre with largest current flow (63% of the current), as highlighted earlier. In the 65%-40 kA RVE model, fibre stress was much more uniform across all fibres. In both cases, the areas of matrix with high stress closely matched the regions with severe matrix damage (Figure 16). Note that the RVE model presented herein consisted of randomly distributed carbon fibres within the epoxy matrix. Thus, the regions of relatively high stress concentration can be varied depending on fibre arrangement and volume fraction.

Figure 17.

Predicted stresses within the matrix and fibres for both 55%-40 kA and 65%-40 kA representative volume element models at t = 30 µs.

Figure 18 shows the predicted fibre-matrix interfacial debonding in both the 55%-40 kA and 65%-40 kA RVE models. As with other damage modes, the predicted fibre-matrix interfacial debonding closely matched the thermal damage profile (Figure 12) and primary current flow paths (Figure 13). While interfacial damage generally followed the thermal damage profile, this damage mode was somewhat independent of the ductile matrix damage profile. This is likely due to the temperature dependence of the interlaminar properties, shown in Table 4. As with other damage modes, interfacial damage varies with V_f. The interfacial damage, taken as the area of the cohesive surface damage criteria (CSDMG) ≥0.99, was 274 μm² for the 55%-40 kA RVE model, but only 72 μm² for the 65%-40 kA RVE model.

Figure 18.

Predicted fibre-matrix interfacial damage in the 55%-40 kA and 65%-40 kA representative volume element models at t = 30 µs.

Debonding between the carbon fibres and epoxy matrix was also observed within the 55%-40 kA RVE model subjected to thermal loading only. Figure 19(a) shows the resulting fibre-matrix debonding in the whole RVE and in a zoomed-in image. In both cases elements with DUCTCRT ≥0.99 have been removed since fibre-matrix debonding can occur due to a combination of matrix damage and interfacial damage. As shown in Figure 19(a), the RVE model can predict considerable fibre-matrix interfacial debonding, consistent with the scanning electron microscope (SEM) images (Figure 19(b) and (c)) of representative carbon/epoxy laminates subjected to a 38 kA peak current amplitude.⁶⁷ This is important for evaluating various damage modes in the lightning attachment location. The RVE model proposed in this work, to the best of the authors’ knowledge, is the first model predicting lightning-induced fibre-matrix interfacial damage.

Figure 19.

Interfacial fibre-matrix debonding prediction at t = 30 µs compared with SEM images from Ref. 67.

Combined pressure-temperature loading

As expected, results for combined pressure-temperature loading show a small difference to the results with temperature loading alone. The matrix damage volume for the 55%-40 kA RVE model with temperature loading only was 17,523 μm³, while this increased slightly to 18,483 μm³ (<5%) with the addition of pressure loading. This finding again confirms a negligible effect of lightning mechanical loading on overall damage assessment.

Results summary

In a CFRP composite, the polymer matrix undergoes severe thermal decomposition when subjected to a lightning strike. The matrix decomposition (occurring due to pyrolysis) results in irreversible thermal damage. The extent of irreversible matrix damage, indicative of lightning thermal damage in a CFRP composite, strongly depends on a thermal decomposition temperature range (300–800°C, herein) and a heating rate. At the lightning attachment location, the rate of temperature increase may locally exceed 10,000°C/min⁷ and the onset and evolution of irreversible thermal damage may require a different decomposition temperature range to assess more accurate lightning thermal damage. Similarly, a quantitative measure of thermo-mechanical damage may depend on a threshold value of DUCTCRT.

Table 6 provides a full summary of the predicted thermal damage, peak temperature, and ductile matrix damage, indicative of thermo-mechanical damage in all the RVE models studied herein. In the table, thermal damage was obtained from the volume of elements in the epoxy matrix region with temperatures ≥300 or 500°C after coupled thermal-electric analysis while ductile matrix damage was calculated from elements in the epoxy matrix region where DUCTCRT ≥0.70 or 0.99 after thermo-mechanical analysis. When the threshold temperature was 300°C, predicted thermal damage increased steadily with increasing peak current amplitude (with a fixed V_f of 55%) and V_f (with a fixed peak current of 40 kA). In contrast, with a threshold temperature of 500°C, predicted thermal damage increased with increasing peak current amplitude (at V_f = 55%), but decreased with V_f. This is reasonable since a higher peak current produces more thermal damage and a higher fibre volume fraction produces more stiff and strong carbon fibres, indicating less relatively weak epoxy matrix in the RVE model, contributing to less thermal damage. This suggests that using a lower threshold temperature may overestimate the overall thermal damage in the RVE model, which makes drawing general conclusions challenging. As expected, peak temperature calculated from each RVE model increased with increasing peak current amplitude and decreased with increasing V_f, regardless of the chosen threshold temperature. Predicted ductile matrix damage after thermo-mechanical analyses decreased with increasing peak current amplitude. As expected, a lower DUCTCRT threshold value increased the ductile matrix damage volume in the RVE model. Another important finding is that, from Table 6, thermal loading associated with thermal strain/expansion has the largest contribution to thermo-mechanical damage. Pressure loading, even applied together with thermal loading, had a negligible effect (<5% difference) on predicted overall thermo-mechanical damage. This matches the previous conclusions at the mesoscale level of analysis.²⁶

Table 6.

Summary of thermal damage, peak temperature, and mechanical damage for representative volume elements.

RVE Model	Thermal-electric analyses			Thermo-mechanical analyses
	Thermal damage (μm³)		Peak Temp. (^oC)	Ductile matrix damage (μm³)
	>300°C	>500°C		Pressure loading	Temp Loading		Combined pressure and temp. Loading
	>300°C	>500°C		Pressure loading	>0.70	>0.99	>0.70	>0.99
55%-10 kA	1057	545	736	—	—	—	—	—
55%-20 kA	976	1540	754	—	—	—	—	—
55%-40 kA	3920	2211	788	0	27,468	17,523	28,336	18,483
60%-40 kA	5121	1401	685	—	—	—	—	—
65%-40 kA	5918	285	590	—	9078	4514	—	—

Note: thermal damage was obtained from the volume of elements in the epoxy matrix region with temperatures ≥300 or 500°C after coupled thermal-electric analyses, and ductile matrix damage was calculated from elements in the epoxy matrix region where DUCTCRT ≥0.70 or 0.99 after thermo-mechanical analyses.

Limitations and future work

The proposed RVE-based models presented previously have, for the first time, been developed to characterize complex lightning damage mechanisms at a fundamental level, i.e., in the fibre, matrix and fibre-matrix interface. In contrast to macroscale models where all loading/boundary conditions are directly monitored by actual lightning tests, our microscale models were developed with some assumptions since these were unknown at a fundamental level.

It was assumed that the loading and grounding of the microscale models could be derived from the macroscale/mesoscale models widely used in literature, with the approach in this paper shown in Figure 6. While lightning damage in mesoscale models is driven by the current distribution^4,68 and the boundary conditions⁴⁹ of the experiment and corresponding model, deriving the exact loading and grounding of an RVE scale model is challenging. Since these are location- and time-specific; these are basically unknown and require some idealisation. Therefore, this work has assumed that (current and mechanical) loading on an RVE surface will be a fraction of the total arc loading. Grounding of the side surfaces would force the current to flow into particular sides of the RVE modelling region. In reality, at a macroscale, the current distribution is driven by impedance and electrical potential. Therefore, this work has assumed that since the RVE loading surface will be totally encapsulated by even the smallest lightning arc channel (R ≥5 mm) the impedance and electrical potential will be uniform in the in-plane directions but will vary in the through-the-thickness. This is assumed to promote current flow through-the-thickness of the RVE, hence the choice of bottom surface grounding. Note that we considered three representative grounding conditions that may provide their upper and lower limits. Future work could further enhance and validate these assumptions by incorporating the RVE in a mesoscale model and completing multi-scale modelling to further understand the loading and grounding required for modelling of microscopic behaviour only i.e., creating a model as shown in Figure 6(b).

The material properties assumed in this work have been derived from the mesoscale models used in literature. For example, the electrical conductivity of the epoxy matrix was derived from the value used for through-thickness conductivity in mesoscale models (i.e. “resin-rich” regions) in preceding research.^2,4,19 Further work is required to establish the temperature dependent electrical and thermal properties of the individual fibres and matrix at conditions suitable for a lightning strike. Percolation effects,^48,69 which can influence through–thickness properties, should be considered in future work. These interfacial effects could be incorporated in the current RVE models using cohesive elements with thermal and electrical gap conductances. Higher and lower interfacial conduction could be simulated by adjusting the cohesive layer thickness and gap conductance. In addition, once the matrix reached a char state it was assumed the conductivity reached 1 × 10⁵ 1/Ω.mm, in line with previous research.^2,4,19 This value was used to represent matrix sublimation without complex element deletion. However, the exact conductivity of the char material may be smaller than the conductivity of carbon fibres in the through-thickness direction. Therefore, further development to incorporate accurate element deletion and matrix charring would enhance the models.

Conclusions

A series of three-dimensional (3D) representative volume element (RVE) models were developed to predict lightning-induced thermal and mechanical damage of unidirectional carbon fibre-reinforced epoxy matrix composite laminates. The RVE models were created with random fibre distribution in epoxy matrix, all measuring 50 μm × 50 μm × 50 μm. The effects of various grounding conditions (one-, two-, and four-side grounding), fibre volume fractions (V_f = 55, 60%, and 65%), and peak current amplitudes (10, 20, and 40 kA) of applied impulse waveforms on microscopic damage in RVE models were characterized to better understand various damage mechanisms associated with lightning strikes. The proposed microscale lightning thermo-mechanical damage models are fundamentally different from the normal mesoscale or macroscale analyses used to predict the overall lightning damage in composite laminates. The primary goal of this work was to characterize complex lightning damage mechanisms at a fundamental level, not to assess the severity of overall lightning damage in composite laminates.

A two-step sequential modelling procedure was applied to decouple (i) thermal damage predicted using a coupled thermal-electrical model and (ii) thermo-mechanical damage obtained from subsequent dynamic pressure loading, thermal loading, and its combined effects. Predicted thermal damage was almost constant up to 20 kA peak current, before increasing significantly at 40 kA peak current. An increase in V_f steadily decreased the extent of severe thermal damage in the RVE models. Thermo-mechanical simulations showed that subsequent dynamic pressure loading was negligible and thermal loading associated with thermal strain/expansion had the largest contribution to thermo-mechanical damage during a lightning strike.

One limitation of the proposed model is that it does not account for carbon fibre damage. This may be acceptable since each fibre is surrounded by widespread matrix damage, thus total lightning damage is nearly equivalent to the matrix damage. Accurate modelling of carbon fibre damage should consider sequential chemical reactions including oxidation, char formation due to matrix decomposition, and ablation/sublimation near their critical temperature (∼3000°C). Although one can model these complex damage mechanisms, carbon fibre damage is highly localised, so cannot represent total lightning damage in carbon/epoxy composites.

Footnotes

Declaration of Conflicting Interests

The author(s) declare no conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Scott L J Millen

Juhyeong Lee

References

Chemartin

Lalande

Peyrou

, et al. Direct effects of lightning on aircraft structure : analysis of the thermal, electrical and mechanical constraints. J Aerosp Lab 2012; 1–15.

Millen

SLJ

Murphy

. Modelling and analysis of simulated lightning strike tests: a review. Compos Struct 2021; 274: 114347. DOI: 10.1016/j.compstruct.2021.114347

Sun

, et al. Experimental study of the damage behaviour of laminated CFRP composites subjected to impulse lightning current. Composites Part B: Engineering 2022; 239: 109949. DOI: 10.1016/j.compositesb.2022.109949

Lee

Lacy Jr

Pittman Jr

. Coupled thermal electrical and mechanical lightning damage predictions to carbon/epoxy composites during arc channel shape expansion. Compos Struct 2021; 255: 112912. DOI: 10.1016/j.compstruct.2020.112912

Harrell

Madsen

Thomsen

, et al. On the effect of dielectric breakdown in UD CFRPs subjected to lightning strike using an experimentally validated model. Appl Compos Mater 2022; 29: 1321–1348. DOI: 10.1007/s10443-022-10014-7

Zhang

Cheng

, et al. Lightning strike damage on the composite laminates with carbon nanotube films: protection effect and damage mechanism. Composites Part B: Engineering 2019; 168: 342–352. DOI: 10.1016/j.compositesb.2019.03.054

Millen

SLJ

Ashworth

Farrell

, et al. Understanding and representing heating and heating rate effects on composite material properties for lightning strike direct effect simulations. Composites Part B: Engineering 2022; 228: 109438. DOI: 10.1016/j.compositesb.2021.109438

Boushab

Gharghabi

Lee

, et al. Lightning arc channel effects on surface damage development on a PRSEUS composite panel : an experimental study. Composites Part B: Engineering 2021; 224: 109217. DOI: 10.1016/j.compositesb.2021.109217

Lee

Lacy

Pittman

. Lightning mechanical damage prediction in carbon/epoxy laminates using equivalent air blast overpressure. Composites Part B: Engineering 2021; 212: 108649. DOI: 10.1016/j.compositesb.2021.108649

10.

Sun

, et al. An experimental study of impulse-current-induced mechanical effects on laminated carbon fibre-reinforced polymer composites. Composites Part B: Engineering 2021; 225: 109245. DOI: 10.1016/j.compositesb.2021.109245

11.

Xia

Mei

Zhang

, et al. Fabrication of the silver modified carbon nanotube film/carbon fiber reinforced polymer composite for the lightning strike protection application. Composites Part B: Engineering 2020; 180: 107563. DOI: 10.1016/j.compositesb.2019.107563

12.

Lee

Gharghabi

Boushab

, et al. Artificial lightning strike tests on PRSEUS panels. Composites Part B: Engineering 2018; 154: 467–477. DOI: 10.1016/j.compositesb.2018.09.016

13.

Hirano

Katsumata

Iwahori

, et al. Artificial lightning testing on graphite/epoxy composite laminate. Composites Part A: Applied Science and Manufacturing 2010; 41: 1461–1470. DOI: 10.1016/j.compositesa.2010.06.008

14.

Feraboli

Miller

. Damage resistance and tolerance of carbon/epoxy composite coupons subjected to simulated lightning strike. Composites Part A: Applied Science and Manufacturing 2009; 40: 954–967. DOI: 10.1016/j.compositesa.2009.04.025

15.

Feraboli

Kawakami

. Damage of carbon/epoxy composite plates subjected to mechanical impact and simulated lightning. J Aircr 2010; 47: 999–1012. DOI: 10.2514/1.46486

16.

Shah

SZH

Lee

. Stochastic lightning damage prediction of carbon/epoxy composites with material uncertainties. Compos Struct 2022; 282: 115014. DOI: 10.1016/j.compstruct.2021.115014

17.

Wang

, et al. Thermal dynamic damage of aircraft composite material suffered from lightning channel attachment based on moving mesh method. Compos Sci Technol 2021; 214: 109003. DOI: 10.1016/j.compscitech.2021.109003

18.

Liu

Wang

. Modeling the lightning continuing current electric arc discharge and material thermal damage: effects of combinations of amplitude and duration. International Journal of Thermal Sciences 2021; 162: 106786. DOI: 10.1016/j.ijthermalsci.2020.106786

19.

Chen

Wang

, et al. The coupling mechanism and damage prediction of carbon fiber/epoxy composites exposed to lightning current. Compos Struct 2018; 203: 436–445. DOI: 10.1016/j.compstruct.2018.07.017

20.

Abdelal

Murphy

. A multiphysics simulation approach for efficient modeling of lightning strike tests on aircraft structures. IEEE Trans Plasma Sci 2017; 45: 725–735. DOI: 10.1109/TPS.2017.2673543

21.

Millen

SLJ

Murphy

. Modelling and understanding edge glow effects on material failure resulting from artificial lightning strike. Compos Struct 2021; 278: 114651. DOI: 10.1016/j.compstruct.2021.114651

22.

Foster

Abdelal

Murphy

. Quantifying the influence of lightning strike pressure loading on composite specimen damage. Appl Compos Mater 2018; 26: 115–137. DOI: 10.1007/s10443-018-9685-1

23.

Karch

Arteiro

Camanho

. Modelling mechanical lightning loads in carbon fibre-reinforced polymers. Int J Solids Struct 2019; 162: 217–243. DOI: 10.1016/j.ijsolstr.2018.12.013

24.

Lee

Shah

SZH

. Predicting stochastic lightning mechanical damage effects on carbon fiber reinforced polymer matrix composites. 37th Am Soc Compos Tech Conf 2022, pp. 1–16. https://dpi-proceedings.com/index.php/asc37/article/view/36374.

25.

Sun

Tian

, et al. Experimental and numerical analysis of damage mechanisms for carbon fiber-reinforced polymer composites subjected to lightning strikes. Eng Fail Anal 2020; 118: 104894. DOI: 10.1016/j.engfailanal.2020.104894

26.

Millen

SLJ

Murphy

Catalanotti

, et al. Coupled thermal-mechanical progressive damage model with strain and heating rate effects for lightning strike damage assessment. Appl Compos Mater 2019; 26: 1437–1459. DOI: 10.1007/s10443-019-09789-z

27.

Wan

Dong

Zhi

, et al. Analysis on electrical and thermal conduction of carbon fiber composites under lightning based on electrical-thermal-chemical coupling and arc heating models. Compos Struct 2019; 229: 111486. DOI: 10.1016/j.compstruct.2019.111486

28.

Dong

Wan

Ping

, et al. Coupled thermal-mechanical damage model of laminated carbon fiber/resin composite subjected to lightning strike. Compos Struct 2018; 206: 185–193. DOI: 10.1016/j.compstruct.2018.08.043

29.

Dong

Wan

Guo

, et al. Damage analysis of carbon fiber composites exposed to combined lightning current components D and C. Compos Sci Technol 2019; 179: 1–9. DOI: 10.1016/j.compscitech.2019.04.030

30.

Sun

Tian

, et al. Lightning strike-induced dynamic conduction characteristics and damage behavior of carbon fiber-reinforced polymer composites. Compos Struct 2021; 275: 114391. DOI: 10.1016/j.compstruct.2021.114391

31.

Ahmed

Zheng

Zhang

, et al. Impact response of carbon/kevlar hybrid 3D woven composite under high velocity impact: experimental and numerical study. Appl Compos Mater 2020; 27: 285–305. DOI: 10.1007/s10443-020-09809-3

32.

Muñoz

Seltzer

Sket

, et al. Influence of hybridisation on energy absorption of 3D woven composites under low-velocity impact loading. Modelling and experimental validation. Int J Impact Eng 2022; 165: 104229. DOI: 10.1016/j.ijimpeng.2022.104229

33.

Wielhorski

Mendoza

Rubino

, et al. Numerical modeling of 3D woven composite reinforcements: a review. Composites Part A: Applied Science and Manufacturing 2022; 154: 106729. DOI: 10.1016/j.compositesa.2021.106729

34.

Sevenois

RDB

Garoz

Gilabert

, et al. Avoiding interpenetrations and the importance of nesting in analytic geometry construction for Representative Unit Cells of woven composite laminates. Compos Sci Technol 2016; 136: 119–132. DOI: 10.1016/j.compscitech.2016.10.010

35.

Ullah

Zhou

Kaczmarczyk

, et al. A unified framework for the multi-scale computational homogenisation of 3D-textile composites. Composites Part B: Engineering 2019; 167: 582–598. DOI: 10.1016/j.compositesb.2019.03.027

36.

Gao

Shakoor

Domel

, et al. Predictive multiscale modeling for unidirectional carbon fiber reinforced polymers. Compos Sci Technol 2020; 186: 107922. DOI: 10.1016/j.compscitech.2019.107922

37.

Liao

Tan

Zhou

, et al. Multi-scale modelling of dynamic progressive failure in composite laminates subjected to low velocity impact. Thin-Walled Structures 2018; 131: 695–707. DOI: 10.1016/j.tws.2018.07.047

38.

Pathan

Tagarielli

Patsias

, et al. A new algorithm to generate representative volume elements of composites with cylindrical or spherical fillers. Composites Part B: Engineering 2017; 110: 267–278. DOI: 10.1016/j.compositesb.2016.10.078

39.

Park

S-M

Lim

Seong

, et al. Efficient generator of random fiber distribution with diverse volume fractions by random fiber removal. Composites Part B: Engineering 2019; 167: 302–316. DOI: 10.1016/j.compositesb.2018.12.042

40.

Cai

, et al. A new multi-scale framework for analyzing the thermo-mechanical coupling mechanism of the composite laminates. Mech Adv Mat Struct 2022: 1–15. DOI: 10.1080/15376494.2022.2122641

41.

Yang

Cao

, et al. Micromechanical modelling and simulation of unidirectional fibre-reinforced composite under shear loading. J Reinf Plast Compos 2015; 34: 72–83. DOI: 10.1177/0731684414562873

42.

Totry

Gonzalez

Llorca

. Influence of the loading path on the strength of fiber-reinforced composites subjected to transverse compression and shear. Int J Solids Struct 2008; 45: 1663–1675. DOI: 10.1016/j.ijsolstr.2007.10.014

43.

Sun

Zhou

Meng

, et al. Failure criteria of unidirectional carbon fiber reinforced polymer composites informed by a computational micromechanics model. Compos Sci Technol 2019; 172: 81–95. DOI: 10.1016/j.compscitech.2019.01.012

44.

Chen

Wan

Ismail

, et al. Micromechanical analysis of UD CFRP composite lamina under multiaxial loading with different loading paths. Compos Struct 2021; 269: 114024. DOI: 10.1016/j.compstruct.2021.114024

45.

González

Llorca

. Mechanical behavior of unidirectional fiber-reinforced polymers under transverse compression: microscopic mechanisms and modeling. Compos Sci Technol 2007; 67: 2795–2806. DOI: 10.1016/j.compscitech.2007.02.001

46.

Guo

Wang

, et al. Enhanced lightning strike protection of carbon fiber composites using expanded foils with anisotropic electrical conductivity. Composites Part A: Applied Science and Manufacturing 2019; 117: 211–218. DOI: 10.1016/j.compositesa.2018.11.022

47.

Wang

Huang

, et al. Thermal-damage suppression of composite material by anti-lightning NCF/CIPCF film. Compos Sci Technol 2022; 227: 109604. DOI: 10.1016/j.compscitech.2022.109604

48.

Chippendale

Golosnoy

. Percolation effects in electrical conductivity of carbon fibre composites. In: IET 8th Int. Conf. Comput. Electromagn. (CEM 2011), 2011, p. 1–2. DOI: 10.1049/cp.2011.0094

49.

Millen

SLJ

Murphy

. Understanding the influence of test specimen boundary conditions on material failure resulting from artificial lightning strike. Eng Fail Anal 2020; 114: 104577. DOI: 10.1016/j.engfailanal.2020.104577

50.

Millen

SLJ

Murphy

Abdelal

, et al. Sequential finite element modelling of lightning arc plasma and composite specimen thermal-electric damage. Comput Struct 2019; 222: 48–62. DOI: 10.1016/j.compstruc.2019.06.005

51.

SAE Aerospace . Aerospace Recommended Practice ARP5412B; Vol. 4970; 1999.

52.

Ogasawara

Hirano

Yoshimura

. Coupled thermal-electrical analysis for carbon fiber/epoxy composites exposed to simulated lightning current. Composites Part A: Applied Science and Manufacturing 2010; 41: 973–981, https://doi.org/10.1016/j.compositesa.2010.04.001

53.

Xiong

Dai

, et al. The development of shock wave overpressure driven by channel expansion of high current impulse discharge arc. Phys Plasmas 2018;25:032115-1–7. DOI: 10.1063/1.5013296

54.

Lee

Lacy

Pittman

, et al. Numerical estimations of lightning-induced mechanical damage in carbon/epoxy composites using shock wave overpressure and equivalent air blast overpressure. Compos Struct 2019; 224: 111039. DOI: 10.1016/j.compstruct.2019.111039

55.

ABAQUS 2016 documentation. ABAQUS Theory Manual. 2017.

56.

Lani

Morelle

Bailly

, et al. Characterization and modelling of the strain-rate, temperature and pressure dependence of the deformation of a highly crosslinked aerospace grade epoxy resin, 2015.

57.

Littell

Ruggeri

Goldberg

, et al. Measurement of epoxy resin tension, compression, and shear stress–strain curves over a wide range of strain rates using small test specimens. J Aerosp Eng 2008; 21. DOI: 10.1061/(ASCE)0893-1321(2008)21:3(162)

58.

Hooputra

Gese

Dell

, et al. A comprehensive failure model for crashworthiness simulation of aluminium extrusions. International Journal of Crashworthiness 2004; 9: 449–464. DOI: 10.1533/ijcr.2004.0289

59.

Khaled

Shyamsunder

Robbins

, et al. Framework for predicting failure in polymeric unidirectional composites through combined experimental and computational mesoscale modeling techniques. Fibers 2021; 9: 50. DOI: 10.3390/fib9080050

60.

Wan

Ismail

Zhu

, et al. Computational micromechanics-based prediction of the failure of unidirectional composite lamina subjected to transverse and in-plane shear stress states. J Compos Mater 2020; 54: 3637–3654. DOI: 10.1177/0021998320918015

61.

Benzeggagh

Kenane

. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-. Modebending Apparatus 1996; 56: 1–11.

62.

Kamiyama

Hirano

Ogasawara

. Delamination analysis of CFRP laminates exposed to lightning strike considering cooling process. Compos Struct 2018; 196: 55–62. DOI: 10.1016/j.compstruct.2018.05.003

63.

Naya

González

Lopes

, et al. Computational micromechanics of the transverse and shear behavior of unidirectional fiber reinforced polymers including environmental effects. Composites Part A: Applied Science and Manufacturing 2017; 92: 146–157. DOI: 10.1016/j.compositesa.2016.06.018

64.

Melro

Camanho

Andrade Pires

, et al. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part I – constitutive modelling. Int J Solids Struct 2013; 50: 1897–1905. DOI: 10.1016/j.ijsolstr.2013.02.009

65.

Varna

Berglund

Ericson

. Transverse single-fibre test for interfacial debonding in composites: 2. Modelling. Composites Part A: Applied Science and Manufacturing 1997; 28: 317–326. DOI: 10.1016/S1359-835X(96)00125-X

66.

Lee

Lacy

Pittman

, et al. Thermal response of carbon fiber epoxy laminates with metallic and nonmetallic protection layers to simulated lightning currents. Polym Compos 2017; 39: E2149–E2166. DOI: 10.1002/pc.24502

67.

, et al. Experimental study of damage characteristics of carbon woven fabric/epoxy laminates subjected to lightning strike. Composites Part A: Applied Science and Manufacturing 2015; 79: 164–175. DOI: 10.1016/j.compositesa.2015.09.019

68.

Foster

Abdelal

Murphy

. Understanding how arc attachment behaviour influences the prediction of composite specimen thermal loading during an artificial lightning strike test. Compos Struct 2018; 192: 671–683. DOI: 10.1016/j.compstruct.2018.03.039

69.

Choi

H-J

Kim

Ahn

, et al. Electrical percolation threshold of carbon black in a polymer matrix and its application to antistatic fibre. Sci Rep 2019; 9: 6338. DOI: 10.1038/s41598-019-42495-1

Microscale modelling of lightning damage in fibre-reinforced composites

Abstract

Keywords

Introduction

Representative volume element modelling

Background summary and paper aim

Methodology

Representative volume element modelling

Thermal-electric representative volume element model

Thermo-mechanical model

Model summary

Thermal damage results

Effects of grounding condition

Effects of peak current amplitude

Effects of fibre volume fraction

Thermo-mechanical damage results

Pressure loading only

Temperature loading only

Combined pressure-temperature loading

Results summary

Limitations and future work

Conclusions

Footnotes

Declaration of Conflicting Interests

Funding

ORCID iDs

References