Nomenclature

3D

three-dimensional

$a$ , $c$

waveform and root radius fitting parameters

$b$ , $f$ , $\nu$

waveform fitting parameters

B-K

Benzeggagh-Kenane

CFRP

carbon fibre-reinforced polymer

$d$

interlaminar damage variable

${d_i}$

intralaminar damage variables

${E_0}$

initial channel energy density

${E_{i{\rm{\;}}}}$

Young’s moduli

ECF

expanded copper foil

$f( \eta )$

dimensionless supersonic shock wave pressure function

${\,f_{fxc{\rm{\;}}}}$

longitudinal compressive strength ratio scaling parameter

${\,f_{\mathcal{G}1 + }}$

fraction of the fracture toughness dissipated in the second arm of the bilinear softening law

${\,f_{XC}}$

longitudinal compressive strength ratio at the inflection point

$f_{XC}^{ef}$

effective longitudinal compressive strength ratio at the inflection point

${\,f_{XT}}$

longitudinal tensile strength ratio at the inflection point

${\mathcal{G}_{1 + }}$

longitudinal tensile fracture toughness

${\mathcal{G}_{1 - }}$

longitudinal compressive fracture toughness

${\mathcal{G}_{2 + }}$

transverse tensile fracture toughness

${\mathcal{G}_6}$

in-plane shear fracture toughness

$\mathcal{G}_6^{ef}$

effective in-plane shear fracture toughness

${\mathcal{G}_{Ic{\rm{\;}}}}$

mode I interlaminar fracture toughness

${\mathcal{G}_{IIc}}$

mode II interlaminar fracture toughness

${G_{ij}}$

shear moduli

${I_3}$

stress invariant of transverse isotropy

${I_{{\rm{max}}}}$

peak current

$J(t)$

transient lightning current amplitude

$k$

spatial distribution parameter of the transient plasma pressure

$K$

cohesive surface stiffness

${K_P}$

shear modulus degradation factor

LSP

lightning strike protection

$p\left( {\eta ,t} \right)$

transient supersonic shock wave pressure

${p_{\rm{e}}}( {r,t} )$

spatial distribution of the transient plasma pressure

${p_{\rm{e}}}(t)$

transient effective plasma pressure

${p_{\rm{m}}}( {r,t} )$

effective magnetic surface pressure

$r$

radial coordinate from the lightning attachment point

$R(t)$

transient shock wave radius

${R_R}(t){\rm{\;}}$

transient arc root radius

${R_{{R_0}}}$

initial arc root radius

${S_L}$

in-plane shear strength

${S_{LP}}$

in-plane shear yield stress

$t$

time variable

${t_{{\rm{max}}}}$

time at the peak current

${t_s}$

root radius fitting parameter

UD

unidirectional

VISAR

velocity interferometer system for any reflector

X, Y, Z

coordinate axes

${x_0}$ , ${y_0}$

lightning attachment (offset) coordinates

${X_C}$

longitudinal compressive strength

${X_T}$

longitudinal tensile strength

$Y_C$

transverse compressive strength

${Y_T}$

transverse tensile strength

Greek symbol

$\alpha $ , $\beta $

root radius fitting parameters

${\alpha _0}$

fracture angle under uniaxial transverse compression

$\beta $

supersonic shock wave pressure parameter for $t \lt 0.5\ \mu \textrm{s}$

$\eta $

dimensionless radius

${\eta _{{\rm{B}} - {\rm{K}}}}$

interlaminar fracture toughness mixed-mode parameter

${\eta _G}$

in-plane shear fracture toughness scaling parameter

${\mu _0}$

permeability of the free space

$\mu_r$

relative permeability

$\nu_{ij}$

Poisson’s ratios

${\rho _0}$

initial mass (air) density

${\sigma _{ij}}$

components of the stress tensor

${\tau _1}$ , ${\tau _2}$

interlaminar shear strengths

${\tau _3}$

interlaminar normal strength

2.0 Mechanical lightning loads

In protected CFRP subjected to a transient lightning current, the shock waves generated by the exploding vaporised materials at the arc attachment area (or within the material volume extremely heated by the generated Joule heat) confined between the composite material and the paint layer are the main contribution to lightning-induced mechanical damage [Reference Karch, Honke, Steinwandel and Dittrich4, Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5, Reference Duval, Lepetit, Revel and Peres27, Reference Wulbrand and Karch28] (see the generated mechanical impulses in Figs 4 and 5). It should be noted that due to the extremely high energy density of a lightning strike at the arc root area, an extremely short evaporation of the material cannot be avoided, which then leads to an explosion and ultimately the generation of shock waves. However, the contributions from the shock waves caused by the supersonic expansion of the hot plasma channel [Reference Karch, Honke, Steinwandel and Dittrich4, Reference Plooster29–Reference Karch, Wulbrand and Müller31] and the magnetic volume forces (magnetic pressures) caused by the impressed current flow in the electrically conducting structures [Reference Karch, Honke, Steinwandel and Dittrich4, Reference Karch and Wulbrand32] cannot be neglected (see Section 5), especially in the unpainted cases. Because the presence of the LSP layer minimises thermal damage and temperature rise in the underlying plies [Reference Millen, Murphy, Abdelal and Catalanotti13], in a first approximation, thermal and thermal-mechanical effects can be neglected [Reference Pedro, Gonçalves, Soares, Arteiro, Honke and Karch33].

2.1 Root radius

The root radius expansion ${R_R}(t)$ is calculated using the finite volume degradation approach [Reference Karch, Arteiro and Camanho25]. The root radius data computed by Karch et al. [Reference Karch, Arteiro and Camanho25] for two protected CFRP with different expanded copper foils (ECFs) – 2Cu4-100FA, with a surface weight of 73.3 g/m ${{\rm{\;}}^2}$ (ECF 73), and 3Cu7-100FA, with a surface weight of 195.3 g/m ${{\rm{\;}}^2}$ (ECF 195) – are fitted using a double exponential function:

(1) \begin{equation}{R_R}(t) = \begin{cases}{{R_{{R_0}}} + a\big( {{e^{\alpha t}} - {e^{\beta t}}} \big),} & {}t \le {t_s}\\{c,} & {}{t \gt {t_s}}\end{cases} \end{equation}

with ${R_{{R_0}}} = $ 0.5 mm ensuring an initial arc root radius of 0.5 mm independently of the parameters $a$ , $\alpha $ and $\beta $ , which are obtained by curve fitting of the finite volume results [Reference Karch, Arteiro and Camanho25].

2.2 Near-surface explosion of the lightning protection layer

The explosion of the LSP layer is assumed to be caused mainly by the rapid vaporisation of the metallic mesh of the LSP layer due to Joule resistive heating and direct heat transfer from the plasma channel at the arc root area. It is assumed that the spatial distribution of the expanding plasma pressure can be approximated by a Gaussian bell-shaped curve:

(2) \begin{equation}{p_e}( {r,t} ) = {p_e}(t){\rm{\;\;}}{k^{ - 2}}{\rm{\;\;exp}}\left[ { - \frac{{{r^2}}}{{{{( {k{\rm{\;\;}}{R_R}(t)} )}^2}}}} \right]\end{equation}

where the effective plasma overpressure ${p_e}(t)$ is numerically determined as described by Karch et al. [Reference Karch, Arteiro and Camanho25]. The parameter $k \lt 1$ in Equation (2) is used to limit the pressure distribution mainly to the arc root area $\pi R_R^2(t)$ by controlling the shape of the Gaussian curve.

2.3 Supersonic plasma expansion

A supersonic shock wave is created by the overpressure caused by the rapid temperature rise of the plasma channel. Following Lin’s approach [Reference Lin34], the shock wave expands radially according to the following equation [Reference Karch, Arteiro and Camanho25]:

(3) \begin{equation}R(t) \cong {\begin{cases}{1.004{{\left( {\dfrac{{{E_0}}}{{{\rho _0}}}} \right)}^{0.25}}\sqrt {0.5}}\, {}, & t \lt 0.5\ \mu s\\[8pt]{1.004{{\left( {\dfrac{{{E_0}}}{{{\rho _0}}}} \right)}^{0.25}}\sqrt t }\, {}, & t \ge 0.5\ \mu s\end{cases}} ,\end{equation}

where ${E_0}$ is the initial channel energy density and ${\rho _0}$ the initial mass (air) density. The following relationship between the initial channel energy density ${E_0}$ (J/m) and the lightning current peak ${I_{{\rm{max}}}}$ (A) can be established [Reference Karch, Wulbrand and Müller31, Reference Karch, Schreiner, Honke and Wolfrum35]:

(4) \begin{equation}{E_0} \cong 0.45 \times {10^{ - 2}}{\rm{\;}}{( {{I_{{\rm{max}}}}} )^{1.25}}\end{equation}

Finally, the shock wave pressure is described by:

(5) \begin{equation}p( {\eta ,t} ) \cong {\begin{cases}{\beta \cdot t} {}, & t \lt 0.5\, \mu s\\[3pt]{0.180{{\left( {{\rho _0}{E_0}} \right)}^{0.5}}f( \eta ){\rm{\;}}{t^{ - 1}} - 1} {}, & t \ge 0.5\, \mu s\\[3pt]0 {}, & t \gt 100\, \mu s\end{cases}}\!,\end{equation}

where $\beta $ ensures a linear evolution of the pressure with time before the radial shock wave detaches from the hot, overpressured lightning channel, $f( \eta )$ is a dimensionless pressure function and $\eta = r/R(t)$ is the dimensionless radius [Reference Karch, Wulbrand and Müller31].

2.4 Magnetic forces

The current flow in the LSP layer creates a magnetic pressure whose amplitude depends on whether it lies inside ( $0 \lt r \lt {R_R}(t)$ ) or outside the arc root radius ( $r \gt {R_R}(t)$ ) [Reference Karch, Arteiro and Camanho25]. Taking this into account, the effective magnetic surface pressure acting on the upper side of the LSP layer is given as:

(6) \begin{equation}{p_{\rm{m}}}( {r,t} ) = {\begin{cases}{\dfrac{{{\mu _r}{\rm{\;}}{\mu _0}{\rm{\;}}J{{(t)}^2}}}{{8{\pi ^2}}}\dfrac{{{r^2}}}{{R_R^4(t)}}} {}, {}&0 \lt r \lt {R_R}(t)\\[9pt]{\dfrac{{{\mu _r}{\rm{\;}}{\mu _0}{\rm{\;}}J{{(t)}^2}}}{{8{\pi ^2}{r^2}}}} {}, {}&r \gt {R_R}(t)\end{cases}} \end{equation}

where ${\mu _0} = 4\pi \times {10^{ - 7}}$ Hm ${{\rm{\;}}^{ - 1}}$ is the permeability of the free space, ${\mu _r} = 1$ is the relative permeability of the non-magnetic LSP layer, and $J(t)$ is the total lightning current amplitude [Reference Karch, Arteiro and Camanho25].

3.0 Composite damage model

A modification [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26] to the continuum damage mechanics model proposed by Maimí et al. [Reference Maimí, Camanho, Mayugo and Dávila36, Reference Maimí, Camanho, Mayugo and Dávila37] is implemented in a user-defined subroutine VUMAT for the commercial explicit finite element solver Abaqus/Explicit [38] to model failure of the composite plies. Each ply, which is assumed transversely isotropic, must be represented explicitly in the finite element model. To account for ply thickness effects, in-situ strengths are defined as a function of ply thickness [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26, Reference Camanho, Dávila, Pinho, Iannucci and Robinson39]. Therefore, adjacent plies with the same fibre orientation are modelled by a single element through the thickness [Reference Arteiro, Gray and Camanho40].

Bilinear softening laws are employed to model longitudinal damage growth in tension and in compression (according to the parameters ${\,f_{\mathcal{G}1 + }}$ , ${\,f_{XT}}$ and ${\,f_{XC}}$ ). To represent the effects of hydrostatic pressure [Reference Camanho, Arteiro, Melro, Catalanotti and Vogler41], the in-plane shear strength ${S_L}$ , the fracture toughness for longitudinal compression ( ${\mathcal{G}_{1 - }}$ ) and the longitudinal compressive strength ratio at the inflection point ( ${\,f_{XC}}$ ) are scaled as a function of the applied hydrostatic pressure [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26]. In addition, the fracture toughness for in-plane shear ( ${\mathcal{G}_6}$ ) is scaled to account for the effect of transverse compressive stresses on the fracture toughness associated with in-plane shear fracture [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26].

To ensure mesh-independent results after damage onset, an energy regularisation approach based on the fracture energy associated to each failure mechanism and based on the characteristic length of the finite elements is used [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26, Reference Bažant and Oh42].

Delamination between plies is modelled using the interaction properties with cohesive behaviour available in Abaqus [38]. The cohesive interaction is defined by a bilinear traction-separation damage law accounting for mode dependency according to the B-K law [Reference Benzeggagh and Kenane43], and the onset of interlaminar damage is predicted by a quadratic stress-based criterion. A friction coefficient of 0.3 is assumed [Reference Arteiro, Gray and Camanho40].

The continuum damage mechanics and cohesive zone models [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26, Reference Arteiro, Gray and Camanho40] employed in this work do not include strain rate effects. Although the authors recognise that strain rate effects would be required for higher fidelity of the mechanical response predictions, and despite the availability of strain rate-dependent formulations in literature [Reference Rodrigues Lopes, Camanho, Andrade Pires and Arteiro44, Reference Cózar, Maimí, González, Camanho and Otero45], in the present work, the authors opted for employing a composite damage model that has been thoroughly validated in the past [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26, Reference Pedro, Gonçalves, Soares, Arteiro, Honke and Karch33, Reference Arteiro, Gray and Camanho40, Reference Ni, Kopp, Kalfon-Cohen, Furtado, Lee, Arteiro, Borstnar, Mavrogordato, Helfen, Sinclair, Spearing, Camanho and Wardle46–Reference Vallmajó, Descamps, Arteiro and Turon48] before moving to newer formulations.

4.0 Composite material, coupon geometry and finite element model

Following Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5], CFRP samples made from eight T700/M21 unidirectional (UD) plies with a quasi-isotropic [45/0/–45/90] $_{\rm{S}}$ lay-up are considered in this study. Each T700/M21 UD ply is 0.262 mm thick. Two sets of samples were tested by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5], protected using ECF 73 and ECF 195, respectively. The samples were subjected to standard arc attachment tests, with a lightning waveform of D type. The experimental setup included a very fast camera (1,000,000 fps) filtered to avoid saturation of the sensors by the light emitted by the electrical arc in the front face of the sample, where the lightning discharge is applied. In the back face, Doppler interferometry (or velocity interferometer system for any reflector – VISAR) was applied to measure deflection of the sample. Five VISAR sensors were used to measure the displacement speed on different points of the sample, located on a square with 49 mm-long sides (illustrated in the finite element mesh in Fig. 2). Sensor 5 (IDF5) was positioned in line with the position of the electrode, on the opposite face. A polyurethane paint system with a thickness of 200 ${\rm{\mu }}$ m was applied to each sample over a surfacing film of 127 ${\rm{\mu }}$ m. Assuming that part of the film fills the surface roughness, a 300 ${\rm{\mu }}$ m dielectric coating (paint) layer is assumed in this study.

Figure 2. Composite plate finite element mesh.

Tables 1 and 2 show, respectively, the UD ply elastic properties and the strength and material degradation properties used in the simulations, including the in-plane yield stress and shear modulus degradation factor ( ${S_{LP}}$ and ${K_P}$ ), the longitudinal tensile and compressive plain strengths ( ${X_T}$ and ${X_C}$ ), the UD transverse tensile and compressive strengths ( ${Y_T}$ and ${Y_C}$ ), the assumed composite fracture angle in uniaxial transverse compression ( ${\alpha _0}$ ), the UD in-plane shear strength ( ${S_L}$ ), the longitudinal tensile fracture toughness and the fraction of the fracture toughness dissipated in the second arm of the bilinear softening law ( ${\mathcal{G}_{1 + }}$ and ${\,f_{\mathcal{G}1 + }}$ ), the longitudinal compressive, transverse tensile and in-plane shear fracture toughness ( ${\mathcal{G}_{1 - }}$ , ${\mathcal{G}_{2 + }}$ and ${\mathcal{G}_6}$ ), the frictional parameter ( ${\,f_{fxc}}$ ) used to scale the crushing (compression) strength factor ( ${\,f_{XC}}$ ) as a function of the hydrostatic pressure [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26]:

(7) \begin{equation}f_{XC}^{ef} = {\rm{min}}\{ {1,{\,f_{XC}}( {1 - {\,f_{fxc}} \langle -{I_3}\rangle} )} \},{\rm{\;\;\;\;}}{I_3} = {\sigma _{22}} + {\sigma _{33}}\end{equation}

and the material dependent empirically derived enhancement factor ( ${\eta _G}$ ) used in the definition of the effective fracture toughness for in-plane shear in the presence of transverse compressive stresses [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26]:

(8) \begin{equation}\mathcal{G}_6^{ef} = {\mathcal{G}_6}( {1 - {\eta _G} \langle- {\sigma _{22}}}\rangle )\end{equation}

Table 1. T700/M21 UD ply density, Young’s moduli, shear moduli and Poisson’s ratios

Table 2. T700/M21 strengths and material degradation properties [Reference Huchette50–Reference Charrier53] used in the continuum damage mechanics model [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho26]

Table 3 shows the interlaminar properties used in the simulations, including the stiffness of the cohesive surfaces before damage onset ( $K$ ) and the mixed-mode parameter of the B-K law ( ${\eta _{B - K}}$ ) [Reference Benzeggagh and Kenane43]. The shear strengths ( ${\tau _1} = {\tau _2}$ ) are determined as a function of the mode I and mode II fracture toughness ( ${\mathcal{G}_{Ic}}$ and ${\mathcal{G}_{IIc}}$ ), and the normal strength ( ${\tau _3}$ ) [Reference Turon, Camanho, Costa and Renart49]:

(9) \begin{equation}{\tau _1} = {\tau _2} = {\tau _3}\sqrt {\frac{{{\mathcal{G}_{IIc}}}}{{{\mathcal{G}_{Ic}}}}} \end{equation}

Table 3. T700/M21 interlaminar properties

The effective properties of the LSP layers ECF 73 and ECF 195 filled with epoxy resin M21 were determined using a micro-mechanical finite element homogenisation approach by Karch et al. [Reference Karch, Arteiro and Camanho25]. The properties used in the simulations are given in Table 4.

Table 4. Effective properties of ECF 73 and ECF 195 filled with epoxy resin M21

The tested samples are square laminated plates (450 mm-long sides) supported on a fixed plate with a circular opening (340 mm diameter), fixed using screws placed along a circumference concentric with the central opening of the support plate (370 mm diameter). The CFRP sample is modelled fixing the out-of-plane displacements on the supported area outside the circular opening and fixing the nodes along the circumference where the screws are placed (Fig. 2). The finite element mesh used in this study was carefully chosen to ensure that the local response below the impact point is accurately captured, especially in the initial nanoseconds corresponding to the peak value of the overpressure caused by the explosion of the LSP layer (see Fig. 3) within a very small region (expanding arc root) [Reference Pedro, Gonçalves, Soares, Arteiro, Honke and Karch33]. The 8-node C3D8R 3D brick elements with reduced integration are used to model each layer of the composite structure under consideration, including the LSP layer. It is assumed that, due to its soft behaviour, the stiffness and strength of the paint is negligible, and therefore it is not explicitly represented in the finite element model, but implicitly accounted for in the computed overpressure caused by the explosion of the surface conducting materials.

Figure 3. Calculated overpressure caused by the explosion of the LSP layers as a function of time.

5.0 Lightning current and lightning loads implementation

To represent the measured lightning transient D current waveform tested by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5], a modified double exponential function was proposed by Karch et al. [Reference Karch, Arteiro and Camanho25]:

(10) \begin{equation}J(t) = {\begin{cases}{c( {{e^{ - at}} - {e^{ - bt}}} )} {}, & {}{t \lt {t_{{\rm{max}}}}}\\[3pt]{c( {{e^{ - af( {t - \nu} )}} - {e^{ - bf( {t - \nu} )}}} )} {},& {}{t \ge {t_{{\rm{max}}}}}\end{cases}} \end{equation}

where ${t_{{\rm{max}}}}$ is the time at the current peak. The parameters for the D current waveform with a time to peak of 26.8 ${\rm{\mu }}$ s, a decay time of 48.5 ${\rm{\mu }}$ s and a peak current of 96.4 kA [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5] are summarised in Table 5.

Table 5. Waveform parameters for the 26.8/48.5 ${\rm{\mu }}$ s, 96.4 kA modified double exponential function

This modified D current waveform [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5] is then used to calculate the root radii, following Karch et al. [Reference Karch, Arteiro and Camanho25]. The parameters of the fitting function in Equation (1) are obtained using the curve fitting application from MATLAB. These parameters are shown in Table 6.

Table 6. Parameters for root radius fitting of the correct transient D current waveform with a time to peak of 26.8 ${\rm{\mu }}$ s, a decay time of 48.5 ${\rm{\mu }}$ s and a peak current of 96.4 kA [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5]

The effective overpressure caused by the explosion of the LSP layer is directly dependent on the root radius. Therefore, the pressure magnitudes for the two ECFs are determined considering the same transient D current waveform [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5]. The effective plasma overpressure is generated considering protected CFRP without paint and with a paint thickness (dielectric coating) of 200, 300 and 400 ${\rm{\mu }}$ m. Figure 3 shows the effective plasma overpressure determined following the approach proposed by Karch et al. [Reference Karch, Arteiro and Camanho25], with the parameter $k$ in Equation (2) set to 0.4 and 0.1 for ECF 73 and ECF 195, respectively. This parameter was selected based on the available experimental data for a paint thickness of 300 ${\rm{\mu }}$ m [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5, Reference Karch, Arteiro and Camanho25] and a sensitivity analysis using the composite damage model and the finite element model presented respectively in Sections 3 and 4. It is important to emphasise that parameter $k$ is fitted for each LSP layer considering the results for just one case (with a defined paint thickness), since $k$ is assumed to be mainly dependent on the properties of the LSP layer, and not the paint thickness; no other parameter in Section 2 is calibrated to fit velocity or deflection test data. Figure 4 shows the impulses generated by the explosion of the LSP layers as a function of time considering the protected CFRPs with ECF 73 and ECF 195, without paint and with a paint thickness of 200, 300 and 400 ${\rm{\mu }}$ m. As can be observed in Figs 3 and 4, independently of the LSP layer, the overpressure and, consequently, the impulse due to the explosion of the LSP layers are mostly negligible if the structure is unpainted. The presence of paint leads to significantly higher loading, whose magnitude increases with increasing paint thickness.

Figure 4. Calculated impulse caused by the explosion of the LSP layers as a function of time.

The supersonic plasma expansion overpressure is determined using the peak current ${I_{{\rm{max}}}}$ of the transient D current waveform [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5] in Equation (4). In this case, ${\rm{\beta }} = 392.2952$ bar/ ${\rm{\mu }}$ s. Because the supersonic plasma expansion overpressure is just a function of the current amplitude, its contribution to mechanical loads is the same independently of the LSP layer used.

To determine the pressure due to the magnetic forces given by Equation (6), the total lightning current amplitude $J(t)$ approximated by the modified double exponential function in Equation (10) is used [Reference Karch, Arteiro and Camanho25]. Figure 5 shows the impulses generated by the supersonic plasma expansion and by the magnetic forces. Compared with the impulses caused by the explosion of the LSP layers in Fig. 4, it is clear that all contributions are comparable in the case of unpainted structures. In the case of painted structures, and for the paint thicknesses considered in this study, the effects of the supersonic plasma expansion and of the magnetic forces (Fig. 5) are lower than the effects of the explosion of the LSP layer (Fig. 4), which has the highest contribution to the mechanical lightning loads. Nevertheless, and even in the case of painted structures, the supersonic plasma expansion and the magnetic forces cannot be neglected.

Figure 5. Calculated impulse caused by the supersonic plasma expansion and by the magnetic forces as a function of time.

Following Karch et al. [Reference Karch, Arteiro and Camanho25], the mechanical pressure fields resulting from the near-surface explosion of the LSP layer, from the supersonic plasma expansion and from the magnetic field caused by the impressed current flow in the electrically conducting structures are implemented in the finite element software Abaqus [38] using the user-defined subroutine VDLOAD and applied to the top surface of the LSP layer in the composite plate. Because the VISAR deflection measurements conducted by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5] indicate that the lightning strike did not initiate from the middle of the CFRP samples, the pressure fields are centred at the average offset coordinates (Table 7), obtained by fitting a 2D symmetric Gauss function to the deflection measurements for different times [Reference Karch, Arteiro and Camanho25].

Table 7. Offset coordinates for ECF 73 and ECF 195

6.0 Finite element results

Using the finite element model described in Section 4, the mechanical response and damage of the composite panels tested by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5] are predicted, including two LSP layers and an effective paint thickness (including surfacing film) of 300 ${\rm{\mu }}$ m. The predictions are then extended to the unpainted case and to a paint thickness of 200 and 400 ${\rm{\mu }}$ m for the two LSP layers, ECF 73 and ECF 195. It is important to note that the overpressure due to the near-surface explosion of the LSP layer is the lightning load that changes with the thickness of the paint layer. Therefore, this is the dominant lightning load contribution for the effect of paint thickness (see Fig. 3).

6.1 ECF 73 LSP layer

Figures 6–9 show the predicted displacement and velocity histories on the five points identified in Fig. 2 for the unpainted and painted cases with different paint thicknesses. For 300 ${\rm{\mu }}$ m paint thickness, the predictions are compared with the VISAR experimental measurements performed by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5]; the displacement and velocity histories correlate well between simulations and experiments, validating the hypothesis that the paint layer does not need to be explicitly modelled in the mechanical finite element analysis. It is also clear that the presence of paint, and consequent enhancement of the contribution of the near-surface explosion of the LSP layer, leads to a significantly higher maximum displacement and velocity for the painted cases when compared with the unpainted case. Moreover, as the paint thickness increases, the maximum displacement and velocity also increases, demonstrating that implicitly accounting for the paint layer in the formulation of the mechanical lightning loads allows capturing the trends observed in experiments [Reference Bigand and Duval2, Reference Murillo, Flourens, Garcia, Duval and Cavaliere7].

Figure 6. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for the unpainted case of ECF 73.

Figure 7. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for a paint thickness of 200 ${\rm{\mu }}$ m and ECF 73.

Figure 8. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for a paint thickness of 300 ${\rm{\mu }}$ m and ECF 73, and comparison with VISAR experimental measurements [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5].

Figure 9. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for a paint thickness of 400 ${\rm{\mu }}$ m and ECF 73.

Figures 10–13 show, respectively, the predictions of longitudinal, transverse, shear and interlaminar damage on ECF 73 for the unpainted case and for different paint thicknesses. No damage is observed in the unpainted panels, while an increasing extent of intralaminar and interlaminar damage with increasing paint thickness is observed in the panels coated with a paint layer.

Figure 10. Longitudinal damage ( ${d_1}$ ) extent considering the effect of paint thickness for ECF 73. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

Figure 11. Transverse damage ( ${d_2}$ ) extent considering the effect of paint thickness for ECF 73. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

Figure 12. Shear damage ( ${d_6}$ ) extent considering the effect of paint thickness for ECF 73. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

Figure 13. Interlaminar damage ( $d$ ) extent considering the effect of paint thickness for ECF 73. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

6.2 ECF 195 LSP layer

Following the study of ECF 73, Figs 14–17 now show the numerical results obtained for ECF 195. Again, the predictions for 300 ${\rm{\mu }}$ m paint thickness are compared with the VISAR experimental measurements performed by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5], whose good correlation reinforces the hypothesis that the paint layer does not need to be explicitly modelled in the mechanical finite element analysis. Also here, the presence of paint, and consequent enhancement of the contribution of the near-surface explosion of the LSP layer, leads to an increase of the displacement and velocity amplitudes, which keep increasing as the thickness of the paint layer increased from 200 to 400 ${\rm{\mu }}$ m. This again demonstrates, now for a higher electrically conducting LSP layer, that implicitly accounting for the paint layer in the formulation of the mechanical lightning loads allows capturing the trends observed in experiments [Reference Bigand and Duval2, Reference Murillo, Flourens, Garcia, Duval and Cavaliere7].

Figure 14. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for the unpainted case of ECF 195.

Figure 15. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for a paint thickness of 200 ${\rm{\mu }}$ m and ECF 195.

Figure 16. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for a paint thickness of 300 ${\rm{\mu }}$ m and ECF 195, and comparison with VISAR experimental measurements [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5].

Figure 17. Finite element predictions of the displacement and velocity histories on the five points identified in Fig. 2 for a paint thickness of 400 ${\rm{\mu }}$ m and ECF 195.

Figures 18–21 show, respectively, the predictions of longitudinal, transverse, shear and interlaminar damage on ECF 195. As in ECF 73 (Figs 8–11), no damage is observed in the unpainted panels, while an increasing extent of intralaminar and interlaminar damage with increasing paint thickness is observed in the panels with paint.

Figure 18. Longitudinal damage ( ${d_1}$ ) extent considering the effect of paint thickness for ECF 195. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

Figure 19. Transverse damage ( ${d_2}$ ) extent considering the effect of paint thickness for ECF 195. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

Figure 20. Shear damage ( ${d_6}$ ) extent considering the effect of paint thickness for ECF 195. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

Figure 21. Interlaminar damage ( $d$ ) extent considering the effect of paint thickness for ECF 195. The 0° direction is parallel to the (horizontal) X-axis. The layer count starts from the bottom of the laminate (Layer 1) to the top of the laminate (Layer 8). The latter is immediately below the LSP layer, where the VDLOAD pressure profiles are applied.

8.0 Concluding remarks

The models proposed in this work predict the faster deflection and the increase of damage extent with the increase of paint thickness so far observed experimentally [Reference Bigand and Duval2, Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5, Reference Murillo, Flourens, Garcia, Duval and Cavaliere7] without fitting parameters calibrated based on test data for different paint thicknesses. Moreover, they satisfactorily capture the effect of paint thickness for LSP layers of different areal weight, following the experimentally observed trends reported in literature [Reference Bigand and Duval2, Reference Lepetit, Escure, Guinard, Revel, Peres and Duval5, Reference Murillo, Flourens, Garcia, Duval and Cavaliere7]. It is important to mention that the possibility of having a tool that is capable of predicting the effect of paint thickness on the damage induced by lightning strikes could support the assessment of the thickness of the paint layer that minimises mechanical damage induced by lightning strikes, at the same time it complies with the environmental protection requirements. Improvements to the models presented in this study include the addition of strain-rate and temperature dependency in the constitutive models, which are expected to increase the level of fidelity of these simulations. Moreover, building upon the development and validation of appropriate models for the thermal, thermal-mechanical and mechanical lightning loads (the latter being addressed in the present work) and for the composite material behaviour could ultimately support the qualification and certification by analysis of new LSP technologies [Reference Kumar, Yokozeki, Karch, Hassen, Hershey, Kim, Lindahl, Barnes, Bandari and Kunc8] and of protected (generally painted) composite aero-structures subjected to lightning strike.