A Study of Whole Joint Model Calibration Using Quasi-Static Modal Analysis

Mechanical assemblies rely on joining technologies to produce complex structures from readily manufacturable piece parts. One of the most common joining technologies is the bolted joint due to its commercial availability, ease of assembly and disassembly, low cost, and load carrying ability. Other types of joints include press fits, welds, tape joints, dovetail joints, and finger joints. All these joints introduce mechanical interfaces that lead to contact and frictional forces to carry loads between parts. Frictional contact is localized to the region of the interface, but it can introduce nonlinear behavior in the dynamic response of the global structure. These interfaces contribute to the overall stiffness and load carrying capability of the structure and can also account for up to 90% of the energy dissipation [1].

The influence of joints and mechanical interfaces on structural dynamic responses has been studied by a vast community of researchers partly due to the importance of understanding the damping mechanisms in high-consequence systems such as turbine engines [2] or aerospace structures [3,4]. Several modeling approaches are documented in the literature that capture the nonlinear physics associated with frictional contact, for example, see the review papers by Gaul and Nitsche [5] and Bograd et al. [6]. One strategy to model friction in joints is to explicitly model the contacting surfaces using the finite element method, and use frictional contact elements to capture the local interactions between element nodes and/or faces [7]. Most commercial, nonlinear finite element packages utilize Coulomb friction elements to capture in-plane tractions tangent to the contacting surfaces whose normal pressures vary with time. Several researchers have used this approach to model joints and steady-state responses to harmonic excitations via the harmonic balance method [8–13]. The harmonic response idealizes the general loading a structure may experience during its service lifetime, and thus a more practical method would allow for general transient response prediction to loads such as impact, random, or sine loading. Due to the excessive cost of direct numerical time integration with iterative contact algorithms, solutions of transient vibrations with a fully discretized interface become intractable.

Many researchers have explored other opportunities to model frictional effects in mechanical joints using low-order, or simplified, representations of the contacting surfaces. The review paper by Mathis et al. [14] provides a historical overview of the various rate independent and dependent constitutive models used to capture nonlinear dissipation in joints. Many of these constitutive elements can be implemented into finite element models (FEMs) using a whole joint approach [15] to relate forces between two contacting surfaces. A whole joint is created by constraining the nodes on a contact surface to a single, virtual node using either averaging or rigid multipoint constraints (MPCs). The constitutive elements are added between the relative degrees-of-freedom (DOFs) of the virtual node pairs to account for stiffness or damping localized to the jointed interface. Each whole joint has six independent DOF and can be populated with linear elements such as elastic springs and dashpots or nonlinear elements derived to capture the hysteretic behavior. The computational cost of the whole joint modeling approach can be significantly reduced by creating a reduced order structural dynamics model with Hurty/Craig-Bampton (HCB) super elements [16,17] that preserve the physical DOFs localized at the joint. This leads to a low-order model form that can be efficiently integrated in time to predict transient solutions to general load cases, while approximately accounting for important nonlinear physics in the joint.

One of the key and practical challenges associated with whole joint modeling is determining the parameters of the constitutive elements used to capture the joint physics. The unknown joint parameters often require some calibration or tuning based on measured data from test hardware, which is not always available during the engineering design phase. Other methods rely on simulated data from a high-fidelity nonlinear finite element model but could be costly to produce the necessary reference data. Several researchers have presented optimization studies to calibrate nonlinear joint parameters in a structural model. Lacayo and Allen [18] presented a methodology to calibrate physical Iwan elements in an HCB model to experimental data of the Brake-Reuß beam using quasi-static modal analysis (QSMA) and postprocessing methods to obtain the instantaneous frequency and damping from ring-down experiments, e.g., see Ref. [19] for an overview of these techniques. Other studies of parameter identification of whole joint models include the work by Charalampakis and Dimou [20] and Charalampakis and Koumousis [21] who developed a method to determine Bouc–Wen element parameters using different global optimization schemes. The research by Wang et al. [22] used high-fidelity finite element simulations of oscillatory and monotonic loading of a single bolted joint to curve-fit the parameters of an extended Iwan element capable of capturing nonzero stiffness in macroslip. Other studies have proposed the use of direct time-integrated responses, such as the work of Wang et al. [23], where a joint model updating scheme was developed using analytical mode decomposition to extract the instantaneous modal characteristics to serve as the reference data for calibration. Alternatively, Oldfield et al. [24] fit the parameters of both Bouc–Wen and parallel Jenkins element models to match the joint hysteresis predicted from a detailed finite element model at various load levels. Several of these existing studies have utilized either static or dynamic simulation data to provide reference solutions to calibrate the joint parameters, but none to date have been capable of capturing the dynamic effects of the modal responses using static analyses only.

This research presents a framework that calibrates an HCB reduced order model (ROM) with whole joints and four-parameter Iwan elements [25] to capture the global nonlinear response from simulated reference data of a high-fidelity model of a bolted assembly. The method calibrates the joint parameters such that the nonlinear modal characteristics of the structure are preserved, i.e., the amplitude-dependent natural frequencies and damping ratios for all modes of interest. This is achieved through the use of the QSMA approach originally developed by Festjens et al. [26] and later extended to ROMs with Iwan elements in Ref. [18] and full-order finite element models in Ref. [27]. The QSMA approach applies a force proportional to the linearized mode shape of interest to the nonlinear, preloaded structural model and uses a quasi-static solver to obtain the response over a range of force levels. From the response, the modal force–displacement and modal hysteresis curves are obtained to provide an estimate of the amplitude-dependent natural frequency and damping ratio as a function of the modal response amplitude. A global optimization algorithm with a multiobjective function is utilized to determine the joint parameters that minimize the mean absolute difference in amplitude-dependent natural frequencies and damping ratios between the reference QSMA results computed from the full-order model and the QSMA response from the ROM with Iwan elements. This approach has the advantage that the simulated data does not require expensive dynamic simulations and can be readily obtained from both high-fidelity and reduced order models. In addition, the QSMA results capture the effective loading on the joint due to the modal deflection shape with only static analyses and reveal how the presence of joints influence the modal characteristics of interest. Performing QSMA on the full-order model helps determine which modes need to be calibrated with nonlinear elements in the whole joint model.

The sections of this paper are as follows: Section 2 presents the whole joint modeling approach combined with the HCB model reduction theory. This section also presents the numerical optimization scheme used to calibrate the four-parameter Iwan elements to capture microslip behavior. The results of the calibration of the C-Beam assembly with two bolted interfaces are investigated in Sec. 3, while validation with direct time-integrated solutions are provided in Sec. 4. Conclusions are provided in Sec. 5.

The model calibration process begins with the development of a detailed finite element model of the jointed structure. The process described here is applicable to most finite element software packages that accommodate nonlinear static solutions with frictional contact, modal analysis linearized about a preloaded state, and Hurty/Craig-Bampton substructure generation. A flowchart describing the general steps of the procedure is shown in Fig. 1. Step 1 in the process uses quasi-static modal analysis to extract the amplitude-dependent natural frequencies and damping ratios for the mode(s) of interest from the full-order model with detailed frictional contact at the interfaces. Step 2 reduces the finite element model using the Hurty/Craig-Bampton method by creating whole joints based on the surfaces determined to be in contact from the preload analysis in step 1. The HCB model with whole joints is populated with constitutive elements to describe the joint behavior, for example, with either linear springs or four-parameter Iwan elements. The final step involves the calibration of the parameters in the selected whole joint model to capture the nonlinear modal characteristics (i.e., amplitude-dependent natural frequency and damping ratio) predicted with QSMA in step 1. A global optimization routine with a multiobjective function explores the parameter space to discover the optimal whole joint model parameters to capture the correct physics in all the nonlinear modes of interest. The following subsections provide details of the theoretical aspects of each step outlined in Fig. 1.

Quasi-static modal analysis is a general procedure that has been developed for either full-order finite element models with frictional contact created within a finite element package [27] or on ROMs with whole joint representations that utilize four-parameter Iwan elements [18]. The method assumes that the physical model for the joint or mechanical interface is described with elements that are rate independent and satisfy Masing's hypothesis. Both the four-parameter Iwan element and Coulomb friction elements used throughout this work satisfy these conditions for the mode shapes studied since they predominantly load the joint in shear and do not significantly alter the normal loads (i.e., in the case of Coulomb frictional with variable normal loads). The reader is referred to Refs. [18,27] for complete details of the QSMA approach; however, a brief overview is provided here for completeness.

This quasi-static equation is computed over a range of modal force amplitudes by incrementally increasing the scale of α. The equilibrium response computed over this range is dependent on the force level and is defined as $x(α)$⁠.

The fact that Masing's rule can be used to construct the hysteresis curve provides a significant computational savings since only the backbone curve of the modal force–displacement needs to be simulated. The alternative to this is to directly calculate each hysteresis loop at the modal response levels of interest, which would add an excessive computational burden particularly when applying this method to large-scale finite element models. The resulting nonlinear natural frequencies and damping ratios from Eqs. (7) and (9) characterize the dynamics of the jointed system and provide insight into how the modes are influenced by frictional contact in the joints. In this research, these solutions serve as the reference solution to calibrate ROMs created with whole joints, as shown in Sec. 3.2. A similar procedure to perform QSMA on the ROM with four-parameter Iwan elements can be found in Ref. [18]. The QSMA algorithm is orders of magnitude faster with the ROMs compared with the full-order model (e.g., on the order of seconds for a ROM with 50DOF and two nonlinear elements) and hence has the advantage of allowing full exploration of the parameter space with global optimization schemes. Additionally, if the ROM captures the nonlinear characteristics with respect to the vibration modes of interest, then it is expected that these ROMs should generally provide accurate predictions of forced transient responses at a fraction of the computational cost of the full-order model. Section 4 examines the accuracy and efficiency of the calibrated ROMs in detail.

A whole joint approach with structural dynamic models approximates the frictional contact in the interfaces with MPCs and one-dimensional constitutive elements describing the joint forces. Whole joint models have the advantage of reducing time and length scales in the finite element model hence speeding up simulation times for forced transient response predictions over long periods. An example of a whole joint model within a bolted joint is shown in Fig. 2 (note the view is exploded for visualization). The nodes determined to be in contact from a nonlinear preload analysis are condensed down to a single virtual node using so-called spidering approaches to create MPCs with either rigid or averaging type elements. The single, vertical bar in Fig. 2 represents a zero-length element connecting coincident virtual nodes with six linearly independent DOFs that can be represented with various types of constitutive elements that describe the joint forces.

A reasonable assumption for jointed structures is that, away from the joint, the model behaves linearly such that deformations are small, and the material strains within its linear elastic range. With these assumptions, a ROM with whole joints can be formulated using component mode synthesis techniques to preserve the physical DOF at the whole joint interface to allow for either linear or nonlinear constitutive elements to describe the joint forces. In this work, the HCB method [16,17] is used to reduce the linear, interior (i-set) DOF with a truncated set of fixed-interface modes, and statically reduce the virtual (b-set) DOF at the whole joint connections with static constraint modes. The theory is briefly reviewed here.

Several constitutive elements exist to capture the physics in the mechanical interface. This research uses the four-parameter Iwan element developed by Segalman [25] to capture frictional microslip which causes nonlinear energy dissipation and stiffness loss in the joint. The Iwan element is a one-dimensional constitutive model that is a parallel arrangement of series spring-and-slider elements with a power-law distribution describing the population density of slider strengths. In this research, the four-parameter Iwan element is used to model the shear forces within the whole joint ROM; the theory is briefly reviewed here but the interested reader should refer to Ref. [25] for full details.

The dimensionless values, x(t, ϕ), are the displacements of the Jenkins elements with a slip displacement, ϕ. The population density of Jenkins elements is ρ(ϕ), and u(t) and f_NL are the joint displacement and internal nonlinear force, respectively. The Iwan element in Eq. (12), assuming a power-law distribution, is completely defined by the parameters provided in Table 1. As discussed in Ref. [25], these parameters are preferred since they are measurable quantities. The slip force has physical meaning since it can be approximated from the normal force calculations with an assumed Coulomb friction coefficient. The joint stiffness is approximated from modal vibration tests performed at low excitation levels to determine the joint stiffness when no slip occurs. The other two coefficients, χ and β, are more elusive and must be calibrated to specific experiments or simulations, such as the method proposed in this work. For this work, the reference data comes from a high-fidelity model with the assumption that the model parameters are known. Additional calibration of the full-order model parameters may be necessary prior to proceeding with the procedure outlined in this paper but is considered beyond the scope of the present study.

A multiobjective optimization scheme is used to determine the four Iwan element parameters in Table 1 within each of the whole joints in the ROM. These parameters are calibrated by matching the amplitude-dependent natural frequency and damping ratio obtained from the QSMA algorithm applied to the ROM, which can be sampled quickly due to the low-order of the system. The results from the full-order finite element model, as discussed in Sec. 2.1, serve as the reference data. To this end, the non-dominated sorting genetic algorithm II (NSGA-II) from the python package DEAP is used to perform the numerical optimization. NSGA-II has been shown to perform better than other traditional multiobjective optimization techniques [29]. Metaheuristic algorithms such as genetic algorithms (GA) are desirable in this application because the Iwan element parameter design space was observed to contain multiple local minima in previous work [30] and GA performs better in these multimodal design spaces than gradient-based methods [31]. The disadvantage of evolutionary algorithms is that they typically require many more iterations than their gradient-based counterparts, which is why using an efficient ROM for optimization is important. Additionally, GA avoids the computational expense of traditional Monte Carlo sampling methods and reduces the number of samples needed to achieve an optimal solution. A Monte Carlo sampling method explores the design space at random; however, a genetic algorithm is focused on populating samples near the Pareto-front. The Pareto-front is the subset of non-dominated solutions, i.e., solutions for which there is no better performing alternative for at least one of the objectives. NSGA-II achieves this goal by taking the best candidates at each iteration and recombining them; this is similar to how biological evolution works. Like evolution, the algorithm introduces random mutations at each iteration to diversify the population of solutions and avoid being trapped in local minima.

The amplitude-dependent natural frequencies, $ω~r(αi,p)$⁠, and damping ratios, $ζ~r(αi,p)$⁠, are computed from the whole joint reduced order model using QSMA for an arbitrary set of whole joint parameters, $p$⁠. The modal data from the full-order model, ω_r(α_i) and ζ_r(α_i), are computed offline. The advantage to using QSMA in the reduced model is the cost efficiency to evaluate one instance of the objective function. Other approaches that rely on dynamic, time-integrated solutions are costlier to simulate, and more difficult to automate with signal processing tools within an optimization scheme [32].

The following subsections describe the details of the models and results used to evaluate the proposed calibration approach to determine the whole joint Iwan parameters within a nonlinear finite element model. The finite element model is introduced in Sec. 3.1 and a discussion of the QSMA results obtained from the full-order model are provided in Sec. 3.2. The results obtained by calibrating the Iwan parameters to match the nonlinear modal characteristics are presented in Sec. 3.3.

The full-order finite element model of the C-Beam bolted assembly was meshed in CUBIT [33]. A mesh convergence study (coarse, medium, and fine) was performed to determine the mesh refinement needed for accurate predictions of the nonlinear preload and linearized vibration modes. It was determined that there was less than 1% difference between the mode frequencies of the medium and fine mesh. The difference in contact forces between the three refinement levels was less than 1%, so the medium mesh was selected for computational efficiency. The final (i.e., medium) mesh had 579,092 eight-noded reduced-integration hexahedral elements and is shown in Fig. 3. The model had a uniform mesh along most of the length of the beam except near the contact area, where a finer mesh density was used to better resolve the contact in the bolted region, as illustrated in Fig. 3. The input location for the transient forces in Sec. 4 is also annotated in this figure.

The contact interfaces were modeled with a Coulomb friction coefficient of 0.2 within all contacting surfaces, which include the bolt head/nut to washer, washer to beam, and beam to beam interfaces. The contact is localized to the bolted regions on each end of the beam. The individual C-Beam and the bolt piece parts were assigned elastic properties corresponding to 440C stainless steel with a Young's modulus of 200,000 MPa, Poisson's ratio of 0.283, and mass density of 7612 kg/m³. The bolt preload was applied by prescribing an axial artificial strain to the shank of the bolts corresponding to a 28,024 N axial preload force. Four nodes (two nodes on either bolted interface) were fixed to prevent rigid body motion during the simulations. All nonlinear quasi-static preload analyses were performed using the implicit solver with Sierra/SM finite element package [34]. The target relative residual for the implicit solver was on the order of 10⁻⁶, while the acceptable relative residual was an order of magnitude higher.

Following the bolt preload simulation for the C-Beam assembly, the finite element model was linearized about the preloaded state to compute the preloaded modes of the assembly. The computational modal analysis was performed using Sierra/SD [35]. With the computed modes, QSMA was performed by applying a body force corresponding to $fr(α)=Mϕrα$⁠, where ϕ_r is the linearized mode shape of interest and α was varied incrementally as the joint underwent microslip and eventually macroslip. The quasi-static force analysis restarted from the preload equilibrium state using the Sierra/SM implicit solver. For the modal force simulations, the target relative residual was set to 10⁻⁴ while the acceptable relative residual was set to one order of magnitude higher. Several ranges of acceptable residuals were explored, and it was found that smaller tolerances did not make a significant difference in the quality of the results but added significant computational burden to the solution. However, after a certain relative residual threshold (10⁻³) the solution began to deteriorate. The heuristics of choosing solver settings has been found to be problem-dependent and required considerable manual tuning before a converged model was found [27]. The process of selecting adequate solver settings for this problem was largely an iterative process that sought to balance runtime, accuracy, and numerical convergence. The solver relative residual was first set to the largest number that resulted in a converged solution. From that point, the relative residual was decreased gradually until sufficient convergence of the quantities of interest (i.e., the displacements) was observed. Convergence was deemed sufficient when there was less than 1% change in the peak displacements.

QSMA was performed for both mode 1 (223.4 Hz) and mode 8 (1142.3 Hz) since these mode shape deformations activated the jointed interface by imparting shear stresses in the x-axis. The deformation shapes of these modes shown in Fig. 4 reveal the strains induced in the bolted region of the assembly. These modes are also symmetric about the mid-plane and have the potential to couple through the nonlinearity in the contact regions. The whole joint models developed to capture the joint physics utilized four-parameter Iwan elements in the relative x-displacement for each of the interfaces. In this scenario, the Iwan elements needed to be tuned to capture the nonlinear behavior of the two modes simultaneously.

The amplitude-dependent natural frequencies and damping ratios for the full-order finite element model are shown in Fig. 5 as a function of the modal displacement magnitude, q₁(α) and q₈(α). The curves for positive and negative loading were obtained to confirm Masing's hypothesis, indicating that the forward and reverse backbone curves produced the same quantitative results. As illustrated in the QSMA analysis, the joint remained stuck for both modes at low displacement amplitudes and the estimated natural frequencies converged to the natural frequencies predicted from the modal analysis linearized about the preloaded state. As the modal displacement amplitude increased, the resonant frequency decreased significantly since the shear loads within the mechanical interface started to overcome the in-plane frictional forces. The joint began to slip in low pressure regions which led to a decrease in joint stiffness. Similarly, the damping increased as the modal displacement amplitude increased due to the energy dissipation from frictional losses in the interface. Each QSMA curve shows a distinct bend or knee in the data, occurring around 3 · 10⁻⁵ and $8⋅10−6m−kg$ for modes 1 and 8, respectively, which occurred around the transition from microslip to macroslip. These curves serve as the reference data to the whole joint model calibration performed in Sec. 3.3. The ranges shown here correspond to the ranges used in the optimization scheme, i.e., 3.5 · 10⁻⁶ to $6.9⋅10−5m−kg$ and 2.3·10⁻⁶ to $2.3⋅10−5m−kg$ for modes 1 and 8, respectively. These ranges were selected to capture the power-law region associated with microslip as well as the transition region.

The amplitude-dependent damping ratios in Fig. 5 were estimated from the backbone curve using Masing's rule to construct the cyclic hysteresis curve. This hypothesis was investigated to confirm the validity of the assumption by directly calculating the modal hysteresis curve with the full-order finite element model. This was performed for mode 1 by loading/unloading the structure with the modal force in Eq. (5) over three full cycles. The peak loading amplitude, α_max, was chosen to be $38,618N/kg$⁠, resulting in a peak modal displacement amplitude of $4.8⋅10−5m−kg$⁠. This level was chosen since it produced a response well into the macroslip regime as evidenced by the plots in Fig. 5. The overlay of the direct hysteresis and the constructed hysteresis from the backbone using Masing's rule are shown in Fig. 6. The direct hysteresis curve shows to be in a steady-state condition as the repeated cycles nearly overlay, which agree well with the Masing's hysteresis constructed from the initially loaded backbone curve. Based on the area within the modal hysteresis curves, the modal dissipated energy calculated using Eq. (8) was 2.15 N − m using Masing's rule and 2.23 N − m using the direct approach. The constructed hysteresis curve, in this case, produced an error of 3.9% for the damping ratio estimate. The use of Masing's rule allows only the backbone of the load–displacement curve to be calculated, offering significant computational savings without the need to directly compute each modal hysteresis curve explicitly. The assumption of Masing's rule needs to be critically evaluated for QSMA applied to finite element models with frictional contact nonlinearities with varying normal loads, as the approximation may breakdown for modes that load the joint in other patterns (e.g., bending or separating). The approximation works well for joints predominately loaded in shear, but an in-depth study of these limitations is beyond the scope of this paper.

Based on the bolt preload simulation, an HCB reduced order model was generated by creating MPC spiders within the bolted interface to reduce the kinematics down to a single, virtual node. The MPC was created to take all nodes determined to be in contact and used rigid bar elements to constrain the nodes on the contact surface down to a single point. The area of the MPC region is shown earlier in Fig. 2, which resulted in two whole joints at each end of the beam assembly. The two interfaces were represented with a four-parameter Iwan element in the relative x-displacement direction, and linear spring elements in the remaining directions (y- and z-displacements, x-, y-, and z-rotations) that were manually tuned to match the linear mode natural frequencies outside of the Iwan element optimization procedure. The HCB ROM was generated with a total of 25 fixed-interface modes, spanning frequencies up to 6200 Hz, and 27 static constraint modes, resulting in a model size of 52DOF. There was a total of 6DOF for each of the four virtual nodes in the whole joint, and 3DOF for the input location of the applied force at the midpoint of the top beam (used for transient response validation studies in Sec. 4). The practice of including the input location DOF in the HCB model may lead to slower convergence of the HCB ROM with number of modes because the input DOF is restrained, resulting in mode shapes that may be different than the system modes [36]. For this reason, convergence was ensured by retaining enough fixed-interface modes.

Due to the symmetry of the system, the size of parameter set within the optimization routine was reduced by assuming a priori that the whole joint parameters were identical on each end, hence resulting in a four-dimensional parametric space. It should be noted that the linear spring elements in the other whole joint directions were calibrated to match other linear mode frequencies of interest (e.g., mode 2 imparts a bending moment about the z-axis, thus calibration of the z-rotational spring was required to capture the natural frequency of the ROM). The calibration of the whole joint elements in the HCB ROM was performed by running the global optimization algorithm using only the nonlinear data as described in Sec. 2.3.

The reference data for the nonlinear calibration were the two sets of amplitude-dependent frequency and damping ratio curves corresponding to mode 1 and mode 8 obtained with QSMA on the full-order model; see data in Fig. 5. NSGA-II minimized the objective function defined in Eq. (14), which consisted of four different objectives to be minimized simultaneously (frequency and damping for the modes 1 and 8). Since a multiobjective scheme was used, there was no need to weigh the objectives differently. Instead, NSGA-II found the Pareto-front corresponding to the subset of non-dominated solutions. The population size was 50 and the algorithm was run for 150 generations. Each generation had 20% probability of crossover and 50% mutation. The results of the optimization search are shown in Fig. 7, where Figs. 7(a) and 7(b) distinguish between mode 1 and mode 8, respectively, and the abscissa and ordinate correspond to the frequency and damping error, respectively. The entire population of the genetic algorithm is shown in black while the samples along the Pareto-front are shown in magenta; the green samples correspond to the best fits that are evaluated further. The Pareto-front resides in a four-dimensional space, so the figures represent a projection onto a two-dimensional space for ease of visualization. The Pareto-front produced by NSGA-II uses the Euclidean distance to compute a cuboid surrounding the neighbors around each point along the front and defines the crowding distance as the average of the cuboid sides [29]. Diversity along the Pareto-front is encouraged by selecting individuals with large crowding distance across the objective space. In addition to the NSGA-II results, a Monte Carlo sampling scheme was used to compare the efficiency of both techniques in reproducing the Pareto-front using the same number of sample evaluations. The samples from the Monte Carlo simulation are shown in gray in the same figure. As can be seen by the data, the Monte Carlo solutions were uniformly spread over wide regions, whereas the NSGA-II solutions were concentrated around the Pareto-front. The Monte Carlo approach failed to populate the Pareto-front with the same number of iterations required for NSGA-II, thus making the latter approach more efficient in terms of sampling the parameter space of the ROM. This challenge worsens when more parameters are included due to the well-known curse of dimensionality.

The five optimal solutions in Fig. 7 (shown in green) are evaluated in more detail. Each of these solutions all reside on the Pareto-front and correspond to the four optimal solutions for each individual objective (i.e., mode 1 frequency and damping, mode 8 frequency and damping), as well as the best overall case that minimized the sum of the norms of each objective function. The resulting amplitude-dependent frequency and damping curves for the five optimal cases are plotted in Fig. 8 alongside the QSMA data obtained from the full-order model. The curves illustrate the tradeoffs that occurred when trying to minimize multiple objectives. For example, the “M8 best damping” solution provided the best agreement with the full-order model QSMA data for mode 8 damping ratio, but overestimated the damping ratio, particularly in the microslip power-law region, for mode 1. The “M1 best frequency” solution provided better agreement over the linear portion of the curve than over the transition region. This behavior may be due to the errors in the linear range being larger than the range within the transition region, ultimately driving the solution to emphasize correlation in the linear portion. The “best overall” solution appears to perform quite well at capturing each nonlinear mode characteristics. The only notable error with this solution is the damping ratio for mode 8, particularly in the lower modal displacement amplitude region where microslip occurred. Regardless, the NSGA-II global optimization algorithm was able to identify a parameter set that effectively produced a model that captured two nonlinear vibration modes simultaneously. This is mentioned to be a significant challenge by previous studies on whole joint model calibration [18].

The Iwan parameters corresponding to the five optimal cases are listed in Table 2. There were some clear differences in the parameter values that were identified to capture the various objective functions. For example, the joint stiffness, K_T, spanned two orders of magnitude for the “M1 best damping” and “M1 best frequency.” The stiffness term could have been fixed and determined based on the linear frequency of modes 1 and 8; however, it was discovered that the global optimizer was able to better capture the full range of frequency and damping by allowing this value to vary. Similarly, in theory the slip force, F_s, could be determined a priori from the normal bolt preload (28,024 N) and the Coulomb friction coefficient (μ = 0.2), which produced an estimated slip force of 5605 N. This value was also allowed to deviate from this theoretical value for the same reason as the joint stiffness. The biggest discrepancy in slip force occurred for the “M1 best frequency” with a slip force 14.6% lower than the theoretical value. The other two parameters, χ and β, have less physical meaning to the joint model and hence the optimization approach was needed to properly tune these values. Note that the optimizer set χ to nearly −1.0 for most of the optimal solutions, which is the known limit of the χ parameter. Likewise, the β parameter was set to a very small number, approaching the limit of 0.0 for β. In other words, the optimizer set these two parameters to values close to their allowed limits to match the reference data. This result indicates that the Iwan element may be reaching its limit to represent the bolted joint under consideration, and that other joint models should be explored in future research. An in-depth analysis of the physical implications of the optimized Iwan parameter values is beyond the scope of this paper.

To validate the calibration approach, the HCB reduced order models with calibrated Iwan parameters were used to predict the transient forced response to an impulse load in Sec. 4.1 and a white-noise random excitation in Sec. 4.2. The accuracy of the simulations was assessed, and the overall cost analysis is provided in Sec. 4.3. For the transient analyses presented here, the only source of dissipation in the models is from friction effects and numerical damping. No other viscous damping models were used to capture other sources of damping in the structure.

The transient response of the full-order and the reduced order model were evaluated for a haversine impulse load. The full-order model response was evaluated using the nonlinear explicit solver in Sierra/SM. Five different HCB models were run with the optimal Iwan parameters listed in Table 2. The plots in Fig. 9 show the transient time histories and Fourier spectra for each input amplitude considered, namely, a 1000 Hz haversine pulse with a peak force of 222, 890, and 2224 N. This range was carefully selected to excite the nonlinear model near the linear regime where no slip occurs in the joint, and into the nonlinear regime to the edge of the microslip regime near the transition region. The input force was applied in the +y-direction at the midpoint node of the top beam.

The data in Table 3 summarizes the TRAC comparison of the drive-point displacements for the various models. The TRAC metric measures the level of correlation between two signals where a value of 1.0 indicates perfect correlation and 0.0 indicates the signals are completely orthogonal. One can observe that the “best damping” models have the worst TRAC correlation, and the “best frequency” and “best overall” all have better correlation with values ranging from 0.76 to 0.88. One should be careful drawing conclusions from such correlation metrics for time-domain analysis since it is known that the TRAC penalizes phase error more so than amplitude error, and hence these results do not provide an absolute error measurement between the ROMs. The results presented later in Fig. 13 visualize the types of errors encountered for the M1 “best frequency” and “best damping” transient responses. The quality of the ROM should be determined by the needs of the modeling efforts.

The time histories of the full-order model and “best overall” HCB ROM are shown in Fig. 10 for visual comparison at all three considered excitation levels. As illustrated, the calibrated “best overall” HCB model matches the full-order response closely, both in phase and magnitude, over a 100 ms period. The best agreement occurs at the higher forcing amplitudes, where the joint dissipates energy due to microslip in the Iwan element. The lower 222 N force case is the least accurate but still captures the transients and peak amplitudes well while responding in the linear regime. The TRAC for the 222 N and 2224 N case produced the same correlation of 0.76, but the visual comparison would suggest that the 2224 N case agrees better with the full FEM, thus confirming the challenges with time history correlation metrics. Note that the time histories evaluated here excite several modes of the structure, and hence lead to complex multimodal behavior within a nonlinear system.

Each term in the summation was used to determine the relative importance of each mode to the measurement point of interest. It was determined that modes 1, 2, and 8 were the dominant modes participating in the drive-point responses in Fig. 10. The mode shape for mode 2 shown in Fig. 11 and the mode shapes for modes 1 and 8 in Fig. 4 reveal that the peak displacement of the modes occurred at the drive-point, or midpoint of the top beam, thus confirming the expectation that these modes contribute significantly to the applied force. The peak modal displacement amplitude for modes 1 and 8 over the three excitation levels is listed in Table 4, along with their normalized response with respect to the peak force amplitude. The peak modal amplitudes can be directly related to the QSMA plots in Fig. 8 to show what level of nonlinear response the given force imparts on the structure. For example, the 2224 N excites the mode 1 coordinate to a peak amplitude of $2.85⋅10−5m−kg$ which occurs just prior to the bend of the QSMA curve transition between micro- and macroslip.

The drive-point modal contributions for modes 1, 2, and 8 are plotted in Fig. 12 for both the full-order model and the “best overall” ROM at the 222 N and 2224 N input levels. The HCB model was able to capture both the damping and the frequency of all the significant modal contributions to the drive-point response. The mode 1 contribution revealed that the decay envelope was significantly steeper for the 2224 N case versus the 222 N case, confirming that the nonlinear, amplitude-dependent damping was accurately captured with the calibrated Iwan elements. Furthermore, the mode 8 damping and phase both agreed well with the full-order model characteristics which confirmed that mode 8 also dissipated more energy at higher force levels. Interestingly, the mode 2 response amplitude remained flat over the entire period of response for both load cases. This mode imparted a bending moment on the joint, as seen in the deformation shape in Fig. 11, and did not induce in-plane shear stresses known to produce frictional losses in lap joints. This mode, as expected, did not have a nonlinear amplitude dependence of the damping ratio based on the loads imparted on the joint, and as a result did not draw energy from the system like modes 1 and 8.

Much of the focus in this subsection has been on the “best overall” HCB ROM and how it compared with the transient response predicted by the full-order model. The other optimal ROMs that were based on the best fit objective for either the frequency or damping are further evaluated here. Specifically, the plots in Fig. 13 show the drive-point contribution of mode 1 for the “M1 best damping” and “M1 best frequency” ROMs in comparison with the full FEM for the two extreme load cases (i.e., 222 N and 2224 N). These plots reveal how the optimizer and the resulting response predictions performed when predicting the ring-down response. As expected, the “M1 best damping” captured the envelope of decay for mode 1 but sacrificed the accuracy of the phase between the two signals. In contrast, the “M1 best frequency” performed well at capturing the correct phase for the low and high force levels but slightly overpredicted the decay envelope amplitude. These results highlight both the challenge and flexibility with optimizing the Iwan parameters and can produce different results depending on what is determined to be the most critical feature of the model. The time histories also reveal the challenge associated with time-domain correlation metrics such as TRAC in quantifying model errors. The “best damping” error associated with phase was clearly penalized far more than the “best frequency” amplitude error, as reported in Table 3. The quality of the ROM will depend on the type of input and metrics of interest for each specific analysis.

Next the full-order model and “best overall” ROM were evaluated via comparison of drive-point responses to a white-noise broadband excitation with flat-line input of 19.8 N²/Hz (1.0 lbf²/Hz) between 10 and 1250 Hz to excite up to the eighth mode. The plots in Fig. 14 show the time history and power spectral density (PSD) of the input with a root-mean-square (RMS) force level of 155.2 N. This force time history was scaled to effectively increase or decrease the input RMS to excite a range of linear and nonlinear responses in the models. All simulations were performed in the time domain over a 1 s period, and the PSDs were all computed with the same method and sampling rate to fairly compare the different models and force levels. The simulations of the full-order model were performed using the Sierra/SM explicit solver. It is possible that small errors in the explicit routine are accumulated over time, especially for long-duration simulations such as those for random vibration and could have an effect on the accuracy of the full-order model results.

The “best overall” ROM was simulated over a range of RMS input force levels to evaluate the influence of nonlinearity on the modes and observe the evolution of the peaks in the PSD at increasing amplitudes. The drive-point response PSDs are shown in Fig. 15 for different RMS force levels listed in the legend. Note the relative differences between the mode 1 and 2 peaks: at low amplitudes, the mode 1 peak was higher than mode 2. However, as the force level increased, mode 2 became the higher peak due to increased damping in mode 1, revealing the effect of energy dissipation in the joint captured by the calibrated Iwan elements. The nonlinear effect was most noticeable for the 1552 N RMS force level as the mode 1 and mode 8 peaks broadened and shifted significantly. This RMS level was the highest forcing level explored with the HCB ROM and was beyond the range of validity for the ROM. Much like the observations in Sec. 4.1, the mode 2 response did not change (i.e., shift or broaden) with increasing force level.

Two RMS force levels were chosen to evaluate the performance of the ROM via direct comparison with the predictions from the full-order model. The random input forces with RMS values of 310.9 N and 3.109 N were applied to the full-order model and integrated using explicit time integration in Sierra/SM. The comparison of the time domain response of the full-order model versus the ROM is shown in Fig. 16 for the higher 310.9 N RMS case. Note that only snap shots of the data are shown to better visualize the correlation between signals. For this higher force level, the ROM time history response matched the full-order model response closely at the start of the simulation, as shown in Fig. 16(a) spanning a period of 0.1 s. The response agreed well up until approximately 0.5 s, where it was observed that the ROM response started to deviate from the reference response, as shown in Fig. 16(b). Interestingly, the phase between the two signals was still in agreement beyond this point, but the amplitudes were quite different. Similar behavior was also observed for the lower amplitude case of 3.109 N RMS but was not shown for brevity.

The drive-point response PSDs are shown in Fig. 17 for both the lower and higher input levels considered in this study. For the PSD of the 310.9 N RMS case, the ROM response matched the peak amplitudes corresponding to modes 1 and 8 response. However, the mode 2 peak was slightly overpredicted by the ROM in comparison with the full-order model. An interesting observation was that the ROM also severely underpredicted the peak just beyond the mode 8 peak, occurring around 1480 Hz. This frequency is beyond the excitation range, so it would not be expected that this mode be strongly excited by the input but rather through modal coupling. Further investigation of the mode shapes revealed that mode 12, with a resonant frequency of 1494 Hz, corresponded to the third out-of-phase bending of the beam assembly. This symmetric shape correlated with mode 2, thus allowing for the exchange of vibrational energy via modal coupling, explaining the discrepancy of the peak around 1480 Hz. Since the ROM only captured the joint physics of these out-of-phase bending modes with a linear spring, it would not be expected to explicitly capture modal coupling. This phenomenon can be confirmed by observing the drive-point PSD at the lower force level, where the mode 12 peak was in much better agreement with the full-order model.

To further investigate the observed discrepancies in the higher force level, a spectrogram was computed using the short-time Fourier transform (STFT), as shown in Fig. 18. The spectrogram revealed that the ROM was not predicting vibrational energy beyond mode 8, whereas the full-order model predicted significant response amplitude in the vicinity of 1500 Hz, and beyond. This observation agreed with the discrepancy in the peaks shown in the drive-point PSDs. The response content beyond the input frequency of 1250 Hz can be explained through nonlinear phenomena known as mode coupling, indicating that mode 2 was exchanging energy with mode 12 when responding in the nonlinear response range. The ROM was developed to capture the energy exchange between modes 1 and 8, which both fall within the excitation bandwidth.

The computational burden of the direct transient solutions is shown in Table 5 to evaluate the efficiency of each model evaluated in Sec. 4. The most expensive model is the full-order model, which was computed on 112 processors on a high-performance computer (HPC) using the Sierra finite element packages. It took approximately 48 h to solve the full-order model using the Sierra/SM explicit solver to obtain 100 ms of response data to the haversine impulse. Direct transient random vibration analysis was even more expensive and took 20 days to obtain 1.0 s of data. The ROMs presented here have an initial offline cost so that their online cost can be drastically reduced. The offline cost of QSMA on the full-order model included the initial preload, the linearized modal analysis about the preloaded configuration, and the QSMA modal force computation, adding up to a total of 30 h of simulation time of the full-order model with 112 processors. Calibration of the Iwan parameters of the HCB model was not a trivial effort and took 5 h on a single processor. It should be noted that each QSMA run only took 1–2 s with the reduced HCB model, however the global optimization searched on the order of tens of thousands of parameter sets to reach the optimal solutions. Once the HCB model had been calibrated, the transient simulations were able to be obtained extremely quickly relative to the full-order model, taking less than 1 min to obtain 1.0 s of random vibration data using an iterative Newmark integration method. The two models used different solvers and codes, so this comparison is not an absolute measure of their performance. This comparison is meant to show the impact of the different methods on the workflow and expected savings given the available simulation tools.

A model calibration procedure has been proposed to effectively and efficiently determine the parameters of a four-parameter Iwan element used to capture the effect of nonlinear joint physics in structural dynamics reduced order models with whole joint representations. A high-fidelity model with frictional contact is developed and the quasi-static modal analysis approach is utilized to determine the nonlinear modal characteristics of the preloaded assembly. The results reveal how the modal frequencies and damping ratio change due to the presence of frictional contact within jointed structures. With the full-order model, a reduced order model is created based on the preload simulation by condensing the contact surface area down to a single, virtual node with multipoint constraints. A Hurty/Craig-Bampton substructure model is created by reducing the interior degrees-of-freedom away from the joint with fixed-interface modes while preserving the nonlinear virtual nodes with static constraint modes. The whole joint model is placed between connecting virtual nodes and is modeled with either linear or nonlinear one-dimensional constitutive elements depending on the physics of the structure. In this study, a four-parameter Iwan element is used to model the in-plane shear loads in the joint known to contribute most to the nonlinear energy dissipation in bolted connections. The parameters in the nonlinear whole joint model are then calibrated by utilizing a genetic algorithm that simulates the model with the quasi-static modal analysis technique over a broad range of the parameter space and modal response amplitudes. A multiobjective function is minimized based on the nonlinear modes of interest that the reduced order model is obligated to capture.

A nonlinear finite element model of a C-Beam assembly was used to demonstrate the calibration method and determine the in-plane Iwan element parameters to simultaneously match two nonlinear modes of interest, namely, the first and third in-phase bending modes. The NSGA-II algorithm focused the parameter search on solutions near the Pareto-front, leading to a more efficient calibration algorithm compared with traditional Monte Carlo simulations. An optimal parameter set was determined to match the two nonlinear modal characteristics of the full-order finite element model. With the calibrated model, the accuracy of the predictions was evaluated for two broadband excitations, namely, an impulsive haversine force and a white-noise broadband excitation with flat-line input. The results for the impulsive loading agreed well with the reference data from the full-order model. A modal decomposition of the time histories revealed the nonlinear characteristics of each mode and the importance of them on the nonlinear response. The random excitation simulations produced good agreement between the reduced and the full-order model; however, some slight discrepancies were observed due to coupling between two mode pairs that were not modeled as nonlinear (i.e., first and third out-of-phase bending). Regardless, the results agreed well, and the reduced order model showed orders of magnitude reduction in the computational cost.

The methodology and results presented here revealed that multiple modes were able to be calibrated with a single Iwan element in the whole joint formulation. This is significant in that the reduced model accurately captures the overall effect that a physics-based model would have on the global response of the structure. In addition, the methodology showed that the reduced model can be calibrated without experimental data in an efficient way using only quasi-static simulations. An important aspect of the approach is that it can help speed up simulations and analysis for rapid design analysis when experimental hardware is not readily available. This method allows for a modeling approach that predicts and captures the nonlinear damping characteristics of a bolted/jointed assembly, given a predictive full-order model, without the requirement of test data to calibrate to. This leads to more of a predictive approach that helps design engineers critically evaluate joint designs under transient dynamics loads and ensure that the assembled structure provides sufficient energy dissipation for the operating environments of the system.

The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. The authors would like to thank Matthew Allen and Robert Lacayo for their insightful discussions throughout the development of this research and for allowing the authors to leverage their QSMA algorithms.

The submitted paper has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science Laboratory, is operated under Contract No. DE-AC02-06CH11357. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doe-public-accessplan

This work was sponsored by Sandia National Laboratories' Advanced Simulation and Computing (ASC) Physics and Engineering Models (PEM) subprogram.

There are no conflict of interest.