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

Authors: C. P. Hsu; C. F. Hung; J. Y. Liao
Journal: Shock and Vibration
 (2018)

---

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


---

## Body

## 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 thechebop 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 Figure1, whose length isl 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(1)∂2∂x2EIx∂2yx,t∂x2+ρAx∂2yx,t∂t2=0,wherey(x,t) denotes the transverse deflection of the beam,A(x) is the cross-sectional area atx,I(x) is the area moment of inertia ofA(x), and ρ andE 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(2)yx,t=ϕxe-iωt,(1) becomes(3)d2dx2EIxd2ϕxdx2-ω2ρAxϕx=0.The corresponding external boundary conditions at both ends are given by(4)EIxd2ϕxdx2+JL+MLeL2ω2dϕxdx-kRLdϕxdx-MLeLω2ϕx=0,(5)ddxEId2ϕxdx2+MLeLω2dϕxdx+kTLϕx-MLω2ϕx=0,atx= 0, and(6)EId2ϕxdx2+kRRdϕxdx-JR+MReR2ω2dϕxdx-MReRω2ϕx=0,(7)ddxEId2ϕxdx2+MReRω2dϕxdx-kTRϕx+MRω2ϕx=0,atx=l.Here,kTR,kRR,MR,JR,eR andkTL,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:(8)X=xl;ΦX=ϕxl;Ω4=ρA0ω2l4EI0;μmL=MLρA0l;μmR=MRρA0l;KRL=kRLlEI0;KTL=kTLl3EI0;KRR=kRRlEIl;KTR=kTRl3EIl;γL=JLMLl2;γR=JRMRl2;δL=eLl;δR=eRl;iX=IxI0;aX=AxA0,where the range ofXis [0,1], andA0 andI0 are the cross-sectional area and the moment of inertia atx = 0, respectively.Substituting (8) into (3), the nondimensional form of (3) is(9)d2dX2iXd2ΦXdX2-Ω4aXΦX=0.Using (8) in (4) to (7), the boundary conditions at both ends in nondimensional form become as follows:(10)Φ′′0+μmLδL2+γL2Ω4-KRLΦ′0-μmLδLΩ4Φ0=0,(11)Φ′′′0-i′0Φ′′0+μmLδLΩ4Φ′0+KTL-μmLΩ4Φ0=0,(12)Φ′′1+KRR-μmRδR2+γR2Ω4Φ′1-μmRδRΩ4Φ1=0,(13)Φ′′′1-i′1i1Φ′′1+μmRδRΩ4Φ′1-KTR-μmRΩ4Φ1=0,where(14)Φ′X=dΦXdX,Φ′′X=d2ΦXdX2,Φ′′′X=d3ΦXdX3,i′X=diXdX.


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

Consider the free vibration of a straight EB beam consisting ofN uniform sections elastically restrained at both ends, as shown in Figure 2. The total length isl, and the ratio of the distance from the right end of theSj section to the left end of the beam is Ri(i=1,2,…,N). Let us define the junction parameters ηj, γj, and μj for the beam sectionsjandj+1 (j= 1, 2, …,N-1) as follows:(15)ηj=Ej+1Ij+1EjIj,γi=ρj+1Aj+1ρjAj,μj=γjηj,forj=1,2,⋯,N-1.HereEj, ρj,Aj, andIj are Young’s modulus, density per length, the cross-section area, and the inertia of cross-section area ofj-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(16)d4ΦiXidXi4-Ωi4ΦiXi=0,Xi∈0,Ri,i=1,2,⋯,N,where Ωi2=ωρiAil4/EiIi,andΩj2=μjΩj+12,j=1,…,N-1.Figure 2
A stepped beam with elastically restrained ends.The interface conditions at the junction of each segment are(17)identical jump:ΦiRj=Φi+1Rj,Φi′Rj=Φi+1′Rj,scaled jump:Φi′′Rj=ηjΦi+1′′Rj,Φi′′′Rj=ηjΦi+1′′′Rj,and the boundary conditions at the both ends are(18)Φ1′′0=KRLΦ1′0,Φ1′′′0=-KTLΦ10,ηN-1ΦN′′1=-KRRΦN′1,ηN-1ΦN′′′1=KTRΦN1.


### 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(19)Mx=0=EIϕ′′x=0=kθLϕ′x=0,Vx=0=EIϕ′′′x=0=-kTLϕx=0,Mx=l=EIϕ′′x=l=-kθRϕ′x=l,Vx=l=EIϕ′′′x=l=kTRϕx=l.Specifically, ϕ is written as a finite sum, the so-called Galerkin approximation:(20)ϕx=∑i=1nqiψix.The assumed mode expansion of ψi(x) can be written in terms of trigonometric and hyperbolic functions as(21)ψix=C1sinκix+C2cos⁡κix+C3sinh⁡κix+C4cosh⁡κix,κi2=ωiρAEI.Let us define the following nondimensional spring constants,(22)β1=kθLlEI,β2=kTLl3EI,β3=kθRlEI,β4=kTRl3EI.The boundary conditions in (19) can be given by the following nondimensional forms.(23)ψi′′0-β1lψi′0=0,ψi′′′0+β2l3ψi0=0,ψi′′l+β3lψi′l=0,ψi′′′l-β4l3ψil=0.According to descriptions found in the Appendix, the analytical modal solution of the beam with partially clamped boundary conditions can be described by(24)ψi, exactx=sin⁡kixC-CHΛi5+S-SHβ3Λi4+2β2SHΛi2+β1CH+β1C+2β3CHβ2Λi+SH+Sβ1β2β3+cos⁡kixSH-SΛi5+2β1CH+β3CH+β3CΛi4+2β2β3SHΛi3-Sh+Sβ1β2Λi-CH+Cβ1β2β3+sinhkixC-CHΛi5+S-SHβ3Λi4+2β2SΛi2-β1C+β1CH+2β3Cβ2Λi-S+SHβ1β2β3+coshkixSH-SΛi5+2β1C+β3CH+β3CΛi4+2β1β3SΛi3+SH+Sβ1β2Λi+CH-Cβ1β2β3,whereΛi =k(ωi)l=kilis a nondimensional natural frequency,S,C,SH, andCH 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:(25)ψi, asympx=sinkix-Λi5-Λi4β3+2β2Λi2+β1β2+2β3β2Λi+β3β2β1+coskixΛi5+β3+2β1Λi4+2β1β3Λi3-β1β2Λi-β3β2β1+e-kixΛi5+Λi4β3+β1β2Λi+β3β2β1+e-kil-xC-SΛi5+β3S+C+2β1CΛi4+2β1β3SΛi3+β2SΛi2+β2β1S-C-2β3CΛi-β1β2β3S+C.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).


## 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 Figure1, whose length isl 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(1)∂2∂x2EIx∂2yx,t∂x2+ρAx∂2yx,t∂t2=0,wherey(x,t) denotes the transverse deflection of the beam,A(x) is the cross-sectional area atx,I(x) is the area moment of inertia ofA(x), and ρ andE 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(2)yx,t=ϕxe-iωt,(1) becomes(3)d2dx2EIxd2ϕxdx2-ω2ρAxϕx=0.The corresponding external boundary conditions at both ends are given by(4)EIxd2ϕxdx2+JL+MLeL2ω2dϕxdx-kRLdϕxdx-MLeLω2ϕx=0,(5)ddxEId2ϕxdx2+MLeLω2dϕxdx+kTLϕx-MLω2ϕx=0,atx= 0, and(6)EId2ϕxdx2+kRRdϕxdx-JR+MReR2ω2dϕxdx-MReRω2ϕx=0,(7)ddxEId2ϕxdx2+MReRω2dϕxdx-kTRϕx+MRω2ϕx=0,atx=l.Here,kTR,kRR,MR,JR,eR andkTL,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:(8)X=xl;ΦX=ϕxl;Ω4=ρA0ω2l4EI0;μmL=MLρA0l;μmR=MRρA0l;KRL=kRLlEI0;KTL=kTLl3EI0;KRR=kRRlEIl;KTR=kTRl3EIl;γL=JLMLl2;γR=JRMRl2;δL=eLl;δR=eRl;iX=IxI0;aX=AxA0,where the range ofXis [0,1], andA0 andI0 are the cross-sectional area and the moment of inertia atx = 0, respectively.Substituting (8) into (3), the nondimensional form of (3) is(9)d2dX2iXd2ΦXdX2-Ω4aXΦX=0.Using (8) in (4) to (7), the boundary conditions at both ends in nondimensional form become as follows:(10)Φ′′0+μmLδL2+γL2Ω4-KRLΦ′0-μmLδLΩ4Φ0=0,(11)Φ′′′0-i′0Φ′′0+μmLδLΩ4Φ′0+KTL-μmLΩ4Φ0=0,(12)Φ′′1+KRR-μmRδR2+γR2Ω4Φ′1-μmRδRΩ4Φ1=0,(13)Φ′′′1-i′1i1Φ′′1+μmRδRΩ4Φ′1-KTR-μmRΩ4Φ1=0,where(14)Φ′X=dΦXdX,Φ′′X=d2ΦXdX2,Φ′′′X=d3ΦXdX3,i′X=diXdX.


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

Consider the free vibration of a straight EB beam consisting ofN uniform sections elastically restrained at both ends, as shown in Figure 2. The total length isl, and the ratio of the distance from the right end of theSj section to the left end of the beam is Ri(i=1,2,…,N). Let us define the junction parameters ηj, γj, and μj for the beam sectionsjandj+1 (j= 1, 2, …,N-1) as follows:(15)ηj=Ej+1Ij+1EjIj,γi=ρj+1Aj+1ρjAj,μj=γjηj,forj=1,2,⋯,N-1.HereEj, ρj,Aj, andIj are Young’s modulus, density per length, the cross-section area, and the inertia of cross-section area ofj-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(16)d4ΦiXidXi4-Ωi4ΦiXi=0,Xi∈0,Ri,i=1,2,⋯,N,where Ωi2=ωρiAil4/EiIi,andΩj2=μjΩj+12,j=1,…,N-1.Figure 2
A stepped beam with elastically restrained ends.The interface conditions at the junction of each segment are(17)identical jump:ΦiRj=Φi+1Rj,Φi′Rj=Φi+1′Rj,scaled jump:Φi′′Rj=ηjΦi+1′′Rj,Φi′′′Rj=ηjΦi+1′′′Rj,and the boundary conditions at the both ends are(18)Φ1′′0=KRLΦ1′0,Φ1′′′0=-KTLΦ10,ηN-1ΦN′′1=-KRRΦN′1,ηN-1ΦN′′′1=KTRΦN1.


## 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(19)Mx=0=EIϕ′′x=0=kθLϕ′x=0,Vx=0=EIϕ′′′x=0=-kTLϕx=0,Mx=l=EIϕ′′x=l=-kθRϕ′x=l,Vx=l=EIϕ′′′x=l=kTRϕx=l.Specifically, ϕ is written as a finite sum, the so-called Galerkin approximation:(20)ϕx=∑i=1nqiψix.The assumed mode expansion of ψi(x) can be written in terms of trigonometric and hyperbolic functions as(21)ψix=C1sinκix+C2cos⁡κix+C3sinh⁡κix+C4cosh⁡κix,κi2=ωiρAEI.Let us define the following nondimensional spring constants,(22)β1=kθLlEI,β2=kTLl3EI,β3=kθRlEI,β4=kTRl3EI.The boundary conditions in (19) can be given by the following nondimensional forms.(23)ψi′′0-β1lψi′0=0,ψi′′′0+β2l3ψi0=0,ψi′′l+β3lψi′l=0,ψi′′′l-β4l3ψil=0.According to descriptions found in the Appendix, the analytical modal solution of the beam with partially clamped boundary conditions can be described by(24)ψi, exactx=sin⁡kixC-CHΛi5+S-SHβ3Λi4+2β2SHΛi2+β1CH+β1C+2β3CHβ2Λi+SH+Sβ1β2β3+cos⁡kixSH-SΛi5+2β1CH+β3CH+β3CΛi4+2β2β3SHΛi3-Sh+Sβ1β2Λi-CH+Cβ1β2β3+sinhkixC-CHΛi5+S-SHβ3Λi4+2β2SΛi2-β1C+β1CH+2β3Cβ2Λi-S+SHβ1β2β3+coshkixSH-SΛi5+2β1C+β3CH+β3CΛi4+2β1β3SΛi3+SH+Sβ1β2Λi+CH-Cβ1β2β3,whereΛi =k(ωi)l=kilis a nondimensional natural frequency,S,C,SH, andCH 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:(25)ψi, asympx=sinkix-Λi5-Λi4β3+2β2Λi2+β1β2+2β3β2Λi+β3β2β1+coskixΛi5+β3+2β1Λi4+2β1β3Λi3-β1β2Λi-β3β2β1+e-kixΛi5+Λi4β3+β1β2Λi+β3β2β1+e-kil-xC-SΛi5+β3S+C+2β1CΛi4+2β1β3SΛi3+β2SΛi2+β2β1S-C-2β3CΛi-β1β2β3S+C.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], thekth-order Chebyshev polynomialTN(s) and Chebyshev-Gauss-Lobatto (CGL) collocation pointssj, in the domain of sj∈[-1,1], can be expressed, respectively, by(26)TNsj=cos⁡Ncos-1sj,sj=cos⁡jπN.The derivative of a smooth g(x) function at the CGL pointssj can then be computed via matrix-vector multiplication, which can be formally represented as(27)g′s0g′s1⋮g′sN=DNgs0gs1⋮gsN,where DN is a (N+1)×(N+1) matrix.The elements within the matrix can be expressed as(28)DN00=2N2+16,DNNN=-DN00,DNjj=-sj2sin2jπ/N,j=1,…,N-1,DNij=-1i+j+1ci2cjsin2π/Ni+jsin2π/Ni-j,i≠j,i,j=0,…,N-1,cj=2,j=0orN1,j=1,…,N-1.Moreover, its higher derivative matrices may be expressed as(29)DNmij=mDNm-1iiDNij-si-sj-1DNm-1ij,i≠j,DNmii=-∑j=0,j≠iNDNmij.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:(30)Ax-1JDN4+2J′DN3+J′′DN2Φ0Φ1⋮ΦN=λΦ0Φ1⋮ΦN,where Xi∈[0,1] is a transformation of CGL collocation pointssj,Xi = (si +1)/2, and(31)J=diagiX0iX1⋮iXN,J′=diagi′X0i′X1⋮i′XN,J′′=diagi′′X0i′′X1⋮i′′XN,Ax=diagaX0aX1⋮aXN,Φi=ΦXi,λ=Ω4.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:(32)DN×NXN×1=λXN×1,where XN×1 denotes the CGL collocation vector and DN×N is the system matrix.VectorX can be decomposed as {Xc}, {Xλ}, and {Xr}. {Xc} is the collection of homogeneous boundary and interface conditions. {Xλ} is the matrix of boundary and interface conditions containing the eigenvalues. {Xr} is the residual matrix of the governing equation of an EB beam. Then, (32) can be partitioned as(33)DrrDrcDrλDcrDccDcλDλrDλcDλλN×NXrXcXλN×1=λXrXcXλN×1.The homogeneous boundary conditions should satisfy the following equation:(34)CcrCccCcλX=0·Xc.When boundary conditions include the eigenvalues, such as in (10) to (13), the solutions should satisfy the following equation:(35)CλrCλcCλλX=λ·Xλ.Substituting (34) and (35) into (33), by the row replacement approach we obtain the following expression:(36)DrrDrcDrλCcrCccCcλDλr-CλrDλc-CλcDλλ-CλλXrXcXλ=λXr00.Simplifying (36) by reducing the order of the matrix, we obtain the final expression with eigenvalues of boundary conditions:(37)Drr-RrcRcc-1RcrXr=λ·Xr,where(38)Rcc=CccCcλDλc-CλcDλλ-Cλλ,Rcr=CcrDλr-Cλr,Rcλ=RrcRcc,Rrc=DrcDrλ.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 largeN. 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 forN+m columns by barycentric resampling discretization and expanded to a rectangular matrix (N× [N+m], wherem is the number of boundary conditions).(b)
Them sets of boundary conditions formed by CSM are then placed above the rectangular matrix and reconstruct a square DN+m×N+m 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 Ds∈RN×N,X∈RN×1, and matrix Cs∈Rm×N represents the imposed boundary and interface conditions. The rank of matrices [Ds] and [Cs] isN andm (N>m), respectively. The eigenproblem of (33) must satisfy the requirement, [Cs]{X}=  {0}. Let the null space of [Cs] be matrix [Nc] and satisfy the identity [Cs][Nc] = [0]. It should be noted that {X} is a vector in physical space, and it may be projected to vector {ξ} in the null space with matrix [Nc]; i.e., {X} = [Nc]{ξ}. As a result, (33) may be projected to the null space of [Cs] as follows:(39)DsNcξ=λNcξ.Premultiplying both sides of (39) with the transposed [Nc], and taking advantage of the orthonormal property in [Nc], we get(40)NcTDsNcξ=λNcTNcξ=λξ,Dsνξ=λξ,Dsν=NcTDsNc.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 {Xλ} = [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 theChebfun Toolbox

In this section, we discuss the solution roadmap by using theChebfun 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, Figure3 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 theChebfun 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: (41)Ω4ΦX=F0X,DΦX,where F0X,D=i′′XD2+2i′XD3+iXD4,DkΦX=dkΦX/dxk.The derivative of (41) with respect toX is(42)Ω4Φ′X=F1X,DΦX,where F1X,D=i′′′XD2+3i′′XD3+3i′XD4+iXD5.For boundary conditions (10)-(13), all the terms containing Ω4 can be removed by substituting the RHSs of (41) and (42), together with a simple manipulation to get(43)Φ′′0-KRLΦ′0+μmLδL2+γL2F10,D-δLF00,DΦ0=0,(44)Φ′′′0-i′0Φ′′0+μmLδLF10,D-F00,D+KTLΦ0=0,(45)Φ′′1+KRRΦ′1-μmRδR2+γR2F11,D+δRF01,DΦ1=0,(46)Φ′′′1-i′1i1Φ′′1+μmRδRF11,D+F01,D-KTRΦ1=0.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 usingN-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 theRj ratio at each segment, (16) may be transformed into a Chebfun object by multiplied with the Heaviside function:(47)1+μj-1·HeavisideXRj-1D4Φ=Ω4Φ.On the other hand, the LHS of (18) can be expressed using a differential operator notation as(48)D2Φ0-KRLDΦ0=0,D3Φ0+KTLΦ0=0,while the LHS of (18) can be rewritten using the differential operator notation by(49)ηN-1D2Φ1+KRRDΦ1=0,ηN-1D3Φ1-KTRΦ1=0.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(Φ,Rj); jump(DΦ,Rj); jump(D2Φ,Rj); jump(D3Φ,Rj)].Since thejump 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 Figure4. 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.


## 4.1. Row Replacement and Rectangular Spectral Collocation Approaches.

From (9) to (13), one may obtain a generalized eigenproblem as follows:(32)DN×NXN×1=λXN×1,where XN×1 denotes the CGL collocation vector and DN×N is the system matrix.VectorX can be decomposed as {Xc}, {Xλ}, and {Xr}. {Xc} is the collection of homogeneous boundary and interface conditions. {Xλ} is the matrix of boundary and interface conditions containing the eigenvalues. {Xr} is the residual matrix of the governing equation of an EB beam. Then, (32) can be partitioned as(33)DrrDrcDrλDcrDccDcλDλrDλcDλλN×NXrXcXλN×1=λXrXcXλN×1.The homogeneous boundary conditions should satisfy the following equation:(34)CcrCccCcλX=0·Xc.When boundary conditions include the eigenvalues, such as in (10) to (13), the solutions should satisfy the following equation:(35)CλrCλcCλλX=λ·Xλ.Substituting (34) and (35) into (33), by the row replacement approach we obtain the following expression:(36)DrrDrcDrλCcrCccCcλDλr-CλrDλc-CλcDλλ-CλλXrXcXλ=λXr00.Simplifying (36) by reducing the order of the matrix, we obtain the final expression with eigenvalues of boundary conditions:(37)Drr-RrcRcc-1RcrXr=λ·Xr,where(38)Rcc=CccCcλDλc-CλcDλλ-Cλλ,Rcr=CcrDλr-Cλr,Rcλ=RrcRcc,Rrc=DrcDrλ.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 largeN. 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 forN+m columns by barycentric resampling discretization and expanded to a rectangular matrix (N× [N+m], wherem is the number of boundary conditions).(b)
Them sets of boundary conditions formed by CSM are then placed above the rectangular matrix and reconstruct a square DN+m×N+m 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 Ds∈RN×N,X∈RN×1, and matrix Cs∈Rm×N represents the imposed boundary and interface conditions. The rank of matrices [Ds] and [Cs] isN andm (N>m), respectively. The eigenproblem of (33) must satisfy the requirement, [Cs]{X}=  {0}. Let the null space of [Cs] be matrix [Nc] and satisfy the identity [Cs][Nc] = [0]. It should be noted that {X} is a vector in physical space, and it may be projected to vector {ξ} in the null space with matrix [Nc]; i.e., {X} = [Nc]{ξ}. As a result, (33) may be projected to the null space of [Cs] as follows:(39)DsNcξ=λNcξ.Premultiplying both sides of (39) with the transposed [Nc], and taking advantage of the orthonormal property in [Nc], we get(40)NcTDsNcξ=λNcTNcξ=λξ,Dsνξ=λξ,Dsν=NcTDsNc.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 {Xλ} = [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 theChebfun Toolbox

In this section, we discuss the solution roadmap by using theChebfun 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, Figure3 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 theChebfun 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: (41)Ω4ΦX=F0X,DΦX,where F0X,D=i′′XD2+2i′XD3+iXD4,DkΦX=dkΦX/dxk.The derivative of (41) with respect toX is(42)Ω4Φ′X=F1X,DΦX,where F1X,D=i′′′XD2+3i′′XD3+3i′XD4+iXD5.For boundary conditions (10)-(13), all the terms containing Ω4 can be removed by substituting the RHSs of (41) and (42), together with a simple manipulation to get(43)Φ′′0-KRLΦ′0+μmLδL2+γL2F10,D-δLF00,DΦ0=0,(44)Φ′′′0-i′0Φ′′0+μmLδLF10,D-F00,D+KTLΦ0=0,(45)Φ′′1+KRRΦ′1-μmRδR2+γR2F11,D+δRF01,DΦ1=0,(46)Φ′′′1-i′1i1Φ′′1+μmRδRF11,D+F01,D-KTRΦ1=0.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 usingN-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 theRj ratio at each segment, (16) may be transformed into a Chebfun object by multiplied with the Heaviside function:(47)1+μj-1·HeavisideXRj-1D4Φ=Ω4Φ.On the other hand, the LHS of (18) can be expressed using a differential operator notation as(48)D2Φ0-KRLDΦ0=0,D3Φ0+KTLΦ0=0,while the LHS of (18) can be rewritten using the differential operator notation by(49)ηN-1D2Φ1+KRRDΦ1=0,ηN-1D3Φ1-KTRΦ1=0.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(Φ,Rj); jump(DΦ,Rj); jump(D2Φ,Rj); jump(D3Φ,Rj)].Since thejump 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 Figure4. 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 andChebfun 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 arenπ, and the corresponding mode shapes are sin(nπ), respectively, wheren is the mode order.With theChebfun 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.Table1 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 areN=30 andN=200, respectively. However, in accordance with the sampling theorem, the maximum number of modes (nM) should not be larger thanN/2 in order to get rid of aliasing. In Table 1, we display the cases ofnM =10 andnM =100 only. As the results indicate, the accuracy of the present approach outperforms the RSC approach in both cases. However, for the RSC approach withN=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:(50)ifnM<11,N=3nM;else,N=2nM.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: (51)ERSC=maxn⁡1-∫0lϕn2xRSCdx∫0lϕn2xexactdx,ENSA=maxn⁡1-∫0lϕn2xNSAdx∫0lϕn2xexactdx,n=1,2,…,5.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 chosennM.Table 1
Comparison of dimensionless natural frequencies fornM =10 and nM =100 by exact solutions, RSC, and present methods.
CGL points n M Method Ω 1 Ω 2 Ω 3 … Ω n M - 2 Ω n M - 1 Ω n M 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.15927Table 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-08Figure 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β1, β2, β3, and β4 in (22). The higher-order asymptotic modal solution, obtained in (25), is denoted by PC1. Moreover, PC2 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 PC0, 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 108, respectively. PC1 is substituted in combination with the asymptotic free and fixed βi into (25). In Figure 3, PC2 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 (PC1) and the NSA approach (PC2). Then, these results are compared to the asymptotic modal solutions as proposed in [8] (PC0). The normalized maxi-max deviations of the natural frequency between PC1, PC2, and PC0 are defined by(52)ePC1=maxn⁡1-Ωn,PC1Ωn,PC0,ePC2=maxn⁡1-Ωn,PC2Ωn,PC0,and the normalized maxi-max deviations of the mode shape between PC1, PC2, and PC0 are defined by(53)EPC1=maxn⁡1-∫0lϕn2xPC1dx∫0lϕn2xPC0dx,EPC2=maxn⁡1-∫0lϕn2xPC2dx∫0lϕn2xPC0dx.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 PC1 and PC2 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 PC1 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 PC1 and PC2 with four classical BCs imposed.
Boundary condition Mode order Maxi-max deviations of natural frequency Maxi-max deviations of mode shape PC1 PC2 PC1 PC2 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.075393Figure 6
Normalized maxi-max deviations of natural frequency by using PC1 and PC2 (ePC1 and ePC2).Figure 7
Normalized maxi-max deviations of mode shape by using PC1 and PC2 (EPC1 and EPC2).


### 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 (β3 →∞, β4 →∞), 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].Table4 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., β1=0, β2=0, β3→∞, β4→∞, μ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, β1=0, β2=0, μR = δR = γR =μL =0, δL = 0.4 and 0.6, β3 = 0.1, 1.0 and 10, and β4= 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 (β1=0, β2=0, β3→∞, β4→∞, and μR = δR = γR = δL =0).
γL α μL Ω 1 Ω 2 Ω 3 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.807008Table 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, β1=0, β2=0, and μR = δR = γR = δL =0).
δL β4 γL β3 Ω 1 Ω 2 Ω 3 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.960296Figure 8
The first three normalized mode shapes of a cantilever cone beam with different left tip-mass ratio (μL).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 (γL) at left end.Figure 11
The first three normalized mode shapes of a cantilever cone beam with different clamped conditions (β3 and β4) 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 Section4.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 segmentsS1 andS2. The length, width, thickness, mass per unit length, and cross-sectional moment of inertia are denoted byli,bi,hi,mi, andIi of segmentSi, respectively. Both segments have the same Young’s modulus (E) and density (ρ). Three configurations of stepped beam are presented and displayed in Figure 12:ConfigurationA:b2 =ςb1,h2 =h1,m2/m1 = ς andI2/I1 = 1Configuration B:b2 =b1,h2 =ςh1,m2/m1 = ς andI2/I1 = ς3Configuration C:b2 =ςb1,h2 =ςh1,m2/m1 = ς2 andI2/I1 = ς4Figure 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, withR1 (R1=l1/l1+l2) denoting the length ratio of the first segment. For ς = 0.5 andR1 = 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 (R1), 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 withR1 = 0.375, while Figure 14 shows the mode shapes of stepped beam A withR1 = 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 withR1= 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 (R1/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.4007Figure 13
The first three normalized mode shapes of clamped-free stepped beams A, B, and C (ς= 0.5 and R1 = 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 R1 = 0.25).Figure 15
The first three normalized mode shapes of simply supported stepped beam A (R1 = 0.375, 0.75, andς= 0.5).Next, as shown in Figure16, 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 isl, with the length and depth of each segmentl1,l2,l3 andh1,h2,h3, respectively. Each segment has the same width, Young’s modulus and density denoted byb,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 ς1 =h2/h1 = 2 and ς2 =h3/h2 = 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 Figure16.Compared to the simply supported prismatic EB beam, doubling the thickness (ς1 = 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.


## 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 arenπ, and the corresponding mode shapes are sin(nπ), respectively, wheren is the mode order.With theChebfun 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.Table1 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 areN=30 andN=200, respectively. However, in accordance with the sampling theorem, the maximum number of modes (nM) should not be larger thanN/2 in order to get rid of aliasing. In Table 1, we display the cases ofnM =10 andnM =100 only. As the results indicate, the accuracy of the present approach outperforms the RSC approach in both cases. However, for the RSC approach withN=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:(50)ifnM<11,N=3nM;else,N=2nM.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: (51)ERSC=maxn⁡1-∫0lϕn2xRSCdx∫0lϕn2xexactdx,ENSA=maxn⁡1-∫0lϕn2xNSAdx∫0lϕn2xexactdx,n=1,2,…,5.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 chosennM.Table 1
Comparison of dimensionless natural frequencies fornM =10 and nM =100 by exact solutions, RSC, and present methods.
CGL points n M Method Ω 1 Ω 2 Ω 3 … Ω n M - 2 Ω n M - 1 Ω n M 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.15927Table 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-08Figure 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β1, β2, β3, and β4 in (22). The higher-order asymptotic modal solution, obtained in (25), is denoted by PC1. Moreover, PC2 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 PC0, 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 108, respectively. PC1 is substituted in combination with the asymptotic free and fixed βi into (25). In Figure 3, PC2 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 (PC1) and the NSA approach (PC2). Then, these results are compared to the asymptotic modal solutions as proposed in [8] (PC0). The normalized maxi-max deviations of the natural frequency between PC1, PC2, and PC0 are defined by(52)ePC1=maxn⁡1-Ωn,PC1Ωn,PC0,ePC2=maxn⁡1-Ωn,PC2Ωn,PC0,and the normalized maxi-max deviations of the mode shape between PC1, PC2, and PC0 are defined by(53)EPC1=maxn⁡1-∫0lϕn2xPC1dx∫0lϕn2xPC0dx,EPC2=maxn⁡1-∫0lϕn2xPC2dx∫0lϕn2xPC0dx.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 PC1 and PC2 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 PC1 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 PC1 and PC2 with four classical BCs imposed.
Boundary condition Mode order Maxi-max deviations of natural frequency Maxi-max deviations of mode shape PC1 PC2 PC1 PC2 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.075393Figure 6
Normalized maxi-max deviations of natural frequency by using PC1 and PC2 (ePC1 and ePC2).Figure 7
Normalized maxi-max deviations of mode shape by using PC1 and PC2 (EPC1 and EPC2).


## 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 (β3 →∞, β4 →∞), 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].Table4 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., β1=0, β2=0, β3→∞, β4→∞, μ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, β1=0, β2=0, μR = δR = γR =μL =0, δL = 0.4 and 0.6, β3 = 0.1, 1.0 and 10, and β4= 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 (β1=0, β2=0, β3→∞, β4→∞, and μR = δR = γR = δL =0).
γL α μL Ω 1 Ω 2 Ω 3 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.807008Table 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, β1=0, β2=0, and μR = δR = γR = δL =0).
δL β4 γL β3 Ω 1 Ω 2 Ω 3 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.960296Figure 8
The first three normalized mode shapes of a cantilever cone beam with different left tip-mass ratio (μL).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 (γL) at left end.Figure 11
The first three normalized mode shapes of a cantilever cone beam with different clamped conditions (β3 and β4) 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 Section4.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 segmentsS1 andS2. The length, width, thickness, mass per unit length, and cross-sectional moment of inertia are denoted byli,bi,hi,mi, andIi of segmentSi, respectively. Both segments have the same Young’s modulus (E) and density (ρ). Three configurations of stepped beam are presented and displayed in Figure 12:ConfigurationA:b2 =ςb1,h2 =h1,m2/m1 = ς andI2/I1 = 1Configuration B:b2 =b1,h2 =ςh1,m2/m1 = ς andI2/I1 = ς3Configuration C:b2 =ςb1,h2 =ςh1,m2/m1 = ς2 andI2/I1 = ς4Figure 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, withR1 (R1=l1/l1+l2) denoting the length ratio of the first segment. For ς = 0.5 andR1 = 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 (R1), 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 withR1 = 0.375, while Figure 14 shows the mode shapes of stepped beam A withR1 = 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 withR1= 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 (R1/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.4007Figure 13
The first three normalized mode shapes of clamped-free stepped beams A, B, and C (ς= 0.5 and R1 = 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 R1 = 0.25).Figure 15
The first three normalized mode shapes of simply supported stepped beam A (R1 = 0.375, 0.75, andς= 0.5).Next, as shown in Figure16, 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 isl, with the length and depth of each segmentl1,l2,l3 andh1,h2,h3, respectively. Each segment has the same width, Young’s modulus and density denoted byb,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 ς1 =h2/h1 = 2 and ς2 =h3/h2 = 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 Figure16.Compared to the simply supported prismatic EB beam, doubling the thickness (ς1 = 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 thechebop 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.



---

*Source: 2487697-2018-07-16.xml*