Transfer Matrix Method for the Determination of the Natural Vibration Characteristics of Realistic Thrusting Launch Vehicle—Part I

Abstract

The feasibility of using the transfer matrix method (TMM) to compute the natural vibration characteristics of a flexible rocket/satellite launch vehicle is explored theoretically. In the approach to the problem, a nonuniform free-free Timoshenko and Euler-Bernoulli beamlike structure is modeled. A provision is made to take into consideration the effects of shear deformation and rotary inertia. Large thrust-to-weight ratio leads to large axial accelerations that result in an axial inertia load distribution from nose to tail which causes the development of significant compressive forces along the length of the launch vehicle. Therefore, it is important to take into account this effect in the transverse vibration model. Once the transfer matrix of a single component has been obtained, the product of all component matrices composes the matrix of the entire structure. The frequency equation and mode shape are formulated in terms of the elements of the structural matrices. Flight test and analytical results validate the present TMM formulas.

1. Introduction

The transient mass and structural characteristics of a typical multistage rocket vehicle require the natural-vibration characteristics of the vehicle known at least for the ignition and burnout time of each stage of flight and frequently for other conditions such as those at Mach number of 1, maximum dynamic pressure, and minimum lift. Thrust-to-weight and length-to-diameter (slenderness) ratios of the launch vehicle are a continuing need in order to launch heavier payload, increasing the range and reducing the aerodynamic drag. As the slenderness ratio increases, the structure becomes more flexible and the natural frequency of transverse vibrations reduces. Knowing the vibration characteristic of the rocket is very important for the computation of its structural dynamic response, aeroelasticity, flight mechanics, and control. In addition, the presence of large thrust/acceleration leads to significant axial forces along the length of the rocket [1], which have a major impact on the transverse vibration frequencies of beams. Since inertial measuring units (IMUs) also sense the local body vibrations, the mass and stiffness nonuniformities plus the thrust action on elastic missiles can potentially influence their measurements and thus must be properly accounted for in an aeroelastic simulation [2].

A highly flexible launch vehicle can be modeled as a free-free beamlike structure with lateral vibration [2–4]. The selection criteria for either Euler-Bernoulli or Timoshenko beam theories are generally determined by some rules involving beam dimensions (slenderness ratio). The Euler-Bernoulli beam theory is well suited to model the behavior of flexure-dominated (or “long”) beams, whereas the Timoshenko theory applies for shear-dominated (or “short”) beams where the effect of shear deformation and rotary inertia has a significant influence on the vibration characteristics [5]. In this report [5], the authors have built a method of determining modal data of a non uniform beam with effects of shear deformation and rotary inertia. They applied the method to a typical space research launch vehicle and its fourth-stage configuration. The results prove that rotary inertia and shear deformation have a significant effect on the frequencies of the fourth stage but only minor effect on the overall vehicle frequencies. Most research on high flexible launch vehicle’s vibration characteristics does not consider the effect of thrust. However, large thrust has a great influence on the natural frequencies and mode shapes and has to be taken into account. The research on vibration characteristics of beams under the axial compressive/tensile force has a long history, where also applications to realistic missile/launch vehicles may be found [1, 2, 6].

In this paper, the transfer matrix method (TMM) will be applied which has been developed for a long time and widely used in structure mechanics of the linear and nonlinear systems [7]. Holzer initially applied TMM to solve problems of torsion vibrations of rods [8]. Myklestad applied TMM to determine the bending-torsion modes of beams [9]. Thomson applied TMM to more general vibration problems [10]. Pestel and Leckie listed transfer matrices for elasto-mechanical elements up to 12th order [11]. Rubin provided a general treatment of transfer matrices and their relation to other forms of frequency response matrices [12, 13]. Transfer matrices have been applied to a wide variety of engineering programs by a number of researchers, including Targoff [14], Lin and McDaniel [15], and McDaniel and Murthy [16]. Many articles studied and improved the finite element transfer matrix for structure dynamics [17, 18]. Recently, Rui et al. [19] developed the transfer matrix method of multibody system (MS-TMM) and discrete time transfer matrix method (DT-TMM) for multibody system dynamics.

The purpose of this paper is to present a solution procedure for computing the natural vibrations of a non uniform beamlike structure, which satisfies the specific needs of a rocket vehicle designer. Because rocket vehicles are exposed to large axial loads caused by the thrust and inertia forces of the vehicle masses, the stability problem is important. Problems of this kind are non conservative. To ensure proper calculation of the critical load, the dynamic stability criterion must be applied which requires a vibration analysis of the loaded system. The analysis developed in this paper is based on the transfer matrix method of multibody system [19] and based on Timoshenko and Euler-Bernoulli beam theory to calculate the frequencies and mode shapes of realistic missiles. The method will be proven to accurately describe not only the mode shape but also slope angle of the elastic curve, rotation angle, mode moments, and mode shears. It is inherently applicable to complex structures such as multistage launch vehicle.

The text is organized as follows. In Section 2, the general theorem of the method is shown. In Section 3, some results calculated by TMM and the other method are given that can validate the proposed method. The conclusions and future works are presented in Section 4.

2. Modeling of a Launch Vehicle as a Nonuniform Beamlike Structure

Figure 1 shows a multistage launch vehicle where 𝑀_𝑚 and 𝐿_𝑚 are the total mass and the total length, respectively. The distributions of axial rigidity, flexural rigidity, mass per unit length, and shear rigidity along the longitudinal 𝑥-axis are given by 𝐸𝐴(𝑥), 𝐸𝐼_𝑧(𝑥), 𝜌𝐴(𝑥), and 𝜅_𝑠𝐺𝐴(𝑥), respectively. Herein, 𝜌 is the mass density, 𝐴 is the beam cross-sectional area, 𝐸 is Young’s modulus of elasticity, 𝐺 is the shear modulus, and 𝜅_s is the Timoshenko shear coefficient depending on the cross section of the beam. The axial inertia moment of the cross section is 𝐼_𝑧. The vehicle is subjected to an engine thrust 𝑇 at its rear end and leads to an axial force 𝑆(𝑥) (positive in tension) in each cross section. In the present study, the structure is divided into 𝑛 segments, where 𝐸𝐴(𝑥), 𝐸𝐼_𝑧(𝑥), 𝜌𝐴(𝑥), and 𝜅_𝑠𝐺𝐴(𝑥) distributions are uniform in each segment. Note that 𝑆(𝑥) will not be uniform due to the accelerated motion of the missile [2]. In the following, the segments will be considered as axially loaded beams, and transfer matrices are derived for this kind of elements.

Figure 1 

Multistage launch vehicle scheme.

2.1. Equations of Motion of Axially Loaded Beams

The Timoshenko model is an extension of the Euler-Bernoulli model by taking into account two additional effects: shear force effect and rotary motion effect [20]. In any beam, except one subject to pure bending only, additional deflection due to the shear stress occurs. The exact solution of the beam vibration problem requires this deflection to be considered. The angle between beam axis and 𝑥-axis, given by slope 𝜕𝑣(𝑥,𝑡)/𝜕𝑥 of the elastic curve, can be decomposed into two parts, namely, the rotation angle 𝜃_𝑧(𝑥,𝑡) of the cross section about the 𝑧-axis due to pure bending and the additional angle 𝛾_𝐺(𝑥,𝑡) caused by shear strain, as shown in Figure 2. It shows that 𝑢(𝑥,𝑡) and 𝑣(𝑥,𝑡) stand for displacements along the coordinate axes 𝑥 and 𝑦, respectively, where 𝑡 is the time.

Figure 2 

Free body diagram of beam element of infinitesimal length 𝑑𝑥 in its deformed state.

For Timoshenko beam, the total slope of the beam, the bending moment 𝑀_𝑧(𝑥,𝑡), and the shearing force 𝑄_𝑦(𝑥,𝑡) are given by 𝜕𝑣(𝑥,𝑡)𝜕𝑥=𝜃_𝑧(𝑥,𝑡)+𝛾_𝐺(),𝑀𝑥,𝑡_𝑧(𝑥,𝑡)=−𝐸𝐼_𝑧𝜕𝜃_𝑧(𝑥,𝑡),𝑄𝜕𝑥_𝑦(𝑥,𝑡)=𝜅_𝑠𝐺𝐴𝛾_𝐺(𝑥,𝑡)=𝜅_𝑠()𝐺𝐴𝜕𝑣𝑥,𝑡𝜕𝑥−𝜃_𝑧(𝑥,𝑡). (1)(2)(3)

Another factor that affects in the lateral vibration of the beam, which is neglected in Euler-Bernoulli’s model, is the fact that each beam section rotates slightly in addition to its lateral motion when the beam deflects. The influence of the beam section rotation is taken into account through the moment of inertia 𝐽𝜕²𝜃_𝑧(𝑥,𝑡)/𝜕²𝑡, which modifies the equation of moments acting on an infinitesimal beam element 𝑑𝑥. Herein, 𝐽=𝜌𝐼_𝑧𝑑𝑥 denotes the rotational inertia.

The equations of motion of the Timoshenko beam under axial loading can be obtained by three approaches: the principle of virtual work, Hamilton’s principle, and the classical dynamic equilibrium method. The derivation of these equations is based on the assumption that the shear force acting on the cross sections is normal to the deflected axis of the beam column and therefore inclined at an angle 𝛾_𝐺(𝑥,𝑡) with respect to the vertical direction. In this section, the linearized equations of motion describing the transverse vibration of a free-free beam will be derived by using the Timoshenko beam theory and the classical dynamic equilibrium method. In order to describe the effect of the axial force, it is assumed that the axial compression force is tangential to the slope of the beam. It could be also assumed that this force is normal to the direction of the shear force, where in [21] both cases have been considered and for both cases the equations of motion have been derived. But in [9] nothing has been said on which method is more accurate. However, in [22] it has been indicated that the equations of motion, which follow from the first assumption, are more accurate. Therefore, also in this paper it is assumed that the axial force is tangential to the slope of the beam. Consider the free body diagram of Figure 2; the equations of motion along the axial and transverse directions and the rotational equation of motion are given by 𝑄_𝑦()()𝑥,𝑡𝜕𝑣𝑥,𝑡−𝑄𝜕𝑥^𝑦(𝑥,𝑡)+𝜕𝑄_𝑦(𝑥,𝑡)×()𝜕𝑥𝑑𝑥𝜕𝑣𝑥,𝑡+𝜕𝜕𝑥²𝑣(𝑥,𝑡)𝜕𝑥²𝑑𝑥−𝑆+𝑆+𝜕𝑆𝜕𝜕𝑥𝑑𝑥=𝜌𝐴𝑑𝑥²𝑢(𝑥,𝑡)𝜕𝑡²,−𝑄_𝑦(𝑥,𝑡)+𝑄_𝑦(𝑥,𝑡)+𝜕𝑄_𝑦(𝑥,𝑡)()𝜕𝑥𝑑𝑥−𝑆𝜕𝑣𝑥,𝑡+𝜕𝑥𝑆+𝜕𝑆()𝜕𝑥𝑑𝑥𝜕𝑣𝑥,𝑡+𝜕𝜕𝑥²𝑣(𝑥,𝑡)𝜕𝑥²𝜕𝑑𝑥=𝜌𝐴𝑑𝑥²𝑣(𝑥,𝑡)𝜕𝑡²,𝑄_𝑦(𝑥,𝑡)𝑑𝑥2+𝑄_𝑦(𝑥,𝑡)+𝜕𝑄_𝑦(𝑥,𝑡)𝜕𝑥𝑑𝑥𝑑𝑥2+𝑀_𝑧()−𝑀𝑥,𝑡_𝑧(𝑥,𝑡)+𝜕𝑀_𝑧(𝑥,𝑡)()𝜕𝑥𝑑𝑥+𝑆𝜕𝑣𝑥,𝑡𝜕𝑥𝑑𝑥2−𝑆+𝜕𝑆()𝜕𝑥𝑑𝑥𝜕𝑣𝑥,𝑡+𝜕𝜕𝑥²𝑣(𝑥,𝑡)𝜕𝑥²𝑑𝑥𝑑𝑥2=𝜌𝐼_𝑧𝜕𝑑𝑥²𝜃_𝑧(𝑥,𝑡)𝜕𝑡². (4)(5)(6)

The high-order term (𝑑𝑥²) appearing in conjunction with the compression term is neglected [22]. If the segment is considered as bar vibrating in axial direction with a uniform strain in every segment, the uniaxial force resulting from stress is 𝑆=𝐸𝐴(𝑥)𝜕𝑢/𝜕𝑥. After simplifying calculations, (4)–(6) reduce to 𝜕𝐸𝐴²𝑢𝜕𝑥²𝜕=𝜌𝐴²𝑢(𝑥,𝑡)𝜕𝑡²,𝜕𝑄_𝑦(𝑥,𝑡)𝜕𝜕𝑥+𝑆²𝑣(𝑥,𝑡)𝜕𝑥²𝜕=𝜌𝐴²𝑣(𝑥,𝑡)𝜕𝑡²,−𝜕𝑀_𝑧(𝑥,𝑡)𝜕𝑥+𝑄_𝑦(𝑥,𝑡)=𝜌𝐼_𝑧𝜕²𝜃_𝑧(𝑥,𝑡)𝜕𝑡². (7)(8)(9)

Substituting the expressions for bending moment (2) and shear force (3) into the two governing equations (8) and (9) gives 𝜅_𝑠𝜕𝐺𝐴²𝑣(𝑥,𝑡)𝜕𝑥²−𝜕𝜃_𝑧(𝑥,𝑡)𝜕𝜕𝑥+𝑆²𝑣(𝑥,𝑡)𝜕𝑥²𝜕=𝜌𝐴²𝑣(𝑥,𝑡)𝜕𝑡²,𝐸𝐼_𝑧𝜕²𝜃_𝑧(𝑥,𝑡)𝜕𝑥²+𝜅_𝑠()𝐺𝐴𝜕𝑣𝑥,𝑡𝜕𝑥−𝜃_𝑧(𝑥,𝑡)=𝜌𝐼_𝑧𝜕²𝜃_𝑧(𝑥,𝑡)𝜕𝑡². (10)(11)

These results are the governing equations for the Timoshenko beam under axial force. There are two modes of deformation coupled by the governing equations. One mode of deformation is simply the transverse deflection of the beam as measured by 𝑣(𝑥,𝑡); the other mode is the transverse shearing deformation 𝛾_𝐺(𝑥,𝑡), as measured by the difference (𝜕𝑣(𝑥,𝑡)/𝜕𝑥−𝜃_𝑧(𝑥,𝑡)). It is more convenient to deal with a higher-order differential equation than with two-coupled equations of lower order. The steps that need to be performed to merge the two equations (10) and (11) into a single equation are the following: (1) differentiate (11) with respect to 𝑥 resulting in partial derivatives 𝜕𝜃_𝑧(𝑥,𝑡)/𝜕𝑥, 𝜕³𝜃_𝑧(𝑥,𝑡)/𝜕𝑥³, and? 𝜕³𝜃_𝑧(𝑥,𝑡)/𝜕𝑥𝜕𝑡²; (2) derive the expression of? 𝜕𝜃_𝑧(𝑥,𝑡)/𝜕𝑥 as a function of 𝑣(𝑥,𝑡) from (10); by differentiating this expression first with respect to 𝑥 two times and later with respect to 𝑡 two times, find expressions for 𝜕³𝜃_𝑧(𝑥,𝑡)/𝜕𝑥³ and 𝜕³𝜃_𝑧(𝑥,𝑡)/𝜕𝑥𝜕𝑡²?as functions of 𝑣(𝑥,𝑡); (3) substitute these three partial derivatives of? 𝜃_𝑧(𝑥,𝑡) into the differentiated equation (11) in order to have a single differential equation in the only unknown dependent variable 𝑣(𝑥,𝑡). After some simple algebra, the latter reads [23] as 𝐸𝐼_𝑧𝑆1+𝜅_𝑠𝜕𝐺𝐴⁴𝑣(𝑥,𝑡)𝜕𝑥⁴−𝐸𝐼_𝑧𝜌𝐴𝜅_𝑠𝐺𝐴+𝜌𝐼_𝑧+𝑆𝜌𝐼_𝑧𝜅_𝑠×𝜕𝐺𝐴₄𝑣(𝑥,𝑡)𝜕𝑥²𝜕𝑡²𝜕−𝑆²𝑣(𝑥,𝑡)𝜕𝑥²𝜕+𝜌𝐴²𝑣(𝑥,𝑡)𝜕𝑡²+𝜌²𝐼_𝑧𝐴𝜅_𝑠𝜕𝐺𝐴⁴𝑣(𝑥,𝑡)𝜕𝑡⁴=0. (12)

Equations (7) and (12) can be transformed into an ordinary differential equation (ODE) by separation of variables. When the beam performs one of its natural modes of oscillation, that is, steady-state vibration with angular frequency 𝜔, the displacements, rotation angle, moment, and shears take the following form: 𝑢()=𝑢(𝑥)(),𝑣()=𝑣(𝑥)(),𝜃𝑥,𝑡cos𝜔𝑡𝑥,𝑡cos𝜔𝑡_𝑧(𝑥,𝑡)=𝜃_𝑧(𝑥)(),𝑀cos𝜔𝑡_𝑧(𝑥,𝑡)=𝑀_𝑧(𝑥)(),𝑄cos𝜔𝑡_𝑥(𝑥,𝑡)=𝑄_𝑥(𝑥)(),𝑄cos𝜔𝑡_𝑦(𝑥,𝑡)=𝑄_𝑦(𝑥)cos(𝜔𝑡). (13)

Herein, 𝑢(𝑥), 𝑣(𝑥), 𝜃_𝑧(𝑥), 𝑀_𝑧(𝑥), 𝑄_𝑥(𝑥), and 𝑄_𝑦(𝑥) are the amplitudes. It should be noted that 𝜅_𝑠𝐺𝐴,𝐸𝐼_𝑧 and 𝐸𝐴 are functions of 𝑥 only; therefore 𝜕𝜅_𝑠𝐺𝐴≡𝑑𝜅𝜕𝑥^𝑠𝐺𝐴,𝜕𝑑𝑥𝐸𝐼_𝑧≡𝑑𝜕𝑥𝐸𝐼_𝑧,𝜕()𝑑𝑥𝐸𝐴≡𝑑()𝜕𝑥𝐸𝐴.𝑑𝑥 (14)

Substituting (13) in (7) and (12) and performing the indicated differentiations yield 𝑑²𝑢(𝑥)𝑑𝑥²+𝑚𝜔²𝑢(𝑥)𝑆𝐸𝐴=0,1+𝜅_𝑠𝑑𝐺𝐴⁴𝑣(𝑥)𝑑𝑥⁴+−𝑆𝐸𝐼_𝑧+𝜌𝐴𝜔²𝜅_𝑠+𝐺𝐴𝜌𝐴𝑟²_𝑔𝜔²𝐸𝐼_𝑧𝑆1+𝜅_𝑠𝑑𝐺𝐴²𝑣(𝑥)𝑑𝑥²,+𝜌𝐴²𝑟²_𝑔𝜔⁴𝐸𝐼_𝑧𝜅_𝑠−𝐺𝐴𝜌𝐴𝜔²𝐸𝐼_𝑧𝑣(𝑥)=0, (15)(16) where 𝑟_𝑔≡𝐼_𝑧/𝐴 is the radius of gyration. From the derivation, it is clear that 𝑢(𝑥) and 𝑣(𝑥), which are found as solutions of (15) and (16), are the longitudinal and transverse vibration mode shapes of the considered beam for a corresponding frequency 𝜔. If the Euler-Bernoulli model is adopted instead, the rotary inertia and shear deformation effects are neglected. The corresponding equation of motion can be derived by setting 𝜌𝐼_𝑧 (represented in rotary inertia) equal to zero in (6) and by letting the modulus 𝐺 tend to infinity (which is equivalent to imposing vanishing shear deformation). With these simplifications, (16) reduces to the equation of an axially loaded Euler-Bernoulli beam: 𝑑⁴𝑣(𝑥)𝑑𝑥⁴−𝑆𝐸𝐼_𝑧𝑑²𝑣(𝑥)𝑑𝑥²−𝑚𝜔²𝐸𝐼_𝑧𝑣(𝑥)=0. (17)

The shear forces 𝑄_𝑥(𝑥,𝑡) and 𝑄_𝑦(𝑥,𝑡) at any section of the beam are related to the deflections 𝑢(𝑥,𝑡)?and?𝑣(𝑥,𝑡) and rotation angle 𝜃_𝑧(𝑥,𝑡) by𝑄_𝑥()()𝑥,𝑡=−𝐸𝐴𝜕𝑢𝑥,𝑡,𝑄𝜕𝑥_𝑦(𝑥,𝑡)=𝜕𝑀_𝑧(𝑥,𝑡)+𝜕𝑥𝑚𝑟²_𝑔𝜕²𝜃_𝑧(𝑥,𝑡)𝜕𝑡²()+𝑆𝜕𝑣𝑥,𝑡.𝜕𝑥 (18a)(18b) Using (13),?(18a) and (18b) can be written as 𝑄_𝑥(𝑥)𝜕𝑢(𝑥)=−𝐸𝐴,𝑄𝜕𝑥_𝑦(𝑥)=𝜕𝑀_𝑧(𝑥)+𝜕𝑥𝑚𝑟²_𝑔𝜕²𝜃_𝑧(𝑥)𝜕𝑡²𝜕𝑣(𝑥)+𝑆.𝜕𝑥 (18c)(18d)

In the model for a launch vehicle (Figure 1), the beam will represent the segments. For 𝑗th segment of the rocket, let 𝑥_𝑗=𝑥/𝐿_𝑚 be the dimensionless length coordinate such that it takes values 0≤𝑥_𝑗≤𝑙_𝑗(𝑙_𝑗=𝑙_𝑗/𝐿_𝑚), where 𝑙_𝑗 is the length of the segment. To express (15) and (16) in dimensionless form, the following quantities are defined: 𝑢𝑢=𝐿_𝑚,𝑣𝑣=𝐿_𝑚,𝐷^𝑛=𝑑^𝑛𝑑𝑥^𝑛_𝑗,𝑘²=𝑆_𝑗𝐿²_𝑚𝐸𝐼_𝑧0,𝑅²=𝑟_𝑔𝑗𝐿_𝑚²,𝑝²=𝜌𝐴₀𝜔²𝐿⁴_𝑚𝐸𝐼_𝑧0,𝑀_𝑧𝑗(𝑥)=𝑀_𝑧𝑗(𝑥)𝐿_𝑚𝐸𝐼_𝑧0𝛼,𝑄_𝑦𝑗(𝑥)=𝑄_𝑦𝑗(𝑥)𝐿²_𝑚𝐸𝐼_𝑧0𝛼,𝑄_𝑥𝑗(𝑥)=𝑄_𝑥𝑗(𝑥)𝐸𝐴₀,𝜃_𝑧𝑗(𝑥)=𝜃_𝑧(𝑥),𝛽=𝜌𝐴₀𝜔²𝐿²_𝑚𝐸𝐴₀,𝛼=1+𝑘²𝜅_𝑠,𝐺𝐴𝐸𝐼_𝑧𝑗=𝐸𝐼_𝑧0𝐸𝐼_𝑧,𝜅_𝑠𝐺𝐴_𝑗=𝐸𝐼_𝑧0𝐿²_𝑚𝜅_𝑠𝐺𝐴,𝐸𝐴_𝑗=𝐸𝐴₀𝐸𝐴,𝜌𝐴_𝑗=𝜌𝐴₀𝑚, (19) where 𝐸𝐼_𝑧0, 𝐸𝐴₀,?and?𝜌𝐴₀ are reference quantities. In order to simplify the notation, 𝐸𝐼_𝑧, 𝐸𝐴,?and?𝜌𝐴 are consistently used for 𝐸𝐼_𝑧, 𝐸𝐴,?and?𝑚, respectively. Inserting (19) into (15) and (16), one gets 𝐷²[𝑢(𝑥)]+𝛽²𝑚𝑢(𝐸𝐴𝑥)𝐷=0,₄[𝑣(𝑥)]+𝑏₀𝐷²[𝑣(𝑥)]+𝑐₀𝑣(𝑥)=0, (20)(21) where 𝑏₀=−𝑘²𝐸𝐼_𝑧𝛼+𝑚𝑝²1𝜅_𝑠+𝐺𝐴𝛼𝑅²𝐸𝐼_𝑧,𝑐₀=𝑝²𝐸𝐼_𝑧𝛼𝑝²𝑅²𝑚²𝜅_𝑠−𝐺𝐴𝑚. (22)

2.2. Solution to the Beam Differential Equations (Longitudinal and Flexural)

The characteristic equation of the flexural beam (21) has the following form: 𝑧²+𝑏₀𝑧+𝑐₀=0, (23) where 𝑧=𝐷²[𝑣(𝑥)], with solutions 𝑧₁=(−𝑏₀+√Δ)/2, and 𝑧₂=(−𝑏₀−√Δ)/2 where Δ=𝑏²₀−4𝑐₀. It can be seen from (22) thatΔ>0 and z₂̃<0 for all 𝜔. However, z₁̃ has two sign possibilities.(1)𝑧₁>0⇔√Δ>𝑏₀⇒𝑏²₀−4𝑐₀>𝑏₀⇒𝑐₀<0. That means that 𝑝²<𝜅_𝑠𝐺𝐴/𝑅²𝑚 or 𝜔²<𝜅_𝑠𝐺𝐴/𝜌𝐼. (2)𝑧₁<0⇔𝜔>𝜅_𝑠𝐺𝐴/𝜌𝐼. In the present study, only 𝑧₁>0 will be considered. Therefore, the roots of (23) are 𝐷[𝑣(𝑥)]∈{±𝑧₁,±𝑖−𝑧₂}. This gives a general solution in the following form 𝑣_𝑗(𝑥)=𝐶^′_1𝑗𝑒𝑧_1𝑥+𝐶^′_2𝑗𝑒⁻𝑧_1𝑥+𝐶^′_3𝑗𝑒^𝑖𝑧_2𝑥+𝐶^′_4𝑗𝑒^−𝑖𝑧_2𝑥, (24)where the constants 𝐶^′₁, 𝐶^′₂, 𝐶^′₃, and 𝐶^′₄ are to be determined from the boundary conditions of the beam. Noting that sinh𝑥=(1/2)(𝑒^𝑥−𝑒^−𝑥), cosh𝑥=(1/2)(𝑒^𝑥+𝑒^−𝑥), sin𝑥=(1/2𝑖)(𝑒^𝑖𝑥−𝑒^−𝑖𝑥), and cos𝑥=(1/2)(𝑒^𝑖𝑥+𝑒^−𝑖𝑥), (24) can be rewritten as 𝑣_𝑗(𝑥)=𝐶_1𝑗cosh𝜆_1𝑗𝑥+𝐶_2𝑗sinh𝜆_1𝑗𝑥+𝐶_3𝑗cos𝜆_2𝑗𝑥+𝐶_4𝑗sin𝜆_2𝑗𝑥, (25) where 𝜆²_1𝑗=𝑧₁=−𝑏₀+√Δ2;𝜆²_2𝑗=𝑧₂=𝑏₀+√Δ2. (26) Substituting (3) into (1) and the result in (2) yields the variable 𝑀_𝑧(𝑥,𝑡). Then substitute (8) into 𝑀_𝑧(𝑥,𝑡), use (13), differentiate it with respect to 𝑥, and substitute it back into (18d). This finally gives 𝑄_𝑦(𝑥)1−𝜌𝐴𝑟²_𝑔𝜔²𝜅_𝑠𝐺𝐴=−𝐸𝐼_𝑧𝑆1+𝜅_𝑠𝜕𝐺𝐴³𝑣(𝑥)𝜕𝑥³−𝐸𝐼_𝑧𝜌𝐴𝜔²𝜅_𝑠𝐺𝐴+𝜌𝐴𝑟²_𝑔𝜔²𝜕𝑣(𝑥)−𝑆.𝜕𝑥 (27) The shear force can be obtained by the substitution of (19) into (27) as 𝑄_𝑦𝑗(𝑥)=𝐴_𝑄𝐷³[𝑣(𝑥)]+𝐵_𝑄𝐷[𝑣(𝑥)], (28) where 𝐴_𝑄=−𝐸𝐼_𝑧1−𝑝²𝑅²𝑚/𝜅_𝑠,𝐵𝐺𝐴_𝑄=−𝑝²𝐸𝐼_𝑧𝑚/𝜅_𝑠𝐺𝐴+𝑝²𝑅²𝑚−𝑘²𝛼1−𝑝²𝑅²𝑚/𝜅_𝑠.𝐺𝐴 (29) Substituting the result (28) into (1) and (2) and after several manipulations, one gets 𝜃_𝑧𝑗(𝑥)=𝐴_𝜃𝐷³[𝑣(𝑥)]+𝐵_𝜃𝐷[𝑣(𝑥,𝑀)]_𝑧𝑗(𝑥)=𝐴_𝑀𝐷⁴[𝑣(𝑥)]+𝐵_𝑀𝐷²[𝑣(𝑥)], (30) where 𝐴_𝜃𝛼=−𝜅_𝑠𝐴𝐺𝐴^𝑄,𝐵_𝜃𝛼=1−𝜅_𝑠𝐵𝐺𝐴^𝑄,𝐴_𝑀=𝐸𝐼_𝑧𝜅_𝑠𝐴𝐺𝐴^𝑄,𝐵_𝑀=𝐸𝐼_𝑧𝜅_𝑠𝐵𝐺𝐴^𝑄−𝐸𝐼_𝑧𝛼. (31) For the longitudinal vibration (20), the characteristic equation has the following form: 𝑧²+𝛽²=0⇒𝑧=±𝑖𝛽, (32) where 𝑧=𝐷[𝑢(𝑥)]. This gives the solution 𝑢_𝑗(𝑥)=𝐶_5𝑗cos𝛽_𝑗𝑥+𝐶_6𝑗sin𝛽_𝑗𝑥. (33) In dimensionless form, the shear force 𝑄_𝑥𝑗(𝑥) at section 𝑗th can be obtained by substituting (33) into (18c) and (19): 𝑄_𝑥𝑗(𝑥)=𝐴_𝑢𝐷[𝑢(𝑥)],where𝐴_𝑢=−𝐸𝐴𝛽_𝑗. (34)

2.3. Transfer Matrix Method (TMM) for Distributed Parameter Models

The introduction of the transfer matrix method into the continuum modeling process provides a very useful tool to facilitate the application of distributed parameter models to more complex configurations. In TMM approach, the continuity relations of displacement, rotation angle, moment, and shear at the interface of adjacent beam segments are performed. For the flexural beam, denoted by superscript 𝑓, (25), and (28)–(31) can be written in matrix form as 𝐙^𝑓(𝑥)=𝐁^𝑓(𝑥)𝐂^𝑓. (35) Herein, 𝐙^𝑓(𝑥) is the state vector at an arbitrary position 𝑥 containing the model variables, that is, 𝐙^𝑓(𝑥)={𝑣,𝜃_𝑧,𝑀_𝑧,𝑄_𝑦}^𝑇. The coefficient vector 𝐂^𝑓={𝐶₁,𝐶₂,𝐶₃,𝐶₄}^𝑇 summarizes unknown constants to be adopted to boundary conditions and𝐁^𝑓(𝑥)=𝐴𝑐ℎ𝑠ℎ𝑐𝑠_𝜃𝜆³₁+𝐵_𝜃𝜆₁𝑠ℎ𝐴_𝜃𝜆³₁+𝐵_𝜃𝜆₁𝑐ℎ𝐴_𝜃𝜆³₂−𝐵_𝜃𝜆₂𝑠−𝐴_𝜃𝜆³₂+𝐵_𝜃𝜆₂𝑐𝐴_𝑀𝜆⁴₁+𝐵_𝑀𝜆²₁𝑐ℎ𝐴_𝑀𝜆⁴₁+𝐵_𝑀𝜆²₁𝑠ℎ𝐴_𝑀𝜆⁴₂−𝐵_𝑀𝜆²₂𝑐𝐴_𝑀𝜆⁴₂−𝐵_𝑀𝜆²₂𝑠𝐴_𝑄𝜆³₁+𝐵_𝑄𝜆₁𝑠ℎ𝐴_𝑄𝜆³₁+𝐵_𝑄𝜆₁𝑐ℎ𝐴_𝑄𝜆³₂−𝐵_𝑄𝜆₂𝑠−𝐴_𝑄𝜆³₂+𝐵_𝑄𝜆₂𝑐, (36)where 𝑐ℎ=cosh(𝜆₁𝑥), 𝑠ℎ=sinh(𝜆₁𝑥), 𝑐=cos(𝜆₂𝑥) and 𝑠=sin(𝜆₂𝑥). Considering a beam segment 𝑗th of length 𝑙_𝑗, the input end 0 may be assigned to 𝑥=0 resulting in 𝐙^𝑓_0,𝑗=𝐁^𝑓(0)_𝑗𝐂^𝑓_𝑗. (37) Therefore, the coefficient vector can be expressed as 𝐂^𝑓_𝑗=𝐁^𝑓(0)_𝑗⁻¹𝐙^𝑓_0,𝑗. (38) For the beam segment output end (𝑗+1) at 𝑥=𝑙_𝑗, (35) and (38) give 𝐙^𝑓_𝑗+1,𝑗=𝐁^𝑓𝑙_𝑗𝐂^𝑓_𝑗=𝐁^𝑓𝑙_𝑗𝐁^𝑓(0)_𝑗⁻¹𝐙^𝑓_0,𝑗. (39) Applying the continuity equations as 𝐙^𝑓_𝑗,𝑗=𝐙^𝑓_{𝑗−1,𝑗}, (40) this results in 𝐙^𝑓_𝑗+1,𝑗=𝐔^𝑓_𝑗𝐙^𝑓_0,𝑗, (41) where 𝐔^𝑓_𝑗=𝐁^𝑓𝑙_𝑗𝐁^𝑓(0)_𝑗⁻¹=𝑣₁₁𝑣₁₂𝑣₁₃𝑣₁₄𝑣₂₁𝑣₂₂𝑣₂₃𝑣₂₄𝑣₃₁𝑣₃₂𝑣₃₃𝑣₃₄𝑣₄₁𝑣₄₂𝑣₄₃𝑣₄₄^𝑓_𝑗. (42) The components 𝑣_𝑙𝑚(𝑙,𝑚=1,2,3,4) of the transfer matrix 𝐔^𝑓_𝑗?are listed in the appendix. For the longitudinal vibration, denoted by superscript 𝑙, the state vector is given as 𝐙^𝑙(𝑥)={𝑢,𝑄_𝑥}^𝑇, the coefficient vector 𝐂^𝑙={𝐶₅,𝐶₆}^𝑇, and 𝐁^𝑙(𝑥)=^𝑐𝑠_{−𝐴𝑢𝑠𝐴𝑢𝑐}, where 𝑐=cos(𝛽𝑥)and𝑠=sin(𝛽𝑥). Following the pervious procedures, the longitudinal transfer matrix 𝐔^𝑙_𝑗 of the 𝑗th segment may be deduced as 𝐔^𝑙_𝑗=𝑐𝑠𝐴_𝑢−𝐴_𝑢𝑠𝑐^𝑙_𝑗=𝑢₁₁𝑢₁₂𝑢₂₁𝑢₂₂^𝑙_𝑗. (43) By combining the state vectors 𝐙^𝑙(𝑥) and 𝐙^𝑓(𝑥) in one column, that is, 𝐙(𝑥)={𝑢,𝑣,𝜃_𝑧,𝑀_𝑧,𝑄_𝑥,𝑄_𝑦}^𝑇, the transfer matrix of a beam vibrating longitudinally in the 𝑥-direction and transversely in the 𝑦-direction is 𝑢𝐔=¹¹000𝑢₁₂00𝑣₁₁𝑣₁₂𝑣₁₃0𝑣₁₄0𝑣₂₁𝑣₂₂𝑣₂₃0𝑣₂₄0𝑣₃₁𝑣₃₂𝑣₃₃0𝑣₃₄𝑢₂₁000𝑢₂₂00𝑣₄₁𝑣₄₂𝑣₄₃0𝑣₄₄. (44) A vibrating beam comprised of 𝑛-segments, see Figure 3, is used as an example to show how to deduce the overall transfer equation and overall transfer matrix of the system. There are 𝑛+1 connection and boundary points where 𝐙_0,1, and 𝐙_𝑛+1,𝑛 are the boundary state vectors.

Figure 3 

Beam with 𝑛 segments.

The transfer equations of the 𝑛 segments, respectively, are 𝐙_2,1=𝐔₁𝐙_0,1,𝐙_3,2=𝐔₂𝐙_1,2,…,𝐙_𝑛+1,𝑛=𝐔_𝑛𝐙_{𝑛−1,𝑛}. (45) The overall transfer equation of system is given by 𝐙_𝑛,𝑛+1=𝐔_𝑛𝐙_{𝑛,𝑛−1}=𝐔_𝑛𝐔_𝑛−1𝐙_{𝑛−1,𝑛−2}=⋯=𝐔𝐙_0,1, (46) where 𝐔=^𝑛_𝑘=1𝐔_𝑘=𝐔_𝑛𝐔_𝑛−1⋯𝐔₁=𝑢₁₁000𝑢₁₂00𝑣₁₁𝑣₁₂𝑣₁₃0𝑣₁₄0𝑣₂₁𝑣₂₂𝑣₂₃0𝑣₂₄0𝑣₃₁𝑣₃₂𝑣₃₃0𝑣₃₄𝑢₂₁000𝑢₂₂00𝑣₄₁𝑣₄₂𝑣₄₃0𝑣₄₄ (47) is 6×6 transfer matrix of the whole beam structure, which is obtained by multiplying the transfer matrices of every segment in sequence.

2.4. Determination of Natural Frequencies and Mode Shapes of a Free-Free Beam

The components of the transfer matrix are all functions of the natural frequency of the system. The overall transfer equation (47) only involves boundary state vectors, and the state vectors at all other connection points do not appear. The state vectors at the boundary are composed of displacements, rotation angles, bending moments, and shear and axial forces which are partly known and partly unknown. For free ends on both sides, the force and moment components are zero. Figure 4 shows an example sketch of an elastic beam divided into 9 segments and free boundaries at sections 0 and 10.

Figure 4 

Beam example with free-free boundary conditions.

From (46), one obtains the relationship between left side 𝐙_0,1={𝑢𝑣𝜃_𝑧𝑀_𝑧𝑄_𝑥𝑄_𝑦}^𝑇_0,1 and right side 𝐙_10,9={𝑢𝑣𝜃_𝑧𝑀_𝑧𝑄_𝑥𝑄_𝑦}^𝑇_10,9 to take the following form:𝑢_10,9=𝑢₁₁𝑢_0,1+𝑢₁₂𝑄_𝑥0,1,𝑣_10,9=𝑣₁₁𝑣_0,1+𝑣₁₂𝜃_𝑧0,1+𝑣₁₃𝑀_𝑧0,1+𝑣₁₄𝑄_𝑦0,1,𝜃_𝑧10,9=𝑣₂₁𝑣_0,1+𝑣₂₂𝜃_𝑧0,1+𝑣₂₃𝑀_𝑧0,1+𝑣₂₄𝑄_𝑦0,1,𝑀_𝑧10,9=𝑣₃₁𝑣_0,1+𝑣₃₂𝜃_𝑧0,1+𝑣₃₃𝑀_𝑧0,1+𝑣₃₄𝑄_𝑦0,1,𝑄_𝑥10,9=𝑢₂₁𝑢_0,1+𝑢₂₂𝑄_𝑥0,1,𝑄_𝑦10,9=𝑣₄₁𝑣_0,1+𝑣₄₂𝜃_𝑧0,1+𝑣₄₃𝑀_𝑧0,1+𝑣₄₄𝑄_𝑦0,1. (48a)(48b)(48c)(48d)(48e)(48f)For a free-free beam, the bending moment and forces at its left and right sides are zero, that is,𝑀_𝑧0,1=𝑄_𝑥0,1=𝑄_𝑦0,1𝑀=0,_𝑧10,9=𝑄_𝑥10,9=𝑄_𝑦10,9=0, (49a)(49b)resulting in lower three equations of (48a)–(48f) which may be written in matrix form as 𝐙^′_10,9=𝐔^′𝐙^′_0,1⇒000_10,9=0𝑣₃₁𝑣₃₂𝑢₂₁000𝑣₄₁𝑣₄₂𝑢𝑣𝜃_𝑧0,1. (50) Non-trivial solutions for (50) require that its coefficient determinant vanishes, that is, Δ(𝜔)=det𝐔^′=0⇒0𝑣₃₁𝑣₃₂𝑢₂₁000𝑣₄₁𝑣₄₂=0, (51) where det𝐔^′ is a function of 𝜔. Equation (51) is called the characteristic equation or frequency equation whose roots are the natural frequencies 𝜔 of the undamped system. These roots are determined by “guess and check.” For a given frequency, the algorithm checks the sign of the determinant Δ(𝜔) in (51). Starting from an initial frequency, it is increased by Δ𝜔 and the sign of Δ(𝜔+Δ𝜔) is rechecked. A root occurs when there is a sign change. If there is no sign change, the incremental Δ𝜔 is increased by a factor (e.g., 1.6) and search is continued. Else 𝜔 is overestimated; Δ𝜔 is decreased by a factor (e.g., 3.2) and the calculation starts from the old frequency before the sign change occurs. The frequency is obtained if Δ𝜔−𝜀≤0 where 𝜀 (e.g., 10^-9) is the required precision for the frequency estimation. It is anticipated that the second root is at least twice the value of the first. Hence, a larger Δ𝜔 in the algorithm is assumed for the iteration of the second root, and so on. Once the converged value of a natural frequency 𝜔_𝑘 (where 𝑘 is a mode counter) is determined, one may generate the corresponding free vibration mode shapes. For a free-free beam, the full initial state variables (or initial parameters) are given by 𝐙_0,1=𝑢(0)=?,𝑣(0)=?,𝜃_𝑧(0)=?,𝑀_𝑧(0)𝑄=0,_𝑥(0)=0,𝑄_𝑦(0)=0^𝑇_0,1. (52) In (52), the values of 𝑀_𝑧0,1, 𝑄_𝑥0,1, and 𝑄_𝑦0,1 are equal to zero (49a), while those of 𝑢_0,1, 𝑣_0,1, and 𝜃_0,1 are unknowns. Using Singular Value Decomposition (SVD), the unknown variables are determined from the solution of linear homogenous equations (50). SVD is based on a theorem from linear algebra, which says that a rectangular 𝑚×𝑛 matrix 𝐔^′ can be broken down into the product of three matrices—an orthogonal 𝑚×𝑚 matrix 𝐔, an 𝑚×𝑛 diagonal matrix 𝐃, and the transpose of an 𝑛×𝑛 orthogonal matrix 𝐕 [25]: 𝐔^′_𝑚×𝑛=𝐔_𝑚×𝑚𝐃_𝑚×𝑛𝐕^𝑇_𝑛×𝑛, (53) where 𝐔^𝑇𝐔=𝐈; 𝐕^𝑇𝐕=𝐈. The columns of 𝐔 are orthogonal eigenvectors of 𝐔^′𝐔^′𝑇, the columns of 𝐕 are orthonormal eigenvectors of 𝐔^′𝑇𝐔^′, and 𝐃 is a diagonal matrix containing the square roots of eigenvalues from 𝐔 or 𝐕 in descending order. Substitution of (53) in (50) yields 𝐙^′_10,9=𝐔𝐃𝐕^𝑇𝐙^′_0,1⋯𝐔𝐃𝐕^𝑇𝐙^′_0,1=𝟎. (54) The solutions 𝐙^′_0,1are given by the columns of? 𝐕 corresponding to zero singular values. By using (45) further, the values of the state vector are found for each of the other segments.

2.5. Axial Force for Constant Thrust Trajectory

The axial force 𝑆 is defined as compressive force acting on the beam for the 𝑗th rocket segment and may be calculated at the centre of gravity of the segment by summing up the inertias of all the segments preceding the current one and taking the average inertia force for the current segment [2, 6]. For the constant thrust trajectory, it can be presented as 𝑆_𝑗(𝑥)=𝑇𝑀_𝑚𝑀_𝑚−^𝑗−1_𝑘=1𝜌𝐴₀𝑚_𝑘𝑙_𝑘−12𝜌𝐴₀𝑚_𝑗𝑙_𝑗. (55)

3. Comparison Examples

3.1. Free-Free Uniform Euler-Bernoulli Beam under Axial Force

In the case where the shear deformation and rotary inertia are disregarded (see (17)), the free-free beam has exact analytical solutions. From Table 1 a perfect concordance [24] and TMM for only one spanwise segment can be observed. Herein, 𝛾(≡𝑇/𝑇_𝑐𝑟≡𝑇𝐿²_𝑚/𝜋²𝐸𝐼_𝑧) is the ratio of axial force to the buckling load 𝛾(≡𝑇/𝑇_𝑐𝑟≡𝑇𝐿²_𝑚/𝜋²𝐸𝐼_𝑧). A negative sign in 𝛾 indicates a compressive force.

Table 1 

Frequency parameter 𝑝² for a free-free Euler-Bernoulli square beam under axial load for 𝐸𝐼_𝑧=200×10⁹?N/m²; 𝐺=1×10⁹⁰?N/m²; 𝜌=7800?kg/m³; 𝑏=ℎ=0.01?m; 𝐿_𝑚=1?m.

3.2. Free-Free Uniform Timoshenko Beam

A good agreement can be also observed from Table 2 between measured [26], ANSYS [26], and the TMM values for a single segment. The geometry and material data are as follows: 𝐸𝐼_𝑧=72.66611×10⁹N/m²,𝐺=27.17481×10⁹N/m²,𝜌=2817.9111kg/m³,𝑏=0.01009m,ℎ=0.01004m,𝐿_𝑚𝑘=0.04004m,_𝑠=0.869877423. (56) Percentage differences between the predicted and the measured values are shown in parentheses.

Table 2 

Natural frequencies (Hz) according to measured values, ANSYS, and TMM of a free-free short Timoshenko aluminum square beam.

3.3. Free-Free Nonuniform Beam: Two-Stage Launch Vehicle

A numerical example has been carried out by using Euler-Bernoulli beams for a two-stage launch vehicle [27]. Figure 5 shows the axial variations of flexural stiffness 𝐸𝐼_𝑧 and mass per unit length 𝜌𝐴 distributions. Table 3 illustrates a comparison of TMM using different equal segments (i.e., uniform 𝐸𝐼_𝑧 and 𝜌𝐴 distributions along each segment) with in-flight frequencies of oscillation obtained from telemetry data with computed first-mode undamped natural frequencies of the free-free two-stage vehicle [27]. A very good concordance can be observed. The number of segments depends mainly on the physical properties and the complexity of the launch vehicle. However, for this example, 400 segments give good convergence as shown in Table 3.

Table 3 

A comparison of TMM using different numbers of segments with flight test and computed data.

(a)

(b)

Figure 5 

Variation of (a) flexural stiffness and (b) mass distribution, for a two-stage launch vehicle prior the launch [27].

3.4. Free-Free Nonuniform Beam: Multistage Launch Vehicle with the Effect of Shear Deformation and Rotary Inertia

A numerical example is also presented for an actual vehicle [5]. The vehicle has four stages. The computation is carried out with both Euler-Bernoulli and Timoshenko approaches for high aspect ratio (𝐿/𝐷≈25, for whole launch vehicle) and low aspect ratio (𝐿/𝐷≈5, for fourth stage only). The basic physical characteristics of the multistage vehicle required as input are illustrated in Figure 6 including the axial variations of flexural stiffness 𝐸𝐼_𝑧, mass per unit length 𝜌𝐴, shear rigidity 𝑘_𝑠𝐺𝐴, and radius of gyration 𝑟_𝑔. The number of segments used is 1000. The first three frequencies computed by TMM are in good agreement with [5] as shown in Table 4. The slope of the elastic curve is computed by using (1). Inspection of the data presented in Table 4 and Figures 7, 8, and 9 shows that rotary inertia and shear deformation have minor effects for 𝐿/𝐷≈25 only, while it has a significant effect on the data for 𝐿/𝐷≈5. Dominant effects are noted in reduction of frequencies, and large effects are pronounced on the moment and shear.

Table 4 

First three frequencies (rad/sec) of high and low aspect ratio for a multistage launch vehicle: comparison between [5] and TMM.

(a)

(b)

(c)

(d)

Figure 6 

Parameters of a multistage launch vehicle [5]: (a) flexural stiffness, (b) shear rigidity, (c) mass, and (d) radius of gyration.

(a)

(b)

(c)

(d)

Figure 7 

Effects of shear deformation and rotary inertia on (a) mode shape, (b) rotation angle, (c) moment, and (d) shear of first mode for 𝐿/𝐷≈25: ♦ Timoshenko model; —Euler-Bernoulli model.

Figure 8 

Effects of shear deformation and rotary inertia on the first three natural vibrations on mode shape ((a), (b), (c)) and rotation angle and slope of elastic curve ((d), (e), (f)) for upper stage of the launch vehicle (𝐿/𝐷≈5).

Figure 9 

Effects of shear deformation and rotary inertia on the first three natural vibrations on moment ((a), (b), (c)) and shear force ((d), (e), (f)) for upper stage of the launch vehicle (𝐿/𝐷≈5): - - - Timoshenko model; —Euler-Bernoulli model.

4. Typical Guided Rocket for Constant Thrust

A numerical example has also been carried out for a typical guided rocket with flexural rigidity, axial compressive force, and mass per unit length distributions as shown in Figures 10(a), 10(b), and 10(c). The structure is divided into 600 equal segments, and two moments of the flight time are selected, that is, initial (𝑡=0) and burn-out (𝑡=50sec) times. The axial compressive force is considered for a nominal value of thrust as shown in Figure 10(b).

(a)

(b)

(c)

Figure 10 

Variation of (a) flexural stiffness, (b) axial force due to constant nominal thrust = 1×106N, and (c) mass per unit length of a typical rocket distributed along the vehicle length for two selected times of flight (—𝑡=0 and - - - 𝑡=50sec).

Figure 11 shows the effect of increased operating thrust on the first natural frequency at 𝑡=0 (maximum mass) and 𝑡=50sec (minimum mass). It highlights the following general notes: (1) increasing the thrust magnitude tends to a reduction in the vehicle’s frequencies; (2) mass depletion (𝑡=50sec) leads to increase in the acceleration and the axial force distribution as shown in Figure 10(b), and consequently it exhibits higher vehicle’s frequency in comparison with maximum mass (𝑡=0) but lower allowable thrust range. This fact can be shown readily from Figure 11. For 𝑡=50sec, the critical buckling thrust (associated with zero frequency) is lower than that at 𝑡=0. In turn, it will limit the maximum operating thrust.

Figure 11 

Effect of operating thrust on the first natural frequency at two of selected flight times (—𝑡=0 and - - - 𝑡=50sec).

Figure 12 illustrates the effect of maximum operating thrust (in the vicinity of critical buckling thrust taken from Figure 11), that is, 𝑇=1.5×10⁶N|_𝑡=0 and 𝑇=1.2×10⁶N|_𝑡=50sec, on the first natural mode shape. For both cases, the deflections of the mode shape with thrust are decreased at the rear end in comparison without thrust, while increasing at the nose. It should be noted that this fact has been also concluded by Pourtakdoust and Assadian [2].

Figure 12 

Effect of two maximum thrust values on the first natural mode shape at two of selected flight times.

5. Conclusions

Starting from an elementary beam formulation for a nonuniform flexible beamlike structure, the free vibration analysis is performed using the Transfer Matrix Method (TMM). The free-free beam is subjected to variable axial compressive force and modeled by both Euler-Bernoulli and Timoshenko theory. The latter solution includes the influences of rotary inertia and shear deformation. For forward applications, the beam is divided into piecewise constant property segments. For the free vibration, the exact solution is obtained using the beam functions. Numerical examples for launch various vehicles are included to demonstrate the validity of TMM.

The methodology proposed in this research with inclusion of rotary inertia, shear deformation, and thrust effect allows for more realistic simulations of a flexible vehicle. The accuracy of this methodology has a great influence on the vibrational characteristics, control/navigation system, vehicle dynamics and aeroelastic analyses.

Ongoing research to be published in part II is devoted to modeling of Closed Structure Systems. This system consists of a multibeam-like structure, rigid joints, springs, and rigid bodies. Special attention is focused on how the transfer equations and transfer matrices of the global system can be developed conveniently (see Figure 13).

Figure 13 

Closed structure system.

Appendix

Considering (42), the components of the transfer matrix 𝐔^𝑓_𝑗 can be stated as follows: 𝑣₁₁=𝜆²₁𝐵_𝑀+𝐴_𝑚𝜆²₁𝑐+𝜆²₂𝐵_𝑀−𝐴_𝑚𝜆²₂𝑐ℎ𝜆²₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂,𝑣₁₂=−𝜆₁𝐵_𝑄+𝐴_𝑄𝜆²₁𝑠+𝜆₂𝐵_𝑄−𝐴_𝑄𝜆²₂𝑠ℎ𝜆₁𝜆₂𝜆²₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃,𝑣₁₃=−𝑐+𝑐ℎ𝜆²₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂,𝑣₁₄=𝜆₁𝐵_𝜃+𝐴_𝜃𝜆²₁𝑠+−𝐵_𝜃𝜆₂+𝐴_𝜃𝜆³₂𝑠ℎ𝜆₁𝜆₂𝜆²₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃,𝑣₂₁=𝜆₁𝜆₂𝜆₁𝐵_𝑀+𝐴_𝑚𝜆²₁−𝐵_𝜃+𝐴_𝜃𝜆²₂𝑠+𝐵_𝜃+𝐴_𝜃𝜆²₁𝜆₂𝐵_𝑀−𝐴_𝑚𝜆²₂×𝜆𝑠ℎ₂₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂⁻¹,𝑣₂₂=−𝐵_𝑄+𝐴_𝑄𝜆²₁𝐵_𝜃−𝐴_𝜃𝜆²₂𝑐+𝐵_𝜃+𝐴_𝜃𝜆²₁𝐵_𝑄−𝐴_𝑄𝜆²₂×𝜆𝑐ℎ₂₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃⁻¹,𝑣₂₃=𝜆₂𝐵_𝜃−𝐴_𝜃𝜆²₂𝑠+𝜆₁𝐵_𝜃+𝐴_𝜃𝜆²₁𝑠ℎ𝜆²₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂,𝑣₂₄=𝐵_𝜃+𝐴_𝜃𝜆²₁𝐵_𝜃−𝐴_𝜃𝜆²₂(𝑐−𝑐ℎ)𝜆²₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃,𝑣₃₁=𝜆²₁𝐵_𝑀+𝐴_𝑚𝜆²₁𝜆²₂−𝐵_𝑀+𝐴_𝑚𝜆²₂(𝑐−𝑐ℎ)𝜆²₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂,𝑣₃₂=𝜆₂𝐵_𝑄+𝐴_𝑄𝜆²₁𝐵_𝑀−𝐴_𝑀𝜆²₂𝑠+𝜆₁𝐵_𝑀+𝐴_𝑀𝜆²₁𝐵_𝑄−𝐴_𝑄𝜆²₂×𝜆𝑠ℎ₂₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃⁻¹,𝑣₃₃=𝜆²₂𝐵_𝑀−𝐴_𝑀𝜆²₂𝑐+𝜆²₁𝐵_𝑀+𝐴_𝑀𝜆²₁𝑐ℎ𝜆²₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂,𝑣₃₄=𝜆₂𝐵_𝜃+𝐴_𝜃𝜆²₁−𝐵_𝑀+𝐴_𝑀𝜆²₂𝑠−𝜆₁𝐵_𝑀+𝐴_𝑀𝜆²₁𝐵_𝜃−𝐴_𝜃𝜆²₂×𝜆𝑠ℎ₂₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃⁻¹,𝑣₄₁=𝜆₁𝜆₂𝜆₁𝐵_𝑀+𝐴_𝑚𝜆²₁−𝐵_𝑄+𝐴_𝑄𝜆²₂𝑠+𝐵_𝑄+𝐴_𝑄𝜆²₁𝜆₂𝐵_𝑀−𝐴_𝑚𝜆²₂×𝜆𝑠ℎ₂₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂⁻¹,𝑣₄₂𝐵=−^𝑄+𝐴_𝑄𝜆²₁𝐵_𝑄−𝐴_𝑄𝜆²₂(𝑐−𝑐ℎ)𝜆²₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃,𝑣₄₃=𝜆₂𝐵_𝑄−𝐴_𝑄𝜆²₂𝑠+𝜆₁𝐵_𝑄+𝐴_𝑄𝜆²₁𝑠ℎ𝜆²₁+𝜆²₂𝐵_𝑀+𝐴_𝑀𝜆₁−𝜆₂𝜆₁+𝜆₂,𝑣₄₄=𝐵_𝜃+𝐴_𝜃𝜆²₁𝐵_𝑄−𝐴_𝑄𝜆²₂𝑐−𝐵_𝑄+𝐴_𝑄𝜆²₁𝐵_𝜃−𝐴_𝜃𝜆²₂×𝜆𝑐ℎ₂₁+𝜆²₂𝐴_𝜃𝐵_𝑄−𝐴_𝑄𝐵_𝜃⁻¹. (A.1)

Acknowledgment

The authors express their indebtedness to Professor Dieter Bestle (Engineering Mechanics and Vehicle Dynamics, Brandenburg University of Technology, Germany) during his 5 months staying at their Institute of Launch Dynamics/NUST for the constructive comments/suggestions that have been contributed to the improvement of the paper.