## Abstract

Stepped beams constitute an important class of engineering structures whose vibration response has been widely studied. Many of the existing methods for studying stepped beams manifest serious numerical difficulties as the number of segments or the frequency of excitation increase. In this article, we focus on the transfer matrix method (TMM), which provides a simple and elegant formulation for multistep beams. The main idea in the TMM is to model each step in the beam as a uniform element whose vibration configurations are spanned by the segment’s local eigenfunctions. Utilizing these local expressions, the boundary conditions at the ends of the multistep beam as well as the continuity and compatibility conditions across each step are used to obtain the nonlinear eigenvalue problem. Also, and perhaps more importantly, we provide a reformulation for multistep Euler–Bernoulli beams that avoids much of the numerical singularity problems that have plagued most of the earlier efforts. When this reformulation is extended to multisegment Timoshenko beams, the numerical difficulties appear to be mitigated, but not solved.

## 1 Introduction

Stepped beams constitute an important class of engineering structures whose vibration response has been widely studied. Some of the methods for analyzing these beams include the transfer matrix method (TMM) [1,2], finite difference method [3], finite element method [4], Green’s function method [5], and composite element method [6]. In this article, we focus on the TMM that provides a simple and elegant formulation for multistep beams. Also, and perhaps more importantly, we provide a reformulation for multistep beams that avoid the numerical singularity problems that have plagued most earlier efforts. Further, this reformulation is extended to accommodate short beam segments better represented by a Timoshenko–Ehrenfest beam model [7].

The direct approach to modeling a beam with N steps as N connected beam segments is often called the numerical assembly method, or the full matrix approach, and it involves writing a system of 4N equations that incorporate the boundary conditions as well as the continuity and compatibility conditions. The continuity conditions ensure that the transverse vibration and the angle of rotation of the cross sections of the beam segments are equal across every step, while the compatibility conditions set the moment and shear forces equal at each interface. However, one of the drawbacks of the direct approach is that the size of the matrix M in the resulting nonlinear eigenvalue problem M(λ)x = 0 grows with the number of segments. Although the matrix M remains sparse, there will be an overall increase in the computational time of the corresponding characteristic equation since the computational complexity for determinants is roughly $O(N3)$ [8]. One approach to prevent the increase in the size of M is to utilize the TMM.

The main idea in the TMM is to model each step in the beam as a uniform element whose vibration configurations are spanned by the segment’s local eigenfunctions. Utilizing these local expressions, the boundary conditions at the ends of the multistep beam as well as the continuity and compatability conditions across each step are used to obtain the nonlinear eigenvalue problem, see Sec. 5 for more details. One of the advantages of this formulation is that the size of the resulting matrix M in the nonlinear eigenvalue problem remains constant regardless of the number of the steps in the beam. In the early formulation of the TMM, the matrix M is a square 4 × 4 matrix of real entries [1,2], but the size of this matrix can be further reduced to 2 × 2 for certain boundary conditions. Specifically, Koplow et al. [9] provided a closed form solution for a discontinuous Euler–Bernoulli (EB) beam with one step for the case of the free-free boundary condition which was studied for a concentrated harmonic load at the boundaries. The solution of Koplow et al. for a single step EB beam was later extended for arbitrary dynamic loads in Ref. [10], while a multistep version for EB was described in Ref. [11].

While many studies on TMM present formulations that theoretically accommodate beams with an arbitrary number of steps, in practice it is only possible to analyze beams with a limited number of steps. There are two reasons for this limitation, one is related to the numerical instability of the prevalent eigenmode expressions, and the second is a consequence of the repeated matrix multiplications in the TMM. Specifically, many of the tools that leverage EB or Timoshenko–Ehrenfest (TB) theory with multistep beams utilize the conventional expressions for the modes which suffer from an inherent numerical instability due to the presence of hyperbolic trigonometric terms [12]. These terms are unbounded as the frequency gets higher or the mode shape approaches the far end of the beam, and the resulting round-off errors lead to numerical instabilities that preclude the analysis for high frequencies. This problem is exacerbated when beam theory is utilized to analyze multistep beams with TMM since the numerical instabilities occur sooner because the determinant in the resulting matrix includes products of the hyperbolic terms. While some prior works have attempted to mitigate this issue [13], this problem remains unsolved even for beams with a small number of segments. In fact, the issue of numerical instabilities in the eigenvalue problem even for a single-segment beam was only recently resolved for the Euler–Bernoulli beam [14,15] and for the Timoshenko–Ehrenfest beam [15]. In the following development, an effort was made to keep the notation of Ref. [15]. The approach developed in these two studies requires that the growth rate of the hyperbolic terms must exactly cancel to obtain a physically feasible, bounded solution. Consequently, Refs. [14,15] rewrite the traditional, numerically unstable expressions so that they are numerically tractable. Unfortunately, the characteristic equations for multistep beams are affected by the number of steps and by each segment’s mechanical properties; therefore, the formulation in Refs. [14,15] cannot be employed directly to stabilize the multistep beam problem. While some of the published works utilize the finite element formulation [16], obtaining analytic expressions is useful for gaining an understanding of the characteristics of the solution, particularly asymptotic behavior. In addition, analytic expressions (when they can be evaluated) are necessary to gain confidence in numerical methods such as finite elements (mesh convergence, code validation, etc.).

In this article, we reformulate the solutions to the eigenvalue problem starting with basis functions introduced by Ref. [17] for Euler–Bernoulli beams and generalized here for Timoshenko–Ehrenfest beam models as well. The resulting novel expressions are then combined with the transfer matrix method to obtain the eigenvalue problem for multistep beams. This approach enables the computation of the eigenvalues and the eigenmodes for beams with long spans and large number of segments. The resulting eigenvalues are verified against finite element computations. Our results show that it is possible to reliably obtain the eigenmodes of multistep Euler–Bernoulli beams. For the Timoshenko–Ehrenfest beam, while our approach enables the computation of the eigenmodes well beyond what current analytical methods allow, it does have some limitations. Specifically, because the eigenmode expressions for the Timoshenko–Ehrenfest beam explicitly contain the square of the frequency, the resulting expressions for multistep beams are ill conditioned, and numerical errors grow with the number of beam segments.

### 1.1 Limitations in Current Analytical Formulations for Multistep Beams.

There are two limitations to current analytical formulations for multistep beams: (1) the analytical expressions for each beam segment contain hyperbolic terms, which make the characteristic function ill conditioned, and (2) obtaining the characteristic function either using the direct method or the TMM involves taking the determinant of a matrix with hyperbolic terms. This speeds up the occurrence of numerical instabilities even at a low number of segments.

To show the source of the numerical difficulties consider the mode shape expressions for the Euler–Bernoulli (Eq. (1)) and the Timoshenko–Ehrenfest beam models (the case λ < α in Table 1). It is the sinh (μx/L) and cosh (μx/L) that cause problems as μ grows large and x/L → 1.

Table 1

The conventional mode shapes’ expressions for the Timoshenko–Ehrenfest beam model [18]

 λ < α $[U(x)Φ(x)]=A1[sinhμxLλ+μ2μcoshμxL]+A2[coshμxLλ+μ2μsinhμxL]+A3[sinωxL−λ−ω2ωcosωxL]+A4[cosωxLλ−ω2ωsinωxL]$⁠. λ = α $[U(x)Φ(x)]=A1[01]+A2[1αxL]+A3[sinωxL−α−ω2ωcosωxL]+A4[cosωxLα−ω2ωsinωxL]$ λ > α $[U(x)Φ(x)]=A1[sinθxL−λ−θ2θcosθxL]+A2[cosθxLλ−θ2θsinθxL]+A3[sinωxL−λ−ω2ωcosωxL]+A4[cosωxLλ−ω2ωsinωxL]$
 λ < α $[U(x)Φ(x)]=A1[sinhμxLλ+μ2μcoshμxL]+A2[coshμxLλ+μ2μsinhμxL]+A3[sinωxL−λ−ω2ωcosωxL]+A4[cosωxLλ−ω2ωsinωxL]$⁠. λ = α $[U(x)Φ(x)]=A1[01]+A2[1αxL]+A3[sinωxL−α−ω2ωcosωxL]+A4[cosωxLα−ω2ωsinωxL]$ λ > α $[U(x)Φ(x)]=A1[sinθxL−λ−θ2θcosθxL]+A2[cosθxLλ−θ2θsinθxL]+A3[sinωxL−λ−ω2ωcosωxL]+A4[cosωxLλ−ω2ωsinωxL]$

Note: Here, α = AL2/I, where A is the cross-sectional area, L is the length of the beam segment, and I is the area moment of inertia of the cross section about the neutral axis, $λ=TT2ω^2$, where $ω^$ is the beam’s natural frequency, $TT=Lρ/Gκ$. See Eq. (3a) for the expressions for ω, μ, and θ.

If we were inclined to use expressions like those of Eq. (1) for Euler–Bernoulli beams or of Table 1 for Timoshenko–Ehrenfest beam in Galerkin approximation for each section of a segmented beam, we would again find ourselves having to deal with the associated numerical difficulties where some of the component terms grow unboundedly with frequency, even though we know some are finite.

## 2 Beam Models

This section describes the conventional mode shape expressions for the Euler–Bernoulli (Sec. 2.1) and the Timoshenko beam models (Sec. 2.2), and introduces some useful notation that we will use in the subsequent presentation.

### 2.1 Euler–Bernoulli Beam.

The general form for the deformations of an Euler–Bernoulli beam in free vibration is
$u(x)=A1sinhμxL+A2coshμxL+A3sinμxL+A4cosμxL$
(1)
where A1A4 are functions of μ, x is the axial position along the beam’s span, L is the length of the beam, μ = (αλ)1/4, α = A L2/I, $λ=TEω^$, $TE=Lρ/E$, $ω^$ is the physical natural frequency of the beam, A is the cross section area, I is the moment of inertia of the cross section, ρ is the mass density, and E is the modulus of elasticity.
Since we are interested in multistep beams, we use Eq. (1) and the local coordinates shown in Fig. 1 to define the matrix $MkEB(xk,ω^),$ which corresponds to the kth beam segment such that
$MkEB(xk,ω^)Ak=[uk(xk)uk′(xk)EkIkuk″(xk)EkIkuk‴(xk)]$
(2)
for all xk ∈ (0, Lk), where $AkT=[A1k,A2k,A3k,A4k]$. The rows of matrix $MkEB(xk,ω^)$ correspond to the displacement, slope, bending moment, and shear force at location xk of beam segment k, respectively.
Fig. 1
Fig. 1
Close modal

### 2.2 Timoshenko–Ehrenfest Beam.

The general form for the deformation of a uniform Timoshenko–Ehrenfest beam in free vibration depends on the ratio λ/α as shown in Table 1. Define the dimensionless constant γ = 2/E. The real quantities ω, μ (for the case λ < α), and θ (for the case λ > α) are given in terms of λ as Ref. [18]
$ω2=12λ(1+γ)(Δ1/2+1)$
(3a)
$μ2=12λ(1+γ)(Δ1/2−1)$
(3b)
$θ2=12λ(1+γ)(1−Δ1/2)$
(3c)
$whereΔ=1−4γ(1+γ)2(1−αλ)$
(3d)
In this case, μ and ω are functions of the dimensionless frequency $λ$, which is related to the physical frequency $ω^$ by $λ=ω^TT$, where $TT=Lρ/(Gκ)$, G is the shear modulus, and κ is the shear correction coefficient.
Similar to Eq. (2), we can define a matrix $MkTB(xk,ω^)$ whose rows contain the displacement, slope, bending moment, and shear force of the Timoshenko–Ehrenfest beam, respectively, according to
$MkTB(xk,ω^)Ak=[Uk(xk)Φ(xk)EkIkΦk′(xk)AkGkκk(Uk′(xk)−Φ(xk))]$
(4)
where $AkT=[A1k,A2k,A3k,A4k]$.

## 3 Reformulation of the Eigenmode Expressions

The enabling idea to reformulating the eigenmode expressions comes from the observation that the term A1sinh μ(x/L) + A2cosh μ(x/L) in Eq. (1) and both the previous term and the term A1cosh μ(x/L) + A2sinh μ(x/L) in Table 1 must be bounded for the displacement to be physically meaningful. This implicitly puts some constraint on the relationship between the coefficients A1 and A2.

The assertion that these problematic terms are bounded can be written as follows:
$A1(μ)sinh(μxL)+A2(μ)cosh(μxL)=O(1)$
(5)
from which we write
$A1(μ)2(eμ(x/L)−e−μ(x/L))+A2(μ)2(eμ(x/L)+e−μ(x/L))=A1(μ)+A2(μ)2eμ(x/L)+−A1(μ)+A2(μ)2e−μ(x/L)=O(1)$
(6)
where the big O notation $O(1)$ means that the absolute value of the left-hand side of the equation is at most some constant. We see that it is only the coefficients of eμ(x/L) that are relevant to the boundedness issue:
$A1(μ)+A2(μ)2eμ(x/L)=O(1)$
(7)
This permits to replace the big O in the aforementioned equation with P(μ), where P(μ) is a bounded function of μ. We can now express A2 in terms of A1 and P(μ) according to
$A2(μ)=−A1(1+P(μ)e−μ)$
(8)
This formulation is used in Sec. 3.1 to derive alternative expressions for the eignemodes of Euler–Bernoulli beams, while Sec. 3.2 shows the corresponding derivation for Timoshenko beams.

### 3.1 The Euler–Bernoulli Beam.

We express Eq. (1) as follows:
$u(x)=−A1e−μ(x/L)−(A1P(μ)/2)(e−μ(1−(x/L))+e−μ(1+(x/L)))+A3sinμxL+A4cosμxL$
(9)
and we define $A^2=−A1P(μ)/2$ and absorb the minus sign in the first term into $A^1$, so
$u(x)=A^1e−μ(x/L)+A^2(e−μ(1−(x/L))+e−μ(1+(x/L)))+A3sinμxL+A4cosμxL$
(10)
which can take on only bounded values.

### 3.2 The Timoshenko Beam.

We may now rewrite the first equation in Table 1 as follows:
$[U(x)Φ(x)]=−A1[e−μ(x/L)−λ+μ2μe−μ(x/L)]−A1P(μ)/2[(e−μ(1−(x/L))+e−μ(1+(x/L)))λ+μ2μ(e−μ(1−(x/L))−e−μ(1+(x/L)))]+A3[sinωxL−λ−ω2ωcosωxL]+A4[cosωxLλ−ω2ωsinωxL]$
(11)
Since P(μ) is yet unspecified other than being a bounded function of μ, we may just replace the term −A1P(μ)/2 with the unknown coefficient $A~2$. Further, we also replace −A1 with $A~1$ to obtain
$[U(x)Φ(x)]=A~1[e−μ(x/L)−λ+μ2μe−μ(x/L)]+A~2[(e−μ(1−(x/L))+e−μ(1+(x/L)))λ+μ2μ(e−μ(1−(x/L))−e−μ(1+(x/L)))]+A3[sinωxL−λ−ω2ωcosωxL]+A4[cosωxLλ−ω2ωsinωxL]$

## 4 Full Matrix Approach

Although the full matrix approach is not as efficient at the transfer matrix solution approach detailed here, it is worth presenting for the purpose of comparison.

The continuity across the interface between segment k and segment k + 1 requires
$Mk(Lk,ω^)Ak=Mk+1(0,ω^)Ak+1$
(12)
where M is the transfer matrix of either Euler–Bernouli (EB) or Timoshenko–Ehrenfest (TB) beams in Eqs. (2) and (4), respectively.

To capture the two homogeneous boundary conditions at the external boundaries of the first beam, we define the 2 × 4 matrix BF consisting of the rows of $M1(0,ω^)$ associated to those two degrees-of-freedom. To capture the two homogeneous boundary conditions at the external boundaries of the last beam, we define the 2 × 4 matrix BL consisting of the rows of $ML(LL,ω^)$ associated to those two degrees-of-freedom.

This imposition of the boundary conditions and continuity across interfaces is illustrated for the case of three segments in Eq. (13). For such a full matrix to admit a solution, the matrix must be singular and code was written to assemble these matrices and evaluate their determinants over a fine distribution of frequencies over the range of interest. In the regions of zero crossing where the sign of the determinant changes, yet finer local meshes were created, the determinants were evaluated over those finer meshes, those values were fit with cubic splines, and zeros of those splines were then found to obtain fine resolution of zeros.
$[BF(ω^)⋅0⋅0M1(L1,ω^)M2(0,ω^)⋅0⋅0M2(L2,ω^)M3(0,ω^)BL(ω^)]{⋅0A1⋅0A2⋅0A3}={012}$
(13)

## 5 Transfer Matrix Method

We now develop the equations for the stepped beam by representing the deformation fields of each step using the expressions in Eq. (1) and Table 1 for the Euler–Bernoulli and the Timoshenko–Ehrenfest models, respectively. To facilitate the following analysis, we define the array $Uk(xk,ω^)$ and matrix $Mk(xk,ω^)$ and use them to describe the TMM without any reference to the beam model used but with the understanding that
$Uk(xk,ω^)=Mk(xk,ω^)Ak={MkEB(xk,ω^)Ak,for EB beamMkTB(xk,ω^)Ak,for TB model$
(14)
where Eqs. (2) and (4) give the specific forms of beam states $Uk(xk,ω^)$ for Euler–Bernouli and Timoshenko–Ehrenfest beam models, respectively. These matrices Mk will be essential in finding the eigen-frequencies of the stepped beam.
Equation (14) can be inverted to express the coefficient array $Ak$ in terms of the state at the beginning of the beam and that in turn is used to provide the state anywhere within each beam:
$Uk(xk,ω^)=Mk(xk,ω^)(Mk−1(0,ω^)Uk(0,ω^))$
(15)
Enforcing the continuity condition, which requires that the displacement and slope are equal across each step, and the compatibility condition, which requires that moment and shear are equal across the steps permits us to observe that the beam states as defined in Eqs. (2) and (4) are continuous across beams:
$Uk(0,ω^)=Uk−1(Lk−1,ω^)$
(16)
which, with Eq. (15), leads us to recursively define matrices $M~k(ω^)$ by
(17)
We may now write
$Uk(xk,ω^)=[Mk(xk,ω^)Mk−1(0,ω^)]M~k−1(ω^)U1(0,ω^)=M^k(xk,0,ω^)M~k−1(ω^)U1(0,ω^)$
(18)
Setting k = N and xk = LN, we obtain
$UN(LN,ω^)=M^N(xN,0,ω^)M~N−1(ω^)U1(0,ω^)=M~N(ω^)U1(0,ω^)$
(19)
Letting $UF=U1(0,ω^)$ and $UL=UN(LN,ω^)$, we now have
$UL=M~N(ω^)UF$
(20)
Let the 2 × 1 vectors $U¯F$ and 02 represent, respectively, the nonzero and the zero rows of UF. Similarly, let $U¯L$ represent the 2 × 1 nonzero rows of UL. By using these definitions, we can describe any homogeneous boundary condition at the two ends of the stepped beam according to
$UF=PFT[U¯F02],UL=PLT[U¯L02]$
(21)
where $PFT$ and $PLT$ are permutation matrices that reorder the entries of UF and UL, respectively, so that the nonzero entries in these vectors float to the top.
Substituting Eq. (21) into Eq. (20) gives
$[U¯L02]=(PLM~N(ω^)PFT)[U¯F02]=M¯(ω^)[U¯F02]$
(22)
We can write the matrix $M¯(ω^)$ using its submatrices $M¯ij$ according to
$[U¯L02]=[M¯11M¯12M¯21M¯22][U¯F02]=[M¯11U¯FM¯21U¯F]$
(23)
where we see that the bottom part of the equation gives rise to the following expression:
$M~21(ω^)U¯F=02$
(24)
which is an eigen-equation for natural frequency $ω^$ and the compatible nonhomogeneous kinematics $U¯F$ at the origin of the stepped beam. We find the beam deformations corresponding to this eigen solution by evaluating UF using the first of Eqs. (21) and employing it in Eq. (18).

### 5.1 Properties of the Transfer Matrix.

To study the properties of the transfer matrix, we investigate the form of the matrix and its determinant for a uniform, homogeneous beam of multiple segments. Specifically, we construct a transfer matrix that describes the compatibility and continuity condition across one interface, and we explore the behavior of this matrix as more segments are added, i.e., as the power of the resulting matrix is increased.

By using these functions, we can infer some of the numerical properties of $M^$ by studying the matrix product $Mk(L,ω^)(Mk(0,ω^))−1$, see Eq. (20). The products $Mk(L,ω^)(Mk(0,ω^))−1$ obtained using the reformulated eigenmodes are listed below for the Euler–Bernoulli (Eq. (25)) and the Timoshenko–Ehrenfest (Eq. (29)) beam models, respectively. A single segment is represented by n = 1 in each of these equations, and n > 1 corresponds to n such segments attached end to end.
$(e−μM^EB(L,0,ω^))n=[M11EBM12EBM21EBM22EB]$
(25)
where the submatrices are given by Eq. (26).
$M11EB=[14(2e−nμcos(nμ)+e−2nμ+1)e−2nμL(2enμsin(nμ)+e2nμ−1)4μe−2nμμ(−2enμsin(nμ)+e2nμ−1)4L14(2e−nμcos(nμ)+e−2nμ+1)]M12EB=[e−2nμL2(−2enμcos(nμ)+e2nμ+1)4EIμ2e−2nμL3(−2enμsin(nμ)+e2nμ−1)4EIμ3e−2nμL(2enμsin(nμ)+e2nμ−1)4EIμe−2nμL2(−2enμcos(nμ)+e2nμ+1)4EIμ2]M21EB=[e−2nμEIμ2(−2enμcos(nμ)+e2nμ+1)4L2e−2nμEIμ(−2enμsin(nμ)+e2nμ−1)4Le−2nμEIμ3(2enμsin(nμ)+e2nμ−1)4L3e−2nμEIμ2(−2enμcos(nμ)+e2nμ+1)4L2]M22EB=[14(2e−nμcos(nμ)+e−2nμ+1)e−2nμL(2enμsin(nμ)+e2nμ−1)4μe−2nμμ(−2enμsin(nμ)+e2nμ−1)4L14(2e−nμcos(nμ)+e−2nμ+1)]$
(26)
Note that all the terms on the right-hand side of Eq. (25) have at most polynomial growth with frequency, regardless of the size of n. Products for dissimilar segments will be similarly well behaved, and the terms (eμ/μ3) can be grouped on the left-hand side. Referring to Eq. (24), these terms factor out of the expression $M^21$ and can be dropped in the eigen-analysis.
Similarly, we can write for a cantilevered Timoshenko–Ehrenfest beam
$PL(M^TB(ω^))nPFT=[M11TBM12TBM21TBM22TB]$
(27)
where each of the submatrices is given by Eq. (29) using the constants
$EIμL=ω^μ^1,AGκμL=ω^μ^2λ+μ2μ=ω^Γ2,λ−ω2ω=ω^Γ1ω=ψμ,whereψ=Δ12+1Δ12−1.$
(28)
$M11TB=[cos(nμψ)−cosh(nμ)ω2(ψΓ1−Γ2)μ1sinh(nμ)Γ1+sin(nμψ)Γ2ωΓ1μ2+ψωΓ2μ2sin(nμψ)Γ1−sinh(nμ)Γ2ψωΓ1μ1−ωΓ2μ1(cosh(nμ)−cos(nμψ))Γ1Γ2(Γ1+ψΓ2)μ2]M12TB=[ψcosh(nμ)Γ1−cos(nμψ)Γ2ψΓ1−Γ2sin(nμψ)(AGκΓ2−μ2)+sinh(nμ)(AGκΓ1+ψμ2)ω(Γ1+ψΓ2)μ2ω(ψsinh(nμ)−sin(nμψ))Γ1Γ2ψΓ1−Γ2cos(nμψ)Γ1(μ2−AGκΓ2)+cosh(nμ)Γ2(AGκΓ1+ψμ2)(Γ1+ψΓ2)μ2]M21TB=[ψcos(nμψ)Γ1−cosh(nμ)Γ2ψΓ1−Γ2ω(ψsin(nμψ)+sinh(nμ))Γ1Γ2μ1(Γ1+ψΓ2)μ2sinh(nμ)(AGκΓ2−μ2)−sin(nμψ)(AGκΓ1+ψμ2)ω(ψΓ1−Γ2)μ1cosh(nμ)Γ1(μ2−AGκΓ2)+cos(nμψ)Γ2(AGκΓ1+ψμ2)(Γ1+ψΓ2)μ2]M22TB=[ψω2(cos(nμψ)−cosh(nμ))Γ1Γ2μ1Γ2−ψΓ1ωμ1(ψsin(nμψ)Γ1(AGκΓ2−μ2)+sinh(nμ)Γ2(AGκΓ1+ψμ2))(Γ1+ψΓ2)μ2ψωsinh(nμ)Γ1(μ2−AGκΓ2)+ωsin(nμψ)Γ2(AGκΓ1+ψμ2)ψΓ1−Γ2(cos(nμψ)−cosh(nμ))(AGκΓ2−μ2)(AGκΓ1+ψμ2)(Γ1+ψΓ2)μ2]$
(29)

Some comments are necessary here about the character of the matrix represented in Eq. (29). Here, we see exponential growth (cosh and sinh) in both n and frequency. This seems to be a feature inherent to the Timoshenko–Ehrenfest beam such that a reformulation of the eigenmode expressions and the direct application of TMM cannot mitigate that growth. This is a topic of continuing effort.

## 6 Verification of Reformulated Method for Euler–Bernoulli Beams Against Finite Element Calculations

Though the cited methods that have been employed up to now appear to have made headway against the computational impediments to this problem, they do not appear to have advanced much beyond problems of six segments. A paradigm is that of Beam V of Ref. [13], where the Euler–Bernoulli beam under discussion had six segments of varying lengths, bending stiffnesses EI, and cross-sectional inertias ρA.

Analysis were done as follows:

• Boundary conditions were selected. (Results of problems involving clamped-free BCs are presented below.)

• The frequency line between 0 and ωmax is divided into many small intervals.

• Using the TMM, the determinant of $M~21$ is evaluated at the frequencies separating intervals.

• The determinants computed earlier are interpolated by cubic splines and zero-crossings are calculated for those splines. The locations of zero crossings are then refined using local meshed in the manner discussed for the full matrix approach.

This is a very inefficient method of root finding, but it is sufficient to assess the computational advantages of the TMM method. Two cases were considered: (1) the beam V geometry of Ref. [13] and (2) a sequence of five Xu Beam V geometries attached end-to-end for a total of 30 segments. In each case, the beam was cantilevered on one end and free at the other end. For comparison, the beam segments were modeled by finely meshed cubic finite elements with lumped masses, e.g., see Ref. [19]. The results for the first ten eigenvalues for the original configuration of Xu et al. and the case of the 30 segment beam are shown in Table 2, where one sees generally good agreement. Better agreement might be expected if a more sophisticated beam element were used in the finite element analysis.

Table 2

Comparison of TMM frequency predictions with those of cubic finite elements for EB beam

Xu 6-segmentXu 6-segment × 5
TMMFETMMFE
3.43773.4640.150070.21716
20.75621.0410.903950.85888
57.93559.9462.38032.3615
118.49115.694.50194.6127
196.11188.367.20087.5972
285.3283.510.89611.398
393.46397.1914.95215.866
520.3527.4219.77121.083
669.7679.225.33827.017
852.93843.5631.54633.483
Xu 6-segmentXu 6-segment × 5
TMMFETMMFE
3.43773.4640.150070.21716
20.75621.0410.903950.85888
57.93559.9462.38032.3615
118.49115.694.50194.6127
196.11188.367.20087.5972
285.3283.510.89611.398
393.46397.1914.95215.866
520.3527.4219.77121.083
669.7679.225.33827.017
852.93843.5631.54633.483

## 7 Verification of Reformulated Method for Timoshenko–Ehrenfest Model Against Finite Element Calculations

The full matrix technique, the transfer matrix technique, and a finite element analysis were performed on a problem of beam consisting of 20 identical segment. The finite element analysis was performed in a matlab code using the beam elements of Friedman and Kosmatka [20] using 100 elements per segment. The problem parameters are provided in Table 3, and the first 13 resonant frequencies are provided in Table 4. We see that for this problem the numerical results of the full matrix method and the transfer matrix method are identical—as they should be since these are mathematically identical. Those results are very similar to the finite element predictions. Differences between the full matrix method and transfer matrix method are illustrated in Fig. 2. The characteristic function obtained numerically from the full matrix method seems to oscillate about zero but with amplitude that grows exponentially. In contrast, the characteristic function obtained numerically from the transfer matrix method seems to oscillate about zero in a linearly growing manner.

Fig. 2
Fig. 2
Close modal
Table 3

Geometric and material properties of Timoshenko–Ehrenfest beam

 Length (L) 0.5 C-S area (A) 1E−4 Moment of inertia (I) 8.33E−10 Young’s modulus (E) 200E9 Shear modulus (G) 75E9 Shear correction coefficient (κ) 5/6 Mass density (ρ) 8050
 Length (L) 0.5 C-S area (A) 1E−4 Moment of inertia (I) 8.33E−10 Young’s modulus (E) 200E9 Shear modulus (G) 75E9 Shear correction coefficient (κ) 5/6 Mass density (ρ) 8050
Table 4

Resonant Frequencies in rad/s for Timoshenko–Ehrenfest beam by three methods

Full matrixTransfer matrixFinite element
1284.01284.01285.8
3526.23526.23534.6
6878.66878.66901.2
11299.911299.911347.2
16755.116755.116839.6
23202.123202.123338.2
30594.630594.630798.1
38882.938882.939170.2
48015.448015.448403.5
57939.557939.558445.1
68602.668602.669242.3
79953.179953.180742.4
91940.791940.792894.2
Full matrixTransfer matrixFinite element
1284.01284.01285.8
3526.23526.23534.6
6878.66878.66901.2
11299.911299.911347.2
16755.116755.116839.6
23202.123202.123338.2
30594.630594.630798.1
38882.938882.939170.2
48015.448015.448403.5
57939.557939.558445.1
68602.668602.669242.3
79953.179953.180742.4
91940.791940.792894.2

Note: Because of the use of lumped mass matrices, one does not expect these finite element frequency calculations to converge consistently from above or from below.

## 8 Contributions and Conclusions

Of greater practical importance, the problem of eigen-analysis of multisegment Euler–Bernoulli beams appears to have been resolved through the combination of reformulation of the eigenmode expressions and application of the transfer matrix method. While previous methods struggled to get accurate solution with a handful of beam steps, this approach has little difficulty in dealing with problems of tens of beam segments.

The second contribution is the observation of a fundamental difference in the behavior of the multisegment TB versus multisegment EB beams: there is a numerical singularity intrinsic to the multisegment Timoshenko beam that is not present in the corresponding Euler–Bernoulli beam. This feature does not appear to be well understood currently and is a topic of continuing research.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## References

1.
Sato
,
H.
,
1983
, “
Free Vibration of Beams With Abrupt Changes of Cross-section
,”
J. Sound. Vib.
,
89
(
1
), pp.
59
64
.
2.
Bapat
,
C.
, and
Bapat
,
C.
,
1987
, “
Natural Frequencies of a Beam With Non-Classical Boundary Conditions and Concentrated Masses
,”
J. Sound. Vib.
,
112
(
1
), pp.
177
182
.
3.
Aydin
,
A. S.
, and
Aksu
,
G.
,
1986
, “
A Finite Difference Method for the Free Vibration Analysis of Stepped Timoshenko Beams and Shafts
,”
Mech. Mach. Theory
,
21
(
1
), pp.
1
12
.
4.
Jaworski
,
J.
, and
Dowell
,
E.
,
2008
, “
Free Vibration of a Cantilevered Beam With Multiple Steps: Comparison of Several Theoretical Methods With Experiment
,”
J. Sound. Vib.
,
312
(
4
), pp.
713
725
.
5.
Lee
,
J.
, and
Bergman
,
L.
,
1994
, “
The Vibration of Stepped Beams and Rectangular Plates by an Elemental Dynamic Flexibility Method
,”
J. Sound. Vib.
,
171
(
5
), pp.
617
640
.
6.
Lu
,
Z.
,
Huang
,
M.
,
Liu
,
J.
,
Chen
,
W.
, and
Liao
,
W.
,
2009
, “
Vibration Analysis of Multiple-Stepped Beams With the Composite Element Model
,”
J. Sound. Vib.
,
322
(
4–5
), pp.
1070
1080
.
7.
Elishakoff
,
I.
,
2019
, “
Who Developed the So-Called Timoshenko Beam Theory?
,”
Math. Mech. Solids
,
25
(
1
), pp.
97
116
.
8.
Aho
,
A. V.
,
Hopcroft
,
J. E.
, and
Ullman
,
J. D.
,
1974
,
The Design and Analysis of Computer Algorithms
,
,
.
9.
Koplow
,
M. A.
,
Bhattacharyya
,
A.
, and
Mann
,
B. P.
,
2006
, “
Closed Form Solutions for the Dynamic Response of Euler–Bernoulli Beams With Step Changes in Cross Section
,”
J. Sound. Vib.
,
295
(
1
), pp.
214
225
.
10.
Yu
,
H.
,
Yang
,
Y.
, and
Yuan
,
Y.
,
2018
, “
Analytical Solution for a Finite Euler–Bernoulli Beam With Single Discontinuity in Section Under Arbitrary Dynamic Loads
,”
Appl. Math. Model.
,
60
, pp.
571
580
.
11.
Stanton
,
S.
, and
Mann
,
B.
,
2010
, “
On the Dynamic Response of Beams With Multiple Geometric or Material Discontinuities
,”
Mech. Syst. Signal. Process.
,
24
(
5
), pp.
1409
1419
. cited By 10.
12.
Young
,
D.
, and
Felgar
,
R. P.
,
1949
, “Tables of Characteristic Functions Representing Nomal Modes of Vibration of a Beam,” University of Texas at Austin, Austin, TX, Technical Report No. 4913.
13.
Xu
,
W.
,
Cao
,
M.
,
Ren
,
Q.
, and
Su
,
Z.
,
2013
, “
Numerical Evaluation of High-Order Modes for Stepped Beam
,”
ASME J. Vib. Acoust.
,
136
(
1
), p.
014503
.
14.
Goncalves
,
P. J. P.
,
Peplow
,
A.
, and
Brennan
,
M. J.
,
2018
, “
Exact Expressions for Numerical Evaluation of High Order Modes of Vibration in Uniform Euler-Bernoulli Beams
,”
Appl. Acoust.
,
141
, pp.
371
373
.
15.
Khasawneh
,
F. A.
, and
Segalman
,
D.
,
2019
, “
Exact and Numerically Stable Expressions for Euler-Bernoulli and Timoshenko Beam Modes
,”
Appl. Acoust.
,
151
, pp.
215
228
.
16.
Yu
,
S. D.
, and
Cleghorn
,
W. L.
,
2000
, “
Free Vibration of a Spinning Stepped Timoshenko Beam
,”
ASME J. Appl. Mech.
,
67
(
4
), pp.
839
841
.
17.
Gartner
,
J. R.
, and
Olgac
,
N.
,
1982
, “
Improved Numerical Computation of Uniform Beam Characteristic Values and Characteristic Functions
,”
J. Sound. Vib.
,
84
(
4
), pp.
481
489
.
18.
Van Rensburg
,
N.
, and
Van der Merwe
,
A.
,
2006
, “
Natural Frequencies and Modes of a TimoShenko Beam
,”
Wave Motion
,
44
(
1
), pp.
58
69
.
19.
Malkus
,
D. S.
,
Cook
,
R. D.
,
Plesha
,
M. E.
, and
Witt
,
R. J.
,
2001
,
Concepts and Applications of Finite Element Analysis
,
Wiley
,
New York
.
20.
Friedman
,
Z.
, and
Kosmatka
,
J.
,
1993
, “
An Improved Two-Node Timoshenko Beam Finite Element
,”
Comput. Struct.
,
47
(
3
), pp.
473
481
.