Design of multi-directional acoustic cloaks using two-dimensional shape
optimization and the boundary element method
Peter RISBY ANDERSEN(1) , Vicente CUTANDA HENRÍQUEZ(2) , Lorenzo SANCHIS(3) , José SÁNCHEZ-DEHESA(4)(1) Centre for Acoustic-Mechanical Micro Systems, Technical University of Denmark, Denmark, prand@elektro.dtu.dk(2) Centre for Acoustic-Mechanical Micro Systems, Technical University of Denmark, Denmark, vcuhe@elektro.dtu.dk(3) Wave Phenomena Group, Dep. Electronic Engineering, Universitat Politècnica de València, Spain, lorenzo.sanchis@uv.es(4) Wave Phenomena Group, Dep. Electronic Engineering, Universitat Politècnica de València, Spain, jsdehesa@upv.esAbstract
Acoustic cloaking is a technique that seeks hiding objects in a sound field by reducing or cancelling theirscattered sound pressure. Incident waves are restored to as close as possible their original undisturbed form afterhitting the cloaked object. One technique for achieving this goal is the design of additional scatterers around theobject, which, properly distributed and shaped, can create cloaking at the design frequency. A newly developed
numerical technique combining the Boundary Element Method (BEM) with shape optimization is applied in thiswork for two-dimensional cloaking of a cylinder. The shapes of the scatterers are optimized for the cloakingof the whole setup with waves impinging on the cylinder from several different directions. The results show ameasure of the amount of cloaking depending on the direction at the range around the design frequency. The
optimization results are compared with existing one-directional cloaks. The impact of visco-thermal acoustic
losses in the cloaking design is also evaluated by means of a BEM implementation with losses.
Keywords: Acoustic cloaking, Boundary Element Method, Shape Optimization

1

INTRODUCTION

Designing acoustic cloaking devices by scattering cancellation posses a challenge. Especially, if the requirement
is to offer a cloak that can operate with incident waves at different angles, i.e., a multi-directional cloak. In
recent years, several examples of cloaking devices that are based on scattering cancellation have been proposedfor one-directional cloaking [1, 2, 3], but also examples of multi-directional acoustic cloaks exist [4]. However,obtaining multi-directional acoustic cloaks with similar performance to its one-directional counterpart is a nontrivial task. In the Refs. [1, 2, 4], the cloaks are based on design methods that utilize simple shapes restricted
to cylinders or Beziér scatters. Therefore, the design space is limited, meaning that finding satisfactory solutions can be cumbersome. In this paper, we will apply a newly developed shape optimization technique to theoptimization of a multi-directional cloak. The shape optimization technique is applicable to both isentropic andlossy optimization and relies on the Boundary Element Method (BEM), which will allow more freedom for thecreation of the scatter shapes [5, 6].
In the following, firstly the BEM for scatting problems is introduced, hereafter, the initial design proposal
and parametrization of the multi-directional cloak is shown. This is followed by a section on the optimization
problem and its implementation. In the end, the final shape-optimized design is presented and its performanceis validated through several benchmarks, including a comparison with existing one-directional acoustic cloaking
designs.

2BOUNDARY ELEMENT METHOD FOR SCATTERING PROBLEMS

The shape optimization technique is based on the BEM by direct collocation. The formulation of the BEM
takes its starting point at the Kirchhoff-Helmholtz integral equation, given by5600

∂ p(Q)∂ G(R)
p(Q)dΓ − G(R)dΓ + pI
(1)
∂n
Γ
Γ ∂n
where p is the acoustic pressure, C is the integral free term, P is a collocation point, Q are integration points
on the generator Γ, G(R) is the two dimensional fundamental solution with R = |P − Q| and pI is an incident
pressure. Discretization using collocation and assuming all surfaces to be rigid, i.e. ∂ p/∂ n = 0, Eq. (1) can be
written in matrix form asZZC(P)p =Ap + pI = 0.

(2)

The assembly of the A matrix is carried out using a MATLAB MEX/C++ compiled boundary element code.
the scattered pressure in the domain can be obtained by using field point calculations
ps = A p p,(3)where ps is the scattered pressure at different field points and A p is the field point matrix.

Control points

SymmetryΩdy [m]Plane waveθ
Objectx [m]Figure 1. The initial design guess of the multi-directional cloak. The cloak design is based on 3 symmetricallydescribed scatterers that are duplicated every 18 degrees around the cylindrical object, resulting in 20 pairs ofidentical scatterers. A detailed illustration of the parametrization is seen in the upper right corner, where theblue dots represent the location of the control points that forms the three scatterer pairs. Incident to the design
is a plane wave oscillating at 340 Hz with an incident angle of θ = 0. The radius of the cylindrical object is 1m.

5601

3INITIAL DESIGN AND PARAMETRIZATION

The initial design guess and the parametrization type are important factors for successful shape optimization.
In the case of multi-directional cloaking it is a non-trivial task to guess on a good initial design. The initial
guess used here is seen in Figure 1. This design is based on an initial optimization process using the analyticaloptimization approach found in Ref. [2], where the cloak is described analytically from Bézier scatters andoptimized using a combination of a genetic algorithm and simulated annealing [7]. It should be noted that theinitial design guess is already an improvement in terms of cloaking of the object.
To conduct the shape optimization the three scatterers are parametrized symmetrically using C2 -continuous
parabolic blending to describe the mesh between control points [8]. The locations of the control points aremarked by blue dots in Figure 1. The total number of control points is 45 with each control point having
the possibility to move in both the x- and y-direction, leading to a total of 90 design variables and a meshconsisting of 5104 quadratic continuous elements. The design of the 3 scatterers is copied for every 18 degrees
around the object, meaning that, due to symmetry, the design has 20 possible directions with the same cloaking
performance.

4OPTIMIZATION PROBLEM

The goal is to minimize the scattered pressure in the domain Ωd (the blue region defined in Figure 1). This isaccomplished by solving the optimization problem stated asZ

min: φ =v

s. t.

|ps |2 dΩΩdκ(v) − κmax ≤ 0(4)Dmin − β D(v) ≤ 0.
Eq. (2),where φ is the objective function, v is a vector containing all the design variables, κ is the curvature at eachcontrol point, κmax is the maximum curvature allowed, D(v) is the distance between control points, Dmin is aminimum allowed distance and β is a self-intersection parameter. A more elaborate description of the curvatureand distance constraints can be found in Refs. [5] and [6]. In the following results, κmax = 250 m−1 andDmin = 1.3 cm. The optimization problem is solved using the gradient-based sequential quadratic programmingalgorithm in the MATLAB function "fmincon".
4.1 Sensitivity analysis
The gradients with respect to the design variables v are calculated using a semi-analytical discrete adjoint approach. Therefore, a gradient with respect to a single design variable v can be obtained by [9] ∂φ∂A
dφ=+ℜ λT+ pI,dv
∂v
∂v

(5)

here λ are the Lagrange multipliers, which can be found by solving the system∂φ∂φ−iA λ =−∂ pr∂ piT

T,

(6)

where the subscripts r and i denote the real and imaginary parts of the complex pressure. In general, the matrix√
derivatives are obtained by finite difference at a matrix level using a step length equal to eps. The sensitivity
analysis is performed in parallel using a combination of the MATLAB "parfor" loop and OpenMP.

5602b)a)
Symmetryx [m]Initialdesignφy [m]

Control points

Object without cloak

Optimized
designIteration number

Figure 2. a) Shape-optimized design of the multi-directional cloak, including a magnified area presenting thelocation of the control points. b) evolution of the objective function φ during optimization (continuous blueline). Additionally, the plot includes the value of the objective function when only the cylindrical object isconsidered (dashed red line).

5RESULTS

In this section, the results of the shape optimization of the multi-directional cloak are presented. This includes
several performance measures studying the frequency dependence, the angle of the incident wave, the effect ofviscothermal losses and the visibility of the cloak compared to two existing one-directional cloaking designs.
The full optimization time is approximately 3 days and is is performed on a 6 core Xeon E5-1650v3 processor
with 128 GB of memory. Figure 2a shows the shapes of the 3 scatterers for the optimized design. The design
is very similar to the initial guess, indicating that the initial design is very important for the shape optimization
outcome. In Figure 2b, the evolution of the objective function for each iteration is shown. Aditionally, thefigure also includes the objective function when only the cylindrical object is considered. The optimized design
is an improvement over the initial guess, with a value of the objective function around φ = 4.9, lower than in
the initial design (φ = 6.8). However, the design does not achieve an objective function close to zero, which
would correspond to perfect cloaking and complete reconstruction of the incident wave.
5.1 Cloak performance as a function of the incident angle and frequency
The objective function with respect to the frequency of the incident wave is plotted in Figure 3a. As expected,
the optimized design reduces the objective function at the optimization frequency. Additionally, the optimized
design shows a better "broadband" performance as compared to the initial design. It should be noted, that thisbroadening of the response is obtained using a single optimization frequency and it is not included as an optimization criterion.
Moreover, Figure 3b shows how the objective function of the initial and optimized designs depended on theincident angle of the impinging plane wave. The incident angle has very little influence on the cloaking performance. Figure 3a and Figure 3b also includes comparative Finite Element Method (FEM) simulations of theoptimized design. The FEM simulations are performed using COMSOL Multiphysics.

5603a)b)φφOptimization
frequency

Frequency [Hz]

θ [degrees]

Figure 3. a) Frequency dependence of the objective function for both the initial (blue) and the optimized (red)design. The figure also includes simulations including viscous and thermal losses (green). b) Dependence of
the incident angle of the impinging plane wave for the initial (blue) and the optimized (red) design. circularsymbols represent equivalent FEM simulations for comparison.
5.2 Effect of viscous and thermal dissipation
The effect of viscous and thermal dissipation on acoustic cloaks based on scattering cancellation was previously
studied in Ref. [10]. In this reference, it is concluded that the effect of viscothermal losses can be neglected
when performing simulations on cloaking devices. However, since the design proposed here is fundamentallydifferent from the designs found in Refs. [1, 2, 6], with much more compactly arranged scatterers, a study ofthe effect of losses on the objective function seems appropriate. Therefore, the frequency dependent behaviorof the objective function is also plotted in Figure 3a for both full viscothermal BEM and FEM simulations(green curves). The lossy BEM simulations are based on the Kirchhoff decomposition of the full linearized
Navier-Stokes (FLNS) equations and implemented following the approach in Ref. [5]. On the other hand, theFEM simulations are carried out using the FLNS implementation in COMSOL Multiphysics. It is observedthat the effect of losses on the objective function is small, and simulations of the objective function are very
similar to the isentropic computations. However, small discrepancies between the two lossy BEM and FEM
simulations can be observed. A possible explanation for the discrepancies might be due to differences in theirmesh densities. The lossy FEM simulations are performed using a mesh with approximately 8 million degrees offreedom. On the contrary, the BEM calculations are based on the same mesh as the isentropic BEM calculations.
Similar small discrepancies between viscothermal FEM and BEM calculations are found in Ref. [10].
5.3 Visibility comparison with existing one-directional cloaks
A common measure to validate the cloaking performance is the so-called average visibility. The definition of
the average visibility and how to calculate it can be found in, e.g., the Refs. [1, 10]. In Figure 4, the averagevisibility of the shape-optimized multi-directional cloak is plotted and compared to two one-directional cloaks
that can be found in Refs. [1] and [6]. The two one-directional cloaks are either created using a distribution of
120 cylinders or using shape optimized scatterers based on the same approach used in this paper. The cloaks
are designed with different operational frequencies in mind. Therefore, the frequency behavior is normalizedwith respect to their individual operational frequency f0 . As expected, the visibility of the multi-directional

5604

Average visibility

cloak is not as low as for the cloaks optimized for a single direction.

f / f0Figure 4. The visibility of the optimized multi-directional cloak compared to two one-dimensional cloaks foundin the literature. The black curve is the visibility of the cloak found in Ref. [1] which consists of 120 cylinders,
and the red curve is the visibility of the shape-optimized one-dimensional cloak found in Ref. [6]. The green
curve represents the visibility of the shape optimized multi-directional cloak as seen in Figure 2a. Here, thefrequency f is normalized with respect to the operational frequency of the different cloaks. The operationalfrequency is denoted by f0 .

6CONCLUSIONSIn this work, a newly developed shape optimization technique has been utilized to improve the design of amulti-directional cloak based on Bézier scatterers. Although the shape-optimized multi-directional cloak is out
performed by the one-directional cloaks at the design frequency, it offer a broader frequency range of operation,
also compared with the initial omni-directional design. Interestingly, the proposed design is very robust withrespect to the angle of the plane wave (Figure 3b). In fact, the cloaking performance does not depend on theincident angle making the design truly omni-directional.
The shape optimization results presented here should be considered an initial step towards creating high-performing
multi-directional acoustic cloaking devices. Several measures can be taken to improve on the optimized designby, e.g., adding more scatterers or altering the initial location of the scatters, and perhaps including the locationof the scatter as additional design variables. Moreover, optimization studies can be performed with less restrictive constraints on, e.g., the curvature κ. This would allow for more design freedom perhaps improving thecloaking performance. The negligible angle dependence of the optimized cloak suggests that the periodic scat-5605

terers could be made broader and fewer. It can also be concluded that, even though the scatterers of the initialand optimized designs are compactly arranged, the dissipative effects can safely be neglected in the optimization
procedure.

ACKNOWLEDGEMENTS
J. Sánchez-Dehesa acknowledges the support by the Spanish Ministerio de Economıá y Competitividad, and theEuropean Union Fondo Europeo de Desarrollo Regional (FEDER) through Project No. TEC2014-53088-C3-1-R.

REFERENCES[1] García-Chocano,V. M.; Sanchis L.; Díaz-Rubio, A.; Martínez-Pastor, J.; Cervena, F.; Llopis-Pontiveros, R.;
Sánchez-Dehesa, J. Acoustic cloak for airborne sound by inverse design. Applied Physics Letters, 99 (7),2011, 074102[2] Lu, Z.; Sanchis L.; Wen J.; Cai L.; Bi,Y.; Sánchez-Dehesa, J. Acoustic cloak based on Bézier scatterers.
Scientific Reports, 8 (1), 2018, 12924[3] Andkjær, J.; Sigmund, O. Topology optimized cloak for airborne sound. ASME J. Vib. Acoust., 135, 2013,041011.
[4] Jo, C.; Jeong, J.; Kwon B.; Park, K.; Oh, I. Omnidirectional two-dimensional acoustic cloak by axisymmetric cylindrical lattices. Wave Motion, 54 , 2015, 157-169.
[5] Risby Andersen, P.; Cutanda Henríquez, V.; Aage, N. Shape optimization of micro-acoustic devices including viscous and thermal losses. Journal of Sound and Vibration, 447, 2019, 120-136.
[6] Risby Andersen, P. Modelling of acoustic viscothermal losses using the Boundary Element Method: Frommethod to optimization. PhD thesis, The Technical University of Denmark, Lyngby, Denmark, 2019[7] Sanchis, L., Cryan, M. J., Pozo, J., Craddock, I. J.; Rarity, J. G. Ultra high Purcell factor in photonic
crystal slab microcavities. Phys. Rev. B, 76 , 2007, 157-169.
[8] Johnston, P. R. C2-continous elements for boundary element analysis. International Journal for NumericalMethods in Engineering, 40 , 1997, 2087-2108.
[9] Dühring, M. B.; Jensen, J. S.; Sigmund, O. Acoustic design by topology optimization. Journal of Soundand Vibration, 317 , 2008, 557-575.
[10] Risby Andersen, P.; Cutanda Henríquez, V.; Aage, N. Sánchez-Dehesa, J.; Visco-thermal effects on an
acoustic cloak based on scattering cancellation, Proceedings of NOise and Vibration Emerging Methods,Ibiza, Spain, 20185606