Robert Lacayo
email: lacayo@wisc.edu

Luca Pesaresi
email: luca.pesaresi12@imperial.ac.uk

Johann Gross

Daniel Fochler

Jason Armand

Loic Salles

Christoph Schwingshackl
email: c.schwingshackl@imperial.ac.uk

Matthew Allen
email: matt.allen@wisc.edu

Matthew Brake
email: brake@rice.edu

Nonlinear modeling of structures with bolted joints: a comparison of two approaches based on a time-domain and frequency-domain solver

Keywords: friction, harmonic balance, modal analysis, model updating, damping, nonlinear vibration

Nonlinear modeling of structures with bolted joints: A comparison of two approaches based on a time-domain and frequency-domain solver

Introduction

Mechanical joints are a fundamental element of many built-up structures ensuring structural integrity, component alignment, and often sealing. The presence of a bolted joint in a structure can also introduce additional energy dissipation and lead to a change in the overall stiffness, which directly influences the damping and natural frequencies of the structure. Given that most structures are expensive to construct and test, it is important to be able to generate models beforehand that can accurately predict the vibrational amplitudes of such structures in order to ensure accurate performance and life predictions. Experience shows that structures with bolted joints exhibit damping that changes with the amplitude of excitation [START_REF] Gaul | Nonlinear dynamics of structures assembled by bolted joints[END_REF][START_REF] Segalman | Handbook on dynamics of jointed structures[END_REF][START_REF] Hartwigsen | Experimental study of non-linear effects in a typical shear lap joint configuration[END_REF][START_REF] Jalali | Identification of nonlinear bolted lap-joint parameters by force-state mapping[END_REF][START_REF] Deaner | Application of viscous and Iwan modal damping models to experimental measurements from bolted structures[END_REF][START_REF] Roettgen | Nonlinear characterization of a bolted, industrial structure using a modal framework[END_REF][START_REF] Brake | The Mechanics of Jointed Structures: Recent Research and Open Challenges for Developing Predictive Models for Structural Dynamics[END_REF]. This adds an additional layer of complexity to their analysis since linear models are insufficient to capture the overall forced response of such a system. The addition of nonlinear constitutive models that give a better description of bolted joint mechanics is therefore needed for such analysis.

However, a bolted joint model that includes too detailed a description of the nonlinear tribological processes occurring at the interface is far too expensive to be used for nonlinear dynamic analysis of a complicated structure. Such models are typically studied by the solid mechanics community, and would involve separate and highly-refined meshes between the interfacing surfaces and other fixture mechanisms (bolts, washers, nuts, etc.), which are coupled with nonlinear friction and contact models [START_REF] Gray | A global bolted joint model for finite element analysis of load distributions in multi-bolt composite joints[END_REF][START_REF] Kim | Finite element analysis and modeling of structure with bolted joints[END_REF][START_REF] Reid | Detailed modeling of bolted joints with slippage[END_REF]. Other works within the same community investigate the fundamental science (friction, adhesion, fretting, elastic deformation, etc.) underlying contact mechanics, and typically involve simplified models such as a punch on a semi-infinite plane [START_REF] Dini | Comprehensive bounded asymptotic solutions for incomplete contacts in partial slip[END_REF][START_REF] Ciavarella | Reduced dependence on loading parameters in almost conforming contacts[END_REF][START_REF] Ahn | Shakedown of coupled two-dimensional discrete frictional systems[END_REF][START_REF] Flicek | Progress in the application of notch asymptotics to the understanding of complete contacts subject to fretting fatigue[END_REF]. These methodologies, however, cannot currently be scaled to three-dimensional dynamic analysis of engineering structures such as that studied here. Further, the length scale of the whole structure (where stresses are of interest to engineers) differs very greatly from that of the joint and so it seems wise to simplify the representation of the joint. Hence, the current state-of-the-art in dynamic joint modeling reflects a compromise between the level of detail in modeling joint kinematics versus the computational effort required to estimate the dynamic response.

The past few decades have seen the development of a number of alternative bolted joint modeling approaches that fall somewhere between low-cost and high fidelity. In the most simple case the assembled structure can be represented as linear mass-spring models, connected by one or a few nonlinear macro-elements composed of a parallel or series arrangement of springs, dashpots, and Coulomb sliders [START_REF] Gaul | Nonlinear dynamics of structures assembled by bolted joints[END_REF][START_REF] Griffin | Friction damping of resonant stresses in gas turbine engine airfoils[END_REF][START_REF] Sanliturk | Harmonic balance vibration analysis of turbine blades with friction dampers[END_REF]. At the other end of the spectrum are the highfidelity joint models produced by the solid mechanics community [START_REF] Gray | A global bolted joint model for finite element analysis of load distributions in multi-bolt composite joints[END_REF][START_REF] Kim | Finite element analysis and modeling of structure with bolted joints[END_REF][START_REF] Reid | Detailed modeling of bolted joints with slippage[END_REF]. A modeling approach that also leans heavily towards the detailed side is one that incorporates zero-thickness elements [START_REF] Goodman | A model for the mechanics of jointed rock[END_REF][START_REF] Kaliakin | Insight into deficiencies associated with commonly used zero-thickness interface elements[END_REF][START_REF] Süß | Investigation of a jointed friction oscillator using the multiharmonic balance method[END_REF] and thin-layer elements [START_REF] Desai | Thin-layer element for interfaces and joints[END_REF][START_REF] Mayer | Segment-to-segment contact elements for modelling joint interfaces in finite element analysis[END_REF][START_REF] Sharma | Analysis and implementation of thin-layer element for interfaces and joints[END_REF] in the contact interface. Originally proposed by the geomechanics community, these elements capture the stick-slip condition as they undergo shear deformation. Most approaches are somewhere in between these extreme cases, and are either based on time [START_REF] Segalman | Handbook on dynamics of jointed structures[END_REF][START_REF] Abad | Design of experiments and energy dissipation analysis for a contact mechanics 3D model of frictional bolted lap joints[END_REF][START_REF] Eriten | Nonlinear system identification of frictional effects in a beam with a bolted joint connection[END_REF] or frequency [START_REF] Bograd | Modeling the dynamics of mechanical joints[END_REF][START_REF] Panning | Optimization of the contact geometry between turbine blades and underplatform dampers with respect to friction damping[END_REF][START_REF] Krack | Vibration prediction of bladed disks coupled by friction joints[END_REF][START_REF] Firrone | The effect of underplatform dampers on the forced response of bladed disks by a coupled static/dynamic harmonic balance method[END_REF] domain solvers that simulate the nonlinear dynamic response of structures containing bolted joints.

Despite the research effort, there are currently no predictive models for joint dynamics and the current benchmark is in error of around 25% in terms of frequency estimation and two orders of magnitude for the damping predictions [START_REF] Segalman | Handbook on dynamics of jointed structures[END_REF][START_REF] Brake | The Mechanics of Jointed Structures: Recent Research and Open Challenges for Developing Predictive Models for Structural Dynamics[END_REF]. In this work, two joint modeling approaches have been selected for a detailed evaluation of their prediction capabilities: the whole-joint approach based on a time-domain solver, and a frequency-domain approach based on the multi-harmonic balance method (MHBM).

The whole-joint approach, introduced by Segalman in [START_REF] Segalman | Handbook on dynamics of jointed structures[END_REF][START_REF] Segalman | Modelling joint friction in structural dynamics[END_REF], represents an area of contact surrounding a bolt as a coupling of two rigid surfaces by a single nonlinear constitutive element, typically an Iwan element [START_REF] Iwan | A distributed-element model for hysteresis and its steady-state dynamic response[END_REF][START_REF] Segalman | A four-parameter Iwan model for lap-type joints[END_REF], which models microslip in a joint interface. Since it is difficult to predict the parameters of an Iwan element for a given geometry and material, an Iwan element is traditionally tuned to capture measurements of the hysteresis of the joint as a function of loading amplitude [START_REF] Segalman | Handbook on dynamics of jointed structures[END_REF][START_REF] Rajaei | Development of generalized Iwan model to simulate frictional contacts with variable normal loads[END_REF]. However, in structures that contain more than one joint, it becomes difficult to isolate the contribution of each joint to the hysteresis observed in the measurements. As a result, subsequent works have measured the net effect of all of the joints on each mode of vibration and have tuned an Iwan element to model each structural mode, which is assumed to be uncoupled from all other modes [START_REF] Deaner | Application of viscous and Iwan modal damping models to experimental measurements from bolted structures[END_REF][START_REF] Roettgen | Nonlinear characterization of a bolted, industrial structure using a modal framework[END_REF][START_REF] Segalman | A modal approach to modeling spatially distributed vibration energy dissipation[END_REF]. This approach is often justified in structures that exhibit weak nonlinearities such that there is negligible change in the mode shapes of the structure [START_REF] Deaner | Application of viscous and Iwan modal damping models to experimental measurements from bolted structures[END_REF][START_REF] Roettgen | Nonlinear characterization of a bolted, industrial structure using a modal framework[END_REF] yet the natural frequency and the damping of each mode vary with amplitude. Hence these modal Iwan models have been tuned to capture this amplitude-change in frequency and damping of each mode.

However, the work of [START_REF] Lacayo | A numerical study on the limitations of modal Iwan models for impulsive excitations[END_REF] has demonstrated that the uncoupled modes assumption can lead to error in the predicted modal damping as amplitudes become large, and a systematic method has yet to be developed that reintroduces coupling to these modal Iwan models. Hence, a different tactic is pursued in this work that returns to the whole-joint formulation using discrete Iwan elements for modeling the physical joint, but still updates those Iwan elements so that the change in modal frequency and damping predicted by simulation matches that which is measured. Even then, despite the model containing only a few Iwan elements, it is still expensive to compute the transient response, and one would have to perform lengthy simulations of the nonlinear ring-down response of a structural model in order to extract these amplitude-dependent properties [START_REF] Lacayo | A numerical study on the limitations of modal Iwan models for impulsive excitations[END_REF]. To address this, a quasi-static technique was recently developed in [START_REF] Festjens | A numerical tool for the design of assembled structures under dynamic loads[END_REF] that was later simplified in [START_REF] Lacayo | Updating structural models containing nonlinear Iwan joints using quasi-static modal analysis[END_REF] to create a highly efficient and accurate algorithm for computing amplitudedependent frequency and damping for models containing Iwan elements. The authors of [START_REF] Lacayo | Updating structural models containing nonlinear Iwan joints using quasi-static modal analysis[END_REF] then used the algorithm to update the Iwan parameters in a finite element model to match the nonlinearity seen in experimental measurements on the transient free-response of a bolted structure. It is not clear, however, whether the algorithm could be used for updating against forced sinusoid response measurements, and this paper seeks to address that.

The second selected joint modeling approach is based on a multi-harmonic balance solver [START_REF] Petrov | Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks[END_REF][START_REF] Petrov | State-of-the-art dynamic analysis for non-linear gas turbine structures[END_REF][START_REF] Petrov | A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces[END_REF]. The forced response of the system is obtained through a frequency-domain solver, which uses a multi-harmonic description of the displacements and nonlinear contact forces. When the MHBM is used, the contact forces are computed in the time domain by following an alternating frequency-time procedure [START_REF] Salles | Dual time stepping algorithms with the high order balance method for contact interfaces with fretting-wear[END_REF][START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF], which allows for a wide variety of contact nonlinearities, including Coulomb friction, to be computed. The linear components are thereby represented via reduced finite element models [START_REF] Besselink | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF] that only maintain the information at the contact nodes, thus significantly reducing the computational time. This modeling approach was successfully applied to a series of industrial test cases, including underplatform dampers [START_REF] Armand | A multiscale approach for nonlinear dynamic response predictions with fretting wear[END_REF][START_REF] Pesaresi | Modelling the nonlinear behavior of an underplatform damper test rig for turbine applications[END_REF] and flange joints [START_REF] Di Maio | Development of a test planning methodology for performing experimental model validation of bolted flanges[END_REF], where a highly-detailed model was able to capture the amplitude-dependent damping and frequency under high-level excitation. One of the major constraints in these previous studies was the effective generation of the nonlinear mesh itself, and a novel approach, based on a Voronoi discretization [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF], is introduced in this study to make the model setup faster while maintaining good accuracy. Thus, the main advantage associated with detailed MHBM modeling is that an accurate description of the interface can be achieved at a reasonable computational cost.

In this work, the proposed changes in the above two approaches were implemented to provide them with better capability for practical updating and validation against forced response measurements. The two approaches are then used independently to create two models of the lap-joint benchmark structure, the Brake-Reuß beam [START_REF] Smith | The effects of boundary conditions, measurement techniques, and excitation type on measurements of the properties of mechanical joints[END_REF][START_REF] Dossogne | Experimental assessment of the influence of interface geometries on structural dynamic response[END_REF], and allow, for the first time, a detailed evaluation of their capabilities and restrictions in capturing the dynamic behavior of a bolted lap joint. A linear finite element model of the Brake-Reuß beam was created initially (Section 2), and its linear modes were analyzed to gain insight into the shapes and natural frequencies of its first three bending modes. Section 3 describes the development of the two joint models, wherein the linear model is augmented with a nonlinear component tailored to each approach. Then, for each joint model, a framework for updating is demonstrated with the goal of matching the forced response measurements of the Brake-Reuß beam around the resonance of the first bending mode. Simulations were then conducted on the updated models to assess their ability to predict the forced response of the second and third bending modes, as Section 4 details. A quantitative comparison between the simulation results for the two approaches is given, followed by a discussion in Section 5 of their advantages and disadvantages in capturing the nonlinear dynamic behavior of a bolted lap joint.

The Brake-Reuß Beam

The Brake-Reuß beam [START_REF] Smith | The effects of boundary conditions, measurement techniques, and excitation type on measurements of the properties of mechanical joints[END_REF][START_REF] Dossogne | Experimental assessment of the influence of interface geometries on structural dynamic response[END_REF] was chosen as the benchmark structure to help facilitate the comparison between the two modeling approaches. The structure consists of two identical, prismatic beams that possess half-height extensions near their ends. These half-height extensions overlap and are bolted to each other with three bolts to form a lap joint, as seen in Fig. 1. The material of both beams is stainless steel alloy 304, and the assembled beam has dimensions of 720.9 × 25.4 × 25.4 mm. Each of the three bolt assemblies consists of a 5/16-24 SAE grade 8 steel hex bolt and nut pair with washers under the bolt head and the nut. One of the main ideas of this design is to keep the two lateral halves of the beam relatively long in comparison to the length of the central lap joint. This ensures that, when the beam is excited, the relatively large kinetic energy seen at the ends of the beam stresses the central lap joint, causing it to activate friction more easily.

Linear Finite Element Model

A linear finite element model of the Brake-Reuß beam, depicted in Fig. 2, was developed in Abaqus to provide the dynamic response of the individual components before assembly. The finite element model consisted of uncoupled meshes for the two beam components as well as for the bolt assemblies. The bolts were included in the model to capture their stiffness and mass effects on the overall model response [START_REF] Schwingshackl | Modeling and validation of the nonlinear dynamic behavior of bolted flange joints[END_REF], and they are the mechanism that clamps the two beams together and applies a pressure on the contact interface. The bolt, nut, and washer assemblies were each lumped into a single mesh since slipping is unlikely to occur between these parts when clamped due to the extremely high contact pressure [START_REF] Schwingshackl | Modeling and validation of the nonlinear dynamic behavior of bolted flange joints[END_REF]. Care was taken to ensure a matching mesh at the contact interface between the two beam components and at the six contacts between the washers and the beam components. Each beam component mesh contains 8240 8-node hexahedron elements in addition to 340 6-node wedge elements on account of the shape irregularity introduced by the circular bolt holes. The elements for both the beam and the bolts use an isotropic, linear elastic material model with a Youngs modulus of 190 GPa, a Poissons ratio 0.32, and a mass density 7800 kg m -3 (stainless steel 304). These material properties were determined using a model updating routine to match the measured natural frequencies of the first three modes of only one of the beam components taken in free-free configuration.

A linear modal analysis was conducted on the finite element model of the assembled beam, with glue interfaces under the bolts, to understand the first three bending modes of the Brake-Reuß beam. Two configurations were considered, (i) an unrestrained and (ii) a fully stuck one, where the focus was on the bending modes orthogonal to the joint interface.

In the unrestrained configuration, there was no coupling between the two surfaces of the contact interface, and the two surfaces were free to separate or penetrate each other. However, the nodes under the washers were rigidly coupled to their coincident counterparts on the beam meshes to prevent rigid motion of the individual components. This configuration represents one extreme linear case of the model in which there is no interaction in the joint interface. In the stick configuration, the 592 coincident node pairs on the contact interface, in addition to those pairs under the bolt meshes, were rigidly coupled. This represented the other extreme case where the contact interface is completely stuck with infinite contact stiffness.

The natural frequencies of the first three bending modes are listed in Table 1 for the two configurations under free-free boundary conditions. As expected, a significant increase in natural frequencies can be observed for the fully coupled interface, which provides some quantitative insight on the significance of the joint on the response. The mode shapes for the unrestrained configuration are shown in Fig. 3. Since these two configurations represent two linear extremes in the model, it can be assumed that any nonlinear forced response analysis will result in frequency response functions (FRFs) having resonances somewhere in between the listed natural frequencies. The unrestrained configuration was used as the base linear model for the nonlinear analyses in Sections 2.2 and 3. Owing to its free interface, this configuration allows the addition of contact elements to model the nonlinear constraint introduced by the frictional interface.

Nonlinear Static Pressure Analysis

For the MHBM approach, it was necessary to have an accurate description of the contact pressure and gap information at the interface because they are required parameters of the 3D nonlinear contact elements in the joint model. This analysis allows the calculation of the static load on each contact element, whereas the dynamic loading is accounted for by the contact element itself during the nonlinear multi-harmonic forced response calculation. For this reason, a nonlinear quasi-static analysis was performed in Abaqus by applying a pre-tension on the bolts. The bolt load, F , corresponding to an applied torque, T was computed via an empirical formula [START_REF]Determination of coefficient of friction of bolt/nut assemblies under specified conditions[END_REF],

T = F (0.159P + 0.578d 2 µ T + 0.5D f µ H ) , (1)
where P is the bolt pitch, d 2 is the nominal bolt diameter, and D f is the average of the bolt cap and shaft diameters. The thread friction coefficient and bolt head friction coefficient µ T = µ H = 0.6 were referenced from [START_REF] Schwingshackl | Measured and estimated friction interface parameters in a nonlinear dynamic analysis[END_REF]. The axial bolt pre-tension obtained from Eq. ( 1) for a 20 N m bolt torque is 11580 N. This preload was applied to the unrestrained configuration model using the bolt-load feature in Abaqus.

A surface-to-surface approach with a direct enforcement method (based on Lagrange multipliers) was used to solve the contact problem in Abaqus. A penalty method formulation with a friction coefficient µ = 0.6 (at room temperature) [START_REF] Schwingshackl | Measured and estimated friction interface parameters in a nonlinear dynamic analysis[END_REF] was used to describe the frictional behavior. The results of the analysis show that the pressure is maximum around the bolt holes and reduces concentrically, leading to an area of lower compression between the holes (see Fig. 4). The pressure tends to zero near both ends of the contact interface in the y direction, and a small gap appears at the free ends due to the Poisson effect. The pressure results of the nonlinear static analysis were used as a starting point to generate the nonlinear dynamic models, as explained in Section 3.

Experimental Setup

The same Brake-Reuß beam tested in [START_REF] Smith | The effects of boundary conditions, measurement techniques, and excitation type on measurements of the properties of mechanical joints[END_REF][START_REF] Dossogne | Experimental assessment of the influence of interface geometries on structural dynamic response[END_REF] was used to collect forced-response data for this work. Before testing, each bolt was tightened first to 14 N m using a torque wrench, with the center bolt tightened first before the outer two bolts, and then all again tightened to 20 N m in the same order. Care was taken to ensure a gap between the vertical surfaces at the ends of the half-width extensions by separating the surfaces with 14-pt cardstock (a business card). This bolting procedure led to the best repeatability in the joint [START_REF] Smith | The effects of boundary conditions, measurement techniques, and excitation type on measurements of the properties of mechanical joints[END_REF], which ensures reliable measurements. As shown in Fig. 5, the Brake-Reuß beam was suspended 50 mm from each end with fishing line hung in series with a twisted bungee cord. The length of the suspension from the top support to the beam is 1.7 m to fully decouple the beam from the support structure. The excitation was provided by a Brüel & Kjaer PM Type 4809 electrodynamic shaker, which sat on a large steel table whose ground supports contained vibration isolators. The threaded hole in the beam (point A in Fig. 1) was carefully aligned with the axis of the shaker rod to avoid placing a torque load on the beam. Although six PCB 35AO3 Triaxial ICP accelerometers were attached to the beam, only the response at the far corner from the shaker (point B in Fig. 1) was considered in this work. A PC-controlled LMS 16-channel spectral analyzer was used as both the signal generator for the shaker and the acquisition system for the accelerometers.

Output Accelerometer

Joint Interface Shaker Driving Point

Fishing Line The force-controlled, stepped-sine response was collected for three frequency bands (swept upwards): 150 to 170 Hz (peak for the first bending mode), 560 to 580 Hz (second bending), and 1150 to 1190 Hz (third bending). For each of these three bands, four tests were conducted at driving force magnitudes of 0.1, 0.5, 1.0, and 2.0 N. The twenty-times increase in force (up to 2 N) was the maximum range of amplitudes used for this study, which already encompasses a significant portion of the micro-slip regime. Higher forces were attempted but abandoned due to difficulties encountered in controlling a constant force near the resonances of the beam at higher excitation levels. The response was measured from the accelerometer at point B, and the acceleration frequency response function (intertance) curves were constructed for each force level in each frequency band. These curves served as the reference for the model updating and validation in this study.

To illustrate the nonlinear behavior of the BRB structure, the inertance curves for the first bending mode are shown in Fig. 6. The 0.1 N forcing case exhibits the most linear behavior as implied by the relatively symmetric peak. As the forcing magnitude increases, the peak resonance point shifts downward in frequency (softening) and the peak becomes shorter and wider (damping increases). Additionally, the shape of the peak loses symmetry and becomes slightly skewed towards the left, which reflects the softening behavior of the beam with increasing response amplitudes. These forced response trends are very typical in structures containing bolted joints and indicate a certain amount of sliding in the contact interface.

Considered Approaches for Modeling Lap Joints

The main challenge in attempting to capture the nonlinear dynamic behavior of a bolted lap joint lies in the representation of the contact interface, which depends on the constitutive element used to provide the nonlinearity. The two considered joint modeling approaches, the time domain whole-joint approach and the MHBM approach, each utilize a different type of nonlinear constitutive element. The MHBM approach uses a nonlinear element that describes a local stick-slip-separation condition on the interface, whereas the nonlinear element for the whole-joint approach effectively models the cumulative effect of many friction elements that slip one after another. Due to the significant differences in the two modeling approaches, a thorough discussion on how to introduce the joint nonlinearity to the linear finite element model is presented in the following subsections.

Time Domain Approach Using the Whole-Joint Model

The intent of the whole-joint approach is to capture the cumulative behavior of the underlying nonlinear contact phenomena occurring in the interface as a whole, and lump it into a single "whole-joint model." This usually involves some manner of physical reduction on the nodes within a contact area underneath a bolt, after which the two interfacing areas are coupled via a nonlinear constitutive element that describes the hysteresis of the entire joint as it flexes [START_REF] Segalman | Modelling joint friction in structural dynamics[END_REF].

A virtual node, representing one end of the nonlinear element, can be made dependent on the motion of the nodes within a contact patch using a multi-point constraint (MPC) defined such that the translational and rotational degrees of freedom of the virtual node represent an average of those in the contact area of the finite element model. The virtual node corresponding to the opposing contact area then serves as the other end point for relative motion in the nonlinear element, as seen in Fig. 7. The forces exerted by the nonlinear element on the virtual nodes are then distributed to the nodes in their respective contact areas using that same constraint mechanism. Further details on the whole-joint model for the Brake-Reuß beam is given in Section 3.1.1.

MPC Spider Contact Patch

Virtual Nodes Nonlinear Element

Model

The selected nonlinear element for the whole-joint model was Segalman's four-parameter Iwan model [START_REF] Segalman | A four-parameter Iwan model for lap-type joints[END_REF]. The Iwan model is one of the most common microslip elements in literature for its relatively simple implementation and its ability to capture power-law damping using a distribution of friction sliders that replicate joint microslip. In general, an Iwan model, as shown in Fig. 8a, is a one-dimensional, parallel implementation of a large (ideally infinite) number of friction slider (Jenkins) elements having different slip thresholds θ i , where the thresholds are defined by a statistical distribution function [START_REF] Iwan | A distributed-element model for hysteresis and its steady-state dynamic response[END_REF]. Segalman's four-parameter formulation defines the slider distribution such that the energy dissipated per cycle by the Iwan element in microslip follows a power-law relationship to the element force (Fig. 8b). The four-parameter Iwan model was originally developed after an extensive testing and modeling campaign [START_REF] Segalman | Handbook on dynamics of jointed structures[END_REF] revealed that the Coulomb model of friction was inadequate for capturing the power-law dissipation exhibited by experimental joints. This power-law relationship continues to be observed in more recent studies [START_REF] Deaner | Application of viscous and Iwan modal damping models to experimental measurements from bolted structures[END_REF][START_REF] Roettgen | Nonlinear characterization of a bolted, industrial structure using a modal framework[END_REF][START_REF] Lacayo | Updating structural models containing nonlinear Iwan joints using quasi-static modal analysis[END_REF], where Iwan models were used successfully to characterize joint behavior.

Despite its merits in modeling the macro-mechanics of a joint, the Iwan model is limited by its inability to accurately represent the contact interface in detail. The model carries a high degree of epistemic uncertainty due to this limitation, and ongoing research is attempting to address this shortcoming (see [START_REF] Brake | Reconciling whole joint models and the preservation of local kinematics[END_REF] for example). Nonetheless, since the dynamics of the structure as a whole is of modeling interest rather than the mechanics occurring in the joint, there is some confidence that the Iwan modal should capture the nonlinearity associated with at least one mode given the achievements of modal Iwan models [START_REF] Deaner | Application of viscous and Iwan modal damping models to experimental measurements from bolted structures[END_REF][START_REF] Roettgen | Nonlinear characterization of a bolted, industrial structure using a modal framework[END_REF]. As a whole-joint model composed of internal sliders, only a single Iwan element is sufficient to represent the hysteresis of the entire contact area. By physical analogy, each slider in the Iwan element could be thought to represent a microscopic region in the interface that undergoes slipping and sticking as the joint flexes. When only a fraction of all the sliders starts slipping, the Iwan element exhibits microslip behavior. When the final slider passes its slip threshold, the Iwan element as a whole enters macroslip and follows a constitutive relationship equivalent to Coulomb friction until load reversal. Since Segalman's four-parameter Iwan element is one-dimensional, the virtual nodes of a whole joint model would typically be coupled with multiple Iwan elements all oriented in the different translational and rotational directions. In addition, these Iwan elements would be supplemented with springs and dashpots placed in parallel to provide linear sources of joint stiffness and damping [START_REF] Deaner | Application of viscous and Iwan modal damping models to experimental measurements from bolted structures[END_REF]. However, it is challenging to determine a unique set of parameters for each of these elements, and, depending on the joint behavior, not all directions my be necessary. For the whole joint models of the Brake-Reuß beam, only a single Iwan element was oriented along the axial (y) direction of the beam since the investigated bending modes were expected to exercise friction along this direction predominantly. A single spring was also added along the surface-normal (x) direction to account for the additional stiffness due to the bolt preload. The glued bolt assemblies were assumed to carry the residual stiffness of the joints, so no additional springs were added in the remaining directions. The linear damping was assumed to be accounted for by the background modal damping, so no dashpots were added.

(a) ) , ( 1  t x k ) , ( 2  t x k ) , ( 3  t x k ) , ( 4  t x k  u F (b) log(Joint Force) ) e l c y C / n o i t a p i s s i D ( g o l microslip macroslip slope = 3 + χ
In a first-attempt arrangement of the Iwan models, the interface of the unrestrained-configuration of the finite element model was divided into three contact patches, one under each bolt, with the virtual nodes located at the center of each bolt hole on the interface. Unfortunately, as will be detailed in Section 3.1.3, it was found that the first bending mode hardly stretched the Iwan elements in this arrangement, causing them to exert too little change in the damping of the structure to be comparable with the measurements. This problem was attributed to the symmetry in the shape of the first bending mode (see Fig. 3), where there is hardly any relative motion between the two interfacing surfaces near the center of the joint, and most relative motion occurs far from the bolts where the beam is thinnest. Hence the outer two contact patches were further split to divide the interface into five patches, as shown in Fig. 9. This setup produced ample damping from the Iwan elements at the edges (joints 1 and 5), so it was used as the final configuration for the upcoming model updating. The four parameters of Segalman's Iwan model are F S , K T , χ, and β [START_REF] Segalman | A four-parameter Iwan model for lap-type joints[END_REF]. Here, F S is the slip force at the inception of macroslip, K T is the tangent stiffness at zero load equilibrium, χ controls the power-law slope of dissipation in microslip, and β relates to the ratio of the tangent stiffness at zero-load equilibrium versus the stiffness near the inception of macroslip. With five Iwan elements and five normal-direction springs, this model contained a total of 25 parameters. As a further means of reducing the problem, the parameters for joint 5 were assumed to be the same as those of joint 1, and likewise joint 4 was made the same as joint 2, based on the symmetry of the bending modes. This leaves 15 unique parameters to determine. Although the slip force, F S , and tangent stiffness, K T , parameters could be deduced based on first principles, the χ and β parameters are not as intuitive. There is currently a lack of literature to help correlate these parameters to design variables such as the joint geometry, materials, and surface properties.

As an alternative, the 15 parameters were selected after tuning them to experimental data. Specifically, the parameters were tuned to reproduce the measured amplitude-dependent frequency and damping changes observed in the first bending mode from Fig. 6. The precedent for tuning Iwan models to modal frequency and damping was set in [START_REF] Deaner | Application of viscous and Iwan modal damping models to experimental measurements from bolted structures[END_REF], which also gave a procedure for tuning such models individually to the response of each mode. In this work, however, the Iwan models are in the physical space, so the projection of their nonlinearity onto each mode is not straightforward enough to employ the same tuning procedure. Instead, a quasi-static method for predicting the modal response of the nonlinear model was used to update the joint models iteratively. See Section 3.1.3 for more details.

Mathematical Theory

The whole-joint approach is amenable to simulations in the time domain since there are only a few Iwan elements to resolve. The effect of an Iwan element on a finite element model manifests itself as a nonlinear force applied equally and oppositely on the coupled virtual nodes in the system equations of motion. In the physical space, the equations of motion are represented as

Mẍ(t) + C ẋ(t) + Kx(t) + f nl (x(t), ẋ(t)) = p(t), (2)
with the displacement, velocity, and acceleration vectors denoted accordingly as x, ẋ, and ẍ. M, C, and K are the linear mass, damping, and stiffness matrices, respectively. The vector p tracks time-varying external forces, and f nl is the vector of nonlinear forces that are applied between contacting nodes by discrete, force-constitutive joint elements.

To reduce the computational time for the nonlinear dynamic simulation, the finite element model (minus the Iwan elements and normal spring elements) was reduced using Hurty/Craig-Bampton component modes [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF]. This uses the Ritz transformation

x = Ψx, (3)
where Ψ is the transformation matrix composed of fixed-interface modes and constraint modes. For the unrestrained configuration finite element model, the 10 virtual nodes in the interface were kept as boundary nodes to form 60 constraint modes (six degrees of freedom per virtual node). The rest of the model was reduced to form 26 fixedinterface modes. The shapes of these fixed-interface modes span up to the fourth cantilivered bending mode from each beam component, which, when combined with the constraint modes, was deemed more than sufficient to capture the first three global bending modes of the full Brake-Reuß beam.

The reduction was completed using the structural dynamics solver, SIERRA/SD [START_REF]Sierra Structural Dynamics Development Team, Sierra structural dynamics-user's notes[END_REF] (a research code developed at Sandia National Laboratories), which includes a reduction algorithm based on the Hurty/Craig-Bampton method to produce the reduced mass matrix M = Ψ T MΨ and the reduced stiffness matrix K = Ψ T KΨ for the beam. These matrices were then imported into MATLAB, where the damping matrix C and the Iwan and spring elements were added. The transformed system equations of motion become

Mẍ + C ẋ + Kx + Ψ T f nl Ψx, Ψ ẋ = Ψ T p, (4)
in which x, ẋ, and ẍ are the generalized displacement, velocity, and acceleration coordinates that govern the fixedinterface modes and constraint modes. The damping matrix C was formed by projecting the diagonal matrix of linear, background modal damping ratios ζ 0r into the reduced space,

C = MΦdiag [2ζ 0r ω 0r ] Φ T M. ( 5
)
Here Φ is the matrix of linear, mass-normalized mode shapes φ r , from the eigenvector problem,

K0 -ω 2 0r M φ r = 0. (6)
The stick stiffness matrix is

K0 = K + Φ T (∂f nl /∂x) | x=0
, where (∂f nl /∂x)| x=0 is the gradient of joint forces (also known as the joint tangent stiffness) near equilibrium. The tangent stiffness near equilibrium for the Iwan element is exactly its K T parameter, so K0 was equivalently formed by substituting all the Iwan elements for linear springs having stiffness K T .

From previous ring-down testing on the experimental beam [START_REF] Dossogne | Experimental assessment of the influence of interface geometries on structural dynamic response[END_REF], the background damping ratios, ζ 0r , for the first, second and third bending modes were determined to be 0.2%, 0.1%, and 0.1%, respectively.

As a final means of reduction, and in preparation for the quasi-static technique for tuning the Iwan elements, the system equations of motion in Eq. ( 4) underwent another basis change to remove the rigid body modes. This basis change stems from the definition for constraining the rigid body modal coordinates,

Φ T R Mx = 0, (7)
where Φ R is a set of the six zero-energy mode shapes (mass-normalized) produced from Eq. ( 6). The non-trivial solutions to Eq. ( 7) lie in the null space of Φ T R M, and the set of all non-trivial solutions,

L = null(Φ T R M), (8)
forms the Ritz basis matrix L that is used to transform Eq. ( 4) in the same manner as Ψ before it. Substituting x = Lq and premultiplying by L T yields the final system equations of motion,

Mq + Ĉ q + Kq + L T Ψ T f nl (ΨLq, ΨL q) = L T Ψ T p. (9)
These equations can be solved in the time domain to obtain the nonlinear dynamic response of the system.

Model Updating

There is currently not a readily available means to intuitively predict the four parameters of Segalman's Iwan model, so the parameters were instead tuned based on the given inertance data. Effectively, the simulated response was made to match the measured change in the natural frequency and modal damping with amplitude when the first bending mode was selectively-excited at forced resonance.

An equivalent modal damping ratio (one that includes nonlinear damping due to the Iwan elements) that would allow the reduced model of the beam to simulate the inertance at resonance point of the rth mode can be calculated as

ζ r = Ψ B L φr φT r LΨ T A 2 |a B (ω r )| , (10)
in which Ψ B is the row in the transformation matrix associated with output point B, and Ψ A refers to the row associated with input point A (see Appendix A for a full derivation). The mass-normalized mode shape φr is computed from the eigenvector problem in Eq. ( 6), except that K and M are used as the stiffness and mass matrices. The value |a B (ω r )| is the magnitude value of the inertance FRF as measured from output point B at the resonance frequency ω r . The amplitude of the modal displacement coordinate governing this mode shape can be calculated as

|H r | = |P A | |a B (ω r )| ω 2 r Ψ B L φr , (11)
where |P A | is the amplitude of the harmonic force at input point A. By applying Eq. ( 10) and Eq. ( 11) on the resonance point for the first bending mode (measured at the four different excitation amplitudes), a set of four values each for ω r , ζ r , and |H r | was computed for that mode. In Fig. 10, the computed ω r and ζ r for the first bending mode are each plotted with respect to |H r | to observe how these measured values for the natural frequency and damping change with increasing amplitude.

Similar curves for amplitude-dependent frequency and damping were calculated for the joint model using quasistatic modal analysis [START_REF] Lacayo | Updating structural models containing nonlinear Iwan joints using quasi-static modal analysis[END_REF]. The technique, described in detail in Appendix A, uses Masing's rules to construct a modal force-deflection hysteresis using the quasi-static response to a monotonically-increasing load distributed over the entire finite element model in the shape of the mode of interest. The slope of the secant line through this hysteresis is the square of the predicted natural frequency ωr , and the equivalent damping ratio ζr is the ratio of the area enclosed by the hysteresis over the total modal kinetic energy for that load level. Similar to the measured properties, ωr and ζr were associated with a modal deflection η r computed for several different load levels.

The predicted ωr and ζr versus η r curves for the first bending mode can then be plotted together with the measured ω r and ζ r versus |H r | curves like those seen in Fig. 10. To tune the Iwan elements, the goal would be to adjust their parameters so that the predicted curves fit reasonably well with the measured curves for the first bending mode. Although a global optimizer -such as the genetic algorithm used in [START_REF] Kuether | Parameter estimation of joint models using global optimization[END_REF] -could achieve this, the parameters for the five joint models were instead adjusted by hand to gain insight on how they affect the predicted frequency and damping curves.

Initially, the parameters for all Iwan elements were set to F S = 1000 N, K T = 7×10 5 N mm -1 , χ = -0.6, and β = 1 based on the trends seen in [START_REF] Lacayo | A numerical study on the limitations of modal Iwan models for impulsive excitations[END_REF]. The resulting predicted ωr and ζr curves achieved a surprisingly close fit with the measured ω r and ζ r , but the damping for the 0.1 N amplitude case was far under-predicted. The 15 unique parameters were then adjusted iteratively in an attempt to achieve a better fit. In the course of this fitting procedure, the following observations were made:

1. Changing the parameters for joint 3 (see Fig. 9) had no effect on ωr and ζr . The curves showed more sensitivity to changes in the parameters of joints 1 and 5 than for joints 2 and 4.

2. Increasing F S , χ, or β for joints 1, 2, 4, and 5 generally decreased the nonlinear contribution to ζr and decreased the shift in ωr , making the model more linear.

3. Increasing K T for joints 1, 2, 4, and 5 also made the model more linear, but it additionally increased the lowamplitude, linear natural frequency up to a limit of convergence. This upper limit corresponds to the linear case where the Iwan model exhibits infinite stiffness.

4. Increasing the stiffness of the linear x-direction spring moved the ωr curve upward in frequency without changing the damping or the frequency shift with amplitude.

Observation 1 shows a trend of increasing sensitivity for joint models further away from the center of the beam. This makes sense given the spatial symmetry in the first bending mode about the center of the beam, which implies that the Iwan elements closer to the center of the beam experience less flexure. Observations 2 and 3 made the fitting a much less straightforward process. Essentially, none of the parameters revealed any distinct effects on the shapes of the predicted ωr and ζr curves that could help determine their value. This made it difficult to fit the model curves to the measured curves because a change in any parameter would upset the fit made by another parameter. This difficulty is exacerbated by the fact that three unique sets of joint parameters (recall that joint pairs 1+5 and 2+4 were assumed equivalent) were being fit to one set of measured data, which is analogous to an under-determined problem.

Through manual tuning, a set of parameters was eventually found that matched the damping well, as seen in Fig. 10, but the relatively small shift in natural frequency could not be matched accurately. While iterating on the parameters, the resulting model curves produced cases where either the frequency or the damping was fit better, which seemed to imply that the model was inhibited by a trade-off between the two properties.

To confirm this trade-off, and to see whether a better parameter set could be found, a Monte Carlo study was conducted over the possible parameter space. As before, the parameters for joints 1 and 5 were taken to be equal, as were those of joints 2 and 4, but the parameters for joint 3 were omitted from the study due to observation 1. The normal direction spring stiffness values were also omitted since they do not influence the damping or the frequency shift (per observation 4), and can be adjusted after the study to raise or lower the frequency curve as necessary. Of the remaining parameters, F S , was allowed to vary between 20 and 2000 N, K T between 5×10 4 and 5×10 6 N mm -1 , χ between -0.99 and -0.01, and β between 0.01 and 100. The parameter values were sampled logarithmically for 50,000 trials.

For each trial, the model ωr and ζr versus η r curves were computed and then interpolated along |H r | to retrieve the estimated values associated with the four measured load cases. The value of ωr associated with the 0.1 N load case was then subtracted from all four interpolated ωr values to compute the frequency shift relative to the 0.1 N load case. Finally, the root-mean-square (RMS) error was computed for the estimated log( ζr ) relative to the measured log(ζ r ) and for the estimated ωr shift relative to the measured ω r shift. All of these calculations were done in MATLAB ® . Figure 11 plots the RMS error in frequency shift versus the RMS error in the natural log of the equivalent damping ratio for all trials. Ideally, a low RMS error for both properties is desired, but Fig. 11 shows a Pareto-optimal limit marking a clear trade-off in ability to match one over the other. On one extreme is Fig. 10, where the damping can be matched, but the frequency shift only reaches half the range of that measured. On the other extreme, the frequency shift can be matched, but the damping grossly overestimates that measured one, as seen in Fig. 12. Even in attempt to achieve a balance in error (Fig. 13), one must accept a substantial amount of error in the damping curve.

The cause for the trade-off between the frequency shift and the damping is not understood, but it is suspected to be a limitation of the four-parameter Iwan model. It may be that the four-parameter Iwan model is too simplistic in its formulation, and a model with more physics built into it may be more appropriate. Of course, other epistemic uncertainties could be important contributors to these errors, such as the choice of using a Iwan element as the friction model, the choice of dividing the interface into five contact patches with a nonlinear element oriented in one direction only, or even the choice of using the whole-joint formulation with tied nodes defined by averaging MPCs. The choices made in the above analysis were based on engineering judgment given the limited data and computational resources available to reasonably capture the shift in the resonance peak of the first bending mode.

Barring a complete redesign of the joint model, any of these Pareto-optimal cases could be acceptable depending on the objectives of the joint model. For this study, priority was given to match the damping because the observed frequency shift was relatively small while the damping had a direct and strong influence on the peak FRF at resonance. The joint parameters (listed in Table 2) that led to the predicted curves in Fig. 10 were deemed sufficient and used as the final configuration for the whole-joint approach.

Node-to-Node Approach Using Multi-Harmonic Balance Method (MHBM)

In contrast to the whole-joint approach from the previous section, the MHBM approach models the interface explicitly by coupling a large number coincident nodes between the two surfaces each with a friction contact element. Specifically, the 3D contact element developed by [START_REF] Petrov | Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks[END_REF] was used for the nonlinear model in this case.

Model

Each 3D contact element (see Fig. 14) includes three springs: two oriented along each of the two tangential directions to describe the tangential contact stiffness affecting the in-plane motion (2D Jenkins element), and a spring oriented in the normal direction to represent the normal contact stiffness of the interface asperities. All springs are connected in to the same Coulomb friction slider to allow stick, slip, separation or permanent gap states during a vibration cycle [START_REF] Petrov | Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks[END_REF]. The spatial variation of the contact pressures predicted from the nonlinear static analysis is paramount for accurately capturing the microslip dynamics of the interface. Recent work by [START_REF] Dossogne | Experimental assessment of the influence of interface geometries on structural dynamic response[END_REF][START_REF] Goyder | Measurement of damping due to bolted joints[END_REF] indicates that most of the dissipation occurs away from the bolted joints where the contact pressure is low. While the contact pressure in these regions is known to fluctuate with time during dynamic response [START_REF] Seeger | In situ measurements of interfacial contact pressure during impact hammer tests[END_REF], a full coupling between the normal and tangential directions ensures that this fluctuation is captured and that the tangential forces are modified accordingly.

The 3D contact element can be described by intuitive physical parameters, which are the friction coefficient µ, the initial local preload N 0 , the per-unit-area stiffnesses in the normal (k n ) and two tangent (k t ) directions, and a gap parameter to help describe the initial separation between the two coupled nodes observed in Fig. 4. Unlike the Iwan element, whose microslip transition is formulated internally, here this transition is captured by a grid of elements that

Joint

Iwan Element (y-direction) Spring (x-direction) Typically, all finite-element node pairs at the contact interface would be coupled with a contact element, but a smaller set was used in this work in order to reduce the computational time for the model updating. The nonlinear nodes were thereby selected so as to cover the interface almost uniformly (see Fig. 15a); a slightly higher density of nonlinear nodes was given in low-pressure regions where dissipation was likely to happen, whereas a lower density of nodes was given around the bolt holes where stick was expected due to high contact pressure. A total of 86 nodes were selected, as shown by the white dots in Fig. 15.

F S [N] K T [N mm -1 ] χ β k x 1,
The values for the normal load N 0 in Fig. 15a and the gap parameter in 15b are different for each contact element, and they were based on the static pressure calculation performed in Abaqus (see Section 2.2). To get N 0 for each node, the pressure in the region nearest to the nonlinear node was reduced to a normal force acting on that node. This required a partitioning of the contact interface into areas whose borders were optimized by proximity to nearby nonlinear nodes, which was achieved via a Voronoi polygon calculation [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF]. The polygons were computed using an algorithm, based on Dirichlet tessellation, implemented in the MATLAB ® function, voronoi. The resulting Voronoi diagram, representing the nonlinear mesh for the MHBM, is shown in Fig. 15.

If computed rigorously, N 0 would then be the integral of the static pressure over the area enclosed by each Voronoi polygon. However, Abaqus only exports nodal forces, so N 0 for each node was computed as the sum of the normal forces acting on the nodes inside each Voronoi polygon. Similarly the gap parameter for each nonlinear node was taken as an average of the gap distance measured between the coincident nodes in each polygon.

Parameters k n and k t both describe per-unit-area stiffness, so the actual stiffnesses associated with the 3D contact element is a multiplication of these values with the area of the Voronoi polygon. The values of k n and k t can generally be measured in a specifically-designed test rig [START_REF] Schwingshackl | Measured and estimated friction interface parameters in a nonlinear dynamic analysis[END_REF]. In this study, a simplification was done by assuming that k n and k t are equal, since the same order of magnitude is expected in elastic axi-symmetric contacts [START_REF] Medina | Analytical and numerical models for tangential stiffness of rough elastic contacts[END_REF]. Unfortunately an initial investigation, based on the measured values, led to a strong deviation from the measured FRF curves, so the stiffness values were instead obtained via an updating process similar to the one described in [START_REF] Medina | Analytical and numerical models for tangential stiffness of rough elastic contacts[END_REF]. The friction coefficient was kept constant at µ = 0.6, which was expected for a dry, room temperature, steel-to-steel contact undergoing dynamic motion [START_REF] Schwingshackl | Measured and estimated friction interface parameters in a nonlinear dynamic analysis[END_REF]. For all dynamic simulations with this model, a single harmonic excitation was then applied to the shaker input location (point A in Fig. 1), and the response was calculated at the accelerometer location (point B) by taking multiple response harmonics into account.

Mathematical Background

To complete the dynamic simulations, a solver was used that combines a model reduction technique with the multi-harmonic balance method (MHBM) [START_REF] Salles | Dual time stepping algorithms with the high order balance method for contact interfaces with fretting-wear[END_REF]. A single-harmonic balance approach [START_REF] Bograd | Modeling the dynamics of mechanical joints[END_REF] was considered for this study, but it was found to be unable to capture the complex behavior occurring at the interface [START_REF] Salles | Numerical round robin for prediction of dissipation in lap joints[END_REF][START_REF] Gross | A numerical round robin for the prediction of the dynamics of jointed structures[END_REF]. Therefore a multi-harmonic balance solver, which forms part of the software code, FORSE (Forced Response Suite) [START_REF] Petrov | Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks[END_REF][START_REF] Petrov | State-of-the-art dynamic analysis for non-linear gas turbine structures[END_REF][START_REF] Petrov | A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces[END_REF], developed at Imperial College London for industrial scale use [START_REF] Pesaresi | Modelling the nonlinear behavior of an underplatform damper test rig for turbine applications[END_REF][START_REF] Sever | Experimental and numerical investigation of rotating bladed disk forced response using underplatform friction dampers[END_REF], was used.

Following the MHBM procedure, the response for each degree of freedom (DOF) in the system, expressed by Eq. ( 2), is represented as a Fourier series including n harmonics,

x(t) = X 0 + n j=1 X c j cos(jω 1 t) + X s j sin(jω 1 t), (12)
where X c j and X s j are the harmonic coefficients for each DOF, and ω 1 is the principal frequency of oscillation. Substituting Eq. ( 12) into Eq. ( 2), and balancing the harmonic terms with a Galerkin projection, a system of equations to determine all harmonic components can be obtained in the frequency domain as follows,

Z(ω)X + F(X) = P, (13)
where

X = X 0 T , X c 1 T , X s 1 T , • • • , X c n T , X s n T T
is the vector of all harmonic coefficients, F and P are the harmonic representations of f nl and p, respectively, and Z is the dynamic stiffness matrix of the linear part of the system. The model reduction described in [START_REF] Petrov | A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces[END_REF] is applied on this system, and the harmonic coefficients are calculated using a Newton-Raphson iterative solver. With each iteration k, the vector of harmonic coefficients is updated as

X (k+1) = X (k) -Z(ω) + ∂F(X (k) ) ∂X -1 Z(ω)X (k) + F(X (k) ) -P , (14)
where ∂F(X (k) )/∂X is the harmonic tangent stiffness matrix. In [START_REF] Petrov | Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks[END_REF], the authors derive an analytical expression for ∂F(X (k) )/∂X for the contact element discussed in Section 3.2.1.

Model Updating

After setting up the model as described in Section 3.2.1, a preliminary study on the first bending mode was carried out to evaluate the performance of the model against experiments and update it where necessary. In an initial study, a simplification was made to the normal load and gap parameters for some of the nonlinear elements. Specifically, the nonlinear nodes in the areas of zero pressure near the left and right edges of the interface (see Fig. 4a) were assumed to be initially in contact but completely unloaded (i.e. zero gap and N 0 = 0). This assumption constitutes a simplification because it represents the standard approach adopted if no gap information was available. As described below, this seemingly small simplification has enormous consequences on the simulated response of the model, and the results are indicative of the potential errors that can arise from the zero-gap assumption.

An initial value of 5000 N mm -3 was used for the contact stiffnesses k n and k t , leading to the experimental and simulated nonlinear FRFs for the 0.1 N excitation level (the most "linear" response) in Fig. 16. It can be observed that the initial model overestimates the resonance frequency and under-predicts the peak amplitude, which implies that much more energy dissipation is predicted as compared to that of the experimental FRF. To understand this behavior better, Fig. 17 shows the average contact static gap and pressure for different "columns" (parallel to the left edge) of nonlinear nodes. The pressure values were computed as the sum of the nodal forces applied on the five nonlinear nodes in one column divided by the total area of the Voronoi polygons surrounding those nodes. The nodes that are located closer to the edge show an increasing average gap, while the average pressure shows 4 MPa for the column furthest from the edge but very quickly drops to almost zero closer to the edge. This shows that the simplifying assumption of initial contact by neglecting the gap in the left and right extremities of the interface is inaccurate; the friction elements in the gap region experience contact and sliding during a vibration cycle, leading to an overestimation of the energy dissipation and the stiffness of the joint. To overcome this issue, a new, updated model was then introduced that included the initial gap for nodes near the edges. The experimental and simulated nonlinear FRFs of this updated model for different excitation levels, are shown in Fig. 18a. The results show much better agreement with the experimental results for both resonance frequency and amplitude. However, a hardening effect can be observed for the model at an excitation of 2 N that is absent in the experiments. This effect is caused by the closing of the gap at high amplitudes, leading to a sudden increase in stiffness for the joint. When all the elements in this gap region were removed from the model (see Fig. 18b), a much better agreement was achieved. This implies that the actual gap might be larger than the predicted one such that this gap does not close during a vibration cycle.

To further improve the model, a study on how the contact stiffnesses (k t , k n ) affect the resonance frequency of the first bending mode was performed for the lowest excitation level (0.1 N), as shown in Fig. 18c. A value of k t = k n = 6000 N mm -3 was selected for the final model since it produced a natural frequency in close agreement with that of the experiments. Unlike in the whole-joint model, where all the excitation levels of the first mode were used for the updating procedure, here only the low level excitation measurement (0.1 N) was used to update the contact stiffness. For this reason, all the higher excitation levels of the first mode can be regarded as predictions.

Simulation Results for the Updated Models

The Whole-Joint Model

Equation ( 9) was integrated numerically using an average-acceleration Newmark algorithm [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF] to compute the response over time to a harmonic input at point A (refer Fig. 2). To produce an FRF, and therefore make the output of the whole-joint model comparable to the MHBM results, the Newmark algorithm was placed inside of a loop that checked for steady-state convergence in the response history every 20 cycles. When convergence was reached, a sine wave with a frequency equal to that of the driving frequency was fit to the time signal. The amplitude and phase shift of the fitted sine wave was then converted into a complex amplitude associated with that driving frequency. The steadystate loop was repeated for each driving frequency of interest (within the three bands measured experimentally), and all again repeated for each forcing amplitude (0.1, 0.5, 1.0, and 2.0 N). The resulting simulated FRFs for the first three bending modes are shown in Fig. 19 along with the experimental FRFs. The joint parameters used in this model are those listed in Table 2.

The simulated curves for the first bending mode show good agreement, though this was expected because the model was fitted to this mode. Since the model was fitted with more weight on the damping (see Fig. 10), the heights of the peaks of each simulated FRF are almost identical to those measured in experiment. In contrast, the shift in the resonance frequency only covers half the measured range, which is consistent with the natural frequency curve in Fig. 10a. However, because the frequency shift is small, this error is considered acceptable.

To explore the predictive capabilities of the developed model, additional bending modes where investigated. For the second bending mode, the simulated FRFs show surprisingly good agreement considering that the model was not tuned to this mode. The maximum amplitude of the simulated FRF for the 0.1 N load case is somewhat underpredicted, but the peak heights agree quite well for the other load cases. Similar to first bending, the discrepancy in the absolute value of the resonance frequency is small, but the resonance frequency of the simulated FRF for the 2.0 N load case appears to align well with the measured one.

In contrast to the first and second bending modes, the simulated FRFs for the third bending mode disagree the most with those measured. The simulation is not able to capture the trend of increasing peak heights as measured from experiments, instead showing an increase in damping at higher amplitudes. The predicted resonance frequencies are somewhat below the measured ones (2.5% underestimation), but the model was able to capture the softening effect. The dashed lines in Fig. 19 represent the curves simulated from the quasi-static modal analysis. For Fig. 19a in particular, this curve uses the exact same data used to construct the model damping curve shown in Fig. 10b, but the data has been recast as peak FRF versus frequency following the derivation in Appendix A. As can be seen, the QSMA curves do well in passing through the peak points of the simulated FRFs, effectively following the backbone curve of the FRFs. Considering that the QSMA curves were computed in a fraction of a second using a desktop computer, these backbone curves can be used reliably to quickly assess the ability for a whole-joint model to capture the measured peak response, which is demonstrative of its efficiency in model updating.

The MHBM Model

The updated MHBM model from Section 3.2.3 was used to simulate the first three bending modes of the Brake-Reuß beam, the response for which was compared to the experimental data (see Fig. 20). The simulations for the first bending mode are in good agreement with the experimental curves since the model was updated with regards to these curves. The model accurately captures both the slight frequency shift and the amplitudes at the different excitation levels.

With regards to the predictive capabilities of the MHBM model, the amount of shift in the resonance frequency for the second bending mode is captured quite accurately, even if a slightly lower resonance frequency is predicted (0.7% underestimation). The amplitude agreement is somewhat poorer for the lightly damped modes, but significantly improves for higher excitation levels.

The frequencies of the third bending mode are slightly underestimated, and, similarly to the whole-joint modeling approach, the amplitude predictions show an increase in damping with the excitation, which is in contrast to the experimental findings.

Comparison Between Approaches

To compare the two approaches, the key properties of each FRF peak -namely the peak value at resonance and its frequency shift with increasing excitation level -were examined. For the MHBM approach, these properties were retrieved by finding the highest response point in each FRFs. The same could be done for the simulated FRFs from the whole-joint model, but the backbone curves computed using quasi-static modal analysis were used instead since they represent the amplitude trends at resonance more continuously. Figure 21 shows the shift in the resonance frequency (relative to the 0.1 N case) and peak response amplitude over different input forces for the first three bending modes. This allowed for a summary comparison between the responses of the two approaches and the experiments.

In addition, Table 3 lists the percent error in absolute value of the resonant frequency, the resonant frequency shift, and the peak FRF amplitude estimated by the modeling approaches relative to the experimentally measured values.

Overall, a very good accuracy in amplitude estimation is obtained for the first bending mode, which reflects the successful model updating procedure, but both approaches tend to lose some accuracy for the second and particularly the third bending mode. With regards to the frequency shift prediction, the MHBM approach is more accurate in the first two modes, with an average difference of 13.2% (first bending) and 8.9% (second bending) compared to the difference of 50% and 63.6% from the whole-joint approach. The latter, however, better predicts the frequency shift for the third mode with an average difference of 16.5% compared to the 38.5% of the MHBM. As for the absolute frequency at the 0.1 N excitation level, the whole-joint model is closer for the second mode, being only 0.22% lower than the experiments compared to the the -0.7% of the MHBM, but the MHBM model is closer for the third bending with a -1.19% difference compared to the -2.47% of the whole-joint model. Both approaches seem to slightly under- predict the absolute frequency of the higher modes, which indicate softer joint models compared to the real one. The whole-joint approach is better at capturing the peak amplitudes, having a difference to the experiments on average 1.5%, 4.9%, and 45.8% for the three modes compared to an average difference of 4.7%, 13.7%, and 59.2% obtained for the MHBM approach. As already discussed, neither approach was able to predict the increasing peak amplitude trend of the third bending mode well. Disregarding the fact that both joint models estimate lighter damping than suggested by the experimental FRFs, the joint models predict a decrease in the peak inertance for third bending as the excitation force increases, whereas the reverse trend is observed in experiment. The cause for this reverse trend is unknown, but it could be attributed to the relatively large flexure of the joint in third bending, as illustrated in Fig. 3c, that may add an additional, as of yet unexplored, mechanism to the joint.

Discussion

In creating and updating the joint models, and in comparing the results between the two approaches, a number of insights on their advantages and disadvantages were gained.

Whole Joint Modeling approach

The main advantage for the whole-joint approach lies in its computational efficiency for model updating using the quasi-static modal analysis. For the joint model discussed in Section 3.1.3, which contains only five Iwan elements, a MATLAB implementation of the quasi-static method could compute 75 load cases on a curve for amplitude-dependent damping in only a fraction of a second using a desktop computer with a 2.93 GHz processor. Such speed makes it practical to run an overnight Monte Carlo analysis that computes tens of thousands of curves associated with different Iwan parameter sets. Depending on the objective function, an optimal parameter set could be found for a model that produces an accurate approximation of the measured forced response. As a time-domain method, the validated model could then be used to predict other types of responses such as shock, ring-down decay, and random vibration.

Although computing the steady-state response to forced input is comparatively expensive (∼2 hours per FRF curve in Fig. 19), quasi-static modal analysis can be used to approximate the backbone curve for an FRF peak in a matter of seconds. This capability allows a time-domain model to be readily comparable to frequency-domain models by examining the forced resonance point of the simulated response. Furthermore, whereas each FRF of a different forcing amplitude contributes one datum to an amplitude curve (like those shown in Fig. 21), quasi-static modal analysis can quickly compute dozens of load steps for more smooth amplitude curves.

The drawback to the whole-joint approach is its lack of detail in the description of the local kinematics and its heavy reliance on experimental data. Due to the lumped nature of a whole joint, many assumptions must be made about the location of the contact, the method of reducing the contact, and the choice of the force-constitutive elements to use. Even the number of, and relative size of, the contact patches is an open question. These assumptions hold many epistemic uncertainties that are difficult to resolve until the model is validated with experimental measurements.

Segalman's four-parameter Iwan element helps to reduce some uncertainties through its ability to reproduce the power-law relationship exhibited by the measured damping. However, its parameters are neither intuitive nor able to be measured directly, and they must be iterated upon to determine a set of values that fits the measurements well. Even then, the energy dissipated by the four-parameter Iwan element is coupled with its tangent stiffness, which, for the whole-joint model examined in this work, resulted in a fit that could be tailored to either the damping or the frequency shift, but not both. This may imply that the four-parameter Iwan model lacks a certain complexity in the physics it attempts to describe, and a model with more physics built into it may be needed.

Multi-Harmonic Balance model

The joint model of the multi-harmonic balance method attempts to represent the friction interface in a more physical manner. The contact elements used resemble a series-spring-Coulomb model whose parameters can be estimated by referencing literature and performing a static pressure analysis on the contact interface, which lessens their reliance on experimental data. By using many friction elements spread over a contact interface, these methods can capture microslip behavior by tracking in detail the locations of stick, slip, and separation on an interface.

The node-to-node approach with 3D contact elements does rely strongly on the correct calculation of the initial static pressure and gap in the interface. The work shown in Section 3.2.3 has shown that this approach is sensitive to changes in the static normal load and gap parameters, and some errors can be introduced if the contact in the real interface differs from the static calculation. This could occur when low tolerances in the manufacturing process lead to a lack of planarity of the surfaces, which can result in a different conformity at the contact interface. For this reason, a high level of control is required during the manufacturing and assembly of the test article in order to produce a configuration that the nonlinear model is more capable of predicting.

If brought into the time domain, the presence of tens, if not hundreds, of friction elements would result in a massive computational effort on the part of the numerical integrator, and possibly strain convergence. In the frequency domain, however, the Jacobian of the nonlinear contact element could have an analytical form, as is the case for the 3D contact element in the MHBM approach. This helps accelerate convergence in solvers, making such contact elements far more efficiently used in harmonic balance methods than for time-domain methods. The drawback is that the simulations are restricted to periodic inputs, though this is not a problem if the harmonic response of a structure is desired.

Lessons Learned in Joints Modeling

In the bigger picture, the two modeling approaches do have their appropriate domains of application. Due to their lumped formulation, whole-joint models are best used if the influence that the joint exerts on the parent structure is of greater interest than the dynamics occurring in the joint itself. If the response of the structure can be measured, then a time-domain model can be tuned quickly using quasi-static modal analysis, and such a model can be readily used to simulate the structural response to any sort of input excitation. In absence of experimental data, the interface should be modeled with a more intuitive, detailed contact formulation, similar to the MHBM model proposed here, for which the response can be computed more efficiently using the harmonic balance method. Node-to-node contact has the potential to reduce the computational error when compared to the whole-joint approach, but some difficulties can arise on the correct selection of the density and distribution of the contact nodes in order to find the optimum compromise between accuracy and computational speed.

From the analysis completed in this work, several recommendations can be made regarding future modeling efforts. For time-domain methods, it is paramount to improve the constitutive models used in the interface. Until geometric, tribological, and material properties can be used to deduce the parameters of a lumped constitutive model from first principles, whole-joint modeling efforts will be limited by the need for calibration using large amounts of experimental data. For frequency-domain methods, the biggest impact of research currently is in developing new methods to reduce the number of interface degrees of freedom necessary to accurately model the nonlinear response of a jointed system. Node-to-node contact has the potential to exhibit lower model-form error than whole-joint approaches, but some difficulties can arise on the correct selection of the density and distribution of the contact nodes in order to find the optimum compromise between accuracy and computational speed.

Conclusions

A comparative study was completed that assessed a time-domain joint modeling approach and a frequency-domain approach on their ability to reproduce the nonlinear forced response behavior of an experimental Brake-Reuß beam. To enable this comparison, current methods for these two approaches were expanded upon to facilitate model updating against forced response measurements. In this work, the model pertaining to each approach was updated to match the FRF curves around the resonance of the first bending mode of a Brake-Reuß beam, and then those updated models were used to predict the forced harmonic response near the measured resonances of the second and third bending modes. To compare the methods, a criterion was established that examined how well the approaches captured the measured shift in both frequency and response amplitude of the peak resonance points of each mode as the excitation magnitude was increased. In formulating and updating the models, and viewing the results of the dynamic simulations, insights were gained on their advantages and disadvantages as well as their performance relative to one another.

The two approaches performed about equally well considering their achievements and trade-offs. The wholejoint approach seems to have the potential to follow the experimental backbone curves accurately given its proficient utilization of experimental data, but the choices made in modeling have limited the present model to either match the frequency shift or the damping shift. A better understanding of the epistemic uncertainties associated with the whole-joint modeling approach could solve this problem. On the other hand, the MHBM model has achieved a good estimate of both the frequency shift and the damping using predictive modeling techniques that do not rely heavily on experimental data. In addition, the discretization of the contact area with a superimposed rough contact mesh based on the Voronoi polygons -together with a strategy to collect local contact parameters -proved to be a very effective way to increase the computational speed without significantly compromising on the quality of the results.

Both modeling techniques demonstrated that the dynamics occurring in a joint interface play a large role in invoking a nonlinear response. It is not sufficient simply to insert nonlinear constitutive models at the interfaces of a finite element model; attention must also be given towards how the motion of the joint structure excites those elements. For example, within the joint interface of the Brake-Reuß beam model studied in this work, the nonlinear elements in the regions of low-pressure close to the free edges of the joint where the main contributors towards the nonlinear dynamic behavior since the interfacing surfaces undergo relatively large deformations, whereas elements the high pressure areas (near the bolts) did not contribute significantly to the nonlinearity.

Lastly, this work demonstrated that joint modeling approaches with dissimilar formulations can still be compared by tracking the change in the forced resonance response of the parent structure with amplitude. To facilitate this comparison, the joint model formulations developed by previous works were improved upon to address some of their weaknesses, so this work was as much an exploratory study to push the current state-of-the-art of joints modeling as it was a comparative study. Even so, the joint models were produced from sound engineering judgment by the specialists behind each technique, and the achieved agreement between the models and the experimental data is quite satisfying despite the complex mechanics occurring in the real joint.

The modal force then becomes F r (α) = α φT r M φr = α, and the modal displacement is η r (α) = φ T r Mq(α). The set of η r (α), and associated α, forms the initial loading (backbone) curve, on which Masing's rules for path symmetry in cyclic loading are applied to get the forward F1 and reverse F2 curves of the full hysteresis [START_REF] Segalman | A four-parameter Iwan model for lap-type joints[END_REF], as seen in The hysteresis can then be used to determine the natural frequency and damping applicable for that level of excitation. The slope of the secant to the initial curve is used to calculate the natural frequency, ωr (α) = α η r (α) .

(A.

3)

The energy dissipated per cycle Dr (α), calculated from the area between F1 and F2 , is related to the modal damping ratio as in [START_REF] Lacayo | A numerical study on the limitations of modal Iwan models for impulsive excitations[END_REF],

ζr (α) = Dr (α) 2π (ω r (α)η r (α)) 2 + ζ 0r (A.4)

Here, ζ 0r accounts for linear, viscous sources of damping prevalent at very low amplitude levels. Both ωr (α) and ζr (α) are typically plotted versus η r (α) to show how the damping and natural frequency of mode r change with modal displacement amplitude for the model. These amplitude curves would then be compared against those measured from experiments or other truth models.

Data Retrieval of Modal Properties from FRFs

It is difficult to ascribe an accurate equivalent viscous damping ratio for experimental structures that behave nonlinearly, especially when their measured FRFs show noticeable softening behavior as the forcing amplitude increases like those of the Brake-Reuß beam. For this work, the focus for damping was instead shifted to the model; the modal damping ratios that caused the model to follow the peak FRF points measured from experiment were sought.

To retrieve modal natural frequency and damping from the forced response of a finite element model, the quasilinearity assumption is applied to the modal equations of motion as transformed from Eq. ( 9) using Eq. (A.2), giving

η + diag [2ζ r ω r ] η + diag ω 2 r η = ΦT L T Ψ T p. (A.5)
Here the modal damping ratio ζ r and natural frequency ω r for each mode are allowed to vary instantaneously with amplitude in order to account for the influence of the joint forces f nl . In another sense, this can be thought of as fitting a localized, linear model within one cycle of the steady-state, modal response. A harmonic force is applied at input point A having the complex form p A = P A e iΩt , where P A is a complex amplitude, Ω is the driving frequency, and i = √ -1. The steady state of the modal response is assumed to be reached, with the displacement, velocity, and acceleration having the forms η = He iΩt , η = iΩHe iΩt , and η = -Ω 2 He iΩt ,

Figure 2 :

 2 Figure 2: Linear finite element model of the Brake-Reuß beam. Points A and B reference the same locations as those in Fig. 1.

Figure 3 :

 3 Figure 3: Mode shapes for (a) first bending, (b) second bending, and (c) third bending in the unrestrained configuration of the Brake-Reuß beam model.

Figure 4 :

 4 Figure 4: Abaqus calculation results for static contact (a) pressure and (b) gap distance.

Figure 5 :

 5 Figure 5: Experimental setup of a Brake-Reuß beam.

Figure 6 :

 6 Figure 6: Experimental FRFs around the resonance for the first bending mode of the Brake-Reuß beam.

Figure 7 :

 7 Figure 7: Exploded view of an implementation of a whole joint in a finite element model.

Figure 8 :

 8 Figure 8: (a) A parallel-series Iwan element, and (b) the power-law relationship between the energy dissipated per cycle and the force exerted by the Iwan element within the microslip regime.

Figure 9 :

 9 Figure 9: Exploded view of the five-joint configuration used in the whole-joint modeling approach. Joint models are labeled 1 through 5.

3 )Figure 10 :

 310 Figure 10: Comparison between measurement and the whole-joint model of (a) the natural frequency and (b) modal damping ratio for the first bending mode of the Brake-Reuß beam.

  Fig. 10Fig. 13

Fig. 13

 13 Fig. 10Fig. 13

Fig. 12 RMS

 12 Fig. 12 RMS Error of log(Damping) RMS Error of Frequency Shift

Figure 11 :

 11 Figure 11: Plot of RMS error in the natural frequency shift versus the RMS error of the log of damping computed in the Monte Carlo study on the whole-joint model.

3 )Figure 12 : 3 )Figure 13 :

 312313 Figure 12: Comparison between the estimated and measured (a) natural frequency shift and (b) damping ratio for an optimal fit of the frequency shift in the whole-joint model. (a)

5

 5 each have a different initial preload.

Figure 14 :

 14 Figure 14: Node-to-node joint model in the MHBM approach. Coincident nodes are coupled with a 3D contact element (shown on the right).

Figure 15 :

 15 Figure 15: (a) Sum of pressure forces and (b) average gap distance computed among all nodes inside each Voronoi polygon.

Figure 16 :

 16 Figure 16: Comparison of FRFs between experiment and the initial model for the MHBM approach. The FRFs correspond to location B and show the resonance of the first bending mode.

Figure 17 :

 17 Figure 17: (a) Illustration of the columns of nonlinear nodes, and (b) the average gap and (c) the average contact pressure for each column. Nonlinear nodes are indicated with a dot.

Figure 18 :

 18 Figure 18: (a) FRF comparison between experiment and MHBM simulation (model includes all contact nodes in gap region) for all excitation levels, (b) the change in simulation FRF (2 N) as certain columns of contact elements in the gap region are removed, and (c) the variation of the first-bending mode resonant frequency with contact stiffness.

Figure 19 :

 19 Figure 19: (a) FRF comparison the whole-joint model and experiment for different excitation levels in (a) the first bending mode, (b) the second bending mode, and (c) the third bending mode. The simulated FRFs for the third bending mode were terminated at 1146 Hz to save on computational time.

Figure 20 :

 20 Figure 20: (a) FRF comparison the MHBM model and experiment for different excitation levels in (a) the first bending mode, (b) the second bending mode, and (c) the third bending mode.

1 ]Figure 21 :

 121 Figure 21: Comparison between experiment and the modeling approaches of (a,c,e) the shift in the resonance frequency and (b,d,f) the peak inertances in the (a,b) first, (c,d) second, and (e,f) third bending modes of the Brake-Reuß beam.

FigFigure A. 1 :

 1 Figure A.1: Construction of a Masing hysteresis curve for quasi-static modal analysis.

Table 1 :

 1 Linear natural frequencies for the unrestrained and stick configurations of the Brake-Reuß beam.

	Mode	Unrestr. Nat. Freq. [Hz] Stick Nat. Freq. [Hz]
	1st Bending	130	201
	2nd Bending	542	619
	3rd Bending	1073	1232

Table 2 :

 2 Parameters for the whole joint models.

Table 3 :

 3 Percent error in FRF resonant frequency, frequency shift, and peak inertance for the two joint modeling approaches.

	Bending Mode	Force Level [N]	Resonance Frequency Error [%] Whole Joint MHBM	Frequency Shift Error [%] Whole Joint MHBM	Peak Inertance Error [%] Whole Joint MHBM
	1st	0.1	-0.13	0.01	-	-	1.3	6.8
	1st	0.5	0.18	-0.06	-51.0	14.6	-1.5	-7.7
	1st	1.0	0.36	0.13	-52.0	-6.7	0.3	1.1
	1st	2.0	0.46	0.25	-46.9	-18.4	2.9	-3.3
	2nd	0.1	-0.22	-0.70	-	-	-12.7	-25.9
	2nd	0.5	-0.11	-0.72	-67.0	10.0	-1.5	-20.8
	2nd	1.0	-0.06	-0.71	-63.9	-4.8	1.7	-6.6
	2nd	2.0	-0.02	-0.66	-60.1	-11.9	3.7	-1.6
	3rd	0.1	-2.47	-1.19	-	-	64.9	91.2
	3rd	0.5	-2.47	-1.23	-2.5	60.5	54.4	66.4
	3rd	1.0	-2.46	-1.23	-17.2	31.0	41.7	47.6
	3rd	2.0	-2.41	-1.24	-30.0	24.0	22.2	31.5

Acknowledgments

This research was conducted at the Nonlinear Mechanics and Dynamics (NOMAD) Research Institute hosted by Sandia National Laboratories in Albuquerque, NM, and funded through the WSEAT program of the NNSA.

Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energys National Nuclear Security Administration under contract DE-NA-0003525.

The authors would also like to thank Christian Swacek, Pascal Reuß, Rob Kuether, and Tim Truster for their contributions to aspects of this research.

Appendix A.

Quasi-static modal analysis is a technique that allows the analyst to quickly compute how the modal damping and natural frequencies change with amplitude in a finite element model of structure. This technique was introduced by the authors of [START_REF] Lacayo | Updating structural models containing nonlinear Iwan joints using quasi-static modal analysis[END_REF], who made two simplifying assumptions in the groundwork theory established in [START_REF] Festjens | A numerical tool for the design of assembled structures under dynamic loads[END_REF]: namely 1) negligible change in the mode shapes, and 2) extension of the "joint region" to cover the entire structure of interest. With selective excitation of a mode, [START_REF] Lacayo | Updating structural models containing nonlinear Iwan joints using quasi-static modal analysis[END_REF] shows that the estimated amplitude-dependent damping and natural frequencies produced from quasi-static modal analysis agreed remarkably well with those obtained using Hilbert transform methods. The caveat is that this technique assumes that the nonlinear elements satisfy the Masing hypothesis, which the Iwan element satisfies [START_REF] Iwan | A distributed-element model for hysteresis and its steady-state dynamic response[END_REF].

Theory

At the heart of the technique is the construction of a modal hysteresis curve after applying a static load on the structure in the shape of the mode of interest. Under the quasi-static assumption, Eq. ( 9) reduces to

Here, α is a scalar multiplier for the force, and φr is the mass-normalized mode shape of the rth mode as calculated from Eq. ( 6) except using K and M as the stiffness and mass matrices. Equation (A.1) is computed using a Newton method to retrieve the displacements q(α) for a number of force levels. The applied load and displacements are then projected onto the modal coordinates η using the transformation

where H is the vector of complex amplitudes for modal displacement. Substituting these forms into Eq. (A.5) and solving for H yields

where Ψ A refers to the row in the Hurty-Craig-Bampton transformation matrix associated with input point A. The complex physical displacement at the measured point B can then be retrieved through X B (Ω) = Ψ B LΦH(Ω), for Ψ B refers to the row in associated with output point B.

The above equations are utilized to determine a damping ratio that would make the model reproduce the peak points of the FRF measured from the experimental Brake-Reuß beam. Given the inertance FRF a B measured at output point B,

the resonant frequency is retrieved from the 90°phase crossing,

At the rth resonance point, the physical response can be considered as being dominated by the rth mode, giving output response

(A.9) Substituting Eq. (A.9) into Eq. (A.7) and rearranging for ζ r yields

Lastly, the magnitude of modal displacement at the resonance of the rth mode was retrieved by rearranging the absolute value form of

For this work, |a B (ω r )| was obtained from the experimental magnitude FRF measured at force magnitude |P A |. The resonance point of each FRF curve contributed one datum to the set of ω r , ζ r , and associated |H r | for the rth FRF peak. These measured properties were considered the target to which the set of ωr (α), ζr (α), and associated η r (α) obtained from quasi-static modal analysis (QSMA) were fit by tuning the joint model parameters.

After the fitting, the FRF backbone curve was retrieved using the QSMA curves. The model peak inertance was computed using the QSMA damping by rearranging Eq. (A.