Induced uncertainties during the filament winding (FW) process may cause a significant stochastic variation in the mechanical behaviour of composite shells. This paper aims to develop a novel and deep uncertainty quantification (UQ), sensitivity and reliability analyses of filament wound shells co nsidering manufacturing uncertainties. Firstly, a progressive damage analysis is performed to estimate their deterministic burst pressure. Then, a signal -to-noise (SNR) approach is employed using the Taguchi method for sensitivity analysis and screening un certainties arising from manufacturing. Initial results reveal that the shells are more sensitive to thickness uncertainties for thinner structures. Then, probabilistic and reliability analyses are carried out using the Boosted Decision Trees Regression (BDTR) approach. Despite the complexity and non -linear relationships in the problem, the developed BDTR -based metamodel shows powerful predictive performance. A comparative study shows that ply thickness uncertainty leads to a significant underestimation of failure probability. For expensive and time-consuming models in that only a few runs can be affordable, a modified approximation method for reliability analysis is proposed. Results indicate a high capability at estimating failure probability with high acc uracy. Keywords: Uncertainty quantification; Taguchi approach; Stochastic modelling; Cylindrical shells; Probabilistic analysis 1. Introduction Due to the well -known advantages and unique characteristics of composite materials, filament -wound cylindrical shells are increasingly being utilised as load -bearing parts in numerous fields, such as aeronautical, aerospace, energy, and marine structures. Fibre -reinforced composites usually have complex microstructures and their manufacturing processes are not triv ial, which naturally makes them susceptible to variations in their microstructure, geometry and material properties . These variations can * Corresponding author. Email address: humberto.almeida@qub.ac.uk 2 be interpreted as uncertainties, which may be the cause of variation in their mechanical properties [1]. Generally, uncertainties can be categori sed into two types: (i) aleator ic, which arises from inherently random effects and (ii) epistemic , which is due to the lack of information and knowledge in any activity or phase of the modelling process because of ignorance of the environment and system variables [2]. If these are not properly accounted for, unexpected failure might take place. Conventional deterministic approaches to analys ing and design ing composite structures may undermine their lightweight potential, which leads to conservative design and the use of higher safety factors. Therefore, the consideration of these uncertainties in the design and analysis of composite structures is vital to take full advantage of their high lightweight potential [3]. It is unfeasible to fully control all parameter s during the manufacturing of composites [2], particularly, in the FW process, in which fibre volume fraction (V f), thickness, and winding angle might vary [4]. As a consequence, these parameters can vary in a stochastic m anner, which makes deterministic approaches non -suitable to accurately describe such structures. For instance, Rafiee et al. [5] consider ed a uniform distribution for winding angle and V f as random variables for composite cylinders to stochastically study their functional failure. Azizian and Almeida [4] proposed a framework to capture various multiscale uncertainties for efficient and fast stochastic, probabilisti c and reliability analyses of filament -wound composite tubes. They concluded that ply thickness and winding angle are the most influential factors on the stochastic burst pressure of composite tubes. Progressive damage analysis of composite structures cons idering uncertainties is a complex and computationally expensive process [2, 4, 6-9]. The majority of studies on this topic [7, 8, 10] uses Monte Carlo (MC) method for sampling purposes. However, when the failure probability is low and for complex systems with large variab ilities in the governing variables [2, 4], MC requires a large number of samples to predict failure probability accurately [11]. Consequently, progressive damage analysis of composite structures along with MC can be expensive and challenging. Towards overcoming this issue, metamodels (or surrogate models) are an efficient alternative [12], where the output is only assessed for a subspace of algorithmically chosen inputs and then an equivalent model is built up to mimic the adjacent mapping of the input -to/output system [13]. Several studies [4, 13-15] on composite structures indicate the great potential of metamodels to decrease the computational cost of analysis . Another way to estimate the input/output relationship and the level of sensitivity of each parameter is through a sensitiv ity analysis (SA), which aims at quantifying the relative importance of input parameter(s). Methods of SA fall into two categories [16, 17]: Local and global SA. Local SA focuses on the local impact of input parameters. It computes the gradient of the response concerning its parameters around a nominal value. Global SA quantifies the output uncertainty due to the uncertainty in the input 3 parameters, either considering them individually or iteratively. In this way, uncertainties can be ranked in order of importance on the targeted output , which makes SA essentia l to save computational efforts . Conceição António and Hoffbauer [18] developed approximation models for the reliability analysis of composite laminates. Artificial neural networks (ANN) have been used along with genetic algorithms (GA) for structural reliability analysis and to study the uncertainty propagation of mechanical properties on the response of composite laminate structures under an imposed reliability level. The ANN was later coupled with a Monte Carlo procedure and the variability of the output was carried out usi ng a global sensitivity analysis based on Sobol indices. Ellul and Camilleri [19] developed probabilistic progressive damage modelling for filament wound cylindrical pressure vessels under internal pressure, in which the material properties were considered as the sources of uncertainties. They c oncluded that more efficient and lighter vessels could be designed with their approach. Zhou et al. [20] proposed an adaptive Kriging metamodeling approach for the failure probability, and local and global sensitivity of compos ite radome structures. There are numerous SA approaches available and each has its benefits and drawbacks [21, 22]. Nevertheless, these studies report issues such as complexity, high computational time , and the absence of synergism among input parameters . In this context, this work aims at carryi ng out an uncertainty quantification for composite tubes considering winding angle and thickness as the sources of uncertainties. Screening and sensitivity analys es are carried out using the Taguchi approach. In addition, efficient probabilistic, stochastic and reliability analyses are performed considering these manufacturing uncertainties that arise from the filament winding process. 2. Deterministic progressive damage model Azizian et al. [23] proposed a p rogressive damage model for estimating the burst pressure of filament - wound tubes using the approach proposed by Lapczyk and Hurtado [24] for composite tubes following the ASTM D1599 standard. The deterministic material properties are presented in Tables 1 -2. The investigated carbon fibre/epoxy tubes are subjected to uniform internal pressure with restrained -end [25] boundary conditions. The structure is 660 mm long, with a 50 -mm radius, and each ply is 0.25 -mm thick. The laminate is made of four plies with windi ng angles at ±55°, which is well -known as the optimum winding angle for these boundary conditions [26]. Abaqus finite element (FE) platform has been used to perform numerical analyses. Table 1. Experimentally -measured material properties utilised in the computational analyses [27]. 4 Elastic constants value Strengths Value 𝐸𝑙𝑎𝑠𝑡𝑖𝑐 𝑚𝑜𝑑𝑢𝑙𝑢𝑠 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑖𝑛𝑎𝑙 129.30 GPa 𝑇𝑒𝑛𝑠𝑖𝑙𝑒 𝑠𝑡𝑟𝑒𝑛𝑔𝑡 ℎ 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑖𝑛𝑎𝑙 1409.9 MPa 𝐸𝑙𝑎𝑠𝑡𝑖𝑐 𝑚𝑜𝑑𝑢𝑙𝑢𝑠 Transverse 9.11 GPa 𝑇𝑒𝑛𝑠𝑖𝑙𝑒 𝑠𝑡𝑟𝑒𝑛𝑔𝑡 ℎ 𝑇𝑟𝑎𝑛𝑠𝑣𝑒𝑟𝑠𝑒 764.1 MPa 𝑃𝑜𝑖𝑠𝑠𝑜𝑛 ’𝑠 𝑟𝑎𝑡𝑖𝑜 Plane 1−2 0.32 𝐶𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 𝑠𝑡𝑟𝑒𝑛𝑔𝑡 ℎ 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑖𝑛𝑎𝑙 42.5 MPa 𝑆ℎ𝑒𝑎𝑟 𝑚𝑜𝑑𝑢𝑙𝑢𝑠 in−plane 5.44 GPa 𝐶𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 𝑠𝑡𝑟𝑒𝑛𝑔𝑡 ℎTransverse 134.5 MPa 𝑆ℎ𝑒𝑎𝑟 𝑚𝑜𝑑𝑢𝑙𝑢𝑠 Transverse 2.10 GPa 𝑆ℎ𝑒𝑎𝑟 𝑠𝑡𝑟𝑒𝑛𝑔𝑡 ℎ in−plane 68.9 MPa Table 2. Damage evolution and stabilization values used as inputs in the numerical models. 𝐹𝑟𝑎𝑐𝑡𝑢𝑟𝑒 𝑒𝑛𝑒𝑟𝑔𝑦 𝑇𝑟𝑎𝑛𝑠𝑣𝑒𝑟𝑠𝑒𝑇𝑒𝑛𝑠𝑖𝑙𝑒 /𝑐𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 1.6 N/mm [28] 𝐹𝑟𝑎𝑐𝑡𝑢𝑟𝑒 𝑒𝑛𝑒𝑟𝑔𝑦 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑖𝑛𝑎𝑙𝑇𝑒𝑛𝑠𝑖𝑙𝑒 /𝑐𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 22.5 N/mm [28] 𝑉𝑖𝑠𝑐𝑜𝑠𝑖𝑡𝑦 𝑐𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 𝑇𝑟𝑎𝑛𝑠𝑣𝑒𝑟𝑠𝑒𝑇𝑒𝑛𝑠𝑖𝑙𝑒 /𝑐𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 0.0005 𝑉𝑖𝑠𝑐𝑜𝑠𝑖𝑡𝑦 𝑐𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑖𝑛𝑎𝑙𝑇𝑒𝑛𝑠𝑖𝑙𝑒 /𝑐𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 0.00005 3. Probabilistic analysis The probabilistic analysis herein proposed aims to consider usual manufacturing inconsistencies during the winding of composite tubes, namely the V f and winding angle [29]. Uncertainties related to Vf reflect on thickness variations [4]. Azizian and Almeida Jr [4] concluded that manufacturing uncertainties are the most influential uncertainties during the filament winding of c omposite tubes. Hence, this study focusses on quantifying the uncertainties associated with the thickness and winding angle of filament wound tubes. The structure under investigation has four plies with winding angles of ±55°. Each macro -layer is composed of two sub -layers consist ing of winding angles +𝜃 and −𝜃, which are here treated as random variables. In addition, the thickness of all plies is randomly selected, and they are sorted in descending order and assigned to plies from the innermost ply to th e outermost one. This strategy is purposeful and intelligent because it is planned in conformity with the FW process, in which excess resin moves from the innermost layer to the surface of filament wound composite cylindrical products [29]. For a normal distribution, 99.73% of data lies in this interval [ 𝜇-3𝜎, 𝜇+3𝜎] [30], in which 𝜎 and 𝜇 are standard deviation and mean, respectively. These intervals are used for the first, second, third and fourth sub -plies, respectively, as follows: [μ -3σ, μ], [μ -2σ, μ+σ], [μ -σ, μ+2σ] and [μ, μ+3σ]. This assumption abides by three important subjects: (i) the sum of the mean of these intervals is equal to the sum of the mean of considering four thickness variables with this interval [μ -3σ, μ+3σ]; (ii) the wall thickness of the shell 5 does not exceed the largest and lowest values; and (iii) all intervals involve the mean /deterministic value (i.e. 0.25 mm), where their statistical properties are shown in Table 3 . Table 3 . Statistical details of the considered input parameters. Random variables Mean COV Probability distribution model Winding angle (+𝜃/−𝜃) +55°/ −55° 2% [4] Normal [4] Thickness (mm) 0.25 5% [19] Normal [4, 31] 3.1 Taguchi -based screening and sensitivity analysis The Taguchi approach [32] is used here for screening and SA of the random variables. Taguchi approach is efficient at providing robust design solutions and at low computational efforts via extracting more quantitative data from fewer samples. The main advantage of the Taguchi appro ach over the other ones, such as the one-factor -at-a-time method [33] (OFAT ) is that numerous factors can be simultaneously considered. Consequently, it can capture synergistic effects of parameters. A few SA [34, 35] have been carried out using Taguchi. They used per cent contribution (PC) from analysis of variance (ANOVA) and signal –to–noise ratio (SNR [36]) methods. The SNR was used to show the importance of factors along with PC. However, the current study presents a new strategy, as described next. Firstly, for all random variables ( Table 3 ), levels of the Taguchi method are calculated according to this interval: [μ, μ+3 𝝈]. The lower limit of interval (μ) is level 1 and μ+3𝝈 is level 2. By using the 2 - level Taguchi approach available in Minitab statistical software, a virtual test is designed with an orthogonal array L4 , that is, 2 levels and 2 factors [32]. The response is Equation 1. This equation demonstrates that with the variation of random parameters from its mean μ to μ+3 𝝈, the relative change to its deterministic value can be provided. This strategy also determines sources of uncertainties with significant impacts on scattering the response and can be used for screening purposes to reduce the number of parameters in UQ. 𝑇𝑎𝑔𝑢𝑐 ℎ𝑖 𝑟𝑒𝑠𝑝𝑜𝑛𝑠𝑒 =|𝑃𝐵𝑢𝑟𝑠𝑡𝑇𝑎𝑔𝑢𝑐 ℎ𝑖−𝑃𝐵𝑢𝑟𝑠𝑡𝐷𝑒𝑡𝑒𝑟𝑚𝑖𝑛𝑖𝑠𝑡𝑖𝑐 | (1) Then, progressive damage modelling is carried out according to suggested trials by Taguchi. After that, the analysis is accomplished with “ the Smaller is better ” [37] or “the Larger is better ” [37] options. Since the smaller scattering from the deterministic response is favourable, the “ Smaller is better ” option is chosen here. Then, the response table for SNR can be obtained. Table 4 shows the SNR response table for a carbon fibre/epoxy composite tube with a 0.25 mm ply thickness and deterministic burst pressure of 14.7 MPa. 6 Table 4 . Response table for SNR ( ply thickness is 0.25 mm). Level Thickness Winding angle 1 −4.802 −3.807 2 −7.088 −7.226 Delta 3.006 3.420 Rank 2 1 SNR combines mean and variance, in which higher values indicate a higher impact on the response. Delta level evaluates the relative effect of each factor on the response by calculating the difference between the highest and lowest SNRs for each factor [38]. Based on these Delta values, Minitab assigns a Rank for each factor. Rank 1 indicates the highest Delta value; consequently, the most effective factor on the response [38], Rank 2 shows the second highest Delta value, and so on. For entering these efficient criteria in SA and screening, a new index is herein introduced, which this study names as Delta contribution index (DCI). For a problem with 𝑛 factors: 𝐷𝐶𝐼 𝑖(%)=𝐷𝑒𝑙𝑡𝑎 𝑖 ∑ 𝐷𝑒𝑙𝑡𝑎 𝑖𝑛 𝑖=1×100 (2) Therefore, for the above problem, DCI for thickness and winding angle are 46.78% and 53.22%, respectively, which were calculated through Equation 2 using data from Table 4 . 3.2. Reliability analysis based on machine learning For a particular problem, the sampling method and surrogate modelling techniques should be selected depending on the computational efficiency, desired level of accuracy, nature and number (dimension) of input parameters, presence of noise in sampling data and complexity of the model [13]. Before using a specific surrogate modelling technique, it is necessary to check thoroughly the fitting quality and prediction capability [39, 40]. To have an efficient and high -accur acy metamodel, a thorough investigation was carried out on available techniques [4, 13, 39, 40]. Figure 1 shows the proposed workflow for BDTR -based Monte Carlo reliability analysis herein utilised. The first step lies in constructing a metamodel based on response surface methodology (RSM) [4]. If the built model does not meet the desired level of accuracy, then the work continues using BDTR -based models, which have shown high capability at constructing surrogates from complex models [4, 13, 14]. Since the RSM technique tends to disperse sample points around the boundaries of the design space and place a few 7 points at its centre, space -filling sampling methods are used to reinforce the dataset. Several space -filling sampling methods can cover a multivariable space, such as generations of Halton and Sobol sequ ences, Latin Hypercube sampling (LHS), and the inverse transform sampling method (ITM) [41]. Here, the ITM method is employed due to its simplicity in covering the whole design space. After an initial assessment, the Boosted Decision Tree Regression (BDTR) [42] was chosen for this study. Despite its powerful predictive performance, the BDTR method has been underexploited for UQ and reliability analysis of composite structures. Figure 1. Workflow for BDTR -based Monte Carlo reliability analysis. The BDTR method has some important advantages leading to superiority over most traditional modelling approaches [42]. There is no need for the elimination of outliers or prior data transformation. They can automatically handle interactive effects among predictors and fit complex non -linear relationships. The BDTR possesses the strengths of two algorithms: (i) boosting approaches that are adaptive methods for combining many simple models to render better and enhanced predictive 8 performance. For this study, least -squares boosting (LS -Boost [43]) is used. (ii) regression decision trees that relate a response to its predictors by recursive binary splits [42]. As explained in Figure 1 , the RSM - based dataset is merged with the I TM ones. After that, these 160 samples are used for BDTR -based metamodelling, where 90% of the data with a 0.1 learning rate was set aside for training along with using random partition data for cross -validation. The seed for the random number generator an d ‘MinLeaf’ is set to 5 (this specifies the least number of data instances that can exist within a terminal leaf node of a decision tree) using the ‘ RegressionTree.template’ available in Matlab. 4. Results and discussion 4.1 SNR -based sensitivity analysis The proposed strategy for SA is implemented for four thickness values 0.15, 0.2, 0.25, and 0.3 -mm with a related deterministic burst pressure of 9, 12.5, 14.7, and 17.3 MPa, respectively. For all cases, the coefficient of variation (COV) is 5%. The COV an d mean of winding angles are 2% and ±55°, respectively. Figure 2 and Table 5 show the obtained results , in which 𝑆𝑁𝑅 𝑡 and 𝑆𝑁𝑅 𝑤𝑎 are the SNR related to thickness 𝑡 and winding angle 𝑤𝑎, respectively. Using Table 5 and Equation 2, 𝐷𝐶𝐼 indices are calculated for all thicknesses ( Figure 2 ). The main advantage of this approach is its simplicity and the consideration of input factors simulta neously, unlike other SA methods, such as the straightforward one -factor -at-a-time method (OFAT) [33], which treats input variables individually. Here, all input variables change in their statistical interval and affect the response synergisticall y. Along with the Taguchi method, which reduces the number of required experiments/analyses, these features render simplicity and low computational cost to the proposed SA strategy, in contrast to other SA methods, such as the Sobol approach [21, 22, 33]. The proposed Taguchi -based strategy can be classified as a global sensitivity method as it assesses the sensitivity of the distribution of the entire parameters . In addition, it determines how input parameters influence the output when all inputs are varied simultaneously. It provides the relative proportion of each input in scattering the response. From a design point of view, composite structures are influenced by many factors; hence, developing a SA able to carry out a comparative study among effective factors to detect the most influential ones is of great importance. This approach can also identify insignificant inputs, saving time and costs for UQ analysis. Therefore, using this approach can be useful for screening uncertainties in high -dimensional problems. 9 Table 5. Response table for SNRs. Thickness (mm) SNR Level 1 Level 2 Delta Rank 0.15 𝑆𝑁𝑅 𝑡 −4.082 −8.136 4.053 1 𝑆𝑁𝑅 𝑤𝑎 −6.444 −6.955 0.510 2 0.2 𝑆𝑁𝑅 𝑡 −0.82785 −4.31364 3.48578 2 𝑆𝑁𝑅 𝑤𝑎 0 −4.72756 4.72756 1 0.25 𝑆𝑁𝑅 𝑡 −4.802 −7.088 3.006 2 𝑆𝑁𝑅 𝑤𝑎 −3.807 −7.226 3.420 1 0.3 𝑆𝑁𝑅 𝑡 −4.082 −7.543 3.461 1.5 𝑆𝑁𝑅 𝑤𝑎 −4.082 −7.543 3.461 1.5 Figure 2 Obtained 𝐷𝐶𝐼 by applying a variation to the mean of ply thickness. Figure 2 shows the DCI results for progressive damage analysis of the tubes. The DCIs for each thickness value, which shows the percentage of contribution to scattering the bu rst pressure concerning input uncertainties, are different. For 0.15 mm, the DCIs are 88.82% and 11.18% for ply thickness and winding angle, respectively. This indicates that the contribution of thickness uncertainties on the 10 scattering of burst pressure i s more relevant for thinner tubes. This means that relying on a deterministic estimation of burst pressure might lead to premature and unexpected failure. Therefore, using probabilistic and reliability -based analyses and designs is essential to ensure safe ty. Figure 2 shows that upon increasing ply thickness from 0.15 to 0.2 mm, the winding angle has a higher contribution. By increasing thickness from 0.2 to 0.25 mm, this value decreases to 53.22%. However, the winding angle is still the dominant factor. Then, the DCI indices become the same for both parameters when the thickness increases to 0.3 mm. Moreover, this shows that ply thickness and winding angle are relevant factors for filament wound tubes. 4.2 BDTR -based Monte Carlo reliability analysis The explained procedure for BDTR (Figure 1 ) was implemented on the dataset , in which the results are presented in Figures 3-5. Figure 3 presents the actual against predicted values of training data points . Despite the complexity of the model, the built BDTR can predict the output very well. Figure 4 shows the mean -squared error (MSE) against the number of trees for both training and test set loss . Both training and test set loss have decreasing trend s and nearly reach zero. For more assurance, the built BDTR is tested on 20 unseen data. The r esults indicate a strong predictive performance of BDTR in metamodelling. Th erefore, this high accuracy model is further used for probabilistic analysis of the composites herein under investigation . Figure 5 shows inputs (predictors) sorted by relative importance. It demonstrates that winding angle inputs a ffect the burst pressure more than thickness. This is in agreement with Figure 2 results (t = 0.25 mm). Two different and distinct studies have similar results , which valid ate the proposed SNR -based sensitivity approach . A comparative discussion on the influence of ply thickness uncertainty on the reliability of the tubes is developed by using the BDTR -based metamodel. Ply thickness uncertainty is considered in three different scenarios: (i) deterministic ply thickness val ue (0.25 mm) disregarding its uncertainty ( S1): (ii) the wall thickness of the cylinder is divided into four thicknesses with equal mean and COV (S2); and (iii) due to production inconsistencies in the FW process, excess resin moves from the innermost laye r to the surface of filament wound composite cylindrical products. Hence, ply thicknesses are sorted in descending order and assigned to plies from the innermost ply to the outermost one. Each ply has a different mean and COV. Considering each ply thicknes s as an independe nt uncertainty and formulating them separately into the problem is categorised into the mesoscale level [44] (S3). After that, using the randtool (interactive random number generation tool) available in Matlab and data from Table 3 , 105 samples are generated. Then, they are simulated by using the developed BDTR -based metamodel. Figure 11 6 shows the ob tained res ults, namely the probability density function (PDF) for these cases and their cumulative density function (CDF). Figure 3. BDTR results against training procedure. Figure 4. Convergence for the BDTR: t est and training error s. 12 Figure 6 shows a stochastic variation of the burst pressure for the tubes considering the three investigated scenarios. The deterministic burst pressure is 14.7 MPa, which was found using progressive damage analysis. Figure 6(a ) shows that despite ignoring thickness uncertainties, the burst pressure scatters from 12.5 to 16.0 MPa. This is in agreement with Figure 5 and shows the importance and high effect of winding angle uncertainties on scattering the burst pressure of the structures. Figure 6(b ) shows the PDF considering thickness uncertainty . The assumption that all plies have the same stochastic variation of ply thickness is common in the vast majority of uncertainty analyses of composite tubular structures [44]. The most interesting aspect of these results is that the burst pressure follows the Weibull distribution (see Equation 3) . A closer analysis of Figure 6(a) -(b) shows that by adding thickness uncertainties, the scattering range of burst pressure increases from 11 to 17 MPa. This means that an increase in the number of uncertainties can result in severe scattering of the burst pressure. Consequently, it requires careful monitoring of the FW process to reduce sources of uncertainties. This careful monitoring leads to assuring that the produced structures reach their maximum load -bearing capacity and be re liable to avoid unexpected failure. Figure 5. Predictors sorted on relative importance. For winding angle: ply 1 = sub -ply 1 at +55° and sub-ply 2 at −55°; ply 2 = sub -ply 3 at +55° and sub -ply 4 at −55°. 13 𝐹𝐵𝐹𝑃=20.44 14.71(𝐹𝐵𝐹𝑃 14.71)20.44−1 𝑒−(𝐹𝐵𝐹𝑃 14.71)20.44 (3) where 𝐹𝐵𝐹𝑃 is the burst pressure. Figure 6(d ) shows that the cumulative probability for all three scenarios at a pressure level of 14.7 MPa is 0.60, 0.63 and 0.75, respectively. This indicates the likelih ood of taking place burst failure for 60%, 63% and 75% out of these 105 samples, which are lower than the deterministic burst pressure. Hence, relying on deterministic analysis and ignoring uncertainties leads to unreliable results. In Figure 6 (c), one can see that the distribution no longer follows Weibull, but instead a non -normal one. In addition, considering thickness results in more scattered burst pressure. For instance, for a pressure of 13 MPa (from Figure 6 (d)), the cumulative distribution function (CDF) increases from S1 to S3. In other words, the scattering of burst pressure increases at pressures lower than the deterministic one, which means that the uncertainties play a key role in the burst pressu re of the tubes. Considering that such manufacturing uncertainties are unavoidable, a reliability analysis is strongly recommended. The BDTR metamodel was used for probabilistic analysis and then a limit state function ( 𝐿𝑆𝐹) is defined (Equation 4) for the reliability analysis. The occurrence of LSF ≤ 0 indicates that the estimation of the BDTR exceeds the threshold value ( 𝑃𝐵𝑢𝑟𝑠𝑡𝑇). Hence, the burst pressure 𝑃𝑓 can be defined as per Equation 5. 𝐿𝑆𝐹 =𝑃𝐵𝑢𝑟𝑠𝑡𝐵𝐷𝑇𝑅−𝑃𝐵𝑢𝑟𝑠𝑡𝑇 (4) 𝑃𝑓=𝑃(𝐿𝑆𝐹 ≤0) (5) By adopting the Monte Carlo sampling method for the six input variables herein considered based on the probability distribution presented in Table 3 , 105 random samples were generated. The corresponding burst pressure for each realisation of random variables was computed through the BDTR metamodel (𝑃𝐵𝑢𝑟𝑠𝑡𝐵𝐷𝑇𝑅). Then, the 𝑃𝑓 was calculated through Equations 4 -5. Monte Carlo probability convergen ce analysis was then carried out and the results show that a probability convergence happens with more than 20,000 random samples ( Figure 7 ). Therefore, due to the high capability of the metamodel to simulate large samples and ensur e the accuracy of the Monte Carlo simulation, 105 random samples were used for reliability analysis. Reliability analysis was accomplished on the three considered scenarios (S1, S2 and S3), and Figure 8 shows the results of this comparative study. Figure 8 shows that ignoring the ply thickness uncertainty for reliability analysis of the tubes can lead to underestimating or overestimating failure probability. The comparative analys is of ply thickness at different scales revealed that considering ply thickness leads to a significant underestimation of failure probability for filament -wound tubes. 14 Figure 6. Probability density function s against burst pressure for (a) S1, (b) S2 , (c) S3 , and (d) cumulative distribution function for all cases. 4.3 Reliability analysis using Monte Carlo : a modified approximation method There are some simulations that only a few runs can be affordable because their high -fidelity models can be time -consuming. Besides, constructing metamodels of complex problems often requires a large number of samples to obtain reasonable accuracy. Having a large number of samples, nonetheless , does not always guarantee the desired level of accuracy. Therefore, another strategy is herein proposed as an alternative to metamodel -based approaches. Shadabfar and Wang [11] proposed a simplified Monte Carlo approximation method for straightforward estimation of 𝑃𝑓. However, their efficient method suffers a drawback. The number of random samples required to implement is variable and different from one problem to another one. Their appro ach is modified here. In contrast to their blind sampling, this paper uses the qua si inverse transform sampling method (ITM), which covers all space, and selected 15 random samples can be representatives of the whole interval. The modified approach for approximation of failure probability is as follow s: Using the Bergman ranking method ( see Equation 6), a dummy probability ( 𝐷𝑃) is calculated, where 𝑛 is the total number of samples (40 for this study) and 𝑖 is the counter. A thorough assessment o f various problems [4, 11, 41] with differ ent inputs/dimensions and diverse degree s of complexity showed that 40 samples are enough to ensure that samples are representative of the whole design space to estimate the 𝑃𝑓 with sufficient accuracy. Using this 𝐷𝑃 as a probability, mean and 𝐶𝑂𝑉 available in Table 3 and statistical data in section 3 , the inverse of the normal cumulative distribution ( 𝐼𝐶𝐷𝐹 ) is used for the 6 random variables herein taken into consideration . This paper calls these generating strategy as the quasi -ITM sampling method. After that, these generating samples juxtapose in a randomised manner. Then, progressive damage modelling is carried out for these points to estimate the burst failure pressure. The 𝐿𝑆𝐹 parameter (Equation 4) for generat ing random samples ( 𝑋𝑖) is calculated and then sorted from small to large. For each 𝑋𝑖, Gumbel distribution (Equation 7) is used to calculate the probability ( 𝑝𝑖). Then, for each 𝑝𝑖, the corresponding inverse of the standard normal cumulative distribution function (𝑍𝑖=Ø−1(𝑝𝑖)) is calculated. The fitted curve is used to estimate the failure probability (Figure 9). The failure probability is calculated using the standard normal distribution function (Ø) at this intercept ion (Equation 8 , in which the mean is zero and the standard deviation is 1 ). The mechanism behind this approach is that when the Monte Carlo method counts the number of less -than-zero samples, it uses the trend that data displays as it approaches the negative border. The intersection shows the initial failure level [11]. DP=𝑖−0.5 𝑛 (6) pi=i n+1 (7) Ø(𝑥)=1 σ√2π𝑒−1 2(𝑥−𝜇 𝜎)2 (8) To test this modified Monte Carlo approximation reliability method, it is implemented on the deterministic burst pressure (14.7 MPa). In addition, it is used for case S3 with 40 samples, thereafter a comparative study is carried out with the BDTR -based Monte Carlo simulation method with 100,000 samples. The calculated failure probability at 14.7 MPa is 74% and 75% for the modified approximation (Figure 9 ) and BDTR methods, respectively. Figure 10 shows the implementation of these two methods for all pressures. The results show that the modified approximation Monte Carlo method with only 40 samples can estimate the failure probability with sufficient accuracy against the need for 100,000 samples from BDT R one. Figure 10 also reveals that these two obtained curves are very similar even though the strategies are different. This indicates the high capability of the modified approximation method in 16 estimating 𝑃𝑓 with significant less computational efforts. Therefore, this simplified method can be used for complex and computationally expensive problems or where building high -accuracy metamodels are unfeasible. Figure 7. Convergence of the Monte Carlo simulation at 13 MPa . Figure 8. The failure probability for S1, S2 and S3 cases using 105 samples. 17 Figure 9. Standard normal cumulative distribution against LSF at a pressure of 14.7 MPa . Figure 10. Failure probability for both BDTR /Monte Carlo reliability analysis and the proposed modified approximation approach . 18 5. Conclusion A progressive damage model was developed to estimate the deterministic burst pressure of filament wound tubes . The n, probabilistic and reliability analyses were carried out to investigate the influence of manufacturing uncertainties on the burst failure of composite tubes . A novel signal -to-ratio strategy was proposed using the Taguchi approach for sensitivity analysis and screening of the uncertainties . A novel index called Delta contribution was introduced to estimate the contribution of each input parameter while they interact simultaneously . The most relevant uncertainties in the filament winding process were considered: ply thickness and winding angle. A working guide was proposed for machine learning -based Monte Carlo simulation reliability analysis, in which a workflow was proposed for constructing high - accuracy metamodels. A hybrid utilisation of the response surface method and inverse transform sampl ing strategy was used to cover the entire design space. Then, to build the metamodel, the boosted decision trees regression method was used. The o btained results show a powerful predictive performance of the BDTR approach . After that, th e developed BDTR -based meta model was used for probabilistic and reliability analysis of the tubes . Key results reveal that the contribution of thickness uncertainties on the burst pressure scattering is more impactful when thinner tubes are considered. However, with increasi ng the wall thickness of cylinders, both uncertainties of thickness and winding angle show equal contribution in scattering the burst pressure. Hence, probabilistic and reliability analyses were carried out considering these manufacturing - induced uncertain ties. The findings of a comparative study indicate that ignoring the ply thickness uncertainty for reliability analysis of composite cylindrical shells can lead to underestimating or overestimating failure probability. At last, a computationally -cheap approach for reliability analysis was proposed. A comparative study was carried out between this method and the developed BDTR -based one. Results indicate the high capability of the modified approximation method in estimating failure probability with suffi cient and acceptable accuracy with low computational efforts . Future work will focus on a thorough experimental investigation to identify more sources of induced aleatory and epistemic uncertainties during the FW process that might affect the probability of failure of filament wound structures , such as composite pressure vessels. Data availability statement The d ataset and Matlab code s will be made available upon request. 19 Acknowledgements HA is supported by the Royal Academy of Engineering under the Research Fellowship scheme [Grant No. RF/201920/19/150].