## Abstract

In this paper, the forced response of a two degrees-of-freedom (DOF) bilinear oscillator with initial gaps involving inelastic collision is discussed. In particular, a focus is placed upon the experimental verification of the generalized bilinear amplitude approximation (BAA) method, which can be used for the accurate estimation of forced responses for bilinear systems with initial gaps. Both experimental and numerical investigations on the system have been carried out. An experimental setup that is capable of representing the dynamics of a 2DOF oscillator has been developed, and forced response tests have been conducted under swept-sine base excitation for different initial gap sizes. The steady-state response of the system under base excitation was computed by both traditional time integration and BAA. It is shown that the results of experiments and numerical predictions are in good agreement especially at resonance. However, slight differences in the responses obtained from both numerical methods are observed. It was found that the time duration where the DOFs are in contact with each other predicted by BAA is longer than that predicted by time integration. Spectral analyses have also been conducted on both experimental and numerical results. It was observed that in a frequency range where intermittent contact between the masses occurs, super-harmonic components of the excitation frequency are present in the spectra. Moreover, as the initial gap size increases, the frequency band where the super-harmonic components are observed decreases.

## 1 Introduction

In many engineering disciplines, it is becoming important to conduct vibration analysis of damaged structures because the vibration responses can be utilized for structural health monitoring and damage detection. However, in many cases, damaged structures show nonlinearities in their dynamic responses, which hinders the application of traditional vibration tools such as modal analysis and frequency response analysis. In particular, structures with breathing cracks [1–4] and delaminated composites [5] are known to behave as piecewise-linear (PWL) nonlinear dynamical systems because they involve repetitive opening and closing of mechanical boundaries. Namely, the dynamics of PWL systems include sudden changes in the system characteristics, especially the stiffness of the oscillators. Therefore, even if the PWL system is simple, it shows strong nonlinearities in its dynamic responses.

Since the pioneering work by Shaw and Holmes [6], many studies have been conducted to understand the fundamental nature of PWL systems and to explore their practical importance [7–11]. For instance, Ing et al. studied the dynamics of PWL systems that are close to glazing incidence [12]. Dyskin et al. studied resonances of impact oscillators [13]. They also studied the dynamics of a chain of bilinear oscillators [10] and showed that there are some patterns of resonance in the chain of oscillators. Andreaus and De Angelis conducted numerical studies on the dynamic response of a single degree-of-freedom (DOF) oscillator excited by a base acceleration and constrained by unilateral constraints [14]. Later they conducted experimental studies on the single DOF oscillator subjected to unilateral constraints under harmonic excitation [15,16]. Furthermore, they investigated the influence of geometrical and mechanical characteristics of isolation and mitigation devices on the vibro-impact system experimentally and numerically [17,18]. Moreover, the computation of nonlinear vibration problems such as contact phenomena using time integration requires a small time-step size to produce accurate results. Recently, efforts for reducing such computational costs have been made by several researchers [1,7,19–22]. In particular, an approximation method to produce forced response of *bilinear* systems has been proposed by Jung et al. [23], which is called the bilinear amplitude approximation (BAA) method. This method was then extended to be able to handle bilinear systems with gaps or prestress between the oscillators [24]. The method is able to produce an accurate approximation of the forced response of bilinear systems, which arise from mechanical systems involving repetitive opening and closing of contacting surfaces, or *inelastic collision* between the DOF. BAA has also been applied to study the dynamics of high dimensional bladed disks with intermittent contact due to cracks [25,26]. However, the method has never been validated with experimental data. Therefore, in this paper, special attention is given to the experimental verification of the BAA method through comparison between the results obtained by BAA, time integration, and most importantly, experimental data. By providing the experimental verification of BAA, it is expected that it be used for predicting dynamic responses of various bilinear systems with high fidelity. Obtaining high-fidelity experimental data for PWL systems is challenging, mostly because the changes in the dynamics of the PWL systems are non-smooth. Among the PWL systems, several experimental studies have been conducted for impact oscillators. For instance, a forced response study has been conducted on a 1DOF impact oscillator consisting of a mass rolling on a parabolic surface and a rigid stop by Todd and Virgin [27], where chaotic behavior and multiple impacts were observed. Virgin et al. also studied an impact oscillator and showed that rich dynamical behavior can be seen experimentally [28]. Bureau et al. conducted an experimental bifurcation analysis of an impact oscillator [29]. Skurativskyi et al. also investigated a 1DOF impact oscillator with a compliant obstacle numerically and experimentally and reported various rich dynamical behavior including grazing bifurcation and chaos [30]. Numerical and experimental studies have also been conducted for impacting pendulums [31,32]. They investigated the bifurcation behaviors of impacting pendulums from both numerical and experimental studies. To date, however, no experimental work has been conducted on PWL systems involving inelastic collision, which arise in systems involving soft contact, or repetitive opening and closing of mechanical interfaces such as cracks, mechanical interfaces, or delamination. This paper is organized as follows. In Sec. 2, the mathematical model for nonlinear time integration is introduced, and the experimental setup that represents the behavior of the bilinear system is presented. In Sec. 3, results of parameter estimation for the experimental setup are presented. In Sec. 4, the results of numerical and experimental forced responses are presented. In particular, the accuracy and applicability of BAA are discussed in detail. Conclusions of the paper are given in Sec. 5.

## 2 Methodology

This section provides a mathematical model of a 2DOF PWL oscillator, as well as an overview of BAA. The experimental setup for the forced response study is then presented.

### 2.1 Mathematical Model.

*open*and

*closed*states. The open state means that the DOFs are not in contact with each other, whereas the closed state means that they are in contact with each other. Figure 1(a) shows the initial state of the 2DOF PWL model of interest where each of the two masses is connected to a base by a spring and a damper. The masses are separated by an initial gap

*g*

_{0}. The contact between the masses and the base is frictionless.

*m*

_{1}and

*m*

_{2}denote the masses of the DOFs,

*c*

_{1}and

*c*

_{2}are the damping coefficients, and

*k*

_{1}and

*k*

_{2}are the spring constants. Now, consider that the system is excited by a base excitation

*x*

_{0}(

*t*). Then the masses start to move, where

*x*

_{1}and

*x*

_{2}are the displacements of the masses measured from their static equilibria, as is shown in Fig. 1(b). The state of the PWL system switches between the open state and closed state depending on the relative displacement between

*x*

_{1}and

*x*

_{2}. Now, consider introducing a

*gap function*,

*g*(

*x*

_{1},

*x*

_{2}), that represents the gap between the masses with the existence of an initial gap

*g*

_{0}, i.e.,

*g*(

*x*

_{1},

*x*

_{2}) > 0, the 2DOF system is in its open state with no contact between the oscillators as is shown in Fig. 1(b). On the other hand, when

*g*(

*x*

_{1},

*x*

_{2}) ⩽ 0, the system is in the closed state where the oscillators are connected to each other with slight allowance of penetration between the masses, as is shown in Fig. 1(c). Since the system is assumed to follow Newton’s third law, there is a contact force between the masses when they are in the closed state. In this study, the contact force is computed by the penalty method. Namely, the contact force is defined as a product between contact stiffness

*k*′ and the gap function. The equation of motion of the 2DOF PWL system is written as

*f*

_{c}(

*x*

_{1},

*x*

_{2}) is the contact force that hinders the masses from penetration, which is written as

*k*′ switches depending on the relative displacement between

*x*

_{1}and

*x*

_{2}. Namely,

*k*′ is equal to

*k** when the oscillator has penetration or

*k*′ = 0, i.e.,

### 2.2 Bilinear Amplitude Approximation.

*n*DOF mechanical system with piecewise linearity, the equations of motion of the sliding and open states are written as

**x**

_{s}and

**x**

_{o}denote the vectors of generalized coordinates for the sliding and the open systems;

**M**

_{s},

**C**

_{s},

**K**

_{s}are the mass, damping, and stiffness matrices for the sliding state;

**M**

_{o},

**C**

_{o},

**K**

_{o}are the mass, damping, and stiffness matrices for the open state;

**F**(

*t*) is the external forcing vector; and

**F**

_{g}is a constant force caused by the contact stiffness when gaps or prestress exist. Denoting the displacement of the

*i*th DOF measured from the static equilibrium position for the open system as

*x*

_{i}and the displacement of the

*i*th DOF measured from the static equilibrium position for the sliding system as $x~i$, a coordinate transformation can be defined, i.e., $xi=x~i+\delta i$ for

*i*= 1, …,

*n*, where

*δ*

_{i}is the static equilibrium displacement of the masses caused by the contact stiffness

*k**. After applying this coordinate transformation to the sliding system only, (5) can be written as

**x**

_{o}=

**Φ**

_{o}

**q**

_{o}on (6) and (7), they can be written as

**Φ**

_{s}and

**Φ**

_{o}are the modal matrices for the sliding and open states, respectively, and $q~s$ and

**q**

_{o}are the corresponding modal coordinates. Equations (8) and (9) can then be decoupled because of the orthonormality of the modes and are expressed as a set of independent equations as

*n*

_{s}and

*n*

_{o}are the number of modes for the sliding and open states, respectively;

*ζ*

_{s,i}and

*ζ*

_{o,i}are the modal damping ratios;

*ω*

_{s,i}and

*ω*

_{o,i}are the undamped natural frequencies; and

*f*

_{s,i}and

*f*

_{o,i}are the modal forces.

*s*

_{1,i},

*s*

_{2,i},

*o*

_{1,i}, and

*o*

_{2,i}are scalar coefficients of the linear transient responses to be determined;

*ω*

_{sd,i}and

*ω*

_{od,i}are the damped natural frequencies corresponding to the natural frequencies

*ω*

_{s,i}and

*ω*

_{o,i}; and $\theta s,i=tan\u22121(2\zeta s,i\omega s,i\omega /(\omega s,i2\u2212\omega 2))$ and $\theta o,i=tan\u22121(2\zeta o,i\omega o,i\omega /(\omega o,i2\u2212\omega 2))$. The angle

*α*is the phase difference between the excitation and the steady-state responses. The coefficients

*α*,

*o*

_{1,i},

*o*

_{2,i},

*s*

_{1,i}, and

*s*

_{2,i}in Eqs. (12) and (13) are the unknown coefficients that are determined by enforcing the compatibility conditions on the physical displacement and velocity at the moment when the system switches its state [24]. The compatibility conditions are enforced by minimizing the differences in displacements and velocities by an optimization solver with respect to those coefficients and the time in the sliding state. This produces an approximate solution of the underlying bilinear system.

### 2.3 Experimental Setup.

This section describes the experimental setup that is capable of simulating and measuring the dynamics of PWL nonlinear oscillators. These experiments involve inelastic collisions of oscillators with an initial gap between the masses. Figure 2 shows the experimental setup. The system consists of two moving components with adjustable masses, which are placed on a linear guide system, and their movements are aligned along a single direction. The entire experimental setup is fixed on an electrodynamic shaker table (SSV-750, San-Esu, Japan) so that the base excitation to the system can be applied. The setup is placed such that the direction of the shaker’s excitation coincides with the direction of the linear guide. Springs are attached to the two sliding components and they are connected to a fixture that is fixed to the shaker table. By adjusting the position where the springs are attached to the fixture, the length of the initial gap between the sliding components can be adjusted. In order to simulate an inelastic collision between the oscillators, viscoelastic rubber sheets (GP-35L, Naigai Rubber Industry, Japan) are attached to the surfaces of the moving components that are in contact with each other during oscillation. By adjusting the mass of the components, the natural frequencies can be adjusted. This is necessary because the masses of the components need to be adjusted so that their natural frequencies differ in order to ensure that the components collide with each other during the oscillation. The displacements of the moving components and the table are measured by laser displacement sensors (IL-300, Keyence, Japan).

## 3 Parameter Estimation

First, experiments have been conducted to estimate the parameters in (2) by studying the forced response of each single DOF system. This process was necessary to account for the masses of the springs themselves and the viscous frictional forces of the linear guide rail and the ones due to the dynamic deformation of the springs.

*x*

_{0}=

*X*

_{0}e

^{jωt}, the resulting response can also be written as a harmonic function,

*x*

_{i}=

*X*

_{i}e

^{jωt}. By substituting these into (2) with

*f*

_{c}(

*x*

_{1},

*x*

_{2}) = 0, the displacement amplitude of the single DOF system for the base excitation can be analytically obtained as (see, e.g., Ref. [33])

*ω*is the excitation frequency of the shaker,

*ω*

_{i}is the natural frequency of the

*i*th oscillator, and

*ζ*

_{i}is the damping ratio. Based on (14), forced responses were computed with the initial sets of parameters shown in Table 1. These responses are plotted in Fig. 3.

Symbols | m_{1}(kg) | c_{1}(N s/m) | k_{1}(N/m) | m_{2}(kg) | c_{2}(N s/m) | k_{2}(N/m) |
---|---|---|---|---|---|---|

Initial value | 1.84 | 25.00 | 386.72 | 2.39 | 25.00 | 386.72 |

Estimated value | 1.80 | 31.14 | 430.00 | 2.46 | 19.13 | 417.78 |

Symbols | m_{1}(kg) | c_{1}(N s/m) | k_{1}(N/m) | m_{2}(kg) | c_{2}(N s/m) | k_{2}(N/m) |
---|---|---|---|---|---|---|

Initial value | 1.84 | 25.00 | 386.72 | 2.39 | 25.00 | 386.72 |

Estimated value | 1.80 | 31.14 | 430.00 | 2.46 | 19.13 | 417.78 |

*p*_{i}= [

*m*

_{i},

*c*

_{i},

*k*

_{i},

*s*,

*r*],

*m*

_{i}is the equivalent mass,

*c*

_{i}is the damping coefficient,

*k*

_{i}is the spring constant, and

*s*and

*r*are variable weighting coefficients. Δ

*X*

_{i}(

*p*_{i}) and Δ

*ω*

_{i}(

*p*_{i}) are defined as

*s*+

*r*= 1. The minimization was conducted using the function “fmincon” in matlab. The estimated values of the parameters based on the minimization are shown in Table 1, and the calculated forced responses with the estimated parameters are shown in Fig. 3. The weighting coefficients converged to

*s*= 0.76 and

*r*= 0.24 for mass 1 and

*s*= 0.49 and

*r*= 0.51 for mass 2. The interpretation of the estimated values is as follows. First, the optimization process is attempted to increase the values of

*k*

_{1}and

*k*

_{2}or decrease the values of

*m*

_{1}and

*m*

_{2}because the resonant frequencies of the oscillators with the initial parameter values are larger than those obtained from the experiments. Therefore, it is reasonable to see that the values of

*k*

_{1}and

*k*

_{2}are both increased. This also corresponds to the additional stiffness due to the frictional force of the linear guide rail. The increase in

*m*

_{2}makes sense because of the additional inertial effect due to the spring motions. The decrease in

*m*

_{1}was unexpected, however, this is due to the fact that the measured damping of the oscillator 1 was higher than that of the oscillator 2, and the optimization solver is needed to increase

*c*

_{1}and decrease

*m*

_{1}because it needed to increase the value of

*k*

_{1}.

m_{i}(kg) | c_{i}(N s/m) | k_{i}(N/m) | s | r | |
---|---|---|---|---|---|

Upper bound | 3.00 | 50.00 | 430.00 | 1 | 1 |

Lower bound | 1.80 | 0.00 | 380.00 | 0 | 0 |

m_{i}(kg) | c_{i}(N s/m) | k_{i}(N/m) | s | r | |
---|---|---|---|---|---|

Upper bound | 3.00 | 50.00 | 430.00 | 1 | 1 |

Lower bound | 1.80 | 0.00 | 380.00 | 0 | 0 |

Slight differences between the measured and computed responses with the estimated parameters have been observed. These differences are attributed to the fact that the viscous friction between the linear guide also creates a slight nonlinearity from a frictional effect. However, the predicted amplitudes of displacement at the resonances with the optimized parameter values are in good agreement with the measured amplitudes. This agreement is important because the amplitude of oscillations of the system is used for discussing the validity of BAA.

## 4 Forced Response Analysis

With the parameters obtained in the previous section, nonlinear forced response analyses of the 2DOF system have been conducted using time integration and BAA. The results are then compared with the experimental results.

### 4.1 Convergence Study.

First, convergence studies were conducted to investigate the effects of penalty coefficient *k** on the response. The 2DOF system represented by (2) is constructed in the matlab Simulink environment and excited by a swept-sine harmonic excitation with an amplitude of 10 mm and a frequency ranging from 1 to 6 Hz. A fourth-order Runge–Kutta method (ode4) was used as the time integration solver with a fixed time-step size of 0.001 s. In Fig. 4, forced responses of the system are shown for various values of *k**. One can see that the value of *k** greatly affects the results of the forced response. The results do converge when *k** reaches 1.0 × 10^{5} N/m. Therefore, in order to reduce the computational cost while maintaining the numerical accuracy, *k** = 1.0 × 10^{5} N/m was used to compute the nonlinear forced responses that are shown in Sec. 4.2.

### 4.2 Nonlinear Forced Response.

Both experimental forced response tests and numerical forced response calculations have been done for the 2DOF PWL system. In particular, the effects of gap size between the masses on the response and the validity of BAA are discussed. Namely, three gap values are chosen, i.e., *g*_{0} = 0 mm, 2 mm, and 11 mm.

First, forced response tests have been conducted with swept-sine base excitation for the 2DOF PWL system. The excitation has been applied for 1–6 Hz.

Next, to obtain the steady-state response of the system, forced response calculations have been conducted by numerical time integration and BAA. At each frequency, the measured displacement at the base is used for evaluating *x*_{0} in (2). Also, the estimated parameters discussed in Sec. 3 have been used for the numerical simulations. The contact parameter value of *k** for time integration is equal to 1.0 × 10^{5} N/m, which was obtained from the convergence study shown in Sec. 4.1.

First, forced response results for *g*_{0} = 0 mm are discussed. Figure 5 shows the displacement amplitude of both masses versus the excitation frequency for *g*_{0} = 0 mm. As can be seen, results produced by time integration and BAA are in excellent agreement for the entire frequency range. Moreover, they qualitatively agree with the experimental result especially at resonance, with at most 3% error in the maximum amplitude and 5% in the resonant frequency. However, there are some differences in the results of time integration or BAA and that of the experiment. This may suggest that the mathematical model with the estimated parameter values is not accurate enough to quantitatively capture not only the resonance but also the dynamics with frequencies that are higher or lower than the resonant frequency. To take a closer look at the behavior of the masses at the resonance, the time histories of the gap function for a period of oscillation computed by time integration and BAA are shown in Fig. 6. Note that the horizontal axes of Fig. 6 are normalized with respect to the period of excitation. When the gap function becomes equal to or slightly less than zero, it means that the system is in the closed state. As can be seen, both time integration and BAA produced a gap function between 0 and 10 mm. However, some differences between time integration and BAA are observed. Namely, time integration produced a shorter closed state than that produced by BAA.

Second, forced response results for *g*_{0} = 2 mm are discussed. Figure 7 shows the displacement amplitude of both masses versus the excitation frequency for *g*_{0} = 2 mm. Again, the results computed by time integration and BAA are in excellent agreement for the entire frequency range. Note that the response computed by BAA does not contain results with the excitation frequency below 1.25 Hz or above 5.5 Hz because the contact between the masses do not occur in these frequency ranges. Time histories of the gap function at the resonance are shown in Fig. 8. As can be seen, the duration of the closed state predicted by time integration is still shorter than that predicted by BAA. However, the duration of the closed state predicted by BAA has become slightly shorter than that predicted by BAA for the case with *g*_{0} = 0 mm.

Finally, the results of the forced responses with a larger gap size, i.e., *g*_{0} = 11 mm, are discussed. In Fig. 9, there is a primary resonance at 2.1 Hz and a secondary resonance at 2.7 Hz in the results from experiments. This phenomenon is also observed in time integration with slight differences in the maximum amplitude and frequencies, as seen in Fig. 9. On the other hand, the response produced by BAA is smooth at the resonance and does not show the multiple resonances. Figure 10 shows the time histories of the gap function at 2.4 Hz. The contact duration predicted by time integration is still shorter than that predicted by BAA. However, the discrepancy between the results of time integration and BAA, especially in terms of the duration of the closed state, is smaller than that observed for the cases with *g*_{0} = 0 and *g*_{0} = 2 mm. This is because the contact duration becomes small as the gap size becomes large, and the contact states predicted by both methods become close to each other. However, the phenomenon of multiple resonances is not observed in the results computed by BAA. This may be due to the fact that BAA does not capture some frequency components in the response, as will be discussed in Sec. 4.3.

It is noted that as the initial gap size increases, the difference between the time history obtained by time integration and that obtained by BAA decreases, as can be seen in Figs. 6, 8, and 10. This is because the duration of contact period decreases as the initial gap size increases. As can be seen in Fig. 8(a), the time integration method predicts that the gap function increases up to 4mm after the contact and it stays below 5mm until the normalized time is 0.45. Then it increases again and reaches its maximum. On the other hand, based on the results of BAA as shown in Fig. 8(b), the gap function is forced to be in the closed state when the normalized time is between 0 and 0.4. This strict enforcement of the closed state is achieved by an artificial penalty damping factor that is applied when the masses are in contact with each other [24]. This damping effect is not included in the model used for time integration. Therefore, since this enforcement of the closed state generates the difference between the results obtained by time integration and BAA, if this time period decreases due to an increase in the initial gap, BAA will approximate the dynamics of the system better, as seen in Fig. 10.

### 4.3 Spectral Analysis.

In order to further examine the dynamics of the system, particularly when the masses are in the frequency ranges where intermittent contact takes place, spectral analyses have been applied to both experimental and numerical forced responses. For the experimental data, the wavelet transform has been applied to observe the instantaneous spectrum of the experimental data in frequency because the experiment involved a relatively fast frequency sweep. The bump wavelet was applied as the mother wavelet in the wavelet analysis. For numerical responses, the fast Fourier transform (FFT) was applied to compute the frequency spectrum at each excitation frequency. Note that both the wavelet transform and FFT produce comparable results. Figure 11 shows the spectral distribution of the measured and numerically computed responses of both masses with *g*_{0} = 2 mm. For numerically computed responses, both time integration and BAA have been considered. In Figs. 11(a) and 11(d), one can observe an order line that represents the frequency component of the excitation frequency for the entire frequency range. Furthermore, one can also see the super-harmonic components, or the frequency components that are integer multiples of the excitation frequency in the frequency range from 1.5 Hz to 4 Hz for both masses. This means that there is intermittent contact between the masses in that frequency range, because intermittent contact of masses is known to produce super-harmonics. Also in Figs. 11(b) and 11(e), one can see that time integration produces the order line representing the frequency component that is identical to the excitation frequency. Furthermore, one can also observe the super-harmonic components that are observed in the experiments. Note that a subharmonic occurs at about 5 Hz and there is also a zero frequency component present in time integration, but these do not occur in the experiments.

The results computed by BAA are shown in Figs. 11(c) and 11(f). Again, the spectral distributions computed by BAA do not contain results when the excitation frequency is below 1.5 Hz or above 5.5 Hz, because the contact between the masses do not occur in these frequency ranges. Namely, the response is linear in these frequency ranges, and these linear responses contain the frequency component that is identical to the excitation frequency, which is denoted by the dashed line. As can be seen in Figs. 11(c) and 11(f), in addition to the line corresponding to the excitation frequency, a frequency component that is twice the excitation frequency is also observed. Other than these frequency components, no other super-harmonic components are observed, even when the intermittent contact between the masses is strong. This means that BAA does not produce higher harmonics in the response, even though it captures the behavior of intermittent contact between the masses. This aligns with the key idea behind BAA, which assumes that there are only two states in the dynamics of the bilinear system, i.e., closed and open states. It is interesting to note that BAA does not have any subharmonic or zero frequency components, which, unlike time integration, matches the experiments. Furthermore, it is noted that there is a possibility that higher order nonlinearity has been induced at the contact. For instance, for *g*_{0} = 2 mm at 2.4 Hz, the gap function produced by BAA is closed from 0 to approximately 0.4 normalized time and open for the rest of the vibration cycle, as can be seen in Fig. 8(b). This means that the contact was assumed to be perfectly inelastic or soft. On the other hand, from 0 to approximately 0.4 normalized time, the gap function produced by time integration is open with a non-zero gap opening of 3–4 mm. This means that the contact between the masses was not perfectly inelastic. This behavior of the masses generated higher order harmonics in the spectra of the response of both masses. Since we observe the higher order harmonics, both in the results of experiments and of time integration as shown in Fig. 11 such phenomenon likely occurred during the experiments.

Figure 12 shows the spectral distributions of the displacement of both masses with a larger initial gap size, *g*_{0} = 11 mm. As seen in Figs. 12(a) and 12(d), measured displacements contain not only the frequency component that is identical to the excitation frequency but also the super-harmonic components especially in the vicinity of the resonances. By comparing Fig. 11(a) or 11(d) with Fig. 12(a) or 12(d), one can see that the width of the frequency bands where the super-harmonic components are observed are narrower in Fig. 12(a) or 12(d), than those observed in Fig. 11(a) or 11(d). This trend is shown even more dramatically in the results for time integration, as shown in Figs. 12(b) and 12(e). Namely, the frequency band where super-harmonic components are observed is much smaller than that observed in Figs. 11(b) and 11(e). Since the occurrence of super-harmonic components in the response indicates that intermittent contact between the masses occurs, this result indicates that when the gap size between the masses increases, the frequency band where intermittent contact occurs decreases. This means that as the gap size increases, the masses tend to make less contact with each other as they vibrate. Indeed, the time history of the gap function at resonance shown in Fig. 10(a) shows that the duration of contact is short, which justifies this argument. Figures 12(c) and 12(f) show the results computed by BAA. As can be seen, the intermittent contact between the masses occurs only in the vicinity of the resonance. This is in agreement with the results given by time integration. This time for BAA, almost no super-harmonic component is observed in the response, which means that even though intermittent contact between the masses occurs in the frequency range, the contact is so mild that it almost has no effect on the dynamics of the masses. This indicates that for the prediction of the dynamics of bilinear systems involving higher harmonic components over a very short time period, such as during impact, further extension of the method may be necessary to fully capture the dynamics. However, it is important to note that BAA does match time integration quite well (see Fig. 9), and there is a bigger difference between the experiment and both computational methods.

## 5 Conclusions

In this paper, the forced response of bilinear systems with initial gaps involving inelastic collision was investigated by numerical simulation and experiments. A focus was placed upon the experimental verification of the generalized BAA, which was developed for the accurate estimation of forced responses of bilinear systems. By conducting spectral analyses on the experimental results and numerical forced responses of the bilinear system, it was shown that super-harmonic components of the excitation frequency can be observed. Moreover, for a small gap size, time integration appeared to match the higher harmonic content of the experiments better than BAA. This is due to the assumptions in the BAA method that lead to some simplifications in the response that remove some of the higher order harmonics. At the same time, BAA was a better match of the experiments than time integration when considering subharmonics. Furthermore, it was shown that as the gap size becomes large, the frequency range that includes super-harmonic components becomes narrow. To further expand the capability of BAA, extension of the method to account for higher-order dynamics would be beneficial.

## Acknowledgment

This paper is based on work partially supported by the National Science Foundation under Grant No. 1902408, program manager Dr. Robert Landers. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the National Science Foundation.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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