Abstract
Damage or failure of a relatively small component can precipitate the failure of a larger part of a structure. The behavior of damaged or worn joints is of particular concern. To address contact/impact in structural systems, this work models a structural lap joint from first principles. A beam with four stops and gaps is used to simulate a loose or damaged lap joint, which also represents designed manufacturing clearances in mechanical systems. The goal is to generate frequency responses to identify the local shock effect due to impact. Spatial and temporal solutions are presented for an example case. Converged time histories were used to generate the impulse as a metric of frequency response. Facilitating mode contribution calculations, the metric of impulse proves to be an excellent indicator of complexities in the beam's motion due to excitation frequency. Noncontact regions, sticking motions, local extrema, grazing impacts, and aperiodicities are identifiable for specific operating parameters. These conditions indicate when harmful impact may occur that can ultimately cause local damage within a structure. Knowledge of dangerous operating conditions can better focus on inspection before propagation occurs.
1. Introduction
Damage or failure of a relatively small component can precipitate the failure of a larger part of a structure. Even the failure of a single connection can cause overload on other structural members, which may consequently fail. Thus, the behavior of damaged or worn joints is of particular concern. In addition to worn parts, many engineering systems with design/manufacturing clearances lead to nonlinearity and contact. The dynamics of loose joints involve nonlinear support conditions and contact/impact phenomena.
Many studies have been performed to model contact/impact, and Gilardi and Sharf [1] provided a good review of different methodologies. Moon and Shaw [2] used a digital simulation of a single mode mathematical model of a cantilever elastic beam impacting a stop at one end. The Runge-Kutta algorithm along with finite difference method was used to provide the integration of bilinear spring equations, which show chaotic vibrations that are qualitatively similar to experimental observations. Shaw [3] considered an elastic beam with a one-sided amplitude constraint subjected to periodic excitation. Experimental results were compared with those given by a theoretical model based upon a single mode analysis of [2]. With just a single mode, the model predicted multiple subharmonic resonances, period doublings, and some chaotic regimes, all of which were found experimentally.
Pun et al. [4, 5] analyzed the free and forced vibration behavior of an L-shaped beam with a limit stop by using the incremental harmonic balance method. The authors used frequency response functions (FRFs) to examine complex internal resonances. In free vibration, their analysis revealed the presence of multiple internal resonances involving interactions among the first five modes. Their results showed that mode interaction may occur if harmonics of lower modes are near to the natural frequency of a higher mode. The conditions for the existence of internal resonance were explored, and it was shown that a prerequisite is the presence of bifurcation points in the system’s frequency response. In one undamped case study, three in-phase and two out-of-phase solution branches were found. The resonance curve is extremely complicated with multiple branches and interactions between the first four modes. The amplitudes of the higher harmonics are heavily influenced by damping, which can effectively attenuate internal resonances.
Metallidis and Natsiavas [6] investigated the response of a continuous system with clearance and motion-limiting constraints. An exact periodic response analysis was presented for a periodically excited deformable rod whose motion was constrained by a flexible obstacle. Direct integration of a rod’s equation of motion for selected parameter combinations revealed chaotic motions.
Chattopadhyay [7] presented an approximate approach for a continuous system with nonlinear boundary conditions. A two-span clamped-free undamped uniform beam with an intermediate two-sided motion limiting constraint (two-sided gap) was considered. The qualitative behavior of this cantilever beam under harmonic base excitation was shown to exhibit disconnected resonance curves and multiple valued responses typical of a nonlinear system. Once the impact becomes steady and periodic, the nonlinear shearing force also becomes periodic.
Ervin and Wickert [8] investigated the forced response dynamics of a clamped-clamped beam with an attached rigid body impacting a compliant base structure. The impacts were modeled by a linear restoring spring force. The method of separation of variables was used to solve the piecewise linear equation of beam motion with the rigid body at each state of in-contact and not-in-contact. Compatibility conditions were applied at their junctions. Increase of contact stiffness shifts the frequency of single peak pairs to higher frequencies as well as decreases peak amplitude. Also, the frequency response for higher contact stiffness is more complicated, especially for frequencies between those of coordinated contact and noncontact modes. Ervin [9] also analytically investigated the repetitive impact dynamics of two orthogonal pinned-pinned beams subjected to base excitation using the same method presented in [8]. The same method used in [8, 9] will be used herein.
To address contact/impact phenomena in civil systems, this work will model a structural lap joint from first principles. In this work, a beam with four stops and gaps is used to simulate a loose or damaged lap joint, which also represents designed manufacturing clearances in mechanical systems. The goal is to generate frequency responses to identify the local shock effect due to impact.
2. Materials and Methods
This section presents Euler-Bernoulli beam theory as applicable to impact/contact. The equation of motion for the forced lateral vibration of an undamped Euler-Bernoulli beam is provided in such references as Rao [10]. The spatial solution for a damped beam is the same, but the temporal solution is different for contact/impact.
2.1. Equation of Motion
Let a beam be subjected to the sinusoidal table motion with amplitude of and frequency of . Also let be the absolute displacement, the relative displacement, and the damping per unit length of the beam. The equation of motion for the forced lateral vibration of a nonuniform damped beam subjected to the base motion is where is Young’s modulus, is the moment of inertia of the beam’s cross section about the -axis [11], is the beam’s mass density, and is the beam’s cross-sectional area. Note that right hand side in (1) indicates that the harmonic base excitation is equivalent to applying a negative harmonic force to the mass.
The relationship between deflection bending moment and shear force can be expressed as
Two initial conditions and four boundary conditions are needed since (1) involves a second order time derivative and fourth order spatial derivative. The solution of this differential equation may be obtained by many methods, Laplace transforms for example. In this work the method of separation of variables will be used because it provides a specific insight into the beam’s physical motion. Denoted as modal decomposition, let where is the relative displacement, is the spatial solution or mode shape, and is the generalized time solution. Substituting (3) into (1) leads to where and the dot superscript represents a single derivative with respect to time. Rearranging (4) provides where is a natural frequency of system and is a positive constant. Equation (6) can be written as two equations: Rearranging provides where The next section will handle spatial considerations (see (8)) and Section 2.3 will address the temporal solution (see (9)).
2.2. Natural Frequencies and Mode Shapes
The spatial solution of (8) can be expressed as where , , , and are constant coefficients. These unknown constants and the value of can be determined from the boundary conditions of the beam.
In order to calculate the natural frequencies and mode shapes, the general beam in Figure 1 is assumed. Representing adjustable boundaries at and , is the lateral spring stiffness, while is the torsional spring stiffness. Due to contact at , any mode shape of the entire beam is divided into two parts: Since there are eight unknowns in (12), eight conditions are needed. Four boundary conditions apply on shear and moment at and (13a), (13b), (13g), and (13h). The other four are compatibility conditions at the contact points (13c), (13d), (13e), and (13f); that is, the displacement, slope, shear and moment just to the left of the contact point must be equal to the same displacement, slope, shear, and moment just to the right of the contact point . These boundary conditions are as follows:in which is shear force, is bending moment, is deflection, is slope, and represents a derivative with respect to . Assuming a prismatic beam, (2a) and (2b) providesUsing (13a), (13b), (13c), (13d), (13e), (13f), (13g), and (13h) with (14a) and (14b), eight equations will be obtained. These equations can be written in a matrix form as where the 8 × 8 -matrix contains terms employing , , , and as well as contact spring stiffnesses, trig functions, and hyperbolic trig functions. This is an eigenproblem as the coefficient matrix is noninvertible. For a nontrivial solution, the determinant of the coefficient matrix must vanish. The eigenvalues are obtained by setting this determinant equal to zero. Each can then be substituted in (15), and all eight coefficients respective to can be obtained. This solves the eigenvector problem and provides a mode shape.
Employing a single contactor at , the model in Figure 1 is insufficient to represent a structural connection. The more complex beam model used in this work is sketched in Figure 2. A cantilever beam with four springs will be modelled to analyze two contact locations in an overlapping joint.
Based upon any potential beam displacement, nine possible states of motion were assumed. In the free state, there is no spring in contact with beam. States 2 and 4 represent the states in which one of the end springs (either the bottom one or the top one at ) is in contact with beam. In states 3 and 5, the beam is in contact with either bottom or top intermediate spring at . In states 6, 7, 8, and 9, two bottom springs, two top springs, or one of each are in contact. These states are defined as shown in Table 1.
For each set of cases, the stiffness of springs are described in Table 2. Note that the upper springs ( and ) have the same stiffness as the lower contactors ( and ). The support at has a translational stiffness and rotational stiffness that are both set to approach infinity. A two-sided displacement limiting contactor is located at both and , and the relative contact stiffnesses are represented by the pair of and and the pair of and , respectively.
Gaps between the springs and the beam are adjustable to represent tight to loose joints due to damage. Note that mode shapes were obtained for the boundary conditions in which the gap between the beam and supports (springs) were not directly considered as shown in Figure 1. Since a beam with 2-sided gap is considered, the gap will be added to the temporal solution later as an initial displacement of each case. This will be discussed in more detail in Section 2.3.
It is desirable to normalize the mode shapes in respect to mass for modal analysis. This normalization makes the subsequent use of normalized mode shapes more convenient. Using a modal mass scheme with the mass normalized mode shapes are defined as The orthogonality of mode shapes for transverse beam vibration was tested and verified in this work. Thus, This principle was even proven for two orthogonal beams by Ervin [9].
In order to perform numerical simulations, the parameters listed in Table 3 are used. The material is steel, is the moment of inertia for a rectangular cross section, and the contact point was chosen at the middle of the beam .
Table 4 provides the first three natural frequencies for different states and two values of spring stiffnesses. The contact stiffnesses are nondimensionalized by the cantilever beam stiffness, so . For the contact is considered as “soft,” and for the contact is considered as “hard.” The active springs column describes which of the springs in Figure 2 are in contact for each state. In the free state, there is no active spring/contact involved, so the stiffness of springs does not affect the natural frequencies. For other states, increasing the stiffness leads to an increase in the natural frequencies. Furthermore, the effect of this stiffness on the lower frequencies is greater than on the higher frequencies [12].
The effect of contact stiffness is also evident in mode shapes. For N/m ( or “soft”), the stiffness of springs is relatively small. Thus, there is little difference between a mode shape of a contact state and free state. The mode shapes show more variation for N/m ( or “hard”). Here, the mode shapes are restricted in deflection at contact points due to the high stiffness. For example, in states 2 and 4 where the end springs are contacting, the displacement of the endpoint of the mode shapes is significantly less than the mode shapes at the same point for state 1, the free state. For example in “hard” contact, the mode shape at in state 2 is 0.32% and 1.67% of the mode shape at free state for the first and for second modes, respectively. The same behavior is shown in states 3 and 5 for the midpoint as well as in states 6, 7, 8, and 9 for midpoint and endpoint.
2.3. Temporal Solution
The forced vibration solution of a beam can be determined using the mode superposition principle. For this, the relative deflection of the beam is assumed as where is the th mass-normalized mode shape for the beam and is the generalized coordinate in the th mode. There are four gaps at each of the four springs; to better represent the deflected shape these gaps are converted to a gap function . This function is equal to the displacement of the beam at just before the change of condition. In other words, the displacement of the beam at each time is equal to the sum of the displacement of the beam at the end of previous condition (, previous deflected shape) and the displacement of the beam in new condition . Equation (19) can be rewritten at each dynamical condition of the beam as where is the change counter on the dynamical condition of the beam. is the th normal mode shape for the beam after changes in condition, and is the generalized coordinate in the th mode after changes in condition. is the displacement of the beam just before the last change in condition (previous deflected shape) such that For example, consider the beam with only one spring at the bottom end. There will be only two possible dynamical states for such a system. The first one is free, and the second one is when the end spring is at contact. Each time the beam contacts the spring (impact) or leaves the spring (rebound), will increase by one; each time the beam leaves the spring (rebound). This index is initially equal to one. Now consider that the beam is initially free, and it is excited. For each impact state, the is equal to the displacement of the beam just before the impact at free state (previous condition) and vice versa for rebound. This system as well as two possible dynamical states and mapping between them are shown in Table 5. For a beam contacting four springs, the same concept is still true with the difference that there are nine possible dynamical states as mentioned in Section 2.2.
By substituting (20) into the equation of motion of the beam subjected to the base motion, (1), Using (8) Multiplying both sides by and using orthogonality of mode shapes, Using (16) and rearranging Rewriting where is the modal mass and it is equal to 1 due to normalization, is the modal damping of , and is the modal stiffness of . This equation is a second order nonhomogeneous differential equation, and thus two initial conditions are required. The solution contains a homogeneous part and two particular parts due to the base excitation and additional deflected shape . The particular parts are generated by the associated family of solutions, and then general solution can be obtained by using two initial conditions. For the base excitation part of the particular solution, the first part of right hand side of (26) with the base’s sinusoidal equation provides where is second derivative of sinusoidal table motion (acceleration) with amplitude of and frequency of . Since the right hand side is a sine function, the solution should be a combination of sine and cosine, so and are determined as
For deflected shape portion of particular solution, the second part of right hand side of (26) is not a function of time. Thus the second part of particular solution is in which the numerator is the equivalent shear force causing the initial condition from the deflected shape. The displacement is derived by dividing the equivalent force by the stiffness .
The th homogeneous solution of (26) is where and are constants which can be determined using initial conditions. is th damped frequency, and it is calculated as The damping ratio for all states is defined as The general solution of (26) is the sum of (29), (31), and (32) such that where and (see (32)) are unknown.
Two initial conditions are needed to calculate constant coefficients and . The displacement and velocity do not change during a switch of dynamic state; in other words, if the contact occurs at , the displacement and the velocity of the beam are equal just before the impact/rebound at and after the impact/rebound at . For displacement compatibility, Substituting (20) provides Further substituting of (21) provides Note that is the time in which the beam is leaving the state and is entering state . Then by cancelling the equal terms and rearranging, Multiplying both sides by and using orthogonality, where is the time of contact or rebound. Substituting in (20) gives the deflected shape at contact which is equal to deflected shape before contact.
For velocity compatibility, Substituting (20) provides Multiplying both sides by and using orthogonality, Substituting (16), is called the mapping or transformation matrix from state to state . By substituting these two initial conditions in the general solution and its derivative and defining the elements of coefficient matrix [a] as Two unknown constants can be determined as Using the calculated , the general solution will be The first term on the right hand side denotes the transient vibration, and the other two terms represent the steady state vibration. Once (49) is solved for , the total solution at each condition can be determined by the summation of (20).
Note that due to numerical instabilities or limitations in computational power, determining mode shapes up to infinity is impossible. In this study the mode shapes are determined up to modes, where represents the minimum number of modes to which the solution converges. Convergence occurs when the relative error between solution for modes and solution for modes is less than a prescribed tolerance.
3. Results and Discussion
Model results for a specific parameter set are presented for the purposes of understanding the complex response. This study analyzes a cantilever beam with four contact springs (Figure 2) experiencing base excitation. As shown in Figure 3, the four points (P2, P3, P4, and P5) have associated four springs (, , , and ) and associated gaps (, , , and ) with respect to the beam. As does point P4 to , point P5 represents the displacement at the top of the beam that coordinates with upward motion through to touch . As does point P2 to , point P3 represents the displacement at the bottom of the beam that coordinates with motion downward through to touch . Note that all springs move with the base, so the relative gap is constant over time. A root-finding algorithm using adaptive time stepping and the bisection method are used to solve for the time of contact at any of the four possible points. For each time step, this is accomplished by testing the displacements at P2, P3, P4, and P5 versus gap (, , , and ).
The basis parameters in Table 6 were selected to generate a flexible cantilever with four identical contactors. The two intermediate springs are located at the halfway point, and the two boundaries are fixed () and displacement limited (). For finding the natural frequencies, a bisection algorithm with maximum machine precision (epsilon = 2.22e−16) was used on intervals of 7.38e−15 rad/s. Also, for finding the time of contact, the same bisection algorithm was used on intervals of 2.2e−11 s. The relative displacement of points P2, P3, P4, and P5 due to sinusoidal base excitations at different frequencies is calculated from initial conditions through the transient to steady state.
A comprehensive convergence study was performed to determine the minimum number of mode shapes to which time history converges. Increasing beyond this will increase the computational time and required hardware with diminishing returns. In short, beam response does not significantly vary beyond eight mode shapes for the free case and five mode shapes for contact cases. The acceptable tolerances on maximum free and contact displacements are 8.0e−6 m and 6.1e−5 m, respectively.
The response of the system is then calculated for a range of base excitation frequencies (). To quantify the response effects directly, comparing time histories is inefficient due to wide variety of transient and steady state responses. An indicator is needed to show the effect of varying system response for different excitation frequencies. In this work, impulse is used as the indicating function: impulse is the change in the momentum of the system and indicates dynamic load applied on the system. Impulse for each cycle is defined herein as the area of relative deflection between impact and rebound multiplied by the respective spring stiffness; thus, traditional force multiplied by time units is reached. The impulse is not a perfect indicator, but it well represents not only deflection but also time in contact. When either parameter slightly changes, the frequency response alters and can even appear as aperiodic response. Impulse can also show sudden jumps in the frequency response when grazing impact occurs; these are an artificial effect [8].
Thus, for each excitation frequency, the relative displacement is run through steady state, and a few repeating cycles are stored. If a time history appears transient after many thousands of impacts, aperiodic response is assumed, and the next 1,000 cycles were captured. Then the relative system responses for the entire beam and the impulses of bottom and top end points (P2 and P4, resp.) were calculated. The impulse frequency response was then plotted as provided in Figure 4 for the example configuration. Note that excitation frequencies lower than 230 rad/s reveal two conditions: no contact (no plot entry) or sticking motion (a zero-value plot entry). At high frequencies the response is simply repeating the excitation with radian phase lag. That is, the beam is moving with the table; while contact does occur, it is more quasistatic rather than shock in nature.
The impulse frequency response of Figure 4 illustrates that the first system response peak occurs at 515 rad/s. This value is just 3 rad/s lower than the first natural frequency of dynamical states in which two springs are in contact (states 6, 7, 8 and 9), indicating the dominance of these modes. The second system peak exists at 820 rad/s, and it is just 6 rad/s higher than the second natural frequency of the states with both springs in contact. The third peak occurs at 2950 rad/s, and it is 20 rad/s lower than the sixth natural frequency of the states with both springs in contact.
The observation that the combined motion incorporates both noncontact and contact frequencies is quite standard for nonlinear dynamics; however, the quantification of each contribution is unique to this study. These peak responses show that the system response is controlled by contact frequencies, but the slight differences show the influence of the free mode set. That is, the first peak is pulled to a 0.58% lower frequency by the second natural frequency of state 1 (free). The second peak is pulled slightly higher (0.78%) due to third natural frequency of state 1. The third peak is pulled to a 0.67% lower frequency by the fifth mode of state 1. Note also that the ratio of impulse at the second peak to the first peak is 0.5 and between the third peak and second peak is 1.5.
There also exist several regions of interest which show abrupt changes of impulse, such as near 1000 rad/s. In general, these occur at some harmonics of states with two springs in contact. Sudden changes of phase lag most often occur at the frequencies of impulse discontinuity. Additionally, the impulses at beam’s top and bottom often become asymmetric. Near zero velocity impact (grazing impact) also occurs, which is observed as near zero impulses in Figure 4.
Unique aperiodic responses are observed at two locations: approximately 1398 rad/s and 2160 rad/s. In the former region, the impulse response does not reach steady state even after 2,500 impacts and rebounds for the excitation frequencies between 1398 rad/s and 1399 rad/s. For the 1,000 cycles shown, the response is not steady but it is bounded to a limit; for example, at 1399 rad/s, the impulse does not exceed 456.88 N·s. This is not the case for the region between excitation frequencies of 2160 rad/s and 2161 rad/s. Here, the impulse appears unbounded. The response does not reach steady state even after 4,000 impacts and rebounds; and examination of the time histories explains the reason: an asymmetrical internal resonance. For the frequency of 2060 rad/s, the deflection of spring increases after each impact while the deflection of spring decreases. As a result, the impulse for point P2 decreases, and it will reach zero because point P2 does not impact the spring after a certain time. Meanwhile, the impulse for point P4 increases. At the frequency of 2061 rad/s, the deflection of spring increases after each impact; as a result the impulse for point P2 starts from a low magnitude and increases after each impact. Meanwhile, the impulse for point P4 decreases to zero. Thus, the contact at these two locations experiences beat phenomenon.
4. Conclusions
Impulse is an excellent measure of the shock input to the structure due to contact/impact. Measuring both contact time and contactor deflection, the frequency response on impulse reveals complexities in the beam’s motion due to excitation frequency. The metric of impulse represents the shock load on a structure; however, it can produce what appears as unusual behavior that requires explanation via time history investigation. Impulse can also identify noncontact and sticking motions, which usually are not a dynamical matter of concern in this work. Peaks of potential damage are identifiable by peaks in the impulse frequency response; these are regions to avoid for any operating frequencies. Due to low amplitude, a displacement-limited frequency region may be a reasonable operating realm despite the aperiodicity. In contrast, a design must avoid the internal resonance near 2160 rad/s for this example.
The aim is to analyze loose joints in order to eventually determine their reliability before structural failure occurs. Affecting repetitive impacts, instabilities occur where the phase lag is in transition. Aperiodic response occurred at just two regions, and at one region it was limited by an upper bound. The system response is most affected by mode shapes of states with two springs in contact; for the major portion of a full cycle, the beam has both springs on the same side in contact. These conditions at which harmful impact occurs can ultimately cause local damage within a structure. The knowledge of operating conditions and any instabilities or peaks can better focus on inspection before propagation occurs.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of the paper. No financial gain is expected from this work. The authors retain the rights to the enclosed materials, and any reproduction of these contents requires the corresponding author’s consent.