Stochastic multi-scale finite element based reliability analysis for laminated composite structures

This paper proposes a novel multi-scale approach for the reliability analysis of composite structures that accounts for both microscopic and macroscopic uncertainties, such as constituent material properties and ply angle. The stochastic structural responses, which es- tablish the relationship between structural responses and random variables, are achieved using a stochastic multi-scale ﬁnite element method, which integrates computational ho- mogenisation with the stochastic ﬁnite element method. This is further combined with the ﬁrst- and second-order reliability methods to create a unique reliability analysis frame- work. To assess this approach, the deterministic computational homogenisation method is combined with the Monte Carlo method as an alternative reliability method. Numerical examples are used to demonstrate the capability of the proposed method in measuring the safety of composite structures. The paper shows that it provides estimates very close to those from Monte Carlo method, but is signiﬁcantly more eﬃcient in terms of computational time. It is advocated that this new method can be a fundamental element in the development of stochastic multi-scale design methods for composite structures. proposed methods inherit the shortcoming of the gradient-based methods, where the accuracy of the results reduces with increases in nonlinearity of the limit state function and coeﬃcients of variations of the random variables. Two multilayer laminates studied in the well-known World Wide Failure Exercise have been investigated using the PSMFEM-FORM method - a quasi-isotropic laminate made of AS4/3501-6 carbon ﬁbre reinforced polymer composite and an E-glass/LY 556 epoxy laminate. Numerical examples illustrate the performance and capability of the proposed method. The proposed method will serve as a fundamental component in the development of stochastic multi-scale design method for composite structures.


Introduction
Typical composite components are laminates comprising layers of fibre reinforced composite laminae, each of which are made of fibres embedded in matrix. The assembly of the fibres and matrix materials to create a lamina, as well as the lay up and curing of laminae, is a complicated process and may involve a lot of uncertainty. Therefore, the material properties of a composite laminate are random in nature. Sources of significant uncertainty include: variations in volume fractions of fibre and matrix, voids in the matrix and between fibres and matrix, imperfect bonding between constituents, cracks, fibre damage, random and/or contiguously packed fibres, misaligned fibres, temperature effects, non-uniform curing of the matrix material, residual stresses, etc. Uncertainties in these factors propagate to a larger scale and are reflected in variability of the stiffness and strength that characterise the overall structural behaviour [1][2][3][4][5][6][7] . Consequently, high safety factors of the order of 8-10 [8] are introduced in current deterministic based structural design, thereby not taking full advantage of composite materials. These issues may be addressed in a probability based design context [9][10][11] equivalent to Eurocode, for example.
The reliability index is one of the widely accepted indicators to measure safety of engineering structures subject to uncertainties [12] . The reliability analysis of composites involves several key issues [13][14][15] . First, many parameters may need to be considered as random variables, but they can be associated with different length scales: micro-scale (constituent level, fibre/matrix), meso-scale (ply level), or macro-scale (component level) [16] . Hence, the choice of uncertainties will directly determine which mechanical model should be used in the structural analysis stage. For instance, a micro-mechanical model is needed when considering micro-scale parameters as random variables. Second, composite materials display a wide variety of failure mechanisms due to their complex structure and manufacturing process, and a range of failure criteria have been proposed, such as maximum stress/strain criteria and polynomial criteria (for instance, Tsai-Wu criteria [17] ). Despite extensive research, such as the well known three stages of World Wide Failure Exercises [18][19][20] , a complete and validated methodology for predicting the behaviour of composite structures including the effects of damage has not yet been fully achieved. Qualitative evaluations [13,21,22] , quantitative experimental comparisons [23] , and numerical comparisons [24] have been made for deterministic failure criteria and quantitative comparison of their performance considering uncertainty was reported in [25] , with broad differences. Hence, the choice of failure criterion to establish the limit state function is critical to conduct reliability analysis. Third, for some relatively simple composite structures, analytical formulations have been developed [26][27][28] , while finite element reliability analysis methods are necessary to handle more complex structures [29] . Finally, the essence of reliability analysis is to calculate the failure probability of structural components or systems, which is expressed by the convolution integral. It is impossible to calculate the integration directly due to its multi-dimensionality, and therefore numerical methods have to be used. Monte Carlo simulation (MCS) is a straightforward option [30] . The accuracy of MCS directly depends on the number of simulations, and it is recommended to conduct at least N simulations [31] , where N = − ln (1 − C) / p f with p f is the expected failure probability and C is a given confidence level (normally 95%), to obtain a sufficiently accurate estimate of p f . However, the failure probability is generally a small value of the order of 10 −5 , e.g. Eurocode requires β ≥ 3.8 or p f ≤ 7 × 10 −5 , and implementing thousands of simulations are thus extremely time consuming, in particular when involving finite element calculations for complex structural systems. Approximation methods such as First-Order Reliability Method (FORM) and Second-Order Reliability Method (SORM) have become popular alternatives due to their efficiency. Recently, new reliability methods have been adopted or proposed to conduct reliability analyses for composite structures, such as Artificial Neuronal Networks [32,33] . Nakayasu et al., [25,34] compared different structural reliability analysis methods for composites and stated a preference of the FORM method.
It has long been recognized that laminate stress analysis and lamina failure criteria are two critical elements in the analysis of composite laminates. In the past years, various homogenization methods have been developed in order to determine the macroscopic material properties of a heterogeneous material from its constituents. However, most of the previous work on reliability analysis for composites had focused only on meso-scale parameters such as ply material properties. It has been widely accepted that the mechanical behaviour of composites are strongly affected by their microscopic variations in material properties [16,35] and homogenization methods have proven to be capable of predicting accurately the effective material properties. Accordingly, the combination of probabilistic modelling and micromechanics seems to be an appropriate approach to achieve the consistent characterisation of composites behaviour [36,37] . The accuracy of structural reliability estimation may be improved by using multi-scale methods [38,39] . Although stochastic multi-scale finite element methods have been developed for many years by various researchers, e.g. [40][41][42] , except in a few cases [43][44][45][46][47][48] , multi-scale modelling of such materials has been limited to purely deterministic analyses.
The objective of the present study is to propose and evaluate a multi-scale finite element based reliability analysis approach in order to address some of the above mentioned challenges. This approach combines a state-of-the-art computational multi-scale homogenization method with composite mechanics and structural reliability analysis. It enables uncertainties in both microscopic and macroscopic parameters to be considered. Stochastic structural responses are obtained from a stochastic multi-scale finite element method, which establishes the relationship between structural responses and microscopic random variables. The commonly used FORM and SORM are coupled with the multi-scale finite element method to conduct reliability analyses. Numerical studies are performed to illustrate the procedure for the reliability analysis of composite structures and to demonstrate the efficiency and accuracy of the proposed approach. The proposed method will serve as a fundamental component in the development of stochastic multi-scale design method for composite structures. Innovative construction and building technologies are usually driven by the developments of new constructional materials and/or structural forms [49,50] .

Stochastic homogenization method for composite materials
First we will summarize the basic assumptions and the final formulae of the probabilistic homogenization method for the estimation of effective elastic moduli developed by the authors in [42] . The class of homogenization-based multi-scale constitutive models employed in the present study is characterised by the assumption that the strain ( ε ) and stress tensor ( σ ) at a point of the so-called macro-continuum are the volume average of their respective microscopic counterpart fields over a pre-specified Representative Volume Element (RVE).
where ε μ , σ μ and u μ denote the strain, stress and displacement fields of the RVE and V μ is the volume of the RVE.
The RVE displacement field can be decomposed as: where the first term is associated with a homogeneous strain field, corresponding to the macroscopic strain u * = ε y ( y is the local RVE position vector), and the second term is the displacement fluctuation field. Similarly, we can decompose the RVE strain field into a homogeneous part and a fluctuation part as: The volume average of this strain field must be equivalent to the macroscopic strain, thus To make the upscaling transition, we assume the well known Hill-Mandel principle holds, requiring that Thus, the RVE equilibrium problem for a linear-elastic material consists of finding, for a given macroscopic deformation ε , a kinematically admissible displacement fluctuation field η ∈ V, such that where C μ is material constitutive tensor. Given uncertainties in the material properties, defined by vector b , the constitutive matrix C μ and displacement fluctuation ˜ u μ are stochastic functions of b . According to the perturbation method [51] , an arbitrary stochastic function, φ( 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 [52,53] , the commonly used second-order approximation is sufficient for the present study. It is therefore adopted here and can be written as: where b is the mean value of the random vector b ; D b i (·) and H b i b j (·) denote the first-and second-order partial derivatives of ( ·), and is a scalar representing a given small perturbation. By expanding the stochastic functions C μ , ˜ u μ and ∇ s ˜ u μ in Eq. (6) to the forms similar to Eq. (7) and substituting into Eq. (6) , and equating terms of equal orders of , we arrive at the following zeroth-, first-and second-order virtual work principles: • The second-order where N denotes shape function, B is the strain-displacement matrix, ˜ a μ is nodal displacement fluctuation vector, δa is virtual nodal displacement fluctuation vector, a * denotes the given nodal displacement vector. The finite element approximation to the zeroth-, first-and second-order variational principles in Eqs. (8) - (10) , respectively are obtained as: • The zeroth-order • The first-order with the RVE finite element stiffness matrix and its first-and second-order partial derivatives defined as: The solution for the microstructure displacement fluctuation field ˜ a μ and its derivatives D b p ˜ a μ and H b p b q ˜ a μ can be found by applying boundary conditions consistent with the applied macro strain field. The three classic boundary conditions are linear displacement, periodic displacement and anti-periodic traction, and constant tractions, the details of which are given in [42] . With these at hand, the effective constitutive relation between the applied macrostrain ε and the homogenized macrostress σ can be calculated given, Similarly, the first-order partial derivative, D b p C , and the second-order partial derivatives, H b p b q C can be calculated.

Macroscale finite element formulation for calculating structural response
A stochastic structural analysis at the macroscale can be realised by using the standard perturbation-based stochastic finite element method, requiring the successive solution to the following set of equations for the nodal displacements u and then its first-and second-order partial derivatives, The effective constitutive matrix for the macroscopic analysis is obtained from the homogenization method in previous section. Hence, the stiffness matrix of macroscale structural elements and its first-and second-order partial derivatives in Eqs. (16) -(18) are given by In the above, C k denotes the transformed constitutive matrix of the ply, orientated by the angle θ : where R is the strain coordinate transformation matrix. It should be noted that the first-and second-order derivatives of C k with respect to the ply orientation angle are obtained as follows (21) and After solving Eqs. (16) -(18) for the nodal displacement field u and its first-and second-order derivatives, D b i u b and , structural responses and their derivatives can be straightforwardly deduced. Here, we give details to compute stress and its first-and second-order derivatives that will be used in the reliability analysis in next section when using stress based failure criterion such as Tsai-Wu criterion in Eq. (27) . The stress in k th layer and its derivatives can be calculated through:

Fundamental of structural reliability analysis
In deterministic analysis of composites, failure of a structural component is typically determined by a particular failure criterion, e.g. maximum stress failure theory. However, in a reliability analysis context, it is the probability of failure, p f , that is considered, where applied loads and structural properties are considered random variables. These variables are collected in a n -dimensional vector x with assigned distributions, F i (x i ) , i = 1 , 2 , . . . , n, reflecting their uncertainties. Structural failure events are specified in terms of limit-state functions (LSF) g ( s ( x )), where s denotes the response quantity, such as displacements, strains, stresses, forces and cumulative response measures. Thus, the structural component or system is considered to have failed when the LSF g ( x ) ≤ 0, e.g., a response quantity exceeding a safe threshold.
A mathematical statement for the reliability problem is that p f is given by a multiple integral over the failure domain where f ( x ) is the joint probability density function of the random vector x . An exact or closed form solution of Eq. (26) is not usually available when involving many random variables and a nonlinear LSF. Hence, various numerical methods are usually employed to solve the problem in an approximate manner such as MCS, FORM, SORM, response surface method, etc.
In the present study, MCS and FORM are adopted to combine with the multi-scale finite element method to derive reliability formulations for composite structures.

Limit state function for composites
The first step in reliability analysis is to establish the LSF. For laminated fibre reinforced composites, this is based on lamina failure theory. There are several failure theories available, such as maximum stress, Tsai-Hill and Tsai-Wu. Among these methods, the Tsai-Wu failure criterion [17] is the most frequently used because it takes into account the interactions between different stress components and provides a relatively accurate prediction of failure compared with other theories [22] . Here, we adopt the Tsai-Wu criterion to illustrate the creation of a LSF for ply failure. The procedure is applicable to other failure criteria, for instance, six widely used failure criteria were compared in [39] on the performance of evaluating reliability of laminated composites by using the multi-scale finite element reliability approaches developed in the present paper. The LSF based on the Tsai-Wu criterion is expressed as: where σ denotes stress in the ply in various directions with 1 representing the direction along fibre and 2 and 3 representing transverse and through thickness directions, respectively, T 11 and C 11 are tensile and compressive strengths of ply in fibre direction, respectively, T 22 and C 22 are transverse tensile and compressive strengths, respectively, S 12 is the in-plane shear strength, S 23 is the transversal shear strength, and F 12 is the interaction parameter. A biaxial tensile test is required to determine F 12 and several empirical formulae have been suggested in the literature. Here we use the Mises-Hencky definition: In addition, the transverse shear strength, S 23 , is particularly difficult to characterise experimentally. [54] attempted to establish a relationship by using micromechanics, whereby the asymptotic relationship is derived as (28) where η has some specific non-dimensional value that is to be determined from where ψ is given by with ν is the matrix material Poisson's ratio, and λ = T /C with T and C the tensile and compressive strength of the isotropic matrix material, respectively.

Multi-scale finite element based Monte Carlo simulation method
A direct way to compute the probability of failure defined in Eq. (26) is by Monte Carlo simulation. Here, we combine the previously described multi-scale finite element method (MFEM) with MCS to form a simulation-based reliability approach (MFEM-MCS). In this approach, primitive random variable distributions are used together with the MFEM to provide the values required for the evaluation of the LSF, Eq. (27) . In the present study, this method is implemented in MoFEM (meshoriented finite element method) program [55] and the primary random variables are sampled according to their probabilistic features by a random number generator. The MFEM is then used to calculate the stresses which, together with the strength variables, are used to evaluate whether the failure condition ( Eq. (27) ) is violated or not. This procedure is repetitively executed for a desired number of trials, N . The number of trials for which g ( x ) ≤ 0 is counted and denoted as, n f . The probability of failure can thus be approximately calculated as

Multi-scale finite element based first-order reliability method
MCS is a computationally demanding option for obtaining acceptable results especially for problems with many random variables, such as, in particular, found in the case of composites. Importance sampling may be used to reduce the number of required samples, but this requires careful selection of the distributions used to create realisations around the failure point [56] . A combination of efficient approximation methods with the multi-scale finite element method is thus required. The approximation methods depend on the calculation of failure probability for quantifying reliability. Normalization recognizes the direct link between the notion of reliability, β, and its measure by a probability, p f , and it is expressed as where is the standard normal cumulative distribution function. By transforming LSF, g ( x ), is transformed into standard normal space as G ( y ), the reliability is geometrically defined as the shortest distance from the origin in the standard normal space to the surface G ( y ) = 0 or highest probability density among all possible realisations in the failure domain G ( y ) ≤ 0. Hence, reliability calculation is actually to find an optimal point, namely the design point, y * , on the surface G ( y ) = 0 to satisfy the constraint equation Consequently, solving Eq. (33) is a key element in estimating the reliability index. Gradient based approximation methods, such as FORM and SORM, are widely used. FORM is considered in the present study. In FORM, the transformed LSF G ( y ) is replaced by its hyperplane tangential to the design point through the first-order Taylor series expansion, which is written as where ∇G ( y ) = [ ∂ G/∂ y 1 , . . . , ∂ G/∂ y n ] denotes the gradient vector. According to the chain rule of differentiation, the gradient vector of G can be obtained from where ∂g ∂ s is the partial derivative of the limit state function with respect to structural response quantity, ∂ s ∂ x is the partial derivative of structural response quantity with respect to the basic random variables, and ∂ x ∂ y is the Jacobian of the probability transformation. From the limit state function g ( x ) such as Eq. (27) , the components in ∂g ∂ s can be easily obtained through simple algebraic operations. For terms in the Jacobian ∂ x ∂ y , if the random variables are independent, the non-diagonal elements are zeros and the diagonal components can be calculated as where f i ( x i ) is the probability density function (PDF) of the random variable x i , and φ( y i ) is the univariate standard normal PDF. For dependent random variables, a decomposition step is needed to convert them into independent random variables prior to projecting them into standard normal space. Two types of decomposition method, namely Nataf transformation and Rosenbaltt transformation, are available depending on whether the marginal distribution are known or not, respectively. The details of these two transformation methods can be easily found in the textbook such as [56] .
It is clear that the only unknown part in Eq. (35) is the derivative of the response quantity with respect to the basic random variable, ∂ s ∂ x . In finite element reliability analysis, these gradients can be obtained using the classic finite difference method, the direct differential method and the perturbation method. For iterative calculations, it is essential such gradients are computed efficiently and accurately. The application of the finite difference method is impractical as it is quite time consuming. Although, the direct differential method is an ideal option due to its efficiency, it is not readily available for the multi-scale finite element case. In the present study, the perturbation method based formulations given in Section 2 is utilized. The gradients of structural response quantities with respect to the basic input random variables in Eq. (35) can be obtained from Eqs. (16) - (18) .
An approach, namely PSMFEM-FORM, that integrates the perturbation based stochastic multi-scale finite element method (PSMFEM) and FORM is proposed to perform reliability analyses for composite structures. The implementation is schematically illustrated in Fig. 1 , and the main steps are summarized as follows: (1) define a limit state function as Eq. (27) for the reliability module; (2) perform the probability transformation to obtain the realization of the original random variables x i from the trial point y i in the standard normal space; (3) pass the trial point x i to the multi-scale finite element module, and separate these realizations into those related to microscale and macroscale. They are then passed to the corresponding finite element sub-module; (4) send the response quantities, s ( x ) and the gradient of the response, ∂ s / ∂ x , which contribute to the evaluation of the limit state function g ( x ) and its gradient ∂ g / ∂ x , from the finite element module to the reliability module; (5) check for convergence of the search to the design point and determine the search direction d and search step size κ if convergence is not achieved; (6) compute the failure probability using Eq. (32) . It should be noted that several iterations are usually required to reach the convergence. The dashed arrows in Fig. 1 illustrate the interaction between the reliability module and the multi-scale finite element module that transfer realizations of the random variables and return the response quantities.

Multi-scale finite element based second-order reliability method
It is well known that FORM is an effective solution strategy if the distance to the limit state has only one minimum, and the function is nearly linear in the neighbourhood of the design point. However, if the failure surface is highly nonlinear, the failure probability estimated by FORM may give unreasonable and inaccurate results. To resolve this problem, higher order approximation has been developed [57][58][59][60] . Here, we integrate second-order reliability method with perturbation based stochastic finite element method to develop a multi-scale finite element based second order reliability method. We adopted the point-fitting SORM developed by Zhao and Ono [60,61] . A brief review of this method and how it links with PSMFEM is given in the following. Similar as FORM, SORM is based on Taylor series expansion to approximate the LSF but with second order approximation. Hence, the second-order Taylor expansion of a LSF in y -space at design y * can be expressed as: Dividing by | ∇G ( y * )| and considering G ( y * ) = 0 , we obtain (38) where α = −∇G ( y * ) / |∇G ( y * ) | is the directional vector at the design point in y -space, H is the scaled second-order derivatives of G ( y * ), known as the scaled Hessian matrix, H = H / |∇G ( y * ) | = ∇ 2 G ( y * ) / |∇G ( y * ) | , and β F = αy * is the first-order reliability index.
The principle of point-fitting methods is to define the second-order surface approximation in terms of a set of fitting points on the limit state surface in the vicinity of the design point. These points, 2 n + 1 in number, are selected along the coordinate axes in standard normal space. The limit state surface is adopted by a rotational parabolic surface of diameter 2 R , where R is the average principal curvature radius expressed as: where κ S is the sum of the principal curvatures of the limit state surfaces at the design point, which can be expressed as: H j j − α T H α (40) By doing so, the performance functions in y -space, can be expressed simply as: A critical part of this calculation is to evaluate the Hessian matrix, which requires second-order derivatives of LSF with respect to the assumed random variables. This requirement may be particularly onerous in the case of a limit-state function that are described by numerical algorithms, such as finite element calculation. The perturbation technique offers a convenient method to obtain first-and second-order derivatives derivatives, for example. As shown in the proceeding section, the second-order derivatives of structural responses in x -space are obtained from the second-order perturbation equation as given in Eq. (13) or Eq. (23) , and then the chain rule is used to transform it from x -space to standard norm space. By doing so, we can get the required Hessian matrix for Eq. (40) to calculate principle curvature κ for the limit state surface.
Afterwards, the second order reliability index can be computed through the empirical formulae given in Eq. (42) . As it links with PSMFEM, we therefore named it as PSMFEM-SORM. The procedure of the SORM approximation according to the above formulation is also schematically illustrated in Fig. 1 .

Numerical examples
In this section, numerical examples are presented to demonstrate the methods described in the preceding sections. A single layer laminate is investigated first to compare MFEM-MCS with PSMFEM-FORM and PSMFEM-SORM in structural reliability prediction. Two multilayer laminates, including one made of carbon fibre reinforced epoxy matrix and another made of glass fibre reinforced epoxy matrix, are further considered to show the capabilities of PSMFEM-FORM and SORM.

Example 1: single layer laminate
In this example, we study a single layer laminate that is made of carbon fibres embedded in an epoxy matrix with material properties indicated in [62] and listed in Table 1 . The dimensions of the structure are given in Fig. 2 (a). The structure is subject to a uniformly distributed load in the fibre direction (or x -direction) to simulate an uniaxial tension test. Considering the symmetry of structural dimensions, boundary conditions and loading, only 1/8 of the structure is needed, and the substructure is shown in Fig. 2 (b). To calculate the effective stiffness matrix, a cubic RVE shown in Fig. 2 (c) is constructed with hexagonal arrangements of unidirectional fibres. A local Cartesian coordinate system (1-2-3) is introduced at the microscale with axis 1 aligned parallel to the direction of the fibres.
Four cases are considered to demonstrate accuracy and efficiency of the proposed method. In each case, different parameters are considered as random variables. For example, in case 1 only the applied force is treated as a random variable, whereas all 13 parameters are random variables in case 4. The reliability calculation is implemented in MoFEM, where the modelled structure is discretized into four-node tetrahedral elements. PSMFEM-FORM is conducted first, and the estimated probability of failure is then used to determine the number of simulations for MFEM-MCS. For PSMFEM-FORM, as described in the flowchart Fig. 1 , once the initial realizations of considered random variables are passed to the multi-scale finite element module, the effective constitutive matrix of the lamina and its first-order derivatives are calculated using Eqs. (11) , (12) and (15) . With the obtained effective constitutive matrix, the stiffness matrix and its first-order partial derivatives of macroscopic element are readily calculated from Eq. (19) . Structural response quantities and their first-order derivatives are then computed using Eqs. (16) and (17) . These are returned to the reliability module to evaluate the LSF and its gradient in Eqs. (27) and (35) , respectively. Then the convergence criterion checks whether the trial point is acceptable. If convergence  is reached, the nominal probability of failure is determined. Otherwise, a new trial point, y i +1 = y i + κd , is generated. Afterwards, SORM is called to provide enhanced results with second order approximation of the LSFs. For MFEM-MCS, once the probability of failure, p f , from PSMFEM-FORM is obtained, the number of simulations can be determined and the procedure described in Section 3.3 is followed.
Results obtained from PSMFEM-FORM and PSMFEM-SORM and MFEM-MCS are compared in terms of computation time and reliability index, as listed in Table 2 . In general, these results show good agreement between the three approaches in all four cases. It has a trend that PSMFEM-SORM becomes close to MFEM-MCS with the increase of the number of random variables, which introduces the nonlinearity of the LSF. As expected, PSMFEM-FORM and PSMFEM-SORM have a significant advantage in computational time. Moreover, as a by-product, PSMFEM-FORM provides sensitivity factors that measure the relative importance of each random variable. Sensitivity factors presented in Fig. 3 show that the scatter in the applied force has the most significant influence on the ply reliability.

Multilayer laminates
As demonstrated in the previous section, PSMFEM-FORM and PSMFEM-SORM can provide a sufficiently accurate estimate of reliability index or probability of failure whilst also providing significant reduction in computational effort. In this section, PSMFEM-FORM and PSMFEM-SORM were used to analyse two multilayer laminates, which were studied in the first World Wide Failure Exercise (WWFE-I) [18] .

Example 2: carbon fibre reinforced composite laminate
A quasi-isotropic (90 °/ ± 45 °/0 °) s AS4/3502 laminated plate used in the WWFE-I is considered in this example as it is an important class of composites, familiar to the aerospace industry [63] . Failure of the quasi-isotropic layup is probably the most important single case that can be studied, as it represents a limiting case of all possible layups. The other limiting case is the unidirectional form itself, and all other cases lie between these two extremes. Thus, the quasi-isotropic case is one of the most severe tests of using a lamina-level failure criterion to predict damage and failure of a laminate [54] . The geometric features of this laminate are shown in Fig. 4 a, and the total thickness of the laminate is 1.1 mm with all layers having identical thickness [63] . Constituent material properties, ply strengths and ply orientation angles were considered as random variables, and their statistical information including distribution type and parameters are given in Table 3 . Material properties were assumed to be random variables with normal distributions. Ply strengths had lognormal distributions, and the four ply orientation angles were considered to be uniformly distributed random variables.
First, ply failure under uniaxial tensile stress in the x-direction (see Fig. 4 b) was studied. Results are shown in Fig. 5 . Both PSMFEM-FORM and PSMFEM-SORM were used to calculate reliability of the laminated under various amplitudes of stress, and results from PSMFEM-FORM are given in solid lines while circles for PSMFEM-SORM. It can be seen that Ply 1 (90 °p ly) with fibres perpendicular to the applied force has a much greater probability to fail compared with Ply 4 (0 °ply) with fibres parallel to the force. The plies, whose orientations are between these two cases, have a chance of failure in between.  In Fig. 5 , only the probabilities of failure for Ply 2 (45 °ply) are given as results for Ply 3 ( −45 • ply) are almost identical to those for Ply 2. Clearly, these results indicate that the probability of failure of the plate depends on the probability of failure of the Ply 1.
Next, a consideration of the laminate under biaxial tension (see Fig. 4 c) was carried out. A tensile uniform pressure loading is applied on all side edges perpendicular to the x -and y -axes. Estimates of reliability indices and probabilities of failure are provided in Table 4 . Again, both FORM and SORM were used. As expected, all plies have quite similar reliability estimates or failure probabilities. This is due to the symmetry of the problem. Interestingly, the probabilities of failure of the plies under biaxial tension are greater than corresponding results for uniaxial tension given in Fig. 5 . For instance, the probability of failure of the 90 °ply under σ x = σ y = 150 MPa is 8 . 86 e − 6 , while the probability is about 2 . 5 e − 6 when the plate is under tension in x-direction of σ x = 150 MPa. Although the increase of the probability of failure is not significant, it implies that this type of plate has higher resistance when subjected to a biaxial load state.

Example 3: glass fibre reinforced composite laminate
Glass fibre reinforced polymer composite is an important material class in engineering applications, such as civil engineering [64] , due to their relatively low cost. A laminate made of E-glass/LY 556 epoxy studied in WWFE-I [18] is considered in this section. The layer orientation of this laminate is illustrated in Fig. 6 a, and the total thickness of the laminate is 2 mm with the ± 30 °layers having 82.8% of the total thickness of the laminate and the 90 °plies having the rest 17.2% [63] . In this reliability analysis, constituent material properties, ply strengths and ply orientation angles were considered as random variables, and their statistical information including distribution type and parameters are given in Table 5 Fig. 6. GFRP laminate (a) geometry, (b) uniaxial tension, (c) biaxial tension.  variables with normal distributions. For macroscopic uncertainties, the lamina strengths ( X T , X C , Y T , Y C , S 12 ) of the material were simulated as independent random variables with lognormal distributions, and the three ply orientation angles ( θ 1 , θ 2 , θ 3 ) were considered as uniformly distributed random variables.
The laminate under uniaxial tensile stress in the x-direction (see Fig. 6 b) was investigated first. The probabilities of failure of each ply calculated by the PSMFEM-FORM and SORM methods are shown in Fig. 7 , where the results from PSMFEM-FORM are depicted with solid line while red circles for PSMFEM-SORD. Various magnitudes of applied stress were investigated. Since the 30 °ply, i.e. Ply 2, and −30 • ply, i.e. Ply 3, have similar failure behaviour, only the results for the Ply 2 are given in Fig. 7 . It can be seen that the ply, i.e. Ply 1, whose fibres perpendicular to the applied force is more prone to failure than the ply with fibres in direction 30 °or −30 • . These results imply that the probability of failure of the plate depends on the probability of failure of the Ply 1. In addition, results from FORM and SORM are very similar for most loading cases, although there are slight biases when probabilities of failure are close to one.
The laminate was further considered under combined loading conditions (see Fig. 6 c). An additional force σ y was applied in y -direction with σ y = σ x . Reliability indices and probabilities of failure for the plate under biaxial tension stress of 40 MPa, 50 MPa, 60 MPa and 70 MPa are calculated and listed in Table 6 with the use of both FORM and SORM. In contrast to the uniaxial tension stress case, the results indicate that ± 30 °plies have greater probabilities of failure compared with 90 °p ly. In other words, the introduction of ± 30 °plies provide protection to the 90 °ply as its probability of ply decreases. Moreover, the plies demonstrate different responses for uniaxial tension and biaxial tension. For 90 °ply, the probabilities

Conclusions
In this paper, a multi-scale finite element based reliability analysis method is proposed for composite structures. The proposed method enables both microscopic uncertainties, such as those in constituent material properties, and macroscopic uncertainties, such as ply orientation angles, to be taken into account. The new reliability analysis framework (PSMFEM-FORM and PSMFEM-SORM) couples a recently developed perturbation-based stochastic multi-scale finite element method with the first-and second-order reliability methods. Formulations of the stochastic homogenization method and its finite element implementation are described. Procedure to perform the proposed multi-scale finite element based reliability analysis methods is thoroughly introduced. In addition, Monte Carlo simulation has been coupled with the deterministic multi-scale finite element method to form the MFEM-MCS to evaluate the performance of PSMFEM-FORM and PSMFEM-SORM.
Numerical examples have been provided to demonstrate the capability of the proposed method in measuring the safety of composite structures induced by variations in microscopic and macroscopic parameters. A single layer laminate is designed to compare the three methods, namely MFEM-MCS, PSMFEM-SORM and PSMFEM-FORM, and metrics, such as reliability index, probability of failure, and computational time, are used to evaluate their relative performance. It is found that the three methods provide very close estimates of reliability indices and probabilities of failure for different cases, but the PSMFEM-FORM and PSMFEM-SORM are much more efficient in terms of computational cost, especially for cases with small probabilities of failure. It should also be kept in mind that the proposed methods inherit the shortcoming of the gradientbased methods, where the accuracy of the results reduces with increases in nonlinearity of the limit state function and coefficients of variations of the random variables. Two multilayer laminates studied in the well-known World Wide Failure Exercise have been investigated using the PSMFEM-FORM method -a quasi-isotropic laminate made of AS4/3501-6 carbon fibre reinforced polymer composite and an E-glass/LY 556 epoxy laminate. Numerical examples illustrate the performance and capability of the proposed method. The proposed method will serve as a fundamental component in the development of stochastic multi-scale design method for composite structures.