Rate dependent compressive failure and delamination growth in multidirectional composite laminates

Abstract

A novel intralaminar model has, for the first time, been applied and validated for the rate-dependent failure of multidirectional carbon/epoxy laminates. Quasi-static compressive failure is evaluated by the growth of intralaminar rate-dependent damage combined with the interaction of cohesive zones for interlaminar delamination. A special feature of the intralaminar model is the homogenised ply response, allowing simultaneous damage-degradation of the polymer matrix combined with the fibres. To model the observed quasi-brittle failure response of the plies under finite deformation, we have used a viscoelastic-viscoplastic matrix combined with damage and isotropic hardening behaviour. Elastic transverse isotropy is used to model the fibre reinforcement of the plies. Standard cohesive surfaces are used to model the initiation and propagation of delamination. Numerical simulations using ABAQUS/Explicit are performed to predict the growth and delamination of intralaminar damage under compression in different laminates with 56 plies of IM7/8552 carbon/epoxy. Predictions of stress versus strain and damage growth are shown to agree well with experimental results for a range of strain rates and stacking sequences.

Keywords

Angle-ply laminate continuum damage viscoelasticity-viscoplasticity isotropic hardening delamination cohesive surface

Introduction

Due to their high stiffness and strength, laminated composite materials are widely used in modern lightweight structures such as aerospace structures.^1,2 In many structural applications, composite materials are exposed to various dynamic loads such as the strike of birds on aerostructures³ and the crash of large vehicles,⁴ producing multiaxial dynamic states of stress. In addition, the strain rate has significant effects on mechanical properties, such as strength and stiffness.^5–7 Several authors investigated the effect of strain rates on unidirectional compounds,^8–10 as reviewed by Singh¹¹ and Sierakowski.¹² In order to accurately predict the initiation and propagation of intra- and interlaminar dynamic damage, the establishment of a strain-rate-dependent constitutive model at the ply level of the laminate is crucial. Some models based on phenomenological invariance exist and describe viscoelastic/viscoplastic effects in polymer composites, but without damage growth.^13–15 Continuum damage-based models were proposed by,^16,17 including fibre kinking theory, for example.^18–20 Larsson et al.²¹ proposed a micromechanics-based model dependent on the strain rate for unidirectional polymer composites at the ply level. This model was later extended to include the continuum damage formulation²² and was applied to different material systems²³ to demonstrate the robustness of the model. In the model, the polymer matrix is assumed to be viscoelastic-viscoplastic with damage, whereas the carbon fibres are transversely isotropic elastic.

Fibre reinforced composites usually demonstrate intra- and interlaminar damage modes that develop and interact with each other in a complicated way. To exploit the damage mechanisms of composite materials under dynamic loading, many failure prediction models have been developed, and the progressive damage method and the cohesive zone method have been commonly adopted^24,25 to evaluate complicated damage patterns in composite structures. Delamination is commonly considered to be an important failure mechanism of composite laminates with weak interfaces that severely decrease the stiffness and strength of the structures. Delamination problems are often represented by the initiation and propagation of ply separation on the basis of cohesive zone models. In finite element (FE) analysis, this is often achieved through cohesive elements/surfaces.^26,27 Cohesive zone models describe the relationship between traction and displacement jump in the cohesive zone ahead of a crack tip. Several shapes have been proposed for the involved traction-separation law, including bilinear, exponential²⁸ and trapezoidal laws.²⁹ The objective of this paper is to apply the strain rate-dependent damage model developed in^22,23 while including isotropic hardening to realistic multidirectional carbon/epoxy laminates with multiple layers subjected to static and dynamic compressive loading rates. The current model version²³ excludes damage growth in the viscoelastic range, included in,²² as it was found to result in excessive softening in disagreement with experimental observations for typical epoxy resins. The intralaminar model is combined with a cohesive model to validate the growth of intra- and interlaminar damage. To describe the compressive failure of the composite, two components are proposed: a rate-dependent model for the quasistatic intralaminar response and a surface-based cohesive zone model for characterising the interlaminar delamination. In this paper, a bilinear cohesive law is used. When applied to large structures, the shape of the cohesive law has been shown not to influence the overall response, as long as the area below the traction separation law is equal to the resistance to fracture and the element size is not larger than the process zone.³⁰ Therefore, for simplicity, the bilinear traction separation laws^31,32 are the most commonly used to describe delamination in composite materials. The laminate model is implemented as a material subroutine in the finite element code ABAQUS/Explicit. The proposed model is validated by detailed comparisons between numerical predictions and experimental results for laminates with 56 plies in different layups.

Strain-rate-dependent damage model

To capture the nonlinear dynamic response of a composite structure, a strain-rate-dependent damage prediction model based on the progressive continuum damage and the cohesive law is established. Relevant modelling details are presented in this section. The following sub-sections briefly present the constitutive laws used to describe the behaviour of the fibre and matrix material of the fibre reinforced angle plies.

Constitutive model for carbon fibres

The rectangular RVE used to develop the micromechanical relations is denoted V_⊡ as shown in Figure 1(a). The fibre region is indicated by $B_{⊡}^{0, f} \in B_{⊡}^{0}$ , and the polymer matrix region is indicated by $B_{⊡}^{0, m}$ in $B_{⊡}^{0}$ . Fibre orientation is defined by the orientation tensor M = N ⊗ N .

Figure 1.

Schematic representation of (a) a representative volume element (RVE) and (b) a rheological model for the rate-dependent effective stress deviator of the polymer matrix.

The fibres are considered transversely elastic isotropic, characterised by the strain energy ψ^f = ψ^f ( M , E ), where E is the Lagrange strain. Total Kirchhoff stress in the fibre is defined by

τ^{f} = \bar{F} \cdot S^{f} \cdot {\bar{F}}^{t} = τ_{d}^{f} + τ_{v}^{f} + τ_{s}^{f} + τ_{a}^{f} with S^{f} = \frac{\partial ψ^{f}}{\partial E}

(1)

where

\bar{F}

is the (macroscopic) deformation gradient and S ^f is the corresponding second Piola Kirchhoff stress. The fibre strain energy ψ^f contains strain energies that correspond to deviatoric stress

τ_{d}^{f}

, volumetric stress

τ_{v}^{f}

, longitudinal (in-plane) fibre shear

τ_{s}^{f}

and uniaxial fibre response

τ_{a}^{f}

related to orientation tensor M , see.²¹

Dynamic model with damage/hardening for polymer matrix

The matrix is assumed viscoelastic-viscoplastic, enhanced by continuum damage to describe the degradation due to shearing in the polymer, as proposed in Ref. 22. In the present paper, the viscoplastic model is extended to include isotropic hardening. Figure 1(b) shows the rheological model for the matrix material. The viscoelastic model is based on a generalised three-parameter Maxwell model defined by quasistatic and dynamic shear moduli G₁, G₀ (via G₂) and a viscoelastic damper with relaxation time t_2*. This is coupled to a viscoplastic model defined by a quasistatic pressure-sensitive yield in the polymer (with parameters σ_t and σ_c representing uniaxial tension and compression), a viscoplastic relaxation time t_* and the constant hardening modulus H. The viscosity parameters act as “delay” parameters, activating the viscoplastic behaviour at “high” stresses. Following,²² the total Kirchhoff stress in the matrix is obtained as

τ^{m} = \bar{F} \cdot S^{m} \cdot {\bar{F}}^{t} = f [α] ({\hat{τ}}_{1, d}^{m} + {\hat{τ}}_{2, d}^{m}) + τ_{v}^{m} 1 with S^{m} = \frac{\partial ψ^{m}}{\partial E}

(2)

where S ^m is the second Piola Kirchhoff stress and

{\hat{τ}}_{d}^{m} = {\hat{τ}}_{1, d}^{m} + {\hat{τ}}_{2, d}^{m}

is the effective stress deviator and

τ_{v}^{m}

is the mean stress of the matrix. In view of the rheology in Figure 1(b), the hypoelastic-inelastic response is postulated in the rate form of the effective Kirchhoff stresses as

\begin{array}{lc} \overset{□}{\hat{τ}}_{1, d}^{m} = 2 G_{1} I_{d}^{sym} : (d - d_{p}) \end{array}

(3a)

\begin{array}{lc} \overset{□}{\hat{τ}}_{2, d}^{m} = 2 G_{2} I_{d}^{sym} : (d - d_{p} - d_{v}) \end{array}

(3b)

\begin{array}{l} τ_{v}^{m} = K \log [J] \end{array}

(3c)

where

\overset{□}{\hat{τ}}^{m} : = \dot{\bar{\bar{F} \cdot S^{m} \cdot {\bar{F}}^{t}}}

is the Oldroyd rate of the effective stress and

I_{d}^{sym}

is the symmetric fourth order deviatoric projection tensor. Stress rates are related to the rate of deformation tensor d in terms of the shear moduli G₁ and G₂ and the bulk modulus K of the matrix material. d _v and d _p are the viscoelastic and viscoplastic counterparts of d . Furthermore, f [α] = (1 − α)² + r is the damage degradation function, where α is the damage variable and r is the residual stress. This only degrades the matrix shear response, as seen in equation (2). Hardening is represented by the microhardening stress κ associated with linear isotropic hardening defined as

\dot{κ} = H \dot{k}

(4)

where k is the isotropic hardening variable.

Evolution rules

Following the developments in,^22,23 the viscoplastic/viscoelastic and hardening rules for the matrix shear response are defined as

\begin{array}{l} d_{p} = λ \frac{\partial ϕ^{*}}{\partial {\hat{τ}}_{d}^{m}} = λ f with f = \frac{3}{2} \frac{{\hat{τ}}_{d}^{m}}{{\hat{τ}}_{e}^{m}} \end{array}

(5a)

\begin{array}{l} d_{v} = \frac{1}{2 G_{2} t_{2 *}} {\hat{τ}}_{2, d}^{m} \end{array}

(5b)

\begin{array}{l} \dot{k} = λ \end{array}

(5c)

The yield and plastic functions ϕ and ϕ* are defined as

ϕ = {\hat{τ}}_{e}^{m} - (c_{y} + γ p + κ), ϕ^{*} = {\hat{τ}}_{e}^{m} - c_{y} with {\hat{τ}}_{e}^{m} = \sqrt{\frac{3}{2}} | {\hat{τ}}_{d}^{m} | and p = - τ_{v}^{m}

(6)

where

{\hat{τ}}_{e}^{m}

is the effective von Mises stress, c_y is the cohesive yield stress parameter and γ is the friction parameter. These parameters are defined by the stresses σ_t and σ_c as defined in.^22,23

A Perzyna-type Bingham model is adopted for viscoplasticity. λ is the viscoplastic multiplier defined as

λ = \frac{1}{t_{*}} η [ϕ] \geq 0 with η [ϕ] : = \frac{< {\hat{τ}}_{e}^{m} - (c_{y} + γ p + κ) >}{3 G_{0}}

(7)

where t_* is the viscoplastic relaxation time and η is the Bingham overstress function.

Following,^22,23 the progression of α for the damage degradation function f [α] in equation (2) is rate-dependent and is driven by viscoplastic dissipation in the matrix. The damage evolution is defined as

l_{c} \dot{α} = v^{*} ⟨ α^{stat} [α] - α ⟩ with α^{stat} = \frac{A_{T} [α] l_{c}}{G_{c}}

(8)

The damage parameters involved are: the fracture energy $G_{c}$ , the internal length parameter l_c, and the damage progression speed parameter v*. α^stat represents the source of damage energy due to inelastic deformation. $A_{T}$ is the damage driving energy consisting of only plastic work $B$ , obtained by reformulating the dissipation rate $D^{m}$ in^22,23 in the total dissipation defined as

D_{T}^{m} : = \int_{0}^{t} D^{m} d t = f [α] {\hat{W}}^{p} + ψ_{h}^{m} + \int_{0}^{α} A_{T} d α \geq 0 with ψ_{h}^{m} = \int_{0}^{k} κ d k

(9)

A_{T} = B with B = - f [α] {\hat{W}}^{p}

(10)

where

{\hat{W}}^{p}

\int_{0}^{t} {\hat{w}}^{p} d t

is the total (integrated) effective plastic work dissipated during the damage evolution process, which initiates as viscoplasticity commences. Also note that part of the total dissipation is stored in the microstructure of the material corresponding to the irreversible stored hardening energy

ψ_{h}^{m}

Homogenized UD-ply stress response

Following the developments in,^22,23 it appears that the homogenized UD-ply stress response is obtained explicitly as

\bar{S} = v^{m} S^{m} + v^{f} S^{f} + v^{m} (a^{m} F + b^{m} G) : (S^{m} - S^{f})

(11)

where v^m and v^f are the matrix and fibre volume fraction. S ^m and S ^f are the piecewise constant stress fields in the matrix and the fibres. a^m and b^m are scalar microfields in the matrix representing the microscopic strain variation. To express the homogenised response in spatial form, the homogenised Kirchhoff stress is obtained by push-forward transformation

\bar{τ} = \bar{F} \cdot \bar{S} \cdot {\bar{F}}^{t}

. The parameters a^m & b^m are identified by solving local stress equilibrium conditions²³ defined as

δ a^{m} v^{m} f : (S^{m} - S^{f}) = 0, δ b^{m} v^{m} g : (S^{m} - S^{f}) = 0

(12)

A major feature of homogenisation is that the homogenised response of a unidirectional composite ply assumes isostrain response along the fibres, whereas isostress response is assumed for the transverse normal and in-plane shear directions of the ply. To improve the response of the elastic model, Halpin-Tsai-based predictions for the numerical model are used to consider the nonuniform stress distribution in the matrix.²³ Furthermore, in²³, the isostress field has been adjusted by separating normal stresses controlled by the parameter a^m and shear stresses controlled by the parameter b^m.

Interlaminar traction-separation model

To capture the delamination of the ply to the ply under dynamic loading conditions, the standard surface-based cohesive zone method in Abaqus is used. Surface-based cohesive behaviour is chosen primarily because the interface thickness is negligible compared to the thickness of the ply. Thus, a cohesive zone model is used to degrade the traction vector of the interlaminar interface. The model consists of a traction-stick-separation law. An elastic penalty is used to handle the stick case, while a stress-based delamination initiation criterion is used along with a fracture mechanics-based delamination propagation to describe traction-separation.^31,33

Material and parameters

Uniaxial compression tests were performed on angle-ply laminates from an IM7/8552 carbon epoxy composite under quasi-static and dynamic loads by Schaefer et al.¹⁰ Angle-ply samples of [±θ]_14s were produced, where θ = 15°, 30°, 45°, 60°, 75°and 90° alternating the orientation of the layer after each ply during the lay-up.

The ±90° laminate is simply a unidirectional plate tested in the direction transverse to the fibres. The IM7/8552 prepreg was cured to produce a consistent composite laminate with a fibre volume fraction, v_f = 0.58 and a layer thickness, t = 0.14 mm. The samples of dimensions 7.82 × 7.82 × 25.4 mm were then cut from the laminate. Strain rates greater than 500/s were achieved in dynamic tests, carried out in split Hopkinson bar setups, while in quasi-static loading, a servohydraulic test machine (Instron 8500) was used to test the strain rate value of the order of 10⁻⁴/s.

The applied boundary conditions are represented in Figure 2(a) and correspond to the experimental testing boundary conditions in compression. The applied velocity is $v = \dot{ϵ} L$ , where $\dot{ϵ}$ is the applied strain rate (Table 1) and L is the specimen length (Figure 2(b)). Table 1 lists the tests used in the calibration and validation of the presented constitutive model. To replicate the experimental tests, a velocity boundary condition is applied on one end of the sample (see Figure 2(a)). With regard to the material properties used in this paper, it is assumed that the carbon fibres are linear elastic and transversely isotropic, and that the matrix follows the viscoelastic-viscoplastic model with damage discussed in Section on page 4. The properties of the elastic material of the constituents are given in Tables 2 and 3. Some of the properties are obtained from the manufacturer’s data sheet,^34,35 while others are estimated.

Figure 2.

Considered angle-ply composite laminate subjected to compression velocity loading. (a) Geometry of the specimen with boundary conditions and (b) 3D FE model with cohesive surfaces between each ply.

Table 1.

Specifications of angle-ply laminates used for numerical validation under quasistatic and dynamic loading.¹⁰

Fibre angle [±θ]	Strain rates ( $\dot{ϵ}$ ) in compression [/s]
30°	1 × 10⁻⁴(QS)	500
60°	1 × 10⁻⁴(QS)	500
90°	1 × 10⁻⁴(QS)	800

Table 2.

Material parameters for the carbon fibre.

$E_{L}^{f}$	$E_{T}^{f}$	$G_{L T}^{f}$	$ν_{L T}^{f}$	$ν_{T T}^{f}$	v ^f
GPa	GPa	GPa	–	–	%
276	25	30	0.25	0.5	58

Table 3.

Material parameters for the epoxy at quasi-static loading.

G ₁	ν ^m	σ _t	σ _c	$G_{c}$
GPa	–	MPa	MPa	N/mm
1.24	0.37	100	130	2.1

In simulations, an effective quasistatic shear modulus $(G_{1}^{⋆} = 1.4 G_{1})$ is used to improve the prediction of elastic stiffness under in-plane shear and transverse loading. Similarly, simulations are also based on an effective toughness $(G_{c}^{⋆} = 3.1 G_{c})$ because intralaminar damage is caused by the shearing of the matrix material.³⁶ $G_{c}$ is the toughness (Mode I) of the neat resin, reported as 0.68 N/mm³⁵. The definition of effective quantities $G_{1}^{⋆}$ and $G_{c}^{⋆}$ has been thoroughly explained by Singh et al.²³

A standard surface-based cohesive zone method in Abaqus is used for delamination modelling. A bilinear elastic traction-separation behaviour together with maximum stress based delamination initiation and energy based delamination evolution is used. Various cohesive strength values have been assumed by research groups under different loading cases for IM7/8552 composite laminates.³⁷ Interface strengths

(τ_{I}^{c}, τ_{I I}^{c}, τ_{I I I}^{c})

are used in delamination initiation criteria based on maximum stress. The assumed values for interfacial strength and interlaminar fracture toughness are reported in Table 4, from.^36–38 The values of the interfacial strengths are assumed to be similar to the values of the intralaminar strength in the plane for the composite laminate IM7/8552. This assumption is based on the fact that the delamination begins and spreads rapidly, so it becomes difficult to accurately measure the opening of the tip of the crack.³⁷ Mode I interlaminar fracture toughness (G_IC) is obtained from³⁸ and Mode II and III interlaminar fracture toughnesses (G_IIC and G_IIIC) are based on the arguments of.³⁶ These parameters are used in the energy-based delamination criterion by Benzeggagh-Kenane (B-K).

Table 4.

Parameters of surface-based cohesive behaviour.

Interface strength	Interlaminar fracture toughness
MPa	N/mm
$τ_{I}^{c} = 60$	G_IC = 0.2
$τ_{I I}^{c} = τ_{I I I}^{c} = 90$	G_IIC = G_IIIC = 0.6

Calibration of model parameters at the ply level

An approach based on numerical optimization, using the least square method, is applied to estimate viscoelastic-viscoplastic parameters for the epoxy matrix material, as explained in.^22,23 In this approach, a calibration routine finds the dynamic model parameters by minimizing the error between homogenized material response at a material point and the uniaxial stress-strain response from the experiments. The identification of the model parameters at the ply level is based on calibration of the 30° off-axis specimen in compression.⁹ For calibration, the constitutive FORTRAN VUMAT driver for the ply response is implemented in MATLAB to handle uniaxial compression at the material point level relative to the compression tests in.⁹ Note that the experimental strains were measured locally using strain gauges. Therefore, calibration is based on linking the strain at the material point to the experimental strain gauge signal, i.e. assuming uniform strains. We emphasise that calibration is considered up to the point (usually peak stress) where strain localisation begins in the test specimen.

The calibrated parameters are summarised in Table 5. The calibrated parameters include dynamic shear modulus G₂, viscoelastic damper t_2* of the viscoelastic model, viscoplastic damper t_* and isotropic hardening H for the viscoplastic model. In addition to them, two more parameters are needed for the initiation and evolution of intralaminar damage. These parameters are: internal length l_c and fracture area propagation velocity v*. The procedure to determine the model parameters has been explained in detail in.^22,23 The calibrated response for the off-axis case 30° is shown in Figure 3, corresponding to the (calibrated) parameters in Table 5. The experiments and model responses in Figure 3 are shown by solid and dashed lines, whereas the quasistatic and dynamic responses are solid blue and red curves. The damage growth predicted by the model is shown by semi-solid lines (dotted dashes).

Table 5.

Estimated model parameters for neat epoxy matrix material.

G ₂	t _2*	t _*	H	l _c	v*
GPa	s	s	GPa	mm	mm/s
1.84	5 × 10⁻⁶	1 × 10⁻⁵	1.48	0.005	$\dot{ϵ} \times 25.4 \times 5$

Figure 3.

Comparison of experiments and uniaxial stress-strain predictions of the material points in compression under quasistatic and dynamic loading for a 30°off-axis specimen. The coloured stars represent the onset of yielding in the matrix.

Numerical validation

Finite element analysis

The numerical results are presented in this section and compared to the experimental results.

Based on published experimental data,^9,10 including laminate specifications and experimental test results, finite element simulation is performed using the Abaqus/Explicit solver. The strain-rate-dependent damage model is implemented through a user-defined material (VUMAT) subroutine. The compressive behaviour of IM7/8552 angle-ply composite laminates is investigated. The laminate specifications and mechanical properties of the constituents are shown in Tables 1–3, respectively.

The explicit integration scheme, based on backward Euler method, element deletion technique, together with enhanced hourglass, is used. Figure 2(b) shows the numerical finite element model. The laminate is modelled with one element per ply along the thickness direction. Reduced integration of 8-node solid elements C3D8R is used to discretise the model. Convergence is achieved with a solid element of 0.3 × 0.3 × 0.3 mm³ which results in about 2200 elements to discretize the single ply as shown in Figure 2(b). Cohesive surfaces are defined for all ply-to-ply contact interfaces. The parameters of the behaviour of the cohesive interface are reported in Table 4.

The model parameters are estimated based on the 30° off-axis experiments using calibration as discussed in Section on page 11. This set of model parameters is used to validate the material model by predicting intra- and interlaminar failure of IM7/8552 angle-ply composite laminates at different compression strain rates. An average uniaxial strain in the xy plane and an average axial stress on the yz plane are used to represent the uniaxial stress-strain curves, as shown in Figures 4(a)–6(a). In Figures 4(b) and 5(b), the view represents the top ply seen from below. The comparison between numerical predictions and experimental observations of a 30° laminate under quasistatic and dynamic loading and under compression is presented in Figure 4. Reasonable agreement is observed between the numerical and experimental curves (Figure 4(a)). Note that the softening and subsequent drop in average stress correspond to the localisation of the strain. Previous studies of mesh convergence for the rate-dependent damage model (8) have shown that the localised damage evolution is mesh objective.^20,23 The stiffness predictions agree well with the experimental observations for the quasi-static and dynamic loadings. An overestimation of the strength prediction is observed under a dynamic load of 500/s. Figures 4(b)–4(c) show the FE predictions for dynamic intra- and interlaminar failure modes of a laminate ±30° at 500/s. A localised intralaminar failure is observed transverse to the loading direction. The localised deformation occurs on the loading surface. Figure 4(c) shows the distribution of out-of-plane intralaminar damage and delamination growth through the thickness of the laminate. The delamination begins and is mainly distributed in the upper plies close to the loading surface because of the lack of out-of-plane constraint at the free surface.

Figure 4.

Experimental and predicted stress-strain responses, FE predictions for dynamic damage growth in ±30° composite laminate in compression. Damage distribution maps α[ x ] are indicated by coloured circles in 4a. (a) Averaged uniaxial stress-strain response under quasi-static and dynamic loading. (b) In-plane (X-Y) intralaminar damage growth at 500/s. (c) Out-of-plane (X-Z) intralaminar damage growth within the plies (top) and delamination in cohesive surfaces (bottom) at 500/s.

Figure 5.

Experimental and predicted stress-strain responses, FE predictions for dynamic damage growth in ±60° composite laminate in compression. The damage distribution maps α are indicated by the coloured circles in 5a. (a) Averaged uniaxial stress-strain response under quasi-static and dynamic loading. (b) In-plane (X-Y) intralaminar damage growth at 500/s. (c) Out-of-plane (X-Z) intralaminar damage growth within the plies (top) and delamination in cohesive surfaces (bottom) at 500/s.

Figure 6.

Experimental and predicted stress-strain responses and FE predictions for dynamic failure in ±90° composite laminate in compression. The strain distribution maps, ϵ_xx, are indicated by the coloured circles in 6a. Experimental images by Koerber et al.⁸ ⓒ Elsevier. Reproduced with permission. (a) Averaged uniaxial stress-strain response under quasi-static and dynamic loading. (b) In-plane (X-Y) strain ϵ_xx distribution predicted for 800/s and experimental strain measurement for 276/s. (c) Out-of-plane (X-Z) inclined shear bands predicted for 800/s and fractographic observation for 276/s

Figure 5(a) shows the averaged uni-axial stress-strain response at strain rates of 10⁻⁴/s and 500/s in compression of ±60° laminate. Under quasistatic and dynamic loading, good agreement is observed for stiffness and strength prediction. The ±60° laminate failed in the plane and through the thickness. In Figure 5(a) and (b) localised intralaminar failure is observed in the plane. Failure in the plane begins in the corner region of the sample and spreads along the loading direction, due to the forced rotation of the ±60° fibres in the upper layer and the lack of support on the free surfaces. Damage to the interlaminar delamination is seen in Figure 5(c). A progression of through-thickness damage is observed on the loading surface due to the lack of support at the free top surface. Figure 5(c) further illustrates the interaction between the intralaminar damage detected by the ply model and the interlaminar delamination detected by the cohesive surface model. In Figure 5(c), it can be seen that the delamination begins on the free surface of the specimen. The interactive nature of the out-of-plane intra- and interlaminar failure is also visible in the xz plane. It is well known that the onset of delamination can be related to matrix cracking and that the extent of delaminations can be related to zones with matrix cracking, that is, intralaminar damage. An early modelling study of how matrix cracks initiate delamination was provided in.³⁹ More recently, Zubillaga et al.⁴⁰ developed the MCID criterion for matrix crack-induced delamination.

Figure 6(a) illustrates the experimental and numerical averaged uniaxial stress-strain curves at strain rates of 10⁻⁴/s and 800/s for a ±90° laminate in compression (corresponding to transverse compression). Good agreement is obtained for the prediction of stiffness and an overprediction of strength for a dynamic loading of 500/s. Under quasi-static loading, good agreement is obtained for stiffness and strength prediction. Patterns of the in-plane shear band can be observed in Figure 6(b). Multiple localised in-plane shear bands are obtained transverse to the loading direction at 800/s in the xy-plane. These patterns are in line with the observations of Koerber et al.⁸ with dynamic experimental tests performed on the IM7/8552 material system under transverse compression. Figure 6(c) illustrates the out-of-plane shear localisation that occurs in the ±90° laminate under dynamic compression. Multiple shear bands circa 45° from the loading direction are observed in the xz-plane. A similar observation, but at a lower loading rate, was made by Koerber et al.⁸

Discussion

Figure 7 shows a comparison of experimental observations and numerical predictions for 15°, 45°, 75°and 90° off-axis specimens with IM7/8552 unidirectional plies under quasi-static and dynamic loading and in compression. The representation of the results is similar to that described in the section on page 11.

Figure 7.

Comparison of experiments and material point uniaxial stress-strain predictions in compression under quasi-static and dynamic loading for different off-axis specimens. The coloured stars represent the beginning of the yielding. (a) Comparison of 15° off-axis response. (b) Comparison of 45° off-axis response. (c) Comparison of 75° off-axis response. (d) Comparison of 90° off-axis response.

Based on the FORTRAN VUMAT implementation of the 3D constitutive ply model, a Matlab code was developed to describe the uniaxial compressive ply tests at the material point level. Predictions at the material point level assume homogeneous straining in the specimens, allowing the calibration of the material point level considered to be based on the prelocalized ply behaviour. Therefore, numerical validations are based on model parameters obtained after calibrating on a sample 30° off-axis with quasistatic and dynamic compression loading (see Section on page 11). Figure 7 shows that good agreement is achieved between compressive experimental observations and numerical predictions for different off-axis specimens with IM7/8552 unidirectional plies under quasistatic and dynamic loading. In particular, the quasi-static results are better predicted than in Ref. 22. This is due to the introduction of an additional model parameter, that is, isotropic hardening, which is estimated by calibration, as explained on page 11 (see Table 5). The formulation is described in the section on page 5. This parameter is used to model the strain hardening behaviour of the matrix caused by the reorientation of the covalent chains during plastic deformation, where the straightening of the chains results in strain hardening.⁴¹ We also observe an increase in initial stiffness for the dynamic case compared to the static model responses because of the rate-sensitive viscoelastic polymer matrix constituent. Note that the experimental curves in Figure 7 reflect the average strain in the region where the strain was measured (a central strain gauge), while the predictions for the material point reflect the response up to the ultimate failure in the region with strain localisation. Therefore, the measured failure strain in the experiments is much lower than the ultimate failure strain of the material (see Figure 7).

The resulting model response of the neat resin for quasistatic loading is shown in Figure 8. In view of the expected behaviour of the polymer, it was observed that the model reasonably well predicts the stiffness up to the onset of the yield in compression. Once the yield stress is reached, the model response is then governed by viscoplasticity coupled to isotropic hardening and continuum damage. In general, a good agreement is obtained between the model response and the experiment (reproduced from⁴¹) in compression.

Figure 8.

Comparison of the predictions of the numerical model with the experimental results of the 8552 resin in compression under quasi-static loading.⁴¹ The beginning of the yielding is shown by a coloured star.

Conclusions

A strain rate-dependent ply level damage model has been successfully applied in FE modelling of quasistatic and dynamic compressive failure processes of angle ply laminates from IM7/8552 carbon/epoxy. The investigation was carried out for various strain rates using Abaqus/Explicit, where the interlaminar behaviour was described by a standard Abaqus traction-separation law. The ply-level strain-rate effects were described by a recently reported homogenised material model involving elastic transversely isotropic carbon fibres embedded in a viscoelastic-viscoplastic polymer matrix. The strain hardening was modelled by including isotropic hardening in the viscoplastic response of the matrix. Model predictions were validated against the experimentally observed rate-dependent damage growth of angle-ply laminates.¹⁰ Relatively good qualitative and quantitative correlation was achieved between the stress-strain response of the numerical models and the experimental results for the angle-ply laminates of IM7/8552. The predicted growth of intra- and interlaminar damage to composite laminates is consistent with experimental observations and confirms the effectiveness of the proposed strain-rate-dependent damage model. Note that strains will be relatively uniform until strain localisation, which is followed by rapid damage growth and subsequent final failure. It appears that the actual delamination growth of each interface is generally triggered by growth of the damage at the ply level.

Footnotes

Declaration of conflicting interests

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

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Swedish Foundation for Strategic Research (SSF) through the contract dnr FID16-0041 and by the internal SK-development funds of RISE.

ORCID iDs

Vivekendra Singh

Robin Olsson

Data Availability Statement

The data given in this article is available on request.

References

Kurşun

Şenel

Enginsoy

. Experimental and numerical analysis of low velocity impact on a preloaded composite plate. Adv Eng Software 2015; 90: 41–52.

Delbariani-Nejad

Farrokhabadi

Jafari

. An energy based approach for reliability analysis of delamination growth under mode I, mode II and mixed mode I/II loading in composite laminates. Int J Mech Sci 2018; 145: 287–298.

Georgiadis

Gunnion

Thomson

, et al. Bird-strike simulation for certification of the boeing 787 composite moveable trailing edge. Compos Struct 2008; 86: 258–268.

Savage

Bomphray

Oxley

. Exploiting the fracture properties of carbon fibre composites to design lightweight energy absorbing structures. Eng Fail Anal 2004; 11: 677–694.

Harding

. Effect of strain rate and specimen geometry on the compressive strength of woven glass-reinforced epoxy laminates. Composites 1993; 24(4): 323–332.

Hsiao

Daniel

Cordes

. Strain rate effects on the transverse compressive and shear behavior of unidirectional composites. J Compos Mater 1999; 33(17): 1620–1642.

Gilat

Goldberg

Roberts

. Experimental study of strain-rate-dependent behavior of carbon/epoxy composite. Compos Sci Technol 2002; 62: 1469–1476.

Koerber

Xavier

Camanho

. High strain rate characterisation of unidirectional carbon-epoxy IM7/8552 in transverse compression and in-plane shear using digital image correlation. Mech Mater 2010; 42(11): 1004–1019.

Schaefer

Werner

Daniel

. Strain-rate-dependent failure of a toughened matrix composite. Exp Mech 2014; 54: 1111–1120.

10.

Schaefer

Daniel

. Strain-rate-dependent yield criteria for progressive failure analysis of composite laminates based on the Northwestern failure theory. Exp Mech 2018; 58: 487–497.

11.

Singh

. Literature survey of strain rate effects on composites. TR18-001 (open), Swerea SICOMP (now RISE). Sweden: Mölndal, 2018.

12.

Sierakowski

. Strain rate effects in composites. Appl Mech Rev 1997; 50(12): 741–761.

13.

Eskandari

Andrade Pires

Camanho

, et al. Analyzing the failure and damage of FRP composite laminates under high strain rates considering visco-plasticity. Eng Fail Anal 2019; 101: 257–273.

14.

Gerbaud

Otero

Bussetta

, et al. An invariant based transversely-isotropic constitutive model for unidirectional fibre reinforced composites considering the matrix viscous effects. Mech Mater 2019; 138: 103146.

15.

Koerber

Kuhn

Ploeckl

, et al. Experimental characterization and constitutive modeling of the non-linear stress–strain behavior of unidirectional carbon–epoxy under high strain rate loading. Adv Model Simul Eng Sci 2018; 5: 17.

16.

Camanho

Bessa

Catalanotti

, et al. Modeling the inelastic deformation and fracture of polymer composites – Part II: Smeared crack model. Mech Mater 2013; 59: 36–49.

17.

Maimí

Camanho

Mayugo

, et al. A continuum damage model for composite laminates: part I – constitutive model. Mech Mater 2007; 39: 897–908.

18.

Costa

Fagerström

Olsson

. Development and validation of a finite deformation kinking model for crushing of composites. Compos Sci Technol. 2020; 197: 108236.

19.

Costa

Gutkin

Olsson

. Mesh objective implementation of a fibre kinking model for damage growth with friction. Compos Struct 2017; 168: 384–391.

20.

Larsson

Gutkin

Rouhi

. Damage growth and strain localization in compressive loaded fiber reinforced composites. Mech Mater 2018; 127: 77–90.

21.

Larsson

Singh

Olsson

, et al. A micromechanically based model for strain rate effects in unidirectional composites. Mech Mater 2020; 148(103491): 193–212.

22.

Larsson

Singh

Olsson

, et al. A micromechanically based model for dynamic damage evolution in unidirectional composites. Int J Solid Struct 2022; 238(111368): 111368.

23.

Singh

Larsson

Olsson

, et al. A micromechanics based model for rate dependent compression loaded unidirectional composites. Compos Sci Technol; 2023: 232: 109821.

24.

Shi

Pinna

Soutis

. Modelling impact damage in composite laminates: A simulation of intra- and inter-laminar cracking. Compos Struct 2014; 114: 10–19.

25.

Tan

Falzon

Chiu

, et al. Predicting low velocity impact damage and compression-after-impact (CAI) behaviour of composite laminates. Composites Part A 2015; 71: 212–226.

26.

Crisfield

Davies

GAO

, et al. Progressive delamination using interface elements. J Compos Mater 1998; 32(14): 1246–1272.

27.

Yang

Cox

. Cohesive models for damage evolution in laminated composites. Int J Fract 2005; 133: 107–137.

28.

Siegmund

Ramani

. Rate-dependent crack growth in adhesives. Int J Adhesion Adhes 2003; 23: 9–13.

29.

Marzi

Hesebeck

Brede

, et al. A rate-dependent cohesive zone model for adhesively bonded joints loaded in mode I. J Adhes Sci Technol 2009; 23: 881–898.

30.

Jousset

Rachik

. Comparison and evaluation of two types of cohesive zone models for the finite element analysis of fracture propagation in industrial bonded structures. Eng Fract Mech 2014; 132: 48–69.

31.

Camanho

Davila

de Moura

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

32.

Johnson

Pickett

Rozycki

. Computational methods for predicting impact damage in composite structures. Compos Sci Technol 2001; 61: 2183–2192.

33.

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.

34.

Hexcel Corporation . Hextow IM7 carbon fibre datasheet. Stamford CT, USA: Global version, 2020. https://www.hexcel.com/Resources/DataSheets/Carbon-Fiber

35.

Hexcel Corporation . Hexply 8552 product data sheet. Stamford CT, USA: EU version, 2020. https://www.hexcel.com/Site/Resources/DataSheets/Prepreg.aspx

36.

Xia

Hutchinson

. Mode II fracture toughness of a brittle adhesive layer. Int J Solid Struct 1994; 31(8): 1133–1148.

37.

Ridha

Chen

, et al. On cohesive element parameters and delamination modelling. Eng Fract Mech 2019; 206: 278–296.

38.

Czabaj

Ratcliffe

. Comparison of intralaminar and interlaminar mode I fracture toughnesses of a unidirectional IM7/8552 carbon/epoxy composite. Compos Sci Technol 2013; 89: 15–23.

39.

Choi

HYT

Chang

. A new approach toward understanding damage mechanisms and mechanics of laminated composites due to low-velocity impact: Part II–Analysis. J Compos Mater 1991; 25: 1012–1038.

40.

Zubillaga

Turon

Maimí

, et al. An energy based failure criterion for matrix crack induced delamination in laminated composite structures. Compos Struct 2014; 112: 339–344.

41.

Naya

. Prediction of mechanical properties of unidirectional FRP plies at different environment conditions by means of computational micromechanics. Spain: Polytechnic University of Madrid, 2017. Doctoral thesis.