Vibration reduction and stability analysis of damping rings on nonlinear free-free beam

Two damping rings (DRs) are added on a nonlinear free-free elastic beam to reduce vibrations. This elastic beam with free-free boundary conditions can be used to simulate the high speed aircraft fuselage vibrations. Both isotropic beam and compound beam with double-section are investigated. The beam is subjected to a distributed load with the harmonic wind force. The method of multiple scales ((MOMS) a perturbation technique) is employed to analyze this nonlinear problem in frequency domain. DRs with various locations, masses, spring constants, and damping coefficients are analyzed and the maximum amplitude counter plots (MACPs) are compiled to identify DRs’ combinations for the optimal damping effects. The Floquet theory is used to make the basin of attraction (BOA) charts and analyze the stability of the system. This research can be applied in a wide range of engineering problems, such as rocket structures, high-speed aircrafts, satellites, missiles, etc.


Introduction
Studies of vibrations have always been a concern for researchers and engineers because it may cause the structural fatigue or failure. Beams are widely applied in engineering problems, such as wings in aerospace engineering, bridges in civil engineering, and train rails in mechanical engineering, etc. Many studies on beam vibrations have been performed. Ö zkaya 1 considered a beam-mass system under simply supported end conditions. The effects of positions, magnitudes, and the number of the masses were investigated. Mundrey 2 considered the 2D Euler-Bernoulli beam resting on an elastic foundation to simulate the railway track. A moving load on the 2D beam was studied. Chang 3 studied nonlinear vibrations in carbon nanobeams under magnetic field. The analytical study was performed by using He's variational method. Wang and Wei 4 used a fluidconveying nonlinear beam and nonlinear springs to simulate the vibration of a fluid-conveying tube placed on an elastic foundation. The internal resonance was found under a certain combination of fundation elastic modulus and flow speed. In most nonlinear beam problems, the internal resonance is a major point of discussion. Due to nonlinearity, internal resonance generally occurs in modes that are not being directly excited by external forces. Exciting higher modes can lead to high amplitude vibrations in lower modes. 5 What is interesting is that in common 3D beams with symmetrical cross sections; one-to-one internal resonance is the most likely to take place among the various degrees of freedom. As its resonant frequency is the same with each other, it is also called primary resonance. Pai 6 analyzed the primary resonance in a 3D nonlinear composite rotating beam. Stoykov and Ribeiro 7 examined the stability of a 3D nonlinear rotating beam based on Timoshenko's theory and took into account the deformation caused by twists and warps. They used Floquet theory to analyze the stability of the system and determined that the primary resonance produces supercritical symmetry-breaking bifurcation in beams with a square cross-section. Chang and Wang 8 studied the primary resonance in a 3D nonlinear free-free beam. Tekin et al. 9 considered the three-to-one internal resonance in 2D multiple stepped beam systems with clampedclamped, pinned-pinned, and clamped-pinned supports. Wang and Chang 10 analyzed the primary and internal resonance in 3D free-free double-section beam.
The way to effectively suppress vibrations is always a problem to engineers. Wang and Kuo 11 examined the vibrations in a nonlinear hinged-free beam resting on a nonlinear elastic foundation. They discovered that with a certain spring constant in the elastic foundation, threeto-one internal resonance occurs in the first and second modes of the system. They prevented internal resonance and reduced vibrations by adding a tuned mass damper (TMD) on the beam. Their results revealed that the best damping effects could be obtained when the tuned mass damper was placed between 0.25l and 0.5l from the fixed end of the beam. Wang and Liang 12 studied the damping effects of a lumped-mass vibration absorber on a hinged-hinged nonlinear beam resting on a nonlinear elastic foundation. Using the maximum amplitude contour plots (MACPs), they obtained the optimal mass and spring constant combination of the lumped-mass vibration absorber for the best vibration reduction. Wang et al. 13 also found that when the multiple value of mass and spring constant (m D k) of the tuned mass damper equals to a certain value (m D k = 0.0475), the best vibration reduction effect can be obtained. Wang and Hsiao 14 studied the off-shore wind turbine tower by using the 3D nonlinear multi-loaded slender beam to simulate the vibration of the turbine tower. They added two damping rings (DRs) to prevent internal resonance and suppress vibration. They concluded that adding one of the damping rings on the top of the turbine tower and the other one on the ocean surface of the tower will produce better damping effects. These studies demonstrate that changing the location, mass, spring constant, and damping coefficient of the tuned mass damper are feasible approaches to prevent internal resonance and mitigate vibrations.
This study considers a nonlinear free-free beam subjected to a distributed load with the harmonic wind force. The isotropic beam and compound beam with double-section are investigated. Two damping rings are added on the beam. Damping rings with various locations, masses, spring constants, and damping coefficients are analyzed and the maximum amplitude counter plots are compiled to identify damping rings' combinations for the optimal damping effects. This study also uses Floquet theory to analyze the stability of this flow-structure coupled system. The results are verified by numerical method (the 4th order Runge-Kutta method) and ANSYS simulations.

Isotropic free-free beam model
This study analyzes the vibrations on a straight 3D free-free nonlinear beam, which can be simulated to high speed aircrafts such as rocket, missiles, etc. The coordinate definitions of the beam and the relationships between various external forces and damping rings are presented in Figure 1. Three consecutive Euler angles are used to relate the deformed and undeformed states. The equations contain structural coupling terms and quadratic and cubic nonlinearities due to curvature and inertia. Based on Newton's second law, Euler's angle transformation, and Taylor series expansion, the dimensionless equations of motion of the nonlinear beam can be expressed as follows 11,12,15 (please see Appendix A for details): € y + m zy y iv =À c y _ y À m xy (g 0 z 00 ) 0 À m zy ½y 0 (y 0 y 00 + z 0 z 00 ) 0 0 À (1 À m zy )½(y 00 g 2 À z 00 g) 0 Figure 1. Schematic model of an isotropic free-free beam with multiple forces and damping rings.
€ z + z iv =À c z _ z + m xy (g 0 y 00 ) 0 À ½z 0 (y 0 y 00 + z 0 z 00 ) 0 0 + (1 À m zy )½(y 00 g + z 00 g 2 ) 0 + y 000 where y and z are beam displacements in the y-and zdirections, respectively. Those terms with a '' . '' on the upper or upper right denotes d/dt, and with a ''#'' denotes d/dx. In which c y,z represent the dimensionless structural damping coefficient in the y-and z-directions, respectively. The mass moments of inertia are denoted as j y,z in the y-and z-directions, respectively. They are non-dimensionalized by beam mass (m) and length (l). F y,z are the dimensionless forces in the y-and z-directions, respectively. The dimensionless flexural rigidities of the elastic beam are defoined as m xy =D xx / D yy , m zy = D zz /D yy . In which D xx =GI x , D yy =EI y and D zz =EI z , respectively for isotropic beam. G is elastic shear modulus, and E is Young's modulus. I x,y,z represent the area moment of inertia in the x-, y-and zdirections, respectively. Assuming the beam is subjected to a distributed load with the harmonic wind force q y, z e iO y, z t and its associated unsteady aerodynamic force F D,A , that is, In which q y,z and O y, z represent the magnitude and frequency of the harmonic wind force in y-and z-direction, respectively. The unsteady aerodynamic force is adopted from van Horssen and Boertjens 16,17 and are non-dimensionalized as: 2mv 2 . are aerodynamic coefficients in the yand z-directions. In which r a is air density, U y, z are wind speeds in the y-and z-directions, d is diameter of the beam (d is the dimensionless value of d), a D0U , a A0U are lifting coefficients, and a D1U , a A1U and a D2U , a A2U are the first-order derivative and second-order derivative modified terms, respectively. It is also noted that both a D0U and a A0U are zero for non-rotating body, both a D1U and a A1U are 2p, and the second-order derivative modified terms (a D2U , a A2U ) are not considered in this approach.
Damping rings can be considered as external forces produced by two tuned mass dampers in the y-and zdirections of the beam. The influence of damping rings will be discussed detaily in later sections. The equation of motion of the damping ring is derived by using Newton's second law as follow: m Dr € y Dr ( t)À f s ( y( x, t)À y Dr ( t))+ g s ( _ y( x, t)À _ y Dr ( t)) The index r = 1,2 represents two damping rings. Replacing y Dr with z Dr in equation (6) The concept of the damping ring is originally designed to reduce the vibration of space rocket in this work. One can use the rubber, silicone or any other suitable materials to make the damping ring. The material properties ( g s and f s ) can be got by Using Newton's second law, the integrated dimensionless equations of motion of the nonlinear beam can thus be written as € y + m zy y iv =À (c y Àâ D1U )_ y À m xy (g 0 z 00 ) 0 À m zy ½y 0 (y 0 y 00 + z 0 z 00 ) 0 0 À (1 À m zy )½(y 00 g 2 À z 00 g) 0 € z + z iv =À (c z Àâ A1U )_ z + m xy (g 0 y 00 ) 0 À ½z 0 (y 0 y 00 + z 0 z 00 ) 0 0 + (1 À m zy )½(y 00 g + z 00 g 2 ) 0 + y 000 where l Dr (r = 1,2) indicates the location of the first and second damping ring on the beam, and dðÞ is Dirac delta function.

Double-section free-free beam model
This section considers the compound structural problem. The schematic of the double-section free-free beam is shown in Figure 2. Two damping rings are added on the double-section beam. Similarly, the equation of motion of damping rings can be obtained by using Newton's second law. The dimensionless equations of motion of this double-section beam can be written as: € y (i) + m zyi y iv (i) =À (c y Àâ D1U )_ y (i) À m xyi (g 0 (i) z 00 (i) ) 0 À ½y 0 (i) (y 0 (i) y 00 (i) + z 0 (i) z 00 (i) ) 0 0 The boundary conditions and the compatibility equations on the joint of two beam sections are where h = l 1 = l 1 + l 2 ð Þ.

Systematic analysis of damping rings
The results from Wang and Chang 10 indicated that when y-direction is excited and with no dampers, the primary resonance would occur in the y-and z-directions. The three-to-one internal resonance will occur in the beam's 1st and 2nd modes. Therefore, two damping rings are equipped on the beam in order to prevent internal resonance and reduce the vibrations. This study encompasses the frequency responses of the beam system with damping rings in the 1st and the 2nd modes. The fixed points plots facilitate the analysis of the influence of damping rings on internal resonance in the system.

Theoretical model of damping rings
The equation of motion of two damping rings in the y-direction is By using the separation of variables and letting y n = j y 0 n u n , equation (18) can be written as substituting the generalized coordinate j y 0 n = B yn (T 1 )e Ài1 yn e iv yn T 0 + cc (n = 1, 2) into equation (18), and then the damping rings' equation becomes: g s m Dr iv yn B yn e Ài1 yn e iv yn T 0 À cc À Á u n + f s m Dr B yn e Ài1 yn e iv yn T 0 + cc Where cc represents complex conjugate. Hence, the solution of equation (20) is H yn e Ài1 yn e iv yn T 0 + cc: ð21Þ Substituting equation (21) into equation (20), the H yn and H yn can be expressed as follows: Adding equations (21) into equation (12), one can get the beam equation with damping rings afterwards. The method to obtain the solution in the z-direction is the same with the y-direction. Furthermore, the solution for the double-section beam in time domain is conformity with the isotropic beam. The analyzed process will not repeated here.

Method of multiple scales and frequency response
This study uses the method of multiple scales to analyze the steady state frequency response (fixed points) of this coupled nonlinear system. The time scale is divided into fast and slow time scales. Supposing that T 0 = t is the fast-time scale, T 1 = e 2 t is the slow-time scale, and the approximation to the solutions of y(x, t) and z(x, t) can be expressed as where e is the time scale of small disturbances and is a minimum value. The influence of high-order terms such as e 5 , e 7 . in the system are neglected. Under the assumptions of small perturbations and weak nonlinear vibrations, the damping coefficient c is scaled as e 2 c, and the forcing term q is scaled as e 3 q, respectively. To facilitate the analysis, only the first two terms are selected of the unsteady aerodynamic force for the external force applied on the beam. The orders ofâ D0U andâ D1U are chosen as e 3 and e 2 , respectively. Similarly, the orders ofâ A0U andâ A1U are set as e 3 and e 2 . Substituting these principles into equations (12) and (13), one can obtain equations of motion for different orders.
The frequency response of the system, for example the case of isotropic beam, can be obtained by substituting equation (21) into equations (12) and (13) of order e 3 terms. Taking into account the frequency responses of the system, and assuming that the forcing frequency is where s is tuned frequency near the beam's linear natural frequency v n . The secular terms must be selected to obtain the solvability conditions. Regarding the y-direction's first mode, for example, the harmonic secular terms include the terms of v 1 . Regarding the z-direction's first mode, the harmonic secular terms also include the harmonic terms of v 1 . The secular terms are designated to zero to derive the solvability conditions. Supposing that the y-direction's first mode is excited (q y1 e iOt = q y1 e isT 1 e iv 1 T 0 ), the secular terms of the y-direction's first mode is pre-multiplied by e i1 y1 . Again, the e i1 z1 is premultiplied to the secular terms of the z-direction's first mode. The G A = 1 y1 À 1 z1 , G B = sT 1 + 1 y1 are defined in this study. The following relationships can also be obtained in steady state: In addition, the aerodynamic functions are set as U y, z = constant (i.e.â D0U = const.,â D1U = const.,â A0U = const.,â A1U = const.). Then the real and imaginary parts of the solvability conditions for the yand z-direction equations are: Next, take the square of equations (27) and (28) and sum them together to eliminate the terms associated with time in G B : Once the wind speedq y1 is chosen, the solvability conditions (equations (29)-(31)) can be solved for the amplitudes (B y1 and B z1 ) by using numerical methods. The frequency response of the 2nd mode can be solved similarly. The various parameters of damping rings can be comprehensively analyzed. The same method is used to obtain the frequency response in double-section beam system and is not detailed here. Figure 3(a)-(d) are the fixed points plots of the double-section beam withq y1 = 10, diameter ratio = 1/0.75, and no DRs are installed. Figure 3 (b) demonstrates that in the case of h = 0:33 and when the 2nd mode was excited, the 1st mode's amplitude (the unexcited mode amplitude B y1 ) is larger than the 2nd mode (B y2 ). The I.R. is triggered in this case.  Figures 3 and 4, the DRs not only reduce the vibration, prevent the I.R., but also avoid the primary resonance (B z s are always smaller than the B y s). Therefore, this study will only consider the y-direction's amplitudes on the following Sections.

Floquet theory and stability analysis
This section discusses the effects of damping rings on system stability. In equations (8), (9), (12), and (13), the magnitude of wind speed is changed and the Floquet theory is used to analyze the stability of the system. In order to analyze the stability of the system, the perturbation technique is employed. Assuming that where j yn represents the periodic solution, andj yn is the disturbance term. By substituting these two expressions into equations (8), (9), (12), and (13), applying orthogonal properties, the dynamic equation for the 1st and the 2nd mode in y-direction can be obtained as: ½( j ya +j ya )( j yb +j yb )( j yc +j yc ) + ( j ya +j ya )( j zb +j zb )( j zc +j zc ) The dynamic equation of the damping ring is The Floquet transition matrix can be obtained by setting the state variables for one period. The differential equation with perodic coefficient (e.g. equations (33) and (34)) is expressed in state variable form as where x =j y1 , _ j y1 ,j y2 , _ j y2 , Á Á Á n o T . It is noted that x 1 (0) = 1 0 0 Á Á Á 0 f g T is set as the first initial condition. By using the fourth-order Runge-Kutta method to solve equation (36), x 1 can be obtained. Again, let x 2 (0) = 0 1 0 Á Á Á 0 f g T and put into equation (36) as the second initial condition, and by using the fourth-order Runge-Kutta method to get the results. Following this process, one can get D is the Floquet transition matrix. The eigenvalues of the Floquet transition matrix are shown to be exponentials occurring in Floquet's form of the initial-value solution. Let y i = a + ib as the eigenvalues of D, and the relationship between system eigenvalues (L) and y i are where T is the period of the system. The stability criteria are defined by the Floquet multipliers (F.M.) and F:M: = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi a 2 + b 2 p . When F.M. \ 1, the system is stable; when F.M. . 1, the system is unstable. There are two types for the stable cases, one is asymptotically stable, which means that the system will converge to a certain value; the other is stable limit cycle, which represents that the system will converge to certain periodic region. If some F.M. \ 1 and some F.M. . 1, it can be called the unstable limit cycle of the saddle type, and it is unstable. The stability analysis will be studied to find the diverge wind speed. The damping effects for damping rings will also be investigated. The results are discussed in the next Section.

Damping effect analysis of damping rings
This section will analyze the vibration reduction effects of damping rings for different mass (m Dr ), location (l Dr ), spring constant (f s ), and damping coefficient (g s ). The 3D maximum amplitude plots (3D MAPs) and maximum amplitude contour plots (MACPs) are made to determine the optimal parameter combinations of damping rings.
For the case of isotropic beam, by solving the solvability conditions (equations (29)-(31)), the fixed points plots are made for different parameter combinations: show the isotropic beam's 1st mode maximum amplitudes in y-direction when the 1st mode was excited. Figure 5(a)-(c) are the 3D MAPs for the cases of damping ring's damping coefficient g s = 0.1, 0.5 and 0.9, respectively. This study fixes the 2nd DR near the beam end (l D2 = 0.9999) and changes the 1st DR's positions (l D1 = 0.0001;0.8). In Figure 5, the x-axis is the mass of the damping ring, the y-axis is the position of the damping ring, and the z-axis is the maximum amplitude corresponding to the combinations of the mass and the position of the damping rings. The different amplitude intervals are divided by different colors: red indicates higher amplitude interval, while blue indicates lower amplitude interval, which represents the better vibration reduction effect. These figures show that the minimum amplitude occurs when f s = 9.0, and the damping rings are placed at 0.0001 and 0.9999.
Although the 3D MAPs can provide the general design concept of damping rings, it is difficult to see the optimal combination of damping rings due to numerous damping rings' parameters. Therefore, this study makes the MACPs and analyzes the influence of damping rings' positions and masses on the system. According to the results of the 3D MAPs, f s = 9.0 has the best damping effect. Different g s case of MACPs are made to examine the g s effects. Figure 6(a)-(c) show that when damping rings are located on 0.0001 and 0.9999, the system has the best vibration reduction effect. The reason is that the maximum amplitudes of the first mode of the elastic beam are located on 0 and 1 (both ends of the beam), respectively. The optimal vibration reduction effect can be achieved by placing damping rings on the maximum amplitude of each mode. It is noteworthy that the larger the g s might not have the better effect of vibration reduction. The best damping effect occurs when g s = 0.5. By comprehensive analyzing the 3D MAPs and MACPs, this study concludes that the optimal damping rings vibration reduction combination for the y-direction's first mode is m Dr = 0.04, l D1 = 0.0001, l D2 = 0.9999, f s = 9.0, g s = 0.5.
Next, damping effects of the double-section beam are analyzed. Figsure 7(a) and (b) are the MACPs for the double-section beam with h = 0:33 and the y-direction's 1st mode was excited. The 2nd damping ring is fixed on the 0.9999 of the beam. The position of the 1st damping ring is altered. Figure 7  Same thing happens in the MACPs of the cases of h = 0:51 when the y-direction's first mode was excited (Figure 8(a) and (b)). These results conclude that the optimal damping effect can be reached when m Dr = 0.01, l D1 = 0.0001, l D2 =0.9999, f s = 9.0, g s = 0.5. Table 1 shows the results of the vibration reduction effects with damping rings. In which Amp. (no DRs) is the maximum amplitude of each mode without damping rings, 10 and Amp. (with DRs) is the maximum amplitude of each mode with damping rings. The ''Effect'' is expressed as a percentage of the vibration reduction effects for each mode under different cases. Table 1 shows that vibration reduction effect is remarkable when the DRs are added. In addition, the vibration reduction effect of the DRs on the first mode is greater than that on the second mode.
To further verify the correctness of the MACPs, the fourth-order Runge-Kutta (RK-4) method is used to solve the perturbation equations in time domain (equations (33)-(35)). Figure 9 is an example of the y-direction's 1st mode MACP of the double-section beam with diameter ratio 1/0.75, h = 0:51, when the 1st mode was excited, where the f s = 9.0 and g s = 0.9. Figure 9(a) is the fixed points plot and time response plot for m Dr = 0.1 and l D1 = 0.8. Figure 9

ANSYS simulation
Several cases are studied by the ANSYS to further verify the results of the damping rings' damping effect from MACPs. It is noted that the beam length is assumed to be 1 and the output amplitudes of the ANSYS are assured to be the same as the dimensionless values. The beam model was created by Solidworks commercial software. The material was set as ''structural steel'' and imported into ANSYS. The boundary conditions were set as free on both ends of the beam. More than 20,000 nodes and 12,000 elements were chosen. The simple harmonic function with constant amplitude was set as normal force and applied on the beam surface. Figure 10 shows the ANSYS simulation of the isotropic beam y-direction's 1st mode amplitude, when the 1st mode was excited. Instead of g s = 0.5, we choose the DR's parameters as g s = 0.9, f s = 9.0, m Dr = 0.04, and l D1 = 0.0001 to verify different cases.  Figure 10(b) is the ANSYS simulation with no DRs. The maximum amplitude is 0.0485, which agrees with the result from Table 1 (0.051088). Figure 10(c) is the ANSYS simulation with the DRs. The maximum amplitude is 0.00288, which agrees with the result from MACP. Figure 10(d) is the numerical result with no DRs by using RK-4 method. The maximum amplitude approaches to 0.051, which agrees with the results from Table 1 and ANSYS. Figure 10(e) is the result with DRs by using RK-4 method. The maximum amplitude approaches to 0.0029, which agrees with the results from MACP and ANSYS. Figure 11 shows the ANSYS simulation of the double-section beam y-direction's 1st mode amplitude with diameter ratio 1/0.75, h=0.33, when the 1st mode was excited. The DR's parameters are the same as in Figure 10. Figure 11(a) shows the minimum amplitude is 0.004. Figure 11 (b) is the ANSYS simulation with no DRs. The maximum amplitude is 0.0604, which agrees with the result from Table 1 (0.061192). Figure 11(c) shows the maximum amplitude is 0.0036 with DRs, which agrees with the result from MACP. Figure 11(d) is the numerical result with no DRs by using RK-4 method. The maximum amplitude approaches to 0.061, which agrees with the results from MACP and ANSYS. Figure 11(e) is the result with DRs by using RK-4 method. The maximum amplitude approaches to 0.0036, which agrees with the results from MACP and ANSYS. Again, Figure 12 shows the amplitude simulation of y-direction's 2nd mode, when the 2nd mode was excited. Same DR's parameters are chosen. Figure 12(a) shows the minimum amplitude is 0.0014. Figure 12 (b) shows the maximum amplitude is 0.0191 with no DRs, which agrees with the result from Table 1 (0.020528). Figure 12(c) is the ANSYS simulation with the DRs. The maximum amplitude is 0.00135, which agrees with the result from MACP. Figure 12(d) is the numerical result with no DRs by using RK-4 method. Figure 12(e) is the result with DRs by using RK-4 method. Their maximum amplitudes agree with the results from MACP and ANSYS.

Analysis of system stability
The stability criteria of Floquet theory is used to verify system stability for different wind speeds. By applying Floquet multipliers, if the system is stable, it is marked with a black dot, and if the system is unstable, it is not marked at all. This study analyzes the stability of the system with the value of wind speed varies from Mach  number (M) = 0.54;1.71 (the correspondingq y1 are 1;11), and draws the corresponding 3D basin of attraction, as shown in Figure 13(a). The basin of attraction is a set of initial conditions leading to longtime dynamic motion that approaches to an attractor. Figure 13(a) shows that the stable region of the system decreases in size with the wind speed, and when the wind speed is M = 1.71, the system has hardly stable region. If M increases further, the damping term (c y Àâ D1U ) (or (c z Àâ A1U )) changes from positive to negative and the divergence may occur. To verify the stability prediction on Figure 13(a), the basin of attraction at M = 1.33 (correspondingq y1 = 6) is chosen as an example (Figure 13(b)). The Poincare´map ( Figure  13(c) and (d)) is used to verify the results from Figure  13 (b). Figure 13 Figure  13(c), it can be found that the Poincare´map of the stable region is in a stable condition, while Figure 13(d) shows a diverging system consistent with predictions of basin of attraction. Next, a basin of attraction of the isotropic beam with damping rings at wind speeds M = 0.54;1.71 is drawn to analyze the effect of additional dampers on the stability of the system, as shown in Figure 14(a). Figures 13(a) and 14(a) show that the stability region of the system increases after adding damping rings regardless of the wind speed. By observing Figure 14(a), it can be found that the system has a bigger stable region even when M = 1.71, meaning that the damping ring can reduce the amplitude, increase the stability of the system, and improve the divergence speed of the system. Similarly, to verify the accuracy of the basin of attraction, two damping rings are added to the system. Taking the basin of attraction at M = 1.33 as an example, Figure 14(b) is verified with Poincareḿ aps, as shown in Figure 14(c) and (d). Figure 14(c) and (d) confirm that the dotted block is indeed stable, while the white block is unstable.
In the double-section beam analysis, the case of h = 0.33 is taken as an example, and its 3D basin of attraction is made, as shown in Figure 15(a). Figure  15(a) shows a basin of attraction with wind speed M = 0.48;1.53. Figure 15(a) indicates that the stable region decreases while the wind speed increases. Furthermore, by comparing Figures 13(a) and 15(a), we can find that the stable region of the double-section beam is smaller than that of the isotropic beam. The 3D basin of attraction with damping rings is shown in Figure 16  rings are added. Figure 16(a) shows that, even when M = 1.53, there still exists stable region. The results confirm that the damping ring can increase the range of the stable region and the divergence speed of the system for a certain extent. In addition, taking the basin of attraction at M = 1.12 as an example, Figure 16

Conclusion
In this study, a nonlinear free-free beam is used as the main model to simulate the vibration of a rocket structure. The harmonic wind force is considered. The method of multiple scales and fixed points plots are used to analyze and verify the theory in this study. Two damping rings are equipped on the elastic beam to prevent internal resonance and reduce vibrations. The maximum amplitude contour plots (MACPs) and time response plots are created to analyze the influence of damping rings' parameters (mass, position, spring constant, and damping coefficient) on the beam system and identify the optimal damping effect. The Floquet theory, basin of attraction, and Poincare´maps are also used to determine the stability of the system at different wind speeds. Finally, the following conclusions are proposed for this study.
1. The optimal placement of the damping ring depends on the maximum amplitude of each mode shape. Placing damping rings near the both ends of the beam will result in best damping effects. 2. Two damping rings can effectively prevent internal resonance and reduce vibration. For the isotropic beam, the optimal combination of vibration reduction parameters is m Dr = 0.04, l D1 = 0.0001, l D2 = 0.9999, f s = 9.0, g s = 0.5. For the double-section beam, the optimal parameters combination is m Dr = 0.01, l D1 = 0.0001, l D2 = 0.9999, f s = 9.0, g s = 0.5 and 0.9.     3. The basin of attraction, using Floquet theory, could predict the system divergence in cases with different wind speeds. In addition, the stability regions increase significantly after damping rings are added. The damping rings play an important role in beam vibration reduction.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by the Ministry of Science and Technology of Taiwan, Republic of China (grant number: MOST 109-2224-E-006-004).
noted that for a circular-cross-section isotropic uniform beam, j y =j z and D zz =D yy , the nonlinear coupling terms with the y-and z-directions are neglected. The torsional equation for g is not considered in this work.
Therefore, the 3D flexural-flexural vibration of uniform isotropic beam equation can be obtained as in equation (1).