## 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 4*N* 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.

λ < α |

$[U(x)\Phi (x)]=A1[sinh\mu xL\lambda +\mu 2\mu cosh\mu xL]+A2[cosh\mu xL\lambda +\mu 2\mu sinh\mu xL]+A3[sin\omega xL\u2212\lambda \u2212\omega 2\omega cos\omega xL]+A4[cos\omega xL\lambda \u2212\omega 2\omega sin\omega xL]$. |

λ = α |

$[U(x)\Phi (x)]=A1[01]+A2[1\alpha xL]+A3[sin\omega xL\u2212\alpha \u2212\omega 2\omega cos\omega xL]+A4[cos\omega xL\alpha \u2212\omega 2\omega sin\omega xL]$ |

λ > α |

$[U(x)\Phi (x)]=A1[sin\theta xL\u2212\lambda \u2212\theta 2\theta cos\theta xL]+A2[cos\theta xL\lambda \u2212\theta 2\theta sin\theta xL]+A3[sin\omega xL\u2212\lambda \u2212\omega 2\omega cos\omega xL]+A4[cos\omega xL\lambda \u2212\omega 2\omega sin\omega xL]$ |

λ < α |

$[U(x)\Phi (x)]=A1[sinh\mu xL\lambda +\mu 2\mu cosh\mu xL]+A2[cosh\mu xL\lambda +\mu 2\mu sinh\mu xL]+A3[sin\omega xL\u2212\lambda \u2212\omega 2\omega cos\omega xL]+A4[cos\omega xL\lambda \u2212\omega 2\omega sin\omega xL]$. |

λ = α |

$[U(x)\Phi (x)]=A1[01]+A2[1\alpha xL]+A3[sin\omega xL\u2212\alpha \u2212\omega 2\omega cos\omega xL]+A4[cos\omega xL\alpha \u2212\omega 2\omega sin\omega xL]$ |

λ > α |

$[U(x)\Phi (x)]=A1[sin\theta xL\u2212\lambda \u2212\theta 2\theta cos\theta xL]+A2[cos\theta xL\lambda \u2212\theta 2\theta sin\theta xL]+A3[sin\omega xL\u2212\lambda \u2212\omega 2\omega cos\omega xL]+A4[cos\omega xL\lambda \u2212\omega 2\omega sin\omega xL]$ |

Note: Here, *α* = *AL*^{2}/*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, $\lambda =TT2\omega ^2$, where $\omega ^$ is the beam’s natural frequency, $TT=L\rho /G\kappa $. See Eq. (3*a*) 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

### 2.1 Euler–Bernoulli Beam.

*A*

_{1}–

*A*

_{4}are functions of

*μ*,

*x*is the axial position along the beam’s span,

*L*is the length of the beam,

*μ*= (

*αλ*)

^{1/4},

*α*=

*A L*

^{2}/

*I*, $\lambda =TE\omega ^$, $TE=L\rho /E$, $\omega ^$ 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.

*k*th beam segment such that

*x*

_{k}∈ (0,

*L*

_{k}), where $AkT=[A1k,A2k,A3k,A4k]$. The rows of matrix $MkEB(xk,\omega ^)$ correspond to the displacement, slope, bending moment, and shear force at location

*x*

_{k}of beam segment

*k*, respectively.

### 2.2 Timoshenko–Ehrenfest Beam.

*λ*/

*α*as shown in Table 1. Define the dimensionless constant

*γ*=

*Gκ*

^{2}/

*E*. The real quantities

*ω*,

*μ*(for the case

*λ*<

*α*), and

*θ*(for the case

*λ*>

*α*) are given in terms of

*λ*as Ref. [18]

*a*)

*b*)

*c*)

*d*)

*μ*and

*ω*are functions of the dimensionless frequency $\lambda $, which is related to the physical frequency $\omega ^$ by $\lambda =\omega ^TT$, where $TT=L\rho /(G\kappa )$,

*G*is the shear modulus, and

*κ*is the shear correction coefficient.

## 3 Reformulation of the Eigenmode Expressions

The enabling idea to reformulating the eigenmode expressions comes from the observation that the term *A*_{1}sinh *μ*(*x*/*L*) + *A*_{2}cosh *μ*(*x*/*L*) in Eq. (1) and both the previous term and the term *A*_{1}cosh *μ*(*x*/*L*) + *A*_{2}sinh *μ*(*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 *A*_{1} and *A*_{2}.

*e*

^{μ(x/L)}that are relevant to the boundedness issue:

*P*(

*μ*), where

*P*(

*μ*) is a bounded function of

*μ*. We can now express

*A*

_{2}in terms of

*A*

_{1}and

*P*(

*μ*) according to

### 3.1 The Euler–Bernoulli Beam.

### 3.2 The Timoshenko Beam.

*P*(

*μ*) is yet unspecified other than being a bounded function of

*μ*, we may just replace the term −

*A*

_{1}

*P*(

*μ*)/2 with the unknown coefficient $A~2$. Further, we also replace −

*A*

_{1}with $A~1$ to obtain

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

To capture the two homogeneous boundary conditions at the external boundaries of the first beam, we define the 2 × 4 matrix *B*_{F} consisting of the rows of $M1(0,\omega ^)$ 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 *B*_{L} consisting of the rows of $ML(LL,\omega ^)$ associated to those two degrees-of-freedom.

## 5 Transfer Matrix Method

*M*

_{k}will be essential in finding the eigen-frequencies of the stepped beam.

*k*=

*N*and

*x*

_{k}=

*L*

_{N}, we obtain

_{2}represent, respectively, the nonzero and the zero rows of

*U*

_{F}. Similarly, let $U\xafL$ represent the 2 × 1 nonzero rows of

*U*

_{L}. By using these definitions, we can describe any homogeneous boundary condition at the two ends of the stepped beam according to

*U*

_{F}and

*U*

_{L}, respectively, so that the nonzero entries in these vectors float to the top.

*U*

_{F}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.

*n*= 1 in each of these equations, and

*n*> 1 corresponds to

*n*such segments attached end to end.

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

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.

Xu 6-segment | Xu 6-segment × 5 | ||
---|---|---|---|

TMM | FE | TMM | FE |

3.4377 | 3.464 | 0.15007 | 0.21716 |

20.756 | 21.041 | 0.90395 | 0.85888 |

57.935 | 59.946 | 2.3803 | 2.3615 |

118.49 | 115.69 | 4.5019 | 4.6127 |

196.11 | 188.36 | 7.2008 | 7.5972 |

285.3 | 283.5 | 10.896 | 11.398 |

393.46 | 397.19 | 14.952 | 15.866 |

520.3 | 527.42 | 19.771 | 21.083 |

669.7 | 679.2 | 25.338 | 27.017 |

852.93 | 843.56 | 31.546 | 33.483 |

Xu 6-segment | Xu 6-segment × 5 | ||
---|---|---|---|

TMM | FE | TMM | FE |

3.4377 | 3.464 | 0.15007 | 0.21716 |

20.756 | 21.041 | 0.90395 | 0.85888 |

57.935 | 59.946 | 2.3803 | 2.3615 |

118.49 | 115.69 | 4.5019 | 4.6127 |

196.11 | 188.36 | 7.2008 | 7.5972 |

285.3 | 283.5 | 10.896 | 11.398 |

393.46 | 397.19 | 14.952 | 15.866 |

520.3 | 527.42 | 19.771 | 21.083 |

669.7 | 679.2 | 25.338 | 27.017 |

852.93 | 843.56 | 31.546 | 33.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.

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 |

Full matrix | Transfer matrix | Finite element |
---|---|---|

1284.0 | 1284.0 | 1285.8 |

3526.2 | 3526.2 | 3534.6 |

6878.6 | 6878.6 | 6901.2 |

11299.9 | 11299.9 | 11347.2 |

16755.1 | 16755.1 | 16839.6 |

23202.1 | 23202.1 | 23338.2 |

30594.6 | 30594.6 | 30798.1 |

38882.9 | 38882.9 | 39170.2 |

48015.4 | 48015.4 | 48403.5 |

57939.5 | 57939.5 | 58445.1 |

68602.6 | 68602.6 | 69242.3 |

79953.1 | 79953.1 | 80742.4 |

91940.7 | 91940.7 | 92894.2 |

Full matrix | Transfer matrix | Finite element |
---|---|---|

1284.0 | 1284.0 | 1285.8 |

3526.2 | 3526.2 | 3534.6 |

6878.6 | 6878.6 | 6901.2 |

11299.9 | 11299.9 | 11347.2 |

16755.1 | 16755.1 | 16839.6 |

23202.1 | 23202.1 | 23338.2 |

30594.6 | 30594.6 | 30798.1 |

38882.9 | 38882.9 | 39170.2 |

48015.4 | 48015.4 | 48403.5 |

57939.5 | 57939.5 | 58445.1 |

68602.6 | 68602.6 | 69242.3 |

79953.1 | 79953.1 | 80742.4 |

91940.7 | 91940.7 | 92894.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.