2.1 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. It is assumed that the spatial distribution of the expanding plasma pressure can be approximated by a Gaussian bell-shaped curve:

(4) \begin{align}{p_{\textrm{e}}}\left( {r,t} \right) = {p_{\textrm{e}}}\left( t \right){\textrm{ }}{k^{ - 2}}{\textrm{ exp}}\left[ { - \frac{{{r^2}}}{{{{\left( {k{\textrm{ }}{R_R}\left( t \right)} \right)}^2}}}} \right]\end{align}

where the effective plasma overpressure ${p_{\textrm{e}}}\left( t \right)$ is numerically determined separately for the two ECFs considering the corresponding root radii and the transient current waveform D with a time to peak of 26.8 $\mu $ s, a decay time of 48.5 $\mu $ s and a peak current of 96.4 kA [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44], as described by Karch et al. [Reference Karch, Arteiro and Camanho5]. The parameter $k \lt 1$ in Equation (4) is used to limit the pressure distribution mainly to the arc root area $\pi R_R^2\left( t \right)$ . A sensitivity study on the influence of this parameter shows that it affects the predicted velocity and displacement fields in the first 80 $\mu $ s, but the displacements after that point, including their peak values, depend mostly on the effective overpressure ${p_{\textrm{e}}}\left( t \right)$ numerically determined using the one-dimensional shock peening model [Reference Karch, Arteiro and Camanho5], namely when $k$ is not too small (i.e., $k\ge 0.1$ ). Based on the available experimental data [Reference Karch, Arteiro and Camanho5, Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44] and the sensitivity analysis, this parameter is set to 0.4 and to 0.1 for ECF 73 and ECF 195, respectively, suggesting differences on the spatial pressure distribution, within the arc root area, caused by the explosion of the LSP layers with 73.3 g/m² and 195.3 g/m² surface weight. This is a limitation of the one-dimensional shock peening model of the overpressure caused by the explosion of the LSP layers [Reference Karch, Arteiro and Camanho5, Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44], highlighting the need to develop a hemispherical shock peening model, currently not available, to account not only for the time evolution of the overpressure, but also for its spatial distribution.

The generated effective plasma overpressure considering protected CFRP with a paint thickness of 300 $\mu $ m (covering the LSP layer) is shown in Fig. 4. The resulting mechanical impulses due to the near-surface explosion of the LSP layer are shown in Fig. 5 for CFRP protected with ECF 73 and in Fig. 6 for CFRP protected with ECF 195. As can be observed, for both LSP layers, the near-surface explosion is the transient mechanical load with the highest contribution.

Figure 4. Calculated overpressure caused by explosion of the LSP layers (ECF 73 and ECF 195) as a function of time considering a transient current waveform D with a time to peak of 26.8 ${{\mu }}$ s, a decay time of 48.5 ${{\mu }}$ s and a peak current of 96.4 kA.

Figure 5. Calculated mechanical impulses for the three main contributions to lightning-induced mechanical loads on CFRP protected with ECF 73 as a function of time considering a transient current waveform D with a time to peak of 26.8 ${{\mu }}$ s, a decay time of 48.5 ${{\mu }}$ s and a peak current of 96.4 kA.

Figure 6. Calculated mechanical impulses for the three main contributions to lightning-induced mechanical loads on CFRP protected with ECF 195 as a function of time considering a transient current waveform D with a time to peak of 26.8 ${{\mu }}$ s, a decay time of 48.5 ${{\mu }}$ s and a peak current of 96.4 kA.

2.2 Supersonic plasma expansion

2.2.1 Strong shock wave approximation

A supersonic shock wave is created by the overpressure caused by the rapid temperature rise of the plasma channel. Following Lin [Reference Lin45], the cylindrical shock wave expands radially according to the following time varying shock wave radius [Reference Karch, Arteiro and Camanho5]:

(5) \begin{align}R\left( t \right) \cong \left\{ {\begin{array}{*{20}{l}}{1.004{{\left( {\frac{{{E_0}}}{{{\rho _0}}}} \right)}^{0.25}}\sqrt {0.5} } {}{,{\textrm{ }}t \lt 0.5 \mu s} \\[8pt]{1.004{{\left( {\frac{{{E_0}}}{{{\rho _0}}}} \right)}^{0.25}}\sqrt t } {}{,{\textrm{ }}t \ge 0.5 \mu s} {}{}\end{array}} \right.\end{align}

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_{{\textrm{max}}}}$ (A) can be established [Reference Karch, Schreiner, Honke and Wolfrum4, Reference Karch, Wulbrand and Müller46]:

(6) \begin{align}{E_0} \cong 0.45 \times {10^{ - 2}}{\textrm{ }}{\left( {{I_{{\textrm{max}}}}} \right)^{1.25}}\end{align}

Finally, the shock wave overpressure, which only depends on the peak current ${I_{{\textrm{max}}}}$ of the waveform, is described by:

(7) \begin{align}p\left( {\eta ,t} \right) \cong \left\{ {\begin{array}{*{20}{l}}{\beta \cdot t} {}{,{\textrm{ }}t \lt 0.5 \mu s} \\[5pt]{0.180{{\left( {{\rho _0}{E_0}} \right)}^{0.5}}f\left( \eta \right){\textrm{ }}{t^{ - 1}} - 1} {}{,{\textrm{ }}t \ge 0.5 \mu s} \\[5pt]0 {}{,{\textrm{ }}t \gt 100 \, \mu s} {}{}\end{array}} \right.\end{align}

where $\beta = 392.2952$ bar/ $\mu $ s is determined to approximately match the spatially uniform pressure distribution given by the strong shock wave approximation when $t = 0.5$ $\mu $ s for a current waveform with a peak current ( ${I_{{\textrm{max}}}}$ ) of 96.4 kA, $f\left( \eta \right)$ is a dimensionless pressure function and $\eta = r/R\left( t \right)$ is the dimensionless radius [Reference Karch, Wulbrand and Müller46]. Figures 7 and 8 show, respectively, the computed shock radius and the pressure profile evolution with time. The resulting mechanical impulse due to the supersonic plasma expansion is shown in Figs. 5 and 6 (noting that the supersonic plasma expansion overpressure is independent of the LSP layer used). As can be observed, although, in this case, the contribution of the supersonic plasma expansion is lower than the contribution of the near-surface explosion, it cannot be neglected.

Figure 7. Supersonic plasma expansion shock radius considering a transient current waveform D with a peak current of 96.4 kA.

Figure 8. Supersonic plasma expansion pressure profile evolution of the strong shock wave approximation and comparison with the ideal gas, Doan-Nickel and Plooster EOS, at 5 ${{\mu }}$ s, 10 ${{\mu }}$ s, 20 ${{\mu }}$ s, 50 ${{\mu }}$ s and 100 ${{\mu }}$ s from the highest (top left) to the lowest (bottom) pressure profiles, respectively.

2.3 Magnetic forces

As explained by Karch et al. [Reference Karch, Arteiro and Camanho5], the lightning current flows from the plasma channel into the LSP layer through the arc root surface. Therefore, the current flow in the LSP layer creates a magnetic pressure that depends on the region of the sample under ( $0 \lt r \lt {R_R}\left( t \right)$ ) or beyond the arc root radius ( $r \gt {R_R}\left( t \right)$ ). Taking this into account, the effective magnetic surface pressure acting on the upper side of the LSP layer is given as:

(8) \begin{align}{p_{\textrm{m}}}\left( {r,t} \right) = \left\{ {\begin{array}{*{20}{l}}{\frac{{{\mu _r}{\textrm{ }}{\mu _0}{\textrm{ }}I{{(t)}^2}}}{{8{\pi ^2}}}\frac{{{r^2}}}{{R_R^4\left( t \right)}}} {}, {}{{\textrm{ }}0 \lt r \lt {R_R}\left( t \right)} \\[5pt]{\frac{{{\mu _r}{\textrm{ }}{\mu _0}{\textrm{ }}I{{(t)}^2}}}{{8{\pi ^2}{r^2}}}} {}, {}{{\textrm{ }}r \gt {R_R}\left( t \right)}\end{array}} \right.\end{align}

where ${\mu _0} = 4\pi \times {10^{ - 7}}$ H $ \cdot $ m ${{\textrm{ }}^{ - 1}}$ is the permeability of the free space, ${\mu _r} = 1$ is the relative permeability of the non-magnetic LSP layer, and $I\left( t \right)$ is the total lightning current amplitude, here approximated by the modified double exponential function in Equation (2) [Reference Karch, Arteiro and Camanho5]. The magnetic pressure distribution is shown in Fig. 9. The resulting mechanical impulses due to the magnetic forces are shown in Fig. 5 for CFRP protected with ECF 73 and in Fig. 6 for CFRP protected with ECF 195. As can be observed, although, in this case, for both LSP layers the contribution of the magnetic forces is the lowest, it cannot be neglected.

Figure 9. Maximum magnetic pressure evolution with time for two LSP layers (ECF 73 and ECF 195) considering a transient current waveform D with a time to peak of 26.8 ${{\mu }}$ s, a decay time of 48.5 ${{\mu }}$ s and a peak current of 96.4 kA.

3.1 Finite element model

Following Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44], CFRP samples made from 8 T700/M21 UD plies with a quasi-isotropic [45/0/–45/90]_S lay-up are considered in this study. Each T700/M21 UD ply is 0.262 mm thick. Samples protected using ECF 73 and ECF 195 were tested, and each sample had a polyurethane paint layer with a thickness of 300 $\mu $ m (including surfacing film).

Square laminated plates, 450 mm $ \times $ 450 mm, were tested, supported on a fixed plate with a circular opening of 340 mm diameter, and fixed using screws placed along a circumference of 370 mm diameter concentric with the central opening of the support plate (Fig. 10a). Here, the full CFRP sample is modelled, fixing the out-of-plane displacement on the supported area outside the circular opening and fixing the nodes along the circumference where the screws are placed (Fig. 10b).

A fine mesh is used at the lightning strike location (Fig. 10b). Due to the sample size, a progressively coarser mesh was used as the distance to the centre of the sample increased. However, due to the small characteristic size (approximately 0.3 mm) of the finite elements in the refined region of the finite element model, the stable time increment of the explicit integration scheme is very small, and it tends to decrease with the simulation due to the collapse of a small number of volumetric elements below the attachment point of the lightning strike. To keep the computational cost of the explicit simulations reasonable, the selective mass scaling procedure available in Abaqus/Explicit [56] is used to limit the stable time increment. Because this procedure applies mass scaling just to the elements whose stable time increment is below a specified value, it was possible to limit the relative mass change to just 2 $ \times $ 10 ${{\textrm{ }}^{ - 9}}$ , while achieving a reasonable computational cost, with negligible changes in the mechanical response of the plates subjected to lightning strike.

A modification [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57] to the continuum damage mechanics model proposed by Maimí et al. [Reference Maimí, Camanho, Mayugo and Dávila58, Reference Maimí, Camanho, Mayugo and Dávila59] is used to represent failure of the composite plies. This modified version of the continuum damage mechanics model has a decoupled longitudinal tensile and compressive elastic behaviour that accounts for different stiffness in tension and in compression parallel to the fibre direction, and it uses bilinear softening laws to model longitudinal damage growth in tension and in compression, including a residual stress plateau in compression to more accurately capture crushing effects [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57]. Tables 3 and 4 show, respectively, the UD ply elastic properties and the strength and material degradation properties used in the simulations. ${S_{LP}}$ and ${K_P}$ are the in-plane shear yield stress and shear modulus degradation factor, ${X_T}$ , ${X_C}$ , ${f_{XT}}$ and ${f_{XC}}$ are the longitudinal tensile and compressive plain strengths and the pull-out (tension) and crushing (compression) strength factors, ${Y_T}$ and ${Y_C}$ are the UD transverse tensile and compressive strengths, ${\alpha _0}$ is the composite fracture angle in uniaxial transverse compression, ${S_L}$ is the UD in-plane shear strength, ${{\mathcal G}_{1 + }}$ and ${f_{{\mathcal G}1 + }}$ are the longitudinal tensile fracture toughness and the fraction of the fracture toughness dissipated due to fibre pull-out, ${{\mathcal G}_{1 - }}$ , ${{\mathcal G}_{2 + }}$ and ${{\mathcal G}_6}$ are the longitudinal compressive, transverse tensile and in-plane shear fracture toughness, ${f_{fxc}}$ is a frictional parameter used to scale ${f_{XC}}$ as a function of the hydrostatic pressure, and ${\eta _G}$ is a material dependent empirically derived enhancement factor used in the definition of an effective fracture toughness for in-plane shear in the presence of transverse compressive stresses [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57]. All material properties are independently determined; the elastic and strength properties follow standard test methods [Reference Charrier60–Reference Vandellos63], while the material degradation properties follow the approaches detailed by Furtado et al. [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57].

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

Figure 10. (a) Configuration of the lightning strike test support and sample. (b) Composite plate finite element mesh.

Table 4. T700/M21 strengths and material degradation properties [Reference Charrier60–Reference Vandellos63] used in the modified continuum damage mechanics model [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57]

The interaction properties with cohesive behaviour available in Abaqus [56] are used to model delamination, considering a bilinear traction-separation damage law. Mode-dependent strengths and fracture energies are considered. The onset of interlaminar damage is predicted by a quadratic stress-based criterion. Table 5 shows the interlaminar properties used in the simulations, where $K$ is the stiffness of the cohesive surfaces before damage onset, ${{\mathcal G}_{Ic}}$ and ${{\mathcal G}_{IIc}}$ are respectively the mode I and mode II fracture toughness, ${\eta _{{\textrm{B}} - {\textrm{K}}}}$ is the mixed-mode parameter of the B-K law [Reference Benzeggagh and Kenane64], and ${\tau _3}$ the normal strength. The shear strengths ${\tau _1} = {\tau _2}$ are determined as a function of the mode I and mode II fracture toughness and the normal strength ${\tau _3}$ as [Reference Turon, Camanho, Costa and Renart65]:

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

Table 5. T700/M21 interlaminar properties

To account for friction between delaminating interfaces subjected to sliding movements, tangential friction is assigned to the cohesive surface interactions. Following previous work with a similar modelling strategy [Reference Arteiro, Gray and Camanho66], a friction coefficient of 0.3 is assumed.

The LSP layers are assumed homogeneous orthotropic linear-elastic solids. Following Karch et al. [Reference Karch, Arteiro and Camanho5], the effective properties of ECF 73 and ECF 195 filled with epoxy resin M21 were determined from a micro-mechanical finite element homogenisation approach using the commercial finite element package ANSYS. Periodic boundary conditions were applied to a representative volume element (RVE) of a unit cell of the metallic mesh filled with resin by coupling nodes of the finite elements on the opposite boundary surfaces of the unit cell. Three uniaxial tensile strain loads and three shear strain loads ( ${{\bf{E}}_j}$ ) are solved to obtain the full stiffness matrix. The macroscopic stresses ( ${{\bf{\Sigma }}_i}$ ) are obtained by averaging the microscopic stress within the RVE, and the $j$ th column of the effective stiffness tensor is obtained from:

(10) \begin{align}C_{ij}^{\textrm{*}} = \frac{{{{\bf{\Sigma }}_i}}}{{{{\bf{E}}_j}}}\end{align}

The properties obtained for both LSP layers are given in Table 6. A similar numerical finite element approach can be used to determine the effective electrical and thermal conductivities of orthotropic hybrid ECF/epoxy resin structures. For more details about the homogenisation approach, the reader is referred to Karch et al. [Reference Karch, Arteiro and Camanho5].

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

Following Karch et al. [Reference Karch, Arteiro and Camanho5], the mechanical forces resulting (i) from the near-surface explosion of the LSP layer, (ii) from the supersonic plasma expansion and (iii) from the magnetic field caused by the impressed current flow in the electrically conducting structures for a transient current waveform D with a time to peak of 26.8 $\mu $ s, a decay time of 48.5 $\mu $ s and a peak current of 96.4 kA [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44], presented in Section 2, are implemented in the finite element software Abaqus/Explicit [56] using the user-defined subroutine VDLOAD. The use of the explicit solver Abaqus/Explicit [56], on one hand, facilitates the solution of the continuum damage mechanics model [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57], implemented in a user-defined subroutine VUMAT, circumventing convergence issues, typical in this type of models, without the need to implement viscous regularisation procedures. On the other hand, the adoption of an explicit dynamics solution is compatible with the highly transient and short time lightning strike problem considered in this study.

To implement the overpressure of the supersonic plasma expansion given by the CFD calculations (Fig. 8), the CFD results are exported to data files as described in Section 2.2.2. For each time increment of the finite element structural simulation of the lighting strike, on each node of the surface where the pressure is applied, the VDLOAD subroutine reads the 2,500 $ \times $ 701 results data file (Section 2.2.2), and the overpressure of the supersonic plasma expansion at the current time and location is determined by bilinear interpolation between the available position and time data points extracted from the CFD analysis.

As discussed by Karch et al. [Reference Karch, Arteiro and Camanho5], the velocity interferometer system for any reflector (VISAR) deflection measurements conducted by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44] strongly indicate that the lightning strike did not initiate from the middle of the CFRP samples. Therefore, Karch et al. [Reference Karch, Arteiro and Camanho5] used a 2D symmetric Gauss function to fit the deflection in the measurement points 1–5 (Fig. 10) for different times. Karch et al. [Reference Karch, Arteiro and Camanho5] then averaged the time dependent offset coordinates, obtaining the values given in Table 7. The mechanical pressures described in Section 2.0 are centred at these points.

Table 7. Offset coordinates for ECF 73 and ECF 195

3.2 Finite element results

The results are presented for the two CFRP samples, protected with ECF 73 and with ECF 195. For each LSP layer, the displacement and velocity results at the five VISAR points (1–5 in Fig. 10) are presented, followed by layer-by-layer damage maps, to assess the damage induced by the simulated lightning strike through the thickness of the laminate.

3.2.1 ECF 73 LSP layer

Figure 11 shows the finite element prediction of the deformed shape obtained with the plasma expansion shock wave from the CFD calculations based on Plooster EOS. Figure 12 compares the predictions using the analytical strong shock wave approximation and the CFD data based on Plooster EOS. Similar results are obtained for the remaining EOS. In all cases, the displacements and velocities are very similar in the initial 100 $\mu $ s. Slightly larger displacements are then predicted for time periods above 100 $\mu $ s when considering the solutions based on CFD data to compute the plasma expansion shock wave.

Figure 11. Deformed shape of the bottom layer (‘Layer 1’) along a centre line of the laminated plate for ECF 73 within a time ${\textrm{t}} \le 500$ ${{\mu }}$ s (colour bar) obtained using the plasma expansion shock wave from the CFD calculations based on Plooster EOS.

Figure 12. Finite element predictions of the displacements and velocities at the five VISAR points (1–5 in Fig. 10) obtained using the plasma expansion shock waves from the analytical approximation (solid lines) and from the CFD calculations based on Plooster EOS (dashed lines) for ECF 73.

Figures 13–16 show the predicted layer-by-layer damage maps for ECF 73. Here, the onset and propagation of damage can be assessed, considering the intralaminar damage variables in the longitudinal direction ( ${d_1}$ ), transverse direction ( ${d_2}$ ) and shear ( ${d_6}$ ), following Furtado et al. [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57], and interlaminar damage ( $d$ ) [Reference Turon, Camanho, Costa and Renart65]. No major differences between the predictions obtained using the plasma expansion shock waves from the strong shock wave approximation and from the CFD data are observed, independently of the EOS used in the CFD calculations.

Figure 13. Longitudinal (fibre) damage maps ( ${{\textrm{d}}_1}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 73. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

Figure 14. Transverse (matrix) damage maps ( ${{\textrm{d}}_2}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 73. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

Figure 15. Shear (matrix) damage maps ( ${{\textrm{d}}_6}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 73. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

Figure 16. Interlaminar damage maps ( ${\textrm{d}}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 73. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

3.2.2 ECF 195 LSP layer

Figure 17 presents the finite element prediction of the deformed shape obtained with the plasma expansion shock wave from the CFD calculations based on Plooster EOS when the CFRP laminate is protected by ECF 195. Figures 18 and 19 compare the predictions using the analytical strong shock wave approximation and the CFD data based, respectively, on the EOS with the ideal gas assumption and on Plooster EOS. The results of the latter are similar to those obained using the CFD data based on Doan-Nickel EOS. While slightly lower displacements are obtained in the initial 300 $\mu $ s with the CFD data based on the EOS with the ideal gas assumption, slightly larger displacements are predicted for time periods above 100 $\mu $ s when employing the plasma expansion shock wave from the CFD calculations based on Doan-Nickel and Plooster EOS, compared with the solution from the analytical strong shock wave approximation. These differences, however, have just a mild impact on the predicted damage footprints (Figs. 20–23).

Figures 20–23 show the predicted layer-by-layer damage maps for ECF 195, where the onset and propagation of damage can be assessed, considering the intralaminar damage variables in the longitudinal direction ( ${d_1}$ ), transverse direction ( ${d_2}$ ) and shear ( ${d_6}$ ), following Furtado et al. [Reference Furtado, Catalanotti, Arteiro, Gray, Wardle and Camanho57], and interlaminar damage ( $d$ ) [Reference Turon, Camanho, Costa and Renart65]. As expected, the damage extent on the CFRP laminate protected by ECF 195 is lower than that protected by the lighter ECF 73 (Figs. 13–16). Yet, no major differences between the predictions obtained using the plasma expansion shock waves from the strong shock wave approximation and from the CFD data are observed, similarly to what was observed for ECF 73 (Figs. 13–16).

Figure 17. Deformed shape of the bottom layer (‘Layer 1’) along a centre line of the laminated plate for ECF 195 within a time ${\textrm{t}} \le 500$ ${{\mu }}$ s (colour bar) obtained using the plasma expansion shock wave from the CFD calculations based on Plooster EOS.

Figure 18. Finite element predictions of the displacements and velocities at the five VISAR points (1–5 in Fig. 10) obtained using the plasma expansion shock waves from the analytical approximation (solid lines) and from the CFD calculations based on the EOS with the ideal gas assumption (dashed lines) for ECF 195.

Figure 19. Finite element predictions of the displacements and velocities at the five VISAR points (1–5 in Fig. 10) obtained using the plasma expansion shock waves from the analytical approximation (solid lines) and from the CFD calculations based on Plooster EOS (dashed lines) for ECF 195.

Figure 20. Longitudinal (fibre) damage maps ( ${{\textrm{d}}_1}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 195. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

Figure 21. Transverse (matrix) damage maps ( ${{\textrm{d}}_2}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 195. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

Figure 22. Shear (matrix) damage maps ( ${{\textrm{d}}_6}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 195. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

Figure 23. Interlaminar damage maps ( ${\textrm{d}}$ ) considering the plasma expansion shock waves from the strong shock wave approximation and from the CFD calculations for ECF 195. The 0 ${{\textrm{ }}^ \circ }$ 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. The whole square layers (450 mm-long sides) are shown. All images are presented on the same scale to facilitate the comparison between the extent of the different damage mechanisms.

3.3 Model validation

The predictions of the mechanical response obtained with the model proposed in this work using the plasma expansion shock waves obtained from the strong shock wave approximation is validated against direct measurements of the out-of-plane displacement and velocity of coupons tested by Lepetit et al. [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44]. The comparison with the experimental VISAR measurements [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44] for ECF 73 and ECF 195 is shown in Figs. 24–25, respectively.

Figure 24. Finite element predictions of the displacements and velocities at the five VISAR points (1–5 in Fig. 10) obtained using the plasma expansion shock waves from the analytical approximation (solid lines) and corresponding experimental VISAR measurements (dashed lines) for ECF 73 [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44].

Figure 25. Finite element predictions of the displacements and velocities at the five VISAR points (1–5 in Fig. 10) obtained using the plasma expansion shock waves from the analytical approximation (solid lines) and corresponding experimental VISAR measurements (dashed lines) for ECF 195 [Reference Lepetit, Escure, Guinard, Revel, Peres and Duval44].

As can be observed, the predictions of the out-of-plane velocity represent very well the experimental results at the different measurement points, in particular the peak velocities. Regarding the predicted deflections, although they exceed the experimental results at the initial 50–300 $\mu $ s, the results then converge at around 500 $\mu $ s.

Remarkable improvements can be observed with respect to what was previously obtained by Karch et al. [Reference Karch, Arteiro and Camanho5] with shell elements in terms of velocity and deflection predictions. On the other hand, the results in Figs. 24–25 also demonstrate that the models of the lightning loads proposed by Karch et al. [Reference Karch, Arteiro and Camanho5] can accurately represent the mechanical contributions from lightning strike, in particular when used together with detailed 3D damage models for CFRPs, as employed in this work.

Figures 13–16 and Figs. 20–23 show the contour plots of the damaged elements (intralaminar damage ${d_1}$ , ${d_2}$ and ${d_6}$ ) and interfaces (interlaminar damage $d$ ) corresponding to damage onset (partially damaged elements/interfaces, with $0\leqslant{d_i}\lt1, i = 1,2,6$ and $0\leqslant{d}\lt1$ , in black) and failure (fully damaged elements, with ${d_i} \approx 1, i = 1,2,6$ and $d \approx 1$ , in red) for the strong shock wave approximation. Due to the progressively coarser mesh from the centre of the sample, the damage onset maps are less accurate in the coarse regions. This can lead to extensive areas of damage onset, e.g., for ECF 73 (Figs. 14–15) closer to the supporting edge. Therefore, any quantitative measure of damage extent would be affected by some level of subjectivity. Nevertheless, the damage patterns are those usually observed on CFRP samples subjected to simulated lightning strikes and attributed to mechanical effects: damage tends to concentrate below the region of the initial attachment, mainly interlaminar and matrix damage at the middle and bottom plies, and to spread away from the initial attachment point predominantly at the top plies. Although the extent of interlaminar damage ( $d$ ) predicted by the models (which is almost negligible) is substantially lower than the total delaminated area identified by C-scan (whose pictures are not publicly available, but were considered for comparison and model validation), the combined extent of matrix-dominated damage ( ${d_2}$ , ${d_6}$ and $d$ ) below the second ply follows the trends observed experimentally: damage extends until the few rear plies, and ECF 73 leads to a total damaged area below the second ply approximately three times larger than ECF 195. It is important to stress that it is often assumed that damage at the interface between the first and second plies in the tests is essentially surface damage, and thus, mostly affected by thermally induced effects. For that reason, it is common practice in industry to neglect the damaged area at the first interface when determining the damaged area by C-scan after testing, thus considering only the damaged area below the second ply.