The natural origin of the fibers combined with random production flaws results in significant uncertainties in the properties of natural fiber reinforced composites. A probabilistic assessment can help to characterize the uncertainties and evaluate the reliability of natural fiber composites , enabling their use in engineering designs . Toward this end, this study aims to quantify the uncertainties in the tensile strength and frequency response of a unidirectional flax/epoxy composite due to the variability of various input parameters, including the fiber material properties and manufacturing flaws. Based on the available data in the l iterature, the non-deter ministic input variables were divided into normal and uniform variables using a statistical test. A computationally efficient response surface approach based on the polynomial chaos expansion was adopted to conduct the uncertainty analysis with multi -type uncertain variables . Moreover, the results were validated by the direct Monte Carlo simulation to demonstrate the accuracy and efficiency of the surrogate model. 2 1. Introduction The use of natural plant fiber reinforcements in polymer composite materials has recently gained attraction due to their promising properties such as high specific stiffness, good acoustic insulation, and v ibration damping, and lower environmental impacts [1,2] . These advantages of natural fibers make them suitabl e for replacing synthetic glass fiber reinforced composites in many nonstructural and semi -structural applications, and such substitutions based on short natural fiber preforms are becoming increasingly common in the transportation and sport goods industri es [3]. Besides, with the aim of utilizing the maximum load -bearing potential of natural fibers, an increasing number of investigations are no w beginning to research into continuous long natural fiber composites [4–6]. However, there are still important issu es that limit natural fiber’s future use in engineering designs, including the di fficulties in predicting structural performance as a result of the natural variability in fiber characteristics [7]. The variations in plant origin, quality of the crop, planting condition, harvesting, extraction and refinement method, and fiber characterizat ion are considered as primary aleatoric uncertainty sources in the fiber properties. The uncertainties in fiber physical and mechanical properties generally lead to significant dispersion in the mechanical and structural performance of the natural fiber re inforced composites. Due to the dependency on a large number of stochastic variables, broad analysis of uncertainty propagation in the material behavior and structural response of composites are quite necessary when natural fiber composites are intended to be used in engineering applications. Uncertainty analysis approaches are useful tools in ord er to probabilistically assess the behavior, evaluate the reliability, and quantify the contributions of specific sources of uncertainty to the overall uncertainty of composite structures . So far, t he uncertainty analysis of synthetic fiber composites has been investig ated by a number of researchers, and various methods have been developed . Probabilistic uncertainty analysis methods, such as traditional Mont Carlo s imulation (MCS) methods, stochastic finite element, and response surface methodology (RSM) have been widel y applied in the analysis and prediction of uncertain response of composite structures [8–14]. The direct MCS is conventionally a robust uncertainty propagation method to evaluate the performance variability of composite structures when there are uncertainties in th e input parameters . Despite the superiority of the MCS in terms of precision, it is practically impossible to use it in computationally intensive problems such as high-fidelity models or nonlinear problems with a large number of uncertai n parameters. Inste ad, surrogate model approaches, 3 which are constructed based on a limited set of actual input/output data points, are a suitable method when dealing with such complex problems. Several methods of application of surrogate models have been reported in the lit erature for evaluation of uncertainties in composite laminates , such as Kriging method [15,16] , radial basis func tion [17], polynomial chaos expansion (PCE) [18–20], and artificial neural network (ANN) [21,22] . State -of-the-art reviews on the surrogate models for evaluating the uncertainty in structural responses of composite laminates can be found in [23]. The response surface methods (RSMs ) are computationally efficient surrogate models that can approximate the response of a system based on the actual response at a limited set of algorithmically chosen points from uncertain inputs. The RSM is ideal for uncertainty propagation of computationally expensive models , such as uncertainties in the structural response of composite s, and evaluating their reliability [18]. It has been shown that the application of the non-intrusive polynomial chaos expansion method is more efficient for dynamic analysis of stochastic systems compared to other RSMs [23,24] . Dey et al. [25] considered fuzzy membership functions to represent the parameter uncertainties and successfully applied the non-intrusive polynomial chaos expansion method to uncertainty propagation in the dynamic character istics of composite structures. Jacquelin et al. [26] extended this method to an uncertain system with mixed uncertainties of random and fuzzy variables and proposed the product of Hermite and Legendre polynomials as a PCE to deal with both types of uncertainties . Recently, Chen et al. [18] proposed a PCE based approach for uncertainty analysis of composite structures considering mixed input uncertainties, includin g normal random and interval variables. In the authors’ w ork, the type of uncertain variables is determined based on the available data for input parameters. In other words, an uncertain parameter is treated as a normal random variable if there exists suff icient data to fit a normal distribution; otherwise, it is considered as an interval variable. A similar approach was also employed in constructing an efficient surrogate model for reliability evaluation of composite structures by considering the uncertain parameters as random or interval variables [27]. A more general mixed uncertainties related to uncertain structural, material and geo metric parameters with multiple uncertainty types were ta king into account for uncertainty quantification of composite laminates by Peng et al. [20]. They divided the uncertain inp ut parameters into different variable types based on the data available for each parameter a nd applied the PCE method to analyze the uncertainties . It can be concluded from the previous studies t he PCE with multi -type uncertainties is an efficient method for uncertainty evaluation of composites . Nevertheless , the validity and 4 effectiveness of this approach when dealing with a large number of input uncertainties from different types have not been fully studied . As described earlier, unlike the synthetic ﬁbre reinforced composites, there are a larger number of uncertainty sources in the constituents and manufacturing process of natural fiber composites. However, there are few studies dedicated to a syst ematic assessm ent of the uncertainty in the properties of natural fiber composites . As one of the few studies in this regard, Blanchard et al. [28] conducted a comparati ve study on the reliabilit y of a flax an d E-glass fiber composite structure using the direct MCS . It was concluded that the structure made of the flax fiber composite needs to be heavier than the E - glass structure to achieve equivalent reliability. Considering the advantages of th e PCE approach, this study aims to provide a computationally efficient framework to quantify the uncertainties in the structural properties of natural fiber composites. Based on the mechanical properties reported in the lite rature for one type of natural fibers, namely flax fibers, and the statistical tests , the uncertain input parameters are divided into two types, strong statistical and interval variables. The rest of this paper is organized as follows. Section 2 describes the numerical procedure for cons tructing the response surface model. Then, the uncertainty analysis problems of the strength and natural frequencies of the unidirectional (UD) flax/epoxy laminate are discussed in Section 3. Finally, the conclusions are pre sented in Section 4. 2. Polynom ial chaos expansion method for uncertainty analysis with multi - type uncertain parameters According to the concept of the polynomial chaos expansion (PCE) introduced by Wiener , the combination of orthogonal polynomial bases is a suitable method for approxim ation of unknown or complex correlations [29]. The PCE is , therefore, a proper choice for estimating the variability of uncertain systems . Based on the simplified form of the PCE, the uncertain response of a system can be expressed as a function of normall y distributed random variables (𝜉𝑖) using the following convergent series [30] 5 𝑋(𝜉1,𝜉2,…,𝜉𝑚) = 𝑐0+∑𝑐𝑖1𝐻1(𝜉𝑖1)𝑚 𝑖1=1 +∑∑𝑐𝑖1𝑖2𝐻2(𝜉𝑖1,𝜉𝑖2)𝑖1 𝑖2=1𝑚 𝑖1=1 +∑∑∑𝑐𝑖1𝑖2𝑖3𝐻3(𝜉𝑖1,𝜉𝑖2,𝜉𝑖3)𝑖2 𝑖3=1+𝑖1 𝑖2=1𝑚 𝑖1=1…. (1) where 𝐻𝑝(𝜉𝑖1,𝜉𝑖2,…,𝜉𝑖𝑝),𝑖1,𝑖2,…,𝑖𝑝∈{1,2,…,𝑚} denotes the Hermite -chaos of order 𝑝, which is a set of multidimensional Hermite polynomials, 𝜉𝑗,𝑗∈{1,2,…,𝑚}, are 𝑚 independent standard normal random input variables, and 𝑐𝑖1,𝑐𝑖2,…,𝑐𝑖𝑝are deterministic constants. The variable 𝜉𝑖 is determined by means of the following change -of-variable 𝜉𝑖=𝑥𝑖−𝜇𝑥𝑖 𝑓𝜎𝑥𝑖 (2) where 𝑥𝑖 is normal variables, 𝜇𝑥𝑖 and 𝜎𝑥𝑖are the corresponding mean and standard deviation, respectively, and 𝑓 is a constant coefficient which is determined based on the confidence level. Here, the constant 𝑓 is set equal to 3 according to the three -sigma rule of the normal distribution . The m ultidimensional Hermite polynomials can be obtained by 𝐻𝑝(𝜉𝑖1,𝜉𝑖2,…,𝜉𝑖𝑝)= (−1)𝑝𝑒 1 2𝝃𝑇𝝃𝜕𝑝𝑒 −1 2𝝃𝑇𝝃 𝜕𝜉𝑖1𝜕𝜉𝑖2… 𝜕𝜉𝑖𝑝 (3) where 𝝃 denotes the 𝑝-dimensional vector of standard normal ran dom variables (𝜉𝑖1,𝜉𝑖2,…,𝜉𝑖𝑝), 𝑖1,𝑖2,…,𝑖𝑝∈{1,2,…,𝑚}. The coefficients 𝑐𝑗 in Eq .(1) are learnable parameters or weights which can be determined using 𝑁 sample points of normal ran dom input variables, i.e. 𝝃𝒊,𝑖∈{1,2,…,𝑁}, and the corresponding actual responses 𝑋(𝝃𝒊),𝑖=1,2,…,𝑁. A suitable sampling strategy, such as the sparse collocation method , can be employed to determine the 𝑁 sample points [31,32] . The response surface model of Eq.(1) can be expressed in the matrix form as 𝑿=𝑯𝑪 (4) Here , 𝑿=𝑿(𝑋1,𝑋2,…,𝑋𝑁)𝑇 is the vector of actual responses for 𝑁 sample points , in which 𝑋𝑖,𝑖∈{1,2,…,𝑁} is an actual outcome obtained through deterministic analysis of the problem with r espect to a sample point of uncertain variables. In Eq. (4), 𝑪 is the polynomial coeﬃcient vector and 𝑯 is the matrix of all Hermite polynomials for 𝑁 sample points which is defined as 6 𝑯= [ 1𝐻1(𝜉11)⋯𝐻2(𝜉𝑖11 ,𝜉𝑖21)…𝐻𝑝(𝜉𝑖11,𝜉𝑖21,…,𝜉𝑖𝑝1)… 1𝐻1(𝜉12)…𝐻2(𝜉𝑖12 ,𝜉𝑖22)…𝐻𝑝(𝜉𝑖12,𝜉𝑖22,…,𝜉𝑖𝑝2)… ⋮⋮⋱ ⋮⋱ ⋱ ⋱ 1𝐻1(𝜉1𝑁)…𝐻2(𝜉𝑖1𝑁 ,𝜉𝑖2𝑁)…𝐻𝑝(𝜉𝑖1𝑁,𝜉𝑖2𝑁,…,𝜉𝑖𝑝𝑁)…] (5) Thus , the polynomial coeﬃcient vector can be determined based on the N input sample points using 𝑪=(𝑯𝑻𝑯)−𝟏𝑯𝑻𝑿 (6) One can make good use of sampling techniques corresponding to the probability distribution of the input variables and numerical matrix inversion to reduce the computational cost of Eq. (6). Besides, f or simplici ty, Eq. (1) can be written as 𝑋(𝝃)=∑𝑐̂𝑗Ψ𝑗(𝝃)𝑀 𝑗=1 (7) where 𝑐̂𝑗 is identical to coefficients 𝑐0,𝑐𝑖1,𝑐𝑖2,…,𝑐𝑖𝑝and there is a one -to-one correspondence between the polynomial functions Ψ𝑗(𝝃) and 𝐻𝑝(𝜉𝑖1,𝜉𝑖2,…,𝜉𝑖𝑝). 𝑀 is the total number of undetermined polynomial coeﬃcients and can be determined , based on the highest order of polynomials ( 𝑝) and number of uncertain variables ( 𝑚), by (𝑚+𝑝)! 𝑚!𝑝!. So far, only norma l random variables have been accounted for in the approximation of the structural response 𝑋 since the random input of Hermite -chaos polynomials of Eq.(1) is limited to normal random variables . The reason is that the Hermite polynomials are uncorrelated when random input ( 𝝃) is normal . This assures the optimal exponential convergence rate of the series of Eq.(7) [30]. However, based on the available data for the uncertain input parameters, there is usually more than one type of random inputs. Here, the random inputs are cate gorized into two types , normal random variables and uniform variables. The classification of input variables is carried out based on the Anderson -Darling (AD) goodness -of-fit test results. An input is considered as a normal random variable if available dat a for that parameter fits a normal distribution ; otherwise, it is considered as a uniform variable. According to the idea of the multi -type random variables PCE [18,26] , the uniform random variables can be incorporated in the polynomial chaos coefficients by expressing them a s a function of the orthogonal Legendre polynomials. This is because the orthogonality weighting functions of t he Legendre polynomials provide a good match with the probability density f unction of the uniform random variables [30]. Considering the Legendre polynomials of order 2, the coefficients 𝑐̂𝑗,𝑗∈{1,2,…,𝑀} in Eq.(7) take the form 7 𝑐̂𝑗=𝑎𝑗0+1 4∑𝑏𝑗𝑙𝛾𝑙𝑛 𝑙=1+1 2∑𝑑𝑗𝑙(3𝛾𝑙2−1)𝑛 𝑙=1 (8) where 𝑎𝑗0,𝑏𝑗𝑙,𝑑𝑗𝑙 are the polynomial coefficients and 𝛾𝑙,𝑙∈{1,2,…,𝑛} are the 𝑛 independent standard uniform random variable, which can be obtained by 𝛾𝑙=𝑦𝑙−𝑦𝑙𝑐 𝑦𝑙𝑟 (9) where 𝑦𝑙𝑐 and 𝑦𝑙𝑟 are the m ean value and the interval radius of the uniform variable, respectively. By substituting Eq. (8) into Eq. (7) , the latter can be expressed as a function of two types of random variables as 𝑋(𝜉1,𝜉2,…,𝜉𝑚,𝛾1,𝛾2,…,𝛾𝑛) =∑(a𝑗0+1 4∑𝑏𝑗𝑙𝛾𝑙𝑛 𝑙=1𝑀 𝑗=1 +1 2∑𝑑𝑗𝑙(3𝛾𝑙2−1)𝑛 𝑙=1)Ψ𝑗(𝜉1,𝜉2,…,𝜉𝑚) (10) To construct the response surface, t he unknown coefficients a𝑗0,𝑏𝑗𝑙 and 𝑑𝑗𝑙 can be determined by solving Eq. (10) by means of an adequate number of sample points . The Latin Hyper cube Sampling (LHS) method [29] and t he sparse Cleanshaw -Curtis (SCC) collocation strategy [31] can be employed to determine the sample points for the standard uniform and normal random variables, respectively . The samp le points should be determined in a two -level sampling strategy. In the inner level, 𝑁 sample points are determined for the normal random variables by means of the SCC collocation strategy as (𝜉1𝑣,𝜉2𝑣,…,𝜉𝑚𝑣),𝑣∈{1,2,…,𝑁}. In the outer level , depending on the accuracy of the polynomial chaos expansion , 𝐾 sample point s are determined for the uniform variables appeared in the polynomial coefficients using the LHS method as (𝛾1𝑤,𝛾2𝑤,…,𝛾𝑚𝑤),𝑤∈{1,2,…,𝐾}. Considering Eq.(4)-(6) and applying the sample points of the uniform variables , the coefficients of the Hermite -chaos polynomials in Eq. (10) can be written as the following matrix form 𝜞𝑬𝑗=𝑪𝑗,𝑗∈{1,2,…,𝑀} (11) where 8 𝜞= [ 1(𝛾1 4)1 ⋯(𝛾𝑛 4)1 (3𝛾12−1 2)1 …(3𝛾𝑛2−1 2)1 1(𝛾1 4)2 …(𝛾𝑛 4)2 (3𝛾12−1 2)2 …(3𝛾𝑛2−1 2)2 ⋮⋮⋱⋮ ⋱ ⋱ ⋱ 1(𝛾1 4)𝐾 …(𝛾𝑛 4)𝐾 (3𝛾12−1 2)𝐾 …(3𝛾𝑛2−1 2)𝐾 ] (12) is the matrix of Legendre polynomials of order 2 for 𝐾 sample points of uniform variables and 𝑬𝑗=[𝑎𝑗0𝑏𝑗1…𝑏𝑗𝑛𝑑𝑗1…𝑑𝑗𝑛]𝑇 (13) 𝑪𝑗=[𝑐̂𝑗1𝑐̂𝑗2…𝑐̂𝑗𝐾]𝑇 (14) Using Eq. (11)-(14), the unknown coefficients 𝑎𝑗0,𝑏𝑗𝑖 and 𝑑𝑗𝑖 of Eq. (10) are calculated by 𝑬𝑗=(𝚪𝑇𝚪)−1𝚪𝑇𝑪𝑗 (15) Consequently, the response surface model of Eq. (10) is obtained. Here, the number of sample points in LHS and the collocation level in t he SCC can be determined based on the accuracy of the developed RSM compared to the actual response of the composite structure. One can make good use of the root mean square errors (RMSE) method to evaluate the precision and the convergence of the surrogat e mode l at any given number of sample points for uniform variables . The RMSE is given by 𝑅𝑀𝑆𝐸=[1 𝑁∑(𝑋𝑖−𝑋𝑖𝑑)2𝑁 𝑖=1]1 2 (16) where 𝑋𝑖 is the predicted response by the RSM and 𝑋𝑖𝑑 is the actual response for an identical input. Clearly, the smaller the value of RMSE , the higher p recision of the RSM prediction . The detailed procedure of solving Eq.(10) can be found in [18]. 3. Uncertainty analysis of mechanical properties of flax fiber com posites by means of PCE -Legendre RS M Flax fiber i s one of the most desirable natural fibers to replace synthetic glass fibers in many applications [33–36]. This is mainly due to the promising characteristics of flax fibers like excellent mechanical pr operties, lower density, and higher aspect ratio. As a result, the mechanical properties of flax fibers and their composites have been extensively reported in the literature. However, their properties have a relatively large variation that arises from their natural variab ilities , such as origin, variety, and maturity of the plant, as well as processing variabilities, such as harvesting, extraction , and refinement of fibers. The uncertainty effects of 9 these factors are generally hard or impossible to control due to their stochastic nature . Thus, it is necessary to consider the uncertainties of properties when modeling and predicting the behavior of flax fiber reinforced composites. The uncertainty information of the elastic properties and density for flax fi bers based on the available experimental data from the literature are summarized in Table 1. One can apply a distributional test to the experimental data obtained from different sources to evaluate data distribution if their probab ility distributions are not specified. Here, t he Anders on-Darling test [37], with a confidence level of 9 5%, is used to test the data for normality. The uncertain parameter is considered as a normal random variable (𝜉) if the calculated observation signiﬁcance level was greater than 0.0 5; otherwise, it is treated as a uniformly distributed variable (𝛾). Besides, when there is insufficient data for an unc ertain parameter, that parameter is assumed as a uniformly dist ributed variable (between its upper and lower limit). Table 1. Statistical information of f lax fiber properties [7,38–52]. Variable Type Mean SD Lower limit Upper limit Young’s modulus, 𝐸𝑓 (GPa) Normal 52.01 17.12 - - Shear modulus, 𝐺𝑓 (GPa) 2.5 0.2 Density, 𝜌𝑓 (g.cm-3) 1.47 0.24 Poisson ’s ratio, 𝜈𝑓 (%) Uniform - - 27 31 The variables related to matrix properties are assumed to be deterministic. Since thermo set resin systems are widely used for flax fiber composites , the average properties of a common epoxy resin with 𝐸𝑚=3.70 GPa, 𝐺𝑚=1.37 GPa, 𝜈𝑚=0.35, and 𝜌𝑚=1.15 g.cm-3 are considered for the analysi s [44]. Considering the uncertainties in the elastic properties and density of flax fibers , the effective elastic properties of the composite can be estimated as a function of the non-deterministic fiber parameters through the rule of mixtures as 𝐸1=𝑉𝑓𝐸𝑓+(1−𝑉𝑓)𝐸𝑚 𝐸2=𝐸𝑚1+0.5𝛿𝑉𝑓 1−𝛿𝑉𝑓, 𝐺12=𝐺𝑚(1+𝑉𝑓)𝐺𝑓+(1−𝑉𝑓)𝐺𝑚 (1−𝑉𝑓)𝐺𝑓+(1+𝑉𝑓)𝐺𝑚 𝑣12=𝑉𝑓𝑣𝑓+(1−𝑉𝑓)𝑣𝑚 (17) 10 where 𝐸1 and 𝐸2 are the effective elastic modulus in the fiber and transverse direction , respectively , 𝑉𝑓 is the fiber volume fraction, 𝐺12 is the in -plane shear modulus , and 𝑣12 is the in-plane Poisson’s ratio. The coefficient 𝛿 is defined as 𝛿=𝐸𝑓−𝐸𝑚 𝐸𝑓+0.5𝐸𝑚 (18) Eq.(17) is based on the assumption that both the materials are linear elastic and the fiber -matrix bonding is perfect [53,54] . Besides, the variability in the local fiber orientation is considered by treating the fiber volume fraction (𝑉𝑓) as a random variable. It is assumed that its distribution is normal , with a mean of 40% and a standard deviation of 11%. In addition to the elastic properties of the fiber, the strength parameters of the UD flax/epoxy laminate are assumed to be nondeterministic and treated as random variables . The statistical information of the strength variables of UD flax/epoxy composite s obtained from the literature is given in Table 2. Table 2. Statistical information of strength parameters for the UD flax/epoxy composite laminate [7,38,39,43 – 48,50–52,55–64]. Variable Type Mean SD Lower limit Upper limit Longitudinal tensil e, 𝑋𝑡 (MPa ) Normal 291.4 84.8 - - Transverse tensil e, 𝑌𝑡 (MPa) 25.9 13.2 In-plane shear, 𝑆 (MPa) 37.7 5.7 Longitudinal compression, 𝑋𝑐 (MPa) Uniform - - 110 136.9 Transverse compression, 𝑌𝑐 (MPa) 76 100 Another uncertainty source is the misalignment of ply o rientation (𝜃) that incurred during the layup process. In the uncertainty analysis of the composite, t he orientation of each ply is considered as an independent normal random variable as 𝜃𝑖 ~ 𝒩(𝜇𝜃𝑖 ,𝜎𝜃𝑖2),𝑖∈{1,2,…,𝑁} (19) where 𝜇𝜃 is the nominal pl y orientation , 𝜎𝜃 is the standard deviation of ply orientation, and 𝑁 is the number of plies. In contrast to the assumption of deterministic values for the standard deviation of material properties, the ply orientation variability ( 𝜎𝜃) is non -deterministic and subject to uncertainties due to uncertainties in the processing of fibers as well as uncertainties during the lay -up. Hence, a parametric sensitivity analysis concerning the dispersion of the ply orientation is required. For this purpose, t he uncertainty analysis is conducted using four 11 different values of 𝜋/360, 𝜋/180, 𝜋/90, and 𝜋/60 for the standard deviation of the ply orientation. In the subsequent sections, the resultant uncertainties in the failure load and natural frequencies of a UD flax/epoxy composite plate are evaluated , and t he validity of the PCE -Legendre method in assessing the uncertainties in the material response of natural fiber composites with high-dimensional uncertain input par ameters are demonstrated. The analysis i s performed on a desktop computer with an Intel CoreTM i7-4790 CPU and 16 GB RAM. 3.1. Prediction of fail ure load of a UD flax /epoxy composite laminate As mentioned in the previous section, there are considerable uncertainties in the primitive variables of flax fiber composites. It is necessary to consider the uncertainty of these variables when simulating the response of a flax fiber composite. In this s ection, the probability distribution of the first-ply failure load of a flax/epoxy composite laminate is predicted using the PCE-Legendre based RSM . The lay-up of the laminate is [0]12 and it is subjected to the tensile load 𝑁𝑥 per unit length, as shown in Figure 1. Figure 1. UD flax/epoxy lamin ate and the loading direction Considering the classical laminat e theory and the plane stress assumption, the deterministic stresses in the 𝑘𝑡ℎ ply with a misalignment of 𝜃𝑘 can be calculated as follows {𝜎𝑥 𝜎𝑦 𝜎𝑥𝑦}(𝑘) =[𝑄̅11𝑄̅12𝑄̅16 𝑄̅12𝑄̅22𝑄̅26 𝑄̅16𝑄̅26𝑄̅66](𝑘) {𝜀𝑥0+𝑧𝜅𝑥 𝜀𝑦0+𝑧𝜅𝑦 𝛾𝑥𝑦0+𝑧𝜅𝑥𝑦} (20) where 𝑄̅𝑖𝑗 are the components of the transformed lamina stiffness matrix of the 𝑘𝑡ℎ ply, and (𝜀𝑥0𝜀𝑦0𝛾𝑥𝑦0)𝑇and (𝜅𝑥𝜅𝑦𝜅𝑥𝑦)𝑇are the in-plane strains and curvatures associated with the middle surfac e of the laminate, respectively . The strain and cur vature components of the middle surface can be obtained based on the tensile load 𝑁𝑥 as 𝜃𝑥 𝑦 𝑧 12 { 𝜀𝑥0 𝜀𝑦0 𝛾𝑥𝑦0 𝜅𝑥 𝜅𝑦 𝜅𝑥𝑦} = [ 𝐴11𝐴12𝐴16𝐵11𝐵12𝐵16 𝐴12𝐴22𝐴26𝐵12𝐵22𝐵26 𝐴16𝐴26𝐴66𝐵16𝐵26𝐵66 𝐵11𝐵12𝐵16𝐷11𝐷12𝐷16 𝐵12𝐵22𝐵26𝐷12𝐷22𝐷26 𝐵16𝐵26𝐵66𝐷16𝐷26𝐷66] −1 { 𝑁𝑥 0 0 0 0 0} (21) where the laminate coefficients 𝐴𝑖𝑗, 𝐵𝑖𝑗, and 𝐷𝑖𝑗 are defined in terms of the transformed lamina stiffnesses 𝑄̅𝑖𝑗 as (𝐴𝑖𝑗,𝐵𝑖𝑗,𝐷𝑖𝑗)=∑∫𝑄̅𝑖𝑗(𝑘)(1,𝑧,𝑧2)𝑑𝑧𝑧𝑘+1 𝑧𝑘𝑁 𝑘=1 (22) where 𝑁 is the number of plies. The components of the transformed lamina stiffness matrix are obtained as a function of the lamina stiffness es matrix and the ply orientation as 𝑄̅11(𝑘)=𝑄11cos4𝜃𝑘+2(𝑄12+2𝑄66)sin2𝜃𝑘cos2𝜃𝑘+𝑄22sin4𝜃𝑘 𝑄̅12(𝑘)=(𝑄11+𝑄22−4𝑄66)sin2𝜃𝑘cos2𝜃𝑘+𝑄12(sin4𝜃𝑘+cos4𝜃𝑘) 𝑄̅22(𝑘)=𝑄11sin4𝜃𝑘+2(𝑄12+2𝑄66)sin2𝜃𝑘cos2𝜃𝑘+𝑄22cos4𝜃𝑘 𝑄̅16(𝑘)=(𝑄11−𝑄12−2𝑄66)sin𝜃𝑘cos3𝜃𝑘 +(𝑄12−𝑄22+2𝑄66)sin3𝜃𝑘cos𝜃𝑘 𝑄̅26(𝑘)=(𝑄11−𝑄12−2𝑄66)sin3𝜃𝑘cos𝜃𝑘 +(𝑄12−𝑄22+2𝑄66)sin𝜃𝑘cos3𝜃𝑘 𝑄̅66(𝑘)=(𝑄11+𝑄22−2𝑄12−2𝑄66)sin2𝜃𝑘cos2𝜃𝑘+𝑄66(sin4𝜃𝑘 +cos4𝜃𝑘) (23) The components of the lamina stiffness matrix , 𝑄𝑖𝑗, are related to the uncertain effective elastic properties of Eq. (17) by 𝑄11= 𝐸1 1−𝑣12𝑣21 𝑄12= 𝑣12𝐸2 1−𝑣12𝑣21=𝑣21𝐸1 1−𝑣12𝑣21 𝑄22= 𝐸2 1−𝑣12𝑣21 𝑄66=𝐺12 𝑣12 𝐸1=𝑣21 𝐸2 (24) One can plug Eq. (21)-(24) into Eq. (20) and use a failure criterion to predict the value of 𝑁𝑥 corresponding to the first -ply failure of the laminate . Here, the Tsai -Wu criterion is applied to 13 determine the failure of a lamina. According to this criterion, ply f ailure occurs when the following inequality is sati sfied 𝐹1𝜎1+𝐹2𝜎2+𝐹11𝜎2 1+𝐹22𝜎2 2+𝐹66𝜎2 12+2𝐹12𝜎1𝜎2≥1 (25) where the strength coeﬃcients are defined in terms of the lamina strengths as 𝐹1=1 𝑋𝑡−1 𝑋𝑐,𝐹2=1 𝑌𝑡−1 𝑌𝑐,𝐹11=1 𝑋𝑡𝑋𝑐,𝐹22=1 𝑌𝑡𝑌𝑐,𝐹12 =−1 2√𝐹11𝐹22,𝐹66=1 𝑆2 (26) in which the lamina strengths are random variables, as given in Table 2. For a given set of input parameters , the nonlinear Eq. (25) can be solved for 𝑁𝑥 to ob tain the failure load of each lamina , where t he least value corresponds to the first -ply fail ure load. The PCE-Legendre method described in section 2 is applied to construct the response surface for 𝑁𝑥 in MATLAB, considering all 23 uncertain uniform an d normal variables, including fiber elastic properties ( Table 1), lamina strengths ( Table 2), and the misalignment of ply orienta tions , i.e. 𝜃𝑖,𝑖∈{1,2,…,12}. Based on the level of accuracy, the number of sample points is determined and applied to solve the Eq. (25) to obtain the corresponding outputs. In this problem, the required accuracy was achieved when the order of Hermite polynomial chaos (𝑝) and the collocation level were equal to 4 and 2, respectively . The distribution s of the first -ply fai lure load when the standard deviation of the ply orientation (𝜎𝜃) is 𝜋/180, obtained from the direct MCS and the PCE-Legendre based RSM are compared in Figure 2. Based on the central limit theorem, it can be concluded that th e variation of 𝑁𝑥 due to the multi -type uncertainties in elastic properties, strengths, and the ply orientation is normal. This enable s us to use the mean value and the standard deviation to describe the 𝑁𝑥. As can be seen, the distribution s of 𝑁𝑥 obtained from the direct MCS and the PCE -Legendre method are in an excellent correlation . It should be noted, h owever, the required number s of solving the nonlinear inequality of Eq. (25) to obtain 𝑁𝑥 are considerabl y different in each method. T he number of sample points required to construct the response surface model using the PCE-Legendre is 3065, while the dir ect MCS needs at least 40000 sample points to reach the desired convergence. Therefore , the computational efficiency of the PCE-Legendre RSM is significantly higher than the direct MCS . 14 Figure 2. Probability density function of the first-ply failure load ( 𝑁𝑥). The number of sample points and the corresponding mean s and standard deviation s (SD) of 𝑁𝑥 for four different values of 𝜎𝜃 are compared in Table 3. Although the difference between the results of the RSM and MCS method s is increased by increasing the standard deviation of the ply orientatio n, the magnitude of the relative errors , particularly for the mean values, is not significant. Consequently, c onsidering the number of uncertain input parameters, it can be said that the PCE-Legendre response surface model is able to represent the uncertai nties of the failure load of the flax fiber composite to a good accuracy , with much less computational cost. Due to the lack of knowledge about the degree of the ply misalignment, one can repeat the uncertainty analysis of 𝑁𝑥 for various values of 𝜎𝜃 to evaluate its effects on the strength of the laminate. The sensitivity of the mean value and standard deviation of 𝑁𝑥 to the ply misalignment dispersion using the PCE-Legendre method is shown in Figure 3a. The results of the direct MCS are also given for validation. Both the mean value and the standard deviation of 𝑁𝑥 show a dependence on 𝜎𝜃, registering a reduction of up to 4.5% and 8.2%, respectively, when 𝜎𝜃 is increased from 0.5 to 3 degrees . The red uction predicted by the direct MCS is slightly less than that of the PCE-Legendre method – 4.1% and 6.4% in the mean and the standard devi ation of 𝑁𝑥, respectively . An interesting fact seen in Figure 3a is that the standard de viation of 𝑁𝑥 decreases by increasing the dispersion of the ply misalignment . This means that the uncertainty of 𝑁𝑥 is inversely related to the uncertainty of the ply orientations. It must be noted that this correlation is only observed when all the uncertainties (i.e . fiber properties, lamina strengths, and misalignment of ply 15 orientations ), are consider ed in the analysis . In other words, an opposite behavior (increasing the standard deviation of 𝑁𝑥by increasing 𝜎𝜃) is obtained when the uncertain inputs are limited to ply orientations. This can be clearly seen from Figure 3b in which all input parameters are deterministic except the ply orientations. The mean values are used for t he deterministic inputs. Table 3. Uncertainty analysis results of the first -ply failure load ( 𝑁𝑥) for different standard deviation of ply orientation using RSM and MCS. 𝜎𝜃 PCE -Legendre based RSM MCS Error (%) No. of simulations Mean (N/mm) SD (N/mm) No. of simulations Mean (N/mm) SD (N/mm) Mean SD 𝜋 360 3065 1105.87 314.64 30000 1106.1 320.23 0.02 1.75 𝜋 180 1092.53 309.21 40000 1097.04 316.03 0.41 2.16 𝜋 90 1076.92 300.1 60000 1080.48 309.66 0.33 3.09 𝜋 60 1055.89 288.75 70000 1060.94 299.75 0.48 3.67 Figure 3. Mean and standard deviation of 𝑁𝑥 against the standard deviation of ply orientation considering : (a) all uncertainties and (b) only uncertainties in ply orientations. 3.2. Uncert ainty analysis of natural frequencies of a UD flax /epoxy composite laminate In the second case study , the effect of uncertainties on the natural frequencies of a [0]12 flax/epoxy laminate is investigated . The fiber material properties ( Table 1), the fiber volume fraction, and the misalignment of ply orientations are considered as the stochastic input (a) (b) 16 variables . In this problem , the only uniform variable is the Poisson’s ratio of the fiber, 𝜈𝑓, and all other inputs are norm ally distributed. The effective elastic properties of the individual lamina are obtained using Eq. (17). Thus , the natural frequency problem can be stated as 𝜔𝑛=𝜔𝑛(𝐸𝑓 ,𝐺𝑓 ,𝑉𝑓 ,𝑣𝑓 ,𝜌𝑓 ,𝜃1 ,𝜃2 ,….,𝜃12) (27) where 𝜔𝑛 is the nth natural frequen cy. To obtain the frequency response of the composite, a n 11mm × 11mm laminate with 12 layers is modeled using the conventional shell element in Abaqus/Explic it. After a mesh convergence study, t he finite element analysis (FEA) results using the mean values as inputs were compared with experimental data from the literature [65] to validate the model. Table 4 compares the first four natural frequencies obtained from the FEA and experimental tests o f a [0]12 flax/epoxy laminate . The comparison of FEA and experimental results indicates that the model is able to capture the frequency response of the composite with acceptable accuracy . Table 4. FEA and experime ntal results of the natural frequency analysis of the flax/epoxy laminate. Mode No. Experimental (Hz) FEM (Hz) Error (%) 1 58.82 59.87 1.79 2 125.7 127.69 1.58 3 377.96 371.67 1.66 4 481.35 488.30 1.44 A Python script in conjunction with a MATLAB co de is developed to construct the RSM following the approach described for the PCE-Legendre method and based on the specified accuracy level, i.e. the order of Hermite polynomial chaos and the collocation level. In this fully automatic procedure, t he sample points of the uniform and normal parameters are generated by the MATLAB code employing the LHS and SCC sampling strategies . Then, t he Python script runs the FEA in Abaqus for each set of the sample points and save the results for post-processing in the MA TLAB. Finally, the RSM is constructed by calculating the coefficients of Eq.(10) and the sensitivity analysis is implemented. This model is then used for the uncertainty quantification. In order to conduct the direct MCS, a Pytho n script is written capable of generat ing stochastic input parameters, running the FEA , and performing statist ical analysis on the outputs . In this uncertainty analysis , the order of Hermite polynomial chaos and the collocation level are considered equal t o 4 and 2 , respectively , to achieve the desired precision. 17 The results of the uncertainty analysis of the firs t four natural frequencies are given in Table 5. The results of the direct MCS are also presented for validation purposes . As can be seen, the results of the PCE-Legendre method are in good agreement with the direct MCS results. The relativ e errors of the mean values and standard deviations calculated by the PCE-Legendre response surface model are less than 0.15% and 1.5%, r espectively. Figure 4 compares t he cumulative distribution functions (CDFs) of the natural frequencies obtained from both the methods, which confirms that the results of the RSM are highly correlated to the direct MCS results. Table 5. Uncertainty analysis results of the first four natural frequencies using the RSM and MCS. Mode No. PCE -Legendre based RSM (Hz) MCS (Hz) Error (%) Mean SD Mean SD Mean SD 1 52.38 9.88 52.46 9.72 0.14 1.51 2 122.4 9.34 122.4 3 9.36 0.02 0.15 3 326.12 60.54 325.74 60.45 0.12 0.14 4 454.68 50.86 455.3 9 50.31 0.15 1.09 Figure 4. Probability of natural frequencies, (a) mode 1 and 2, (b) mode 3 and 4. To demonstrate the efficiency of the RSM, the computation time of the uncertainty analysis using the direct MCS and the PCE -Legendre response surface are calculated and compared in Figure 5. The initial computation time of the RSM indicates the total time of running th e FEA (a) (b) 18 of the sample points to construct the surrogate model. Based on the chosen order of Hermite polynomial chaos, the collocation level , and the LHS in this problem, the number of sample points is 1635 and the corre sponding initial computation time is 369.7 minutes. The difference between the computation time of the MCS and PCE-Legendre method is more noticeable as the number of simulations increases. Therefore , the PCE-Legendre method is significantly more efficient than the direct MCS in terms of the co mputational cost, particularly for high dimension uncertainties where the convergence rate of the MCS is very slow. Figure 5. Computation time of the uncertainty analysis of frequency response using the MCS and the RSM. As men tioned in the previous section, the sensitivity of the uncertainty analysis to the dispersion of the ply misalignment needs to be evaluated . Thus, the influence of the standard deviation of the ply orientation on the frequency response is investigated by m eans of the PCE - Legendre RSM . Figure 6 illustrates the variation of mean value s and stand ard deviation s of the first four natural frequencies with respect to various standard deviation s of the ply orientation (𝜎𝑓), ranging fro m 0.5 to 3 degree s. The other uncertainties (i.e. material properties) are the same as discussed earlier. The uncertainties analysis results show that the mean and standard deviation of all the first four natural frequencies remain almost unchanged with increasing the standard deviation of the ply orientation up to 3 degree s. In other words, considering all the so urces of the uncertainty (i.e. fiber material properties, fiber volume fraction, and ply orientation) , the effect of small misalignment of ply or ientations on the frequency response is negligible. This can be attributed to the high uncertainty in the fiber properties compared to the uncertainty considered for the 19 ply orientations. Therefore, it can be concluded that the frequency response of the UD flax/epoxy laminate is not sensitive to small misalignments of the ply orientation. Figure 6. Mean and standard deviation of the first four natural frequencies against the standard deviation of ply orientation. 4. Conclusio n The natural origin of fibers and unavoidable defects during the extraction stage , combined with random manufacturing flaws , have always given rise to considerable uncertainties in the properties of natural fiber reinforced composites. Therefore, the pre diction of the uncertainty of the mechanical and structural behavior of natural fiber composites is crucial for the design process. Mode No. 1 Mode No. 2 Mode No. 3 Mode No. 4 20 This work presents an uncertainty analysis of tensile strength and frequency response of a UD flax/epoxy composite, as one of the widely studied natural fiber composites . The variability of the fiber elastic properties, fiber density, lamina strengths, fiber volume fraction, and the misalignment of ply orientations are considered as the uncertainty sources . The uncertain input variables are divided into normal and uniform variables by means of the Anderson -Darling test , based on the various experimental data acquired from the literature. Due to the high number of uncertain input variables, a computationally efficient RSM based on the PCE is adopted and modified accordingly to take the multi -type stochasticity of the inputs into account. For each case study, a direct MCS is performed to validate the accuracy of the stochastic model. In both cases, with a different number of non -deterministic input var iables, the results of PCE- Legendre RSM demonstrate a strong agreement with the direct MCS results. Therefore, it is concluded that the modified RSM is capable of capturing the stochastic mechanical behavior of natural fiber composite s to good accuracy, wi th significantly less computational effort in comparison with the direct MCS. Besides, a parametric sensitivity analysis with respect to the dispersion of the ply misalignment is conducted to identify the effect of uncertainties intro duced during the lay -up. It is found that the mean and standard deviation of the maximum tensile force decrease with increasing the standard dev iation of the ply orientation . However, the frequency response shows almost no dependence on the standard deviat ion of the ply orientation up to 3 degrees . The uncertainty analysis of natural fiber composites can help to characterize the tails of the distribution of mechanical properties, allowing evaluation of the ir reliability and enabling their use in engin eering designs.