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

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.


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. 1 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][5][6][7][8][9] Some failure criteria require expensive and timeconsuming 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][12][13][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. 13 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 5 and Hashin 7 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. 17 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, 25 combined longitudinal tension and in-plane shear 1 and triaxial loadings. 26 Sun et al. 1 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. 27 Six out of 10 representative biaxial loading conditions were considered due to the transverse isotropy of the crosssection 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. 27 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 5 and Hashin 7 ) and micromechanical simulations. The experimentally validated computational approach provides fast screening capabilities to improve material down-selection as part of a design cycle.

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. 21 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 28 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 firstorder 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 À6 s. 21,26 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, 29 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. 30 Constitutive models, reported in the literature, to describe the mechanical behaviour of polymeric materials under multiaxial loading 31 include the extended Drucker-Prager (D-P) yield model, associated with a ductile damage criterion, 21 the modified Drucker-Prager plastic damage model, 19 and the elastoplastic model with an isotropic damage constitutive model. 32 In this study, the polymer matrix was modelled as an isotropic elastoplastic solid. The modified Drucker-Prager plastic damage model 33 was used to model the mechanical behaviour of epoxy under multiaxial stress states. The yield surface of the epoxy is expressed by where I 1 stands for the first invariant of the stress tensor, J 2 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 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. 19 Figure 2 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.
where t 0 n , t 0 s and t 0 t 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 δ 0 s ) to 1 (at final failure δ f s ), which is shown in Figure 2(b). The energy-based Benzeggath-Kenane damage propagation criterion 35 is adopted to consider the dependence of the fracture energy dissipation on fracture modes during the damage evolution of the cohesive elements, which reads where η C BK is the BK power exponent, G C n and G C s are the normal and shear fracture energies, respectively. G n and G s are the reciprocal work under mixed mode damage propagation.
The calibrated interface strengths in normal and shear directions from our previous study 26 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 19 and experimental studies. 36 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 2 . 19 In this study, the fracture energy of 2 J/m 2 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 2 . 19 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. 39 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 δ 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 τ 12 σ 2 ¼ 0.24, 0.49, 0.9, 1.5, 3.5.

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.

Tsai-Wu failure criterion
The Tsai-Wu failure criterion 5 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 where F i , F ij are the coefficients associated with the material strengths determined using experiments. Considering the random distribution of fibres within the matrix at the crosssection 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 where the coefficients can be calculated as shown below 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. 40

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 7 : Computational results and discussions

Transverse dominated loadings
Failure envelope comparison.  13 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, 41 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. 27 It can be found that both failure criteria agree reasonably well with experimental data and numerical 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 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°f or 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 τ 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 sheardominated loads. 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 τ 12 σ 2 ¼ 1.5, one main shear fracture plane is observed close to the   left edge, while the shear fracture plane becomes random when τ 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 19 with a fracture angle of 56°and by experimental observations with an angle of 53°, 40 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.

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.
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 crosssection in composites, which changes the direction of the fracture propagation.

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 6 and experimentally. 42 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.
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 outof-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, 44 IM7/8552 8 and E-glass/ RP528. 45 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 inplane 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, 6 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 σ 0  is characterised by the equivalent plastic strain PEEQ). τ 0 where the recommended range of p ðÀÞ 'k is 0.25-0.30 for carbon fiber reinforced composites. 6 Here, p ðÀÞ 'k = 0.25 is selected for the calculation of p ðÀÞ '' according to equation (13).

Conclusion
This study introduced a 3D high-fidelity micromechanicsbased 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) highfidelity 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 19 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.

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.