Shock and Vibration
Volume 2018 (2018), Article ID 2487697, 21 pages
https://doi.org/10.1155/2018/2487697

Research Article

A Chebyshev Spectral Method with Null Space Approach for Boundary-Value Problems of Euler-Bernoulli Beam

C. P. Hsu,1,2 C. F. Hung,¹ and J. Y. Liao²

¹Department of Engineering Science and Ocean Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd., Taipei, Taiwan
²National Chung Shang Institute of Science and Technology, PB 90008-6-16, Longtang, Taoyuan, Taiwan

Correspondence should be addressed to C. F. Hung; hungef@ntu.edu.tw

Received 7 March 2018; Accepted 21 June 2018; Published 16 July 2018

Academic Editor: José J. Rangel-Magdaleno

Copyright © 2018 C. P. Hsu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We proposed a Chebyshev spectral method with a null space approach for investigating the boundary-value problem of a nonprismatic Euler-Bernoulli beam with generalized boundary or interface conditions. It is shown here that, with few vital improvements, a Chebyshev spectral collocation approach can be systematically applied to modeling nonprismatic Euler-Bernoulli beams with eigenvalue embedded tip-massed boundary conditions as well as the jump conditions that appear at the stepped interfaces. This study also presents a numerical stable asymptotic modal solution for the higher-order modes of a partially clamped beam and show that the proposed approach validates the robust higher-order modal solutions. Through a sequence of four increasingly complicated examples, using the proposed approach with higher-order modes, generalized boundary conditions, and interface jump conditions of nonprismatic beams, the results are in excellent agreement with those reported in the literature using various other approaches.

1. Introduction

Vibration analyses of Euler-Bernoulli (EB) beams are widely used in mechanical, civil, and aerospace engineering [1, 2]. According to a variety of in situ applications, modal solutions of nonprismatic beams with generalized boundary conditions (BCs) are prerequisite in preliminary design. In addition to classical boundary conditions in many vibration textbooks, generalized boundary conditions such as damaged support, tip-mass, multiple-beams connected with elastic stiffness [3], and imperfect beams [4] have been adopted in many applications. For example, Olgac and Jalili [5] design a delayed resonator vibration absorber with a partially clamped beam. Moreover, beam models composed of a double layer carbon nanotube with elastically supported foundations and partially clamped boundary conditions are widely used to emulate the natural frequencies of integrated circuit assemblies under different external forces [6]. Recently, based on Jeong et al. [7], we developed a mechanically excited resonant beam for a pyroshock source simulator by using a partially clamped EB, and a configuration of a two-step beam is surveyed for the evaluation of the effects on adapter plates and test items.

Although there are a lot of papers in the literature that utilize nonclassical boundary conditions to analyze the free vibration of beams, they often tend to restrict their analysis to a low-order mode shapes, usually less than a dozen, while higher-order modes (more than 30) are rarely discussed. Nevertheless, Goncalves, Brennan and Elliott [8] investigate the numerical evaluation of higher-order modes of vibrations in uniform EB beams with various classical boundary conditions, and some asymptotic expressions are obtained. In this paper, a more general approach for partially clamped EB beams is presented and we obtain a unique set of asymptotic modal solution for higher-order modes. Cha, Wang, and Liao [9] analyzed chatter stability of surface grinder machines, modeled by a moving lumped grinding wheel element over an elastic worktable, with the last assumed to be a partially clamped EB beam. The assumed mode expansion is used for the EB beam, and asymptotic higher-order modes are required to facilitate the analysis.

The exact analytical solutions for free vibration of EB beams have been given only for prismatic beams with simply supported boundary conditions. Lai, Hsu, and Chen [10, 11] use a modified Adomian decomposition method (ADM) to analyze the modal solutions of a tip-massed tapered beam. However, the tip-masses at the ends of the beam introduce eigenvalue embedded boundary conditions that further complicate the modal solutions of the tip-massed EB beam. In addition, Mao and Pietrzko [12] discuss the modal solutions of a stepped beam using ADM. For a beam with a discontinuous cross-section, together with the boundary conditions at both ends, there are also interface conditions at the step location of the two beam segments. However, the interface conditions must obey displacement compatibility and force equilibrium relationships. Again, Mao discusses the modal solutions of nonprismatic EB beams for piezoelectric modal sensors by using the transformation matrix method of [13]. Moreover, the modal solutions of a stepped beam with multiple cracks and different boundary conditions are discussed in Attar [14]. Lastly, Sarkar and Ganguli [15] provide a modal tailoring and closed-form solutions for a rotating nonprismatic beam by using an inverse problem approach.

Besides the above-mentioned approach, the spectral method is another numerical approach for EB beams, which can potentially provide superior accuracy and domain flexibility [16, 17]. Because of the superior convergence of the Chebyshev spectral method (CSM), it has been widely used for the modal solution of various Timoshenko beams [18–21]. However, in the Mathematical Institute at University of Oxford, Professor Trefethen has initiated and held the Chebfun project since 2002. The Chebfun toolbox under the MATLAB environment enables us to formulate the spectral collocation matrix automatically for the integral and integrodifferential equations [22–27]. Driscoll, Bornemann, and Trefethen [25] solve linear differential equations based on the developed Chebfun toolbox with the object-oriented MATLAB symbolic scripts and introduce the “chebop” system, where the ODE and boundary conditions may be represented with symbolic chebop commands.

For boundary and interface conditions issues that arise in CSM-type approaches, Driscoll and Hale [26] suggest a rectangular spectral collocation method (RSC) without any rows deletions to formulate the corresponding CSM. Moreover, Smith, Laoulache, and Heryudono [28] focus on the boundary conditions of CSM and discuss the row replacement method, the fictitious point method, and the rectangular collocation method, all implemented as velocity and vorticity conditions for solving the unsteady Navier-Stokes equations of an annular Couette flow.

As in the classical boundary conditions and low-order modal solutions, the eigensolutions of CSM may be easily obtained by using the overloaded “eigs” command in the Chebfun toolbox. However, in some cases such as with nonclassical boundary conditions and higher-order modes, it is not always working well. These nonclassical boundary conditions may result in an ill-conditioned matrix and numerical instability in the eigensolutions. Wilkinson [29] explains the reasons of an ill-condition matrix and Saad [30] provides a typical preconditioned iterative approach for computing large eigenvalue problems. Since in the null space approach we project the original solution space to its constrained space instead, many null space applications such as solving nonlinear system of equations [31], data reconciliation, and variable classification have been applied in of Mitsas [32]. Considering that the rigid-body modes are in the null space of the matrix of the structure system, Felippa and Park [33] apply the null space approach to the free-free flexibility matrix of the structures and explore the relationship between the rigid-body modes and boundary conditions.

The main purpose of this paper is to propose a systematic procedure using the Chebyshev spectral method with the null space approach (NSA) for the ill-conditioned system matrix caused by modeling the nonprismatic EB beams with generalized boundary and interface conditions. Another contribution of this paper is to introduce the chebop system in the Chebfun toolbox for modeling the governing equation with boundary or interface conditions of EB beams. The symbolic operator commands in the chebop are dramatically simple and one may represent the governing equation and boundary or interface conditions of an EB beam in just a few lines long executable MATLAB code. After constructing the CSM collocated matrix, one may use the NSA to project the ill-conditioned matrix to the null space of the boundary or interface conditions, thereby significantly improving the condition of the spectral collocation matrix. For validation purposes, the exact modal solutions of the simply supported EB beam are used to benchmark the accuracy of higher-order modal solutions by CSM with NSA. This paper also derives higher-order asymptotic modal solutions for the partially clamped EB beam and benchmarks the accuracy of the derived asymptotic solution with the proposed approach. Finally, several examples of a nonprismatic EB beam with generalized boundary and interface conditions, such as tip-massed and stepped beams, are used to demonstrate the proposed CSM with NSA approach. Those examples show the straightforward symbolic formulation in accordance with the governing equation and boundary conditions of an EB beam, in obtaining the robust eigensolutions of boundary-value problems.

2. Governing Equation and Boundary Conditions of EB Beam

2.1. The Governing Equation of a Nonprismatic EB Beam with Tip-Masses of Rotary Inertia and Nonclassical Boundary Conditions

Consider a nonprismatic EB beam as shown in Figure 1, whose length is l with both ends restrained by the translation and torsion springs. At the both ends of the beam, there are concentrated masses that comprise the rotatory inertia of the mass.

Figure 1: A nonprismatic EB beam with elastically restrained ends, which supports tip-masses of rotatory inertia at both ends.

The equation of motion for the transverse vibration of the beam is given by

where y(x, t) denotes the transverse deflection of the beam, A(x) is the cross-sectional area at x, I(x) is the area moment of inertia of A(x), and ρ and E are the mass density (mass per unit volume) and Young’s modulus, respectively.

Assuming the beam is oscillating harmonically at a frequency (rad/s) and the solutions of transverse deflection to (1) are

(1) becomes

The corresponding external boundary conditions at both ends are given by

at x = 0, and

at x = l.

Here, kTR, kRR, MR, JR, eR and kTL, kRL, ML, JL, eL, borrowed from the nomenclatures of [9], are the translational springs, the torsion springs, the concentrated masses attached at beam tip, the moments of inertia of the tip-masses, and the eccentricities which are the distances between the beam tip and the center of the tip-mass at the right end and left end of the beam, respectively.

Equations (3) to (7) can be converted to a nondimensional form by introducing the following quantities:

where the range of X is , and A₀ and I₀ are the cross-sectional area and the moment of inertia at x = 0, respectively.

Substituting (8) into (3), the nondimensional form of (3) is

Using (8) in (4) to (7), the boundary conditions at both ends in nondimensional form become as follows:

where

2.2. The Governing Equation of a Stepped EB Beam with Nonclassical Boundary Conditions

Consider the free vibration of a straight EB beam consisting of N uniform sections elastically restrained at both ends, as shown in Figure 2. The total length is l, and the ratio of the distance from the right end of the Sj section to the left end of the beam is . Let us define the junction parameters ηj, γj, and μj for the beam sections j and j+1 (j = 1, 2, …, N-1) as follows:

Here Ej, , Aj, and Ij are Young’s modulus, density per length, the cross-section area, and the inertia of cross-section area of j-th section, respectively. Assuming the stepped beam is oscillating harmonically at a frequency ω (rad/s), the nondimensional governing equation for each segment can be written as

where

Figure 2: A stepped beam with elastically restrained ends.

The interface conditions at the junction of each segment are

and the boundary conditions at the both ends are

2.3. Numerical Instability in Higher-Order Modes of an EB Beam for Partially Clamped Boundary Conditions

This section derives an asymptotic solution for the higher-order modes of a partially clamped EB beam by using trigonometric and modified hyperbolic functions. The evaluation of hyperbolic functions [8] in the analytical solutions of higher-order modes of a partially clamped EB beam will cause numerical instability issues. The numerical ill-conditions occur in the argument of the hyperbolic functions and lead to a quantity larger than that can be represented by the machine. We derive the analytical solution of a partially clamped EB beam and its asymptotic forms for higher-order modes. A detailed derivation is given in the Appendix for the interested reader.

The governing equation of a partially clamped uniform EB beam is shown in (1). The partially clamped condition may be represented by elastically restrained translation and torsion springs. The corresponding boundary conditions at both ends of the beam can be expressed as

Specifically, ϕ is written as a finite sum, the so-called Galerkin approximation:

The assumed mode expansion of can be written in terms of trigonometric and hyperbolic functions as

Let us define the following nondimensional spring constants,

The boundary conditions in (19) can be given by the following nondimensional forms.

According to descriptions found in the Appendix, the analytical modal solution of the beam with partially clamped boundary conditions can be described by

where Λi = k(ωi)l = kil is a nondimensional natural frequency, S, C, SH, and CH represent sin(Λi), cos(Λi), sinh(Λi), and cosh(Λi), respectively.

However, (24) is numerically unstable using double precision floating point arithmetic with MATLAB for mode order greater than 13, due to the presence of hyperbolic functions (e.g., MATLAB cannot tell the difference between sinh13π−1 and sinh13π). The key step in approximating (24) for higher-order mode shapes is to use the modified exponent functions replacing the original hyperbolic functions. This step can solve the problem of divergence of the analytical solution while preserving the orthogonality of the higher-order modes. Refer to the Appendix for the detailed higher-order asymptotic modal solutions. The asymptotic solution for higher-order mode shapes of a uniform EB beam with partially clamped at both ends is shown below:

To tackle the numerical error issue indicated in [8], we apply (24) with Λi ≤ 6π instead of Λi ≤ 13π as a breakpoint to compute the lower-order analytical solutions and the higher-order asymptotic solutions. Therefore, for Λi ≤ 6π, use (24) to obtain the mode shapes; otherwise use (25).

3. The Chebyshev Spectral Method and the Differential Operator in ODE

From [16], the kth-order Chebyshev polynomial TN(s) and Chebyshev-Gauss-Lobatto (CGL) collocation points sj, in the domain of , can be expressed, respectively, by

The derivative of a smooth function at the CGL points sj can then be computed via matrix-vector multiplication, which can be formally represented as

where is a (N+1)×(N+1) matrix.

The elements within the matrix can be expressed as

Moreover, its higher derivative matrices may be expressed as

Thus, the differential operators in the governing equation and the BCs can be approximated by the differential matrix with the Chebyshev spectral approach. As an example, for a prismatic EB beam, the differential governing equation in (9) can be easily transformed to the following matrix eigenvalue problem:

where is a transformation of CGL collocation points sj, Xi = (si +1)/2, and

It should be noted that the LHS summed up matrices in (30) form the so-called system matrix of a nonprismatic EB beam.

4. The Chebyshev Spectral Method for Boundary-Value Problems of an EB Beam

There are two kinds of boundary conditions for the nonprismatic EB beam problem: one is the homogeneous Dirichlet boundary condition and the other is the boundary condition embedded with an eigenvalue, for example, tip-masses at both ends of the beam. The system matrix for the first kind boundary-value problem can be solved straightforwardly with a standard eigensolver. However, it is not easy to formulate the eigenvalue problem for the system matrix with the second kind boundary-value conditions imposed.

For the second kind boundary-value condition, we will briefly introduce the row replacement and the rectangular spectral collocation approaches in [21, 26] and present a detail null space approach to the corresponding boundary-value problem of an EB beam.

4.1. Row Replacement and Rectangular Spectral Collocation Approaches.

From (9) to (13), one may obtain a generalized eigenproblem as follows:

where denotes the CGL collocation vector and is the system matrix.

Vector can be decomposed as , , and . is the collection of homogeneous boundary and interface conditions. is the matrix of boundary and interface conditions containing the eigenvalues. is the residual matrix of the governing equation of an EB beam. Then, (32) can be partitioned as

The homogeneous boundary conditions should satisfy the following equation:

When boundary conditions include the eigenvalues, such as in (10) to (13), the solutions should satisfy the following equation:

Substituting (34) and (35) into (33), by the row replacement approach we obtain the following expression:

Simplifying (36) by reducing the order of the matrix, we obtain the final expression with eigenvalues of boundary conditions:

where

The row replacement method can be used to obtain the feasible eigensolutions for the well-conditioned system matrix after being reduced from the boundary conditions. However, for the classical two-point boundary-value problem, the eigensolutions can be achieved easier by removing the top and bottom rows of the differential matrix. However, for EB beams with interface jump conditions, row replacement can become rather complicated in building up the spectral collocation matrix in the LHS of (36). Moreover, the row replacement method does not always produce eigenvalues that have the same qualitative properties, such as symmetry, as the original operator. Thus, the matrix in (37), even for the simplest differential equation, is full and ill-conditioned, so it is in general not advisable to solve (37) using a direct method for large N. Instead, some appropriate preconditioning approach should be used in [16].

Driscoll and Hale [26] propose a rectangular spectral approach, based upon resampling the differentiated polynomials into a lower-degree subspace to make the differentiation matrices. The rectangular collocation approach is described as follows:

(a)The square system matrix (N×N) originally constructed according to the CSM is allocated for N+m columns by barycentric resampling discretization and expanded to a rectangular matrix ( N× [N+m], where m is the number of boundary conditions).(b)The m sets of boundary conditions formed by CSM are then placed above the rectangular matrix and reconstruct a square matrix for eigensolver.

This approach avoids ambiguities that arise when applying the classical row deletion method in boundary-value problems. However, as indicated in [28], it shows several orders of magnitude less accuracy when evaluated using the resampling points, but higher accuracy when evaluated directly at the original collocation points.

4.2. Projection by Using the Null Space Approach

In this paper, we propose a novel null space approach to solve the eigenvalue problems of the system matrix by preserving its original dimension and characteristics. According to the definition, by (33) we have ,, and matrix represents the imposed boundary and interface conditions. The rank of matrices [Ds] and [] is N and m (N > m), respectively. The eigenproblem of (33) must satisfy the requirement, []=  . Let the null space of [] be matrix [Nc] and satisfy the identity [][Nc] = . It should be noted that is a vector in physical space, and it may be projected to vector in the null space with matrix [Nc]; i.e., = [Nc] As a result, (33) may be projected to the null space of [] as follows:

Premultiplying both sides of (39) with the transposed [Nc], and taking advantage of the orthonormal property in [Nc], we get

Consequently, (33) is reduced to (40) by projecting the system matrix of the physical space to the null space which is imposed by the boundary and interface conditions. For , an eigenvector in (40), one may restore the eigenvector in physical space by setting = [Nc].

The null space approach may actually improve the matrix condition of EB beams. For problems with large variations or consisting of independent rows or columns in the elements of matrix [Ds], a robust null space preconditioning approach will improve the potential ill-conditioned problems, such as the system matrix formulated from a partially clamped EB beam, and overcome the numerical issues in the modal solutions that include both high and low order modes.

4.3. Roadmap for EB Beam Eigenvalue Problems with the Aid of the Chebfun Toolbox

In this section, we discuss the solution roadmap by using the Chebfun toolbox and deal with EB beams with various boundary and interface conditions. The Chebfun toolbox in MATLAB environment is developed by Professor L. N. Trefethen in the Cambridge University [24, 25]. This toolbox is a collection of commands based on a Chebyshev expansion algorithm, where the number of expansion terms is controlled by mechanical precision. The chebop utility in the Chebfun toolbox enables a symbolic operator syntax to be applied to numerical objects and implements the recasting process of governing equations of differential equations. In addition, this toolbox is constructed in MATLAB by matrices and vectors and uses automatic Chebyshev polynomial interpolation to represent functions and automatic spectral collocation methods to approximate operators. With the MATLAB built-in numerical functions of the toolbox, one may quickly and accurately get the asymptotic solutions of various problems.

As an illustrated example, Figure 3 is a solution roadmap for a partially clamped EB beam (see (16) and (17)). In the first step, the governing equation of the EB beam with partially clamped boundary conditions is constructed by using the chebop syntax. Through the chebop system, one may easily define the corresponding EB beam governing equation as well as the left and right boundary conditions (L.op, L.lbc, and L.rbc) and convert them into the Chebyshev spectral matrix in the MATLAB working space. Next, it shows the null space approach for building the orthonormal space of the boundary conditions (nuA) and the projection of the governing differential matrix into the corresponding null space (Aend) thus obtaining the solution with the MATLAB built-in eigensolver. The last step is projecting the eigenvectors backward into the physical space (Vx), which represents the physical modal shapes of the EB beam.

Figure 3: The MATLAB scripts and procedures of CSM with NSA.

4.4. Modal Analysis of a Nonprismatic Tip-Massed EB Beam with Partially Clamped Boundary Conditions

We demonstrate how to use CSM together with NSA through the Chebfun toolbox to solve the nonprismatic EB beam that tip-massed and elastically restrained at both ends.

The governing equation in (9) may be expressed with differential operators as follows:

where

The derivative of (41) with respect to X is

where

For boundary conditions (10)-(13), all the terms containing can be removed by substituting the RHSs of (41) and (42), together with a simple manipulation to get

After that, one may follow the roadmap in Figure 3, where the boundary conditions (43) to (46) can be converted to the [Cs] (4×N) matrix by using N-points CSM, and we get the corresponding (N-4) ×N matrix of the null space. With the NSA (see (39)-(40)) applied, the modal solutions can then be obtained by using a standard eigensolver.

4.5. Modal Analysis of Stepped EB Beams

For a stepped EB beam, (16) and (17) form the governing equation and interface conditions for each segment, and (18) is the boundary condition at both ends. Equations (16) and (18) can be defined with the help of the chebop system as shown in the previous section. According to the Rj ratio at each segment, (16) may be transformed into a Chebfun object by multiplied with the Heaviside function:

On the other hand, the LHS of (18) can be expressed using a differential operator notation as

while the LHS of (18) can be rewritten using the differential operator notation by

Furthermore, in (17), the displacement compatibility (identical jump conditions) and force equilibrium (scaled jump conditions) at the interface of each section can be established by using the jump command in chebop system, as follows:

 [jump; jump; jump; jump].

Since the jump command in the chebop system can only handle the identical jump conditions but is unable to exactly convert the scaled jump conditions of the stepped beam, extra manipulations on the rows of the interface conditions in the differential matrix are needed to define the scaled jump conditions.

For example, a two-step EB beam is shown in Figure 4. With ten CGL collocated points for each beam segment, its differential matrix should consist of 30 rows and columns. In that matrix, the first four rows are formulated using the boundary conditions of both ends, while the fifth to twelfth rows are formulated using the displacement compatibility (identical jump) conditions and force equilibrium (scaled jump) conditions at the interface of the steps. It should be noted that the scaled jump conditions, appearing on rows 7, 8, 11, and 12, should be modified with the scale ηj in (17). The last 18 rows (the DM part in Figure 4) are the governing equation of each beam segment stacked up diagonally with each corresponding differential matrix. After that, the corresponding null space matrix is obtained by using these first 12-row partition (the BIC part in Figure 4). With the NSA (see (39) and (40)) applied, the modal solutions can be obtained by using a standard eigensolver.

Figure 4: The CSM differential matrix with ten CGL collocated points for three-steps EB beam.

5. Numerical Calculations

In order to verify the proposed CSM with NSA approach for analyzing the free vibration of EB beam step by step, four numerical examples with different beam styles and boundary conditions are discussed in this section. All results were obtained using MATLAB R2013B and Chebfun 5.1 toolbox on a HP EliteDesk with an Intel Core i5-6500MQ @ 3.30 GHz CPU and 4 GB RAM.

5.1. Simply Supported Prismatic EB Beam

The exact modal solution of a simply supported prismatic EB beam is used for validating the accuracy of the proposed CSM with NSA approach as well as the CSM with RSC (resampling spectral collocation) approach. For a simply supported prismatic EB beam, the exact nondimensional modal frequencies are nπ, and the corresponding mode shapes are sin(nπ), respectively, where n is the mode order.

With the Chebfun toolbox, one may use the chebop system to define the governing equation and boundary conditions and can obtain the eigensolution using the overloaded eigs command in the Chebfun toolbox. It should be noted that this eigensolution is based on the RSC method built in the Chebfun toolbox. The proposed CSM with NSA approach is illustrated in Figure 3, which projects the differential matrix from the EB beam to the null space of the boundary conditions. In order to benchmark the accuracy of the eigensolutions with these two approaches, the natural frequencies of the 15th and 100th mode are solved separately.

Table 1 compares the exact natural frequencies with the RSC and present approach (CSM with NSA) for a simply supported EB beam. The selected CGL collocated points for the differential matrices are N=30 and N=200, respectively. However, in accordance with the sampling theorem, the maximum number of modes (nM) should not be larger than N/2 in order to get rid of aliasing. In Table 1, we display the cases of nM =10 and nM =100 only. As the results indicate, the accuracy of the present approach outperforms the RSC approach in both cases. However, for the RSC approach with N=200, results are unclear since the lower-order modes are distorted. We conjecture that the resampling algorithm of RSC goes unstable and oversample the lower-order modes as claimed in [28]. Consequently, we suggest the following rule of thumb for the present approach to select the number of CGL points:

The eigenvectors (mode shapes) solved by these two approaches can be compared to the exact solution, and we obtain the normalized maxi-max error as follows:

Following the criteria in (50) for selecting number of CGL collocated points and the maximum number of modes (nM), Figure 5 and Table 2 show the normalized maxi-max errors of the natural frequencies and mode shapes in the first five modes, analyzed by using the RSC and NSA approaches. It is clear that the accuracy of the modal solutions by using the NSA approach outperforms the RSC approach at all chosen nM.

Table 1: Comparison of dimensionless natural frequencies for nM =10 and nM =100 by exact solutions, RSC, and present methods.

CGL points
nM
Method

…

Exact
3.14159
6.28319
9.42478
…
40.84070
43.98230
47.12389

30
10
RSC
3.14076
6.29299
9.42456
…
40.84072
43.98240
47.12414

Present
3.14159
6.28319
9.42478
…
40.84070
43.98225
47.12390

Exact
3.14159
6.28319
9.42478
…
307.87608
311.01767
314.15927

200
100
RSC
3.709444+ (3.70944i)
8.96197
9.05268
…
307.87394
311.04696
314.16051

Present
3.14159
6.28319
9.42478
…
307.87608
311.01767
314.15927

Table 2: Maxi-max errors of natural frequencies and mode shapes by using RSC and NSA approaches.

Mode order
Maxi-max errors of natural frequencies
Maxi-max errors of mode shapes

RSC
NSA
RSC
NSA

5
2.73E-06
2.09E-11
1.39E-06
6.77E-11

10
1.69E-03
2.79E-11
1.44E-02
3.33E-11

15 to 20
0.001685759
< 1.56E-09
0.014375
< 7.33E-10

25 to 45
0.556364694
< 1.53E-08
0.332986
< 5.01E-08

30 to 45
0.556364694
2.09E-09
0.332986
2.17E-08

50 to 110
1.046748562
3.86E-10
0.865779
< 4.59E-07

115 to 170
1.755796169
< 1.35E-09
0.931819
< 2.96E-08

175 to 200
0.909712404
< 3.59E-09
0.994279
< 5.22E-08

Figure 5: Normalized maxi-max errors of frequencies and mode shapes by using RSC and NSA approaches.

5.2. Higher-Order Modal Solutions for a Partially Clamped Prismatic EB Beam

The partially clamped conditions of an EB beam are defined by four elastic restrained coefficients β₁, β₂, β₃, and β₄ in (22). The higher-order asymptotic modal solution, obtained in (25), is denoted by PC₁. Moreover, PC₂ denotes the null space approach (NSA) as shown in Figure 3 for dealing with the eigenproblem of a partially clamped prismatic EB beam. The exact modal solution of a partially clamped uniform EB beam is not found in the literature, to the best of our knowledge. The modal solutions with classical boundary conditions (clamped (C), free (F), pinned (P), slider (S), and simply supported (S-S)) are adopted as a baseline and denoted by PC₀, such as the asymptotic solution proposed in [8] for higher-order modal solutions.

The asymptotic magnitude of the elastic coefficients for free and fixed conditions of a partially clamped beam is around 0 and 10⁸, respectively. PC₁ is substituted in combination with the asymptotic free and fixed βi into (25). In Figure 3, PC₂ uses the Chebfun toolbox to solve the modal solutions with the asymptotic magnitudes of βi. For PC0, we use the classical boundary conditions as proposed in [8]. With the classical C-C, C-F, C-S, and C-P boundary conditions, the modal solutions for the first 200 modes of the EB beam are obtained by higher-order asymptotic analysis of the partial-clamped solutions (PC₁) and the NSA approach (PC₂). Then, these results are compared to the asymptotic modal solutions as proposed in [8] (PC₀). The normalized maxi-max deviations of the natural frequency between PC₁, PC₂, and PC₀ are defined by

and the normalized maxi-max deviations of the mode shape between PC₁, PC₂, and PC₀ are defined by

The maxi-max deviations of the modal solutions at the specified modes of order (10:10:100 and 125:25:200) are shown in Figures 6 and 7 and Table 3, respectively. As we can see, both PC₁ and PC₂ approaches attain around the same accuracy level. In this case, they are in excellent agreement with the modal solution of both approaches and this encourages us to conclude the following:

(a)The result from PC₁ validates the asymptotic higher-order modal solutions of a partial-clamped beam in the Appendix and avoids the numerical problems that arise at higher-order modes.(b)The accuracy of the proposed CSM with NSA approach is nearly equal to the asymptotic analytical solution. This approach offers a robust and systematic procedure without going through sophisticated symbolic derivations needed for the EB beam with a variety of boundary conditions.

Table 3: Normalized maxi-max deviations of modal solutions for PC₁ and PC₂ with four classical BCs imposed.

Boundary condition
Mode order
Maxi-max deviations of natural frequency
Maxi-max deviations of mode shape

PC₁
PC₂
PC₁
PC₂

C-C
10 to 200
2.44E-06
2.44E-06
0.001926
0.001926

C-F
10 to 90
3.70E-06
3.70E-06
0.001014
3.70E-06

100 to 200
3.70E-06
3.70E-06
0.001014
< 2.69E-04

C-P
10 to 50
2.34E-06
2.34E-06
2.84E-06
2.84E-06

60 to 80
2.34E-06
2.34E-06
7.35E-06
7.35E-06

90 to 100
2.34E-06
2.34E-06
8.49E-06
8.49E-06

125 to 200
2.34E-06
2.34E-06
< 1.04E-04
< 1.04E-04

C-S
10 to 200
2.34E-06
2.34E-06
0.075393
0.075393

Figure 6: Normalized maxi-max deviations of natural frequency by using PC1 and PC2 ( and ).

Figure 7: Normalized maxi-max deviations of mode shape by using PC1 and PC2 ( and ).

5.3. Modal Solution for a Partially Clamped and Tip-Massed Nonprismatic EB Beam

In this example, we attempt to validate the robustness of our proposed approach for a tip-massed nonprismatic EB beam with a variety of BCs. It should be noted that, for a tip-massed EB beam, the boundary conditions are embedded in the eigenvalue to be solved. The governing equation (9) and BCs (see (10)-(12)) are formulated with the Chebfun toolbox and the chebop system through (42)-(46) in Section 4.4. Then the modal solutions are obtained with the NSA projection and solved by the MATLAB eig command. However, for comparison with [11], the BCs of the beam are tip-massed at the left end and fixed at the right end (β₃ → ∞, β₄ → ∞), while the other relevant beam parameters have the same value as those in [11] (γL = 0, 0.3 and 0.6, α = 1.2 and 2; μL = 0.2 and 2). Finally, the lowest three natural frequencies and associated mode shapes of a tip-massed nonprismatic EB beam with various BCs are evaluated and compared with [11].

Table 4 shows the first three nondimensional natural frequencies of a left-end tip-massed cantilever cone beam by using the proposed approach and compares the results with [11]. Using different taper ratio, left-end tip-mass, and tip-mass rotatory inertia (i.e., β₁=0, β₂=0, β₃→∞, β₄→∞, μR = δR = γR = δL =0, γL = 0 and 0.6, α= 1.2 and 2, and μL = 0.2 and 2), the accuracy of the natural frequencies competes with [11] up to 3 decimal points, with the corresponding mode shapes shown in Figure 8. Moreover, Table 5 and Figures 9–11 display the results for a cantilever cone beam with the effects of tip-mass rotatory inertia and eccentricity at the left end, (i.e., αb=αh=α=1.1, β₁=0, β₂=0, μR = δR = γR = μL =0, δL = 0.4 and 0.6, β₃ = 0.1, 1.0 and 10, and β₄= 1 and ∞), and we obtain the same accuracy level as in the previous case. Consequently, for this example of a nonprismatic EB beam with eigenvalue embedded BCs, the robustness and accuracy of the proposed approach can be verified by the numerical results that compare favorably to [11].

Table 4: The first three dimensionless natural frequencies of a cantilever cone beam with tip-mass, the rotatory inertia of mass, and its eccentricity at the left end (β₁=0, β₂=0, β₃→∞, β₄→∞, and μR = δR = γR = δL =0).

γL
α
μL

Present
(Ref. [7])
Present
(Ref. [7])
Present
(Ref. [7])

0.0
1.2
0.2
1.805113
1.805116
4.531398
4.531399
7.682832
7.682833

2.0
1.180606
1.180605
4.242896
4.242897
7.478641
7.478641

2.0
0.2
2.392504
2.392498
5.375541
5.375538
8.914098
8.914100

2.0
1.485717
1.485704
5.108411
5.108415
8.760282
8.760286

0.3
1.2
0.2
1.759466
1.759470
3.493511
3.493510
5.666955
5.666955

2.0
1.122790
1.122787
2.228089
2.228088
5.079316
5.079316

2.0
0.2
2.268671
2.268663
3.648412
3.648385
6.290178
6.290136

2.0
1.363150
1.363162
2.253308
2.253207
5.818125
5.818075

0.6
1.2
0.2
1.634321
1.634324
2.779848
2.779846
5.472033
5.472032

2.0
0.991768
0.991773
1.795564
1.795562
5.051224
5.051225

2.0
0.2
1.954460
1.954460
3.049005
3.048977
6.199510
6.199473

2.0
1.124889
1.124886
1.935342
1.935209
5.807056
5.807008

Table 5: The first three dimensionless natural frequencies of a cone beam with tip-mass, the rotatory inertia of mass, and its eccentricity at the left end (αb=αh=α=1.1, β₁=0, β₂=0, and μR = δR = γR = δL =0).

δL
β₄
γL
β₃

Present
(Ref. [7])
Present
(Ref. [7])
Present
(Ref. [7])

0.4
∞
0.6
0.1
0.467442
0.467440
1.854357
1.854354
4.358113
4.358113

1.0
0.755242
0.755247
1.948854
1.948854
4.454283
4.454283

10.0
0.934564
0.934568
2.182581
2.182581
4.834589
4.834589

1.0
0.1
0.441666
0.441668
1.583464
1.583463
4.256809
4.256810

1.0
0.703336
0.703347
1.685049
1.685046
4.358454
4.358455

10.0
0.850072
0.850076
1.921416
1.921416
4.751734
4.751734

1.0
0.6
0.1
0.461904
0.461893
1.165267
1.165266
2.480518
2.480518

1.0
0.700856
0.700855
1.197688
1.197686
2.634190
2.634190

10.0
0.801466
0.801468
1.238414
1.238412
2.891271
2.891271

0.6
∞
0.6
0.1
0.443601
0.443598
1.873050
1.873048
4.436802
4.436802

1.0
0.712292
0.712291
1.983124
1.983125
4.530416
4.530416

10.0
0.873798
0.873804
2.248240
2.248241
4.905603
4.905603

1.0
0.1
0.423210
0.423220
1.621962
1.621961
4.295470
4.295470

1.0
0.671928
0.671929
1.732488
1.732486
4.395387
4.395387

10.0
0.809967
0.809974
1.984642
1.984642
4.784972
4.784972

1.0
0.6
0.1
0.439333
0.439333
1.162506
1.162500
2.520126
2.520126

1.0
0.671066
0.671065
1.184763
1.184762
2.684893
2.684893

10.0
0.773950
0.773954
1.212627
1.212624
2.960297
2.960296

Figure 8: The first three normalized mode shapes of a cantilever cone beam with different left tip-mass ratio ().

Figure 9: The first three normalized mode shapes of a cantilever cone beam with different taper ratio of the beam (α).

Figure 10: The first three normalized mode shapes of a cantilever cone beam with different rotatory inertias () at left end.

Figure 11: The first three normalized mode shapes of a cantilever cone beam with different clamped conditions ( and ) at right end.

5.4. Modal Analysis of a Uniform Stepped EB Beam with Classical Boundary Conditions

In order to validate the proposed method for the modal analysis of a stepped beam, here we focus on the interface conditions between the two segments of the beam. For the displacement compatibility and force equilibrium conditions on the interface, one may go through Section 4.5 to handle the identical jump and scaled jump conditions with a modified jump command in the chebop system.

Several numerical examples with different classical boundary conditions and step locations are discussed. Consider a stepped beam consisting of two uniform segments S₁ and S₂. The length, width, thickness, mass per unit length, and cross-sectional moment of inertia are denoted by li, bi, hi, mi, and Ii of segment Si, respectively. Both segments have the same Young’s modulus (E) and density (ρ). Three configurations of stepped beam are presented and displayed in Figure 12:

 Configuration A: b₂ = ςb₁, h₂ = h₁, m₂/m₁ = ς and I₂/I₁ = 1 Configuration B: b₂ = b₁, h₂ = ςh₁, m₂/m₁ = ς and I₂/I₁ = ς³ Configuration C: b₂ = ςb₁, h₂ = ςh₁, m₂/m₁ = ς² and I₂/I₁ = ς⁴

Figure 12: Three cases for a stepped beam.

In this case, all the configurations of the stepped beam are compared to the results in [12]. Stepped beams with classical boundary conditions (clamped-free and simply supported) and various segment combinations are analyzed, with R₁ () denoting the length ratio of the first segment. For ς = 0.5 and R₁ = 0.25, 0.375, 0.625, and 0.75, the modal frequencies of each stepped beam are shown in Table 6. There are three types of stepped beams with two different boundary conditions and four different length ratio (R₁), and the results of the first four natural frequencies are compared with those in [12]. We can see that the accuracy of the natural frequencies directly compete with [12] up to 4 decimal points. Figures 13–15 show the first three mode shapes of the stepped beams with ς = 0.5. In more detail, Figure 13 shows the mode shapes of clamped-free stepped beams A, B, and C with R₁ = 0.375, while Figure 14 shows the mode shapes of stepped beam A with R₁ = 0.25, as well as the clamped-free and simply supported boundary conditions. Figure 15 shows the mode shapes of simply supported stepped beam A with R₁= 0.375 and 0.75. Consequently, for this two-segment EB beam example, the robustness and accuracy of the proposed approach can be verified by numerical results that directly compete with [12].

Table 6: The first four nondimensional natural frequencies for clamped-free and simply supported beams with a step ratio ς = 0.5.

BCs
Mode index
Beam type
Methods
Length ratio of the first segment (R₁/l)

0.25
0.375
0.625
0.75

Clamped-free
1
Beam-A
Present
2.0849
2.1526
2.1526
2.0849

Ref. [9]
2.0849
2.1526
2.1526
2.0849

Beam-B
Present
1.6687
1.8696
2.1193
2.0813

Ref. [9]
1.6687
1.8696
2.1193
2.0813

Beam-C
Present
1.7151
1.9777
2.3651
2.2454

Ref. [9]
1.7151
1.9777
2.3651
2.2454

2
Beam-A
Present
4.9153
4.7950
4.7950
4.9153

Ref. [9]
4.9153
4.7950
4.7950
4.9153

Beam-B
Present
3.9528
3.9387
4.0976
4.6760

Ref. [9]
3.9528
3.9387
4.0976
4.6760

Beam-C
Present
4.1395
4.1638
3.9885
4.8251

Ref. [9]
4.1395
4.1638
3.9885
4.8521

3
Beam-A
Present
7.9045
7.8343
7.8343
7.9045

Ref. [9]
7.9045
7.8343
7.8343
7.9045

Beam-B
Present
6.1635
6.1399
6.8479
6.9380

Ref. [9]
6.1635
6.1399
6.8479
6.9380

Beam-C
Present
6.3841
6.0002
7.0045
6.7812

Ref. [9]
6.8341
6.0002
7.0045
6.7812

4
Beam-A
Present
10.9711
11.0414
11.0414
10.9711

Ref. [9]
10.9711
11.0414
11.0414
10.9711

Beam-B
Present
8.3054
8.8732
9.4841
9.9945

Ref. [9]
8.3054
8.8733
9.4841
9.9945

Beam-C
Present
8.2156
8.8775
9.3342
10.1310

Ref. [9]
8.2156
8.8775
9.3342
10.1310

Simply supported
1
Beam-A
Present
3.1160
3.0890
3.0890
3.1160

Ref. [9]
3.1160
3.0890
3.0890
3.1160

Beam-B
Present
2.2270
2.2637
2.5443
2.8238

Ref. [9]
2.2270
2.2637
2.5443
2.8238

Beam-C
Present
2.1602
2.1257
2.3111
2.6160

Ref. [9]
2.1602
2.1257
2.3111
2.6160

2
Beam-A
Present
6.2218
6.2831
6.2831
6.2218

Ref. [9]
6.2218
6.2831
6.2831
6.2218

Beam-B
Present
4.5706
4.9160
5.5403
5.4805

Ref. [9]
4.5706
4.9160
5.5403
5.4805

Beam-C
Present
4.3995
4.7788
5.7147
5.4039

Ref. [9]
4.3995
4.7788
5.7147
5.4039

3
Beam-A
Present
9.4248
9.4667
9.4667
9.4248

Ref. [9]
9.4248
9.4667
9.4667
9.4248

Beam-B
Present
7.0984
7.6854
8.1049
8.6897

Ref. [9]
7.0984
7.6854
8.1049
8.6897

Beam-C
Present
6.9661
7.7817
7.9441
8.8739

Ref. [9]
6.9661
7.7817
7.9441
8.8739

4
Beam-A
Present
12.6251
12.5075
12.5075
12.6251

Ref. [9]
12.6253
12.5076
12.5076
12.6253

Beam-B
Present
9.6906
9.9389
10.8726
11.4763

Ref. [9]
9.6908
9.9390
10.8727
11.4765

Beam-C
Present
9.6683
10.0408
10.9903
11.4005

Ref. [9]
9.6684
10.0409
10.9905
11.4007

Figure 13: The first three normalized mode shapes of clamped-free stepped beams A, B, and C (ς = 0.5 and = 0.375).

Figure 14: The first three normalized mode shapes of stepped beam A with clamped-free and simply supported boundary conditions (ς = 0.5 and = 0.25).

Figure 15: The first three normalized mode shapes of simply supported stepped beam A ( = 0.375, 0.75, and ς = 0.5).

Next, as shown in Figure 16, we consider a simply supported uniform stepped EB beam with three segments; we illustrate a piezoelectric modal sensor as in [13]. The total length of the two-step beam is l, with the length and depth of each segment l₁, l₂, l₃ and h₁, h₂, h₃, respectively. Each segment has the same width, Young’s modulus and density denoted by b, E, and ρ, respectively. The boundary condition of the stepped beam is simply supported and the scaling ratio from the left to the right of each segment is set at ς₁ = h₂/h₁ = 2 and ς₂ = h₃/h₂ = 0.5. The position of the first step is located at 17l/38 while the second is located at 21l/38. The modal analysis of this two-step beam is carried out by our proposed approach. The first five natural frequencies are 3.2778, 6.2920, 9.8657, 12.6356, and 16.4596, respectively, very close to the results in [13]. Moreover, Figure 17 shows the first five normalized mode shapes of this two-step beam.

Figure 16: Uniform three stepped EB beams with S-S BCs.

Figure 17: The first five normalized mode shapes for the beam shown in Figure 16.

Compared to the simply supported prismatic EB beam, doubling the thickness (ς₁ = 2) of the second segment will increase the nondimensional fundamental natural frequency from 3.14159 to 3.2778 (about 4%). Also, this two-step beam example can be used to evaluate the influence of the natural frequency for applying adapter plate or fixture to the tunable resonant beam of a pyroshock simulator. Otherwise, the high-order mode shapes (17th-20th) plots of last two examples are provided in Supplementary Materials section (available here) for the interested reader.

6. Conclusion and Discussion

The Chebyshev spectral method with the null space approach for the modal solutions of a nonprismatic EB beam obeying general boundary and interface conditions has been presented. A unique and more general asymptotic modal solution for the higher-order modes of partially clamped EB beam has also been derived. With the chebop system in the Chebfun toolbox, the symbolic formulation for the governing equation and the boundary as well as interface conditions of an EB beam has been converted to a differential system matrix preconditioned with the null space approach. Results at different boundary or interface conditions for nonprismatic EB beams are benchmarked with those reported in the literature using alternative methods and validate the general modal solutions for EB beams with eigenvalue embedded boundary conditions, higher-order modes, and nonprismatic cross-sections. For further applications, it is recommended that tunable resonant beams for simulating pyroshock should be evaluated by the proposed approach, thereby minimizing the required trials and errors for tuning the shock spectrum knee frequency of the resonant beam.

Appendix

Higher-Order Asymptotic Modal Solutions for an EB Beam with Partially Clamped Boundary Conditions

Substituting the characteristic function in (22) into the boundary conditions in (21), with , yields the following homogeneous matrix equation (A.1) in terms of to ,

where , , , and represent sin(Λi), cos(Λi), sinh(Λi), and cosh(Λi), respectively.

The determinant of the LHS coefficient matrix is set to zero to obtain the characteristic equation, from which the nondimensional natural frequencies Λi = kil can be determined from (A.2) as follows:

where

Then, one may obtain the solutions of C₁ to C₄ with (A.1) as

Introducing C₁ to C₄ into (21), we obtain the following analytical solution of the i-th mode shape for partially clamped EB beam:

For the first 12 modes, using (A.2) and (A.3), a numerical root-finder may search the nondimensional natural frequency for the corresponding mode and obtain the mode shape by using (A.4) and (A.5). However, (A.2) and (A.5) permit the numerical evaluation for the first 12 modes (e.g., Λi ≤12π) or so; otherwise both equations are numerically ill-conditioned due to the double precision floating point arithmetic of the MATLAB environment. For Λi≥13π, an asymptotic simplification approach is used as follows:

Substituting (A.6) into (A.2) and (A.5), with the aid of symbolic computation software, we obtain the asymptotic higher-order modal solutions for the partially clamped EB beam as

which reduces (A.3) to the following P₀ to P₄ asymptotic coefficients in the characteristic equation (A.2):

With (A.2) and (A.8), one may use a numerical root-finder in searching the nondimensional natural frequency Λi for the corresponding higher-order mode and obtain the mode shape using (A.7).

It should be noted that, in Table 2 of [8], there are ten formulas for the natural frequencies and mode shapes for ten corresponding classical boundary conditions. However, by using (A.2) and (A.8) as well as feasible values for βi, this appendix provides a unique set of generalized modal solution. For the case of a simply supported beam (β₁= β₃ = 0 and β₂ = β₄ = ∞), with (A.2) and (A.8), the higher-order nondimensional natural frequency is Λi = nπ (P₂ = 0), and the corresponding mode shape is sin(kix) by (A.7).

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank Professor Michalopoulos for providing helpful suggestions that improved the presentation and content of this paper.

Supplementary Materials

File name: Supp_Mat.zip. S1: uniformeuler_review.m. Matlab M-file demonstrates eigensolution of an EB beam under generalized boundary conditions by using proposed Chebyshev spectral method and null space approach. S2: FIG 8_S-17_S.pdf. FIG 8_S-FIG_11S are plots of the higher-order mode shapes (17th-20th) in numerical example of Section 5.3. FIG 13_S-FIG_15S and FIG_17S are plots of the higher-order mode shapes (17th-20th) in numerical example of Section 5.4. (Supplementary Materials)

References

Ping Lou and Qing-yuan Zeng, “Formulation of Equations of Motion for a Simply Supported Bridge under a Moving Railway Freight Vehicle,” Shock and Vibration, vol. 14, pp. 1–18, 2007.

Ki-Young Koo and Jin-Hak Yi, “Substructural Identification of Flexural Rigidity for Beam-Like Structures,” Shock and Vibration, vol. 2015, pp. 1–15, 2015.

Q. Mao, “Free vibration analysis of elastically connected multiple-beams by using the Adomian modified decomposition method,” Journal of Sound and Vibration, vol. 331, no. 11, pp. 2532–2542, 2012.

P. F. Pai, L. Huang, S. H. Gopalakrishnamurthy, and J. H. Chung, “Identification and applications of boundary effects in beams,” International Journal of Solids and Structures, vol. 41, no. 11-12, pp. 3053–3080, 2004.

N. Olgac and N. Jalili, “Modal analysis of flexible beams with delayed resonator vibration absorber: Theory and experiments,” Journal of Sound and Vibration, vol. 218, no. 2, pp. 307–331, 1998.

E. B. Magrab, Vibrations of elastic systems: with applications to MEMS and NEMS, vol. 184, Springer Science and Business Media, 2012.

J. W. Jeong, J. H. Lim, K. W. Kim, and J. J. Lee, “Development of a Point Pyroshock Source Simulator,” Shock and Vibration, vol. 2017, pp. 1–16, 2017.

P. J. P. Gonçalves, M. J. Brennan, and S. J. Elliott, “Numerical evaluation of high-order modes of vibration in uniform Euler-Bernoulli beams,” Journal of Sound and Vibration, vol. 301, no. 3–5, pp. 1035–1039, 2007.

K. C. Cha, N. Wang, and J. Y. Liao, “Dynamics and cutting stability of the dynamically loaded worktable subjected to simply supported conditions,” The International Journal of Advanced Manufacturing Technology, vol. 71, no. 1–4, pp. 605–620, 2014.

J.-C. Hsu, H.-Y. Lai, and C.-K. Chen, “An innovative eigenvalue problem solver for free vibration of uniform Timoshenko beams by using the Adomian modified decomposition method,” Journal of Sound and Vibration, vol. 325, no. 1-2, pp. 451–470, 2009.

H.-Y. Lai, C. K. Chen, and J.-C. Hsu, “Free vibration of non-uniform Euler-Bernoulli beams by the Adomian modified decomposition method,” Computer Modeling in Engineering & Sciences, vol. 34, no. 1, pp. 87–113, 2008.

Q. Mao and S. Pietrzko, “Free vibration analysis of stepped beams by using Adomian decomposition method,” Applied Mathematics and Computation, vol. 217, no. 7, pp. 3429–3441, 2010.

Q. Mao, “Design of piezoelectric mode sensor for non-uniform Euler–Bernoulli beams with rectangular cross-section by using differential transformation method,” Mechanical Systems and Signal Processing, vol. 33, pp. 142–154, 2012.

M. Attar, “A transfer matrix method for free vibration analysis and crack identification of stepped beams with multiple edge cracks and different boundary conditions,” International Journal of Mechanical Sciences, vol. 57, no. 1, pp. 19–33, 2012.

K. Sarkar and R. Ganguli, “Modal tailoring and closed-form solutions for rotating non-uniform Euler-Bernoulli beams,” International Journal of Mechanical Sciences, vol. 88, pp. 208–220, 2014.

J. Shen, T. Tang, and L. Wang, “Applications in Multi-Dimensional Domains,” in Spectral Methods, vol. 41 of Springer Series in Computational Mathematics, pp. 367–413, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.

L. Fox and I. B. Parker, Chebyshev Polynomials in Numerical Analyis, Oxford University Press, 1st edition, 1968.

J. Lee and W. W. Schultz, “Eigenvalue analysis of Timoshenko beams and axisymmetric Mindlin plates by the pseudospectral method,” Journal of Sound and Vibration, vol. 269, no. 3-5, pp. 609–621, 2004.

M. S. Sari and E. A. Butcher, “Free vibration analysis of non-rotating and rotating Timoshenko beams with damaged boundaries using the Chebyshev collocation method,” International Journal of Mechanical Sciences, vol. 60, no. 1, pp. 1–11, 2012.

N. Wattanasakulpong and Q. Mao, “Dynamic response of Timoshenko functionally graded beams with classical and non-classical boundary conditions using Chebyshev collocation method,” Composite Structures, vol. 119, pp. 346–354, 2014.

W. Y. Wang, J. Y. Liao, and L. W. Hourng, “Matrix-based Chebyshev spectral approach to dynamic analysis of non-uniform timoshenko beams,” Structural Engineering and Mechanics, vol. 39, no. 5, pp. 669–682, 2011.

L. N. Trefethen, Spectra Methods in MATLAB, SIAM, Philadelphia, Pennsylvania, 1st edition, 2001.

T. A. Driscoll, “Automatic spectral collocation for integral, integro-differential, and integrally reformulated differential equations,” Journal of Computational Physics, vol. 229, no. 17, pp. 5980–5998, 2010.

“The Chebfun Team, University of Oxford,” http://www.chebfun.org/.

T. A. Driscoll, F. Bornemann, and L. N. Trefethen, “The chebop system for automatic solution of differential equations,” BIT Numerical Mathematics, vol. 48, no. 4, pp. 701–723, 2008.

T. A. Driscoll and N. Hale, “Rectangular spectral collocation,” IMA Journal of Numerical Analysis (IMAJNA), vol. 36, no. 1, pp. 108–132, 2016.

J. L. Aurentz and L. N. Trefethen, “Block operators and spectral discretizations,” SIAM Review, vol. 59, no. 2, pp. 423–446, 2017.

B. Smith, R. Laoulache, and A. Heryudono, “Implementation of Neumann boundary condition with influence matrix method for viscous annular flow using pseudospectral collocation,” Journal of Computational and Applied Mathematics, vol. 285, pp. 100–115, 2015.

J. H. Wilkinson, The Algebraic Problem, Oxford University Press, Oxford, UK, 1965.

Y. Saad, Numerical Methods for Large Eigenvalue Problems, Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, USA, 2nd edition, 2011.

P.-y. Nie, “A null space method for solving system of equations,” Applied Mathematics and Computation, vol. 149, no. 1, pp. 215–226, 2004.

C. L. Mitsas, “Data reconciliation and variable classification by null space methods,” Measurement, vol. 43, no. 5, pp. 702–707, 2010.

C. Felippa and K. Park, “The construction of free–free flexibility matrices for multilevel structural analysis,” Computer Methods Applied Mechanics and Engineering, vol. 191, no. 19-20, pp. 2139–2168, 2002.