## Abstract

The Wentzel–Kramers–Brillouin (WKB) approximation is applied for asymptotic analysis of time-harmonic dynamics of corrugated elastic rods. A hierarchy of three models, namely, the Rayleigh and Timoshenko models of a straight beam and the Timoshenko model of a curved rod is considered. In the latter two cases, the WKB approximation is applied for solving systems of two and three linear differential equations with varying coefficients, whereas the former case is concerned with a single equation of the same type. For each model, explicit formulations of wavenumbers and amplitudes are given. The equivalence between the formal derivation of the amplitude and the conservation of energy flux is demonstrated. A criterion of the validity range of the WKB approximation is proposed and its applicability is proved by inspection of eigenfrequencies of beams of finite length with clamped–clamped and clamped-free boundary conditions. It is shown that there is an appreciable overlap between the validity ranges of the Timoshenko beam/rod models and the WKB approximation.

## 1. Introduction

The Wentzel–Kramers–Brillouin (WKB) approximation is a well-established and recognized tool in theoretical physics in general, and in quantum mechanics in particular. Its application for solving the Schrödinger equation provides simple-structured solutions that describe the motion of particles and phenomena such as tunnelling. Examples can be readily found in classical compilations and textbooks such as Fröman & Fröman [1] or Bender & Orszag [2]. The fundamental feature of the WKB method is its ability to approximate the solutions of linear differential equations with varying coefficients. This unique characteristic of the method immediately suggests its suitability for describing wave motion in inhomogeneous media and consequently the method has found its way to the field of acoustics. For problems of wave propagation in non-uniform acoustic ducts the method has been applied by for instance Rienstra [3], Cooper & Peake [4] and Brambley & Peake [5]. The method has also been applied in elastodynamics and, more specifically, for analysis of wave motion of planar curved elastic strips by Scott & Woodhouse [6], guided waves in a circular elastic annulus by Gridin *et al.* [7], and straight non-uniform rods in plane stress by Rosenfeld & Keller [8].

Another class of problems in elasticity where the WKB method is applicable, but, to the best of the authors knowledge has not yet been used profoundly, is the wave propagation in non-uniform rods in the framework of the classical ‘plane cross-sectional’ model. The pioneering work along these lines has been done by Pierce [9], who used the WKB method to determine the wavenumbers in straight inhomogeneous Bernoulli–Euler and Timoshenko beams. However, the amplitude variation was determined from the energy conservation condition, rather than by means of the formal WKB approximation. An expression for the eigenfrequencies of a cantilevered Bernoulli–Euler beam is presented in the paper, with no numerical results being reported. In a paper by Langley [10], the propagation of flexural waves in a straight Bernoulli–Euler beam was studied by deriving a short-wave perturbation solution similar to a WKB solution. The short-wave perturbation solution found in Langley [10] was derived after recasting the governing equation to a first-order matrix vector system and transform this to ‘wave coordinates’ based on the eigenvectors and eigenvalues of the system. More specifically the case analysed was the reflection and transmission coefficients for a non-uniform pre-compressed beam resting on an elastic foundation with non-uniform properties.

A relatively low interest in the WKB approximations in the classical rod theories can be explained by an obvious contradiction between the limitations of the two. Namely, the WKB approximation is, in essence, a short-wave or high-frequency method, whereas the classical rod theories are formulated for long-wave or low-frequency wave motion. More specifically, the Timoshenko beam model may be regarded as accurate up to and slightly above the shear wave cut-on [11], whereas the validity range for the Bernoulli–Euler beam model is much narrower. This naturally raises the question on whether there is an appreciable overlap of frequencies, where both the plane cross-sectional hypothesis and the WKB approximation are valid. To answer this question, results determined by the WKB method should be compared with results determined by a substantially different method, for example, by numerical solutions. In the above-mentioned references on rods, such a comparison has not been done and one of the goals of this paper is to accomplish this. This is done by comparing eigenfrequencies determined by the WKB method and numerical methods for various rods when specialized to different boundary conditions. Determining eigenfrequencies of non-uniform beams as done in this paper is, however, not the main motivation of using the WKB method, but rather eigenfrequencies are regarded as a convenient measure to assess the accuracy of the WKB method when applied to plane cross-sectional problems.

Analysis of elastic waves in non-uniform beams by the WKB approximation is useful in many practical aspects. Asymptotic analysis is, in general, known as a powerful tool to reveal the underlying physics of complicated problems and to provide some approximate yet simple and useful relationships between the parameters of a given problem and its solution. These insights can greatly facilitate the engineer during the design process of, say, a suspension system composed of non-uniform beams. The important information on wave propagation, cut-on/cut-off frequencies, etc., of a non-uniform waveguide, available from the WKB method is quite unique as no standard numerical method provides that.

The classical application of the WKB approximation in quantum mechanics (or in linear acoustics) is concerned with solving a single partial differential equation of the second order with varying coefficients. In elastodynamics, the approximation is applied to solve a system of two differential equations, which are of the second order each. The motion of curved rods in the framework of Timoshenko formulation are generally governed by six coupled differential equations with varying coefficients and therefore the WKB solution of three coupled equations, representing a planar curved corrugated Timoshenko rod, reported in this paper contains aspects of novelty in its own right.

The study presented in this paper is devoted to demonstrate how the WKB method is generalized to, in a formal manner, solve systems of differential equations with slowly varying coefficients and to demonstrate that there is an appreciable overlap between the high-frequency WKB method and the low-frequency theory of rods, which employs the plane cross-sectional hypothesis. The application of the WKB approximation will be exemplified on the following non-uniform elastic members: (i) *straight Rayleigh rod*, (ii) *straight Timoshenko rod* and (iii) *curved Timoshenko rod*. This choice of increasingly advanced models will accomplish the demonstration of how simple and systematic the WKB approximation in fact is. In §4, the accuracy of the WKB solutions when applied to elasticity problems based on the plane cross-sectional hypothesis is gauged. This is done by the aforementioned comparison of eigenfrequencies.

## 2. Models of non-uniform beams

Consider a rod with a circular cross section with a varying diameter and define the axial length scale of inhomogeneity as *l** (* denotes dimensional quantities), i.e. this is the length over which the diameter changes significantly. More specifically, the diameter may be defined as
*d*′* is the change in diameter and *s** is the axial coordinate. By scaling lengths using *ϵ*≪1. Thus, the diameter is dependent on a slow coordinate *d**(*ϵs*), and likewise are the area and moment of inertia: *A**(*ϵs*)=(*π*/4)*d**^{2}(*ϵs*), *I**(*ϵs*)=(*π*/64)*d**^{4}(*ϵs*). This is the typical way of stating the requirement of slow variation of a waveguide even though nothing, at this stage, is said about the wavelength. It will be shown later on that stating *ϵ*≪1 is not by itself sufficient to ensure accuracy of the WKB approximation.

The key assumption in classical rod theories is the assumption that a cross section remains plane in the course of deformation. The simplest of these models is the classical Bernoulli–Euler model, in which the cross section is additionally assumed to remain perpendicular to the middle axis during deformation [12]. This has the immediate consequence that the deformation of the beam is completely given by the lateral deflection *u**(*s**,*t**) of the centreline. Owing to this assumption, the elastic energy stems only from the deflection caused by the bending moment, which is given by *Q**(*s**,*t**)=−*E***I**(*s*)(∂^{2}*u**(*s**,*t**)/∂*s**^{2}), where *E** is Young's modulus, thus the elastic energy density is *L** is the length of the beam. The kinetic energy of a beam under deformation is
*ρ** is the mass density. The first term in (2.3) is the conventional Bernoulli–Euler inertia term and the second term is the rotary inertia correction term introduced by Rayleigh [13]. For that reason, this model will henceforth be referred to as the Rayleigh beam model. The governing equation of such a beam can be found by setting up the action integral and using Hamilton's principle
*e*^{−iω*t*} and then scale lengths by *β*(*s**,*t**). Thus, the shear force in a Timoshenko beam cannot be determined only from the deflection of the centreline, but rather by the difference between the slope of the centreline and the rotation of a cross section
*G** is the shear modulus, which is rewritten by *G**=*E**/(2(1+*ν*)) where *ν* is Poisson's ratio which is also dimensionless. The internal moment is given by
*u**(*s**,*t**) and *β**(*s**,*t**), so contrary to the Bernoulli–Euler and Rayleigh beam models the Timoshenko beam model is governed by two equations. Again by adopting the temporal dependence *e*^{−iω*t*} and perform non-dimensional scalings the following equations are found:
*p*=λ/(2(1+*ν*)) has been introduced for convenience. A complete derivation of these equations is available in Shames & Dym [12]. The curved Timoshenko beam is naturally based on the same kinematic hypothesis as the straight Timoshenko beam, but due to the curvature it couples to axial motion denoted by *w*(*s*,*t*). A formal derivation requires the involvement of differential geometry and the use of Hamilton's principle. In non-dimensional form, the governing equations in the frequency domain are given by
*κ* is the non-dimensional curvature found from the radius of curvature *R** indicated in figure 1.

In the limit of *e*^{iks−iωt} are;

Rayleigh rod

A preliminary justification of whether an overlap exists between the short-wave WKB method applied to long-wave beam models is now given. In Chapman & Sorokin [14], it is demonstrated that the beam model formulation for uniform beams effectively corresponds to neglecting *O*(*h*^{2}/λ^{2}) terms, where *h* is the local cross-wise dimension and λ is the wavelength. Denoting *μ*=*h*/λ and by introducing the relation *μ*=*ϵ*^{q}, it can be shown that for 1/2<*q*<1 both the long-wave beam assumption, *μ*≪1, and short-wave WKB assumption, λ/*l*≪1, are fulfilled when *ϵ*≪1. Meanwhile does it ensure that the largest neglected term in the beam model is smaller than the smallest term accounted for in the WKB method.

## 3. Asymptotic solutions

The WKB solutions to the above-mentioned problems are derived by the formal procedure of asymptotics. It should be mentioned that when solving systems of coupled equations the transport equation is found by imposing orthogonality between the resonant terms of the *O*(*ϵ*) terms and the eigenvector of the *O*(1) terms. Details are given in the following sections.

### (a) The non-uniform Rayleigh beam

The WKB approximation takes its point of departure in a scaling of the axial coordinate
*O*(1) terms yields
*A*(*S*) and *I*(*S*). Collecting the *O*(*ϵ*) terms yields

### (b) The non-uniform Timoshenko beam

The results in this section can either be found by following the procedure demonstrated in §3*c* or alternatively by a deduction of the results presented in §3*c* involving the limit

### (c) The curved non-uniform Timoshenko beam

The derivation of this problem is shown in some detail in the following. Firstly, the axial coordinate *s* is scaled to *S* and the WKB ansatz can be put down as
*O*(1) terms leads to
*k*(*S*), with no odd powers and therefore has a closed form solution, that in the interest of brevity is omitted here. The wavenumbers come in pairs of left and right going waves denoted by ±*k*_{u}(*S*), ±*k*_{β}(*S*) and ±*k*_{w}(*S*), where the indices designate the dominant type of deformation. This can be analysed by inspection of the components of the eigenvector *v*_{j} for a certain wavenumber *k*_{j}(*S*). In the limit of small *κ*, i.e. when the rod is straightened out, the dispersion relation factorizes itself into those for plane axial waves in a corrugated rod and the flexure and shear waves in a corrugated Timoshenko rod
*B*_{0}(*S*) and *W*_{0}(*S*) to *U*_{0}(*S*) determined from (3.12) are
*O*(*ϵ*) in the governing equation after inserting the WKB ansatz yields an inhomogeneous version of (3.12)
*B*_{0,j}(*S*) and *W*_{0,j}(*S*) in (3.19) and to construct an eigenvector *F*_{j}⋅*v*_{j}=0 the transport equation for the present problem emerges. Solving it provides the following amplitude:

### (d) Energy flux

It is remarkable that the transport equation in all of the above cases turns out as
*C*(*S*)*U*_{0}(*S*) must vanish. This has the immediate consequence that an explicit expression for *U*_{0}(*S*) can be determined with no further problems and, furthermore, that the quantity *C*(*S*)*U*_{0}(*S*) is conserved over *S*, i.e. this is an adiabatic invariant. This was recognized by Pierce [9] as it can be proved that conservation of energy flux to the leading order leads to exactly the same requirement for the amplitude function. This can be exemplified by considering the straight Timoshenko beam. The energy flux is
^{☆} denotes complex conjugate. The above expression can be rewritten to
*B*_{0}(*S*) to *U*_{0}(*S*) is

## 4. Determination and comparison of eigenfrequencies

In this section, it is tested how the solutions obtained in the previous section perform numerically when compared to reference solutions. The reference solutions are obtained by the finite element (FE) method as this can be regarded as a reliable and easily applied tool as far as free and forced vibrations of beams of finite length are concerned and therefore is ideal for the comparison. Another advantage is that the FE method is most accurate at low frequencies and therefore is expected to be accurate where the WKB method is expected to be least accurate. For the WKB solutions this requires, however, that these are elaborated to predict standing waves in a beam and thereby match an FE model. Eigenfrequencies, which can be determine via the FE method with no further difficulty, are chosen as the result for the actual comparison. Thus, the WKB method must be specialized to fulfil certain boundary conditions and thereafter the eigenfrequencies can be determined and compared to those of the FE method. Here, clamped–clamped and clamped-free boundary conditions are considered.

A clamped boundary condition is stated as *u*(*s*)=*u*′(*s*)=0 for a Rayleigh beam and as *u*(*s*)=*β*(*s*)=*w*(*s*)=0 for the Timoshenko beam. At a free end, the Rayleigh boundary conditions are *u*′′(*s*)=*u*′′′(*s*)=0 corresponding to zero bending moment and zero shear force, respectively. For a curved Timoshenko beam, the conditions at a free end are
*κ* is simply set to zero in (4.1). It is stressed that taking the derivatives related to the boundary equations produces a series of *O*(1), *O*(*ϵ*), *O*(*ϵ*^{2}),… terms and to be consistent with the WKB solutions found only the *O*(1) terms are retained.

The cross section is defined as circular which immediately suggests a standard choice of shear correction factor
*L** is the length of the beam. These parameters correspond to a length scale of inhomogeneity of *l**=10 *m*. The first six eigenfrequencies can be found in tables 1 and 2 for the Rayleigh beam model, tables 3 and 4 for the Timoshenko beam model, and tables 5 and 6 for the curved Timoshenko beam model. Along with the FE and WKB eigenfrequencies, their mutual deviations determined as *BEAM189* that allows for tapering. This way of modelling ensures that the plane cross-sectional hypothesis is preserved and consequently it is the governing differential equations presented in (2.6), () and () that are approximated numerically. When modelling the Rayleigh beam the shear deformation has simply been suppressed. The results reported in the tables have been generated with 60 elements in the axial direction. Eigenfrequencies of the same models have also been conducted with 20, 30, 40 and 50 elements and their convergence towards the results reported has been observed.

A common way of defining the wave regime in which the WKB approximation is valid is by
*et al.* [7]. This basically requires the wavelength to be short compared with the length scale of inhomogeneity. Remark here the identity *kl*=*k*/*ϵ*. An alternative to the above criteria, equation (4.3), is now proposed based on the fact that the WKB approximation is concerned with the change of the waveguide over one wavelength. Firstly, the change of the waveguide can be stated as the spatial derivative of the diameter
*rϵ* that is the controlling factor in the change of waveguide and not *ϵ* alone. Thus, an alternative version of the failure criteria is likely related to the quantity *k*/*rϵ*. Define
*ξ*∈]0,*L*[ and *x*=*u*, *β*, or *w*. These subscripts again designate the dominant type of deformation for the given wavenumber. It is stressed here that none of the wavenumbers corresponds solely to either flexural, axial or shear motion, but by inspection of the modal coefficients it can be seen that they do correspond dominantly to the indicated type of deformation. The expression in (4.5) is a convenient quantity to track for each eigenmode to assess how much larger than unity *Φ*_{x} is required to be for the WKB method to work well. This criteria can be interpreted as consisting of two conditions: a mild slope condition of the waveguide walls, i.e. *ϵr* must be small, and a short-wave condition, i.e. *k*/*ϵ* must be large, and together these form a criteria that requires only small changes of the waveguide over one wavelength. The values of *Φ*_{x} are given in tables 1–6.

In the following, some comments and observations pertaining to the results are given. From tables 1–4, it is clear, first of all, that the results are in general much better captured for the clamped–clamped model than the clamped-free model of the same length which naturally follows from the fact that the length of the standing wave is shorter in the clamped–clamped model. As a general trend, it is clear that the results improve for large values of *Φ*_{u}. The deviations for the WKB method compared to the FE method are 10%, the value of *ϵ*, when *Φ*_{u} is about 4 for the Rayleigh beam and a threshold *Φ*_{u}<2 confines the errors of the Timoshenko model by as little as 4%.

From tables 5 and 6, it broadly holds that the results agree for large values of *Φ*_{x}, with the exception of a minor increase around the fifth eigenfrequency. In table 6, it is seen that in the transition between the fourth and the fifth mode *Φ*_{w} has changed from being zero to be larger than zero, but smaller than one. Thus, this eigenfrequency is the first one where the axial wave has cut on, and therefore can have a global impact in the motion of the beam, but the value of *Φ*_{w} is not large enough for the asymptotics to work very well and therefore the result is determined with some deviation. As *Φ*_{w} is determined from the minimum value of *k*(*S*), there is no cut-on/cut-off location at this frequency in the beam, but there will be so at a slightly lower frequency.

It is worth considering if an eigenmode that is associated with several different valued *Φ*_{x} is reasonably accurate and how this can be assessed. Thus, the question is whether an eigenfrequency compares well to the exact result even though its associated values of *Φ*_{x} do not collectively indicate accuracy of the result. The following dispute demonstrate the relevance of such a consideration. Mode one in table 6 is found with a large deviation and has one component, *Φ*_{u}>2, which suggests that it is probably accurate, whereas the two others, *Φ*_{β}<2 and *Φ*_{w}<2, suggest that it is not very accurate. Contrary, mode three in table 5 is found with a small deviation, but *Φ*_{w} suggests that it is wrong while *Φ*_{u} suggests that it is correct. The disagreement has a simple explanation that becomes evident by plotting the actual mode shapes. These are seen in figures 2 and 3. Here, it is seen in figure 3 that mode three is dominated by flexural motion, *u*(*S*), and the fact that the value of *Φ*_{w} suggests the mode is outside the WKB validity range is therefore not relevant. From figure 2, it is clear that axial motion cannot be neglected when compared with flexural motion and therefore the approximation fails. In other words, the measure of ‘asymptoticness’, *Φ*_{x}, needs to be accompanied by an assessment of the involvement of modes for multi-modal problems. This involvement of modes can possibly be further quantified by the modal participation factor known from forced vibration [16]. Finally, it is seen that clamping both ends of the curved Timoshenko beam drastically increases the first eigenfrequency when compared with its cantilevered counterpart. This has a simple physical explanation: the bending couples to stretching due to the curvature, but stretching resistance is high since the beam is clamped at both ends.

A comparison of dispersion relationships is shown in appendix A in which the Rayleigh and Timoshenko rod theories are compared to the Pochhammer–Chree elastodynamics solution for waves in a straight elastic rod with circular cross section. The comparison demonstrates that Timoshenko theory for the parameters chosen can be expected to work well until about 8 kHz, whereas the Rayleigh theory breaks down at a much earlier stage. In fact, it can be observed from the results in tables 1 and 3 that already at the sixth eigenfrequency the Rayleigh theory diverges by 17% from the Timoshenko theory.

The attention is now turned to an inspection of how the results compare in the *ϵ*−*r* parameter space—the two parameters free to choose in the modulation of the waveguide with the given functional form of *d*(*s*). In table 7 is the deviation between WKB method and FE method shown for the first eigenfrequency for a clamped–clamped Rayleigh beam at various combination of *ϵ* and *r*. Level curves are drawn in the table that bounds *Φ*_{u} less than *π*, 2, and 1. It can be observed from the results in table 7 that none of the level curves perfectly captures the incorrect results. Thus, it is not possible to state a minimum value of *Φ*_{u}, but only require that *Φ*_{u}≫1. It is observed, however, that in many cases *Φ*_{u}>3 is actually sufficient. Additionally, it is pointed out that the measure for validity proposed in this paper, equation (4.5), correctly predicts that, results based on a large value of *ϵ*, but a small value of *r*, and vice versa, can be determined correctly. This is an improvement compared to the more traditional measure in (4.3).

## 5. Conclusion

In the work presented in this paper, the WKB method is applied to approximate the solution of the wave propagation problem in non-uniform beams. It is demonstrated that the leading order WKB approximation corresponds to conservation of energy flux. The solutions are specialized to clamped–clamped and clamped-free boundary conditions, and a number of eigenfrequencies are determined and compared to those determined by the FE method. A parameter study is conducted and the validity range of the WKB method and the plane cross-sectional hypothesis are discussed. It is found that

It is found that there is an appreciable overlap of validity between the WKB method and Timoshenko beam theory, the overlap with Rayleigh (or Bernoulli–Euler) beam theory is found to be much more narrow.

The simplicity of the WKB approximation, and the fact that the dispersion relation matches that of the unperturbed waveguide, illustrates the power of the method in the analysis of non-uniform rods. This is fortunate as the wavenumbers often will be known *a priori*; these are just the same as the wavenumbers found for the corresponding uniform waveguide. It seems natural to extend the analysis to waveguides with constant cross section, but varying curvature. These kind of structures are of practical interest in the analysis of the linear time harmonic dynamics of suspension systems, pipes, curved rails, etc., where one would be concerned with noise transmission and energy flux. Since the ‘standard’ WKB approximations reported in this paper, as is shown, conserves the energy flux it seems obvious to propose to analyse turning point problems. Doing so allows to study problems where the net energy flux is zero because of the wave reflection at a turning point. This will give detailed information on the internal reflection of waves in a non-uniform beam structure and can potentially be used to analyse asymptotically stop and pass bands in rods with periodic curvature.

## Funding statement

The research is funded by *The Danish Council for Independent Research, Technology and Production Sciences*.

## Acknowledgements

The authors address their sincere gratitude to Prof. C. J. Chapman for the many discussions on the topic.

## Appendix A. Comparison of rod theories and elastodynamics

Here, the validity ranges of the Rayleigh beam theory and the Timoshenko beam theory are demonstrated and discussed briefly. In figure 4, the dispersion diagram for a beam in terms of the Rayleigh beam model, the Timoshenko beam model and the Pochhammer–Chree exact solution of elastodynamics is shown. For details on the Pochhammer–Chree solution, see Chree [17]. The plot illustrates how well the first branch of the Timoshenko model compares to the Pochhammer–Chree solution in a large frequency range (until *ω**≈14 *kHz* corresponding to *ω*≈3), whereas the Rayleigh model compares well only at much lower frequencies. There is an obvious mismatch between Timoshenko theory and Pochhammer–Chree when it comes to the second branch (the shear wave). This is not of any great importance for frequencies lower than this cut-on frequency as the shear wave in this regions is evanescent. It should be mentioned that the cut-on can be perfectly matched with a proper selection of λ. For the non-uniform Timoshenko beams, the lowest cut-on is found to be *ω**≈8 *kHz* for both the straight and curved beam model with the parameters in §4. The disagreement between Rayleigh and Timoshenko model can even be observed by comparing the results in, for instance, tables 1 and 3. Here, the deviation between the Rayleigh results and the Timoshenko results is about 17% already at the sixth eigenfrequency.

- Received October 29, 2013.
- Accepted March 25, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.