Structures whose thickness follow a power law profile exhibit the “acoustic black hole” (ABH) effect and can be used for effective vibration reduction. However, it is difficult to know a priori what constitutes the best design. A new block matrix formulation of the transfer matrix method is developed for use in the optimization of an ABH vibration absorber at the end of a cantilever beam. Results indicate that introduction of the ABH significantly alters the dynamics of the beam, which must be considered in determining the optimal design for a given vibration reduction problem.

Mironov1 was the first to demonstrate the theoretical possibility of zero reflection of bending waves in a wedge whose thickness decreases to zero according to the relation h(x)=ϵxm for m2, and others2 have since shown this possibility for other thickness profiles. As the waves approach the zero-thickness end of such a so-called acoustic black hole (ABH), the wavespeed approaches zero and the transit time to the end approaches infinity. However, in practice there will always be a nonzero truncation at the tip of the wedge, so that the thickness profile will be h(x)=ϵxm+h0. The result of this modification is that the transit time to the end becomes finite, which can result in reflection of as much as 70% of the incident energy even for small values of h0.3 

However, the shortcomings of real wedges can be drastically improved with the addition of a thin layer of viscoelastic damping material near the tip. In such a case, the ABH works by reducing the bending wavelength and increasing the transverse amplitude, thereby focusing strain in a small region where damping material can more effectively absorb and dissipate energy. Such ABH vibration absorbers have received increased attention in recent years as a means to effectively reduce structural vibrations without adding weight.4–10 This paper explores the optimal design of a one-dimensional ABH vibration absorber located at the end of a cantilever beam, using a block transfer matrix approach.

The dynamic equation for an Euler-Bernoulli beam can be written as

μ2wt2+2x2(EI2wx2)q(x,t)=0,
(1)

where q(x,t) is the external force applied to the beam, E is the Young's modulus of the beam material, μ is its density per unit length, and I is the second moment of inertia about the beam's neutral axis. In the case that E and I are independent of position, there are no external forces, and assuming a steady-state, time-harmonic solution, w(x,t)=W(x)ejωt, one form of the solution to Eq. (1) is

W(x)=W0c1(x)+θ0kc2(x)M0EIk2c3(x)Q0EIk3c4(x),
(2)

where k4=μω2/EI, W0 is the vertical displacement of the beam's neutral axis at x=0, θ0 is its slope, M0 is the moment about the neutral axis, and Q0 is the shear force. Additionally, c1(x)=(cosh(kx)+cos(kx))/2, c2(x)=(sinh(kx)+sin(kx))/2, c3(x)=(cosh(kx)cos(kx))/2, and c4(x)=(sinh(kx)sin(kx))/2. For a finite beam of length L, it is possible to relate the state variables at x=0, u1=[W0θ0M0Q0]T, to the state variables at x=L, u2=[WLθLMLQL]T, through the transfer matrix

Z=[c1(L)c2(L)kc3(L)EIk2c4(L)EIk3kc4(L)c1(L)c2(L)EIkc3(L)EIk2EIk2c3(L)EIkc4(L)c1(L)c2(L)kEIk3c2(L)EIk2c3(L)kc4(L)c1(L)].
(3)

This relation may be written succinctly as u2=Zu1. This transfer matrix method can also be used to relate the state variables at any intermediate point along the beam, so long as each segment is uniform and homogenous. In this way, the dynamics of complex beam geometries can be analyzed by partitioning them into segments that are approximately uniform and homogenous. For a beam partitioned into N segments with N+1 nodes, the two ends of the beam can be related

uN+1=(i=1NZi)u1=Au1.
(4)

The Π notation is used in Eq. (4) and subsequently to denote sequential left multiplication of matrices. Given a set of known boundary conditions, the unknown state variables in u1 and uN+1 can be solved by explicitly calculating A and carrying out the appropriate algebra. State variables at intermediate points can then be determined using sequential application of transfer matrices.

If the beam is excited by a harmonic point force, then the transfer matrix formulation is modified as follows. Suppose a force with magnitude q0 is applied at the node 1<j<N+1. The state vector at this node would then be calculated as uj=Zj1uj1+q, where q=[000q0]T. The two ends of the beam are then related

uN+1=(i=1NZi)u1+(i=jNZi)q=Au1+Bq
(5)

and the unknown state variables can be solved just as before.

For reasons of numerical stability, calculations were implemented using the Riccati transfer matrix method developed by Horner and Pilkey.11 The advantage of the Riccati transfer matrix method is that, like the generalized Riccati transformation, it converts a numerically unstable two-point boundary value problem into a stable initial value problem. Moreover, the formulation of Horner and Pilkey requires no numerical integration, and the matrix components can be determined analytically. The reader is referred to Ref. 11 for details of the derivation, but in short the state vector and transfer matrix are modified to be of the form

[fe]i+1=[Z11Z12Z21Z22]i[fe]i+[dfde]i,
(6)

where fi are the state variables that are homogeneous at x=0 and ei are the complimentary forces or displacements. For the cantilever beam with a harmonic point force at node j, fi=[QiMi]T, ei=[Wiθi]T, dif=[q0δi(j1)0]T, and die=0, with Zi defined appropriately. The generalized Riccati transformation relating fi to ei is given by fi=Siei+pi, where Si is a 2×2 matrix and pi is a 2×1 vector.12 The corresponding equation for ei is ei=Tiei+1+gi, where Ti is similarly a 2×2 matrix and gi is a 2×1 vector. Given these two equations and the modified transfer matrix equation, the following recursive relations can be deduced:

Si+1=(Zi11Si+Zi12)(Zi21Si+Zi22)1,
(7)
Ti=(Zi21Si+Zi22)1,
(8)
pi+1=(Zi11pi+dif)Si+1(Zi21pi+die),
(9)
gi+1=Ti(Zi21pi+die).
(10)

Equations (7)–(10) are solved sequentially from node i=1 to node i=N+1, the boundary conditions at x=L are applied, and finally fi and ei are solved in reverse from node i=N+1 to node i=1.

Because Zi is frequency-dependent, fi and ei would normally need to be calculated separately for each frequency of interest. For additional speed-up, a block matrix approach was developed to handle the calculation of multiple frequencies at once. For K analysis frequencies, 2×1 vectors are converted to 2K×1 vectors and 2×2 matrices are converted to 2K×2K block-diagonal matrices. For example, fi becomes [fiT(ω1)fiT(ω2)fiT(ωK)]T, while Si becomes

[Si(ω1)000Si(ω2)000Si(ωK)].
(11)

Analysis can then be carried out exactly as in the single-frequency case. Although this block approach can theoretically handle all analysis frequencies simultaneously, for reasons of memory limitations the size of the matrices were restricted in this study so that 4K2(N+1)106.

The objective function chosen for this problem is

J(h)=CLω1ω20L|jωW|2dxdω=Cω1ω2Ẇ2dω,
(12)

where W(x,ω,h) is the solution to the steady-state equation

μω2W+2x2(EI2Wx2)q0δ(xxq)=0.
(13)

That is, W is the (complex) displacement amplitude of a beam under a harmonic point force, |jωW|2 is the square magnitude of its velocity, and Ẇ2 is the spatially averaged square velocity. Because μ and I are both functions of h(x), W will also depend upon h(x). The normalization factor, C, is chosen such that J(h)=1 for an unmodified beam with thickness profile h(x)=h1. The objective function thus represents the percent reduction in the total, spatially averaged square velocity response across the frequency band [ω1,ω2], as compared to an unmodified, uniform beam. To limit the design space, the set of possible thickness profiles is defined by the set V,

V={hL(0,L)|h(x)={ϵxm+h00xLABHh1LABH<xL,ϵ=h1h0(LABH)m,0LABHLABHLABH+,0h0h0h0+,0<mmm+},
(14)

which defines the geometry of a beam with an ABH at one end, and includes the ad hoc bounds set upon the design variables used in the optimization procedure. These design variables are the length of the ABH taper, LABH, the minimum thickness of the taper, h0, and the power-law exponent of the taper, m. Having defined the objective function and the set of possible solutions, the optimization problem is formulated as

minhVJ(h) such that ω[ω1,ω2],
EI2Wx2|x=0=EI3Wx3|x=0=0,
(15)
W(L)=0,Wx|x=L=0,
(16)

where the four constraints represent the free-clamped boundary conditions of a cantilever.

Table 1 lists the parameters and material properties used in this optimization study. A free damping layer was added to the thinnest 25% of the ABH taper model, which was determined by the authors to be the best trade-off between additional mass and damping for the most extreme taper profile considered in this study. This added damping was incorporated into the model using a complex Young's modulus and additional loss factor.13 

A routine was written to calculate J(h) given the design variables LABH, h0, and m. The geometry is first divided into a minimum of 50 segments per wavelength, which was shown to produce a relative error less than 10−5 without significantly increasing calculation time. The displacement, W, due to the point force is calculated using the block Riccati transfer matrix method described in Sec. 2 with appropriate boundary conditions. Finally, the square velocity, |jωW|2, is calculated and integrated across x and ω using the trapezoidal rule, and multiplied by the appropriate factors to obtain J(h). Analysis frequencies were chosen to range from 50 to 2000 Hz with eleven points per 3 dB bandwidth, assuming a structural Q of 238.5. This routine was used with the Borg multi-objective evolutionary algorithm14 for 10000 function evaluations to determine the optimal design.

Table 2 shows the optimal design variables, as well as the value of J(h) for this design compared to the unmodified beam. A graphical depiction of the optimal taper profile is shown in Fig. 1. From the results in Table 2, it is clear that the ABH reduces the objective function by a factor of 3.35 within the frequency band considered. As predicted by classic ABH theory, a longer taper is optimal for vibration reduction. However, a larger taper power was not found to be optimal, which is contrary to classic ABH theory. This may be due to the frequency range used in the objective function, which would penalize a higher taper power if it resulted in a greater number of structural modes in the frequency range [ω1,ω2].

Looking at Fig. 2, which shows the optimal design's spatially averaged square velocity response, Ẇ2, as a function of frequency, there is clearly a greater number of resonances compared to an unmodified beam, although it should be noted that the average response is still lower for the optimal design. Indeed, classic ABH theory predicts that for increasing taper power, ABH modes should move closer together in frequency. The inclusion of an ABH vibration absorber therefore significantly alters the dynamics of the system. This may be detrimental if trying to control a discrete set of resonances without affecting other frequency regions, but is not necessarily so for broadband vibration reduction, since the overall structural losses also increase with the inclusion of an ABH vibration absorber. Future work will investigate the effect of frequency range on the optimal ABH design.

As the results of this study corroborate, ABH vibration absorbers are highly effective at reducing bending vibrations in structures. However, the inclusion of an ABH vibration absorber significantly alters the dynamics of the structure, and the frequency range of interest must be considered. The optimization of ABH vibration absorbers can be a valuable tool in practical applications, so long as the objective function is tailored to the particular application and design outcomes.

1.
M. A.
Mironov
, “
Propagation of a flexural wave in a plate whose thickness decreases smoothly to zero in a finite interval
,”
Akust. Zh.
34
(
3
),
318
319
(
1988
)
M. A.
Mironov
, [
Sov. Phys. Acoust.
34
(
3
),
318
319
(
1988
)].
2.
A.
Karlos
,
S. J.
Elliott
, and
J.
Cheer
, “
Higher-order WKB analysis of reflection from tapered elastic wedges
,”
J. Sound Vib.
449
,
368
388
(
2019
).
3.
V. V.
Krylov
, “
New type of vibration dampers utilising the effect of acoustic ‘black holes
,’”
Acta Acust. Acust.
90
(
5
),
830
837
(
2004
).
4.
V.
Denis
,
F.
Gautier
,
A.
Pelat
, and
J.
Poittevin
, “
Measurement and modelling of the reflection coefficient of an acoustic black hole termination
,”
J. Sound Vib.
349
,
67
79
(
2015
).
5.
M. R.
Shepherd
,
P. A.
Feurtado
, and
S. C.
Conlon
, “
Multi-objective optimization of acoustic black hole vibration absorbers
,”
J. Acoust. Soc. Am.
140
(
3
),
EL227
EL230
(
2016
).
6.
J. G.
Ih
,
M.
Kim
,
I. J.
Lee
, and
J. S.
Jensen
, “
Truncated acoustic black hole structure with the optimized tapering shape and damping coating
,” in
Proceedings of the 45th International Congress and Exposition on Noise Control Engineering
(
2016
).
7.
L.
Tang
and
L.
Cheng
, “
Enhanced acoustic black hole effect in beams with a modified thickness profile and extended platform
,”
J. Sound Vib.
391
,
116
126
(
2017
).
8.
C. A.
McCormick
and
M. R.
Shepherd
, “
Optimal design and position of an embedded one-dimensional acoustic black hole
,” in
Proceedings of the 47th International Congress and Exposition on Noise Control Engineering
(
2018
).
9.
X.
Li
and
Q.
Ding
, “
Analysis on vibration energy concentration of the one-dimensional wedge-shaped acoustic black hole structure
,”
J. Intel. Mat. Syst. Str.
29
(
10
),
2137
2148
(
2018
).
10.
W.
Huang
,
H.
Zhang
,
D. J.
Inman
,
J.
Qiu
,
C. E.
Cesnik
, and
H.
Ji
, “
Low reflection effect by 3D printed functionally graded acoustic black holes
,”
J. Sound Vib.
450
,
96
108
(
2019
).
11.
G. C
Horner
and
W. D.
Pilkey
, “
The Riccati transfer matrix method
,”
J. Mech. Design
100
(
2
),
297
302
(
1978
).
12.
G. B.
Rybicki
and
P. D.
Usher
, “
The generalized Riccati transformation as a simple alternative to invariant imbedding
,”
Astrophys. J.
146
,
871
879
(
1966
).
13.
L.
Cremer
and
M.
Heckl
,
Structure-Borne Sound
, 2nd ed., translated by
E. E.
Ungar
(
Springer-Verlag Berlin
,
Heidelberg
,
1988
), Chap. 3.
14.
D.
Hadka
and
P.
Reed
, “
Borg: An auto-adaptive many-objective evolutionary computing framework
,”
Evol. Comput.
21
(
2
),
231
259
(
2013
).