Reference Value Selection in a Perturbation Theory Applied to Nonuniform Beams

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


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][2][3][4][5][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][13][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][8][9][10][11][12][13][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][17][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   -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][2][3][4][5][6].

Problem Formulation
2.1.Perturbation Theory.In this work, the free transverse vibration of an Euler-Bernoulli (EB) beam model with spatially dependent bending stiffness () and mass per unit length () is considered.The functions () and () are assumed to have derivatives of all orders except on the finite sets of points   ⊂ [0, ] and   ⊂ [0, ], respectively, where  is the total length of the system.That is, () and () are  ∞ -functions on [0, ] \   and [0, ] \   , respectively.The partial differential equation (PDE) for the free transverse vibration of the system is for 0 <  <  and  > 0. The boundary conditions at the endpoints  * = 0 and  are The 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  * and  * 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  and the time scale  2 √ * / * , and let , , and  denote dimensionless quantities.Then, assume a separable solution of the form (, ) = ()  .Assuming the variations of the spatially dependent functions () and () from the respective reference values are small, the ODE can be written as where  is a small parameter and  Ê() = (()− * )/ * and  ρ() = (() −  * )/ * .Applying the same steps to the boundary conditions in (2a) and (2b) yields In 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 () = ∑ ∞ =0     () and  = ∑ ∞ =0     , 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   () and   .

Mass Normalization Condition.
In a vibrations analysis, it is typical to mass normalize the mode shapes.Denote the mass normalized mode shapes by () and the dimensionless mass normalized mode shapes by ().In the current work, the nondimensionalized mass normalization condition is given by 1 = ∫ 1 0 {1+ ρ()} 2 ()d.Using the perturbation theory expansion for the mode shape the zeroth-order, O(1), and first-order, O(), conditions are

O(1)
Problem.From the perturbation theory, the ODE for the O(1) problem is d 4  0 /d 4 −  2 0  0 = 0 for 0 <  < 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 0 (0) = 0, d 0 (0) d = 0, The 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 th mode of the system, a subscript  is added to the terms in the pertur-Shock and Vibration bation theory expansion of the mode shape and frequency.

O(𝜖)
Problem.The ODE for the O() problem is As is typically done in perturbation theory, a solvability condition is applied to ensure that a solution to (10) exists, Integration by parts and using the fact that CC, CF, and FF boundary conditions are considered in this work, a solution for  1 can be found as Next the ODE in (10) is solved for  1 ().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  =   ∪   ⊂ [0, ].Define P as the set of all the points in  divided by ; thus, P ⊂ [0, 1].Suppose the size of the set P is  − 1 and order the points in P such that 0 <   1 <   2 < ⋅ ⋅ ⋅ <   −1 < 1.Then each of the  intervals on which both Ê and ρ are of class  ∞ can be denoted using a superscript (), i.e., using the notation Ê () and ρ () .The th interval is defined as [  −1 ,    ] with   0 = 0 and    = 1.The solution for  1 () in (10) can be determined over each of the  intervals.The general solution for the th interval, where  () 1, () is a particular solution, is Boundary 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 solution  1 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    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.
The boundary conditions for the O() 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  (1)  1 (0) = 0, d 1 (0) d = 0,  () 1 (1) = 0, (1)  1 (0) = 0, 1 (0) First, the boundary conditions at  = 0 are applied, then the continuity conditions applied sequentially beginning from the left end of the system,   1 , up to the right end of the system,   −1 .After these have been applied, a single boundary condition at  = 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  1 , the fourth and final boundary condition will automatically be satisfied.The final expression for the th section of the mode shape, with  ≤ 1 ≤ , is given by where  ũ() for CC, CF, and FF boundary conditions is In (17) the function ũ() 1 () is calculated in terms of the O(1) solution  0 () and the functions  () 1, () with 1 ≤  ≤ .The complete expression for ũ() 1 () is given in Appendix A. Applying the mass normalization condition in (6), the final constant  (1)  1 is determined and the solution on the th interval is denoted by  ()  1 ().The final coefficient is given by A solution to (20) exists when the solvability condition of the integral over the dimensionless domain of the product of  0 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  2 can be found as

Dimensional Solution and Dependence on Reference
Values.In (3), a small parameter was introduced for the perturbation theory via  Ê() = (() −  * )/ * and  ρ() = (() −  * )/ * .Using these, expressions for the frequency and mode shape perturbation solution can be found in terms of the reference values  * and  * .In addition, these expressions will depend on known quantities such as the bending stiffness () and the mass per unit length (). where Shock and Vibration The complete expressions of  1 () and  2 () 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 are 1/√kg.For the choice of length and time scales used in the current work, the mass scale is taken to be  * .To reintroduce dimensions, the mode shapes obtained from the perturbation theory must be divided by the square root of the mass scale.

Reference Selection Method 1:
Minimizing   Norm of Difference.The first proposed approach for determining the reference values is to minimize the   norm of the difference between the reference values and its corresponding spatial function.Mathematically, where  is a positive integer.Minimizing the   norm of the difference results in a single value for  * and  * for all the modes.In general the minimization procedure listed in ( 26) is performed numerically.

Reference Selection Method 2:
Minimizing   Norm of Perturbation.The second proposed approach for determining the reference values is to minimize the   norm of perturbation terms in (3).Mathematically, where  is a positive integer.Minimizing the   norm of the perturbation results in a single value for  * and  * 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, ].In (27) the domain [0, 1] is used since this is the domain in consideration when applying the perturbation theory.

Reference Selection Method 3: Gradient Norm Minimization.
The final proposed approach to determine the reference values is minimizing the  2 (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  * and  * .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  * and  * 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 minimize: ∇     2 subject to: min and the  2 norm of the gradient for the perturbation theory frequency of ( 25) is given by For 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  *  and  *  .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  *  = 2 1 / 0 and  *  = 2 2 / 0 .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].

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 constant () and () 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 both () and () are not constant over the domain.It is useful to rewrite the dimensional frequency in (25) [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.

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 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.
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 -value is determined and presented in the figures.This is achieved for each system setup by considering multiple -values and calculating the errors for the variable parameter being analyzed.Then, the 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.
When 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 .

Shock and Vibration
Consider the results in Figures 2 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  0 / = 1/4 and  0 / = 3/4.The results for  0 / = 1/4 given variable  are identical to the results for  0 / = 3/4 given 1/, and vice versa.For example, the errors reported for  = 3/2 given  0 / = 1/4 are identical to the errors reported for  = 2/3 given  0 / = 3/4.The same results for symmetry hold true for variable .
Finally, consider the -values that produce the smallest area under the error curves for the norm minimization techniques in Figures 2 to 3. For variable  the -value used always approaches infinity with the exception of the CF results.For variable  the -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 -value should be considered, whereas when the perturbation in the system is due to a change in the mass per unit length, a -value of less than 10 should be employed.
Next, dimensional mode shapes associated with the fundamental frequency of the system are presented in Figure 4.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  0 / = 1/4 is used and the values of  and  are 1 unless otherwise specified.The -values used are taken from the corresponding previous results for variables  and .
For the results presented in Figure 4, 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  0 / = 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  in the system is given by () = {1 +  Ê()}(d 2 ()/d 2 ).In the current work the first-order correction to the mode shapes is found.Therefore, it is assumed that () =  0 () +  1 ().Choose an arbitrary point    ∈ P, consider the difference in the value of the moment to the left and right of this point, and use the small parameter  definition.

𝑀 (𝑋
It 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 point    ∈ P can be determined.The shear in the system at a given point  is () = (d/d)() and the result of the difference is given by It 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.
Table 1 presents the absolute errors in predicting the moment and shear for the various boundary conditions. 1 that the norm minimization methods provide significantly better estimations for the moment and shear when compared to the firstorder 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.

It is clearly demonstrated in Table
The final results presented for the stepped beam system are the frequency error results for a variable step location.Figures 5 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 -value is determined for the norm minimization methods for each system setup.Note that since now a different parameter is being varied,  0 /, the optimal value is not expected to be the same as the case for variable  and .
In Figures 5 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 firstorder 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 Figures 5 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 Figures 5 and 6, it is seen that the largest errors occur near  0 / = 1/2 and the metric value is 0 when  0 / = 0 and  0 / = 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  0 / = 0 or  0 / = 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.

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 Figure 7.In this example, at the middle point of the groove the system reaches a thickness value of ℎ, where ℎ is the thickness of the beam where there is no groove.Thus, when  > 1, the beam becomes thicker in the grooved section as in Figure 7(a); when 0 <  < 1, the beam becomes thinner in the grooved section as in Figure 7(b).If  = 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.
Consider a single grooved section shown in Figures 7(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  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  As 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  3 , 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 values were considered and the errors calculated for the Shock and Vibration variable parameter are being analyzed.Then, the -value corresponding to the smallest area under the error curves was determined and plotted in the figures.The smallest  3 value considered produces a 50% and 20.6% decrease in the bending stiffness and mass per unit length of the system, respectively.The largest  3 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  3 , 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 .For this reason  3 value is used as the independent variable in the figures.
In  it is seen that when  3 = 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. Tese 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 firstorder 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 Figures 8-11, the results for the errors of a given system are not symmetric about  3 = 1 and there are larger errors when  3 = 0.5 compared to when  3 = 1.5.To analyze this, consider a  3 value a distance Δ from  3 = 1.Then the points to the left and right of  3 = 1 that are equidistant have a value of  = 3 √ 1 ± Δ.For a given Δ, the percentage change in the  value from  = 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  3 = 1 the effect on the mass per unit length is greater than for an equidistant point to the right of  3 = 1.Due to this, there are higher errors for the cases when  3 < 1 than when  3 > 1 and the results for the error are not symmetric about  3 = 1.
Most importantly, the results of Figures 8-11 demonstrate the accuracy for a large range of  3 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  3 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 Figures 8 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  3 > 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 Figures 10 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  * and  * 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  3 = 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).Finally, the errors for a system with a fixed  3 value and variable number of grooves in a system with a fixed length are presented in Figures 13-16.Multiple  3 values are considered.As before, for the norm minimization methods multiple -values were considered and the errors calculated for the variable parameter are being studied.Then, the value corresponding to the smallest area under the error curves was determined and plotted in the figures.It is not expected that the optimal -values for systems considering variable number of grooves are the same as the case for variable  3 values.

Shock and Vibration
In Figures 13-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  3 is near one.The largest sum of absolute error in the frequencies across all boundary conditions for the first-order and secondorder gradient methods and the difference and perturbation norm minimization methods is 23.436%, 5.709%, 0.884%, and 0.888%, respectively.
Overall in Figures 13-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  3 = 0.5, the next largest are for  3 = 1.5, the second smallest are for  3 = 0.9, and the smallest errors are for  3 = 1.1.It was previously shown that for equidistant points from  3 = 1 that the  3 < 1 value will report higher errors than the  3 > 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  3 = 1 is not considered.
It is seen in Figures 13 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  3 .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 Figures 15 and 16  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  * and  * 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  3 = 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  3 = 0.5 and  3 = 1.5.It should be noted that the pattern observed is due to the manner in which the -value was determined.If, for example, the optimal -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  3 are discussed, which will result in a different pattern for Figures 10 and 11.Regardless of which -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 Figures 13 to 16 for variable number of grooves considering multiple  3 values in comparison to the results observed in Figures 8 to 11 for variable  3 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  3 value such that it is further from  3 = 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 -values being determined for the cases of variable  3 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 -value will have an effect on the behaviour of the errors.Choosing two sets of optimal -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  was used throughout all the simulations.

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

Figure 2 :
Figure 2: Frequency error results for variable .

Figure 3 :
Figure 3: Frequency error results for variable .

Figure 4 :
Figure 4: Mode shape displacement, moment, and shear for CC boundary conditions.

Figure 5 :
Figure 5: Frequency error results with a decrease in  or  for variable  0 /.

Figure 6 :
Figure 6: Frequency error results with an increase in  or  for variable  0 /.

Figure 9 :
Figure 9: Frequency error results for variable  using the secondorder gradient minimization.

Figure 10 :
Figure 10: Frequency error results for variable  using the norm minimization of the difference.

Figure 11 :
Figure 11: Frequency error results for variable  using the norm minimization of the perturbation.

Figure 12 :
Figure 12: Mode shape displacement, moment, and shear for clamped-free boundary conditions.

Figure 13 :
Figure 13: Frequency errors for variable groove number and firstorder gradient minimization.

Figure 14 :
Figure 14: Frequency errors for variable groove number and second-order gradient minimization.

Figure 15 :
Figure 15: Frequency errors for variable groove number and norm minimization of difference.

Figure 16 :
Figure 16: Frequency errors for variable groove number and norm minimization of perturbation.

Table 1 :
Absolute percentage of error in moment and shear prediction of perturbation methods compared to analytical results.