Acoustic topology optimization of sound absorbing materials directly from subdivision surfaces with isogeometric boundary element methods

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 362 (2020) 112806 www.elsevier.com/locate/cma

Acoustic topology optimization of sound absorbing materials directly from subdivision surfaces with isogeometric boundary element methods Leilei Chena , Chuang Lua,d , Haojie Lianb , Zhaowei Liuc , Wenchang Zhaod , Shengze Lig , Haibo Chend , Stéphane P.A. Bordasb,e,f ,∗ b

a College of Architecture and Civil Engineering, Xinyang Normal University, PR China Institute for Computational Engineering, Faculty of Science, Technology and Communication, University of Luxembourg, Luxembourg c Glasgow Computational Engineering Centre, School of Engineering, University of Glasgow, UK d CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, University of Science and Technology of China, PR China e School of Engineering, Cardiff University, The Parade, CF24 3AA, Cardiff, UK f China Medical University Hospital, China Medical University, Taichung, Taiwan, ROC g National Innovation Institute of Defense Technology, Beijing, PR China

Received 25 August 2019; received in revised form 1 November 2019; accepted 21 December 2019 Available online xxxx

Abstract This paper presents an acoustic topology optimization approach using isogeometric boundary element methods based on subdivision surfaces to optimize the distribution of sound adsorption materials adhering to structural surfaces. The geometries are constructed from triangular control meshes through Loop subdivision scheme, and the associated Box-spline functions that generate limit smooth subdivision surfaces are employed to discretize the acoustic boundary integral equations. The effect of sound-absorbing materials on the acoustic response is characterized by acoustic impedance boundary conditions. The optimization problem is formulated in the framework of Solid Isotropic Material with Penalization methods and the sound absorption coefficients on elements are selected as design variables. The potential of the proposed topology optimization approach for engineering prototyping is illustrated by numerical examples. c 2020 Elsevier B.V. All rights reserved. ⃝ Keywords: Acoustics; Topology optimization; Isogeometric analysis; Subdivision surface; Boundary element methods; Adjoint variable method

1. Introduction Noise control is an important issue in numerous engineering and musical applications. Adhering sound absorbing materials to the surfaces of structures can effectively reduce noise level. Subject to the constraint of weights and costs, full coverage of sound absorption materials on structural surfaces is impractical. Topology optimization offers ∗ Corresponding author at: Institute for Computational Engineering, Faculty of Science, Technology and Communication, University of

Luxembourg, Luxembourg. E-mail address: [email protected] (S.P.A. Bordas). https://doi.org/10.1016/j.cma.2019.112806 c 2020 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

a valuable design automation tool to optimize the material distribution [1–3]. In structural topology optimization, the finite element method (FEM) is most popular [4–7], but for infinite domain problems like acoustics, the Boundary Element Method (BEM) is rather appealing, as only surfaces need to be discretized and because the domain need not be truncated [8,9]. Traditional FEM and BEM with Lagrange polynomial basis functions both rely on polygonal meshes, resulting in cumbersome preprocessing procedures, geometrical errors, and low continuity of the field variables [3]. These drawbacks can be overcome thanks to the concept of isogeometric analysis (IGA) proposed by Hughes et al. [10]. The essence of IGA is to employ the same basis functions as those used in computer-aided design (CAD) to represent the domain geometry to discretize the unknown physical fields when solving partial differential equations. IGA is particularly attractive for topology optimization [11–13] as (1) the design variables are linked with exact geometries; (2) any geometry update is directly inherited by the field approximation; (3) the final geometry need not be postprocessed for visualization purpose; (4) using Geometry Independent Field approximaTion (GIFT) space can be locally refined whilst keeping the geometry untouched [14]. IGA was initially developed in the context of Finite Element Methods (IGAFEM), but in some circumstances it is desirable to combine IGA with Boundary Element Methods (IGABEM) [15–18]. The advantages of IGABEM are threefold: (1) it is capable of utilizing CAD data immediately without needing volume parameterization from geometric surfaces [19–21], since both of them are boundary-represented [22–26]; it inherits the merits of conventional BEM in addressing infinite domain problems [27–29]; (3) IGABEM is ideal for free boundary problems such as crack growth as no volume parameterization is needed [18,30]. IGA develops with the advance of CAD modeling techniques. As the industrial standard in CAD, Non-Uniform Rational B-splines (NURBS) are widely used in IGA [31–34]. Nevertheless, NURBS do not naturally permit local adaptive refinements and often lead to gaps or overlaps across patches in complicated geometries [35–37]. Sederberg et al. introduced T-splines to address these issues [38]. T-splines enable T-junctions in control meshes, allow for local refinement and yield watertight geometries with one single control mesh. Despite making a significant step forward, T-splines are not without limitations. T-spline basis functions are not guaranteed to be linearly independent; only a subset of T-spline models under some restrictions on the T-mesh is analysis-suitable in IGA [39]. In addition, the complexity of knot insertion algorithms for adaptive refinements poses implementation difficulties. The review article [40] provides a detailed survey to the development of local refinable splines, including PHT-splines [41–46], Hierarchical B-splines [47,48], and LR B-splines [49,50], and identifies the strengths and weaknesses of these methods. Another promising alternative is taking subdivision surfaces as CAD models in IGA. Compared with NURBS or other competitors, a salient feature of subdivision surface is that it is able to address extraordinary vertex. In addition, its recent advance also allows exact description of conic sections [51]. Subdivision surfaces were first introduced in [52,53] by generalizing knot insertion rules for B-splines on regular grids, which were then extended to triangular B-splines using the so-called Loop scheme [54]. Subdivision surfaces are ubiquitous in animation and computer games, and are now also available in various CAD packages such as Rhinoceros 3D, CATIA and SolidWorks. Subdivision surfaces exhibit flexibility for representing complicated topologies [55–61] while preserving watertight geometries and supporting local refinement. Furthermore, subdivision surfaces possess salient features for analysis such as linearly independent basis functions and the capability of addressing extraordinary vertices. They are also easy to implement, which make them amenable to numerical analysis and topology optimization. Numerical analysis using subdivisions surfaces has been applied in many engineering applications, such as shell analysis [62,63], shell-acoustic coupling problems [64], electrostatics [65,66], and shape optimization [67,68]. This paper presents an acoustic topology optimization paradigm using isogeometric boundary element methods based on subdivision surfaces (SS-IGABEM). The topology optimization algorithm adopts the Solid Isotropic Material with Penalization (SIMP) method for its versatility and efficiency. The Burton–Miller formulation that combines conventional boundary integral equation with its normal derivative is used to eliminate fictitious eigenfrequencies [69]. The sensitivity analysis is conducted through the adjoint variable method [70–73] for its advantage over direct differentiation methods [74,75] in handling a large number of design variables. Furthermore, we use a bespoke Fast Multipole Methods (FMM) to accelerate both the acoustic analysis and the sensitivity analysis [74]. We organize the paper as follows. Section 2 reviews the preliminary concepts of subdivision surfaces. Section 3 introduces IGABEM discretization for exterior acoustic analysis with subdivision surface basis functions. Section 4 discusses the topology optimization procedure including the optimization formulations, sensitivity analysis based on adjoint variable methods and the admittance model. Section 5 provides numerical examples to demonstrate the validity of the proposed approach, followed by conclusions in Section 6.

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

3

Fig. 1. Subdivision refinement of an initial control polygons (left column). The three polygons in the right most column are generated by repeated subdivision of the initial polygon. S 1 represents one level of subdivision, S 2 stands for two level of subdivision. Notice that the polygon gradually becomes smooth as the number of refinements increases. In the limit, the curves tend to cubic B-spline curves.

Fig. 2. Vertex and edge points stencils for segment (red circles represent modified vertices and green circles represent new edge points). According to the weights, the solid line on the right indicates the next level control polygon, and the dashed line is consistent with the segments on left. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2. Loop subdivision surface For the sake of completeness, we briefly review the concept of Loop subdivision surfaces. The interested reader can refer to [62,63] for more details. 2.1. Loop subdivision scheme Geometric modeling with subdivision surfaces consists in subdividing the coarse mesh recursively until the desired limit surface is reached. The Loop refinement scheme falls into the category of approximate refinement schemes. Not only new points are introduced in the refinement procedure, but existing points are also relocated. To illustrate the refinement scheme, let us start from subdivision curves as a simple example. Fig. 1 shows how the boundary of a polygon approaches a smoothed closed curve with increasing subdivision steps. As shown in Fig. 2, k each edge [xik , xi+1 ] is first split into two segments by inserting a mid-point, where k is level number and i is vertex number. For a given coarse polygon of level k with vertex coordinates xik , a refined polygon of level k + 1 with vertex coordinates x2ik+1 is computed as x2ik+1 =

1 k 3 1 k x + x k + xi+1 8 i−1 4 i 8

(1)

k+1 and edge points x2i+1 are calculated as mid-point:

1 k 1 k x + xi+1 (2) 2 i 2 Loop subdivision surfaces are constructed in a similar manner as subdivision curves. Fig. 3 demonstrates how a triangular mesh converges to a smooth surface through a number of subdivisions using the Loop scheme. The refinement of a triangular mesh by quadrisection is illustrated in Fig. 4. The number of edges connected to a vertex is referred to as the valence of this vertex. A vertex is regular when its valence N = 6 and irregular when N ̸= 6. k+1 x2i+1 =

4

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

Fig. 3. Subdivision refinement of a complex model, the control mesh is generated with Loop subdivision. The number of elements in initial control mesh is 488, and the number of elements subdivided once is 1952, increased fourfold.

Fig. 4. Vertex and edge points stencils for regular and irregular vertices and edges (red circles represent modified vertex and green circles represent new edge points, The dashed lines emphasize that the vertex applies to arbitrary valence N ). Refinement of mesh by twice quadrisection, the valence of an irregular vertex does not increase as the number of refinements increases, but the number of regular points increases exponentially. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

By introducing a new vertex at the midpoint of each edge, each triangle of the control mesh is subdivided into four sub-triangles. The position of new vertices can be calculated from the vertices of previous level as k+1 xiv =

N 5 k 3 ∑ k xiv + x 8 8N i=1 i

(3)

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

5

Fig. 5. A single regular triangular patch defined by 12 control vertices (reft) and an irregular triangular patch defined by N + 6 control vertices (right).

where iv denotes the ith vertex point and the coordinates of new edge points are 3 1 3 1 k+1 xie = x1k + x2k + x3k + x4k 8 8 8 8 where ie means the ith edge point.

(4)

2.2. Subdivision surface basis functions In practice, one cannot conduct unlimited number of subdivisions. An important technical advance in subdivision surfaces is to generate limiting surfaces by constructing an elementwise mapping from the subdivided meshes. In this way, the subdivided meshes play the same role as control grids of NURBS in geometry modeling, and the limiting surface resembles the physical surface in NURBS. Such mapping of subdivision surfaces is formulated as xe (θ1 , θ2 ) =

Na ∑

Ba (θ1 , θ2 )X ae

(5)

a=1

where xe is the Cartesian coordinates of a point in element e on limiting surface, (θ1 , θ2 ) ∈ [0, 1]2 the local coordinates within element e, Ba the box-spline basis functions, and X a the nodal coordinates of the patch in which the element e is located. Here a patch is defined as collection of elements which share at least one vertex with the target element, as shown in Fig. 5. We remark that the term of element in subdivision surfaces is analogous to that in classical FEM/BEM in the sense it provides a local parameterized domain for piecewise mapping to the physical surface. However, the vertices in control meshes are not located on the limiting surface since the box spline basis functions are not interpolatory. Moreover, it is worth noting that the basis functions which have support over the element e are associated with the nodes of its patch. A patch is regular if it is determined by 12 control vertices, i.e Na = 12, and the corresponding basis functions for evaluating Loop subdivision surfaces use bivariate quartic box splines: ) 1 ( 4 B1 (u, v, w) = u + 2u 3 v 12 1 4 B2 (u, v, w) = (u + 2u 3 w) 12 1 4 B3 (u, v, w) = (u + 2u 3 w + 6u 3 v + 6u 2 vw + 12u 2 v 2 + 6uv 2 w + 6uv 3 + 2v 3 w + v 4 ) 12 1 B4 (u, v, w) = (6u 4 + 24u 3 w + 24u 2 w2 + 8uw3 + w4 + 24u 3 v + 60u 2 vw + 36uvw 2 12 + 6vw3 + 24u 2 v 2 + 36uv 2 w + 12v 2 w 2 + 8uv 3 + 6v 3 w + v 4 ) 1 4 B5 (u, v, w) = (u + 6u 3 w + 12u 2 w 2 + 6uw3 + w4 + 2u 3 v + 6u 2 vw + 6uvw 2 + 2vw 3 ) 12 1 B6 (u, v, w) = (2uv 3 + v 4 ) 12

6

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

1 4 (u + 6uw 3 + 12u 2 w 2 + 6uw3 + w 4 + 8u 3 v + 36u 2 vw + 36uvw 2 + 8vw3 12 + 24u 2 v 2 + 60uv 2 w + 24v 2 w 2 + 24uv 3 + 24v 3 w + 6v 4 ) 1 4 (u + 8u 3 w + 24u 2 w 2 + 24uw 3 + 6w4 + 6u 3 v + 36u 2 vw + 60uvw2 B8 (u, v, w) = 12 + 24vw 3 + 12u 2 v 2 + 36uv 2 w + 24v 2 w 2 + 6uv 3 + 8v 3 w + v 4 ) 1 B9 (u, v, w) = (2uw 3 + w 4 ) 12 1 B10 (u, v, w) = (2v 3 w + v 4 ) 12 1 B11 (u, v, w) = (2uw 3 + w 4 + 6uvw 2 + 6vw3 + 6uv 2 w + 12v 2 w2 + 2uv 3 + 6v 3 w + v 4 ) 12 1 4 (w + 2vw 3 ) (6) B12 (u, v, w) = 12 where the barycentric coordinates (u, v, w) are subject to the constraint u + v + w = 1, and (u, v) = (θ1 , θ2 ). The basis functions in Eq. (6) are plotted in Fig. 6. Compared to NURBS basis functions that are defined recursively, the subdivision surface basis functions are polynomials and thus can be evaluated more efficiently. For irregular patches, the initial mesh is subdivided at least once to isolate the extraordinary vertices so that each patch contains one irregular vertex only. In a triangle with irregular vertices, Stam [76] proposed an elegant solution based on the eigendecomposition of the subdivision matrix. It needs several steps of refinement until the point under consideration is located in a patch of regular elements (Fig. 5). B7 (u, v, w) =

3. Isogeometric boundary element method with subdivision surfaces 3.1. Acoustic boundary integral equation Consider a problem of acoustic scattering of an incident wave as shown in Fig. 7. The Helmholtz equation governing the acoustic field in the exterior infinite domain Ω can be transformed to the conventional boundary integral equation (CBIE): ∫ ∫ ∂G(x, y) 1 (7) ∀x ∈ S, p(x) + p(y)dS(y) = G(x, y)q(y)dS(y) + pinc (x) 2 S ∂n(y) S where S denotes the structural surface which is the boundary of the acoustic domain Ω ∈ R3 , and x, y ∈ S are known as source point and field point, respectively. n(y) is the normal at the field point y. p(y) is the acoustic pressure field and its normal flux is q(y) = ∂ p(y)/∂n(y). The total acoustic pressure p(x) consists of the incident acoustic pressure pinc and scattered acoustic pressure pscat , i.e. p(x) = pinc (x) + pscat (x). G(x, y) and ∂G(x,y) are ∂n(y) kernel functions which are written as: eikr 4πr ∂G(x, y) eikr ∂r =− (1 − ikr ) 2 ∂n(y) 4πr ∂n(y) G(x, y) =

(8) (9)

√ where i is the imaginary unit defined as i := −1, k denotes the wave number, and r = |x − y|. It is worth noting that Eq. (7) alone may result in instability of the solution at some fictitious eigen-frequencies. The Burton–Miller method that combines the conventional boundary integral equation (Eq. (7)) and its normal derivative formulation overcomes this difficulty. By differentiating Eq. (7) with respect to the normal direction n(x) at the source point x, we arrive at the following hypersingular boundary integral equation (HBIE) ∫ ∫ 1 ∂ 2 G(x, y) ∂G(x, y) ∂ pinc (x) q(x) + p(y)dS(y) = q(y)dS(y) + (10) 2 ∂n(x) S ∂n(x)∂n(y) S ∂n(x)

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

Fig. 6. Box spline basis functions.

Fig. 7. The exterior acoustic scattering problem.

7

8

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

The linear combination of Eqs. (7) and (10) is expressed as follows ∫ ∫ 1 1 ∂G(x, y) ∂ 2 G(x, y) p(x) + α q(x) + p(y)dS(y) + α p(y)dS(y) 2 2 S ∂n(y) S ∂n(x)∂n(y) ∫ ∫ ∂ pinc (x) ∂G(x, y) q(y)dS(y) + pinc (x) + α (11) = G(x, y)q(y)dS(y) + α ∂n(x) S S ∂n(x) where α is the coupling constant: α = i/k for k > 1, and α = i otherwise. To simulate the sound absorption characteristics of absorbing materials on the structural surface, an impedance boundary condition is introduced as q(y) = ikβ(y) p(y)

(12)

where β(y) denotes the normalized acoustic admittance at the field point y. Substituting Eq. (12) into Eq. (11) yields ∫ ∫ 1 + ikαβ(x) ∂G(x, y) ∂ 2 G(x, y) p(x) + p(y)dS(y) + α p(y)dS(y) 2 S ∂n(y) S ∂n(x)∂n(y) ∫ ∫ ∂ pinc (x) ∂G(x, y) ikβ(y) p(y)dS(y) + pinc (x) + α (13) = G(x, y)ikβ(y) p(y)dS(y) + α ∂n(x) S S ∂n(x) After solving the boundary integral equation, the boundary acoustic field can obtained. Then the pressure at any inspection point x f within the acoustic domain which is away from the structural surface can be calculated as ∫ ∫ ∂G(x f , y) p(x f ) = G(x f , y)q(y)dS(y) − p(y)dS(y) + pinc (x f ) (14) ∂n(y) S S 3.2. Discretization A key step in the present work is to employ the Loop subdivision surfaces basis functions that were used for limit surface generation to discretize the acoustic pressure and the normal flux field in boundary integral equation (Eq. (13)). pe (θ1 , θ2 ) = qe (θ1 , θ2 ) =

Na ∑ a=1 Na ∑

Ba (θ1 , θ2 ) p˜ ae

(15)

Ba (θ1 , θ2 )q˜ae

(16)

a=1

where pe and qe denote the sound pressure and its normal flux at a point with parametric coordinates (θ1 , θ2 ) on the boundary element Se , Ba indicates the box-spline basis functions, and p˜ ae and q˜ae stand for the nodal parameters associated with the sound pressure and normal flux, respectively. We adopt collocation-based approach to generate the linear system of equations. In NURBS based IGA, various collocation schemes have been considered [77], such as Greville abscissae, Demko points and the maxima of B-splines. In the present work, the collocation points are chosen as a series of points on the physical surface mapped from the nodes of the control mesh (Fig. 8), which correspond to the maxima of the subdivision surface basis functions. Such collocation scheme has proven to be successful in SS-IGABEM for acoustic analysis [64]. It would be interesting to make a detailed comparison between different collocation schemes in subdivision surfaces in the future. By making Eq. (13) hold at the selected collocation points, the following system of equations are given as ( ) ∫ Ne ∑ Na ∑ 1 + ikαβ ∂G (xc , y) Ba (y (θ1 , θ2 )) G (xc , y) + α p(xc ) = ikβ e p˜ ae (y (θ1 , θ2 )) dSe 2 ∂n(x) Se e=1 a=1 ( ) ∫ Ne ∑ Na ∑ ∂G(xc , y) ∂ 2 G(xc , y) ∂ pinc (xc ) e − p˜ a (y (θ1 , θ2 )) Ba (y) +α dSe + pinc (x) + α (17) ∂n(y) ∂n(x)∂n(y) ∂n(xc ) Se e=1 a=1 where the subscript c denotes the index of collocation points and Ne the number of elements. The above equation can be written in matrix form as: [H − GY] p˜ = pinc

(18)

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

9

Fig. 8. A regular patch (left) and its associated limit surface (right). Hollow black circles are control vertices and solid black circles are collocation points on the limit surface.

inc ˜ is the column vector where p f = p(x f ) and p inc f = p (x f ). H and G are the coefficient matrices of BEM, p containing nodal parameters discretizing the boundary sound pressures field. pinc is a vector collecting the acoustic incident pressure at the collocation points, and Y is the admittance matrix: ⎡ ⎤ β1 0 ⎢ ⎥ .. (19) Y = ik ⎣ ⎦ .

β Ne

0

After Eq. (18) is solved, the acoustic pressure at an inspection point in the infinite acoustic domain can be calculated by combining Eqs. (12), (14) and (16): ∫ Na Ne ∑ ∑ ( ) e e ikβ p˜ a (y (θ1 , θ2 )) Ba (y (θ1 , θ2 )) G x f , y dSe p(x f ) = −

e=1 a=1 ∫ Ne ∑ Na ∑ p˜ ae (y (θ1 , θ2 )) Se e=1 a=1

Se

Ba (y)

∂G(x f , y) dSe + pinc (x f ) ∂n(y)

(20)

In matrix form, Eq. (20) can be expressed by p f = −[H f − G f Y]p + p inc f

(21)

inc where p f = p(x f ), p inc f = p (x f ), H f is a row vector collected from the second term of the right-hand side of Eq. (20), and G f Y is constructed from the coefficients of the first term on the right-hand side.

4. Topology optimization model with SS-IGABEM 4.1. Definition of the optimization and sensitivity analysis The acoustic topology optimization problem can be formulated as ⎧ 2 ⎪ ⎪min Π = | p f | = p f p f ⎨ ∑ Ne

s.t. ⎪ ⎪ ⎩

e=1 ρe Ae t

≤ V∗ 0 ≤ ρe ≤ 1 V

(22)

where p f in objective function is the sound pressure at the inspection point x f and p f denotes its complex conjugate. The variables in the constraint functions are related to impedance boundary conditions. t is the thickness of the absorbing material layer which is assumed to be a constant around the whole surface. ρe and Ae are the absorbing material density and area of the boundary element S∑ e (e = 1, 2, . . . , Ne ), respectively. V is the total volume of the Ne absorbing materials which can be calculated as V = e=1 Ae t. V ∗ is the constraint of the material volume fraction.

10

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

Since we use a gradient-based optimization algorithm to solve the optimization problem, the sensitivities of the objective function with respect to design variables should be evaluated. We use the adjoint variable method to compute the sensitivity. From Eqs. (18) and (21), the objective function can be expressed as ] [ ] [ (23) Π = Π ( p f ) + λT1 (H − GY)p˜ − pinc +λ2 p f + (H f − G f Y)p˜ − p inc f where Π is a real-valued function, λ1 and λ2 are arbitrarily chosen adjoint variables. By differentiating the above equation with respect to design variable ρe , we arrive at the following equation: [ ] ∂pf ∂pf ∂Π ∂ p˜ ∂Y = pf + pf + λT1 (H − GY) −G p˜ ∂ρe ∂ρe ∂ρe ∂ρe ∂ρe [ ] ∂pf ∂ p˜ ∂Y + λ2 + (H f − G f Y) − Gf p˜ (24) ∂ρe ∂ρe ∂ρe After rearrangement of terms, Eq. (24) can be rewritten as ∂pf ∂Y ∂Y ∂Π = −λT1 G p˜ − λ2 G f p˜ + (λ2 + 2 p f ) ∂ρe ∂ρe ∂ρe ∂ρe [ T ] ∂ p˜ + λ1 (H − GY) + λ2 (H f − G f Y) (25) ∂ρe Because the adjoint variables λ1 and λ2 can be arbitrarily chosen, the derivative of the objective function with respect to design variable ρe can be expressed by ∂Π ∂Y ∂Y = −λT1 G p˜ + λ2 G f p˜ (26) ∂ρe ∂ρe ∂ρe when the following adjoint equations are satisfied { λ2 + 2 p f = 0 (27) λT1 (H − GY) + λ2 (H f − G f Y) = 0 The above equation can be solved iteratively using GMRES. After solving the adjoint equation, λ1 and λ2 can be substituted back into Eq. (26) to evaluate the sensitivities and the terms of λT1 G and λ2 G f can be computed by the fast multipole method (FMM). 4.2. Material interpolation scheme and update of design variables According to the Delany–Bazley–Miki model, the normalized impedance of porous materials is given as ( )−0.632 ( )−0.632 f f + 0.1071i z = 1 + 0.0699 σ σ

(28)

where σ is the flow resistance rate of the material (N s/m4 ), f is the frequency (Hz), and the corresponding normalized admittance value is expressed as 1 β0 = (29) z According to the Solid Isotropic Material with Penalization (SIMP) method, the element admittance can be interpolated based on the following equation: βe = β0 ρeη

(30)

where η is the penalty factor which plays a role in making the intermediate density approach zero or one. The determination of penalty factors is discussed in detail in [78]. The penalty factor η should be greater than 1 so that intermediate densities are unfavorable in optimization. If one specifies the larger value of penalty factor, the intermediate densities decrease and the final optimized model is closer to 0–1 design. However, the increase of η may give a solution which is local minima and highly sensitive to the initial design. In the present work, η is set to be 3 as commonly used in topology optimization literature [2,8]. After the convergence of optimization, we use a Heaviside function as a second filter to eliminate intermediate densities.

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

11

Fig. 9. Flowchart of topology optimization.

In this paper, we employ the method of moving asymptotes (MMA) as the optimizer to update the design variables with convergence criterion ⏐ i+1 ⏐ ⏐Π − Π i⏐ <τ (31) Πi where Π i represents the value of the objective function at the ith iteration. The complete optimization process is shown in Fig. 9.

5. Numerical examples In order to investigate the effectiveness and applicability of the topology optimization approach using SS-IGABEM, several numerical examples are presented in this section. The algorithm has been accelerated with Fast Multipole Method. In the following numerical examples, the geometries contain extraordinary points, and the incident plane acoustic wave is prescribed as pinc = Aeikx·d , where A is the amplitude of incident wave that travels along with direction d with |d| = 1. The iterative convergence criterion τ is set to 1.0 × 10−4 .

12

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

Fig. 10. Control meshes of a sphere model with different subdivision levels.

Fig. 11. The error of acoustic pressure with the increase of degree of freedoms.

5.1. Plane wave scattered by a rigid sphere Firstly we test current algorithm in acoustic analysis by analyzing a problem of acoustic scattering by a rigid sphere with the radius of 1. We begin with a coarse polygon mesh and apply Loop subdivision method to recursively subdivide the control mesh (Fig. 10). It is noted that the repeated subdivision increases the number of degree of freedoms, but the limiting surface remains the same which can be constructed using box spline basis functions from any level of refined control mesh. Hence, the subdivision surfaces enable us to perform h-refinement in exact geometries like in NURBS-based IGA. We investigate the acoustic pressure at point (0, 1, 1) and make a comparison between the numerical results and analytical solutions. As shown in Fig. 11, the relative error of numerical acoustic pressure decreases as the number of elements increases, Fig. 12 shows the solution of sound pressure at the computing point (2, 0, 0). From this figure, we can find that conventional single boundary integral equation (Eq. (7)) results in a deviation from the analytical solution at some fictitious eigen-frequencies. However, the combination of conventional boundary integral equation and hypersingular boundary integral equations with Burton–Miller formulation (Eq. (11)) can overcome this problem and keep numerical results agree well with analytical results at all frequencies. It is underlined that the high accuracy of the numerical solution is achieved even in the presence of extraordinary points.

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

13

Fig. 12. The real part and imaginary part of acoustic pressure at a computing point (2, 0, 0) in terms of frequency by SS-IGABEM using different integral equations.

Fig. 13. Optimization problem of rigid rounded edge cube.

5.2. Plane wave scattered by a rigid rounded edge cube The distribution of sound absorbing materials on the surface of a rounded edge cube under the action of plane waves is investigated, as shown in Fig. 13. Its initial mesh has 386 elements and 768 vertices. The optimization goal is to minimize the sound pressure at a sample point A with coordinate (5, 0, 0). The plane wave propagates forward along the x direction. The origin of the Cartesian coordinate system is located at the center of the cube. The design variable ρe varies from 0 to 1, with zero standing for no sound absorbing material coverage and one denoting full material coverage. The intermediate value between zero and one will be filtered in the final step of the optimization to facilitate industrial manufacturing. The convergence of the objective function and material volume fraction during the optimization process is depicted in Fig. 14. The initial value of the objective function corresponds to the full coverage design. Due to the influence of volume constraint, the volume fraction decreases sharply at the first iterative steps, which leads to the rapid increase of the objective function. After that, the volume ratio becomes stable at 0.5 that satisfies the volume constraint, and the optimization algorithm keeps minimizing the objective function until convergence according to the sensitivity of the objective function. The final optimized distribution of sound absorbing material on the subdivision surfaces of the cube model is shown in Fig. 15. The black elements indicate that sound absorbing materials are applied onto their surfaces, and the white elements denote that the surface of these elements is free of sound absorbing materials. Some gray elements with an intermediate density appear in these optimized distribution figures. The use of sparse triangular grids leads to asymmetric material distribution in the optimized design (Fig. 15a). In comparison, further subdivision can result in symmetrical material distribution in the optimized model (Fig. 15b) with a smooth acoustic pressure and flux field (Fig. 16). The noise reduction level of sound absorbing material is affected by the frequency of noise source, i.e. incident acoustic wave, so the optimized material distribution varies with different frequencies of incidence waves. Fig. 17 shows the optimized material distribution contour at four different frequencies. It can be seen from the figure that

14

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

Fig. 14. Iteration history of the cube model after subdivision once during the optimization process. The frequency of the incident wave is 200 Hz. The cube model after one level of subdivision has 3072 triangular elements and 1538 vertices. The constraint of the volume fraction V ∗ is set to 0.5. The initial values of the design variables ρe are all set to 1.

Fig. 15. The optimized distribution of sound absorbing materials on the cube surfaces with different levels of subdivisions. The figures on the top are viewed from the direction of wave propagation, and that on the bottom are viewed from the opposite direction.

when the excitation frequency is different, the final optimal material distribution are also different, exhibiting the dependence of the acoustic topology optimization results on frequencies. Moreover, when the frequency is relatively low, the sound absorbing material will aggregate over a large area, and with the increase of the frequency, the distribution of the sound absorbing material becomes more dispersed. This is because as the frequency increases, the wavelength becomes smaller, the period of the phase change of the corresponding acoustic wave interference becomes shorter, and the sound absorbing material is distributed in accordance with the interference condition. In addition, the use of subdivision technology makes the optimization results smoother and more symmetric, which is beneficial to industrial manufacturing. 5.3. A block combination model In this section, the model mentioned in introduction section is used for the optimization analysis. The objective function and the material volume fraction in terms of iteration steps are shown in Fig. 18. Similar to Fig. 14, the

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

15

Fig. 16. The contour of the acoustic pressure and flux on the cube surface with optimized distribution of sound absorbing materials. The figures in the first row show the contour of sound pressure, and the second row the flux. The grid with two levels of subdivision is used for the numerical solution.

Fig. 17. The optimized distribution of sound absorbing materials of cube surface at different frequencies. The figures in the first row show the results corresponding to the control meshes with one level of subdivision, and the second row the meshes with two levels of subdivision.

objective function in the initial step is the smallest because of the full coverage of sound absorbing materials. Then, due to the constraint of material volume fraction, the coverage material decreases sharply, resulting in the rapid increase of the objective function. Finally, the use of optimization algorithm makes the objective function approach to the optimized solution stably. The final optimized distribution of sound absorbing materials on the subdivision surface of this initial grid and its acoustic pressure distribution map are shown in Fig. 19. There exist some intermediate density elements in these optimized distribution figures because of the complexity of the structural model. In fact, for complex models, it is difficult to eliminate intermediate density elements completely, but the use of subdivision schemes can effectively reduce the number of intermediate density elements.

16

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

Fig. 18. Iteration history of elastic submarine shell for the optimization problem.

Fig. 19. Material distributions and acoustic pressure at front or behind the plane wave.

Fig. 20. The initial control mesh of submarine model with 9510 elements and 19 016 vertices (left). The smooth limit surface of the submarine model (right).

5.4. Submarine model Now we consider the optimization of the sound absorbing material distributions on the surface of a submarine model (Fig. 20). This example exhibits the capability of the proposed algorithm in handing complicated geometry.

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

17

Fig. 21. Iteration history of elastic submarine shell for the optimization problem.

Fig. 22. The distribution of sound absorbing materials under different iteration steps.

The submarine is immersed in water, and the speed of sound propagation is 1482 m/s. The plane wave propagates along positive z direction. The optimization goal is to minimize the sound pressure magnitude at a point (0,0,6). The volume fraction constraint is set to be V ∗ = 0.5, and the initial values of the design variables are specified as one. The iteration history of the objective function and material volume fraction is shown in Fig. 21. Herein, four different incident wave frequencies, f = 50, 80, 100, and 150 Hz are considered. Similar to Figs. 14 and 18, the objective function value increases rapidly in the initial iteration steps because of the material volume constraints,

18

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

Fig. 23. Variation of acoustic pressure amplitude under different iteration steps.

and then steadily decreases until it converges. The optimization process at four frequencies converges almost at 30 iteration steps. The material volume fraction oscillates slightly in the first 10 iterations and then approximates to 0.5 smoothly. The material distribution on the structural surface in the iteration procedure is presented in Fig. 22, and the corresponding acoustic pressure contours are shown in Fig. 23. From Fig. 22, we can find that a large number of gray elements with intermediate density exist in the middle and early stages of the optimization process, such as for n iter = 5 and n iter = 10. More areas on structural surface for n iter = 10 are filled with red and green than that for n iter = 5, indicating that more elements with artificial density of 0 or 1 appear. It demonstrates that the use of density filter enables the optimization to approach to 0–1 design. The optimization process converges at n iter = 31, and we can hardly observe intermediate density elements. Fig. 23 shows a slight change in acoustic pressure distribution, although accompanied by a remarkable change in material density. Additionally, the acoustic pressure on the structural surface behind the plane wave is much smaller than that in front of the plane wave because the transmission of sound wave is scattered by the structure. Three different frequencies, f = 50, 80, and 150 Hz are considered for optimization analysis, and the final optimized material distribution on the structural surface is shown in Fig. 24. The material distributions at different frequencies are also different, so the optimized results are frequency dependent. In order to satisfy the accuracy requirement, the topology optimization at frequency band requires different scales of discrete meshes. The higher the frequency is, the denser the mesh is required. The use of subdivision technology can refine meshes to meet discretization requirements at high frequencies, but maintain the same geometry. It eliminates the costly mesh regeneration from the original CAD models.

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

19

Fig. 24. The optimal distribution of sound absorbing materials under different frequencies.

6. Conclusion Acoustic topology optimization of sound absorbing materials with isogeometric boundary element methods is presented. The subdivision surfaces are used to construct watertight CAD models with complicated geometries, and the associated box-splines basis functions are adopted to construct the physical surface and discretize acoustic boundary integral equations. The present method enables acoustic topology optimization to be conducted directly from CAD models without any meshing procedures and meanwhile eliminating geometric errors. The work can be extended to topology optimization of acoustic-structure interaction problems in future. In addition, other optimization algorithms for example level-set methods and deformable simplicial complex methods [79–81] can be incorporated. Because excitation frequency is commonly in a frequency band rather than a single value, it is meaningful to investigate the average optimized distribution results by performing acoustic topology optimization for incident waves in a frequency range. Our method is currently limited to dealing with C 1 continuous surfaces. We will incorporate the techniques which are able to handle sharp edges and corners (C 0 continuities) in the future work. Acknowledgments Financial supports from the National Natural Science Foundation of China (11702238, 51904202, 11901578), and Nanhu Scholars Program for Young Scholars of XYNU are acknowledged. We thank the financial support of Intuitive modeling and SIMulation platform (IntuiSIM) (PoC17/12253887) grant by Luxembourg National Research Fund. Haibo Chen appreciates the support of the National Natural Science Foundation of China (NSFC) under Grant No. 11772322. References [1] [2] [3] [4]

S. Marburg, Developments in structural-acoustic optimization for passive noise control, Arch. Comput. Methods Eng. 27 (2002) 291–370. M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods and Applications, Springer, 2003. M.B. Dühring, J.S. Jensen, O. Sigmund, Acoustic design by topology optimization, J. Sound Vib. 317 (3–5) (2008) 557–575. J. Du, N. Olhoff, Minimization of sound radiation from vibrating bi-material structures using topology optimization, Struct. Multidiscip. Optim. 33 (2007) 305–321.

20

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

[5] G.H. Yoon, J.S. Jensen, O. Sigmund, Topology optimization of acoustic-structure interaction problems using a mixed finite element formulation, Internat. J. Numer. Methods Engrg. 70 (2007) 1049–1075. [6] Z. Kang, X. Zhang, S. Jiang, G. Cheng, On topology optimization of damping layer in shell structures under harmonic excitations, Struct. Multidiscip. Optim. 46 (2012) 51–67. [7] J. Zhu, F. He, T. Liu, W. Zhang, Q. Liu, C. Yang, Structural topology optimization under harmonic base acceleration excitations, Struct. Multidiscip. Optim. 57 (2018) 1061–1078. [8] W. Zhao, L. Chen, H. Chen, S. Marburg, Topology optimization of exterior acoustic-structure interaction systems using the coupled FEM-BEM method, Internat. J. Numer. Methods Engrg. 119 (2019) 1–28. [9] L. Chen, C. Zheng, H. Chen, S. Marburg, Structural-acoustic sensitivity analysis of radiated sound power using a finite element/ discontinuous fast multipole boundary element scheme, Internat. J. Numer. Methods Fluids 82 (2016) 858–878. [10] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (39–41) (2005) 4135–4195. [11] H.K. Seo, Yudeok, S. Youn, Isogeometric topology optimization using trimmed spline surfaces, Comput. Methods Appl. Mech. Engrg. 199 (49–52) (2010) 3270–3296. [12] L. Dedè, M.J. Borden, T.J.R. Hughes, Isogeometric analysis for topology optimization with a phase field model, Arch. Comput. Methods Eng. 19 (3) (2012) 427–465. [13] X. Qian, Topology optimization in B-spline space, Comput. Methods Appl. Mech. Engrg. 265 (2013) 15–35. [14] E. Atroshchenko, S. Tomar, G. Xu, S.P.A. Bordas, Weakening the tight coupling between geometry and simulation in isogeometric analysis: From sub-and super-geometric analysis to geometry-independent field approximation (GIFT), Internat. J. Numer. Methods Engrg. 114 (10) (2018) 1131–1159. [15] R.N. Simpson, S.P.A. Bordas, J. Trevelyan, T. Rabczuk, A two-dimensional isogeometric boundary element method for elastostatic analysis, Comput. Methods Appl. Mech. Engrg. 209 (2012) 87–100. [16] R.N. Simpson, S.P.A. Bordas, H. Lian, J. Trevelyan, An isogeometric boundary element method for elastostatic analysis: 2D implementation aspects, Comput. Struct. 118 (2013) 2–12. [17] M.A. Scott, R.N. Simpson, J.A. Evans, S. Lipton, S.P.A. Bordas, T.J.R. Hughes, T.W. Sederberg, Isogeometric boundary element analysis using unstructured T-splines, Comput. Methods Appl. Mech. Engrg. 254 (2013) 197–221. [18] X. Peng, E. Atroshchenko, P. Kerfriden, S.P.A. Bordas, Linear elastic fracture simulation directly from CAD: 2D NURBS-based implementation and role of tip enrichment, Int. J. Fract. 204 (2017) 55–78. [19] G. Xu, B. Mourrain, R. Duvigneau, A. Galligo, Parameterization of computational domain in isogeometric analysis: methods and comparison, Comput. Methods Appl. Mech. Engrg. 200 (23–24) (2011) 2021–2031. [20] G. Xu, B. Mourrain, R. Duvigneau, A. Galligo, Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications, Comput. Aided Des. 45 (2) (2013) 395–404. [21] G. Xu, M. Li, B. Mourrain, T. Rabczuk, J. Xu, S.P.A. Bordas, Constructing iga-suitable planar parameterization from complex cad boundary by domain partition and global/local optimization, Comput. Methods Appl. Mech. Engrg. 328 (2018) 175–200. [22] K. Li, X. Qian, Isogeometric analysis and shape optimization via boundary integral, Comput. Aided Des. 43 (11) (2011) 1427–1437. [23] K. Kostas, A. Ginnis, C. Politis, P. Kaklis, Ship-hull shape optimization with a T-spline based BEM–isogeometric solver, Comput. Methods Appl. Mech. Engrg. 284 (2015) 611–622. [24] H. Lian, P. Kerfriden, S.P.A. Bordas, Implementation of regularized isogeometric boundary element methods for gradient-based shape optimization in two-dimensional linear elasticity, Internat. J. Numer. Methods Engrg. 106 (12) (2016) 972–1017. [25] H. Lian, P. Kerfriden, S.P.A. Bordas, Shape optimization directly from cad: an isogeometric boundary element approach using T-splines, Comput. Methods Appl. Mech. Engrg. 317 (2017) 1–41. [26] S. Li, J. Trevelyan, Z. Wu, H. Lian, D. Wang, W. Zhang, An adaptive SVD–Krylov reduced order model for surrogate based structural shape optimization through isogeometric boundary element method, Comput. Methods Appl. Mech. Engrg. 349 (2019) 312–338. [27] R.N. Simpson, M.A. Scott, M. Taus, D.C. Thomas, H. Lian, Acoustic isogeometric boundary element analysis, Comput. Methods Appl. Mech. Engrg. 269 (2014) 265–290. [28] L. Chen, L. Liu, W. Zhao, C. Liu, An isogeometric approach of two dimensional acoustic design sensitivity analysis and topology optimization analysis for absorbing material distribution, Comput. Methods Appl. Mech. Engrg. 336 (2018) 507–532. [29] R.N. Simpson, Z. Liu, Acceleration of isogeometric boundary element analysis through a black-box fast multipole method, Eng. Anal. Bound. Elem. 66 (2016) 168–182. [30] X. Peng, E. Atroshchenko, P. Kerfriden, S.P.A. Bordas, Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth, Comput. Methods Appl. Mech. Engrg. 316 (2017) 151–185. [31] D. Schillinger, J.A. Evans, A. Reali, M.A. Scott, T.J.R. Hughes, Isogeometric collocation: Cost comparison with galerkin methods and extension to adaptive hierarchical NURBS discretizations, Comput. Methods Appl. Mech. Engrg. 267 (2013) 170–232. [32] C. Liu, L. Chen, W. Zhao, H. Chen, Shape optimization of sound barrier using an isogeometric fast multipole boundary element method in two dimensions, Eng. Anal. Bound. Elem. 85 (2017) 142–157. [33] L. Coox, F. Greco, O. Atak, D. Vandepitte, W. Desmet, A robust patch coupling method for NURBS-based isogeometric analysis of non-conforming multipatch surfaces, Comput. Methods Appl. Mech. Engrg. 316 (2017) 235–260. [34] A. Vuong, C. Giannelli, B. Juttler, B. Simeon, A hierarchical approach to adaptive local refinement in isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 200 (49) (2011) 3554–3567. [35] V.P. Nguyen, P. Kerfriden, M. Brino, S.P.A. Bordas, E. Bonisoli, Nitsche’s method for two and three dimensional NURBS patch coupling, Comput. Mech. 53 (6) (2014) 1163–1182. [36] V.P. Nguyen, P. Kerfriden, S.P.A. Bordas, Two-and three-dimensional isogeometric cohesive elements for composite delamination analysis, Composites B 60 (2014) 193–212.

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

21

[37] Q. Hu, F. Chouly, P. Hu, G. Cheng, S.P.A. Bordas, Skew-symmetric nitsche’s formulation in isogeometric analysis: Dirichlet and symmetry conditions, patch coupling and frictionless contact, Comput. Methods Appl. Mech. Engrg. 341 (2018) 188–220. [38] T.W. Sederberg, J. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCs, ACM Trans. Graph. 22 (3) (2003) 477–484. [39] Y. Bazilevs, V.M. Calo, J. Cottrell, J.A. Evans, T.J.R. Hughes, S. Lipton, M.A. Scott, T.W. Sederberg, Isogeometric analysis using T-splines, Comput. Methods Appl. Mech. Engrg. 199 (5) (2010) 229–263. [40] X. Li, F. Chen, H. Kang, J. Deng, A survey on the local refinable splines, Sci. China Math. 59 (4) (2016) 617–644. [41] X. Li, J. Deng, F. Chen, Surface modeling with polynomial splines over hierarchical T-meshes, Vis. Comput. 23 (12) (2007) 1027–1033. [42] J. Deng, F. Chen, X. Li, C. Hu, W. Tong, Z. Yang, Y. Feng, Polynomial splines over hierarchical T-meshes, Graph. Models 70 (4) (2008) 76–86. [43] X. Li, J. Deng, F. Chen, Polynomial splines over general T-meshes, Vis. Comput. 26 (4) (2010) 277–286. [44] N. Nguyenthanh, H. Nguyenxuan, S.P.A. Bordas, T. Rabczuk, Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids, Comput. Methods Appl. Mech. Engrg. 200 (21) (2011) 1892–1908. [45] N. Nguyenthanh, J. Kiendl, H. Nguyenxuan, R. Wuchner, K.U. Bletzinger, Y. Bazilevs, T. Rabczuk, Rotation free isogeometric thin shell analysis using PHT-splines, Comput. Methods Appl. Mech. Engrg. 200 (47) (2011) 3410–3424. [46] C. Anitescu, N. Hossain, T. Rabczuk, Recovery-based error estimation and adaptivity using high-order splines over hierarchical T-meshes, Comput. Methods Appl. Mech. Engrg. 328 (2018) 638–662. [47] D.R. Forsey, R.H. Bartels, Hierarchical B-spline refinement, ACM Siggraph Comput. Graph. 22 (4) (1988) 205–212. [48] C. Hofreither, B. Jüttler, G. Kiss, W. Zulehner, Multigrid methods for isogeometric analysis with THB-splines, Comput. Methods Appl. Mech. Engrg. 308 (2016) 96–112. [49] T. Dokken, T. Lyche, K.F. Pettersen, Polynomial splines over locally refined box-partitions, Comput. Aided Geom. Design 30 (3) (2013) 331–356. [50] K.A. Johannessen, T. Kvamsdal, T. Dokken, Isogeometric analysis using LR B-splines, Comput. Methods Appl. Mech. Engrg. 269 (2014) 471–514. [51] F. Cirak, Q. Long, Subdivision shells with exact boundary control and non-manifold geometry, Internat. J. Numer. Methods Engrg. 88 (9) (2011) 897–923. [52] E.E. Catmull, J.H. Clark, Recursively generated B-spline surfaces on arbitrary topological meshes, Comput. Aided Des. 10 (1978) 350–355. [53] D. Doo, M.A. Sabin, Behaviour of recursive division surfaces near extraordinary points, Comput. Aided Des. 10 (6) (1978) 356–360. [54] C. Loop, Smooth Subdivision Surfaces Based on Triangles (Master’s thesis), Department of Mathematics, University of Utah, 1987. [55] Z. Huang, J. Deng, G. Wang, A bound on the approximation of a Catmull-Clark subdivision surface by its limit mesh, Comput. Aided Geom. Design 25 (2008) 457–569. [56] D. Toshniwal, H. Speleers, T.J.R. Hughes, Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations, Comput. Methods Appl. Mech. Engrg. 327 (2017) 411–458. [57] M. Wu, B. Mourrain, A. Galligo, B. Nkonga, Hermite type spline spaces over rectangular meshes with arbitrary topology, Commun. Comput. Phys. 21 (2017) 835–866. [58] C. Lee, Automatic metric 3D surface mesh generation using subdivision surface geometrical model. Part 2: Mesh generation algorithm and examples, Int. J. Numer. Methods Eng. 56 (2003) 1615–1646. [59] Q. Pan, G. Xu, Y. Zhang, A unified method for hybrid subdivision surface design using geometric partial differential equations, Comput. Aided Des. 46 (2014) 110–119. [60] X. Wei, Y. Zhang, T.J.R. Hughes, M.A. Scott, Truncated hierarchical Catmull-Clark subdivision with local refinement, Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20. [61] Q. Pan, G. Xu, G. Xu, Y. Zhang, Isogeometric analysis based on extended Catmull-Clark subdivision, Comput. Math. Appl. 71 (2016) 105–119. [62] F. Cirak, M. Ortiz, P. Schröder, Subdivision surfaces: a new paradigm for thin-shell finite-element analysis, Internat. J. Numer. Methods Engrg. 47 (12) (2000) 2039–2072. [63] F. Cirak, M.J. Scott, E.K. Antonsson, M. Ortiz, P. Schröder, Integrated modeling, finite-element analysis, and engineering design for thin-shell structures using subdivision, Comput. Aided Des. 34 (2) (2002) 137–148. [64] Z. Liu, M. Majeed, F. Cirak, R.N. Simpson, Isogeometric FEM-BEM coupled structural-acoustic analysis of shells using subdivision surfaces, Internat. J. Numer. Methods Engrg. 113 (9) (2018) 1507–1530. [65] K. Bandara, F. Cirak, O. Steinbach, J. Zapletal, Boundary element based multiresolution shape optimisation in electrostatics, J. Comput. Phys. 297 (2015) 584–598. [66] J. Li, D. Dault, B. Liu, Y. Tong, B. Shanker, Subdivision based isogeometric analysis technique for electric field integral equations for simply connected structures, J. Comput. Phys. 319 (2016) 145–162. [67] K. Bandara, F. Cirak, Isogeometric shape optimisation of shell structures using multiresolution subdivision surfaces, J. Phys. Conf. Ser. 734 (3) (2016) 032142. [68] K. Bandara, T. Rüberg, F. Cirak, Shape optimisation with multiresolution subdivision surfaces and immersed finite elements, Comput. Methods Appl. Mech. Engrg. 300 (2016) 510–539. [69] A. Burton, G. Miller, The application of integral equation methods to the numerical solution of some exterior boundary-value problems, Proc. R. Soc. A 323 (1553) (1971) 201–210. [70] R. Troian, F. Gillot, S. Besset, Adjoint sensitivity related to geometric parameters for mid-high frequency range vibroacoustics, Struct. Multidiscip. Optim. 52 (4) (2015) 803–811.

22

L. Chen, C. Lu, H. Lian et al. / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112806

[71] H. Denli, J. Sun, Structural-acoustic optimization of sandwich cylindrical shells for minimum interior sound transmission, J. Sound Vib. 316 (1) (2008) 32–49. [72] K. Koo, B. Pluymers, W. Desmet, S. Wang, Vibro-acoustic design sensitivity analysis using the wave-based method, J. Sound Vib. 330 (17) (2011) 4340–4351. [73] L. Chen, S. Marburg, H. Chen, H. Zhang, H. Gao, An adjoint operator approach for sensitivity analysis of radiated sound power in fully coupled structural-acoustic systems, J. Comput. Acoust. 25 (01) (2017) 1750003. [74] C. Zheng, T. Matsumoto, T. Takahashi, H. Chen, A wideband fast multipole boundary element method for three dimensional acoustic shape sensitivity analysis based on direct differentiation method, Eng. Anal. Bound. Elem. 36 (3) (2012) 361–371. [75] L. Chen, H. Lian, Z. Liu, H. Chen, E. Atroshchenko, S.P.A. Bordas, Structural shape optimization of three dimensional acoustic problems with isogeometric boundary element methods, Comput. Methods Appl. Mech. Engrg. 355 (2019) 926–951. [76] J. Stam, Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values, in: Siggraph, Vol. 98, Citeseer, 1998, pp. 395–404. [77] D. Schillinger, J.A. Evans, A. Reali, M.A. Scott, T.J.R. Hughes, Isogeometric collocation: Cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations, Comput. Methods Appl. Mech. Engrg. 267 (2013) 170–232. [78] W. Zhao, C. Zheng, C. Liu, H. Chen, Minimization of sound radiation in fully coupled structural-acoustic systems using FEM-BEM based topology optimization, Struct. Multidiscip. Optim. 58 (1) (2017) 115–128. [79] A.N. Christiansen, M. Nobel-Jørgensen, N. Aage, O. Sigmund, J.A. Bærentzen, Topology optimization using an explicit interface representation, Struct. Multidiscip. Optim. 49 (3) (2014) 387–399. [80] H. Lian, A.N. Christiansen, D.A. Tortorelli, O. Sigmund, N. Aage, Combined shape and topology optimization for minimization of maximal von Mises stress, Struct. Multidiscip. Optim. 55 (5) (2017) 1541–1557. [81] M. Zhou, H. Lian, O. Sigmund, N. Aage, Shape morphing and topology optimization of fluid channels by explicit boundary tracking, Internat. J. Numer. Methods Fluids 88 (6) (2018) 296–313.