Optimizing the Shape of a Horn

This example shows how to apply boundary shape optimization to a simple axisymmetric horn. For the sake of simplicity, the far-field sound pressure level is maximized for a single frequency and in a single direction. The focus is on the optimization procedure, which involves choice of objective function, and optimization solver settings. The deformation of the geometry is introduced using the Shape Optimization functionality of COMSOL Multiphysics.

The model was inspired by the work of Erik Bängtsson, Daniel Noreland, and Martin Berggren (Ref. 1).

This application requires the Acoustics Module and the Optimization Module.

Figure 1: The initial configuration is a simple cone (the z<0 part of the gray area) in an infinite baffle.

Model Definition

A plane-wave mode feeds an axisymmetric horn radiating from an infinite baffle toward an open half space; see Figure 1. The radius of the feeding waveguide is assumed to be fixed, as well as the depth of the horn and the size of the hole where the horn is attached to the baffle. By varying the curvature of the initially conical surface of the horn, its directivity and impedance can be changed.

The surface is parameterized assuming that the radius and the z position of the horn deviates from the simple cone by a set of 8th-order Bernstein polynomials. The maximum displacement (in the two directions) is given by the dmax parameter. The number of optimization variables is determined by the order of the polynomial. Using a higher order gives more freedom and potentially a better final value of the objective function, but it also makes the optimization process more sensitive and can generate a shape that is less suitable for production.

Optimization can only be applied to real-valued functions because the minimum of a complex-valued function is not well-defined. But the raw result from a frequency-domain acoustics simulation is a complex-valued pressure field. From this you generate a scalar, real-valued quantity to be used as the objective function in the optimization process. However, any operation which converts a complex number to a real value is necessarily non-analytical, which means that its derivative is not uniquely defined.

The gradient-based optimization solver in the Optimization Module by default evaluates derivatives of the objective function via the solution of an adjoint equation. This procedure requires that the symbolic derivative of any non-analytic function is selected in a special way. The default behavior of the composite functions abs(z) and conj(z), which are most commonly used to obtain a real-valued objective function, is to return a derivative parallel to the real axis. However, this behavior is not appropriate for the adjoint method, where you instead need the definitions

(1)

It is indeed possible to redefine the symbolic derivatives of built-in functions in COMSOL Multiphysics, but in this case it is more convenient to use the special function realdot(z1, z2), which evaluates as real(z1·conj(z2)) but differentiates according to Equation 1. In particular, as a measure of the transmission properties of the horn, use an expression of the form realdot(pm, pm)/p02, where pm is the pressure measured at a specific point in front of the horn and p0 is the (real-valued and constant) amplitude of the incoming wave.

If you choose to evaluate pm in the near-field, or can afford to include a sufficiently large domain in front of the horn to effectively measure a far-field value at a point in the model, you can simply measure pm as the local pressure in a geometry vertex. However, in order to optimize the directivity pattern in an efficient way, pm should be defined using an integral representation of the exterior-field pressure as a function of the angle from the axis. Typically, for speaker it will be the pressure or SPL evaluated at a 1 m distance.

COMSOL Multiphysics contains optimized code for evaluating such exterior field integrals. This is, however, a pure postprocessing feature that does not support the automatic differentiation required by the adjoint method. Therefore, return to the definition of the Helmholtz-Kirchhoff integral as given in its axisymmetric form. See the Theory Background for the Pressure Acoustics Module section in the Acoustics Module User’s Guide for further details. To simplify things we will use the expression of the exterior field in the far-field limit.

(2)

If the infinite baffle is placed at z = 0, its effect is the same as that of adding a mirror image of the horn and at the same time removing the baffle. If, in addition, the integration surface is taken to be the wide end of the horn, in the plane of the baffle, most of the terms in Equation 2 cancel out, and all that is left is

(3)

where J0 is the Bessel function of the first kind of order 0,

is the far-field evaluation distance, and the angle θ from the axis has been introduced as a parameter. This integral is easily implemented in COMSOL Multiphysics using an integration operator.

In this way the metric for the optimization is based on the far-field limit of the exterior field, which is valid for distances larger than the Rayleigh radius given by

Optimization as a rule implies many evaluations of the model for different designs, which can be very time consuming. In addition, the solver can be asked to evaluate each design at a number of frequencies and optimize with respect to the sum of the objective function values for each frequency. In this tutorial, a single frequency of 5000 Hz has been selected in order to make it possible to experiment with other aspects of the model. For example, changing the parameter θ, you can easily study the effect on the horn shape of optimizing the output at a specified angle from the axis.

Results and Discussion

By changing the shape of the horn within the limits of the selected parameterization, the on-axis sound pressure level can be raised by about 1 dB compared to the simple cone in Figure 1.

Figure 2: The final shape of the horn, optimized for on-axis SPL at 5000 Hz.

The improvement is rather small, because the initial design also shows a marked directivity, as can be seen from Figure 3. Obviously, the optimal shape with respect to on-axis SPL leads to deep undesirable minima in other directions.

Figure 3: Radiation plot of the original (dashed blue) and final (solid black) designs.

Optimizing with respect to a slight off-axis direction can give you a more uniform far-field pattern, but may also result in a deep minimum on the axis. Try for example to set the off-axis angle θ to 22° (in the parameters list set theta equal to 22[deg]).

Varying the frequency reveals that the on-axis improvement is robust over a wide frequency band; see Figure 4.

Figure 4: The objective function is plotted as a function of the frequency for the initial as well as the optimized design.

To search for a stable and practically useful horn design, you might instead create a composite objective function as a weighted sum of transmission values evaluated for a number of discrete directions, or choose to minimize the deviation from the mean SPL over a range of angles. In addition, you would also want to optimize with respect to more than one frequency, and experiment with different parameterizations.

Notes About the COMSOL Implementation

COMSOL Multiphysics implements the parameterization as a prescribed boundary displacement in the Moving Mesh feature. The mesh is allowed to move freely in the conical part of the horn, but otherwise kept fix. Some measures must be taken to avoid inverted elements when the shape of the cone is changed. Firstly, a quad mesh is used, since quads are less likely to become inverted, compared to triangles.

Secondly, the amplitude of the boundary displacement is restricted by the dmax parameter and the polynomial order. These constraints are intended to keep the mesh element volumes positive at all times and must not be active at the optimum point. You perform this sanity check as a final postprocessing step.

A time-harmonic Pressure Acoustics, Frequency Domain interface solves for the pressure field inside the horn and in a small spherical domain surrounding its opening. The air domain is terminated by a spherical PML layer which absorbs outgoing waves in such a way that the artificial termination of the domain has no influence on the near field. An accurate near field is sufficient, since the far-field result is based on an integral representation evaluated in the plane of the baffle. A plane wave radiation boundary condition on the waveguide attached to the narrow end of the horn feeds the horn with a plane wave of amplitude 1 Pa.

Set up the optimization problem in the study node. The pressure measured by an integration coupling variable, according to Equation 3, is inserted into a scalar objective function equal to the SPL value, or 10·log10(0.5·realdot(pm, pm)/(20·10−6)2).

Reference

1. E. Bängtsson, D. Noreland, and M. Berggren, “Shape Optimization of an Acoustic Horn,” Technical Report 2002-019, Department of Information Technology, Uppsala University, May 2002.

Acoustics_Module/Optimization/horn_shape_optimization

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.025.

Locate the section. In the text field, type -0.2.

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.3.

In the text field, type 90.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.1

Click

Click the

button in the toolbar.

Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

In the text field, type 0 0.1 0.1 0.025 0.025 0.

In the text field, type 0 0 0 -0.175 -0.175 -0.175.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
theta	0[deg]	0 rad	Polar angle
f0	5[kHz]	5000 Hz	Optimization frequency
df	50[Hz]	50 Hz	Optimization frequency bandwidth
Lfar	1[m]	1 m	Far-field evaluation distance

Definitions

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 6 only.

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
pm	intop1(besselj(0,acpr.krsin(theta))*pz/2/pi/Lfar)	N/m²	Measured pressure

In the definition of pm, intop1() is the name of your integration operator, r is the radial coordinate, acpr.k is the local wave number, theta is the observation angle, and pz is the derivative of the pressure with respect to the z-coordinate.

Global Variable Probe 1 (var1)

In the toolbar, click