A Bayesian Multi-Fidelity Neural Network to Predict Nonlinear Frequency Backbone Curves

Digital engineering leads to more reliance on simulation models to evaluate the design performance of engineered structures. Computational models can be leveraged to evaluate several design options to mature the product definition prior to manufacturing and testing the design. By learning and failing in the digital engineering space, better design choices are realized prior to entering into the expensive and iterative build-test-redesign cycle. Early in a design process, designs change rapidly and many of the details are ill-defined. For this reason, low-fidelity, parametric simulation models offer a valuable design tool to evaluate quantities of interest for a wide parameter space. Low-fidelity models are computationally efficient but are limited in their predictive accuracy, thus they are most useful for relative comparisons. As the design cycle progresses and more definition is provided, the model fidelity tends to increase as well. High-fidelity simulation models require significant time investment to generate, due to the complexity needed to emulate the physical hardware. These models are developed in preparation for more detailed analysis, e.g., stress analysis during mechanical environments, as well as validation efforts. The high-fidelity models can be used for qualification activities to provide a technical basis to demonstrate that the designed structure will meet its requirements and survive the operational and environmental loads.

The challenge with this simulation model framework is that often the high-fidelity model is accurate, but inefficient, thus limiting the number of queries to explore the parameter space. The low-fidelity simulation models offer the contrary, where their efficiency allows the evaluations of many parameter samples but lack the resolution needed to capture the physics accurately. For these reasons, multifidelity, parametric surrogate models have been developed to be both efficient at sampling the parameter space, and as accurate as the high-fidelity simulations. Multifidelity sampling has been used to improve the performance of Monte Carlo methods in estimating statistics from data [1–4]. Beyond the estimation of statistics from processes with varying fidelity, researchers have leveraged these data to train surrogate models. The surrogate model is often formulated as a nonlinear autoregressive model based on the linear fusion scheme proposed by Kennedy and O'Hagan [5]. This concept has been extended to capture nonlinear correlations where the problem is then formulated as finding the unknown coefficients of nonlinear autoregressive model through optimization or training of machine learning models such as Gaussian processes [6,7] or neural networks [7–11]. Multifidelity surrogate modeling has been successfully applied to a wide range of applications as reported in Ref. [12] where the authors used the pyapprox package [13] to fuze data from various sources for modeling complex systems. Multifidelity methods have even been extended to deal with dissimilar parameters [14]. The availability of an efficient multifidelity model enables the application of Bayesian methods for calibration of structural dynamics models [15–17] due to the ability of these methods to consider uncertainty. In general, the large body of research into multifidelity modeling has demonstrated the method's versatility and robustness as documented in the comprehensive review provided in Ref. [18].

In 2021, the Tribomechadynamics (TMD) research challenge (TRC) was announced to the research community [19]. The purpose of the challenge was to invite participants to create blind predictions of the nonlinear dynamic response of a mechanical structure that had not been fabricated. The goal is to evaluate the predictive capability of the communities physics-based modeling approaches in the absence of calibration data and physical experiments. As part of the challenge, the organizers released the cad geometry along with mechanical drawings, material/surface specifications, and assembly procedures. Each participant was requested to deliver the amplitude-dependent natural frequency and damping ratio of the first elastic mode of the assembly, in order to characterize the nonlinear structural dynamic performance of the design.

The nonlinear frequency and damping backbone curves requested as deliverables from the research challenge are useful in describing the change in resonant behavior of the structure at varying levels of vibration response. Classical structural dynamics theory is rooted in linear system theory and describes the dynamics of the system using modal analysis, which is the characterization of the invariant mode shapes, natural frequencies and damping ratios. In reality, there are several nonlinear physics common to mechanical systems that violate the assumption of linearity, including frictional contact, geometric nonlinearity due to large deformations, and material constitutive laws. The nonlinear frequency and damping backbones are rooted in nonlinear normal mode (NNM) theory [20,21], which provides a theoretical foundation to describe the amplitude, or energy, dependence of the modal characteristics. Several methods are available to estimate the backbone curves both from computational mechanics models [22] or from experiments [23–25]. The goal of the TRC is to evaluate the current methods from the research community to simulate the nonlinear physics associated with the design as well as the physical models used to represent the nonlinear physics.

In this research, we develop a parametric, multifidelity surrogate model that efficiently and accurately predicts the nonlinear frequency backbone curve for the TRC benchmark structure. The authors participated in the research challenge previously mentioned, and through the course of the effort, the team developed both a low- and high-fidelity model to capture the geometry of the structure and nonlinear physics expected to be present during the vibration response. The multifidelity surrogate model is trained from these two models of varying fidelity and is readily simulated to predict and compare the response to experimental measurements. This paper extends multifidelity surrogate models to operate with nonlinear frequency-dependent response amplitude data while previous work had focused on linear models, or scalar functions. Furthermore, in contrast with previous work, this research employs a Bayesian neural network for direct probabilistic inference with the surrogate model whereas to the best of the authors' knowledge, had not been done previously in this context. This research describes the development of this model and demonstrates its accuracy and usefulness in comparing to experimental data measured on the manufactured hardware.

The remainder of the paper is organized as follows: Sec. 2 provides a brief background of the Tribomechadynamics Research Challenge and the associated structural design and experiments. Section 3 introduces the multifidelity physics-based models developed to perform multiquery simulations needed to train the surrogate model. Section 4 describes how the multifidelity model results are leveraged to develop a surrogate model and it's uses for Bayesian calibration and uncertainty quantification. Numerical results demonstrating the accuracy of the surrogate model are presented in Sec. 5, followed by the conclusion in Sec. 6.

The work presented in this paper is based on the benchmark structure introduced during the Tribomechadynamics (TMD) Research Challenge (TRC). A brief summary of the challenge and experimental campaign is described in Secs. 2.1 and 2.2.

The TRC, first announced to the research community in 2021, invited researchers from different institutions to make blind simulation predictions of the nonlinear resonant behavior of a proposed structural assembly that simultaneously exhibits nonlinearity from both large deformations and bolted joints. The structure consists of four main parts: a thin panel, a support structure, and two blades to compress the ends of the panel with six bolts on each end. A schematic of the proposed design is shown in Fig. 1 and the released data can be found in the “TRC Challenge—design documents” project repository [19]. The thin panel is nominally a 300 mm × 100 mm flat plate with a 1.5 mm thickness and made of stainless steel (AISI 301 material number 1.4310). The panel is bolted between the blades and pillars on the support structure with six M6 bolts on each end and washers installed underneath the bolt heads. The rear support is designed to be bolted onto a slip table of a large shaker with a total of ten M10 bolts and washers. The support structure and the two blades are made of hardened steel (material number 1.7147). The thin plate is designed to be nominally flat but a slight angle of 2.2 deg at the support pillar surface produces an arch in the panel when installed, resulting in a nonflat panel with some residual stress when in its assembled state. The structure setup is meant to introduce contact nonlinearity via the frictional contact between the panel, blades, and support, and a geometric nonlinearity due to coupling in the axial stretching and bending of the thin panel. It was expected that the two nonlinearities coupled through the axial stretching of the panel and frictional contact resisting the in-plane deformations.

A total of eight research institutions participated in the TRC, where the main objective was to predict the nonlinear frequency and damping backbone curves of the lowest frequency (elastic) mode as a function of response amplitude at the center of the thin panel. Different methodologies were pursued by different research groups, including the authors of this article, and the summary of the predictions are documented in detail in Ref. [26]. In addition to the collection of blind predictions, the results are compared against experimental measurements in effort to assess the accuracy of the different proposed techniques and developed models. A brief summary of the experimental efforts are described in Sec. 2.

The benchmark structure was manufactured, assembled and tested during the Tribomechadynamics Research Camp during the summer of 2022. The experimental setup is shown in Fig. 2, which was designed to excite the lowest-frequency mode via base excitation applied to the support structure by bolting it onto a slip table of a large shaker. The amplitude-dependent frequency and damping ratio of the targeted mode were obtained by performing nonlinear phase-resonance tests using the approach in Ref. [27]. The full details of the experiment along with other measurements are provided in Ref. [28].

The researchers who conducted the experiments noted that the measurements were sensitive to temperature variations which introduced additional variability in the response. It was also observed that the nominally flat thin plate had an initial curvature prior to installation. They measured a 1.2 mm deviation from the flat plate and they attributed the cause to the rolling manufacturing process. While these details are not included explicitly in the models that will be described in the next section, they will be helpful in understanding the results presented later in Sec. 5.

The authors developed two different physics-based models for the research challenge, as discussed in Ref. [29]. A low- and high-fidelity model were created to account for the geometric and frictional contact nonlinearity in the system, however only the high-fidelity model was officially submitted as the blind prediction for the research challenge [26]. The high-fidelity finite element model (FEM) is presumed to be highly accurate due to the explicit and detailed representation of the geometry and the physics. On the other hand, the low-fidelity model is a simplification of the geometry and physics in the structure, such as contact pressure distribution and the exclusion of three-dimensional stress states. An initial sensitivity analysis was performed using the low-fidelity model, which revealed that the angle and thickness of the thin panel, as well as the coefficient of static friction, were the most sensitive parameters to the nonlinear backbone curves. More detail on this analysis is provided in Sec. 4.2. As a result of this finding, these three variables were chosen to parameterize the high-fidelity FEM.

where M is the mass matrix, $ϕ$ is the mass normalized mode shape vector of interest, and α is a scalar to define the amplitude of the body force. The modal body force in Eq. 1 is applied to all degrees-of-freedom in the finite element model, which approximately loads the nonlinear model into the shape of the mode of interest. After the release of the modal body force, a transient dynamics simulation then simulates the ring-down behavior and the resultant backbone curves can be estimated with time-frequency analysis techniques [25]. Note that it is assumed that the mode is sufficiently isolated from other modes (and that no modal coupling occurs). In this research, backbone curves were extracted using the short-time Fourier transform (STFT) of the response, as described in Ref. [33]. Since the modal body force approximately isolates the ring-down response with a monoharmonic signal, the use of time-frequency analysis techniques are appropriate to estimate the amplitude-dependent frequency and damping of the structure [34].

The following steps summarize the simulation workflow:

1. Quasi-static preload: the bolted interfaces are preloaded to solve for the initial preloaded equilibrium state.

2. Modal analysis: an eigen-analysis of the linearized preloaded model is performed.

3. Quasi-static modal force: the preloaded model is subjected to the application of a modal body force in order to deform the structure in the shape of the mode of interest. This is the modal body force defined with the QSMA technique in Refs. [30–32].

4. Dynamic ring-down: after the modal force application, all external forces are removed to allow the model to respond in free vibration.

5. Post-processing: the nodal displacement time histories from the free response simulation are post-processed to obtain the amplitude-dependent frequency and damping ratio for the mode of interest using the STFT method in Ref. [33].

The high-fidelity finite element mesh of the TMD structure was developed in CUBIT [35] to capture the full three-dimensional complexity of the assembly. The close-up image in Fig. 3 shows the specific detail of the blade, panel and pillar interaction. Sierra solid mechanics (Sierra/SM) [36] and Sierra structural dynamics (Sierra/SD) [37] finite element software was used to perform the high-fidelity simulations. A selective deviatoric element formulation was used for all the analyses with the high-fidelity model. This element formulation corresponds to a fully integrated element with respect to the deviatoric stress response and under-integrated with respect to the hydrostatic pressure response [36]. Initial preloading due to bolt tightening was performed by applying an artificial strain to the bolt shanks required to reach the nominal preload force using a quasi-static analysis in Sierra/SM. Following this approach, the preloaded model was linearized about the equilibrium state to perform a linearized modal analysis in the Sierra/SD code. The mass normalized modes were computed and used to output the body forces in Eq. (1), which were then applied to the preloaded system in Sierra/SM. The body force corresponding to the first panel bending mode was applied rapidly via a cosine ramp function and then released to initiate the free transient response. The frequency of the linearized model varied between 112 and 125 Hz depending on the value of the model parameters that were varied. The last step in the simulation was to obtain the nonlinear ring-down response using an explicit time integration scheme with a critical time-step of roughly 50 nanoseconds.

The low-fidelity model of the curved plate was developed in matlab using a two-dimensional corotational beam element formulation [38,39] to capture the geometric nonlinearity. Two-dimensional Jenkins elements [40] with variable normal load were added to the beam elements at the location of the bolted interface in order to capture the effects of frictional contact. These friction elements were able to be preloaded to capture the initial preload force in the model. It was assumed that the normal force from the bolt preload was uniformly distributed, which is a known limitation due to the more complex contact pressure distribution observed in the high-fidelity model. A schematic of the low-fidelity model is shown in Fig. 4. The low-fidelity model consists of only 136 beam elements, making it efficient to perform the quasi-static preload, linearized eigen-analysis, and nonlinear transient ring-down analysis. One distinction of the approach developed for the low-fidelity model is that the preloaded model was loaded with the modal body force using a quasi-static analysis, then these initial conditions were used for the transient dynamic analysis without any external load.

The preceding section described the two multifidelity physics-based models used to simulate the nonlinear modal characteristics of the structure. Due to the multiple steps involved in the simulation sequence, the workflow can become computationally expensive, in particular when running the high-fidelity FEM at the resolution required to obtain accurate predictions. To this end, a parametric, multifidelity surrogate model is developed to obtain fast predictions of the nonlinear frequency backbone curve. A fast surrogate model enables parametric uncertainty analysis and Bayesian calibration to existing test data. The following sections describe this process.

There are many relevant physical parameters that influence the response of the TRC structure. However, generating training data with the high-fidelity model quickly becomes very expensive with the addition of new model parameters. To reduce the number of parameters considered, an initial sensitivity study was performed with the low-fidelity model to determine the top three most influential parameters. The sensitivity was assessed by fitting a random forest regression model [41] to 150 Latin Hypercube [42] samples generated with the low-fidelity model and then inspecting the normalized feature importance of the random forest model on a held-out test set (i.e., permutation feature importance). This approach was computationally efficient and informative in down-selecting the most influential features to select for the parameterized models. It should be noted that this approach does not perform high-order sensitivity analysis but it is still useful for feature selection. The parameters considered in this study were the following: material density (rho), Young's modulus, preload force, friction (mu), thickness, and plate angle. Figure 5 illustrates the ranking of feature importance extracted from the random forest algorithm. As shown, the angle, thickness, and friction were the most influential parameters according to this study, so these were the parameters that were selected for the multifidelity surrogate model.

Here, N_lf and N_hf are the sizes of the low- and high-fidelity sample sets, respectively. The variational posterior distribution is represented by $q(w(j)|π)$⁠, which is chosen to be Gaussian in this work and is parameterized by π, which is a vector containing the mean and standard deviation of the distribution. The weights prior distribution, $P(w(j))$⁠, is assumed to be Gaussian as well, and $w(j)$ represents the jth Monte Carlo sample drawn from the variational posterior $q(w(j)|π)$⁠. In Eq. (6), the first two terms represent the Kullback–Leibler (KL) divergence between the variational distribution and the prior, and can be interpreted as a complexity cost (i.e., the difference between the proposed prior and the learned variational distribution). The third term of Eq. (6) represents the likelihood cost, or the goodness of fit to the data. This loss term is an approximation to the variational free energy cost function [45,46]. The total loss function is minimized using the Adam optimizer [47] with a learning rate set to $10−4$⁠.

The training data were generated by varying the blade angle θ (which controls the initial curvature as shown in Fig. 4), the panel thickness t, and the friction coefficient μ. The ranges are provided in Table 1 and were determined based on the tolerances specified in the drawings, except for friction which was determined based on engineering judgment. A central composite design (CCD) sample technique was used for the high-fidelity FEM, resulting in 15 samples which took roughly 100 h on 48 processors for each simulation case. In contrast, a Latin Hypercube sample (LHS) was generated with the low-fidelity model. This LHS set consisted of 150 samples and it took 25 min on one processor for each simulation case. These training samples were combined to train the neural network. An additional five random samples were run with the high-fidelity FEM to define a test set. The distribution of these samples is shown in Fig. 7.

As described in Sec. 2.2, the experiment exhibited variability due to temperature effects as well as due to assembly-to-assembly variations. Although the model did not account for these effects explicitly, the chosen parameters are able to serve as a proxy to capture some of the effects that temperature and re-assembly have on the structure. Namely, temperature variations may affect the curvature and bolt preload which are effectively captured by varying the friction coefficient and blade angle. For these reasons, the distribution of the parameters was updated to better reflect the variability observed in the test data.

Because the relationship between the parameters λ and the amplitude-dependent frequency f given by $F(x,λ)$ is nonlinear, the likelihood function cannot be easily expressed analytically. This fact poses a challenge for obtaining the posterior analytically which is often the case for nonlinear systems. In this work, the posterior is approximated with Markov chain Monte Carlo (MCMC) sampling of the trained multifidelity surrogate model. The Bayesian inference algorithms were implemented in Python using the tensorflow probability (TFP) package. The MCMC samplers in TFP were used to obtain a chain of dependent samples to approximate the posterior distributions of the input parameters (i.e., thickness, friction, and angle). Specifically, the no U-turn (NUT) sampler was used [48], which is part of the Hamiltonian Monte Carlo sampling algorithms.

The high-fidelity training data samples generated with CCD sampling (Fig. 7) showed excellent agreement with the experimental data, as shown in Fig. 8. These results demonstrate that the parameter bounds chosen span the range of the experimental data prior to any model updating. The subsequent Bayesian model updates are meant to refine the assumed uniform distributions but are not expected to deviate dramatically from the assumed uniform distribution.

The multifidelity surrogate model was trained for 30,000 epochs. A range of hyper-parameters were explored and it was found that the following architecture had the best performance: $Gl$⁠: six layers, six units and $Gh,nl$⁠: four layers, six units, with a Swish activation function for all hidden layers, and a linear function for the output layer. The predictions of the test set are shown in Fig. 9. These results show generally good agreement between the predicted frequency and the curves obtained from the high-fidelity model. It was observed that in some cases the actual backbone is just outside the predictive interval. It is expected that adding more high-fidelity training samples would improve the prediction accuracy, but, as explained previously, each training sample comes with additional high computational cost.

Once trained, the surrogate model can be used to visualize the variability of the amplitude-dependent frequency as a function of the normalized input parameters. To illustrate this sensitivity, each parameter was varied one at a time ranging from its lower bound value to its upper bound value (see Table 1) while setting the other two parameters to either $0.1,0.5$⁠, or 0.9 to explore the potential interactions between these parameters. Figure 10 illustrates the results of this sensitivity exploration where the first row corresponds to setting the other parameters to 0.1, the second to 0.5, and the third to 0.9.

Each individual column in Fig. 10 is analyzed to understand the effect of the different parameters on the response. The first column focuses on the thickness variation. For example, it is observed that when the other two parameters are low (0.1), the change from low to high thickness values exhibits the opposite trend than it does when the parameters are set to either medium (0.5) or high (0.9). This result is likely driven by the fact that if the friction coefficient is low, then the joint is allowed to slip earlier and the nonlinearity is driven by joint slip instead of the geometric effects. For the second and third rows, increasing the thickness then increases the bending stiffness of the mode. The second column illustrates high variability due to variations in the friction coefficient. As expected, decreasing the friction shifts the curves downwards which is associated with softening due to the earlier onset of joint slip. The third column shows the effect of the angle. The response variability due to the angle variation seems to be a bit more pronounced when the other parameters are set to low values. As the thickness and friction increase, the angle is less influential to the response, however, it is observed that increasing the angle value shifts the curves upwards which is associated with a more pronounced effect of the geometric nonlinearity.

The sensitivity analysis provides intuition into the effects of each of the parameters and demonstrates how they are inherently coupled to produce the resultant shape of the frequency backbone curve. As described in Sec. 4.3, Bayesian inference can be used to better understand what the statistical distribution of input parameters should be to reproduce the observed experimental data, given the existing model. To this end, the MCMC algorithm was run for 30,000 steps, with 2,000 burn-in steps discarded. The effective number of samples was over 1200 for all parameters which is more than enough independent samples based on the recommendation by Gelman et al. [49]. The posterior parameter distributions are shown in Fig. 11. As shown, the model update was minimal, which may typically imply that the response was not sensitive to variations in the parameters. However, given the sensitivity analysis performed and the results that will be presented, it is thought that this result is indicative of good coverage over the experimental distribution obtained with the parameters uniform distributions. This result is consistent with the good qualitative agreement between the high-fidelity CCD samples and the experimental data shown in Fig. 8.

It was observed that two distinct clusters could be distinguished in the experimental data based on visual inspection of the backbone curves. To explore this behavior, a K-Means algorithm was used to cluster the frequency curves into two groups, and then Bayesian inference was performed on each cluster separately. Performing the inference this way, the algorithm revealed certain trends that can be observed in the histograms of the parameter posterior distributions in Figs. 12 and 13 for cluster 1 (bottom group) and 2 (top group), respectively. The posterior samples for the two clusters are shown in Fig. 14. The posterior distributions suggest that the friction coefficient and the thickness are influencing the response, with a lower friction coefficient shifting the curves down (i.e., more pronounced softening). Similarly, the curves on the top (cluster 2) seem to be associated with a slightly larger friction coefficient and a concentration on larger thickness values which may indicate that the geometric nonlinearity overcomes the joint nonlinearity earlier. These statistical distributions, combined with the intuition drawn from Fig. 10, suggest that separation of these two clusters is driven by the type of nonlinearity that dominates the response. Cluster 1, which is associated with lower friction coefficient values might be dominated by joint slip (i.e., it slips earlier) while cluster 2, which is characterized by a higher friction coefficient and large thickness values, may be dominated by the geometric nonlinearity. Variations in both of these parameters could be related to the temperature or assembly-to-assembly variations mentioned in Ref. [26]. The friction value may act as a proxy for preload force as well, since the preload may also be variable in the real experimental setup and could be influenced by the temperature. The statistical distribution for the angle was hardly updated by the Bayesian inference, suggesting that the angle variability helps explain the response variability within each cluster but not necessarily the cluster identity. It was observed that the posterior predictive distributions obtained do not envelope the experimental data perfectly. This result is likely due to the limited parameter range explored. It is possible that expanding the parameter ranges would help capture the dataset extremes better. However, this expansion would require retraining of the model and additional high-fidelity training data, which is too expensive to obtain.

A multifidelity surrogate model was developed to efficiently compute the nonlinear frequency backbone curves for the TRC benchmark structure. The multifidelity model combined a high-fidelity and low-fidelity finite element model to produce fast but accurate predictions of the backbone curves. The multifidelity model was constructed using an ensemble of neural networks including a component with Bayesian layers to quantify the network's prediction uncertainty. In addition, the network architecture included explicit dependence of the correction functions on the low-fidelity network output while preserving the additive skip connection, similar to the original nonlinear autoregressive formulation. The resulting surrogate model was used to explore the sensitivity of the backbone curves to the input parameters as well as performing Bayesian inference to gain insight into the trends observed in the experimental data.

This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns all right, title and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan². This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.