Contour Method with Uncertainty Quantification: A Robust and Optimised Framework via Gaussian Process Regression

Over the past 20 years, the Contour Method (CM) has been extensively implemented to evaluate residual stress at the macro scale, especially in products where material processing is involved. Despite this, insufficient attention has been devoted to addressing the problems of input data filtering and residual stress uncertainties quantification. The present research aims to tackle this fundamental issue by combining Gaussian Process Regression (GPR) with the CM. Thanks to its stochastic nature, GPR associates a Gaussian distribution with every subset of data, thus holding the potential to model the inherent uncertainty of the input data set and to link it to the measurements and the surface roughness. The conventional and unrobust spline smoothing process is effectively replaced by the GPR which is capable of providing uncertainties over the fitting. Indeed, the GPR stochastically and automatically identifies the fitting parameter, thus making the experimental data post-processing practically unaffected by the user’s experience. Moreover, the final residual stress uncertainty is efficiently evaluated through an optimised Monte Carlo Finite Element simulation, by appropriately perturbing the input dataset according to the GPR predictions. The simulation is globally optimised exploiting numerical techniques, such as LU-factorisation, and developing an on-line convergence criterion. In order to show the capability of the presented approach, a Friction Stir Welded plate is considered as a case study. For this problem, it was shown how residual stress and its uncertainty can be accurately evaluated in approximately 15 minutes using a low-budget personal computer. The method developed herein overcomes the key limitation of the standard spline smoothing approach and this provides a robust and optimised computational framework for routinely evaluating the residual stress and its associated uncertainty. The implications are very significant as the evaluation accuracy of the CM is now taken to a higher level.


Introduction
The Contour Method (CM) is a destructive evaluation technique for assessing residual stress (RS) at the macro scale. Essentially, the CM entails sectioning a specimen that is meant to contain residual stress. As a result of the sectioning operation, two new surfaces are created, which are naturally subjected to an elastic stress relaxation phenomenon. This relaxation triggers a deformation over the cut surfaces, and the correspondent out-of-plane displacement is experimentally measured. Eventually, these measurements are used as an input to numerical models to back-calculate the RS over the cut surface, which was present prior to the sectioning operation [1].
Due to the relatively young age of the CM, however, little research has been devoted to the problem of uncertainties arising throughout the RS evaluation process. Correct evaluation of errors and uncertainties is of fundamental importance to both appreciate the goodness of the evaluated RS and to achieve a superior degree of reliability of structures then employed during the design process. Modern advanced structural design approaches require this additional information to provide a probabilistic assessment of structural integrity. These approaches have turned out to be of vital importance to design lightweight and high performance components, by providing an appraisal of the survival likelihood under operating conditions, instead of using overabundant safety factors [15].
At present, the problem of uncertainty related to the numerical manipulation of the experimentally measured displacement data has been only preliminary addressed by Olson et al. [16]. The authors of [16] accounted for two crucial sources of error: (i) the uncertainty in the experimentally acquired displacement data; (ii) the uncertainty arising as a result of the displacement data interpolation and smoothing. These uncertainties propagate throughout the operational steps of the CM and eventually influence the results in terms of RS. Specifically, the source of error (i) can be separated as the sum of two contributions: the material surface roughness and the intrinsic measurement error of the used measuring machine (usually a Coordinate Measuring Machine (CMM)), and thus being of stochastic nature. As such, this error can be quantified using a Monte Carlo approach, for a certain number of trials (numerical calculations), by adding a normally distributed random noise to the measured data. The error induced by (i) is then given by the standard deviation of RS from the Monte Carlo simulations. Regarding the source of error labelled as (ii), the state-of-the-art of the CM recommends to use tensor product surfaces, e.g. bivariate splines, to fit the out-of-plane displacement data for smoothing and interpolation purposes [16][17][18]. Some identification strategies of the bivariate splines are proposed by the literature, based upon a trade-off between appropriately capturing the gradients of the measured out-of-plane deformation and filtering out the noise already expected in the experimental measurement. It is evident how the process is massively influenced by the operator's experience when selecting the fitting parameters involved, i.e. the degree of the tensor product surface and the number of nodes of the interpolation grid.
For this reason, a robust evaluation of the error due to the interpolation process becomes unrealistic in the present state-of-the-art. Assuming that the interpolation process was sufficiently accurate, possibly owing to the user's experience in the use of the CM, Olson et al. proposed to assess the model error by retrieving the standard deviation of the RS computed by the CM for different values of the spline fitting parameters [16]. Finally, the total error due to the contribution of errors (i) & (ii) was essentially computed by taking their quadrature. Indeed, the method was applied on a set of five Al alloy samples. More recently, this approach has been successfully applied also to an aluminium T-section, a stainless steel dissimilar plate, a titanium electron beam welded plate, stainless steel and a nickel based alloy forged specimen [19].
As outlined before, the Olson et al. state-of-the-art approach to uncertainty quantification may be effective in some cases. Nevertheless, it would seem to lack robustness, thus highly affected by user's experience and confidence. Earlier research proposed a more sophisticated identification strategy of the fitting parameters [20]. Still, it appeared deterministic and influenced by user's experience and ability.
An advanced and powerful solution to tackle the central issue of the uncertainty quantification involved in the CM is through Gaussian Processes (GPs) [21]. GPs have recently gained popularity in the machine learning community in regression and classification tasks. Specifically, GPR is a probabilistic supervised learning technique that calculates uncertainties over predictions by assuming that the training data are jointly Gaussian [22]. This GPR can also account for the inherent uncertainty of the input dataset. Only in recent years, GPR has been applied to tackle engineering problems, demonstrating its effectiveness to model random phenomena pertinent to this field. For instance, the engineering research community applied GPR to predict the water inflow into tunnels [23], the maximum vertical displacement of a bridge subjected to uncertain load conditions [24], the mechanical response of marine structures [25], pile load bearing capacity [26], blast-induced ground vibration [27], and in control theory [28]. The application of GPR in solid mechanics is still extremely scarce, but its potential is incredibly attractive, particularly if applied to uncertainty determination problems or material response behaviour.
In this manuscript, a GPR-based approach is developed to address the outstanding issue related to the uncertainty quantification of RS when employing the CM. Dedicated optimisation tools are used to minimise the influence of user's confidence on the numerical processing of the input displacement data [29], and to devise a robust and accurate calculation framework. Firstly, GPR is applied to fit the outof-plane displacement of each half-plate obtained from the CMM measurement after the material cut, and to quantify the uncertainty of the fitting predictions. Following, each half-plate is independently analysed, thus deliberately differing from the traditional CM procedure; effects of this choice will be discussed. In this respect, each fitted displacement is used as a boundary condition in a separate numerical simulation to carry out the RS evaluation of the correspondent half-plate. In particular, a Monte Carlo approach is pursued to numerically perturb the boundary condition (prescribed displacement) according to the GPR-estimated uncertainty. This allowed for the appraisal of the uncertainty associated with the local perturbation. Eventually, the results of the Monte Carlo trials are statistically post-processed to obtain both the expected value of RS and its associated uncertainty. Alongside, an optimised open-source computational framework is set up to reduce the computation effort required to run the Monte Carlo simulations. In order to demonstrate the effectiveness of the presented methodology, a friction stir welded (FSW) aluminium plate is considered as a case study. Eventually, the performance of the proposed method is critically discussed, and the major improvements with respect to the previously used evaluation framework are highlighted.

Practical Implementation of the Contour Method
One of the prerequisites for a correct evaluation of RS is the limited invasiveness of the cutting process, which may itself introduce additional RS within the probed sample. To help achieve this ideal condition, wire electrical discharge machines (WEDM) are strongly recommended as cutting tools [30], due to the shallow affected layer generated [31]. As far as the surface displacement measurement is concerned, high-resolution CMM with contact probing system, have been identified as the preferred equipment to measure the out-of-plane displacement relaxation. The experimental data obtained from CMM is then processed to serve as an input to the FEM simulation that follows. Differently from other works seen in the literature, in the present work the data analysis of experimental measurements was performed using GPR (see "Gaussian Process Regression in the Present Study").

Case Study: Friction Stir Welded Al Alloy
RS is well-known to affect any type of welded structure, mainly due to large heat gradients involved during the manufacturing process. For this reason, in the last decades a considerable amount of effort has been put by manufacturers and researchers to mitigate this issue, primarily by reducing the amount of heat involved. Solid-state welding methods serve this purpose since the joining process occurs below the melting point of the material involved [32]. One of the most consolidated solid-state welding techniques is FSW, which allows abutting edges of workpieces, usually plates, to be joined by the mechanical action of a non-consumable rotating pin [33]. Despite the lower temperatures involved with respect to more conventional fusion welding processes, FSW is still affected by RS, although less severely. For this reason, researchers have been investigating its presence by using several experimental techniques, amongst which the CM. For instance, the CM was used to determine RS in a 2 mm thick, dissimilar aluminium-copper FSW butt-weld [34]. With regard to aluminium joints, the CM was used to examine RS in a 4 mm thick AA6061T6 FSW butt-weld [35] and in 4 mm thick AA2024-T3 FSW butt-weld [36]. Additionally, the CM was exploited to analyse RS in a 25.4 mm thick dissimilar FSW butt-weld, made of AA7050-T451 and AA2024-T351 [37].
In the present case study, the FSW technique was employed to butt-weld a AA6082-T6 4 mm thick plate. The experimental results shown here were obtained and presented in a previous publication of the same authors, in which the goal was the comparison of different solid state welding techniques [38]. Figure 1 shows the characteristics of the FSW buttwelded joint studied here such as, the dimensions of the analysed plate, the position of the advancing & retreating sides, and the welding direction. The advancing and retreating sides are the locations where the rotation speed of the pin and the welding feed speed sum and subtract, respectively. The chemical composition of the parent Al alloy is summarised in Table 1.
A CDM Rovella 650 © WEDM was used to section the FSW plate along its cross section. The reference position of Fig. 1 The FSW plate analysed in the present work (dimensions in mm). According to the CM procedure the plate was cut in half in correspondence of the WEDM cut surface S, which is indicated in blue the cut surface is shown in Fig. 1. The cut was performed using a wire of 0.25 mm diameter. The core of the WEDM wire was composed of CuZn36, whereas its external surface was brass-coated. The cut provided two half-plates, namely F1 and F2 ( Fig. 1) of equal length and two correspondent cut surfaces, S 1 and S 2 . Low-roughness of the cut surfaces was ensured by selecting a WEDM cutting speed of 5 mm/min.
The obtained surfaces S 1 and S 2 , exhibited the expected out-of-plane deformation promoted by this stress relaxation, and its measurement was carried out by means of a Hexagon Global S © CMM equipped with a 1 mm diameter ruby spherical probe. Two identical raster-scan patterns were defined over both S 1 and S 2 , which consisted of a regular grid whose nodal spacing was 0.75 × 0.25 mm. Following the defined raster-scan pattern, the CMM probed the cut surfaces S 1 and S 2 giving the displacement maps shown in Fig. 2(a) and (b). In Fig. 2, the couple (x i , y i ) refers to the coordinates of the grid node where the displacement value z 1 (x i , y i ) and z 2 (x i , y i ) were sampled. The index h assumes either 1 or 2 to refer to the half-plate F1 or F2, respectively. Therefore, z h will indicate the i-th displacement map, z 1 or z 2 . Accordingly, S h will denote one of the two cut surfaces, S 1 or S 2 .

Definition of the Finite Element Model
Aiming at performing the CM evaluation of RS in the FSW plate, a three-dimensional FE model was realised, Fig. 3. The dimensions and geometry of the model were chosen to comply with those of the obtained half-plates ( Fig. 1). The 3D domain was then discretised with 70200 tetrahedral elements opportunely distributed over the domain, using Gmsh [39]. In this instance, m indicates the mesh node lying on S h , whose coordinates are (x m , y m ) . Additional boundary conditions were prescribed on the corners of the cut surface for restraining rigid body motions but avoiding a structurally indeterminate model. Since for the present problem a purely elastic relaxation is hypothesised, the material was assumed as homogeneous and linear elastic. Therefore, E = 70000 MPa and = 0.3 were adopted respectively as Young's modulus and Poisson's ratio.
The current application of the CM differed from its standard protocol [1] as the half-plates F1 and F2 were independently analysed. If the standard CM procedure had been pursued, the displacement maps would have been aligned, interpolated over a common grid (even using GPR) and averaged. This would have provided a single displacement map that would have eventually been transferred into a FE model. Conversely, the method devised herein may overcome the possible arbitrariness of the averaging operation, provided that the presence of shear residual stress is negligible. Indeed, the method would compensate for the mismatch in terms of stress that the cut could have been induced since, in principle, each half-plate had shared the same RS state before the cut occurred.
For the sake of clarity, the following explanation of the used procedure will refer to a single half-plate. Nevertheless, in this case study the procedure is replicated identically to the counterpart.

Gaussian Process Regression in the Present Study
As far as z h (x i , y i ) is concerned, the GPR starts by defining the following dataset Specifically, z h is modelled as follows: where N(⋅ | 0, 2 n ) is the Gaussian distribution with zero mean and variance 2 n . Provided the dataset D, the GPR aims to build a probabilistic model that predicts f h at a new input m , i.e mesh nodes laying on S h . The standard assumption made by GPR is that f h is a trajectory generated by a Gaussian Process (GP) with mean M( ) and kernel K( , � ): This choice allows one to inject prior information in the learning process coming from, e.g., mechanistic models, by selecting a particular functional form for M and K (in case of no prior knowledge available it is common practice to simply select M( ) = 0 ∀ ∈ ℝ 2 [21]). Note that the above assumption is not limiting, as most of the kernels used in practice are universal approximators: the space of functions that can be represented with a GP with a universal kernel is dense with respect to the space of continuous functions. Therefore, this confers GPR the ability to approximate arbitrarily well any continuous function over a compact set [40].
Although various kernels are available in the literature [21], in this paper, because of its universality and wide use [40], the squared exponential kernel was considered: where C and l are the so-called hyperparameters. In particular, C represents how much the function f h can span. By contrast, given two generic samples z h ( i ) and z h ( j ) , l is the length scale that quantifies how much the correlation between these samples decreases as their relative distance For the sake of convenience, the hyperparameters are collected in the vector = C l ⊤ , namely the vector of the hyperparameters. (1) The computational implementation of GPR was performed using the Python module scikit-learn [29]. Upon incorporating the additive noise n = 1 µm for the input dataset, the optimised hyperparameters for f 1 and f 2 were identified according to the computational algorithm Algorithm 1. Alongside, the hyperparameters for both f 1 and f 2 are given in Table 2.
It should be noted that n quantifies the noise embedded in the experimental measurements. This assumption aimed to estimate the uncertainty due to the CMM resolution and the surface roughness of the cut surface. A potential strategy to estimate this value more accurately would be the roughness measurements of a WEDM cut on a stress-free component, which is directly linked with the WEDM process parameters. Conversely, the CMM uncertainty is provided by the manufacturer. Through the newly proposed approach the user has the capability to promptly set different values of n when this parameter in not known to evaluate its impact on the evaluated RS field. Nevertheless, this is out of the scope of the present paper.
As common in the literature [21], the identification of the hyperparameter was carried out automatically by maximising the logarithm of the likelihood p( h | , ): where h is the vector containing all the output data z h ( i ) , is the identity matrix, and is the covariance matrix defined as: Once the hyperparameters were computed, the prior mean and the kernel function were determined. Therefore, it was possible to make predictions at a new input m . In GPR this Exemplification of the FE model used for both the half-plates F1 and F2. The discretisation was performed using tetrahedral elements. The interpolated displacement (not to scale) is indicated in purple and it was applied on the cut surface S as a boundary condition. The additional boundary conditions prescribed on the corners of the cut surface are used to suppress rigid body motions. Note that the mesh is not representative of that employed in the present analysis is obtained by computing the conditional distribution of the prior GP with the observed data. In particular, in the case of GPR it is possible to use the closure property of Gaussian random variables with respect to conditional distributions and obtain that the posterior distribution of f h ( m ) is still Gaussian and with mean and variance respectively given by: Note that the computation of equations (6) and (7) only requires matrix operations and the overall complexity is cubic in the size of the dataset (due to the need to compute the inverse of a matrix of size × ). It is interesting to note that equation (7) models the uncertainty regarding the prediction at m and allows one to readily obtain confidence intervals on the predictions of the model.

Computation of Residual Stress
Exploiting the hyperparameters in Table 2, the displacement map was forecasted at every mesh node m belonging to the cut surface of the related half-plate. The predicted value of displacement at where is the model stiffness matrix, is the unknown nodal displacement vector and is the known force vector. The expansion of equation (8) leads to: where the i-th row of the system is associated to a specific degree of freedom (DOF) of each node of the model, in total N. Note that the following explanation still holds regardless of the order of the rows in equation (9). For the sake of clarity, the DOFs of the problem are sorted as shown in equation (9), i.e. the non-zero b i ∀ i = 1, 2, … , s at those DOFs where a displacement boundary condition was imposed. Conversely, the null entries b i ∀ i = r, … , N correspond to the unconstrained DOFs. After the application of the boundary conditions, i.e. the displacement boundary condition and the additional boundary conditions to prevent rigid body motions (Fig. 3), equation (9) is conveniently transformed as: The N × N matrix on the left-hand side of equation (10) is the modified stiffness matrix, and it is named as ′ . The righthand side of this equation consists of the modified load vector, which is labelled as ′ . Herein, load is intended in its broad sense as it embodies the knowledge of the imposed nodal displacements [u m ] . For the particular problem presented here, the first three rows of equation (10) are associated to the constrained DOFs by the nodal boundary conditions employed to cancel rigid body motions, i.e. u i = 0 ∀ i = 1, 2, 3 . The rows from 4 to s are related to the z-DOFs of the mesh nodes m belonging to the cut surface S h . Therefore, the corresponding equations give u i = [u i ] ∀ i = 4, … , s.
The numerical rearrangement just outlined enabled all the known quantities involved in the model to be gathered in a single vector, i.e. that on the right-hand side of equation (10). Such a unique characteristic of the system in equation (10) permitted a random disturbance vector to be introduced on the right-hand side of equation (10), thus perturbing the imposed boundary conditions: The non-zero entries of the disturbance vector are given by a Gaussian random noise having variance equal to [u m ] , over the prediction [u m ] . By defining p as the disturbance vector, equation (11) is contracted as follows: The system in equation (12) was repeatedly solved for both the half-plates F1 and F2 as in a standard Monte Carlo simulation. Aiming at optimising the overall computational process of the simulation, the LU-factorisation of ′ was computed only once and then employed for the computation of the nodal displacements in each simulation of the Monte Carlo approach. Subsequently, through equation (8) the nodal reactions were computed. Lastly, the RS was evaluated. The total number of simulations was arbitrarily set to 1000. Each i-th simulation provided the full-field RS over the correspondent cut surface S 1 and S 2 , namely (i) zz,1 (x, y) and (i) zz,2 (x, y) . Following, (i) zz,1 (x, y) and (i) zz,2 (x, y) were collated in the set R (i) (x, y) (the dependence on the coordinates is omitted for each (i) zz,k to lighten the notation, whereas it is kept for R(x, y) to explicitly indicate the dependence of R (i) on the spatial coordinates): Before the WEDM cut occurred, the cut surfaces S 1 and S 2 had been the counterpart of each other and, in principle, had shared exactly the same RS state. Hereafter, the cut surface will be univocally indicated as S. Furthermore, the related RS state of S at the i-th iteration will be denoted by (i) zz (x, y) and its associated uncertainty by U (i) zz (x, y) . According to the set of results R (i) (x, y) (equation (13) Given that U (i) zz (x, y) is the standard deviation of the expected value of RS, it represents a confidence interval of approximately 68%.
The number of trials for the Monte Carlo simulation was set to 1000. However, the convergence of the simulation was periodically monitored through the following residual indicator: zz (x, y) is the expected value in equation (14) calculated at the j-th iteration, and ‖ ⋅ ‖ 1 is the 1-norm pointwisely computed with respect to each mesh node (x m , y m ) ∈ S: Therefore, equation (17) gives: In this regard, the convergence could be considered as reached when i is less than a threshold value th opportunely selected.
The developed computational framework is briefly illustrated in Algorithm 2. Additionally, Algorithm 2 allows for certain flexibility in computing RS. In this instance, if the uncertainty quantification were of marginal interest for the users, they would perform the GPR on the input displacement data, impose p = and run Algorithm 2 once. In this case, the convergence check would not be required. Alternatively, the user would perform GPR only to smooth the experimental dataset disregarding the uncertainty quantification and pursue the traditional CM procedure.
It is important to mention that the computation of the standard deviation as in equation (15) may overestimate the uncertainty compared to the standard CM procedure, in particular instances. Moreover, probed samples containing relevant amount of shear residual stress may induce an out-ofplane displacements mismatch between the two cut surfaces. In the present study, the residual shear stress is thought to play a negligible role, although not fully confirmed, so a first estimate of the sample repeatability error was sought. This uncertainty could be reduced if more samples of the same plate were available. Unfortunately, in the present work this aspect could not be verified, given the limited number of data (only two cut surfaces). Therefore, a reliable estimate of the error due to the analysis of multiple samples, i.e. the repeatability error [19,43], could not be accomplished. In this instance, the overestimation of the uncertainty provided by equation (15) could be seen as a means to compensate for: i) the averaging uncertainty that is often underestimated or neglected when pursuing the standard CM procedure; ii) the mismatch in terms of RS that may arise over the cut surface of each half-plate as a result of the WEDM cut. Nevertheless, the authors' aim is to propose an alternative approach that is sufficiently flexible to be applied under different circumstances. Figure 4 shows the GPR predictions of the measured outof-plane displacement (z 1 and z 2 in Fig. 2). In this figure, the z-axis reports the predicted value of the displacement, namely [u m ] , in correspondence of each mesh node (x m , y m ) belonging to the cut surface S h . The colour bar indicates the uncertainty of the prediction [u m ] , namely √ [u m ] . Globally, the GPR-estimated uncertainty reaches values around 1.30 µm. However, at x = 0 and x = 168 mm for both S 1 and S 2 , higher values of uncertainty can be noticed. This deviation from the general trend of √ [u m ] could be attributed to a lack of data at the surface perimeter. Thus, the GPR was able to fit the edge-data but resulting in a higher level of uncertainty, as high as 2.2 µm.

Results and Discussion
Referring to a plane section of the data set shown in Fig. 4(a), taken at y m = 2 , Fig. 5 aims to show more exactly the effectiveness of the GPR. In particular, the red dots refer to the measured experimental data of z 1 (Fig. 2(a)), while the black dots represent the predicted value [u m ] . The error bars are the correspondent standard deviations of each mesh point, i.e. ± √ [u m ] . From the comparison reported in Fig. 5 it is evident that the GPR successfully fitted the data on the mesh nodes. Besides, the GPR thoroughly reproduced the data trend and filtered out the high frequency noise that affected the experimental measurement, while at the same time accounting for such a noise through the point-by-point estimate. Furthermore, it is interesting to note that the prescription of n = 1 µm as the noise level allowed for the automatic detection of the outliers, i.e. the points resulted in being outside the error bars. These points were not meaningful for the regression task, and without GPR they should have been manually removed before the CM evaluation process.
The FE Monte Carlo simulation provided the set of results R ( ) (x, y) (equation (13)). The expected value of the RS estimator ( ) zz (x, y) and its associated uncertainty U ( ) zz (x, y) were given by equations (14) and (15), respectively. Hereafter, the superscripted ( ) will be dropped to lighten the notation. The contour of zz (x, y) and U zz (x, y) are depicted in Fig. 6(a)-(b). Although the discussion of the evaluated RS distribution was already presented in a previous manuscript [38], it is worth commenting that the FSW process gave rise to a higher tensile RS within the weld affected region.
To help visualise the stress distribution within the plate cross-section and appreciate the relevant gradients, three horizontal lines scans (L) were defined over the cut surface S (Fig. 6(a)). Specifically, L = T , L = M and L = B refer to the top, middle-thickness and bottom line, located at y = 4 mm, y = 2 mm, y = 0 mm, respectively. Accordingly, zz (L) indicates the expected value of the RS extracted along the path L. These paths allowed the "M-like" shape RS distribution to be unveiled, see Fig. 6(c). The figure illustrates the comparison between zz (L). Such a characteristic distribution of RS is aligned with earlier findings [35,36], and matched the author's previously reported results [38], despite the different interpolation methodology adopted.
As far as the propagated uncertainty is concerned, Fig. 6(b) shows that U zz (x, y) is less than 16 MPa over most of the cut surface, apart from localised areas near the perimeter of S. Most likely, these exceptions are due to edge cutting artefacts which are known to considerably affect measurements in thin parts [6]. Lastly, a considerably high uncertainty was encountered at the top-left corner of S. Within this restricted area, U zz (x, y) is about 47 MPa. This could be interpreted as being a result of a lack of data within this localised area. In particular, GPR extrapolates the displacement data at the perimeter of the cut surface, especially at corners, more than it would do over its interior, where a large amount of data is present. Consequently, GPR predicts displacements with higher uncertainty in such regions, giving rise to higher uncertainty in terms of stress as well. Despite this, the presented approach provided reliable results.
The uncertainty data were also extracted along the same line scans L and displayed in Fig. 7(a)-(c) along with the expected value of the RS. In this figure the colour-filled bands stand for the interval zz (L) ± U zz (L) , thus representing a confidence interval of about 68%.
From a structural integrity viewpoint, much attention should be devoted to the weld region. In particular, the line scan zz (M) revealed that zz (x, y) is characterised by comparable peak values at the advancing and retreating sides; around 110 MPa. The developed strategy enabled the RS evaluation, i.e. zz (x, y) , to be strengthened with a measure of uncertainty. Specifically, U zz (x, y) was estimated to be 20 MPa at both advancing and retreating sides. Thus, at such regions it is possible to assert that the RS is 110 ± 20 MPa, with a confidence level of 68%.
From a computational perspective, the proposed evaluation approach appears to outperform spline interpolation commonly adopted in the CM. In particular, the GPR allows the optimal fitting parameters to be automatically identified. Furthermore, the GPR enables the user to effectively take into account measurement uncertainties and the influence of the roughness due to the WEDM cut. These features makes the devised GPR-supported approach way faster and less unwieldy in comparison with the consolidated spline-based one. With particular regard to the CM uncertainty quantification, the GPR-supported approach presented herein seems to introduce several advantages over that discussed in [16]. For instance, the uncertainty due to the CMM measurement and interpolation were assumed as separated contributions in [16]. Therefore, the assessment of their correspondent errors required two independent sets of simulations which were eventually combined after. By contrast, the GPR combined  (Fig. 4(a)) taken at y m = 2 and experimental data extracted at y = 2 . The black dots represent the GPR prediction both these uncertainties before the FE simulation, and therefore a single set of Monte Carlo FE simulation was necessary to quantify the uncertainties, allowing for a considerable reduction of the computational cost.
Moreover, the definition of the model error in [16] implies that several combinations of fitting parameters need to be surveyed to find those more suitable for the specific case study. Specifically, the users should choose, based upon their experience, the most appropriate degree of the spline as well as the most suitable number of spline nodes (four parameters in total) to perform the interpolation of the measured data. Conversely, the GPR overcomes this main limitation by reducing the number of parameters involved in the fitting process (C and l, see Table 2) and determining them optimally and automatically via probabilistic modelling (see "Gaussian Process Regression in the Present Study"). Therefore, the probabilistic nature of the GPR makes the entire evaluation more robust against the user's experience and confidence as far as the data analysis is concerned.
An additional shortcoming, linked with the use of splineinterpolation, is the numerical instability produced when attempting to fit the data in the proximity of the contour of the surface [6]. This instability is dramatically reduced when the GPR is employed since the GPR is less prone to overfitting and numerical instabilities due to its probabilistic nature [21]. For this reason, it was possible to extract the data very close to the upper and lower edge of S (line scans in Fig. 6(a)). Although the interpolation near the surface contour is always affected by the scarce experimental data available in these regions, the GPR approach can effectively quantify its uncertainty, given that GPR provides the uncertainty of the interpolation, i.e. √ [u m ], differently from spline fitting.
The entire simulation was carried out on a PC equipped with an Intel ® Core ™ i7-7500U CPU (@ 2.70 GHz) and 8 GB RAM. The simulation lasted approximately 2 and a half hours for 1000 iterations, whereas each of them took about 8 seconds. The aforementioned LU-factorisation (Computational Procedure) of the modified FE stiffness matrix ′ (equation (12)) led to a substantial reduction of the computational time for each Monte Carlo trial. Preparatory tests showed that if the LU-factorisation had not been adopted, the solution of each trial would have lasted 50 times longer.
Although the total number of simulations was set to 1000, the convergence indicator i (equation (17)) suggests that the simulation could have been interrupted approximately after 100 trials, i.e. after about 15 minutes. In this instance, i exhibited a drop of two orders of magnitude from the first trial to the hundredth, whereas i decreased only about one order of magnitude until the thousandth trial ( Fig. 8(a)). Practically, this denotes a rapid convergence of the RS estimator (i) zz (x, y) . In order to further optimise the computational cost, a convergence control can be easily implemented on-line, i.e. during the Monte Carlo simulation. To this end, it suffices to define a convergence criterion for i , such as i ≤ th , where th is a threshold. For instance, the convergence can be rationally considered to be achieved when i has decreased by two orders of magnitude. Therefore, for the particular case of the present study, th can be reasonably fixed to 0.1 (the red horizontal line in Fig. 8(a)). Aiming to display the rapid convergence of the RS estimator of (equation (14), (i) zz (x, y) extracted at L = M , i.e. (i) zz (M) , was stored at the iterations i = 1, 10, 50, 100, 500, 1000 and eventually plotted on the same graph, (Fig. 8(b)). As it can be seen, after the hundredth simulation the fluctuations become negligible, making the RS profiles overlapping. Henceforth, the convergence of the RS estimator can be assumed as attained.

Conclusions
A comprehensive methodology to tackle the uncertainty quantification of RS when utilising the CM has been developed and illustrated in the present paper. Using GPR, it was possible to circumvent the lack of robustness implied by the well-established spline-based interpolation and smoothing strategy of the experimental data of the CMM out-of-plane displacements. Specifically, the GPR was adopted given that the optimal fitting parameters can be stochastically and automatically determined. Therefore, this approach allowed for the minimisation of user's intervention which is inevitably reflected in the lack of uniqueness of the fitting process. Besides providing the data fitting, the GPR estimated the uncertainty that arose as a consequence of this procedure. Furthermore, the resulting uncertainty also embedded the sources of uncertainties associated with the intrinsic CMM measurement error, and the cut surface roughness inherited from the WEDM cut.
Given that the GPR condensed these sources of uncertainty along with the smoothed dataset, the latter was appropriately perturbed and fed into a single Monte Carlo FE simulation to estimate the associated uncertainty as concerned the RS. An efficient strategy to compute the Monte Carlo simulation was developed by exploiting the LU-factorisation of the modified FE stiffness matrix. In particular, this matrix was computed only once and then employed for the solution of each Monte Carlo iteration.
Finally, a convergence criterion was developed to stop the simulation when a satisfactory result was achieved. If this criterion is implemented on-line, then the evaluation of uncertainties becomes rapid and accurate, sufficiently rapid to be routinely employed.
Aiming to assess the performance and the effectiveness of this method, a 4 mm thick AA6082-T6 FSW was considered as a case study. The simulation provided the RS field along with the associated uncertainty quantification, thus providing additional information and also giving and appraisal of the reliability of the results. For this particular case study the method revealed, with a confidence level of 68%, that the RS in the FSW joint reached magnitudes as high as 110 ± 20 MPa. The entire simulation can be performed, by implementing the proposed convergence criterion, within approximately 15 minutes with a standard low-budget personal computer.
Finally, it is important to note that the developed method has deviated from the standard CM procedure. In particular, the displacement input data from two surfaces from the same cut were not averaged before performing the Monte Carlo FE simulation. Therefore, two full-field RS maps were provided and eventually post-processed at the end of each trial. Additional errors could affect this strategy due to possible shear stress present over the cut surface, usually reduced by the averaging process. However, this effect was thought to be negligible in the analysed case-study. Without averaging the input data, the presented post-processing may, in turn, lead to an overestimation of the uncertainty. Nonetheless, the arbitrariness of the averaging operation envisaged by the standard protocol induces errors as well that are not accounted for.
The authors believe that the proposed strategy may represent a practical methodology that can help CM-users to routinely evaluate RS along with its associated uncertainty.
Funding Open access funding provided by Università degli Studi di Udine within the CRUI-CARE Agreement. No funding was received for conducting this study.