Comparison of Monte Carlo and bootstrap analyses for residual life and confidence interval

In the rapidly evolving innovative and industrial world, the high-tech production of sophisticated components is critical to contemporary civilisation. A primary concern has always been about durability: how long and substantial the product will serve in an application. Components such as aircraft [1], bridges [2], pressure vessels [3], and oil and gas components [4] are products that serve in zones of extreme loading. The engineering field faces significant difficulties in preventing these problems. However, analysing the causes of the failure can help prevent future material failures.

When a material fails due to surface crack growth, propagation and fracture, the crack is linked to material disintegration [5], which begins from within and becomes visible. Imperfection [6] in a material can cause the weakening of links between molecules, resulting in disintegration and the formation of cracks inside the material when subjected to intensely repetitive stress from static or cyclic loading. The surface crack propagates as the loading continues, contributing to material failure through the fracture. Several engineering products do not last the intended duration of ideal service as predicted by the manufacturer under conditions of surface cracks. This issue threatens engineering reliability and puts people’s lives at risk. One parameter that is crucial for determining reliability is the stress intensity factor (SIF) [7].

The SIF was calculated using the virtual crack closure method (VCCM) [8] as VCCM does not require a special mesh arrangement around the crack front [9]. The SIF was selected due to its advantages of computing in three different patterns of yield stress distribution such as constant, linearly and quadratically varying [10]. When the SIF was greater than the critical SIF, catastrophic failure could occur. The SIF of three-dimensional fatigue crack propagation must be computed end to end of the crack front. A three-dimensional surface crack was evaluated throughout the fatigue crack growth process. The value of SIF was affected by random parameters in the material properties— surface crack size and location [11]. Thus, a robust technique was needed to analyse a model’s crack growth, life span and SIF.

To analyse the surface crack growth based on fracture mechanism, FEM software was used. Advanced FEM was needed to update the dynamic domain during crack propagation. Thus, researchers enhanced the FEM into X-FEM [12], boundary element method [13], extended FEM [14] and S-version FEM [15]. S-version FEM has vast applications since the method has an advantage of rebuilding the crack domain. The overall domain remains the same, but the crack domain was updated separately, thus reducing the simulation time. Even though S-version FEM was used for pressure vessels [16], internal flaws [17], the welded material [18] and many more applications, only a few researchers updated the S-version FEM with the probabilistic analysis. Akramin et al. [19] simulated surface cracks using the Monte Carlo sampling method. Husnain et al. [20] used bootstrap analysis for crack growth prediction. A comparison of Monte Carlo sampling with the bootstrap method for surface crack analysis was not performed in any other study.

Since the S-version FEM is unable to handle uncertainties in parameters, numerical analysis with advanced computational analysis is required, such as a probabilistic approach. A probabilistic approach embedded with S-version FEM was used to compute uncertainties during numerical calculation. Randomness in material characteristics, crack size, component geometry and applied load all impact the results. Even though the probabilistic approach involves various techniques, two approaches, namely, Monte Carlo and bootstrap, are discussed in this article.

The Monte Carlo technique was selected for the fundamentals and benchmark purposes as it can lead to accurate size distribution [21], and the bootstrap technique estimates results by averaging various data samples. Unlike Monte Carlo analysis, bootstrap analysis estimates confidence intervals more straightforwardly [22]. Since bootstrap analysis does not simulate the whole sample but represents the whole sample in an averaging of sample data, it will reduce the repeating iteration and simulation time.

The objectives of this study focus on the modelling of surface cracked structures. The calculation of SIF and crack growth is presented in this article. Additionally, the cracked structure lifespan was predicted under fatigue loading using Monte Carlo and bootstrap approaches.

The methodology used in this study involves calculating the SIF, surface crack growth and fatigue life using the S-version FEM for non-random input or deterministic parameters. Meanwhile, the S-version FEM was embedded with Monte Carlo and bootstrap analyses for random input parameters.

The S-version FEM was derived from [23] the hp-FEM technique [24]. Further theoretical background of the S-version FEM is given in Ref [16]. The finite element formulation for the S-version FEM can be expressed as follows: 1 [KGGKGLKLGKLL]{uGuG}={FGFL} $$\left[ {\matrix{ {{K_{GG}}} & {{K_{GL}}} \cr {{K_{LG}}} & {{K_{LL}}} \cr } } \right]\left\{ {\matrix{ {{u^G}} \cr {{u^G}} \cr } } \right\} = \left\{ {\matrix{ {{F_G}} \cr {{F_L}} \cr } } \right\}$$

where

L and G in Eq. 2 represent local and global elements, respectively. ξ^G, η^G, ζ^G and ξ^L, η^L, ζ^L are the Gaussian points for global and local elements, respectively. Figure 1 shows the global and local elements with Gaussian points in the S-version FEM. Global and local elements have different coordinate systems. Thus, the Newton–Raphson method was implemented to calculate the Gaussian point in the global coordinate system iteratively. Once both global and local coordinates were calculated, the strain–displacement matrix could be determined. Displacement in each element would lead to the calculation of SIF at the crack front.

Fig. 1.

The SIF was calculated using the VCCM [25]. Figure 2 shows the VCCM concept. The VCCM was based on the displacement of the crack front. S1J$S_1^J$ and S2J$S_2^J$ (darker colour in Figure 2) are the finite faces at the crack front segment. The opening displacement, u_i, on the five nodes at the crack front was measured for the energy release rate calculation, G. Based on the energy release rate, SIF can be computed to calculate the crack growth rate.

Fig. 2.

Crack growth was calculated based on Paris’ law.

Eq. 3 illustrates the relationship between the equivalent SIF, ΔK_eq, and the crack growth rate da/dN. n and C are the material constants and are typically introduced as Paris n and Paris C, respectively. a is the crack depth or length, and N is the number of cycles. The equivalent SIF, ΔK_eq, is computed based on the Richard criterion [26]. The Richard criterion produces 3-dimensional crack growth direction and is proven by series of deterministic analyses using the S-version FEM [15, 16]. Thus, it will be appropriate to use the same criterion in this research study. In the future, other methods might be used to perform the analysis. The angle of the crack propagation, φ₀, was computed as follows:

where K_I, K_II and K_III are SIF for fracture modes I, II and III, respectively. Once the crack growth rate and angle were determined, the probabilistic method was embedded using the S-version FEM.

A probabilistic approach is suitable to determine the randomness of the input parameters as this frequently happens in fatigue surface crack analysis. Two probabilistic approaches are discussed here, namely, Monte Carlo and bootstrap. Monte Carlo is frequently used as a benchmark since it represents the actual occurrences in a distribution. The major drawback is time consumption. Furthermore, the Monte Carlo consumes much time to represent the whole distribution as the occurrences at the tail and head of a distribution rarely happen. Thus, bootstrap is considered an alternative to troubleshoot the disadvantage. Nonetheless, Monte Carlo still has been used as a benchmark.

The Monte Carlo method requires four significant steps. It starts with random sampling parameters, then generates data from the selected distribution, generates random parameters to produce an output and computes the probabilistic information.

The sampling of the random parameter, u_i, is equated with the cumulative distribution function (CDF), F_x, of the random variable x_i. The equation is presented as follows:

If X is lognormal distributed, N(μ_X, σ_X ), then s = (x_i – μ_X )/σ_X is a standard lognormal variate. The random parameter, u_i, can be shown as follows: 6 ui=FX(xi)=Φ(si)=Φ(xi−μX)/σX $${u_i} = {F_X}\left( {{x_i}} \right) = \Phi \left( {{s_i}} \right) = \Phi \left( {{x_i} - {\mu _X}} \right)/{\sigma _X}$$

Thus,

where Φ⁻¹ is the inverse of CDF of the lognormal distribution.

Once the random parameter is generated, the output from the S-version FEM can be produced.

If 100 random parameters are generated, 100 outputs are computed. Based on the 100 outputs, probabilistic information could be extracted. For example, the distribution of results and probability of failure could be computed from the 100 outputs. Since the Monte Carlo method takes a long time to cover the whole distribution completely, the bootstrap method can be used for comparison.

Bootstrap is used to generate sample parameters. Unlike Monte Carlo, the sample was generated based on the whole distribution. Bootstrap selects a few samples and calculates the mean. Based on the mean value, the random parameter was generated. The random parameter was used to compute the outputs. For instance, if Monte Carlo simulates a 100 samples, which means 100 times, bootstrap might reduce the simulation into 10 simulations only. A 100 samples were classified into 10 groups, each with 10 samples. The mean from every group will be calculated and fed into the S-version FEM. A comparison of the concept is shown in Figure 3. Bootstrap was used to calculate the mean parameter in a smaller number of samples and let the mean value represent the random parameter. The random parameters vary according to the material and geometry model.

Fig. 3.

A four-point bending model was constructed for the comparison between Monte Carlo, bootstrap and deterministic (non-probabilistic) approaches. The four-point bending model was 160 mm × 60 mm × 20 mm. Figure 4 shows the model with a surface crack at the centre and the surface crack plane. Material properties for four-point bending are presented in Table 1. The material was aluminium 7075-T6 with the chemical composition of Zn (5.6 wt%), Mg (2.5 wt%), Cu (1.6 wt%) and Cr (0.3 wt%). The specimen was prepared in a log traverse direction. The four-point model parameters and type of distribution were used, as described in Ref [27]. All random parameters were computed based on lognormal distribution as the lognormal distribution has a small coefficient of variation. The initial crack length and depth were replicated in the experimental results.

Fig. 4.

The four-point bending specimen was tested on a servo-hydraulic fatigue testing machine, as shown in Figure 5. The stress ratio was 0.1. The maximum load was 40 kN and 22 kN for mean stress. The experiment was carried out with a frequency of 20 Hertz. The load was cyclically loaded until 10,000 cycles. Then, the stress ratio was changed to 0.8 for 100,000 cycles. The stress ratio was 0.8 for benchmark development on the crack surface. The stress ratio was between 0.1 and 0.8 until the specimen was fractured. Once the specimen was fractured, the crack length and depth were measured using the Baty Vision System, as shown in Figure 6. The experimental work was repeated until the distribution of results was constructed. Based on the statistical data, bounds of results were produced to show the range of the results. Lower and upper bounds were computed.

Fig. 5.

Fig. 6.

The data for lower and upper bounds were calculated based on the mean and standard deviation from the simulation results. The minimum result denotes the lower bound, and the maximum result denotes the upper bound. Furthermore, the bounds of the probabilistic information in the S-version FEM results were constructed from the random parameters. Therefore, there were lower and upper bounds in the results, which represent the ranges of the variation of the results. In addition, fatigue life was predicted to be in between bounds. The bounds allow for a reliable and efficient solution to the problems associated with a certain percentage of confidence level.

The 95% confidence bounds are the most used, but the percentage can be set at any level between 0% and 100%. The confidence bound denotes the probability that a result will fall between a range of values. Confidence bounds measure the degree of confidence in the probabilistic method. Thus, in this research, 95% confidence bounds were used since it is most likely the actual range of the samples. It is reasonable to expect 95% confidence bounds to fall within a specific distribution range. For the fatigue crack growth rate prediction, 95% confidence bounds were used, which are used to justify the prediction with experimental results. The 95% bound was calculated using the following equation: 8 95% Confidence Bound=μX∓1.96σXN $${\rm{95\% ConfidenceBound}} = {\mu _X} \mp 1.96{{{\sigma _X}} \over {\sqrt N }}$$

where N is the total number of samples, μ_X is the sample mean and σ_X is the sample standard deviation.

In this study, an algorithm for probabilistic surface crack analysis was developed. Validating the results with other deterministic numerical solutions was essential. Therefore, the probabilistic and deterministic approaches were compared. Then, the reliability of the developed model was judged. The verification was focussed on SIF, surface crack growth and fatigue life.

A model was selected to calculate the SIF. Table 2 shows the details of the tension model. First, tension model A was chosen to demonstrate the capability of the developed probabilistic S-version FEM. Then, to validate the accuracy of SIF prediction, the developed probabilistic S-version FEM was used to generate SIF values in lognormal distributions.

The SIF assesses the stress intensity at the crack front and numerically evaluates the possibility of crack enlargement of an existing crack. The practical solutions [28] for crack propagation SIF were also plotted in the various graphs used for benchmarking to compare the closeness of the simulation results. Figure 7 depicts the threedimensional tension model [29]. The Newman– Raju solution is widely used to benchmark SIF in three-dimensional finite crack bodies. As a result, the Newman–Raju solution was used in this evaluation. In addition, many past studies have used it for validation purposes. The tension model was created with a surface crack in the rectangular model’s centre. For verification purposes, the mode I SIF at the crack front is presented for comparison purposes.

Fig. 7.

The tension model A had a crack of 4.5 mm length and 4.87 mm depth. Therefore, the crack size aspect ratio, crack shape aspect ratio and width aspect ratio were 0.8, 1.00 and 0.5, respectively, as shown in Table 2. The tension model A was simulated to demonstrate the capability of the developed probabilistic S-version FEM.

Figure 8 shows the normalised SIFs from 0 angles until 1 radian at the crack front. The result shows zero angle at the crack length axis and one radian at the crack depth axis. Figure 8 is computed from Ref [29], Monte Carlo, bootstrap and deterministic techniques. The deterministic is a nonprobabilistic approach where the standard deviation of each parameter is set as zero. For the probabilistic approach, the bounds were computed based on Eq. (8). Monte Carlo and bootstrap were simulated for 20 samples of SIFs with the lognormal distribution.

Fig. 8.

According to Ref [27], the crack shape is fraught with uncertainty, complicating the prediction of a probabilistic lifetime. The uncertainty is the cause of the crack shape aspect ratio’s randomness. Then, the bounds were calculated using the minimum and maximum values from the generated SIF samples. In addition, the SIF was considered acceptable since the values were greater at 0 than at 1 angle. Therefore, this demonstrates that numerical calculations produce good results. The behaviour of SIFs from the Monte Carlo and bootstrap approaches, which were embedded in the Sversion FEM, agrees well with the Newman–Raju solution since the trend is close to each other. Even though the Newman–Raju deviated from the maximum bound, the trend of all methods agrees.

Figures 9A–9C show the SIFs by deterministic, Monte Carlo and bootstrap approaches. Every beach mark (BM) of fatigue crack is presented by BM 1 until BM 6. All beach marks were compared with critical SIF (K_IC). The SIF value of SIFs was acceptable since its values were far below the K_IC of 29 MPa.√m. If the SIF value was more than K_IC, catastrophic failure will occur. It showed that the load applied does not produce a high SIF value until it reaches the critical limit.

Fig. 9.

Meanwhile, surface crack growth was verified using a four-point bending model under the mode I condition. The crack length measured between the beach marks by previous researchers was compared to that developed by the probabilistic S-version FEM analysis. The crack growth was predicted using the bending model, as depicted in Figure 10. Two supports and a two-point load were computed using the model. The zero-degree fatigue crack surface was used for crack growth validation. The analysis demonstrated the adaptability of the developed approach to various models. The deceleration effect on crack growth caused by bending was a unique feature of the bending model. The material properties were simulated using the specifications listed in Table 1. Aluminium alloy 7075-T6 was chosen as the material due to its significance in the engineering industry.

Fig. 10.

Figure 11 shows the boundary condition of the four-point bending model. Constraints at the bottom and load were applied at the top. A semielliptical crack shape was introduced in the middle of the model. The semi-elliptical crack was surrounded by a denser mesh. The domain of the semielliptical crack surrounded by a denser mesh was introduced as a local mesh. The initial crack length and depth were referred to in a previous experimental study [30]. The four-point bending model was set up as the actual experimental setup. This way, the experiment results were examined by using the four-point bending model, especially in the local mesh area. The pre-crack was modelled as precisely as per the specimen in the experiment. Since fatigue load was applied in the experimental setup, the minimum load was 4.5 kN, and the maximum load was 45 kN. More details of the precrack size can be referred to in Table 2. The standard deviation for critical SIF, K_IC and yield tensile strength were set to zero because both parameters were used as a comparison value. For instance, the critical SIF was used as a comparison value for fatigue crack to reach the unstable crack growth. In addition, both comparison values were not affected by the calculation of SIF, surface crack growth and fatigue life.

Fig. 11.

The deterministic results were compared with the experimental results, showing a deviation. Troubleshooting was needed to improve the analysis. A probabilistic approach was embedded in the S-version FEM to improve the outputs and demonstrate the capability of the developed software. Hence, Monte Carlo and bootstrap provided a probabilistic view of the outcomes. The probabilistic crack growth model was created using 20 samples in total. As the samples increase, the accuracy will improve.

Therefore, Figure 12 illustrates the beach marks formed at the semi-elliptical surface crack. The focus was on the crack front since the propagation started from here. As the load was repeatedly applied, the crack front slowly propagates. The phenomena keep repeating as the fatigue phenomena were replicated. It was described in six beach marks, as shown in Figure 12. Two beach marks were used for the validation process to avoid complex figures and graphs. As a result, the graph turned out to be unreadable. The deceleration effect on crack growth caused by bending was a distinct feature of the bending model. A higher SIF value at the 0 and 180 crack angles than the 90-degree angle shows the acceleration of crack propagation— higher crack propagation at the non-90°.

Fig. 12.

The initial crack length simulated by Monte Carlo and bootstrap was 14.00 mm, and the initial crack depth was 8.00 mm, as shown in Figure 13. Using Eq. (3), the fatigue power parameter, n, was set as deterministic (mean zero) to control the acceleration of fatigue crack growth. The maximum applied load was 45 kN, and the minimum applied load was 4.5 kN. The fatigue region or reliable crack growth was the only area where surface crack propagation was considered. As a result, the developed probabilistic in the S-version FEM is a reliable simulation tool for predicting crack growth under fatigue loading, as shown in Figure 13.

Fig. 13.

Figure 13 shows the comparison for two beach marks between the experimental approach by Kikuchi et al. [30], deterministic, Monte Carlo and bootstrap approaches, respectively. The lognormal distribution was used to generate the mean Monte Carlo and bootstrap because it impossibly produces negative values in a random number. Thus, it is fit for computing a real fatigue life and SIF value. The prediction of the probabilistic S-version FEM of surface crack growth remains valid, and the trend was acceptable compared to the experimental results. Figure 13 illustrates that the 95% confidence bound for the third beach mark and sixth beach mark was compared with the experimental results. For the third beach mark, results for deterministic, bootstrap and Monte Carlo approaches were very close to the experimental result. The results also lie inside the range of the bound for 95% confidence bounds. The 95% confidence bounds were calculated using Eq. (8). Meanwhile, for the sixth beach mark, a little deviation in Figure 13 was due to the randomness of occurrences computed during the simulation. Nonetheless, Monte Carlo and bootstrap show approximate solutions with the experimental results.

Figure 14 shows the comparison for six beach marks, and the red rectangle was focused on the sixth beach mark. The result for the first beach mark until fifth beach marks was closer to the experimental result. Meanwhile, for the sixth beach mark, the mean for Monte Carlo slightly deviated from the experimental result. Unlikely, for the mean bootstrap, it remained close to the experimental result, as shown in the red rectangle in Figure 14. This was because the Monte Carlo was capable of simulating the rare occurrences in a parameter distribution. Thus, the deviation might happen when a high number of samples were conducted.

Fig. 14.

Meanwhile, since only 20 samples were generated for bootstrap, it cannot replicate the experimental results closely. Therefore, more samples are suggested for the bootstrap method. The number of samples for the bootstrap method is also a good topic for a new investigation.

The developed probabilistic S-version FEM was then validated by predicting fatigue life. As a result, the experiment in Ref [31], deterministic, Monte Carlo and bootstrap approaches were compared to predict the fatigue life in Figure 15. A four-point bending model based on the crack growth prediction was used for verification. The same material was used, which was aluminium alloy (7075-T6). The maximum load applied was 45 kN, and the minimum load applied was 4.5 kN. The prediction of fatigue life was simulated using Monte Carlo and bootstrap with 20 samples. Then, the average of 20 samples was used to calculate the mean fatigue life. Monte Carlo was close to the experimental results when the comparison was made. However, deterministic and bootstrap methods showed deviation from the experimental results. The divergence in the results showed that the randomness and uncertainty in the analysis were modelled thoroughly. However, the Monte Carlo was the best method to predict fatigue life. Twenty samples for bootstrap are insufficient for generating a thorough fatigue life distribution.

Fig. 15.

Based on the research objective for this study, the probabilistic S-version FEM successfully predicted the fatigue crack growth and fatigue life. The SIF was determined and validated with the primary reference for analytical SIF calculation produced by Newman and Raju. A 10% deviation between the analytical solution and developed probabilistic S-version FEM is acceptable.

Fatigue crack growth determined using the bootstrap technique has a small advantage. The fatigue crack growth produced by bootstrap was 20% closer to the experimental results than the Monte Carlo method. However, the randomness of random parameters during the generating process by Monte Carlo led to the deviation.

Meanwhile, for fatigue life prediction, the Monte Carlo technique is the best method to predict fatigue life. The randomness of generated parameters was scattered in the whole distribution, and Monte Carlo can still produce a viable fatigue life prediction. Monte Carlo prediction was the closest to the experimental result.

Thus, the Monte Carlo and bootstrap were superior to another based on the specific aim and targets of the simulation.