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

Authors: Laith K. Abbas; Li Min-Jiao; Rui Xiao-Ting
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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 Section2, 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

Figure1 shows a multistage launch vehicle where Mm and Lm 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 x-axis are given by EA(x), EIz(x), ρA(x), and κsGA(x), respectively. Herein, ρ is the mass density, A is the beam cross-sectional area, E is Young’s modulus of elasticity, G 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 Iz. The vehicle is subjected to an engine thrust T at its rear end and leads to an axial force S(x) (positive in tension) in each cross section. In the present study, the structure is divided into n segments, where EA(x), EIz(x), ρA(x), and κsGA(x) distributions are uniform in each segment. Note that S(x) 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 x-axis, given by slope ∂v(x,t)/∂x of the elastic curve, can be decomposed into two parts, namely, the rotation angle θz(x,t) of the cross section about the z-axis due to pure bending and the additional angle γG(x,t) caused by shear strain, as shown in Figure 2. It shows that u(x,t) and v(x,t) stand for displacements along the coordinate axes x and y, respectively, where t is the time.Figure 2
Free body diagram of beam element of infinitesimal lengthdx in its deformed state.For Timoshenko beam, the total slope of the beam, the bending momentMz(x,t), and the shearing force Qy(x,t) are given by
(1)∂v(x,t)∂x=θz(x,t)+γG(x,t),(2)Mz(x,t)=-EIz∂θz(x,t)∂x,(3)Qy(x,t)=κsGAγG(x,t)=κsGA(∂v(x,t)∂x-θz(x,t)).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 inertiaJ∂2θz(x,t)/∂2t, which modifies the equation of moments acting on an infinitesimal beam element dx. Herein, J=ρIzdx 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γG(x,t) 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
(4)Qy(x,t)∂v(x,t)∂x-(Qy(x,t)+∂Qy(x,t)∂xdx)×(∂v(x,t)∂x+∂2v(x,t)∂x2dx)-S+(S+∂S∂xdx)=ρAdx∂2u(x,t)∂t2,(5)-Qy(x,t)+(Qy(x,t)+∂Qy(x,t)∂xdx)-S∂v(x,t)∂x+(S+∂S∂xdx)(∂v(x,t)∂x+∂2v(x,t)∂x2dx)=ρAdx∂2v(x,t)∂t2,(6)Qy(x,t)dx2+(Qy(x,t)+∂Qy(x,t)∂xdx)dx2+Mz(x,t)-(Mz(x,t)+∂Mz(x,t)∂xdx)+S∂v(x,t)∂xdx2-(S+∂S∂xdx)(∂v(x,t)∂x+∂2v(x,t)∂x2dx)dx2=ρIzdx∂2θz(x,t)∂t2.The high-order term (dx2) 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 S=EA(x)∂u/∂x. After simplifying calculations, (4)–(6) reduce to
(7)EA∂2u∂x2=ρA∂2u(x,t)∂t2,(8)∂Qy(x,t)∂x+S∂2v(x,t)∂x2=ρA∂2v(x,t)∂t2,(9)-∂Mz(x,t)∂x+Qy(x,t)=ρIz∂2θz(x,t)∂t2.Substituting the expressions for bending moment (2) and shear force (3) into the two governing equations (8) and (9) gives
(10)κsGA(∂2v(x,t)∂x2-∂θz(x,t)∂x)+S∂2v(x,t)∂x2=ρA∂2v(x,t)∂t2,(11)EIz∂2θz(x,t)∂x2+κsGA(∂v(x,t)∂x-θz(x,t))=ρIz∂2θz(x,t)∂t2.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 byv(x,t); the other mode is the transverse shearing deformation γG(x,t), as measured by the difference  (∂v(x,t)/∂x-θz(x,t)). 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 x resulting in partial derivatives ∂θz(x,t)/∂x, ∂3θz(x,t)/∂x3, and  ∂3θz(x,t)/∂x∂t2; (2) derive the expression of  ∂θz(x,t)/∂x as a function of  v(x,t) from (10); by differentiating this expression first with respect to x two times and later with respect to t two times, find expressions for ∂3θz(x,t)/∂x3 and ∂3θz(x,t)/∂x∂t2 as functions of v(x,t); (3) substitute these three partial derivatives of  θz(x,t) into the differentiated equation (11) in order to have a single differential equation in the only unknown dependent variable v(x,t). After some simple algebra, the latter reads [23] as
(12)EIz(1+SκsGA)∂4v(x,t)∂x4-(EIzρAκsGA+ρIz+SρIzκsGA)×∂4v(x,t)∂x2∂t2-S∂2v(x,t)∂x2+ρA∂2v(x,t)∂t2+ρ2IzAκsGA∂4v(x,t)∂t4=0.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:
(13)u(x,t)=u^(x)cos(ωt),v(x,t)=v^(x)cos(ωt),θz(x,t)=θ^z(x)cos(ωt),Mz(x,t)=M^z(x)cos(ωt),Qx(x,t)=Q^x(x)cos(ωt),Qy(x,t)=Q^y(x)cos(ωt).Herein,u^(x), v^(x), θ^z(x), M^z(x), Q^x(x), and Q^y(x) are the amplitudes. It should be noted that  κsGA,EIz and EA are functions of  x only; therefore
(14)∂(κsGA)∂x≡d(κsGA)dx,∂(EIz)∂x≡d(EIz)dx,∂(EA)∂x≡d(EA)dx.Substituting (13) in (7) and (12) and performing the indicated differentiations yield
(15)d2u^(x)dx2+m-ω2EAu^(x)=0,(16)(1+SκsGA)d4v^(x)dx4+[-SEIz+ρAω2κsGA+ρArg2ω2EIz(1+SκsGA)]d2v^(x)dx2,+((ρA)2rg2ω4EIzκsGA-ρAω2EIz)v^(x)=0,
where rg≡Iz/A is the radius of gyration. From the derivation, it is clear that u^(x) and v^(x), 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 ρIz (represented in rotary inertia) equal to zero in (6) and by letting the modulus G 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:
(17)d4v^(x)dx4-SEIzd2v^(x)dx2-m-ω2EIzv^(x)=0.The shear forcesQx(x,t) and Qy(x,t) at any section of the beam are related to the deflections  u(x,t) and v(x,t) and rotation angle θz(x,t) by(18a)Qx(x,t)=-EA∂u(x,t)∂x,(18b)Qy(x,t)=∂Mz(x,t)∂x+m-rg2∂2θz(x,t)∂t2+S∂v(x,t)∂x.
Using (13), (18a) and (18b) can be written as
(18c)Q^x(x)=-EA∂u^(x)∂x,(18d)Q^y(x)=∂M^z(x)∂x+m-rg2∂2θ^z(x)∂t2+S∂v^(x)∂x.In the model for a launch vehicle (Figure1), the beam will represent the segments. For jth segment of the rocket, let  x-j=x/Lm be the dimensionless length coordinate such that it takes values 0≤x-j≤l-j(l-j=lj/Lm), where l-j is the length of the segment. To express (15) and (16) in dimensionless form, the following quantities are defined:
(19)u~=u^Lm,v~=v^Lm,Dn=dndx-jn,k-2=SjLm2EIz0,R-2=(rgjLm)2,p-2=ρA0ω2Lm4EIz0,M~zj(x-)=M^zj(x)Lm(EIz0α),Q~yj(x-)=Q^yj(x)Lm2(EIz0α),Q~xj(x-)=Q^xj(x)EA0,θ~zj(x-)=θ^z(x),β=ρA0ω2Lm2EA0,α=1+(k-2κsGA¯),EIzj=EIz0EIz¯,κsGAj=(EIz0Lm2)κsGA¯,EAj=EA0EA¯,ρAj=ρA0m¯,
where EIz0, EA0, and ρA0 are reference quantities. In order to simplify the notation, EIz, EA, and ρA are consistently used for EIz¯, EA¯, and m¯, respectively. Inserting (19) into (15) and (16), one gets
(20)D2[u~(x-)]+β2m¯EA¯u~(x-)=0,(21)D4[v~(x-)]+b0D2[v~(x-)]+c0v~(x-)=0,
where
(22)b0={(-k-2EIz¯α)+m¯p-2[1κsGA¯α+R-2EIz¯]},c0=p-2EIz¯α[p-2R-2m¯2κsGA¯-m¯].


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

The characteristic equation of the flexural beam (21) has the following form:
(23)z~2+b0z~+c0=0,
where z~=^D2[v~(x-)], with solutions z~1=(-b0+Δ)/2, and z~2=(-b0-Δ)/2 where Δ=b02-4c0. It can be seen from (22) thatΔ>0 and z~2<0 for all ω. However, z~1 has two sign possibilities.(1)
z~1>0⇔Δ>b0⇒b02-4c0>b0⇒c0<0. That means that p-2<κsGA/R-2m¯ or ω2<κsGA/ρI.(2)
z~1<0⇔ω>κsGA/ρI. In the present study, only z~1>0 will be considered. Therefore, the roots of (23) are D[v~(x-)]∈{±z~1,±i-z~2}. This gives a general solution in the following form
(24)v~j(x-)=C1′jez~1x-+C2′je-z~1x-+C3′jeiz~2x-+C4′je-iz~2x-,where the constants C1′, C2′, C3′, and C4′ are to be determined from the boundary conditions of the beam. Noting that sinhx=(1/2)(ex-e-x), coshx=(1/2)(ex+e-x), sinx=(1/2i)(eix-e-ix), and cosx=(1/2)(eix+e-ix), (24) can be rewritten as
(25)v~j(x-)=C-1jcosh(λ1jx-)+C-2jsinh(λ1jx-)+C-3jcos(λ2jx-)+C-4jsin(λ2jx-),
where
(26)λ12j=|z~1|=(-b0+Δ)2;λ22j=|z~2|=(b0+Δ)2.
Substituting (3) into (1) and the result in (2) yields the variable Mz(x,t). Then substitute (8) into Mz(x,t), use (13), differentiate it with respect to x, and substitute it back into (18d). This finally gives
(27)Q^y(x)(1-ρArg2ω2κsGA)=-EIz(1+SκsGA)∂3v^(x)∂x3-(EIzρAω2κsGA+ρArg2ω2-S)∂v^(x)∂x.
The shear force can be obtained by the substitution of (19) into (27) as
(28)Q~yj(x-)=AQD3[v~(x-)]+BQD[v~(x-)],
where
(29)AQ=-EIz¯(1-p-2R-2(m¯/κsGA¯)),BQ=-(p-2(EIz¯m¯/κsGA¯)+p-2R-2m¯-k-2)α(1-p-2R-2(m¯/κsGA¯)).
Substituting the result (28) into (1) and (2) and after several manipulations, one gets
(30)θ~zj(x-)=AθD3[v~(x-)]+BθD[v~(x-)],M~zj(x-)=AMD4[v~(x-)]+BMD2[v~(x-)],
where
(31)Aθ=-ακsGA¯AQ,Bθ=1-ακsGA¯BQ,AM=EIz¯κsGA¯AQ,BM=EIz¯κsGA¯BQ-EIz¯α.
For the longitudinal vibration (20), the characteristic equation has the following form:
(32)z2+β2=0⇒z=±iβ,
where z~=^D[u(x-)]. This gives the solution
(33)u~j(x-)=C-5jcos(βjx-)+C-6jsin(βjx-).
In dimensionless form, the shear force Q~xj(x-) at section jth can be obtained by substituting (33) into (18c) and (19):
(34)Q~xj(x-)=AuD[u~(x-)],whereAu=-EA¯βj.


### 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 superscriptf, (25), and (28)–(31) can be written in matrix form as
(35)Zf(x-)=Bf(x-)C-f.
Herein, Zf(x-) is the state vector at an arbitrary position x- containing the model variables, that is, Zf(x-)={v~,θ~z,M~z,Q~y}T. The coefficient vector C-f={C-1,C-2,C-3,C-4}T summarizes unknown constants to be adopted to boundary conditions and(36)Bf(x-)=[chshcs(Aθλ13+Bθλ1)sh(Aθλ13+Bθλ1)ch(Aθλ23-Bθλ2)s(-Aθλ23+Bθλ2)c(AMλ14+BMλ12)ch(AMλ14+BMλ12)sh(AMλ24-BMλ22)c(AMλ24-BMλ22)s(AQλ13+BQλ1)sh(AQλ13+BQλ1)ch(AQλ23-BQλ2)s(-AQλ23+BQλ2)c],where ch=cosh(λ1x-),  sh=sinh(λ1x-),  c=cos(λ2x-)  and  s=sin(λ2x-). Considering a beam segment jth of length l-j, the input end 0 may be assigned to x-=0 resulting in
(37)Z0,jf=[Bf(0)]jC-jf.
Therefore, the coefficient vector can be expressed as
(38)C-jf=[Bf(0)]j-1Z0,jf.
For the beam segment output end (j+1) at  x-=l-j, (35) and (38) give
(39)Zj+1,jf=[Bf(l-)]jC-jf=[Bf(l-)]j[Bf(0)]j-1Z0,jf.
Applying the continuity equations as
(40)Zj,jf=Zj-1,jf,
this results in
(41)Zj+1,jf=UjfZ0,jf,
where
(42)Ujf=[Bf(l-)]j[Bf(0)]j-1=[v11v12v13v14v21v22v23v24v31v32v33v34v41v42v43v44]jf.
The components vlm(l,m=1,2,3,4) of the transfer matrix Ujf are listed in the appendix. For the longitudinal vibration, denoted by superscript l, the state vector is given as Zl(x-)={u~,Q~x}T, the coefficient vector C-l={C-5,C-6}T, and Bl(x-)=[cs-AusAuc], where c=cos(βx-)ands=sin(βx-). Following the pervious procedures, the longitudinal transfer matrix Ujl of the jth segment may be deduced as
(43)Ujl=[csAu-Ausc]jl=[u11u12u21u22]jl.
By combining the state vectors Zl(x-) and Zf(x-) in one column, that is, Z(x-)={u~,v~,θ~z,M~z,Q~x,Q~y}T, the transfer matrix of a beam vibrating longitudinally in the x-direction and transversely in the y-direction is
(44)U=[u11000u1200v11v12v130v140v21v22v230v240v31v32v330v34u21000u2200v41v42v430v44].
A vibrating beam comprised of n-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 n+1 connection and boundary points where Z0,1, and Zn+1,n are the boundary state vectors.Figure 3
Beam withn segments.The transfer equations of then segments, respectively, are
(45)Z2,1=U1Z0,1,Z3,2=U2Z1,2,…,Zn+1,n=UnZn-1,n.
The overall transfer equation of system is given by
(46)Zn,n+1=UnZn,n-1=UnUn-1Zn-1,n-2=⋯=UZ0,1,
where
(47)U=∏k=1nUk=UnUn-1⋯U1=[u11000u1200v11v12v130v140v21v22v230v240v31v32v330v34u21000u2200v41v42v430v44]
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 Z0,1={u~v~θ~zM~zQ~xQ~y}0,1T and right side Z10,9={u~v~θ~zM~zQ~xQ~y}10,9T to take the following form:(48a)u~10,9=u11u~0,1+u12Q~x0,1,(48b)v~10,9=v11v~0,1+v12θ~z0,1+v13M~z0,1+v14Q~y0,1,(48c)θ~z10,9=v21v~0,1+v22θ~z0,1+v23M~z0,1+v24Q~y0,1,(48d)M~z10,9=v31v~0,1+v32θ~z0,1+v33M~z0,1+v34Q~y0,1,(48e)Q~x10,9=u21u~0,1+u22Q~x0,1,(48f)Q~y10,9=v41v~0,1+v42θ~z0,1+v43M~z0,1+v44Q~y0,1.For a free-free beam, the bending moment and forces at its left and right sides are zero, that is,(49a)M~z0,1=Q~x0,1=Q~y0,1=0,(49b)M~z10,9=Q~x10,9=Q~y10,9=0,resulting in lower three equations of (48a)–(48f) which may be written in matrix form as
(50)Z10,9′=U′Z0,1′⇒{000}10,9=[0v31v32u21000v41v42]{u~v~θ~z}0,1.
Non-trivial solutions for (50) require that its coefficient determinant vanishes, that is,
(51)Δ(ω)=detU′=0⇒|0v31v32u21000v41v42|=0,
where detU′ 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 ωk (where k 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
(52)Z0,1={u~(0)=?,v~(0)=?,θ~z(0)=?,M~z(0)=0,Q~x(0)=0,Q~y(0)=0}0,1T.
In (52), the values of M~z0,1, Q~x0,1, and Q~y0,1 are equal to zero (49a), while those of u~0,1, v~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 m×n matrix U′ can be broken down into the product of three matrices—an orthogonal m×m matrix U-, an m×n diagonal matrix D-, and the transpose of an n×n orthogonal matrix V- [25]:
(53)U′m×n=U-m×mD-m×nV-n×nT,
where U-TU-=I; V-TV-=I. The columns of U- are orthogonal eigenvectors of U′U′T, the columns of V- are orthonormal eigenvectors of U′TU′, 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
(54)Z10,9′=U-D-V-TZ0,1′⋯U-D-V-TZ0,1′=0.
The solutions Z0,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.


### 2.5. Axial Force for Constant Thrust Trajectory

The axial forceS is defined as compressive force acting on the beam for the jth 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
(55)Sj(x)=TMm[Mm-∑k=1j-1ρA0m¯klk-12ρA0m¯jlj].


## 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 x-axis, given by slope ∂v(x,t)/∂x of the elastic curve, can be decomposed into two parts, namely, the rotation angle θz(x,t) of the cross section about the z-axis due to pure bending and the additional angle γG(x,t) caused by shear strain, as shown in Figure 2. It shows that u(x,t) and v(x,t) stand for displacements along the coordinate axes x and y, respectively, where t is the time.Figure 2
Free body diagram of beam element of infinitesimal lengthdx in its deformed state.For Timoshenko beam, the total slope of the beam, the bending momentMz(x,t), and the shearing force Qy(x,t) are given by
(1)∂v(x,t)∂x=θz(x,t)+γG(x,t),(2)Mz(x,t)=-EIz∂θz(x,t)∂x,(3)Qy(x,t)=κsGAγG(x,t)=κsGA(∂v(x,t)∂x-θz(x,t)).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 inertiaJ∂2θz(x,t)/∂2t, which modifies the equation of moments acting on an infinitesimal beam element dx. Herein, J=ρIzdx 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γG(x,t) 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
(4)Qy(x,t)∂v(x,t)∂x-(Qy(x,t)+∂Qy(x,t)∂xdx)×(∂v(x,t)∂x+∂2v(x,t)∂x2dx)-S+(S+∂S∂xdx)=ρAdx∂2u(x,t)∂t2,(5)-Qy(x,t)+(Qy(x,t)+∂Qy(x,t)∂xdx)-S∂v(x,t)∂x+(S+∂S∂xdx)(∂v(x,t)∂x+∂2v(x,t)∂x2dx)=ρAdx∂2v(x,t)∂t2,(6)Qy(x,t)dx2+(Qy(x,t)+∂Qy(x,t)∂xdx)dx2+Mz(x,t)-(Mz(x,t)+∂Mz(x,t)∂xdx)+S∂v(x,t)∂xdx2-(S+∂S∂xdx)(∂v(x,t)∂x+∂2v(x,t)∂x2dx)dx2=ρIzdx∂2θz(x,t)∂t2.The high-order term (dx2) 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 S=EA(x)∂u/∂x. After simplifying calculations, (4)–(6) reduce to
(7)EA∂2u∂x2=ρA∂2u(x,t)∂t2,(8)∂Qy(x,t)∂x+S∂2v(x,t)∂x2=ρA∂2v(x,t)∂t2,(9)-∂Mz(x,t)∂x+Qy(x,t)=ρIz∂2θz(x,t)∂t2.Substituting the expressions for bending moment (2) and shear force (3) into the two governing equations (8) and (9) gives
(10)κsGA(∂2v(x,t)∂x2-∂θz(x,t)∂x)+S∂2v(x,t)∂x2=ρA∂2v(x,t)∂t2,(11)EIz∂2θz(x,t)∂x2+κsGA(∂v(x,t)∂x-θz(x,t))=ρIz∂2θz(x,t)∂t2.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 byv(x,t); the other mode is the transverse shearing deformation γG(x,t), as measured by the difference  (∂v(x,t)/∂x-θz(x,t)). 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 x resulting in partial derivatives ∂θz(x,t)/∂x, ∂3θz(x,t)/∂x3, and  ∂3θz(x,t)/∂x∂t2; (2) derive the expression of  ∂θz(x,t)/∂x as a function of  v(x,t) from (10); by differentiating this expression first with respect to x two times and later with respect to t two times, find expressions for ∂3θz(x,t)/∂x3 and ∂3θz(x,t)/∂x∂t2 as functions of v(x,t); (3) substitute these three partial derivatives of  θz(x,t) into the differentiated equation (11) in order to have a single differential equation in the only unknown dependent variable v(x,t). After some simple algebra, the latter reads [23] as
(12)EIz(1+SκsGA)∂4v(x,t)∂x4-(EIzρAκsGA+ρIz+SρIzκsGA)×∂4v(x,t)∂x2∂t2-S∂2v(x,t)∂x2+ρA∂2v(x,t)∂t2+ρ2IzAκsGA∂4v(x,t)∂t4=0.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:
(13)u(x,t)=u^(x)cos(ωt),v(x,t)=v^(x)cos(ωt),θz(x,t)=θ^z(x)cos(ωt),Mz(x,t)=M^z(x)cos(ωt),Qx(x,t)=Q^x(x)cos(ωt),Qy(x,t)=Q^y(x)cos(ωt).Herein,u^(x), v^(x), θ^z(x), M^z(x), Q^x(x), and Q^y(x) are the amplitudes. It should be noted that  κsGA,EIz and EA are functions of  x only; therefore
(14)∂(κsGA)∂x≡d(κsGA)dx,∂(EIz)∂x≡d(EIz)dx,∂(EA)∂x≡d(EA)dx.Substituting (13) in (7) and (12) and performing the indicated differentiations yield
(15)d2u^(x)dx2+m-ω2EAu^(x)=0,(16)(1+SκsGA)d4v^(x)dx4+[-SEIz+ρAω2κsGA+ρArg2ω2EIz(1+SκsGA)]d2v^(x)dx2,+((ρA)2rg2ω4EIzκsGA-ρAω2EIz)v^(x)=0,
where rg≡Iz/A is the radius of gyration. From the derivation, it is clear that u^(x) and v^(x), 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 ρIz (represented in rotary inertia) equal to zero in (6) and by letting the modulus G 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:
(17)d4v^(x)dx4-SEIzd2v^(x)dx2-m-ω2EIzv^(x)=0.The shear forcesQx(x,t) and Qy(x,t) at any section of the beam are related to the deflections  u(x,t) and v(x,t) and rotation angle θz(x,t) by(18a)Qx(x,t)=-EA∂u(x,t)∂x,(18b)Qy(x,t)=∂Mz(x,t)∂x+m-rg2∂2θz(x,t)∂t2+S∂v(x,t)∂x.
Using (13), (18a) and (18b) can be written as
(18c)Q^x(x)=-EA∂u^(x)∂x,(18d)Q^y(x)=∂M^z(x)∂x+m-rg2∂2θ^z(x)∂t2+S∂v^(x)∂x.In the model for a launch vehicle (Figure1), the beam will represent the segments. For jth segment of the rocket, let  x-j=x/Lm be the dimensionless length coordinate such that it takes values 0≤x-j≤l-j(l-j=lj/Lm), where l-j is the length of the segment. To express (15) and (16) in dimensionless form, the following quantities are defined:
(19)u~=u^Lm,v~=v^Lm,Dn=dndx-jn,k-2=SjLm2EIz0,R-2=(rgjLm)2,p-2=ρA0ω2Lm4EIz0,M~zj(x-)=M^zj(x)Lm(EIz0α),Q~yj(x-)=Q^yj(x)Lm2(EIz0α),Q~xj(x-)=Q^xj(x)EA0,θ~zj(x-)=θ^z(x),β=ρA0ω2Lm2EA0,α=1+(k-2κsGA¯),EIzj=EIz0EIz¯,κsGAj=(EIz0Lm2)κsGA¯,EAj=EA0EA¯,ρAj=ρA0m¯,
where EIz0, EA0, and ρA0 are reference quantities. In order to simplify the notation, EIz, EA, and ρA are consistently used for EIz¯, EA¯, and m¯, respectively. Inserting (19) into (15) and (16), one gets
(20)D2[u~(x-)]+β2m¯EA¯u~(x-)=0,(21)D4[v~(x-)]+b0D2[v~(x-)]+c0v~(x-)=0,
where
(22)b0={(-k-2EIz¯α)+m¯p-2[1κsGA¯α+R-2EIz¯]},c0=p-2EIz¯α[p-2R-2m¯2κsGA¯-m¯].


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

The characteristic equation of the flexural beam (21) has the following form:
(23)z~2+b0z~+c0=0,
where z~=^D2[v~(x-)], with solutions z~1=(-b0+Δ)/2, and z~2=(-b0-Δ)/2 where Δ=b02-4c0. It can be seen from (22) thatΔ>0 and z~2<0 for all ω. However, z~1 has two sign possibilities.(1)
z~1>0⇔Δ>b0⇒b02-4c0>b0⇒c0<0. That means that p-2<κsGA/R-2m¯ or ω2<κsGA/ρI.(2)
z~1<0⇔ω>κsGA/ρI. In the present study, only z~1>0 will be considered. Therefore, the roots of (23) are D[v~(x-)]∈{±z~1,±i-z~2}. This gives a general solution in the following form
(24)v~j(x-)=C1′jez~1x-+C2′je-z~1x-+C3′jeiz~2x-+C4′je-iz~2x-,where the constants C1′, C2′, C3′, and C4′ are to be determined from the boundary conditions of the beam. Noting that sinhx=(1/2)(ex-e-x), coshx=(1/2)(ex+e-x), sinx=(1/2i)(eix-e-ix), and cosx=(1/2)(eix+e-ix), (24) can be rewritten as
(25)v~j(x-)=C-1jcosh(λ1jx-)+C-2jsinh(λ1jx-)+C-3jcos(λ2jx-)+C-4jsin(λ2jx-),
where
(26)λ12j=|z~1|=(-b0+Δ)2;λ22j=|z~2|=(b0+Δ)2.
Substituting (3) into (1) and the result in (2) yields the variable Mz(x,t). Then substitute (8) into Mz(x,t), use (13), differentiate it with respect to x, and substitute it back into (18d). This finally gives
(27)Q^y(x)(1-ρArg2ω2κsGA)=-EIz(1+SκsGA)∂3v^(x)∂x3-(EIzρAω2κsGA+ρArg2ω2-S)∂v^(x)∂x.
The shear force can be obtained by the substitution of (19) into (27) as
(28)Q~yj(x-)=AQD3[v~(x-)]+BQD[v~(x-)],
where
(29)AQ=-EIz¯(1-p-2R-2(m¯/κsGA¯)),BQ=-(p-2(EIz¯m¯/κsGA¯)+p-2R-2m¯-k-2)α(1-p-2R-2(m¯/κsGA¯)).
Substituting the result (28) into (1) and (2) and after several manipulations, one gets
(30)θ~zj(x-)=AθD3[v~(x-)]+BθD[v~(x-)],M~zj(x-)=AMD4[v~(x-)]+BMD2[v~(x-)],
where
(31)Aθ=-ακsGA¯AQ,Bθ=1-ακsGA¯BQ,AM=EIz¯κsGA¯AQ,BM=EIz¯κsGA¯BQ-EIz¯α.
For the longitudinal vibration (20), the characteristic equation has the following form:
(32)z2+β2=0⇒z=±iβ,
where z~=^D[u(x-)]. This gives the solution
(33)u~j(x-)=C-5jcos(βjx-)+C-6jsin(βjx-).
In dimensionless form, the shear force Q~xj(x-) at section jth can be obtained by substituting (33) into (18c) and (19):
(34)Q~xj(x-)=AuD[u~(x-)],whereAu=-EA¯βj.


## 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 superscriptf, (25), and (28)–(31) can be written in matrix form as
(35)Zf(x-)=Bf(x-)C-f.
Herein, Zf(x-) is the state vector at an arbitrary position x- containing the model variables, that is, Zf(x-)={v~,θ~z,M~z,Q~y}T. The coefficient vector C-f={C-1,C-2,C-3,C-4}T summarizes unknown constants to be adopted to boundary conditions and(36)Bf(x-)=[chshcs(Aθλ13+Bθλ1)sh(Aθλ13+Bθλ1)ch(Aθλ23-Bθλ2)s(-Aθλ23+Bθλ2)c(AMλ14+BMλ12)ch(AMλ14+BMλ12)sh(AMλ24-BMλ22)c(AMλ24-BMλ22)s(AQλ13+BQλ1)sh(AQλ13+BQλ1)ch(AQλ23-BQλ2)s(-AQλ23+BQλ2)c],where ch=cosh(λ1x-),  sh=sinh(λ1x-),  c=cos(λ2x-)  and  s=sin(λ2x-). Considering a beam segment jth of length l-j, the input end 0 may be assigned to x-=0 resulting in
(37)Z0,jf=[Bf(0)]jC-jf.
Therefore, the coefficient vector can be expressed as
(38)C-jf=[Bf(0)]j-1Z0,jf.
For the beam segment output end (j+1) at  x-=l-j, (35) and (38) give
(39)Zj+1,jf=[Bf(l-)]jC-jf=[Bf(l-)]j[Bf(0)]j-1Z0,jf.
Applying the continuity equations as
(40)Zj,jf=Zj-1,jf,
this results in
(41)Zj+1,jf=UjfZ0,jf,
where
(42)Ujf=[Bf(l-)]j[Bf(0)]j-1=[v11v12v13v14v21v22v23v24v31v32v33v34v41v42v43v44]jf.
The components vlm(l,m=1,2,3,4) of the transfer matrix Ujf are listed in the appendix. For the longitudinal vibration, denoted by superscript l, the state vector is given as Zl(x-)={u~,Q~x}T, the coefficient vector C-l={C-5,C-6}T, and Bl(x-)=[cs-AusAuc], where c=cos(βx-)ands=sin(βx-). Following the pervious procedures, the longitudinal transfer matrix Ujl of the jth segment may be deduced as
(43)Ujl=[csAu-Ausc]jl=[u11u12u21u22]jl.
By combining the state vectors Zl(x-) and Zf(x-) in one column, that is, Z(x-)={u~,v~,θ~z,M~z,Q~x,Q~y}T, the transfer matrix of a beam vibrating longitudinally in the x-direction and transversely in the y-direction is
(44)U=[u11000u1200v11v12v130v140v21v22v230v240v31v32v330v34u21000u2200v41v42v430v44].
A vibrating beam comprised of n-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 n+1 connection and boundary points where Z0,1, and Zn+1,n are the boundary state vectors.Figure 3
Beam withn segments.The transfer equations of then segments, respectively, are
(45)Z2,1=U1Z0,1,Z3,2=U2Z1,2,…,Zn+1,n=UnZn-1,n.
The overall transfer equation of system is given by
(46)Zn,n+1=UnZn,n-1=UnUn-1Zn-1,n-2=⋯=UZ0,1,
where
(47)U=∏k=1nUk=UnUn-1⋯U1=[u11000u1200v11v12v130v140v21v22v230v240v31v32v330v34u21000u2200v41v42v430v44]
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 Z0,1={u~v~θ~zM~zQ~xQ~y}0,1T and right side Z10,9={u~v~θ~zM~zQ~xQ~y}10,9T to take the following form:(48a)u~10,9=u11u~0,1+u12Q~x0,1,(48b)v~10,9=v11v~0,1+v12θ~z0,1+v13M~z0,1+v14Q~y0,1,(48c)θ~z10,9=v21v~0,1+v22θ~z0,1+v23M~z0,1+v24Q~y0,1,(48d)M~z10,9=v31v~0,1+v32θ~z0,1+v33M~z0,1+v34Q~y0,1,(48e)Q~x10,9=u21u~0,1+u22Q~x0,1,(48f)Q~y10,9=v41v~0,1+v42θ~z0,1+v43M~z0,1+v44Q~y0,1.For a free-free beam, the bending moment and forces at its left and right sides are zero, that is,(49a)M~z0,1=Q~x0,1=Q~y0,1=0,(49b)M~z10,9=Q~x10,9=Q~y10,9=0,resulting in lower three equations of (48a)–(48f) which may be written in matrix form as
(50)Z10,9′=U′Z0,1′⇒{000}10,9=[0v31v32u21000v41v42]{u~v~θ~z}0,1.
Non-trivial solutions for (50) require that its coefficient determinant vanishes, that is,
(51)Δ(ω)=detU′=0⇒|0v31v32u21000v41v42|=0,
where detU′ 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 ωk (where k 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
(52)Z0,1={u~(0)=?,v~(0)=?,θ~z(0)=?,M~z(0)=0,Q~x(0)=0,Q~y(0)=0}0,1T.
In (52), the values of M~z0,1, Q~x0,1, and Q~y0,1 are equal to zero (49a), while those of u~0,1, v~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 m×n matrix U′ can be broken down into the product of three matrices—an orthogonal m×m matrix U-, an m×n diagonal matrix D-, and the transpose of an n×n orthogonal matrix V- [25]:
(53)U′m×n=U-m×mD-m×nV-n×nT,
where U-TU-=I; V-TV-=I. The columns of U- are orthogonal eigenvectors of U′U′T, the columns of V- are orthonormal eigenvectors of U′TU′, 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
(54)Z10,9′=U-D-V-TZ0,1′⋯U-D-V-TZ0,1′=0.
The solutions Z0,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.


## 2.5. Axial Force for Constant Thrust Trajectory

The axial forceS is defined as compressive force acting on the beam for the jth 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
(55)Sj(x)=TMm[Mm-∑k=1j-1ρA0m¯klk-12ρA0m¯jlj].


## 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, γ(≡T/Tcr≡TLm2/π2EIz) is the ratio of axial force to the buckling load γ(≡T/Tcr≡TLm2/π2EIz). A negative sign in γ indicates a compressive force.Table 1
Frequency parameterp-2 for a free-free Euler-Bernoulli square beam under axial load for EIz=200×109 N/m2; G=1×1090 N/m2; ρ=7800 kg/m3; b=h=0.01 m; Lm=1 m.
Axial loadγ
Mode
p
-
2 [24]
p
-
2 [TMM]
−0.999
1
0.84680644
0.84680659
2
7.22386030
7.223860352
3
10.6308452
10.630845189
0.0
1
0.00000000
0.000000000
2
4.73004074
4.7300407450
3
7.85320462
7.853204624
4
10.9956078
10.995607838
1.0
1
3.26787488
3.267874878
2
5.59376731
5.593767306
3
8.35398454
8.353984538
4
11.3264204
11.326420425


### 3.2. Free-Free Uniform Timoshenko Beam

A good agreement can be also observed from Table2 between measured [26], ANSYS [26], and the TMM values for a single segment. The geometry and material data are as follows:
(56)EIz=72.66611×109N/m2,G=27.17481×109N/m2,ρ=2817.9111kg/m3,b=0.01009m,h=0.01004m,Lm=0.04004m,ks=0.869877423.
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.
Mode
Measured [26]
ANSYS [26]
TMM
1
27359.6
27417.6 (+0.212%)
27407.233 (+0.174%)
2
60862.0
60882.0 (+0.0328%)
60851.321 (+0.0175%)
3
97609.5
97734.4 (+0.128%)
97796.367 (+0.191%)
4
131494
131658.0 (+0.124%)
132277.502 (+0.595%)


### 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 EIz and mass per unit length ρA distributions. Table 3 illustrates a comparison of TMM using different equal segments (i.e., uniform EIz and ρA 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.
TMM
Flight test [27] time from launch(0.25 to 0.75 sec)
Computed data [27] time from launch(0.0 to 1.0 sec)
Mode
Number of segments
50
100
400
1000
1
12.95
37.18
10.58
10.58
≈10.5
10.5 to 10.9
2
11.55
35.75
35.03
35.02
—
—
3
10.58
66.25
62.41
62.54
—
—Variation of (a) flexural stiffness and (b) mass distribution, for a two-stage launch vehicle prior the launch [27].
(a)
(b)


### 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 (L/D≈25, for whole launch vehicle) and low aspect ratio (L/D≈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 EIz, mass per unit length ρA, shear rigidity ksGA, and radius of gyration rg. 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 L/D≈25 only, while it has a significant effect on the data for L/D≈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.
Mode
Timoshenko beam
Euler-Bernoulli beam
With rotary inertia and shear deformation
With rotary inertia only
With shear deformation only
(without rotary inertia and shear deformation)
Reference [5]
TMM
Reference [5]
TMM
Reference [5]
TMM
Reference [5]
TMM
Whole multistage launch vehicle (L/D≈25)
1
16.017
16.005
16.054
16.042
16.083
16.069
16.120
16.107
2
46.521
46.492
47.195
47.172
46.770
46.735
47.457
47.428
3
93.943
93.710
97.566
97.393
94.727
94.471
98.686
98.322
Upper stage of the launch vehicle (L/D≈5)
1
334.776
335.019
400.968
401.468
348.945
349.146
438.557
439.003
2
760.445
760.596
1104.970
1106.305
803.663
803.364
1367.317
1369.986
3
960.819
959.593
1655.463
1659.978
1049.897
1048.378
2292.381
2290.001Parameters of a multistage launch vehicle [5]: (a) flexural stiffness, (b) shear rigidity, (c) mass, and (d) radius of gyration.
(a)
(b)
(c)
(d)Effects of shear deformation and rotary inertia on (a) mode shape, (b) rotation angle, (c) moment, and (d) shear of first mode forL/D≈25: ♦ Timoshenko model; —Euler-Bernoulli model.
(a)
(b)
(c)
(d)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 (L/D≈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 (L/D≈5): - - - Timoshenko model; —Euler-Bernoulli model.


## 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, γ(≡T/Tcr≡TLm2/π2EIz) is the ratio of axial force to the buckling load γ(≡T/Tcr≡TLm2/π2EIz). A negative sign in γ indicates a compressive force.Table 1
Frequency parameterp-2 for a free-free Euler-Bernoulli square beam under axial load for EIz=200×109 N/m2; G=1×1090 N/m2; ρ=7800 kg/m3; b=h=0.01 m; Lm=1 m.
Axial loadγ
Mode
p
-
2 [24]
p
-
2 [TMM]
−0.999
1
0.84680644
0.84680659
2
7.22386030
7.223860352
3
10.6308452
10.630845189
0.0
1
0.00000000
0.000000000
2
4.73004074
4.7300407450
3
7.85320462
7.853204624
4
10.9956078
10.995607838
1.0
1
3.26787488
3.267874878
2
5.59376731
5.593767306
3
8.35398454
8.353984538
4
11.3264204
11.326420425


## 3.2. Free-Free Uniform Timoshenko Beam

A good agreement can be also observed from Table2 between measured [26], ANSYS [26], and the TMM values for a single segment. The geometry and material data are as follows:
(56)EIz=72.66611×109N/m2,G=27.17481×109N/m2,ρ=2817.9111kg/m3,b=0.01009m,h=0.01004m,Lm=0.04004m,ks=0.869877423.
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.
Mode
Measured [26]
ANSYS [26]
TMM
1
27359.6
27417.6 (+0.212%)
27407.233 (+0.174%)
2
60862.0
60882.0 (+0.0328%)
60851.321 (+0.0175%)
3
97609.5
97734.4 (+0.128%)
97796.367 (+0.191%)
4
131494
131658.0 (+0.124%)
132277.502 (+0.595%)


## 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 EIz and mass per unit length ρA distributions. Table 3 illustrates a comparison of TMM using different equal segments (i.e., uniform EIz and ρA 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.
TMM
Flight test [27] time from launch(0.25 to 0.75 sec)
Computed data [27] time from launch(0.0 to 1.0 sec)
Mode
Number of segments
50
100
400
1000
1
12.95
37.18
10.58
10.58
≈10.5
10.5 to 10.9
2
11.55
35.75
35.03
35.02
—
—
3
10.58
66.25
62.41
62.54
—
—Variation of (a) flexural stiffness and (b) mass distribution, for a two-stage launch vehicle prior the launch [27].
(a)
(b)


## 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 (L/D≈25, for whole launch vehicle) and low aspect ratio (L/D≈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 EIz, mass per unit length ρA, shear rigidity ksGA, and radius of gyration rg. 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 L/D≈25 only, while it has a significant effect on the data for L/D≈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.
Mode
Timoshenko beam
Euler-Bernoulli beam
With rotary inertia and shear deformation
With rotary inertia only
With shear deformation only
(without rotary inertia and shear deformation)
Reference [5]
TMM
Reference [5]
TMM
Reference [5]
TMM
Reference [5]
TMM
Whole multistage launch vehicle (L/D≈25)
1
16.017
16.005
16.054
16.042
16.083
16.069
16.120
16.107
2
46.521
46.492
47.195
47.172
46.770
46.735
47.457
47.428
3
93.943
93.710
97.566
97.393
94.727
94.471
98.686
98.322
Upper stage of the launch vehicle (L/D≈5)
1
334.776
335.019
400.968
401.468
348.945
349.146
438.557
439.003
2
760.445
760.596
1104.970
1106.305
803.663
803.364
1367.317
1369.986
3
960.819
959.593
1655.463
1659.978
1049.897
1048.378
2292.381
2290.001Parameters of a multistage launch vehicle [5]: (a) flexural stiffness, (b) shear rigidity, (c) mass, and (d) radius of gyration.
(a)
(b)
(c)
(d)Effects of shear deformation and rotary inertia on (a) mode shape, (b) rotation angle, (c) moment, and (d) shear of first mode forL/D≈25: ♦ Timoshenko model; —Euler-Bernoulli model.
(a)
(b)
(c)
(d)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 (L/D≈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 (L/D≈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 Figures10(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 (t=0) and burn-out (t=50sec) times. The axial compressive force is considered for a nominal value of thrust as shown in Figure 10(b).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 (—t=0 and - - - t=50sec).
(a)
(b)
(c)Figure11 shows the effect of increased operating thrust on the first natural frequency at t=0 (maximum mass) and t=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 (t=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 (t=0) but lower allowable thrust range. This fact can be shown readily from Figure 11. For t=50sec, the critical buckling thrust (associated with zero frequency) is lower than that at t=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 (—t=0 and - - - t=50sec).Figure12 illustrates the effect of maximum operating thrust (in the vicinity of critical buckling thrust taken from Figure 11), that is, T=1.5×106N|t=0 and T=1.2×106N|t=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 Figure13).Figure 13
Closed structure system.



---

*Source: 764673-2013-03-14.xml*