Efficient strategies for reliability analysis and uncertainty quantification for filament-wound cylinders under internal pressure

Abstract

Induced uncertainties during the filament winding 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 considering 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 uncertainties 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 from machine learning algorithms. 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 accuracy.

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.^1,2 Fibre-reinforced composites usually have complex microstructures, and their manufacturing processes are not trivial, making them susceptible to variations in their microstructure, geometry, and material properties. These variations can be interpreted as uncertainties, which may be the cause of variation in their mechanical properties. ³ Generally, uncertainties can be categorised into two types: (i) aleatoric, 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. ⁴ If these are not properly accounted for, unexpected failure might take place. Conventional deterministic approaches to analysing and designing 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. ⁵

It is unfeasible to fully control all parameters during the manufacturing of composites, ⁴ particularly, in the filament winding (FW) process, in which fibre volume fraction ( $V_f$ ), thickness, and winding angle might vary. ⁶ Consequently, these parameters can vary in a stochastic manner, which makes deterministic approaches non-suitable to accurately describe such structures. For instance, Rafiee et al. ⁷ considered a uniform distribution for winding angle and $V_f$ as random variables for composite cylinders to stochastically study their functional failure. Azizian and Almeida Jr ⁶ proposed a framework to capture various multiscale uncertainties for efficient and fast stochastic, probabilistic and reliability analyses of filament-wound composite tubes. They investigated uncertainties that may arise from material microstructure and complexity and inconsistencies in their manufacturing process through a screening approach. The results, however, revealed that ply thickness and winding angle are the most influential factors on the stochastic burst pressure of composite tubes. Therefore, this paper as a step forward will focus more on these manufacturing-induced uncertainties, which have significant effects on the mechanical behaviour of composite cylinders.

Progressive damage analysis of composite structures considering uncertainties is a complex and computationally expensive process.^4,6,8–11 The majority of studies on this topic^9,10,12 use Monte Carlo (MC) method for sampling purposes. However, when the failure probability is low and for complex systems with large variabilities in the governing variables,^4,6 MC requires a large number of samples to predict failure probability accurately. ¹³ 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, ¹⁴ 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. ¹⁵ Several studies^6,15–17 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 sensitivity analysis (SA), which aims at quantifying the relative importance of input parameter(s). Methods of SAs fall into two categories:^18,19 Local and global SAs. 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 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 essential to save computational efforts. Conceição António and Hoffbauer ²⁰ developed approximation models for the reliability analysis of composite laminates. Artificial neural networks (ANNs) have been used along with genetic algorithms (GAs) 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 using a global sensitivity analysis based on Sobol indices. Ellul and Camilleri ²¹ 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 concluded that more efficient and lighter vessels could be designed with their approach. Zhou et al. ²² proposed an adaptive Kriging metamodeling approach for the failure probability, and local and global sensitivity of composite radome structures. There are numerous SA approaches available and each has its benefits and drawbacks.^23,24 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 to carry out an uncertainty quantification for composite tubes considering winding angle and thickness as the sources of uncertainties. Screening and sensitivity analyses are carried out using the Taguchi approach. In addition, probabilistic, stochastic and reliability analyses are performed considering these manufacturing uncertainties that arise from the FW process.

Deterministic progressive damage analysis

The progressive damage model used here is based on the principles of continuum damage mechanics (CDM) ²⁵ using the formulation developed by Lapczyk and Hurtado. ²⁶ The damage model uses the Hashin failure criterion to set the initiation of fibre failure and the damage evolution laws can be found in Ref. 27. The deterministic material properties are shown in Tables 1 and 2. The investigated carbon fibre/epoxy tubes are subjected to uniform internal pressure with restrained-end ³⁰ boundary conditions, as shown in Figure 1. The structure is 660 mm long, with a radius of 50 mm, and each ply has a nominal thickness ( $t_n$ ) of 0.25 mm. The laminate is made of four plies with winding angles of ±55°, which is well-known as the optimum winding angle for these boundary conditions. ³¹ The simulations are carried out in Abaqus finite element (FE) package. A convergence analysis was performed, and the converged FE mesh is depicted in Figure 1, which has 2112 elements and 2261 nodes. Shell elements (type S4R) with reduced integration and hourglass control have been used in all simulations.

Table 1.

Experimentally-measured material properties utilised in the computational analyses. ²⁸

Elastic constants	Value	Strengths	Value
${Elastic\hspace{0.17em}modulus}^{Longitudinal}$	129.30 GPa	${Tensile\hspace{0.17em}strength}^{Longitudinal}$	1409.9 MPa
${Elastic\hspace{0.17em}modulus}^{\mathit{T}\mathit{r}\mathit{a}\mathit{n}\mathit{s}\mathit{v}\mathit{e}\mathit{r}\mathit{s}\mathit{e}}$	9.11 GPa	${Tensile\hspace{0.17em}strength}^{Transverse}$	42.5 MPa
${Poisson’s\hspace{0.17em}ratio}^{\mathit{P}\mathit{l}\mathit{a}\mathit{n}\mathit{e}\hspace{0.17em}1-2}$	0.32	${Compressive\hspace{0.17em}strength}^{Longitudinal}$	764.1 MPa
${Shear\hspace{0.17em}modulus}^{\mathit{i}\mathit{n}-\mathit{p}\mathit{l}\mathit{a}\mathit{n}\mathit{e}}$	5.44 GPa	${Compressive\hspace{0.17em}strength}^{\mathit{T}\mathit{r}\mathit{a}\mathit{n}\mathit{s}\mathit{v}\mathit{e}\mathit{r}\mathit{s}\mathit{e}}$	134.5 MPa
${Shear\hspace{0.17em}modulus}^{\mathit{T}\mathit{r}\mathit{a}\mathit{n}\mathit{s}\mathit{v}\mathit{e}\mathit{r}\mathit{s}\mathit{e}}$	2.10 GPa	${Shear\hspace{0.17em}strength}^{\mathit{i}\mathit{n}-\mathit{p}\mathit{l}\mathit{a}\mathit{n}\mathit{e}}$	68.9 MPa

Table 2.

Damage evolution and stabilization values used as inputs in the numerical models

${Fracture\hspace{0.17em}energy}_{Transverse}^{Tensile/compressive}$	1.6 N/mm 29
${Fracture\hspace{0.17em}energy}_{Longitudinal}^{Tensile/compressive}$	22.5 N/mm 29
${Viscosity\hspace{0.17em}coefficient}_{Transverse}^{Tensile/compressive\hspace{0.17em}direction}$	0.0005
${Viscosity\hspace{0.17em}coefficient}_{Longitudinal}^{Tensile/compressive\hspace{0.17em}direction}$	0.00005

Figure 1.

The model with ¼ of symmetry highlighting boundary conditions and the converged FE mesh.

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. ³² Uncertainties related to $V_f$ reflect on thickness variations because of variations in the volume of the matrix during the winding process. ⁶ Azizian and Almeida Jr ⁶ concluded that manufacturing uncertainties are the most influential uncertainties during the filament winding of composite 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 consisting of winding angles $+\theta$ and $-\theta$ , 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 the 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. ³² For a normal distribution, 99.73% of data lies in this interval $\left[\mu \pm 3\sigma \right]$ ³³ , in which $\sigma$ and $\mu$ 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σ] (see Figure 2(b) ). These intervals will be used in the third scenario of this paper. 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 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. In all probabilistic and reliability analyses of this paper, ply thickness and winding angles as induced uncertainties during the production of filament wound composite cylinders are inputs of study and the burst failure pressure is the output.

Figure 2.

Schematic representation of the cylinder (a) disregarding thickness uncertainty and with the same thickness for every layer (cases S1 and S2) and (b) case S3, where the thickness of every layer is uncertain and follows a normal distribution with different functions.

Table 3.

Statistical details of the input parameters considered.

Random variables	Mean	COV (%)	Probability distribution model
Winding angle $\left(+\theta /-\theta \right)$	$+55°/-55°$	2 6	Normal 6
Thickness (mm)	0.25	5 21	Normal6,34

Taguchi-Based screening and sensitivity analysis

The Taguchi approach ³⁵ is used here for screening and sensitivity analysis of the random variables. Taguchi approach is efficient at providing robust design solutions at low computational efforts via extracting more quantitative data from fewer samples. The main advantage of Taguchi approach over others, such as the one-factor-at-a-time method ³⁶ (OFAT), is that numerous factors can be simultaneously considered. Hence, it can capture the synergistic effects of parameters. A few SAs^37,38 have been carried out using Taguchi. They used per cent contribution (PC) from analysis of variance (ANOVA) and signal–to–noise ratio (SNR ³⁹ ) 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, two levels and two factors. ³⁵ The Taguchi response is presented in equation (1), where the variation of random parameters from its mean μ to μ + 3σ, the relative change to its deterministic value can be provided. This determines sources of uncertainties with a significant impact on scattering the Taguchi response and can be used for screening purposes to reduce the number of parameters in the UQ.

$Taguchi\hspace{0.17em}response=\left|P_{Burst}^{Taguchi}-P_{Burst}^{Deterministic}\right|$

(1)

Then, the progressive damage analysis is carried out according to the suggested trials by Taguchi. After that, the analysis is accomplished with “the smaller is better” or “the larger is better” ⁴⁰ options. Since the smaller scattering from the deterministic response is favourable, the “smaller is better” option is chosen here. Next, the response table for SNR can be obtained. Table 4 shows the SNR response table using a 0.25 mm ply thickness and deterministic burst pressure of 14.7 MPa for the tubes in this study.

Table 4.

Response table obtained from the SNR analysis.

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. ⁴¹ 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, ⁴¹ Rank 2 shows the second highest Delta value, and so on. Due to the insertion of these efficient criteria in SA and screening, a new index is herein introduced, which this study names as Delta contribution index ( $\mathrm{D}\mathrm{C}\mathrm{I}$ ). For a problem with $n$ factors, we have

${DCI}_i(\%)=\frac{{Delta}_i}{{\displaystyle \sum }_{i=1}^n{Delta}_i}\times 100$

(2)

The Delta contribution index thus describes the relative weight/contribution of each Delta considered. Therefore, for the above problem, $\mathrm{D}\mathrm{C}\mathrm{I}$ for thickness and winding angle is 46.78% and 53.22%, respectively, which were calculated through equation (2) using data from Table 4.

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. ¹⁵ Before using a specific surrogate modelling technique, it is necessary to check thoroughly the fitting quality and prediction capability.^42,43 To have an efficient and high-accuracy metamodel, a thorough investigation was carried out on available techniques.^6,15,42,43 Figure 3 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). ⁶ 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.^6,15,16 Since the RSM technique tends to disperse sample points around the boundaries of the design space and place a few 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 sequences, Latin Hypercube sampling (LHS), and the inverse transform sampling method (ITM). ⁴⁴ 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) ⁴⁵ was chosen for this study.

Figure 3.

Workflow for BDTR-based Monte Carlo reliability analysis.

The BDTR method has some important advantages leading to superiority over most traditional modelling approaches. ⁴⁵ There is no need for the elimination of outliers or prior data transformation. It 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 performance. For this study, least-squares boosting (LS-Boost) ⁴⁶ is used. (ii) Regression decision trees that relate a response to its predictors by recursive binary splits. ⁴⁵ As shown in Figure 3, the RSM-based dataset is merged with the ITM 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 learning rate determines the speed at that the learner converges to the optimal solution and defines the step size during the learning process. The seed for the random number generator and ‘MinLeaf’ is set to five (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.

Results and discussion

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.0, 12.5, 14.7, and 17.3 MPa, respectively (statistical input values shown in Table 3). Statistical details of the input parameters considered.). Figure 4 and Table 5 show the obtained results, in which ${SNR}_t$ and ${SNR}_{wa}$ are the SNR related to thickness $t$ and winding angle $wa$ , respectively. Using Table 5 and equation (2), $DCI$ indices are calculated for all thicknesses (Figure 4).

Figure 4.

Obtained $\mathrm{D}\mathrm{C}\mathrm{I}$ by applying a variation to the mean of ply thickness.

Table 5.

Response table for the SNR analysis.

Thickness (mm)	SNR	Level 1	Level 2	Delta	Rank
0.15	${SNR}_t$	−4.082	−8.136	4.053	1
	${SNR}_{wa}$	−6.444	−6.955	0.510	2
0.2	${SNR}_t$	−0.82785	−4.31364	3.48578	2
	${SNR}_{wa}$	0	−4.72756	4.72756	1
0.25	${SNR}_t$	−4.802	−7.088	3.006	2
	${SNR}_{wa}$	−3.807	−7.226	3.420	1
0.3	${SNR}_t$	−4.082	−7.543	3.461	1.5
	${SNR}_{wa}$	−4.082	−7.543	3.461	1.5

Figure 4 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 burst 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 scattering of burst pressure is 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 safety. Figure 4 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.2%. However, the winding angle is still the dominant factor. Then, the $\mathrm{D}\mathrm{C}\mathrm{I}$ 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 under internal pressure.

Thin-walled composite shells are very thickness-dependent given that this geometric parameter drives the determination of whether the shell is thin or thick, for instance. Regarding the winding angle, due to the high anisotropy of long carbon fibres, even low variations around the optimum winding angle of ±55° lead to drastic changes in the burst pressure of composite shells. These reasons explain why the ply thickness has very high relative importance over the winding angle for a mean ply thickness of 0.15 mm, and it has lower importance when the mean ply thicknesses increase.

Thus, for a thickness of 0.15 mm, the ply thickness plays a major role in the burst pressure of the cylinder, whilst the winding angle becomes more influential as the cylinders get thicker given that the variations in the winding angle affect the burst pressure than thickness variations.

BDTR-based Monte Carlo reliability analysis

The explained procedure for BDTR (Figure 3) was implemented on the dataset and the results are presented. Figure 5 presents the actual against predicted values of training data points. Figure 6 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 trends and nearly reach zero. Figure 7 shows inputs (predictors) sorted by relative importance. It demonstrates that winding angle inputs affect the burst pressure more than thickness. Two different and distinct studies have similar results, which validate the proposed SNR-based sensitivity approach. Figures 5 and 7 show that thinner structures in which thickness plays a role as their weak point are more sensitive to thickness uncertainties.

Figure 5.

BDTR results against training procedure.

Figure 6.

Convergence for the BDTR: test and training errors.

Figure 7.

Predictors sorted by 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°.

A comparative discussion on the influence of 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 (see Figure 2): S1: deterministic ply thickness value (0.25 mm) disregarding its uncertainty while the winding angle uncertainty is considered; S2: the wall thickness of the cylinder is divided into four thicknesses with equal mean and coefficient of variation (COV) (see Table 3); and S3: in the FW process, excess resin moves from the innermost layer 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 following the introduced intervals in Section 3. After that, using the randtool (interactive random number generation tool) available in Matlab and data from Table 3, 100,000 samples are generated. Then, they are simulated by using the developed BDTR-based metamodel. Figure 8 shows the obtained results, namely the probability density function (PDF) and cumulative density function (CDF).

Figure 8.

Probability density functions against burst pressure for (a) S1, (b) S2, (c) S3, and (d) cumulative distribution function for all cases.

Figure 8 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 8(a) shows that despite ignoring thickness uncertainties, the burst pressure scatters from 12.5 to 16.0 MPa. This agrees with Figure 7 and shows the importance and high effect of winding angle uncertainties on scattering the burst pressure of the structures. Figure 8(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. ⁴⁷ The most interesting aspect of these results is that the burst pressure follows the Weibull distribution (see equation (3)). A closer analysis of Figures 8(a) and (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 reliable to avoid unexpected failure.

$f\left(BFP\right)=\frac{20.44}{14.71}{\left(\frac{BFP}{14.71}\right)}^{20.44-1}{e}^{-{\left(\frac{BFP}{14.71}\right)}^{20.44}}$

(3)

where $BFP$ is the burst failure pressure.

Figure 8c shows 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 8(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 pressure of the tubes. Considering that such manufacturing uncertainties are unavoidable, a reliability analysis is strongly recommended. Figure 8(d) presents 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 likelihood of taking place burst failure for 60%, 63% and 75% out of these 100,000 samples, which is lower than the deterministic burst pressure. Hence, relying on deterministic analysis and ignoring uncertainties leads to unreliable results.

The BDTR metamodel was used for probabilistic analysis and then a limit state function ( $\mathrm{L}\mathrm{S}\mathrm{F}$ ) is defined (equation (4)) for the reliability analysis. The occurrence of LSF ≤0 indicates that the estimation of the BDTR exceeds the threshold value ( $P_{Burst}^T)$ . Hence, the burst pressure $P_f$ can be defined as per equation (5).

${LSF=P}_{Burst}^{BDTR}-P_{Burst}^T$

(4)

$P_f=P\left(LSF\le 0\right)$

(5)

By adopting the Monte Carlo sampling method for the six input variables herein considered based on the probability distribution presented in Table 3, ${10}^5$ random samples were generated. The corresponding burst pressure for each realisation of random variables was computed through the BDTR metamodel ( $P_{Burst}^{BDTR}$ ). Then, the $P_f$ was calculated through Equations (4) and (5). Monte Carlo probability convergence analysis was then carried out and the results show that a probability convergence happens with more than 20,000 random samples (Figure 9). Therefore, due to the high capability of the metamodel to simulate large samples and ensure the accuracy of the Monte Carlo simulation, 100,000 random samples were used for reliability analysis. Reliability analysis was accomplished on the three considered scenarios (S1, S2 and S3). Figure 10 reveals that ignoring the ply thickness uncertainty for reliability analysis of the tubes can underestimate the probability of failure. The comparative analysis of ply thickness at different scales revealed that considering ply thickness and ignoring manufacturing-induced uncertainties lead to a significant underestimation of failure probability for filament-wound tubes.

Figure 9.

Convergence of the Monte Carlo simulation at 13 MPa.

Figure 10.

The failure probability for S1, S2 and S3 cases using 100,000 samples.

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 many samples, nonetheless, does not always guarantee the desired level of accuracy. Shadabfar and Wang ¹³ proposed a simplified Monte Carlo approximation method to estimate $P_f$ . In contrast to their blind sampling, we propose here a quasi-inverse transform sampling method (ITM), which covers all design space, and selected random samples can be representative of the whole interval. The modified approach for approximation of failure probability is as follows: Using the Bergman ranking method (equation (6)), a dummy probability ( $\mathrm{D}\mathrm{P}$ ) is calculated, where $n$ is the total number of samples (40 for this study) and $i$ is the counter. A thorough assessment of various problems^6,13,44 with different inputs/dimensions and diverse degrees of complexity showed that 40 samples are enough to ensure that samples are representative of the whole design space to estimate the $P_f$ with sufficient accuracy. Using this $\mathrm{D}\mathrm{P}$ as a probability, mean and $\mathrm{C}\mathrm{O}\mathrm{V}$ available in Table 3 and statistical data in Section 3, the inverse of the normal cumulative distribution function ( $\mathrm{I}\mathrm{C}\mathrm{D}\mathrm{F}$ ) is used for the six random variables herein taken into consideration. After that, these samples juxtapose randomly. Then, the progressive damage analysis is carried out for these points to estimate the burst failure pressure. The $\mathrm{L}\mathrm{S}\mathrm{F}$ parameter (equation (4)) for generating random samples ( $X_i$ ) is calculated and then sorted from small to large. For each $X_i$ , Gumbel distribution (equation (7)) is used to calculate the probability ( $p_i$ ). Then, for each $p_i$ , the corresponding inverse of the standard normal cumulative distribution function ( $Z_i={Ø}^{-1}\left(p_i\right)$ ) is calculated. The fitted curve is used to estimate the failure probability (Figure 11). The failure probability is calculated using the standard normal distribution function ( $Ø$ ) at this interception (equation (8), in which the mean is zero and the standard deviation is 1). The mechanism behind this approach is that when the MC 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. ¹³

$\mathrm{D}\mathrm{P}=\frac{i-0.5}{n}$

(6)

${\mathrm{p}}_{\mathrm{i}}=\frac{\mathrm{i}}{\mathrm{n}+1}$

(7)

$Ø\left(x\right)=\frac{1}{\mathrm{\sigma }\surd 2\mathrm{\pi }}{e}^{-\frac{1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2}$

(8)

Figure 11.

Standard normal cumulative distribution against LSF at a pressure of 14.7 MPa.

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 11) and BDTR methods, respectively. Figure 12 presents 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 BDTR one. Figure 12 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 estimating $P_f$ with significant less computational efforts. Therefore, this simplified method can be used for complex and computationally expensive problems or when building high-accuracy metamodels is unfeasible.

Figure 12.

Failure probability for both BDTR/Monte Carlo reliability analysis and the proposed modified approximation approach.

Conclusion

A progressive damage model was developed to estimate the deterministic burst pressure of filament wound tubes. Then, 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 sampling strategy was used to cover the entire design space. Then, to build the metamodel, the boosted decision trees regression method (BDTR) was used.

The results showed a powerful predictive performance of the BDTR approach. The developed BDTR-based metamodel 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 increasing 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 uncertainties. 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 sufficient and acceptable accuracy with low computational efforts. Future work will focus on a thorough experimental investigation to identify more sources of uncertainties in the FW process and take thickness and winding angle local measurements to develop FE models considering intra-ply variations of these parameters.