Uncertainty quantification of elastic material responses: testing, stochastic calibration and Bayesian model selection

Motivated by the need to quantify uncertainties in the mechanical behaviour of solid materials, we perform simple uniaxial tensile tests on a manufactured rubber-like material that provide critical information regarding the variability in the constitutive responses between different specimens. Based on the experimental data, we construct stochastic homogeneous hyperelastic models where the parameters are described by spatially independent probability density functions at a macroscopic level. As more than one parametrised model is capable of capturing the observed material behaviour, we apply Baye theorem to select the model that is most likely to reproduce the data. Our analysis is fully tractable mathematically and builds directly on knowledge from deterministic finite elasticity. The proposed stochastic calibration and Bayesian model selection are generally applicable to more complex tests and materials.


Introduction
The study of material elastic properties has traditionally used deterministic approaches, based on ensemble averages, to quantify constitutive parameters [36]. In practice, these parameters can meaningfully take on different values corresponding to possible outcomes of the experiments. The art and challenge of experimental setup is to get as close as possible to the ideal situations that can be analysed mathematically. From the mathematical modelling point of view, stochastic representations accounting for data dispersion are needed to improve assessment and predictions [9,14,18,26,45,50,61,71].
Recently, stochastic models (described by a strain-energy density) were proposed for nonlinear elastic materials, where the parameters are characterised by probability distributions at a continuum level [37,[65][66][67][68][69]. These are advanced phenomenological models that rely on the finite elasticity theory [15,48,76] and on the maximum entropy principle to enable the propagation of uncertainties from input data to output quantities [62]. The principle of maximum entropy was introduced by   [19,20] (see also [21]) and states that "The probability distribution which best represents the current state of knowledge is the one with largest entropy (or uncertainty) in the context of precisely stated prior data (or testable information)." The notion of entropy (or uncertainty) of a probability distribution was first defined by Shannon (1948) [57,64] in the context of information theory. These models can further be incorporated into Bayesian approaches [3,30] for model selection and updates [37,45,46,56].
For non-deterministic material models, two important questions arise, namely "how do material constitutive laws influence possible equilibrium states and their stability?" and "what is the effect of the probabilistic parameters on the predicted elastic responses?" Presently, theoretical approaches have been able to successfully contend with cases of simple geometry such as cubes, spheres, shells and tubes. Specifically, within the stochastic framework, the classic problem of the Rivlin cube was considered in [38], the static and dynamic inflation of cylindrical and spherical shells was treated in [34,35], the static and dynamic cavitation of a sphere was analysed in [33,40], and the stretch and twist of anisotropic cylindrical tubes was examined in [39]. These problems were drawn from the analytical finite elasticity literature and incorporate random variables as basic concepts along with mechanical stresses and strains. Such problems are interesting in their own right and may also inspire further thinking about how stochastic extensions can be formulated before they can be applied to more complex systems.
To investigate the effect of probabilistic parameters in the case of more realistic geometries and loading conditions, computational approaches have been proposed in [68,69]. However, for real materials, most available data consist of mean values from which deterministic models are usually derived, and there is a lack of experimental data reported in the literature that are directly suitable for stochastic modelling.
For rubber-like materials, the first experimental data for load-deformation responses of different material samples under large strains were provided by Rivlin and Saunders [55]. Following the phenomenological models for vulcanised rubber developed by Treloar (2005) [75] (read also [13]), many deterministic hyperelastic models calibrated to the mean data values were described (for example, in [7,8,17,36,49,70,77]). To capture the variability of Rivlin and Saunders' data (see also [75, p. 224], or [76, p. 181]), probability distributions for the random shear modulus under relatively small strains were obtained in [34]. In recognition of the fact that a crucial part in assessing the elasticity of materials is to quantify the uncertainties in their mechanical responses, for rubber and soft tissues under large strain deformations, explicit stochastic hyperelastic models based on datasets consisting of mean values and standard deviations were developed in [37], while statistical models derived from numerically generated data were proposed in [6,44]. In all these developments, the aim was to find the most appropriate general method of characterising the properties of nonlinear elastic materials at a macroscopic (non-molecular) level.
In this study, we first report on simple experimental results for different samples of silicone rubber material in finite uniaxial tension (Section 2). Full details of the manufactured material are presented in Section 2.1, while specific descriptions of the experimental setup and techniques used during testing are outlined in Sections 2.2 and 2.3. Employing the stochastic calibration method and Bayesian selection criterion proposed in [37], we then obtain stochastic homogeneous hyperelastic models, where the parameters are random variables that are constant in space, and are calibrated to our experimental measurements (Section 3). The assumptions and ideas that we adopt for the stochastic modelling are discussed in Section 3.1, while calibration results and details on how Bayes' theorem can be employed explicitly to select the best performing model are given in Sections 3.3 and 3.4. In Section 4, we draw concluding remarks.

Experimental testing
In this section, we describe the experimental setup and techniques used to measure large elastic deformations of manufactured silicone rubber specimens. Our experimental results capture the inherent variability in the acquired data between different specimens during tensile tests.

Specimen manufacture
Two batches of silicone were manufactured in order to investigate the consistency of the specimen behaviour within the same batch and within different batches. Each sample was a simple rectangular shape with rounded edges. The approximate (a) (b) Fig. 1 Marked undeformed specimens, with height a of 100 mm, width of 10 mm, and depth of 4 mm: a single specimen where the longer double arrow shows the grip length, which was 30 mm as standard and reduced to 20 mm for some tests, while the shorter arrow indicates the gauge length, which was 20 mm as standard and 10 mm when using the reduced grip length; b multiple specimens geometric parameters of each sample before deformation were height of 100 mm, width of 10mm, and a depth 4 mm. Examples of undeformed samples are illustrated in Fig. 1.
• For Batch 1, tensile testing specimens were cast using Tech-Sil 25 Silicone (Technovent). This is a two part silicone, with a standard mixture ratio of 9:1 for Part A:Part B, respectively, as per the manufacturer's recommendation, and is generally allowed to cure at room temperature for at least 24 h. The silicone was mixed and de-gassed prior to casting, ensuring an even mixture, and no air bubbles were present in the tensile specimens. Testing specimens of equal dimensions were made within a mould that is typically used for this purpose. The silicone was removed from the mould for testing after 4 weeks. • For Batch 2, the same make of silicone as in Batch 1 was used, with slight variations in the mixture components in order to simulate an error that would be within a realistic experimental range. Namely, the mixture ratio was 8.96:1. Tensile testing specimens with the specified geometry were created using this mixture. The silicone was also left to cure at room temperature, and taken out of the mould for testing after 2 weeks.
The full details of the silicone mixture for the two batches are recorded in Table 1.

Experimental set-up
Two testing sessions were conducted, to allow for variability in the experimental results due to degradation of silicone properties over different testing days. Uniaxial tests were conducted on the Zwick Roel Z050 testing machine, with a 1 KN load cell to measure tensile force. Specimens were mounted using a set of roll clamps, and the general experiment setup can be observed in Fig. 2. The standard test method used consisted of a pre-load of 2 N, a 30-mm grip length of the specimens, and loading at a speed of 30 mm/min. However, testing parameters, such as the grip length and testing speed, were varied between specimens, as detailed in Table 2. Tests were stopped once specimens reached approximately 100% strain measured optically (see Section 2.3 for details). All specimens were tested with the same technician loading the test specimens and using the test machine.

Optical strain measurement
A video strain gauge system (Imetrum) was employed to capture the global deformation of the specimen during tensile tests. The system works by imaging the specimen in an unloaded state, and then tracking the position of markers on the surface of the specimen throughout loading to measure displacements [10,72,78,79]. The system was used with a single camera with a general purpose lens, and calibrated using markers of a known distance apart within the field of view, as per  the manufacturer's instructions. When processing the captured video, a digital strain gauge was placed between two marker points. Specimens were marked using a permanent marker in the positions indicated in Fig. 3, so the digital gauge length could vary. An example of a marked silicone specimen is shown in Fig. 1. When processing the data, the target size for each marker was varied. The target area defines the area in which the software locates the surface marker. Additionally, within the software, the user has the option to control certain elements of the tracking algorithm, for example, in order to account for the stretch and rotation in specimens subjected to large deformations. The post-processing variations are outlined in Table 3, and were allowed for the purpose of accounting for varying user preferences.
The data values measured and recorded during our experimental tests capture the inherent variation between the constitutive behaviour of the different specimens listed in Table 2. The values for the applied force versus maximum vertical displacement in tensile loading for the individual specimens are shown in Fig. 4. Note that, for the rubber-like material, the force needed to be increased to a critical level before non-zero macroscopic deformations could be observed. For the specimens in each batch, the data are represented together in Fig. 5 where the average values, which are typically used for the calibration of deterministic models, are marked by red lines. In this figure, different quantities of interest are also illustrated. Specifically, in addition to the applied force versus maximum vertical displacement observed experimentally, we have calculated the following nonlinear quantities: • The first Piola-Kirchhoff (PK) tensile stress, representing the force per unit area in the reference configuration, where F is the applied tensile force and A = 40 mm 2 is the cross-sectional area.  Table 2 • The nonlinear stretch modulus [36]  where a is the stretch ratio and T = aP is the Cauchy stress (representing the force per unit area in the current configuration), with P the first Piola-Kirchhoff stress given by Eq. 1. • The nonlinear shear modulus, given by the universal formula [36] where E is the stretch modulus defined by Eq. 2.
For formal definitions and explicit derivations of key nonlinear elastic moduli in homogeneous isotropic finite elasticity, and their universal relations under large strains, we refer to [36].

Stochastic modelling
Next, we construct specific stochastic homogeneous hyperelastic models where the parameters are characterised by spatially independent probability distributions and optimised to the collected data.

Stochastic isotropic incompressible hyperelastic models
We recall that a homogeneous hyperelastic material is defined by a strain-energy function W (F), depending on the deformation gradient, F, with respect to a reference configuration [15,48,76]. Traditionally, this is characterised by a set of deterministic model parameters, which contribute to the constant elastic moduli under small strains, or to the nonlinear elastic moduli, which are functions of the deformation under large strains [36]. In contrast, a stochastic hyperelastic model has parameters which are defined by probability density functions [37,[65][66][67]. Typically, each model parameter is described in terms of its mean value and its variance, which hold information regarding the range of values about the mean value [5,6,18,29,44]. Here, we construct stochastic homogeneous isotropic incompressible hyperelastic models that rely on the following assumptions [33-35, 37, 38, 40]: where {λ i } i=1,2,3 and {T i } i=1,2,3 denote the principal stretches and principal Cauchy stresses, respectively, such that In Eq. 4, the strict inequality '>' is replaced by "≥" if any two principal stretches are equal. (A4) For any given finite deformation, at any point in the material, the shear modulus, μ, and its inverse, 1/μ, are second-order random variables, i.e., they have finite mean value and finite variance [65][66][67][68][69].
In practice, elastic moduli can take on different values, corresponding to possible outcomes of experimental tests. The maximum entropy principle allows us to explicitly construct the probability laws of the model parameters, given the available information. Approaches for the explicit derivation of probability distributions for the elastic parameters of stochastic homogeneous isotropic hyperelastic models calibrated to experimental data for rubber-like material and soft tissues were proposed in [37,67].

Hypothesis testing
Before, we attempted to construct models based on the collected data, we applied standard statistical tests [11] to verify whether we can treat the entire dataset as one, or we must model the dataset of each batch separately. First, we used an unpaired t test [52] to compare the shear modulus data at small strain for the two batches. This test provided a p value of approximately 10 −12 , which was less than the 'standard' 0.05 significance level, and therefore allowed us to reject the null hypothesis that the moduli of the two batches come from the same distribution. The significant differences between the data corresponding to the two batches are illustrated in Fig. 6 . Further, we used a χ 2 (chi-square) goodness-of-fit test [51] to check that we cannot reject the null hypothesis that the shear modulus data of each batch come from a gamma distribution. The p values for the batch 1 and 2 fits were 0.1 and 0.07, respectively. Thus, we could not reject the hypothesis that the shear moduli are gamma-distributed. Figure 7 illustrates the gamma and normal (Gaussian) probability distributions fitted to the shear modulus data at small strain. We note the similarity between the represented gamma and normal distributions. This is not accidental because, when ρ 1 → ∞, the gamma probability distribution Eq. 8 is approximated by a normal distribution (see Appendix A for a detailed proof). Table 4 provides the fitted variables for both distributions. As we see, ρ 1 is very large compared to ρ 2 , justifying the   Table 4 similarity between the two distributions. However, elastic moduli cannot be characterised by the normal distribution since this distribution is defined on the entire real line whereas elastic moduli are typically positive. Thus, the gamma distribution is chosen for its consistency with the mechanical properties of elastic materials.

Stochastic calibration
Applying the stochastic method developed in [37], we specialise to three different stochastic homogeneous isotropic strainenergy functions which we calibrate to the collected data. These strain-energy functions are presented in Table 5, while their constitutive parameters fitted to the experimental measurements for batch 1 and batch 2 are recorded in Table 6. Note that the calculated shear moduli are comparable to those recorded in Table 4. The load-deformation results are plotted in Figs. 8 and 9, demonstrating that all three models perform well, but exhibit different levels of accuracy when compared to the actual data.

Bayesian model selection
In this section, we employ Bayes' theorem [3,30] to select the best perfoming model. We denote by P (M) the prior probability of a model M before the data values D are taken into account, and by P (D|M) the likelihood, or the probability of obtaining the data values D from the model M. Bayes' theorem states that where P (M|D) is the posterior probability of the model M, and P (D) is the normalisation value, known as the marginal likelihood. This theorem provides also a methodology for estimating the odds for a model M (i) to another model M (j ) in light of the data D, i.e., where is known as the Bayes factor. Formula (11) states that the posterior odds O ij for the model M (i) against the model M (j ) , given the data D, are equal to the prior odds multiplied by the Bayes factor. In particular, if the models have equal prior

Mooney-Rivlin
Cp αp 1−a 3 [47] C p independent of deformation; probabilities, P (M (i) ) = P (M (j ) ), i.e., there is no prior favourite, then, by Eq. 11, the posterior odds are equal to the Bayes factor, i.e., O ij = B ij . However, if the Bayes factor is equal to 1, then Occam's razor [22][23][24]73] would imply that a larger prior probability should be assigned to the simpler model than to the more complex one for reasons of parsimony.
Maintaining a general framework, we assume P (D|M) to be an arbitrary probability that is symmetric about the mean value D = 0 and decreasing in the absolute value of D (see also the imaginative cartoon by D. M. Titterington (1982) [74]). In this case, the Bayes factor B ij satisfies the inequality as follows [4]: where D (i) and D (j ) designate the standard deviation that the predicted quantity of interest computed with the model M (i) and M (j ) , respectively, deviates from the observed data value D. The formula for calculating D (i) is as follows, where, for the quantity of interest, Q (i) is the expected value computed with the model M (i) , and D and D represent the experimentally observed mean value and standard deviation, respectively. The expression on the right-hand side of Eq. 13 provides an explicit lower bound on the Bayes factor B ij . Thus, taking equal prior probabilities, the lower bound on the Bayes factor given by Eq. 13 represents a lower bound on the posterior odds. This lower bound is an estimate of the amount of evidence against the model M (i) , i.e., the maximum support for the model M (j ) provided by the data. Table 6 Parameters of stochastic constitutive models given in Table 5 Table 6, showing: a the first Piola-Kirchhoff tensile stress, b the nonlinear stretch modulus, c the nonlinear shear modulus, and d the relative error for the shear modulus mean values We now apply these bounds to select the best performing model among the models calibrated to our experimental data. For example, at a = 1.15, for each model recorded in Table 6, we calculate the standard deviations that the mean shear modulus μ deviates from the experimental mean data value 0.2909, given that the experimental standard deviation is 0.0170. Applying formula (14), we obtain the following: (b) (c) (d) Fig. 9 Stochastic models calibrated to the batch 2 data, with the parameters recorded in Table 6, showing: a the first Piola-Kirchhoff tensile stress, (b) the nonlinear stretch modulus, c the nonlinear shear modulus, and d the relative error for the shear modulus mean values taking P (D|M (1) ) = 1 − P (D|M (3) ), by the lower bound on the Bayes factor B 31 , the likelihood of obtaining the data with the Ogden model is P (D|M (3) ) ≥ 0.6634. Similarly, assuming equal prior probabilities and taking P (D|M (2) ) = 1 − P (D|M (3) ), the lower bound on B 32 implies P (D|M (3) ) ≥ 0.6179. Therefore, the data at a = 1.15 are more likely to be reproduced by the Ogden model than by the other two models. In Fig. 10, we illustrate the lower bounds on the Bayes factor B ij for each model M (i) against another model M (j ) at various stretch ratios. From these bounds, we infer that the Gent-Gent model is generally more likely to represent the data than the Mooney-Rivlin model, and the Ogden model is more likely than both the Mooney-Rivlin and the Gent-Gent models. These results are consistent with those in Figs. 8d and 9d where the relative errors for the mean values of the shear modulus are shown. The relative error was calculated in the usual way, i.e., |μ model − μ data |/μ data , where μ model is the mean value of the shear modulus for the respective model and μ data is the mean value of the shear modulus calculated from the experimental data.

Conclusion
We report on experimental tests on different samples of a manufactured rubber-like material under large tensile loading, and employ a stochastic strategy to derive constitutive models that take into account the variability in the collected data. Specifically, we construct isotropic incompressible hyperelastic models with model parameters defined as spatially independent random variables characterised by probability density functions at a continuum level. In addition, we provide a methodology where an explicit lower bound on the Bayes factor is used to compare different models, and then applied to select the model that is most likely to reproduce the data. Ideally, these models are calibrated and validated on multiaxial test data [31,32,41,63]. In this case also, our stochastic calibration and Bayesian model selection can be employed to obtain suitable models.
Our analysis is fully tractable mathematically and builds directly on knowledge from deterministic finite elasticity. This study highlights the need for continuum models to consider the variability in the elastic behaviour of materials at large strains, and complements our previous theoretical investigations of how elastic solutions of fundamental problems in nonlinear elasticity can be extended to stochastic hyperelastic models [33][34][35][38][39][40].
Thus, the limiting moment generating function of Y when ρ 1 → ∞ takes the form as follows: The above limit can be calculated as follows: Therefore, which is the moment generating function of a normal-distributed random variable. Thus, the limiting distribution of a gamma distribution with ρ 1 → ∞ is the normal distribution. This completes the proof.