- Email: [email protected]
A wide class of analytical solutions of the Webster equation
Journal Pre-proof A wide class of analytical solutions of the Webster equation M. Bednarik, M. Cervenka
PII:
S0022-460X(19)30732-1
DOI:
https://doi.org/10.1016/j.jsv.2019.115169
Reference:
YJSVI 115169
To appear in:
Journal of Sound and Vibration
Received Date: 15 April 2019 Revised Date:
20 December 2019
Accepted Date: 25 December 2019
Please cite this article as: M. Bednarik, M. Cervenka, A wide class of analytical solutions of the Webster equation, Journal of Sound and Vibration (2020), doi: https://doi.org/10.1016/j.jsv.2019.115169. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Credit Author Statement
Michal Bednarik conceived the presented analytical method, wrote the text of the manuscript and performed all the analytical calculations. Milan Cervenka verified all presented results based on numerical calculations, drew the figures in the manuscript and helped with some statements in the text of the manuscript. Both authors discussed the results and commented on the manuscript.
A wide class of analytical solutions of the Webster equation M. Bednarika,∗, M. Cervenkaa a Czech
Technical University in Prague, Faculty of Electrical Engineering, Technicka 2, 166 27 Prague 6, Czech Republic
Abstract This paper aims at presenting closed-form general analytical solutions of the Webster equation describing plane elastic or acoustic waves. The considered radius functions of non-uniform cross-sectioned rods or ducts are based on the triconfluent Heun functions and contain some optional parameters enabling us to set various profiles of the radius functions in a relatively wide range, while it is possible to employ the presented exact general analytical solution of the Webster equation for all selected profiles. If the radius functions are predetermined, then the derived general analytical solution can also be employed for their triconfluent Heun approximations, including certain polynomial ones. The applicability and correctness of the derived analytical solutions are demonstrated by calculations of natural frequencies and mode shapes for representative radius functions while the results based on approximate analytical solutions are verified numerically. Keywords: Triconfluent Heun equation, Triconfluent Heun function, Webster equation, non-uniform cross-sectioned waveguide, horn
1. Introduction There are different technological devices for the transmission of acoustic energy into an acoustic load which are based on waveguiding systems. We can meet them, for example, in heat and mass transfer processes, plastic welding, ∗ Corresponding
author Email address: [email protected] (M. Bednarik)
Preprint submitted to Journal of LATEX Templates
December 26, 2019
metal treatment, ultrasonic sonochemistry, acoustic transmission lines, phononic and sonic crystals etc. Under simplifying assumptions, we can describe the propagation of sound (or elastic waves) in waveguiding systems with variable cross-sections as pipes, horns, concentrators, ultrasonic rod waveguide-radiators, by the Webster equation (see e.g. [1, 2, 3]). The solution of the Webster equation enables us to optimally design the waveguiding systems with respect to the characteristics of our interest. The performance of many ultrasonic devices depends on a proper design of a sonotrode shape (concentrator). For this reason, it is important to known the dependence of fundamental modal properties (natural frequencies, mode shapes) of sonotrode shapes on their geometrical parameters (see e.g. [4]). Also in the design of ultrasonic rod waveguide-radiators, it is desirable to know the solution of the Webster equation (see e.g. [5]). Locally periodic structures (phononic or sonic crystals) composed of periodically repeated elements are used in many fields of science and technology (vibration isolators and mechanical filters to building frames, bridge trusses, railway tracks, compound pipes etc.) (see e.g. [6]). Therefore, the determination of widths and locations of band-gaps and the corresponding attenuation levels associated with these systems represents an important practical problem. Furthermore, we may encounter the need to know the solution of Webster equation in musical acoustics in exploring and designing wind instruments, where the main function of the horn (bell) is to form a more gradual impedance transition between the high inner tube impedance and the low outside air impedance (see e.g. [3]). The Webster equation has its use in the design of acoustical horns applied to loudspeakers (see e.g. [2]) and it is helpful in designing very small (even miniaturized) components (e.g. small transducers, earplugs) (see e.g. [7]). The solution of the Webster equation also finds application in acoustic filters and acoustic resonators (see e.g. [8, 9, 10, 11]). However, to the best of the authors’ knowledge, only comparatively few publications are available that present the exact analytical solutions of the Webster equation for various cross-section profiles of waveguides (see e.g. [12, 13, 14, 15, 2
16, 17]). Therefore, the primary motivation of this work is to present both new exact analytical solutions of the Webster equation and a method that enable us to find approximate analytical solutions of this equation for the cases when the exact ones are not known. To accomplish this intention, we used cross-section profiles based on the triconfluent Heun functions whose evaluation is already a standard part of mathematical software such as Maple. The functions have attracted the attention of scientists mainly during the last decade. Nowadays, we can observe a growing interest in their utilization in various fields of physics even if they are still waiting for their broader application particularly in acoustics and elastodynamics. The advantage of these functions is that they contain a sufficient number of parameters for modification of the cross-section profiles, without the need to employ many analytical solutions of the Webster equation. As shown below, it is possible in this way to vary shapes and sizes of the cross-section profiles in a comparatively wide range. For example, using these parameters it is possible to set the cross-section profiles of rods that are changing smoothly and continuously (e.g. with zero derivatives at the endpoints) along a desired dimension to avoid failures which are likely to initiate at the places where changes in material properties are not gradual. This paper is organized as follows. The Webster equation for elastic and acoustic waves is presented in Section 2. Based on mapping the Schrödinger-like differential equation into the Webster equation, its general solution is derived in Section 3. The radius-function profiles expressed using the triconfluent Heun functions for which exact analytical solutions of the Webster equation exist are given in Section 4. The possibilities of expressing the radius-functions based on polynomial solutions of the triconfluent Heun differential equation are also presented in this section. Section 5 is devoted to approximate analytical solutions of the Webster equation for predetermined radius-functions. The applicability and correctness of the derived analytical solutions are demonstrated by calculations of natural frequencies and mode shapes for representative radius-functions in Section 6. Finally, Section 7 provides our conclusions. In order to make this paper sufficiently self-contained and more clear, some features concerning the 3
triconfluent Heun differential equation and its solutions including polynomial ones are outlined in Appendix.
2. Webster equation We consider plane elastic wave propagation through an isotropic non-uniform cross-sectioned slender rod or tube, see Fig. 1 where the geometry is specified. Assuming small vibration amplitudes and small transverse dimensions then longitudinal propagation of elastic waves can be described by means of the following model equation (in the absence of body force) (see e.g. [5]) ∂ ∂u(z, t) ∂ 2 u(z, t) , EA(z) = ρs A(z) ∂z ∂z ∂t2
(1)
where u(z, t) is the displacement of the rod along the length direction from its equilibrium position at point z and time t. Further, E, A(z), and ρs are Young’s modulus, cross-sectional area function and density of the rod, respectively. If we assume linear plane wave acoustic propagation along a narrow fluid-filled tube of gradually varying cross-section then the corresponding model equation can be expressed as (see e.g. [18]) ∂ ∂p(z, t) ∂ 2 p(z, t) BA(z) = ρA(z) , ∂z ∂z ∂t2
(2)
where p(z, t), ρ and B are the acoustic pressure, fluid density and the adiabatic bulk modulus, respectively. Limiting ourselves to the harmonic solutions of Eqs. (1) and (2), i.e. u(z, t) = √ U (z) exp(−jωt) and p(z, t) = P (z) exp(−jωt), where j = −1 is the imaginary unit, and ω is the angular frequency, then the model equations (1) and (2) reduce to the Webster equation: d2 Φ(z) 1 dA(z) dΦ(z) + + k 2 Φ(z) = 0 . 2 dz A(z) dz dz
(3)
Here Φ(z) stands for either U (z) or P (z), k = ω/c is the wave number and c is p the wave speed. The wave speed for elastic waves is c = ρs /E and for acoustic p waves c = ρ/B. 4
r(z)
A(z) = πr2(z)
0
H
z
Figure 1: Rod/tube with variable cross-section.
Using the Webster equation for plane elastic and acoustic waves has its limitations that are discussed in, e.g. [5, 19, 20, 21]. Employing the Webster equation enables us to treat propagation of elastic or acoustic waves as an one-dimensional problem if the characteristic width of a rod or tube is small compared to the wavelength and the cross section A(z) varies relatively slowly. Considering an axially symmetrical waveguide, i.e. A(z) = πr2 (z) (r(z) is a radius-function) then, for convenience, we can rewrite Eq. (3) into the following form:
d2 Φ(s) 2 dR(s) dΦ(s) + + K 2 Φ(s) = 0 , ds2 R(s) ds ds
(4)
where R(s) = r(s)/r0 (dimensionless radius-function), s = z/`, K = k`, ` is a characteristic length, and r0 = r(0).
3. General analytical solution of the Webster equation We can solve the Webster equation (4) based on a time-independent onedimensional Schrödinger-like differential equation d2 ψ(s) 2 + K + G(s) ψ(s) = 0 , 2 ds 5
(5)
that can be mapped into Eq. (4) by the following transformation (6)
ψ(s) = R(s)Φ(s) . Substituting this transformation relation into Eq. (4) we obtain d2 Φ(s) 2 dR(s) dΦ(s) Φ(s) d2 R(s) 2 + + K + G(s) Φ(s) = 0 . + 2 2 ds R(s) ds ds R(s) ds
(7)
By identifying this equation with Eq. (4) we obtain the following condition on the transformation
d2 R(s) + G(s)R(s) = 0 . (8) ds2 Thus, searching for the general solution of the Webster equation (4) is reduced to finding a general solution of the Schrödinger-like equation (5). 3.1. General solution of the Schrödinger-like equation We can see from Eq. (8) that the (generating) function G(s) determines the radius function R(s). If we were able to solve Eq. (5) for a given function G(s), then using the transformation relation (6) we could solve the Webster equation (4) for the cross-sectional area function A(s) = πr02 R(s)2 . It is obvious that for G = a ≥ 0 one can very easily solve Eqs. (5) and (8). The corresponding radius-functions are presented, e.g. in the works [9, 12, 5]. Employing more complex generating functions other exact analytical solutions can be found, e.g. in the works [14, 16, 17]. In this work, we assume that the generating function G(s) is a quartic polynomial, which enables us to solve Eq. (8) exactly and in addition, such a generating function contains a sufficient number of coefficients that allow us to determine an appropriate radius profile. The generating function can be expressed as G(s) = a0 + a1 (s − s0 ) + 2a2 (s − s0 )2 + a3 (s − s0 )3 − a4 (s − s0 )4 ,
a4 6= 0 . (9)
Substituting this generating function into Eq. (5) we obtain d2 ψ(s) 2 + K + a0 + a1 (s − s0 ) + 2a2 (s − s0 )2 + ds2 a3 (s − s0 )3 − a4 (s − s0 )4 ψ(s) = 0 . (10) 6
In order to find the solution of Eq. (10) we first employ the following transformation of the dependent variable ψ(s) = exp
Z 1 F (s)ds ϕ(s) 2
(11)
for Eq. (10) and we attain d2 ϕ dϕ F (s)2 1 dF (s) + F (s) + + + K 2 + a0 + a1 (s − s0 )+ 2 ds ds 4 2 ds 2a2 (s − s0 )2 + a3 (s − s0 )3 − a4 (s − s0 )4 ϕ(s) = 0 . (12) Substituting 2 √ a3 16a2 a4 + 3a23 + F (s) = −2 a4 (s − s0 ) − 3 4a4 8a42
(13)
into Eq. (12) we obtain the form of the triconfluent Heun differential equation ( ) 2 √ d2 ϕ a3 16a2 a4 + 3a23 dϕ − 2 a4 (s − s0 ) + + − 3 ds2 4a4 ds 8a42 " 5 8a1 a24 + 8a2 a3 a4 + a33 − 16a42 (s − s0 )+ 8a24 # 5 64(K 2 + a0 )a34 + 64a22 a24 + 16a2 a23 a4 + 32a3 a42 + a43 ϕ(s) = 0 . (14) 64a34 Using the following transformation for an independent variable: a3 σ = Q (s − s0 ) − , 4a4 where
√ 13 2 a4 Q= , 3
(15)
(16)
we can rewrite Eq. (14) into its canonical form: ! 2 d2 ϕ(σ) 16a a + 3a dϕ(σ) 2 4 3 + − 3σ 2 − 3 dσ 2 dσ 2 8Qa4 64(K 2 + a0 )a34 + 64a22 a24 + 32a2 a33 a4 + 16a1 a3 a24 + 3a43 + 64Q2 a34 # 5 8a1 a24 + 8a2 a3 a4 + a33 − 16a42 σ ϕ(σ) = 0 . (17) 8Q3 a24 7
By comparing Eq. (17) with the canonical triconfluent Heun equation (A.3) in the Appendix, we obtain α=
64(K 2 + a0 )a34 + 64a22 a24 + 32a2 a23 a4 + 16a1 a3 a24 + 3a43 , 64Q2 a34 24a4 (a1 a4 + a2 a3 ) + 3a33 16a2 a4 + 3a23 β= , γ=− . (18) 5 3 16a42 8Qa42
Employing Eq. (A.16) we can write the general solution of Eq. (17) as ϕ(σ) = A1 THF (α, β, γ; σ) + A2 exp σ 3 + γσ THF (α, −β, γ; −σ) , (19) where A1 and A2 are integration constants and THF stands for the triconfluent Heun function. Using the relation (A.18) and (15) we can write the general solution (integration constants are arbitrary) of Eq. (10) in the following form:
3 Q s− ψ(s) = A1 exp −
a3 4a4
− s0
3
+ γQ s − 2
a3 4a4
− s0 ×
a3 THF α, β, γ; Q s − − s0 + 4a4 3 a3 a3 3 Q s − 4a4 − s0 + γQ s − 4a4 − s0 A2 exp × 2 a3 − s0 THF α, −β, γ; −Q s − . (20) 4a4 3.2. General solution of the Webster equation for the specific class of radius functions Employing the transformation relation (6) and the solution (20) it is possible to write the general closed-form solution of the Webster equation (4) for the
8
specific class of the radius functions that is given in the Sec. 4, as 3 3 Q a3 a3 exp − Q2 s − 4a − s − γ s − − s 0 0 2 4a4 4 Φ(s) = A1 × R(s) a3 THF α, β, γ; Q s − − s0 + 4a4 3 Q3 Q a3 a3 exp 2 s − 4a4 − s0 + γ 2 s − 4a4 − s0 × A2 R(s) a3 THF α, −β, γ; −Q s − − s0 4a4 = A1 Φ1 (s) + A2 Φ2 (s) . (21) Based on what has been shown above, there are two approaches to solving the Webster equation. The first approach consists in solving Eq. (8) for a predetermined generating function G(s), which leads to obtaining the radiusfunction R(s). It is obvious that if we know the solution of Eq. (5) at the same time, then analytical solution (21) is exact. The second approach, which is often convenient from practical point of view, is based on a predetermined radius-function R(s) when the function G(s) is calculated by using Eq. (8). However, to make the analytical solution (21) exact, it is necessary to solve Eq. (5) exactly as well. Unfortunately, there are only a few radius-functions for which we can solve the Webster equation exactly. To avoid this fact we can employ a quartic polynomial that approximates the function G(s). In this case, the Webster equation solution (21) is no longer exact, though it is an analytical one. 4. Radius functions for the exact solution of the Webster equation One of the intended aims of this work is to extend the group of radius profiles for which exact solutions of the Webster equation are known. Based on the results from Sec. 3 we can accomplish this task by employing the generating function G(s) in the form the quartic polynomial (9). However, from a practical point of view, it is better to reduce the number of degrees of freedom in the 9
search for the radius function by omitting some of lower terms of the quartic polynomial (a4 > 0). For this reason, we will no longer consider the cubic term in this polynomial to reduce the number of optional parameters. With respect to the identity (A.17)) we assume that a2 6= 0. Substituting the generating function (9) with a3 = 0 into Eq. (8) we obtain an equation that is of the same form as Eq. (10). This fact enables us to directly utilize the solution (20) where we substitute for K 2 and a3 the zero value. After small algebraic modifications we can write the solution of Eq. (8) as 2 3 s a2 R(s) = C1 exp − Q3 s − ss0 + s20 − × 2 3 a4 2 a2 + a0 a4 3a1 3Q2 a2 THF , , − ; Q(s − s ) + √ 0 a4 Q2 2 a4 a4 2 s a2 3 − ss0 + s20 − × C2 exp Q3 s 2 3 a4 2 a2 + a0 a4 3a1 3Q2 a2 THF , − , − ; −Q(s − s ) = √ 0 a4 Q2 2 a4 a4 C1 R1 (a0 , a1 , a2 , a4 , s0 ; s) + C2 R2 (a0 , a1 , a2 , a4 , s0 ; s) . (22) We can see that the derived radius-function (22) can be adjusted based on 7 optional parameters, i.e. two integration constants C1,2 and a group of five internal parameters (a0 , a1 , a2 , a4 , s0 ). This fact enables us to find an entire family of radius-functions of various profiles for which we have one exact solution of the Webster equation (4) by substituting the relation (22) for R(s) in the solution (21). If we want to determine the radius at the point s = 0 in order to correspond to the radius r0 then we choose values of the integration constants C1,2 in such a way that the radius-function R(s = 0) = 1. In this case, required values of the radius-function R(s) at another point si ∈ (0, sH ], where sH = H/` is a dimensionless length of the rod or tube, are given only by the choice of the optional parameters. To do that, we can change one of the optional parameters, and we fix the remaining ones until the radius function reaches the required value at the given point si . It is worth noting that if the argument of the triconfluent Heun functions is 10
equal to zero then these functions equal unity for each of the parameters α, β, γ. Additionally, the first derivative of the Heun functions are always equal to zero if their arguments are zero. Examples of radius-function profiles for various optional parameters are shown in Fig. (2)(a). However, if we require to determine a radius r0 at the point s = 0 and at
(a)
(b)
1.25
1.25 R(s)
1.5
R(s)
1.5 1
1
0.75
0.75
0.5
0.5 0
0.5
1 s
1.5
2
0
0.5
1 s
1.5
2
Figure 2: (a) Examples of radius-function profiles R(s) based on the formula (22) for a given sequence of optional parameters (C1 , C2 , a0 , a1 , a2 , a4 , s0 ): (1,0,-0.1,0,0.01,0.08,0) – black line, (0,1,-1.99,0,-0.04,0.001,0) – blue line, (1,0,-3,0,-0.05,0.001,0) – red line, (1,0,0,0,0.01,0.001,0) – brown line; (b) Examples of radius-function profiles R(s) based on the formula (23) for a given sequence of optional parameters (a0 , a1 , a2 , a4 , s0 , si , ri /r0 ): (-0.7,2,0.2,0.66,0,2,3) – black line, (0,0,2.32856,0.85,1,1,3) – red line, (0,0,-0.01,0.045,1,2,0.5) – blue line, (0,0,-1.5,0.55,1,2,0.5) – brown line.
the same time a radius ri at the point si ∈ (0, sH ] with fixed optional internal parameters then we need to employ both the constants C1,2 that depend on a value of both the radii r0 and ri . In this case, we can write the radius-function (22) as R(s) =
ri r0 R2 (0)
− R2 (si )
W
R1 (s) −
ri r0 R1 (0)
− R1 (si )
W
R2 (s) ,
(23)
where W = R2 (0)R1 (si ) − R2 (si )R1 (0) . Here for brevity, we do not write the parameters in functions R1,2 (s). 11
(24)
4.1. Some polynomial solutions for the radius function Due to the fact that it is relatively difficult to estimate the profiles of the triconfluent Heun functions based on their parameters α, β, γ, it seems to be useful to express them as terminating power series, i.e., polynomials, if certain conditions are satisfied, see Appendix A.2. Based on the theorem presented in Appendix A.2 where the meaning of the used symbols below is also given, it is possible to write the first five polynomial solutions of the triconfluent equation: 1. For N = 0 (β = 3), we can directly find that D1 ≡ α = 0 and P0 = p0 and the triconfluent function equals to unit, i.e. THF(0, 3, γ; σ) = 1 ;
(25)
2. For N = 1 (β = 6), D2 ≡ α2 + 3γ = 0, P1 (σ) = σ −
α ; 3
3. For N = 2 (β = 9), D3 ≡ α3 + 12αγ + 36 = 0, P2 (σ) = σ 2 −
α α2 1 σ+ − ; 3 36 α
4. For N = 3 (β = 15), D4 ≡ α5 + 60α3 γ + 756α2 + 576αγ 2 + 5184γ = 0, γ α2 7 2 α α3 + − αγ − ; P3 (σ) = σ 3 − σ 2 + σ− 3 2 18 162 54 3 5. For N = 4 (β = 12), D5 ≡ α4 + 30α2 γ + 216α + 81γ 2 = 0, α P4 (σ) = σ − σ 3 + 3 4
where α=
a22 + a0 a4 , a4 Q2
3 α2 2 α 5 4 2 + γ σ − + αγ + σ+ 18 3 162 27 3 α4 2 5 γ2 + α2 γ + α + , 1944 81 18 9
γ=−
3Q2 a2 a4
and σ = Q(s − s0 ) .
(26)
It should be noted that the polynomial solutions for N ≥ 1 do not represent triconfluent Heun functions according to the commonly used definition because they do not satisfy the condition b1 = 0, as it is supposed in (A.11). 12
For N > 4, we can continue finding the couples (α, γ) that enable us to express the polynomial solutions but higher order polynomials lose their usefulness in the sense of estimation of their profiles and therefore we do not consider them here. Using the first two polynomial solutions above, we obtain the following expressions for the radius function R(s): 1. N = 0;
3 R(s) = C1 exp − Q3 s 2
s2 a2 − s0 s + s20 − 3 a4
(27)
.
2. N = 1; 2 s a2 a22 + a0 a4 3 3 2 − s0 s + s0 − . Q(s − s0 ) − R(s) = C1 exp − Q s 2 3 a4 3a4 Q2 (28)
(a)
(b)
1.5
1.75 1.5 1.25
1
R(s)
R(s)
1.25
1
0.75
0.75
0.5
0.5 0
0.5
1 s
1.5
2
0
0.5
1 s
1.5
2
Figure 3: (a) Examples of radius-function profiles R(s) based on the formula (27) for a given sequence of optional parameters (C1 , a2 , a4 , s0 ): (1,0.161,0.12075,0) – black line, (1,0.005,0.0011,0) – blue line, (1,0.005,0.008,1) – red line, (1,1,6,1) – brown line; (b) Examples of radius-function profiles R(s) based on the formula (28) for a given sequence of optional parameters (C1 , a0 , a2 , a4 , s0 ): (0.5510,-2,-0.1,0.0947,0) – black line, (-0.3712,3,0.1,0.3,1) – red line, (0.5166,-3,0.1,0.3,0) – blue line, (0.9927,-2,0,0.15,1) – brown line.
The polynomial solutions of Eq. (8) for N ≥ 1 may be expressed as linear combinations of the two standard solutions (22). 13
It is necessary to note that it is not possible to express the independent function R2 (s) in the polynomial form for the same optional parameters. Thus, for the second independent solution, it is necessary to use the corresponding triconfluent Heun function or to employ the relation (A.22) if we are able to integrate it. For this reason, it is more convenient to employ the polynomial solutions mainly for the method based on setting suitable optional parameters. The material functions based on the first two polynomial solutions enable us to calculate the corresponding optional parameters analytically, some results being shown in Fig. 3. For the remaining polynomials, it is easier to seek the values of the respective parameters numerically.
5. Approximate analytical solutions of the Webster equation If the radius function R(s) is predetermined then the function G(s) can be calculated based on Eq. (8), i.e. G(s) = −
1 d2 R(s) . R(s) ds2
(29)
However, we generally cannot solve Eq. (5) for the function G(s) given by Eq. (29). There are predetermined radius functions for which we can approximate this function by a quartic polynomial with sufficient accuracy, i.e. G(s) ≡ −
1 d2 R(s) u W (s) = R(s) ds2 a0 + a1 (s − s0 ) + 2a2 (s − s0 )2 + a3 (s − s0 )2 − a4 (s − s0 )4 .
(30)
The coefficients of the quartic polynomial W (s) can be calculated for instance based on an interpolation or by a least-squares approximation. Replacing the function G(s) by the polynomial W (s) in Eq. (5) we can employ the function (21) as the analytical solution of Eq. (4). In other words, we substitute coefficients a0 , a1 , a2 , a3 and a4 from the quartic polynomial (30) into the solution (21).
14
6. Natural frequencies and mode shapes for non-uniform rods Without loss of generality, we limit ourselves only to elastic waves (Φ(s) ≡ U (s)) in this section, where the presented procedures are also directly applicable to acoustic waves in fluids. Based on the following representative examples, we demonstrate the applicability and correctness of the derived solutions. 6.1. Representative examples based on the exact solutions We can calculate natural frequencies based on the solution (21) for the following three classical boundary conditions: 1. Clamped-clamped rod (CCR): U (0) = A1 U1 (0) + A2 U2 (0) = 0 , U (sH ) = A1 U1 (sH ) + A2 U2 (sH ) = 0 ; (31) 2. Clamped-free rod (CFR): U (0) = A1 U1 (0) + A2 U2 (0) = 0 , dU (s) dU1 (s) dU2 (s) = A + A =0; 1 2 ds s=sH ds s=sH ds s=sH
(32)
3. Free-free rod (FFR): dU2 (s) dU1 (s) dU (s) = A1 + A2 =0, ds s=0 ds s=0 ds s=sH dU2 (s) dU (s) dU1 (s) + A2 = A1 ds ds ds s=sH
s=0
= 0 . (33) s=sH
For all the reported radius-functions below, which are found based on the predetermined generating functions (9), we assume the internal parameters a3 = s0 = 0 and the endpoint sH = 1. In addition, we select the radiusfunctions with the zero derivatives at the endpoints and R(0) = 1. To meet these requirements for the endpoint s = 0 we choose the integration constants C1 = C2 = 1/2 in the relation (22). The zero derivative of the radius function at the endpoint s = sH is accomplished by numerical calculation of the internal
15
parameter a1 by fixing the predetermined parameters a0 , a2 6= 0, and a4 > 0. Thus, the radius-function reads 2 1 3 a2 s R(s) = exp − Q3 s − THF (α, β, γ; Qs) + 2 2 3 a4 2 1 3 3 a2 s exp Q s − THF (α, −β, γ; −Qs) , (34) 2 2 3 a4 dR(s) dR(s) =0, = R(0) = 1 , ds s=0 ds s=sH where Q=
√ 13 2 a4 a2 + a0 a4 3a1 3Q2 a2 , α= 2 , β= √ , γ=− . 2 3 a4 Q 2 a4 a4
(35)
(36)
Based on the solution (21) and the radius function (34) we can, after some algebra, express the general solution for the radius-function (34) as U (s) = A1 U1 (s) + A2 U2 (s) = 2 2THF α + K , β, γ; Qs 2 Q h i A1 + 2 3 THF (α, β, γ; Qs) + exp 3Q s s3 − aa42 THF (α, −β, γ; −Qs) 2 2THF α + K Q2 , −β, γ; −Qs h i A2 . (37) 2 exp −3Q3 s s3 − aa24 THF (α, β, γ; Qs) + THF (α, −β, γ; −Qs) Non-trivial solutions of pairs of homogeneous equations (31)-(33) result in the vanishing of determinants. For this purpose, we can use the following relations that are obtained based on the general solution (37): dU1 (s) dU2 (s) a2 a2 U1 (0) = 1 , U2 (0) = 1 , = √ , = − √ . (38) ds s=0 a4 ds s=0 a4 To meet the condition that the determinants of the corresponding pairs of homogeneous equations are zero, then taking into account the relations (38) we have to solve the following transcendental equations for eigenvalues K: 1. Clamped-clamped rod (CCR): 1 a2 K2 exp 3Q3 − THF α + 2 , −β, γ; −Q − 3 a4 Q K2 THF α + 2 , β, γ; Q = 0 , (39) Q 16
2. Clamped-free rod (CFR): 1 a2 K2 exp 3Q3 − THF α + 2 , −β, γ; −Q − 3 a4 Q 2 K K2 THF0 α + 2 , −β, γ; −Q − THF α + 2 , β, γ; Q − Q Q K2 THF0 α + 2 , β, γ; Q = 0 , (40) Q 3. Free-free rod (FFR): 1 a2 K2 exp 3Q3 − THF α + 2 , −β, γ; −Q − 3 a4 Q 2 K K2 THF0 α + 2 , −β, γ; −Q + THF α + 2 , β, γ; Q + Q Q K2 THF0 α + 2 , β, γ; Q = 0 , (41) Q where the prime stands for the derivative with respect to the variable s, see (A.14). To demonstrate the applicability of the presented results we choose four representative radius-functions that meet the above-mentioned requirements. The chosen radius-functions are shown in Fig. (4) and are specified by the optional parameters that are arranged in parameter vectors v = (C1 , C2 , a0 , a1 , a2 , a3 , a4 , s0 ). First (fundamental) eigenvalues are calculated from the transcendental equations (39)-(41) for the four considered rod radius-functions, see Tab. 1 (the mathematical software Maple (version 17) was used for the evaluations of the triconfluent Heun functions). The eigenvalues K are related to natural frequencies via the relation ω = Kc/`. Three normalized mode shapes for the radius function (see Fig. 5(a)) that has the zero derivatives at the endpoints and takes the same value at these points, i.e. R(0) = R(1), are shown in Fig. 5(b). 6.2. Representative examples based on the approximate solutions As representative examples, we assume the following two predetermined radius functions: R(s) = exp 17
s4 2
,
(42)
1.2
R(s)
1
0.8
0.6
0.4 0
0.2
0.4
0.6
0.8
1
s Figure
4:
R(s)
Radius-functions
(C1 , C2 , a0 , a1 , a2 , a3 , a4 , s0 ):
for
various
vectors
of
parameters
(0.5, 0.5, −7, 25.74836, 10, 0, 70, 0)
=
v1
-
v
black
= line;
v2 = (0.5, 0.5, −7, 34.68600, 2, 0, 70, 0) - blue line; v3 = (0.5, 0.5, 0, 9.16938, 0.1, 0, 30, 0) - red line; v4 = (0.5, 0.5, 0, 6.62284, 0.1, 0, 20, 0) - brown line.
CCR
CFR
FFR
v1
2.4832
1.7573
3.9262
v2
2.4766
1.9611
4.0974
v3
2.9169
2.2463
3.7469
v4
2.9539
2.0434
3.5166
Table 1: First (fundamental) eigenvalues for four rod radius functions represented by corresponding parameter vectors v for the boundary conditions (39)-(41)
R(s) =
1 , 1 + s2 − 0.5s4
(43)
for which we can find approximate analytical solutions employing the solution (21) as it is described in Sec. 5. Considering the radius function (42) we can interpolate the corresponding func-
18
(b)
(a) 1 0.5 U (s)
R(s)
1.3 1.2
0
1.1
-0.5
1
-1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
s Figure 5:
0.6
0.8
1
s
(a) Radius-function R(s) for CFR with the vector of parameters v
=
(0.5, 0.5, −7, 25.1795, 2, 0, 36.9072, 0); (b) The corresponding three mode shapes for the first three eigenvalues: K = 1.4917 - black line; K = 4.8762 - blue line; K = 7.9167 - red line.
tion G(s) from Eq. (29) by the following quartic polynomial G(s) = −
1 d2 R(s) u W (s) = −16.3856s4 +18.6067s3 −13.00832s2 +0.7872s . R(s) ds2 (44)
The radius-function is depicted in Fig. 6(a). Substituting the coefficients of the polynomial (44) into the solution (21) we obtain the approximate analytical solution for longitudinal elastic waves in the slender non-uniform cross-sectioned rod. Employing this approximate solution then we can calculate e.g., the corresponding mode shapes of CCR for the first three eigenvalues that are shown and compared with numerical results in Fig. 6(b) (all numerical calculations were conducted by means of numerical software COMSOL Multiphysics). As can be seen from the comparison in Fig. 6(b), the differences between analytical and numerical results are indistinguishable. The radius function (43) is shown in Fig. 7(a). The corresponding polynomial G(s) can be approximated as G(s) = −
1 d2 R(s) u W (s) = −16.1372s4 +34.5048s3 −23.4950s2 +0.4608s+2 . R(s) ds2 (45)
19
(a)
(b) 1
1.6 U (s)
R(s)
0.5 1.4
0
1.2
-0.5
1
-1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
s
0.6
0.8
1
s
Figure 6: (a) Radius functions R(s) (42) for CCR; (b) The corresponding three normalized mode shapes for the first three eigenvalues: K = 3.4237 - black line; K = 6.4686 - blue line; K = 9.5555 - red line. The corresponding numerical results are represented by dots.
After substituting the coefficients of the quartic polynomial into the solution (21) we obtain the approximate analytical solution that can be used for calculation of eigenvalues and mode shapes as in the previous example, Fig. see 7(b). By comparison with numerical results, we can see that the differences are again indistinguishable.
(b)
1
1
0.9
0.5 U (s)
R(s)
(a)
0.8
0
-0.5 0.7
-1 0
0.2
0.4
0.6
0.8
1
0
s
0.2
0.4
0.6
0.8
1
s
Figure 7: (a) Radius functions R(s) (42) for CCR; (b) The corresponding three normalized mode shapes for the first three eigenvalues: K = 3.1763 - black line; K = 6.2972 - blue line; K = 9.4347 - red line. The corresponding numerical results are represented by dots.
20
7. Conclusion In this work, we have presented a closed-form analytical solution of the Webster equation for both elastic and acoustic waves propagating through waveguides with variable cross-sections. If a radius-function is given by the triconfluent Heun functions, then the presented analytical solution is exact. The use of triconfluent Heun functions allows to set the profile by seven optional parameters that affect the radius function profile, both in shape and size over a relatively wide range. We can use the presented analytical solution for all radius-function profiles designed in this way. To facilitate the search for suitable parameters for determining the profile shapes of radius-functions, we have included polynomial forms of the triconfluent Heun functions that can be found under certain conditions specified in the Appendix. In the case of a predetermined radius-function, then after its approximation by a polynomial solution of the triconfluent Heun equation, we can also employ the derived analytical solution of the Webster equation but only as an approximate one. We have demonstrated the applicability and correctness of the derived analytical solutions in calculations of natural frequencies and mode shapes for representative radius-functions. The approximate analytical solutions have been compared with numerical results to show that acceptable accuracy is achieved. The presented analytical solution of the Webster equation extends the group of already known exact solutions, and this solution is applicable to a wide range of radius-functions. In addition, the derived solution can be relatively straightforwardly extended even for the case of longitudinal viscoelastic waves or for the case of thermo-viscous fluids. Acknowledgement This work was supported by the the Grant Agency of the Czech Republic (GACR) grant No. 18-24954S.
21
Appendix A. The triconfluent Heun equation and its solution Here, we outline some facts concerning the triconfluent Heun equation and its solutions. The Heun equation is known to be a standard form of the second-order Fuchsian differential equation with four regular singularities {0, 1, a, ∞}. The canonical form of the general Heun equation is (see, e.g., [22, 23, 24]): γ d2 ϕ δ dϕ αβσ − f + + + + ϕ(σ) = 0 , dσ 2 σ σ − 1 σ − a dσ σ(σ − 1)(σ − a)
(A.1)
where a 6= 0, 1 is generally a complex parameter and the complex-valued exponent parameters α, β, γ, δ, and satisfy the (Fuchsian) condition (see, e.g., [22, 24]): 1+α+β =γ+δ+,
(A.2)
and f is the accessory parameter (it does not affect the exponent parameters). The Heun equation is invariant under α ←→ β.
Appendix A.1. General solution of the triconfluent Heun equation The triconfluent Heun equation is derived from the Heun equation by the coalescence of the three finite regular singular points with infinity ([24]). The canonical triconfluent Heun equation has the form (see, e.g., [24],[23]): dϕ(σ) d2 ϕ(σ) − (3σ 2 + γ) + [α + (β − 3)σ]ϕ(σ) = 0 , 2 dσ dσ and its symmetric (self-adjoint) used form is: d2 Φ(σ) γ2 3 2 9 4 + α− + βσ − γσ − σ Φ(σ) = 0 , dσ 2 4 2 4 where
σ 3 + γσ Φ(σ) = exp − 2
ϕ(σ) .
(A.3)
(A.4)
(A.5)
The quartic polynomial in Eq. (A.4) does not contain σ 3 term, and the coefficient of σ 4 is set equal to 9/4. This normalization can be arranged by subjecting the independent variable σ to an innocuous transformation of the form 22
σ 0 = d1 σ + d0 (d0 and d1 are constants), but other normalizations could equally well be used. The normalization herein is commonly used, see e.g., [24]. An analytical solution ϕ(σ) of Eq. (A.3) can be expressed using a power series of the form ϕ(σ) =
∞ X
(A.6)
bn σ n .
n=0
The derivatives of ϕ with respect to σ are: ∞ ∞ X dϕ(σ) X = nbn σ n−1 = (n + 1)bn+1 σ n , dσ n=0 n=0 ∞ ∞ X d2 ϕ(σ) X n−2 = n(n − 1)bn σ = (n + 2)(n + 1)bn+2 σ n . dσ 2 n=0 n=0
(A.7)
(A.8)
By substituting the power series (A.6), (A.7), and (A.8) into Eq. (A.3), we obtain the recurrence relation: n(n − 1)bn − (n − 1)γbn−1 + αbn−2 + (β + 6 − 3n)bn−3 = 0 .
(A.9)
Shifting the summation index leads to: (n + 2)(n + 3)bn+3 − (n + 2)γbn+2 + αbn+1 + (β − 3 − 3n)bn = 0 , (A.10) where b−1 = 0 , b0 = 1 , b1 = 0 , b2 = −
α . 2
(A.11)
From the recurrence relation (A.9), we can write: bn =
γ(n − 1)bn−1 − αbn−2 − (β + 6 − 3n)bn−3 , n(n − 1)
n≥3.
(A.12)
Additional information regarding this recursive relation is provided in the article [25] where the authors deal with an investigation of the asymptotics of a linearly independent set of solutions for a four-term recurrence relation controlling the coefficients appearing in the power series expansion for the triconfluent Heun functions.
23
For (α, β, γ) ∈ C, the triconfluent Heun function represents a solution of the triconfluent Heun equation (A.3) and can be written as THF(α, β, γ; σ) =
∞ X
bn σ n ;
(A.13)
|σ| < 1 ,
n=0
where bn are given by (A.11) and (A.12). The derivative of the triconfluent Heun function (the prime triconfluent Heun function) can be expressed as ∞ X d THF(α, β, γ; σ) ≡ THF0 (α, β, γ; σ) = nbn σ n−1 ; dσ n=0
|σ| < 1 .
(A.14)
We often need the second derivative of the triconfluent function. For this purpose, we use Eq. (A.1) and we obtain: d2 THF(α, β, γ; σ) = (3σ 2 + γ)THF0 (α, β, γ; σ) dσ 2 − [α + (β − 3)σ]THF(α, β, γ; σ) . (A.15) The general solution of the triconfluent Heun equation (A.3) is (see e.g., [24]): ϕ(σ) = C1 THF(α, β, γ; σ) + C2 exp(σ 3 + γσ)THF(α, −β, γ; −σ) ,
(A.16)
where C1 and C2 are integration constants. There is the following identity: THF(α, β, 0; σ) = exp(σ 3 )THF(α, −β, 0; −σ) ,
(A.17)
which means that if γ = 0 then the presented solutions in Eq. (A.16) are not linearly independent. We can obtain the general solution of the symmetric form of the triconfluent equation (A.4) using the solution (A.16) and the transformation relation (A.5): σ 3 + γσ Φ(σ) = C1 exp − THF(α, β, γ; σ)+ 2 3 σ + γσ THF(α, −β, γ; −σ) . C2 exp 2 (A.18)
24
Thus, in general, the triconfluent Heun differential equation is, by definition, a second-order ordinary differential equation of the form d2 Φ + Π(σ)Φ(σ) = 0 , dσ 2
(A.19)
where the coefficient function Π(σ) is some quartic polynomial in the independent variable σ. In principle, the coefficient function Π(σ) could be considered to be a polynomial of a higher degree than four. Solutions Φ(σ) of Eq. (A.19) when the function Π(σ) is a quintic, sextic or higher polynomial, would be special functions of new types. Since the classification scheme has not yet been extended to the higher degrees of polynomials than four, and thus is not part of the mathematical software such as Maple, we have limited ourselves in this work to the quartic polynomials. Appendix A.2. Polynomial solution of the triconfluent Heun equation We can see that the evaluation of successive terms of a series solution (A.6) or (A.13) to the triconfluent Heun equation (A.3) is carried out by means of a recurrence relation (A.12) where the coefficient bn depends upon n, the previous values of br (r < n), and the three parameters (α, β, γ). Based on this recursive relation, it is possible to express a solution of the triconfluent equation by a finite polynomial if certain conditions are met. The conditions for the existence of polynomial solutions can be formulated by the following theorem (see e.g., [24, 25]): Theorem 1. Suppose that in the triconfluent Heun equation (A.3) the parameters α, β and γ satisfy the following two conditions: (a) β = 3(N + 1),
N = 0, 1, 2 . . . ,
(b) The determinant DN +1 (α, γ) = 0 of the matrix
25
MN +1
α
3N 0 0 = .. . 0
−γ
2·1
0
0
···
α
−2γ
2·3
0
···
3(N − 1)
α
−3γ
3·4
···
0 .. .
3(N − 2) .. .
α .. .
−4γ .. .
··· ..
0
···
0 0 . .. . −(N − 1)γ N (N − 1) α −N γ 3·1 α (A.20) 0
. 3·3
α
0
3·2
0
0
Then the triconfluent Heun equation (A.3) has a polynomial (Liouvillian) solution ϕ(σ) = P β −1 (σ) , 3
where P β −1 (σ) denotes a polynomial of degree 3
β 3
− 1, whose coefficients pn
(n = 0, 1, . . . , β3 − 1) are the solutions of the following linear system of equations
p0
Mβ 3
p1 .. .
=0.
p β −1
(A.21)
3
If ϕ1 (σ) is the polynomial solution of the triconfluent Heun equation (A.3) then a second independent solution ϕ2 (σ) can be calculated by the Wronskian method (see e.g., [26]) Z
σ
ϕ2 (σ) = ϕ1 (σ)
exp
R w
(3v 2 + γ)dv dw . ϕ1 (w)2
(A.22)
References [1] A. G. Webster, Acoustical impedance, and the theory of horns and of the phonograph, Proceedings of the National Academy of Sciences 5 (1919) 275–282. doi:10.2307/84052. [2] L. L. Beranek, T. J. Mellow, Acoustics: Sound Fields and Transducers, Academic Press, 2012.
26
[3] A. Chaigne, J. Kergomard, Acoustics of Musical Instruments, SpringerVerlag, 2016. [4] M. V. Guiman, I. C. Rosca, A new approach on vibrating horns design, Shock and Vibration 8532021 (2017) 1–12. doi:doi.org/10.1155/2017/ 8532021. [5] K. Hornisova, P. Billik, Some properties of horn equation model of ultrasonic system vibration and of transfer matrix and equivalent circuit methods of its solution, Ultrasonics 54 (2014) 330 – 342. doi:https: //doi.org/10.1016/j.ultras.2013.05.003. [6] V. S. Sorokin, Effects of corrugation shape on frequency band-gaps for longitudinal wave motion in a periodic elastic layer, The Journal of the Acoustical Society of America 139 (2016) 1898–1908. doi:10.1121/1.4945988. [7] P. Honzík, S. Durand, N. Joly, M. Bruneau, On the acoustic transfer function of slowly tapered small horns filled with thermo-viscous fluid, Acta Acustica united with Acustica 99 (2013) 694–702.
doi:10.3813/aaa.
918648. [8] Y. A. Ilinskii, Nonlinear standing waves in an acoustical resonator, The Journal of the Acoustical Society of America 104 (1998) 2664–2674. doi: 10.1121/1.423850. [9] O. Rudenko, A. Shvartsburg, Nonlinear and linear wave phenomena in narrow pipes, Acoustical Physics 56 (4) (2010) 429–434. doi:10.1134/ S1063771010040044. [10] M. Cervenka, M. Soltes, M. Bednarik, Optimal shaping of acoustic resonators for the generation of high-amplitude standing waves, The Journal of the Acoustical Society of America 136 (2014) 1003–1012. doi: 10.1121/1.4892751.
27
[11] M. Cervenka, M. Bednarik, Acoustic bandpass filters employing shaped resonators, Journal of Sound and Vibration (2016) 76–88doi:10.1016/j. jsv.2016.06.045. [12] E. Eisner, Complete solutions of the "Webster" horn equation, The Journal of the Acoustical Society of America 41 (1967) 1126–1146. doi:10.1121/ 1.1910444. [13] S. Abrate, Vibration of non-uniform rods and beams, Journal of Sound and Vibration 185 (4) (1995) 703 – 716. doi:doi.org/10.1006/jsvi.1995. 0410. [14] B. Kumar, R. Sujith, Exact solutions for the longitudinal vibration of nonuniform rods, Journal of Sound and Vibration 207 (1997) 721 – 729. doi: doi.org/10.1006/jsvi.1997.1146. [15] Q. Li, Exact solutions for free longitudinal vibration of stepped non-uniform rods, Applied Acoustics 60 (1) (2000) 13 – 28. doi:doi.org/10.1016/ S0003-682X(99)00048-1. [16] A. Raj, R. Sujith, Closed-form solutions for the free longitudinal vibration of inhomogeneous rods, Journal of Sound and Vibration 283 (3) (2005) 1015 – 1030. doi:https://doi.org/10.1016/j.jsv.2004.06.003. [17] B. Yardimoglu, L. Aydin, Exact longitudinal vibration characteristics of rods with variable cross-sections, Shock and Vibration 18 (2011) 555–562. doi:10.3233/SAV-2010-0561. [18] B. A. Donskoy, Dimitri M.; Cray, Acoustic particle velocity horns, The Journal of the Acoustical Society of America 131 (2012) 3883–3890. doi: 10.1121/1.3702432. [19] A. D. Pierce, Acoustics: an Introduction to its Physical Principles and Applications, McGraw-Hill Book Company, Inc., 1981.
28
[20] S. W. Rienstra,
Webster’s horn equation revisited,
nal on Applied Mathematics 65 (2005) 1981–2004.
SIAM Jourdoi:10.1137/
s0036139902413040. [21] P. A. Martin, On Webster’s horn equation and some generalizations, The Journal of the Acoustical Society of America 116 (2004) 1381. doi:10. 1121/1.1775272. [22] G. Kristensson, Second Order Differential Equations:
Special Func-
tions and Their Classification, Springer New York, 2010. doi:10.1007/ 978-1-4419-7020-6. [23] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, Handbook of Mathematical Functions, Cambridge University Press, 2010. [24] A. Ronveaux (Ed.), Heun’s Differential Equations, Oxford Science Publications, 1995. [25] D. N. M. Batic, D.; Mills-Howell, Potentials of the Heun class: The triconfluent case, Journal of Mathematical Physics 56 (2015) 052106. doi: 10.1063/1.4921344. [26] K. Riley, M. Hobson, S. J. Bence, Mathematical Methods for Physics and Engineering, Cambridge University Press, 2006.
29
Declaration of interests ⊠ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
PII:
S0022-460X(19)30732-1
DOI:
https://doi.org/10.1016/j.jsv.2019.115169
Reference:
YJSVI 115169
To appear in:
Journal of Sound and Vibration
Received Date: 15 April 2019 Revised Date:
20 December 2019
Accepted Date: 25 December 2019
Please cite this article as: M. Bednarik, M. Cervenka, A wide class of analytical solutions of the Webster equation, Journal of Sound and Vibration (2020), doi: https://doi.org/10.1016/j.jsv.2019.115169. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Credit Author Statement
Michal Bednarik conceived the presented analytical method, wrote the text of the manuscript and performed all the analytical calculations. Milan Cervenka verified all presented results based on numerical calculations, drew the figures in the manuscript and helped with some statements in the text of the manuscript. Both authors discussed the results and commented on the manuscript.
A wide class of analytical solutions of the Webster equation M. Bednarika,∗, M. Cervenkaa a Czech
Technical University in Prague, Faculty of Electrical Engineering, Technicka 2, 166 27 Prague 6, Czech Republic
Abstract This paper aims at presenting closed-form general analytical solutions of the Webster equation describing plane elastic or acoustic waves. The considered radius functions of non-uniform cross-sectioned rods or ducts are based on the triconfluent Heun functions and contain some optional parameters enabling us to set various profiles of the radius functions in a relatively wide range, while it is possible to employ the presented exact general analytical solution of the Webster equation for all selected profiles. If the radius functions are predetermined, then the derived general analytical solution can also be employed for their triconfluent Heun approximations, including certain polynomial ones. The applicability and correctness of the derived analytical solutions are demonstrated by calculations of natural frequencies and mode shapes for representative radius functions while the results based on approximate analytical solutions are verified numerically. Keywords: Triconfluent Heun equation, Triconfluent Heun function, Webster equation, non-uniform cross-sectioned waveguide, horn
1. Introduction There are different technological devices for the transmission of acoustic energy into an acoustic load which are based on waveguiding systems. We can meet them, for example, in heat and mass transfer processes, plastic welding, ∗ Corresponding
author Email address: [email protected] (M. Bednarik)
Preprint submitted to Journal of LATEX Templates
December 26, 2019
metal treatment, ultrasonic sonochemistry, acoustic transmission lines, phononic and sonic crystals etc. Under simplifying assumptions, we can describe the propagation of sound (or elastic waves) in waveguiding systems with variable cross-sections as pipes, horns, concentrators, ultrasonic rod waveguide-radiators, by the Webster equation (see e.g. [1, 2, 3]). The solution of the Webster equation enables us to optimally design the waveguiding systems with respect to the characteristics of our interest. The performance of many ultrasonic devices depends on a proper design of a sonotrode shape (concentrator). For this reason, it is important to known the dependence of fundamental modal properties (natural frequencies, mode shapes) of sonotrode shapes on their geometrical parameters (see e.g. [4]). Also in the design of ultrasonic rod waveguide-radiators, it is desirable to know the solution of the Webster equation (see e.g. [5]). Locally periodic structures (phononic or sonic crystals) composed of periodically repeated elements are used in many fields of science and technology (vibration isolators and mechanical filters to building frames, bridge trusses, railway tracks, compound pipes etc.) (see e.g. [6]). Therefore, the determination of widths and locations of band-gaps and the corresponding attenuation levels associated with these systems represents an important practical problem. Furthermore, we may encounter the need to know the solution of Webster equation in musical acoustics in exploring and designing wind instruments, where the main function of the horn (bell) is to form a more gradual impedance transition between the high inner tube impedance and the low outside air impedance (see e.g. [3]). The Webster equation has its use in the design of acoustical horns applied to loudspeakers (see e.g. [2]) and it is helpful in designing very small (even miniaturized) components (e.g. small transducers, earplugs) (see e.g. [7]). The solution of the Webster equation also finds application in acoustic filters and acoustic resonators (see e.g. [8, 9, 10, 11]). However, to the best of the authors’ knowledge, only comparatively few publications are available that present the exact analytical solutions of the Webster equation for various cross-section profiles of waveguides (see e.g. [12, 13, 14, 15, 2
16, 17]). Therefore, the primary motivation of this work is to present both new exact analytical solutions of the Webster equation and a method that enable us to find approximate analytical solutions of this equation for the cases when the exact ones are not known. To accomplish this intention, we used cross-section profiles based on the triconfluent Heun functions whose evaluation is already a standard part of mathematical software such as Maple. The functions have attracted the attention of scientists mainly during the last decade. Nowadays, we can observe a growing interest in their utilization in various fields of physics even if they are still waiting for their broader application particularly in acoustics and elastodynamics. The advantage of these functions is that they contain a sufficient number of parameters for modification of the cross-section profiles, without the need to employ many analytical solutions of the Webster equation. As shown below, it is possible in this way to vary shapes and sizes of the cross-section profiles in a comparatively wide range. For example, using these parameters it is possible to set the cross-section profiles of rods that are changing smoothly and continuously (e.g. with zero derivatives at the endpoints) along a desired dimension to avoid failures which are likely to initiate at the places where changes in material properties are not gradual. This paper is organized as follows. The Webster equation for elastic and acoustic waves is presented in Section 2. Based on mapping the Schrödinger-like differential equation into the Webster equation, its general solution is derived in Section 3. The radius-function profiles expressed using the triconfluent Heun functions for which exact analytical solutions of the Webster equation exist are given in Section 4. The possibilities of expressing the radius-functions based on polynomial solutions of the triconfluent Heun differential equation are also presented in this section. Section 5 is devoted to approximate analytical solutions of the Webster equation for predetermined radius-functions. The applicability and correctness of the derived analytical solutions are demonstrated by calculations of natural frequencies and mode shapes for representative radius-functions in Section 6. Finally, Section 7 provides our conclusions. In order to make this paper sufficiently self-contained and more clear, some features concerning the 3
triconfluent Heun differential equation and its solutions including polynomial ones are outlined in Appendix.
2. Webster equation We consider plane elastic wave propagation through an isotropic non-uniform cross-sectioned slender rod or tube, see Fig. 1 where the geometry is specified. Assuming small vibration amplitudes and small transverse dimensions then longitudinal propagation of elastic waves can be described by means of the following model equation (in the absence of body force) (see e.g. [5]) ∂ ∂u(z, t) ∂ 2 u(z, t) , EA(z) = ρs A(z) ∂z ∂z ∂t2
(1)
where u(z, t) is the displacement of the rod along the length direction from its equilibrium position at point z and time t. Further, E, A(z), and ρs are Young’s modulus, cross-sectional area function and density of the rod, respectively. If we assume linear plane wave acoustic propagation along a narrow fluid-filled tube of gradually varying cross-section then the corresponding model equation can be expressed as (see e.g. [18]) ∂ ∂p(z, t) ∂ 2 p(z, t) BA(z) = ρA(z) , ∂z ∂z ∂t2
(2)
where p(z, t), ρ and B are the acoustic pressure, fluid density and the adiabatic bulk modulus, respectively. Limiting ourselves to the harmonic solutions of Eqs. (1) and (2), i.e. u(z, t) = √ U (z) exp(−jωt) and p(z, t) = P (z) exp(−jωt), where j = −1 is the imaginary unit, and ω is the angular frequency, then the model equations (1) and (2) reduce to the Webster equation: d2 Φ(z) 1 dA(z) dΦ(z) + + k 2 Φ(z) = 0 . 2 dz A(z) dz dz
(3)
Here Φ(z) stands for either U (z) or P (z), k = ω/c is the wave number and c is p the wave speed. The wave speed for elastic waves is c = ρs /E and for acoustic p waves c = ρ/B. 4
r(z)
A(z) = πr2(z)
0
H
z
Figure 1: Rod/tube with variable cross-section.
Using the Webster equation for plane elastic and acoustic waves has its limitations that are discussed in, e.g. [5, 19, 20, 21]. Employing the Webster equation enables us to treat propagation of elastic or acoustic waves as an one-dimensional problem if the characteristic width of a rod or tube is small compared to the wavelength and the cross section A(z) varies relatively slowly. Considering an axially symmetrical waveguide, i.e. A(z) = πr2 (z) (r(z) is a radius-function) then, for convenience, we can rewrite Eq. (3) into the following form:
d2 Φ(s) 2 dR(s) dΦ(s) + + K 2 Φ(s) = 0 , ds2 R(s) ds ds
(4)
where R(s) = r(s)/r0 (dimensionless radius-function), s = z/`, K = k`, ` is a characteristic length, and r0 = r(0).
3. General analytical solution of the Webster equation We can solve the Webster equation (4) based on a time-independent onedimensional Schrödinger-like differential equation d2 ψ(s) 2 + K + G(s) ψ(s) = 0 , 2 ds 5
(5)
that can be mapped into Eq. (4) by the following transformation (6)
ψ(s) = R(s)Φ(s) . Substituting this transformation relation into Eq. (4) we obtain d2 Φ(s) 2 dR(s) dΦ(s) Φ(s) d2 R(s) 2 + + K + G(s) Φ(s) = 0 . + 2 2 ds R(s) ds ds R(s) ds
(7)
By identifying this equation with Eq. (4) we obtain the following condition on the transformation
d2 R(s) + G(s)R(s) = 0 . (8) ds2 Thus, searching for the general solution of the Webster equation (4) is reduced to finding a general solution of the Schrödinger-like equation (5). 3.1. General solution of the Schrödinger-like equation We can see from Eq. (8) that the (generating) function G(s) determines the radius function R(s). If we were able to solve Eq. (5) for a given function G(s), then using the transformation relation (6) we could solve the Webster equation (4) for the cross-sectional area function A(s) = πr02 R(s)2 . It is obvious that for G = a ≥ 0 one can very easily solve Eqs. (5) and (8). The corresponding radius-functions are presented, e.g. in the works [9, 12, 5]. Employing more complex generating functions other exact analytical solutions can be found, e.g. in the works [14, 16, 17]. In this work, we assume that the generating function G(s) is a quartic polynomial, which enables us to solve Eq. (8) exactly and in addition, such a generating function contains a sufficient number of coefficients that allow us to determine an appropriate radius profile. The generating function can be expressed as G(s) = a0 + a1 (s − s0 ) + 2a2 (s − s0 )2 + a3 (s − s0 )3 − a4 (s − s0 )4 ,
a4 6= 0 . (9)
Substituting this generating function into Eq. (5) we obtain d2 ψ(s) 2 + K + a0 + a1 (s − s0 ) + 2a2 (s − s0 )2 + ds2 a3 (s − s0 )3 − a4 (s − s0 )4 ψ(s) = 0 . (10) 6
In order to find the solution of Eq. (10) we first employ the following transformation of the dependent variable ψ(s) = exp
Z 1 F (s)ds ϕ(s) 2
(11)
for Eq. (10) and we attain d2 ϕ dϕ F (s)2 1 dF (s) + F (s) + + + K 2 + a0 + a1 (s − s0 )+ 2 ds ds 4 2 ds 2a2 (s − s0 )2 + a3 (s − s0 )3 − a4 (s − s0 )4 ϕ(s) = 0 . (12) Substituting 2 √ a3 16a2 a4 + 3a23 + F (s) = −2 a4 (s − s0 ) − 3 4a4 8a42
(13)
into Eq. (12) we obtain the form of the triconfluent Heun differential equation ( ) 2 √ d2 ϕ a3 16a2 a4 + 3a23 dϕ − 2 a4 (s − s0 ) + + − 3 ds2 4a4 ds 8a42 " 5 8a1 a24 + 8a2 a3 a4 + a33 − 16a42 (s − s0 )+ 8a24 # 5 64(K 2 + a0 )a34 + 64a22 a24 + 16a2 a23 a4 + 32a3 a42 + a43 ϕ(s) = 0 . (14) 64a34 Using the following transformation for an independent variable: a3 σ = Q (s − s0 ) − , 4a4 where
√ 13 2 a4 Q= , 3
(15)
(16)
we can rewrite Eq. (14) into its canonical form: ! 2 d2 ϕ(σ) 16a a + 3a dϕ(σ) 2 4 3 + − 3σ 2 − 3 dσ 2 dσ 2 8Qa4 64(K 2 + a0 )a34 + 64a22 a24 + 32a2 a33 a4 + 16a1 a3 a24 + 3a43 + 64Q2 a34 # 5 8a1 a24 + 8a2 a3 a4 + a33 − 16a42 σ ϕ(σ) = 0 . (17) 8Q3 a24 7
By comparing Eq. (17) with the canonical triconfluent Heun equation (A.3) in the Appendix, we obtain α=
64(K 2 + a0 )a34 + 64a22 a24 + 32a2 a23 a4 + 16a1 a3 a24 + 3a43 , 64Q2 a34 24a4 (a1 a4 + a2 a3 ) + 3a33 16a2 a4 + 3a23 β= , γ=− . (18) 5 3 16a42 8Qa42
Employing Eq. (A.16) we can write the general solution of Eq. (17) as ϕ(σ) = A1 THF (α, β, γ; σ) + A2 exp σ 3 + γσ THF (α, −β, γ; −σ) , (19) where A1 and A2 are integration constants and THF stands for the triconfluent Heun function. Using the relation (A.18) and (15) we can write the general solution (integration constants are arbitrary) of Eq. (10) in the following form:
3 Q s− ψ(s) = A1 exp −
a3 4a4
− s0
3
+ γQ s − 2
a3 4a4
− s0 ×
a3 THF α, β, γ; Q s − − s0 + 4a4 3 a3 a3 3 Q s − 4a4 − s0 + γQ s − 4a4 − s0 A2 exp × 2 a3 − s0 THF α, −β, γ; −Q s − . (20) 4a4 3.2. General solution of the Webster equation for the specific class of radius functions Employing the transformation relation (6) and the solution (20) it is possible to write the general closed-form solution of the Webster equation (4) for the
8
specific class of the radius functions that is given in the Sec. 4, as 3 3 Q a3 a3 exp − Q2 s − 4a − s − γ s − − s 0 0 2 4a4 4 Φ(s) = A1 × R(s) a3 THF α, β, γ; Q s − − s0 + 4a4 3 Q3 Q a3 a3 exp 2 s − 4a4 − s0 + γ 2 s − 4a4 − s0 × A2 R(s) a3 THF α, −β, γ; −Q s − − s0 4a4 = A1 Φ1 (s) + A2 Φ2 (s) . (21) Based on what has been shown above, there are two approaches to solving the Webster equation. The first approach consists in solving Eq. (8) for a predetermined generating function G(s), which leads to obtaining the radiusfunction R(s). It is obvious that if we know the solution of Eq. (5) at the same time, then analytical solution (21) is exact. The second approach, which is often convenient from practical point of view, is based on a predetermined radius-function R(s) when the function G(s) is calculated by using Eq. (8). However, to make the analytical solution (21) exact, it is necessary to solve Eq. (5) exactly as well. Unfortunately, there are only a few radius-functions for which we can solve the Webster equation exactly. To avoid this fact we can employ a quartic polynomial that approximates the function G(s). In this case, the Webster equation solution (21) is no longer exact, though it is an analytical one. 4. Radius functions for the exact solution of the Webster equation One of the intended aims of this work is to extend the group of radius profiles for which exact solutions of the Webster equation are known. Based on the results from Sec. 3 we can accomplish this task by employing the generating function G(s) in the form the quartic polynomial (9). However, from a practical point of view, it is better to reduce the number of degrees of freedom in the 9
search for the radius function by omitting some of lower terms of the quartic polynomial (a4 > 0). For this reason, we will no longer consider the cubic term in this polynomial to reduce the number of optional parameters. With respect to the identity (A.17)) we assume that a2 6= 0. Substituting the generating function (9) with a3 = 0 into Eq. (8) we obtain an equation that is of the same form as Eq. (10). This fact enables us to directly utilize the solution (20) where we substitute for K 2 and a3 the zero value. After small algebraic modifications we can write the solution of Eq. (8) as 2 3 s a2 R(s) = C1 exp − Q3 s − ss0 + s20 − × 2 3 a4 2 a2 + a0 a4 3a1 3Q2 a2 THF , , − ; Q(s − s ) + √ 0 a4 Q2 2 a4 a4 2 s a2 3 − ss0 + s20 − × C2 exp Q3 s 2 3 a4 2 a2 + a0 a4 3a1 3Q2 a2 THF , − , − ; −Q(s − s ) = √ 0 a4 Q2 2 a4 a4 C1 R1 (a0 , a1 , a2 , a4 , s0 ; s) + C2 R2 (a0 , a1 , a2 , a4 , s0 ; s) . (22) We can see that the derived radius-function (22) can be adjusted based on 7 optional parameters, i.e. two integration constants C1,2 and a group of five internal parameters (a0 , a1 , a2 , a4 , s0 ). This fact enables us to find an entire family of radius-functions of various profiles for which we have one exact solution of the Webster equation (4) by substituting the relation (22) for R(s) in the solution (21). If we want to determine the radius at the point s = 0 in order to correspond to the radius r0 then we choose values of the integration constants C1,2 in such a way that the radius-function R(s = 0) = 1. In this case, required values of the radius-function R(s) at another point si ∈ (0, sH ], where sH = H/` is a dimensionless length of the rod or tube, are given only by the choice of the optional parameters. To do that, we can change one of the optional parameters, and we fix the remaining ones until the radius function reaches the required value at the given point si . It is worth noting that if the argument of the triconfluent Heun functions is 10
equal to zero then these functions equal unity for each of the parameters α, β, γ. Additionally, the first derivative of the Heun functions are always equal to zero if their arguments are zero. Examples of radius-function profiles for various optional parameters are shown in Fig. (2)(a). However, if we require to determine a radius r0 at the point s = 0 and at
(a)
(b)
1.25
1.25 R(s)
1.5
R(s)
1.5 1
1
0.75
0.75
0.5
0.5 0
0.5
1 s
1.5
2
0
0.5
1 s
1.5
2
Figure 2: (a) Examples of radius-function profiles R(s) based on the formula (22) for a given sequence of optional parameters (C1 , C2 , a0 , a1 , a2 , a4 , s0 ): (1,0,-0.1,0,0.01,0.08,0) – black line, (0,1,-1.99,0,-0.04,0.001,0) – blue line, (1,0,-3,0,-0.05,0.001,0) – red line, (1,0,0,0,0.01,0.001,0) – brown line; (b) Examples of radius-function profiles R(s) based on the formula (23) for a given sequence of optional parameters (a0 , a1 , a2 , a4 , s0 , si , ri /r0 ): (-0.7,2,0.2,0.66,0,2,3) – black line, (0,0,2.32856,0.85,1,1,3) – red line, (0,0,-0.01,0.045,1,2,0.5) – blue line, (0,0,-1.5,0.55,1,2,0.5) – brown line.
the same time a radius ri at the point si ∈ (0, sH ] with fixed optional internal parameters then we need to employ both the constants C1,2 that depend on a value of both the radii r0 and ri . In this case, we can write the radius-function (22) as R(s) =
ri r0 R2 (0)
− R2 (si )
W
R1 (s) −
ri r0 R1 (0)
− R1 (si )
W
R2 (s) ,
(23)
where W = R2 (0)R1 (si ) − R2 (si )R1 (0) . Here for brevity, we do not write the parameters in functions R1,2 (s). 11
(24)
4.1. Some polynomial solutions for the radius function Due to the fact that it is relatively difficult to estimate the profiles of the triconfluent Heun functions based on their parameters α, β, γ, it seems to be useful to express them as terminating power series, i.e., polynomials, if certain conditions are satisfied, see Appendix A.2. Based on the theorem presented in Appendix A.2 where the meaning of the used symbols below is also given, it is possible to write the first five polynomial solutions of the triconfluent equation: 1. For N = 0 (β = 3), we can directly find that D1 ≡ α = 0 and P0 = p0 and the triconfluent function equals to unit, i.e. THF(0, 3, γ; σ) = 1 ;
(25)
2. For N = 1 (β = 6), D2 ≡ α2 + 3γ = 0, P1 (σ) = σ −
α ; 3
3. For N = 2 (β = 9), D3 ≡ α3 + 12αγ + 36 = 0, P2 (σ) = σ 2 −
α α2 1 σ+ − ; 3 36 α
4. For N = 3 (β = 15), D4 ≡ α5 + 60α3 γ + 756α2 + 576αγ 2 + 5184γ = 0, γ α2 7 2 α α3 + − αγ − ; P3 (σ) = σ 3 − σ 2 + σ− 3 2 18 162 54 3 5. For N = 4 (β = 12), D5 ≡ α4 + 30α2 γ + 216α + 81γ 2 = 0, α P4 (σ) = σ − σ 3 + 3 4
where α=
a22 + a0 a4 , a4 Q2
3 α2 2 α 5 4 2 + γ σ − + αγ + σ+ 18 3 162 27 3 α4 2 5 γ2 + α2 γ + α + , 1944 81 18 9
γ=−
3Q2 a2 a4
and σ = Q(s − s0 ) .
(26)
It should be noted that the polynomial solutions for N ≥ 1 do not represent triconfluent Heun functions according to the commonly used definition because they do not satisfy the condition b1 = 0, as it is supposed in (A.11). 12
For N > 4, we can continue finding the couples (α, γ) that enable us to express the polynomial solutions but higher order polynomials lose their usefulness in the sense of estimation of their profiles and therefore we do not consider them here. Using the first two polynomial solutions above, we obtain the following expressions for the radius function R(s): 1. N = 0;
3 R(s) = C1 exp − Q3 s 2
s2 a2 − s0 s + s20 − 3 a4
(27)
.
2. N = 1; 2 s a2 a22 + a0 a4 3 3 2 − s0 s + s0 − . Q(s − s0 ) − R(s) = C1 exp − Q s 2 3 a4 3a4 Q2 (28)
(a)
(b)
1.5
1.75 1.5 1.25
1
R(s)
R(s)
1.25
1
0.75
0.75
0.5
0.5 0
0.5
1 s
1.5
2
0
0.5
1 s
1.5
2
Figure 3: (a) Examples of radius-function profiles R(s) based on the formula (27) for a given sequence of optional parameters (C1 , a2 , a4 , s0 ): (1,0.161,0.12075,0) – black line, (1,0.005,0.0011,0) – blue line, (1,0.005,0.008,1) – red line, (1,1,6,1) – brown line; (b) Examples of radius-function profiles R(s) based on the formula (28) for a given sequence of optional parameters (C1 , a0 , a2 , a4 , s0 ): (0.5510,-2,-0.1,0.0947,0) – black line, (-0.3712,3,0.1,0.3,1) – red line, (0.5166,-3,0.1,0.3,0) – blue line, (0.9927,-2,0,0.15,1) – brown line.
The polynomial solutions of Eq. (8) for N ≥ 1 may be expressed as linear combinations of the two standard solutions (22). 13
It is necessary to note that it is not possible to express the independent function R2 (s) in the polynomial form for the same optional parameters. Thus, for the second independent solution, it is necessary to use the corresponding triconfluent Heun function or to employ the relation (A.22) if we are able to integrate it. For this reason, it is more convenient to employ the polynomial solutions mainly for the method based on setting suitable optional parameters. The material functions based on the first two polynomial solutions enable us to calculate the corresponding optional parameters analytically, some results being shown in Fig. 3. For the remaining polynomials, it is easier to seek the values of the respective parameters numerically.
5. Approximate analytical solutions of the Webster equation If the radius function R(s) is predetermined then the function G(s) can be calculated based on Eq. (8), i.e. G(s) = −
1 d2 R(s) . R(s) ds2
(29)
However, we generally cannot solve Eq. (5) for the function G(s) given by Eq. (29). There are predetermined radius functions for which we can approximate this function by a quartic polynomial with sufficient accuracy, i.e. G(s) ≡ −
1 d2 R(s) u W (s) = R(s) ds2 a0 + a1 (s − s0 ) + 2a2 (s − s0 )2 + a3 (s − s0 )2 − a4 (s − s0 )4 .
(30)
The coefficients of the quartic polynomial W (s) can be calculated for instance based on an interpolation or by a least-squares approximation. Replacing the function G(s) by the polynomial W (s) in Eq. (5) we can employ the function (21) as the analytical solution of Eq. (4). In other words, we substitute coefficients a0 , a1 , a2 , a3 and a4 from the quartic polynomial (30) into the solution (21).
14
6. Natural frequencies and mode shapes for non-uniform rods Without loss of generality, we limit ourselves only to elastic waves (Φ(s) ≡ U (s)) in this section, where the presented procedures are also directly applicable to acoustic waves in fluids. Based on the following representative examples, we demonstrate the applicability and correctness of the derived solutions. 6.1. Representative examples based on the exact solutions We can calculate natural frequencies based on the solution (21) for the following three classical boundary conditions: 1. Clamped-clamped rod (CCR): U (0) = A1 U1 (0) + A2 U2 (0) = 0 , U (sH ) = A1 U1 (sH ) + A2 U2 (sH ) = 0 ; (31) 2. Clamped-free rod (CFR): U (0) = A1 U1 (0) + A2 U2 (0) = 0 , dU (s) dU1 (s) dU2 (s) = A + A =0; 1 2 ds s=sH ds s=sH ds s=sH
(32)
3. Free-free rod (FFR): dU2 (s) dU1 (s) dU (s) = A1 + A2 =0, ds s=0 ds s=0 ds s=sH dU2 (s) dU (s) dU1 (s) + A2 = A1 ds ds ds s=sH
s=0
= 0 . (33) s=sH
For all the reported radius-functions below, which are found based on the predetermined generating functions (9), we assume the internal parameters a3 = s0 = 0 and the endpoint sH = 1. In addition, we select the radiusfunctions with the zero derivatives at the endpoints and R(0) = 1. To meet these requirements for the endpoint s = 0 we choose the integration constants C1 = C2 = 1/2 in the relation (22). The zero derivative of the radius function at the endpoint s = sH is accomplished by numerical calculation of the internal
15
parameter a1 by fixing the predetermined parameters a0 , a2 6= 0, and a4 > 0. Thus, the radius-function reads 2 1 3 a2 s R(s) = exp − Q3 s − THF (α, β, γ; Qs) + 2 2 3 a4 2 1 3 3 a2 s exp Q s − THF (α, −β, γ; −Qs) , (34) 2 2 3 a4 dR(s) dR(s) =0, = R(0) = 1 , ds s=0 ds s=sH where Q=
√ 13 2 a4 a2 + a0 a4 3a1 3Q2 a2 , α= 2 , β= √ , γ=− . 2 3 a4 Q 2 a4 a4
(35)
(36)
Based on the solution (21) and the radius function (34) we can, after some algebra, express the general solution for the radius-function (34) as U (s) = A1 U1 (s) + A2 U2 (s) = 2 2THF α + K , β, γ; Qs 2 Q h i A1 + 2 3 THF (α, β, γ; Qs) + exp 3Q s s3 − aa42 THF (α, −β, γ; −Qs) 2 2THF α + K Q2 , −β, γ; −Qs h i A2 . (37) 2 exp −3Q3 s s3 − aa24 THF (α, β, γ; Qs) + THF (α, −β, γ; −Qs) Non-trivial solutions of pairs of homogeneous equations (31)-(33) result in the vanishing of determinants. For this purpose, we can use the following relations that are obtained based on the general solution (37): dU1 (s) dU2 (s) a2 a2 U1 (0) = 1 , U2 (0) = 1 , = √ , = − √ . (38) ds s=0 a4 ds s=0 a4 To meet the condition that the determinants of the corresponding pairs of homogeneous equations are zero, then taking into account the relations (38) we have to solve the following transcendental equations for eigenvalues K: 1. Clamped-clamped rod (CCR): 1 a2 K2 exp 3Q3 − THF α + 2 , −β, γ; −Q − 3 a4 Q K2 THF α + 2 , β, γ; Q = 0 , (39) Q 16
2. Clamped-free rod (CFR): 1 a2 K2 exp 3Q3 − THF α + 2 , −β, γ; −Q − 3 a4 Q 2 K K2 THF0 α + 2 , −β, γ; −Q − THF α + 2 , β, γ; Q − Q Q K2 THF0 α + 2 , β, γ; Q = 0 , (40) Q 3. Free-free rod (FFR): 1 a2 K2 exp 3Q3 − THF α + 2 , −β, γ; −Q − 3 a4 Q 2 K K2 THF0 α + 2 , −β, γ; −Q + THF α + 2 , β, γ; Q + Q Q K2 THF0 α + 2 , β, γ; Q = 0 , (41) Q where the prime stands for the derivative with respect to the variable s, see (A.14). To demonstrate the applicability of the presented results we choose four representative radius-functions that meet the above-mentioned requirements. The chosen radius-functions are shown in Fig. (4) and are specified by the optional parameters that are arranged in parameter vectors v = (C1 , C2 , a0 , a1 , a2 , a3 , a4 , s0 ). First (fundamental) eigenvalues are calculated from the transcendental equations (39)-(41) for the four considered rod radius-functions, see Tab. 1 (the mathematical software Maple (version 17) was used for the evaluations of the triconfluent Heun functions). The eigenvalues K are related to natural frequencies via the relation ω = Kc/`. Three normalized mode shapes for the radius function (see Fig. 5(a)) that has the zero derivatives at the endpoints and takes the same value at these points, i.e. R(0) = R(1), are shown in Fig. 5(b). 6.2. Representative examples based on the approximate solutions As representative examples, we assume the following two predetermined radius functions: R(s) = exp 17
s4 2
,
(42)
1.2
R(s)
1
0.8
0.6
0.4 0
0.2
0.4
0.6
0.8
1
s Figure
4:
R(s)
Radius-functions
(C1 , C2 , a0 , a1 , a2 , a3 , a4 , s0 ):
for
various
vectors
of
parameters
(0.5, 0.5, −7, 25.74836, 10, 0, 70, 0)
=
v1
-
v
black
= line;
v2 = (0.5, 0.5, −7, 34.68600, 2, 0, 70, 0) - blue line; v3 = (0.5, 0.5, 0, 9.16938, 0.1, 0, 30, 0) - red line; v4 = (0.5, 0.5, 0, 6.62284, 0.1, 0, 20, 0) - brown line.
CCR
CFR
FFR
v1
2.4832
1.7573
3.9262
v2
2.4766
1.9611
4.0974
v3
2.9169
2.2463
3.7469
v4
2.9539
2.0434
3.5166
Table 1: First (fundamental) eigenvalues for four rod radius functions represented by corresponding parameter vectors v for the boundary conditions (39)-(41)
R(s) =
1 , 1 + s2 − 0.5s4
(43)
for which we can find approximate analytical solutions employing the solution (21) as it is described in Sec. 5. Considering the radius function (42) we can interpolate the corresponding func-
18
(b)
(a) 1 0.5 U (s)
R(s)
1.3 1.2
0
1.1
-0.5
1
-1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
s Figure 5:
0.6
0.8
1
s
(a) Radius-function R(s) for CFR with the vector of parameters v
=
(0.5, 0.5, −7, 25.1795, 2, 0, 36.9072, 0); (b) The corresponding three mode shapes for the first three eigenvalues: K = 1.4917 - black line; K = 4.8762 - blue line; K = 7.9167 - red line.
tion G(s) from Eq. (29) by the following quartic polynomial G(s) = −
1 d2 R(s) u W (s) = −16.3856s4 +18.6067s3 −13.00832s2 +0.7872s . R(s) ds2 (44)
The radius-function is depicted in Fig. 6(a). Substituting the coefficients of the polynomial (44) into the solution (21) we obtain the approximate analytical solution for longitudinal elastic waves in the slender non-uniform cross-sectioned rod. Employing this approximate solution then we can calculate e.g., the corresponding mode shapes of CCR for the first three eigenvalues that are shown and compared with numerical results in Fig. 6(b) (all numerical calculations were conducted by means of numerical software COMSOL Multiphysics). As can be seen from the comparison in Fig. 6(b), the differences between analytical and numerical results are indistinguishable. The radius function (43) is shown in Fig. 7(a). The corresponding polynomial G(s) can be approximated as G(s) = −
1 d2 R(s) u W (s) = −16.1372s4 +34.5048s3 −23.4950s2 +0.4608s+2 . R(s) ds2 (45)
19
(a)
(b) 1
1.6 U (s)
R(s)
0.5 1.4
0
1.2
-0.5
1
-1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
s
0.6
0.8
1
s
Figure 6: (a) Radius functions R(s) (42) for CCR; (b) The corresponding three normalized mode shapes for the first three eigenvalues: K = 3.4237 - black line; K = 6.4686 - blue line; K = 9.5555 - red line. The corresponding numerical results are represented by dots.
After substituting the coefficients of the quartic polynomial into the solution (21) we obtain the approximate analytical solution that can be used for calculation of eigenvalues and mode shapes as in the previous example, Fig. see 7(b). By comparison with numerical results, we can see that the differences are again indistinguishable.
(b)
1
1
0.9
0.5 U (s)
R(s)
(a)
0.8
0
-0.5 0.7
-1 0
0.2
0.4
0.6
0.8
1
0
s
0.2
0.4
0.6
0.8
1
s
Figure 7: (a) Radius functions R(s) (42) for CCR; (b) The corresponding three normalized mode shapes for the first three eigenvalues: K = 3.1763 - black line; K = 6.2972 - blue line; K = 9.4347 - red line. The corresponding numerical results are represented by dots.
20
7. Conclusion In this work, we have presented a closed-form analytical solution of the Webster equation for both elastic and acoustic waves propagating through waveguides with variable cross-sections. If a radius-function is given by the triconfluent Heun functions, then the presented analytical solution is exact. The use of triconfluent Heun functions allows to set the profile by seven optional parameters that affect the radius function profile, both in shape and size over a relatively wide range. We can use the presented analytical solution for all radius-function profiles designed in this way. To facilitate the search for suitable parameters for determining the profile shapes of radius-functions, we have included polynomial forms of the triconfluent Heun functions that can be found under certain conditions specified in the Appendix. In the case of a predetermined radius-function, then after its approximation by a polynomial solution of the triconfluent Heun equation, we can also employ the derived analytical solution of the Webster equation but only as an approximate one. We have demonstrated the applicability and correctness of the derived analytical solutions in calculations of natural frequencies and mode shapes for representative radius-functions. The approximate analytical solutions have been compared with numerical results to show that acceptable accuracy is achieved. The presented analytical solution of the Webster equation extends the group of already known exact solutions, and this solution is applicable to a wide range of radius-functions. In addition, the derived solution can be relatively straightforwardly extended even for the case of longitudinal viscoelastic waves or for the case of thermo-viscous fluids. Acknowledgement This work was supported by the the Grant Agency of the Czech Republic (GACR) grant No. 18-24954S.
21
Appendix A. The triconfluent Heun equation and its solution Here, we outline some facts concerning the triconfluent Heun equation and its solutions. The Heun equation is known to be a standard form of the second-order Fuchsian differential equation with four regular singularities {0, 1, a, ∞}. The canonical form of the general Heun equation is (see, e.g., [22, 23, 24]): γ d2 ϕ δ dϕ αβσ − f + + + + ϕ(σ) = 0 , dσ 2 σ σ − 1 σ − a dσ σ(σ − 1)(σ − a)
(A.1)
where a 6= 0, 1 is generally a complex parameter and the complex-valued exponent parameters α, β, γ, δ, and satisfy the (Fuchsian) condition (see, e.g., [22, 24]): 1+α+β =γ+δ+,
(A.2)
and f is the accessory parameter (it does not affect the exponent parameters). The Heun equation is invariant under α ←→ β.
Appendix A.1. General solution of the triconfluent Heun equation The triconfluent Heun equation is derived from the Heun equation by the coalescence of the three finite regular singular points with infinity ([24]). The canonical triconfluent Heun equation has the form (see, e.g., [24],[23]): dϕ(σ) d2 ϕ(σ) − (3σ 2 + γ) + [α + (β − 3)σ]ϕ(σ) = 0 , 2 dσ dσ and its symmetric (self-adjoint) used form is: d2 Φ(σ) γ2 3 2 9 4 + α− + βσ − γσ − σ Φ(σ) = 0 , dσ 2 4 2 4 where
σ 3 + γσ Φ(σ) = exp − 2
ϕ(σ) .
(A.3)
(A.4)
(A.5)
The quartic polynomial in Eq. (A.4) does not contain σ 3 term, and the coefficient of σ 4 is set equal to 9/4. This normalization can be arranged by subjecting the independent variable σ to an innocuous transformation of the form 22
σ 0 = d1 σ + d0 (d0 and d1 are constants), but other normalizations could equally well be used. The normalization herein is commonly used, see e.g., [24]. An analytical solution ϕ(σ) of Eq. (A.3) can be expressed using a power series of the form ϕ(σ) =
∞ X
(A.6)
bn σ n .
n=0
The derivatives of ϕ with respect to σ are: ∞ ∞ X dϕ(σ) X = nbn σ n−1 = (n + 1)bn+1 σ n , dσ n=0 n=0 ∞ ∞ X d2 ϕ(σ) X n−2 = n(n − 1)bn σ = (n + 2)(n + 1)bn+2 σ n . dσ 2 n=0 n=0
(A.7)
(A.8)
By substituting the power series (A.6), (A.7), and (A.8) into Eq. (A.3), we obtain the recurrence relation: n(n − 1)bn − (n − 1)γbn−1 + αbn−2 + (β + 6 − 3n)bn−3 = 0 .
(A.9)
Shifting the summation index leads to: (n + 2)(n + 3)bn+3 − (n + 2)γbn+2 + αbn+1 + (β − 3 − 3n)bn = 0 , (A.10) where b−1 = 0 , b0 = 1 , b1 = 0 , b2 = −
α . 2
(A.11)
From the recurrence relation (A.9), we can write: bn =
γ(n − 1)bn−1 − αbn−2 − (β + 6 − 3n)bn−3 , n(n − 1)
n≥3.
(A.12)
Additional information regarding this recursive relation is provided in the article [25] where the authors deal with an investigation of the asymptotics of a linearly independent set of solutions for a four-term recurrence relation controlling the coefficients appearing in the power series expansion for the triconfluent Heun functions.
23
For (α, β, γ) ∈ C, the triconfluent Heun function represents a solution of the triconfluent Heun equation (A.3) and can be written as THF(α, β, γ; σ) =
∞ X
bn σ n ;
(A.13)
|σ| < 1 ,
n=0
where bn are given by (A.11) and (A.12). The derivative of the triconfluent Heun function (the prime triconfluent Heun function) can be expressed as ∞ X d THF(α, β, γ; σ) ≡ THF0 (α, β, γ; σ) = nbn σ n−1 ; dσ n=0
|σ| < 1 .
(A.14)
We often need the second derivative of the triconfluent function. For this purpose, we use Eq. (A.1) and we obtain: d2 THF(α, β, γ; σ) = (3σ 2 + γ)THF0 (α, β, γ; σ) dσ 2 − [α + (β − 3)σ]THF(α, β, γ; σ) . (A.15) The general solution of the triconfluent Heun equation (A.3) is (see e.g., [24]): ϕ(σ) = C1 THF(α, β, γ; σ) + C2 exp(σ 3 + γσ)THF(α, −β, γ; −σ) ,
(A.16)
where C1 and C2 are integration constants. There is the following identity: THF(α, β, 0; σ) = exp(σ 3 )THF(α, −β, 0; −σ) ,
(A.17)
which means that if γ = 0 then the presented solutions in Eq. (A.16) are not linearly independent. We can obtain the general solution of the symmetric form of the triconfluent equation (A.4) using the solution (A.16) and the transformation relation (A.5): σ 3 + γσ Φ(σ) = C1 exp − THF(α, β, γ; σ)+ 2 3 σ + γσ THF(α, −β, γ; −σ) . C2 exp 2 (A.18)
24
Thus, in general, the triconfluent Heun differential equation is, by definition, a second-order ordinary differential equation of the form d2 Φ + Π(σ)Φ(σ) = 0 , dσ 2
(A.19)
where the coefficient function Π(σ) is some quartic polynomial in the independent variable σ. In principle, the coefficient function Π(σ) could be considered to be a polynomial of a higher degree than four. Solutions Φ(σ) of Eq. (A.19) when the function Π(σ) is a quintic, sextic or higher polynomial, would be special functions of new types. Since the classification scheme has not yet been extended to the higher degrees of polynomials than four, and thus is not part of the mathematical software such as Maple, we have limited ourselves in this work to the quartic polynomials. Appendix A.2. Polynomial solution of the triconfluent Heun equation We can see that the evaluation of successive terms of a series solution (A.6) or (A.13) to the triconfluent Heun equation (A.3) is carried out by means of a recurrence relation (A.12) where the coefficient bn depends upon n, the previous values of br (r < n), and the three parameters (α, β, γ). Based on this recursive relation, it is possible to express a solution of the triconfluent equation by a finite polynomial if certain conditions are met. The conditions for the existence of polynomial solutions can be formulated by the following theorem (see e.g., [24, 25]): Theorem 1. Suppose that in the triconfluent Heun equation (A.3) the parameters α, β and γ satisfy the following two conditions: (a) β = 3(N + 1),
N = 0, 1, 2 . . . ,
(b) The determinant DN +1 (α, γ) = 0 of the matrix
25
MN +1
α
3N 0 0 = .. . 0
−γ
2·1
0
0
···
α
−2γ
2·3
0
···
3(N − 1)
α
−3γ
3·4
···
0 .. .
3(N − 2) .. .
α .. .
−4γ .. .
··· ..
0
···
0 0 . .. . −(N − 1)γ N (N − 1) α −N γ 3·1 α (A.20) 0
. 3·3
α
0
3·2
0
0
Then the triconfluent Heun equation (A.3) has a polynomial (Liouvillian) solution ϕ(σ) = P β −1 (σ) , 3
where P β −1 (σ) denotes a polynomial of degree 3
β 3
− 1, whose coefficients pn
(n = 0, 1, . . . , β3 − 1) are the solutions of the following linear system of equations
p0
Mβ 3
p1 .. .
=0.
p β −1
(A.21)
3
If ϕ1 (σ) is the polynomial solution of the triconfluent Heun equation (A.3) then a second independent solution ϕ2 (σ) can be calculated by the Wronskian method (see e.g., [26]) Z
σ
ϕ2 (σ) = ϕ1 (σ)
exp
R w
(3v 2 + γ)dv dw . ϕ1 (w)2
(A.22)
References [1] A. G. Webster, Acoustical impedance, and the theory of horns and of the phonograph, Proceedings of the National Academy of Sciences 5 (1919) 275–282. doi:10.2307/84052. [2] L. L. Beranek, T. J. Mellow, Acoustics: Sound Fields and Transducers, Academic Press, 2012.
26
[3] A. Chaigne, J. Kergomard, Acoustics of Musical Instruments, SpringerVerlag, 2016. [4] M. V. Guiman, I. C. Rosca, A new approach on vibrating horns design, Shock and Vibration 8532021 (2017) 1–12. doi:doi.org/10.1155/2017/ 8532021. [5] K. Hornisova, P. Billik, Some properties of horn equation model of ultrasonic system vibration and of transfer matrix and equivalent circuit methods of its solution, Ultrasonics 54 (2014) 330 – 342. doi:https: //doi.org/10.1016/j.ultras.2013.05.003. [6] V. S. Sorokin, Effects of corrugation shape on frequency band-gaps for longitudinal wave motion in a periodic elastic layer, The Journal of the Acoustical Society of America 139 (2016) 1898–1908. doi:10.1121/1.4945988. [7] P. Honzík, S. Durand, N. Joly, M. Bruneau, On the acoustic transfer function of slowly tapered small horns filled with thermo-viscous fluid, Acta Acustica united with Acustica 99 (2013) 694–702.
doi:10.3813/aaa.
918648. [8] Y. A. Ilinskii, Nonlinear standing waves in an acoustical resonator, The Journal of the Acoustical Society of America 104 (1998) 2664–2674. doi: 10.1121/1.423850. [9] O. Rudenko, A. Shvartsburg, Nonlinear and linear wave phenomena in narrow pipes, Acoustical Physics 56 (4) (2010) 429–434. doi:10.1134/ S1063771010040044. [10] M. Cervenka, M. Soltes, M. Bednarik, Optimal shaping of acoustic resonators for the generation of high-amplitude standing waves, The Journal of the Acoustical Society of America 136 (2014) 1003–1012. doi: 10.1121/1.4892751.
27
[11] M. Cervenka, M. Bednarik, Acoustic bandpass filters employing shaped resonators, Journal of Sound and Vibration (2016) 76–88doi:10.1016/j. jsv.2016.06.045. [12] E. Eisner, Complete solutions of the "Webster" horn equation, The Journal of the Acoustical Society of America 41 (1967) 1126–1146. doi:10.1121/ 1.1910444. [13] S. Abrate, Vibration of non-uniform rods and beams, Journal of Sound and Vibration 185 (4) (1995) 703 – 716. doi:doi.org/10.1006/jsvi.1995. 0410. [14] B. Kumar, R. Sujith, Exact solutions for the longitudinal vibration of nonuniform rods, Journal of Sound and Vibration 207 (1997) 721 – 729. doi: doi.org/10.1006/jsvi.1997.1146. [15] Q. Li, Exact solutions for free longitudinal vibration of stepped non-uniform rods, Applied Acoustics 60 (1) (2000) 13 – 28. doi:doi.org/10.1016/ S0003-682X(99)00048-1. [16] A. Raj, R. Sujith, Closed-form solutions for the free longitudinal vibration of inhomogeneous rods, Journal of Sound and Vibration 283 (3) (2005) 1015 – 1030. doi:https://doi.org/10.1016/j.jsv.2004.06.003. [17] B. Yardimoglu, L. Aydin, Exact longitudinal vibration characteristics of rods with variable cross-sections, Shock and Vibration 18 (2011) 555–562. doi:10.3233/SAV-2010-0561. [18] B. A. Donskoy, Dimitri M.; Cray, Acoustic particle velocity horns, The Journal of the Acoustical Society of America 131 (2012) 3883–3890. doi: 10.1121/1.3702432. [19] A. D. Pierce, Acoustics: an Introduction to its Physical Principles and Applications, McGraw-Hill Book Company, Inc., 1981.
28
[20] S. W. Rienstra,
Webster’s horn equation revisited,
nal on Applied Mathematics 65 (2005) 1981–2004.
SIAM Jourdoi:10.1137/
s0036139902413040. [21] P. A. Martin, On Webster’s horn equation and some generalizations, The Journal of the Acoustical Society of America 116 (2004) 1381. doi:10. 1121/1.1775272. [22] G. Kristensson, Second Order Differential Equations:
Special Func-
tions and Their Classification, Springer New York, 2010. doi:10.1007/ 978-1-4419-7020-6. [23] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, Handbook of Mathematical Functions, Cambridge University Press, 2010. [24] A. Ronveaux (Ed.), Heun’s Differential Equations, Oxford Science Publications, 1995. [25] D. N. M. Batic, D.; Mills-Howell, Potentials of the Heun class: The triconfluent case, Journal of Mathematical Physics 56 (2015) 052106. doi: 10.1063/1.4921344. [26] K. Riley, M. Hobson, S. J. Bence, Mathematical Methods for Physics and Engineering, Cambridge University Press, 2006.
29
Declaration of interests ⊠ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: