1 Introduction

Due to their high strength-to-weight ratio, carbon-fibre reinforced composite laminates, using a thermosetting epoxy polymer matrix, have been widely used in the manufacture of structural components for commercial aircraft and vehicles [1,2,3]. However, their wider application in primary load-carrying structures has been hindered by their reduced impact performance compared to that of more traditional metallic materials [4,5,6]. In particular, impact loading can lead to damage which may not be necessarily visible but which may, nonetheless, lead to a significant strength reduction [7,8,9,10]. To enhance confidence when using composite materials in the latest generation of aircraft and vehicles a better understanding of the impact behaviour of composite materials is needed. Driven by this requirement, considerable efforts have been made by both academia and industry to study and improve the impact performance of such composite laminates e.g. [11,12,13,14,15,16,17,18,19,20,21,22,23,24].

Experiments are an essential way to understand the impact behaviour of composite laminates but only a relatively small number of comparative studies have been reported on the effects of the shape of the impactor head. Mitrevski et al. [17] have reported the effect of impactor shape on the impact behaviour of woven carbon-fibre/epoxy-matrix composites. Their studies involved the low-velocity impact of the composites using 12 mm diameter hemispherical, ogival and conical head-shaped impactors at impact energies of 4 J and 6 J. They found that the energy absorbed by the specimen was the highest for the conical impactor which also produced the largest penetration depth. However, the peak load was greatest for the hemispherical impactor which also produced the shortest contact duration. In a subsequent paper [18] it was reported that different impactor head-shapes produced different levels of the various types of damage that were observed, such as fibre breakage, matrix cracking and delamination, which affected the residual properties of the composite. It was found that the relatively blunt hemispherical impactor produced the largest damage area and the type of damage was dominated by delamination. Robinson and Davies [19] have reported results from drop-weight tests using a hemispherically-shaped impactor at impact velocities of up 6.0 m/s striking woven glass- and carbon-fibre composites based on a matrix of an acrylic polymer. They found that the impact damage caused was mainly a function of the impact energy, and not the velocity or mass, separately, of the impactor. Further, the measured peak load was also shown to correlate with the impact energy independent of the impactor mass used, indicating that the low-velocity impact of these specimens was a quasi-static process. Their test results also indicated that the peak load versus damage size relationship was independent of the diameter of the impactor. More recently, Sevkat et al. [20] have undertaken drop-weight tests at an impact velocity of 4.1 m/s on woven glass-fibre/carbon-fibre hybrid epoxy composites. They reported that the contact area between the impactor and the composite test specimen was a critical factor in the response of the composites. Impactors with relatively large contact surfaces produced a higher initial load peak, a higher maximum load, a shorter contact duration for the time-scale of the impact event and more extensive delaminations between the hybrid layers.

Another effective approach to investigate the impact response of composite laminates is by developing analytical and numerical models. Analytical models are more suited to understanding the physics of the impact effects and have better computational efficiency, since they tend to be relatively tractable using simple geometries and boundary conditions. Numerical models are more suitable when trying to represent reality more closely, which should also lead to a better understanding of the impact behaviour of composites. For example, to simulate the behaviour of composite laminates under a large mass, low-velocity impact, Olsson proposed a simplified analytical model to predict the initiation and propagation of impact damage in composites [21, 22]. Analytical results were compared with experimental results obtained from composites with different lay-ups, geometries and boundary conditions. It was found that the critical load for delamination growth was insensitive to the specimen geometry and boundary conditions. A relatively high computational efficiency, for the level of accuracy achieved, was demonstrated by the analytical model when it was employed to predict the impact behaviour of composite laminates which did not experience extensive fibre breakage. On the other hand, Bouvet et al. [23, 24] have used a numerical model to capture the different types of damage within composite laminates subjected to a relatively low-velocity/low-energy impact. In their numerical model, three types of damage, namely fibre failure, matrix cracking and delamination, were considered. Delamination was conventionally modelled using interface elements based on a fracture mechanics approach. Fibre failure was simulated using degradation of the volumetric elements. For simulating intralaminar matrix cracking, a series of longitudinal strips, with one volume element through the ply thickness, were modelled and connected using zero-thickness interface elements.

In the present work, a detailed study of the impact behaviour of composite laminates was conducted using both experimental and numerical methods. For the experimental studies, panels of carbon-fibre/epoxy-matrix composite laminates, with a lay-up of [+453/03/−453/03]s, were manufactured. The composite specimens, produced from the composite panels, were then employed to perform the drop-weight impact experiments. The specimen dimensions and experimental set-up described in ASTM D-7136 standard [25] were adopted. Hemispherical steel impactors or flat-ended steel impactors were used to strike the composite specimens with 15 J of energy. Subsequently, all the tested specimens were examined using ultrasonic C-scanning to evaluate the extent of damage due the different impactors. For the numerical studies, a detailed finite-element (FE) model was developed within the commercial software code ‘Abaqus’ to simulate the drop-weight impact tests. This approach follows on from previous research performed at Queen’s University Belfast, in the Advanced Composites Research Group, and a novel aspect of the current research is to enhance the previous numerical analyses with an elastic-plastic model. The validated FE model, incorporating the elastic-plastic material response, was employed (a) to determine the damage distribution and deformation in the composite specimen and (b) to gain further insights into the effects of the shape of the impacting projectile on the impact behaviour of composite laminates.

2 Impactors and Specimens

2.1 Impactors

Two types of impactors with different shaped heads were used: a hemispherical steel (i.e. a round-nosed steel (RNS)) impactor and a flat-ended (i.e. a flat-faced steel (FFS)) steel impactor. The hemispherically-headed steel impactor and the flat-ended steel impactor were manufactured from stainless steel and are illustrated in Fig. 1. They had masses of 5.20 and 5.25 kg, respectively, and the impact velocities used were 2.40 m/s and 2.39 m/s respectively. The impact energy was 15 J in both cases. The impactor heads had a diameter of 16 mm, with the head of the flat-ended impactor having a chamfer of about 45o around its periphery, as shown in Fig. 1. The diameter of the main body of the impactor was 20 mm and the load cell, with a data sampling rate of 500 kHz, was located in the forward section of the impactor.

Fig. 1
figure 1

Illustrations of: (a) the hemispherically-headed (RNS) steel and (b) the flat-ended steel (FFS) impactor with a 45o chamfer on the edge to the flat-ended cylindrical section

2.2 Composite Specimens

Composite laminates were manufactured from unidirectional (T700) carbon-fibre reinforced epoxy-matrix pre-pregs. The stacked prepregs, provided by Beian Ltd., China, were cured in an autoclave. A schematic of the vacuum-bagging procedure for the autoclave and the curing cycle for the T700/epoxy pre-pregs (0.125 mm in thickness) are presented in Fig. 2. Composite test specimens for the drop-weight impact experiments were produced from the manufactured composite panels and the cured composite laminates had a thickness of 3 mm. All the composite specimens possessed the same lay-up of [+453/03/−453/03]s. (The reason for selecting this particular composite lay-up, with blocking of the plies, was to give more detectable delaminations, rather than to absorb more impact energy.) The dimensions of the composite specimens were 150 mm × 100 mm, as defined in the ASTM D7136 standard [25].

Fig. 2
figure 2

Diagrams of (a) the vacuumn-bag assembly and (b) the curing cycle for the T700/epoxy prepregs

3 Experimental Procedures

The drop-weight impact experiments were performed using a drop-tower system provided by Instron CEAST, Milan, Italy, as shown in Fig. 3. In the experiments the hemispherically-headed impactor and the flat-ended steel impactor were employed to strike the composite specimens with 15 J of energy. Each target composite specimen was placed on a steel picture-frame which had outer dimensions that matched those of the composite specimens and with a 125 mm × 75 mm cut-out. This assembly was clamped to the base of the drop-weight tower using four toggle clamps with rubber tips, which prevented slippage of the composite specimen during the impact test [25]. For each test, three replicate composite specimens were impacted and the signal from the load cell was not subjected to any software filtering. The contact time was determined from the load versus time trace, as recorded via the load cell in the impactor. The velocity at impact, to give an impact energy of 15 J, was pre-determined in the drop-weight tests, via the software supplied by CEAST, to an accuracy of two decimal places and a catching system was employed to arrest the impactor after a single impact. The displacement of the test specimen was calculated, again via the software supplied by CEAST, by double integration of the load versus time response from the load cell in the impactor using Newton’s second law. All the tested composite specimens were inspected using ultrasonic C-scanning, as shown in Fig. 4, to assess the extent of any impact-induced delamination damage. The C-scan equipment had a 5 MHz probe with a scanning width of 34 mm [26].

Fig. 3
figure 3

The set-up for the drop-weight experiment (left) and the target composite specimen (right)

Fig. 4
figure 4

The ultrasonic C-scan equipment for impact damage assessment

4 Experimental Results

4.1 Loading Response

The loading responses of the impacted composite specimens when struck by the hemispherical steel and the flat-ended steel are shown in Fig. 5. At the same impact energy level of 15 J, the average contact time, i.e. 5.8 ms, for the specimens impacted by flat-ended steel impactor is 0.9 ms (i.e. about 13.5%) shorter than when using the hemispherical steel impactor, where the average contact time is of 6.7 ms. An average maximum impact load of 7.1 kN was measured on the flat-ended steel impactor upon striking the composite specimens. This is about 50% higher than that for the hemispherical steel impactor, where the average maximum impact load was 4.7 kN. The maximum displacement measured for the composite test specimen struck with the hemispherical impactor is 6.0 mm. For the flat-ended impactor, the maximum displacement is reduced to 4.7 mm. The reason for these differences in the values of the displacement is that the load introduced by the flat-ended impactor was distributed over a larger area. As a result, compared to the hemispherically-headed impactor, the flat-ended impactor is arrested more quickly and the initial effective stiffness of the specimen is higher, which leads to a shorter impact duration and lower specimen displacement.

Fig. 5
figure 5

Experimental results for (a) the load versus time and (b) the load versus displacement curves obtained from the impact experiments performed using the hemispherical steel (RNS) and the flat-ended steel (FFS) impactors. (Results from the three replicate tests are shown)

4.2 Impact Damage

For all the composite test specimens, no damage could be detected from a simple visual inspection of the specimens after the impact test. However, a slight indentation mark was apparent on the front face of all the test specimens at the impact site. From the C-scan tests for both types of impactor it was found that, since the composite has a lay-up of [+453/03/−453/03]s, there was a propensity for the delamination to grow at the interfaces between the 45o and 0o plies, with the orientation of the lower ply determining the direction of delamination propagation. Figure 6 shows the delamination maps, measured via the C-scan tests, obtained from the composite test specimens impacted by the hemispherical steel and flat-ended steel impactors. The average areas of the damage maps produced by the hemispherical steel impactor and flat-ended steel impactors are similar. However, the specimens struck with the hemispherical impactor show a continuous damage area, which is centred around the point of impact, whilst the specimens impacted with the flat-ended impactor exhibit two separate damage areas, with a central zone exhibiting no damage. These differences in the extent and shape of the damage are explained below from the results of the numerical modelling studies.

Fig. 6
figure 6

Interlaminar C-scan damage maps obtained from the specimens impacted with the hemispherical steel (RNS) and the flat-ended steel (FFS) impactors. (Results from the three replicate test specimens are shown. The 0o fibre direction is indicated. The right-hand side scale indicates the location of the delamination as a function of the depth through the thickness of the specimen, where the red colour is representative of the front (impacted) surface and the blue colour is representative of the rear (non-impacted) surface of the composite specimen)

5 Finite-Element (FE) Numerical Studies

5.1 The Damage Model

The implementation of the FE model is shown schematically in the flow charts in Fig. 7, for one computational time-step and for a single element, for modelling the interlaminar and intralaminar damage. Both the flowchart of the main model and the flow-chart for the elastic-plastic damage model are shown and these two models run simultaneously. The FE model would typically run from 0 to about 8 ms, with about 100 time-steps being employed and a computational time-step was performed for every appropriate single element in the FE model. The numerical model runs were stopped when the defined total time for the impact event had expired.

Fig. 7
figure 7

The implementation of the FE numerical model showing schematically the flow charts, for one computational time-step and for a single element, for modelling interlaminar and intralaminar damage. With both the flowchart for the main model and the flow-chart for the elastic-plastic damage model shown. The impact FE model would run typically from 0 to ca. 8 ms, with ca. 100 time-steps being employed

5.1.1 The Intralaminar Damage Model

The types of intralaminar damage that are typically observed in composite laminates consist of permanent deformation of the matrix, matrix cracking, fibre debonding and localised fibre failure.

The Initiation of Intralaminar Damage

In the intralaminar damage model that is proposed an extended elastic-plastic model was used to capture the material response prior to the initiation of matrix cracking, etc., since this enables a more accurate prediction of the impact behaviour of the composite laminate. The constitutive relation for the elastic-plastic model can be obtained by combining the classic elastic model with the extended plastic model and is given in [27, 28] as:

$$ \left\{\begin{array}{c}d{\varepsilon}_{11}\\ {}d{\varepsilon}_{22}\\ {}d{\varepsilon}_{33}\\ {}d{\varepsilon}_{12}\\ {}{d\varepsilon}_{13}\\ {}{d\varepsilon}_{23}\end{array}\right\}=\left[\begin{array}{cccccc}1/{E}_{11}& -{\nu}_{21}/{E}_{11}& -{\nu}_{31}/{E}_{11}& 0& 0& 0\\ {}-{\nu}_{12}/{E}_{22}& 1/{E}_{22}& -{\nu}_{32}/{E}_{22}& 0& 0& 0\\ {}-{\nu}_{13}/{E}_{33}& -{\nu}_{23}/{E}_{33}& 1/{E}_{33}& 0& 0& 0\\ {}0& 0& 0& 1/{G}_{12}& 0& 0\\ {}0& 0& 0& 0& 1/{G}_{13}& 0\\ {}0& 0& 0& 0& 0& 1/{G}_{23}\end{array}\right]\left\{\begin{array}{c}d{\sigma}_{11}\\ {}d{\sigma}_{22}\\ {}d{\sigma}_{33}\\ {}d{\sigma}_{12}\\ {}d{\sigma}_{13}\\ {}d{\sigma}_{23}\end{array}\right\}+\left\{\begin{array}{c}d{\varepsilon}_{11}^p\\ {}d{\varepsilon}_{22}^p\\ {}d{\varepsilon}_{33}^p\\ {}d{\varepsilon}_{12}^p\\ {}d{\varepsilon}_{13}^p\\ {}d{\varepsilon}_{23}^p\end{array}\right\} $$
(1)

where ij (i, j = 1, 2, 3) are the incremental total strain tensors and ij (i, j = 1, 2, 3) are the incremental stress tensors. The parameters νij (i, j = 1, 2, 3, i ≠ j) are the Poisson’s ratios, Eii (i, j = 1, 2, 3) are the Young’s moduli either for tension or compression loading, which are generally considered to be similar for composite laminates [29], and Gij (i, j = 1, 2, 3,  i ≠ j) are the shear moduli. The parameters \( d{\varepsilon}_{ij}^p\ \left(i,j=1,2,3\right) \) represent the incremental plastic strain tensors, which are related to the effective stress, σeff, and effective plastic strain,\( {\varepsilon}_{eff}^p \). The effective stress, σeff, is given by [30, 31]:

$$ {\sigma}_{eff}=\sqrt{\frac{3}{2}\left({\sigma}_{22}^2+{\sigma}_{33}^2\right)-3{\sigma}_{22}{\sigma}_{33}+3{a}_{66}\left({\tau}_{13}^2+{\tau}_{12}^2+{\tau}_{23}^2\right)} $$
(2)

where the value of a66 can be readily determined from off-axis experiments conducted at different values of the off-axis angle using a unidirectional composite [27, 31,32,33]. The relationship between the effective stress, σeff, and the effective plastic strain,\( {\varepsilon}_{eff}^p \), can be expressed as a power-law function, given by [31,32,33]:

$$ {\varepsilon}_{eff}^p=A{\sigma_{eff}}^n $$
(3)

where A and n are nonlinear coefficients, which are determined by fitting to the \( {\sigma}_{eff}\ \mathrm{versus}\ {\varepsilon}_{eff}^p \) data, again as obtained from the off-axis experiments when different off-axis angles are employed using a unidirectional composite [27, 31,32,33]. The determination of the parameter, a66, and the nonlinear coefficients, A and n, allow the calculation of the incremental plastic strain tensors, \( d{\varepsilon}_{ij}^p\left(i,j=1,2,3\right) \), as given by:

$$ \left\{\begin{array}{c}d{\varepsilon}_{11}^p\\ {}d{\varepsilon}_{22}^p\\ {}d{\varepsilon}_{33}^p\\ {}d{\varepsilon}_{12}^p\\ {}d{\varepsilon}_{13}^p\\ {}d{\varepsilon}_{23}^p\end{array}\right\}=\frac{An}{{\sigma_{eff}}^{n-1}}\left\{\begin{array}{c}0\\ {}3\left({\sigma}_{22}-{\sigma}_{33}\right)/2{\sigma}_{eff}\\ {}3\left({\sigma}_{33}-{\sigma}_{22}\right)/2{\sigma}_{eff}\\ {}3{a}_{66}{\tau}_{12}/2{\sigma}_{eff}\\ {}3{a}_{66}{\tau}_{13}/2{\sigma}_{eff}\\ {}3{a}_{66}{\tau}_{23}/2{\sigma}_{eff}\end{array}\right\}d{\sigma}_{eff} $$
(4)

The three-dimensional (3D) Northwestern University (NU) damage criteria were employed to capture the initiation of intralaminar damage such as matrix cracking, etc. These damage initiation criteria proposed by Daniel et al. [34,35,36] are partially interactive damage criteria, which means that more than one stress component may be used in them. Damage criteria were defined which represent the intralaminar damage that may be induced in a composite ply via: longitudinal tensile stresses, longitudinal compressive stresses, transverse tensile stresses, transverse compressive stresses, transverse shear stresses, through-thickness tensile stresses, through-thickness compressive stresses and through-thickness shear stresses. The mathematical details of the damage initiation model for a composite ply are given in Table 1, where σij and τ13 are the normal and shear stress components, respectively, and FiT (i = 1, 2, 3), FiC (i = 1, 2, 3) and FiS (i = 2, 3) are the tensile, compressive and shear damage criteria in the three material directions, respectively. The terms XT, YT and ZT denote the allowable tensile strengths in the three material directions. Similarly, XC, YC and ZC are the allowable compressive strengths in the three material directions. Finally, S12, S13 and S23 represent the allowable shear strengths in the corresponding material directions.

Table 1 The mathematical definitions for the intralaminar damage initiation model

The Propagation of Intralaminar Damage

Corresponding to the damage initiation mechanisms defined in Hashin’s damage criteria, eight damage parameters, d1t, d1c, d2t, d2cd2s, d3t, d3c and d3s may be defined to indicate the propagation of intralaminar damage in a composite ply, see Table 2. A general form of the damage variable, for a particular type of damage initiation, is given by:

$$ d=\frac{\varepsilon^f\left(\varepsilon -{\varepsilon}^0\right)}{\varepsilon \left({\varepsilon}^f-{\varepsilon}^0\right)} $$
(5)

where d denotes the general damage parameter. The strain, ε, is the equivalent strain in the composite ply. The strain values, ε0 and εf, are the equivalent strains corresponding to initial failure and final failure, respectively. The final failure strain tensor, \( {\varepsilon}_{ij}^f\ \left(i,j=1,2,3\right) \), can be determined by the following equation:

$$ {\varepsilon}_{ij}^f=2{G}_{cij}/\left({\sigma}_{ij}^0{l}_c\right) $$
(6)

where Gcij is the fracture toughness corresponding to the different principal materials directions, \( {\sigma}_{ij}^0 \) is the stress tensor corresponding to damage initiation and lc is the characteristic length which can be determined based on the volume of the elements. During the evolution of intralaminar damage, the elasticity matrix needs to be degraded to compute the values of the degraded stresses. To achieve this, four damage variables, d1, d2, d3 and ds, which reflect longitudinal damage, transverse damage and through-thickness damage and shear damage, respectively, were derived from the damage parameters, d1t, d1c, d2t, d2cd2sd3t, d3c and d3s, corresponding to the types of damage previously discussed, as follows:

$$ \mathrm{Longitudinal}\ \mathrm{damage}:{d}_1=\left\{\begin{array}{c}{d}_{1t},\kern0.5em \mathrm{and}\ {\upsigma}_{11}\ge 0\\ {}{d}_{1c},\kern0.5em \mathrm{and}\ {\sigma}_{11}\le 0\end{array}\right. $$
(7)
$$ \mathrm{Transverse}\ \mathrm{damage}:{d}_2=\left\{\begin{array}{c}{d}_{2t},\kern0.5em \left|{\sigma}_{22}\right|\ge \left|{\tau}_{12}\left({\tau}_{23}\right)\right|\ \mathrm{and}\ {\sigma}_{22}\ge 0\\ {}{d}_{2c},\kern0.5em \left|{\sigma}_{22}\right|\ge \left|{\tau}_{12}\left({\tau}_{23}\right)\right|\ \mathrm{and}\ {\sigma}_{22}\le 0\\ {}{d}_{2s},\kern0.5em \left|{\upsigma}_{22}\right|\le \left|{\tau}_{12}\left({\uptau}_{23}\right)\right|\kern6em \end{array}\right. $$
(8)
$$ \mathrm{Through}-\mathrm{thickness}\ \mathrm{damage}:{d}_3=\left\{\begin{array}{c}{d}_{3t},\kern0.5em \left|{\sigma}_{33}\right|\ge \left|{\tau}_{13}\left({\tau}_{23}\right)\right|\ \mathrm{and}\ {\sigma}_{33}\ge 0\\ {}{d}_{3c},\kern0.5em \left|{\sigma}_{33}\right|\ge \left|{\tau}_{13}\left({\tau}_{23}\right)\right|\ \mathrm{and}\ {\sigma}_{33}\le 0\\ {}{d}_{3s},\kern0.5em \left|{\upsigma}_{33}\right|\le \left|{\tau}_{13}\left({\uptau}_{23}\right)\right|\kern6em \end{array}\right. $$
(9)
$$ \mathrm{Shear}\ \mathrm{damage}:{d}_s=1-\left(1-{d}_1\right)\left(1-{d}_2\right)\left(1-{d}_3\right) $$
(10)
Table 2 The mathematical definitions for the intralaminar damage propagation model

The derived damage variables, d1, d2, d3 and ds, have the value of 0 when the element is undamaged and 1 when fully damaged and were used to degrade the elasticity matrix to form the damaged elasticity matrix, Cd. Then, the degraded stresses were calculated from σ = Cdεe, in which σ is the resulting stress tensor and εe is the current elastic strain tensor. These degraded stresses and strains were then read as being the new ‘material state’, for a given element, in the next step of the run of the FE model, see Fig. 7. For the simulations of the extent of intralaminar damage as a function of the time-scale of the impact event that are deduced from the model then, following earlier work [37], a value for the damage parameter of equal to, or greater, than 0.9 is used to define the relatively intense intralaminar damage, i.e. to calculate the areas indicated by a red colour in the figure shown later.

5.1.2 The Interlaminar Damage Model

Interlaminar damage involves the delamination, i.e. interlaminar cracking, of layers of the composite laminate. In the interlaminar damage model an elastic constitutive relation was employed to capture the response prior to delamination initiation and a quadratic-stress criterion was employed to govern damage initiation at the composite interface [29, 30]. The evolution of interlaminar damage during the impact event was modelled using a linear-softening material model embedded into a bilinear cohesive (i.e. damage) surface law, where the traction (i.e. stress), t, is plotted versus the displacement, δ; and the area under the t versus δ bilinear law represents the fracture toughness [29, 38]. The initial stiffness of the bilinear cohesive law is defined as k. This cohesive surface law was implemented as a sub-routine in the FE code [29] via an interface element. For such delamination propagation, an interlaminar damage parameter, dinter, was defined to degrade the cohesive stiffness as the interlaminar damage progressed. The details of the mathematical definitions are given in Table 3. In Eq. (27), δ33, δ31 and δ32 are the displacements and t33, t31 and t32 are the cohesive stresses. The equation, k = αE33/tp, as reported in [38], was employed to determine the initial cohesive stiffness, ki(i = 33, 31, 32); where α is a constant much larger than unity, i.e. α ≫ 1, E33 is the transverse Young’s modulus of the composite and tp is the thickness of two adjacent three-ply block layers (i.e. 0.750 mm). The values of the cohesive stiffness along the normal, k33, the first shear, k31, and the second shear, k32, directions were assumed to be equal. In Eq. (28) the term ti(i = 33, 31, 32) represents the current normal or shear stresses and \( {t}_i^0\ \left(i=33,31,32\right) \) represents the normal and shear cohesive strengths. The value of the cohesive strength of \( {t}_{33}^0=43\ \mathrm{MPa} \) was determined based on the theory proposed in [38, 39], which maintains computation accuracy whilst avoiding a very fine mesh and a commensurate increase in computational cost. The shear cohesive strengths, \( {t}_{31}^0={t}_{32}^0 \), may be determined from \( {t}_{31}^0\left({t}_{32}^0\right)={t}_{33}^0\sqrt{\frac{G_{IIc}}{G_{Ic}}} \), as reported in [39], where GIc and GIIc are the interlaminar Mode I and Mode II critical energy release rates, respectively, and the energy-based Benzeggagh-Kenane (B-K) exponent, ηBK, was invoked to allow for mode-mixity effects [39]. In Eq. (29) the term \( \delta =\sqrt{{\delta_{33}}^2+{\delta_{31}}^2+{\delta_{32}}^2} \) is the equivalent displacement at a delamination interface. The displacement values \( {\delta}^0=\sqrt{{\delta_{33}^0}^2+{\delta_{31}^0}^2+{\delta_{32}^0}^2} \) and \( {\delta}^f=\sqrt{{\delta_{33}^f}^2+{\delta_{31}^f}^2+{\delta_{32}^f}^2} \) are the equivalent displacements corresponding to initial failure and final failure, respectively [29, 39]. The calculated interlaminar damage parameter, dinter, see Eq. (29), was used to degrade the cohesive stiffness matrix once the initiation criterion had been reached [29, 39]. The cohesive law properties were used as the ‘material state’ for the cohesive surface model, for a given element, in each step of the FE model, see Fig. 7. For the simulations of the location and extent of interlaminar cracking, i.e. delamination, as a function of the impact velocity and energy, and of the time-scale of the impact event, that are deduced from the model then, following earlier work [37], a value of interlaminar damage parameter, dinter, being equal to, or greater, than 0.9 is used, i.e. to calculate the areas shown by a red colour in the figures shown later.

Table 3 Mathematical definitions for the interlaminar damage model

5.2 Finite-Element Model

The FE model was created in Abaqus 2018 for the two different impactor geometries, as shown schematically in Fig. 8. Each composite ply was modelled as a three-dimensional deformable part and meshed using eight-noded linear-reduced integration (C3D8R) solid-elements with a size of 1 mm × 1 mm. This mesh size was found to give mesh-independent results. In the computational model, the hemispherical steel impactor and the flat-ended steel impactor were modelled as analytical rigid surfaces. In the FE model, 112,076 elements were used for the composite specimen impacted by the hemispherical impactor and 113,024 when impacted by the flat-ended impactor. Cohesive surfaces were defined at the 0o/45o interfaces to capture the interlaminar damage [40, 41]. The general contact algorithm was invoked to govern the global interaction between the impactor and the composite specimen [42, 43]. Friction coefficients of 0.2 and 0.25 were defined in the global contact and cohesive contact, respectively [44, 45]. The computation accuracy was set as ‘double precession’ to reduce the accumulation of errors when running the simulation. Computations were completed using 16 CPUs on a Linux Cluster with a run time of 12 to 15 h, depending on the geometry of the impactors. The off-axis tension tests reported in [27] were employed to characterise the nonlinear shear behaviour of the composite laminates. The input parameters [27, 46, 47] required for the FE model are presented in Table 4. The numerical model runs were stopped when the defined total time for the impact event had expired.

Fig. 8
figure 8

Simulation schematics for investigating the effects of impactor geometry for: (a) the hemispherical steel (RNS) impactor and (b) the flat-ended steel (FFS) impactor

Table 4 Input properties for the FE modelling of the carbon-fibre/epoxy-matrix composite ply [27, 46, 47]

5.3 Results and Discussions

5.3.1 Predicted Loading Response and Damage

The loading responses obtained from the simulation modelling studies for the hemispherical steel impactor and the flat-ended steel impactor are compared with the corresponding experimental results in Fig. 9. As may be clearly seen, the numerical predictions are in very good agreement with the experimental results. Indeed, the model predicts accurately the effects of the shape of the impactor on the effective stiffness of the composite specimens and duration times for the impact event, as well as the maximum loads and displacements that were experimentally measured.

Fig. 9
figure 9

The experimental and FE numerically-predicted results compared for (a) the load as a function of time curves obtained from the hemispherical steel (RNS) impactor case and (b) the flat-ended steel (FFS) impactor case; and the load as a function of displacement curves obtained from (c) the hemispherical steel (RNS) impactor case and (d) the flat-ended steel (FFS) impactor case. (The experimental results from the three replicate tests are shown)

Figure 10 shows the interlaminar damage suffered by the composite specimens from the experimental studies, via the C-scan tests, and the numerical FE simulations. For the C-scan results, the white-dashed line represents a map of the damage that was induced. For the numerically-predicted results, a red-coloured zone indicates the relative extent of delamination induced by the impact event. As may be seen, again a very good correlation is obtained between the experimental and numerical results. For both types of impactor used on the composite, with a lay-up of [+453/03/−453/03]s, the simulation clearly shows the delamination growing along the interfaces between 45o and 0o plies, with the orientation of the lower ply determining the direction of delamination growth, as was observed in the experiments. Considering further the extent of interlaminar damage that was induced by the impact, then Fig. 11 shows cross-sectional images obtained from modelling the composite specimens, with respect to the maximum out-of-plane displacement, with the interlaminar damage induced by the impact event being indicated by the red coloured zone. Figure 11a and b are for the composite specimens which have been struck by the hemispherical steel impactor and the flat-ended steel impactor, respectively. As in Fig. 10, the composite specimen struck by the hemispherical shaped impactor presents a continuous damage area, whilst the specimen struck by the flat-ended steel impactor shows two separate damage areas, as was indeed observed experimentally.

Fig. 10
figure 10

Interlaminar damage maps obtained from the experimental C-scan studies and the FE numerical modelling simulations for the hemispherical steel (RNS) impactor and the flat-ended steel (FFS) impactor. (The 0o fibre direction is indicated. The right-hand side scale for the experimental C-scan results indicates the location of the delamination as a function of the depth through the thickness of the specimen, where the red colour is representative of the front (impacted) surface and the blue colour is representative of the rear (non-impacted) surface of the composite specimen)

Fig. 11
figure 11

Interlaminar damage, at the maximum displacement, obtained from the FE numerical modelling simulations, with the central region magnified. For: (a) the hemispherical steel (RNS) impactor and (b) the flat-ended steel (FFS) impactor

Further, to study the effects of the impactor shape on the intralaminar damage suffered by the composite specimens, the extent of intralaminar matrix damage was also extracted from the simulation results and is shown in Fig. 12. (Shown as cross-sectional images at the point of maximum displacement of the composite specimen, with the intralaminar damage induced by the impact event being indicated by the red coloured zone.) Figure 12a and b are for the composite specimens struck by the hemispherical steel impactor and by the flat-ended steel impactor, respectively. Figure 12a shows that the hemispherical steel impactor causes localised damage throughout the thickness of the composite, directly below the impactor. This through-thickness intralaminar damage steadily increases in extent towards the rear face of the specimen. In contrast, for the flat-ended steel impactor the through-thickness intralaminar damage is very much concentrated adjacent to the rear face of the specimen and extends below the impact site, as may be seen in Fig. 12b, unlike the delamination damage which forms a ring around the impact site, see Fig. 11b. (It should be noted that the C-scan test only detects delamination damage and this intralaminar damage was not visible from the C-scans, or with the naked eye.)

Fig. 12
figure 12

Intralaminar matrix damage, at the maximum displacement, obtained from the FE numerical modelling simulations, with the central region magnified. For: (a) the hemispherical steel (RNS) impactor and (b) the flat-ended steel (FFS) impactor

Finally, some of the key results obtained from the experimental studies and the modelling predictions, i.e. maximum load, maximum displacement, impactor contact duration time and damage area, are summarised in Table 5. The results clearly confirm the ability of the model to capture accurately the impact response of, and the damage area in, the composite specimens.

Table 5 Comparison of some key results obtained from the experimental studies and the FE numerical simulations for the impact of the composite specimens struck with the two different impactors

5.3.2 Predicted Deformation of the Composite Specimen

The deformation of the composite specimens, at the maximum displacement of the specimens, may be obtained from the modelling simulation of the hemispherical steel impactor and the flat-ended steel impactor. The results are presented in Fig. 13 for the distribution of the through-thickness strain, ε33, in the composite specimens, at the maximum displacement, obtained from the simulations, with the central region magnified, for: (a) the hemispherical steel impactor and (b) the flat-ended steel impactor. For both the hemispherical steel impactor and the flat-ended steel impactor there is a good correlation between (a) the strain, ε33, distribution maps and (b) the extent and location of the intralaminar matrix damage in the composite specimens, e.g. compare Figs. 12 and 13. For the hemispherical steel impactor, the deformation, which was initially localised and located in the central region below the impactor, propagated to the peripheral regions during the interaction between the impactor and the specimen. For the flat-ended steel impactor, the damage initiated in the composite specimen in the region contacted by the circular edge of the impactor. Also, it is notable that less curvature was observed in the central region below the flat-ended steel impactor, compared with when the hemispherical impactor was used. Interestingly, these observations also explain why very little delamination damage occurred in the central region of the composite specimen when it was impacted with the flat-ended steel impactor, see Figs. 10 and 11b, since delamination damage tends to initiate in regions of high curvature or high-strain gradient.

Fig. 13
figure 13

The distribution of through-thickness strain, ε33, at the maximum displacement, obtained from the FE numerical modelling simulations, with the central region magnified. For: (a) the hemispherical steel (RNS) impactor and (b) the flat-ended steel (FFS) impactor. (The level of strain is indicated in the scale bar)

From Figs. 10 to 13 it is clear that appreciable strain and damage have occurred in the composite specimens at the maximum displacement. These effects arise, of course, in bringing either impactor to rest. At this point the initial kinetic energy of the impactor has been expended in intralaminar and interlaminar damage in the composite specimen, along with some limited slip and friction at the boundary supports. Also, there is elastic strain-energy stored in the specimen. Indeed, this stored elastic strain-energy is mostly responsible for rebounding the impactor. It is interesting to note that the flat-ended steel impactor is arrested more rapidly, and the maximum displacement is reduced, compared with using the hemispherical impactor. Therefore, the effective stiffness of the composite specimen is lower when using the hemispherical steel impactor and, hence, a greater overall displacement of the specimen is observed. However, the higher localised strains generated by the hemispherical impactor produced somewhat more damage in the composite specimen, again as was accurately predicted by the FE model, see Table 5.

6 Conclusions

A hemispherical steel impactor and a flat-ended steel impactor have been used to impact carbon fibre-reinforced/epoxy-matrix composite laminates with an energy level of 15 J. All the composite specimens were then inspected using an ultrasonic C-scan test to assess the impact damage. It was found that the composite specimens struck with a flat-ended steel impactor suffered a higher maximum load upon impact, but experienced a smaller out-of-plane displacement, than those struck with the hemispherical steel impactor. The comparison of the experimental C-scan images showed that the delamination maps of the composite specimens impacted with the hemispherical steel and the flat-ended steel impactors were similar in area but with somewhat more areal damage occurring when the hemispherical steel impactor was used. However, the composite specimens impacted with the hemispherical steel impactor showed a continuous delamination area, which was centred around the point of impact, whilst the flat-ended steel impacted specimens exhibited two separate delamination areas, with a central zone exhibiting no delamination. To further understand the effects of the geometry of the impactor head, a three-dimensional, finite-element (FE) computational model was developed to simulate the impact response of the composite laminates when struck with the two different types of impactor. It was found that the FE numerical simulations, with an embedded elastic-plastic damage model for the material response, very accurately predicted the key parameters of the impact event, such as (a) the shape and area of the delamination, (b) the maximum load and displacement, and (c) the duration time of the impact event. Furthermore, these numerical simulations have explained the underlying reasons for the effects of the shape of the impactor head on the impact behaviour of the composite specimens. These studies have, therefore, increased our understanding of the behaviour of composite laminates when subjected to a relatively low-velocity impact and have proposed, and validated, an elastic-plastic three-dimensional FE numerical model for predicting the response of the laminate. This model will play a significant role in future academic research and industrial applications with respect to predicting and improving the impact behaviour of composite laminates.