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

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.


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][3][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 elastomechanical 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.

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.

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 V(, )/ 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 V(, ) stand for displacements along the coordinate axes  and , respectively, where  is the time.For Timoshenko beam, the total slope of the beam, the bending moment   (, ), and the shearing force   (, ) are given by V (, )  =   (, ) +   (, ) , (, ) =     (, ) =    ( V (, )  −   (, )) . ( 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  2   (, )/ 2 , 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 The high-order term ( 2 ) 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 (, ) Substituting the expressions for bending moment (2) and shear force (3) into the two governing equations ( 8) and ( 9) gives 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 V(, ); the other mode is the transverse shearing deformation   (, ), as measured by the difference (V(, )/ −   (, )).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   (, )/,  3   (, )/ 3 , and  3   (, )/ 2 ; (2) derive the expression of   (, )/ as a function of V(, ) from (10); by differentiating this expression first with respect to  two times and later with respect to  two times, find expressions for  3   (, )/ 3 and  3   (, )/ 2 as functions of V(, ); (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 V(, ).After some simple algebra, the latter reads [23] as 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: Herein, û(), V(), θ (), M (), Q (), and Q () are the amplitudes.It should be noted that   ,   and  are functions of  only; therefore Substituting ( 13) in ( 7) and ( 12) and performing the indicated differentiations yield where   ≡ √  / is the radius of gyration.From the derivation, it is clear that û() and V(), 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: The shear forces   (, ) and   (, ) at any section of the beam are related to the deflections (, ) and V(, ) and rotation angle   (, ) by Using ( 13), (18a) and (18b) can be written as 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: where   0 ,  0 , and  0 are reference quantities.In order to simplify the notation,   , , and  are consistently used for   , , and , respectively.Inserting ( 19) into ( 15) and ( 16), one gets where

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 where ℎ = cosh( 1 ), ℎ = sinh( 1 ),  = cos( 2 ) and  = sin( 2 ).Considering a beam segment th of length   , the input end 0 may be assigned to  = 0 resulting in Therefore, the coefficient vector can be expressed as For the beam segment output end ( + 1) at  =   , (35) and (38) give Applying the continuity equations as this results in where  The components V  (,  = 1, 2, 3, 4) of the transfer matrix U   are listed in the appendix.For the longitudinal vibration, denoted by superscript , the state vector is given as Z  () = {ũ, Q }  , the coefficient vector C  = { 5 ,  6 }  , and where  = cos( ) and  = sin( ).Following the pervious procedures, the longitudinal transfer matrix U   of the th segment may be deduced as By combining the state vectors Z  () and Z  () in one column, that is, Z() = {ũ, Ṽ, θ , M , Q , Q }  , the transfer matrix of a beam vibrating longitudinally in the -direction and transversely in the -direction is 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 Z 0,1 , and Z +1, are the boundary state vectors.
The transfer equations of the  segments, respectively, are The overall transfer equation of system is given by where is 6 × 6 transfer matrix of the whole beam structure, which is obtained by multiplying the transfer matrices of every segment in sequence.

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.
Q 10,9 =  21 ũ0,1 +  22 Q 0,1 , (48e) For a free-free beam, the bending moment and forces at its left and right sides are zero, that is, M 10,9 = Q 10,9 = Q 10,9 = 0, resulting in lower three equations of (48a)-(48f) which may be written in matrix form as Non-trivial solutions for (50) require that its coefficient determinant vanishes, that is, where det U  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 Z 0,1 = {ũ (0) =?, Ṽ (0) =?, θ (0) =?, M (0) = 0, In ( 52), the values of M 0,1 , Q 0,1 , and Q 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 U  can be broken down into the product of three matrices-an orthogonal  ×  matrix U, an  ×  diagonal matrix D, and the transpose of an  ×  orthogonal matrix V [25]: where U  U = I; V  V = I.The columns of U are orthogonal eigenvectors of U  U   , the columns of V are orthonormal eigenvectors of U   U  , and D is a diagonal matrix containing the square roots of eigenvalues from U or V in descending order.Substitution of (53) in (50) yields The solutions Z  0,1 are given by the columns of V corresponding to zero singular values.By using (45) further, the values of the state vector are found for each of the other segments.

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

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, (≡ /  ≡  2  / 2   ) is the ratio of axial force to the buckling load (≡ /  ≡  2  / 2   ).A negative sign in  indicates a compressive force.

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: Percentage differences between the predicted and the measured values are shown in parentheses.

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.

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.

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 Figure 11 shows the effect of increased operating thrust on the first natural frequency at  = 0 (maximum mass) and  = 50 sec (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 ( = 50 sec) 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  = 50 sec, the critical buckling thrust (associated with zero frequency) is lower than that at  = 0.In turn, it will limit the maximum operating thrust.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 6 N| =0 and  = 1.2 × 10 6 N| =50 sec , 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].

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).for the constructive comments/suggestions that have been contributed to the improvement of the paper.

Figure 2 :
Figure 2: Free body diagram of beam element of infinitesimal length  in its deformed state.

Figure 4 :
Figure 4: Beam example with free-free boundary conditions.
), 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 ( = 50 sec) times.The axial compressive force is considered for a nominal value of thrust as shown in Figure10(b).

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

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

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

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

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