Probability embedded failure prediction of unidirectional composites under biaxial loadings combining machine learning and micromechanical modelling

This study presents a data-driven, probability embedded approach for the failure prediction of IM7/8552 unidirectional carbon fibre reinforced polymer (CFRP) composite materials under biaxial stress states based on micromechanical modelling and artificial neural networks (ANNs). High-fidelity 3D representative volume element (RVE) finite element models were used for the generation of failure points. Fibre failure and the friction between fibres and matrix after fibre/matrix debonding were taken into consideration and implemented as VUMAT subroutines, respectively. Uncertainty quantification was conducted based on a coupled experimental– numerical approach and failure probabilities were inserted into the failure points to generate the database for the training of ANNs. A total of 15 biaxial stress combinations were considered for the generation of datasets. Two strategies were considered for the construction of form-free failure criteria based on the ANNs for regression and classification problems. It is found that for the regression problems, an ANN model with 2 hidden layers and 64 neurons can achieve a mean square error (MSE) of 0.027% and a mean absolute error (MAE) of 0.78%. For the classification problems, an ANN model with 3 hidden layers and 32 neurons, presents an excellent performance in the prediction with a probability of 98.1%. A good agreement was observed between the failure strength of composites under transverse and in-plane shear predicted by these ANNs and failure envelopes theoretically predicted by Tsai–Wu and Hashin failure criteria.


Introduction
Carbon Fibre Reinforced Polymer (CFRP) composites are finding increasing applications in aerospace, automotive, marine and many other industries due to their excellent strength and stiffness per unit weight, as well as the capability of tailoring to achieve desired stiffness and strength in a specific direction.However, despite years of extensive research on the failure analysis of CFRP composites, a complete and validated methodology for predicting their failure, including the progressive failure process, has yet to be fully realised.This is largely due to their complex failure mechanisms compared to isotropic materials.CFRP composites are hierarchical, spanning three different length scales (i.e.micro-, meso-and macro-), and it is this 'tyranny of scales' that makes failure prediction challenging, especially under multiaxial loading conditions.
A large number of failure theories, criteria and models have been proposed to predict the ultimate strength of composite lamina/ and accurate modelling of matrix and fibre/matrix interface.The detailed stress or strain field generated can lead to accurate estimates of the onset and propagation of damage as well as the prediction of failure strength [9][10][11].Moreover, computational micromechanics modelling provides an alternative to the experimental approaches when it comes to complicated loading conditions such as biaxial or triaxial loads.Recently, micromechanics-based modelling was successfully applied to investigate the mechanical behaviour of composites under different combined loading conditions, such as combined transverse tension/compression and out-of-plane shear [12,13], combined transverse tension/compression and in-plane shear [13][14][15][16][17], and combined transverse compression and axial tension [18], combined longitudinal compression and in-plane shear [13].During the combined transverse compression and in-plane shear loadings, the friction between the fibres and matrix, after the fibre/matrix interface fails, plays an important role in the failure mode transition between interface-dominated failure in pure shear and matrix-dominated failure in pure compression [15].This leads to an increase of shear strength under a moderate transversal compressive stress state, which has been observed experimentally [19] and predicted analytically [20].Such effects under moderate compression cannot be predicted by computational micromechanics-based modelling if the fibre/matrix friction is omitted, as reported in [21,22].The increase in shear strength can also be captured by assuming arbitrarily large fibre/matrix interface shear strength [16].This, however, results in unrealistic predictions of the transverse tension and in-plane shear strengths [15].The combined effect of friction with cohesive behaviour has been addressed by other authors [23,24] when using traditional cohesive elements, which eventually led to the development of a cohesive element formulation that takes both mechanisms into account.
In recent years, artificial neural networks (ANNs), within the framework of machine learning, have been explored in the failure prediction of composite materials [21,25,26].An ANN consists of three kinds of layers: an input layer, hidden layer(s) and an output layer.ANN can discover intricate patterns from a large database by using a backpropagation algorithm to indicate how a network should change its internal parameters that are used to compute the representation in each layer from the representation in the previous layer to match the target [27].Thus, it is commonly used either to construct a surrogate model that can improve the computational efficiency, or to approximate a function which has no explicit mathematical form or is difficult to obtain.Yan et al. [28] employed ANNs to construct a surrogate model using regression for the constitutive law and classification for the damage information based on unit cell modelling.Failure analysis was performed on open-hole composite laminates under uniaxial tension with the ANN-based surrogate model and it was found that the proposed model offered substantive computational benefits over conventional FE models while maintaining sufficient accuracy.Liu et al. [29] proposed a new failure criterion for fibre tows within textile composites based on micromechanical modelling with a mechanics of structure genome (MSG) and a deep neural network model.Failure analysis of a plain weave composite example was performed with the proposed yarn failure criterion and other traditional failure criteria (i.e.Hashin, Tsai-Wu).It was found that the proposed failure criterion performs better than traditional ones with regards to accuracy.In addition, a neural network enhanced system, containing a subsystem with multiple neural networks was proposed to learn nonlinear inplane shear constitutive law and failure initiation criterion of CFRP composite lamina.Lee et al. [30] conducted biaxial tests on crossply CFRP composite tube under combined axial tension/compression and torsion with a range of biaxial stress ratios.An ANN was used to construct the failure criterion for the composite subjected to certain biaxial loadings and compared to the ones obtained from the Tsai-Wu failure criterion and a combined optimised tensor polynomial theory.It was found that the ANN performs best among these three methods, yielding the smallest root-mean-square error.In previous work [21], a high-fidelity micromechanics-based RVE model was utilised to generate the failure points under triaxial transverse, out-of-plane and in-plane shear loadings, and an ANN was used to construct the failure prediction of IM7/8552 unidirectional composites specifically for those loadings.This study extended the previous work to a more general failure prediction which covers various biaxial loading conditions and takes the failure probability into account.The more data the ANN receives from other loadings, the more powerful and accurate it becomes.
Computational micromechanics analysis is usually adopted to assess conventional failure criteria through a few biaxial loading conditions [12,13,[15][16][17].The objective of this study is to construct failure criteria for IM7/8552 UD composites under biaxial loadings via ANN models based on computational micromechanical modelling considering failure probability.The training database contains the failure data associated with failure probability from all possible biaxial loading conditions.Failure probability is obtained through coupled numericalexperimental data under transverse and in-plane shear, which is assumed to follow a Weibull distribution.In this study, computational micromechanics-based RVE models are developed to provide off-line training datasets of UD CFRP composite laminae under biaxial loadings, described in Section 2. The RVE model is built with three different phases (i.e.fibres, matrix and fibre/matrix interfaces).Periodic boundary conditions are used to apply biaxial loadings.The friction between the fibres and matrix, after interface failure, is enabled by implementing a cohesive frictional damage model via a VUMAT in ABAQUS.The fibre failure is predicted using the maximum stress failure criterion, and the effects of hydrostatic stress on the yielding and damage behaviour of the epoxy matrix is modelled with a Drucker-Prager plastic damage model.Five RVE models, with different arrangements of fibres, are used to compute the failure strengths under biaxial stress states.In Section 3, densification and pre-processing of data obtained from the five 3D high-fidelity RVE models are performed using failure probability with curve fitting.One database with failure probability is used for regression problems and another one without failure probability is used for classification problems.Section 4 illustrates the setup, training, validation and testing processes of ANNs for regression and classification problems.Results and discussion for both regression and classification problems of composite failure are given in Section 5.

3D micromechanics-based RVE models
RVEs have been shown to be a reliable way of predicting the elastic mechanical properties of the UD composite lamina.In our previous works [16,21] and that of others [9,13], it has been shown that 50 fibres are sufficient to capture the failure of unidirectional composites with a fibre volume fraction of 60%.By setting the average fibre radius to 3.5 μm, based on experimental measurements of IM7/8552 UD lamina [31], an RVE of 50 μm × 50 μm is obtained.It was reported in [32] that the depth of the RVE has an insignificant influence on the failure strength prediction under transverse loads.Therefore, an RVE with 5 μm depth was used for the numerical simulation under transverse loads [21], while an RVE with 50 μm depth was adopted to capture fibre failure, as shown in Fig. 1.The microstructure is generated by a discrete element method-based approach combining experimental data and an initial periodic shaking algorithm [33].This algorithm overcomes the jamming limit appeared in previous methods [34] and is capable of generating high volume fractions of fibres with random distributions and controllable inter-fibre distances by shaking an initial hexagonal packing of the fibres.Two-dimensional RVE model with periodic fibre distributions was generated and the final RVE was achieved by extruding the model along the fibre direction.
Five RVE models with random fibre distributions were created for the generation of failure data, in order to provide training and validation datasets for the ANN, see Fig. 1(a).These models were developed using ABAQUS/Explicit, in which the fibres and the matrix were discretised with hexahedral solid elements with reduced integration scheme (C3D8R) and wedge elements (C3D6).A convergence study was conducted to determine the size of mesh with a good compromise between result accuracy and computational cost.The interface was meshed using first-order cohesive elements (COH3D8), as shown in Fig. 1(b).A periodic mesh was used on the RVE to allow the application of periodic boundary conditions (PBCs).Mass scaling is sometimes applied to accelerate numerical simulations in ABAQUS/Explicit.A ratio of the kinetic energy over the internal energy of the system below 10% suggests that the mass scaling has insignificant effects on the numerical results [35].Consequently, the stable time increment was set to 5 × 10 −6 s in this study.

Constitutive model of fibres
The carbon fibres were modelled as linearly elastic, transversely isotropic solids.For the fibre failure cases where longitudinal tension and compression are involved, the maximum stress failure criterion was used to model the brittle fracture of fibres.The criterion was implemented into ABAQUS/Explicit via a VUMAT subroutine.The maximum stress failure criteria compare the measured strengths of a given material with the calculated stresses.Therefore, tensile and compressive failure of fibres can be predicted by: where   and   are the failure indices for longitudinal tensile and compression failure, respectively. 11 is the stress of fibre in the longitudinal direction, and the fibre fails immediately when the indices reach one.  and   are the uniaxial longitudinal tensile and compressive failure strengths of fibres, respectively.Under longitudinal compression, composite failure is triggered by the mechanism of fibre kinking, which is strongly influenced by elastic properties and geometry of the fibres [36] as well as manufacturing-induced void content [37] and other imperfections [38].A 3D RVE model containing 50 fibres resulted in 1.2 million DOFs with a fibre length of 500 μm and took over 72 h to complete a single simulation under uniaxial longitudinal compression using an HPC cluster and GPU acceleration [39].It is impractical to use such a model to generate failure points under biaxial loadings.
Although the single fibre unit cell RVE can predict a similar ultimate load as the one obtained from the multi-fibre RVE, it cannot capture kinking banding features [39].On the other hand, in the biaxial loading condition, two different failure modes compete at the transition point and one failure mode dominates when it is far from the transition point [15,16].It is still challenging to conduct the biaxial longitudinal related simulations with the fibre kinking feature taken into account.
When such a model is developed, failure data under biaxial longitudinal related loadings can be easily replaced for an improved ANN.Therefore, to keep the loading approach consistent with periodic boundary conditions, a simplified RVE model was adopted with degraded fibre strengths to represent the aforementioned uncertainties.It should be noted that although this RVE model cannot capture the kinking band, it can predict the failure strength under longitudinal loadings and is still capable of capturing other failure modes under biaxial loads.The initially degraded failure strength of fibres calibrated from experiments are 3.95 GPa and 1.74 GPa under longitudinal tension and compression, respectively, which are listed in Table 1.

Constitutive model of matrix
The polymeric matrix is generally modelled using an isotropic elasto-plastic criterion with damage.In order to capture matrix failure, the Mohr-Coulomb criterion [14], the extended Drucker-Prager yield model associated with a ductile damage criterion [9,16], the modified Drucker-Prager plastic damage model [15,21,42] and the paraboloidal elasto-plastic yielding criterion with an isotropic damage constitutive model [13,22] are commonly adopted.In this paper, the built-in Drucker-Prager (DP) plastic damage model in ABAQUS was adopted for predicting the mechanical performance of 8552 epoxy subjected to multiaxial stress states.The DP plastic damage model is based on the yield function proposed by Lubliner et al. [43] associated with modifications accounting for damage evolution subjected to tensile and compressive loads [44].The yield function is described in terms of a von Mises equivalent stress and hydrostatic pressure stress as [45]: wherein  is the pressure sensitivity factor of the yield function, ranging from 0 to 0.5,  is the von Mises equivalent stress,  is the hydrostatic pressure stress, σ is the maximum principal stress, ⟨⟩ represents the Macaulay brackets to return the argument if positive and zero otherwise. 0 and  0 represent biaxial and uniaxial compressive yield strengths, respectively.  and   are effective tensile and compressive yield stress, respectively.  is the ratio of the second stress invariant on the tensile meridian to that on the compressive meridian [46].The plasticity model assumes non-associated potential plastic flow.The plastic flow of the material is controlled by the flow potential , and is expressed as where  is a parameter referred to as the eccentricity,  0 is the uniaxial tensile stress at failure,  is the dilation angle in the  −  plane.The internal friction angle controls the dependence of the plastic behaviour of the material on hydrostatic pressure.After onset of damage under tensile loads, in order to ensure the correct energy dissipation of the matrix, an exponential cohesive law is adopted to control the quasibrittle behaviour of the material, which is characterised by a scalar damage variable [46].For the behaviour of matrix under compressive loads, perfect plastic behaviour is assumed based on experimental observations [47].More details about the Drucker-Prager plastic damage model and its numerical implementation can be found in [21,46].

Constitutive model of interface
The interface between fibres and matrix is modelled as cohesive elements governed by the bilinear traction-separation law [46].The thickness and stiffness of the fibre/matrix interface were identified in our previous work [21].The damage onset of this interface is determined by a quadratic interaction criterion, and damage propagation is controlled by the energy-based mix-mode Benzeggagh-Kenane fracture criterion [48].The interfacial strengths are calibrated from uniaxial experiments (i.e.tension, compression and shear) [41], while the fracture energies are set at 2 J∕m 2 [15] and 100 J∕m 2 [15,21] for Mode I and Mode II, respectively.It was experimentally observed that the shear strength increases in UD composite laminae under a moderate transverse compressive stress when subjected to biaxial loads [19], compared to pure shear stresses.This is difficult to be captured adequately with current built-in cohesive zone element in ABAQUS [46] since when modelling with the conventional cohesive element, the friction can only be considered when the cohesive element is totally damaged and removed from the finite element mesh.This phenomenon was also theoretically predicted by Puck's failure theory and numerically simulated with a cohesive surface together with a pure Coulomb model [15] or a cohesive-friction damage approach implemented into ABAQUS via a VUMAT subroutine [17,23,49].Here the cohesive-frictional model for the fibre/matrix interface was used and its main part was briefly reported below, and its implementation and one element verification can be found in Appendix.
The constitutive equation of the standard cohesive law is defined as where the    ,   and   are the traction, penalty stiffness and displacement for Mode I, respectively;    ,   and   (i=1,2) are the traction, stiffness and displacement for Mode II and Mode III, respectively.The mix-mode equivalent displacement jump  and mode-dependent penalty stiffness   can be defined as [49]: Following Turon et al. [24], a local mixed mode ratio  is defined to describe the evolution of the damage variable when subjected to mixed-mode loading, hence the ratio  reads: Here assuming   =  1 =  2 , the mode-dependent penalty stiffness   can be expressed as   = (1 − )  +   .The damage activation function  () for the interface damage initiation is defined [49]: where () is a function of equivalent displacement which increases monotonically and   is the threshold function.The interface damage is initiated when the function  () reaches zero.These functions are updated at every time increment  during the analysis where  indicates the actual time, and defined as: where  0 and   are the displacement jumps for the interface damage onset and final failure, respectively.The Benzeggagh and Kenane (BK) criterion [48] is adopted for the determination of these parameters: where  0  and    are the displacement jumps at damage onset and final failure in pure Mode I, respectively. 0  and    are the equivalent pure shear mode onset and final displacement jumps (  = √  2 1 +  2 2 ), respectively.In addition,  represents the BK law coefficient and here is taken as 1.45 [15,21] These parameters can be directly calculated from pure mode stiffnesses and critical energy release rates since the standard cohesive model is bilinear: where  0  and   represent the critical traction at damage onset and critical energy release rates for pure mode, respectively.
The damage variable  can be expressed as a function of the threshold function (  ), mix-mode onset ( 0 ) and failure displacement jumps (  ) [24], reads: Therefore, the tractions obtained from the cohesive damage model for mixed-mode loadings are listed as: The post-failure behaviour of cohesive element can be governed by the Coulomb friction law over the cracked region, in which the displacement is normally divided into two parts, namely an elastic part (  ) and inelastic (  ) part, where   (  =  −   ).The traction in the cracked part is described: The Coulomb friction function is introduced as: where  is the friction coefficient and the equivalent shear traction is expressed by the single shear mode as: The evolution of inelastic displacement jump    is defined by a non-associative relationship as [50]: with  the displacement jump and the multiplier δ is uniquely defined in the Kuhn-Tucker conditions: Finally, the cohesive damage law associated with the Coulomb friction for the post-failure behaviour is given by:

Periodic boundary and loading conditions
Periodic boundary conditions (PBCs) are normally used to ensure the periodic displacement field and traction continuity within the RVE once it is constructed with periodic mesh.In the commercial FE code ABAQUS/ Explicit [46], the PBCs are implemented by generating the relations between the nodes of surfaces and their counterparts on the opposite side.The unified PBCs equations read in terms of displacement vectors   ,   and   : where   ( = 1, 2, 3) and   ( = 1, 2, 3) are the length of the RVE and imposed displacement, respectively.Three dummy nodes on three principal axes are introduced to apply loads and origin of the RVE is fixed to avoid rigid body motion.The imposed strains were computed from the imposed displacements divided by the corresponding lengths, while the predicted normal and shear stresses were computed from the resultant normal and tangential forces acting on the RVE faces divided by the cross-sectional area.
Although some of the biaxial loadings conditions are unlikely to occur in practice, it is still worth obtaining all the possible failure points under biaxial loadings to ensure the generality of the proposed failure criteria.Meanwhile, the failure data obtained under biaxial loadings can lay the foundation for the failure data under triaxial or other multiaxial stress states.Fig. 2 presents an example of biaxial loading conditions and the combinations of biaxial loads in 3D stress states.Take biaxial longitudinal and transverse loads as an example, failure strengths of the RVE under longitudinal tension/compression and transverse tension/compression should be computed.However, due to the transversal isotropy of the cross section of RVEs, the transverse stresses in two directions (i.e. 2 and  3 ) and in-plane shear stresses (i.e. 6 =  12 ) and  5 =  13 are equivalent.Therefore, a total 15 stress combinations is reduced to 9 combinations, which are shown in Fig. 2(b).There are 5 to 20 different representative loading cases in each stress combination, depending on the combination with transversal and/or shear stresses.One cell in the matrix shown in Fig. 2(b) represents a biaxial stress combination, and same colours in the matrix represent equivalent loading conditions.

Probability-embedded failure prediction framework
Composite materials have wider scatter in their mechanical properties compared to homogeneous materials.Therefore, it is important to take the scattering of the properties as well as their average value into account when designing composite structures [30].Variation in fibre distribution, voids in the matrix and fibre waviness exist in composites even though the manufacturing process is strictly controlled.Such microstructural uncertainties introduced by the manufacturing process result in the stochastic composite strengths under multiaxial stress state [51].Imtiaz and Liu [52] proposed a framework to predict the failure envelopes of composite materials under biaxial loadings considering failure probability with plane stress assumption.They considered the scattering in strength parameters and quantified the failure probability in composite lamina.They proposed three probability quantification methods, namely: (a) a pure data driven approach, (b) a probability-embedded failure criteria with an assumption of self-similar failure surface, and (c) probability-embedded failure criteria with an assumption of self-similar failure surface and Weibull distribution.Here in this work, a framework is extended from plane stress problems into three-dimensional stress problems by using the third approach which quantifies the failure probability via Weibull distribution with an assumption of self-similar failure surface, thus taking a step further towards the data-driven failure prediction of composites under general loading conditions.Fig. 3 shows a schematic diagram for failure prediction of composites under biaxial stress states (  −   , ,  = 1, 2 … 6,  ≠ ).In the context of composite failure theory,  1 ,  2 …  6 represent the stresses in longitudinal, transverse and shear directions.Mean failure surface with average failure probability of 50% can be determined with curve fitting techniques.The failure surface with probability of    corresponds to maximum failure probability since it covers all possible experimental data, while the region within the failure surface with the probability of    represents safe stress states.The experimental data in biaxial stress plane can be divided into different small segments along angular coordinate in which each segment represents a specific range of loading conditions (i.e.different applied strain or stress ratios), such as region  in Fig. 3. Ideally, each segment should be considered independently so that a specific distribution of scattering in failure strength can be obtained.However, the independent consideration of scattering distribution for each segment may require a large number of biaxial/multiaxial experiments which is practically impossible and some biaxial/multiaxial loading conditions may not even be achieved experimentally, which should rely on virtual experimentation from numerical simulations.
Based on the assumption of self-similar failure surface, the mean failure surface with average probability (  = 50%) can be expanded or contracted with different failure probabilities.The failure criteria with mean failure probability reads: Thus failure surfaces with other failure probabilities can be determined with where (  ) is the scattering parameter, which is used to determine the relationship between the mean failure surface with the failure probability of 0.5 and other failure surfaces.Here the scattering parameter Z is defined as: where  is the distance from origin to the mean failure surface for a specific loading condition, and Here the mean failure surface is determined using curve fitting techniques, and then it is used to compute the (  ) for experimental and/or numerical data.
The scattering distribution of failure strength of composite materials is usually characterised by normal, log-normal or Weibull distribution [30,53].Here in this study, the normal and Weibull distributions are considered to describe the scattering of strengths for different loading conditions.Typical probability density functions (PDFs) are expressed as where   is the scattering parameter of th sample,  and  are the mean value and standard deviation of the normal distribution, respectively. and  are the shape and scale parameters in the probability density function of Weibull distribution.The cumulative distribution functions for both methods can be written in terms of the scattering parameter (  ) as where   is failure probability of   ,  is the shape parameter to describe the scattering degree and  0 is a scale parameter of the distribution, also termed as the mean value in the Weibull cumulative distribution function.The mean rank method is used to determine the value of   (  ) at   [52].The Benard's approximation for the median rank method reads where  is the total number of the samples.The shape parameter and mean value of the Weibull distribution can be determined from It should be noted that  <  0 means lower failure strength compared to the strength on the mean failure surface and vice versa.Thus, the number of failure cases increase as  increases.

Coupled experimental-numerical datasets with failure probability
A mesh convergence study was performed to ensure that the RVE models capture accurately the mechanical behaviour and corresponding progressive damage process of each constituent while reducing  computational costs.A good compromise should be made between the accuracy of numerical results and computational time via a comparison with experimental results.The RVE1 model from Fig. 1(b) was selected for the convergence study under various uniaxial loading conditions.Fig. 4 presents the comparison of simulation time, predicted elastic properties and failure strengths under various transverse and shear related uniaxial loadings from the convergence study.Three types of RVE models were generated with a different number of elements, ranging from 12 (low mesh density), 22 (medium mesh density) to 42 (high mesh density) elements, see Fig. 4(a).A minimum number of 12 elements was adopted due to the complexity of the RVE models.It can be seen in Fig. 4(c-d) that the numerical results obtained from the RVE models with three mesh densities are in good agreement with the experimental data.However, it is very challenging to accurately capture the stress distribution of the matrix in fibre-rich regions with the low mesh density.The computing time of the numerical simulations with low mesh density takes 1-1.5 h, compared to 5-7 h with medium mesh density and 27-45 h with high mesh density, see Fig. 4(b).Therefore, the RVE model with 22 elements was adopted for the rest of the simulations under biaxial loadings.Fig. 5 presents a comparison between numerical predictions and experimental results of a UD composite lamina under transverse and in-plane shear loading conditions.The microscale uncertainties, such as fibre distribution within the RVE model, were taken into consideration.The numerical results were obtained from the RVE models with five different fibre distributions, which are denoted by differently-shaped black data in Fig. 5.This numerical data can enrich the dataset obtained from expensive and time-consuming biaxial experiments [19], which is represented by the blue squares.A black line was fitted based on both experimental and numerical datasets, using a univariate spline method in Python due to the smooth surface it can generate with a limited number of datasets, and is defined as the mean failure surface.
Based on the definition of the scattering parameter,   was computed for each test in a total of 116 experimental and numerical failure cases.Numerical and experimental data within similar scattering parameter ranges were grouped together.Then   (  ) for scattering data were computed using the median rank method, and the shape parameter and mean value of Weibull distribution can be determined.Fig. 6 shows how the two parameters of the Weibull distribution were determined.A comparison between the fitted Weibull distribution   function and the Bernard approximation is also shown. and  0 are calculated to be 17.771 and 1.029, respectively.The failure probability curve obtained through the Weibull distribution function agrees well with the failure probability from the Bernard approximation, Fig. 6(b).In addition, it was found that the two-parameter Weibull distribution function is sufficient to characterise the failure probability of composite strength under biaxial stress states as it can be represented by a straight line, Fig. 6(a).Recently, the lower and upper bounds on the probability of failure were calculated by means of random set theory, in which an efficient probability-based reliability method, called subset simulation, was employed [54].The lower and upper bounds were used in the failure analysis of composite in [16], in which the bounds were determined with the lowest and highest biaxial failure strengths obtained from numerical simulations to assess the conventional failure criteria.In this study, due to the limitation of experimental and numerical data, a 95% confidence interval was considered to describe the failure probability, in which the lower bound    and the upper bound    are set to 0.025 and 0.975, respectively.

Database generation for regression problem
Due to the transverse isotropy of the cross section of the UD composite lamina, a total of 15 biaxial stress combinations can be reduced to 9 combinations, shown in Fig. 2. To rationalise the computational cost of micromechanics-based high fidelity modelling, the first RVE model was selected for the rest of the loading cases to determine the mean failure surface.In addition, because of the lack of experimental data of the other biaxial loading cases, an assumption was made that the strengths share the same scattering pattern with the combined transverse and in-plane shear loading as shown in Fig. 5. Fig. 7 presents the failure envelopes of an IM7/8552 UD composite lamina under different biaxial loadings with failure probabilities.Failure points were collected from biaxial loading cases to determine the mean failure surface.Here 2.5% lower and 97.5% upper bounds were selected for the demonstration considering a 95% confidence interval.It should be noted that it is still a matter of debate whether the failure envelope for composite lamina under biaxial transverse compression should be open (i.e.no failure) or not [8], therefore from the perspective of numerical modelling, an extra artificial failure point was added to close the failure envelope, see Fig. 7(d).The resolution of the database depends on the interval selection of failure strength on the failure surface and the failure probability of the surface.In this study, 1% interval of failure probability between the lower and upper bounds and 300 data points on the failure envelope with failure probabilities in each stress combination were selected, resulting in a database containing half a million samples.

Database generation for classification problem
Material failure can also be treated as a classification problem.With a proper failure criterion defined, the state of safety or failure can be classified.In most engineering applications, a credibility variable is commonly used as a safety factor to measure the confidence of the numerical simulation.One of the advantages of ANN is that the output of a classification problem (i.e.accuracy) is actually a probability of the categorisation based on the failure criterion.In this section, a failure classification problem was introduced based on the previously generated database for regression problems.The mean failure surface with the failure probability of 50% was used as the failure envelope in the classification problem.Thus, the region, where points are inside the failure surfaces with probabilities between 2.75% and 50%, is regarded as a safe region and their probabilities were replaced with 0; while the region inside the failure surfaces in which points with probabilities between 50% and 97.5%, is regarded as a failure region and their probabilities were replaced with 1. Then a new database was generated for the failure classification problem.The failure envelopes, and safe and failure regions for IM7/8552 composite lamina in different biaxial loading conditions is shown in Fig. 8.

Artificial neural networks
An artificial neural network (ANN) is a machine learning algorithm which mimics the neurological processes of a human brain, and is used to tackle various problems in industry, such as image recognition, machine translation, novel material discovery, and more recently, failure prediction of materials.ANN generally consists of three basic layers: an input layer, hidden layers, and an output layer.The input layer, where data is entered, and the output layer are connected through the hidden layer in a neural network.Neurons in each layer are fully or partially connected to the neuron in preceding and subsequent layer with unique weights.The input signal propagates through the hidden layers on a layer-by-layer basis towards the output layer.The backpropagation training algorithm [55] is usually adopted to determine the weights (  ) to iteratively minimise a chosen cost function.The training process is terminated when the loss function (i.e.mean-squareerror (MSE), root-mean-square-error (RMSE) or mean-absolute-error (MAE)) reaches a pre-specified threshold, or a pre-specified number of learning epochs is completed.
These weights between neighbouring neurons help to determine the importance of any given variable, with larger weights contributing more significantly to the output.A bias term   is used to cover a wider range.All inputs from the last layer are then multiplied by their respective weights and then summed to the bias, so the input   obtained by a given neuron , where  is the number of neurons at the current layer.An activation function is used to compute a scalar value for each neuron with   in the current layer to decide whether the information carried by this neuron should pass to the next layer.There are four commonly used activation functions, namely Tanh, ReLu, Sigmoid and Softplus.In this study, a ReLu activation function is used, which has been proved to have a better training performance for deep neural networks compared to others [56].During the training process, the weights and bias of the network are updated to minimise the loss function.The loss function measures the difference between the true and predicted values.The loss function can be defined as the commonly used mean square error (MSE) and mean absolute error (MAE), in Eq. ( 29).The MSE and MAE represent the average of the squared difference and the absolute difference between the original and predicted values in the dataset, which measures the variance of the residuals and the average of the residuals in the dataset, respectively.In most regression problems, the MSE was used for the training and validation of the ANNs [28], here the MAE was used for the evaluation of the network; while the binary cross-entropy loss function (Eq.( 30)) was used for the classification problems.
where  is the number of training samples,   is the true output in the training datasets while ŷ is the prediction of output using the ANN architecture.
where   and   are the score and the ground truth label for the class   .Since the classification problem only has two classes, so  = 2 and  2 = 1 −  1 ,  2 = 1 −  1 .In general, the training process can be considered as a typical application of a gradient-based optimisation algorithm and statistical estimation in a back-propagation manner.In the current study, the Adam [57], an algorithm for first-order gradientbased optimisation of stochastic objective functions, was used based on adaptive estimates of lower-order moments for both regression and classification problems.In this optimisation process, the weights and bias of all the network are updated iteratively until the desired error tolerance is met or the maximum number of iterations (epoch) is reached.Due to the paired input and output data, both training for regression and classification problems are categorised as supervised learning.
Under biaxial loading, there are a total of 15 possible stress combinations, or biaxial loading conditions, as shown in Fig. 2(b), where one grid entry represents one biaxial loading condition.The regression problems of composite failure can be dealt with using two strategies.In the first strategy, only one ANN is used to cover all 15 biaxial loadings.To more accurately predict the failure probability under biaxial loadings, and to consider the failure mode in each loading, the second strategy uses 15 ANN models which are used for the prediction of failure probability under 15 different biaxial loading conditions, Fig. 9.
In the context of neural networks, the probability embedded failure criterion for the regression problem is the mapping of failure probability of composite materials from different stress state combinations.The ANN can provide a form-free criterion for the mathematical form defined in Eq. (22).A database, containing 505,000 samples generated from 3D high-fidelity micromechanical failure analysis and biaxial experimental data, was used to feed the ANNs.Each sample accepts six stress components as input and provide a failure probability as output.The database consists of three parts: training data, validation data and test data, and was split in a 80:10:10 ratio.The validation data is used for the optimisation of hyperparameters of the ANN, while the test data is to test the performance of the ANN with unseen data.
For the regression problem, in the first strategy, the database, with the 15 loading cases, was split as a whole in the 80:10:10 ratio for training, validation and test, respectively.While in the second strategy, each dataset of the database were split in the same ratio for the training of each ANN separately.For the classification problem, the database transformed from the regression problem was also split in the 80:10:10 ratio for the training.The focus of this study was on the biaxial loading conditions, and more data from the micromechanical failure analysis under triaxial loading conditions will be generated in future to enrich the database.

Probability embedded failure criteria with one ANN for regression problems
The ANN structure and hyperparameters (i.e.batch size and number of epochs) have significant influences on its performance, thus manual parametric studies for the hyperparameters optimisation were conducted on the number of hidden layers and neurons as well as the batch size and the number of epochs.Fig. 10 shows the comparison of MSE and MAE in training and validation processes, and predicted and true values with different ANN structures.With three hidden layers, the MSE decreases as the number of neurons increases.A better prediction can be obtained with the ANN model with 128 neurons in each hidden layer compared to the ones with 32 and 64 neurons.In addition, when one more hidden layer was added to the ANN model with 128 neurons, the MSE obtained in validation process decreases by 26.5% and the prediction capability improves as well.Therefore, the six-layer ANN model with 128 neurons in each hidden layer was used for the parametric study of hyperparameters.The batch size is a hyperparameter of gradient descent that controls the number of training samples to update the model's internal parameters, while the number of epochs is another hyperparameter of gradient descent that controls the number of complete passes through the entire training dataset.Fig. 11 shows the comparison of MSE and MAE in training and validation processes, and predicted versus true values with different batch sizes.It can be seen that when the batch size is 500, the predicted values obtained from the ANN model are underestimated and the MAE has large fluctuations.When the batch size increases to 1000, better performance of the ANN model is achieved.However, the prediction was not improved significantly when the batch size increases to 2000 with the predicted values being slightly overestimated.This is probably because the structure of the current ANN model has limited capability to take 2000 sample at a time.Therefore, the selected ANN model has four hidden layers, and 128 neurons in each layer.With the optimal batch size, the MSE and MAE obtained from predictions can reach 0.094% and 1.6%, respectively, and the predicted values have excellent agreement with the true values.

Probability embedded failure criteria with multiple ANNs for regression problems
It is found in Figs.10-11 that the MSE, as an indicator of accuracy in the ANN prediction, can reach 0.094%, however, due to the structure of the datasets, this error cannot be improved significantly by adjusting the ANN structure or hyperparameter.When multiple ANNs are used, corresponding to specific loading conditions, this error can be significantly improved to 0.027%, which ensures the accuracy of the prediction in biaxial loadings.An ANN model, containing 2 hidden layers with 64 neurons in each layer, was used for the construction of ANNs in the second strategy, which consists of 15 ANNs.Different datasets were fed to the ANN model to construct failure criterion for each loading case with failure probabilities.The training, validation and test processes of the ANN model for the biaxial longitudinal and transverse loading combination is selected for demonstration, which is shown in Fig. 12. Fast convergence of MAE and MSE can be observed in training and validation processes.Excellent agreement is found between the predicted and true values, with a MSE value of 0.027% and a MAE value of 0.75%, which suggests an accurate and robust ANN model.

Probability embedded failure criteria with one ANN for classification problems
For the classification problem, a four-layer ANN was initially used for training, in which the hidden layer has 6 neurons.This is similar to the ANN architecture used in our previous work [21] for training the failure points under triaxial loading conditions.However, it can be found in Fig. 13 that only 77.9% prediction accuracy can be achieved.When the number of the neurons in hidden layers increase from 6 to 16, the prediction accuracy of the ANN models increase by 21%.Further increases of neurons do not contribute significantly to the improvement of the prediction accuracy.Comparing the loss and accuracy obtained from training and validation processes, the ANN model with 32 neurons in the hidden layers converges faster than the one with only 16 neurons.In addition, it can be seen from the ANN model with three hidden     layers that one more hidden layer added to the ANN model can only improve the ANN convergence without much improvement in prediction accuracy.The computing time for a whole training, validation and testing cycle obtained from all ANN models was less than 4 min on a workstation with Intel(R) Xeon(R) Gold 6226R processor containing 16 cores.stresses were generated within the lower and upper bounds following Weibull distribution, see Fig. 5. Regarding the regression problem of composite failure, the overall ANN, aimed for all loading conditions, with the structure of 6-128-128-128-128-1 and a single ANN, aimed specifically for transverse and in-plane shear loading, with the structure of 6-64-64-1, were selected for the failure prediction of composites, which can be found in Fig. 14(a) and (b), respectively.The colour of dots represents different failure probability, ranging from 2.5%-97.5%,at the biaxial stress state, and lighter dot means the composite is highly likely to fail at this stress state.The black line is fitted mean failure envelope, representing the failure probability of 50%, obtained from coupled numerical and experimental data, see Fig. 5.It can be found that few lighter dots predicted by the overall ANN fall within the black line, while the single ANN performs better with the comparison of the fitted envelope, due to its smaller MSE.Both of failure criteria cannot capture the increase of shear strength under moderate compressive stresses in the transverse compression and in-plane shear quadrant, and Hashin failure criterion overestimates the biaxial failure strength when shear stress is dominant in the transverse tension and in-plane shear quadrant.

Comparison between biaxial failure strengths predicted by trained ANN models and failure criteria
Regarding the classification problem of composite failure, a shallow ANN with the structure of 6-32-32-32-1 was selected for the failure prediction of composites under combined transverse and in-plane shear loadings.The Fig. 14(c) represents the comparison of failure points predicted by the ANN and failure envelopes, in which the red star represents the failure point while the black dot represents the safe point.Overall, it can be clearly seen that the ANN performs well compared to the fitted mean failure envelope.Still, 3 black dots, which represent safe stress states, fall beyond the fitted envelope due to the prediction accuracy of 98.1%.In addition, comparison between the failure prediction from the ANN for the classification problems and the overall ANN for the regression problem where the threshold of the failure probability is set to 50%, are conducted in Fig. 14(d).Red and black colours represent failure and safe states, respectively; circle and star represent the ANN for classification problems and the overall ANN for regression problems, respectively.It is found that both predictions have excellent accuracy regarding the fitted envelope, while few wrong prediction can also be found in both predictions, this discrepancy mainly comes from the fitting nature of artificial neural network.

Conclusions and future work
This study explored the application of artificial neural networks (ANNs) in the failure prediction of unidirectional carbon fibre reinforced polymer composites based on micromechanics modelling.3D high-fidelity representative volume element (RVE) models were set up to provide the failure points under biaxial loading conditions.Five RVE models with different fibre distributions were used to provide failure points under transverse and in-plane shear loadings for the determination of the mean failure surface.Two ABAQUS VUMAT subroutines were developed to predict fibre failure and capture the friction under transverse compression and in-plane shear loadings.A coupled experimental-numerical framework was applied to conduct uncertainty quantification and generate datasets for transverse and in-plane shear loadings with probabilities.A total of 15 biaxial stress combinations was considered for the generation of the datasets.In each combination, a minimum of 10 representative biaxial loading cases were computed to provide the failure points.Three data-driven probability embedded strategies were considered to construct the failure criteria for unidirectional CFRP composites in terms of regression and classification problems.Parametric studies were conducted for the selection of ANN models.
It was found that all three strategies have excellent performance in the probabilistic failure prediction of composite materials.For the failure regression problems, in the first strategy, the selected ANN model with 4 hidden layers, and 128 neurons per layer had a mean square error (MSE) of 0.094% and a mean absolute error (MAE) of 1.6%; while in the second strategy, the single ANN model within the framework achieved an MSE of 0.027% and an MAE of 0.75%.Both strategies are consequently considered to be promising in the probabilistic failure prediction of UD CFRP composites.For the failure classification problems, a shallow ANN model with only 3 hidden layers, in which 32 neurons were used, demonstrated excellent performance in the failure prediction of composites under biaxial loadings, with a prediction accuracy of 98.1%.The comparison between the failure strength of composite materials predicted by these three ANNs and failure envelopes predicted by Tsai-Wu and Hashin failure criteria was conducted and it was found that all three networks have good performance regarding the comparison with the fitted envelope.Both of failure criteria cannot capture the increase of shear strength under moderate compressive stresses in the transverse compression and inplane shear quadrant, and Hashin failure criterion overestimates the biaxial failure strength when shear stress is dominant in the transverse tension and in-plane shear quadrant.The computing time for a single numerical simulation of biaxial loadings took 6-11 h depending on different biaxial stress combinations, while a complete ANN training, validation and testing cycle took less than 4 min on a workstation with Intel(R) Xeon(R) Gold 6226R processor containing 16 cores.
The novel aspects of this study include: • The introduction of failure probability into the database used to train ANN models for the failure prediction of IM7/8552 UD composites under biaxial loading conditions.• The comparison of failure prediction between ANN models for regression and classification problems and conventional Hashin and Tsai-Wu failure criteria.• The completeness of failure analysis of IM7/8552 UD composites under all possible biaxial loading conditions using coupled computational micromechanics and artificial neural networks.
It should be noted that this ANN-based failure approach for failure prediction based on 3D high-fidelity micromechanics modelling is demonstrated for IM7/8552 unidirectional composite materials.Extra effort should be made for the failure prediction of composites under biaxial loading conditions if a new material system is introduced, resulting in the necessity of the generation of general failure criteria.There are two limitations of this study.On the one hand, the fibre failure was simulated with a maximum stress failure criterion, neglecting fibre kinking mechanisms in longitudinal compression and stochastic fibre strength in longitudinal tension, therefore improved RVE models need to be considered to capture such phenomena.On the other hand, the parameters of Weibull distributions used for the rest of the biaxial loading conditions are assumed to be the same as the ones calculated from the coupled experimental-numerical data under transverse and in-plane shear loadings due to the limitation of experimental data.However, this is not always true due to the uncertainties in composites, resulting in the variety of Weibull parameters for different loadings.This can be improved by advanced experimental approaches to generate more specific data or advanced numerical approaches with uncertainties considered.Currently, the machine learning-based failure criteria are being implemented into ABAQUS via a VUMAT subroutine to conduct the progressive failure analysis of unidirectional composites with an open hole and woven composites under complex loadings.

Fig. 1 .
Fig. 1. 3D geometrical RVE models with three phases: (a) RVE models with different fibre distributions, (b) RVE1 with 5 μm depth for inter-fibre failure cases and RVE6 with 50 μm depth for fibre failure cases.

Fig. 2 .
Fig. 2. (a) Biaxial loading conditions imposed on the RVE with  1 and  2 as an example (solid arrows represent tension and dash lines represent compression); (b) Biaxial stress combinations under 3D stress states (Here  1 represents longitudinal tensile/compressive stress,  2 and  3 represent transversal tensile/compressive stresses,  6 ( 12 ) and  5 ( 13 ) represent in-plane shear stresses and  4 ( 23 ) represents out-of-plane shear stress, same colour in different grids represent equivalent loading conditions, while different colours represent different biaxial loadings).

Fig. 4 .
Fig. 4. (a) Comparison between 3D RVE models with different numbers of elements: low mesh density with 12, medium mesh density with 22 and high mesh density with 42 elements; (b) Comparison of computational time obtained from the three types of mesh density; Comparisons of (c) predicted elastic properties and (d) failure strengths with experimental data [40,41] under various uniaxial loadings.

L
.Wan et al.

Fig. 5 .
Fig. 5. Comparison between numerical predictions and experimental data [19], and the upper and lower bounds of normal and Weibull distributions.

Fig. 6 .
Fig. 6.(a) Determination of shape parameter and mean value of the Weibull distribution with median rank method, (b) Comparison between the Bernard approximation and fitted Weibull distribution.

L
.Wan et al.

Fig. 7 .
Fig. 7. Failure envelopes of IM7/8552 UD composite lamina with failure probabilities ( 1 is longitudinal stress,  2 and  3 are transverse stresses,  6 and  5 are in-plane shear stresses and  4 is the transverse shear stress).

Fig. 8 .
Fig. 8. Failure envelopes, and safe and failure regions of IM7/8552 UD composite lamina under different biaxial loading conditions.

Fig. 9 .
Fig. 9.The computational framework containing a set of ANNs for the construction of failure criterion.

L
.Wan et al.

Fig. 10 .
Fig. 10.Comparison of MSE, MAE and prediction obtained from the ANN with different structures.(The batch size and number of epochs used here are 1000 and 300, respectively).

Fig. 11 .
Fig. 11.Comparison of MSE, MAE and prediction from the ANN with different batch sizes.

Fig. 14
Fig. 14 presents comparisons between the failure prediction of composites under transverse and in-plane shear loadings with trained ANNs and Tsai-Wu and Hashin failure criteria.100 random biaxial

Fig. 14 .
Fig. 14.Comparison between failure strengths of composites predicted by trained ANN models, and Tsai-Wu and Hashin failure criteria under biaxial transverse and in-plane shear loadings: (a) overall ANN for regression problem; (b) single ANN for regression problem; (c) overall ANN for classification problem; (d) overall ANN for regression problem and overall ANN for classification problem when set the failure probability of 50% as the threshold.

Fig. A. 1 .
Fig. A.1.Cohesive-frictional damage model for the fibre/matrix interface: (a) Representative Elementary Area and decomposition of displacement jumps [23]; (b) traction-separation laws of standard cohesive zone model (black lines) and cohesive-frictional model (blue lines) ( 0  ,  0  and  0  ,  0  represent the displacements and tractions at the damage onset of standard cohesive and cohesive-frictional damage models, respectively); (c) graphical representation of standard cohesive damage model (black lines) and cohesive-frictional damage model (blue lines) under mixed mode loadings when compression is involved.