Microsoft Word - 6-AITI#5290 270-291.docx Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 Study of Primary and Internal Resonance on 3D Free-Free Double-Section Beam Yi-Ren Wang*, Yun-Shuo Chang Department of Aerospace Engineering, Tamkang University, New Taipei City, Taiwan Received 17 February 2020; received in revised form 25 April 2020; accepted 23 June 2020 DOI: https://doi.org/10.46604/aiti.2020.5290 Abstract This work investigates the primary resonance and internal resonance of a double-section beam with cubic nonlinearities. This model can be applied in a wide range of engineering problems, such as rocket and missile structures. Even space technology has been developed for decades; several nonlinear properties deserve further study, especially, for the internal resonance. The method of multiple scales (a perturbation technique) is employed to analyze this nonlinear problem. This study focuses on finding the forcing conditions of this 3D double-section beam to trigger the often-ignored internal resonance or prime resonance in rocket structures. A primary resonance is found on a uniform free-free beam at certain flight speed. The three-to-one internal resonance of the double-section beam occurs within the first and the second modes in the diameter ratio of 1/0.75 with the length ratio of 0.33 or 0.51. The semi-analytical results are verified by the time marching numerical method. Keywords: nonlinear vibration, internal resonance, method of multiple scales 1. Introduction Studies of vibration have always been a concern for researchers and engineers because it may cause the structural fatigue or failure. Beams are widely applied in engineering, such as wings in aerospace, bridges in civil engineering, and train rails in mechanical engineering. Many studies on beam vibrations have been performed previously. Özkaya [1] researched on a beam-mass system under simply supported end conditions, and the effect of positions, magnitudes, and the number of the masses was investigated. Mundrey [2] considered the 2D Bernoulli-Euler beam resting on an elastic foundation to simulate the railway track and a moving load on the 2D beam was studied. Phuoc Nguyen et al. [3] considered the dynamic response of the Euler-Bernoulli beam subjected to moving oscillators. Chang [4] used a variational method to study nonlinear vibrations in carbon nanobeams under the magnetic field. Zhang et al. [5] provided a controllable active torque actuator model for wind-turbine tower vibrations and demonstrated that the proposed module can effectively mitigate the vibrations in wind turbines during operation. These studies demonstrated a wide application of nonlinear models; yet, a 3D nonlinear model and some unique nonlinear properties should be discussed further. In most nonlinear beam problems, an internal resonance is a major point of discussion. The internal resonance is unique to nonlinear systems in which integer relationships exist among the natural frequencies with various modes. Due to nonlinearity, the internal resonance generally occurs in modes that are not being directly excited by external forces. Exciting the higher modes can lead to high amplitude vibrations in lower modes [6]. Large vibration amplitude of the unexcited mode of nonlinear beams could be triggered due to internal resonance. Therefore, the lower modes should not be overlooked. However, it is interesting that in common 3D beams with symmetrical cross sections; one-to-one internal resonance is the most likely to take * Corresponding author. E-mail address: 090730@mail.tku.edu.tw Tel.: +886-2-26215656; Fax: +886-2-26209746 Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 271 place among the various degrees of freedom. As its resonant frequency is the same as each other, it is also called primary resonance. For instance, Pai [7] analyzed the primary resonance in a 3D nonlinear composite rotating beam. Stoykov and Ribeiro [8] 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. Research has also been conducted on the internal resonance of beams caused by external forces or additional structures, such as those associated with the suspension or support. For instance, van Horssen et al. [9-10] considered the use of nonlinear aerodynamics to create three-to-one internal resonance in systems with elastic foundations or suspension springs. Wang et al. [11] analyzed a hinged-hinged nonlinear Bernoulli-Euler beam with linear and nonlinear tuned dampers. The internal resonance was investigated for the system. The nonlinear properties of the nonlinear beam and the nonlinear damper were studied in-depth. Wang and Hsiao [12] studied the 3D nonlinear multi-loaded slender beam. The nonlinear primary resonance was first found in the wind turbine tower. Tekin et al. [13] considered the three-to-one internal resonance in the multiple stepped beam systems. These previous researches investigated the nonlinear properties of internal resonance and primary resonance. However, the internal and primary resonance conditions might happen simultaneously for a 3D beam. They should be taken into consideration at the same time. In addressing general analytical methods to nonlinear vibrations, Nayfeh and Pai [14] investigated vibrations in nonlinear Euler-Bernoulli beams. They formulated a number of useful nonlinear beams equations in accordance with Newton's laws, Euler angles transformation, and the Karman-type strain-displacement relationship. Nayfeh and Mook [6] also proposed a number of perturbation methods by which to solve nonlinear systems, including the Lindstedt-Poincaré method, the method of averaging, and the method of multiple scales (MOMS). Perturbation methods allow the researchers to get good approximations for systems that the exact solutions are not all easy to be solved. Ji and Zu [15] studied the rotating shaft system of a Timoshenko beam, using MOMS to analyze the natural frequency responses of nonlinear systems. Nayfeh and Nayfeh [16] employed MOMS to identify nonlinear modes and nonlinear frequencies, whereupon they applied the Galerkin method to the analysis of dynamic responses in a nonlinear beam. Mao [17] used the ADM method to analyze the vibration of beams consisting of an arbitrary number of steps through a recursive method. He showed that the ADM offers an accurate and effective method of free vibration analysis of multiple-stepped beams with arbitrary boundary conditions. The 3D free-free beam model can be used to simulate space rocket vibration in high-speed motion. The double-section 3D beam model can also be utilized to model a high-speed double-stage rocket motion. The most concerns on the multi-stage rocket motion are the stage separation dynamic control [18] or the high speed flows over the joint of the space rocket stages [19]. According to the references aforementioned with 3D beams or rocket motions, the nonlinear properties of both primary resonance and internal resonance are less discussed simultaneously in the 3D beams or rocket structures. The present study considers a nonlinear double section free-free beam subjected to the distributed load with the wind force and its associated unsteady aerodynamic force. Since the primary resonance or internal resonance is unique in nonlinear problems, the unexcited modes usually have larger amplitudes than the excited modes. However, these conditions cannot be predicted by using linear theories. It is worth having a deep study on this model. The present work focuses on finding the forcing conditions of this 3D double-section beam to trigger the often-ignored internal resonance or prime resonance in rocket structures. This research uses Newton’s second law to derive the nonlinear equation of motion. The method of multiple scales (MOMS) is used to obtain the steady state frequency response (fixed point). A primary resonance is found for a uniform free-free beam at certain flight speed. The three-to-one internal resonance is also triggered within the first and the second modes in certain combinations of beam length and diameter ratios of a double-section beam. The fixed points plot is used to observe the nonlinear phenomenon in the system. The fourth-order Runge-Kutta method is also performed to verify the results of frequency response and the internal resonance. Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 272 2. Theoretical Model 2.1. Equations of motion for uniform free-free beam This study analyzes the vibrations in a straight 3D free-free nonlinear beam, which could be used to simulate the high speed aircrafts such as rocket, missiles, etc. The schematic of the space rocket, the coordinate definitions of the simulated beam, and the relationships between the various external forces are demonstrated in Fig. 1. Fig. 1(a) shows the space rocket model similar to the Satellite and Rocket Propulsion Lab. at NCKU. It is true that the proposed beam model can be applied to simulate a missile or a rocket vibration motion. Even though the rocket model is hollow, the bulkheads, formers, longerons, stringers and struts strengthen the fuselage structure of the rocket. The beam model is used to simulate the rocket motion for a preliminary study. In this research, the internal resonance and the primary resonance of the rocket nonlinear properties are studied. Fig. 1(b) represents the schematic of the coordinate definitions and the multiple forces applied to a uniform free-free beam. The initial status of the beam is considered straight under the assumption that each cross-section is a plane. The deformation of the beam also follows the stress-strain laws. Three consecutive Euler angles are used to relate the deformed and unreformed states. The equations contain structural coupling terms, quadratic, and cubic nonlinearities due to curvature and inertia. Based on Newton’s 2nd law, Euler angle transformation, and Taylor series expansion, the equations of motion of the nonlinear beam can be expressed as follows [14]: y y y mv c v G F′+ = +ɺɺ ɺ (1) z z zmw c w G F′+ = +ɺɺ ɺ (2) x x x xj c G Mγ γ+ = +ɺɺ ɺ (3) (a) Schematic of NCKU HTPB/N2O hybrid rocket model (b) Schematic of multiple forces applied to a uniform free-free beam Fig. 1 The rocket model and a uniform free-free beam approximated model where m is the beam mass per unit length, ��̅ is the mass moment of inertia in the x-direction, ��̅ , ��̅, and ��̅ represent the damping coefficient of the beam in the y-, z-, and x-axis, respectively. The � � , � � , and �� represent the unsteady harmonic wind force and aerodynamic moment in the y-, z-, and x-directions, respectively. The v, w, and � represent the beam displacement (or twist angle) in the y-, z- and x-directions, respectively. ⋅� denotes �/��̅ , and �� denotes �/��̅ . Furthermore, Gy, Gz, and Gx are defined below for a uniform and isotropic beam [14]: 2 0 2 2 0 ( ) ( ) ( ) ( )[( ) ] (1/ 2) [ ( ) ] x y zz xx zz yy zz x x l G D v D w D v v v w w D D w v w v w dx v m v w dx dx γ γ γ′′ ′ ′ ′′ ′ ′ ′′ ′ ′′ ′ ′′ ′′ ′ ′′′ ′′ ′= − − − + + − − − ′ ′ ′− + ∫ ∫ ∫ ɺɺ (4) Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 273 2 0 2 2 0 ( ) ( ) ( ) ( )[( ) ] (1/ 2) [ ( ) ] x z yy xx yy yy zz x x l G D w D v D w v v w w D D v w v w v dx w m v w dx dx γ γ γ′′ ′ ′ ′′ ′ ′ ′′ ′ ′′ ′ ′′ ′′ ′ ′′′ ′′ ′= − + − + + − + + ′ ′ ′− + ∫ ∫ ∫ ɺɺ (5) 2 2 2 2 0 ( ) ( )[( ) ] [( ) ( ) ] ( )[( ) ] x x xx zz yy x y z G D D D v w v w j v w dx v w j j v w w vγ γ γ′ ′ ′′ ′′ ′′ ′′ ′′ ′ ′ ′ ′ ′ ′ ′= + − − − + − + − − −∫ ɺɺ ɺ ɺ ɺ ɺ ɺ ɺ (6) where �̅ is the position coordinate along the beam in the axial direction, � ̅ is the beam length, ��� = ���̅ , ��� = ���̅ , ��� = ���̅, and G is the elastic shear modulus and E is Young’s modulus. The ��̅,�,� represent the area moment of inertia in the x-, y-, and z-directions, respectively. The dimensionless coefficients are defined as follows: 4 2 2 2 2 2 xy zy / 1, / , / , / / , / , / / , / , / , / yy y y yy z z yy x x yy x x x y y z z xx yy zz yy l l l x x l t t D ml c c l mD c c l mD c c mD j j j ml j j ml j j ml D D D Dµ µ  = = = = =   = = =  = = = = (7) To simplify the notation, the symbols for the dimensionless displacement of the beam were used as the same as the three axes. In other words, y and z are respectively used to represent the transverse (transverse dir., � = �/�)̅ and lateral (side dir., � = �/�)̅ dimensionless displacement functions to avoid influencing the derivation of the theoretical model and to maintain consistency in the overall equations of motion. Therefore, the 3D flexural-flexural-torsional vibration of uniform isotropic beam equation can be written as: 2 0 2 2 1 0 ( ) [ ( ) ] (1 )[( ) ] 1 { [ ( ) ] } 2 x iv zy y xy zy zy x x y y y c y z y y y z z y z z y z dx y y z dx dx F µ µ γ µ µ γ γ′ ′′ ′ ′ ′ ′′ ′ ′′ ′ ′ ′′ ′′ ′ ′′′ ′′ ′ ′+ = − − − + − − − + ′ ′ ′ ′− + + ∫ ∫ ∫ ɺɺ ɺ ɺɺ (8) 2 0 2 2 1 0 ( ) [ ( ) ] (1 )[( ) ] 1 { [ ( ) ] } 2 x iv z xy zy x x z z z c z y z y y z z y z y y z dx z y z dx dx F µ γ µ γ γ′ ′′ ′ ′ ′ ′′ ′ ′′ ′ ′ ′′ ′′ ′ ′′′ ′ ′′ ′+ = − + − + + − + + ′ ′ ′ ′− + + ∫ ∫ ∫ ɺɺ ɺ ɺɺ (9) 2 2 2 2 0 1 ( ) ( ) ( ) ( ) xxy zy y z x x x x x j j c y z y z y z dx y z y z y z M j j j µ µ γ γ γ γ γ γ γ − − ′′ ′′ ′′ ′′ ′′ ′′ ′ ′ ′ ′ ′ ′ ′− = − + − − + − + − − +∫ɺɺ ɺ ɺɺ ɺ ɺ ɺ ɺɺ ɺ (10) In addition, the cross-section of the beam in this study is a circle, and the ratio of the lateral area moments of inertia is equal to 1, which means that �� = 1. The structural damping coefficients are also the same, which means that cy= cz. It is noted that the structural damping is always used to simulate a solid material when it flexes. As for the simulation of beam vibrations, the damping coefficient is rather small, and cy is taken as 0.1 in this study. Assuming the beam is subjected to a distributed load with the harmonic wind force " �,�# $%�&,'(̅. Its associate unsteady aerodynamic force � ),*, that is, � �,� can be expressed as � �,� = " �,�# $%�&,'(̅ + � ),*. The wind force can be normalized as: , , 2 4 , , / / y z y z yy y z y z yy q q D l l D ml = Ω = Ω (11) The windward aerodynamic forces in the y-direction of the beam are given as follows (please refer the research from van Horssen [20]), Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 274 2 2 3 1 2 3 0 2 3 2 2 3 1 2 3 0 2 3 ( , ) ( , ) ( , ) ... 2 ( , ) ( , ) ( , ) ... 2 a y D U D U D U D D U y y y a z A U A U A U A A U z z z U d a a ay x t y x t y x t F a U t t tU U U d a a az x t z x t z x t F a U t t tU U ρ ρ   ∂ ∂ ∂ = + + + +   ∂ ∂ ∂     ∂ ∂ ∂ = + + + +  ∂ ∂ ∂  (12) Normalizing Eq. (12) yields: 2 3 0 1 2 3 2 3 0 1 2 3 ˆ ˆ ˆ ˆ ..., ˆ ˆ ˆ ˆ ..., D D U D U D U D U A A U A U A U A U F a a y a y a y F a a z a z a z  = + + + +  = + + + + ɺ ɺ ɺ ɺ ɺ ɺ (13) where ,-)./ is the lifting coefficients and ,-)0/ is the first-order derivative modified term. Both of them are defined as: 2 0 0 2 ˆ 2 a y D U D U U da a m ρ ω = (14) 1 1 2 ˆ 2 a y D U D U U da a m ρ ω = (15) 2 2 2 ˆ 2 a D U D U da a m ρ ω = (16) 2 3 2 ˆ 2 a D U D U y da a m U ρ ω = (17) 2 0 0 2 ˆ 2 a z A U A U U da a m ρ ω = (18) 1 1 2 ˆ 2 a z A U A U U da a m ρ ω = (19) 2 2 2 ˆ 2 a A U A U da a m ρ ω = (20) 3 3 2 ˆ 2 a A U A U z da a m U ρ ω = (21) The aerodynamic coefficients of the y and z directions, and 1��,� denotes the wind speed in the y and z directions. If 1��,� is also in unsteady form, then 2 0 0ˆ i t D U D Ua a e Ω= ɶ (22) 1 1ˆ i t D U D Ua a e Ω= ɶ (23) 2 0 0ˆ i t A U A Ua a e Ω= ɶ (24) 1 1ˆ i t A U A Ua a e Ω= ɶ (25) It is also noted that ,-)./ is zero for non-rotating body, and ,-)0/ is 22. This study investigates the internal resonance of a 3D nonlinear free-free beam. The main goal is focusing on the structural vibration analysis due to unsteady loads. The fluid-solid interaction aeroelastic problem is not considered in this Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 275 research. The boundary conditions for the free-free beam are simply the structural dynamic model. The dimensionless boundary conditions of the beam are: (0, ) (0, ) (1, ) (1, ) 0 (0, ) (0, ) (1, ) (1, ) 0 (0, ) (1, ) 0 y t y t y t y t z t z t z t z t t tγ γ ′′ ′′′ ′′ ′′′= = = =  ′′ ′′′ ′′ ′′′= = = =  ′ ′= = (26) 2.2. Double section beam model The equations of motion of a double section free-free beam (Fig. 2) can be re-written by using Eqs. (1)-(6) as follows: Fig. 2 Schematic model of double section uniform free-free beam 1 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 ( ) ( ) ( ) ( ) ( ) [ ( ) ] ( )[( ) ] 1 { [ ( ) ] 2 i i i i i i i i i i i y i yy i xx i i zz i i i i i x yy zz i i i i i i i i l x i i i i i l m v c v D v D w D v v v w w D D w v w v w dx v m v w dx γ γ γ − − ′′ ′′ ′ ′′ ′ ′ ′ ′′ ′ ′′ ′ ′+ = − + − + ′′ ′′ ′ ′′′ ′′ ′ ′+ − − − ′ ′ ′− + ∫ ∫ ɺɺ ɺ ɺɺ } i i x i y l dx F′ +∫ (27) 1 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 ( ) ( ) ( ) ( ) ( ) [ ( ) ] ( )[( ) ] 1 { [ ( ) ] 2 i i i i i i i i i i i z i yy i xx i i yy i i i i i x yy zz i i i i i i i i l x i i i i i l m w c w D w D v D w v v w w D D v w v w v dx w m v w dx γ γ γ − − ′′ ′′ ′ ′′ ′ ′ ′ ′′ ′ ′′ ′ ′+ = − + − + ′′ ′′ ′ ′′′ ′′ ′ ′+ − + + ′ ′ ′− + ∫ ∫ ɺɺ ɺ ɺɺ } i i x i z l dx F′ +∫ (28) 1 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( )[( ) ] [( ) ( ) ] ( )[( ) ] i i i i i i i i i x i x i xx i zz yy i i i i i x x i i i i i l y z i i i i i x j c D D D v w v w j v w dx v w j j v w w v M γ γ γ γ γ − ′ ′ ′′ ′′ ′′ ′′+ = + − − − ′′ ′ ′ ′+ − ′ ′ ′ ′+ − − − + ∫ ɺɺ ɺ ɺɺ ɺ ɺ ɺ ɺ ɺ ɺ (29) The index i =1, 2 and represents the section number for this beam. For a circular cross-section beam, the moment of inertia 3�4 = 3�4 , and the side moment of inertia ratio equals 1 ( ��4 = 1), the structural dampings in the y- and z-dir are the same, cy= cz, the two sections of the beam are using the same material E1= E2. The wind force applied on the double-section beam is uniform and expressed as " �,�# $%�&,'(̅ and is the same as in the case of uniform beam (Fig.1). This beam is subjected to unsteady aerodynamic forces in the y and z-dir. These forces are expressed as � ),*, where , , 1 2 4 , , 1 / / y z y z yy y z y z yy q q D l l D m l = Ω = Ω (30) 2 2 3 1 ( ) ( ) ( )1 2 3 0 2 3 2 32 ( ) ( ) ( )1 1 2 3 0 2 3 ( , ) ( , ) ( , ) ... 2 ( , ) ( , ) ( , ) ... 2 a y i i iD U D U D U D D U y y y i i ia z A U A U A U A A U z z z U d y x t y x t y x ta a a F a U t t tU U z x t z x t z x tU d a a a F a U t t tU U ρ ρ   ∂ ∂ ∂ = + + + +   ∂ ∂ ∂     ∂ ∂ ∂ = + + + +   ∂ ∂ ∂  (31) Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 276 The dimensionless forms are: 2 3 0 1 2 3 2 3 0 1 2 3 ˆ ˆ ˆ ˆ ..., ˆ ˆ ˆ ˆ ..., D D U D U D U D U A A U A U A U A U F a a y a y a y F a a z a z a z  = + + + +  = + + + + ɺ ɺ ɺ ɺ ɺ ɺ (32) where 2 1 0 0 2 1 ˆ 2 a y D U D U U d a a m ρ ω = (33) 1 1 1 2 1 ˆ 2 a y D U D U U d a a m ρ ω = (34) 1 2 2 2 1 ˆ 2 a D U D U d a a m ρ ω = (35) 1 3 3 2 1 ˆ 2 a D U D U y d a a m U ρ ω = (36) 2 1 0 0 2 1 ˆ 2 a z A U A U U d a a m ρ ω = (37) 1 1 1 2 1 ˆ 2 a z A U A U U d a a m ρ ω = (38) 1 2 2 2 1 ˆ 2 a A U A U d a a m ρ ω = (39) 1 3 3 2 1 ˆ 2 a A U A U z d a a m U ρ ω = (40) It is noted that Eqs. (30)-(32) are used for the two-stage beam. Instead of using “m” and “d ” in Eqs. (11)-(21), Eqs. (30)-(32) use the dimensions of the first stage beam’s mass “m1” and the characteristic length “�̅0” to express the aerodynamic forces. In order to make dimensionless beam equations, the following definitions are introduced. 1 1 1 1 4 1 2 2 1 1 1 1 2 2 2 1 1 1 / , / , / , / , / / , / , / / , / , / / , / i i i i i i i i i i i i i i i i i i i i yy y y yy z z yy x x yy x x x y y z z xy xx yy zy zz yy x x l l l l y v l z w l t t D m l c c l m D c c l m D c c m D j j j m l j j m l j j m l D D D Dµ µ  = = = = =   = = =   = = =  = = (41) The dimensionless equations of motion of this double section beam can be obtained as: ( ) 1 ( ) ( ) 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 ( ) ( ) ( ) 0 ˆ( ) ( ) [ ( ) ] 1 ˆ{ [ ( ) ] , 1, 2} 2 i i i y i i iv i i y D U i xy i i i i i i i x x i t i i i i i y D U l l y y c a y z y y y z z y y z dx dx q e a i µ γ − Ω ′ ′′ ′ ′ ′ ′′ ′ ′′ ′ ′+ = − − − − + ′ ′ ′ ′ =− + + +∫ ∫ ɺɺ ɺ ɺɺ 　 (42) ( ) 1 ( ) ( ) 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 ( ) ( ) ( ) 0 ˆ( ) ( ) [ ( ) ] 1 ˆ{ [ ( ) ] , 1, 2} 2 i i i z i i iv i i z A U i xy i i i i i i i x x i t i i i i i z A U l l z z c a z y z y y z z z y z dx dx q e a i µ γ − Ω ′ ′′ ′ ′ ′ ′′ ′ ′′ ′ ′+ = − − + − + ′ ′ ′ ′ =− + + +∫ ∫ ɺɺ ɺ ɺɺ 　 (43) Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 277 ( ) 1 ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) , ( 1 2) + , i i i i xxy i i x i i i i i x l x c y z dx y z M j i µ γ γ γ − ′′ ′′ ′ ′ ′− = − + − =∫ɺɺ ɺ ɺɺ ɺ ɺ 　 (44) The boundary conditions (Eq. (45)) and the compatibility equations (Eq. (46)) at the joint of the two beam sections are: (1) (1) (2) (2) (1) (1) (2) (2) (1) (2) (0, ) (0, ) (1, ) (1, ) 0 (0, ) (0, ) (1, ) (1, ) 0 (0, ) (1, ) 0 y t y t y t y t z t z t z t z t t tγ γ ′′ ′′′ ′′ ′′′ = = = =  ′′ ′′′ ′′ ′′′= = = =  ′ ′= = (45) (1) (2) (1) (2) (1) (2) (1) (2) (1) (2) ( , ) (0, ), ( , ) (0, ) ( , ) (0, ), ( , ) (0, ) ( , ) (0, ) y t y t y t y t z t z t z t z t t t η η η η γ η γ ′ ′ = =  ′ ′= =  = (46) where 5 = �0/ �0 + �6�. 2.3. Preturbation analysis Perturbation methods allow the researchers to get good approximations for systems where the exact solutions are not all easy to be solved. This study adopted the method of multiple scales (MOMS) to analyze the frequency response and make the fixed point plots of this nonlinear system. The time scale was divided into fast and slow time scales. Suppose that 7. = � is the fast-time term, 70 = 86� is the slow-time terms, and the expansions of each direction are: 1 3 ( ) ( )0 0 1 ( )1 0 1 1 3 ( ) ( )0 0 1 ( )1 0 1 1 3 ( ) ( )0 0 1 ( )1 0 1 ( , ; ) ( , , ) ( , , ) ..., ( , ; ) ( , , ) ( , , ) ..., ( , ; ) ( , , ) ( , , ) ..., i i i i i i i i i y x t y x T T y x T T z x t z x T T z x T T x t x T T x T T ε ε ε ε ε ε γ ε ε γ ε γ  = + +  = + +  = + + (47) where 8 is the time scale of small disturbances and is a minimum value. For the sake of simplicity, the influence of high-order terms such as 89, 8:... are neglected in the system. Two time scales (T0 and T1) are considered in this study. Under the assumption of nonlinear vibrations, the damping coefficient c is scaled as 86�, the forcing term q is scaled as 8;", respectively. To facilitate this analysis, only the first two terms of the unsteady aerodynamic force are extracted for the external force. The orders of ,-)./ and ,-)0/ in the �) equation are chosen as 8; and 86. Similarly, the orders of ,-*./ and ,-*0/ in the �* equation are set as 8; and 86. These principles are substituted into Eqs. (42)-(44) and obtain the expansion of the equation of the 80 order in the y-direction. ( )2 0 ( )0 ( )0 0, 1, 2iv i i D y y i =+ = 　 (48) The equation of order 8; is: ( ) 1 2 0 ( )1 ( )1 0 1 ( )0 1 0 ( )0 ( )0 ( )0 ( )0 ( )0 ( )0 2 2 ( )0 ( )0 ( )0 0 ˆ2 ( ) [ ( ) ] 1 ˆ{ [ ( ) ] , } 2 1, 2 i i y i i iv i i i y D U i i i i i i x x i t i i i i i y D U l l D y y D D y c a D y y y y z z y y z dx dx q e a i − Ω ′ ′ ′′ ′ ′′ ′ ′+ = − − − − + ′ ′ =′ ′− + + +∫ ∫ ɺɺ 　 (49) The equation of the 80 order in the z direction is: ( )2 0 ( )0 ( )0 0, 1, 2iv i i D z z i =+ = 　 (50) The equation of order 8; is: Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 278 ( ) 1 2 0 ( )1 ( )1 0 1 ( )0 1 0 ( )0 ( )0 ( )0 ( )0 ( )0 ( )0 2 2 ( )0 ( )0 ( )0 0 ˆ2 ( ) [ ( ) ] 1 ˆ{ [ ( ) ] , } 2 1, 2 i i z i i iv i i i z A U i i i i i i x x i t i i i i i z A U l l D z z D D z c a D z z y y z z z y z dx dx q e a i − Ω ′ ′ ′′ ′ ′′ ′ ′+ = − − − − + ′ ′ =′ ′− + + +∫ ∫ ɺɺ 　 (51) The equation of the 80 order in the γ direction is: ( )2 0 ( )0 ( )0 , 1, 20i i xy i i x iD j µ γ γ ′′ =− = 　 (52) The equation of order 8; is: ( )2 0 ( )1 ( )1 0 1 ( )0 0 ( )0 , 1 22 ,i i xy i i i x i x x D D D c D M i j µ γ γ γ γ =′′− = − − + 　 (53) The boundary conditions and the compatibility equations at the joint of the two beam sections are = 80�: (1)0 (1)0 (2)0 (2)0 (1)0 (1)0 (2)0 (2)0 (1)0 (2)0 (0, ) (0, ) (1, ) (1, ) 0, (0, ) (0, ) (1, ) (1, ) 0, (0, ) (1, ) 0. y t y t y t y t z t z t z t z t t tγ γ ′′ ′′′ ′′ ′′′ = = = =  ′′ ′′′ ′′ ′′′= = = =  ′ ′= = (54) (1)0 (2)0 (1)0 (2)0 (1)0 (2)0 (1)0 (2)0 (1)0 (2)0 ( , ) (0, ), ( , ) (0, ), ( , ) (0, ), ( , ) (0, ), ( , ) (0, ). y t y t y t y t z t z t z t z t t t η η η η γ η γ ′ ′ = =  ′ ′= =  = (55) The equation of order 8; is: (1)1 (1)1 (2)1 (2)1 (1)1 (1)1 (2)1 (2)1 (1)0 (2)0 (0, ) (0, ) (1, ) (1, ) 0, (0, ) (0, ) (1, ) (1, ) 0, (0, ) (1, ) 0. y t y t y t y t z t z t z t z t t tγ γ ′′ ′′′ ′′ ′′′ = = = =  ′′ ′′′ ′′ ′′′= = = =  ′ ′= = (56) (1)1 (2)1 (1)1 (2)1 (1)1 (2)1 (1)1 (2)1 (1)1 (2)1 ( , ) (0, ), ( , ) (0, ), ( , ) (0, ), ( , ) (0, ), ( , ) (0, ). y t y t y t y t z t z t z t z t t t η η η η γ η γ ′ ′ = =  ′ ′= =  = (57) 2.4. 3D free-free beam free vibration analysis The purpose of this section is to find the mode shapes of this vibration beam by the free vibration analysis. By using the MOMS (Section 2.3), the 80 order in the y direction is expressed in Eq. (48). The � $�. is divided into time and space by using the separation of variables, defined as � $�. = > ��7 ��, and substituted in Eq. (48) to obtain: 2 40 ( )( ) ( ) ( ) iv D T tX x X x T t α= − = (58) where ? is the eigenvalue of the system. The general solution of > �� is assumed as: 1 2 3 4( ) (sin sinh ) (sin sinh ) (cos cosh ) (cos cosh )X x E x x E x x E x x E x xα α α α α α α α= + + − + + + − (59) Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 279 and substituted into the boundary conditions, which reveals that E1 = E3 = 0 and that the characteristic equation is cos cosh 1α α = (60) then the mode shape of the i th mode can be expressed as the following: )cosh(cos coscosh sinhsin )(sinhsin)()( xxxxx ii ii ii iii αα αα αα ααϕ + − − ++= (61) 3. Frequency Response 3.1. Conditions for internal resonance The natural frequency will be changed if the diameter and length of the beam are altered. The diameter ratios and the sectional length ratios of the compound beam are the main properties to trigger internal resonance. References [21 and 22] showed that the diameter ratio range is from 1/0.6 to 1/0.9 for the space rockets. The case of �̅0/�̅6 = 1/0.75 is chosen as the test case. The Transfer Matrix Method [17, 23] is used to find the relationship between modal frequency and the sectional length ratio (5) for the beam. The conditions for internal resonance can be obtained by a certain combination of 5. Taking an example of the y-direction motion and using the separation of variables, the following expression for y(i)0 can be obtained by: ( ) 0( )0 ( ) ( ) , 1,) 2(i i yy x itϕ ξ == 　 (62) Substituting Eq. (62) into Eq. (48) and dividing by D $� ��E�. ��, it yields: ( )0 0 2 0( ) 4 ( ) ( )( ) , ) 1, ( ( 2 ) iv yi i i y tx t i D x ξϕ λ ϕ ξ = == − 　 (63) where D $� �� can be expressed as: (1) 1 1 1 1 1 1 1( ) cosh( ) sinh( ) cos( ) sin( ) i x A x B x C x D xϕ λ λ λ λ= + + + (64) (2) 2 2 2 2 2 2 2 2( ) cosh( ( )) sinh( ( )) cos( ( )) cos( ( ))x A x B x C x D xϕ η λ η λ η λ η λ η− = − + − + − + − (65) where F0 is the eigen value. From [23], F0 and F6 satisfy the relationship F6 = G6�0̅/G0�6̅� 0/HF0, where A is the area of beam’s cross-section and I is the moment of inertia. Since circular cross-sections are considered, F6 = �̅0/�̅6� 0/6F0 = 1.1547F0. The mode shape is expressed as: ( )(1) (2) (1)( ) ( ) ( ) ( )x x x x H xϕ ϕ ϕ η ϕ η = + − − −  (66) where J � is the Heaviside function. The researchers define matrix K 0� and L 0� as 1 1(1) 2 2 1 1 3 3 1 1 1 0 1 0 0 1.1547 0 1.1547 0.4219 0 0.4219 0 0 0.4871 0 0.4871 λ λ λ λ λ λ      =  −   −  P (67) Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 280 1 1 1 1 1 1 1 1 1 1 1 1(1) 2 2 2 2 1 1 1 1 1 1 1 1 3 3 3 3 1 1 1 1 1 1 1 1 cosh( ) sinh( ) cos( ) sin( ) sinh( ) cosh( ) sin( ) cos( ) cosh( ) sinh( ) cos( ) sin( ) sinh( ) cosh( ) sin( ) cos( ) ηλ ηλ ηλ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ λ ηλ    − =  − −   −  Q (68) K 0� represents the moment of inertia of the 2nd section to the 1st section of the beam. L 0� represents the displacement, slope, bending moment, and shear force shape function. The Transfer Matrix is expressed as: 1(1) (1) (1)−  =  T P Q (69) and also satisfies the following condition, { } { }(1) 2 2 2 2 1 1 1 1 T T A B C D A B C D= T (70) From the boundary condition D 0� �� 0� = D 0� ��� 0�, it can be obtained by: 1 1 1 1 1 0 1 0 0 0 1 0 1 0 A B C D    −     =    −        (71) The external boundary conditions are: cosh sinh cos sin sinh cosh sin cosl S S S S S S S S − −  =  −  B (72) where M = 1.1547 1 − 5�F0 and { }2 2 2 2 0 0 T l A B C D   =     B (73) Substituting Eq. (70) into Eq. (73), it yields: { } { }(1) 1 1 1 1 1 1 1 1 0 0 T T l A B C D A B C D   = =     H B T (74) where O = PQ ̅R 0� 11 12 13 14 21 22 23 24 H H H H H H H H   =     H (75) and is the total transfer matrix. Eqs. (71) and (74) can be solved simultaneously, 1 1 11 12 13 14 1 21 22 23 24 1 1 0 1 0 0 0 1 0 1 0 0 0 A B H H H H C H H H H D −           −      =                 (76) The characteristic equation for the eigen values is: Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 281 11 12 13 14 21 22 23 24 1 0 1 0 0 1 0 1 det( ) 0 H H H H H H H H −   −  =       (77) let 00 1 ( , ) ( ) ( )y n n n y x t t xξ ϕ ∞ = =∑ (78) 11 1 ( , ) ( ) ( )y n n n y x t t xξ ϕ ∞ = =∑ (79) Substituting into Eqs. (78) and (79) into Eqs. (48) and (49) and using orthogonal properties, the following dynamic equations can be obtained, = 80�: 0 0 2 4 0 1 0y m m y mD ξ λ ξ+ = (80) = 8;�: 1 1 0 0 0 0 0 0 0 0 2 4 0 1 0 1 1 0 1 1 1 1 0 0 0 0 , , 1 1 0 ˆ2 ( ) ( 3 ) ( 3 y m m y m y m y D U y m iv m y i y j y k m i j k m i j k m i j k m i j k i j k m y i z j z k m i j k m i j D D D c a D G dx dx dx dx G dx ξ λ ξ ξ ξ ξ ξ ξ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ξ ξ ξ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ∞ = + = − − − ′′ ′′ ′′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′− + + + ′′ ′′ ′′ ′ ′′ ′′′− + ∑ ∫ ∫ ∫ ∫ ∫ ɶ ɶ 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 , , 1 , , 1 ) 1 ( 2 2 ) 2 ( i i i i i i i iv k m i j k m i j k i j k m y i y j y k y j y k y j y k z j z k z j z k z j z k i j k x x x x m i j k i i i j k i i l l l l dx dx dx G dx dx dx dx ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ϕ ϕ ϕ ϕ ϕ ϕ ϕ − − ∞ = ∞ = ′′ ′ ′′′ ′ ′+ + − + + + + + ′′ ′ ′ ′ ′ ′⋅ + ∑ ∫ ∫ ∫ ∑ ∫ ∫ ∫ ɶ ɺɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ ɺɺ 21 1 00 0 1 ˆ) i yi t ym m D U m i dx q e G a dxϕ Ω =   ′ + +     ∑∫ ∫ ∫ɶɶ (81) The frequency as ST = F0T 6 , m=1, 2, 3…can be obtained from Eq. (80). The relationship between different frequency ratios and the beam sectional length ratios (5) is shown in Fig. 3. The integer multiple of frequency ratios has the potential to trigger the internal resonance (I.R.). It is noted that higher orders of structural modes do not exist in linear cases. The energy transferring between lower and higher modes does not happen in the linear structural models. Since there are so many combinations of frequency ratios in this nonlinear system, the fixed point’s plots are used to examine if the internal resonance really happens in these cases. Besides, for the cases with many combinations of the ratios between higher and lower modes frequencies, the energy may not be possible to transfer between these modes. The possibility of higher modes frequency ratios to trigger I.R. is much less than the cases of small frequency ratios. This study also reveals that when 5 = 0.33 & 0.51, S0: S6 = 1: 3, 5 = 0.68 & 0.86, S0: S6 = 1: 5 and 5 = 0.25 & 0.42, S0: S6 = 1: 6 are possible to trigger I.R. According to Eq. (53), only cubic order of E exists in this system, which means only 1:3 I.R. exists in the 1st and 2nd modes. Therefore, only the case of 1:3 I.R. in the 1st and 2nd modes will be investigated in this study. 5 = 0.33 is chosen and substituted into Eqs. (68), (72). By using Eqs. (67)-(76), this study can obtain transfer matrix R 0� and find the first three eigenvalues: 4.03, 6.98, 10.04. The 1st section beam mode shape coefficients can be obtained by Eq. (76). The 2nd section beam mode shape coefficients can be obtained by Eq. (70). The mode shape functions of the 1st and 2nd beam sections can be obtained by Eqs. (65)-(66). Similarly, this study can gain other mode shapes for different 5. Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 282 Fig. 3 Frequency ratios for different beam sectional length ratios 3.2. Frequency analysis Since 1:3 I.R. happens in the 1st and 2nd modes, the 1st and 2nd modes will be discussed in the following analysis. The solutions to the generalized coordinate are introduced: 0 0 1( ) m mi i T m mB T e e cc ς ωξ −= + (82) 0 0 00 1( ) m mi i T m m m mD i B T e e cc ς ωξ ξ ω −= = +ɺ (83) 0 0 00 1 1 1( ( ) ) ( ( ) )m m m mi i T i i T m m m m m mD D i B T e e cc B T e e cc ς ω ς ωξ ω ω ς− −′ ′= + + + (84) where �� ≡ �0 � and cc represents complex conjugate. In Eq. (81), the distributed load is assumed as "T# $%\ = "T# $ ]^_`ab�cd = "T# $`abcd#$]^cd = "T# $bce#$]^cd . If the solving procedure continues, terms containing the factors of the system frequencies (or harmonics) appear on the right-hand of Eq. (81). Terms such as these are called secular terms. Because of the secular terms, the solution of Eq. (81) increases without bound as t increases. The time scale of E0 does not provide a small correction to time scale of E.. The secular terms for S0, S6 − 2S0 harmonics of the y-dir. 1st mode (m = 1) can be selected. Likewise, the secular terms due to harmonic numbers of S6, 3S0for the y-dir 2nd mode, due to harmonic numbers of S0, S6 − 2S0 for the z-dir 1st mode, and due to harmonic numbers of S6, 3S0 for the z-dir 2nd mode can be obtained, respectively. Refer to Appendix for the secular terms. Next is to let the secular terms as 0 to get the solvability conditions for solving the mode amplitudes in the frequency domain. The fixed point’s plots can be obtained to analyze the frequency response for the 1st and 2nd modes of the y- and z-dir. Fig. 4 Fixed point plot of diameter ratio = 1/0.75, 5 = 0.33, exciting the y-dir., 1st mode Fig. 5 Fixed point plot of diameter ratio = 1/0.75, 5 = 0.33, exciting the y-dir., 2nd mode Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 283 Figs. 4-5 are the fixed point plot of the y-dir., 5 = 0.33 and for the 1st and 2nd mode, individually. Even though the 2nd mode is excited, the 1st mode amplitude is larger than the 2nd mode (the excited mode), and thus an I.R. happens in the 1st and 2nd modes. The same thing happens in the case of 5 = 0.51. The fixed point’s plots are shown in Figs. 6-7 for the 1st and 2nd mode are excited. Fig. 6 Fixed point plot of diameter ratio = 1/0.75, 5 = 0.51, exciting the y-dir., 1st mode Fig. 7 Fixed point plot of diameter ratio = 1/0.75, 5 = 0.51, exciting the y-dir., 2nd mode To verify the frequency domain results, the fourth order Runge-Kutta method (RK-4) is used to get the numerical results. The fourth order Runge–Kutta method is a family of implicit and explicit iterative methods in numerical analysis. It is used in temporal discretization for the approximate solutions of linear or nonlinear differential equations. In order to get the dynamic equations for the RK-4 method, the following expressions for y and z are chosen: 1 yn n n y ξ ϕ ∞ = =∑ (85) 1 zn n n z ξ ϕ ∞ = =∑ (86) Substituting Eqs. (85) and (86) into Eqs. (42) and (43) and by using the orthogonal property, the following dynamic equations can be obtained: 2 3 2 2 2 2 1 1 1 1 1 1 1 1 1 1 2 1 1 2 4 1 2 2 1 5 2 2 2 3 2 2 2 1 2 1 2 6 1 2 2 1 2 7 2 2 2 13 1 1 1 1 1 1 1 16 1 ˆ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ( )] [ ( y y D U y y y y z y y y z z y y y z y y y z y y y z z y y z y y y y z z z y c a C C C C C C C ξ ξ ω ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + − + + + + + + + + + + + + + + + + + + ɺɺ ɺ ɶ ɶ ɶ ɶ ɶ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺ 2 2 1 2 1 2 1 2 1 2 1 2 1 2 17 2 1 1 1 1 1 1 15 2 2 1 2 2 2 2 2 2 25 2 1 2 1 2 1 2 1 2 1 2 1 2 2 )] [ ( )] [ ( )] [ ( 2 2 y y y y y y z z z z z z y y y y z z z y y y y z z z y y y y y y y z z z z z C C C ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + + + + + + + + + + + + + + + + + + + ɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ ɺɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ ɺɺ 2 19 2 2 2 2 2 2 2 2 2 25 1 1 1 0 10 )] [ ( )] ˆy z y y y y z z z i t y D U C C q e G a dx ξ ξ ξ ξ ξ ξ ξ ξ ϕ Ω + + + + = + ∫ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɶɶ (87) 2 3 2 2 2 2 2 1 2 2 2 1 1 1 8 1 2 1 1 2 11 1 2 2 1 12 2 2 2 3 2 2 2 1 2 1 2 14 1 2 2 1 2 15 2 2 2 10 1 1 1 1 1 1 1 21 ˆ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ( )] [ y y D U y y y y z y y y z z y y y z y y y z y y y z z y y z y y y y z z z c a C C C C C C C ξ ξ ω ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + − + + + + + + + + + + + + + + + + + + ɺɺ ɺ ɶ ɶ ɶ ɶ ɶ ɶ ɺɺ ɺ ɺɺ ɺ ɶ 2 2 1 1 2 1 2 1 2 1 2 1 2 1 2 22 2 1 1 1 1 1 1 23 2 2 1 2 2 2 2 2 2 24 2 1 2 1 2 1 2 1 2 1 2 ( 2 2 )] [ ( )] [ ( )] [ ( 2 2 y y y y y y y z z z z z z y y y y z z z y y y y z z z y y y y y y y z z z z C C C ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + + + + + + + + + + + + + + + + + + + ɺɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ ɺɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ 1 2 26 2 2 2 2 2 2 2 2 2 24 1 2 2 0 20 )] [ ( )] ˆy z z y y y y z z z i t y D U C C q e G a dx ξ ξ ξ ξ ξ ξ ξ ξ ϕ Ω + + + + = + ∫ ɺɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɶɶ (88) Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 284 2 3 2 2 2 2 1 1 1 1 1 1 1 1 1 1 2 1 1 2 4 1 2 2 1 5 2 2 2 3 2 2 2 1 2 1 2 6 1 2 2 1 2 7 2 2 2 13 1 1 1 1 1 1 1 16 1 ˆ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ( )] [ ( z z A U z z z z y z z z y y z z z y z z z y z z z y y z z y z z z z y y y z c a C C C C C C C ξ ξ ω ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + − + + + + + + + + + + + + + + + + + + ɺɺ ɺ ɶ ɶ ɶ ɶ ɶ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺ 2 2 1 2 1 2 1 2 1 2 1 2 1 2 17 1 2 1 1 1 1 1 1 15 2 2 1 2 2 2 2 2 2 25 2 1 2 1 2 1 2 1 2 1 2 2 2 )] [ ( )] [ ( )] [ ( 2 2 z z z z z z y y y z z z z z z z y y y z z z z y y y z z z z z z z y y y y y C G C C ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + + + + + + + + + + + + + + + + + + + ɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ ɺɺ ɶ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ 1 2 19 2 2 2 2 2 2 2 2 2 25 1 1 1 0 10 )] [ ( )] ˆz y z z z z y y y i t z A U C C q e G a dx ξ ξ ξ ξ ξ ξ ξ ξ ϕΩ + + + + = + ∫ ɺɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɶɶ (89) 2 3 2 2 2 2 2 1 2 2 2 1 1 1 8 1 2 1 1 2 11 1 2 2 1 12 2 2 2 3 2 2 2 1 2 1 2 14 1 2 2 1 2 15 2 2 2 10 1 1 1 1 1 1 1 21 ˆ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ( )] [ z z A U z z z z y z z z y y z z z y z z z y z z z y y z z y z z z z y y y c a C C C C C C C ξ ξ ω ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + − + + + + + + + + + + + + + + + + + + ɺɺ ɺ ɶ ɶ ɶ ɶ ɶ ɶ ɺɺ ɺ ɺɺ ɺ ɶ 2 2 1 1 2 1 2 1 2 1 2 1 2 1 2 22 2 1 1 1 1 1 1 23 2 2 1 2 2 2 2 2 2 24 2 1 2 1 2 1 2 1 2 1 2 ( 2 2 )] [ ( )] [ ( )] [ ( 2 2 z z z z z z z y y y y y y z z z z y y y z z z z y y y z z z z z z z y y y y C C C ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + + + + + + + + + + + + + + + + + + + ɺɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ ɺɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɺɺ ɺ ɺ ɺɺ ɺɺ ɺ ɺ 1 2 26 2 2 2 2 2 2 2 2 2 24 1 2 2 0 20 )] [ ( )] ˆz y y z z z z y y y i t z A U C C q e G a dx ξ ξ ξ ξ ξ ξ ξ ξ ϕΩ + + + + = + ∫ ɺɺ ɶ ɺɺ ɺ ɺɺ ɺ ɶ ɶɶ (90) where the fgh are the mode shape integrations and are shown in the Appendix. Eqs. (87)-(90) can be solved simultaneously by using RK-4. (a) y-dir., 1st mode (E�0 and Ei�0) (b) y-dir., 2nd mode (E�6 and Ei�6) Fig. 8 Phase plot of diameter ratio = 1/0.75, j = k. ll, exciting the y-dir., 1st mode (a) y-dir., 1st mode (E�0 and Ei�0) (b) y-dir., 2nd mode (E�6 and Ei�6) Fig. 9 Phase plot of diameter ratio = 1/0.75, j = k. mn, exciting the y-dir., 1st mode The phase plots are shown in Figs. 8-11. It is recorded that Figs. 8-9 are the cases for the 1st mode excited (the excitation frequency is the 1st mode’s natural frequency). Figs. 10-11 are the cases for the 2nd mode excited (the excitation frequency is the 1st mode’s natural frequency). The Displ. and Vel. represent the E�0,�6,�0 and Ei�0,�6,�0 , respectively. The converged displacements in time domain in Figs. 8-11 agree with the amplitudes form fixed point plots and show the I.R. happens in these cases. Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 285 (a) y-dir., 1st mode (E�0 and Ei�0) (b) y-dir., 2nd mode (E�6 and Ei�6) (c) z-dir., 1st mode (E�0 and Ei�0) Fig. 10 Phase plot of diameter ratio = 1/0.75, j = k. ll, y-dir., exciting the 2nd mode It is noted that the frequency analysis of the uniform beam can also be done by using Eqs. (8)-(10) and apply MOMS, select the secular terms for different harmonics of different modes and solve the solvability conditions to make the fixed point plots, since the procedure is similar to the case of double-section beam, and will not detail here. (a) y-dir., 1st mode (E�0 and Ei�0) (b) y-dir., 2nd mode (E�6 and Ei�6) (c) z-dir., 1st mode (E�0 and Ei�0) Fig. 11 Phase plot of diameter ratio = 1/0.75, j = k. mn, y-dir., exciting the 2nd mode 4. Wind Pressure Effects 4.1. Diameter ratio and length ratio Section 3.2 showed that when the diameter ratio �̅0/�̅6 = 1/0.75 and the length ratio 5 = 0.33&0.51, the 1:3 I.R. internal resonance occurs in the 1st and 2nd modes. Various diameter ratios from 1/0.6 to 1/0.9 will be discussed in this Section. Fig. 12 shows that when the diameter ratios are less than 1/0.78, there will be no I.R. When the diameter ratios are greater than 1/0.78, there will be two possible combinations for the length ratios to trigger I.R. Since there is no integer multiple of the frequency ratio for the first 3 modes in the uniform beam, the I.R. will not occur. However, the primary resonance will be examined for the uniform beam in the following Section. 4.2. Wind pressure and the beam amplitude Section 4.1 finds that the I.R. could be triggered in some certain combinations of diameter ratios and beam sectional length ratios. The I.R. cannot happen in the uniform beam, because there is no integer multiple of frequency ratios. Next step is to study the effect of the dimensionless wind pressure "oT (Refer to Appendix for the definition) on this beam amplitude. For the case of a uniform beam, the value of the dimensionless wind pressure "o�0 varies from 1-10. The responding amplitudes in both y- and z-D.O.F. can be found from Eqs. (8)-(10) by using MOMS method. The aerodynamic force is treated as the external load in this study. Instead of Reynolds number, the dimensionless wind pressure "o�0 is considered in Figs. 13-17. Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 286 Fig. 13 demonstrates the primary resonance phenomenon in an isotropic uniform beam when the 1st mode’s natural frequency in the y-dir. is excited. Fig. 13 presents that the amplitudes in the y-dir. increases linearly when "o�0 increases when "o�0 p 7. The amplitudes in the z-dir. seem to be unexcited. However, in the cases of "o�0 q 7, the amplitudes in the z-dir. increase notably. This implies that the energy is transferring from one D.O.F. to another D.O.F. under these conditions. When "o�0 = 10, the amplitude in the z-dir. is larger than the y-dir. This is a typical primary resonance. The result shows the primary resonance can happen under such conditions. Fig. 12 Frequency ratio plot of a double section beam, exciting the y-dir. 1st mode Fig. 13 Beam amplitude-wind pressure plot of the uniform isotropic beam, exciting the y-dir. 1st mode’s natural frequency Fig. 14 indicates the beam amplitude-wind pressure plot of the double section beam �̅0/�̅6 = 1/0.75, 5 = 0.33, when the y-dir. 1st mode is excited. Fig. 15 reveals the beam amplitude-wind pressure plot of the double section beam�̅0/�̅6 = 1/0.75, 5 = 0.51, when the y-dir. 1st mode is excited. The amplitude of By1 is always larger than the other cases and increases linearly as the wind pressure increases. The case of 5 = 0.33 represents a more flexible beam than the case of 5 = 0.51. This is the reason that By1 has the largest amplitude than the other cases. Fig. 14 Beam amplitude-wind pressure plot of the double section beam �̅0/�̅6 = 1/0.75, 5 = 0.33, exciting the y-dir. 1st mode Fig. 15 Beam amplitude-wind pressure plot of the double section beam �̅0/�̅6 = 1/0.75, 5 = 0.51, exciting the y-dir. 1st mode Figs. 16-17 are the beam amplitude-wind pressure plots of the double section beam �̅0/�̅6 = 1/0.75 when the y-dir. 2nd mode is excited and when 5 = 0.33 and 5 = 0.51, respectively. Even if the 2nd mode is excited and By2 increases as the wind pressure increases, the unexcited D.O.F. amplitude By1 is still larger than that of the excited D.O.F. This implies the energy is transferring in between these two modes. The internal resonance is triggered. In comparing with the primary resonance, the 3:1 I.R. happens earlier (less wind pressure) than the primary resonance (evidenced by the growth of Bz1). The internal resonance has greater effects than the primary resonance and should be noticed under this forcing condition. Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 287 Fig. 16 Beam amplitude-wind pressure plot of the double section beam �̅0/�̅6 = 1/0.75, 5 = 0.33, exciting the y-dir. 2nd mode Fig. 17 Beam amplitude-wind pressure plot of the double section beam �̅0/�̅6 = 1/0.75, 5 = 0.51, exciting the y-dir. 2nd mode 5. Conclusions This study investigates the internal resonance of a double-section beam with cubic nonlinearities. This model can be applied in a wide range of engineering problems, such as rocket and missile structures. Different ratios of beam section lengths and diameters in a variety of external loads (aerodynamic forces or flight speeds of the rockets or missiles) are studied. A primary resonance occurs on a uniform free-free beam at certain flight speed. The three-to-one internal resonance is also triggered within the 1st and the 2nd modes in certain combinations of beam length and diameter ratios. The phase plots and the time marching numerical method are used to verify the semi-analytical results. The findings are concluded as follows: (1) For a uniform beam, and when wind pressure " �0 r 10, primary resonance occurs. (2) The internal resonance does not happen in the uniform beam case, because the frequency ratios of different modes in the uniform beam are not equal to integer multiples. (3) For the double sectional beam and the diameter ratio �̅0/�̅6 = 1/0.75, the length ratio 5 = 0.33&0.51, the 1:3 internal resonance in the 1st and 2nd modes occurs. (4) When the diameter ratios are less than 1/0.78, there will be no I.R. While the diameter ratios are greater than 1/0.78, there will be two possible combinations for the length ratios to trigger I.R. (5) For the case of the double section beam, the internal resonance has greater effects than the primary resonance and should be noticed in this forcing condition. This study reveals the fight conditions to trigger primary resonance or internal resonance. However, the vibration of the beam structure still exists. The investigations of vibration reduction on the 3D beam or rocket model deserve an extension study on this research. The tuned mass damper (TMD) or the damping ring may be added to see the damping effects on the vibration beam. Since the nonlinear phenomenon is difficult to be observed experimentally, the experimental setup for verification of the present predictions is suggested in future research. Acknowledgement This research was supported by the Ministry of Science and Technology of Taiwan, Republic of China (grant number: MOST 108-2218-E-006-021). Conflicts of Interest The authors declare no conflict of interest Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 288 Appendix The secular terms for S0, S6 − 2S0 harmonics of the y-dir. 1st mode (m = 1) are: 1 1 1 1 1 2 1 1 1 1 1 1 2 ( 2 )2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 ( 2 ) ( )2 1 2 2 3 1 1 1 1 1 1 1 1 2 4 2 2 1 ˆ2( ) ( ) 3 2 (2 ) y y y y y y y y y z y z z i i i i i y y y y D U y y y y y i i y y y y z z y z y z z y z i B e B e c a i B e B B e C B B e C B B B e C B B B e B B e C B B B e C B B e ς ς ς ς ς ς ς ς ς ς ς ς ς ω ω ς ω − − − − − − + − − − − + − − − + ′ ′− + − − − − − − + − − ɶ ɶ ɶ ɶ ɶ 2 1 1 2 1 2 2 1 2 1 1 2 1 1 2 ( 2 ) ( ) ( ) 5 1 2 2 6 2 1 2 2 1 2 7 ( 2 ) ( 2 )2 2 2 2 2 1 1 1 16 1 2 17 18 1 2 2 19 2 1 2 17 2 ( ) [3 ( ) 2 ] [ y z y y z z y z z y y y y y y i i i y z z y z z y z z i i i i y y y y y y y y y C B B B e C B B B e B B B e C B B e C B B e C C B B B e C B B e C ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω ω − − − − + − − − + + − − − + − − − + − − + + + + + + ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ 1 1 1 2 1 2 1 1 1 1 1 ( 2 ) ( 2 )2 2 2 2 2 1 2 2 20 19 1 1 1 16 1 2 1 2 17 1 1 2 18 ( 2 ) (2 2 2 2 1 2 2 20 1 1 1 1 1 1 16 1 1 2 2 ( )] 2 2 [(2 ) y y y y y y y y y z y i i i i y y y y y y y y y i i y y y y z z y z y z z B B B e C C B B e C B B e C B B e C B B B e C B B B e B B e C B B B e ς ς ς ς ς ς ς ς ς ς ς ω ω ω ω ω ω − − − − + − − + − − − − + − − − + + − − + − + + + ɶ ɶ ɶ ɶ ɶ ɶ ɶ 1 2 2 1 2 1 2 2 1 2 1 1 2 1 2 1 2 2 1 ) 17 ( 2 ) ( ) ( ) ( )2 2 2 1 18 2 1 2 2 1 2 19 2 1 1 2 17 ( ) ( 1 2 2 20 2 1 2 2 1 2 ( ) ] [ 2 ( z z y z y z z y z z y z z y y z z y z i i y z y z z y z z y z z i i i y z z y z z y z z C B B e C B B B e B B B e C B B B e C B B B e C B B B e B B B e ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω + − − − + − − − + + − − − + − − + − − − + + + + + + + + ɶ ɶ ɶ ɶ ɶ 2 1 1 1 1 1 1 2 2 1 2 1 2 2 1 ) 2 19 1 1 1 1 ( 2 ) ( ) ( 2 )2 2 2 2 1 1 16 2 1 2 2 20 1 2 1 1 2 17 1 2 1 18 ( ) ( 1 2 2 1 2 2 1 2 ) ] ( 2 ) 2 2 2 ( z y y z y y z z y z y z z y z i y z z i y z y z z y z z y z i i y z z y z z C B B B e B B e C B B B e C B B B e C B B e C B B B e B B B e ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω ω ω ω ω ω ω + − − − + − − − − + − − − + − − − + + − + − − + − − + ɶ ɶ ɶ ɶ ɶ 2 1 ) 19 1)z i T y C q e ς σ+ +ɶ ɶ The secular terms for S6, 3S0 harmonics of the y-dir. 2nd mode are: 2 2 2 1 2 2 1 1 1 1 2 1 1 2 33 2 2 2 2 2 2 1 2 2 1 8 1 1 2 9 2 2 10 ( 2 ) ( ) ( )2 1 1 8 1 1 2 1 1 2 11 2 ˆ2( ) ( ) 2 3 ( ) 2 y y y y y y y z y z z y z z i i i i i i y y y y D U y y y y y y y i i i y z y z z y z z y z i B e B e c a i B e B e C B B B e C B B e C B B e C B B B e B B B e C B B ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω ω ς ω − − − − − − − + − − + − − + + ′ ′− + − − − − − − − + − ɶ ɶ ɶ ɶ ɶ 2 2 2 2 1 2 2 2 1 2 1 1 12 ( 2 ) 32 2 3 2 2 2 2 2 10 1 1 21 1 1 2 22 23 32 2 2 3 2 2 1 1 2 22 2 2 24 2 1 1 21 1 1 1 2 23 (2 ) [ 2 ( )] [2 3 ] 2 y y y z y y y y y y i z i i i i y z z y z y y y y i i i i y y y y y y y y y B e C B B B e B B e C B e C B B B e C C B B B e C B B e C G B e C B B B e C ς ς ς ς ς ς ς ς ς ς ω ω ω ω − − − − + − − − − − − − + + + + + + + − − ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ 2 1 1 1 1 2 1 1 2 2 2 1 1 2 2 2 1 2 2 24 ( 2 ) ( ) ( )2 2 1 1 1 21 1 1 2 1 1 2 22 2 1 1 23 ( ) (2 2 1 2 1 1 23 2 1 1 2 1 1 2 [ ( ) 2 ] 2 [( y y z y z z y z z y y y z z i y y i i i i y z y z z y z z y z z i i i y z z y z z y z z B B e C B B e C B B B e B B B e C B B B e C B B B e C B B B e B B B e ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω ω ω ω − − + − − + − − + + − − − − + − − + + + + − + + ɶ ɶ ɶ ɶ ɶ ( ) 1 1 2 1 12 2 2 1 1 2 1 1 2 2 2 2 ) 22 2( 2 )2 2 2 2 2 2 2 2 24 1 1 1 21 ( ) ( ) 1 2 1 1 2 1 1 2 22 ( 2 )2 2 2 2 2 2 2 2 ) (2 ) ] 2 ( ) ( 2 ) y z z y zy y z y z z y z z y y z ii i y z z y z y z i i y z z y z z i i y z z y z C B B B e B B e C B B e C B B B e B B B e C B B B e B B e ς ς ς ς ςς ς ς ς ς ς ς ς ς ς ς ς ω ω ω ω + + − +− − − + − − + − − + + − − − + + + + + − + + − + ɶ ɶ ɶ ɶ ɶ 1 24 2 i T yC q e σ+ ɶ The secular terms for S0, S6 − 2S0 harmonics of the z-dir and 1st mode are: 1 1 1 1 1 2 1 1 1 1 1 21 ( 2 )2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 2 3 ( 2 ) ( )2 2 1 1 1 1 1 1 1 1 2 4 2 1 ˆ2( ) ( ) 3 2 (2 ) z z z z z z z z y z y yz i i i i i i z z z z A U z z z z z z z z i z y y z y z y y z y i B e B e c a i B e B B e C B B e C B B B e C B B B e B B e C B B B e C B B e ς ς ς ς ς ς ς ς ς ς ς ςς ω ω ς ω− − − − − − + − − − + − − − +− ′ ′− + − − − − − − + − − ɶ ɶ ɶ ɶ ɶ 2 1 1 2 1 2 2 1 2 1 1 2 1 1 2 ( 2 ) 5 1 2 2 6 ( ) ( ) ( 2 )2 2 2 2 1 2 2 1 2 7 1 1 1 16 1 2 17 18 ( 2 )2 2 1 2 2 19 2 1 2 17 2 ( ) [3 ( ) 2 ] [ z y z z y y z y y z z z z z z i z y y i i i i z y y z y y z z z z i i z z z z z C B B B e C B B B e B B B e C B B e C B B e C C B B B e C B B e C ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω ω − − − − + − − − + + − − − + − − − + − − + + + + + + ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ 1 1 1 2 1 2 1 1 1 11 2 2 1 2 2 20 19 1 1 1 16 ( 2 ) ( 2 )2 2 2 2 1 2 1 2 17 1 1 2 18 2 1 2 2 20 ( 2 ) (2 2 1 1 1 1 1 1 16 1 1 2 2 ( )] 2 2 [(2 ) z z z z z z z z y zz i i z z z z z i i i z z z z z z z i z y y z y z y y B B B e C C B B e C B B e C B B e C B B B e C B B B e B B e C B B B e ς ς ς ς ς ς ς ς ς ςς ω ω ω ω ω ω − − − − + − − + − − − + − − −− + + − − + − + + + ɶ ɶ ɶ ɶ ɶ ɶ ɶ 1 2 2 1 2 1 2 2 1 2 1 1 2 1 2 1 2 2 1 ) ( 2 )2 17 2 1 18 ( ) ( ) ( )2 2 1 2 2 1 2 19 2 1 1 2 17 1 2 2 20 ( ) ( 2 1 2 2 1 2 ( ) ] [ 2 ( y y z y z y y z y y z y y z z y y z y z y i i i z y y z y y z y y z y y i i z y y z y y C B B e C B B B e B B B e C B B B e C B B B e C B B B e B B B e ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω + − − − + − − − + + − − − + − − + − − − + + + + + + + + ɶ ɶ ɶ ɶ ɶ 2 1 11 1 1 2 2 11 2 1 2 2 ) ( 2 )2 2 19 1 1 1 1 1 1 16 ( ) ( 2 )2 2 2 1 2 1 2 2 20 1 2 1 1 2 17 1 2 1 18 ( ) ( 1 2 2 1 2 2 1 2 ) ] ( 2 ) 2 2 2 ( y z yz z y y z yz z y y z i z y y z y i z y y z y y z y i i z y y z y y C B B B e B B e C G B B B e C B B B e C B B e C B B B e B B B e ς ς ςς ς ς ς ς ςς ς ς ς ς ω ω ω ω ω ω ω + − − +− − − − + − −− − + − − − + + − + − − + − − + ɶ ɶ ɶ ɶ ɶ ɶ 1 2 1 ) 19 1)y y i T z C q e ς ς σ+ +ɶ ɶ The secular terms for S6, 3S0 harmonics of the z-dir and 2nd mode are: Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 289 2 2 2 1 2 2 1 1 1 1 2 1 1 2 33 2 2 2 2 2 2 1 2 2 1 8 1 1 2 9 2 2 10 ( 2 ) ( ) ( )2 1 1 8 1 1 2 1 1 2 11 2 ˆ2( ) ( ) 2 3 ( ) 2 z z z z z z z y z y y z y y i i i i i i z z z z A U z z z z z z z i i i z y z y y z y y z y i B e B e c a i B e B e C B B B e C B B e C B B e C B B B e B B B e C B B − − − − − − − + − − + − − + + ′ ′− + − − − − − − − + − ɶ ɶ ɶ ɶ ɶ ς ς ς ς ς ς ς ς ς ς ς ς ς ς ω ω ς ω 2 2 22 1 2 2 2 1 2 1 1 12 ( 2 ) 32 2 3 2 2 2 2 2 10 1 1 21 1 1 2 22 23 32 2 2 3 2 2 2 1 1 2 22 2 2 24 1 1 21 1 1 1 2 23 1 (2 ) [ 2 ( )] [2 3 ] 2 z z yz z z z z z z i y ii i i z y y z y z z z z i i i i z z z z z z z z z B e C B B B e B B e C B e C B B B e C C B B B e C B B e C B e C B B B e C − − − +− − − − − − − − + + + + + + + − − ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ς ς ςς ς ς ς ς ς ς ω ω ω ω ω 2 1 1 1 1 2 1 1 2 2 1 1 2 12 2 2 2 24 ( 2 ) ( ) ( )2 2 1 1 1 21 1 1 2 1 1 2 22 2 1 1 23 ( ) (2 2 1 2 1 1 23 2 1 1 2 1 1 2 [ ( ) 2 ] 2 [( z z y z y y z y y z z y y zz i z z i i i i z y z y y z y y z y y i ii z y y z y y z y y B B e C B B e C B B B e B B B e C B B B e C B B B e C B B B e B B B e − − + − − + − − + + − − − + − −− + + + + − + + ɶ ɶ ɶ ɶ ɶ ς ς ς ς ς ς ς ς ς ς ς ς ς ςς ω ω ω 1 2 2 2 1 12 1 1 2 1 1 2 2 22 ) 22 ( 2 ) ( 2 )2 2 2 2 2 2 2 2 24 1 1 1 21 ( ) ( ) 1 2 1 1 2 1 1 2 22 ( 2 )2 2 2 2 2 2 2 2 ) (2 ) ] 2 ( ) ( 2 ) y y z y z yz z y y z y y z yz i ii z y y z y z y i i z y y z y y ii z y y z y C B B B e B B e C B B e C B B B e B B B e C B B B e B B e C + + − − + − +− − − + − − + + − − +− + + + + − + + − + ɶ ɶ ɶ ɶ ɶ ς ς ς ς ς ςς ς ς ς ς ς ς ς ςς ω ω ω ω 1 24 2 i T zq e+ ɶ σ where ( )1 2 0 2 1 , 1,m m G dx m ϕ = = ∫ ɶ 　 (A1) ( ) 1 0 1 2 0 , 1, 2 m m m m q dx q dx m ϕ ϕ == ∫ ∫ ɶ 　 (A2) 1 3 2 1 1 1 1 1 1 1 1 10 ( 4 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′ ′′ ′′′ ′= + +∫ɶ ɶ (A3) 1 2 2 2 1 1 1 2 1 1 2 1 2 1 2 1 1 1 2 1 2 10 (3 4 4 4 2 )iv iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′ ′′ ′′′ ′ ′′ ′′′ ′ ′ ′= + + + + +∫ɶ ɶ (A4) 1 2 2 3 1 1 1 2 1 2 2 2 1 2 2 2 1 1 2 2 2 10 (3 4 4 4 2 )iv iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′ ′′ ′′′ ′ ′′ ′′′ ′ ′ ′= + + + + +∫ɶ ɶ (A5) 1 2 2 4 1 1 1 2 1 1 2 1 2 1 2 1 1 2 1 1 2 10 (2 4 3 )iv iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′= + + + + +∫ɶ ɶ (A6) 1 2 5 1 1 1 2 2 1 1 2 1 1 1 2 10 ( 3 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′= + + +∫ɶ ɶ (A7) 1 2 6 1 1 1 2 1 2 2 1 2 2 1 2 20 ( 3 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′= + + +∫ɶ ɶ (A8) 1 2 2 7 1 1 1 2 2 1 2 2 1 2 1 2 2 2 2 1 2 10 (2 3 4 )iv iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′ ′ ′′ ′′′ ′= + + + + +∫ɶ ɶ (A9) 1 3 2 8 2 2 1 1 1 1 1 10 ( 4 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′ ′′ ′′′ ′= + +∫ɶ ɶ (A10) 1 2 2 9 2 2 1 2 1 1 2 1 2 1 2 1 1 1 2 1 2 10 (3 4 4 4 2 )iv iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′ ′′ ′′′ ′ ′′ ′′′ ′ ′ ′= + + + + +∫ɶ ɶ (A11) 1 3 2 10 2 2 2 2 2 2 2 20 ( 4 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′ ′′ ′′′ ′= + +∫ɶ ɶ (A12) 1 2 2 11 2 2 1 2 1 1 2 1 2 1 2 1 1 2 1 1 2 10 (2 4 3 )iv iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′= + + + + +∫ɶ ɶ (A13) 1 2 12 2 2 1 2 2 1 1 2 1 1 1 2 10 ( 3 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′= + + +∫ɶ ɶ (A14) Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 290 1 3 2 13 1 1 2 2 2 2 2 20 ( 4 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′ ′′ ′′′ ′= + +∫ɶ ɶ (A15) 1 2 14 2 2 1 2 1 2 2 1 2 2 1 2 20 ( 3 )iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′= + + +∫ɶ ɶ (A16) 1 2 2 15 2 2 1 2 2 1 2 2 1 2 2 2 1 1 2 2 2 10 ( 3 4 )iv iv C G dxϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ′′ ′′ ′ ′′ ′′′ ′′ ′ ′′′ ′ ′′ ′′′ ′ ′ ′= + + + + +∫ɶ ɶ (A17) 1 1 21 2 2 16 1 1 1 1 1 10 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A18) 1 1 21 17 1 1 1 1 2 1 1 20 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A19) 1 1 21 2 2 18 1 1 2 1 2 10 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A20) 1 1 21 19 1 1 2 1 2 2 1 20 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A21) 1 1 21 2 2 20 1 1 1 2 1 20 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A22) 1 1 21 2 2 21 2 2 1 1 1 10 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A23) 1 1 21 22 2 2 1 1 2 1 1 20 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A24) 1 1 21 2 2 23 2 2 2 1 2 10 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A25) 1 1 21 2 2 24 2 2 2 2 2 20 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A26) 1 1 21 2 2 25 2 1 2 2 2 20 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A27) 1 1 21 26 2 2 2 1 2 2 1 20 1 ( ) i i i i i i i i x x x x i i i i l l l l i C G dx dx dx dx dxϕ ϕ ϕ ϕ ϕ ϕ ϕ − −=  ′′ ′ ′ ′ ′ ′ ′= +   ∑∫ ∫ ∫ ∫ ∫ɶ ɶ (A28) References [1] E. Özkaya, “Non-linear transverse vibrations of a simply supported beam carrying concentrated masses,” Journal of Sound and Vibration, vol. 257, no. 3, pp. 413-424, October 2002. [2] J. S. Mundrey, Railway track engineering, Tata McGraw-Hill, 2000. [3] T. Phuoc Nguyen, D. Trung Pham, and P. Hoa Hoang, “Effects of foundation mass on dynamic responses of beams subjected to moving oscillators,” Journal of Vibroengineering, Vol. 22, no. 2, pp. 280-297, March 2020. [4] T. P. Chang, “Nonlinear free vibration analysis of nano-beams under magnetic field based on non-local elasticity theory,” Journal of Vibroengineering, vol. 18, no. 3, pp. 1912-1919, May 2016. [5] Z. Zhang, S. R. Nielsen, F. Blaabjerg, and D. Zhou, “Dynamics and control of lateral tower vibrations in offshore wind turbines by means of active generator torque,” Energies, vol. 7, no. 11, pp. 7746-7772, November 2014. [6] A. H. Nayfeh and D. T. Mook, Nonlinear oscillations, Wiley-Interscience Publication, New York, pp.95-160, 1995. Advances in Technology Innovation, vol. 5, no. 4, 2020, pp. 270-291 291 [7] P. F. Pai, “Nonlinear flexural-flexural-torsional dynamics of metallic and composite beams,” Ph.D. dissertation, Department of Engineering Science and Mechanics, Virginia Polytechnic Institute and State University, 1990. [8] S. Stoykov and P. Ribeiro, “Stability of nonlinear periodic vibrations of 3D beams,” Nonlinear Dynamics, vol. 66, no. 3, pp. 335, August 2011. [9] W. T. van Horssen and G. J. Boertjens, “On mode interactions for a weakly nonlinear beam equation,” Nonlinear Dynamics, vol. 17, no. 1, pp. 23-40, September 1998. [10] W. T. van Horssen and G. J. Boertjens, “An asymptotic theory for a weakly nonlinear beam equation with a quadratic perturbation,” SIAM Journal on Applied Mathematics, vol. 60, no. 2, pp. 602-632, 2000. [11] Y. R. Wang, C. K. Feng, and S. Y. Chen, “Damping effects of linear and nonlinear tuned mass dampers on nonlinear hinged-hinged beam,” Journal of Sound and Vibration, vol. 430, pp. 150-173, 2018. [12] Y. R. Wang and W. C. Hsiao, “Vibration reduction of damping rings on 3D nonlinear multi-loaded slender beams,” Journal of Chinese Society of Mechanical Engineers, vol. 40, no.4, pp. 327-339, June 2019. [13] A. Tekin, E. Özkaya, and S. M. Bağdatlı, “Three-to-one internal resonance in multiple stepped beam systems,” Applied Mathematics and Mechanics, vol. 30, no. 9, pp. 1131-1142, September 2009. [14] A. H. Nayfeh and P. F. Pai, Linear and nonlinear structural mechanics, New York: John Wiley & Sons, 2004. [15] Z. Ji and J. W. Zu, “Method of multiple scales for vibration analysis of rotor-shaft systems with non-linear bearing pedestal model,” Journal of Sound and Vibration, vol. 218, no. 2, pp. 293-305, November 1998. [16] A. H. Nayfeh and S. A. Nayfeh, “On nonlinear modes of continuous systems,” Journal of Vibration and Acoustics, vol. 116, no. 1, pp. 129-136, January 1994. [17] Q. Mao, “Free vibration analysis of multiple-stepped beams by using Adomian decomposition method,” Mathematical and Computer Modeling, vol. 54, no.1-2, pp. 756-764, March 2011. [18] D. Jeyakumar, K. K. Blswas, and B. Nageswara Rao, “Stage separation dynamic analysis of upper stage of a multistage launch,” Mathematical and Computer Modeling, vol. 41, no. 8-9, pp. 849-866, April-May 2005. [19] G. Srinivas and M. V. S. Prakash, “Aerodynamics and flow characterization of multistage rockets,” IOP Conference Series: Materials Science and Engineering, vol. 197, no. 1, article number 012077, July 2017. [20] W. T. van Horssen, “An asymptotic theory for a class of initial-boundary value problems for weakly nonlinear wave equations with an application to a model of the galloping oscillations of overhead transmission lines,” SIAM Journal of Applied Mathematics, vol. 48, no. 6, pp. 1227-1243, 1988. [21] “WayBack Machine,” https://web.archive.org/web/20150824115334/. [22] “European Space Agency,” http://wsn.spaceflight.esa.int/docs/EUG2LGPr3/EUG2LGPr3-6-SoundingRockets.pdf/. [23] K. Torabi, H. Afshari, and H. Najafi, “Vibration analysis of multi-step Bernoulli-Euler and Timoshenko beams carrying concentrated masses,” Journal of Solid Mechanics, vol. 5, no. 4, pp. 336-349, 2013. Copyright© by the authors. Licensee TAETI, Taiwan. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC) license (https://creativecommons.org/licenses/by-nc/4.0/).