Topology Optimization of Broadband Acoustic Transition Section: A Comparison between Deterministic and Stochastic Approaches

This paper focuses on the topology optimization of a broadband acoustic transition section that connects two cylindrical waveguides with different radii. The primary objective is to design a transition section such that it maximizes the transmission of a planar acoustic wave while ensuring the planarity of the transmitted wave. Helmholtz equation is used to model linear wave propagation in the device.We utilize the finite element method to solve the state equation on a structured mesh of square elements. Subsequently, a material distribution topology optimization problem is formulated to optimize the distribution of sound-hard material in the transition section. We employ two different gradient-based approaches to solve the optimization problem: namely, a deterministic approach using the method of moving asymptotes (MMA), and a stochastic approach utilizing both stochastic gradient (SG) and continuous stochastic gradient (CSG) methods. A comparative analysis is provided among these methodologies concerning the design feasibility and the transmission performance of the optimized designs, and the computational efficiency. The outcomes highlight the effectiveness of stochastic techniques in achieving enhanced broadband acoustic performance with reduced computational demands and improved design practicality. The insights from this investigation demonstrate the potential of stochastic approaches in acoustic applications, especially when broadband acoustic performance is desired.


Introduction
In many acoustic applications, acoustic waveguides connect various parts of a device from the transducer to the radiating aperture or receiver [1,2].This includes the transmission of acoustic waves between parts with different geometry, dimensions, and acoustic impedances.In most cases, a minimum reflection and alteration of the incoming signal is desired, which can be achieved by using an impedance matching section [3,4].
One of the earliest works on acoustic transition sections is the study by Kirby [5], who used numerical simulations to model acoustic wave propagation in two waveguides connected through a transition section.Wadbro [3] used a material distribution topology optimization method to design a transition section between two cylindrical waveguides with different radii to achieve impedance matching.Robertson et al. [4] considered a similar problem of impedance matching between two cylindrical waveguides, although they used a twoneck Helmholtz resonator as the transition section.They showed that perfect impedance matching can be achieved by tuning the dimensions of the Helmholtz resonator.However, this approach has bandwidth limitations with impedance matching achieved only at a narrow range of frequencies close to the resonance frequency of the resonator.Cao et al. [6] took a different approach by considering a two-dimensional acoustic transformation section with an impedance-tunable transformed medium.They showed that desirable broadband impedance matching can be achieved in this way, though in practice, it is very difficult to set a transformation medium with acoustic properties changing according to a given function [6].
Here, we use the material distribution topology optimization method, also known as density-based topology optimization, to design an acoustic transition section for impedance matching.This is a common method in computational design optimization for acoustic problems [7][8][9][10].To obtain a final topology that is favorable for a broad range of frequencies, there are two fundamentally different approaches: First, the problem can be viewed as a topology optimization task for infinitely many load cases.Using an appropriate discretization scheme, the problem can be transformed into a deterministic multi-load formulation, where the objective function admits a finite sum structure [11][12][13], for which various deterministic optimization schemes have been explored.In this contribution, the resulting optimization problems are solved via the method of moving asymptotes (MMA) [14].We investigate how the solution depends on the number of frequencies considered in the optimization.
Changing perspective, the original objective can also be considered as a robust optimization problem [15][16][17], where the broadband frequency range models an underlying uncertainty of load cases.While many approaches for robust optimization again rely on discretizations or series expansions of the full objective, we specifically choose two methods from stochastic optimization which do not follow this philosophy.Namely, we optimize the transition section using the stochastic gradient descent [18] and the continuous stochastic gradient descent [19][20][21] method.Both of these approaches represent probabilistic, sample-based optimization schemes.
The remainder of this paper is structured as follows.In Section 2, we introduce the problem setup, including the geometric configuration and the governing equations.We also discuss the discretization of the state equations using the finite element method.Next, we introduce the objective function for the design optimization problem aiming for a broadband transition.For the different optimization approaches considered in this work, we present the corresponding results of our numerical experiments in Section 3. Finally, in Section 4, we present a concluding discussion that includes a comprehensive comparison of the results obtained by the different approaches.We summarize the main findings and provide insights into the strengths and limitations of each method.

Problem Statement
Consider the cylindrical setup illustrated in Fig. 1, consisting of two semi-infinite pipes connected by a transition section.Assume a planar acoustic wave propagating from left to right in the left pipe.As this incoming wave propagates in the transition section Ω D , highlighted in grey, a part of the wave will propagate through the transition section to the right pipe, while another part will be reflected back to the left pipe.By optimizing the distribution of sound-hard material in the transition section Ω D , this study aims to ensure that the planar incoming wave in the left pipe continues to propagate as a planar wave to the right pipe, despite the change in the diameter of the waveguide.A similar problem was previously considered by Wadbro [3] for a narrower range of frequencies.In this study, we aim to extend the analysis to a broader range of frequencies, which is of practical importance for a wide range of acoustic applications.Additionally, we compare two distinct approaches to solving the optimization problem: the deterministic approach and the stochastic approach.Here and throughout this article, the waves that propagate away from the transition section are termed outgoing waves, and the waves that propagate towards the transition section are termed incoming waves.We note that the outgoing waves are further divided into reflected waves, traveling to the left in the left pipe, and transmitted waves, traveling to the right in the right pipe.

Mathematical model
In this study, we consider linear wave propagation in the cylindrically symmetric setup, illustrated in Fig. 1(a).Specifically, we let P (x, t) = Re p(x)e iωt denote the time-harmonic pressure in the air-filled region.Under these assumptions, the complex pressure amplitude p(x) satisfies the Helmholtz equation in cylindrical coordinates.That is, which describes the behavior of sound waves in the system.Therein, Ω is the air-filled region of the setup, and the wave number k = ω/c is determined by the speed of sound in air c and the angular frequency ω.As mentioned earlier, the design domain Ω D (transition section) may partly be filled with sound-hard material.Fig. 2 shows the air-filled domain Ω and an arbitrary distribution of sound-hard material Ω s in the design domain.By varying the distribution of sound-hard material within the design domain, we can explore how this affects the propagation of sound waves through the transition section and identify optimal designs that minimize wave reflection and ensure planar wave transmission.
Fig. 2 The computational domain Ω comprising the design domain Ω D and two truncated waveguides.An arbitrary sound-hard material distribution inside Ω D is given as Ω s .The remainder of the setup Ω = Ω \ Ω s is filled with air.The boundaries are divided into an axi-symmetric axis Γ sym (dash-dotted line), the sound-hard walls Γ s (solid line) and the boundaries Γ L and Γ R (dashed lines) of the artificially truncated waveguides.
To numerically approximate the infinitely long pipes on both sides of the design domain, we truncate them and use DtN (Dirichlet-to-Neumann) non-reflecting boundary conditions at the artificial boundaries Γ L and Γ R , as illustrated in Fig. 2. Regarding more details about this type of artificial boundary conditions, we refer the reader to the book by Ihlenburg [22] and the appendix of the article by Wadbro [3].Considering these artificial boundary conditions, we obtain the boundary-value problem Conditions (2c) and (2d) ensure that all the outgoing waves are perfectly absorbed.Furthermore, condition (2c) also specifies an incoming planar wave with unit amplitude at Γ L .By multiplying equation (2a) with a test function q and integrating over the domain Ω, the variational form of boundary-value problem (2) can be written as follows.
Find p ∈ H 1 ( Ω) such that: For a given distribution of solid material in the design domain and a given shape of region Ω s , the solution p to equation (3) shows the distribution of the complex pressure in Ω.Following a standard approach in topology optimization, we define a material indicator function α such that α ≡ 0 in Ω s and α ≡ 1 in Ω.Using this function, we extend the integration domain of the domain integrals in variational formulation (3) from Ω to Ω = Ω ∪ Ω s .We note that the computational domain Ω consists of both the air and the solid regions.The resulting reformulation of variational formulation (3) is then given by The solution p to the variational formulation (4) represents the distribution of complex pressure in the waveguide, given a design of solid scatter in Ω D , described by a material indicator function α.

Discretization
We use the finite element method to discretize and numerically solve problem (4) on a structured grid of square elements.Let V be a finite element functional space consisting of continuous and bi-quadratic functions on each element, and let φ j , j = 1, 2, . . ., N be bi-quadratic shape functions, where N is the number of degrees of freedom.Thus, V = span {φ 1 , φ 2 , . . ., φ N }.We approximate the complex pressure p and the test function q by p h ∈ V and q h ∈ V , respectively.Additionally, we approximate the material indicator function α with an element-wise constant function α h .Using the above definitions and approximations, we obtain the discretized version of problem (4) below.
where DtN h represents semi-discrete Dirichlet-to-Neumann type boundary operators on Γ L and Γ R [3, Appendix A].The algebraic or matrix formulation of problem (5) reads T is the vector of nodal values of the complex acoustic pressure amplitude, T is the vector that holds the element values of α h (with N D denoting the number of elements in Ω D ) and 1 N = [1, 1, . . ., 1] T is a vector of length N .Also, the N ×N stiffness K, mass M, and boundary mass M L matrices have components respectively.The boundary matrices B L and B R represent the non-reflecting boundary conditions at Γ L and Γ R , respectively.A detailed derivation of B L and B R is provided in a previous work by Wadbro [3, Appendix A].

Power of outgoing waves
Let Helmholtz equation ( 1) govern the distribution of the complex pressure p in the two semi-infinite pipes on the left and right side of the transmission section illustrated in Fig. 1(a) with sound-hard boundary condition (2b) on the solid walls.Using the separation of variables, the general solution for p in the left and right pipes reads respectively, where z L and z R are the position of z-axis on Γ L and Γ R , the functions f m (r) are the modes at left and right waveguide, e is the base of the natural logarithm, i is the imaginary unit and A L m and B L m are complex constants that determine the amplitude of incoming and outgoing waves at the left waveguide, respectively.Similarly, A R m and B R m are complex constants that determine the amplitude of incoming and outgoing waves at the right waveguide, respectively.Lastly, the constants k m are the so-called reduced wave numbers.In the continuous case, we have an infinite number of modes m = 0, 1, 2, . .., but only the modes with real-values of k m are propagating modes and the ones with imaginary k m will decay exponentially according to the equations (8).These modes are known as evanescent modes.
The mode functions f m (r) should satisfy the following one-dimensional eigenvalue problem in the radial direction on the boundaries Γ L and Γ R : where W is the radius of the pipe.In the continuous case, it is well known that the functions f m are so-called Bessel functions.For our numerical treatment, we can extrapolate the numerical solution on Γ L and Γ R to any point in the waveguides using expansion (8) and viewing the problem continuously in the lengthwise direction and discretely in the radial direction.Note that for a given finite element discretization of equation ( 9), the number of modes that are representable in the discretized case M equals the number of basis functions with support on the boundary.Let f h m be an eigenfunction (mode) corresponding to eigenvalue λ m , where m = 0, 1, . . ., M .Since the complex pressure p satisfies Helmholtz equation ( 1), for the reduced wave number we have Recall that f h m is a propagating mode if its corresponding reduced wave number is real.The smallest eigenvalue in solving problem (9) is 0, which corresponds to the planar wave.So λ 0 = 0 and k 0 = k, and thus, the planar wave mode is always a propagating mode.Moreover, for a given frequency f , there is a finite number M p of propagating modes satisfying the condition λ m ≤ k 2 .The number of propagating modes depends on the frequency of the wave and the radius of the pipe.Thus, for the different radii of the left and right pipes, we may have a different number of propagating modes at Γ L and Γ R , which we denote by M p L and M p R , respectively.In the discretized case, the solution at the boundaries Γ L and Γ R , where z L − z = 0 and From now onward, we occasionally use the superscript X in an expression to represent either L for the left waveguide or R for the right waveguide.Note that the corresponding statement holds in both cases of replacing X by either L or R referring to the left and right waveguides, respectively.Let f h n be the nth propagating mode.Then, we have where the first equality follows from substituting p Γ X from equation (11) into the first expression and the second equality follows from the orthonormality of modes.Considering v X n to be the N × 1 vector representing the nodal values of the discrete mode f h n on the boundary nodes at Γ X (note that all other entries of v X n corresponding to the internal nodes are zero), equation (12) where M X is the boundary mass matrix as defined in equation (7c) at Γ X , and p is the nodal values of the complex pressure.Thus, for a given solution p to problem (6), we can recover the complex amplitudes of the incoming and outgoing waves for each propagating mode at Γ X , using equation (13).
In this study, we only consider the case where we have a planar incoming wave with unit amplitude at Γ L .Therefore, we can rewrite equation (13) as Recall that M p L and M p R are the highest propagating modes in the left and right pipes, respectively.
The power of a propagating wave is proportional to the square of its amplitude and its corresponding reduced wave number.Defining the normalized power of outgoing waves as the power of outgoing wave divided by the power of the unit-amplitude incoming wave and considering equation (14) to compute the amplitude of outgoing waves, we have and ) Here P L m and P R n are the normalized power of the outgoing waves of mode m and n at Γ L and Γ R , respectively.A detailed derivation of the power of outgoing waves for a given amplitude in the discretized case is provided by Wadbro [3,Appendix B].Here and throughout this article, whenever the term power of outgoing wave is used, it means the normalized power of the outgoing wave by dividing the power of the outgoing wave with the power of a unit-amplitude incoming wave imposed at Γ L .

Objective function
As mentioned earlier, the aim of this study is to design the transmission section in Fig. 1 to (i) maximize the transmission and (ii) ensure that the transmitted wave is planar as illustrated in Fig. 1(c).To achieve this, we minimize the sum of the power of all outgoing waves except for the planar wave to the right over the targeted range of frequencies F .Thus, the primary objective function can be written as (17) Note that we have normalized the objective by the length of the targeted frequency range |F |.Note that for each frequency, we need to calculate M p L (f ) and M p R (f ), the number of propagating modes at Γ L and Γ R , respectively.
For binary values of α h = {0, 1}, the optimization problem with objective function ( 17) is a large-scale non-linear integer optimization problem.Also, if α h = 0 in some of the elements, then the system matrix in equation ( 6) becomes singular.To solve the numerical and mathematical issues that arise when solving this problem, a standard approach in topology optimization is to relax the binary value constraint and let α h take values in the range [ϵ, 1], where ϵ is a small number [3,7,9,23,24].Moreover, we aim for a pure solid (α h = ϵ) or air (α h = 1) final design.Thus, we use a combination of filtering and penalty methods to suppress the intermediate values.The non-linear density filters used in the numerical experiments also ensure a size control on the solid region in the design [9,[25][26][27][28] T be the vector of design variables before filtering.Thus, we define the N D × 1 vector α := F (d), where F is a filter operator.To further suppress the intermediate values of the design variables, we add a standard quadratic penalty term [3,9,29,30] to the primary objective function (17) and thus, we define the objective function for the numerical experiments as follows: where γ is the penalty parameter, |Ω D | denotes the size of the design domain.

Numerical experiments
In our numerical experiments, we consider the setup illustrated in Fig. 2 with the following dimensions: The radius and length of the design region Ω D is r D = 50 mm and l D = 50 mm, respectively.The radius and length of the truncated waveguides are r L = 30 mm, r R = 40 mm, and l W = 20 mm, respectively.We aim to maximize the transmission of the planar incoming wave in the frequency range of 4-16 kHz, ensuring that the transmitted wave is also planar.We discretize the computational domain into a structured grid of square elements with a uniform mesh size of h = 0.25 mm, resulting in 250,721 degrees of freedom for the finite element discretization.To solve the optimization problem, we employ three different optimization algorithms: the MMA method [14], the stochastic gradient (SG) method [18], and the continuous stochastic gradient (CSG) method [19][20][21].
We define the performance of a given design at frequency f as the normalized power of the outgoing planar wave, computed using expression ( 16), as follows: To evaluate the performance of the optimized designs, we use a boundary-fitted mesh for the final designs in the Acoustics Modules in COM-SOL Multiphysics.The performance of different designs are compared over a range of frequencies from 4 kHz to 16 kHz with a step size of 20 Hz.

MMA approach
To solve the optimization problem considering objective function (18) using the MMA method, we approximate the integral over the range of targeted frequencies using the function values of the integrand at just a few frequencies.Thus, we discretize the optimization problem as where Q is the number of frequencies subject to optimization, 1 N D is N D ×1 vector with all entries equal to 1, and is the set of admissible designs.The scaling constant Q −1 is neglected in expression (20).This is done because the scaling between the primary objective function and the quadratic penalty term can be also tuned using the penalty parameter γ.
Note that by increasing the number of frequencies subject to optimization Q, we can improve the approximation used to discretize objective function (18).To solve optimization problem (20), we utilize the least squares formulation of the MMA approach described by Svangberg [14,31].Thus, we need the sensitivity information for each part of the objective function in optimization problem (20).The computation of sensitivities for the quadratic penalty term can be readily performed for a given filter F. However, the task of computing the gradient of the power of outgoing modes with respect to the design variables poses a challenge.This process involves determining the gradient of the amplitudes of the outgoing modes with respect to the design variables, as per equations ( 15) and ( 14).Notably, the power of a propagating mode is proportional to the square of its amplitude.The sensitivity computations are done analytically using the adjoint variable method, which is a powerful technique for computing the gradient of a function that depends on the solution of a partial differential equation with respect to the design variables.A detailed derivation of the design sensitivities is provided in the appendix.
The design problem has 40,000 design variables, that is, the number of elements in the design domain Ω D .We set the lower bound for the design variable as ϵ = 10 −8 , the filter radius to 1 mm, and use a so-called continuation approach for the penalty parameter.That is, we solve problem (20) for a sequence of increasing penalty parameters γ i = 10 i , i = 0, 1, . . ., 5, using the previously computed solution as the initial design.The aim is to gradually move the optimizer's focus from the acoustic performance of the device towards obtaining a black-and-white final design.This approach ensures an optimized layout with sharp solid-fluid boundaries, free of any intermediate values of the material indicator functions [7].To solve optimization problem (20), we need to consider a sufficiently large number of frequencies Q to get broadband acoustic performance.However, increasing the number of frequencies subject to optimization will increase the number of times we need to solve the state equation (6).Note that the finite element solver is the primary contributor to the computational costs in the optimizer.Consequently, increasing the number of frequencies will result in a significant increase in computational
Note that here, the convergence criteria is based on the residual norm of the KKT (or the first-order optimality) condition, together with a limitation on the number of iterations in each penalty step.We will use the number of evaluations, defined as the number of times we need to solve the state equation ( 6), as a metric to compare the computational cost between different cases.Fig. 3 shows the optimized design for all three cases.Fig. 4 shows the performance of each of these designs, computed using expression (19).Fig. 5 shows the convergence history for all three MMA cases, where two numerical approximations of the objective function ( 17) are plotted versus the number of evaluations, one is the naive approximation using only the frequencies subject to optimization in each case and the other one is achieved using 150 equidistant frequencies in the range 4 kHz to 16 kHz.
Considering the observations from Fig. 4 and Fig. 5, the following conclusions can be drawn: • In Case I, where only a few frequencies within the targeted range are considered in the optimization, the overall broadband performance of the final layout is poor, with only a few peaks observed at the considered frequencies.The total number of evaluations required for the convergence is approximately 500 evaluations.• Case II, by considering more frequencies to approximate objective function (17), a better overall performance is achieved for the optimized design.However, there are still deeps in the frequency response, indicating weak broadband performance.In this case, Fig. 5 shows that the approximation of the objective function is still poor throughout the optimization.Moreover, approximately 2200 evaluations were required for convergence.• In Case III, a desirable broadband performance is achieved for the optimized layout.However, this improvement comes at the expense of an increase in computational cost, as evident from the approximately 2700 evaluations required to obtain the optimized

Evaluations
Fig. 5 History of convergence in all three cases.The black line indicates the approximation of the objective function (17) using numerical integration utilizing the evaluation of the integrand only at the targeted frequencies in each case.The dash-dotted colored lines demonstrate the true objective function approximated by numerical integration using 150 equidistant frequencies.
layout.As illustrated in Fig. 5, the approximation of the objective function is significantly improved throughout the optimization compared to Case I and Case II.These observations highlight two disadvantages of using a deterministic approach in this problem which can be manifested as follows: 1. To ensure an acceptable broadband performance, an increased number of frequencies needs to be included in the optimization process.This results in a higher number of evaluations and computational costs.2. Evaluating objective function (17) only at specific frequencies leads to designs that are tailored for those particular frequencies.
As the number of frequencies considered in the optimization increases, the resulting designs exhibit numerous free-hanging parts, as depicted in Fig. 3.This phenomenon is characterized by an increasing number of small inclusions from Case I to Case III in Fig. 3.In summary, while the MMA approach demonstrates advantages in achieving desirable broadband performance, it is pivotal to consider the associated computational cost and the design's specificity to the frequencies included in the optimization process.

SG approach
In contrast to MMA, SG does not require a computation of the integral over frequencies appearing in (18).Instead, in each iteration, we draw a random frequency f n ∈ F and evaluate the objective function gradient only for this specific choice.This gradient sample is then used as a search direction for the current optimization step.As a result, compared to MMA, an SG iteration consumes significantly less time, but generally provides a smaller improvement of the objective function.
While SG is typically used with a diminishing learning rate, we fix a constant learning rate and impose a shrinkage in step length directly through move limits.To be precise, in every iteration, the search direction is multiplied by the constant learning rate and afterward projected onto the set where C n = O(1/ √ n).To end up with a blackand-white design, the final result is rounded.This setup yielded the best performance for SG in our numerical experiments.
Due to the stochastic nature of SG, we performed 50 independent optimization runs with 500 iterations each.An overview of the observed performance range is found in Fig. 10.For the best final design obtained, the objective function evolution is shown in Fig. 6, while the final design is illustrated in Fig. 7.

CSG approach
Similar to SG, the CSG search direction is based on stochastic samples of the full objective function gradient.However, since evaluating such a  8 Performance of the optimized designs in the SG approach, the dash-dotted orange line, and in the CSG approach, the blue line.
gradient sample is still computationally expensive, discarding all information after each iteration is rather wasteful.Therefore, in CSG, gradient samples (g i ) i=1,...,n from past iterations are stored.By calculating design-dependent integration weights (β i ) i=1,...,n , the full objective function gradient is then approximated by the continuous stochastic model For our numerical experiments, we choose socalled exact hybrid weights, since they are easily computable due to dim(F ) = 1.More details on how these weights are computed in practice is given by Grieshammer et al. [20].Therein, it was also shown that the approximation error vanishes during the optimization process.That is, As a consequence, CSG inherits strong convergence results from full gradient schemes, like convergence for constant learning rates and line search techniques, while retaining a low cost per iteration, since the integration weight computation is negligible compared to the numerical solution of the state equation.
For a better comparison, we also choose a combination of constant learning rates and move limits for CSG in our experiments.This time, however,

Evaluations
Fig. 9 CSG approximation to the objective function value (17) in each iteration (black).The dash-dotted blue line gives a different approximation to the objective function value, obtained by using numerical integration with 150 equidistant frequencies.
it is not necessary to pick diminishing move limits.Instead, we can adaptively choose these limits based on the progress achieved in the internal CSG model for the objective function.
Again, 50 independent optimization runs with 500 evaluations each were performed.The full overview of results is shown in Fig. 10.The objective function evolution of the run corresponding to the best final result obtained is given in Fig. 9. Therein, we also included the history of objective function approximations by CSG.These approximations indicate the quality of the underlying continuous-stochastic model, which is used internally in CSG.An illustration of the final design can be found in Fig. 7.

Impact of Stochasticity
To better capture the probabilistic nature of SG and CSG, the transmission spectra of all 50 final designs are calculated.Afterwards, for each individual frequency f ∈ F , we determined the respective transmission quantiles P 0.1,0.9(f ) and P 0.25,0.75(f ).Here, P 0.1,0.9(f ) denotes the range of transmission values achieved by all runs, where the highest 10 % and lowest 10 % of values are omitted.Likewise, P 0.25,0.75(f ) indicates the range of transmission values obtained by 50 % of all runs, where the best 25 % and worst 25 % of results are neglected.Thus, the resulting quantile plots in Fig. 10 give a good impression concerning both the average performance of a design obtained by SG or CSG as well as the variance in results, depending on the random sequence of sample frequencies.Note, however, that this form of representation results in smoother-looking spectra, since sharp peaks and other resonance effects are averaged out in the process.

Discussion
To better compare the achieved performances for all methods, we introduce the cumulative performance density CPD : By construction, CPD(p) provides the relative amount of frequencies in our frequency range F , for which the performance is lower than the given threshold p ∈ [0, 1].For example, if a design satisfies CPD(0.5)= 0.25, its performance is less or equal to 0.5 for 25 % of all considered frequencies.Thus, an ideal design, which has a perfect performance of 1 for all f ∈ F satisfies Furthermore, the objective function value (17) 11 Cumulative performance density curves for all final designs in Fig. 3 and Fig. 7 as well as the empty design domain.For a given final design, the CPD curve indicates two measures of quality: the objective function value (area under the graph) and the median performance (intersection with horizontal line at y = 0.5).
interval [0, 1], that is, For each of the optimization approaches, the CPD of the corresponding final design is given in Fig. 11.Therein, we also included the CPD for the empty design domain (α h ≡ 1).As we can see, the final design obtained with MMA in Case I yields the worst performance overall, even falling behind the empty design region.The final designs of SG and MMA Case II have comparable median performances.However, while the final SG design performs rather similar for all frequencies, indicated by the sharp increase of CPD at ∼ 0.9, the MMA Case II design performs poorly for a wide range of frequencies, resulting in a much better final objective function value of SG.The best overall performance is achieved by the final design of MMA Case III, with the final CSG design performing slightly worse.However, recall that the associated numerical effort for MMA Case III is much higher (2704 state equation solutions) when compared to CSG (500 state equation solutions).

Conclusion
In this study, we presented the results of topology optimization applied to a broadband acoustic transition section.We compared the outcomes of a deterministic approach using the method of moving asymptotes (MMA) with two stochastic approaches: stochastic gradient (SG) and continuous stochastic gradient (CSG) methods.
In the case of the MMA approach, we found that achieving optimal broadband performance requires optimizing over an increased number of frequencies (Fig. 4).However, this comes at the cost of a significant increase in computational costs and results in designs with a higher number of free-hanging inclusions, which can negatively impact manufacturability (Fig. 3).On the other hand, the stochastic approaches, SG and CSG, offer a more computationally efficient alternative while still producing optimized designs with improved manufacturability.In particular, we observed that the CSG method outperforms the SG method, as evidenced by the median transmission spectrum of the final designs and the overall frequency response shown in the quantile plots (Fig. 10).
These findings highlight the potential of stochastic approaches in acoustic applications, especially when broadband acoustic performance is desired.By reducing computational costs and improving manufacturability, stochastic methods offer a promising alternative for optimizing acoustic systems in such applications.element E as where p E and (z L m ) E are vectors containing the components of p and (z L m ) corresponding to nodes in element E and M E and K E are the element mass and stiffness matrices for element E, respectively.
Considering z R n to be the solution to the problem (6) with the right-hand side replaced by M R v R n , following the same argumentation as for z L m , we obtain the partial derivatives of B R n with respect to the element values of α E in element E as Note that z L m and z R n are known as adjoint variables.

Fig. 1
Fig. 1 (a) Axi-symmetric setup with cylindrical design domain Ω D in the middle and one cylindrical waveguide on both sides.(b) A 3D visualization of the setup.(c) The targeted wave propagation characteristics in the waveguides, where the planar incoming wave in the left pipe continues to propagate as a planar wave to the right pipe.

Fig. 3
Fig. 3 On the left is an axisymmetric and on the right is a 3D visualization of the optimized designs in (a) Case I, (b) Case II, and (c) Case III

Fig. 4
Fig. 4 Performance of the optimized designs in (a) Case I, the red line (b) Case II, the dash-dotted teal line and (c) Case III, the dashed violet line.

Fig. 6 Fig. 7
Fig.6Value of the objective function(17) for all intermediate SG iterates.Each value was approximated by numerical integration using 150 equidistant frequencies.
Fig.8Performance of the optimized designs in the SG approach, the dash-dotted orange line, and in the CSG approach, the blue line.