Nonlinear Vibrations and Chaotic Behavior of Curved Two-Layer Beam-Like Structures Used in Automotive Industry

The presence of dry friction in structures and machines often leads to undesirable side effects. The continuous changes between static and dynamic friction in these systems result in the nonlinear stick-slip phenomenon. The presence of the stick-slip phenomenon, in addition to making the structure's vibrations unpredictable, causes noise during motion, generates heat, and ultimately leads to fatigue in multi-layer structures. Therefore, the proper and appropriate design of mechanisms that utilize such structures requires accurate modeling of interlayer friction, taking into account the nonlinear stick-slip phenomenon [1].

Multi-layer structures are commonly used in industry as vibration dampers in structures and machines [2]. With technological advancements, the application of such mechanisms is increasing due to their high stiffness and suitable damping capacity. The mentioned benefits have led to the use of layered beams as vibration dampers, such as leaf springs in the suspension systems of rail and road vehicles, and multi-layer structures with bolted and riveted connections or retaining clips. Alongside these benefits, due to the presence of friction in the structure, interlayer friction during vibrational motion constantly transitions from static to dynamic and vice versa, which is referred to as the stick-slip phenomenon [3].

Research on multi-layer beams began in the mid-1960s. Over the years, numerous researches have been conducted in this field. It can be said that the first research on two-layer beams was conducted by Goodman and Klumpp [4]. They considered a two-layer beam with one end fixed, which was a simplified model of turbine blade roots, assuming uniform pressure between the layers and the presence of dry friction between the layers. They modeled the full-slip friction, derived expressions for energy loss and optimal pressure to minimize the stress applied to the layers and increase the fatigue life of turbine blades, and validated their results using empirical data. However, it should be noted that the amount of energy loss due to dry friction between the layers in turbine blades is not only controllable by the pressure between the layers but also depends on the number of layers of the structure.

Adam et al. [5] investigated the transverse vibrations of a two-layer beam with interlayer slip. In this analysis, each layer of the beam was considered separately as the Euler-Bernoulli beam, and the relationship between interlayer slip and the shear force resulting from the presence of clamping bolts was modeled as a linear spring. Finally, for transverse vibrations, a sixth-order differential equation was derived and analyzed under different boundary conditions using modal analysis. Then, assuming low damping, they examined the governing equations for a two-layer beam with pinned-pinned boundary conditions and plotted the response for different loading conditions, considering various modes. Damisa et al. [6], considering a non-uniform pressure profile between the layers of a clamped two-layer beam with slip friction, calculated the small-amplitude transverse vibration response analytically. They then derived the energy loss in one cycle and calculated the interlayer pressure to optimize the dissipated energy. They also obtained their results for harmonic and non-harmonic end forces and concluded that frequency ratio changes play a more significant role than pressure profile changes in optimizing the energy loss in the system. Shokrieh and Rezaei [7] investigated a four-layer steel leaf spring under a specific static load, obtained the maximum displacement and stress using the finite element method, and showed its agreement with experimental results. They then examined the problem of replacing this steel leaf spring with a composite leaf spring, aiming to minimize the spring's geometry (weight of the composite leaf spring) under constraints such as maximum stress and maximum spring displacement. They concluded that the optimal transverse profile of the spring is a second-degree curve, and the spring thickness is linear, while the weight of the composite spring is 80% less than the steel spring, and the applied stress is also significantly less than that of the steel spring.

Li and Hua [8] conducted studies on the free vibrations of a composite two-layer beam. They considered the shear displacement theory (in this hypothesis, the horizontal displacement of the beam is considered as a sinusoidal function along the beam's cross-section), derived the governing equations for the free vibrations of a two-layer beam, and then, using the obtained equations and a method called Wittrick-Williams, calculated the six natural frequencies and mode shapes of the composite two-layer beam.

As mentioned in the previous section, one of the consequences of interlayer friction in two-layer structures is the alternating transition between static and dynamic friction. The presence of the stick-slip phenomenon complicates the governing equations and distorts the phase space in the dynamics of nonlinear systems. Awrejcewicz and Pyryev [9] in their research, investigated the governing equation for the motion of an elastic shaft inside a solid bushing with dry friction between them. In this system, they observed the stick-slip phenomenon between the shaft and the bushing. They analytically examined this system and found conditions for predicting chaotic phenomena, validating their results numerically and defining conditions for controlling chaotic phenomena. Kenmoé et al. [10] considered a model of a system with dry friction and a nonlinear spring force and investigated the nonlinear stick-slip phenomenon and the effect of the nonlinear spring parameter on this phenomenon. Their research included analytical and numerical results, and in their investigations, they found conditions for the existence of periodic, quasi-periodic, and chaotic behaviors of the system. Finally, they studied the dependence of the dynamic behaviors of the mentioned system on the defined nonlinear spring parameter. van de Vrande et al. [11], by examining an independent dynamic system, predicted the stable and unstable regions of stick-slip vibrations in the presence of dry friction. In this research, they approximated the two-parameter friction force using a single-parameter function for analytical examination of the governing equations. The estimation of the largest Lyapunov exponent using the phenomenon of chaos synchronization in a dynamic system in the presence of dry friction was performed by Fu and Wang [12]. They used an innovative method by drawing bifurcation diagrams and the largest Lyapunov exponent to investigate chaotic regions in the system. Stefanski and Kapitaniak [13] introduced a new method for estimating the largest Lyapunov exponent in nonlinear dynamical systems in the presence of discontinuous terms. To validate the proposed algorithm, they considered a Duffing system with a discontinuous friction term and showed the accuracy of the proposed method by comparing bifurcation diagrams and the Lyapunov exponent. Kang et al. [14], using numerical and analytical methods, investigated the nonlinear dynamic behavior of an oscillator with two frictional contact surfaces. They found two response patterns in the stable regime of the mentioned model. Awrejcewicz et al. [15], by considering a two-degree-of-freedom system in the presence of friction and harmonic driving force, calculated the largest Lyapunov exponent. By drawing bifurcation diagrams and the system's phase space, they also studied the dynamic behavior of the system and identified chaotic and regular regions.

In some dampers, multi-layer springs alongside a lubricant are used as vibration dampers in industry. The elastic properties of multi-layer beams, along with their energy dissipation capacity due to interlayer friction in the presence of a lubricant, can be used as a vibration damper (Figure 1). Also, a type of layered structure known as multi-layer leaf springs is used in road and rail vehicles as suspension systems. These types of springs have high strength to withstand heavy loads and suitable damping capacity to dissipate energy resulting from vehicle collisions with road obstacles. An example of these springs is shown in Figure 2.

The study of nonlinear dynamics and the slip behavior of beam-like structures has advanced tremendously in the past several years. The complexity of the vibration responses is influenced by many different types of conditions, including friction, joint conditions, and layering of materials, which has allowed for significant research to be explored. There has been research on nonlinear oscillators with irrational nonlinearities [16], antisymmetric discontinuities with elastic foundations [17], and interfacial slippage in laminated beams conforming to the rules of Coulomb friction [18], [19]. In all cases, it has been shown that complex forms of nonlinear behavior are present and significant to include in proper analytical and numerical models. In addition, researchers have shown us that slip can affect structural integrity, a degradation of stiffness, and energy storage and dissipation in layered systems. Through the use of experimental and computational methods, including frequency response analysis and artificial neural networks, researchers have conducted work to develop and optimize damping characteristics for multi-layered riveted cantilever beams [20]. Some researchers have been able to find ways to account for interfacial properties in complicated systems, consisting of many layers, in analytical studies as well and demonstrate that interface stiffness would show a nonlinear dependence on the behavior of the multi-layered system [21]. One area of research has also included various methods of stick-slipping and how it would affect the resonance, stiffness, and stability of the overall system in their models, as they applied to dynamic models of multi-beam structures with variable joint conditions [22]. In summary, work done by researchers studying systems, or incorporating systems with frictional interfaces and nonlinear dynamic analysis, will ultimately assist engineers in designing and optimizing advanced engineered structures.

Given the mentioned points and a review of previous research in the field of the dynamics of layered structures, it is clear that in previous studies, the effect of stick-slip phenomenon in the presence of interlayer friction in curved two-layered beams has not been considered. The current research, in line with completing previous studies on the vibrations of two-layer beams, investigates the nonlinear vibrations of curved two-layer beams in the presence of interlayer slip, taking into account the stick-slip effect. After deriving the governing vibration equations for the dynamic behavior of curved two-layer beams, different vibration regimes are identified by changing system parameters. One of the main functions of layered structures is to absorb vibrations caused by external excitations. The results of numerical simulations are compared and validated with finite element results in ANSYS software.

In this section, before deriving the differential equations governing the vibrations of curved two-layer beams, due to differences in compatibility relations, natural frequencies, characteristic equations, boundary conditions, and generally the vibration equations of curved beams compared to straight beams, the motion and compatibility relations of single-layer curved beams are first developed. One of the main challenges in examining the dynamic behavior of curved beams is the calculation of natural frequencies and their mode shapes. For this purpose, by deriving the vibration equations and boundary conditions, the natural frequencies and mode shapes of a single-layer curved beam are numerically calculated. Due to the lack of analytical solutions for calculating mode shapes and natural frequencies, for each beam with specific parameters, the mode shapes and natural frequencies should be first calculated [23], [24].

Consider a single-layer curved beam with cross-sectional area $A$, radius of curvature $R$, moment of inertia $I$, and modulus of elasticity $E$. To derive the vibration equations in the radial and circumferential directions, the element of the curved beam shown in Figure 3 is considered.

Assuming the curved beam as Euler-Bernoulli beam, the compatibility equations for bending moment and circumferential force in terms of radial $w$ and tangential $u$ displacements are expressed as follows [25]:

Considering the element shown in Figure 3 and neglecting the rotational inertia effect of the curved beam, one can write:

Using Newton's second law, the vibration equation of the element in the radial direction can be expressed as:

By substituting Eqs. (1) to (3) into the above relation, the governing equation for the vibration of the curved beam in the radial direction is derived as:

By defining dimensionless parameters as follows:

where, $\omega$ is the natural frequency of the system, the dimensionless governing equation for radial vibrations is then expressed as:

Similarly, by writing the governing equation for the vibration of the element shown in the tangential direction, we have:

Finally, by defining dimensionless parameters, the above governing equation is derived as:

To calculate the mode shapes and natural frequencies of curved beams, the differential Eqs. (7) and (10) must be solved alongside the boundary conditions of the beam as an eigenvalue problem. Consider a single-layer curved beam with one clamped end, the boundary conditions can be expressed as follows [23]:

Using the relations defined in Eqs. (1) to (3), one can write:

The problem in calculating the mode shapes and natural frequencies is to solve an eigenvalue problem. The free vibration response of the dimensionless equations is considered as follows:

By substituting the assumed response forms into the governing differential equations, the sixth-order characteristic equation is obtained as:

By solving the above characteristic equation, the natural frequencies of the curved beams can be calculated. The eigenvalues $s_i(\lambda, k)$ that satisfy the sixth-order equation above must also satisfy the boundary conditions expressed in relations (11) and (12). Explicit parametric relations for the eigenvalues of curved beams do not exist; however, numerical analysis shows that the solution of the third-order equation in terms of the parameter $s^2$ has one negative root $\left(-v_{12}\right)$ and two positive roots $\left(v_{22}, v_{32}\right)$ [23]. With the explanations given, the mode shapes of curved beams can be considered as follows:

It should be noted that the coefficients $C_i, D_i$ are not independent of each other, and their dependence can be expressed by substituting the assumed mode shapes into the differential equations. To determine the natural frequencies and mode shapes of curved beams, the problem must be solved numerically. For this purpose, for different values of $k$ and the assumed angle of the curved beam, the eigenvalues and eigenvectors can be calculated. By placing the assumed mode shapes in the bloundary conditions, to calculate non-trivial responses for the coefficients $C_i, D_i$, the determinant of the coefficient matrix must be set to zero. The calculations related to the eigenvalue problem of curved beams are plotted for a curved beam taken from the standard of the Society of Automotive Engineers (SAE), for values $\alpha=1 \times 10^{-2}$ and the curved beam angle $\beta=0.4 \,\mathrm{rad}$. The changes in the determinant of the boundary condition coefficient matrix according to the dimensionless frequency $\lambda$ are plotted in Figure 4. To validate the results obtained with this procedure, the finite element analysis is also performed using ANSYS software. The mode shapes obtained from the numerical results and the findings using finite element method are shown in Figure 5. According to Table 1, it is observed that the results obtained from the present analysis are in very good agreement with the simulated results using ANSYS.

After calculating the mode shapes of curved beams, the dynamic motion of curved two-layer beams can be studied by deriving the governing vibration equations. Consider the element of a two-layer beam shown in Figure 6, the equilibrium equation for the bending moment of the curved two-layer beam is expressed as:

By writing the vibration equation in the radial direction, one has:

By simplifying the above relation, the following equation is obtained:

Using the compatibility relations of curved beams, the governing vibration equation in the radial direction is expressed as:

By defining the following dimensionless parameters:

and simplifying the above relation, the dimensionless vibration equation in the radial direction of the curved two-layer beam is derived as:

To derive the equations of motion in the axial direction, considering interlayer friction, the equilibrium equation for the bending moment of the upper layer of the curved beam is first written as:

which results in:

Finally, the relationship between the bending moment and the shear force is derived as:

Similarly, for the lower layer of the curved beam, one can write:

The vibrational equation in the axial (tangential) direction for the upper beam, by writing Newton's law and using the relations derived for the single-layer beam, is expressed as:

By substituting Eq. (27) into the obtained equation, one can write:

Finally, using the defined dimensionless parameters, the dimensionless dynamic equation for the upper layer is derived as:

Similarly, for the lower beam, the vibration equation in the longitudinal direction, can be expressed as:

The above PDEs will be handled by Galerkin decomposition method [26] and then will be numerically solved by Runge-Kutta $4^{\text{th }}$ order algorithm. The interlayer frictional force is defined by static friction when there is no interlayer slip. When interlayer slip occurs, it should be defined by dynamic friction which is dependent upon the interlayer pressure and the coefficient of kinetic friction between the layers; therefore, one can write it as:

In the next section, using the obtained governing equations in this section, the different phases of dynamics of two-layer beams in the presence of interlayer friction are identified by studying various parameters.

To validate the results of curved two-layer beams, before qualitative studies, the results obtained using the governing equations and those from the finite element analysis are compared. The finite element model of the curved two-layer beam, as shown in Figure 7, consists of 1620 three-dimensional elements and 10296 nodes. Interlayer friction in the two-layer beam is modeled using 174 three-dimensional contact elements and 170 target elements. In the presence of contact elements, the selection of the normal stiffness parameter starts from higher values until the convergence of the problem is achieved. In case of divergence, smaller values are used.

To validate the numerical results for the vibrations of curved two-layer beams in the present modeling, comparisons are made between the results of radial and circumferential vibrations using the present modeling and finite element results. The time history of circumferential and radial vibrations of the curved two-layer beam with one and three mode shapes and the results related to finite element analysis for are shown in Figure 8. It is clearly revealed that the results of the present analysis using three mode shapes has good agreement with the finite element outcomes. A closer look at the diagrams demonstrates that there is very good agreement in the response trends and the maximum and minimum values for both analyses. The time history of lateral and axial vibrations in subgraphs (a) and (b) of Figure 8 shows that the present modeling results in the axial direction have less agreement compared to the radial displacement. This reduction can be due to factors such as the small amplitude of vibrations in this direction, the definition of time intervals for integration in achieving convergence, and greater sensitivity to interlayer friction.

In previous research, interlayer friction has always been considered as full-slip. Although full-slip modeling simplifies the governing equations for the vibrations of two-layer beams, the results obtained in this case can lead to a misunderstanding of the dynamics of two-layer beams. To demonstrate the importance of accurate modeling of interlayer friction considering the stick-slip effect, the results of the simplified full-slip and the stick-slip models are compared in two cases. In the first case, according to subgraphs (a) and (b) of Figure 9, it is observed that although both types of modeling predict periodic behaviors with a limit cycle, the structure of the stable limit cycle in the mentioned models is completely different. By comparing the results presented in subgraphs (c) and (d) of Figure 9, the difference in predicting the dynamics of the two assumptions becomes more apparent. In this case, modeling the interlayer friction with the assumption of the stick-slip phenomenon predicts chaotic behaviors on two irregular tori, while with the assumption of full-slip, at the same parameter values, the trajectories in the transverse phase plane move towards an attracting limit cycle.

In this section, the dynamic behavior of curved two-layer beams with harmonic excitation is investigated, and different regimes in the dynamics of the system are identified. First, the effect of the excitation amplitude is studied by examining the system's behavior using the Poincaré map. The dynamics of the two-layer beam before and after initial slip shows completely different behaviors. The Poincaré map for $f_0=0.5$ is plotted in subgraph (a) of Figure 10. The regular behavior of the system in this case shows that before the initial slip phenomenon in the curved two-layer beam, the trajectories move on a stable torus, and so the behavior is quasi-periodic. However, this type of dynamic regime has a completely fragile structure, such that with an increase in the external amplitude to $f_0=1$, initial slip occurs, and the trajectories in the phase plane move towards a chaotic attractor. The Poincaré map plotted in subgraph (b) of Figure 10 shows the trajectories after the initial slip. The image of the chaotic attractor in this case is shown more prominently in subgraph (c) of Figure 10. Increasing the excitation amplitude further to $f_0=5$, changes the system's behavior. As seen in subgraph (d) of Figure 10, the trajectories in the Poincaré map move towards a stable fixed point. The stable limit cycle in the system's phase plane in this case is shown in subgraph (e) of Figure 10. With a further increase in the excitation amplitude, according to subgraph (f) of Figure 10, the trajectories in the Poincaré map are damped at a slower rate and eventually move towards a stable toros in the phase plane.

The identification of chaotic regions in the nonlinear vibrations of curved two-layer beams is of particular importance. In chaotic regimes, there are frequency ranges where if the excitation frequency falls within this range, it can lead to an increase in the amplitude of vibrations. The results presented in this research indicate that after initial slip, the dynamic behavior of two-layer beams experiences chaotic regimes. Figure 11 shows the phase portrait and Poincaré map for another set of system parameters after initial slip. The irregularity of the trajectories in the phase space in this case indicates the chaotic nature of the dynamics of curved two-layer beams after initial slip.

When bodies are in contact in the lubricated case, the friction force is not only dependent upon the slip velocity between the surfaces but also on the interaction of the two surfaces at boundary or full film lubrication conditions. In the case of the lubricated surface contact, friction decreases with increasing sliding velocity until the full-film lubrication conditions. After that, friction may remain constant, decrease, or increase depending upon viscosity of the lubricant and temperature conditions. Stribeck model is defined as follows:

where, $F_c$ is the sliding Coulomb friction, $F_s$ denotes the maximum static Coulomb friction, $v_s$ represents the slip velocity coefficient, and $k_v$ is the viscous friction coefficient. The Stribeck friction model diagram is shown in Figure 12.

In the presence of a lubricant between the layers of a curved two-layer beam or moisture in the environment, dry friction does not accurately model the slip behavior of contact surfaces. In this case, by defining a new friction model called Stribeck friction, the characteristics of dry friction and viscous friction in interlayer slip are considered. The effect of dry and Stribeck frictions on the vibrational behavior of a two-layer beam with varying excitation amplitude is investigated in Figure 13. According to the presented results, in this set of parameters, the trajectories move towards stable limit cycle attractors. The characteristic of this dynamic behavior in the Poincaré map is a stable fixed point. The appearance of the limit cycle attractors with dry and Stribeck frictions at smaller values of the parameter $f_0$ has a significant difference. With an increase in the excitation amplitude, the effect of a lubricant on the dynamic behavior of the two-layer beam system decreases.

The results of this research demonstrated that the dynamic behavior of curved two-layer beams experiences chaotic regimes after initial slip. The accurate modeling of interlayer friction, considering the stick-slip phenomenon, is crucial for understanding the dynamic behavior of these structures. The findings of this study provided valuable insights into the design and analysis of multi-layer structures subjected to vibrational loads, particularly in applications where damping and energy dissipation are of great importance.