Vibrational energy transfer in coupled mechanical systems with 1 nonlinear joints

10 This study investigates vibrational transfer and energy flow in nonlinearly coupled systems, each subjected 11 to a harmonic force with different excitation frequency. A nonlinear joint having either smooth or non-12 smooth stiffness characteristics at the coupling interface is considered. The steady-state dynamic responses 13 are obtained by a method of harmonic balance with alternating frequency and time and by a direct numerical 14 integration. The time-averaged transmitted power is used to assess the direction of energy flow and the 15 power transfer between the systems. It is shown that as the excitation frequency ratio increases, the point 16 of zero net power transmission between subsystems move to lower frequencies. The cubic stiffness 17 nonlinearity mainly affects the power transfer in the vicinity of the second resonant frequencies. It is also 18 shown that the second resonant frequencies of both subsystems and the point of zero net power transmission 19 shift to higher frequencies when the bilinear stiffness ratio increases. For the power transfer curves, the 20 bilinear stiffness ratio controls the location of the second resonant frequencies. Findings from this study 21 can provide insights for the design of the joint interfacial properties with regards to vibration transfer in 22 coupled systems under multi-frequency excitations.


Introduction
There is a growing interest in comprehensive understanding nonlinear dynamics of engineering systems and a wide range of nonlinear models has been developed.One important model, the Duffing oscillator, having cubic stiffness term and nonlinear restoring force in the governing equation, is widely used to describe different nonlinear dynamical systems including pendulums [1], beam with permanent magnets [2], cables [3], electric circuits [4] and nonlinear isolators [5][6][7].It was also reported that bolted joints can exhibit nonlinear stiffness property with the nonlinear force being a cubic function of the relative displacement [8,9].Different combinations of the coefficients of linear and nonlinear stiffness terms in the Duffing equation will lead to hardening, softening and double-well potential nonlinearities, causing complex nonlinear phenomena such as super-/ sub-harmonic resonances, internal resonances, multiple response states, bifurcation and chaos.
While many nonlinear systems are characterized by a smooth nonlinear function of the displacement or velocity in their governing equations, some systems behave non-smoothly in terms of restoring force and displacement relationship.A typical example is the so-called smooth and discontinuous oscillator (SD oscillator), which was originally proposed to describe a transition from smooth to discontinuous behaviour [10].The SD oscillator was studied and shown to exhibit complex dynamical behaviour including bifurcations and chaos [11,12].Another example is piecewise linear systems, which can be used to represent typical nonlinear systems with motion constraints [13][14][15], dry friction [16,17], asymmetrical stiffness or damping [18,19], and bolted flange joints [21][22][23].
Many recent studies been devoted to exploiting various types of nonlinearities for performance benefits in vibration suppression [24].For example, nonlinear vibration isolators with geometric nonlinearity can have a high-static-low-dynamic characteristic, providing performance enhancement compared with that of conventional linear isolators [5,6,25,26].The use of nonlinear elements in energy harvesting systems [27] and nonlinear energy sinks (NES) [28][29][30][31] has been studied extensively, for the design objective of achieving optimal output power and targeted energy transfer, respectively.Quinn et al. [27] showed that a nonlinear energy harvester outperforms a tuned linear one with a higher efficiency across a broader frequency range.Much less work has been reported on the use of suppression systems with non-smooth nonlinear characteristics.Wang et al. [33] showed superior suppression performance of a piecewise linear NES compared to a linear vibration attenuation system.
While there are many investigations on the dynamic analysis of nonlinear systems, most of them have primarily considered a single frequency excitation.As a result, the comprehensive understanding of the nonlinear dynamics of coupled systems subjected to multi-frequency excitations remains limited.However, in various engineering scenarios, it is common to encounter multiple excitation frequencies simultaneously.
To illustrate, there is in fact a prevalence of multi-frequency excitations in various engineering applications.
In turbomachinery, the vibration of rotating blades and airfoils can result in an unsteady flow subject to two distinct excitation frequencies [34,35].Additionally, an axially transporting beam can experience twofrequency parametric excitation [36], and a dual-rotor system can exhibit two fundamental excitation frequencies induced by a low-pressure rotor and a high-pressure rotor [37].Multi-frequency excitations are also encountered in microelectromechanical systems, such as microbeams and micromirrors [38,39].These examples highlight the presence of multiple excitation frequencies in various domains of engineering.
Unlike linear systems, nonlinear systems with multi-frequency excitations may exhibit super-/subharmonics and combined resonances [40][41][42].In addition, for nonlinear systems, the principle of superposition cannot be used.In view of this, some research attempts have been made.Guskov et al. [43] explored the multi-frequency dynamical behaviour of a modified Jeffcott rotor system using the multidimensional harmonic balance method (MHBM), alternating frequency-time (AFT) and arc-length continuation.Didier et al. [44] used stochastic-MHBM and polynomial chaos expansion method to investigate the nonlinear vibration of a mechanical system with uncertain material and geometrical parameters.The considered system was subjected to unbalanced forces with incommensurable frequencies, leading to quasi-periodic dynamic response.Zhao et al. [45] studied the nonlinear cable vibration forced by two external periodic excitations.The Galerkin method was used to discretize the governing partial differential equations into ordinary differential equations, and the multiple scale method is further applied to obtain the frequency-response functions.
It is noted that previous studies have focused on the dynamic response of systems, and there have been few attempts made on vibration power and energy transfer analysis of nonlinear systems under multifrequency excitations.Power flow analysis (PFA) is a widely accepted method for assessing vibration and energy transmission level in complex dynamical systems.Its concept was first proposed by Goyder and White [46] and has been further developed to study various linear and nonlinear systems [47][48][49][50][51][52].Zhao et al. [53] studied the power flow transfer in space truss structures using a Timoshenko theory.The active control of the minimum power transmission was found to be more effective and achievable than the control of the minimum acceleration.Xie et al. [54] investigated the vibration transfer and power flow characteristic of a propulsion shaft system in underwater vehicles.The vibration attenuation in a shaft-hull system is quantified and evaluated by power flow and mean square velocity level.In recent years, time-averaged power flow quantities, e.g., input, dissipated, and transmitted powers, have also been used to assess the vibration transmission level in the Duffing oscillator [55] and coupled oscillators with smooth or nonsmooth connections [56][57][58][59].
This study investigates the vibration transfer and energy flow in a coupled system with a nonlinear smooth or a non-smooth interface under multi-frequency excitations.Two external harmonic forces with different excitation frequencies are applied to two subsystems.The smooth joint is characterized by a cubic stiffness spring and the non-smooth connection interface is modelled by a spring of piecewise linear restoring force and displacement relationship.The first-order HB and HB-AFT techniques [60][61][62][63] are used to obtain analytical solutions for dynamic response and related power transmission, and the fourth order Runge-Kutta (RK4) method is used as a numerical approach in the time domain.The rest of this article is organized as follows.Section 2 introduces the physical and mathematical model.Section 3 shows the analytical first-order HB, HB-AFT method, and PFA formulations.Two case studies with the smooth and non-smooth interfaces are demonstrated in Sections 4 and 5, respectively.The conclusions are presented in the last section of this article.

Physical and mathematical modelling
Many engineering systems, such as the engine blade-disk dovetail structure, transmission shaft in ship propulsion system, and satellite separation system with clamp-band-joint comprise jointed components that are subjected to dynamic loading, see in Fig. 1(a).Understanding of the nonlinear dynamics including the vibration transmission within of the jointed structures is important to achieve enhanced design.It was reported that bolted joint will cause nonlinear behaviour and can be approximately described by a smooth cubic stiffness [8,9] or a non-smooth bilinear stiffness model [21][22][23].Fig. 1(b) shows a nonlinear joint with a smooth stiffness nonlinearity characterized by a cubic restoring force term: where   () is the restoring force with a smooth joint,  = where   () is the nonlinear restoring force of the bilinear spring,   and  ℎ are the constant spring stiffness coefficients corresponding to positive and negative relative displacement, respectively.In Fig. 1(d), each structure can be further characterized by a SDOF linear subsystem representing the dominant mode.Therefore, the original physical model is simplified as a coupled oscillator system with smooth or non-smooth joint.Subsystem one (S1) consists of a mass  1 subjected to a harmonic force  1 cos  1 , a linear spring with stiffness coefficient  1 , and a viscous damper with damping coefficient  1 .
Subsystem two (S2) comprises a mass  2 with another external force  2 cos  2  attached to a linear spring  2 and a viscous damper  2 .Both masses oscillate horizontally, and the static equilibrium position is taken as reference at which the displacements,  1 =  2 = 0.The equation of motion of the integrated system is where () is the coupling force at the interface, replaced by   () in the case of the smooth nonlinear joint, and by   () or   () for the cases of non-smooth joint.New parameters and variables are introduced below to facilitate dynamic analysis ,  =  10 , where  10 and  20 are the undamped natural frequencies of subsystems one and two, respectively,  is the frequency ratio between them,  is the mass ratio,  0 is the un-stretched length of the spring on the left,  1 and  2 are the non-dimensional displacements of masses  1 and  2 , respectively, ∆ is the nondimensional relative displacement between the masses,  1 and  2 are the non-dimensional damping coefficients,  1 and  2 are the non-dimensional forcing amplitudes, Ω 1 and Ω 2 are the non-dimensional fundamental excitation frequencies,  is the non-dimensional time.By using these dimensionless parameters, the governing equation ( 4) can be written into a non-dimensional form where  = { 1 (),  2 ()} T is the displacement response vector, the primes (′) denote differentiation operations with respect to the non-dimensional time , the symbol "T" denotes taking the transpose of a matrix,   () = { 1 e iΩ 1  ,  2 e iΩ 2  } T denoting the external load vector, , , and  represent the mass, damping and stiffness matrices of the system with and   (,  ′ , ) = {−(∆, ), (∆, )} T represents the force vector generated at the nonlinear joint.Here the excitation frequency ratio is defined as the ratio of the two excitation frequencies Ω 1 and Ω 2 : For the joint with a cubic stiffness nonlinearity, we have where  =   / 1 and  =    0 2 / 1 , representing linear and nonlinear stiffness ratios of the smooth joint, respectively.When the joint is characterized by a piecewise linear spring, the corresponding dimensionless restoring force is where  =   / 1 is the stiffness ratio,  =   /  is the piecewise linear stiffness ratio and  = / 0 is the non-dimensional offset.When the joint is characterized by a bilinear spring, the corresponding dimensionless restoring force is where  =  ℎ / 1 is the stiffness ratio and  =    ℎ ⁄ is the bilinear stiffness ratio.
To examine the vibration transmission and energy flow through the nonlinear joint of the coupled system, it is necessary to solve the nonlinear governing equations.In this study, two different approaches will be adopted.One is the harmonic balance (HB) method based on analytical derivations and the alternating frequency time (AFT) scheme.The other is based on a fourth-order Runge-Kutta (RK) method.
The HB-AFT method has been a widely accepted tool to obtain the periodic responses of a dynamical system and it can provide physical insights into the dynamics of nonlinear systems.The RK method can be used to obtain both periodic or non-periodic responses with high accuracy but at higher computational cost.

HB-based vibration energy flow analysis
In this section, HB-based vibration energy flow analysis is presented.A general approach employing the HB-AFT is introduced to obtain the steady-state response solution of Eq. ( 4).Analytical method using firstorder HB approximations of smooth joint case is also presented.Multiple performance indices, such as time-averaged input and transmitted power, are defined and formulated.

HB-AFT for multi-frequency excitations
The HB-AFT technique is used to obtain the periodic responses of the coupled systems with a nonlinear joint [60][61][62][63].For its implementation, the general solution of Eq. ( 4) can be truncated into -th order Fourier series where =1 or 2 represents the subsystem S1 or S2;  ̃(,) and  ̃(,) are the complex Fourier coefficients of the dimensionless displacement for the -th harmonics associated with excitation frequencies Ω 1 and Ω 2 , respectively; ℜ denotes the operation of taking the real part of a complex number.The corresponding velocity and acceleration are respectively.The nonlinear force applied to the nonlinear joint can be expressed as where  ̃ and  ̃ are the complex variables of nonlinear force with -th harmonics.By inserting Eqs ( 11)-( 14) into Eq.( 4) and balancing the harmonic coefficients at the n-th order, one obtains where  ̃ = { ̃(1,) ,  ̃(2,) } T ,  ̃ = { ̃(1,) ,  ̃(2,) } T ,  ̃ = { ̃, − ̃} T ,  ̃ = { ̃, − ̃} T ,  ̃1 = { 1 , 0} T , and  ̃2 = {0,  2 } T .It is noted that Eqs ( 15) and ( 16) are two nonlinear equations with complex numbers, which can be transformed into four real algebraic equations.Therefore, for the coupled two-DOF system with N-th order harmonics, the total equations will be 2(4 + 2).The solutions of these nonlinear algebraic equations can be obtained by the Newton-Raphson based pseudo arc-length continuation techniques [64][65][66].

Analytical HB approximation
The previous section provides a general procedure to obtain the dynamic response and power flow variables for nonlinear systems with smooth or non-smooth joint based on the HB method.This approach is mainly based on the Fourier Transform and numerical continuation technique, which has sufficient accuracy but relatively large amount of calculation.For a smooth joint, e.g., cubic stiffness nonlinearity, the analytical first-order harmonic balance (HB) approximation can also be used to obtain the dynamic response effectively and efficiently.The steady-state dimensionless displacement of S1 and the relative displacement of the subsystems are expressed by respectively, where , , , and  are the response amplitudes,  1 ,  2 ,  1 and  2 are the corresponding phase angles.Based on Eqs ( 17) and ( 18), the first and second derivatives of the displacements with respect to time can be calculated.By substituting related displacements, velocities and accelerations into governing Eq. ( 4), ignoring high-order terms, and balancing the coefficients of terms cos(Ω 1 ), sin(Ω 1 ), cos(Ω 2 ), and sin(Ω 2 ) , we can obtain eight nonlinear algebraic equations with eight unknowns of response amplitudes and phase angles.See details in the Appendix.They can be solved by a standard Newton-Raphson technique together with the numerical continuation algorithm scheme.
For later analysis, the natural frequencies of the corresponding linear undamped system are determined.
By setting  1 =  2 =  1 =  2 =  = 0, Eq. ( 4) becomes where first-order approximations are used, | 1 | and | 2 | are the response amplitudes of the S1and S2, respectively.The natural frequencies are determined by setting the determinant of the matrix to be zero Eq. ( 20) is the characteristic equation to which the corresponding solutions are the linearized natural frequencies Ω 1 and Ω 2 (assuming Ω n1 < Ω n2 ) .By solving the quadratic equation of Ω 2 , we have

Vibration energy flow quantities
Vibration power flow and energy variables are widely used to evaluate the level of vibration transmission for dynamical systems.In this study, the input and transmitted powers are of interest.
The instantaneous input power is the sum of power injection from external sources in each subsystem, that is, the total energy or power consumption within the system due to the viscous damping according to the law of energy conservation, which can be expressed as where  in1 and  in2 are the instantaneous input power of subsystem one and two, respectively;  1 ′ and  2 ′ are the velocities based on Eq. ( 12) of S1 and S2, respectively.For steady-state motion, the dimensionless time-averaged input power is where  0 is the starting time of integration and   is the averaging time; ( * ) denotes the complex conjugate of a complex number.Starting time is set as  0 = 800 to remove the transient motion, where  = 2/Ω 1 .
Averaging time is   = 1000.The expression of the time-averaged input power obtained using the firstorder HB approximation is provided in the Appendix.
The instantaneous transmitted power is defined as the product of the nonlinear transmitted force and the velocity of subsystem two, representing the power transmission between the two subsystems through the nonlinear smooth/non-smooth joint.Therefore, the non-dimensional instantaneous transmitted power can be expressed as According to the energy conservation law, the sum of the transmitted power  t and the dissipated power  d2 in S2 equals the total power input  in2 .Therefore, the time-averaged transmitted power in steady-state can be written as It should be mentioned that the positive value of the time-averaged transmitted power represents the vibrational power flow and energy transmission from subsystem one to two, which means that subsystem one has higher energy potential than subsystem two, and vice versa.Note that the expression of the timeaveraged transmitted power using the first-order HB method is shown in the Appendix.

Results and discussion
In this section, the dynamic response and vibrational energy transfer of coupled systems with a smooth or a non-smooth joint is presented in section 4.1 and 4.2, respectively.The effects of the excitation frequency ratio  = Ω 2 /Ω 1 , the piecewise linear stiffness ratio  =   /  and the bilinear stiffness ratio  =    ℎ ⁄ on the response amplitude and power flow quantities are examined.

Vibration transmission through smooth joint
Here, the two subsystems are connected by a nonlinear smooth joint with cubic stiffness nonlinearity described by Eq. ( 8).The effects of the excitation frequency ratio on the dynamic response and vibration transmission are investigated and analyzed.
Figure 2 shows the impact of the excitation frequency ratio  on the resonant peaks.Based on the free vibration analysis of the undamped system, the natural frequencies Ω 1 and Ω 2 , as derived from Eqs (21) and ( 22) respectively, are denoted as peaks M and N in the figure.Furthermore, when subsystem two is excited by a frequency Ω 2 (i.e., Ω 1 ), two additional resonances emerge, represented by Ω 1 / and Ω 2 /, and labeled as peaks P and Q respectively.Consequently, the frequency-response curves exhibit a total of four resonant peaks, with peaks M and P corresponding to in-phase motions, while peaks N and Q correspond to out-of-phase motions.It is noted that the frequencies of peaks M and N do not change despite of the variations of the excitation frequency ratio .In the case of system parameters set as  =  =  = 1, peaks M and N are located at Ω 1 = 1 and Ω 2 = √3.In comparison, the other two peak frequencies vary with .Therefore, from low to high frequencies, there are four possible orders of the appearance of the peaks: Type-1: MNPQ, Type-2: MPNQ, Type-3: PMQN, and Type-4: PQMN (sequence PMNQ is not applicable for the current case) depending on the fixed excitation frequency ratio, as shown in Figs 2(a), (b), (c) and (d), respectively.For example, in Fig. 2(a), Type-1 with the frequency ratio  = 1/3, the resonant peaks M, N, P, and Q are located at Ω 1 = 1, Ω 1 = √3, Ω 1 = 3 and Ω 1 = 3√3, respectively.Fig.
2 also shows the influence of the cubic stiffness nonlinearity  on the dynamic response as compared to the reference linear joint case.The response curves of two cases almost coincide at the peaks M and P, meaning that these two resonant peaks remain unchanged with the variations in  value.Due to the hardening stiffness nonlinearity  = 0.5, two peaks N and Q are bent to the high-frequency range, and the jump phenomena and multiple solutions also occur.Because of the large displacement motion near the resonance, the linearization fails in accurate prediction of the dynamic response.The stiffness nonlinearity has major effects in the vicinity of the second resonant peaks.In the high-and low-frequency ranges, the response curves for difference cases merge, indicating that the effects of the frequency ratio and the stiffness nonlinearity are negligible here.This is because that the relative displacement between the subsystems is relatively small in these regions, so that the nonlinear restoring force due to the nonlinearity of the joint is low compared to the linear term.Due to the multi-frequency excitations, there are two primary frequency components.Apart from the primary response component at Ω  = Ω 1 , the other frequency component is related to the excitation frequency ratio, e.g., Ω  = Ω 2 = Ω 1 /3 when  is 1/3 and Ω  = Ω 2 = 3Ω 1 with  being 3.In other words, the instantaneous dynamic response only contains first-order frequency components, and there are no obvious super-or sub-harmonic components.However, at point B', the response only shows a major frequency component at Ω  = Ω 1 while the expected frequency component of Ω  = 3Ω 1 disappears, as shown in Fig. 2(d).The reason is that the excitation frequency of Ω 2 at point B' is away from the resonant peaks P and Q, and the corresponding influence on the dynamic response is negligible.Therefore, the fundamental excitation frequency Ω 1 is dominant at high frequencies for peak Type-4 (i.e.peak sequence PQMN).It demonstrates that the frequency spectra results are highly related to the excitation frequencies Ω 1 , Ω 2 and their relative ratio.Similar phenomena can also be observed in Fig. 2(a) for peak Type-1 (i.e.peak sequence MNPQ), the spectrum shows a stronger frequency component at Ω  = Ω 1 for point A as it is close to the resonant peak M. As for point B, it is near the resonant peak Q and away from the resonant M, therefore, the frequency component at Ω  = Ω 1 /3 is higher than Ω  = Ω 1 .In Figs 3(a) and (b), the dynamic responses associated with points A and A′ are obtained from the HB and RK methods and shown in the time domain.The figure shows that the analytical results using the first-order HB approximations agree well with the direct numerical integration results.Hence, with a balanced consideration of computational efficiency and accuracy, the analytical first-order approximation is used in this section.Fig.
3 also shows that with the frequency ratio of 1/3 and 3, the displacement responses are periodic having periods  0 = 3 and  0 = , respectively, where  0 is one oscillations cycle and  = 2/Ω 1 .
Figures 4(a) and (b) show the effects of the excitation frequency ratio on the relative displacement amplitude of  and the time-averaged input power  ̅ in , respectively.The appearance of the peaks is in sequence of PQMN, i.e., Type-4.Fig. 4(a) shows that there are only two right-bending peaks N and Q in the response curve of , while there are no primary resonance peaks M and P.This behaviour is related to the fact that the coupled subsystems exhibit in-phase motion at the two peak frequencies M and P, and outof-phase motion at N and Q.Four peaks of similar heights are observed in the curves of  ̅ in .Fig. 4 also shows that the two resonances related to excitation frequency Ω 2 (i.e., peaks P and Q) shift to lower frequencies as the frequency ratio  increases.When  increases from 3, to 5, and then to 7, the peak P moves from Ω 1 =1/3 to 1/5 and then to 1/7.In comparison, the peaks M and N remain unchanged regardless of the variations in the excitation frequency ratio.Fig. 4(b) shows that in the low-frequency range, there is a higher level of total power input when the system has a larger frequency ratio.This is because the excitation frequency Ω 2 dominates at low frequencies for peak Type-4 (peak sequence PQMN).Away from the low frequency range, e.g., resonant area around peaks M and N as well as the high-frequency range, the influence of the frequency Ω 2 is weakened, and the excitation frequency Ω 1 plays a major role, with the lines for different cases merge.The behaviour of vibrational energy transfer within the coupled system under different excitation frequency ratios is investigated and the results are shown in Fig. 5.The results indicate that when the excitation frequency is high or close to peaks M and N, the time-averaged transmitted power  ̅  is positive.
This signifies a net power flow from subsystem one to subsystem two through the smooth nonlinear interface.Importantly, there exists a critical frequency at which the power transmission curve changes sign, resulting in zero net energy transfer.Beyond this critical frequency, power starts flowing in the opposite direction, indicating that subsystem two possesses a higher energy potential, and power transfers from subsystem two to subsystem one.In Figure 5, this critical frequency is approximately Ω 1 ≈0.687, 0.468, and 0.341 for the cases where  equals 3, 5, and 7, respectively.As the frequency ratio increases, the peaks P and Q move to the low-frequency range, allowing the critical frequency to reach a new equilibrium point.For peak Type-4 shown in Fig. 5, the fundamental excitation frequency Ω 1 of subsystem one has a major influence around resonances M and N as well as at high frequencies, so subsystem one has higher energy potential in these areas.In comparison, the fundamental excitation frequency Ω 2 of subsystem two controls the power transmission in the lowfrequency range, that is, subsystem two has higher energy potential in this region.Figs 5(a)-(d) further show the dynamic response behaviour of two positions (Ω 1 = 0.2 and Ω 1 = 8) for different cases using phase portrait diagrams.The results in Fig. 5(a) suggest that in the high-frequency region, the excitation frequency ratio has little effect on the transmitted power.This is reflected in the fact that the phase portrait shows only one periodic solution, regardless of the value of .This indicates that the system behaviour is relatively insensitive to the changes in the excitation frequency ratio in this region.In contrast, the results in the lowfrequency region show that the system behaviour is much more sensitive to the changes in the excitation frequency ratio.The phase portrait can show multiple types of solutions, such as two periodic solutions, quasi-periodic solutions, and multi-periodic solutions.This means that changes in the excitation frequency ratio can significantly impact the transmitted power, leading to different system behaviour.
In Fig. 6 the mechanism of power transmission through the nonlinear smooth interface under multifrequency excitations is further investigated.The dotted line represents the time-averaged transmitted power with frequency ratio  = 1/3 (peak sequence MNPQ, i.e., Type-1), and the solid line denote the peak Type-3 (peak sequence PMQN) with frequency ratio  = 3/2.For Type-1, the two resonant peaks of the excitation frequency Ω 2 (peaks P and Q) are in the high-frequency with negative value.The corresponding equilibrium point of power transmission is around Ω 1 ≈ 2.061.For Type-3 case shown in Fig. 6, the downward arrows represent the net power flows from subsystem one to two, and the upward arrow indicates the energy transmission in the opposite direction.In the two frequency ranges Ω 1 ≈ 0.10 to 0.861 and Ω 1 ≈ 1.102 to 1.253, the time-averaged transmitted powers are negative.Combined with Fig. 6, it shows that  ̅ t has a negative value in the vicinity of the resonant peaks P and Q, and  ̅ t is positive near the resonance areas of peaks M and N. It indicates that the frequency ratio has a significant effect on the location of peaks P and Q and the direction of energy transfer in coupled vibration system, while other regions have negligible effects.Fig. 6 provides a potential control method for power transmission by using different excitation frequencies.In addition, the influence of the stiffness nonlinearity of the joint and the forcing amplitude is also considered, and results are presented in the Appendix.

Vibration transmission through non-smooth joint
This section explores the dynamic responses and vibrational energy transfer within the coupled systems, facilitated by a non-smooth joint.Two models are considered: one featuring piecewise linear stiffness (illustrated in Fig. 1c), and the other employing bilinear stiffness (depicted in Fig. 1d).In order to obtain accurate and efficient results for the dynamic response and power flow variables, the seventh-order HB-AFT method is utilized.This section also aims to assess the influence of the stiffness ratio and excitation frequency ratio on the vibration transmission.
Figure 7 examines the impact of the stiffness ratio, represented by the piecewise linear joint  =   /  , on the relative response amplitude  and the time-averaged input power  ̅ in .The findings highlight that the ratio of the two slopes primarily influences the dynamic response and power transmission in the secondary resonances.An increase in the  value, indicating a higher stiffness, results in the secondary resonant peaks shifting towards higher frequencies.Additionally, the study reveals that bending and discontinuity occur in the response and power transmission curves when the relative displacement amplitude surpasses the offset deformation .When the  value is small, i.e.,   <   , a right-bending is observed in the secondary resonant peaks, reminiscent of hardening behaviour.Conversely, when  > 1, indicating   >   , the secondary resonant peaks bend towards the low-frequency range, similar to softening behaviour.employs an analytical first-order HB method (represented by solid lines).The study demonstrates that the piecewise linear stiffness joint can be effectively approximated by a smooth polynomial function with a cubic term, yielding a satisfactory agreement between the two approaches.Furthermore, the study reveals that increasing the stiffness ratio  (where  =   / 1 ) causes the secondary resonant peaks to shift towards higher frequencies.Additionally, a larger value of  leads to a lower level of the relative response amplitude in the low-frequency range.However, in the high-frequency range, the impact of the stiffness ratio  becomes negligible as the frequency-response curves for each case coincide with each other.The results shown in Fig. 9 demonstrate the influence of the bilinear stiffness ratio  on the response amplitude of the relative displacement Y, considering the bilinear stiffness model.The figure compares the HB-AFT results with the numerical integration results obtained from the Runge-Kutta method.As the bilinear stiffness ratio increases, the two resonant peaks shift to higher frequencies, indicating that the natural frequencies of the system increase with the stiffness.In the low-frequency range, the response amplitude is higher for a smaller bilinear stiffness ratio, while in the high-frequency range, the effect of the bilinear stiffness ratio on the response amplitude is not significant.It is noted that the equivalent stiffness using the linearization method can provide a good estimation of the dynamic response, especially in the resonant area and high-frequency range.However, for the low-frequency range, it may lead to an underestimation of the response amplitude, and important dynamic information such as super-harmonics and quasi-periodic motions will not be captured.In comparison, the use of the seventh-order HB-AFT method enables the detection of super-harmonics at low frequencies and the corresponding results are in good agreement with the numerical integration results.at  1 = 1 and  1 ≈ 1.557, respectively.The results also show that in the low-frequency range, the response amplitude and time-averaged input power increase with the excitation frequency ratio.However, in the high-frequency range, the effect of the frequency ratio is insignificant.Combined with the previous findings in Fig. 9, it suggests that the bilinear stiffness ratio  primarily affects the dynamic response and energy transmission curves around the second resonant frequencies, while the frequency excitation ratio ε is responsible for the resonant peaks of  2 .This provides potential methods for controlling and mitigating vibrations and power in nonlinear systems subjected to multi-frequency excitation.12(a) show that the critical frequency of zero net power transmission between the subsystems increases with the bilinear stiffness ratio .It is observed that the two primary resonant peaks of the timeaveraged transmitted power remain unchanged, but the peaks N and Q shift to higher frequencies as the bilinear stiffness ratio increases.Furthermore, the results suggest that a greater value of bilinear stiffness ratio leads to a higher level of power transmission in the high-frequency range, but has little effect on  ̅ t in the low-frequency range.The results in Fig. 12(b) indicate that the frequency bandwidth of positive  ̅ t increases with the frequency ratio .It is found that the point of zero net power transmission shifts to lower frequencies as the frequency ratio increases.These results show that changing the bilinear stiffness and frequency ratios can be a potential way of energy mitigation and targeted energy transmission in nonlinear dual-excitation systems.

Conclusions
In this study we investigated the dynamic response and the vibrational energy transfer between coupled systems subjected to different excitation frequencies.The two subsystems are connected by a nonlinear cubic stiffness, a non-smooth piecewise linear or a bilinear stiffness joint.The first-order HB and the seventh-order HB-AFT methods were used as analytical approximations.The numerical fourth-order Runge-Kutta method was also employed for validation and comparison.The time-averaged input and transmitted powers were used to assess the energy transmission performance.
For the system with the smooth joint, the resonant peaks, caused by fundamental frequency Ω 1 , do not change with the variation of excitation frequency ratio ε.However, the other two peak frequencies change with ε, leading to four possible orders of appearance of the peaks depending on the value of ε.It was found that the dynamic response of each subsystem only contains fundamental excitation frequencies without obvious sub-/super-harmonics.It was also demonstrated that the cubic stiffness nonlinearity mainly affects the vicinity of the second resonant peaks of the energy transmission curves.
For the system featuring a non-smooth joint, the piecewise linear stiffness ratio induces a bending behaviour in the second resonant peaks, resembling either hardening or softening characteristics.In the case of the bilinear stiffness joint, it has been observed that increasing the bilinear stiffness ratio leads to a shifting of the second resonant peaks and the point of zero net energy transmission towards higher frequencies.The excitation frequency ratio exerts a significant influence on the peak frequencies of subsystem two in terms of energy transmission, while the bilinear stiffness ratio primarily affects the characteristics of the second resonant peaks.Moreover, various nonlinear phenomena such as periodic-1, periodic-2, periodic-3, and super-/sub-harmonic resonances have been identified in the system's response.
For both two cases, it is shown that the direction and the amount of time-averaged transmitted power through the nonlinear joint can be tuned by adjusting the excitation frequency and bilinear stiffness ratios, and thus achieving better dynamic performance.The resonances of subsystem two move to lower frequencies as the increase of the excitation frequency ratio.In the low-frequency range, the dynamic response and energy transmission level decrease with the excitation frequency ratio, which is beneficial for vibration suppression.

Fig. 1
Fig.1 Schematic diagram of coupled structures with bolted joint, subject to different excitations  1 e i 1  and  2 e i 2  .(a) Physical model of a bolted joint with dynamic loading.It may exhibit following force-deformation characteristics: (b) a smooth joint with linear stiffness coefficient   and nonlinear stiffness coefficient   ; (c) a non-smooth joint with piecewise linear stiffness coefficients   and   ,  is the offset deformation; (d) a non-smooth joint with bilinear stiffness coefficients  ℎ and   .

Fig. 2
Fig. 2 Effects of the excitation frequency ratio  on the sequence of the resonant peaks computed for (a)  = 1/3, (b)  = 3/5, (c)  = 3/2, (d)  = 3. Peaks M and N are the resonances due to excitation frequency Ω 1 , and peaks P and Q are resonances due to excitation frequency Ω 2 .Frequency spectra diagrams are for the points located at Ω 1 = 0.2 (A and A') and Ω 1 = 8 (B and B').Solid lines: nonlinear stiffness at the joint ( = 0.5).Dashed lines: linear stiffness

Figure 2
Figure 2 also contains the frequency spectra information of four points as marked by A, A', B and B'.

Figure 8
Figure 8 depicts the frequency-response curve of the relative response amplitude Y obtained using two different methods.The first method employs a seventh-order HB-AFT approach (represented by dashed lines) to characterize the dynamic response of the system featuring a piecewise linear stiffness joint.The second method approximates the non-smooth joint by utilizing a smooth cubic stiffness function and

Figures 11 (
Figures 11(a) and (b) show the effects of the excitation frequency ratio  on the relative displacement amplitude and the time-averaged input power of the coupled system with a bilinear stiffness joint.The results indicate that as the frequency ratio increases, peaks P and Q, which are resonances caused by the fundamental frequency  2 , shift toward lower frequencies.In contrast, peaks M and N remain nearly fixed

Figures 12 (
Figures 12(a) and (b) analyze the influence of two parameters, the bilinear stiffness ratio  and the excitation frequency ratio , on the time-averaged transmitted power  ̅ t between two subsystems.The results in Fig. 12(a) show that the critical frequency of zero net power transmission between the subsystems
2 −  1 is the relative displacement between two subsystems,   and   are the linear and nonlinear stiffness coefficients of the smooth joint, respectively.It is worth noting that the first derivative of the nonlinear force   () with respect to relative displacement  leads to a linear stiffness at the original equilibrium position, i.e., 1  } + ℜ{ 2 ′}ℜ{ 2 e iΩ 2  } = (ℜ{∑ iΩ 1  ̃(1,) e iΩ 1