Exploiting the benefits of multi-scale analysis in reliability analysis for composite structures

This paper investigates two critical issues, namely propagation of multi-scale uncertainty, and selection of failure criteria, related to reliability analysis of composites by using multi-scale methods. Due to the multi-scale architecture of composites, uncertainties exist in both microscale and macroscale parameters. It is necessary, therefore, to consider random variables at various length scales to ensure accurate estimates of the reliability of composites. Three types of homogenization methods, namely rule of mixtures, Mori–Tanaka and computational homogenization, are adopted to link these two scales, and to propagate uncertainty from micro to macro scales. By integrating these homogenization methods with the stochastic ﬁnite element method and structural reliability methods, the reliability of composites can be investigated with a limit state function based on a chosen failure criterion. This multi-scale reliability analysis procedure has been applied to analyse laminated ﬁbre reinforced composites made of AS4/3501 carbon/ epoxy. Firstly, a comparative study has been conducted to evaluate the performance of the assumed homogenization methods for the reliability of composites, and to identify advantages compared with a single scale analysis. The results show that multi-scale analysis can provide more accurate reliability estimates. Secondly, several popularly used failure criteria for composites have been compared using multi-scale reliability analysis. (cid:1) 2016 The Authors. Published by Elsevier Ltd.


Introduction
The trend towards the increasing use of composites is being seen in diverse industries including aerospace, automotive, marine and construction.Advances in design and structural analysis approaches are eagerly required to fully exploit the benefits brought by these materials as other new construction materials [1,2].Typical fibre reinforced polymer composite structures are made up of laminates which in turn are obtained by stacking individual plies with different fibre orientation.This leads to three different entities including ply, laminate, and component, whose mechanical behaviour is characterised by three different scales, namely fibre diameter, ply and laminate thickness, respectively.Design performance is often related to the probability of failure and requires an understanding of the interaction of uncertain characteristics at both material and structure levels.It is for these rea-sons that a formal probabilistic analysis is proposed to study the performance of composite materials, components, and structures.
Probability-based methods are powerful tools to design structures under uncertainty.In structural engineering, reliability is one of the most used indicators that interprets response information for design, maintenance, repair, etc.In recent decades, a large number of articles have contributed to the understanding of the probabilistic failure and reliability of composites.For example, [3][4][5].However, most of the previous research consider uncertainties at ply-or/and structural-level parameters only, while uncertainties at micro-level parameters are ignored.Owing to the multi-scale architecture of composite materials, variability in structural responses is affected by the variations in parameters at various length scales.As has long been recognised, laminate stress analysis and lamina failure criterion are two critical elements in failure analysis of laminated composite structures.Multi-scale modelling methods are ideal tools to link micro-scale parameters with macro-scale parameters and to propagate uncertainties from micro-scale to macro-scale [6,7].They have demonstrated their capability to provide sufficiently accurate structural performance simulations due to the fact that it does not require any assumption on the constitutive model at ply level.For instance, unidirectional fibre reinforced composites are commonly assumed to be transversely isotropic materials.Consequently, estimation of reliability can be improved using multi-scale method [8].Several researchers have attempted to conduct multi-scale reliability analyses that consider uncertainties at both micro-scale and macro-scale [9][10][11][12].Concerning the second issue, despite years of extensive research around the world, a complete and validated methodology for predicting the behaviour of composite structures including the effects of damage has not yet been fully achieved with the exception of a few recent research contributions [13,14].A brief review of failure criteria for composites is provided in the next section.Although qualitative evaluations [15][16][17], quantitatively experimental comparisons [18], and numerical comparisons [19] have been made for deterministic failure criteria, quantitative comparisons of their performance considering uncertainty are much less prevalent [20].
The objective of this paper is to exploit benefits provided by the adoption of multi-scale analysis for the reliability analysis of composites.Different homogenization methods including rule of mixtures, Mori-Tanaka and computational homogenization have been adopted to link micro-scale parameters with macro-scale parameters and propagated uncertainties from micro to macro.These methods have then been integrated with stochastic finite element method and structural reliability method to conduct reliability analysis for composites.Benefits of multi-scale analysis have been investigated by a series of comparative studies on the effective elastic properties and their statistics, and reliability estimates.In addition, a comparative study of some frequently used failure criteria has been performed from a structural reliability analysis perspective.This has been conducted by using computational homogenization method based multi-scale reliability analysis.Using this approach for the reliability analysis enables us to investigate the influences of variations in micro-scale and macroscale parameters on the reliability estimation, and to quantify the relative importance of various uncertainties.Consequently, the sources of differences among various failure criteria can be identified in a new and unique way.In comparison with a deterministic failure analysis, a reliability analysis requires to conduct multiple structural analyses.For instance, several iterations are required to obtain a sufficiently accurate estimate of reliability when using the first order reliability method.In practice, finite element methods have become standard tools for numerical simulation of structural response.Hence, failure criteria incorporated in a finite element-based reliability analysis should be efficient in terms of implementation.In accordance with these overarching considerations, six failure criteria were finally selected in the present study: (1) maximum stress, (2) Hashin failure theory, (3) Tsai-Hill failure theory, (4) Tsai-Wu theory, (5) Christensen failure theory, and (6) Hoffman failure theory.Numerical analyses for a lamina and a quasi-isotropic laminate ð0 = AE 45 =90 Þ s made of AS4/3501 carbon/epoxy were adopted for the comparative studies.

Brief literature review on failure criteria for unidirectional fibre reinforced polymer composites
Composite materials display a wide variety of failure mechanisms including fibre failure, matrix cracking, buckling at several scales, and delamination as a result of their complex structural behaviour.Hence, various limit state functions (LSF) can be established depending on the specific problem under consideration.In this study, we limit our focus to laminate failure under in-plane loading conditions, which are largely based on the stress components of an individual ply within the laminate.Here, a brief review of theories for in-plane failure is presented.Readers are referred to find more comprehensive reviews, such as the state-of-the art composite failure theories are included in [17].Generally, failure criteria can be broadly classified into two groups according to whether failure modes are separated or not.

Failure theories without failure modes
Firstly, we focus on failure criteria where the ply failure modes are not considered, but the failure of the entire ply is predicted.From the application point of view, this group is easier and more convenient to use, especially in reliability analysis which involves iterative calculations.However, they are often criticised due to their lack of phenomenological basis.This group includes criteria from papers in which the difference between fibre and matrix failure is either unclear or not specified and the so-called fully interactive criteria.Hence, a single quadratic or higher order polynomial equation containing all stress (or strain) components is used to predict the failure.Tsai-Hill failure theory [21] is one of the criteria in this group that modifies the Hill's anisotropic failure theory derived from the Von Mises yield criterion for metals [22].To account for different strengths in tension and compression, Hoffman added linear terms to Hill's equation, defining the Hoffman failure criterion [23].Tsai and Wu further developed these criteria to improve their performance in the representation of experimental data that results in the well-known curve fitting based Tsai-Wu failure criterion [24].As one of the 19 leading failure criteria used in the worldwide failure exercises, Tsai-Wu failure criterion presented in [25,26] shows good performance [18].

Failure theories with consideration of failure modes
Hashin and Rotem [27], for the first time, proposed failure of laminated composites attributed to different physical phenomena including fibre controlled and matrix controlled failure modes.Many failure criteria have been developed following this idea, and new failure mode based criteria, which are also called physicallybased failure criteria in the literature, are increasingly being proposed.We only mentioned few of them in this brief review.

Fibre failure
Fibre failure occurs due to the accumulation of individual fibre failures within plies, which become critical when there are insufficient intact fibres remaining to carry the required loads.For fibre tensile failure, it is generally acknowledged that the maximum stress or maximum strain theory is sufficient to predict this failure mechanism, such as stated in [28].To consider the interaction between different stress components, more sophisticated models, have been developed.Hashin [29] used a quadratic interaction criterion involving in-plane shear.Chang and Chang [30] improved Hashin's criterion [29] by incorporating nonlinear shear behaviour.Puck and Schürmann [31] modified the maximum strain criterion to include the transverse normal stress through a stress magnification factor.Research in fibre compressive failure is still being performed as the failure mechanism is complex.Depending on the materials, microbuckling and kinking are often recognised as the main failure mode under compressive loading in longitudinal direction [32].A number of approaches have been developed for incorporating the effects of microbuckling and kinking, e.g.[31,33,34].

Matrix failure
Due to the unique feature of fibre-reinforced composites offering relatively lower properties in the transverse direction compared with the high strength and stiffness properties in the longitudinal direction, matrix cracks are usually the first observed form of damage in fibre reinforced composites [35].The matrix microcracks usually initiate at defects or fibre-matrix interfaces, accumulate throughout the laminate, and coalesce, leading to failure across a critical failure plane.Hence, numerous efforts have been dedicated to analysing these complex phenomenon in matrix failure, and a number of approaches have been developed to predict the formation and growth of matrix cracks.As fibre failure, the matrix failure mode is divided into failure in tension and failure in compression.
The classical maximum stress and maximum strain criteria are the basic types of rules used to examine the initiation of matrix failure.Due to the fact that matrix failure in tension generally involves an interaction between the tensile normal and in-plane shear stresses, many authors proposed failure theories to consider these phenomena by assuming a critical fracture plane in the transverse tension direction.Hashin and Rotem [27] proposed the first quadratic interaction model under in-plane loading condition, and this model was further developed in several studies, e.g.[29,36], to include nonlinear shear terms, in situ transverse tensile and shear strengths, crack density, etc.In contrast, Cuntze and Freund [37] proposed a different model that involves the transverse tensile stress and strength and through-thickness shear stress.
For matrix failure in compression, Hashin and Rotem [27] assumed the fracture plane was in the transverse direction and proposed a simple quadratic interaction criterion using the transverse normal and in-plane shear components.This was then modified by Hashin [29] to include the through-thickness strength and Chang et al. [36] by incorporating a nonlinear shear formulation.The criterion proposed by Cuntze and Freund [37] uses only the transverse normal strength, with a combination of several stress invariants.For the criteria considering a non-zero fracture plane angle, this angle must be either assumed or determined by checking all possible angles, though Puck and Schürmann [31] proposed an analytical form for the case of plane stress.Davila et al. [38] developed a set of phenomenological failure criteria denoted as LaRC03 based on the ideas of Hashin and Puck, where no empirical and non-physical fitting parameters are required.Subsequently, improvements on LaRC03 criterion were conducted to extend it to three-dimensional stress states, and these led to the so-called LaRC04 criterion [39].To address the general three-dimensional stress states confronted in the design and analysis of thick composite laminates, 3D failure criteria based on Puck's ideas were developed in [40], and it was subsequently improved to distinct transverse and longitudinal failure mechanisms and that led to a new 3D failure criterion [41].

Reliability formulation for unidirectional fibre reinforced composites based on various failure criteria
Let us consider an unidirectional fibre reinforced composite lamina as shown in Fig. 1, a coordinate system (1-2-3) with the 1-axis aligned parallel to the axis of the fibres is used, and the 2-3 plane is isotropy.Therefore, the stress term such as r 11 represents the stress in the longitudinal direction.X; Y and Z are defined as the strengths in the longitudinal, transversal and lateral (or through thickness) directions.T and C are used to distinguish tensile strength and compressive strength, and subscripts are used to identify their direction.X ¼ T 11 ; Y ¼ T 22 , and Z ¼ T 33 if r 11 P 0; r 22 P 0, and r 33 P 0, respectively, and and Z ¼ C 33 , in the same order.S is used to represent shear strength with S 12 ; S 13 and S 23 being in-plane shear strength, transverse shear strength.

Maximum stress failure theory (MS)
The lamina is considered to have failed if Based on the distortion theory and von Mises' yield criterion, Hill [22] proposed a yield criterion for orthotropic materials The components F; G; H; L; M, and N of the strength criterion depend on the material strengths X; Y and S, and they can be expressed as [21]: For a unidirectional lamina in Fig. 1 with isotropy plane in 2-3, Eq. ( 6) can be simplified, and LSF derived from it can be written as

Hoffman failure criterion (HO)
To account for different strengths in tension and compression, Hoffman [23] added linear terms to Hill's failure criterion in Eq. ( 6) that leads to: Similarly, the failure parameters a ij can be related to the material strengths as ; With T 33 ¼ T 22 ; C 33 ¼ C 22 , and S 12 ¼ S 31 for the concerned lamina (Fig. 1), the failure criterion for a unidirectional lamina is simplified to Therefore, the LSF based on Hoffman failure theory can be established as

Tsai-Wu failure criterion (TW)
To improve the correlation between failure theory and experiment, Tsai and Wu [24] proposed a new interaction equation that uses various strengths in tensor form.This failure criterion is a specialization of the general quadratic failure criterion and expressed in the form where i; j ¼ 1; Á Á Á ; 6 and repeated indices indicate summation, and F i ; F ij are experimentally determined material parameters.The stresses r i are expressed in Voigt notation.
For the studied unidirectional lamina with the 2-3 plane as the isotropy plane, the following equalities hold.Then the general form of the Tsai-Wu failure criterion reduces to where Most of the parameters, F i and F ij , can be directly determined from the uniaxial strengths of laminate using However, the remaining parameters, F 12 and F 23 , are the interactive strength terms.These must be determined experimentally using combined stress tests in order to apply the theory to either the full three-dimensional or the reduced plane stress cases.However, very few materials have been characterized for these interactive parameters.Consequently, several experimentally empirical estimates for their values have been adopted when using this theory, for instance DeTeresa and Larsen [42] proposed simple formulae to calculate the two interactive strength parameters F 12 and F 23 for materials that do not fail under hydrostatic compression or equal transverse compression.Accordingly, the LSF using Tsai-Wu failure criterion is 3.5.Hashin failure theory (HA) As described in Section 2.2, several criteria have been developed to distinguish failure mode from various physical phenomena.Hashin criterion [29] is one of the commonly used, and it is composed of four separate modes including matrix failure in tension and matrix failure in compression, and fibre failure in tension and fibre failure in compression.Accordingly, four LSFs are derived from these failure modes.For tensile matrix mode, For tensile fibre mode, r 11 > 0 For compressive fibre mode, r 11 > 0 3.6.Christensen failure criterion (RC) Christensen [28] derived a macroscopic failure criterion for composite materials using micromechanics, and the failure modes in composites are decomposed into matrix dominant failure and fibre dominant failure.LSFs based on Christensen failure criterion thus have two cases.For matrix controlled failure For fibre controlled failure, the LSF is 3.7.Transverse shear strength for aligned fibre composites Many of the failure criteria discussed in the preceding section require the transverse shear strength, S 23 , to determine failure parameters.Compared with other strengths, S 23 is particularly difficult to characterize experimentally.Inspired by the relationship between shear strength S and tensile and compressive strengths T and C for isotropic material, Christensen [43] developed a similar relationship for transversely isotropic material by using micromechanics.It will be used in the following studies, and it is expressed as: where g has some specific non-dimensional value that is to be determined.
where H is given by where m m is the matrix material Poisson's ratio, and k ¼ T m =C m with T m and C m being the tensile strength and compressive strength of the isotropic matrix material, respectively.

Multi-scale finite element based reliability analysis method for composites
The aim of a reliability analysis is to compute the probability of an event when the LSF reaches g 6 0. For composite materials, failure criterion is integrated with stochastic structural responses obtained from probabilistic structural analysis to form such a limit state function, gðsðxÞÞ, where x are considered random variables, such as material properties.In principle, reliability problem can be stated in the following mathematical expression as the multiple integrals of the joint probability density function f ðxÞ over the failure domain where p f is failure probability of interest.An exact solution of Eq. ( 28) is usually unavailable because of the multi-dimensional convolutions and the nonlinear nature of the LSFs.Hence, numerical approaches, such as Monte Carlo simulation, importance sampling, first/second-order reliability method (FORM/SORM) and Response Surface Methods (RSMs) [44], are usually employed to solve the problem in an approximate manner.For reliability analysis of unidirectional composites, FORM can provide sufficiently accurate relia-bility estimate [20,45].In FORM, the LSF g ¼ 0 is approximated by the first-order Taylor series expansion, and the essence of the approximation is to find an optimal point, namely design point, on g ¼ 0 with highest probability density in the transformed space of uncorrelated standard normal random variables y ¼ yðxÞ.Therefore, the searching for the design point is equivalent to solving a constrained optimization problem where GðyÞ ¼ gðxðyÞÞ denotes the LSF in the standard normal space.
The distance from the origin in the standard normal space to the design point is the reliability index b, and b is related to the failure probability by p f % UðÀbÞ, where U is standard normal cumulative distribution function.Several gradient based algorithms have been proposed in the literature to solve Eq. ( 29) including the HLRL algorithm, improved HLRL algorithm, etc.In the present study, the improved HLRL algorithm [46] is adopted and implemented due to its efficiency and accuracy.According to the chain rule of differentiation, the gradient vector of G can be obtained from where @x @y is Jacobian of the probability transformation, @g @s is straightforwardly obtained since the limit state function g is function of components of s in a simple algebraic form.Hence, the only unknown is @s @x .The classic finite difference method, direct differential method, or perturbation method is usually adopted to calculate @s @x .They are combined with finite element analysis when involving complex structure systems.Due to their multi-scale architecture, variabilities in composites originate from uncertainties at the micro-and macro-scale parameters.Hence, it needs to link micro-scale with macro-scale [6,7].Homogenization methods were adopted to propagate uncertainties from micro to macro, and they were integrated with perturbation method and finite element implementation to get the required derivatives of structural responses with respect to considered variables.The procedure is schematically presented in the flowchart shown in Fig. 2.

Linking microscale with macroscale using homogenization methods
The properties of composites are clearly determined by their constituent materials, their structure and distributions.Several approximate methods ranging from the classic rule of mixtures to leading homogenization theories can be used to predict the effective elastic properties of composites.In this study, three homogenization methods including the rule of mixtures (RM), Mori-Tanaka method (MT) and computational homogenization method (CH) were selected and their performance on reliability analysis was investigated.We considered a composite consisting of fibre with a transversely isotropic material and matrix with an isotropic material.This type of composite is generally considered as a transversely isotropic material, and its elastic properties can be described by five parameters E 1 ; E 2 ; m 12 ; m 23 and G 12 when aligning fibre direction in 1-direction in a 1-2-3 coordination system as shown in Fig. 1.

Rule of mixtures
RM is the simplest method to estimate the effective elastic properties of unidirectional fibre reinforced composites.They can be derived from constituent material properties using the following equations: where E and G are Young's and shear moduli of a material, respectively, m is a Poisson's ratio, and V is volume fraction, with superscripts f and m referring to the fibre and the matrix.With these effective elastic properties, it is straightforward to calculate the constitutive matrix, C RM .

Mori-Tanaka method
An important limitation of RM is the fact that it does not consider the mechanical interaction among different phases.More sophisticated asymptotic homogenization methods have been developed.A popular model is the mean field method proposed by Mori and Tanaka [47], which is based on Eshelby's strainconcentration tensor around an ellipsoidal inclusion in the matrix.The estimate of the overall composite stiffness tensor with where concentration tensor, I is the identity tensor, and E is the Eshelby tensor, which depends only on the inclusion aspect ratio and the matrix elastic constants.S and C are the compliance and stiffness tensors, respectively.For a typical unidirectional fibre reinforced composites consisting of fibre with transversely isotropic material and matrix with isotropic material, the Mori-Tanaka's tensor, A MT , can be found as where the non-zero elements are given by Therefore, the five effective elastic moduli can be recovered from Eq. (32).

Computational homogenization method
Although analytical homogenization approaches, e.g.Mori-Tananka, significantly improve the accuracy in predicting the overall properties of the composite material, limitations on the geometry of the representative microstructure and its constitutive response, which is often assumed to be linearly elastic, hinder their applications on composites with complex architecture, for instance, 3D composites.Multi-scale computational homogenization methods, which enable to analyse general geometries and nonlinear materials, have been developed, and they are becoming standard approaches for composites [48].The essence of computational homogenization method is to discrete fine scale fields on the representative microstructure and to obtain the material response at the macro-scale by solving a boundary value problem defined on a socalled representative volume element (RVE) of the microstructure.Details of the principle of the computation homogenization scheme adopted in this work can be found in e.g.[49].
In creating a numerical approximation, a key feature in solving the RVE boundary value problem is the application of appropriate boundary conditions.Three types of boundary conditions are commonly considered in the literature, namely (1) linear displacement, (2) periodic boundary condition and (3) uniform traction boundary condition.Periodic boundary condition generally provides the most accurate estimation among them.Hence, it will be considered in this study.The generalized RVE boundary condition enforcement approach proposed by Kaczmarczyk et al. [50] and its extension on 3D finite element method [51] was adopted.In a matrix form, CH is to solve the following equation: where k is the Lagrangian vector, e is the prescribed or given strain, and P and D are constraint matrix and global coordinate matrix, which are defined by where H is a matrix associated with boundary condition with details in [50,51], N is the standard shape function matrix and X is a position matrix evaluated at the integration points on the RVE boundary @X.Solving Eq. ( 36), the nodal displacement vector u f g and Lagrangian vector k f g can be obtained.These are then used to numerically calculate the effective tangent matrix using the following formula (see [51] for details of the derivation): Therefore, the effective elastic properties can be obtained from Eq. (38).

Propagating uncertainties from microscale to macroscale
Considering randomness in the material properties of the constituents and defining b ¼ fb 1 ; b 2 ; Á Á Á ; b n g T as an n-dimensional random vector, that, in the present case, comprises Young's modulus, Poisson's ratio, and shear modulus.Using the perturbation technique [52], an arbitrary stochastic function, uðbÞ, can be approximated via a Taylor series expansion.Although the general form of the perturbation-based stochastic finite element method enables higher-order approximations to be assumed, as proposed in [53], the commonly used second-order approximation is adopted here.It can be written as: where b is the mean value of the random vector b; db i denotes the variation around mean value of the i-th random variable, and H b i b j ðuÞ h i denote the first-and second-order partial derivatives of uðbÞ with respect to b i , and is a scalar representing a given small perturbation.
For the stiffness matrices C RM and C MT obtained from RM and MT, respectively, their stochastic functions can be directly expanded to obtained the corresponding approximation functions following Eq.( 39).This is due to the fact that they are explicitly expressed as a function of underlying constituent material properties.Thus, they can be approximated as For the CH-based analysis, the constitutive matrix C CH is numerically calculated.It is implicitly expressed with respect to underlying parameters, it is thus impossible to be approximated directly using Eq. ( 39).Various perturbation technique based probabilistic homogenization schemes have been proposed in the literature, e.g.[54,51].The probabilistic homogenization method developed by the authors in [51] is adopted in the present study, and the method is briefly summarised here.For convenience, Eq. ( 36) is written in a compact form as Considering material properties as random variables, the stiffness matrix K h i and the nodal displacement vector f ûg are thus stochastic functions.By expanding them in the form of Eq. ( 39), substituting into Eq.( 42) and equating terms of equal orders of , we arrive at the following zeroth-, first-and second-order equations: The zeroth-order The first-order The second-order From Eq. ( 36), the block related to stiffness matrix, K l Â Ã , of the microstructure in the compact matrix, K h i , is function of material properties.It can be expressed as and its first-and second-order partial derivatives are where B is the strain-displacement matrix, C l is the constitutive matrix of constituent material, and D bp C l and H bpbq C l are the first-and second-order partial derivatives of the material constitutive matrix.Hence, the expression of K h i and its first-and secondorder partial derivatives can be written as: Computing Eqs. ( 43)-( 45) successively, the compact displacement vector û f g and its first-and second-order partial derivative û f g and H bp bq û È É can be derived.These are then used to calculated the effective elastic moduli and its first-and second-order derivatives using Eq. ( 38) to form its approximation: Given the approximation for CðbÞ in Eqs. ( 40), ( 41) and ( 49), their statistics in terms of mean value and variance can be easily calculated from the following equations for the second-order approximation.Mean value is and the variance is

Stochastic structural responses from stochastic finite element method
Once the effective constitutive matrix for the composite is obtained, the finite element implementation to conduct stochastic structural analysis at the macro level can be realised by using the standard perturbation-based stochastic finite element method.In the present study, only the first order expansion is needed as the first order derivatives of structural responses are required in Eq. ( 30).This can be written as: The zeroth order The first order Solving the Eqs.( 52) and ( 53) consecutively, the nodal displacements u and its first-order partial derivatives, D b i u Â Ã can be obtained.Then other structural response terms can be calculated straightforwardly.It should be noted that the constitutive matrix for macroscopic structural elements is obtained from the homogenization method in the present study.Hence, in Eqs. ( 52) and (53) the stiffness matrix of macro-scale structural elements and its first-order partial derivatives are given by where C k H is the constitutive matrix obtained from homogenization method for the k-th ply.

Procedure of multi-scale reliability analysis
In summary, the procedure of multi-scale finite element reliability analysis starts with the identification of primitive variables and their respective distribution type and corresponding parameters at different scales including constituent material properties, geometrical parameters at the microscale, and variables at lamina level such as ply stacking orientation and lamina thickness.Probabilistic homogenization methods that combines the micromechanics based homogenization with the stochastic finite element method is used to propagate the scatter in primitive variables at micro-scale to uncertainties in the ply level mechanical properties with Eqs. ( 40), ( 41) and (49).Next, another perturbation-based stochastic finite element method is introduced to determine the stochastic structural responses corresponding to the selected primitive variables.Finally, the gradient of structural responses with respect to both micro-scale and macro-scale variables can be obtained, and these are supplied to the Eq. ( 30) to calculate the gradient of the LSF.Hence, the PSMFE and the FORM produce a multi-scale finite element based reliability analysis method.In addition, the commonly used Monte Carlo Importance Sampling (MCIS) is combined with the described multi-scale finite element method to estimate the probability of failure, which is served as a benchmark to verify the FORM based reliability calculation.

Results and discussions
In this research, a single-layer laminate and an 8-layer quasiisotropic laminate studied in the first worldwide failure exercise (WWFE-I) [55] have been analysed.Both of them are made of AS4 carbon fibres and 3501-6 epoxy matrix with a microstructure shown in Fig. 3.The lay-up of the quasi-isotropic laminate is ð90 = AE 45 =0Þ s .It has a total thickness of 1.1 mm with all plies having an identical thickness.Uniaxial and biaxial loading conditions have been considered for the single layer and the quasiisotropic multi-layer laminates.Results of reliability index and the corresponding probability of failure have been calculated by the previously introduced multi-scale finite element based reliability methods, which have been implemented on the Mesh Oriented Finite Element Method (MoFEM) program [56].

Estimating statistics of effective engineering properties from homogenization methods
Before conducting the reliability analysis, a comparative study was performed to show differences among various homogeniza- tion methods on estimating the effective elastic properties.Results for the five independent material properties of the composite are listed in Table 2 and measured data indicated in [55] are also provided for reference.It can be seen that there are some differences between calculated results from homogenization methods and measured data.For instance, the measured longitudinal Young's modulus E 1 is 126 GPa, while the calculated value is around 136 GPa from all the three methods.Gotsis et al. [57] had similar observation when participating WWFE-I with their micromechanics based method.In order to make the calculated properties match the measured data, Gotsis et al. [57] modified the constituent material properties provided in Table 1 to the following values: Although using the modified constituent material properties could make the estimated effective elastic properties match the measured data, the original constituent material properties were used for the following numerical studies.Variations in constituent material properties result in variability in the effective elastic properties.When treating constituent material properties as random variables with statistics given in Table 1, uncertainty can be propagated from micro to macro using Eqs.( 40), ( 41) and (49).To investigate the difference among homogenization methods on propagating uncertainties, the statistics of the effective elastic properties derived from Eqs. ( 50) and ( 51) are compared.The results are given in Table 3. Except for statistics, i.e. mean value and standard deviation or coefficient of variation, the probabilistic homogenization method enables to investigate correlations between the derived material properties at the macroscale.Tables 4-6 give the derived correlation coefficients for the five engineering properties.It is of interest to observe the correlations that exist between these effective engineering parameters.In general, the strength of relationship between variables is measured by the coefficient of correlation.Value in the range from À1.0 to À0.5 or from 0.5 to 1.0 indicates strong correlation, value from À0.5 to À0.3 or from 0.3 to 0.5 for moderate correlation, value from À0.3 to À0.1 or from 0.1 to 0.3 for weak correlation, and value from À0.1 to 0.1 for very weak or none correlation.From results in Table 4 derived from CH, it can be identified that there are strong correlation between G 12 and E 2 or m 12 and E 1 as their coefficients of correlation are larger than 0.5.In contrast, when using experiment test to measure these mechanical properties, they are usually assumed to be independent or it is very difficult to measure these correlations.Hence, it may lead to inaccurate results when assuming independence among mesoscale parameters.This will be further discussed in reliability analysis.

Comparison of multi-scale reliability analyses
As found in the previous section, there are differences among the three homogenization methods in terms of statistics shown in Table 3.A quasi-isotropic ð90 = AE 45 =0 Þ s AS4/3502 laminated plate studied in the WWFE-I was considered in the following study to further investigate the influence of these differences on the reliability of composite.Two loading conditions with one for uniaxial tension and another for uniaxial compression were tested.Reliability indices from these three homogenization methods are quite similar when the plate is under uniaxial tension condition.However, the calculated reliability indices differ from each other when the plate is under uniaxial compression.Results for the À45 ply from CH, RM and MT are shown in Fig. 4 when a series of compression is applied.The limit state function based on Tsai-Wu failure criterion as given in Eq. ( 18) was used.It can be seen that the reliability indices from RM are larger than those from the rest two methods, and CH provides the lowest estimates among these three approaches.Comparing with CH and MT, the results obtained from RM are significantly larger which may indicate more aggressive prediction as larger reliability index refers to the smaller probability of failure.Calculated reliability indices are also summarised in Table 7.The estimated effective elastic properties in Table 3 imply that those provided by MT are closer to CH than corresponding values obtained from RM. Furthermore, sensitivity factors, which are direction cosine to measure the relative importance of uncertainties in each variable, shown in Fig. 5 indicate that random variables have different contributions in various homogenization methods.This comparison may demonstrate that the homogenization methods with consideration of the interaction between various phases of materials provide more rational estimation in reliability analysis.

Comparison of multi-scale and single scale reliability analyses
A further investigation was conducted here to compare multiscale and single scale finite element reliability analysis.For the single scale finite element reliability analysis, only one instance of the finite element analysis, i.e. structural level, is required.In such case, the ply is treated as homogeneous material.Hence, all the random variables are at ply level.In practice, ply material properties are assumed to be independent when obtained from experimental tests, while the results derived from the probabilistic homgoenization methods show that they are correlated.Hence, the first analysis is to investigate the influence of correlation on reliability analysis.In single scale reliability analysis, two groups of calculations were carried out.Ply material properties and their statistics obtained from homogenization methods are used.For the first group, the ply material properties were treated as independent random variables, while they were considered as correlated variables with correlation matrix given in Tables 4-6.The reliability results are shown in Figs.6-8 for properties obtained from CH, RM and MT, respectively.These results indicate that the correlation has certain influences on the reliability calculation.Ignoring the correlation or assuming independences among variables may give a conservative estimate.Furthermore, another comparison was conducted that is between multi-scale analysis and single scale analysis.In single scale analysis, the derived material properties and their statistics from RVE analysis are used.Results in Fig. 9 show that there is a good agreement between these two methods.It indicates that the probabilistic homogenization method can accurately propagate uncertainties from micro to macro.In general, these two sets of comparative studies show that multi-scale reliability analysis can provide an efficient estimation of reliability for composites by including both microscale and macroscale uncertainties.More importantly, ignoring dependence between material properties at macroscale or ply level may result in inaccurate results for the reliability of composite.

Comparison of failure criteria in multi-scale reliability analysis
In this section, we focus on the issue about the selection of failure criterion for composites in establishing limit state function.The performance of six commonly used failure criteria was compared in terms of reliability by using the introduced multi-scale reliability analysis method.In the preceding section, CH based multi-scale reliability analysis method show better performance than the rest two methods.Hence, it was used in the following studies.

Single-layer laminate
A single layer composite plate was studied first.Constituent material properties and lamina strengths are considered as random variables, and their statistical information including the type of distribution, mean value and coefficient of variation is listed in Table 1.In total, twelve random variables with seven representing micro-scale parameters and five for mesoscale parameters were considered.Four analysis cases were chosen: two uniaxial stress cases of (1) tension in the longitudinal direction and (2) compression in the longitudinal direction, and two biaxial stress cases of (3)   For the two uniaxial stress cases, the plate is subjected to a tension of 0:7X T and a compression of À0:7X C , respectively.The estimates of reliability from different failure criteria for these two cases show that the six failure criteria almost provide the same values of the reliability index or failure probability.Thus, these results are not given to keep the article concise.Table 8 shows the results for the lamina under biaxial loading conditions.One is subjected to combined longitudinal and transverse tension with magnitudes of 0:55X T and 0:5Y T , and the other is under combined in-plane shear and longitudinal tension with magnitudes of s xy ¼ 0:25S 12 and r x ¼ 0:6X T .For the biaxial tension case, the three failure mode based criteria provide the most conservative reliability estimates CH MT RM among the six, and the Hoffman criterion provides the smallest estimate of reliability index.The calculations based on the maximum stress, Hashin and Christensen failure theories imply that the failure mechanism is tensile fibre failure.Hence, the probabilities of failure were calculated by using Eqs.( 5), ( 21) and (24), which are essentially based on the maximum stress failure.The sensitivity information shown in Fig. 10 provided by FORM indicates that all the six reliability indices are most sensitive to uncertainty at the longitudinal ply tensile strength.However, the reliability indices provided by the three failure theories without failure mode are negatively correlated to uncertainty in the transverse ply tensile strength.This leads to smaller reliability indices provided by the Tsai-Wu, Tsai-Hill and Hoffman failure theories.
For the tension and shear case, the reliability estimates from the maximum stress and the Christensen criteria are larger than those from the rest four criteria.The failure mode determined by the maximum stress and the Christensen criteria is again fibre tensile failure, while the Hashin theory indicated that the laminate failure is caused by matrix failure.Under the considered loading condition, there has shear stress in the structure.From the LSFs in Eqs. ( 5) and (24), only longitudinal stress and ply strengths are considered in the fibre tensile failure mode.The uncertainties in ply shear strength and the influences of shear modules on resulting structural responses are ignored.It thus leads to differences in the reli-ability estimates obtained from fibre tensile mode based analyses with those from matrix failure mode.It is further confirmed by the sensitivity information shown in Fig. 11, which indicates the uncertainty in shear strength has a significant influence on the reliability of the structure under the considered loading condition.
In addition, a stochastic failure envelope for r 1 versus r 2 based on the Tsai-Wu failure theory is produced.Since there is no analytical solution for Eq. ( 18), it is very time consuming to numerically evaluate it.Hence, only the envelope for reliability index, b ¼ 3:8, which is specified by the Eurocode 1 for ultimate limit states ''appropriate for most cases", is given for illustration purpose as shown in Fig. 12.A set of combinations of r 1 and r 2 were searched to get the reliability index of b ¼ 3:8, and these found r 1 and r 2 (see red circle points in the figure) were then used to establish an optimal envelope (see blue curve in the figure) by curve fitting.In accordance with the deterministic failure envelope, the area inside the envelope is the safe zone that has a smaller probability of failure or higher reliability index of b > 3:8.

Multilayer quasi-isotropic laminate
For the multilayer laminate, two uniaxial stress cases were investigated first.For uniaxial tension case, the results were obtained for the laminate subjected to a tension of 160 MPa, while a compression stress of 350 MPa was applied to the laminate for the uniaxial compression case.In order to simplify the understanding, only results for the first failure ply are reported in Table 9.When the laminate is under uniaxial tensile loading, the 90 ply has the largest probability of initiating damage under the considered loading condition, and the À45 and þ45 plies fail after the 90 ply but before 0 ply.The failure mode based criteria including maximum stress, Hashin and Christensen indicate that the damage initiates from the matrix failure.Comparing the estimates listed in Table 9 from different failure criteria, the Tsai-Wu criterion gives the smallest reliability index or the largest failure probability for the 90 ply, while the three failure mode based criteria provide the largest reliability index.It can be seen from the sensitivity information shown in Fig. 13 for the 90 ply that uncertainties in both micro-scale elastic properties and ply strengths affect the reliability of the structure under uniaxial tension.In addition, the slight differences in the reliabilities between the failure mode based group and the group without failure mode may result from the uncertainties in the longitudinal tensile and compression strengths.
For the uniaxial compression case, the À45 and þ45 plies fail prior to the 0 ply, and the 90 ply fails at last.The failure mode dependent criteria including Hashin and Christensen imply that the failure initiates from the matrix damage for the À45 and þ45 as the criteria of matrix controlled failure gives a larger probability of failure.Furthermore, the results from maximum stress failure theory further indicate that the failure mode is due to shear failure, which agrees with finding from the experiment reported in [58].Comparing the estimates from different failure criteria, the Hashin criterion gives the smallest reliability index or the largest  failure probability for the À45 and þ45 plies, while the Christensen criterion provides the largest reliability index.By observing the sensitivity information shown in Fig. 14, the positive contribution from the uncertainties in the longitudinal ply strength results in slight larger reliability estimates for the Tsai-Wu, Hoffman and Christensen theories based LSFs comparing with those from the rest three failure theories.For the 0 ply, the smallest reliability estimate is given by the Tsai-Wu failure criterion, but all the six failure criteria give similar reliability estimate.Importantly, compared with the judgement from deterministic failure analysis, which concludes that the damage initiated simultaneously in AE45 plies and the 0 ply [43], the reliability analysis finds that the damage initiates at the AE45 plies prior to the 0 ply.To create a biaxial tension loading condition for the quasiisotropic laminate, a tensile uniform pressure loading was applied on all side edges perpendicular to the x-and y-axes.Here, equal loads were applied along the x-and y-directions with r x ¼ r y ¼ 150 MPa.Results for reliability calculations by FORM and MCIS are listed in Table 9.As expected, all layers have quite similar reliability estimates or failure probabilities when using the same failure criterion.This is due to the fact that the laminate is quasi-isotropic.Comparing the estimates from different failure criteria, they generally have some extent scatter.The Tsai-Wu criterion implies that the highest failure load will be predicted as it provides the largest reliability or the smallest probability of failure.The other five theories provide very similar results with the lowest from Tsai-Hill, which implies the most conservative estimation.Furthermore, the three failure mode based criteria show that the damage initiated at the matrix.It can be seen from Fig. 15 that all the six failure theories have similar extents of sensitivity to uncertainties in matrix material properties E m and m m , fibre material properties m f 23 ; m f 12 ; G f 12 ; E f 1 and E f 2 , and ply strength Y T .The reliability estimates obtained from the Tsai-Wu and Hoffman theories based LSFs slightly differ from those from the rest four failure theories based LSFs due to the fact that they are affected by the uncertainties in X T and X C .

Conclusions
In this paper, the two critical issues, propagation of multi-scale uncertainty and selection of failure criterion, for the reliability analysis of composites were investigated.For the multi-scale uncertainty issue, three types of homogenization methods -rule of mixtures, Mori-Tanaka, and computational homogenization, were adopted to link micro and macro scales.They were further combined with the perturbation technique to propagate uncertainty at multiple scales.This enables statistical analysis of material properties at the macroscale when uncertainties are present in the constituent material properties at the microscale.By integrating these homogenization methods with the stochastic finite element method and structural reliability methods, such as the first order reliability method and Monte Carlo importance sampling used in the present study, the structural reliability of composites can be investigated with a selected probability failure function constructed on failure criterion for composite.This multi-scale reliability analysis procedure was applied to analyse laminated fibre reinforced composites made of AS4/3501 carbon/ epoxy.Six popularly used failure criteria in reliability analysis of laminated fibre reinforced polymer composites were compared and their performance were evaluated by using multi-scale reliability analysis.Some conclusions can be drawn from studies conducted in this paper:  By combining homogenization methods with perturbation technique to form probabilistic homogenization methods, it is possible to address the multi-scale uncertainty issue.These probabilistic homogenization methods enable to determine statistics of the effective elastic properties of composite materials with uncertainties in constituent material properties.Various probabilistic homogenization methods produce different estimates of statistics for material properties at the macroscale.These differences may result in different prediction in reliability analysis for composites in certain loading conditions.Among the examined three homogenization methods, the two advanced homogenization methods, i.e.Mori-Tanaka and computational homogenization, provide more rational results in reliability analysis.Probabilistic homogenization methods reveal that there are correlations between material properties at the macroscale, and ignorance of these correlations may lead to in accurate estimate of reliability index for the reliability analysis of composites.The comparison of failure criteria in reliability analysis demonstrates that the difference between the failure mode based criteria and the failure criteria without failure modes becomes more evident when the structures are in complex stress status, such as under biaxial or off-axis loadings.Differences among the examined failure criteria were investigated using multiscale analysis.In general, failure theories without failure modes are preferred in reliability analysis for composites due to their efficiency.

Fig. 1 .
Fig. 1.Coordinate system and profile of a unidirectional fibre-reinforced lamina.

Fig. 13 .Fig. 14 .
Fig. 13.Sensitivity factors of reliability index for the 90 ply of the quasi-isotropic laminate under uniaxial tension.

Fig. 15 .
Fig.15.Sensitivity factors of reliability index for the 90 ply of the quasi-isotropic laminate under biaxial tension.
ÀC 11 6 r 11 6 T 11 ; ÀC 22 6 r 22 6 T 22 ; ÀC 33 6 r 33 6 T 33 ; jr 12 j 6 S 12 ð1Þ is violated.Note that all seven strength parameters are treated as positive numbers, and the normal stresses are positive if tensile and negative if compressive.LSFs derived from these failure conditions are written as: g MS;1 ðxÞ ¼ T 11 ðxÞ À r 11 ðxÞ; or g MS;1 ðxÞ ¼ C 11 ðxÞ þ r 11 ðxÞ ð 2Þ g MS;2 ðxÞ ¼ T 22 ðxÞ À r 22 ðxÞ; or g MS;2 ðxÞ ¼ C 22 ðxÞ þ r 22 ðxÞ ð 3Þ g MS;3 ðxÞ ¼ T 33 ðxÞ À r 33 ðxÞ; or g MS;3 ðxÞ ¼ C 33 ðxÞ þ r 33 ðxÞ 34 GPa, E f 2 ¼ 14:96 GPa, G f 12 ¼ 14:96 Gpa, E m ¼ 5:8 GPa and m m ¼ 0:38.The calculated properties based on these modified materials properties are also given in Table 2.It can be observed that results from MT and CH are close to most of the measured data, while RM gives the least accurate estimates, especially for through-thickness components, i.e.E 2 and m 23 .In addition, the calculated properties from [57] by micromechanics are also given in the table to compare with those from the three homogenizations.Again, results from CH and MT are quite close to those from [57].These comparisons demonstrate that the two advanced homogenization methods, i.e.MT and CH can provide more accurate estimates of the effective elastic properties.

Table 1
Statistical properties of random variables.

Table 2
Comparison of effective elastic properties estimated from various methods (moduli in GPa).

Table 4
Correlations of effective elastic properties derived from computational homogenization method.

Table 5
Correlations of effective elastic properties derived from rule of mixture method.

Table 6
Correlations of effective elastic properties derived from Mori-Tanaaka method.

Table 3
Comparison of the statistics of the effective elastic properties estimated from various methods (moduli in GPa).

Table 7
Reliability index from various situations.
Comparison of reliability index calculated using RVE and Macro.
Fig. 11.Sensitivity factors of reliability index for the single-layer laminate under tension and shear.

Table 8
Reliability indices and probabilities of failure obtained from various failure criteria for the single layer laminate under biaxial loading.

Table 9
Reliability indices and probabilities of failure obtained from various failure criteria for the quasi-isotropic laminate.