## Abstract

While Iwan elements have been used to effectively model the stiffness and energy dissipation in bolted joints, integrating the equations of motion of these elements is fairly expensive since implicit schemes, such as Newmark’s methods, need to be used. This paper presents a method of simulating dynamic systems containing nonlinear Iwan elements that significantly reduce the computation cost by using closed-form expressions for stiffness and damping in the microslip regime and an averaging method for regions of time in which no external force is applied. The proposed algorithm is demonstrated on a single degree-of-freedom (SDOF) system to evaluate the range over which it retains accuracy and the improvement in performance it offers. Although the current implementation is limited to SDOF systems, it can be used to simulate the response of each mode in structures exhibiting weak nonlinearity that can be modeled using the modal Iwan approach. To verify this, the dynamic response of a finite element model of a beam assembly, integrated using the Newmark-$\beta $ method, has been compared with its equivalent modal model integrated using the proposed algorithm. The results show that the algorithm accurately predicts the response in a fraction of the time taken by implicit integration schemes, so long as the modes remain uncoupled and weakly nonlinear.

## 1 Introduction

Mechanical fasteners have long been known to be a source of stiffness and energy dissipation in built-up structures [1,2]. Bolted joints allow slip between contact interfaces, which leads to frictional energy dissipation and changes in stiffness. Both stiffness and energy dissipation (and hence damping) have been found to show amplitude-dependent behavior [3], that is the amount of slip that occurs is dependent on the amplitude of the force applied. At lower amplitudes, the edges of the joint surfaces slip while a majority of the joint remains clamped due to the bolt pre-load. This is known as microslip. In the microslip regime, the stiffness of the joint decreases only slightly with increase in vibration amplitude, but there is significant energy loss which leads to a large increase in damping. This behavior has been observed in numerous experimental studies [4–6]. As the force amplitude increases, the slip region gradually expands until macroslip occurs, which is characterized by relative motion between the surfaces and significant decrease in joint stiffness. Joints are typically designed to maintain their integrity and are hence expected to predominantly show microslip behavior in most realistic structures. In fact, Deaner et al. [5] experimentally observed that very high force levels are required to cause a bolted structure to go into macroslip, making it difficult to even fully characterize this behavior, unless the bolt torques were unrealistically low.

Modeling the above-explained phenomena in detail would make dynamic analysis highly computationally expensive. An alternative is to replace the contact interface with a lumped, hysteretic model capable of simulating the microslip and macroslip behavior observed. One such model is the Iwan model [7], initially used for metal elasto-plasticity. The Iwan model consists of a parallel system of spring-slider units known as Jenkins elements. There have been many adaptations of the Iwan model to capture joint mechanics. The most widely used among them is Segalman’s four-parameter Iwan model [8], with the four parameters accounting for the joint stiffness, the macroslip force, the transition to macroslip, and the power law energy dissipation that many joints have been found to exhibit in microslip. While this is less expensive than modeling the contact in detail, the computational burden is significant in structures with multiple joints or when performing parametric studies. As an alternative, Segalman [9] proposed a modal approach stating that each nonlinear mode can be represented as an single degree-of-freedom (SDOF) system with an Iwan element to account for the nonlinearity, provided the modes are uncoupled and weakly nonlinear. Deaner et al. [5] and Roettgen and Allen [6] showed that the modal Iwan model is capable of describing the nonlinearity in various bolted structures. Further, Lacayo et al. [10] found that this approach can accurately capture the nonlinear response to an impulse-type excitation (that excites different modes of the structure to varying extents) if it consists of one dominant mode and may still be fairly acceptable in cases of more than one dominant mode, provided the mode being studied is a dominant one.

While these developments have helped speed up nonlinear dynamic response analysis, one bottleneck that remains is the time integration. Currently, implicit numerical integration techniques such as Newmark’s methods need to be used with Newton–Raphson iteration schemes to account for the nonlinear force in the Iwan model. These methods are robust and effective but require a small time-step, making the integration computationally expensive. In an effort to reduce integration costs, Brake [11] presented a reduced formulation of the Iwan model that makes it possible to derive an analytical expression for the nonlinear Iwan force. This paper presents another alternative that exploits the weakly nonlinear behavior of bolted joints in the microslip regime to significantly speed up the integration. Using the modal Iwan approach [9], the response of the uncoupled modes can be computed using closed-form expressions for the energy dissipation and joint stiffness, applicable in the microslip regime. Additionally, the averaging method can be used, taking advantage of the fact that the amplitude and phase of the decaying response vary slowly with time in comparison to the response itself [12]. The simulation method presented in this paper also provides the benefit of calculating the amplitude-dependent damping and natural frequency of the system in the course of the time integration, without requiring further post-processing (like the Hilbert transform [13]). Krack et al. [14] presented a similar approach to compute the steady-state and unsteady dynamics of nonlinear modes in the absence of modal interactions. They used a multi-harmonic analysis to numerically compute the amplitude-dependent characteristics of the nonlinear mode of interest, applying the method of averaging to compute the slow dynamics of the system. The present paper proposes a quasi-linear approach in which the response is assumed to be monoharmonic, and the amplitude-dependent damping and natural frequency are calculated using analytical expressions specifically applicable to Iwan elements.

The following section of the paper gives an overview of the four-parameter Iwan model, highlights the drawbacks of the Newmark-$\beta $ method [15] and explains the theory behind the alternative integration algorithm proposed. This is followed by examining its applicability to an SDOF system with three types of input force—an impulse, a sine beat, and a bandlimited random input. Next, to test the algorithm on a more realistic structure, a test case of a finite element (FE) model of two beams bolted together, commonly referred to as the Sumali beam, is presented. The response obtained using the proposed algorithm with a modal Iwan model is compared with that obtained by numerically integrating the finite element model. The algorithm is shown to be fairly accurate and much faster than the Newmark-$\beta $ method.

## 2 Theoretical Background

### 2.1 Overview of the Four-Parameter Iwan Model.

*x*(

*t*) is the imposed displacement,

*u*(

*t*) is the displacement of the Jenkins elements that constitute the Iwan joint, and $\varphi $ is the displacement at which $\rho (\varphi )$ number of sliders slip. Thus, $\rho (\varphi )$ can be understood as the population density of the sliders. Segalman [8] proposed a power-law population distribution, given by Eq. (2)

*R*and

*S*can be understood as the stiffness of the power-law portion of the distribution and the delta function portion of the distribution, respectively. Since the parameters

*R*and

*S*have fractional dimension, Segalman proposed using another set of more intuitive parameters, $[FS,KT,\chi ,\beta ]$, with

*F*

_{S}being the force required to cause macroslip,

*K*

_{T}being the tangential stiffness of the joint at small applied loads, and $\chi $ and $\beta $ being dimensionless parameters that measure the energy dissipation. It must be noted here that the parameter $\chi $ is identical in both parameter sets. Further details about the model formulation and conversion from one set of parameters to another can be found in Ref. [3].

### 2.2 Numerical Integration Using the Newmark-$\beta $ Method.

*C*

_{lin}is the linear damping coefficient,

*K*

_{∞}is the stiffness of the system when the response amplitude is high enough to cause macroslip in the Iwan element, and

*F*

_{nl}(

*t*) is the nonlinear force due to the Iwan element, given by Eq. (1). As seen in Eq. (1), the distribution of Jenkins elements is theoretically continuous. However, to numerically evaluate the same, the equation is discretized, with 100 sliders typically being a good approximation. The discretization procedure has been explained in Ref. [8].

*n*+ 1) represents the current time-step,

*n*represents the previous time-step and $\u25b3t$ is the step size. This is an implicit technique, i.e., the state variables,

*x*and $x\u02d9$, depend not only on historical information but also on the current estimate of the variable $x\xa8$. Hence, a Newton–Raphson iteration scheme is used to converge on a solution for each time-step. The residue function, $\eta NR$, for the Newton–Raphson method, given by Eq. (7), is obtained by substituting Eqs. (4)–(6) in Eq. (3).

*K*

_{nl}, which depend on the current displacement and state of the sliders, are calculated for every iteration. Since $x\u02d9n+1$ and

*x*

_{n+1}in Eq. (7) depend on $x\xa8n+1$, the gradient, $\theta NR$, obtained by taking the partial derivative of Eq. (7) with respect to $x\xa8n+1$, is given by

### 2.3 Proposed Alternative: Using Closed-Form Expressions for the Nonlinear Parameters.

Even though the Iwan element is complicated and can be rigorously modeled as given above, in the microslip regime, the weakly nonlinear behavior of the system can be taken advantage of. The damping and stiffness change slowly with response amplitude in a power-law fashion, as observed experimentally [16,17].

*X*.

*D*(

*X*), and the joint stiffness,

*K*

_{j}(

*X*), (i.e., Eqs. (11) and (12)) that are applicable in the microslip regime, are used,

*K*

_{∞}is the linear, macroslip stiffness.

What remains is a nonlinear, second-order ordinary differential equation (ODE) with inhomogeneity due to the forcing function, which can be integrated using the fourth-order adaptive Runge–Kutta (RK) method. This requires that the ODE be defined in the form $x\u02d9=f(x,t)$ with the initial conditions provided as input. It is important to note that the traditional Iwan element cannot be integrated using Runge–Kutta, because it is hysteretic by definition and hence cannot be written in the form given by Eq. (10) (without some kind of approximation). For that reason, the Newmark-$\beta $ algorithm with a fixed time-step has been used in most cases in the literature. The following subsections discuss two ways to perform this integration, the difference between them being the method used to obtain the response amplitude.

#### 2.3.1 Integrating the Inhomogeneous Ordinary Differential Equation.

*x*(

*t*), and velocity, $x\u02d9(t)$, as state variables if the amplitude

*X*can be expressed as a function of

*x*(

*t*) and $x\u02d9(t)$. Assuming the response to be harmonic at each time instant, the displacement can be given by Eq. (15). It must be noted here that the assumption of a single harmonic response may not be applicable for every type of excitation but is reasonable for weakly nonlinear vibratory systems.

*X*, can then be estimated as

*X*

_{j}corresponds to the amplitude calculated in the previous loop and

*X*

_{j+1}is the amplitude calculated in the current loop. A convergence tolerance of 10

^{−13}was found to be sufficient for the cases presented here.

In this way, all the terms required to perform the integration of the inhomogeneous ODE can be obtained.

#### 2.3.2 Integrating the Homogeneous Ordinary Differential Equation.

The two integration methods explained in the above subsections can be combined to obtain the ring-down response of an SDOF system with an Iwan element. The forced response (inhomogeneous ODE) is obtained using the technique described in Sec. 2.3.1. The amplitude and phase at the end of this integration are calculated and provided as initial conditions to the next integration method, where the free response (homogeneous ODE) is obtained using the method of averaging described in Sec. 2.3.2. Moving forward, the algorithm obtained by combining these two techniques has been referred to as the averaging algorithm. It must be noted, however, that the method of averaging is only applied in the homogeneous ODE integration and not throughout the numerical simulation.

## 3 Test Case 1: SDOF System

The previous described algorithm was tested against the Newmark-$\beta $ method for speed and accuracy. First, an impulsive force of varying amplitude was applied (Sec. 3.1). To further verify the results, it was also tested with a sine beat input of varying bandwidth (Sec. 3.2), and a multi-amplitude, bandlimited random input (Sec. 3.3). The linear and nonlinear system parameters are given in Table 1. These parameters result in a frequency at macroslip, *f*_{n,∞} of 30 Hz, and a stuck natural frequency (frequency at low-vibration amplitude), *f*_{n,0} of 50 Hz, and are representative of the nonlinear response of a typical mode of a structure with bolted joints.

Parameter | Value |
---|---|

Mass (m) | 1 kg |

Macroslip stiffness (K_{∞}) | 3.55 × 10^{4} N/m |

Material damping ($\zeta lin$) | 1 × 10^{−4} |

Iwan joint $[Fs,KT,\chi ,\beta ]$ | [100 N, 6.32 × 10^{4} N/m, − 0.75, 5] |

Parameter | Value |
---|---|

Mass (m) | 1 kg |

Macroslip stiffness (K_{∞}) | 3.55 × 10^{4} N/m |

Material damping ($\zeta lin$) | 1 × 10^{−4} |

Iwan joint $[Fs,KT,\chi ,\beta ]$ | [100 N, 6.32 × 10^{4} N/m, − 0.75, 5] |

### 3.1 Input: Impulse.

For the first input case, the force applied to the SDOF system was one half cycle of a sinusoid with width of 0.02 s. This was applied at various amplitudes, gradually increasing the amplitude such that the Iwan element progressed from microslip to macroslip, which was verified by checking the displacement of the sliders in the Iwan joint. The time responses for a simulation period of 35 s were obtained by each method. These were then processed using the Hilbert transform, as described in Ref. [13], to estimate the instantaneous damping and natural frequency.

The system response can be divided into two parts—the initial period (when the external force is present) and the ring-down period (when the external force is zero). The accuracy in capturing the ring-down behavior can be measured in terms of the amplitude-dependent damping and natural frequency obtained using the Hilbert transform, as shown in Sec. 3.1.1, while the initial response accuracy can be measured in terms of the displacement amplitude just after the impulse ends (i.e., the response amplitude at *t* = *T*_{pulse}), as shown in Sec. 3.1.2.

#### 3.1.1 Accuracy: Transient Behavior.

To test how accurately the Averaging algorithm predicts the transient behavior, the response of the SDOF system was simulated for an impulse of amplitude 50 N using both integration schemes. This force amplitude is low enough for the Iwan element to remain in the microslip regime. The Hilbert transform of the simulated time response was computed and trimmed to mitigate the end effects. To maintain consistency, the same trim points were used for both simulation methods. A piece-wise linear function was fit to the Hilbert transform to minimize noise when computing its derivative. The instantaneous damping and natural frequency obtained from the Hilbert transform were then plotted against the response velocity amplitude, resulting in the plots shown in Fig. 4. To generalize the results, the non-dimensional velocity was defined as the ratio of the velocity to the product ($\omega n,0\varphi max$), where $\omega n,0$ is the stuck natural frequency and $\varphi max$ is the displacement of the joint at macroslip. The error in damping ratio and natural frequency was then calculated as a percentage of the Newmark-$\beta $ solution, which was considered to be the true solution.

It was initially observed that the damping ratio (and natural frequency) estimated by the Newmark-$\beta $ solution deviated from the averaging algorithm at low response amplitudes, as shown in Fig. 3. The Newmark-$\beta $ solution predicted a damping lower than the linear damping of the system, which is unrealistic. It was then found that the solution predicted by Newmark-$\beta $ algorithm was sensitive to the convergence tolerance of the Newton–Raphson iteration loop,$\u03f5NR$. This was surprising, because the tolerance of 10^{−10} had been used in many previous studies [5,6,10] without difficulty. The tolerance was reduced gradually until further reduction did not have a significant effect on the output. Finally, a tolerance of 10^{−15} was found suitable for the range of forces analyzed in this work. The arbitrary nature of this selection highlights a possible drawback of the existing integration method and further warrants development of a more robust technique.

Figures 4(a) and 4(b) show that the transient behavior predicted by the averaging algorithm agrees well with the Newmark-$\beta $ algorithm after reducing the convergence tolerance for the Newton–Raphson iteration loop. It must be noted here that these graphs progress from right to left (i.e., from higher velocity amplitudes to lower velocity amplitudes) since the response velocity decreases with time. A maximum error of $1.86%$ in the damping estimate and $0.07%$ in the natural frequency estimate is observed. The estimations improve as the velocity amplitude reduces with the percentage error being nearly zero at low velocities.

The averaging algorithm also estimates the amplitude-dependent natural frequency and damping ratio for the free decay using the closed-form expressions, given by Eqs. (11)–(14), without requiring any Hilbert processing. The states used by the Newmark-$\beta $ method, on the other hand, are the displacement *x* and velocity, $x\u02d9$, and so the Hilbert transform is required to estimate the response amplitude, *X*, and the instantaneous natural frequency and damping. Figure 5 shows that the amplitude-dependent damping curve obtained from the closed-form expressions agrees well with the one obtained using the Hilbert transform applied to the Newmark-$\beta $ method. The same result could also be drawn from the amplitude-dependent frequency curve, which has not been shown here for brevity. Thus, the averaging algorithm requires only algebraic processing in the form of evaluating Eqs. (11)–(12) to obtain the amplitude-dependent frequency and damping as opposed to Hilbert processing required in the conventional Newmark-$\beta $ algorithm, which is discussed in Ref. [13].

#### 3.1.2 Accuracy: Initial Response.

Figure 6 shows the time response from both integration schemes for low-amplitude impulsive forces. Both the amplitude and the phase of the response predicted by the averaging algorithm are in good agreement with the Newmark-$\beta $ method. However, at higher force amplitudes, when the system is well into macroslip, the averaging solution starts to deviate from the true response, as can be seen in Fig. 7. This is likely due to the fact that Eqs. (11)–(14) are applicable solely in the microslip regime. Another source of error is the approximate amplitude equation(i.e., Eq. (18)) that assumes small changes in amplitude with time. This approximation could fail at high impulsive forces. These inaccuracies, however, only affect the forced part of the response. If the correct amplitude and phase are given to the homogeneous ODE solver in Sec. 2.3.2, and the response is well approximated as harmonic, it will give the correct response from that point forward.

*E*

_{r}=

*E*

_{impulse}/

*PE*

_{max}, where the numerator is the energy imparted by the impulse (Eq. (29)) and the denominator is the approximate potential energy for which the joint fully slips (Eq. (30)). Macroslip occurs when

*E*

_{impulse}≈

*PE*

_{max}, i.e., when

*E*

_{r}≈ 1.

The response amplitude just after the impulse has ended was then plotted against this energy ratio *E*_{r}. Figure 8 shows that for *E*_{r} < 1, there is negligible difference in the velocity amplitude predicted by the two integration schemes. As *E*_{r} increases to values above 9.3, significant deviations in the amplitude predicted by the two methods are observed. We would expect the averaging algorithm to lose accuracy in the macroslip regime (*E*_{r} > 1) since the closed-form expressions used in this algorithm are applicable only in the microslip regime. Surprisingly, however, the averaging algorithm does not significantly deviate from the truth solution until *E*_{r} reaches values nearly ten times greater than the macroslip value. It must be noted here that one could limit the applicability of the algorithm to the microslip regime to avoid the possibility of inaccurate response estimation at very high force values. This can be done by enforcing *E*_{r} < 1, or equivalently, $X/\varphi max<1$. However, in the case studies shown here the algorithm was actually quite accurate up until *E*_{r} = 10.

#### 3.1.3 Speed.

Table 2 provides the simulation time, in seconds, for each of the two methods when computed on an Intel^{®} Core™ i7-950 (3.07 GHz) processor for the case where an impulsive force having an amplitude of 50 N is applied on the SDOF system under consideration.

Method | No. of time-steps | Simulation time (s) |
---|---|---|

Newmark-$\beta $ | 35001 | 72.96 |

Averaging | 1092 + 241 = 1333 | 0.1193 |

Method | No. of time-steps | Simulation time (s) |
---|---|---|

Newmark-$\beta $ | 35001 | 72.96 |

Averaging | 1092 + 241 = 1333 | 0.1193 |

It can be seen that the averaging algorithm significantly speeds up the integration process. This was anticipated as applying the concept of averaging allows for a larger time-step, thus lowering the number of integration steps. 35,001 time-steps were required for the Newmark-$\beta $ approach (200 samples per period) whereas for the same duration of interest, with the method of averaging, the RK integrator used 1092 time-steps to solve the inhomogeneous ODE and another 241 time-steps to solve the homogeneous part, resulting in a total of 1333 time-steps. The homogeneous ODE requires a much smaller number of time-steps than the inhomogeneous section, due to the use of the averaging method. It must also be noted that the time taken by Newmark-$\beta $ scheme depends on the number of sliders used in the integration. The averaging algorithm has no such issue since it uses closed-form expressions instead of calculating the slider positions.

Table 3 shows the time taken to evaluate the nonlinear force (discretized form of Eq. (1)) in the Newmark-$\beta $ algorithm for 30 versus 100 sliders. The times listed are that for a single evaluation, obtained by performing an average over 500,000 calculations to reduce the effect of variability in computation time. A larger number of sliders are required for higher accuracy but also significantly increase the computational burden. On the other hand, the time taken to evaluate the closed-form expressions for amplitude-dependant damping and frequency (Eqs. (11)–(14)) used in the averaging algorithm is orders of magnitude lower.

### 3.2 Input: Sine Beat.

For a sine beat, the force input is created by applying a Blackman–Harris window to a monoharmonic sine wave, resulting in a signal that starts with zero amplitude, increases in amplitude to a peak value, dwells at this peak value and then decreases in amplitude back to zero. The frequency, *f*_{n}, and peak amplitude, *F*, of the sine wave, ramp time (or frequency band, *dF*) and dwell time, *t*_{d}, can be varied. A wider frequency band (or higher value of *dF*) results in a shorter force signal and can be interpreted as getting “closer” to an impulse. For the purpose of this analysis, *dF* was varied to excite gradually narrowing frequency bands, keeping all other parameters constant. A narrower frequency band (or larger ramp time) means the force is applied over a longer duration, as seen in Fig. 9. Ideally, the sine beat would focus on just the resonant frequency. However, for the nonlinear SDOF system being studied here, the effective natural frequency decreases as amplitude increases, requiring a bandwidth wide enough to account for this variation. The parameters of the input force are provided in Table 4.

Parameter | Value |
---|---|

Frequency (f_{n}) | 50 Hz |

Peak amplitude (F) | 35 N |

Dwell time (t_{d}) | 0 s |

Bandwidth (dF) | varies |

Parameter | Value |
---|---|

Frequency (f_{n}) | 50 Hz |

Peak amplitude (F) | 35 N |

Dwell time (t_{d}) | 0 s |

Bandwidth (dF) | varies |

The response was simulated using the two integration schemes for different ramp times and the amplitude of the response at time greater than *T*_{ramp} was tracked. Figure 10 plots the response amplitude for different ramp times on the left *Y*-axis, and the corresponding maximum slider displacement of the last slider, obtained from the Newmark-$\beta $ method, on the right *Y*-axis. If this slider’s displacement is greater than zero, the system has gone into macroslip. It can be observed that within the microslip regime, the Averaging algorithm accurately predicts the response irrespective of the ramp time (or how long the external force lasts). It is only when the system is well into macroslip that significant deviations are observed. This shows that the algorithm can be used for more persistent forces than an impulse, provided the system remains in microslip.

### 3.3 Input: Bandlimited Random Force.

To test the algorithm’s efficiency in simulating the response of randomly excited systems, a third input type was considered. A uniformly distributed random force was generated, which was then filtered to focus on the frequency range of interest. This was done using a low-pass Butterworth filter with a cutoff frequency equal to 1.25 times the stuck frequency (*f*_{n,0}). Three different peak forcing amplitudes were considered, resulting in increasing levels of nonlinearity of the system. The response for each forcing level was computed over a time interval of 3200 s using the Newmark-$\beta $ and the averaging algorithms, respectively, and the power spectral density (PSD) of the response was calculated.

Figure 11(a) compares the PSDs obtained for the three different force levels using the two algorithms. The response is nonlinear at the higher force levels, with a frequency shift of nearly 1.5 Hz for the force level of ±150 N. The closely matching PSDs obtained from the two integration approaches indicate that the averaging algorithm is able to accurately simulate the response of a nonlinear system to a random input. The time-domain displacement, plotted in Fig. 11(b), shows that the response is sinusoidal in nature at each time instant, even though the applied input is random. This explains why the averaging algorithm is still able to predict the response fairly accurately. Table 5 lists the time takenby the Newmark-$\beta $ and the averaging algorithms to compute the response for the three different force amplitudes considered. It must be noted that since the external force is always active, the ODE to be integrated is inhomogeneous throughout the time interval and hence, only the part of the averaging algorithm derived in Sec. 2.3.1 is being implemented. The evaluation times show that the computation time can be reduced by more than an order of magnitude by implementing the amplitude approximation equation (Eq. (18)) along with the recursive scheme described in Sec. 2.3.1, even without the method of averaging being used.

## 4 Test Case 2: Sumali Beam Finite Element Model

To test the algorithm’s accuracy on a more realistic structure, an assembly of two beams commonly referred to as the Sumali beam, was considered. The Sumali beam structure consists of two thin, identical, stainless steel beams of length 508.00 mm, width 50.80 mm, and thickness 6.35 mm, that overlap and are connected using four bolts. The area of overlap extends 355.60 mm along the length of each beam.

Lacayo et al. [10] performed model correlation and updating on a FE model of this beam. The Hurty/Craig-Bampton reduction technique [19,20] was used to obtain a reduced FE model, consisting of 49-degrees-of-freedom (DOF). Each of the four bolts was modeled using the whole-joint approach [3], i.e., they were represented by a single Iwan element in the *x*-translation direction (axial to the beam). All the Iwan elements were assumed to have the same parameters which were iterated upon to match the experimental results obtained by Deaner et al. [5]. The structure was modeled to have free-free boundary conditions. Hence, the first six modes are rigid body modes. Modes 7, 8, and 9 are the first three elastic bending modes, as shown in Fig. 12. These modes exhibit nonlinearity. Further information on the FE model and the updating routine used can be found in Ref. [10]. Apart from this, [10] also implements the modal Iwan modeling approach [9], representing each of the nonlinear modes of the Sumali beam as a SDOF system containing an Iwan element. To do so, first an impulsive force was applied on the FE model to excite only the mode of interest. The response of the 49DOFs to this force was simulated using the Newmark-$\beta $ integration method. The simulated response was then modally filtered, and Hilbert analysis was performed on the velocity of the nonlinear mode of interest (*r*) to obtain its amplitude-dependent natural frequency,$\omega r,meas$, and damping ratio, $\zeta r,meas$. The closed-form expressions for damping and frequency (Eqs. (11)–(14)) were then used, treating the Iwan parameters as variables that were optimized to match $\omega r,meas$ and $\zeta r,meas$. Table 6 lists the Iwan parameters calculated in Ref. [10] for the three nonlinear modes of the Sumali beam.

Parameter | Mode 7 | Mode 8 | Mode 9 |
---|---|---|---|

K_{∞} [s^{−2}] | 4.235 × 10^{5} | 1.673 × 10^{6} | 7.321 × 10^{6} |

$\zeta lin$ | 1 × 10^{−4} | 1 × 10^{−4} | 1 × 10^{−4} |

F_{s} [kg^{1/2} ms^{−2}] | 1.802 | 1.518 | 6.934 |

K_{T} [s^{−2}] | 1.537 × 10^{5} | 1.340 × 10^{5} | 2.822 × 10^{6} |

$\chi $ | −0.1369 | −0.1620 | −0.1008 |

$\beta $ | 2.561 | 2.000 | 1.529 |

Parameter | Mode 7 | Mode 8 | Mode 9 |
---|---|---|---|

K_{∞} [s^{−2}] | 4.235 × 10^{5} | 1.673 × 10^{6} | 7.321 × 10^{6} |

$\zeta lin$ | 1 × 10^{−4} | 1 × 10^{−4} | 1 × 10^{−4} |

F_{s} [kg^{1/2} ms^{−2}] | 1.802 | 1.518 | 6.934 |

K_{T} [s^{−2}] | 1.537 × 10^{5} | 1.340 × 10^{5} | 2.822 × 10^{6} |

$\chi $ | −0.1369 | −0.1620 | −0.1008 |

$\beta $ | 2.561 | 2.000 | 1.529 |

*Z*-axis at the end of the beam. This force level was chosen to ensure that all the Iwan elements are within the microslip region, and the nonlinearity is weak enough for the modal Iwan modeling approach to hold. The equation of motion for the FE model is given by Eq. (31).

*q*

_{0}is the initial displacement, and $q\u02d90$ is the initial velocity. The initial conditions are equal to the displacement and velocity at the last time-step of the numerical integration solution.

Figure 13 shows the response velocity for mode 8 of the Sumali beam obtained by the two methods—integrating the whole-joint model and using the averaging algorithm. The two time histories are in fairly good agreement. The responses can be further processed by performing the Hilbert transform described in Sec. 3.1.1 to obtain the amplitude-dependent natural frequency and damping, shown in Fig. 14. The modal Iwan model can also be integrated using the Newmark-$\beta $ method, the results of which have also been shown in Fig. 14. The natural frequency predicted by the averaging algorithm is in close agreement with the whole-joint model, with the maximum percentage error being $0.021%$. The damping predicted by the averaging algorithm shows some deviation, with the maximum percentage error being $15.58%$, treating the whole-joint model as the truth solution. However, it can be seen in Fig. 14(b) that the damping curves obtained by integrating the modal Iwan model using the Newmark-$\beta $ method and the averaging algorithm (i.e. the dotted and dashed curves in Fig. 14(b)) overlay. Hence, the difference in damping observed between the whole-joint model and the averaging algorithm is a consequence of the modal Iwan approach and not an integration error caused by the method of averaging.

Figure 15 shows the FFT of the response velocity at the end of the beam. Since modes up to 1050 Hz were used to update the FE model, only the modes in that frequency range were compared. It can be seen that there is good agreement between the two integration methods for modes 7, 8, and 9. At higher frequencies, the whole-joint model deviates from the natural frequency of the mode, possibly due to integration errors introduced by the Newmark-$\beta $ method. The modal approach using the averaging algorithm, on the other hand, aligns well with the system natural frequencies.

Table 7 lists the computation time for the three approaches: integrating the whole-joint model using the Newmark-$\beta $ method, integrating the modal model using the Newmark-$\beta $ method, and integrating the modal model using the averaging algorithm. In the modal approach, the response of the linear modes is obtained by direct numerical integration coupled with evaluation of analytical expressions (i.e., Eqs. (33) and (34)). The table highlights the significant improvement in speed the averaging algorithm offers; it is nearly 300 times faster than integrating the whole-joint model and over 80 times faster than using the Newmark-$\beta $ method for the modal model.

Method | Nonlinear modes | Linear modes | Total time (s) | |
---|---|---|---|---|

(a) | Newmark-$\beta $ whole-joint model | – | – | 997.97 |

(b) | Newmark-$\beta $ modal model | 257.19 s | 2.91 s | 260.10 |

(c) | Averaging modal model | 0.36 s | 2.91 s | 3.27 |

Method | Nonlinear modes | Linear modes | Total time (s) | |
---|---|---|---|---|

(a) | Newmark-$\beta $ whole-joint model | – | – | 997.97 |

(b) | Newmark-$\beta $ modal model | 257.19 s | 2.91 s | 260.10 |

(c) | Averaging modal model | 0.36 s | 2.91 s | 3.27 |

This shows that applying the proposed averaging algorithm to realistic structures can be computationally more efficient than integrating the whole FE model, without significant loss of accuracy.

## 5 Conclusions

This paper presents an algorithm that can be effectively used to integrate a single degree-of-freedom system with an Iwan element in order to simulate a ring-down response. The method of averaging can be used when the external force goes to zero to speed up the integration. The proposed algorithm is much more computationally efficient than the prevalent implicit integration approaches with minimal loss in accuracy, even for impulsive forces large enough to cause the joints to go into macroslip. For example, in the SDOF system considered, accurate response is estimated for energy ratios as high as *E*_{r} = 9.3. The sine beat and random input case studies show that the method of averaging gives accurate results irrespective of how long the input force lasts, as long as the system is in or near microslip. In the microslip regime, the amplitude-dependent natural frequency and damping ratio were found to be comparable to the Newmark-$\beta $ algorithm. In fact, the Newmark-$\beta $ algorithm required adjusting the ad-hoc convergence tolerance in the Newton–Raphson iteration loop to accurately capture the damping and frequency at low-vibration amplitudes, making the Averaging solution more reliable in the cases considered. The Averaging algorithm is also preferable if the final goal is to estimate damping and frequency versus time, because no further processing (i.e., with the Hilbert transform) is required to obtain these. A drawback of the integration technique is that the closed-form expressions used are fundamentally limited to microslip only, even though results show that the approach is fairly accurate even for forces larger than macroslip. Another limitation is that this algorithm is applicable to SDOF systems or modal Iwan models; the latter are only valid for MDOF systems if the modes are uncoupled and the nonlinearity is weak. The Sumali beam test case shows that if the modal Iwan approach holds, the averaging algorithm offers great computational efficiency without much loss of accuracy.

## Acknowledgment

This material is based in part upon work supported by the National Science Foundation under Grant No. CMMI-1561810. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.