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.

## 1. Introduction

Mironov^{1} 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)=\u03f5xm$ for $m\u22652$, and others^{2} 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)=\u03f5xm+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.

## 2. Theoretical model

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

where $q(x,t)$ is the external force applied to the beam, $E$ is the Young's modulus of the beam material, $\mu $ 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\omega t$, one form of the solution to Eq. (1) is

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

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

The $\Pi $ 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=Zj\u22121uj\u22121+q$, where $q=[000q0]T$. The two ends of the beam are then related

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

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\theta i]T$, $dif=[q0\delta i(j\u22121)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\xd72$ matrix and $pi$ is a $2\xd71$ vector.^{12} The corresponding equation for $ei$ is $ei=Tiei+1+gi$, where $Ti$ is similarly a $2\xd72$ matrix and $gi$ is a $2\xd71$ vector. Given these two equations and the modified transfer matrix equation, the following recursive relations can be deduced:

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\xd71$ vectors are converted to $2K\xd71$ vectors and $2\xd72$ matrices are converted to $2K\xd72K$ block-diagonal matrices. For example, $fi$ becomes $[fiT(\omega 1)fiT(\omega 2)\cdots fiT(\omega K)]T$, while $Si$ becomes

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)\u2264106$.

## 3. Objective function

The objective function chosen for this problem is

where $W(x,\omega ,h)$ is the solution to the steady-state equation

That is, $W$ is the (complex) displacement amplitude of a beam under a harmonic point force, $|j\omega W|2$ is the square magnitude of its velocity, and $\u27e8W\u03072\u27e9$ is the spatially averaged square velocity. Because $\mu $ 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 $[\omega 1,\omega 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$,

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

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

## 4. Parameters

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}

Parameter . | . | Material Property . | . |
---|---|---|---|

$L$ | 30 cm | $E$ | 70 GPa |

$[\omega 1,\omega 2]$ | $2\pi \u22c5[50,\u20092000]$ rad/s | $\rho $ | 2700 kg/m^{3} |

$h1$ | 6.35 mm | $\eta $ | 0.0001 |

$[LABH\u2212,LABH+]$ | [0, 22.26] cm | $Edamp$ | 9 MPa |

$[h0\u2212,h0+]$ | [0.635, 6.35] mm | $\rho damp$ | 1812 kg/m^{3} |

$[m\u2212,m+]$ | [1, 12] | $\eta damp$ | 0.2 |

$q0$ | 1 N | ||

$xq$ | 22.26 cm |

Parameter . | . | Material Property . | . |
---|---|---|---|

$L$ | 30 cm | $E$ | 70 GPa |

$[\omega 1,\omega 2]$ | $2\pi \u22c5[50,\u20092000]$ rad/s | $\rho $ | 2700 kg/m^{3} |

$h1$ | 6.35 mm | $\eta $ | 0.0001 |

$[LABH\u2212,LABH+]$ | [0, 22.26] cm | $Edamp$ | 9 MPa |

$[h0\u2212,h0+]$ | [0.635, 6.35] mm | $\rho damp$ | 1812 kg/m^{3} |

$[m\u2212,m+]$ | [1, 12] | $\eta damp$ | 0.2 |

$q0$ | 1 N | ||

$xq$ | 22.26 cm |

## 5. Results

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\omega W|2$, is calculated and integrated across $x$ and $\omega $ 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 algorithm^{14} 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 $[\omega 1,\omega 2]$.

. | $LABH$ . | $h0$ . | $m$ . | $J(h)$ . |
---|---|---|---|---|

Unmodified beam | 0 cm | 6.35 mm | 1 | 1.000 |

Optimal design | 22.26 cm | 0.635 mm | 3.06 | 0.2987 |

. | $LABH$ . | $h0$ . | $m$ . | $J(h)$ . |
---|---|---|---|---|

Unmodified beam | 0 cm | 6.35 mm | 1 | 1.000 |

Optimal design | 22.26 cm | 0.635 mm | 3.06 | 0.2987 |

Looking at Fig. 2, which shows the optimal design's spatially averaged square velocity response, $\u27e8W\u03072\u27e9$, 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.