Damage development in stitch bonded GFRP composite plates under low velocity impact: Experimental and numerical results

Abstract

Damage development due to impact needs to be understood to evaluate the consequences of impact on composite structures. This study concentrates on modelling and measuring damage development due to low velocity impact on thick industrial composites made from glass fibre epoxy by vacuum-assisted resin infusion. Cross-plied laminates were tested with different impact energy and different number of interfaces (clustering). Results were compared to a 3D finite element analysis. Interfaces and their damage development were modelled with cohesive elements. Intra ply properties were modelled by progressive failure analysis. Many elements and large memory use were needed to obtain sufficient modelling accuracy. However, all input parameters of the model were based on widely available and independently obtained material properties. Impact force and time to initiate damage and maximum force were measured and related to impact energy and clustering. Damage development was monitored optically in the translucent material for all test cases. The results show that the numerical model using only simple and independently measured material data was able to predict the impact behaviour for the different energies and different stacking sequences.

Keywords

Composite Glass Fibre Reinforced Polymer(GFRP)low velocity impact experimental tests finite element method delamination cohesive elements puck failure criterion

Introduction

Many composite applications may experience impact and the consequence of such impact needs to be evaluated. This study concentrates on low velocity impact caused by, for example, dropped tools, dropped components and minor collisions. Pipes, pressure vessels, ships, wind turbine blades, etc. may suffer from impact damage. For this reason industry needs readily available tools to design for impact resistance or to access the severeness of impact damage.

Impact is a complex interaction of an object with a structural component. It typically causes matrix cracking, fibre failure and delamination.¹ Impact was initially mainly experimentally investigated for thin carbon laminates produced by autoclave and used in aerospace applications.¹ Only few studies have been published about impact on thicker laminates based on glass fibres and produced by vacuum infusion or hand layup, as used in the marine industry and for wind turbine blades. Sutherland and Guedes Soares^2–5 performed an extensive test program to evaluate impact events on marine composites. E-glass woven-roving polyester samples were impacted with an instrumented impact machine with various impactor shapes and specimen geometries.³ Their studies of impact loads on thick samples showed that shear effects, delamination and indentation were important parameters to describe impact.⁵ The effects of composite constituents and laminate thickness were also investigated.⁵

Many numerical models have been developed to characterize impact events on composites, ranging from simple 2D analytical models^6–12 to complex 3D finite element (FE) simulations.^13–18 Two-dimensional analytical models^6–12 are usually applicable to limited impact scenarios and specimen geometries, because generally they cannot characterize the 3D nature of damage properly. The complexity of the impact event due to the highly dynamic behaviour, multiple loads, contact, damage and variable boundary conditions can lead to oversimplifications by the analytical approaches.

Currently, the finite element method (FEM) is most suitable for simulating impact events on composites due to new advances in numerical models and the increase of the available computational power. The FEM approach gives the possibility to simulate different material behaviour simultaneously and to reproduce the complexity of the problem. Several FE modelling methodologies have been presented in recent years to analyse the structural behaviour of composite materials under low velocity impacts.

Different strategies have been proposed for modelling interlaminar damage of composites based on either the continuum damage mechanics approach^19,20 or the fracture mechanics approach.^21,22 More recently, the cohesive approach, based on a combination of strength criteria and fracture mechanics, has gained much interest.^23–29 The quasi-brittle process of intralaminar damage has been mainly modelled with continuum damage mechanics. In these methods behaviour after damage is modelled by a material degradation model.

González et al.¹³ performed a broad study, comparing numerical results with experimental tests. An advanced interface model, based on the work by Turon et al.,³⁰ was used to model delamination. The intralaminar effects were evaluated using the LaRC04 failure theory developed by Maimi et al.^31–33 An advanced and innovative approach was presented by Bouvet et al.,¹⁴ where discrete interface elements were used to model both the interlaminar and intralaminar damages. Both methods gave good agreement between theory and experiments, but they need an extensive set of input parameters requiring large and complicated test programs.

In safety critical applications impact modelling is generally not accepted for evaluating damage today. For example, the DNV offshore standards^34,35 for composite components and composite risers require testing of components with impact damage. Generally available models with clearly defined modelling parameters, together with evidence that experimental results can be reproduced, will be needed before the expensive experimental test requirements in the design codes can be replaced with a modelling approach.

Impact evaluation requires three steps: defining the impact scenario, characterizing the damage in the material caused by the impact and finally evaluating the consequence of this damage on the structural performance of the component. This paper addresses the first two aspects and shows a FE analysis method that can characterize the impact damage well based on readily available material parameters. This characterization can be used for further analysis of the damaged structure.

The impact model is based on 3D FE modelling with cohesive elements for the interlaminar damage and a strength-based failure criterion (based on Puck^36–38 and Hashin³⁹) for the intralaminar damage. This approach is fairly complicated to model and requires long computation times, but it is available to most engineers. Modelling results are compared to experimental tests on simple flat, cross-plied laminates with different stacking sequences (to evaluate the influence of ply clustering). This is seen as a step to build up confidence in the modelling methods, eventually allowing engineers to design impact resistance of low cost composites without much laboratory testing.

Experimental setup

Materials

All tests were done on an R-glass epoxy composite produced by vacuum-assisted resin infusion, as commonly used in wind turbines, ships and offshore applications. The reinforcement was a Saertex unidirectional (UD) stitch bonded layer (a small amount of fibres, less than 0.1% of the total weight, were used in the transverse direction to keep the fibre tows together) made from 3B’s HiPer-tex W2020 R-glass fibre with a weight of 1150 g/m². Momentive EPIKOTE MGS 135 epoxy resin with EPIKURE MGS 137 curing agent was used as resin. A mix ratio of 100:30 was selected with a curing time of 24 h at room temperature and post curing at 80℃ for 15 h.

UD ply properties were measured for this material and results are summarized in Table 1.⁴⁰ Table 1 also shows through thickness ply properties. The values were chosen theoretically based on symmetry considerations for orthotropic plies (same properties in direction 2 and 3), since testing is difficult especially for coarse thick laminates used here. This set of data is the main input for the FE analysis used in this study. All test results were obtained by standard test methods and are usually also readily available for other materials.

Table 1.

Summary of the tested material properties for the used HiPer-Tex R-glass fibre/Momentive EPIKOTE MGS 135 resin and EPIKURE MGS 137 composite produced by vacuum infusion.

Properties	Value
Ply properties
Ply thickness	0.86 mm
Fibre volume fraction	54%
Density	1230 kg/m³
Ply elastic properties	E₁ = 44.87 GPa; E₂ = E₃^a = 12.13 GPa; G₁₂ = G₁₃^a = G₂₃^a = 3.38 GPa; ν₁₂ = ν₁₃^a = 0.30; ν₂₃^a = 0.5
Ply strength (MPa)^b	X_t = 1006.30; X_c = 487.00; Y_t = Z_t^a = 45.95; Y_c = Z_c^a = 131.90; S₁₂ = S₁₃^a = S₂₃^a = 49.51
Interface properties
Elastic properties	K_nn^c = 12.13 GPa; K_ss^c = K_tt^c = 3.38 GPa
Strength (MPa)	t_n^c = 45.95; t_s^c = t_t^c = 49.51
Fracture toughness (N/mm)	G_Ic = 0.98; G_IIc = G_IIIc^c = 3.71
Mode interaction – BK	η = 1.40

Assumed based on symmetry considerations.

X is in the fibre direction, Y is normal to the fibre direction and Z is in the through thickness direction; t is for tension and c for compression; S₁₂ is for shear.

Selected input for cohesive elements, see ‘Numerical simulations – Interlaminar damage model’ section.

Samples for impact testing

Three different cross plied laminates were tested. All laminates had eight plies and a thickness of 6.8 ± 0.2 mm. The laminates are called here L1 [0,0,90,90]_s, L2 [0,90,90,0]_s and L3 [0,90,0,90]_s. The three laminates have the same in-plane laminate properties due to the same amount of plies in 0 and 90 directions. The important difference between the laminates is the number of interfaces (n_d) where delamination can propagate (L1 n_d = 2, L2 n_d = 4 and L3 n_d = 6) due to the different plies clustering. This can be defined as the tendency to stack together plies with the same fibre orientation. For this reason, delamination shape and dimensions are expected to vary for each laminate (for the same impact energy).

Drop weight test

Impact tests were performed according to the standard test method ASTM D713⁴¹ in an Instron CEAST 9350 drop tower machine with a CEAST DAS 64K data acquisition system. A cylindrical impactor of 5.02 kg with 20 mm diameter spherical end was used.

The test specimens had the dimensions 150 ± 0.1 × 100 ± 0.1 mm. They were supported on all sides in such a way that only a central region of 125 × 75 mm was free to deform under the impact. A special clamping system was adopted to prevent out-of-plane movements of the sample during the impact event. Four toggle clamps with a minimum capacity of 1100 N were installed on the longer specimen side, as suggested by the standard.⁴¹ Details about the test configuration are shown in Figure 1. The impact force was measured by an accelerometer embedded in the impactor head via a high-speed data acquisition system triggered on the impact. An acquisition frequency of 1000 kHz was used for all the tests.

Impact energy

The samples were tested with different impact energies by changing the drop height of the impactor and with it the velocity, when the impactor hits the sample. The test standard ASTM D7136⁴¹ suggests an impact energy of 45.5 J for the given specimen thickness. All laminates were tested with this energy. The impact speed was 4.25 m/s in that case. Two other energies were applied to laminate L2: 22.70 J (velocity of 3.00 m/s) and 68.25 J (velocity of 5.2 m/s), which represent +50% and −50% of the reference energy. A summary of the laminates and the corresponding impact energies is given in Table 2 together with the experimental results.

Table 2.

Impact data from experimental tests.

Laminate		Number of interfaces	Impact energies [J]		F_d [N]		t_d [ms]		F_Max [N]		δ_Max [mm]
L2	[0,90,90,0]_s	4	22.70	(–50%)	7222	(–10%)	0.69	(+15%)	12,749	(–15%)	3.66	(–31%)
			45.50	(Ref)	8000	(Ref)	0.60	(Ref)	15,010	(Ref)	5.37	(Ref)
			68.25	(+50%)	10008	(+25%)	0.59	(–2%)	19,899	(+33%)	6.69	(+25%)

L1	[0,0,90,90]_s	2	45.50		8000		0.60		13850		5.77
L2	[0,90,90,0]_s	4							15010		5.37
L3	[0,90,0,90]_s	6							14667		5.51

Numerical simulations

Intralaminar damage model

Matrix cracking was modelled by the Puck criterion^36–38 and fibre failure by the Hashin criteria.³⁹ Both methods were applied on the ply level using a user defined material model coded in a VUMAT subroutine, written in FORTRAN, for Abaqus/Explicit.⁴² Ply damage is evaluated in each integration point at every step of the simulation. The Puck criterion is well suited for modelling impact due to its capability to consider also the out-of-plane components of the stresses in the failure evaluation. Following a short explanation of the Hashin and Puck criteria is reported.

For the Hashin criterion,³⁹ the fibre failure occurs when the parameter $f_{ff} \geq 1$ . $f_{ff}$ is evaluated by equation (1) for tension and compression loads

f_{ff} = {(\frac{σ_{1}}{X_{t}}) 2 + \frac{1}{S_{12}^{2}} (τ_{12}^{2} + τ_{13}^{2}) for σ_{1} \geq 0 (\frac{σ_{1}}{X_{c}}) 2 for σ_{1} < 0

(1)

where

σ_{1}

is the stress in the direction of fibre, and

τ_{12}

and

τ_{13}

are the in-plane and out-of-plane shear stresses. X_t and X_c are the tension and compression strength in the fibre direction. S₁₂ is the in-plane shear strength. More details about this failure theory can be found in Ref. 39.

The Puck failure criterion³⁸ is an interactive stress-based criterion applicable for UD composite lamina. The matrix failure criterion is based on the assumption that failure is created only by the stresses that act on the fracture plane (σ_n, τ_nl and τ_nt) inclined at θ_fp to the material plane. The normal and shear stresses acting on this plane are calculated by rotating the 3D stress tensor from the material coordinate system to the fracture plane using classical tensor transformations. The matrix failure is now described by an interfibre failure criterion (see equation (2)) that is only a function of the stresses acting on the fracture plane. Failure is reached when $f_{iff} \geq 1$

where $R_{⊥}$ is the failure resistance normal to the fibres, $R_{⊥ ψ}$ and $R_{t ∥}$ are the shear resistances and $P_{⊥ ψ}^{+}$ and $P_{⊥ ψ}^{-}$ are the slope parameters representing the internal friction effects in the Mohr–Coulomb failure criterion (more details in Refs. 38,43). The parameters presented before can be calculated using the following expressions

f_{iff} = {\begin{matrix} \sqrt{[(\frac{1}{R_{⊥}} - \frac{P_{⊥ ψ}^{+}}{R_{⊥ ψ}}) σ_{n} (θ)] 2 + (\frac{τ_{nt} (θ)}{R_{⊥ ⊥}}) 2 + (\frac{τ_{nl} (θ)}{R_{⊥ ∥}}) 2} + \frac{P_{⊥ ψ}^{+}}{R_{⊥ ψ}} σ_{n} (θ) for σ_{n} \geq 0 \\ \sqrt{(\frac{τ_{nt} (θ)}{R_{⊥ ⊥}}) 2 + (\frac{τ_{nl} (θ)}{R_{⊥ ∥}}) 2 + [(\frac{P_{⊥ ψ}^{-}}{R_{⊥ ψ}}) σ_{n} (θ)] 2} + \frac{P_{⊥ ψ}^{-}}{R_{⊥ ψ}} σ_{n} (θ) for σ_{n} \leq 0 \end{matrix}

(2)

\frac{P_{⊥ ψ}^{+}}{R_{⊥ ψ}} = \frac{P_{⊥ ⊥}^{+}}{R_{⊥ ⊥}} (\cos ψ) 2 + \frac{P_{⊥ ∥}^{+}}{R_{⊥ ∥}} (\sin ψ) 2 \frac{P_{⊥ ψ}^{-}}{R_{⊥ ψ}} = \frac{P_{⊥ ⊥}^{-}}{R_{⊥ ⊥}} (\cos ψ) 2 + \frac{P_{⊥ ∥}^{-}}{R_{⊥ ∥}} (\sin ψ) 2 (\cos ψ) 2 = \frac{τ_{nt}^{2} (θ)}{τ_{nt}^{2} (θ) + τ_{nl}^{2} (θ)} (\sin ψ) 2 = \frac{τ_{nl}^{2} (θ)}{τ_{nt}^{2} (θ) + τ_{nl}^{2} (θ)} R_{⊥} = Y_{t} R_{⊥ ∥} = S_{21} R_{⊥ ⊥} = \frac{Y_{c}}{2 (1 + P_{⊥ ⊥}^{-})}

(3)

where all properties can be calculated from Y_t,Y_c and S₂₁ as given in Table 1. In addition, the special Puck parameters (

P_{⊥ ∥}^{+} = 0.30, P_{⊥ ∥}^{-} = 0.25, P_{⊥ ⊥}^{+} = 0.225

and

P_{⊥ ⊥}^{-} = 0, 225

) are needed. The measurement of these parameters requires complicated multi-axial testing. However, Puck recommends use of the default parameters reported in Ref. 38. The default parameters were used in this study. A sensitivity study of the effect of the choice of the Puck parameters on impact modelling results was not carried out, because generally these parameters seem to give good results for modelling ply damage in many different cases. If very accurate damage geometries are needed, a more detailed analysis of the choice of the Puck parameters should be done.

The Hashin criterion needs the fibre dominated ply strengths X_t and X_c as well as the shear strength S₁₂ as input parameters. The puck criterion needs the matrix dominated strengths Y_t and Y_c together with the shear strength S₁₂ as input. All these strengths were measured and are reported in Table 1.

Modelling of the propagation of damage was done by progressive failure analysis in this study. When failure is predicted in a particular position of the composite, the values of the stiffness matrix of the FEs are reduced locally according to the failure type predicted by the failure criterion. In the case of fibre failure, all of the components of the stiffness matrix are reduced to zero. In the case of matrix cracking, only the transverse and shear components are reduced. In both cases, the out-of-plane component of the stiffness matrix is kept equal to the undamaged one to avoid any unnatural penetration. This degradation modelling is a simple and reasonable approach that does not need any further testing to establish properties after first damage.

Interlaminar damage model

Delamination onset and propagation was modelled by the use of a cohesive zone model (CZM). The model is based on the assumption that there is a thin layer of a separate material with an independent constitutive law between the composite layers. This interface connects the different composite layers, and it is simulated by thin cohesive elements placed between each layer where the delamination is expected to propagate. Delamination can propagate only between differently oriented layers,¹ and thus, the delamination path is known a priori.

The classical energy based bilinear traction-separation law was used in the cohesive zone elements. The cohesive behaviour was defined directly in terms of a traction-separation law based on the following assumptions:

Linear elastic traction-separation law for the undamaged material;

Damage initiation predicted using the quadratic nominal stress criterion;

Linear degradation law function of the dissipated energy;

The Benzeggagh–Kenane⁴⁴ (BK) law for the mixed opening mode;

The interfacial material properties are matrix dominated. For this reason the following assumptions were made:

Normal modulus K_nn = E₂;

Normal strength t_n = Y_t;

Equal shear moduli in both directions K_ss = K_tt = G₁₂;

Equal shear strengths in both directions t_s = t_t = S₁₂.

The parameters K_nn, K_ss, K_tt, t_n, t_s and t_t are defining strength, stiffness and fracture energy of the CZM model, as described in detail in the ABAQUS manual.⁴⁵ In the approach chosen here, these parameters are defined by the readily available ply properties from Table 1.

Numerical implementation

The impact tests reported earlier were numerically modelled using the FE software Abaqus Explicit 6.11.⁴⁵ A quarter of the model (just for graphical reason) is shown in Figure 2 and it reflects the experimental setup shown in Figure 1.

Experimental test results and simulations

In order to be able to model damage development in each individual ply of the laminate, each layer (both composite layers and cohesive interfaces) needed to be modelled separately. Two C3D8R elements were used through the thickness in each ply. This element is the standard continuum solid hexahedral element with eight nodes and the reduced integration scheme.⁴⁵ The intralaminar damage model, implemented by means of a VUMAT subroutine, was integrated in the elements and evaluated for each integration point at each time increment. The hourglass, by means of the enhanced method, and the distortion controls were used for all of the elements.⁴⁶

The interfaces between differently oriented plies were discretized using the cohesive element COH3D8⁴⁵ by means of the classical bilinear traction-separation law. Each of the interfaces was modelled with a finite thickness of 0.1 mm, which was deducted from the thickness of the adjacent layers to obtain the same total laminate thickness. The total cohesive thickness, compared to the global thickness, represented 2.9% for laminate L1, 5.8% for laminate L2 and 8.7% for laminate L3. The cohesive thickness of 0.1 mm was chosen after several trials to obtain the best compromise between the computational time and the reduction of the global laminate stiffness related to the presence of the cohesive interface.

A mesh sensitivity analysis was conducted to evaluate the best compromise between the computational time and the results accuracy. The results clearly showed that the size of the cohesive elements is of crucial importance for accuracy. It was confirmed that, according to Turon et al.,⁴⁷ a sufficient number of elements needs to be used in the cohesive zone length to obtain the proper evaluation of delamination onset and propagation. For this reason, all three laminates were discretized with a very fine mesh: cohesive elements of 0.5 × 0.5 mm² (one element in the thickness direction) and solid elements of 1 × 1 mm² (two elements in the thickness direction for each ply). A constant in-plane ratio of one to four was used between the solid and cohesive elements. An average of 640 K elements was obtained for the different models. The FE simulations were computationally expensive; approximately 25 GB of memory was needed with a total computational time of about 48 h on two cluster nodes with a total of 32 cores at 2.60 GHz. The long computational time is directly related to the very small stable increment time required for this type of simulation. Moreover, the simulations were conducted without the use of mass scaling in any of the different parts of the model to prevent an effect of mass scaling on the accuracy of the results. This decision is in contrast to what other authors have done in the past,^13,48 but due to the lower number of elements and the available computational power, there was no reason to use the mass scaling with the consequent possibility of affecting the results. It is also expected that computer power will keep increasing with time, making the required memory more widely available and allow shorter calculation times in the near future.

Tie constrains were used to connect each composite layer with the adjacent cohesive interface. This method is quite complex but it gives the possibility to obtain a better control of the element sizes between the cohesive and composite layers. A finer mesh was used in the cohesive layers without leading to the same increase in the remaining elements of the laminate, optimizing the total required computational time.

The plate support was modelled as completely rigid and placed directly in contact with the samples. The four rubber pins holding the sample in place were modelled as a deformable body with a spring load of 250 N on each pin. They were placed directly in contact with the upper part of the sample. The impactor was modelled as a rigid shell body with a concentrated mass of 5.02 kg. An initial vertical velocity, according to the test setup, was prescribed to the impactor. Moreover, a force of 49 N in the vertical direction was applied to the impactor centre of mass to simulate the gravitational force (5.02 kg × 9.8 m/s² = 49.2 N).

The boundary conditions are schematically shown in Figure 2; they should represent the experimental conditions reasonably well.¹³ In addition to the physical conditions, symmetry constraints were applied on the sample centrelines to obtain a more consistent damage development.

The contacts between the impactor, the specimen, the support and the pins were simulated by the use of the general contact algorithm with a penalty enforcement method.⁴⁹ Moreover, due to the possible deletion of the cohesive element by reaching full damage, the same contact formulation was also applied between all the layers to prevent any unnatural penetration. The Coulomb model was selected for the friction formulation. The friction coefficient is strictly related to the material in contact and the surface quality. Several studies^50,51 were conducted to evaluate the friction between the delaminated surfaces showing variation of friction due to the interface angle. Reported friction coefficients are in the range between 0.2 for a 0°/0° interface and 0.8 for a 90°/90° interface. The friction coefficient between the metal part and composite is stated as 0.4.⁴⁸ In the current work, an average friction coefficient, μ = 0.3, was used for all of the implemented contacts with no variation for the different interfaces/parts.

Due to the progressive damage development during the simulation a localized reduction of element stiffness may happen. This can result in excessive distortions for both the cohesive and brick elements, eventually causing numerical instability. To avoid these problems, the fully damaged elements were removed from the calculation.

Impact loads

The impact test program described in Table 2 was carried out with six parallel tests for each sample/energy combination. The experimental results, force versus time or force versus displacement curves, obtained from the instrumented impactor are reported in Figures 4 and 5.

Figure 3 shows test results for laminate L2 impacted with 45.5 J as a typical example. The impact curves (load versus time and load versus displacement) are shown for all six parallel samples together with the averaged curve. The measurements show large vibrations that are typical for impact tests, but it can also be seen that the six parallel tests had very reproducible results. This was the case for all measurements; future plots (Figures 4 and 5) will only show average results.

When damage initiates, the kinetic energy starts to dissipate as irreversible material damage (matrix cracking, fibre failure and delamination). The threshold force F_d and time t_d to initiate damage (mostly delamination) are also shown in Figure 3. Due to the dynamic vibrations it is not always easy to identify the damage threshold values. The damage threshold is defined as the point where the measured structural stiffness is reduced,¹³ i.e. a reduction of the slope of the force–time and force–displacement curves can be seen (Figure 3). The impactor velocity gradually decreases to zero when the deformation is at its maximum δ_Max. At this point the impact force reaches its maximum F_Max.

The influence of the impact energy was tested for three impact energies on laminate L2 as shown in Figure 4. The influence of clustering was measured for laminates L1, L2 and L3, which were impacted with a constant energy of 45.50 J. The results are shown in Figure 5.

Figures 4 and 5 also show the numerical simulations. Oscillations are shown, especially in the impactor force curves, which are related to the numerical instability during the damage propagation. To better compare the data, the numerical curves were filtered by removing the higher frequency components of the signal. Original and filtered curved are reported for all of the impact configurations.

Generally good agreement between experimental data and simulations can be seen for the initial part of the curves and damage initiation (F_d). The maximum impact force and total impact time are slightly mispredicted by the simulations, but it should be noted that these values are highly affected by the complex damage development predicted by the progressive failure analysis. The error of maximal 10% should be acceptable for most practical applications.

Damage development was recorded for all specimens. Typical examples are shown in Figure 9 together with numerical simulations. This will be discussed in more detail in the following sections.

Impact force and damage initiation

Focusing on the force–displacement curves in Figures 4 and 5 it is clearly visible that the initial elastic response of the samples (gradient of the initial part of the curves) was not affected by the impact energy nor by the related change in strain rate. It might be that the change in impact velocity from 3.0 to 5.2 m/s was too small to create any significant changes in material response. It confirms that using of the rate independent materials’ properties in the numerical simulations was a reasonable choice for low velocity impact modelling.

The threshold force for the creation of damage F_d increased with increasing impact energy, while the damage initiation time t_d decreased slightly as given in Table 2.

The threshold value for the damage initiation should be a material/thickness parameter, not affected by the impact energy, as reported by Davies et al. in Refs. 52 –54. They used a simplified approach to evaluate the critical load for the delamination onset based on the Young’s modulus, mode II strain energy release rate, Poisson’s ratio and laminate thickness. The same approach was used by Sutherland and Guedes Soares⁴ for thin and thick glass fibre composite obtaining a good fit with experimental data. This study showed different results. Moreover, the initiation force was also affected by the impact energy, as reported for carbon fibre laminates in Ref. 13. Future studies are needed to better understand this behaviour.

Table 2 presents that clustering had nearly no effect on damage initiation. F_d and t_d were located at almost the same position for all the three laminates. Even if the dimensions and shapes of the delaminations differ due to the different number of interfaces, the damage starts at the same value of force. The FE simulations show the same results for F_d and t_d.

The maximum impact force F_max increased with the impact energy. Typically F_max increases proportional to the impact energy for thin carbon fibre composites.^1,9,13 For the thick glass fibre laminates tested here, the increase was not linear, opposite to what was reported by Sutherland and Guedes Soares.³ The impactor’s maximum displacement δ_Max also increased with the impact energy displaying a fairly linear relationship between change of impact energy and change of δ_Max. Clustering had little effect on the maximum force, as shown in the experimental results and numerical simulations.

Use of the numerical simulation enables calculating the total impact energy and the energy used for damage development. The energies are shown in Figure 6 together with the experimentally measured absorbed impact energies. The agreement between experiments and simulations is good, especially in the initial part of the curves. It can be seen that for laminates L2 exposed to increasing impact energies, damage initiation always happens after about 0.5 ms and about 20% of the total impact energy (Figure 6, top). This corresponds fairly well to the damage initiation times, t_d, given in Table 2.

The numerical results reported in Figure 6 (top) clearly show the increase of the total dissipated energy with the impact energy for both the intralaminar and interlaminar damage. The energy evaluations are even more important to investigate the effect of plies clustering, Figure 6 (bottom). Layup 1, with the highest plies clustering, showed the lowest amount of energy dissipated by damage. Moreover, the damage initiation is predicted after 0.8 ms at about 40% of the impact energy for L1 and the previously found 0.5 ms at about 20% of the impact energy for laminates L2 and L3.

Damage development

Once the damage threshold energy is exceeded, the first damage is matrix cracking. A rough evaluation of matrix cracking was performed by optical observations on the impacted samples. These clearly show a proportional increase of the matrix cracking density with the increase of the impact energy. Due to the thick laminates used in this work, the typical pine tree damage pattern¹ was observed. The matrix cracks started at the layer directly in contact with the impactor due to the high contact stress. The matrix crack density increased from the top ply (closest to the impactor) to the lower plies.

The numerically predicted matrix cracking envelopes for each ply are reported in Figure 8 for the three laminates investigated. The predicted matrix cracking envelopes followed the expected behaviour, developing in the direction parallel to the fibre. The upper layer also presented a high density of matrix cracks, which is related to the higher contact stresses that act locally between the sample and the impactor.

Delaminations were clearly visible in the translucent glass epoxy laminates. Their dimensions and shapes were determined from pictures taken using the backlight technique.^55,56 Figure 9 shows photos and numerical predicted delamination areas for laminates L1, L2 and L3. The backlight technique measures the overall area of delamination, but it cannot clearly measure the delamination size of each interface. In order to precisely evaluate size of delamination at each interface, ultrasound C-Scan and thermography tests, commonly used for aeronautical high quality carbon fibre composite, were performed on the impacted plates. Due to the nature of the composite (thick layer of UD R-glass fibre), both techniques were not able to give any better results. For this reason only the images obtained by the backlight technique were reported here.

The overall measured delamination areas (Figures 7 and 9) increased with both the impact energy and ply clustering. By increasing the impact energy by 50% from 45.5 to 68.25 J, the overall damage increased more than twice compared to the variation between 22.70 and 45.50 J (see Figure 7). For ply clustering the overall delamination area increased proportional to the number of interfaces (Figure 7).

In the FE simulation, delamination is evaluated at each interface by the damage sustained in the cohesive elements. Once an element shows full damage, the element is deleted and its original position is shown as delamination. In Figure 9, the numerical and experimental delamination areas are shown for each configuration. The classical elongated peanut shapes¹ were well represented in the simulations. The main direction of the delamination followed the fibre orientation of the lower layer, as shown in Figures 8 and 9. The shapes of the overall delamination, in both sample directions, were comparable with the experimental tests for all of the tested configurations, as shown in Figure 9. It is clearly shown that the numerical model is able to predict the general trend of the overall delamination area for the different configurations.

The samples showed a small amount of fibre failure for the impact energy of 45.50 and 68.25 J. Fibre failure was limited to the outer layer in the proximity of the contact area with the impactor. The sudden drop in the force versus time curves from Figure 4 at about the maximum impact force for the energy of 45.50 and 68.25 J is also an indication of fibre damage. The drop would be due to the sudden loss of the local material stiffness caused by the fracture of the fibres. The numerical simulations did not show fibre failure for any of the tested configurations in contrast to the experimental tests. This discrepancy can be related to the contact formulation used in the FE model, which did not allow the transition of such high local stresses.

Quantitatively the simulations overestimate the size of the damage zone by approximately 10–15% for all cases. This may be due to variations in material strength, especially the matrix dominated properties. This can be related to the numerical approach used here. The matrix cracking caused the local stiffness to drop to zero in the progressive failure analysis. A continuum damage model^31–33 based on the dissipated energy would have reduced the stiffness more gradually. This would have created slightly less matrix cracking and possibly less overall damage. However, these models^31–33 require more material parameters than the standard properties used here. These extra material measurements are often difficult to obtain, especially for low cost applications. For this reason a simpler approach was applied here, using an advanced failure criterion but with no continuum damage modelling. Even if the matrix cracking distribution is slightly overestimated, the simulated impact curves and the delaminations were well predicted.

Conclusions

Experimental impact tests were conducted on flat cross-plied composite laminates produced with R-glass fibres and epoxy by vacuum infusion. The experimental impact curves were compared with numerical simulations to evaluate the accuracy of the implemented FE model. The influence of the impact energy and clustering effect on the impact-induced damage was also investigated.

Experimental test results

Several tests using the same stacking sequence but with different impact energy were conducted. The shape of the load displacement curves was not affected by the impact energy, but the maximum impact forces, F_max, showed a non-linear increase. The damage onset and propagation could be identified clearly in the force–displacement curves by a variation in the curve's slopes. The total damage (delamination, matrix cracking and fibre failure) increased with the impact energy as expected. Moreover, the increase in impact energy caused a more than proportional increase of the overall delamination area.

Three different layups were impacted at the same energy to investigate the clustering effect, showing almost no variation in the impact load versus displacement curves. The overall projected delamination area and shape was significantly influenced by the number of interfaces available as consequence of the difference in available interface for the delamination onset and propagation. The overall delamination area was almost linearly dependent on the number of available interfaces. No influences on the threshold damage initiation points (F_d, t_d) were observed for the different laminates.

Numerical simulations

The experimental impact tests were numerically modelled using an advanced FE model based on the combination of cohesive elements for the interlaminar damage and strength-based failure criteria for the intralaminar damages. All material input parameters were based on readily available and independently measured data.

The numerical results showed good accuracy predicting: the total impact time, the maximum impact force, the impactor displacement, the energy history, the delamination shape and position and matrix cracking. Although the numerical tool was capable of predicting the impact behaviour, the computational cost was still very expensive. The large number of elements in conjunction with their small characteristic length led to a very large simulation that required a special cluster to be solved in a reasonable time. Considering the steady increase of computer power with time, it is expected that these calculations can be done on standard computers in a few years.

The model was capable to evaluate the influence of the ply clustering on the overall projected delamination shape and position. Moreover, the classical delamination peanut shapes oriented in the fibre direction of the lower layer were well reproduced for the different laminates. Even the increase of the number of interfaces, with a consequent increase of the model complexity, did not affect the result accuracy. For the different impact energies, the model was also able to predict the impact response quite well, but a reduction of the model accuracy was shown by increasing the impact energy. The total absorbed energies for the different configurations were compared with the experimental ones; good agreement was shown for the peak and the initial part of the curves. The dissipated energy for the damage was also reported, giving a powerful tool to directly compare the different configurations. Surprisingly, the total energy dissipated for the damage onset and propagation was lower for the laminate that presented the highest ply clustering. Further research must be performed to better understand these phenomena.

The results clearly show that the numerical model using only simple and typically available measured material data was able to predict the impact behaviour for the different energies and different stacking sequences. A powerful and simple tool to analyse the impact behaviour of thick R-glass fibre composite using the classical in-plane material properties was presented and verified against experimental tests.

Footnotes

Acknowledgements

This work is part of the collaborative project ‘Composite structures under impact loading’ with the industrial partners Flowtite Technology AS, Nammo Raufoss AS and Ragasco AS and the research institutes Norwegian University of Science and Technology (NTNU), SINTEF Materials and Chemistry and SINTEF Raufoss Manufacturing. The authors would also like to thank all partners in the project for constructive discussions.

Funding

The authors would like to express their thanks for the financial support by the Norwegian Research Council (grant 193238/i40) and the industrial partners.

Conflict of interest

None declared.

References

Abrate

. Impact on composite structures, New York, USA: Cambridge University Press, 1998, p.289.

Sutherland

Guedes Soares

. The effects of test parameters on the impact response of glass reinforced plastic using an experimental design approach. Compos Sci Technol 2003; 63: 1–18.

Sutherland

Guedes Soares

. Impact behaviour of typical marine composite laminates. Compos Part B Eng 2005; 37: 89–100.

Sutherland

Guedes Soares

. Impact characterisation of low fibre-volume glass reinforced polyester circular laminated plates. Int J Impact Eng 2005; 31: 1–23.

Sutherland

Guedes Soares

. Impact on low fibre-volume, glass/polyester rectangular plates. Compos Struct 2005; 68: 13–22.

Olsson

. Analytical prediction of large mass impact damage in composite laminates. Appl Sci Manuf 2001; 32(9): 1207–1215.

Naik

Chandra Sekher

Meduri

. Damage in woven-fabric composites subjected to low velocity impact. Compos Sci Technol 2000; 60(5): 731–744.

Zhou

. The use of experimentally-determined impact force as a damage measure in impact damage resistance and tolerance of composite structures. Compos Struct 1998; 42(4): 375–382.

Abrate

. Modeling of impacts on composite structures. Compos Struct 2001; 51(2): 129–138.

10.

Olsson

Falzon

. Delamination threshold load for dynamic impact on plates. Int J Solids Struct 2006; 43(10): 3124–3141.

11.

Gong

Toh

Shim

VPW

. The elastic response of orthotropic laminated cylindrical shells to low-velocity impact. Compos Eng 1994; 4(2): 247–266.

12.

Gong

Shim

VPW

Toh

. Impact response of laminated shells with orthogonal curvatures. Compos Eng 1995; 5(3): 257–275.

13.

González

Maimí

Camanho

. Simulation of drop-weight impact and compression after impact tests on composite laminates. Compos Struct 2012; 94: 3364–3378.

14.

Bouvet

Castanié

Bizeul

. Low velocity impact modelling in laminate composite panels with discrete interface elements. Int J Solids Struct 2009; 46: 2809–2821.

15.

Zemba

Okabe

. A new cohesive model for simulating delamination propagation in composite laminates under transverse loads. Mech Mater 2008; 40: 920–935.

16.

de Moura

MFSF

Gonçalves

JPM

. Modelling the interaction between matrix cracking and delamination in carbon–epoxy laminates under low velocity impact. Compos Sci Technol 2004; 64: 1021–1027.

17.

Reid

Zou

. Modelling damage of multiple delaminations and transverse matrix cracking in laminated composites due to low velocity lateral impact. Compos Sci Technol 2006; 66: 827–836.

18.

Zhang

Zhu

Lai

. Finite element analysis of low-velocity impact damage in composite laminated plates. Mater Des 2006; 27: 513–519.

19.

Choi

Chang

F-K

. A model for predicting damage in graphite/epoxy laminated composites resulting from low-velocity point impact. J Compos Mater 1992; 26: 2134–2169.

20.

Lee

Soutis

. Prediction of impact-induced fibre damage in circular composite plate. Appl Compos Mater 2005; 12: 109–131.

21.

Cheng

. Low-velocity impact-induced damage of continuous fiber-reinforced composite laminates. Part II. Verification and numerical investigation. Compos Part A Appl Sci Manuf 2002; 33: 1063–1072.

22.

Yin

. Low-velocity impact-induced damage of continuous fiber-reinforced composite laminates. Part I. An FEM numerical model. Compos Part A Appl Sci Manuf 2002; 33: 1055–1062.

23.

Camanho

Davila

de Moura

. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003; 37(16): 1415–1438.

24.

Nishikawa

Okabe

Takeda

. Numerical simulation of interlaminar damage propagation in CFRP cross-ply laminates under transverse loading. Int J Solids Struct 2007; 44: 3101–3113.

25.

Borg

Nilsson

Simonsson

. Simulation of low velocity impact on fiber laminates using a cohesive zone based delamination model. Compos Sci Technol 2004; 64: 279–288.

26.

Aoki

Suemasu

Ishikawa

. Damage propagation in CFRP laminates subjected to low velocity impact and static indentation. Adv Compos Mater 2007; 16: 45–61.

27.

Geubelle

Baylor

. Impact-induced delamination of composites: A 2D simulation. Compos Part B Eng 1998; 29: 589–602.

28.

Aymerich

Dore

Priolo

. Simulation of multiple delaminations in impacted cross-ply laminates using a finite element model based on cohesive interface elements. Compos Sci Technol 2009; 69: 1699–1709.

29.

Aymerich

Dore

Priolo

. Prediction of impact-induced delamination in cross-ply composite laminates using cohesive interface elements. Compos Sci Technol 2008; 68: 2383–2390.

30.

Turon

Camanho

Costa

. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006; 38: 1072–1089.

31.

Maimí

Camanho

Mayugo

. A continuum damage model for composite laminates: Part I – Constitutive model. Mech Mater 2007; 39: 897–908.

32.

Maimí

Camanho

Mayugo

. A continuum damage model for composite laminates: Part II – Computational implementation and validation. Mech Mater 2007; 39: 909–919.

33.

Pinho

. Failure models and criteria for FRP under in-plane or three-dimensional stress states including shear non-linearity. 2005. NASA Technical Memorandum 213530 (2005): 18.

34.

DNV. Composite Components. Offshore Standard. 2010, p. 164.

35.

DNV. Composite Risers. Offshore Standard. 2010, p. 46.

36.

Puck

Kopp

Knops

. Guidelines for the determination of the parameters in Puck’s action plane strength criterion. Compos Sci Technol 2001; 62: 371–378.

37.

Puck

Schürmann

. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 1998; 58: 1045–1067.

38.

Puck

Schürmann

. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 2001; 62: 1633–1662.

39.

Hashin

. Failure criteria for unidirectional composites. J Appl Mech 1980; 47: 329–334.

40.

Perillo

Shchebetov

Echtermeyer

. Static ply properties: HiPer-tex™/Epoxy composite. Composites and polymers report series, Trondheim: NTNU – Department of Engineering Design and Materials, 2012, pp. 12–12.

41.

ASTM D7136/D7136M. Standard test method for measuring the damage resistance of a fiber-reinforced polymer matrix composite to a drop-weight impact event.

42.

Perillo

Vedvik

Echtermeyer

. Numerical analyses of low velocity impacts on composite. Advanced modelling techniques. In: Proceedings of the Simulia customer conference – 2012. 2012 2012.

43.

Perillo

Vedvik

Echtermeyer

. Numerical application of three-dimensional failure criteria for laminated composite materials. In: Proceedings of the Simulia customer conference – 2011. 2011 2011.

44.

Benzeggagh

Kenane

. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 1996; 56: 439–449.

45.

Abaqus I. Abaqus Version 6.11 User’s manual, 2011.

46.

Abaqus I. Abaqus Version 6.11 Theory manual, 2011.

47.

Turon

Dávila

Camanho

. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng Fract Mech 2007; 74: 1665–1682.

48.

Lopes

Camanho

Gürdal

. Low-velocity impact damage on dispersed stacking sequence laminates. Part II: Numerical simulations. Compos Sci Technol 2009; 69: 937–947.

49.

Abaqus I. Abaqus 6.11 User’s manual, 2011.

50.

Lancaster

. The effect of carbon fibre reinforcement on the friction and wear of polymers. J Phys D Appl Phys 1968; 1(5): 549–549.

51.

Schön

. Coefficient of friction of composite delamination surfaces. Wear 2000; 237: 77–89.

52.

Zhou

Davies

GAO

. Impact response of thick glass fibre reinforced polyester laminates. Int J Impact Eng 1995; 16: 357–374.

53.

Davies

GAO

Zhang

. Impact damage prediction in carbon composite structures. Int J Impact Eng 1995; 16: 149–170.

54.

Davies

GAO

Hitchings

Zhou

. Impact damage and residual strengths of woven fabric glass/polyester laminates. Compos Part A Appl Sci Manuf 1996; 27: 1147–1156.

55.

Sutherland

Guedes Soares

. Effects of laminate thickness and reinforcement type on the impact behaviour of E-glass/polyester laminates. Compos Sci Technol 1999; 59: 2243–2260.

56.

Sutherland

Guedes Soares

. Impact tests on woven-roving E-glass/polyester laminates. Compos Sci Technol 1999; 59: 1553–1567.