Evaluation of global sensitivity analysis methods for computational structural mechanics problems

Abstract

The curse of dimensionality confounds the comprehensive evaluation of computational structural mechanics problems. Adequately capturing complex material behavior and interacting physics phenomenon in models can lead to long run times and memory requirements resulting in the need for substantial computational resources to analyze one scenario for a single set of input parameters. The computational requirements are then compounded when considering the number and range of input parameters spanning material properties, loading, boundary conditions, and model geometry that must be evaluated to characterize behavior, identify dominant parameters, perform uncertainty quantification, and optimize performance. To reduce model dimensionality, global sensitivity analysis (GSA) enables the identification of dominant input parameters for a specific structural performance output. However, many distinct types of GSA methods are available, presenting a challenge when selecting the optimal approach for a specific problem. While substantial documentation is available in the literature providing details on the methodology and derivation of GSA methods, application-based case studies focus on fields such as finance, chemistry, and environmental science. To inform the selection and implementation of a GSA method for structural mechanics problems for a nonexpert user, this article investigates five of the most widespread GSA methods with commonly used structural mechanics methods and models of varying dimensionality and complexity. It is concluded that all methods can identify the most dominant parameters, although with significantly different computational costs and quantitative capabilities. Therefore, method selection is dependent on computational resources, information required from the GSA, and available data.

Impact Statement

The lack of practical guidance for nonexperts and general users when selecting candidate global sensitivity analysis (GSA) methods promotes the use of inadequate or inefficient approaches for a given workflow, such as the implementation of local methods on nonlinear problems or a higher order effects study when only identification of dominant parameters is required. Using an inadequate GSA method can result in the introduction of error, inconclusive results, and the use of unnecessarily high or inaccurately low datasets. Our objective is to inform a nonexpert user in GSA method selection and best practices for structural mechanics based on modeling and analysis attributes. This article is written from a practical viewpoint.

1. Introduction

The rapid growth of accessible computing power has enabled the formulation and analysis of complex, physics-based computational models in the field of structural mechanics. Researchers can now realistically capture the physics, chemistry, and mechanics necessary to simulate phenomena ranging from the atomic scale (such as exploring the effects of chemical composition on mechanical properties) to the macro scale (such as investigating the progressive failure of layered or hybrid structures). Techniques and methods used to formulate these computational models vary substantially with respect to fidelity, computational demands, and complexity where the selection is based on the physical process to be simulated, data quality, and required accuracy. Simulating a specific scenario in a reasonable amount of time, even one with a high dimensionality and complex behavior, is achievable with today’s computational capabilities. However, comprehensive exploration of the design and analysis space to understand the change in structural performance caused by geometry, loading, boundary conditions, and material properties can quickly become intractable, even when utilizing high-performance computing (HPC) resources.

The curse of dimensionality must therefore be overcome by reducing the parameter space to include only the most influential parameters on the outputs under evaluation and identifying the significant interactions between parameters. Such information then informs uncertainty quantification (UQ) and improved model formulation to ensure that the most dominant parameters and behaviors are fully characterized and captured. Additionally, experimental testing and parameter characterization can then focus on limited data collection to fully quantify the influential input parameters and to better understand significant parameter behavior and interactions that impact model performance. Identification of the most significant input parameters also informs the development of fast-running, yet accurate, lower fidelity models, such as a reduced order surrogate model. Therefore, an approach to efficiently sample the input space so that the model is evaluated as few times as possible while capturing the important statistical properties of the response is essential to inform model development, data collection, and interpretation of the results.

Simply put, sensitivity analysis is a method which links the uncertainty in the output of a computational model to the input parameters. A local sensitivity analysis captures sensitivity relative to an individual parameter and global sensitivity analysis (GSA) evaluates sensitivity with regards to the entire parameter space. Thus, results from a GSA can be used to rank the influence of each parameter on a particular response. While GSA provides the necessary information to identify dominant input parameters, many distinct types of GSA methods are available, presenting a challenge when selecting the optimal approach for a specific problem. This ambiguity is exacerbated as direct comparisons between methods and discussion of method selection are lacking in the GSA literature for structural mechanics. Most of the literature pertaining to structural mechanics is comprised of papers describing the complex mathematics of new or enhanced GSA techniques using structural mechanics problems to demonstrate the method or cases of a single GSA method applied for a specific structural mechanics problem with the objective of identifying the most influential parameters. An example of a new method demonstrated using structural mechanics problems is the investigation of out of plane unidirectional composite laminate properties using Bayesian multimodel inference and importance sampling to calculate probabilistic Sobol’ indices (Zhang et al., 2021a, 2021b).

Representative examples of applied GSA in structural mechanics include the analysis of cable-stayed bridges using the method of Sobol’ (Nariman, 2017), structural vibration of a nuclear reactor using the global Fourier amplitude sensitivity method (Banyay et al., 2020), thrust force ripple analysis of linear switched reluctant motors using the method of Sobol’ (Chen et al., 2021), and a demonstrative article where finite element methods and GSA are used to evaluate models pertaining to vehicle impact, blast containment vessels, spinal impact injuries, and space shuttle flow liners (Thacker et al., 2005).

Comparative studies of GSA methods are essential to evaluate the differences in the results from the varying methods and to inform method selection for a specific problem given that the many techniques are based on distinct methodology, sampling, efficiency, and output. Existing studies are distributed across disparate fields and examples include using models for watershed and snow accumulation with 13 and 5 parameters respectively evaluated using local methods and GSA methods including Sobol’ (Tang et al., 2007), food safety with 5 parameters including the methods of regression, analysis of variance, and graphical methods (Frey et al., 2003), flood forecasting with 13 parameters evaluated using 10 available methods contained in the software package PSUADE and 8 parameters evaluated with local, global screening, and quantitative methods (Gan et al., 2014; Iooss and Lemaître, 2015), radioactive pollution with 21 parameters evaluated with 9 local and GSA methods (Hamby, 1995), water quality with 6 parameters evaluated with local methods including the Morris method (Sun et al., 2012), GSA methods applied to a 5 parameter thermal runaway model of a batch reactor (Saltelli et al., 2005), and building energy assessment with 11 parameters being evaluated with local and GSA methods including extended Fourier amplitude sensitivity test (FAST; Yang et al., 2016). Publications specifically comparing GSA methods for structural mechanics models are much scarcer. An example is a theoretical study of variance-based and derivative-based GSA applied to metal forming where two parameters, yield stress and hardening modulus, are investigated (Arnst and Ponthot, 2014).

Conclusions in the literature are generally that each of the methods is sufficient for the models being evaluated and that the quantitative methods are superior in the evaluation of the parameters. While a method may be sufficient, it may not be efficient, and the computational burden of a quantitative method may not be necessary if only parameter screening is needed. The model with the greatest dimensionality evaluated in the cited comparative studies contained 21 parameters, and the other models were relatively simple and often comprised of deterministic, equation-based models. Structural mechanics problems can grow exponentially with an increase in model complexity due to the inherent nature of defining material behavior in popular commercial simulation programs, for example, finite element (FE) methods, thus requiring the need for GSA comparisons of structural mechanics problems to include models with a large number of varying parameters.

The lack of practical guidance for nonexperts and general users when selecting candidate GSA methods promotes the use of inadequate or inefficient approaches for a given workflow, such as the implementation of local methods on nonlinear problems (Saltelli et al., 2019) or performance of a higher order effects study when only identification of dominant parameters is required. Using an inadequate GSA method for a computational structural mechanics problem leads to the introduction of error, inconclusive results, inefficiency, and the potential for obtaining unnecessarily large datasets. Critically, using a GSA method lacking the required functionality also deprives the user of potentially notable discoveries pertaining to parameter behavior and interactions. Our objective is to inform a nonexpert user in selecting a GSA method for structural mechanics based on modeling and analysis attributes. Along with encouraging best practices for GSA method selection, this article demonstrates the process of applying five popular GSA methods to three computational structural mechanics problems where the process and results are discussed from a practical viewpoint for nonexperts and general users. The GSA methods selected for evaluation are available from the same open-source package and readily implemented using typical desktop computational hardware/resources. Optimal performance is evaluated based on a variety of criteria including convergence, required resources, and relative accuracy.

2. Overview of GSA methods

GSA method selection was based on popularity, documentation, open-source availability, and computational resource requirements. A total of five GSA methods were evaluated representing a range of mathematical approaches, namely Sobol’ Indices, Morris Method, Extended Fourier Amplitude Sensitivity Test (EFAST), Random Sampling-High Dimensional Model Representation (RS-HDMR), and Derivative-Based Global Sensitivity Measure (DGSM). Each of these methods was analyzed using sampling methods considered to be optimal by the literature, and a brief mathematical/conceptual overview of each method is included. Table 1 compares various functionalities of each method along with a reference for additional information, where the ability to estimate the first (Si), total (ST), and higher order (n-order) sensitivity measures is noted. The first-order sensitivity measure is a parameter’s estimated contribution to the variance of the model output independent of possible interaction with other parameters. The total order sensitivity measure is a parameter’s estimated contribution to a model’s variance inclusive of interactions between all other model parameters. A method’s ability to estimate a parameter’s n-order sensitivity measures allows for the inspection of the parameter’s contribution to the model variance while including interactive effects from n number of parameters, where a second-order interaction would yield a parameter’s estimated contribution inclusive of interaction from a second parameter. Although the Morris method and the DGSM do not calculate first-order measures, these methods were included due to their value as grouping/screening procedures where they are commonly used as an initial step to reduce model dimensionality prior to employing more computationally expensive variance-based methods (Saltelli et al., 2006, 2008). All investigated methods are readily implemented using consumer-grade workstations (do not require HPC access) and are available in the open-source Python library, SALib (Jon Herman, 2017).

Table 1. Comparison of the GSA methods

2.1. Sobol’ indices: variance-based

This GSA method, originally formulated by Sobol’ (1990), applies analysis of variance to estimate the significance of an input parameter based on first and total order indices. A significant advantage of the Sobol’ method is that higher-order indices can be calculated to further investigate parameter interactions. This method measures the contribution of an input parameter by first determining the total variance of the model throughout the sample space. Then by fixing a single parameter and evaluating the model throughout the sample space again, a comparison can be made between the total variance of the model and the variance with the single input being fixed. The significance of the fixed parameter is then based on the difference between these two variances. The original Sobol’ method used a Monte Carlo algorithm to determine the first-order sensitivity measure estimates of an arbitrary group of parameters with respect to a function $ f(x) $ decomposed into $ {2}^k $ summands, where k is the number of parameters in the function and $ x=\left[{x}_1,{x}_2,\dots, {x}_k\right] $ . First-order indexes $ {S}_i $ quantify the significance of a single, independent parameter on the uncertainty of the model output.

Multiple modifications and contributions were made to the original algorithm by exploring various sampling methods, estimators, capabilities, and limitations (Jansen et al., 1994; Jansen, 1999; Sobol’, 2001, 2005; Saltelli et al., 2010). An improved sampling method expanded on Sobol’s quasi-random LPt-sequence (Sobol’ and Shukhman, 1995). Of these contributions, the Jansen (1999) and Saltelli et al. (2010) publications are considered particularly noteworthy, where the Saltelli et al. (2010) approach introduces an algorithm where first and total order indices can be calculated simultaneously while utilizing the Jansen (1999) estimator. The total order index, $ {S}_T $ , highlights the significance of a parameter while also considering all k-number interactions with the other parameters with respect to the model output. The distinction between a parameter’s first and total indices is significant in evaluating and identifying interactions with other parameters. If the first and total order indices are similar in value, then minimal interaction is present. If the first and total order indices are significantly different, then the parameter is interacting with other parameters. While not explicitly available in the SALib library, it is also possible for the Sobol’ method to estimate a parameter’s significance based on its $ n\mathrm{th} $ -order index, where the index is based on the parameter’s own contribution and its interaction with all the other parameters at the nth level, where $ n<k $ , with respect to the model output, while if $ n=k $ the resulting estimation is simply the total order index.

The Sobol’ method was chosen for investigation based on its widespread use, the ability to explore nth-order interactions, and its robustness. The specific implementation by Saltelli et al. (2010) with the total-order effects determined using Jansen’s (1999) reversed estimator scheme was evaluated (Puy et al., 2022). A brief mathematical overview begins with the decomposition form of the function in question $ f(x) $ as shown in equation ( $ 1 $ ), where each term is square integrable over $ \varOmega $ ,

(1) $$ f={f}_0+{\sum}_i\hskip0.45em {f}_i+{\sum}_i{\sum}_{i<j}\hskip0.45em {f}_{ij}+\dots +{f}_{1,2\dots k}. $$

Each term in equation ( $ 1 $ ) is obtained from the series of conditionals shown in equations (2a)–(2c),

(2a) $$ {f}_0=E(Y), $$

(2b) $$ {f}_i={E}_{X\sim i}\left({X}_i\right)-E(Y), $$

(2c) $$ {f}_{ij}={E}_{X\sim ij}\left({X}_i,{X}_j\right)-{f}_i-{f}_j-E(Y) $$

and so on for higher-order terms.

The partial variance is shown in equations (3a) and (3b),

(3a) $$ {V}_i={V}_{X_i}\left[{E}_{X\sim i}\left({X}_i\right)\right], $$

(3b) $$ {V}_{ij}={V}_{X_i{X}_j}\left({E}_{X\sim ij}\left({X}_i,{X}_j\right)\right)-{V}_{X_i}\left({E}_{X\sim i}\left({X}_i\right)\right)-{V}_{X_j}\left({E}_{X\sim j}\left(Y|{X}_j\right)\right) \dots $$

and so on for higher-order terms.

The output variance of the function is given in equation (4),

(4) $$ V(Y)={\sum}_i\hskip0.45em {V}_i+{\sum}_i{\sum}_{i<j}\hskip0.45em {V}_{ij}+\dots +{V}_{1,2\dots k}. $$

The first-order sensitivity index is given by equation (5), a ratio between the first-order partial variance and the output variance,

(5) $$ {S}_i=\frac{V_i}{V(Y)}. $$

The total variance is shown in equation (6), where N is the number of samples, A and B are sample matrixes of size $ Nxk $ . $ {A}_B^{(i)} $ is a matrix where the $ i\mathrm{th} $ column is taken from matrix B, and all other columns are that of matrix A,

(6) $$ {V}_T={E}_{X\sim i}\left({V}_{X_i}\left({X}_{\sim i}\right)\right)=\frac{1}{2N}{\sum}_{j=1}^N{\left(f{(A)}_j-f{\left({A}_B^{(i)}\right)}_j\right)}^2. $$

The total sensitivity index for variable $ i $ is shown in equation (7), as a ratio between the total variance and the output variance,

(7) $$ {S}_{T_i}=\frac{V_T}{V(Y)}. $$

Further reading and detailed derivations are provided by Saltelli et al. (2008, 2010).

2.2. The Morris method: derivative-based

The Morris method is a one factor at a time (OFAAT) derivative-based approach to determine parameter influence measures called elementary effects (Morris, 1991). Advancing beyond a simple OFAAT method, the Morris method uses an average of these elementary effects evaluated at multiple points to remove the localizing effects of a single sample method, where derivatives are taken at a single point in the parameter space. The method is found to be efficient when examining datasets with many parameters The Morris method exhibits a minimized computational cost when compared to variance-based methods (Saltelli et al., 2008), but lacks the quantitative ability of the aforementioned Sobol’ method and is therefore considered a qualitative method. The Morris method was chosen based on the reduced computational cost of generating qualitative total order sensitivity measures and the method is commonly used as a screening application to reduce model dimensionality prior to more detailed analysis with an advanced GSA method.

Sampling methods for the Morris method create a trajectory on which the parameters are varied given pseudo-random starting points. Improvements made on Morris’ original sampling method were formulated by Campolongo et al. (2007), where the samples/trajectories are chosen to uniformly stratify the dataset in p-levels, improving the sample space coverage and minimizing overlapping sequences. As an addition to the originally proposed method, Campolongo et al. (2007) suggested taking the absolute value of the calculated elementary effect. This method of using the absolute value is shown to improve the accuracy of the predicted elementary effect when the model exhibits nonmonotonic behavior. Equation (8)) shows the method in which the mean of the elementary effects of each variable is calculated as proposed by Campolongo et al. (2007). The mean elementary effect is often compared to the total order GSA measures estimated by variance-based methods (Campolongo et al., 2007; Saltelli et al., 2008; Sobol’ and Kucherenkob, 2009). Where $ r $ in equation (8) is the number of trajectories or sample sets, and $ \mathrm{E}{\mathrm{E}}_i^j $ is the elementary effect determined from each trajectory for each input variable $ {x}_i=\left[{x}_1,{x}_2,\dots, {x}_k\right] $ where $ k $ is the number of input parameters in $ f(x) $ ,

(8) $$ {\mu}_i^{\ast }=\frac{1}{r}{\sum}_{j=1}^r\left|{\mathrm{EE}}_i^j\right|. $$

Equation (9) shows the method of calculating the standard deviation of the mean elementary effect. This represents the magnitude of interactive effects between variable $ {x}_i $ and all other parameters $ {x}_{\sim i} $ ,

(9) $$ {\sigma}_i^2=\frac{1}{\left(r-1\right)}{\sum}_{j=1}^r{\left({\mathrm{EE}}_i^j-\mu \right)}^2. $$

The analytical derivative in equation (10) is the method used to calculate the elementary effects. Where $ \varDelta $ is a value contained in [1/ $ \left(p-1\right) $ ,…, $ 1-1 $ / $ \left(p-1\right) $ ], where p is the level of discretization of the sample set. Conceptually this is the partial derivative of the model where an input variable $ {x}_i $ is changed by $ \varDelta $ ,

(10) $$ \mathrm{E}{\mathrm{E}}_i=\frac{\left[Y\left({X}_1,{X}_2,\dots, {X}_{i-1},{X}_i+\Delta , \dots, {X}_k\right)-Y\left({X}_1,{X}_2,\dots, {X}_k\right)\right]}{\Delta }. $$

2.3. EFAST: Fourier transformation based

The FAST, originally proposed by Cukier et al. (1978)), uses Fourier transformation coefficients to determine first-order sensitivity indices. Later the FAST method was expanded by Saltelli and Bolado (1998) by implementing an improved sampling method and estimator of total order indices, referred to as the EFAST method. For the first-order indices, it is suggested that the EFAST sensitivity measure is comparable to that of Sobol’s first order, with a reduction in computational cost, while the total order may differ based on the interdependencies of parameters (Saltelli et al., 1999). The EFAST method represents a unique GSA method when compared to variance and derivative-based methods and the potential of decreased computational time compared to Sobol’s method given the more deterministic sampling scheme, while allowing for higher order terms to be included. The k-dimensional sample space is explored using design points taken over a search curve to explore the input space with various frequencies $ \left({\omega}_1,{\omega}_2,\dots, {\omega}_k\right) $ . These frequencies are selected to prohibit interference with one another, allowing for the frequencies to respond uniquely to each input. A Fourier transformation is used on the search curve and the resulting Fourier coefficients are used to calculate the variances required for the sensitivity indices. Equation (11) shows the general form of the search curve proposed by Saltelli et al. (1999), where $ {\omega}_i $ are the frequencies corresponding to each k-input, s is a scaler evaluated throughout the range $ \left(-\pi <s<\pi \right) $ , and $ \varphi $ is a phase-shift coefficient randomly chosen between $ \left(0<\varphi <2\pi \right) $ ,

(11) $$ {x}_i=\frac{1}{2}+\frac{1}{\pi}\left(\sin \left({\omega}_is+{\varphi}_i\right)\right). $$

Equation (12) shows that the function $ f(s) $ is a function of k-numbers of search curves, each evaluated through the range of s, where k is the number of inputs,

(12) $$ f(s)=f\left({x}_1(s),{x}_2(s),\dots, {x}_k(s)\right). $$

Equation (13) shows the method used to calculate the output variance of the model,

(13) $$ \hat{D}=\frac{1}{2\pi }{\int}_{-\pi}^{\pi }{f}^2(s) ds-{\left[\frac{1}{2\pi }{\int}_{-\pi}^{\pi }f(s) ds\right]}^2. $$

Equations (14) and (15) show how the Fourier coefficients are calculated to represent $ f(s) $ as a Fourier series expansion,

(14) $$ {A}_j=\frac{1}{2\pi }{\int}_{-\pi}^{\pi }f(s)\mathrm{cos}\ \mathrm{cos}\;(js)\; ds, $$

(15) $$ {B}_j=\frac{1}{2\pi }{\int}_{-\pi}^{\pi }f(s)\mathrm{sin}\ \mathrm{sin}\;(js) ds, $$

and

(16) $$ {\varLambda}_j={A}_j^2+{B}_j^2. $$

Equation (17) shows how the partial variance of each k-input is calculated using $ \varLambda $ , which is calculated using equation (16) and where p is the number of higher harmonics,

(17) $$ {D}_i=2{\sum}_{p=1}^{+\infty }{\varLambda}_{p{\omega}_i}. $$

Equation (18) shows the method of determining the total variance (Homma and Saltelli, 1996; Saltelli et al., 1999), which is the difference between the model’s output variance and the partial variance of all parameters except the ith term,

(18) $$ {D}_{T_i}=D-{D}_{\sim i}. $$

The total and first-order indices are then calculated using equations (19) and (20), respectively,

(19) $$ {S}_{T_i}=\frac{D_{T_i}}{D} $$

and

(20) $$ {S}_i=\frac{D_i}{D}. $$

2.4. Random sampling-high dimensional model representation: meta-model approach

This method utilizes a meta-model approach where the randomly sampled data is represented by a component function allowing the RS-HDMR method to be used on unstructured datasets. This is particularly useful when the data has already been obtained and the particular sampling method/sequence for the other methods was not utilized, whereas the other methods require specific sampling methods of the inputs. This method first approximates a function representing the dataset using an expansion of a bias function in the form of cubic B-splines followed by improved estimations of the component function using backfitting. Upon convergence of the component function, sensitivity measures are calculated using the equations shown below. A complete derivation and explanation of the RS-HDMR method is provided by Li et al. (2010). The RS-HDMR method was selected given the level of versatility the method offers when applied to existing datasets or datasets that would be too computationally expensive to obtain using specific sampling methods. The sensitivity indices are represented here as total $ \left({S}_{p_j}\right) $ , equation (2 1), structural first order $ \left({S}_{p_j}^a\right) $ , equation (22), and correlated $ \left({S}_{p_j}^b\right) $ , equation (23) contributions.

(21) $$ {S}_{p_j}=\frac{\sum_{s=1}^N{f}_{p_j}\left({x}_{P_j}^{(s)}\right)\left({y}^{(s)}-\underset{\_}{y}\right)}{\sum_{s=1}^N{\left({y}^{(s)}-\underset{\_}{y}\right)}^2}, $$

(22) $$ {S}_{p_j}^a=\frac{\sum_{s=1}^N{\left({f}_{p_j}\left({x}_{p_j}^{(s)}\right)\right)}^2}{\sum_{s=1}^N{\left({y}^{(s)}-\underset{\_}{y}\right)}^2}, $$

(23) $$ {S}_{p_j}^b={S}_{p_j}-{S}_{p_j}^a. $$

2.5. Derivative-based global sensitivity measure

The DGSM method was proposed by Sobol’ and Kucherenkob (2009) as an advancement to the original Morris method. Their work provides a link between the importance measure and the total sensitivity index. The DGSM method is selected here based on the prospect of improving the quantitative ability of the Morris method for total order sensitivity measures leading to improved screening results, where it is most utilized while retaining the relatively low computational cost compared to the more complex methods. Differentiating the DGSM from the improved Morris method proposed by Campolongo et al. (2007), the partial derivatives of a function’s local variance are squared, as opposed to taking the absolute value. Similar to equation (8) in the Morris method, equation (24) shows the method of calculating the importance measure, $ {\nu}_i $ for each variable. The normalization of this importance measure is shown in equation (2 5), where the link between the importance measure and Sobol’s total sensitivity index is provided in Sobol’ and Kucherenkob (2009),

(24) $$ {v}_i={\int}_{H^n}{\left(\frac{\partial f}{\partial {x}_i}\right)}^2 dx. $$

Equation (2 5) shows the link between the importance measure and the total sensitivity index for variable i, where D is the total output variance shown in equation (6),

(25) $$ {S}_i^{\mathrm{tot}}\le \frac{\nu_i}{\pi^2D}. $$

3. Technical approach

The convergence, performance, and computational resource requirements for the GSA methods were evaluated by applying each method to three case studies following the procedure in Figure 1. These case studies consisted of two widespread structural mechanics analysis methods, FE and peridynamics (PD), with varying model fidelity, dimensionality, expected parameter interactions, and complexity. Case 1 was a high-fidelity model with complex, multiple interacting damage mechanics, and sufficient dimensionality relative to computational structural mechanics problems that could be evaluated on consumer-grade hardware. Case 2 was known to have interacting parameters and case 3 was chosen to evaluate the effectiveness of GSA on PD models which included modeling parameters in addition to material property parameters. Surrogate models were formulated for each case study using datasets calculated with the high-fidelity models. Surrogate modeling allows for large sample sizes to be generated for the GSA at a comparably minimal computational cost relative to collecting the necessary data directly from the high-fidelity models.

Figure 1. The technical approach used for evaluating the GSA methods.

3.1. Data acquisition from case study models

The study datasets were obtained from validated macro and microstructural mechanics models created using two different software packages, Abaqus 2019 by Dassault Systèmes and Peridigm (Parks et al., 2012), a PD code developed at Sandia National Laboratories (Silling and Lehoucq, 2010). Each case was chosen to represent a unique analysis condition for GSA.

3.1.1. Case 1: layered structure under four point bend loading (FE FPB 41)

The most complex case study is the 41-parameter FE macroscale model analyzed using Abaqus. This model represents a layered metal/composite Four-Point Bend (FPB) specimen, fabricated from four different Eglass fabric/epoxy composite lamina cocured to a 5456 aluminum substrate (Figure 2a). The simulation exhibits complex, progressive failure with nonlinear behavior (Figure 2b). This high-fidelity model captures five different damage mechanisms that are each explicitly modeled: matrix cracking, fiber fracture, delamination within the composite, disbond at the composite/metal interface, and plastic deformation. All material properties included in the model are GSA parameters and most are described by sparse and/or disparate data. The goal of GSA was to determine the most influential parameters on damage tolerance to inform focused experimental testing to fully characterize all dominant parameters. Further details of the model development and validation are available in Arndt et al. (2022). Initial simulation indicated that plasticity in the metal significantly dominated damage tolerance, followed by shear in the composite. Therefore, it is expected that GSA will indicate that the most influential parameters represent these two damage mechanisms. This model was chosen due to its complex structural behavior and large number of parameters. This case represents the limit of computational time per single analysis and a number of parameters to feasibly generate an adequately populated dataset for surrogate modeling on a commercially available desktop computer.

Figure 2. (a) Solid model for the damage tolerance of a metal/composite cocured structure. (b) Representative example of simulated composite failure using the finite element model.

Table 2 provides the parameter number, description, mean value, and sampling range for each parameter used to generate the dataset for this case. Given the sparseness of the data, the mean values were used to validate the model and a range of $ \pm 20\% $ of the mean values, representing a z-score of 2 for an assumed 10% standard deviation representing typical material variability (Mead et al., 2012), was assigned to each parameter to create the sampling range for both surrogate model development and GSA sampling, This z-score was assumed to capture adequate parameter variation which is critical for GSA. The aluminum substrate is defined using an elastic-plastic Johnson Cook constitutive material model (X0–X4). All interfaces between lamina (interlaminar) and the composite/metal (bimaterial) are explicitly modeled with cohesive elements using a Cohesive Zone Model (CZM). Cohesive elements are defined with maximum principal stress for damage initiation, mixed-mode critical strain energy release rate for damage propagation, and element deletion to represent progressive failure (X29–X40). Each of the lamina types is defined using a Continuum Damage Model (CDM) to capture the in-plane failure (X5–X28). The total energy absorbed by both damage and plastic deformation was cumulative throughout the analysis and chosen as the output value to be evaluated with GSA due to its ease of automated extraction and correlation to damage magnitude.

Table 2. Parameter descriptions and ranges used in the FE FPB model study

3.1.2. Case 2: double cantilever beam with epoxy adhesive (FE DCB 9)

This nine parameter macroscale FE model created in Abaqus represents a double cantilever beam (DCB) specimen fabricated from 5456 aluminum adherends bonded by a thermoset epoxy (Figure 3). DCB testing is a standard method for obtaining mode I fracture properties and DCB simulation is used to explore the effects of varying material properties on disbond between the adherends. This model has the lowest dimensionality of the case studies and the fracture behavior is straightforward and well-characterized. Full model details and validation are provided in Smith (2021). This model was chosen to explore the possibility that a screening method may be adequate for lower dimensionality and less complex GSA scenarios, allowing reduction in computational resource requirements if accuracy is maintained. The parameter information is provided in Table 3. To test the capability of the GSA methods to encompass a larger parameter space relative to values, a range of $ \pm 30\% $ of the means (z-score of 3) was applied for this case. Based on engineering judgment, the mode I energy release rate and adherend stiffness are expected to dominate the parameter space.

Figure 3. FE model of a DCB test configuration.

Table 3. Parameter descriptions and ranges used in the FE DCB 9 model study

The aluminum adherends were defined using linear-elastic constituent properties as plastic deformation was negligible during experimental testing. The cohesive elements used to represent the thermoset epoxy were defined using CZM as in Case 1 and propagation of damage was defined using mode-I critical strain energy release rate. The model output was again the total damage energy absorbed.

3.1.3. PD uniaxial tensile test (PD uniaxial 18)

This 18 parameter microscale model predicts the Young’s modulus of a zirconium diboride (ZrB2) tensile specimen as a function of porosity using PD to test the GSA methods on a distinctly different methodology and size scale than the FE case studies. PD has been successfully demonstrated to predict complex microcrack behavior in composites and metals, as well as unit cells containing voids and cracks (Askari et al., 2008; Madenci et al., 2018). In PD, damage is incorporated into the force function by allowing bonds to break after a certain specified elongation allowing for the inherent inclusion of crack networks in analysis. PD is becoming more widespread for problems displaying multiple crack initiation and propagation, therefore, this analysis method was chosen as an alternative approach to FE. The model was recreated based on a problem from the literature (Guo et al., 2008) and was validated using all four of the reported specimens with varied void density (Figure 4).

Figure 4. PD model and validation results.

The ZrB2 was modeled as a unit cube to minimize computational cost with minimal error given the specimens’ linear-elastic behavior and the quasistatic nature of the loading condition, as the original specimen dimensions were unknown. The voids were introduced in the gauge region of the model specimen based on percent volume where the void’s position and radius were chosen based on a defined Gaussian distribution. The respective void parameters and distribution parameters are shown in Table 3 (X1–X11). GSA was performed to evaluate the effect of the voids and their geometric parameters on the effective stiffness of the uniaxial tensile ZrB2 specimen. The mean parameter values provided in Table 4 are the arithmetic means of the input to the four validation models and the data range was $ \pm 20\% $ of this mean value. PD modeling parameters were included in the GSA for this case. The parameters reported as integers comply with the required format of the PD code. The horizon value was based on Peridigm analysis requirements (Parks et al., 2012), and the maximum value was limited to 0.2 due to the rapidly increasing computational cost above this value. Note that some of the reported values in Table 4 lack units as the values were normalized to the scale of the model. Because the model was developed as a unit cube, normalization was needed to preserve the ratios between void parameters and the model dimensions as well as eliminate the use of exceedingly small values. For this case, a priori prediction of the dominant parameters was challenging as the influence of the modeling parameters was uncertain.

Table 4. Parameter descriptions and ranges used in the PD uniaxial 18 model study

3.2. Surrogate modeling

The computational time required for the parameter sampling needed for GSA using high-fidelity models is typically prohibitive, therefore surrogate models are developed to obtain the large data set needed for GSA. A surrogate model is an approximation of a complex system, based on a limited number of data points, that predicts the relationship between the system inputs and outputs using a mathematical formulation. With proper construction, surrogate models accurately mimic the behavior of the high-fidelity simulation at a much lower computational cost (milliseconds versus minutes/hours per simulation for the three cases).

To formulate the surrogate models, output data was generated at points identified using stratified Latin hypercube sampling (Shields et al., 2015; Shields and Zhang, 2016). This process was automated through Distribution-based Input for Computational Evaluations (DICE) developed by the Naval Surface Warfare Center Carderock Division (NSWCCD) (Nahshon and Reynolds, 2016). DICE is a standalone executable that generates an array of input vectors using a specified sampling method for the provided parameter ranges and distributions.

The surrogate models were trained on these datasets using forward-feed neural networks by means of TensorFlow V2.11 (Abadi et al., 2016). The datasets were broken into three subsets, the training dataset, the validation dataset, and the testing dataset. The training dataset was used to train the neural network (NN) and the validation dataset acted as a benchmark from which accuracy metrics were determined to evaluate the current training iteration of the surrogate model. The validation dataset is determined from random sampling during each training iteration where the data is shuffled and resampled before each iteration. The training/validation dataset were allocated such that the training and validation datasets were 80 and 20% of the training/validation (parent) dataset, respectively. The choice of the training/validation split percentages was based on a traditional ratio when using NN’s, where an 80/20 split was known to result in abundant training data while using a large enough portion of the sample space for validation to prevent overfitting.

This process was repeated until the root-mean-square error (RMSE) of the NN’s estimation of the entire validation dataset was less than 5%. RMSE was chosen as it is more sensitive to higher magnitude errors than mean absolute error. Upon the convergence of the NN validation metric the testing dataset that was withheld from training of the NN was used to predict values for each model, where the accuracy of the predictions made by the NN versus the actual were shown to agree well with the validation accuracy found during training. This agreement between the validation metric and the accuracy of the predictions made by the NN from the testing data versus the actual values refute overtraining of the NN model. Parameters of each of the NN surrogates are shown in Table 5, where the number of hidden layers and the constant number of neurons is specified. The input layer was equal to the number of model parameters. An output layer of a single neuron completed the NN’s architecture, determined by utilizing a grid search method, where the best-performing NN for each model was the resulting surrogate.

Table 5. Architecture of the Forward Feed Neural Network used for surrogate modeling, all of which have the same number of neurons for each layer where the number of hidden layers is specified

Table 6 shows the number of runs used to train the NN surrogate models, the training dataset sizes, the testing dataset sizes, and the RMSEs of the NN’s estimations resulting from the testing datasets. The size of the dataset generated from each parent model was dictated by the model complexity, where the parent dataset used to train and validate the NN increased relative to model dimensionality. The models with a higher number of parameters require considerably more parent runs to achieve an accurate surrogate model. The computational time to generate an adequate dataset is further impacted when considering that the models with a larger number of parameters are more complex for the cases investigated, leading to a much higher computational time per single analysis of the original models.

Table 6. The number of parent model runs required to train the surrogate models, the test/train split, and the RMSE error of the final surrogate model

3.3. Sampling for GSA

Once the surrogate model for each case was formulated, these surrogate models were used to collect datasets for GSA according to the sampling scheme specific to each GSA method in SALib. The sampling strategies for the Sobol’, Morris, EFAST, and DGSM methods are identified in their primary citations (Table 1). The HDMR method allows for any appropriate sampling method, and the stratified Latin hypercube method was applied for consistency. Parameters were normalized within the range of 0–1 to be evaluated by the NN, where normalized parameters increase the rate of training convergence. To determine the convergence rate, the sample size was increased between 300 and 500,000 for each GSA method when possible. Due to numerical or hardware limitations, some sample sizes were not possible when implementing certain methods. Specifically, the EFAST method has limitations at smaller sample sizes, while some sample sizes were too large for the HDMR method when the available RAM was insufficient (<128 GB).

4. Case studies results

4.1. Convergence study

The number of samples at convergence was determined for each method and summarized in Table 7. Convergence is determined as the number of samples required for each method to produce constant estimations (within 5% for two consecutive increases in sample size) for all parameters with a sensitivity measure contributing more than 10% to the model output. The rate of convergence for each method is visualized by heatmaps that report the absolute difference as a percentage between consecutive sample sizes (Figures 5 and 6). The heatmaps present convergence rates for the individual parameters estimated to contribute more than 10%, an average rate for all parameters contributing greater than 10% (GT 10%), an average rate for parameters that were identified to be nonnegligible by inspection of the indices (Top 10 or Top 3), and an average rate for all parameters (All). In general, more parameters resulted in larger datasets to achieve convergence. This finding is important as GSA can provide relevant information without convergence for some objectives. For example, if the objective is to identify the most dominant parameters, this information will be evident at a much smaller computational cost than if the specific ordering or converged influence values for all parameters are needed.

Table 7. The number of samples for convergence required by each method to estimate first and total order GSA measures for parameters contributing more than 10% to the model output

Note. + indicates that the convergence criteria were not achieved.

Figure 5. Heat maps for total order convergence.

Figure 6. Heat maps for first-order convergence.

4.2. Case 1: FE FPB 41

The converged first and total order sensitivity measures calculated by each GSA method for Case 1 are shown in Figure 7a,b, and quantitative sensitivity measures are provided in Table 8. Each of the GSA methods produced similar estimations for first-order indices when available. Total indices differed with HDMR and DGSM predicting higher measures than Sobol and EFAST for the most influential parameter (X2). Meanwhile, the Morris method differentiated the dominant parameters but generally underpredicted parameters of greater significance and overpredicted parameters with less significance as compared to the other methods.

Figure 7. (a) Values of total order SA measures estimated by each method for the FPB A-41 model using a sample size of 50,000 for the HDMR method and converged values at 500,000 for all other methods. (b) Values of first-order SA measures estimated by each method for the FPB A-41 model using a sample size of 50,000 for the HDMR method and converged values at 500,000 for all other methods.

Table 8. Top 5 influential parameters for Case 1 ranked in order from most (I-1) to less (I-5) significant by each method at converged values

Parameters X2, X0, X4, and X31 have the largest overall influence on the model’s output as indicated by both total and first-order indices. All methods ranked these as the top four influential parameters, although the order of importance determined by the Morris method deviated from the other methods for I-3 and I-4 (the third and fourth most influential parameters, respectively). The Morris method switched the importance of X4 and X31; however, there is a small difference in sensitivity measure between the two, indicating that they are similar in importance. Similarly, the I-5 (the fifth most influential parameter) ranking deviated with DGSM and EFAST predicting X33 and the other methods selecting X28. As seen in Table 8, there is a small difference between the two indices, again indicating that these two parameters are similar in importance.

The HDMR method predicted some interaction between parameters as indicated by the difference between the total and first-order measures. Meanwhile, the Sobol’ and EFAST methods calculated nearly identical first and total values, indicating minimal, if any, parameter interaction.

4.3. Case 2: FE DCB 9

The converged first and total order sensitivity indices for the FE DCB 9 model are shown in Figure 8a,b, with the quantitative values provided in Table 9. As with Case 1, all methods predicted similar results for first-order measures, when available, with differences noted in the total measures. The Morris method continued a trend of overestimating the less significant parameters indices while underestimating the significant parameter indices relative to the other methods, while in general predicting the same set of most influential parameters. Meanwhile, the DGSM method again predicted sensitivity measures 20–40% higher than the other methods. However, for Case 2, the HDMR predictions aligned with Sobol’ and EFAST sensitivity measures rather than the DGSM predictions as in Case 1.

Figure 8. (a) Converged values of total order SA measures estimated by each method for the FE DCB 9 model using a sample size of 75,000 for the HDMR method and 200,000 for all other methods. (b) Converged values of first-order SA measures estimated by each method for the FE DCB 9 model using a sample size of 75,000 for the HDMR method and 200,000 for all other methods.

Table 9. The most influential parameters for Case 2 ranked in order from most (I-1) to less significant (I-3) by each method at converged values

For Case 2, parameter X0 was determined to be the most significant parameter with little interaction with other parameters, shown by the large and similar magnitudes for both total and first-order measures. All methods identified the same first three influential parameters, X0, X1, and X8, with the DGSM switching the order of X1 and X8. Given the extremely small difference in sensitivity measures, indicating nearly equal importance between these two parameters, this difference is considered insignificant.

4.4. Case 3: PD uniaxial 18

Total and first-order sensitivity measures for Case 3 are shown in Figure 9, with quantitative values provided in Table 10. As with Case 2, it was observed that both the Morris and the DGSM varied from the other methods for the total order estimations while all other predictions were nearly identical. Parameters X0, X15, and X14 were the most dominant with minimal parameter interaction indicated by the relatively large and similar estimations for both total and first-order measures. The DGSM method alone ranked X15 over X0 as the most influential parameter. All methods except for the DGSM predicted X13 as the fifth most influential parameter, which instead predicted X16.

Figure 9. (a) Converged values of total order SA measures estimated by each method for the PD uniaxial18 model using a sample size of 50,000 for the HDMR method and 200,000 for all other methods. (b) Converged values of first-order SA measures estimated by each method for the PD uniaxial18 model using a sample size of 50,000 for the HDMR method and 200,000 for all other methods.

Table 10. Dominant parameters for PD uniaxial 18 ranked in order from most (I-1) to least (I-5) significant by each method

5. Discussion

5.1. Convergence

As expected, the screening methods converged faster while identifying the same most influential parameters as the quantitative methods, highlighting their primary advantage. Both methods were significantly less computationally demanding than the Sobol’ and EFAST methods for all cases. When comparing the screening methods to the HDMR method, both methods were more efficient for case 1 and equally efficient for case 2, however, both required double the amount of data for case 3 convergence. It should be noted that the input to the HDMR method is the dataset obtained directly from the FE analysis, and the code creates its own surrogate model for GSA. Meanwhile, all other methods use the same ANN surrogate model for each case formulated according to Section 3.2. When directly comparing the two screening methods, the Morris method is significantly more efficient for case 1, the same for case 2, and slightly less efficient for case 3.

When comparing the quantitative methods, Sobol’ and EFAST were the most computationally demanding of all cases. This result was expected for Sobol’, and EFAST was expected to be more computationally efficient than Sobol’ for first-order results with total order efficiency related to interdependencies. However, the EFAST method required the same number of samples for first-order convergence in case 1 as the Sobol method, double the convergence requirement for case 2, and was only more efficient for case 3, even though minimal parameter interactions were noted for all cases. Regarding total indices, EFAST was substantially higher than both other quantitative methods for cases 1 and 2, and lower than Sobol’ but significantly higher than HDMR for case 3. HDMR required the lowest number of samples in all cases, except for first-order case 1 convergence which was equivalent to Sobol’. As previously noted, HDMR used a different surrogate model than the other methods.

Comparing the first and total order convergence results, it is shown that the first order measures converged equal to or faster than the total order measures, except for the Sobol’ method. Because the total order measure is a sum of the first order measure and all interdependent measures, The SALib library version of the Sobol’ method calculates the first and total order measures simultaneously using identical sample matrices, eliminating implementation as the source of this finding. It was observed that the EFAST method required the same number of samples to converge for both first and total order for each of the three cases, where the HDMR method converged more rapidly for the first order than the total order GSA measures for case 1 and 3 and equally for case 2.

The trend for CoC increased as a function of the number of parameters, with two exceptions caused by numerical instability. For Case 3, the Morris method exhibited a single deviation for one of the dominant parameters with an error >10% from the converged sensitivity value in the 7500–10,000 sample range (Figure 2), resulting in the Case 3 CoC to be greater than the Case 1 CoC as a result. Instability in convergence was also noted in the Sobol method for cases 1 and 3 where the dominant parameters oscillated about the converged value after 5000 samples until consistently predicting the converged value at 50,000 samples for both cases.

5.2. Parameter ordering

Tables 8–10 show that all methods predicted the same most dominant parameters for all cases. For case 1 (Table 8), all quantitative methods predicted the same ordering and nearly identical measures. The HDMR method suggested stronger interactions between the parameters indicated by the increased total order measures relative to first-order measures. The screening methods exhibited minor discrepancies from the quantitative methods, namely X4 and X31 were switched by the Morris method and X33 was ranked more influential than X28 by DGSM. These discrepancies were considered to be insignificant as both X4 and X31 had nearly identical measures and X28 and X33 had relatively small measures compared to the top four. As previously stated, a prior UQ study by Arndt et al. (2022) concluded that only the top four parameters significantly influenced the prediction of damage tolerance. Additionally, the predicted parameters align with the Arndt et al findings that aluminum and composite shear plasticity are the dominant failure mechanisms. The top three parameters represent aluminum properties and the next three are shear plasticity properties of the resin.

All methods predicted a single dominant parameter, X0 (the Young’s modulus of the aluminum) for case 2 (Table 9) followed by the Poisson’s ratio for the aluminum and the critical strain energy release rate for the adhesive resin (X1 and X8, respectively). From an engineering perspective, the deformation in the aluminum caused by the effective bending stiffness of the adherends (X0) greatly affected the stress concentration at the delamination front overwhelming the importance of the Mode I fracture energy of adhesive (X8). Surprisingly the Poisson’s ratio (X1) and X8 were ranked as similarly important with all methods except DGSM predicting X1 as slightly more influential. However, all parameters after X0 had a nearly identical and significantly low measure relative to X0.

For case 3, all methods predicted identical dominant parameters except for the DGSM method (Table 10). DGSM results for this case indicated the largest discrepancy in the entire study. Parameters X15 and X0 were switched and assigned nearly identical measures while all other methods ranked X15 first by an appreciable margin over X0. Additionally, DGSM predicted X16 in place of X13 which is considered insignificant given the low measures after the third most influential parameter. As this case was chosen to investigate a problem with modeling parameters in addition to material properties, it is interesting to note that the top two parameters are modeling parameters, the hourglass coefficient (X15) and point volume (X0). It was not surprising that the point volume would have a significant effect on the stiffness of the specimen as this parameter effects the computed density in the energy-based formulation. The sensitivity to X0 was also noted during model construction.

5.3. Method selection

As methods predicted similar results when identifying the most dominant parameters, method selection is dependent on the information required from the GSA and available computational resources. The Sobol’ method requires the most computational resources in general while providing the option for higher-order investigation. When this capability is needed, one option to significantly reduce computational demand is to first use a screening method to eliminate parameters with little influence on the output. The Sobol’ method can then be performed on the reduced parameter space to obtain more detailed, higher-order information on the remaining parameters. The Morris method exhibited a much faster rate of convergence for all test cases while identifying the same dominant parameters as the quantitative methods. Therefore, this method would be efficient at reducing the parameter space prior to quantitative GSA and especially effective when long-running, high-fidelity models are used to generate the dataset. Other screening applications are identification of the significant parameters to inform focused design, analysis, optimization, and UQ.

Both the Morris method and the Sobol’ method require a specific sampling method to be optimally efficient. Therefore, method selection must be determined prior to data generation directly from the high-fidelity models or the sample data determined from a surrogate model following the sample scheme. For existing data sets, the HDMR method uses any random or quasi-random sampling method and converges relatively quickly compared to the other quantitative methods. Therefore, the HDMR method can be applied to experimental test data or other cases with alternative data sampling. However, the HDMR method is limited by RAM capacity, prohibiting large sample sizes without access to HPC. Both the DGSM and the EFAST method performed adequately but with no clear advantage and being limited by sample size, respectively.

5.4. Limitations

The present study is limited to three test cases with simulation times that could be feasibly sampled using a typical desktop computer. Two common analysis methods were investigated, however, the field of structural mechanics encompasses methodology ranging from analytical solutions to complex high-fidelity modeling. Different mathematical formulations, higher dimensionalities and parameter interactions, and output requirements may necessitate a more advanced or specialized GSA method. An open-source code was utilized due to its availability but limits the study as other GSA packages were not compared. For example, noteworthy methods not included in this study are the Bayesian method (Oakley and O’Hagan, 2004), the generalized Sobol’ method for models with multiple outputs (Gamboa et al., 2014), and the Cramér von Mises method (Gamboa et al., 2018), where the last two of these methods can be explored using the open-source Python package UQpy (Olivier et al., 2020). Additionally, the DAKOTA software developed and maintained by Sandia National Laboratory (Adams et al., 2021) provides many options for sensitivity analysis, UQ, and optimization. Hardware limitations were noted with the HDMR method, where limited RAM capacity hindered the convergence study. As the goal of this study was to investigate analysis that could be performed on a typical desktop, HPC was not utilized and would remove this limitation.

6. Conclusions

The full process to perform GSA was presented in a general manner for nonexperts to inform method selection by compiling resources for additional study, highlighting a user-friendly Python package that can be implemented on a typical desktop, and encouraging best practices for computational structural mechanics problems. It was demonstrated that method selection is dependent on computational resources, information required from the GSA, and available data. All methods provided comparable results as a screening method to identify the most dominant parameters, and the more advanced, quantitative methods provided insight into parameter interactions by providing both first and total indices. The Morris method was advantageous as a screening tool using the estimated total order GSA measures. This method predicted the same dominant parameters as the others at a much smaller computational cost, which is especially beneficial when data acquisition is computationally intensive. The DGSM screening method also predicted the same dominant parameters as the other methods with some deviances in the ordering, which would be irrelevant if the objective was to only identify the most influential parameters. DGSM did require a much higher number of samples for Case 1 convergence, the most complex model with the highest number of parameters. Sobol’, EFAST, and HDMR all provided comparable indices for both first and total order GSA measures, with Sobol’ the most computationally expensive in general. However, Sobol’ offers the capability to explore higher-order interactions when necessary for more complex parameter spaces. The HDMR method is advantageous for existing datasets. It did exhibit a RAM limitation which could be overcome by performing GSA using HPC. While the EFAST method did provide comparable ordering and indices to the other quantitative methods, it did not exhibit the expected reduction in computational cost.