Comprehensive inter-fibre failure analysis and failure criteria comparison for composite materials using micromechanical modelling under biaxial loading

Abstract

Inter-fibre failure analysis of carbon fibre-reinforced polymer (CFRP) composites, under biaxial loading conditions, has been a longstanding challenge and is addressed in this study. Biaxial failure analysis of IM7/8552 CFRP unidirectional (UD) composites is conducted under various stress states. Two widely accepted failure criteria, the interactive Tsai-Wu and non-interactive Hashin failure criteria, are comprehensively assessed with finite element-based micromechanical analysis. High-fidelity three-dimensional representative volume elements (RVEs) are subjected to biaxial loadings with imposed periodic boundary conditions. Carbon fibres are assumed to be transversely isotropic and linearly elastic. The Drucker-Prager plastic damage constitutive model and cohesive zone model are utilised to simulate the mechanical response of the matrix and fibre-matrix interface, respectively. Coulomb friction is assumed between the fibres and matrix after interface failure. Two sets of biaxial loading scenarios (i.e. transverse stress dominated and shear stress dominated) with the associated failure modes are selected for the failure analysis and assessment of these failure criteria. A data-driven failure envelope for the composites under biaxial loadings is developed using a univariate cubic spline function. Failure mode transition points are determined under biaxial loadings. It is found that the micromechanics-based numerical model is effective in assessing these two existing criteria.

Keywords

Micromechanics finite element analysis representative volume element biaxial loadings failure criteria failure modes

Introduction

Carbon fibre-reinforced polymer (CFRP) composites are increasingly used in aerospace, automotive and marine industries as lightweight structures, due to their excellent strength- and stiffness-to-weight ratios. To mitigate safety concerns that arise in the use of composite materials, efficient material characterisation techniques are required at the early stages of structural design, especially the failure behaviour under multiaxial loading conditions. In contrast to metallics, CFRP composites are hierarchical with three different length scales (i.e. micro-, meso- and macro-). Damage initiation and propagation mechanisms are associated with specific scales, resulting in difficulties in the prediction of composites failure. In addition, the coexistence of various damage modes and failure mechanisms requires various failure criteria, which are not fully validated, especially under a multiaxial stress state.¹ Therefore, failure analysis of composites under multiaxial loading conditions remains challenging, considering the failure mode interaction and the progressive nature of the failure.

Over the years, a large number of failure theories, criteria and models have been proposed to predict the failure of composite lamina/laminates. These failure criteria are predominantly either strain-based^2,3 or stress-based.^1,4–9 Some failure criteria require expensive and time-consuming experimental campaigns to determine the input parameters for different material systems.^10,11 These widely used failure criteria and theories were assessed in the ‘World Wide Failure Exercises (WWFEs)’,^11–14 according to their predictive capabilities of the failure strength of composite lamina/laminates under various loading conditions. However, significant discrepancies were found between the predictions and experimental results in some of the test cases, and no failure criteria/models were deemed to be able to handle all load cases, especially multiaxial cases.¹³ Despite continued research effort, a consensus on selecting the best criteria to predict composite failure has not been reached. Moreover, due to the limitations of experimental fixtures and specimen geometries, multiaxial loading conditions are still challenging to be achieved, resulting in many existing failure criteria remaining unvalidated. Nonetheless, some of these failure criteria/models such as Tsai-Wu⁵ and Hashin⁷ are still widely used in complex loading conditions, e.g. impact.^15,16 Therefore, it is necessary to assess these failure criteria under multiaxial loadings before they are applied to predict the failure of composite structures under complex loadings.

The availability of high-performance computing and high-fidelity numerical approaches has facilitated the virtual testing and failure analysis of composites by means of computational micromechanics. Computational analysis under multiaxial loads avoids the complexity of physical experimental tests. A representative volume element (RVE) is widely used within the framework of the finite element method by considering the response of the constituents (i.e. fibres, matrix and fibre/matrix interface) to infer the mechanical behaviour of composites.¹⁷ Recently, micromechanics-based modelling has been successfully applied to investigate the mechanical behaviour of composites and assess the abovementioned failure criteria and/or models under different biaxial loading conditions, such as combined transverse and out-of-plane shear,^1,18 combined transverse and in-plane shear,^1,19–24 combined transverse compression and axial tension,²⁵ combined longitudinal tension and in-plane shear¹ and triaxial loadings.²⁶ Sun et al.¹ conducted a computational failure analysis of composites under three sets of biaxial loading conditions and assessed several widely used failure criteria. They found discrepancies between computational results and conventional failure criteria, and the reason lies in the transition of dominant failure mechanisms observed in computational results that have not yet been considered in existing failure criteria. Thus a comprehensive assessment considering failure mechanisms transition is necessary to gain confidence in the proposed failure criteria by exploiting the potential of computational analysis.

In this study, computational micromechanics-based RVE models with randomly distributed fibres were developed to predict the inter-fibre failure envelopes and failure modes of IM7/8552 unidirectional (UD) CFRP composite laminae under biaxial loading conditions. Each RVE model was built with three constituents: the fibres were modelled as linearly elastic and transversely isotropic material; the mechanical behaviour of the matrix was modelled with the Drucker Prager plastic damage model considering the influences of hydrostatic pressure, and the interface was modelled using a bilinear cohesive law and the Benzeggath-Kenane damage propagation law for mix-mode energy dissipation. The post-failure behaviour of the interface between fibres and matrix was governed by Coulomb friction and implemented into ABAQUS via a subroutine.²⁷ Six out of 10 representative biaxial loading conditions were considered due to the transverse isotropy of the cross-section of unidirectional composites. Periodic boundary conditions were applied to achieve the biaxial loading conditions. Five RVE models with dimensions 50 μm × 50 μm × 5 μm were validated by experimental results from literature under combined transverse and in-plane shear stress states. Data-driven failure predictions were generated based on representative failure points in each biaxial loading case. This was extended by machine learning techniques considering failure probability in our previous study.²⁷ A comprehensive failure analysis was conducted based on failure envelopes and failure modes. Representative failure modes and failure mechanisms and the transition between different modes and mechanisms are discussed. A detailed comparison is presented between the failure predictions of the two classical failure criteria (Tsai-Wu⁵ and Hashin⁷) and micromechanical simulations. The experimentally validated computational approach provides fast screening capabilities to improve material down-selection as part of a design cycle.

Micromechanical modelling of composites

3D RVE model set up

Computational micromechanics is based on the analysis of a statistically representative volume element (RVE) of the material under tension, compression, shear or combined stresses. The microstructure of the RVE of the unidirectional composite contains randomly distributed parallel and circular fibres embedded in a polymer matrix. 50 fibres were shown to adequately capture the essential microstructural features of the material during the progressive failure process while maintaining reasonable computational effort.²¹ Five RVEs with the same number of fibres and volume fraction, but different fibre distributions, were generated and results were compared to demonstrate that the RVE is statistically representative.

Random fibre distributions were generated using the discrete element method²⁸ with an average fibre diameter of 7 μm and fibre volume fraction of 60%. As shown in Figure 1(a), the microstructures of the generated five RVEs reveal well-distributed fibres without significant fibre clustering or matrix-rich regions. The final 3D RVE model of the unidirectional composite material can be obtained by extruding the 2D model with periodic fibre distributions along the longitudinal direction. The thickness of RVE1 to RVE5 was set to 5 μm considering the balance between numerical results accuracy and computational costs. The size of the cross-section of the RVEs was 50 μm × 50 μm. Fibres and matrix in these RVEs were discretised with first-order hexahedral elements under a reduced integration scheme (C3D8R) and few tetrahedral elements (C3D6), while the fibre/matrix interface was meshed with first-order cohesive elements (COH3D8). Figure 1(b) illustrates the microstructure of the 3D RVE model of a UD composite with three constituents. Typically, the RVE1 model was discretised with around 20,000 elements in ABAQUS/Explicit to capture the stress concentration between neighbouring fibres and ensure the balance between result accuracy and computational cost. To accelerate the simulation process, mass scaling is normally utilised in the ABAQUS/Explicit and the stable time increment was set at 5 × 10⁻⁶ s.^21,26

Figure 1.

(a) Microstructures of the RVEs with different periodic fibre distributions, (b) 3D micromechanics-based RVE model with three constituents.

Constitutive models of constituents

Carbon fibres were modelled to be linearly elastic, brittle and transversely isotropic and were assumed not to contribute to the failure of the composite under transverse and shear dominated loadings. Previous work has shown that the hydrostatic stress has significant influences on the mechanical behaviour of the polymer,²⁹ and exhibits a completely different behaviour when subjected to various simple uniaxial loading conditions, such as brittle in tension while plastic in compression and shear.³⁰ Constitutive models, reported in the literature, to describe the mechanical behaviour of polymeric materials under multiaxial loading³¹ include the extended Drucker-Prager (D-P) yield model, associated with a ductile damage criterion,²¹ the modified Drucker-Prager plastic damage model,¹⁹ and the elasto-plastic model with an isotropic damage constitutive model.³² In this study, the polymer matrix was modelled as an isotropic elastoplastic solid. The modified Drucker-Prager plastic damage model³³ was used to model the mechanical behaviour of epoxy under multiaxial stress states. The yield surface of the epoxy is expressed by

Φ (I_{1}, J_{2}, σ_{I}, β, α) = \frac{1}{1 - α} (\sqrt{3 J_{2}} + α I_{1} + B 〈 σ_{I} 〉) - σ_{m y c} = 0

(1)

where

I_{1}

stands for the first invariant of the stress tensor, J₂ is the second invariant of the deviatoric stress tensor, α is the pressure-sensitivity parameter of the Drucker-Prager yield criterion, σ_I is the maximum principal stress, <> are Macaulay brackets, and B is a function of the tensile and compressive yield stresses (σ_myt and σ_myc), which is defined as

B = \frac{σ_{m y t}}{σ_{m y c}} (1 - α) - (1 + α)

(2)

wherein α can be determined according to tanβ = 3α from the internal friction angle of the material (β), which controls the hydrostatic pressure dependence on the plastic behaviour. For the matrix behaviour under uniaxial tension, the quasi-brittle behaviour is controlled by an exponential cohesive law after damage onset, characterised by a single normalized scalar damage variable, to ensure the correct energy dissipation of the matrix G_m. For matrix behaviour under uniaxial compression, perfect plasticity is assumed based on the experimental findings.¹⁹ Figure 2(a) shows the mechanical response of the epoxy under uniaxial tension and compression. More details about the constitutive model and the numerical implementation are given in.^19,26,33

Figure 2.

(a) Constitutive model of the epoxy matrix, (b) Constitutive model and damage variable D of fibre/matrix interface under pure shear load (black line) and combined compression and shear loads (blue line).

Regarding the modelling of the fibre/matrix interface, a bilinear cohesive model was used to predict its mechanical behaviour and progressive failure by means of cohesive elements. The initial response of the cohesive model is assumed to be linear elastic governed by penalty stiffnesses (K_nn—normal, K_ss—shear longitudinal and K_tt—shear transversal), which should be large enough to ensure displacement continuity. Here a machine learning based approach was used to determine the stiffnesses,^26,34 which can be found in Table 1. Damage onset is controlled by a quadratic interaction criterion, and the damage occurs when the criterion involving the sum of nominal stress ratios reaches one. The criterion is represented as.

{\frac{〈 t_{n} 〉}{t_{n}^{0}}}^{2} + {\frac{t_{s}}{t_{s}^{0}}}^{2} + {\frac{t_{t}}{t_{t}^{0}}}^{2} = 1,

(3)

where

t_{n}^{0}, t_{s}^{0} a n d t_{t}^{0}

represent the interface strengths in the normal and two shear directions. Damage evolution is defined based on the traction separation law, and the dissipated fracture energy is used to determine the separation displacement at failure. When the damage is initiated, the traction stress

t_{i} (i = n, s, t)

starts to decrease according to a single normalised scalar damage parameter, which monotonically increases from 0 (at damage onset

δ_{s}^{0}

) to 1 (at final failure

δ_{s}^{f}

), which is shown in Figure 2(b). The energy-based Benzeggath-Kenane damage propagation criterion³⁵ is adopted to consider the dependence of the fracture energy dissipation on fracture modes during the damage evolution of the cohesive elements, which reads

G^{C} = G_{n}^{C} + (G_{s}^{C} - G_{n}^{C}) {(\frac{2 G_{s}}{G_{n} + 2 G_{s}})}^{η_{B K}^{c}},

(4)

where

η_{B K}^{C}

is the BK power exponent,

G_{n}^{C}

and

G_{s}^{C}

are the normal and shear fracture energies, respectively.

G_{n}

and

G_{s}

are the reciprocal work under mixed mode damage propagation.

Table 1.

Material properties of the constituents.²⁶

IM7 fibre properties
E₁ (GPa)	E₂ = E₃ (GPa)	$υ_{12} = υ_{13}$	$υ_{23}$	G₁₂ = G₁₃ (GPa)	G₂₃ (GPa)
287	13.4	0.29	0.48	23.8	7
8552 epoxy properties
E (GPa)	$υ_{m}$	$σ_{m y t}$ (MPa)	$σ_{m y c}$ (MPa)	G_m (J/m²)
4.08	0.38	99	130	100
Interface properties
Thickness (mm)	K_nn (GPa/mm)	K_ss = K_tt (GPa/mm)	$t_{n}^{0}$ (MPa)	$t_{s}^{0} = t_{t}^{0}$ (MPa)	G_IC (J/m²)	G_IIC = G_IIIC (J/m²)
0.0001	253	682	58	92	2	100

The calibrated interface strengths in normal and shear directions from our previous study²⁶ were 58 MPa and 92 MPa, respectively. These are in good agreement with the assumption that the interface normal strength is 2/3 of shear strength according to numerical¹⁹ and experimental studies.³⁶ The interface fracture energy in mode I, G_IC, could not be measured experimentally so it is usually assumed to be small and is in the range of 2–5 J/m².¹⁹ In this study, the fracture energy of 2 J/m² in Mode I is adopted in the simulations, and due to the lack of experimental data, the interface fracture energy in shear is assumed to be equal to the matrix cracking fracture energy, 100 J/m².¹⁹ The friction between fibres and matrix after interface failure was considered and implemented in an ABAQUS/Explicit subroutine.^27,37,38 Figure 2(b) shows the constitutive model of the interface under pure shear (black line) and combined compression and shear loads (blue line). The material properties can be found in Table 1.

Periodic boundary conditions and loading cases

Periodic boundary conditions (PBCs) are imposed on the corresponding surfaces of the RVE using linear equations between the periodic nodes at opposite faces to guarantee the periodicity of the displacement and traction. The detailed implementation of PBCs can be found in.³⁹ The strains were computed from the imposed displacement divided by the corresponding lengths, while the normal and shear stresses were computed from the resultant normal and tangential forces acting on the RVE faces divided by the cross-section area. Six loading conditions out of the full complement of 10 are necessary for the failure analysis under biaxial loading conditions due to the transverse isotropy in the cross-section of unidirectional composites. In a loading condition with two stress components, four quadrants are used to describe the biaxial stress space. The calculated stress ratios are based on the resultant stresses obtained from numerical simulations under biaxial loading conditions. It was assumed that the failure envelopes under most of biaxial loadings are polynomial curves, which need a minimum of three points. However, for the loading scenario of transverse compression and in-plane shear, where more complex failure mechanisms are involved (i.e. friction between fibre and matrix), more than five loading cases are required to capture representative failure modes and to provide enough data for the curve fitting of failure envelopes. Figure 3 shows the total number of biaxial loading conditions and an example of the stress ratio determination under combined transverse tension and in-plane shear. In Figure 3(a), the grid represents a biaxial loading condition, and the same colour in different grid locations represents an equivalent loading condition. Here two sets of loading conditions are grouped: transverse dominated loading conditions, including ( $σ_{2}, σ_{3}$ ), ( $σ_{2}, τ_{12}$ ), ( $σ_{2}, τ_{13}$ ), ( $σ_{2}, τ_{23}$ ); and shear dominated loading conditions, including ( $τ_{12}, τ_{13}$ ) and ( $τ_{12}, τ_{23}$ ). Here seven points are considered in the combined transverse tension and in-plane shear (see Figure 3(b)), including one point on each axis, representing uniaxial failure strength, and five points in a quadrant, representing five different biaxial stress ratios at failure. Since displacement load was imposed on the RVE, the resultant reaction forces were obtained to calculate the stresses. Five initial displacement ratios were selected, namely load cases 1–5 (see Figure 3(c)), which represent $\frac{δ_{12}}{δ_{2}} =$ 0.3, 0.6, 1.2, 2.5, 5 for transverse tension and in-plane shear loading cases for the calculation of failure stresses in the biaxial loading conditions. The failure of the whole RVE is determined based on the biaxial stress curve, when either stress starts to decrease. Take the combined transverse tension and in-plane shear case as an example, see Figure 3(c). Five biaxial loading cases were selected and biaxial stress curves were plotted. Both stresses increase linearly from zero, then the failure of the RVE is determined when transverse tensile stress starts to decrease, where the red cross is marked on the curve. The failure stress ratios can be calculated at the failure points accordingly as $\frac{τ_{12}}{σ_{2}} =$ 0.24, 0.49, 0.9, 1.5, 3.5.

Figure 3.

(a) Biaxial loading conditions, (b) An example of combined transverse tension and in-plane shear loading and (c) Stress curves in the ( $+ σ_{2}$ , ${+ τ}_{12})$ quradrant. (Red cross represents the failure point under biaxial loadings and the stress ratio besides the failure point represents the resultant failure stress ratio.)

Determination of the transition point in biaxial loading cases

It is challenging to determine the transition point on the failure envelopes, especially when there is no obvious change on the curve. Two main strategies were used for the determination of the transition points. The first one is based on the change of the failure envelopes. The Puck theory, which is used to determine the transition point, takes two loading conditions into consideration, namely transverse compression and in-plane shear, transverse compression and out-of-plane shear. The failure envelopes obtained from both loading conditions have an obvious peak point, which is used for the determination of transition point. While for other loading conditions, which cannot be performed experimentally and have no corresponding theory for the prediction of transition points, the second strategy was used, which is based on the change of failure modes. For example, the failure mode under pure out-of-plane shear is interface debonding and matrix cracking, while under pure in-plane shear is matrix yielding. When these damage modes coexist in the RVE, the transition point is determined when there is no dominant or less obvious mode at final failure. It should be noted that the second strategy is microstructure-depdendant since failure modes with different microstructure could result in slightly different transition points.

Conventional failure theories

Tsai-Wu failure criterion

The Tsai-Wu failure criterion⁵ is a phenomenological interactive failure criterion for anisotropic composite materials. It is a highly interactive equation, which integrates all stress components in one equation. The general expression of the Tsai-Wu failure criterion is

F_{i} σ_{i} + F_{i j} σ_{i j} = 1, (i, j = 1, 2, 3, 4, 5, 6)

(5)

where

F_{i}, F_{i j}

are the coefficients associated with the material strengths determined using experiments. Considering the random distribution of fibres within the matrix at the cross-section of most UD composites, transverse isotropy is usually assumed to describe the mechanical behaviour of UD composites. Consequently, the Tsai-Wu failure criterion for inter-fibre failure under 3D stress states reads with the fibre direction of 1

F_{2} (σ_{2} + σ_{3}) + F_{22} (σ_{2}^{2} + σ_{3}^{2}) + F_{66} (τ_{12}^{2} + τ_{13}^{2}) + F_{44} τ_{23}^{2} + (2 F_{22} - F_{44}) σ_{2} σ_{3} = 1

(6)

where the coefficients can be calculated as shown below

F_{2} = \frac{1}{Y_{T}} - \frac{1}{Y_{C}}, F_{22} = \frac{1}{Y_{T} Y_{C}}, F_{44} = \frac{1}{S_{23}^{2}}, F_{66} = \frac{1}{S_{12}^{2}}

(7)

with

X_{T}

and

X_{C}

being the tensile and compressive strengths of the composite material along the fibres,

Y_{T}

and

Y_{C}

the tensile and compressive strengths of the composite material in the transverse direction, and

S_{12}

and

S_{23}

the shear strengths along and transverse to the fibres, respectively. In this study, the strengths of IM7/8552 UD composite were obtained from standard experiments under quasi-static uniaxial or pure shear stress states.⁴⁰

Hashin failure criteria

Hashin proposed separate failure criteria for fibre and matrix by assuming a quadratic interaction between the tractions on the failure plane. These criteria can distinguish the failure of fibre and matrix in tension and compression⁷:

Matrix tension $(σ_{2} + σ_{3} \geq 0)$

{(\frac{σ_{2} + σ_{3}}{Y_{T}})}^{2} + \frac{1}{S_{23}^{2}} (τ_{23}^{2} - σ_{2} σ_{3}) + \frac{τ_{12}^{2} + τ_{13}^{2}}{S_{12}^{2}} = 1

(8)

Matrix compression $(σ_{2} + σ_{3} < 0)$

{(\frac{σ_{2} + σ_{3}}{{2 S}_{23}})}^{2} + [{(\frac{Y_{C}}{{2 S}_{23}})}^{2} - 1] \frac{σ_{2} + σ_{3}}{Y_{C}} + \frac{1}{S_{23}^{2}} (τ_{23}^{2} - σ_{2} σ_{3}) + \frac{τ_{12}^{2} + τ_{13}^{2}}{S_{12}^{2}} = 1

(9)

Computational results and discussions

Transverse dominated loadings

Failure envelope comparison

Figure 4 presents the comparison between the numerical results of the UD composite and Tsai-Wu and Hashin failure criteria under biaxial transverse dominated loadings. It can be found in Figure 4(a) (orange square box) that both Tsai-Wu and Hashin failure criteria underestimate the failure strength of the composite under biaxial transverse and out-of-plane tension loadings. Failure strengths predicted by Tsai-Wu and Hashin failure criteria agree with numerical results from the loading case where the stress ratio is −0.28. The comparison suggests that the predicted envelope from the Hashin criterion agrees well with the fitted envelope from numerical failure points to the infinite strength. Due to its closed envelope, the prediction from the Tsai-Wu criterion only agrees with the fitted envelope until the loading condition with a stress ratio of 3.13. It is still an open question whether the failure surface should be open or closed, especially under hydrostatic pressure.¹³ According to the numerical results, the fitted envelope under transverse compression and out-of-plane compression is open since the matrix suffers from hydrostatic pressure due to the biaxial compressive loads and constraints from fibres, resulting in a higher failure strength of the matrix. The comparison between the numerical results and the failure envelopes predicted by the conventional failure criteria under biaxial transverse compression and out-of-plane shear (see Figure 4(b)) shows excellent agreement. However, both failure criteria overestimate the failure strengths when transverse tension was involved. This is mainly because under the transverse tension and out-of-plane shear loadings, the interface failure dominates which is not considered in these conventional failure criteria. Figure 4(c) shows the comparison between the experimental results,⁴¹ numerical results and failure envelopes predicted by the failure criteria under combined transverse and in-plane shear loadings ( $σ_{2}, τ_{12}$ ). Scattered failure points are observed in numerical simulations from five RVEs and experiments, and the increase of shear strength under moderate transverse compression due to friction between fibres and matrix was captured by the numerical simulations via a VUMAT subroutine.²⁷ It can be found that both failure criteria agree reasonably well with experimental data and numerical results except for the part where the transverse compressive stress is around −200 MPa. This is mainly due to the limitation of quadratic mathematical formulation in the Tsai-Wu failure criteria and the inability to determine the actual fracture plane under transverse compression in the Hashin failure criteria. In the comparison of numerical results and the envelopes plotted from failure criteria under combined transverse and in-plane shear ( $σ_{2}, τ_{13}$ ) in Figure 8(d), overestimation and underestimation of failure strength predicted by both failure criteria are found in ( ${- σ}_{2}, τ_{13}$ ) and ( ${+ σ}_{2}, τ_{13}$ ) quadrants, respectively.

Figure 4.

Comparison between failure envelopes of unidirectional composites fitted from numerical failure points, Tsai-Wu and Hashin failure criteria under transverse dominated loadings: (a) transverse and transverse stresses, (b) transverse and out-of-plane shear stresses, (c) transverse and in-plane shear stresses (same loading plane) and (d) transverse and in-plane shear stresses (different loading plane).

Failure modes

Figure 5 shows the comparison of failure modes obtained from different stress ratios under biaxial transverse and out-of-plane tensile/compressive loadings. These stress ratios correspond to the ones labelled in Figure 4(a). It can be found that in the loading of ( ${+ σ}_{2}, + σ_{3}$ ), the failure modes are fibre/matrix interface debonding and matrix tensile failure, and the fracture plane which is prominent at a load ratio of 0.3 becomes less apparent at a ratio of 0.6 until it disappears at a ratio of 1.06 due to the similar value of transverse tensile stress and out-of-plane tensile stress. In the loading regime of ( ${+ σ}_{2}, - σ_{3}$ ), the failure mode of the composite transferred from interface debonding and matrix tensile failure to matrix shear failure, by out-of-plane compression, at a ratio of −1.83, where these failure modes coexisted. The matrix shear failure under out-of-plane compression and the fracture planes are not detected clearly in the ratio of −3.2. The same phenomenon was found in the loading regime of ( ${- σ}_{2}, - σ_{3}$ ) with a ratio of 6.95, while beyond this point, no reduction was detected in both stress curves, indicating no failure was manifested under such loading conditions. This can also be shown in the equivalent plastic strain distribution which becomes random and severe due to the hydrostatic pressure on the matrix.

Figure 5.

Failure modes of the IM7/8552 composite under biaxial transverse and out-of-plane tensile/compressive loadings. (For interpretation purposes, The matrix tensile failure ‘T’ and shear failure ‘SC’ are characterised by the damage variables DAMAGET and PEEQ, respectively. ‘No’ means no failure found in the loading case).

Figure 6 shows the comparison of failure modes obtained from different stress ratios under biaxial transverse tension/compression and out-of-plane shear loadings, and these stress ratios correspond to the ones labelled in Figure 4(b). Under biaxial transverse tension and out-of-plane shear, due to the smaller normal strength of the fibre/matrix interface compared to the tensile strength of the matrix, interface debonding occurs first, followed by matrix cracking. Both damages join together to form a fracture plane either perpendicular to the transverse direction (for transverse tension-dominated loadings) or inclined at an angle of 40° for out-of-plane shear-dominated loadings. Since there is no failure mode transition for this loading case, no transition point can be obtained. For the sake of the discussion on the progressive failure at different stress ratios, three ratios were still used. At the stress ratio $\frac{τ_{23}}{σ_{2}} =$ 0.53, the fibre/matrix interface experiences mixed-mode stresses under transverse tension and out-of-plane shear, initiating faster damage at the interface. Under combined transverse compression and out-of-plane shear loadings, the fracture plane, with an angle of 50° with respect to the plane perpendicular to the compression axis, disappears at the transition point and severe tensile damage is observed under out-of-plane shear-dominated loads.

Figure 6.

Failure modes of the IM7/8552 composite under combined transverse and out-of-plane shear. (For interpretation purposes, The matrix tensile failure ‘T’ and shear failure ‘SC’ are characterised by the damage variables DAMAGET and PEEQ, respectively).

Figures 7 and 8 show the comparison of failure modes obtained from different stress ratios under biaxial transverse and in-plane shear loadings. These stress ratios correspond to the ones labelled in Figures 4(c) and (d), respectively. Similar failure modes can be found under ( $σ_{2}, τ_{12}$ ) and ( $σ_{2}, τ_{13}$ ), though with different stress ratios. Under the combined transverse tension and in-plane shear with $\frac{τ_{12}}{σ_{2}} =$ 1.5, one main shear fracture plane is observed close to the left edge, while the shear fracture plane becomes random when $\frac{τ_{13}}{σ_{2}} =$ 2.43. The difference comes from the fact that the fibre/matrix interface experiences mixed-mode stress from normal and shear directions simultaneously when the RVE1 is subjected to ( $+ σ_{2}, τ_{12}$ ), resulting in its earlier failure which triggers the matrix shear failure; while under ( $+ σ_{2}, τ_{13}$ ), the interface suffers single-mode stress, either in the normal or shear direction, which delays matrix failure. Under combined transverse compression and in-plane shear loadings with dominated compressive stress, a shear band with a fracture angle of 50° is predicted in both loadings. A remarkably similar trend is suggested by the computational micromechanical results¹⁹ with a fracture angle of 56° and by experimental observations with an angle of 53°,⁴⁰ with some differences. These differences are most likely related to the discreteness of the microstructure. At the transition point under combined transverse compression and in-plane shear, both fracture planes with angles of 0° and 53° disappear during the competition of compressive and shear stresses. While beyond the transition point, clear shear bands parallel to the shear directions ( $τ_{12}$ and $τ_{13}$ ) are formed close to the edges.

Figure 7.

Failure modes of the IM7/8552 composite under combined transverse and in-plane shear loadings ( $σ_{2}, τ_{12}$ ). (For interpretation purposes, The matrix tensile failure ‘T’ is characterised by the damage variable DAMAGET. The matrix shear failure under compression ‘SC’ and shear ‘SS’ are characterised by the equivalent plastic strain PEEQ).

Figure 8.

Failure modes of the IM7/8552 composite under combined transverse and in-plane shear loadings $(σ_{2}, τ_{13})$ . (For interpretation purposes, The matrix tensile failure ‘T’ is characterised by the damage variable DAMAGET. The matrix shear failure under compression ‘SC’ and shear ‘SS’ are characterised by the equivalent plastic strain PEEQ).

Shear dominated loadings

Failure envelope comparison

Figure 9 compares the numerical results and the failure envelopes plotted from conventional Hashin and Tsai-Wu failure criteria. The red squares in Figure 9 represent the numerical data for curve-fitting reasons, in which a negative $τ_{12}$ means the opposite direction of its counterpart. It is found in Figure 9(a) that the fitted curve of the failure envelope under biaxial in-plane shear loadings is in excellent agreement with the envelopes plotted with Tsai-Wu and Hashin failure criteria. Comparing the failure points obtained from (+ $τ_{12}, {+ τ}_{13}$ ) and (+ $τ_{12}, {- τ}_{13}$ ), a slight difference was detected regarding the failure points. Different stress ratios in loading cases of (+ $τ_{12}, {+ τ}_{13}$ ) and (+ $τ_{12}, {- τ}_{13}$ ) are most likely due to the microstructure of the selected RVE1, which changes the damage propagation path. Figure 9(b) compares the envelopes fitted from numerical results and plotted from Tsai-Wu and Hashin failure criteria under biaxial in-plane and out-of-plane shear loadings. Excellent agreement is found between these envelopes, which verifies the capability of failure prediction of the unidirectional composite lamina with both failure criteria under biaxial shear stresses.

Figure 9.

Comparison between failure envelopes of unidirectional composites from fitted numerical failure points, Tsai-Wu and Hashin failure criteria under shear-dominated loadings: (a) in-plane shear and in-plane shear stresses, (b) in-plane shear and out-of-plane shear stresses.

Failure modes

Figure 10 shows the comparison of failure modes of the UD composite under biaxial in-plane shear loadings with different directions. It can be found that when the shear stresses are in the same direction and similar, only one main shear plane with an inclination angle of 45° is formed, while two or more shear planes with angles less than 45° are formed when either shear stress is larger than the other. However, several shear planes are formed when the shear stresses are in the opposite direction and similar. The number of planes decreases when either shear stress is larger than the other and two main planes with an inclination angle of 45° can be observed. Figure 11 shows the comparison of failure modes under biaxial in-plane and out-of-plane shear loadings with different directions. The same failure modes and a similar trend of shear bands can be found in both loading conditions. The difference in the fracture angle in both out-of-plane stresses dominated loadings probably comes from the microstructure of the cross-section in composites, which changes the direction of the fracture propagation.

Figure 10.

Failure modes of the IM7/8552 composite under biaxial in-plane shear loadings ( $τ_{12}, τ_{13}$ ).

Figure 11.

Failure modes of the IM7/8552 composite under biaxial in-plane and out-of-plane shear loadings ( $τ_{12}, τ_{23}$ ). (For interpretation purposes, The matrix tensile failure ‘T’ is characterised by the damage variable DAMAGET. The matrix shear failure ‘SS’ is characterised by the equivalent plastic strain PEEQ).

Transition point determination under biaxial loading conditions

The transition point plays an important role in distinguishing failure mode under biaxial loadings, which can be determined analytically⁶ and experimentally.⁴² While no failure transition point is detected if there is only one failure mode existing under biaxial loadings, such as the tensile failure mode (‘T’) in biaxial transverse tension in Figure 5. Table 2 collects transition points of UD IM7/8552 composites from computational analysis under biaxial loading conditions, which are grouped in three sets and compared with numerical and experimental results from the literature.

Table 2.

Stress ratio of transition points of UD IM7/8552 composites under biaxial loadings.

In the transverse dominated biaxial loadings, since the failure modes of the composite are interface and matrix tensile failure under transverse tension and out-of-plane tension/shear loads, there is no failure mode transition point in such biaxial loadings. While under combined transverse tension and out-of-plane compression loads, only matrix shear failure under compression was found and no failure was found under biaxial transverse and out-of-plane compression. Under transverse compression and out-of-plane shear, good agreement in the stress ratio of transition points is found between the one obtained from the IM7/8552 composites and glass FRP composites.^9,43 Off-axial tests have been conducted on different material systems to determine the failure strengths under biaxial loadings, such as4/3501–6,⁴⁴ IM7/8552⁸ and E-glass/RP528.⁴⁵ An excellent agreement can be found between experimental results, even across different material systems. In the biaxial in-plane shear loadings, only matrix shear failure due to shear stress is found thus no transition point is detected, while in the combined in-plane shear and out-of-plane shear loadings, the stress ratio ( ${τ_{23} / τ}_{12}$ ) of −0.77 is determined, beyond which the failure mode of the matrix transforms from shear dominated failure to tension dominated failure.

Based on Puck’s theory,⁶ the corresponding values at transition points under combined transverse and in-plane shear, and transverse and out-of-plane shear loadings can be calculated for comparison with numerical and experimental results

| σ_{2}^{'} | = \frac{Y_{c}}{2 (1 + p_{⊥ ⊥}^{(-)})}

(10)

| τ_{12}^{'} | = S_{12} \sqrt{1 + {2 p}_{⊥ ⊥}^{(-)}}

(11)

| τ_{23}^{'} | = S_{23} \sqrt{1 + {2 p}_{⊥ ⊥}^{(-)}}

(12)

p_{⊥ ⊥}^{(-)} = \frac{1}{2} \sqrt{1 + {2 p}_{⊥ ∥}^{(-)} \frac{Y_{c}}{S_{12 (23)}} - 1}

(13)

where the recommended range of

p_{⊥ ∥}^{(-)}

is 0.25–0.30 for carbon fiber reinforced composites.⁶ Here,

p_{⊥ ∥}^{(-)}

= 0.25 is selected for the calculation of

p_{⊥ ⊥}^{(-)}

according to equation (13).

Conclusion

This study introduced a 3D high-fidelity micromechanics-based RVE model with random fibre distributions to predict the inter-fibre failure envelopes and failure modes of IM7/8552 UD CFRP composite laminae under biaxial loadings. The RVE model was built with three constituents, namely the fibre, matrix and fibre/matrix interface. Fibres were modelled to be transversely isotropic and elastic; the interface was modelled with a cohesive zone model with the consideration of friction via a subroutine; the matrix was modelled by a Drucker-Prager plastic damage model. Only six out of the 10 possible biaxial stress combinations were taken into account due to the transversely isotropic characteristics of the cross-section and symmetry of shear stresses. Periodic boundary conditions were utilised for the application of biaxial loadings. The novelty of the current study lies in the use of (i) high-fidelity 3D micromechanics-based finite element modelling to assess failure criteria and (ii) the determination of transition points and progressive inter-fibre failure analysis of UD composties, under various possible biaxial loading conditions, some of which are not easy to determine experimentally.

Comprehensive biaxial inter-fibre failure analysis of UD IM7/8552 composites and the comparison of failure envelopes obtained from computational analysis and classical failure criteria (Tsai-Wu and Hashin) were conducted considering the failure mode transition. The main findings and conclusions are given below:

• Both failure criteria agree well with the computational results in $({+ σ}_{2}, {- σ}_{3})$ and $({- σ}_{2}, {- σ}_{3})$ while underestimating the failure strength in $({+ σ}_{2}, {+ σ}_{3}) .$ Both criteria generally agree well with the computational and experimental results in $(σ_{2}, τ_{12})$ but neither can capture the increase of shear strength under moderate compressive stress due to friction between fibre and matrix, which was predicted by the RVE model with the frictional-cohesive model. This phenomenon was not found in the numerical results under $(- σ_{2}, τ_{13}) .$ Thus this failure mechanism should be taken into account in future failure criteria modifications or new ones. Excellent agreement was observed between these criteria in $({- σ}_{2}, τ_{23})$ but overestimated the results in $({+ σ}_{2}, τ_{23})$ .

• Both failure criteria agree well with the computational results for shear-dominated loadings. The difference between the results obtained from $(τ_{12}, τ_{13})$ and $(τ_{12}, {- τ}_{13})$ is likely related to the discreetness of the microstructure at the cross-section.

• The stress ratio at the transition points describes the failure transition mechanisms which could be taken into account in the failure criteria for better prediction of composite failure under biaxial loading conditions.

These virtual tests can provide full control of the composite microstructural and material properties, allowing the optimisation and exploration of the microstructure in the improvement of the mechanical performance of composites. Also, this work suggests a need to improve existing mesoscale failure criteria that only rely on ply properties.^4,6 Some important microstructural parameters and/or information are suggested to be taken into account in the existing or newly proposed failure criteria, such as fibre/matrix interfacial debonding,^21,24 fibre/matrix friction¹⁹ and failure transition points.^1,8 Furthermore, 3D high-fidelity RVE modelling of composites under longitudinal loadings is under investigation, which considers fibre kinking and stochastic fibre strength phenomena. The failure points obtained from these simulations can be used to construct fibre dominated failure criteria.

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 study was conducted as part of the Belfast Maritime Consortium UKRI Strength in Places project, “Decarbonisation of Maritime Transportation: A return to Commercial Sailing” led by artemis Technologies, Project no. 107138.

ORCID iDs

Lei Wan

Zahur Ullah

Dongmin Yang

References

Sun

Zhou

Meng

, et al. Failure criteria of unidirectional carbon fiber reinforced polymer composites informed by a computational micromechanics model. Compos Sci Technol 2019; 172: 81–95. DOI: 10.1016/j.compscitech.2019.01.012

Christensen

. Tensor ransformations and ailure riteria for the nalysis of iber omposite aterials. J Compos Mater 1988; 22: 874–897. DOI: 10.1177/002199838802200906

Hart-Smith

. Predictions of the original and truncated maximum-strain failure models for certain fibrous composite laminates. Compos Sci Technol 1998; 58: 1151–1178. DOI: 10.1016/S0266-3538(97)00192-9

Dávila

Camanho

Rose

. Failure criteria for FRP laminates. J Compos Mater 2005; 39: 323–345. DOI: 10.1177/0021998305046452

Tsai

. A eneral heory of trength for nisotropic aterials. J Compos Mater 1971; 5: 58–80.

Puck

Schürmann

. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 2002; 62: 1633–1662. DOI: 10.1016/S0266-3538(01)00208-1

Hashin

. Failure riteria for nidirectional iber omposites. J Appl Mech 1980; 47: 329–334. DOI: 10.1115/1.3153664

Daniel

Fenner

. A new yield and failure theory for composite materials under static and dynamic loading. Int J Solids Struct 2018; 148–149: 79–93. DOI: 10.1016/j.ijsolstr.2017.08.036

Camanho

Arteiro

Melro

, et al. Three-dimensional invariant-based failure criteria for fibre-reinforced composites. Int J Solids Struct 2015; 55: 92–107. DOI: 10.1016/J.IJSOLSTR.2014.03.038

10.

Daniel

Luo

J-J

Schubel

, et al. Interfiber/interlaminar failure of composites under multi-axial states of stress. Compos Sci Technol 2009; 69: 764–771. DOI: 10.1016/j.compscitech.2008.04.016

11.

Hinton

Kaddour

Soden

. A further assessment of the predictive capabilities of current failure theories for composite laminates: comparison with experimental evidence. Compos Sci Technol 2004; 64: 549–588. DOI: 10.1016/S0266-3538(03)00227-6

12.

Soden

Hinton

Kaddour

. Biaxial test results for strength and deformation of a range of E-glass and carbon fibre reinforced composite laminates. Fail Criteria Fibre-Reinforced-Polymer Compos 2004; 62: 52–96. DOI: 10.1016/B978-008044475-8/50004-4

13.

Kaddour

Hinton

. Maturity of 3D failure criteria for fibre-reinforced composites: omparison between theories and experiments: Part B of WWFE-II. J Compos Mater 2013; 47: 925–966. DOI: 10.1177/0021998313478710

14.

Kaddour

Hinton

Smith

, et al. Mechanical properties and details of composite laminates for the test cases used in the third world-wide failure exercise. J Compos Mater 2013; 47: 2427–2442. DOI: 10.1177/0021998313499477

15.

Abir

Tay

Ridha

, et al. Modelling damage growth in composites subjected to impact and compression after impact. Compos Struct 2017; 168: 13–25. DOI: 10.1016/j.compstruct.2017.02.018

16.

Zhang

. An efficient approach for predicting low-velocity impact force and damage in composite laminates. Compos Struct 2015; 130: 85–94. DOI: 10.1016/j.compstruct.2015.04.023

17.

González

Lorca

. Mechanical behavior of unidirectional fiber-reinforced polymers under transverse compression: icroscopic mechanisms and modeling. Compos Sci Technol 2007; 67: 2795–2806. DOI: 10.1016/j.compscitech.2007.02.001

18.

Canal

Segurado

Lorca

. Failure surface of epoxy-modified fiber-reinforced composites under transverse tension and out-of-plane shear. Int J Solids Struct 2009; 46: 2265–2274. DOI: 10.1016/j.ijsolstr.2009.01.014

19.

Naya

González

Lopes

, et al. Computational micromechanics of the transverse and shear behavior of unidirectional fiber reinforced polymers including environmental effects. Compos Part A Appl Sci Manuf 2017; 92: 146–157. DOI: 10.1016/j.compositesa.2016.06.018

20.

Sebaey

Catalanotti

Lopes

, et al. Computational micromechanics of the effect of fibre misalignment on the longitudinal compression and shear properties of UD fibre-reinforced plastics. Compos Struct 2020; 248: 112487. DOI: 10.1016/j.compstruct.2020.112487

21.

Wan

Ismail

Zhu

, et al. Computational micromechanics-based prediction of the failure of unidirectional composite lamina subjected to transverse and in-plane shear stress states. J Compos Mater 2020; 54: 3637–3654. DOI: 10.1177/0021998320918015

22.

Sharma

Daggumati

Gupta

, et al. On the prediction of the bi-axial failure envelope of a UD CFRP composite lamina using computational micromechanics: ffect of microscale parameters on macroscale stress–strain behavior. Compos Struct 2020; 251: 112605. DOI: 10.1016/j.compstruct.2020.112605

23.

Elnekhaily

Talreja

. Effect of axial shear and transverse tension on early failure events in unidirectional polymer matrix composites. Compos Part A Appl Sci Manuf 2019; 119: 275–282. DOI: 10.1016/j.compositesa.2019.01.031

24.

Totry

González

Lorca

. Prediction of the failure locus of C/PEEK composites under transverse compression and longitudinal shear through computational micromechanics. Compos Sci Technol 2008; 68: 3128–3136. DOI: 10.1016/j.compscitech.2008.07.011

25.

Romanowicz

. A numerical approach for predicting the failure locus of fiber reinforced composites under combined transverse compression and axial tension. Comput Mater Sci 2012; 51: 7–12. DOI: 10.1016/j.commatsci.2011.07.039

26.

Chen

Wan

Ismail

, et al. A micromechanics and machine learning coupled approach for failure prediction of unidirectional CFRP composites under triaxial loading: preliminary study. Compos Struct 2021; 267: 113876. DOI: 10.1016/j.compstruct.2021.113876

27.

Wan

Ullah

Yang

, et al. Probability embedded failure prediction of unidirectional composites under biaxial loadings combining machine learning and micromechanical modelling. Compos Struct 2023; 312: 116837.

28.

Ismail

Yang

. Discrete element method for generating random fibre distributions in micromechanical models of fibre reinforced composite laminates. Compos Part B Eng 2016; 90: 485–492. DOI: 10.1016/j.compositesb.2016.01.037

29.

Ward

. Review: The yield behaviour of polymers. J Mater Sci 1971; 6: 1397–1417. DOI: 10.1007/BF00549685

30.

Fiedler

Hojo

Ochiai

, et al. Failure behavior of an epoxy matrix under different kinds of static loading. Compos Sci Technol 2001; 61: 1615–1624. DOI: 10.1016/S0266-3538(01)00057-4

31.

Wan

Ismail

Sheng

, et al. A review on micromechanical modelling of progressive failure in unidirectional fibre-reinforced composites. Compos Part C Open Access 2023; 10: 100348. DOI: 10.1016/j.jcomc.2023.100348

32.

Melro

Camanho

Andrade Pires

, et al. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part I – onstitutive modelling. Int J Solids Struct 2013; 50: 1897–1905. DOI: 10.1016/j.ijsolstr.2013.02.009

33.

Dassault Systèmes . Abaqus 6.13 documentation. 2013.

34.

Ismail

Wan

Chen

, et al. An ABAQUS® plug-in for generating virtual data required for inverse analysis of unidirectional composites using artificial neural networks. Eng Comput 2022; 38: 4323–4335.

35.

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. DOI: 10.1016/0266-3538(96)00005-X

36.

Ogihara

Koyanagi

. Investigation of combined stress state failure criterion for glass fiber/epoxy interface by the cruciform specimen method. Compos Sci Technol 2010; 70: 143–150. DOI: 10.1016/j.compscitech.2009.10.002

37.

Alfano

Sacco

. Combining interface damage and friction in a cohesive-zone model. Int J Numer Methods Eng 2006; 68: 542–582. DOI: 10.1002/nme.1728

38.

Turon

González

Sarrado

, et al. Accurate simulation of delamination under mixed-mode loading using a cohesive model with a mode-dependent penalty stiffness. Compos Struct 2018; 184: 506–511. DOI: 10.1016/j.compstruct.2017.10.017

39.

Garoz

Gilabert

Sevenois

RDB

, et al. Consistent application of periodic boundary conditions in implicit and explicit finite element simulations of damage in composites. Compos Part B Eng 2019; 168: 254–266. DOI: 10.1016/j.compositesb.2018.12.023

40.

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: 1004–1019.

41.

Körber

. Mechanical response of advanced composites under high strain rates. Portugal: Universidade do Porto, 2010.

42.

Vogler

Kyriakides

. Inelastic behavior of an AS4/PEEK composite under combined transverse compression and shear. Part I: experiments. Int J Plast 1999; 15: 783–806. DOI: 10.1016/S0749-6419(99)00011-X

43.

Ashouri Vajari

González

Llorca

, et al. A numerical study of the influence of microvoids in the transverse mechanical response of unidirectional composites. Compos Sci Technol 2014; 97: 46–54. DOI: 10.1016/j.compscitech.2014.04.004

44.

Daniel

Werner

Fenner

. Strain-rate-dependent failure criteria for composites. Compos Sci Technol 2011; 71: 357–364. DOI: 10.1016/j.compscitech.2010.11.028

45.

Gan

Laux

Taher

, et al. A novel fixture for determining the tension/compression-shear failure envelope of multidirectional composite laminates. Compos Struct 2018; 184: 662–673. DOI: 10.1016/j.compstruct.2017.10.030