Fast Simulation of a Single Degree-of-Freedom System Consisting of An Iwan Element Using the Method of Averaging

Shetty, Drithi; Allen, Matthew

doi:10.1115/1.4047070

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- $β$ 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- $β$ 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- $β$ method.

2 Theoretical Background

2.1 Overview of the Four-Parameter Iwan Model.

The Iwan model [7] consists of a parallel arrangement of springs and sliders connected in series, known as Jenkins elements. The constitutive form of the model is given by Eq. ,

F_{nl} (t) = \int_{0}^{\infty} ρ (ϕ) [x (t) - u (t, ϕ)] d ϕ

(1)

where x(t) is the imposed displacement, u(t) is the displacement of the Jenkins elements that constitute the Iwan joint, and

ϕ

is the displacement at which

ρ (ϕ)

number of sliders slip. Thus,

ρ (ϕ)

can be understood as the population density of the sliders. Segalman [8] proposed a power-law population distribution, given by Eq.

ρ (ϕ) = R ϕ^{χ} [H (ϕ) - H (ϕ - ϕ_{max})] + S δ (ϕ - ϕ_{max})

(2)

Segalman’s four-parameter model can thus be represented by the parameter set

[ϕ_{max}, χ, R, S]

⁠,where

ϕ_{max}

is the displacement at which all sliders slip (i.e., macroslip occurs),

χ

is a dimensionless quantity that measures the power-law energy dissipation, and 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,

[F_{S}, K_{T}, χ, β]

⁠, with F_S being the force required to cause macroslip, K_T being the tangential stiffness of the joint at small applied loads, and

χ

and

β

being dimensionless parameters that measure the energy dissipation. It must be noted here that the parameter

χ

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- $β$ Method.

An SDOF system with a single Iwan element has an equation of motion given by Eq. ,

m \ddot{x} + C_{lin} \dot{x} + K_{\infty} x + F_{nl} (t) = F_{ext} (t)

(3)

where 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. . As seen in Eq. , 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].

Equation can be integrated using the Newmark-

β

method derived in Ref. [15]. In this method, at every time-step, the initial guess for acceleration equals the acceleration corresponding to the previous time-step, as shown in Eq. . An initial estimate for velocity and displacement can then be obtained using Eqs. and ,

{\ddot{x}}_{n + 1} = {\ddot{x}}_{n}

(4)

{\dot{x}}_{n + 1} = {\dot{x}}_{n} + △ t [γ_{nb} {\ddot{x}}_{n + 1} + (1 - γ_{nb}) {\ddot{x}}_{n}]

(5)

x_{n + 1} = x_{n} + △ t {\dot{x}}_{n} + \frac{△ t^{2}}{2} [2 β_{nb} {\ddot{x}}_{n + 1} + (1 - 2 β_{nb}) {\ddot{x}}_{n}]

(6)

where (n + 1) represents the current time-step, n represents the previous time-step and

△ t

is the step size. This is an implicit technique, i.e., the state variables, x and

\dot{x}

⁠, depend not only on historical information but also on the current estimate of the variable

\ddot{x}

⁠. Hence, a Newton–Raphson iteration scheme is used to converge on a solution for each time-step. The residue function,

η_{NR}

⁠, for the Newton–Raphson method, given by Eq. , is obtained by substituting Eqs. – in Eq. .

\begin{aligned} η_{NR} = m {\ddot{x}}_{n + 1} + C_{lin} {\dot{x}}_{n + 1} + K_{\infty} x_{n + 1} + F_{nl} (x_{n + 1}, ϕ_{n + 1}) - F_{ext} (t_{n + 1}) \end{aligned}

(7)

The nonlinear joint force and the corresponding nonlinear tangential stiffness K_nl, which depend on the current displacement and state of the sliders, are calculated for every iteration. Since

{\dot{x}}_{n + 1}

and x_n+1 in Eq. depend on

{\ddot{x}}_{n + 1}

⁠, the gradient,

θ_{NR}

⁠, obtained by taking the partial derivative of Eq. with respect to

{\ddot{x}}_{n + 1}

⁠, is given by

θ_{NR} = m + C_{lin} γ_{nb} △ t + [K_{\infty} + K_{nl}] β_{nb} (△ t)^{2}

(8)

The variable

{\ddot{x}}_{n + 1}

is iterated upon until the convergence tolerance

ϵ_{NR}

⁠, given by Eq. , is below a predefined value. A very low convergence tolerance of

ϵ_{NR} = 10^{- 15}

was chosen for the presented studies, the reason behind which will be discussed in Sec. 3.1.1. Once the Newton–Raphson algorithm converges and the corresponding velocity and displacement for a particular time-step are calculated using Eqs. and , the algorithm marches forward in time.

ϵ_{NR} = \frac{η_{NR}^{2}}{| C_{lin} {\dot{x}}_{n + 1} + K_{\infty} x_{n + 1} |^{2}}

(9)

Either a constant or linear acceleration can be assumed between time-steps. The value of

β_{nb}

and

γ_{nb}

depends on the acceleration method used. In the constant acceleration method,

β_{nb} = 1 / 4

and

γ_{nb} = 1 / 2

⁠, whereas in the linear acceleration method,

β_{nb} = 1 / 6

and

γ_{nb} = 1 / 2

⁠. Since the constant acceleration method assumes that the acceleration remains constant within an integration step, a small step size is required to obtain an accurate solution [15]. A step size 200 times smaller than the shortest period was required in prior works [5,6]. Thus, this method can be time-consuming and warrants a more efficient approach.

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].

The alternative approach presented here is to modify the equation of motion to be of the form as in Eq. , where the natural frequency,

ω_{n}

⁠, and the damping,

ζ

⁠, of the system are expressed as functions of the amplitude of the response, referred to as X.

\ddot{x} + 2 ζ (X) ω_{n} (X) \dot{x} + ω_{n} {(X)}^{2} x = \frac{F_{ext}}{m}

(10)

To solve this modified equation of motion, closed-form expressions derived by Segalman [8] for the dissipation per cycle, D(X), and the joint stiffness, K_j(X), (i.e., Eqs. and ) that are applicable in the microslip regime, are used,

D (X) = \frac{4 R X^{χ + 3}}{(χ + 3) (χ + 2)}

(11)

K_{j} (X) = K_{T} [1 - \frac{r^{χ + 1}}{(χ + 2) (β + 1)}]

(12)

where

r = X / ϕ_{max}

⁠. The damping and natural frequency can then be calculated using Eqs. and ,

ζ (X) = ζ_{lin} + \frac{D (X)}{2 π m [ω_{n} (X)]^{2} X^{2}}

(13)

ω_{n} (X) = \sqrt{\frac{K_{j} (X) + K_{\infty}}{m}}

(14)

where

ζ_{lin}

is the linear damping ratio and 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 $\dot{x} = 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- $β$ 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.

The inhomogeneous ODE can be integrated with displacement, x(t), and velocity,

\dot{x} (t)

⁠, as state variables if the amplitude X can be expressed as a function of x(t) and

\dot{x} (t)

⁠. Assuming the response to be harmonic at each time instant, the displacement can be given by Eq. . 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 (t) = X \sin (ω_{d} (X) t + ϕ_{0})

(15)

Differentiating Eq. with respect to time gives

\begin{aligned} \dot{x} (t) & = [ω_{d} (X) + t \frac{d ω_{d} (X)}{d X} \dot{X}] X \cos [ω_{d} (X) t + ϕ_{0}] \\ + \dot{X} \sin [ω_{d} (X) t + ϕ_{0}] \end{aligned}

(16)

If the change in amplitude with time is small,

\dot{X}

is negligible, and Eq. can be simplified to

\dot{x} (t) \approx X ω_{d} (X) \cos (ω_{d} (X) t + ϕ_{0})

(17)

The amplitude, X, can then be estimated as

X \approx \sqrt{x {(t)}^{2} + {[\frac{\dot{x} (t)}{ω_{d} (X)}]}^{2}}

(18)

where

ω_{d} (X) = ω_{n} (X) \sqrt{1 - ζ {(X)}^{2}}

⁠. It can be seen that the amplitude, natural frequency, and damping are interdependent variables. Hence, an iterative scheme, shown in Fig. 1, needs to be implemented within the RK integrator function definition. The amplitude and corresponding

ω_{n}

and

ζ

are calculated iteratively until the absolute difference between the amplitude calculated in successive iterations falls within a pre-defined tolerance

ϵ

⁠, given by Eq. ,

ϵ = | X_{j + 1} - X_{j} |

(19)

where 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⁻¹³ was found to be sufficient for the cases presented here.

Fig. 1

View large Download slide

Depiction of the recursive scheme implemented to estimate the response amplitude within the inhomogeneous ODE function definition

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.

In many cases, the external force lasts only for a short period and one is concerned with the free response of the system, which is governed by Eq. .

\ddot{x} + 2 ζ (X) ω_{n} (X) \dot{x} + ω_{n} {(X)}^{2} x = 0

(20)

While the previous method could be used to integrate the homogeneous ODE as well, a quicker solution can be obtained by changing the state variables. Figure 2 shows an example of the transient response of an SDOF system. It can be seen that although the response itself changes quickly with time, the amplitude and phase of the response change relatively slowly with time. Hence, changing the state variables from displacement and velocity to amplitude and phase allows a larger time-step to be used, thus speeding up the integration further.

Fig. 2

View large Download slide

Transient response showing the amplitude and phase changing slowly with time

The equations for change in amplitude,

\dot{X}

⁠, and phase,

\dot{γ}

⁠, can be found by applying the averaging method [12]. This was done, for example, in Ref. [18]. The displacement and velocity are written considering amplitude and phase to be time-dependent rather than constant, as shown in Eqs. and .

x = X (t) e^{i γ (t)}

(21)

\dot{x} = \dot{X} (t) e^{i γ (t)} + X (t) \dot{γ} (t) i e^{i γ (t)}

(22)

The acceleration can be obtained by differentiating Eq. with respect to time. Since bolted joints exhibit weak nonlinearity in the microslip regime, the frequency (and hence change in phase) varies slowly with time. Therefore,

\ddot{γ} (t)

can be neglected and the acceleration equation can be given as

\ddot{x} = [\ddot{X} (t) + 2 \dot{X} (t) \dot{γ} (t) i - X \dot{γ} {(t)}^{2}] e^{i γ (t)}

(23)

The equations for

x (t), \dot{x} (t)

⁠, and

\ddot{x} (t)

are then substituted in the equation of motion, i.e., Eq. . The real and imaginary parts are equated, resulting in Eqs. and .

\dot{X} (t) = - X (t) ω_{n} (X) ζ (X)

(24)

\dot{γ} {(t)}^{2} - \frac{\ddot{X} (t)}{X (t)} = ω_{n} {(X)}^{2} (1 - 2 ζ {(X)}^{2})

(25)

Equation is a first-order ODE for the instantaneous amplitude. Exploiting the weakly nonlinear behavior of bolted joints,

ω_{n} (X)

and

ζ (X)

in Eq. can be considered constant locally. Thus, differentiating Eq. with respect to time results in the equation,

\ddot{X} (t) = - \dot{X} (t) ω_{n} (X) ζ (X)

(26)

Substituting Eq. for

\dot{X} (t)

in Eq. gives,

\ddot{X} (t) = - X (t) ω_{n} {(X)}^{2} ζ {(X)}^{2}

(27)

which can then be substituted in Eq. to obtain an ODE for the instantaneous phase, as given by Eq. . It must be noted here that

ω_{n}

and

ζ

are locally constant but still need to be calculated for different amplitude values using Eqs. –.

\dot{γ} (t) = ω_{n} (X) \sqrt{1 - ζ {(X)}^{2}}

(28)

In this way, the state variables can be changed from displacement and velocity to amplitude and phase. The ODEs for amplitude and phase (Eqs. and ) can then be integrated using the Runge–Kutta integration scheme, and the response displacement and velocity can be retrieved using Eqs. and .

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- $β$ 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.

Table 1

Properties of the nonlinear SDOF system used in Test case 1

Parameter	Value
Mass (m)	1 kg
Macroslip stiffness (K_∞)	3.55 × 10⁴ N/m
Material damping (⁠ $ζ_{lin}$ ⁠)	1 × 10⁻⁴
Iwan joint $[F_{s}, K_{T}, χ, β]$	[100 N, 6.32 × 10⁴ 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 (⁠ $ω_{n, 0} ϕ_{max}$ ⁠), where $ω_{n, 0}$ is the stuck natural frequency and $ϕ_{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- $β$ solution, which was considered to be the true solution.

It was initially observed that the damping ratio (and natural frequency) estimated by the Newmark- $β$ solution deviated from the averaging algorithm at low response amplitudes, as shown in Fig. 3. The Newmark- $β$ 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- $β$ algorithm was sensitive to the convergence tolerance of the Newton–Raphson iteration loop, $ϵ_{NR}$ ⁠. This was surprising, because the tolerance of 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⁻¹⁵ 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.

Fig. 3

View large Download slide

Comparison of amplitude-dependent damping obtained from the Newmark- $β$ and averaging algorithms. The Newmark- $β$ algorithm has a Newton–Raphson iteration tolerance of $ϵ_{NR} = 10^{- 10}$ ⁠.

Figures 4(a) and 4(b) show that the transient behavior predicted by the averaging algorithm agrees well with the Newmark- $β$ 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.

Fig. 4

View large Download slide

Comparison of (a) damping and (b) frequency obtained from the Newmark- $β$ (with Newton–Raphson tolerance of $ϵ_{NR} = 10^{- 15}$ ⁠) and the averaging algorithms. The % error is calculated treating Newmark- $β$ as the true solution: (a) damping (top) and % error in damping (bottom) and (b) natural frequency (top) and % error in frequency (bottom).

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- $β$ method, on the other hand, are the displacement x and velocity, $\dot{x}$ ⁠, 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- $β$ 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- $β$ algorithm, which is discussed in Ref. [13].

Fig. 5

View large Download slide

Comparison of the damping obtained using Eqs. (11)–(14) in the averaging algorithm with the damping obtained using the Hilbert transform method

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- $β$ 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.

Fig. 6

View large Download slide

Response velocity simulated by the Newmark- $β$ and the averaging algorithms in the microslip regime (E_r = 0.05)

Fig. 7

View large Download slide

Response velocity simulated by the Newmark- $β$ and the averaging algorithms in the macroslip regime (E_r = 11.07)

To further quantify this, the response was simulated for multiple force amplitudes of increasing magnitude using both integration methods. To generalize the results, the non-dimensional force amplitude was expressed as an energy ratio E_r = E_impulse/PE_max, where the numerator is the energy imparted by the impulse (Eq. ) and the denominator is the approximate potential energy for which the joint fully slips (Eq. ). Macroslip occurs when E_impulse ≈ PE_max, i.e., when E_r ≈ 1.

E_{impulse} = \frac{1}{2 m} (\int_{0}^{T_{pulse}} f (t) d t)^{2}

(29)

P E_{max} \approx \frac{1}{2} (K_{T} + K_{\infty}) ϕ_{max}^{2}

(30)

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 / ϕ_{max} < 1$ ⁠. However, in the case studies shown here the algorithm was actually quite accurate up until E_r = 10.

Fig. 8

View large Download slide

Response amplitude after the impulse has ended plotted against energy ratio E_r, where E_r < 1 indicates microslip and E_r ≥ 1 indicates macroslip

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.

Table 2

Simulation time of the Newmark- $β$ and the averaging methods when a 50 N impulsive force is applied on the SDOF system

Method	No. of time-steps	Simulation time (s)
Newmark- $β$	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- $β$ 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- $β$ 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- $β$ 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.

Table 3

Time taken to evaluate the discretized, nonlinear force in the Newmark- $β$ algorithm compared against its averaging counterpart, the time taken to evaluate the closed-form damping and frequency expressions

Method	No. of sliders	Function evaluation time (s)
Newmark- $β$	30	7.04 × 10⁻⁵
Newmark- $β$	100	1.04 × 10⁻⁴
Averaging	–	5.98 × 10⁻⁷

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.

Fig. 9

View large Download slide

Force versus time plot for two different ramp times, showing that a larger ramp time implies a longer duration of force

Table 4

Parameters of the sine beat applied to the nonlinear SDOF system

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- $β$ 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.

Fig. 10

View large Download slide

Response amplitude when t > T_ramp plotted against varying ramp time. The maximum slider displacement of the last slider to slip is tracked on the right axis, with non-zero values corresponding to macroslip.

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- $β$ 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- $β$ 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.

Fig. 11

View large Download slide

(a) PSDs for bandlimited random input with three different peak force amplitudes (Am) and (b) time-domain response to random input with force level ±150 N, obtained by the Nemark- $β$ and averaging algorithms

Table 5

Computation time for the two integration algorithms for a bandlimited random input when the peak force amplitude is (a) 1 N, (b) 50 N, and (c) 150 N

	Peak force (N)	Newmark- $β$ evaluation time (s)	Averaging evaluation time (s)
(a)	±1	8908	386
(b)	±50	7221	380
(c)	±150	7223	385

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- $β$ 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, $ω_{r, meas}$ ⁠, and damping ratio, $ζ_{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 $ω_{r, meas}$ and $ζ_{r, meas}$ ⁠. Table 6 lists the Iwan parameters calculated in Ref. [10] for the three nonlinear modes of the Sumali beam.

Fig. 12

View large Download slide

Linear mode shapes of modes 7, 8, and 9, the nonlinear modes of the Sumali beam: (a) first bending mode (mode 7), (b) second bending mode (mode 8), and (c) third bending mode (mode 9)

Table 6

Modal Iwan parameters for the Sumali beam reduced-order model, adapted from Ref. [10]

Parameter	Mode 7	Mode 8	Mode 9
K_∞ [s⁻²]	4.235 × 10⁵	1.673 × 10⁶	7.321 × 10⁶
$ζ_{lin}$	1 × 10⁻⁴	1 × 10⁻⁴	1 × 10⁻⁴
F_s [kg^1/2 ms⁻²]	1.802	1.518	6.934
K_T [s⁻²]	1.537 × 10⁵	1.340 × 10⁵	2.822 × 10⁶
$χ$	−0.1369	−0.1620	−0.1008
$β$	2.561	2.000	1.529

Parameter	Mode 7	Mode 8	Mode 9
K_∞ [s⁻²]	4.235 × 10⁵	1.673 × 10⁶	7.321 × 10⁶
$ζ_{lin}$	1 × 10⁻⁴	1 × 10⁻⁴	1 × 10⁻⁴
F_s [kg^1/2 ms⁻²]	1.802	1.518	6.934
K_T [s⁻²]	1.537 × 10⁵	1.340 × 10⁵	2.822 × 10⁶
$χ$	−0.1369	−0.1620	−0.1008
$β$	2.561	2.000	1.529

The previously-described modal model has been used here to test the applicability of the Averaging algorithm to a realistic structure. An impulsive force of 0.5 N was applied along the positive 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. .

M \ddot{x} + C_{lin} \dot{x} + K_{lin} x + F_{joint} (x) = F_{ext}

(31)

Time integration of the reduced FE model, henceforth referred to as the whole-joint model, was performed using the Newmark-

β

method, and the response obtained in the physical domain was modally filtered to get the displacement and velocity of the nonlinear modes of the system. This was then compared with the response obtained using the modal Iwan model of the beam, integrated using the averaging algorithm. The equation of motion for an multi-degrees-of-freedom (MDOF) system can be converted from physical to modal coordinates by pre-multiplying the equation by the transpose of the linear, mass-normalized mode shape matrix Φ, assuming the mode shapes remain constant and equal to the linear modes of the system. This is a reasonable assumption for the weak nonlinearity exhibited by bolted joints in the microslip regime, with the results obtained in Refs. [5,6] proving the same. The equation of motion for any of the modes of the system can then be written in the form that is compatible with the averaging algorithm, as shown in Eq. .

{\ddot{q}}_{r} + 2 ζ_{r} ω_{r} {\dot{q}}_{r} + ω_{r}^{2} q_{r} = Φ_{0 r}^{T} F_{ext}

(32)

where

ζ_{r}

and

ω_{r}

correspond to the linear damping ratio and natural frequency, respectively, for linear modes and to the amplitude-dependent damping and natural frequency for nonlinear modes.

Each of the three nonlinear modes can then be treated as an SDOF system which can be integrated in time using the averaging algorithm. The remaining 46 modes can be approximated as linear and integrated directly using any numerical integration method. In fact, for the linear modes, once the external force goes to zero, the classical expressions for displacement and velocity for free vibration of an under-damped SDOF system, given by Eqs. and , respectively, and derived in Ref. [21], can be used,

q_{r} = e^{- ζ_{r} ω_{r} t} [q_{0} \cos (ω_{d} t) + \frac{{\dot{q}}_{0} + ζ_{r} ω_{r} q_{0}}{ω_{d}} \sin (ω_{d} t)]

(33)

{\dot{q}}_{r} = e^{- ζ_{r} ω_{r} t} [{\dot{q}}_{0} \cos (ω_{d} t) - \frac{ω_{r}}{ω_{d}} (ζ_{r} {\dot{q}}_{0} + ω_{r} q_{0}) \sin (ω_{d} t)]

(34)

where

ω_{d}

is the damped natural frequency, q₀ is the initial displacement, and

{\dot{q}}_{0}

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- $β$ 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- $β$ 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.

Fig. 13

View large Download slide

Modal response of the nonlinear mode 8, with the modal model, integrated using the Averaging algorithm, being compared against the modally-filtered response from the whole-joint model

Fig. 14

View large Download slide

Comparison of the amplitude-dependent (a) frequency and (b) damping obtained from three integration approaches: integrating the whole-joint model using the Newmark- $β$ method, integrating the modal model using the averaging algorithm, and integrating the modal model using the Newmark- $β$ method: (a) frequency versus velocity amplitude and (b) damping versus velocity amplitude.

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- $β$ method. The modal approach using the averaging algorithm, on the other hand, aligns well with the system natural frequencies.

Fig. 15

View large Download slide

Fast Fourier transform (FFT) of the response velocity at the end of the beam

Table 7 lists the computation time for the three approaches: integrating the whole-joint model using the Newmark- $β$ method, integrating the modal model using the Newmark- $β$ 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- $β$ method for the modal model.

Table 7

Computation time for the three integration approaches: using the Newmark- $β$ method for (a) the whole-joint model, (b) the nonlinear modes in the modal model, and (c) using the averaging algorithm for the modal model

	Method	Nonlinear modes	Linear modes	Total time (s)
(a)	Newmark- $β$ whole-joint model	–	–	997.97
(b)	Newmark- $β$ 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- $β$ algorithm. In fact, the Newmark- $β$ 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.

References

1.

Gaul

,

L.

, and

Lenz

,

J.

,

1997

, “

Nonlinear Dynamics of Structures Assembled by Bolted Joints

,”

Acta Mech.

,

125

(

1

), pp.

169

–

181

. 10.1007/BF01177306

2.

Ungar

,

E.

,

1973

, “

The Status of Engineering Knowledge Concerning the Damping of Built-Up Structures

,”

J. Sound Vib.

,

26

(

1

), pp.

141

–

154

. 10.1016/S0022-460X(73)80210-X

Google Scholar

Crossref

3.

Segalman

,

D. J.

,

2006

, “

Modelling Joint Friction in Structural Dynamics

,”

Struct. Control Health Monitoring

,

13

(

1

), pp.

430

–

453

. 10.1002/stc.119

Google Scholar

Crossref

4.

Smallwood

,

D. O.

,

Gregory

,

D. L.

, and

Coleman

,

R. G.

,

2000

, “

Damping Investigations of a Simplified Frictional Shear Joint

,”

Sandia National Lab. (SNL-NM)

,

Albuquerque, NM (US)

;

Sandia National Labs., Livermore, CA (US), Techical Report No. SAND2000-1929C

.

5.

Deaner

,

B. J.

,

Allen

,

M. S.

,

Starr

,

M. J.

,

Segalman

,

D. J.

, and

Sumali

,

H.

,

2015

, “

Application of Viscous and Iwan Modal Damping Models to Experimental Measurements From Bolted Structures

,”

ASME J. Vib. Acoust.

,

137

(

2

), p.

021012

. 10.1115/1.4029074

Google Scholar

Crossref

6.

Roettgen

,

D. R.

, and

Allen

,

M. S.

,

2017

, “

Nonlinear Characterization of a Bolted, Industrial Structure Using a Modal Framework

,”

Mech. Syst. Signal Process.

,

84

(

B

), pp.

152

–

170

. 10.1016/j.ymssp.2015.11.010

7.

Iwan

,

W. D.

,

1966

, “

A Distributed-Element Model for Hysteresis and Its Steady-State Dynamic Response

,”

ASME J. Appl. Mech.

,

33

(

4

), pp.

893

–

900

. 10.1115/1.3625199

Google Scholar

Crossref

8.

Segalman

,

D. J.

,

2005

, “

A Four-Parameter Iwan Model for Lap-Type Joints

,”

ASME J. Appl. Mech.

,

72

(

5

), pp.

752

–

760

. 10.1115/1.1989354

Google Scholar

Crossref

9.

Segalman

,

D. J.

,

2010

, “

A Modal Approach to Modeling Spatially Distributed Vibration Energy Dissipation

,”

Sandia National Lab

.,

Albuquerque, NM

,

Technical Report No. SAND2010-4763, 993326

.

10.

Lacayo

,

R. M.

,

Deaner

,

B. J.

, and

Allen

,

M. S.

,

2017

, “

A Numerical Study on the Limitations of Modal Iwan Models for Impulsive Excitations

,”

J. Sound Vib.

,

390

(

C

), pp.

118

–

140

. 10.1016/j.jsv.2016.11.038

11.

Brake

,

M. R. W.

,

2017

, “

A Reduced Iwan Model that Includes Pinning for Bolted Joint Mechanics

,”

Nonlinear Dyn.

,

87

(

2

), pp.

1335

–

1349

. 10.1007/s11071-016-3117-2

Google Scholar

Crossref

12.

Nayfeh

,

A. H.

,

2011

,

Introduction to Perturbation Techniques

,

John Wiley & Sons

,

New York

.

13.

Sumali

,

H.

, and

Kellogg

,

R. A.

,

2011

, “

Calculating Damping From Ring-Down Using Hilbert Transform and Curve Fitting

,”

Sandia National Lab

.,

Albuquerque, NM

,

Technical Report No. SAND2011-1960C

.

14.

Krack

,

M.

,

Panning-von Scheidt

,

L.

, and

Wallaschek

,

J.

,

2014

, “

On the Computation of the Slow Dynamics of Nonlinear Modes of Mechanical Systems

,”

Mech. Syst. Signal Process.

,

42

(

1–2

), pp.

71

–

87

. 10.1016/j.ymssp.2013.08.031

15.

Cook

,

R. D.

,

Malkus

,

D. S.

,

Plesha

,

M. E.

, and

Witt

,

R. J.

,

2007

,

Concepts and Applications of Finite Element Analysis

,

John Wiley & Sons

,

New York

.

16.

Ungar

,

E. E.

,

1964

, “

Energy Dissipation at Structural Joints; Mechanisms and Magnitudes

,”

Defense Technical Information Center

,

Fort Belvoir, VA

,

Technical Report

.

17.

Ames

,

N. M.

,

Lauffer

,

J. P.

,

Jew

,

M. D.

,

Segalman

,

D. J.

,

Gregory

,

D. L.

,

Starr

,

M. J.

, and

Resor

,

B. R.

,

2009

, “

Handbook on Dynamics of Jointed Structures

,”

Sandia National Lab

.,

Albuquerque, NM

,

Technical Report No. SAND2009-4164, 1028891

.

18.

Haslam

,

A. H.

,

Chauda

,

G.

,

Kenia

,

N.

,

Rufat-Meix

,

E. S.

,

Allen

,

M. S.

,

Lacayo

,

R. M.

,

Krack

,

M.

, and

Brake

,

M. R. W.

,

2019

,

Nonlinear System Identification for Joints Including Modal Interactions

,

G.

Kerschen

, ed.,

Springer International Publishing

,

Cham

, pp.

79

–

99

.

19.

Hurty

,

W. C.

,

1965

, “

Dynamic Analysis of Structural Systems Using Component Modes

,”

AIAA J.

,

3

(

4

), pp.

678

–

685

. 10.2514/3.2947

Google Scholar

Crossref

20.

Craig

,

R. R.

, and

Bampton

,

M. C. C.

,

1968

, “

Coupling of Substructures for Dynamic Analyses

,”

AIAA J.

,

6

(

7

), pp.

1313

–

1319

. 10.2514/3.4741

Google Scholar

Crossref

21.

Ginsberg

,

J.

,

2001

,

Mechanical and Structural Vibrations: Theory and Applications

,

John Wiley and Sons, Inc.

,

New York

.

Fast Simulation of a Single Degree-of-Freedom System Consisting of An Iwan Element Using the Method of Averaging

Abstract

1 Introduction

2 Theoretical Background

2.1 Overview of the Four-Parameter Iwan Model.

2.2 Numerical Integration Using the Newmark- $β$ Method.

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

2.3.1 Integrating the Inhomogeneous Ordinary Differential Equation.

2.3.2 Integrating the Homogeneous Ordinary Differential Equation.

3 Test Case 1: SDOF System

3.1 Input: Impulse.

3.1.1 Accuracy: Transient Behavior.

3.1.2 Accuracy: Initial Response.

3.1.3 Speed.

3.2 Input: Sine Beat.

3.3 Input: Bandlimited Random Force.

4 Test Case 2: Sumali Beam Finite Element Model

5 Conclusions

Acknowledgment

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Fast Simulation of a Single Degree-of-Freedom System Consisting of An Iwan Element Using the Method of Averaging

Abstract

1 Introduction

2 Theoretical Background

2.1 Overview of the Four-Parameter Iwan Model.

2.2 Numerical Integration Using the Newmark-β Method.

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

2.3.1 Integrating the Inhomogeneous Ordinary Differential Equation.

2.3.2 Integrating the Homogeneous Ordinary Differential Equation.

3 Test Case 1: SDOF System

3.1 Input: Impulse.

3.1.1 Accuracy: Transient Behavior.

3.1.2 Accuracy: Initial Response.

3.1.3 Speed.

3.2 Input: Sine Beat.

3.3 Input: Bandlimited Random Force.

4 Test Case 2: Sumali Beam Finite Element Model

5 Conclusions

Acknowledgment

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only

2.2 Numerical Integration Using the Newmark- $β$ Method.