Topology Optimization and Verification of an Acoustic Mode in a 2D Room

This tutorial introduces the use of topology optimization in acoustics. The goal of the optimization is to find the optimal material distribution (solid or air) in a given design domain, here the ceiling of a 2D room, that minimizes the average sound pressure level in an objective region. The optimization is carried out for a single frequency. The topology optimized design is further transformed into a geometry and the results of the optimization are verified using sound hard boundaries. See the tutorial Optimizing the Shape of a Horn for instructions on how to perform the optimization for a frequency band.

The model was inspired by the work of M. B. Dühring, O. Sigmund, and J. Søndergaard Jensen (Ref. 1).

This application requires the Acoustics Module and the Optimization Module.

Figure 1: The initial design has a pressure node at the edge of the objective domain.

Model Definition

The room is exited with a monopole point source. The initial sound pressure level is shown in Figure 1. The objective is to minimize the average sound pressure level in the circular domain illustrated in the figure. The optimization is allowed to distribute material in the top of the room (a layer of the ceiling). The top domain is constrained to consists of at least 85% air, but there are no other constraints, so the material does not have to connect with the top boundary.

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 nonanalytic 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.

The model uses the feature to solve the topology optimization problem using the density method. This means that the geometry is defined implicitly using nonphysical material parameters. The local value of a given material parameter is determined by the local value of a control variable field, θc, which is bounded between zero and one, with zero corresponding to some solid material and one to air.

Within the field of acoustic topology optimization it is common to interpolate the inverse density and bulk modulus linearly, that is,

(2)

where ρ1 and K1 are the density and bulk modulus of aluminum, respectively, while ρ2 and K2 are the properties of air. Wherever θ is equal to one, the properties of air are thus used. θ is the material volume factor, which is related to the control variable field, θc, through a filtering and projection step; see the Topology Optimization of an MBB Beam model in the Optimization Module Application Library for details. The main difference is that this model uses milling constraints to reduce the design freedom. Acoustic topology optimization can be sensitive to small design changes, which can easily occur during postprocessing of the topology optimization results. To reduce these issues, a continuation strategy is applied in the projection slope β. The model thus solves a sequence of optimization problems starting from a low projection slope and initializing the optimization with higher values using the previous optimization results.

Results and Discussion

The topology optimization achieves a large reduction of the average sound pressure level in the objective domain by moving a pressure node into the domain. The topology optimization result is further transferred to a new component using a Filter dataset. In this step a solid geometry is generated from the optimized solution. This causes the node to shift a bit, but it still goes through the domain as shown in Figure 2.

Figure 2: The sound pressure level is plotted after the optimization result has been transferred to a new component utilizing sound hard boundaries.

The value of the objective function is computed for a range of frequencies for the initial design as well as the topology optimized design (both before and after the verification step). The graph is plotted in Figure 3. It shows that the objective is only improved in a narrow frequency band around the optimization frequency. The verification step does not seem to cause the objective function to change noticeably.

Figure 3: The objective function is plotted as a function of frequency for the initial and topology optimized designs.

Notes About the COMSOL Implementation

The model uses milling constraints to reduce the design freedom, but this also causes the objective function to be worse than it would be more design freedom.

Reference

1. M.B. Dühring, O. Sigmund, and J.S. Jensen, Optimization of acoustic, optical and optoelastic devices, Kgs. Lyngby, Denmark: Technical University of Denmark (DTU), DCAMM Special Report No. S109, 2009.

Acoustics_Module/Optimization/topology_optimization_2d_room

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

Global Definitions

Parameters 1

You can type in the parameters listed below or load them from topology_optimization_2d_room_parameters.txt.

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
f0	34.5[Hz]	34.5 Hz	Frequency for optimization
W	18[m]	18 m	Room width
H	9[m]	9 m	Room height
dH	1[m]	1 m	Design domain height
Rob	1[m]	1 m	Objective domain radius
xob	16[m]	16 m	Objective domain x-coordinate
yob	2[m]	2 m	Objective domain y-coordinate
rho1	1.204[kg/m^3]	1.204 kg/m³	Air density
K1	141.921e3[Pa]	1.4192E5 Pa	Air bulk modulus
rho2	2643.0[kg/m^3]	2643 kg/m³	Aluminum density
K2	68.7[GPa]	6.87E10 Pa	Aluminum bulk modulus
volfrac	0.85	0.85	Design domain volume fraction
alpha_K	0.001	0.001	Damping coefficient
hmax	0.3	0.3	Maximum element size
rmin	0.1*hmax	0.03	Minimum element size
beta	32	32	Projection slope

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type W.

In the text field, type H.

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


Layer name	Thickness (m)
Layer 1	dH

Clear the check box.

Select the check box.

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type Rob.

Locate the section. In the text field, type xob.

In the text field, type yob.

Locate the section. Select the check box.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.

Locate the section. Select the check box.

Form Union (fin)

In the window, click .

In the window for , click

Design Domain

In the toolbar, click

and choose .

In the window for , type Design Domain in the text field.

Locate the section. In the text field, type H-dH/2.

Pressure Acoustics, Frequency Domain (acpr)

Pressure Acoustics 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

From the K list, choose . In the associated text field, type K1*(1+alpha_K*i).

From the ρ list, choose . In the associated text field, type rho1.

Monopole Point Source 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the QS text field, type 0.02.

Mesh 1

Free Triangular 1

In the toolbar, click

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Click to expand the section. Locate the section. Click the button.

Locate the section.

Select the check box. In the associated text field, type hmax/2.

Size 2

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Click the button.

Locate the section.

Select the check box. In the associated text field, type 0.1*hmax.

Size

In the window, under click .

In the window for , click to expand the section.

In the text field, type hmax.

Click

Definitions

Objective Function

In the toolbar, click

and choose .

In the window for , type Objective Function in the text field.

In the text field, type obj.

Locate the section. From the list, choose .

The realdot operator is analytic, so if we use this to define the probe, it can be used as objective function without enabling splitting of complex variables.

Locate the section. In the text field, type 0.5*realdot(p,p)/acpr.pref_SPL^2.

Study 1

Step 1: Frequency Domain

In the window, under click .

In the window for , locate the section.

Click

In the dialog box, type 24 in the text field.

In the text field, type 0.1.

In the text field, type 42.

Click .

In the window for , click to expand the section.

From the list, choose .

In the window, click .

In the window for , type Study 1 - Initial Design in the text field.

In the toolbar, click

Results

Acoustic Pressure (acpr)

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Sound Pressure Level (acpr)

In the window, click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Click the

button in the toolbar.

Acoustic Pressure (acpr), Sound Pressure Level (acpr)

In the window, under , Ctrl-click to select and .

Right-click and choose .

Study 1 - Initial Design

In the window for , type Study 1 - Initial Design in the text field.

Component 1 (comp1)

Set up the topology optimization by adding a and interpolate the inverse density and inverse bulk modulus linearly.

Density Model 1 (dtopo1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the Rmin list, choose .

In the text field, type hmax.

Use milling constraints to restrict the design freedom.

Click to expand the section. From the list, choose .

In the table, enter the following settings:


X	Y
0	-1

Locate the section. From the list, choose .

In the β text field, type beta.

Locate the section. From the list, choose .

Locate the section. In the θ0 text field, type 1.

Prescribed Material Boundary 1

In the toolbar, click

and choose .

Select Boundary 4 only.

Definitions

Design Domain Variables

In the window, under right-click and choose .

In the window for , type Design Domain Variables in the text field.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
rhod_inv	1/(1/rho2+dtopo1.theta_p*(1/rho1-1/rho2))	kg/m³	Design domain density
Kd_inv	1/(1/K2+dtopo1.theta_p*(1/K1-1/K2))	Pa	Design domain bulk modulus