Computational micromechanics-based prediction of the failure of unidirectional composite lamina subjected to transverse and in-plane shear stress states

Abstract

This paper presents a micromechanics-based 3D finite element model for predicting the damage initiation, propagation, and failure strength of TC33/Epoxy carbon fiber reinforced polymer (CFRP) unidirectional lamina under biaxial loadings. The finite element model is generated by introducing representative volume element (RVE) with a random distribution of fibers and a non-zero thickness, numerically identified interface phase via cohesive elements. In the finite element model, the carbon fibers are considered as elastic, while the elasto-plastic behavior and damage of the matrix are governed by extended Drucker–Prager plastic yielding model and ductile damage criterion. By imposing periodic boundary conditions to the RVEs, various cases subjected to uniaxial and biaxial loading conditions are carried out. During the combined transverse and in-plane shear stress states, a failure transition from compression- or tension-dominated to shear-dominated is captured, and the effects of the interfacial strength on the transition damage mechanisms are discussed. The corresponding failure locus is compared with the upper bound and lower bound predictions of three phenomenological failure criteria (Hashin, Tasi–Wu, and Puck failure criteria) for composites. It was found that in the interface-dominated failure of a CFRP lamina with a weak interface, the Hashin failure criterion performs best among the currently popular failure criteria. However, in the matrix-dominated failure with a strong interface, the Puck failure criterion performs best. Comparing these three criteria, it can be seen the Tsai–Wu may be generally better than both of others as it presents more neutral predictions in both of the examined cases.

Keywords

Carbon fiber reinforced polymer micromechanics finite element method representative volume element failure envelopes damage progression

Introduction

Carbon fiber reinforced polymer (CFRP) composites have been widely used in aerospace and mechanical industries for several decades due to their outstanding specific stiffness and strength. Moreover, the composite laminates with desired modulus and strength in different directions can be achieved by proper design and optimization of an individual lamina. However, in most cases, composite components are usually over-designed considering the safety factors, resulting in the larger and heavier parts. This is mainly due to lack of reliable predictive models and theories to fully understand the mechanical response and failure/damage mechanisms of the composites under various loading conditions.¹ Unlike most homogeneous materials, CFRP composite materials, with heterogeneous nature between different phases, tend to present multiple damage modes depending on the loading conditions, stress states, and possible manufacturing defects.² The coexistence of these various damage modes in CFRP composites implies the necessity of the combination of different failure criteria depending on the loading mode. Besides, a failure locus is usually formed by the intersection of various smooth surfaces in the stress space, and each one represents the critical condition for a given fracture.³ Therefore, it is still challenging to accurately predict failure envelopes of CFRP composites subjected to multi-axial stress states.

Over the last several decades, a large number of failure criteria have been proposed to predict the failure envelopes of the composite materials from limited experimental data, and some of them have been successfully applied to design and calculate the safety factors for the composite structures. They include but not limited to strain-based,⁴ strain-energy-based,⁵ stress-based,^6–9 and phenomenological failure criteria.^10–15 Although many failure criteria have been proposed and the modification of some criteria is still ongoing, the validation of these criteria remains challenging. That is due to the inability of these criteria regarding the accurate prediction of the progressive failure process in composites. An objective assessment of the currently available failure criteria for fiber-reinforced composites was conducted under complex 3D stress state to predict the failure strength and to describe the failure envelop in a series of three World Wide Failure Exercises.^16–22 The Second World Wide Failure Exercise (WWFE-II), of assessing some existing failure criteria for FRP composite laminates, has shown the satisfactory performance of each criterion to various degrees. However, there still exist considerable variations in the accuracy of the predictions by these criteria. More recently, the third World Wide Failure Exercise (WWFE-III) was conducted to highlight the degree of maturity of 12 internationally recognized approaches (some of them are different from the criteria mentioned in WWFE-II) considering their capabilities of detecting the various damages within the composite materials when subjected to multi-axial loading.^20–22 It was found that any two models cannot give identical predictions for any of the 13 test cases. In a few cases, the ratio between the highest and lowest strength predictions can even reach a factor of 20. In addition, the monitoring and visualization of in-situ damage progression during mechanical tests is no doubt challenging and expensive, especially for the multi-axial stress states. Therefore, precise conclusions have not been reached regarding which criterion can best reproduce the physical failure mechanisms and the mechanical strength because of the scarcity of the experimental data, especially under multi-axial stress states. Thus, many criteria have still not been validated for the prediction of the strength of the fiber-reinforced composites. Meanwhile, the input parameters of these aforementioned models were obtained through costly and time-consuming experiments for different material system. However, the results obtained from a given unidirectional (UD) composite material system cannot be extrapolated to other configurations with different fiber volume fraction or constituent properties, leading to a huge amount of investment in their physical characterization.²³

Thanks to the increasing computation power, many of these difficulties can be overcome by taking advantage of computational micromechanics. Computational micromechanics offers a novel approach for a better understanding of the deformation and damage mechanisms by employing the representative volume element (RVE) modeling. Compared to the classic homogenization techniques, computational micromechanics presents two main advantages. On the one hand, this method takes the influences of the geometry and spatial distribution of the three phases into consideration. For example, the size of fiber, the fiber clustering, and the interface connectivity between fiber and matrix are included. On the other hand, the details of the stress and strain distribution under different loading conditions can be captured, leading to more accurate estimation of the onset and progressive process of damage, and the final failure strength.²⁴ Recently, computational analysis was successfully employed to investigate the mechanical response of the different material systems for the fiber-reinforced composite lamina subjected to the different combined loading conditions, such as transverse tension and out-of-plane shear,³ transverse compression and out-of-plane shear,²⁵ transverse compression and in-plane shear,^23,26 transverse tension and in-plane shear,²³ and transverse compression and axial tension.²⁷ Interface modeling is a crucial part of RVE modeling with the finite element method (FEM), and usually, a cohesive crack model is implemented to simulate the mechanical response of the interface between the fiber and matrix. In the linear behavior before the onset of damage, an initial stiffness K_i (10⁵GPa/mm) is used in most research^26–32 to simulate the elastic behavior of the RVE model. The choice of this parameter is based on that it should be large enough to ensure the displacement continuity at the interface and to avoid any modification of the stress fields around the fiber before damage.²⁶ However, it was found that the average Young’s modulus and the strength of the interphase are around five and nine times larger than those of the bulk resin matrix.³² Therefore, the interphase was modeled as a separate zone with the same constitutive and damage models of matrix, which makes the model more complicated. The parameter identification of the interphase of the carbon fiber reinforced composite was conducted by the inverse strategy based on the experimental data, microstructural modeling method, and Kriging metamodel, including the identification of thickness, normal stiffness, and tangent stiffness.³³ This set of parameters was applied to predict the elastic and strength properties of carbon FRP composite yarn, and good agreement was found with other regular simulation models.³⁴ In order to obtain accurate results from FE simulations, the constitutive model used for the matrix material modeling plays an important role. There are several constitutive models proposed for modeling matrix behavior, namely Mohr–Coulomb model (M–C), Drucker–Prager model (D–P), and a new model proposed by Melro et al.,³⁵ which are widely used to conduct failure analysis of composite materials within the framework of micromechanics, see Canal et al.,³ Gonźalez and LLorca,²⁴ Totry et al.,²⁶ Melro et al.,²⁹ Tan et al.,³¹ and Sun et al.³²

In this paper, a micromechanics-based FE model employing RVE modeling is developed for the investigation of the failure/damage mechanisms and assessment of the failure envelopes of UD CFRP composites subjected to transverse and in-plane shear stress states. A novel approach for generating random fiber distributions is applied using the discrete element method (DEM) with high volume fractions and any specified inter-fiber distances.³⁶ Five RVE models with different random fiber distributions and the same volume fraction are built for the investigation of the effects of the different fiber distributions on the mechanical response of the composite. The linear extended D–P model and the ductile damage model are adopted in the modeling of plastic behavior and damage of the matrix, the damage mode transition, and associated change of stress-based yielding law in the D–P model are discussed. In addition, the cohesive crack model is adopted for simulating the mechanical behavior of the interface with the identified elastic parameters from Lu et al. and Zhu et al.^33,34 The predictions with upper and lower bounds from three popular failure criteria are compared with the numerical results in terms of a weak and a strong interface for the evaluation of these criteria under different transverse and in-plane shear stress states.

Computational micromechanics modeling

Computational micromechanics modeling is performed on the RVEs of the UD CFRP subjected to homogeneous stress states, such as tension, compression, and shear. A total number of around 50 fibers is enough to capture adequately the essential features of the microstructure of the material²⁴ while maintaining reasonable computing efforts. In addition, the average radius value and average volume fraction are obtained experimentally to be 3.115 μm and 65.12% in Zhu et al.,³⁴ respectively. The microstructure of each RVE is idealized as the dispersion of 50 circular fibers with an average radius of 3.2 μm randomly embedded in the matrix, making a volume fraction of 65%. The size of RVEs is chosen to be 50 μm × 50 μm × 6 μm, considering the insignificant effects of depth on the transverse properties. The novel approach for generating the 2D random fiber distribution to overcome the jamming limit is adopted from the previous study,³⁶ and extruded along the fiber direction to achieve the final RVE configuration of the UD composite material. FEM models are generated in ABAQUS/Explicit to overcome the convergence difficulty of numerical analysis. The fibers and matrix are meshed using six-node linear triangular prism elements (C3D6) with hourglass control, and the interface is meshed with eight-node cohesive elements (COH3D8). Typically, around 52,000 elements are adopted in each RVE to capture the large stress gradients between neighboring fibers; see Figure 1 for an illustration of the RVE model with random fiber distribution. In order to accelerate the simulation process, mass scaling is normally utilized in the ABAQUS/Explicit to artificially increase the mass of elements, resulting in an increase of the time increment. However, this technique can influence the results largely, especially during an analysis of dynamic study where the inertia effects become dominant. A common way to check the influences of the mass scaling on the numerical results is to compare the kinetic energy with the internal energy, and a ratio below 10% can be regarded as an insignificant effect.^37,38 Therefore, the stable time increment for the mass scaling is selected as 6 × $10^{- 6}$ s to ensure a negligible influence on the results. The initial time step is around 5.37 × $10^{- 12}$ s, and the cost of simulation time is 6–8 h, depending on the loading conditions and the fiber spatial distribution within RVE, with the 8-core Intel E5-2670 (2.6 GHz) processors and 16 GB memory on ARC2 HPC operating system provided by the University of Leeds.

Figure 1.

3D RVE FE model with random fiber distribution and its microstructure distribution.

The response of the UD composite subjected to various loading conditions can be obtained from the homogenized stress σ_ij and the homogenized strain $ɛ_{ij}$ . While the simulated results obtained from the microscale model are only local stress and strain distributions within the RVE, the homogenized stress can be calculated at each timestep from the equation below

σ_{ij} = \frac{1}{V} \int_{V} (σ_{ij}^{loc}) d V

(1)

where

σ_{ij}^{loc}

is the local stress tensor and V is the volume of the RVE. Besides, within the frame of the first-order computational homogenization, the homogenized strain is equal to the macro-strain of the dummy node. It is worth noting that the nonlinear displacements within the RVE, such as the interface debonding when it fails, are taken into consideration within the macro-strain or homogenized strain. Meanwhile, it is important to note that the application of periodic boundary conditions (PBCs) introduces some constraints and limitations regarding the periodicity of the results, which means the damage created within the RVE represents a periodic damage instead of a local defect.³⁹

Constitutive laws of fiber and matrix

Carbon fibers are modeled as linear, elastic, and transversely isotropic solids, and their anisotropies are taken into consideration by defining five independent elastic constants, which are given in Table 1. While for the isotropic epoxy resin matrix in this study, it was assumed to behave as an elastic perfectly plastic solid. Since the behavior of polymers is found to be sensitive to the hydrostatic stress,⁴⁰ the extended linear D–P yield criterion is adopted in conjunction with the ductile damage criterion to simulate the plastic deformation and capture the damage process of the polymer matrix. The extended linear D–P criterion can be expressed as below⁴¹

F = t - p tan β - d = 0, t = \frac{1}{2} q [1 + \frac{1}{k} - (1 - \frac{1}{k}) (\frac{r}{q})^{3}]

(2)

where p is the equivalent pressure stress, β is the slope of the linear yield surface in the p–t plane (also referred to as the friction angle of the material in the D–P criterion), d is the cohesion of the material and related to the yielding stress when the yielding behavior is defined by the uniaxial compression, tension, and shear, q is the Mises equivalent stress, k is the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression and controls the difference in yielding behavior between tension and compression. According to the M–C model and the relation between M–C and D–P models, the internal friction angle β and cohesion of the matrix are

18 . 5^{\circ}

and 62.4 MPa, respectively. The matrix is assumed as perfect plasticity, and no hardening is considered in this study. The plastic flow of the material is controlled by the flow potential G, and is expressed as

G = t - p tan ψ

(3)

where ψ is the dilation angle in the p–t plane. Experimental findings⁴² showed that the use of the associated flow rule for polymers overestimates the extent of plastic dilatancy, thus a non-associated flow rule is utilized to compute the direction of the plastic flow.

Table 1.

Mechanical properties of the fiber and matrix.^23,34

Mechanical properties	TC33 fibers	Epoxy resin matrix
Longitudinal modulus E₁ (GPa)	230	3.08
Transverse modulus E₂ = E₃ (GPa)	106	–
In-plane shear modulus G₁₂ (GPa)	20.5	–
Major Poisson’s ratio $υ_{12}$	0.255	0.35
Transverse Poisson’s ratio $υ_{23}$	0.405	–
Tensile strength (MPa)	–	60.96
Compressive strength (MPa)	–	90.54
Critical fracture energy ( ${J / m}^{2}$ )	–	100

After the yielding process of the polymeric material, an additional criterion is needed for the prediction of the onset of damage. Experimental findings⁴³ show that the polymer exhibits a rather brittle fracture behavior under the uniaxial tension while a large plastic deformation under uniaxial compression and pure shear. This behavior can be governed by the ductile criterion, which assumes that the equivalent plastic strain at the onset of damage ${\bar{ɛ}}_{D}^{pl}$ is a function of stress triaxiality η and strain rate $\overset{\cdot}{\bar{}} ɛ^{pl}$ , where $η = - p / q$ . The criterion for damage initiation (ω_D) is met when the following condition is satisfied

ω_{D} = \int \frac{d {\bar{ɛ}}^{pl}}{{\bar{ɛ}}_{D}^{pl} (η, \overset{\cdot}{\bar{}} ɛ^{pl})} = 1

(4)

Once the ductile damage initiation criterion is met, the following damage propagation is controlled by a damage evolution law, and the whole process from elastic behavior until final failure is illustrated in Figure 2. The solid curve represents the damaged stress–strain response, which manifests in two different forms: softening of the bold yielding stress and degradation of the elasticity, while the dashed curve is the material response in the absence of damage. When material damage occurs, the stress–strain relationship can no longer represent the material behavior accurately. That is because the stress–strain relationship can introduce a strong mesh dependency based on strain localization, in which the energy dissipation decreases with a finer mesh. Therefore, the softening response after damage initiation is characterized by a stress-displacement response rather than a stress–strain response to alleviate mesh dependency of the results. This can be achieved in the FE model by introducing the critical fracture energy, defined as the energy required to open a unit area of crack, and a characteristic length L. The fracture energy can be expressed as

G_{f} = \int {\bar{ɛ}}_{0}^{pl} {\bar{ɛ}}_{f}^{pl} L σ_{y} d {\bar{ɛ}}_{f}^{pl} = \int_{0}^{{\bar{u}}_{f}^{pl}} σ_{y} d {\bar{u}}^{pl}

(5)

Figure 2.

Stress–strain curve of the polymeric matrix with progressive damage.

The characteristic length is calculated based on the element geometry. For 3D geometry, the length is equal to the cube root of the integration point volume for solid elements. The initial plastic strain for the onset of damage is chosen as 0.05 for tension and 0.5 for compression and shear considering their experimental mechanical performances.⁴⁴ $σ_{y 0}$ is the yield stress and same as the strength in the uniaxial loading conditions. ${\bar{ɛ}}_{0}^{pl}$ is the equivalent plastic strain at the onset of the damage and capable of controlling the behavior of the polymer matrix after yielding, and ${\bar{ɛ}}_{f}^{pl}$ is the equivalent plastic strain at the failure when the overall damage variable D = 1. The equivalent plastic displacement at failure is defined as ${\bar{u}}_{f}^{pl} = 2 G_{f} / σ_{y 0}$ . The overall damage variable can capture the combined effect of all active mechanisms and is computed from the individual scalar damage variables d_i, which are modeled for the degradation of the stiffness associated with each active damage mechanism. Hence, at a given time, the stress tensor in the material can be calculated as $σ_{y} = (1 - D) \bar{σ}$ . More details about the damage model and the numerical implementation can be found in Abaqus Users Manual, Version 6.13⁴¹ and Yang et al.,⁴⁴ and the matrix properties can be found in Table 1.

Cohesive element model of interface

The fiber–matrix interface is modeled using cohesive element, which is controlled by the bilinear traction-separation law. The elastic behavior is written in terms of an elastic constitutive matrix that relates the nominal stresses to the nominal strains across the interface. The nominal traction stress vector t consists of three components, t_n, t_s, t_t, which represent the normal and two shear tractions, respectively. The corresponding separations are denoted by δ_n, δ_s, and δ_t, and the original thickness of the cohesive element is denoted by T₀. Then the nominal strains can be defined as

ɛ_{n} = \frac{δ_{n}}{T_{0}}, ɛ_{s} = \frac{δ_{s}}{T_{0}}, ɛ_{t} = \frac{δ_{t}}{T_{0}}

(6)

Therefore, the elastic behavior for the cohesive element can be written in equation (7). For simplicity of computation, uncoupled behavior between the normal and shear components is desired, so the off-diagonal terms in the elasticity matrix are set to be zero and the stiffness in two shear directions is assumed to be equal

\begin{matrix} t = {\begin{matrix} t_{n} \\ t_{s} \\ t_{t} \end{matrix}} = [\begin{matrix} K_{nn} & K_{ns} & K_{nt} \\ K_{ns} & K_{ss} & K_{st} \\ K_{nt} & K_{st} & K_{tt} \end{matrix}] {\begin{matrix} δ_{n} \\ δ_{s} \\ δ_{t} \end{matrix}} = K ɛ \end{matrix}

(7)

Damage is initiated when a quadratic interaction function involving the nominal stress ratios reaches a value of one. This criterion can be represented as

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

(8)

where

〈 〉

is the Macaulay brackets, which return the argument if positive and zero otherwise, and the bracket is also used to signify that a pure compressive deformation or stress state does not initiate damage.

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

represent the peak values of the nominal stress when the deformation is purely normal to the interface and purely in the first and second shear directions, respectively. Damage evolution is defined based on the dissipated fracture energy during the damage propagation. Once the damage initiates, the traction stress t⁰ is reduced depending on the interfacial damage parameter, which monotonically evolves from 0 (in the absence of damage

δ^{0}

) to 1 (at the final failure δ^f), as shown in Figure 3. The energy-based Benzeggagh–Kenane (BK law) damage propagation criterion⁴⁵ is adopted during the damage evolution of the cohesive elements, which is illustrated as below

\begin{matrix} G_{n}^{C} + (G_{s}^{C} - G_{n}^{C}) (\frac{G_{S}}{G_{T}})^{η} = G^{C}, G_{S} = G_{s} \\ + G_{t}, G_{T} = G_{n} + G_{s} \end{matrix}

(9)

where

G_{n}^{C}, G_{s}^{C}

, and

G_{t}^{C}

refer to the critical fracture energies required to cause failure in the normal, the first, and the second shear directions, respectively. Here

G_{s}^{C} = G_{t}^{C}

is set and G^C is the critical total interfacial fracture energy and the η is a cohesive property parameter. The properties of the interface are shown in Table 2, where the elastic properties are obtained from the parameter identification.

Figure 3.

Traction-separation law of the interface.

Table 2.

Mechanical properties of the interface.^23,34

Mechanical properties	Value
Thickness T₀ (µm)	0.139
Normal stiffness K_nn (GPa/mm)	846.8
Shear stiffness K_ss (K_tt) (GPa/mm)	535.6
Normal strength $t_{n}^{0}$ (MPa)	27.5
Shear strength $t_{s}^{0}$ ( $t_{t}^{0}$ ) (MPa)	45.0
Normal critical fracture energy $G_{n}^{C}$ ( ${J / m}^{2}$ )	5
Shear critical fracture energy $G_{s}^{C}$ ( $G_{t}^{C}$ ) ( ${J / m}^{2}$ )	100

Periodic boundary conditions

PBC is imposed on the corresponding surfaces of the RVE by means of introducing the equations between the periodic nodes, in order to guarantee the periodicity of the displacement and traction fields as well as to ensure the continuity between the neighboring RVEs. The unified PBCs are expressed in terms of the displacement vectors $\vec{U_{1}}, \vec{U_{2}}$ , and $\vec{U_{3}}$ which are related to the displacements between the opposite surfaces

\begin{matrix} {\begin{matrix} \vec{u} (0, x_{2}, x_{3}) - \vec{u} (L_{1}, x_{2}, x_{3}) = \vec{U_{1}} \\ \vec{u} (x_{1}, 0, x_{3}) - \vec{u} (x_{1}, L_{2}, x_{3}) = \vec{U_{2}} \\ \vec{u} (x_{1}, x_{2}, 0) - \vec{u} (x_{1}, x_{2}, L_{3}) = \vec{U_{3}} \end{matrix} \end{matrix}

(10)

where L₁, L₂, and L₃ are the lengths of the RVE along with three orthogonal directions, respectively. The absolute formation is used to impose the linear constraint functions on the nodes at parallel opposite pairs of faces, and the edges and the vertices are extracted from the face which they belong to for the equations. The dummy nodes are introduced as reference points to apply the load in three directions, and then different loading conditions can be achieved by applying the displacement loads on the dummy nodes. When a displacement component of the dummy node is set free, this displacement can be computed by the FE solver under stress-free conditions. Therefore, the Poisson effect is permitted in a specific direction. For the uniaxial loading conditions, the load is applied to the axial direction, and then the axial components in the other two directions are set free for the consideration of the Poisson effects. While for combined uniaxial transverse loading perpendicular to the fiber direction and in-plane shear along the fiber, the loading is imposed with

\vec{U_{2}} = (δ_{s}, \pm δ, 0)

, where the

\pm δ

and δ_s represent the transverse tension/compression and shear displacements, respectively. Same as the uniaxial loading conditions, the axial components in the other two directions are set free. More details about the implication of the PBC on the RVE is found in Garoz et al.³⁹

Assessment of the mechanical response of the RVE under uniaxial loading conditions

Before the prediction of the damage mode and envelop of the UD RVE under biaxial stress states is conducted, an assessment of the RVE model, with a reasonable fine mesh, appropriate boundary condition and material properties, is required. In order to retain the accuracy of the numerical results with a reasonable computational cost, five RVE models with different fiber distributions are generated with around 52,000 elements, and the primary results are summarized in Table 3. It is worth noting that the transverse modulus was obtained by averaging Young’s modulus from transverse tension and compression cases, and the average value is almost identical with the value obtained from the idealized single fiber model.³⁴ All of the predicted properties are in good agreement with the single fiber model within a difference of 7%. The details of the stress–strain curve and ductile damage distribution of RVE1 can be found in Figure 4, where all the simulations were practically superposed in the elastic regime, divergences arose at the onset of matrix plastic deformation and interface debonding, then increased/decreased when the composite strength reached its plateau.

Figure 4.

Predicted stress–strain curve of five RVE models with different fiber distributions (left) and predicted damage variable contour of RVE1 (right) under different uniaxial loading conditions: (a) transverse tension; (b) transverse compression; (c) in-plane shear.

Table 3.

Comparison of predicted elastic and plastic properties of five random RVE models with the idealized single fiber model.³⁴

	RVE1	RVE2	RVE3	RVE4	RVE5	Average	Single fiber model³⁴
Transverse modulus
$E_{22} (GPa)$	9.322	10.1	9.765	9.51	9.741	9.688	9.68
In-plane shear modulus
$G_{12} (GPa)$	4.365	4.6	4.44	4.3	4.34	4.41	4.18
Transverse tensile
strength $Y_{T} (MPa)$	38.6	38.5	37.5	39.3	37.75	38.33	41.14
Transverse compressive
strength $Y_{C} (MPa)$	169.4	176.4	177.7	191.4	181	179.18	169.87
In-plane shear strength
$S_{12} (MPa)$	49.2	46	47.8	46.1	47.3	47.28	51.12

RVE: representative volume element.

During the pure transverse tension loading condition, the damage process is mainly dominated by the interface debonding. Non-linearity of normal stress appeared in the stress–strain curve when it meets the interface strength, which is shown in Figure 4(a) (left) and the mechanism of deformation and damage in tension is observed in Figure 4(a) (right), which shows the contour plot of the ductile damage distribution in matrix. The crack always starts at the interface close to fiber cluster, where the fibers are more closer to each other, resulting in the stress concentration in these regions. The damage in the interface initiates firstly by the nucleation of interface cracks, responsible for the non-linearity behavior at small strains, and grows along the interfaces. After damage initiates from the interface, the matrix experiences severe plastic deformation due to the stress concentrations at the interface crack tip. Then the damage propagates along the weakest path at the relatively matrix poor region determined by the spatial distribution of interface cracks, perpendicular to the loading direction. Final failure occurs by linking up the interface cracks through the RVE.

During the pure transverse compression, the damage initiates from the interface in the form of debonding, then the matrix bears the load afterward. After the damage initiation, matrix plastic damage starts to develop at the vicinity of the debonding positions, and the curves in Figure 4(b) (left) become stable, which indicate that the compressive strength is dominated by the matrix. Finally, the matrix cracks at different locations are linked together to form the main shear band passing through the interface cracks with an orientation of 52° relative to the loading direction, which is very close to the published findings^24,44 and is shown in Figure 4(b) (right). All of the five simulations under transverse compression are superposed in the elastic regime and started to diverge at the onset of the matrix plastic deformation. The largest strength among the five RVE models is found in RVE4, which is probably due to boundary effects and the rotation of the fiber as one fiber was divided into a same part at four corners of the RVE model.

During the pure in-plane shear loading, the damage of the UD composite lamina is dominated by interfacial decohesion or by matrix yielding, depending on the interface strength.⁴⁶ Here in this study, the interface strength is smaller than the matrix shear limit, thus, the crack initiates from the interface in a form of interfacial debonding and then the matrix holds the shear loads in the most time. Therefore, the damage starts from the interface and propagates as a result of the interactions between the interfacial debonding and the matrix plastic deformation by linking them together. Several different damage surfaces, passing through the RVE perpendicular to the loading condition, can be found in the RVE at different locations, which is shown in Figure 4(c) (right).

FEM prediction of the damage initiation and propagation and failure strength under transverse and in-plane shear loadings

Once the validation of the RVE model is completed, it could be used to predict damage modes and failure envelop of the UD composite materials subjected to biaxial loading conditions. The topics about the transition damage modes between transverse and in-plane shear stresses, and the failure loci for the whole range of this combined stress state are covered. The scheme of the combined transverse tension/compression and in-plane shear stress state is depicted in Figure 5. Two different loading paths can be applied to the RVE for the combined stress state: the first is to apply a transverse load up to a prescribed stress state and then a shear load is applied until the final failure while the total transverse load acting on the RVE is held constant; the second is to apply the transverse and shear displacement simultaneously to the RVE by a proportional amount, which is designated as radial loading path. It is worth noting that no significant influence of the loading path on the failure loci was found from the biaxial experimental results^47,48 and the combined transverse compression and shear numerical simulations.^47,49 Therefore, here in this study, the radial loading path was applied to the RVE for the combined transverse tension/compression and in-plane shear stress states.

Figure 5.

Schematic of the RVE of UD composite subjected to combined in-plane shear and transverse stresses.

Transition of the damage mechanism

One of the most interesting phenomena in the investigation of progressive damage of composite lamina under transverse and in-plane shear is the finding of transverse point. Two different damage modes exist at the same time around the transition point, and only one damage mode can be observed when it is far away from the transition point. Here in this study, six different loading cases were selected to reveal the transition modes and transition of damage mechanisms. The mechanical behavior of the composite lamina under different loading paths characterized by $γ_{12} / ɛ_{22}$ is shown in Figure 6 in the $σ_{22} - τ_{12}$ stress space for a lamina with a weak interface. Three points on the curve for each case represent the progressive damage process at different strain stages, and the corresponding damage modes and mechanisms can be found in Figures 7 and 8, which show the contour plots of the plastic strain at different stages. It can be found in Figure 6 that normal and shear stresses increase proportionally during the elastic deformation, and this proportionality disappeared with the onset damage in the interface. During the transverse tension and in-plane shear loadings, the normal or shear stress falls very sharply, depending on if it is tension-dominated or shear-dominated damage. While during the transverse compression and in-plane shear, the curves differ from each other and all of them experience a slightly hardening in shear stress and then decrease until failure.

Figure 6.

Mechanical response of the composite lamina subjected to transverse and in-plane shear loadings in the $σ_{22} - τ_{12}$ stress space with a weak interface. The $γ_{12} / ɛ_{22}$ for each curve represents the ratio of shear strain to the transverse tension/compression strain.

Figure 7.

Contour plot of the accumulated plastic strain in the composites subjected to biaxial transverse tension and in-plane shear: (a) $γ_{12} / ɛ_{22} = 0.83 (ɛ_{22} = 1.2 %, γ_{12} = 1 %)$ , (b) $γ_{12} / ɛ_{22} = 12.5 (ɛ_{22} = 0.8 %, γ_{12} = 10 %)$ , and (c) $γ_{12} / ɛ_{22} = 50 (ɛ_{22} = 0.3 %, γ_{12} = 15 %)$ .

Figure 8.

Contour plot of the accumulated plastic strain in the composites subjected to biaxial transverse compression and in-plane shear: (a) $γ_{12} / ɛ_{22} = - 1.25 (ɛ_{22} = - 4 %, γ_{12} = 5 %)$ , (b) $γ_{12} / ɛ_{22} = - 20 (ɛ_{22} = - 1.6 %, γ_{12} = 32 %)$ , and (c) $γ_{12} / ɛ_{22} = - 200 (ɛ_{22} = - 0.16 %, γ_{12} = 32 %)$ .

Transition of the damage mechanism under transverse tension and in-plane shear loadings

At the low ratio ( $γ_{12} / ɛ_{22} = 0.83$ ), the damage is dominated by the interface cracks due to tension traction, and the damage surface is formed by linking up the interface cracks at the relatively matrix poor regions. And the damage plane is always generated perpendicular to the loading direction, which is very similar to the pure transverse tension case. However, when the ratio increases to a moderate value ( $γ_{12} / ɛ_{22} = 12.5$ ), the shear load starts to play a role in the combined stress state. In this case, more cracks can be found at the interface near the neighboring fibers (Figure 7(b1)) and as the damage propagated, large plastic strain accumulation was detected near the interface crack (Figure 7(b2)). At the final damage stage (Figure 7(b3)), two damage planes were found existing at the same time, one is same as the damage at low ratio (Figure 7(a3)) due to the interfacial debonding, and the other one is almost formed due to the shear stress. The same fracture plane was also captured in the high ratio (Figure 7(c3)) where the shear failure is dominant and the position of the other fracture plane is also different from the one found in the moderate ratio. In this case, the damage mode is similar to the pure in-plane shear case where the shear band was formed by linking the interface cracks with matrix plastic deformation regions together. At the moderate ratio where two damage modes exist at the same time, this point can be called as transition point and near this point, two different damage modes can be captured. It was found in this study that the transition point is around $γ_{12} / ɛ_{22} = 12.5$ when a lamina is subjected to transverse tension and in-plane shear. Therefore, beyond this point, the perfect plastic definition of the Drucker–Prager for the matrix should be changed from tension to shear.

Transition of the damage mechanism under transverse compression and in-plane shear loadings

For the transverse compression and in-plane shear cases, the simulation results from the lower ratio ( $γ_{12} / ɛ_{22} = - 1.25$ ) suggest that the cracks initiated from the interface and propagated by the formation of plastic shear band of matrix. This can happen by linking up the interface cracks, shown as Figure 8(a3), which is very similar to the pure transverse compression cases. When the ratio decreases to $γ_{12} / ɛ_{22} = - 20$ , more cracks at the neighboring fibers were captured from the onset of damage (Figure 8(b1)), and shear band begins to vanish due to the shear stress becoming a dominant factor controlling the damage mode. However, the compression-dominated and shear-dominated damage modes were found existing at the same time, compared to the low ratio case. Beyond this transition point, the perfect plasticity definition for the matrix in the Drucker–Prager should be switched to shear from compression. The transition of fracture angle was found from 52° ( $γ_{12} / ɛ_{22} = - 1.25$ ) to 0° ( $γ_{12} / ɛ_{22} = - 20$ ) which is in agreement with the numerical findings²³ considering friction effects. In the high shear ratio ( $γ_{12} / ɛ_{22} = - 200$ ), the composite fails in shear-dominated mode and the final fracture angle/shear plane is similar to the pure in-plane shear.

Comparison between the classical failure criteria and numerical results with a weak interface

Currently, there exist a large number of failure criteria which are mostly stress-based, including fully interactive criteria such as Tsai and Wu¹⁰ and damage mode based criteria like Hashin¹² and Puck and Schürmann.^13,50 The Tsai and Wu¹⁰ criterion predicts failure with a highly integrated equation involving all stress components by combining different damage mechanisms. Different from the Tsai and Wu criterion, the Hashin criterion is capable of distinguishing the fiber and matrix damage initiation in composite materials, and either one is further subdivided into two damage mechanisms such as tensile and compressive modes. Despite the capability of this criterion for predicting the damage in the lamina under uniaxial loading, numerous studies showed that it does not always agree with the experiments accurately, especially under the combined transverse compression and in-plane shear. This disadvantage is due to neglecting the shear hardening effects with the presence of transverse compression.¹⁴ Hashin’s criterion has been further developed by Puck and Schürmann^13,50 by addressing matrix compression failure with a model based on Mohr–Coulomb hypothesis, which assumes that damage is triggered due to the normal stress and tangential stress, acting on the fracture plane with a specific inclination angle to the material plane.

The predictions of those models are based on the mechanical tests, which should provide the lamina strength under transverse tension Y_T, transverse compression Y_C, and in-plane shear S₁₂. As there was no experimental result available, these magnitudes were obtained from the computational micromechanics simulations,³⁴ which is presented in Table 3. The failure envelopes under transverse loads and in-plane shear ( $σ_{22} - τ_{12}$ ) provided by these criteria are compared to numerical simulations, as shown in Figure 9, which indicates large variances among different criteria, especially under transverse compression range. For $σ_{22} > 0$ , the predictions by these three criteria are almost identical to the computational results when the ratio $γ_{12} / ɛ_{22} < 8$ , which is tensile dominant damage triggered by the interface decohesion. However, when $γ_{12} / ɛ_{22} > 12.5$ , the Puck criterion is always consistent with the numerical results, when the plastic yielding definition in Drucker–Prager is set to shear, while the Hashin criterion overestimates. When the ratio is beyond 10, the Puck criterion predicted the failure strength accurately, with the matrix yielding setting to shear; however, when it was set to tension yielding, all of the criteria overestimated the failure strength. That is because the shear failure in the interface is dominant in the high shear load, while these models assume that the main deformation is shear yielding of the matrix instead of the interface decohesion which was found in the numerical results.

Figure 9.

Comparison of the predicted failure surface and Puck, Hashin, and Tsai–Wu failure criteria subjected to transverse load and in-plane shear. The numbers next to points represent the ratio of shear strain (γ₁₂) to transverse strain ( $ɛ_{22}$ for either tension or compression).

For $σ_{22} < 0$ , there existed big differences between the numerical results and the failure criteria, especially when $γ_{12} / ɛ_{22} < - 5$ , while for $γ_{12} / ɛ_{22} > - 2.5$ , Hashin’s criterion is consistent with the numerical results; however, the Puck and Tsai–Wu criteria slightly overestimated the failure strength. Under compression-dominant failure, the damage is caused by the nucleation of interface cracks at the equator of the fibers, which is perpendicular to the loading direction, and continued by the plastic deformation of the matrix in the form of the shear band which links up these cracks. Shear failure is dominant when $γ_{12} / ɛ_{22} < - 20$ , when the interface damage occurs rather than matrix shear yielding. As the interface shear strength does not depend on the compressive normal stress, thus the shear strength of the composite under these loading conditions is independent of compressive stress. Hence, the strength was found to be practically constant in the range of $γ_{12} / ɛ_{22} < - 20$ , no matter the hardening behavior of matrix is set to compression or shear. Nevertheless, the strength predicted by computational micromechanics with shear hardening is slightly higher than the one with compression hardening. In addition, the largest difference between the failure criteria and the numerical results starts from the transition point where $γ_{12} / ɛ_{22} = - 20$ in the shear-dominated region, as criteria predict an increase of the composite strength when the compressive stress increases, which is not realistic for the case of interface-dominated failure.²⁶

Once the transition point subjected to the biaxial loading was found, an appropriate yielding behavior of the matrix described by the Drucker–Prager should be chosen after the point. The lowest and highest strengths of the lamina under different uniaxial loadings from Table 3 were selected as the lower bound and upper bound to plot the failure envelopes from the Puck, Tasi–Wu, and Hashin criteria. These envelopes were compared to the numerical simulations obtained from the five RVE models under biaxial loadings, which can be found in Figures 10 to 12. It was found that during the tensile dominant failure period when $γ_{12} / ɛ_{22} < 8$ , all of the theoretical predictions are in good agreement with the results obtained from numerical simulations. However, when it is beyond the transition point ( $γ_{12} / ɛ_{22} = 12.5$ ) from tensile dominant failure to shear one, the difference occurred. Only the predictions of lower and upper bounds by Hashin criteria can cover all the five numerical simulation results, while the other two slightly underestimate the strength beyond the transition point. When it comes to transverse compression cases, large difference occurred around the transition point. Shear hardening under moderate transverse compressive stress was observed experimentally⁵¹ and is predicted by Puck and Tsai–Wu criteria, as shown in Figures 10 and 11. However, the prediction by the Puck criteria is based on the assumption that the matrix-dominated damage occurs through the surface parallel to the fibers, and it is caused by the normal and tangential stresses acting on the fracture plane, which is not realistic in the case of interface-dominated failure. The Hashin criterion compares well to the numerical simulations when the effects of shear hardening with the presence of transverse compression is neglected, which is suitable for most cases of interface-dominated failure. Below this transition point in the transverse compression side, the prediction by Hashin is better than the other two criteria in most cases.

Figure 10.

Comparison of numerical simulations of UD composite and failure envelopes predicted by Puck criterion under combined transverse and in-plane shear stress states. The numbers next to points represent the ratio of shear strain (γ₁₂) to transverse strain ( $ɛ_{22}$ for either tension or compression).

Figure 11.

Comparison of numerical simulations of UD composite and failure envelopes predicted by Tasi–Wu criterion under combined transverse and in-plane shear stress states. The numbers next to points represent the ratio of shear strain (γ₁₂) to transverse strain ( $ɛ_{22}$ for either tension or compression).

Figure 12.

Comparison of numerical simulations of UD composite and failure envelopes predicted by Hashin criterion under combined transverse and in-plane shear stress states. The numbers next to points represent the ratio of shear strain (γ₁₂) to transverse strain ( $ɛ_{22}$ for either tension or compression).

Comparison between the classical failure criteria and numerical results with a strong interface

As the interface strength can influence the damage modes of UD composite largely subjected to biaxial stress states, it is worth conducting another set of simulations with a strong interface subjected to the same biaxial stress states. The influences of the interface strength was investigated in Canal et al. and Totry et al.^3,26 in the form of a weak and a strong interface subjected to different loading conditions. It was found that damage of composites is controlled by the matrix with a strong interface, and controlled by the interface with a weak interface under transverse compression and longitudinal shear,²⁶ as well as transverse tension and out-of-plane shear.³ Here in this study, six different biaxial loading cases, covering the transition point from compression-dominated failure to shear-dominated shear, were conducted on the five RVE models with different fiber distributions, and 112 MPa was chosen as the interface shear strength to represent a good fiber/matrix bonding.²⁶

The comparison of numerical results and different failure criteria is shown in Figure 13, where the lower and upper bounds for plotting the failure envelopes come from the minimum and maximum strength values of five RVE models under uniaxial loadings. It is worth noting that the numerical simulations can explain the transition from shear failure to compressive failure with different damage modes. The shear failure is characterized by a linear increases of the shear strength until the transition point $γ_{12} / ɛ_{22} = - 20$ , when the applied compressive stress increased. The end of this regime occurred at the transition point, beyond which the damage mode was changed to compressive failure, resulting in a gentle reduction in the maximum shear strength at failure. It should be noted that the agreement is excellent between the numerical results and the whole shape of the failure locus of the Puck and Tasi–Wu failure criteria. Both of these criteria can predict the transition point ( $γ_{12} / ɛ_{22} = - 20$ ) and the failure surface trend accurately, while Tasi–Wu failure criterion slightly overestimated the shear strength from $γ_{12} / ɛ_{22} = - 40$ and the Puck failure criterion can give a better prediction in a higher ratio due to its matrix-dominated damage mode assumption that failure under transverse compression and in-plane shear occurs through the surfaces parallel to the fibers and the damage is triggered by the normal and tangential stresses acting on the fracture plane. However, the Hashin failure criterion always underestimated the strength of UD composite subjected to transverse compression and in-plane shear, and it cannot predict the increase in shear strength in the presence of moderate compressive stresses. This limitation is due to the fact that the model only assumes a quadratic interaction between the stress invariants instead of determining the actual fracture plane.

Figure 13.

Comparison of numerical results and failure envelopes of UD composite predicted by (a) Puck, (b) Tasi–Wu, and (c) Hashin subjected to combined transverse compressive and in-plane shear stress states. The numbers next to points represent the ratio of shear strain (γ₁₂) to transverse compression strain ( $ɛ_{22}$ ).

Conclusions

In this study, a comprehensive set of RVE based on computational micromechanics is developed in commercial FEM software ABAQUS/Explicit to investigate the failure/damage mechanisms and failure envelopes of UD TC33 carbon fiber UD composite lamina subjected to biaxial stress states. Random fiber distributions are generated from 2D DEM models and the RVEs are obtained in forms of three phases including identified interface properties. Taking the statistic into account, five RVEs are adopted for each case. The transition from compression or tension to shear-dominant damage and the orientation of the fracture plane are adequately predicted by the FE simulations, and the key findings are shown as below:

The transition from compression-dominated failure to shear-dominated failure arises when shear strain is larger than 20 times the absolute compressive strain. Thus, during the shear-dominated failure cases, the hardening of matrix should be switched to shear in the D–P model in ABAQUS/Explicit. This also applies to the transition from tension to shear when shear strain is larger than 12.5 times the tension strain.

During the combined transverse tension and in-plane shear loading conditions, the failure of composites with a weak interface is controlled by nucleation and growth of interface tensile cracks in the tension-dominated damage mode. However, when it is beyond the transition point, the failure is controlled by the interface shear cracks in the shear-dominated damage mode. During the transverse compression and in-plane shear loading conditions, the damage is controlled by the nucleation of the interface crack and propagates by the formation of shear band in the matrix under compression-dominated failure. However, when it is beyond the transition point and enters into the shear-dominated failure, the damage is controlled by the interface shear fracture.

The effects of the interface strength on the mechanical response of the composite lamina are investigated, and the upper bound and lower bound predictions of the failure criteria are introduced to eliminate their dependence on the uniaxial strength data. The results clearly revealed that the Hashin failure criterion provides a better estimation for the failure locus of the UD lamina subjected to transverse and in-plane shear loadings, when the mechanical behavior of the composite with a weak interface is controlled by the interface decohesion. However, when it comes to matrix-dominated failure cases where the interface strength is larger than the matrix shear strength, the Puck failure criterion can provide a better estimation for the failure locus under the same biaxial loading conditions. It is hard to say which one is better because in a specific situation, depending on the interface properties, different failure criterion has its own strengths and drawbacks. Comparing these three criteria, it can be seen the Tsai–Wu may be generally better than both of others as it presents more neutral predictions in both of the examined cases. This could be due to the fact that Tsai–Wu criterion is a generalized stress-based criterion but this does not in any way suggest that Tsai–Wu criterion is always more accurate for other cases.

The FE modeling is useful for validation of current composites failure criteria, especially for those loading conditions where experiments are extremely difficult to conduct. These results showed the potential of computational micromechanics simulations to assist the modification of existing failure criteria and the development of new failure criteria for composites in general. With the help of our current models, the complex loading conditions such as combined longitudinal and shear stress states and triaxial loading conditions will be considered. It is important to note that the friction condition between the fiber and matrix is not considered in these cases, which is the next topic in our future work.

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: Lei Wan would like to thank China Scholarship Council (CSC) and the School of Civil Engineering at University of Leeds for their financial support for his PhD studies.

ORCID iDs

Ping Zhu

Dongmin Yang

References

Bunsell

Gorbatikh

Morton

, et al. Benchmarking of strength models for unidirectional composites under longitudinal tension. Compos Part A: Appl Sci Manuf 2018; 111: 138–150.

Hinton

Kaddour

Soden

. Failure criteria in fibre reinforced polymer composites: The world-wide failure exercise, New York: Elsevier, 2004.

Canal

Segurado

LLorca

. Failure surface of epoxy-modified fiber-reinforced composites under transverse tension and out-of-plane shear. Int J Solids Struct 2009; 46(11): 2265–2274.

Hart-Smith

. Predictions of the original and truncated maximum-strain failure models for certain fibrous composite laminates. Compos Sci Technol 1998; 58(7): 1151–1178.

Wolfe

Butalia

. A strain-energy based failure criterion for non-linear analysis of composite laminates subjected to biaxial loading. Compos Sci Technol 1998; 58(7): 1107–1124.

Christensen

. Stress based yield/failure criteria for fiber composites. Int J Solids Struct 1997; 34(5): 529–543.

Cui

Wisnom

. A combined stress-based and fracture-mechanics-based model for predicting delamination in composites. Composites 1993; 24(6): 467–474.

Whitney

Nuismer

. Stress fracture criteria for laminated composites containing stress concentrations. J Compos Mater 1974; 8(3): 253–265.

Sun

Tao

. Prediction of failure envelopes and stress/strain behaviour of composite laminates. Compos Sci Technol 1998; 58(7): 1125–1136.

10.

Tsai

. A general theory of strength for anisotropic materials. J Compos Mater 1971; 5(1): 58–80.

11.

Hoffman

. The brittle strength of orthotropic materials. J Compos Mater 1967; 1(2): 200–206.

12.

Hashin

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

13.

Puck

Schürmann

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

14.

Pinho

Iannucci

Robinson

. Physically-based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: Part I: Development. Compos Part A: Appl Sci Manuf 2006; 37(1): 63–73.

15.

Davila

Camanho

Rose

. Failure criteria for FRP laminates. J Compos Mater 2005; 39(4): 323–345.

16.

Hinton

Soden

. Predicting failure in composite laminates: The background to the exercise. Compos Sci Technol 1998; 58(7): 1001–1010.

17.

Hinton

Kaddour

Soden

. Evaluation of failure prediction in composite laminates: Background to ‘Part C’ of the exercise. Compos Sci Technol 2004; 64(3): 321–327.

18.

Kaddour

Hinton

. Maturity of 3D failure criteria for fibre-reinforced composites: Comparison between theories and experiments: Part B of WWFE-II. J Compos Mater 2013; 47(6–7): 925–966.

19.

Christensen

. The world wide failure exercise II examination of results. J Reinf Plast Compos 2013; 32(21): 1668–1672.

20.

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(20–21): 2427–2442.

21.

Kaddour

Hinton

Smith

, et al. The background to the third world-wide failure exercise. J Compos Mater 2013; 47(20–21): 2417–2426.

22.

Kaddour

Hinton

Smith

, et al. A comparison between the predictive capability of matrix cracking, damage and failure criteria for fibre reinforced composite laminates: Part A of the third world-wide failure exercise. J Compos Mater 2013; 47(20–21): 2749–2779.

23.

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.

24.

González

LLorca

. Mechanical behavior of unidirectional fiber-reinforced polymers under transverse compression: Microscopic mechanisms and modeling. Compos Sci Technol 2007; 67(13): 2795–2806.

25.

Totry

González

LLorca

. Failure locus of fiber-reinforced composites under transverse compression and out-of-plane shear. Compos Sci Technol 2008; 68(3): 829–839.

26.

Totry

González

LLorca

. Prediction of the failure locus of C/PEEK composites under transverse compression and longitudinal shear through computational micromechanics. Compos Sci Technol 2008; 68(15): 3128–3136.

27.

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(1): 7–12.

28.

Arteiro

Catalanotti

Melro

, et al. Micro-mechanical analysis of the in situ effect in polymer composite laminates. Compos Struct 2014; 116: 827–840.

29.

Melro

Camanho

Pires

, et al. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part II - Micromechanical analyses. Int J Solids Struct 2013; 50(11): 1906–1915.

30.

Arteiro

Catalanotti

Melro

, et al. Micro-mechanical analysis of the effect of ply thickness on the transverse compressive strength of polymer composites. Compos Part A: Appl Sci Manuf 2015; 79: 127–137.

31.

Tan

Naya

Yang

, et al. The role of interfacial properties on the intralaminar and interlaminar damage behaviour of unidirectional composite laminates: Experimental characterization and multiscale modelling. Compos Part B: Eng 2018; 138: 206–221.

32.

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.

33.

Zhu

, et al. Identification of the mechanical properties of the carbon fiber and the interphase region based on computational micromechanics and Kriging metamodel. Comput Mater Sci 2014; 95: 172–180.

34.

Zhu

Liu

, et al. Numerical investigation of fiber random distribution on the mechanical properties of yarn in-plain woven carbon fiber-reinforced composite based on a new perturbation algorithm. J Compos Mater 2018; 52(6): 755–771.

35.

Melro

Camanho

Pires

, et al. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part I - Constitutive modelling. Int J Solids Struct 2013; 50(11): 1897–1905.

36.

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.

37.

Wang

Yang

Sun

, et al. Dynamic explicit FE modeling of hot ring rolling process. T Nonferr Metal Soc 2006; 16(6): 1274–1280.

38.

Wang

Long

. Investigation of material deformation in multi-pass conventional metal spinning. Mater Des 2011; 32(5): 2891–2899.

39.

Garoz

Gilabert

Sevenois

, 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.

40.

Ward

. Review: The yield behaviour of polymers. J Mater Sci 1971; 6(11): 1397–1417.

41.

Abaqus Users Manual, Version 6.13, Dassault Systémes Simulia Corp., Providence, Rhode Island, USA.

42.

Subramanya

Viswanath

Narasimhan

. A three-dimensional numerical study of mode I crack tip fields in pressure sensitive plastic solids. Int J Solids Struct 2007; 44(6): 1863–1879.

43.

Fiedler

Hojo

Ochiai

, et al. Failure behavior of an epoxy matrix under different kinds of static loading. Compos Sci Technol 2001; 61(11): 1615–1624.

44.

Yang

Yan

Liu

, et al. Microscopic failure mechanisms of fiber-reinforced polymer composites under transverse tension and compression. Compos Sci Technol 2012; 72(15): 1818–1825.

45.

Benzeggagh

Kenane

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

46.

Totry

Molina-Aldareguía

González

, et al. Effect of fiber, matrix and interface properties on the in-plane shear deformation of carbon-fiber reinforced composites. Compos Sci Technol 2010; 70(6): 970–980.

47.

Vogler

Kyriakides

. Inelastic behavior of an AS4/PEEK composite under combined transverse compression and shear. Part I: Experiments. Int J Plast 1999; 15(8): 783–806.

48.

Vogler

Hsu

S-Y

Kyriakides

. Composite failure under combined compression and shear. Int J Solids Struct 2000; 37(12): 1765–1791.

49.

Totry

González

LLorca

. Influence of the loading path on the strength of fiber-reinforced composites subjected to transverse compression and shear. Int J Solids Struct 2008; 45(6): 1663–1675.

50.

Puck

Schürmann

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

51.

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.