# Reference Value Selection in a Perturbation Theory Applied to Nonuniform Beams

Authors: Blake Martin; Armaghan Salehian
Journal: Shock and Vibration
 (2018)

---

## Abstract

The Lindstedt-Poincaré method is applied to a nonuniform Euler-Bernoulli beam model for the free transverse vibrations of the system. The nonuniformities in the system include spatially varying and piecewise continuous bending stiffness and mass per unit length. The expression for the natural frequencies is obtained up to second-order and the expression for the mode shapes is obtained up to first-order. The explicit dependence of the natural frequencies and mode shapes on reference values for the bending stiffness and the mass per unit length of the system is determined. Multiple methods for choosing these reference values are presented and are compared using numerical examples.


---

## Body

## 1. Introduction

Nonuniformities in a mechanical system can be caused by sudden or gradual changes in geometry or material properties or a combination of multiple components into the system. An example for these nonuniformities includes the addition of components such as electronic cables and power cords used commonly for control of space structures. A simple way to account for the dynamic effects of these additional components is to use ad hoc models that account for the additional mass and stiffness of the cables on the host structure. Previous works by the authors include development of more complex mathematical techniques to include the local stiffening effects of these additive components using homogenization techniques [1–6]. The main underlying assumption behind these works includes the geometric periodicity and a repeated pattern for the cables harnessing the host structure, which is necessary for employing a homogenization technique. In an attempt to improve these mathematical models, the assumption of periodicity of the cable pattern should be corrected to obtain a general solution for systems with additive components that do not obey a repetitive pattern. Therefore, an ultimate goal is to study the dynamic effects of these nonperiodic additive components on a host space structure. Since these additional cables are treated as mechanical elements with stiffness and mass, their local dynamic effects may be approximated by variable mass per unit length and stiffness specific to the point of attachment to the host structure. As a result, differential equations with spatially variable coefficients will be developed in this work for which a perturbation theory or other approximation modelling techniques should be used in order to obtain a solution. Therefore, the current work entails the next attempt in making this modelling a possibility by considering structures that are modelled as nonuniform beams.Some of the previous works on nonuniform structures consider analytical modelling techniques. For example, [7, 8] considered a beam with an exponentially varying width. It is shown that for this type of geometry the spatial ordinary differential equation (ODE), which is used to determine the natural frequencies and mode shapes, reduces to one with constant coefficients. In [9] it was demonstrated that when the area and moment of inertia of a beam structure are of a specific polynomial form, the equation of motion for the transverse bending could be transformed to that of a homogeneous beam for which an analytical solution is easily obtainable. With regard to a system with a step discontinuity in the mass per unit length and bending stiffness, [10, 11] demonstrated that an analytical solution can be found for the transverse vibrations by applying continuity conditions at the location of the step. Other times, the solution to a nonuniform, Euler-Bernoulli beam structure can be expressed in terms of special functions such as Bessel functions and Chebyshev polynomials [12–14]. In [12, 13] it is shown that tapered beams admit a solution that can be expressed using Bessel functions. Chebyshev polynomials are used in [14] for the solution to a beam of exponentially varying width clamped at one end with the other end attached to a translational or torsional spring. Unfortunately, it is not always possible to determine an exact analytical solution or express the solution in terms of these special functions with exception of very specific types of nonuniformities some of which are discussed in [7–14]. This is a major difficulty with obtaining a dynamic solution for spatially varying beams.In the absence of an analytical solution, approximation techniques may be applied in the analysis of nonuniform beams. The finite element method [15] and Rayleigh-Ritz method, [16–18] are common approximation techniques for estimating the natural frequencies and mode shapes of a system. Perturbation theory is also a powerful tool for obtaining approximate solutions for systems with complex geometries for which closed form analytical solutions are not readily available. With respect to a nonuniform Euler-Bernoulli beam, the recently developed method of varying amplitudes, a perturbation theory approach, was applied first in [19] to a system with continuous, and periodically varying, mass, and stiffness properties. While quite a powerful method, this approach primarily builds on the assumption of periodic variations in the system parameters over the length of the structure. The transverse vibration of continuous nonuniform Euler-Bernoulli beams with an axial force was studied in [20, 21], and no additional assumptions on the system, such as periodicity, were imposed. The WKB approximation was employed in [20] and the small parameter multiplying the highest derivative was related to the inverse of the natural frequencies of the system. Although this proves to be an advantage for higher frequencies, the method becomes less effective for smaller frequencies. The Lindstedt-Poincaré method was first applied to a nonuniform Euler-Bernoulli beam in [21]. However, the work in [21] was restricted to continuous systems due to lack of continuity conditions at points of discontinuity. In the current work, the Lindstedt-Poincaré method is applied to a nonuniform Euler-Bernoulli beam with piecewise continuous properties.The assumption of the mass per unit length and bending stiffness being discontinuous at a set of points extends the work of [21] and constitutes one of the contributions of the current work. Accounting for a system with piecewise continuous properties will allow for future studies of more complex systems, such as cable-harnessed structures motivated by the space applications [22, 23]. Next, applying the perturbation theory requires a set of reference values for the bending stiffness and mass per unit length from which the spatial functions are perturbed. Additionally, it is shown that the choice of the reference values will impact the accuracy of the perturbation theory results. Therefore, another contribution of the presented work is exploring three methods, one gradient minimization and two Lp-norm minimization methods, by which the reference values in the perturbation theory can be chosen. One method for finding these reference values is considered in [21]. Further, the proposed gradient minimization method in the presented work extends the results in [21] in a more general form to allow for an approach for systems where more than two reference values are required. Providing a general framework for choosing more reference values will prove beneficial for complex systems such as those including axial forces, damping, or an elastic foundation. In these more complex systems, it may not always be intuitive how the multiple reference values should be selected and the presented method in this work will help with this selection. The final main contribution of the current work is the consideration of higher order corrections for the frequencies and mode shapes. With respect to the mode shapes in particular, the higher order correction proves important in the calculation of moment and shear that otherwise will result in significant errors for these calculations.Several beam geometries for a stepped beam and a periodically grooved beam are considered to compare the accuracy of various reference value selection methods when additional corrections in the perturbation theory are considered. Additionally, the moment and shear in the system are investigated and compared for the cases with and without a first-order correction to the mode shapes. For the case of a stepped beam, an analytical solution can be found [10, 11], which serves as an exact solution to which the perturbation theory results can be compared. Determining a solution for a system with piecewise continuous spatial properties, for which in some sections an analytical solution cannot be found, is the goal of the current work. The result is a method that overcomes the need for specific geometries and does not require periodicity or continuity in the system. Future work includes applying the developed approximation technique to more complex geometries, such as cable-harnessed structures discussed in [1–6].


## 2. Problem Formulation

### 2.1. Perturbation Theory

In this work, the free transverse vibration of an Euler-Bernoulli (EB) beam model with spatially dependent bending stiffnessEIx and mass per unit length ρAx is considered. The functions EIx and ρAx are assumed to have derivatives of all orders except on the finite sets of points PEI⊂0,l and PρA⊂0,l, respectively, where l is the total length of the system. That is, EIx and ρAx are C∞-functions on 0,l∖PEI and [0,l]∖PρA, respectively. The partial differential equation (PDE) for the free transverse vibration of the system is (1)∂2∂x2EIx∂2wx,t∂x2+ρAx∂2wx,t∂t2=0for 0<x<l and t>0. The boundary conditions at the endpoints x∗=0 and l are (2a)wx∗,t=0or∂∂xEIx∂2wx,t∂x2x=x∗=0(2b)∂wx∗,t∂x=0orEIx∗∂2wx∗,t∂x2=0The boundary conditions in (2a) can be interpreted physically as the system having zero displacement or zero shear at the endpoints. The boundary conditions in (2b) correspond physically as the system having either zero slope or zero moment at the endpoints.The problem given in (1)-(2a) and (2b) is then nondimensionalized using the reference values EI∗ and ρA∗ for the bending stiffness and mass per unit length of the system, respectively. These values represent a reference point from which the system parameters are considered to be perturbed. The choice of these reference values is studied in the current work and up to this point remains unfixed. Introduce the length scale l and the time scale l2ρA∗/EI∗, and let X, T, and W denote dimensionless quantities. Then, assume a separable solution of the form W(X,T)=u(X)eiωT. Assuming the variations of the spatially dependent functions EI(X) and ρA(X) from the respective reference values are small, the ODE can be written as (3)d2dX21+ϵEI^Xd2uXdX2-ω21+ϵρA^XuX=0,where ϵ is a small parameter and ϵEI^(X)=(EI(X)-EI∗)/EI∗ and ϵρA^(X)=(ρA(X)-ρA∗)/ρA∗. Applying the same steps to the boundary conditions in (2a) and (2b) yields (4a)uX∗=0orddX1+ϵEI^Xd2uXdX2X=X∗=0(4b)duX∗dX=0or1+ϵEI^X∗d2uX∗dX2=0In the current work, the Lindstedt-Poincaré method outlined in [24] is applied to the problem given by (3) and (4a) and (4b). Using this approach, the mode shape and natural frequency are expanded in powers of the small parameter ϵ as u(X)=∑i=0∞‍ϵiui(X) and ω=∑i=0∞‍ϵiωi, respectively. These are substituted into (3) and (4a) and (4b) and terms with similar powers of ϵ are grouped together and equated to 0. This results in a sequence of problems that are solved for the unknown ui(X) and ωi.


### 2.2. Mass Normalization Condition

In a vibrations analysis, it is typical to mass normalize the mode shapes. Denote the mass normalized mode shapes byϕ¯(x) and the dimensionless mass normalized mode shapes by ϕ(X). In the current work, the nondimensionalized mass normalization condition is given by 1=∫01‍1+ϵρA^Xϕ2(X)dX. Using the perturbation theory expansion for the mode shape the zeroth-order, O(1), and first-order, O(ϵ), conditions are(5)O1:1=∫01ϕ02XdX(6)Oϵ:0=∫012ϕ0Xϕ1X+ρA^Xϕ02XdX


## 2.1. Perturbation Theory

In this work, the free transverse vibration of an Euler-Bernoulli (EB) beam model with spatially dependent bending stiffnessEIx and mass per unit length ρAx is considered. The functions EIx and ρAx are assumed to have derivatives of all orders except on the finite sets of points PEI⊂0,l and PρA⊂0,l, respectively, where l is the total length of the system. That is, EIx and ρAx are C∞-functions on 0,l∖PEI and [0,l]∖PρA, respectively. The partial differential equation (PDE) for the free transverse vibration of the system is (1)∂2∂x2EIx∂2wx,t∂x2+ρAx∂2wx,t∂t2=0for 0<x<l and t>0. The boundary conditions at the endpoints x∗=0 and l are (2a)wx∗,t=0or∂∂xEIx∂2wx,t∂x2x=x∗=0(2b)∂wx∗,t∂x=0orEIx∗∂2wx∗,t∂x2=0The boundary conditions in (2a) can be interpreted physically as the system having zero displacement or zero shear at the endpoints. The boundary conditions in (2b) correspond physically as the system having either zero slope or zero moment at the endpoints.The problem given in (1)-(2a) and (2b) is then nondimensionalized using the reference values EI∗ and ρA∗ for the bending stiffness and mass per unit length of the system, respectively. These values represent a reference point from which the system parameters are considered to be perturbed. The choice of these reference values is studied in the current work and up to this point remains unfixed. Introduce the length scale l and the time scale l2ρA∗/EI∗, and let X, T, and W denote dimensionless quantities. Then, assume a separable solution of the form W(X,T)=u(X)eiωT. Assuming the variations of the spatially dependent functions EI(X) and ρA(X) from the respective reference values are small, the ODE can be written as (3)d2dX21+ϵEI^Xd2uXdX2-ω21+ϵρA^XuX=0,where ϵ is a small parameter and ϵEI^(X)=(EI(X)-EI∗)/EI∗ and ϵρA^(X)=(ρA(X)-ρA∗)/ρA∗. Applying the same steps to the boundary conditions in (2a) and (2b) yields (4a)uX∗=0orddX1+ϵEI^Xd2uXdX2X=X∗=0(4b)duX∗dX=0or1+ϵEI^X∗d2uX∗dX2=0In the current work, the Lindstedt-Poincaré method outlined in [24] is applied to the problem given by (3) and (4a) and (4b). Using this approach, the mode shape and natural frequency are expanded in powers of the small parameter ϵ as u(X)=∑i=0∞‍ϵiui(X) and ω=∑i=0∞‍ϵiωi, respectively. These are substituted into (3) and (4a) and (4b) and terms with similar powers of ϵ are grouped together and equated to 0. This results in a sequence of problems that are solved for the unknown ui(X) and ωi.


## 2.2. Mass Normalization Condition

In a vibrations analysis, it is typical to mass normalize the mode shapes. Denote the mass normalized mode shapes byϕ¯(x) and the dimensionless mass normalized mode shapes by ϕ(X). In the current work, the nondimensionalized mass normalization condition is given by 1=∫01‍1+ϵρA^Xϕ2(X)dX. Using the perturbation theory expansion for the mode shape the zeroth-order, O(1), and first-order, O(ϵ), conditions are(5)O1:1=∫01ϕ02XdX(6)Oϵ:0=∫012ϕ0Xϕ1X+ρA^Xϕ02XdX


## 3. Perturbation Theory Solution

### 3.1.O(1) Problem

From the perturbation theory, the ODE for theO(1) problem is d4u0/dX4-ω02u0=0 for 0<X<1. The ODE is identical to that of a homogeneous EB beam model. In the current work, attention is given to clamped-clamped (CC), clamped-free (CF), and free-free (FF) boundary conditions. The boundary conditions for a CC, CF, and FF system are (7)u00=0,du00dX=0,u01=0,du01dX=0(8)u00=0,du00dX=0,d2u01dX2=0,d3u01dX3=0(9)d2u00dX2=0,d3u00dX3=0,d2u01dX2=0,d3u01dX3=0The solution to the ODE given either CC, CF, or FF boundary conditions is readily available and can be found in many references, for example, [15]. To denote the nth mode of the system, a subscript n is added to the terms in the perturbation theory expansion of the mode shape and frequency.


### 3.2.O(ϵ) Problem

The ODE for theO(ϵ) problem is (10)d4un1dX4-ωn02un1=2ωn0ωn1ϕn0+ωn02ρA^ϕn0-d2dX2EI^d2ϕn0dX2.As is typically done in perturbation theory, a solvability condition is applied to ensure that a solution to (10) exists, 0=∫01‍2ωn0ωn1ϕn0+ωn02ρA^ϕn0-d2/dX2EI^(d2ϕn0/dX2)ϕn0dX. Integration by parts and using the fact that CC, CF, and FF boundary conditions are considered in this work, a solution for ωn1 can be found as (11)ωn1=ωn0∫01EI^d2ϕn0/dX22dX2∫01d2ϕn0/dX22dX-ωn0∫01ρA^ϕn02dX2∫01ϕn02dXNext the ODE in (10) is solved for un1(X). Recall that the bending stiffness and mass per unit length are assumed to have derivatives of all orders except on a finite set of points. Denote P=PEI∪PρA⊂[0,l]. Define P^ as the set of all the points in P divided by l; thus, P^⊂[0,1]. Suppose the size of the set P^ is m-1 and order the points in P^ such that 0<XP1<XP2<⋯<XPm-1<1. Then each of the m intervals on which both EI^ and ρA^ are of class C∞ can be denoted using a superscript (i), i.e., using the notation EI^(i) and ρA^(i). The ith interval is defined as [XPi-1,XPi] with XP0=0 and XPm=1. The solution for un1(X) in (10) can be determined over each of the m intervals. The general solution for the ith interval, where un1,p(i)(X) is a particular solution, is (12)un1iX=An1icos⁡ωn0X+Bn1isin⁡ωn0X+Cn1icosh⁡ωn0X+Dn1isinh⁡ωn0X+ωn12ωn0Xdϕn0dX+un1,piXBoundary conditions and continuity conditions are then applied to allow the solution over each of the intervals to be “connected” to each other. Physically it is expected that the solutionun1 will be continuous and have continuous derivatives. Additional continuity conditions can be determined mathematically by integrating the ODE in (10). The continuity conditions at an arbitrary point XPi are listed in (13a), (13b), (13c), and (13d). The continuity conditions of (13c) and (13d) relate physically to the moment and shear in the solution, respectively. (13a)un1iXPi=un1i+1XPi(13b)dun1iXPidX=dun1i+1XPidX(13c)d2un1iXPidX2=d2un1i+1XPidX2+EI^i+1XPi-EI^iXPid2ϕn0XPidX2(13d)d3un1iXPidX3=d3un1i+1XPidX3+ddXEI^i+1X-EI^iXd2ϕn0XdX2X=XPiThe boundary conditions for theO(ϵ) problem can be simplified using the boundary conditions from the O(1) problem. After simplification, the boundary conditions for a CC, CF, and FF system are (14)un110=0,dun110dX=0,un1m1=0,dun1m1dX=0(15)un110=0,dun110dX=0,d2un1m1dX2=0,d3un1m1dX3=0(16)d2un110dX2=0,d3un110dX3=0,d2un1m1dX2=0,d3un1m1dX3=0First, the boundary conditions atX=0 are applied, then the continuity conditions applied sequentially beginning from the left end of the system, XP1, up to the right end of the system, XPm-1. After these have been applied, a single boundary condition at X=1 is applied. In the current work, the third boundary condition listed in (14) to (16) is imposed. Only a single right hand boundary condition needs to be applied since it follows that, for the value of ωn1, the fourth and final boundary condition will automatically be satisfied. The final expression for the ith section of the mode shape, with i≤1≤m, is given by(17)un1iX=An11ϕn0X+ωn12ωn0Xdϕn0dX+u~n1iX+τu~X,where τu~(X) for CC, CF, and FF boundary conditions is (18)τu~X=-u~n1m1sin⁡ωn0X-sinh⁡ωn0Xsin⁡ωn0-sinh⁡ωn0τu~X=d2u~n1m1/dX2sin⁡ωn0X-sinh⁡ωn0Xωn0sin⁡ωn0+sinh⁡ωn0τu~X=d2u~n1m1/dX2sin⁡ωn0X+sinh⁡ωn0Xωn0sin⁡ωn0-sinh⁡ωn0In (17) the function u~n1(i)(X) is calculated in terms of the O(1) solution ϕn0(X) and the functions un1,p(j)(X) with 1≤j≤i. The complete expression for u~n1(i)(X) is given in Appendix A. Applying the mass normalization condition in (6), the final constant An1(1) is determined and the solution on the ith interval is denoted by ϕn1(i)(X). The final coefficient is given by (19)An11=-∫01ωn12ωn0Xϕn0dϕn0dX+ϕn0u~n1+ϕn0τu~X+12ρA^ϕn02dX


### 3.3.O(ϵ2) Problem

The ODE for theO(ϵ2) problem is (20)d4un2dX4-ωn02un2=2ωn0ωn2ϕn0+2ωn0ωn1ϕn1+ωn12ϕn0+2ωn0ωn1ρA^ϕn0+ωn02ρA^ϕn1-d2dX2EI^d2ϕn1dX2A solution to (20) exists when the solvability condition of the integral over the dimensionless domain of the product of ϕn0 and the right hand side of (20) equal to 0 is satisfied. Integration by parts and considering CC, CF, and FF boundary conditions are used in this work; a solution for ωn2 can be found as (21)ωn2=-ωn122ωn0+∫01EI^d2ϕn0/dX2d2ϕn1/dX2-ωn02ρA^ϕn0ϕn1-ωn0ωn1ρA^ϕn02dX2ωn0∫01ϕn02dX


### 3.4. Dimensional Solution and Dependence on Reference Values

In (3), a small parameter was introduced for the perturbation theory via ϵEI^(X)=EI(X)-EI∗/EI∗ and ϵρA^(X)=ρA(X)-ρA∗/ρA∗. Using these, expressions for the frequency and mode shape perturbation solution can be found in terms of the reference values EI∗ and ρA∗. In addition, these expressions will depend on known quantities such as the bending stiffness EI(X) and the mass per unit length ρA(X). (22)ϕnX=32ϕn0X+fn1XρA∗-fn2XEI∗,ωn=ωn0+2αn1EI∗-2αn2ρA∗+βn1EI∗2+βn2EI∗ρA∗+βn3ρA∗2,where (23)αn1=ωn0∫01EIXd2ϕn0/dX22dX2∫01d2ϕn0/dX22dXβn2=ωn0∫01EIXd2fn1/dX2d2ϕn0/dX2dX2∫01d2ϕn0/dX22dX+ωn0∫01ρAXfn2ϕn0dX2∫01ϕn02dXαn2=ωn0∫01ρAXϕn02dX2∫01ϕn02dXβn3=αn222ωn0-ωn0∫01ρAXfn1ϕn0dX2∫01ϕn02dXβn1=-αn122ωn0-ωn0∫01EIXd2fn2/dX2d2ϕn0/dX2dX2∫01d2ϕn0/dX22dXThe complete expressions of fn1(X) and fn2(X) are listed in Appendix B.The final step is to reintroduce dimensionality. For the frequencies of the system, the nondimensionalized frequencies are divided by the time scale. The dimensions of mass normalized mode shapes are1/kg. For the choice of length and time scales used in the current work, the mass scale is taken to be ρA∗l. To reintroduce dimensions, the mode shapes obtained from the perturbation theory must be divided by the square root of the mass scale. (24)ϕ¯nx=1ρA∗l32ϕn0xl+fn1x/lρA∗-fn2x/lEI∗(25)ω¯nEI∗,ρA∗=1l2EI∗ρA∗ωn0+2αn1EI∗-2αn2ρA∗+βn1EI∗2+βn2EI∗ρA∗+βn3ρA∗2


## 3.1.O(1) Problem

From the perturbation theory, the ODE for theO(1) problem is d4u0/dX4-ω02u0=0 for 0<X<1. The ODE is identical to that of a homogeneous EB beam model. In the current work, attention is given to clamped-clamped (CC), clamped-free (CF), and free-free (FF) boundary conditions. The boundary conditions for a CC, CF, and FF system are (7)u00=0,du00dX=0,u01=0,du01dX=0(8)u00=0,du00dX=0,d2u01dX2=0,d3u01dX3=0(9)d2u00dX2=0,d3u00dX3=0,d2u01dX2=0,d3u01dX3=0The solution to the ODE given either CC, CF, or FF boundary conditions is readily available and can be found in many references, for example, [15]. To denote the nth mode of the system, a subscript n is added to the terms in the perturbation theory expansion of the mode shape and frequency.


## 3.2.O(ϵ) Problem

The ODE for theO(ϵ) problem is (10)d4un1dX4-ωn02un1=2ωn0ωn1ϕn0+ωn02ρA^ϕn0-d2dX2EI^d2ϕn0dX2.As is typically done in perturbation theory, a solvability condition is applied to ensure that a solution to (10) exists, 0=∫01‍2ωn0ωn1ϕn0+ωn02ρA^ϕn0-d2/dX2EI^(d2ϕn0/dX2)ϕn0dX. Integration by parts and using the fact that CC, CF, and FF boundary conditions are considered in this work, a solution for ωn1 can be found as (11)ωn1=ωn0∫01EI^d2ϕn0/dX22dX2∫01d2ϕn0/dX22dX-ωn0∫01ρA^ϕn02dX2∫01ϕn02dXNext the ODE in (10) is solved for un1(X). Recall that the bending stiffness and mass per unit length are assumed to have derivatives of all orders except on a finite set of points. Denote P=PEI∪PρA⊂[0,l]. Define P^ as the set of all the points in P divided by l; thus, P^⊂[0,1]. Suppose the size of the set P^ is m-1 and order the points in P^ such that 0<XP1<XP2<⋯<XPm-1<1. Then each of the m intervals on which both EI^ and ρA^ are of class C∞ can be denoted using a superscript (i), i.e., using the notation EI^(i) and ρA^(i). The ith interval is defined as [XPi-1,XPi] with XP0=0 and XPm=1. The solution for un1(X) in (10) can be determined over each of the m intervals. The general solution for the ith interval, where un1,p(i)(X) is a particular solution, is (12)un1iX=An1icos⁡ωn0X+Bn1isin⁡ωn0X+Cn1icosh⁡ωn0X+Dn1isinh⁡ωn0X+ωn12ωn0Xdϕn0dX+un1,piXBoundary conditions and continuity conditions are then applied to allow the solution over each of the intervals to be “connected” to each other. Physically it is expected that the solutionun1 will be continuous and have continuous derivatives. Additional continuity conditions can be determined mathematically by integrating the ODE in (10). The continuity conditions at an arbitrary point XPi are listed in (13a), (13b), (13c), and (13d). The continuity conditions of (13c) and (13d) relate physically to the moment and shear in the solution, respectively. (13a)un1iXPi=un1i+1XPi(13b)dun1iXPidX=dun1i+1XPidX(13c)d2un1iXPidX2=d2un1i+1XPidX2+EI^i+1XPi-EI^iXPid2ϕn0XPidX2(13d)d3un1iXPidX3=d3un1i+1XPidX3+ddXEI^i+1X-EI^iXd2ϕn0XdX2X=XPiThe boundary conditions for theO(ϵ) problem can be simplified using the boundary conditions from the O(1) problem. After simplification, the boundary conditions for a CC, CF, and FF system are (14)un110=0,dun110dX=0,un1m1=0,dun1m1dX=0(15)un110=0,dun110dX=0,d2un1m1dX2=0,d3un1m1dX3=0(16)d2un110dX2=0,d3un110dX3=0,d2un1m1dX2=0,d3un1m1dX3=0First, the boundary conditions atX=0 are applied, then the continuity conditions applied sequentially beginning from the left end of the system, XP1, up to the right end of the system, XPm-1. After these have been applied, a single boundary condition at X=1 is applied. In the current work, the third boundary condition listed in (14) to (16) is imposed. Only a single right hand boundary condition needs to be applied since it follows that, for the value of ωn1, the fourth and final boundary condition will automatically be satisfied. The final expression for the ith section of the mode shape, with i≤1≤m, is given by(17)un1iX=An11ϕn0X+ωn12ωn0Xdϕn0dX+u~n1iX+τu~X,where τu~(X) for CC, CF, and FF boundary conditions is (18)τu~X=-u~n1m1sin⁡ωn0X-sinh⁡ωn0Xsin⁡ωn0-sinh⁡ωn0τu~X=d2u~n1m1/dX2sin⁡ωn0X-sinh⁡ωn0Xωn0sin⁡ωn0+sinh⁡ωn0τu~X=d2u~n1m1/dX2sin⁡ωn0X+sinh⁡ωn0Xωn0sin⁡ωn0-sinh⁡ωn0In (17) the function u~n1(i)(X) is calculated in terms of the O(1) solution ϕn0(X) and the functions un1,p(j)(X) with 1≤j≤i. The complete expression for u~n1(i)(X) is given in Appendix A. Applying the mass normalization condition in (6), the final constant An1(1) is determined and the solution on the ith interval is denoted by ϕn1(i)(X). The final coefficient is given by (19)An11=-∫01ωn12ωn0Xϕn0dϕn0dX+ϕn0u~n1+ϕn0τu~X+12ρA^ϕn02dX


## 3.3.O(ϵ2) Problem

The ODE for theO(ϵ2) problem is (20)d4un2dX4-ωn02un2=2ωn0ωn2ϕn0+2ωn0ωn1ϕn1+ωn12ϕn0+2ωn0ωn1ρA^ϕn0+ωn02ρA^ϕn1-d2dX2EI^d2ϕn1dX2A solution to (20) exists when the solvability condition of the integral over the dimensionless domain of the product of ϕn0 and the right hand side of (20) equal to 0 is satisfied. Integration by parts and considering CC, CF, and FF boundary conditions are used in this work; a solution for ωn2 can be found as (21)ωn2=-ωn122ωn0+∫01EI^d2ϕn0/dX2d2ϕn1/dX2-ωn02ρA^ϕn0ϕn1-ωn0ωn1ρA^ϕn02dX2ωn0∫01ϕn02dX


## 3.4. Dimensional Solution and Dependence on Reference Values

In (3), a small parameter was introduced for the perturbation theory via ϵEI^(X)=EI(X)-EI∗/EI∗ and ϵρA^(X)=ρA(X)-ρA∗/ρA∗. Using these, expressions for the frequency and mode shape perturbation solution can be found in terms of the reference values EI∗ and ρA∗. In addition, these expressions will depend on known quantities such as the bending stiffness EI(X) and the mass per unit length ρA(X). (22)ϕnX=32ϕn0X+fn1XρA∗-fn2XEI∗,ωn=ωn0+2αn1EI∗-2αn2ρA∗+βn1EI∗2+βn2EI∗ρA∗+βn3ρA∗2,where (23)αn1=ωn0∫01EIXd2ϕn0/dX22dX2∫01d2ϕn0/dX22dXβn2=ωn0∫01EIXd2fn1/dX2d2ϕn0/dX2dX2∫01d2ϕn0/dX22dX+ωn0∫01ρAXfn2ϕn0dX2∫01ϕn02dXαn2=ωn0∫01ρAXϕn02dX2∫01ϕn02dXβn3=αn222ωn0-ωn0∫01ρAXfn1ϕn0dX2∫01ϕn02dXβn1=-αn122ωn0-ωn0∫01EIXd2fn2/dX2d2ϕn0/dX2dX2∫01d2ϕn0/dX22dXThe complete expressions of fn1(X) and fn2(X) are listed in Appendix B.The final step is to reintroduce dimensionality. For the frequencies of the system, the nondimensionalized frequencies are divided by the time scale. The dimensions of mass normalized mode shapes are1/kg. For the choice of length and time scales used in the current work, the mass scale is taken to be ρA∗l. To reintroduce dimensions, the mode shapes obtained from the perturbation theory must be divided by the square root of the mass scale. (24)ϕ¯nx=1ρA∗l32ϕn0xl+fn1x/lρA∗-fn2x/lEI∗(25)ω¯nEI∗,ρA∗=1l2EI∗ρA∗ωn0+2αn1EI∗-2αn2ρA∗+βn1EI∗2+βn2EI∗ρA∗+βn3ρA∗2


## 4. Reference Values Selection

### 4.1. Reference Selection Method 1: MinimizingLp Norm of Difference

The first proposed approach for determining the reference values is to minimize theLp norm of the difference between the reference values and its corresponding spatial function. Mathematically, (26)EI∗=minEI∗∈R∫0lEIx-EI∗pdx1/p,ρA∗=minρA∗∈R∫0lρAx-ρA∗pdx1/p,where p is a positive integer. Minimizing the Lp norm of the difference results in a single value for EI∗ and ρA∗ for all the modes. In general the minimization procedure listed in (26) is performed numerically.


### 4.2. Reference Selection Method 2: MinimizingLp Norm of Perturbation

The second proposed approach for determining the reference values is to minimize theLp norm of perturbation terms in (3). Mathematically, (27)EI∗=minEI∗∈R∫01EIX-EI∗EI∗pdX1/p,ρA∗=minρA∗∈R∫01ρAX-ρA∗ρA∗pdX1/p,where p is a positive integer. Minimizing the Lp norm of the perturbation results in a single value for EI∗ and ρA∗ for all the modes. It should be mentioned that the minimization procedure in (27) is performed with respect to the nondimensional domain [0,1], while the minimization procedure in (26) is performed with respect to the original domain [0,l]. In (27) the domain [0,1] is used since this is the domain in consideration when applying the perturbation theory.


### 4.3. Reference Selection Method 3: Gradient Norm Minimization

The final proposed approach to determine the reference values is minimizing theL2 (Euclidean) norm of the gradient of the expression for the natural frequencies given in (25). The gradient is taken with respect to the reference values EI∗ and ρA∗. The gradient minimization procedure relies on the idea that if the reference values are chosen in a manner such that the perturbation theory frequency value is close to the actual value, then any small change in either of the reference values will not greatly affect the predicted frequency. On the other hand, if one of the reference values was chosen to be either very large or very small compared to the optimal reference value, then it is expected that a small change in its value would entail a large change in the predicted frequency.Additionally, constraints are placed on the minimization. The first constraint that is placed on the gradient norm minimization is that the sum of the first- and second-order corrections to the frequencies is 0. This will provide a simple formula for the frequencies in (25). The second constraints EI∗ and ρA∗ are restricted to be within the minimum and maximum values of the function to which they are related. Mathematically the gradient norm minimization method can be described as (28a)minimize:∇ω¯n2subject  to:minx∈0,lEIx≤EI∗≤maxx∈0,lEIx(28b)minx∈0,lρAx≤ρA∗≤maxx∈0,lρAx(28c)2αn1EI∗-2αn2ρA∗+βn1EI∗2+βn2EI∗ρA∗+βn3ρA∗2=0and the L2 norm of the gradient for the perturbation theory frequency of (25) is given by (29)∇ω¯n2=12l21EI∗ρA∗ωn0-2αn1EI∗-2αn2ρA∗-3βn1EI∗2-βn2EI∗ρA∗+βn3ρA∗22+EI∗ρA∗3ωn0+2αn1EI∗-6αn2ρA∗+βn1EI∗2+3βn2EI∗ρA∗+5βn3ρA∗221/2For each mode of the system that is of interest, a solution to the problem outlined in (28a), (28b), and (28c) is determined. This results in a different reference value for each of the modes and these modally dependent reference values are denoted by EIn∗ and ρAn∗. Obtaining references values that depend on the mode number is in contrast to the previously proposed methods in which a single reference value was produced for all the modes.Finally, it should be noted that if the gradient minimization procedure outlined in (28a), (28b), and (28c) is applied to a system with only a single correction to the frequencies in the perturbation theory, then it can be shown that the reference values are EIn∗=2αn1/ωn0 and ρAn∗=2αn2/ωn0. For this choice of reference values, the gradient norm is 0 and they are also exactly the reference values reported in [21]. Therefore, the gradient norm minimization method outlined in (28a), (28b), and (28c) can be seen as a generalization of the reference selection method detailed in [21].


#### 4.3.1. Proof of Solution Using Gradient Minimization Method

It is required to show that there is at least a single point satisfying the various constraints placed on the gradient minimization problem. In the simple case of a constantEIx and ρAx on the domain, it can be shown by direct substitution that the reference values should be taken as these constants. In this case all the constraints placed by the gradient method are satisfied.Without loss of generality, assume now that bothEI(x) and ρA(x) are not constant over the domain. It is useful to rewrite the dimensional frequency in (25) as ω¯nEI∗,ρA∗=1/l2EI∗/ρA∗ωn0+ΩnEI∗,ρA∗, where ΩnEI∗,ρA∗ represents the sum of the corrections to the frequencies.Consider the set of reference valuesEI1∗=maxx∈[0,l]EIx and ρA1∗=minx∈[0,l]ρAx. Choosing these values, ωn0EI1∗/ρA1∗/l2 will overestimate the nth frequency of the system. Since these reference values overestimate the frequency, the total correction from the perturbation theory will be negative, i.e., ΩnEI1∗,ρA1∗<0. Similarly, consider the set of reference values EI2∗=minx∈[0,l]EIx and ρA2∗=maxx∈[0,l]ρAx. Choosing these values, ωn0EI2∗/ρA2∗/l2 will underestimate the nth frequency of the system. Since these reference values underestimate the frequency, the total correction from the perturbation theory will be positive, i.e., ΩnEI2∗,ρA2∗>0.Referring back to (28a) and (28b), these constraints on the possible values for EI∗ and ρA∗ form a convex set. Since the set is convex, there exists a path in the set from EI1∗,ρA1∗ to EI2∗,ρA2∗. Since it was shown that ΩnEI1∗,ρA1∗<0<ΩnEI2∗,ρA2∗, it follows from a corollary to the intermediate value theorem that there must exist a point EI3∗,ρA3∗ in the convex set such that ΩnEI3∗,ρA3∗=0 [25]. In other words, there is at least one point satisfying the constraint that the sum of the corrections to the natural frequency is equal to 0. This shows that there is at least a single point that satisfies the constraints listed in the gradient minimization method and thus a solution for gradient norm minimization can always be found.


## 4.1. Reference Selection Method 1: MinimizingLp Norm of Difference

The first proposed approach for determining the reference values is to minimize theLp norm of the difference between the reference values and its corresponding spatial function. Mathematically, (26)EI∗=minEI∗∈R∫0lEIx-EI∗pdx1/p,ρA∗=minρA∗∈R∫0lρAx-ρA∗pdx1/p,where p is a positive integer. Minimizing the Lp norm of the difference results in a single value for EI∗ and ρA∗ for all the modes. In general the minimization procedure listed in (26) is performed numerically.


## 4.2. Reference Selection Method 2: MinimizingLp Norm of Perturbation

The second proposed approach for determining the reference values is to minimize theLp norm of perturbation terms in (3). Mathematically, (27)EI∗=minEI∗∈R∫01EIX-EI∗EI∗pdX1/p,ρA∗=minρA∗∈R∫01ρAX-ρA∗ρA∗pdX1/p,where p is a positive integer. Minimizing the Lp norm of the perturbation results in a single value for EI∗ and ρA∗ for all the modes. It should be mentioned that the minimization procedure in (27) is performed with respect to the nondimensional domain [0,1], while the minimization procedure in (26) is performed with respect to the original domain [0,l]. In (27) the domain [0,1] is used since this is the domain in consideration when applying the perturbation theory.


## 4.3. Reference Selection Method 3: Gradient Norm Minimization

The final proposed approach to determine the reference values is minimizing theL2 (Euclidean) norm of the gradient of the expression for the natural frequencies given in (25). The gradient is taken with respect to the reference values EI∗ and ρA∗. The gradient minimization procedure relies on the idea that if the reference values are chosen in a manner such that the perturbation theory frequency value is close to the actual value, then any small change in either of the reference values will not greatly affect the predicted frequency. On the other hand, if one of the reference values was chosen to be either very large or very small compared to the optimal reference value, then it is expected that a small change in its value would entail a large change in the predicted frequency.Additionally, constraints are placed on the minimization. The first constraint that is placed on the gradient norm minimization is that the sum of the first- and second-order corrections to the frequencies is 0. This will provide a simple formula for the frequencies in (25). The second constraints EI∗ and ρA∗ are restricted to be within the minimum and maximum values of the function to which they are related. Mathematically the gradient norm minimization method can be described as (28a)minimize:∇ω¯n2subject  to:minx∈0,lEIx≤EI∗≤maxx∈0,lEIx(28b)minx∈0,lρAx≤ρA∗≤maxx∈0,lρAx(28c)2αn1EI∗-2αn2ρA∗+βn1EI∗2+βn2EI∗ρA∗+βn3ρA∗2=0and the L2 norm of the gradient for the perturbation theory frequency of (25) is given by (29)∇ω¯n2=12l21EI∗ρA∗ωn0-2αn1EI∗-2αn2ρA∗-3βn1EI∗2-βn2EI∗ρA∗+βn3ρA∗22+EI∗ρA∗3ωn0+2αn1EI∗-6αn2ρA∗+βn1EI∗2+3βn2EI∗ρA∗+5βn3ρA∗221/2For each mode of the system that is of interest, a solution to the problem outlined in (28a), (28b), and (28c) is determined. This results in a different reference value for each of the modes and these modally dependent reference values are denoted by EIn∗ and ρAn∗. Obtaining references values that depend on the mode number is in contrast to the previously proposed methods in which a single reference value was produced for all the modes.Finally, it should be noted that if the gradient minimization procedure outlined in (28a), (28b), and (28c) is applied to a system with only a single correction to the frequencies in the perturbation theory, then it can be shown that the reference values are EIn∗=2αn1/ωn0 and ρAn∗=2αn2/ωn0. For this choice of reference values, the gradient norm is 0 and they are also exactly the reference values reported in [21]. Therefore, the gradient norm minimization method outlined in (28a), (28b), and (28c) can be seen as a generalization of the reference selection method detailed in [21].


### 4.3.1. Proof of Solution Using Gradient Minimization Method

It is required to show that there is at least a single point satisfying the various constraints placed on the gradient minimization problem. In the simple case of a constantEIx and ρAx on the domain, it can be shown by direct substitution that the reference values should be taken as these constants. In this case all the constraints placed by the gradient method are satisfied.Without loss of generality, assume now that bothEI(x) and ρA(x) are not constant over the domain. It is useful to rewrite the dimensional frequency in (25) as ω¯nEI∗,ρA∗=1/l2EI∗/ρA∗ωn0+ΩnEI∗,ρA∗, where ΩnEI∗,ρA∗ represents the sum of the corrections to the frequencies.Consider the set of reference valuesEI1∗=maxx∈[0,l]EIx and ρA1∗=minx∈[0,l]ρAx. Choosing these values, ωn0EI1∗/ρA1∗/l2 will overestimate the nth frequency of the system. Since these reference values overestimate the frequency, the total correction from the perturbation theory will be negative, i.e., ΩnEI1∗,ρA1∗<0. Similarly, consider the set of reference values EI2∗=minx∈[0,l]EIx and ρA2∗=maxx∈[0,l]ρAx. Choosing these values, ωn0EI2∗/ρA2∗/l2 will underestimate the nth frequency of the system. Since these reference values underestimate the frequency, the total correction from the perturbation theory will be positive, i.e., ΩnEI2∗,ρA2∗>0.Referring back to (28a) and (28b), these constraints on the possible values for EI∗ and ρA∗ form a convex set. Since the set is convex, there exists a path in the set from EI1∗,ρA1∗ to EI2∗,ρA2∗. Since it was shown that ΩnEI1∗,ρA1∗<0<ΩnEI2∗,ρA2∗, it follows from a corollary to the intermediate value theorem that there must exist a point EI3∗,ρA3∗ in the convex set such that ΩnEI3∗,ρA3∗=0 [25]. In other words, there is at least one point satisfying the constraint that the sum of the corrections to the natural frequency is equal to 0. This shows that there is at least a single point that satisfies the constraints listed in the gradient minimization method and thus a solution for gradient norm minimization can always be found.


## 4.3.1. Proof of Solution Using Gradient Minimization Method

It is required to show that there is at least a single point satisfying the various constraints placed on the gradient minimization problem. In the simple case of a constantEIx and ρAx on the domain, it can be shown by direct substitution that the reference values should be taken as these constants. In this case all the constraints placed by the gradient method are satisfied.Without loss of generality, assume now that bothEI(x) and ρA(x) are not constant over the domain. It is useful to rewrite the dimensional frequency in (25) as ω¯nEI∗,ρA∗=1/l2EI∗/ρA∗ωn0+ΩnEI∗,ρA∗, where ΩnEI∗,ρA∗ represents the sum of the corrections to the frequencies.Consider the set of reference valuesEI1∗=maxx∈[0,l]EIx and ρA1∗=minx∈[0,l]ρAx. Choosing these values, ωn0EI1∗/ρA1∗/l2 will overestimate the nth frequency of the system. Since these reference values overestimate the frequency, the total correction from the perturbation theory will be negative, i.e., ΩnEI1∗,ρA1∗<0. Similarly, consider the set of reference values EI2∗=minx∈[0,l]EIx and ρA2∗=maxx∈[0,l]ρAx. Choosing these values, ωn0EI2∗/ρA2∗/l2 will underestimate the nth frequency of the system. Since these reference values underestimate the frequency, the total correction from the perturbation theory will be positive, i.e., ΩnEI2∗,ρA2∗>0.Referring back to (28a) and (28b), these constraints on the possible values for EI∗ and ρA∗ form a convex set. Since the set is convex, there exists a path in the set from EI1∗,ρA1∗ to EI2∗,ρA2∗. Since it was shown that ΩnEI1∗,ρA1∗<0<ΩnEI2∗,ρA2∗, it follows from a corollary to the intermediate value theorem that there must exist a point EI3∗,ρA3∗ in the convex set such that ΩnEI3∗,ρA3∗=0 [25]. In other words, there is at least one point satisfying the constraint that the sum of the corrections to the natural frequency is equal to 0. This shows that there is at least a single point that satisfies the constraints listed in the gradient minimization method and thus a solution for gradient norm minimization can always be found.


## 5. Numerical Examples

### 5.1. Stepped Beam with Homogeneous Sections

Consider the simple numerical example of a stepped beam with homogeneous sections. Assume that the bending stiffness and mass per unit length of the system are given by(30)EIx=1,0≤x≤l0τ,l0≤x≤l,ρAx=1,0≤x≤l0γ,l0≤x≤l,where τ,γ>0. For this simple test case, there is a single point where both the bending stiffness and mass per unit length are discontinuous. A schematic of the system is shown in Figure 1.Figure 1
Stepped beam system schematic.The error between an analytical and the perturbation theory solutions is used to measure how well the perturbation theory estimates the natural frequencies of the system. Consider the first 10 nonrigid body modes of the system. The error results with variablesτ and γ for multiple boundary conditions are presented in Figures 2 to 3. In each figure, four results are presented: the gradient minimization method with 1 and 2 corrections to the natural frequencies, the norm minimization of the difference, and the norm minimization of the perturbation. As previously noted, the results for the gradient minimization with 1 correction to the frequencies correspond to the work outlined in [21]. For the norm minimization methods an optimal p-value is determined and presented in the figures. This is achieved for each system setup by considering multiple p-values and calculating the errors for the variable parameter being analyzed. Then, the p-value corresponding to the smallest area under the sum of the frequency error curves was determined. This process is also performed for each of the boundary conditions individually.Figure 2
Frequency error results for variableτ.
(a)
Step atl0/l=1/4, CC (b)
Step atl0/l=1/4, CF (c)
Step atl0/l=1/4, FF (d)
Step atl0/l=3/4, CC (e)
Step atl0/l=3/4, CF (f)
Step atl0/l=3/4, FFFigure 3
Frequency error results for variableγ.
(a)
Step atl0/l=1/4, CC (b)
Step atl0/l=1/4, CF (c)
Step atl0/l=1/4, FF (d)
Step atl0/l=3/4, CC (e)
Step atl0/l=3/4, CF (f)
Step atl0/l=3/4, FFWhen bothτ and γ are unity, the system is homogeneous across the entire length. As such, it is expected that the perturbation theory results would exactly match those of the analytical solution. This is observed in Figures 2 to 3 since the error is 0 in this case. Additionally, in each of the figures, the first-order gradient minimization is seen to produce the largest sum of error for a set of system parameters and the norm minimization methods always produce the smallest sum of error. The second-order gradient minimization method produces a sum of error that is between the first-order gradient minimization and norm minimization approaches. The results for the norm minimization methods do not have a consistent pattern as to which of the two approaches will produce a larger or smaller sum of error. It can also be seen that the results for CF boundary conditions are typically smaller than for CC and FF boundary conditions given the same system parameters. In addition, the errors in the frequencies are larger for variable τ compared to variable γ.Consider the results in Figures2 to 3, as the value of τ or γ becomes further from unity the sum of errors increases. This is expected since as these values become increasingly further from 1, the step change in one of the system properties becomes more significant. This causes larger perturbations from the reference value and thus a larger error in the perturbation theory results is anticipated. The largest sum of error occurred when τ or γ was 0.5. Furthermore, it is seen in each of the figures that the results for the error are not symmetric about τ=1 and γ=1. There is also no reason that this type of symmetry would exist as the effect of increasing or decreasing either bending stiffness or mass per unit length will affect the system, and hence perturbation theory results, differently. However, a form of symmetry does exist when comparing the cases of l0/l=1/4 and l0/l=3/4. The results for l0/l=1/4 given variable τ are identical to the results for l0/l=3/4 given 1/τ, and vice versa. For example, the errors reported for τ=3/2 given l0/l=1/4 are identical to the errors reported for τ=2/3 given l0/l=3/4. The same results for symmetry hold true for variable γ.Finally, consider thep-values that produce the smallest area under the error curves for the norm minimization techniques in Figures 2 to 3. For variable τ the p-value used always approaches infinity with the exception of the CF results. For variable γ the p-values are between 3 and 9, inclusively. These results indicate that when the perturbation in the system is due to a change in bending stiffness, a large p-value should be considered, whereas when the perturbation in the system is due to a change in the mass per unit length, a p-value of less than 10 should be employed.Next, dimensional mode shapes associated with the fundamental frequency of the system are presented in Figure4. The results for CC boundary conditions only are presented as they demonstrate the largest errors between perturbation theory and analytical results. In this figure, the moment and shear of the system due to the given mode shape are also presented. A value of l0/l=1/4 is used and the values of γ and τ are 1 unless otherwise specified. The p-values used are taken from the corresponding previous results for variables τ and γ.Figure 4
Mode shape displacement, moment, and shear for CC boundary conditions.
(a)
Displacement,τ=1.5 (b)
Displacement,γ=1.5 (c)
Moment,τ=1.5 (d)
Moment,γ=1.5 (e)
Shear,τ=1.5 (f)
Shear,γ=1.5For the results presented in Figure4, the case of γ=1.5 and τ=1 matches better overall for the displacement, moment, and shear when compared to the case of γ=1 and τ=1.5. A point of interest is that the moment and shear when γ=1 and τ=1.5 are discontinuous for the perturbation theory results, with the point of discontinuity at the location of the step l0/l=1/4, whereas the moment and shear when γ=1.5 and τ=1 are both continuous. This difference is investigated by considering an expression for the continuity of moment and shear for the system.From the right hand side of (4b), the moment at a point x in the system is given by M(X)=1+ϵEI^Xd2ϕ(X)/dX2. In the current work the first-order correction to the mode shapes is found. Therefore, it is assumed that ϕ(X)=ϕ0(X)+ϵϕ1(X). Choose an arbitrary point XPi∈P^, consider the difference in the value of the moment to the left and right of this point, and use the small parameter ϵ definition. (31)MXPi--MXPi+=EIXPi+-EIXPi-EIXPi--EI∗EI∗2d2ϕ0XPidX2+1EI∗d2ϵϕ1XPi+dX2It is seen in (31) that if the bending stiffness of the system is continuous, then the moment will also be continuous. However, if the bending stiffness of the system is not continuous, then the moment of the system may or may not be continuous. In (31) it is also seen that the mass per unit length does not play a role in the continuity of the moment. In the stepped beam case when τ is not equal to one, the bending stiffness is not continuous and this causes the discontinuous moment in Figure 4(c). In the simulations where τ=1 and γ=1.5, the bending stiffness is continuous and therefore the moment in Figure 4(d) is continuous.Using the same approach, the jump in shear at a pointXPi∈P^ can be determined. The shear in the system at a given point X is S(X)=d/dXM(X) and the result of the difference is given by (32)SXPi--SXPi+=dEIXPi+dX-dEIXPi-dXEIXPi--EI∗EI∗2d2ϕ0XPidX2-1EI∗d2ϵϕ1XPi+dX2+EIXPi--EIXPi+1EI∗d3ϵϕ1XPi+dX3-EIXPi--EI∗EI∗2d3ϕ0XPidX3-dEIXPi-/dXEI∗2d2ϕ0XPidX2It can be seen from (32) that if the bending stiffness of the system and its first derivative are continuous, then the shear will also be continuous. However, if the bending stiffness of the system is not continuous or the first derivative is not continuous, then the shear of the system may or may not be continuous. In (32) it is seen that the mass per unit length does not play a role in the continuity of the shear. In the stepped beam case when τ is not equal to one, the bending stiffness and the first derivative are both not continuous and this is the reason for the discontinuous shear in Figure 4(e). In the simulations where τ=1 and γ=1.5 the bending stiffness and its first derivative are continuous and therefore the shear in Figure 4(f) is continuous.Table1 presents the absolute errors in predicting the moment and shear for the various boundary conditions. It is clearly demonstrated in Table 1 that the norm minimization methods provide significantly better estimations for the moment and shear when compared to the first-order gradient method. This is expected for the first-order gradient method since there are fewer corrections to the mode shapes. The error in using the first-order gradient method is typically about 5 times larger when compared to the norm minimization methods, with the errors being up to 30 times higher in the case of the shear for the CF boundary conditions. This is a major contribution of the current work and clearly shows the importance of including the first-order correction to the mode shapes.Table 1
Absolute percentage of error in moment and shear prediction of perturbation methods compared to analytical results.
Method Moment atl0/l Shear atl0/l Moment atl0/l Shear atl0/l Moment atl0/l Shear atl0/l 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + Clamped-Clamped Clamped-Free Free-Free 1st Grad. 62.94% 44.40% 21.44% 17.84% 9.90% 35.02% 11.05% 33.43% 32.74% 0.89% 31.64% 2.54% 2nd Grad. 15.99% 5.32% 4.55% 0.67% 2.33% 7.73% 2.42% 6.87% 10.24% 0.32% 9.53% 0.82% L 6-Diff. - - - - 0.64% 1.89% 0.58% 1.36% - - - - L 6-Pert. - - - - 0.56% 1.60% 0.49% 1.09% - - - - L ∞-Diff. 12.48% 9.04% 3.28% 2.06% - - - - 6.28% 0.22% 5.64% 0.52% L ∞-Pert. 12.48% 9.04% 3.28% 2.06% - - - - 6.28% 0.22% 5.64% 0.52%The final results presented for the stepped beam system are the frequency error results for a variable step location. Figures5 and 6 consider the cases when either γ or τ is equal to 0.9 or 1.1, respectively. CC, CF, and FF boundary conditions are all considered. As was the case for variable τ and γ, an optimal p-value is determined for the norm minimization methods for each system setup. Note that since now a different parameter is being varied, l0/l, the optimal p-value is not expected to be the same as the case for variable τ and γ.Figure 5
Frequency error results with a decrease inγ or τ for variable l0/l.
(a)
τ = 0.9, γ=1, CC (b)
τ = 0.9, γ=1, CF (c)
τ = 0.9, γ=1, FF (d)
τ = 1, γ=0.9, CC (e)
τ = 1, γ=0.9, CF (f)
τ = 1, γ=0.9, FFFigure 6
Frequency error results with an increase inγ or τ for variable l0/l.
(a)
τ = 1.1, γ=1, CC (b)
τ = 1.1, γ=1, CF (c)
τ = 1.1, γ=1, FF (d)
τ = 1, γ=1.1, CC (e)
τ = 1, γ=1.1, CF (f)
τ = 1, γ=1.1, FFIn Figures5 and 6, it can be seen that for a 10% increase or decrease in the bending stiffness or mass per unit length the proposed methods for choosing the reference values in the perturbation theory lead to very small error values. The first-order gradient method produced significantly larger error values when compared to the other proposed methods. For all the presented simulations, the largest sum of absolute error of the frequencies for the first-order gradient method is 0.921%, for the second-order gradient method is 0.050%, for the difference norm minimization is 0.046%, and for the perturbation norm minimization is 0.046%. This indicates that the perturbation theory detailed in the current work provides very accurate estimations for the natural frequencies of the system, even in the case of fairly large changes to the bending stiffness or the mass per unit length.More specifically from Figures5 and 6, the cases when the bending stiffness is varied produce typically larger errors in the frequencies when compared to the cases when the mass per unit length is varied. This is similar to the observations in the variables τ and γ results of Figures 2 and 3, as well as the moment and shear results. Additionally, the error results for CF boundary conditions are typically smaller than for CC and FF boundary conditions given the same system parameters. A similar observation occurred in the previous variables τ and γ results.Finally, in Figures5 and 6, it is seen that the largest errors occur near l0/l=1/2 and the metric value is 0 when l0/l=0 and l0/l=1. When the step location in the beam is near the middle of the system, it becomes more difficult to determine a single reference value that accurately represents the bending stiffness or mass per unit length of the system. Therefore there will be a larger error associated with the perturbation theory methods. When l0/l=0 or l0/l=1, the entire system has a single value for the bending stiffness and mass per unit length. In this case, the system is a homogeneous beam and the perturbation theory solution is exactly the analytical solution. Since the perturbation theory solution is exactly the analytical solution, the error in the frequencies is 0.


### 5.2. Periodically Grooved Beam

The second system to which the perturbation theory developed in the current work is applied is a grooved beam, as shown in Figure7. In this example, at the middle point of the groove the system reaches a thickness value of ch, where h is the thickness of the beam where there is no groove. Thus, when c>1, the beam becomes thicker in the grooved section as in Figure 7(a); when 0<c<1, the beam becomes thinner in the grooved section as in Figure 7(b). If c=1 then there are no grooves in the system and the beam is homogeneous. A system with multiple grooves consists of multiple single groove elements from either Figure 7(a) or Figure 7(b) assembled side by side.Figure 7
Grooved beam system schematics.
(a)
Grooved element withc>1 (b)
Grooved element with0<c<1Consider a single grooved section shown in Figures7(a) and 7(b). In this case, the bending stiffness and mass per unit length of the system are continuous and the derivatives are discontinuous. Let L denote the total length of a single grooved element. The bending stiffness and mass per unit length of a single grooved element can be expressed as (33)EIx=1,0≤x≤L41+c-1x-L/4L/43,L4≤x≤L2c+1-cx-L/2L/43,L2≤x≤3L41,3L4≤x≤L,ρAx=1,0≤x≤L41+c-1x-L/4L/4,L4≤x≤L2c+1-cx-L/2L/4,L2≤x≤3L41,3L4≤x≤LAs in the previous numerical example, the sum of the absolute error in the frequencies is used to provide a measure of how well the perturbation theory compares to an exact solution. The first 10 nonrigid body modes are considered and compared to a finite elements analysis (FEA). The finite element analysis employed in the current work follows the approach outlined in [15]. The error results for a system with variable c3, multiple boundary conditions applied, and varying number of total grooves in a fixed length system are presented in Figures 8–11. As in the stepped beam example, for the norm minimization methods multiple p-values were considered and the errors calculated for the variable parameter are being analyzed. Then, the p-value corresponding to the smallest area under the error curves was determined and plotted in the figures.Figure 8
Frequency error results for variablec using the first-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 9
Frequency error results for variablec using the second-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 10
Frequency error results for variablec using the norm minimization of the difference.
(a)
CC (b)
CF (c)
FFFigure 11
Frequency error results for variablec using the norm minimization of the perturbation.
(a)
CC (b)
CF (c)
FFThe smallestc3 value considered produces a 50% and 20.6% decrease in the bending stiffness and mass per unit length of the system, respectively. The largest c3 value considered produces a 50% and 14.5% increase in the bending stiffness and mass per unit length of the system, respectively. The ratio of the change in the bending stiffness from the centre of the groove to a nongrooved section, which is exactly c3, is larger than the ratio of the change in the mass per unit length from the centre of the groove to a nongrooved section, which is exactly c. For this reason the c3 value is used as the independent variable in the figures.In Figures8–11 it is seen that when c3=1 the errors are 0. This is expected since in this case the beam is homogeneous and the perturbation theory and FEA produce identical results. Additionally in Figures 8–11, it is observed that the first-order gradient minimization method always produces the largest error of the four methods. This is expected since in this case there is only a single correction to the natural frequencies, whereas in the other methods there are two corrections to the natural frequencies and a correction to the mode shapes. These additional corrections provide an increased accuracy. The next largest error is the second-order gradient results in Figure 9, then the norm minimization results in Figures 10 and 11. The largest sum of absolute error in the frequencies across all boundary conditions for the first-order and second-order gradient methods and the difference and perturbation norm minimization methods is 19.835%, 4.501%, 0.921%, and 0.9353%, respectively.As well in Figures8–11, the results for the errors of a given system are not symmetric about c3=1 and there are larger errors when c3=0.5 compared to when c3=1.5. To analyze this, consider a c3 value a distance Δ from c3=1. Then the points to the left and right of c3=1 that are equidistant have a value of c=1±Δ3. For a given Δ, the percentage change in the c value from c=1 is larger when Δ is subtracted in the radicand than when it is added. This indicates that in the case that we are to the left of c3=1 the effect on the mass per unit length is greater than for an equidistant point to the right of c3=1. Due to this, there are higher errors for the cases when c3<1 than when c3>1 and the results for the error are not symmetric about c3=1.Most importantly, the results of Figures8–11 demonstrate the accuracy for a large range of c3 values when using two corrections to the frequencies and a single correction to the mode shapes. In particular, using the norm minimization methods to select the reference values in the perturbation theory, a very small error is produced even when there is up to a 50% increase or decrease in the bending stiffness of the system. Furthermore, the error curve for these methods is quite flat for c3 values near 1 when compared to the gradient minimization methods. These properties highlight the robustness of the proposed methods, in particular the norm minimization approaches.Additionally, it is seen in Figures8 and 9 that the errors for the frequencies are smallest for a system with 1 groove and increase as the number of grooves increases for all the boundary conditions. The one exception to this was for the second-order gradient minimization method with FF boundary conditions for c3>1 where a system with 3 grooves produced slightly larger errors than a system with 5 grooves. This behaviour of increasing frequency error for increasing number of grooves indicates that a system with more rapidly varying bending stiffness and mass per unit length will present high errors when using a perturbation theory employing the gradient minimization method.Next, it is observed in Figures10 and 11 for the difference and perturbation norm minimization methods that there is no consistent pattern with respect to errors in the frequencies and number of grooves. For the norm minimization methods, a single value for EI∗ and ρA∗ is found, whereas for the gradient minimization method a reference value is found for each mode of the system. This is likely the main contributing factor for the lack of pattern in the errors for the norm minimization methods compared to the gradient minimization method. It can be seen, however, that there is a consistent behaviour amongst the results for the difference and perturbation norm minimization.Consider the expressions for the moment and shear using the perturbation theory presented in (31) and (32). For the grooved system in the current example, the bending stiffness of the system is continuous over the domain and discontinuous at the start, centre, and end of each groove. Due to this, it is expected that the perturbation theory results will produce a moment that is continuous and a shear that is discontinuous. As an example, Figure 12 presents the results for CF boundary conditions with c3=1.1 and 1 groove in the system. The expected behaviours for the moment and shear are clearly seen in these results and the advantage of including additional corrections to the mode shape is particularly highlighted in the shear results of Figure 12(c).Figure 12
Mode shape displacement, moment, and shear for clamped-free boundary conditions.
(a)
Displacement (b)
Moment (c)
ShearFinally, the errors for a system with a fixedc3 value and variable number of grooves in a system with a fixed length are presented in Figures 13–16. Multiple c3 values are considered. As before, for the norm minimization methods multiple p-values were considered and the errors calculated for the variable parameter are being studied. Then, the p-value corresponding to the smallest area under the error curves was determined and plotted in the figures. It is not expected that the optimal p-values for systems considering variable number of grooves are the same as the case for variable c3 values.Figure 13
Frequency errors for variable groove number and first-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 14
Frequency errors for variable groove number and second-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 15
Frequency errors for variable groove number and norm minimization of difference.
(a)
CC (b)
CF (c)
FFFigure 16
Frequency errors for variable groove number and norm minimization of perturbation.
(a)
CC (b)
CF (c)
FFIn Figures13–16, it is observed that the first-order gradient minimization method always produces the largest error of the four presented methods. This is expected since in this case there is only a single correction to the natural frequencies, whereas in the other methods there are two corrections to the natural frequencies and a correction to the mode shapes. These additional corrections provide an increased accuracy. Also the figures show that the norm minimization methods produce smaller errors when compared to the second-order gradient method. In addition, the norm minimization methods are seen to provide very accurate results for all the number of grooves considered when c3 is near one. The largest sum of absolute error in the frequencies across all boundary conditions for the first-order and second-order gradient methods and the difference and perturbation norm minimization methods is 23.436%, 5.709%, 0.884%, and 0.888%, respectively.Overall in Figures13–16, there is a clear pattern in the levels of error reported for all the boundary conditions and all the proposed methods for determining the reference value. For a given system, the errors reported are largest for c3=0.5, the next largest are for c3=1.5, the second smallest are for c3=0.9, and the smallest errors are for c3=1.1. It was previously shown that for equidistant points from c3=1 that the c3<1 value will report higher errors than the c3>1 value. Indeed this is observed throughout the presented results for variable number of grooves. In Figures 13–16, it is not expected that the error is 0 for any of the simulations since the case of c3=1 is not considered.It is seen in Figures13 and 14 for the gradient minimization methods that the errors for the frequencies increase as the number of grooves in the system increases and then reaches an upper bound. This upper bound is achieved for a smaller number of grooves by the second-order gradient minimization compared to the first-order method due to the additional corrections considered in the perturbation theory. This is an advantage of the second-order method. This overall increase in the errors up to a given number of grooves in the system is consistent with the behaviour that was observed in the results considering variable c3. For both the first- and second-order gradient minimization methods, the behaviour observed was consistent across all the boundary conditions.Next, it is witnessed for the difference and perturbation norm results in Figures15 and 16 that the behaviour is quite different than for the gradient minimization approach. This is due to the manner in which the reference values are determined in each of these cases. For the gradient minimization a reference value EI∗ and ρA∗ is determined for each mode of the system, whereas the norm minimization methods predict only a single value for all the modes of the system. It can be seen, however, that similar behaviours are observed between the results of the difference and perturbation norm minimization methods. For CC boundary conditions the frequency errors are fairly constant, with a slight overall increase for the case of c3=0.5. For CF and FF boundary conditions, the error in the frequencies noticeably decreases as the number of grooves increases for the cases of c3=0.5 and c3=1.5. It should be noted that the pattern observed is due to the manner in which the p-value was determined. If, for example, the optimal p-value was determined for each data point in the figures by minimizing the errors for the given data point, then the pattern of the errors for variable number of grooves would be different. A similar approach for minimizing the errors for each given data point can be applied for the figures where the effects of a variable c3 are discussed, which will result in a different pattern for Figures 10 and 11. Regardless of which p-value is chosen for the system optimization, the norm minimization method consistently gives better error values compared to the gradient minimization method.A final note is made with regard to the results observed in Figures13 to 16 for variable number of grooves considering multiple c3 values in comparison to the results observed in Figures 8 to 11 for variable c3 value considering multiple number of grooves in the system. For the gradient minimization methods, the method for determining the reference values is prescribed and the results are seen to be consistent across all the simulations. That is, increasing the number of grooves in the system from 1 to 7 produces larger errors, as well as increasing the c3 value such that it is further from c3=1. However, it is seen that slightly different behaviours occur for the results of the difference and perturbation norm minimization. This is due to different optimal p-values being determined for the cases of variable c3 and variable number of grooves. As such, it is not expected that the pattern of the errors is consistent across all the numerical simulations. As previously mentioned, the choice of p-value will have an effect on the behaviour of the errors. Choosing two sets of optimal p-values that are based on which type of variation is allowed in the system demonstrates how the results can be made more accurate than if a single value for p was used throughout all the simulations.


## 5.1. Stepped Beam with Homogeneous Sections

Consider the simple numerical example of a stepped beam with homogeneous sections. Assume that the bending stiffness and mass per unit length of the system are given by(30)EIx=1,0≤x≤l0τ,l0≤x≤l,ρAx=1,0≤x≤l0γ,l0≤x≤l,where τ,γ>0. For this simple test case, there is a single point where both the bending stiffness and mass per unit length are discontinuous. A schematic of the system is shown in Figure 1.Figure 1
Stepped beam system schematic.The error between an analytical and the perturbation theory solutions is used to measure how well the perturbation theory estimates the natural frequencies of the system. Consider the first 10 nonrigid body modes of the system. The error results with variablesτ and γ for multiple boundary conditions are presented in Figures 2 to 3. In each figure, four results are presented: the gradient minimization method with 1 and 2 corrections to the natural frequencies, the norm minimization of the difference, and the norm minimization of the perturbation. As previously noted, the results for the gradient minimization with 1 correction to the frequencies correspond to the work outlined in [21]. For the norm minimization methods an optimal p-value is determined and presented in the figures. This is achieved for each system setup by considering multiple p-values and calculating the errors for the variable parameter being analyzed. Then, the p-value corresponding to the smallest area under the sum of the frequency error curves was determined. This process is also performed for each of the boundary conditions individually.Figure 2
Frequency error results for variableτ.
(a)
Step atl0/l=1/4, CC (b)
Step atl0/l=1/4, CF (c)
Step atl0/l=1/4, FF (d)
Step atl0/l=3/4, CC (e)
Step atl0/l=3/4, CF (f)
Step atl0/l=3/4, FFFigure 3
Frequency error results for variableγ.
(a)
Step atl0/l=1/4, CC (b)
Step atl0/l=1/4, CF (c)
Step atl0/l=1/4, FF (d)
Step atl0/l=3/4, CC (e)
Step atl0/l=3/4, CF (f)
Step atl0/l=3/4, FFWhen bothτ and γ are unity, the system is homogeneous across the entire length. As such, it is expected that the perturbation theory results would exactly match those of the analytical solution. This is observed in Figures 2 to 3 since the error is 0 in this case. Additionally, in each of the figures, the first-order gradient minimization is seen to produce the largest sum of error for a set of system parameters and the norm minimization methods always produce the smallest sum of error. The second-order gradient minimization method produces a sum of error that is between the first-order gradient minimization and norm minimization approaches. The results for the norm minimization methods do not have a consistent pattern as to which of the two approaches will produce a larger or smaller sum of error. It can also be seen that the results for CF boundary conditions are typically smaller than for CC and FF boundary conditions given the same system parameters. In addition, the errors in the frequencies are larger for variable τ compared to variable γ.Consider the results in Figures2 to 3, as the value of τ or γ becomes further from unity the sum of errors increases. This is expected since as these values become increasingly further from 1, the step change in one of the system properties becomes more significant. This causes larger perturbations from the reference value and thus a larger error in the perturbation theory results is anticipated. The largest sum of error occurred when τ or γ was 0.5. Furthermore, it is seen in each of the figures that the results for the error are not symmetric about τ=1 and γ=1. There is also no reason that this type of symmetry would exist as the effect of increasing or decreasing either bending stiffness or mass per unit length will affect the system, and hence perturbation theory results, differently. However, a form of symmetry does exist when comparing the cases of l0/l=1/4 and l0/l=3/4. The results for l0/l=1/4 given variable τ are identical to the results for l0/l=3/4 given 1/τ, and vice versa. For example, the errors reported for τ=3/2 given l0/l=1/4 are identical to the errors reported for τ=2/3 given l0/l=3/4. The same results for symmetry hold true for variable γ.Finally, consider thep-values that produce the smallest area under the error curves for the norm minimization techniques in Figures 2 to 3. For variable τ the p-value used always approaches infinity with the exception of the CF results. For variable γ the p-values are between 3 and 9, inclusively. These results indicate that when the perturbation in the system is due to a change in bending stiffness, a large p-value should be considered, whereas when the perturbation in the system is due to a change in the mass per unit length, a p-value of less than 10 should be employed.Next, dimensional mode shapes associated with the fundamental frequency of the system are presented in Figure4. The results for CC boundary conditions only are presented as they demonstrate the largest errors between perturbation theory and analytical results. In this figure, the moment and shear of the system due to the given mode shape are also presented. A value of l0/l=1/4 is used and the values of γ and τ are 1 unless otherwise specified. The p-values used are taken from the corresponding previous results for variables τ and γ.Figure 4
Mode shape displacement, moment, and shear for CC boundary conditions.
(a)
Displacement,τ=1.5 (b)
Displacement,γ=1.5 (c)
Moment,τ=1.5 (d)
Moment,γ=1.5 (e)
Shear,τ=1.5 (f)
Shear,γ=1.5For the results presented in Figure4, the case of γ=1.5 and τ=1 matches better overall for the displacement, moment, and shear when compared to the case of γ=1 and τ=1.5. A point of interest is that the moment and shear when γ=1 and τ=1.5 are discontinuous for the perturbation theory results, with the point of discontinuity at the location of the step l0/l=1/4, whereas the moment and shear when γ=1.5 and τ=1 are both continuous. This difference is investigated by considering an expression for the continuity of moment and shear for the system.From the right hand side of (4b), the moment at a point x in the system is given by M(X)=1+ϵEI^Xd2ϕ(X)/dX2. In the current work the first-order correction to the mode shapes is found. Therefore, it is assumed that ϕ(X)=ϕ0(X)+ϵϕ1(X). Choose an arbitrary point XPi∈P^, consider the difference in the value of the moment to the left and right of this point, and use the small parameter ϵ definition. (31)MXPi--MXPi+=EIXPi+-EIXPi-EIXPi--EI∗EI∗2d2ϕ0XPidX2+1EI∗d2ϵϕ1XPi+dX2It is seen in (31) that if the bending stiffness of the system is continuous, then the moment will also be continuous. However, if the bending stiffness of the system is not continuous, then the moment of the system may or may not be continuous. In (31) it is also seen that the mass per unit length does not play a role in the continuity of the moment. In the stepped beam case when τ is not equal to one, the bending stiffness is not continuous and this causes the discontinuous moment in Figure 4(c). In the simulations where τ=1 and γ=1.5, the bending stiffness is continuous and therefore the moment in Figure 4(d) is continuous.Using the same approach, the jump in shear at a pointXPi∈P^ can be determined. The shear in the system at a given point X is S(X)=d/dXM(X) and the result of the difference is given by (32)SXPi--SXPi+=dEIXPi+dX-dEIXPi-dXEIXPi--EI∗EI∗2d2ϕ0XPidX2-1EI∗d2ϵϕ1XPi+dX2+EIXPi--EIXPi+1EI∗d3ϵϕ1XPi+dX3-EIXPi--EI∗EI∗2d3ϕ0XPidX3-dEIXPi-/dXEI∗2d2ϕ0XPidX2It can be seen from (32) that if the bending stiffness of the system and its first derivative are continuous, then the shear will also be continuous. However, if the bending stiffness of the system is not continuous or the first derivative is not continuous, then the shear of the system may or may not be continuous. In (32) it is seen that the mass per unit length does not play a role in the continuity of the shear. In the stepped beam case when τ is not equal to one, the bending stiffness and the first derivative are both not continuous and this is the reason for the discontinuous shear in Figure 4(e). In the simulations where τ=1 and γ=1.5 the bending stiffness and its first derivative are continuous and therefore the shear in Figure 4(f) is continuous.Table1 presents the absolute errors in predicting the moment and shear for the various boundary conditions. It is clearly demonstrated in Table 1 that the norm minimization methods provide significantly better estimations for the moment and shear when compared to the first-order gradient method. This is expected for the first-order gradient method since there are fewer corrections to the mode shapes. The error in using the first-order gradient method is typically about 5 times larger when compared to the norm minimization methods, with the errors being up to 30 times higher in the case of the shear for the CF boundary conditions. This is a major contribution of the current work and clearly shows the importance of including the first-order correction to the mode shapes.Table 1
Absolute percentage of error in moment and shear prediction of perturbation methods compared to analytical results.
Method Moment atl0/l Shear atl0/l Moment atl0/l Shear atl0/l Moment atl0/l Shear atl0/l 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + 1 / 4 - 1 / 4 + Clamped-Clamped Clamped-Free Free-Free 1st Grad. 62.94% 44.40% 21.44% 17.84% 9.90% 35.02% 11.05% 33.43% 32.74% 0.89% 31.64% 2.54% 2nd Grad. 15.99% 5.32% 4.55% 0.67% 2.33% 7.73% 2.42% 6.87% 10.24% 0.32% 9.53% 0.82% L 6-Diff. - - - - 0.64% 1.89% 0.58% 1.36% - - - - L 6-Pert. - - - - 0.56% 1.60% 0.49% 1.09% - - - - L ∞-Diff. 12.48% 9.04% 3.28% 2.06% - - - - 6.28% 0.22% 5.64% 0.52% L ∞-Pert. 12.48% 9.04% 3.28% 2.06% - - - - 6.28% 0.22% 5.64% 0.52%The final results presented for the stepped beam system are the frequency error results for a variable step location. Figures5 and 6 consider the cases when either γ or τ is equal to 0.9 or 1.1, respectively. CC, CF, and FF boundary conditions are all considered. As was the case for variable τ and γ, an optimal p-value is determined for the norm minimization methods for each system setup. Note that since now a different parameter is being varied, l0/l, the optimal p-value is not expected to be the same as the case for variable τ and γ.Figure 5
Frequency error results with a decrease inγ or τ for variable l0/l.
(a)
τ = 0.9, γ=1, CC (b)
τ = 0.9, γ=1, CF (c)
τ = 0.9, γ=1, FF (d)
τ = 1, γ=0.9, CC (e)
τ = 1, γ=0.9, CF (f)
τ = 1, γ=0.9, FFFigure 6
Frequency error results with an increase inγ or τ for variable l0/l.
(a)
τ = 1.1, γ=1, CC (b)
τ = 1.1, γ=1, CF (c)
τ = 1.1, γ=1, FF (d)
τ = 1, γ=1.1, CC (e)
τ = 1, γ=1.1, CF (f)
τ = 1, γ=1.1, FFIn Figures5 and 6, it can be seen that for a 10% increase or decrease in the bending stiffness or mass per unit length the proposed methods for choosing the reference values in the perturbation theory lead to very small error values. The first-order gradient method produced significantly larger error values when compared to the other proposed methods. For all the presented simulations, the largest sum of absolute error of the frequencies for the first-order gradient method is 0.921%, for the second-order gradient method is 0.050%, for the difference norm minimization is 0.046%, and for the perturbation norm minimization is 0.046%. This indicates that the perturbation theory detailed in the current work provides very accurate estimations for the natural frequencies of the system, even in the case of fairly large changes to the bending stiffness or the mass per unit length.More specifically from Figures5 and 6, the cases when the bending stiffness is varied produce typically larger errors in the frequencies when compared to the cases when the mass per unit length is varied. This is similar to the observations in the variables τ and γ results of Figures 2 and 3, as well as the moment and shear results. Additionally, the error results for CF boundary conditions are typically smaller than for CC and FF boundary conditions given the same system parameters. A similar observation occurred in the previous variables τ and γ results.Finally, in Figures5 and 6, it is seen that the largest errors occur near l0/l=1/2 and the metric value is 0 when l0/l=0 and l0/l=1. When the step location in the beam is near the middle of the system, it becomes more difficult to determine a single reference value that accurately represents the bending stiffness or mass per unit length of the system. Therefore there will be a larger error associated with the perturbation theory methods. When l0/l=0 or l0/l=1, the entire system has a single value for the bending stiffness and mass per unit length. In this case, the system is a homogeneous beam and the perturbation theory solution is exactly the analytical solution. Since the perturbation theory solution is exactly the analytical solution, the error in the frequencies is 0.


## 5.2. Periodically Grooved Beam

The second system to which the perturbation theory developed in the current work is applied is a grooved beam, as shown in Figure7. In this example, at the middle point of the groove the system reaches a thickness value of ch, where h is the thickness of the beam where there is no groove. Thus, when c>1, the beam becomes thicker in the grooved section as in Figure 7(a); when 0<c<1, the beam becomes thinner in the grooved section as in Figure 7(b). If c=1 then there are no grooves in the system and the beam is homogeneous. A system with multiple grooves consists of multiple single groove elements from either Figure 7(a) or Figure 7(b) assembled side by side.Figure 7
Grooved beam system schematics.
(a)
Grooved element withc>1 (b)
Grooved element with0<c<1Consider a single grooved section shown in Figures7(a) and 7(b). In this case, the bending stiffness and mass per unit length of the system are continuous and the derivatives are discontinuous. Let L denote the total length of a single grooved element. The bending stiffness and mass per unit length of a single grooved element can be expressed as (33)EIx=1,0≤x≤L41+c-1x-L/4L/43,L4≤x≤L2c+1-cx-L/2L/43,L2≤x≤3L41,3L4≤x≤L,ρAx=1,0≤x≤L41+c-1x-L/4L/4,L4≤x≤L2c+1-cx-L/2L/4,L2≤x≤3L41,3L4≤x≤LAs in the previous numerical example, the sum of the absolute error in the frequencies is used to provide a measure of how well the perturbation theory compares to an exact solution. The first 10 nonrigid body modes are considered and compared to a finite elements analysis (FEA). The finite element analysis employed in the current work follows the approach outlined in [15]. The error results for a system with variable c3, multiple boundary conditions applied, and varying number of total grooves in a fixed length system are presented in Figures 8–11. As in the stepped beam example, for the norm minimization methods multiple p-values were considered and the errors calculated for the variable parameter are being analyzed. Then, the p-value corresponding to the smallest area under the error curves was determined and plotted in the figures.Figure 8
Frequency error results for variablec using the first-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 9
Frequency error results for variablec using the second-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 10
Frequency error results for variablec using the norm minimization of the difference.
(a)
CC (b)
CF (c)
FFFigure 11
Frequency error results for variablec using the norm minimization of the perturbation.
(a)
CC (b)
CF (c)
FFThe smallestc3 value considered produces a 50% and 20.6% decrease in the bending stiffness and mass per unit length of the system, respectively. The largest c3 value considered produces a 50% and 14.5% increase in the bending stiffness and mass per unit length of the system, respectively. The ratio of the change in the bending stiffness from the centre of the groove to a nongrooved section, which is exactly c3, is larger than the ratio of the change in the mass per unit length from the centre of the groove to a nongrooved section, which is exactly c. For this reason the c3 value is used as the independent variable in the figures.In Figures8–11 it is seen that when c3=1 the errors are 0. This is expected since in this case the beam is homogeneous and the perturbation theory and FEA produce identical results. Additionally in Figures 8–11, it is observed that the first-order gradient minimization method always produces the largest error of the four methods. This is expected since in this case there is only a single correction to the natural frequencies, whereas in the other methods there are two corrections to the natural frequencies and a correction to the mode shapes. These additional corrections provide an increased accuracy. The next largest error is the second-order gradient results in Figure 9, then the norm minimization results in Figures 10 and 11. The largest sum of absolute error in the frequencies across all boundary conditions for the first-order and second-order gradient methods and the difference and perturbation norm minimization methods is 19.835%, 4.501%, 0.921%, and 0.9353%, respectively.As well in Figures8–11, the results for the errors of a given system are not symmetric about c3=1 and there are larger errors when c3=0.5 compared to when c3=1.5. To analyze this, consider a c3 value a distance Δ from c3=1. Then the points to the left and right of c3=1 that are equidistant have a value of c=1±Δ3. For a given Δ, the percentage change in the c value from c=1 is larger when Δ is subtracted in the radicand than when it is added. This indicates that in the case that we are to the left of c3=1 the effect on the mass per unit length is greater than for an equidistant point to the right of c3=1. Due to this, there are higher errors for the cases when c3<1 than when c3>1 and the results for the error are not symmetric about c3=1.Most importantly, the results of Figures8–11 demonstrate the accuracy for a large range of c3 values when using two corrections to the frequencies and a single correction to the mode shapes. In particular, using the norm minimization methods to select the reference values in the perturbation theory, a very small error is produced even when there is up to a 50% increase or decrease in the bending stiffness of the system. Furthermore, the error curve for these methods is quite flat for c3 values near 1 when compared to the gradient minimization methods. These properties highlight the robustness of the proposed methods, in particular the norm minimization approaches.Additionally, it is seen in Figures8 and 9 that the errors for the frequencies are smallest for a system with 1 groove and increase as the number of grooves increases for all the boundary conditions. The one exception to this was for the second-order gradient minimization method with FF boundary conditions for c3>1 where a system with 3 grooves produced slightly larger errors than a system with 5 grooves. This behaviour of increasing frequency error for increasing number of grooves indicates that a system with more rapidly varying bending stiffness and mass per unit length will present high errors when using a perturbation theory employing the gradient minimization method.Next, it is observed in Figures10 and 11 for the difference and perturbation norm minimization methods that there is no consistent pattern with respect to errors in the frequencies and number of grooves. For the norm minimization methods, a single value for EI∗ and ρA∗ is found, whereas for the gradient minimization method a reference value is found for each mode of the system. This is likely the main contributing factor for the lack of pattern in the errors for the norm minimization methods compared to the gradient minimization method. It can be seen, however, that there is a consistent behaviour amongst the results for the difference and perturbation norm minimization.Consider the expressions for the moment and shear using the perturbation theory presented in (31) and (32). For the grooved system in the current example, the bending stiffness of the system is continuous over the domain and discontinuous at the start, centre, and end of each groove. Due to this, it is expected that the perturbation theory results will produce a moment that is continuous and a shear that is discontinuous. As an example, Figure 12 presents the results for CF boundary conditions with c3=1.1 and 1 groove in the system. The expected behaviours for the moment and shear are clearly seen in these results and the advantage of including additional corrections to the mode shape is particularly highlighted in the shear results of Figure 12(c).Figure 12
Mode shape displacement, moment, and shear for clamped-free boundary conditions.
(a)
Displacement (b)
Moment (c)
ShearFinally, the errors for a system with a fixedc3 value and variable number of grooves in a system with a fixed length are presented in Figures 13–16. Multiple c3 values are considered. As before, for the norm minimization methods multiple p-values were considered and the errors calculated for the variable parameter are being studied. Then, the p-value corresponding to the smallest area under the error curves was determined and plotted in the figures. It is not expected that the optimal p-values for systems considering variable number of grooves are the same as the case for variable c3 values.Figure 13
Frequency errors for variable groove number and first-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 14
Frequency errors for variable groove number and second-order gradient minimization.
(a)
CC (b)
CF (c)
FFFigure 15
Frequency errors for variable groove number and norm minimization of difference.
(a)
CC (b)
CF (c)
FFFigure 16
Frequency errors for variable groove number and norm minimization of perturbation.
(a)
CC (b)
CF (c)
FFIn Figures13–16, it is observed that the first-order gradient minimization method always produces the largest error of the four presented methods. This is expected since in this case there is only a single correction to the natural frequencies, whereas in the other methods there are two corrections to the natural frequencies and a correction to the mode shapes. These additional corrections provide an increased accuracy. Also the figures show that the norm minimization methods produce smaller errors when compared to the second-order gradient method. In addition, the norm minimization methods are seen to provide very accurate results for all the number of grooves considered when c3 is near one. The largest sum of absolute error in the frequencies across all boundary conditions for the first-order and second-order gradient methods and the difference and perturbation norm minimization methods is 23.436%, 5.709%, 0.884%, and 0.888%, respectively.Overall in Figures13–16, there is a clear pattern in the levels of error reported for all the boundary conditions and all the proposed methods for determining the reference value. For a given system, the errors reported are largest for c3=0.5, the next largest are for c3=1.5, the second smallest are for c3=0.9, and the smallest errors are for c3=1.1. It was previously shown that for equidistant points from c3=1 that the c3<1 value will report higher errors than the c3>1 value. Indeed this is observed throughout the presented results for variable number of grooves. In Figures 13–16, it is not expected that the error is 0 for any of the simulations since the case of c3=1 is not considered.It is seen in Figures13 and 14 for the gradient minimization methods that the errors for the frequencies increase as the number of grooves in the system increases and then reaches an upper bound. This upper bound is achieved for a smaller number of grooves by the second-order gradient minimization compared to the first-order method due to the additional corrections considered in the perturbation theory. This is an advantage of the second-order method. This overall increase in the errors up to a given number of grooves in the system is consistent with the behaviour that was observed in the results considering variable c3. For both the first- and second-order gradient minimization methods, the behaviour observed was consistent across all the boundary conditions.Next, it is witnessed for the difference and perturbation norm results in Figures15 and 16 that the behaviour is quite different than for the gradient minimization approach. This is due to the manner in which the reference values are determined in each of these cases. For the gradient minimization a reference value EI∗ and ρA∗ is determined for each mode of the system, whereas the norm minimization methods predict only a single value for all the modes of the system. It can be seen, however, that similar behaviours are observed between the results of the difference and perturbation norm minimization methods. For CC boundary conditions the frequency errors are fairly constant, with a slight overall increase for the case of c3=0.5. For CF and FF boundary conditions, the error in the frequencies noticeably decreases as the number of grooves increases for the cases of c3=0.5 and c3=1.5. It should be noted that the pattern observed is due to the manner in which the p-value was determined. If, for example, the optimal p-value was determined for each data point in the figures by minimizing the errors for the given data point, then the pattern of the errors for variable number of grooves would be different. A similar approach for minimizing the errors for each given data point can be applied for the figures where the effects of a variable c3 are discussed, which will result in a different pattern for Figures 10 and 11. Regardless of which p-value is chosen for the system optimization, the norm minimization method consistently gives better error values compared to the gradient minimization method.A final note is made with regard to the results observed in Figures13 to 16 for variable number of grooves considering multiple c3 values in comparison to the results observed in Figures 8 to 11 for variable c3 value considering multiple number of grooves in the system. For the gradient minimization methods, the method for determining the reference values is prescribed and the results are seen to be consistent across all the simulations. That is, increasing the number of grooves in the system from 1 to 7 produces larger errors, as well as increasing the c3 value such that it is further from c3=1. However, it is seen that slightly different behaviours occur for the results of the difference and perturbation norm minimization. This is due to different optimal p-values being determined for the cases of variable c3 and variable number of grooves. As such, it is not expected that the pattern of the errors is consistent across all the numerical simulations. As previously mentioned, the choice of p-value will have an effect on the behaviour of the errors. Choosing two sets of optimal p-values that are based on which type of variation is allowed in the system demonstrates how the results can be made more accurate than if a single value for p was used throughout all the simulations.


## 6. Conclusion

This paper presents a perturbation theory with two corrections to the natural frequencies and one correction to the mode shapes applied to an Euler-Bernoulli PDE. A contribution of the current paper is extending the work in [21] to allow for piecewise continuous spatial properties in the system. Multiple methods for choosing the reference values from which the bending stiffness and mass per unit length are perturbed are considered. It was shown that the method in which these reference values are chosen plays a significant role in the accuracy of the perturbation theory results. The first reference value selection method proposed in the current work involves minimizing the gradient of the expression for the natural frequencies. The additional methods involve minimizing the Lp norm of the difference and the perturbation of the system from the reference value.Two numerical examples were considered to investigate the accuracy of the proposed methods and are compared to a perturbation theory with only a single correction to the natural frequencies. This was done to demonstrate the significant increase in accuracy of the perturbation theory by including additional corrections to the frequencies and mode shapes. The first numerical example is a system with a single step in the bending stiffness and mass per unit length and the second system has periodic grooves. Multiple boundary conditions and various system parameters were considered. It was shown that the norm minimization methods produced the most accurate results when compared to an analytical solution for the stepped beam and a finite element analysis for the grooved beam. Additionally, the correction to the mode shape was shown to greatly increase the accuracy of the estimates for the moment and shear in the system. The results showed that the norm minimization methods provides a manner in which to choose the reference values of the system such that the perturbation theory solution is also quite robust to variable system parameters.



---

*Source: 4627865-2018-06-27.xml*