An isogeometric approach of two dimensional acoustic design sensitivity analysis and topology optimization analysis for absorbing material distribution

Accepted Manuscript An isogeometric approach of two dimensional acoustic design sensitivity analysis and topology optimization analysis for absorbing material distribution Leilei Chen, Cheng Liu, Wenchang Zhao, Linchao Liu

PII: DOI: Reference:

S0045-7825(18)30149-X https://doi.org/10.1016/j.cma.2018.03.025 CMA 11831

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 6 December 2017 Revised date : 15 March 2018 Accepted date : 16 March 2018 Please cite this article as: L. Chen, C. Liu, W. Zhao, L. Liu, An isogeometric approach of two dimensional acoustic design sensitivity analysis and topology optimization analysis for absorbing material distribution, Comput. Methods Appl. Mech. Engrg. (2018), https://doi.org/10.1016/j.cma.2018.03.025 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights (for review)

Research Highlights 1. Applying wideband fast multipole method to isogeometric BEM based on Burton-Miller with non-singular boundary integral. 2. Structural shape optimization based on isogeometric BEM and OC method for acoustic problems. 3. Absorbing material distribution optimization based on BEM and SIMP method.

*Manuscript Click here to download Manuscript: manuscript.pdf

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked References

An isogeometric approach of two dimensional acoustic design sensitivity analysis and topology optimization analysis for absorbing material distribution Leilei Chena,∗, Cheng Liub , Wenchang Zhaob , Linchao Liua a College

of Architecture and Civil Engineering, Xinyang Normal University, P.R.China of Modern Mechanics, University of Science and Technology of China

b Department

Abstract A study of structural shape optimization and absorbing material distribution optimization of noise barrier structures based on the recently proposed isogeometric analysis method with exact geometric definitions is presented. The acoustic scattering is approximated using the basis functions that represent the geometry. A fast multipole method is applied to accelerate the solution of the boundary element method (BEM). The Burton–Miller formulation is used to overcome the fictitious frequency problem when a single Helmholtz boundary integral equation is used for the exterior boundary–value problem. The strongly singular integrals in the Burton–Miller formulation using the isogeometric BEM are evaluated explicitly and directly, particularly for the sensitivity formulation with a hyper–singular integral. The optimality criteria method is used for two types of optimization analysis, shape optimization and material distribution topology optimization. For the shape optimization, the design variables can be set to the locations of the control points because the control points determine the shape of structure. For the second optimization, a new material interpolation scheme for acoustic problems based on the solid isotropic material with penalization (SIMP) method is given, where the interpolation variable is not the real structural density used in a conventional SIMP, but a fictitious material density that determines the normalized surface admittance. Several examples are presented to demonstrate the validity and efficiency of the proposed algorithm. Keywords: Isogeometric analysis, Fast multipole boundary element method, Shape sensitivity analysis, Topology optimization analysis

1. Introduction The concept of isogeometric analysis (IGA) was proposed in a seminal paper by Hughes et al. (1). It successfully integrates computer–aided design (CAD) tools with numerical methods to increase the efficiency of current engineering workflows. In a conventional engineering workflow, an initial geometric model with material properties is designed using commercial CAD ∗ Corresponding

author Email address: [email protected] (Leilei Chen) Preprint submitted to CMAME

March 15, 2018

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

software. Then, the engineers need to generate an analysis mesh according to the designed CAD model. The mesh is used as a discretization method for numerical analysis. The numerical results guide design updates of the material properties and geometric shape. In this conventional process, engineers waste large amounts of time on data conversion. The mesh–generating stage normally occupies 80% of the time in the entire engineering workflow (1). IGA directly adopts the CAD design models for analysis. The geometry of the design model is always described by non-uniform rational B-splines (NURBS), whereas conventional analytical models use polynomial functions to interpolate the geometry and unknown variables. Use of the NURBS function as an analysis function can completely eliminate the meshing stage. Because it integrates design and analysis into one unified model, huge amounts of time will be saved. After the concept was introduced, IGA expanded rapidly to a number of engineering applications, including potential problems (2; 3; 4), elasticity (5; 6; 7; 8; 9; 10), electromagnetics (11; 12), fluid flow (13; 14; 15), and vibrations (16). IGA is focused mainly on the use of the finite element method because it is the most commonly used numerical method. However, finite element analysis requires volumetric discretization, whereas commercial CAD software preferentially uses a surface representation approach to model the geometry. The boundary element method (BEM) has gradually drawn attention because it requires only boundary discretization, making it a perfect match for IGA. Researchers (7) have developed a two–dimensional (2D) isogeometric BEM to solve elastostatic problems. As the BEM has advantages for problems with infinite domains, isogeometric BEM has also been extended to acoustic applications for both 2D (17) and three–dimensional (3D) (18; 19; 20) problems. Design sensitivity analysis can provide the derivatives of quantities of interest with respect to the design variables. It is a critical step in gradient–based optimization methods and provides guidance regarding the variable changes required to achieve higher acoustic performance. Thus, design sensitivity analysis is also important for acoustic optimization problems (21; 22; 23; 24; 25; 26; 27). We expect IGA to have great potential for application in the field of optimization because a unified model is used for both design and analysis. In conventional optimization schemes, a mesh is generated from a CAD model, and after each analysis step, a new geometric model must be generated; thus, a new model is required for each iteration of the optimization process, and the procedure is time-consuming and inefficient. Instead, by using the concept of IGA, in which the same model is used for CAD and analysis, significant gains are made in time and efficiency. The IGA approach eliminates the expensive mesh generation step, thus offering significant efficiency gains over traditional optimization approaches. In sum, IGA has shown better performance, particularly for gradient-based shape optimization analysis, it has several advantages. First, it provides more accurate sensitivity of objective function for complex geometry because of higher order geometric approximation, such as normal vector and curvature. Second, it makes the design and modification of the complex geometry more flexible without communicating with the CAD description of geometry during optimization process. Many studies of IGA optimization have used the finite element method because it is the most commonly used numerical method for optimization. Few works have adopted the isogeometric BEM for optimization, including (28; 29; 30; 31; 32). Previous studies have focused on application to the elastostatic problem. No works have studied application to acoustics. In this paper, we propose an isogeometric approach to the optimization of 2D acoustic problems. Two types of optimization analysis are presented in this paper. One is shape optimization of the geometry, and the other is topology optimization of the distribution of the absorbing material on a noise barrier. For the shape optimization, NURBS control points or weights are used as design variables to flexibly control the boundary shape, and the optimality crite2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

ria (OC) method is applied to update the design variable. The solid isotropic material with penalization (SIMP) method has been widely applied to a series of structural-acoustic topology optimization problems because it has high computing efficiency and can be simply implemented (33; 34; 35; 36; 37; 38; 39; 40; 41; 42). In fact, the IGA method has been extensively used in a phase field method and SIMP method for spatial domain discretization. Also, a method of density based topology optimization was presented using the design space limited to the B-spline space (43). In conventional element based methods, the material density function is constant in each finite element or boundary element. The jogged boundaries and checkerboarding patterns are the well-known drawbacks of such methods especially when a coarse mesh and lower order elements are used. However, the IGA method has been successfully used to overcome this difficulty; see (43), where density function determines the distribution of material throughout the domain of interest and is approximated by NURBS basis functions. In this paper, similar method is applied for topology optimization analysis of absorbing material distribution, where the interpolation variable is not the real structural density used in a conventional SIMP, but a fictitious material density that determines the normalized surface admittance. For the topology optimization analysis, an optimization method generated by using the BEM with a fast multipole method (FMM) and SIMP is proposed for optimization of the distribution of the absorbing material on the noise barrier’s edge. Finally, several examples validate the efficiency of the algorithm proposed in this paper. 2. 2D acoustic state analysis based on BEM 2.1. 2D acoustic state analysis based on conventional BEM The boundary integral equation (BIE) and normal derivative boundary integral equation (NDBIE) of the Helmholtz half-space problem with an absorbing material can be expressed as (44) # Z " ∂G(r s , r) − ikβ(r s )G(r s , r) p(r s , r0 )ds(r s ) = G(r, r0 ), (1) c(r)p(r, r0 ) + ∂n(r s ) S and ikβc(r)p(r, r0 ) +

Z " S

# ∂G(r s , r) ∂G(r, r0 ) ∂2G(r s , r) − ikβ(r s ) p(r s , r0 )ds(r s ) = , ∂n(r s )∂n(r) ∂n(r) ∂n(r)

(2)

where c(r) is a jump term related to the geometry at the source point r, r s denotes the field point on the boundary, r0 = (x0 , y0 ) is the incident acoustic source point outside the boundary, and β is the acoustic admittance. The Green function G(r, r0 ) represents the incident acoustic pressure and can be expressed as i i h (1) H0 (k|r0 − r|) + H0(1) (k|r¯0 − r|) 4

(3)

CBIE + αNDBIE = 0,

(4)

G(r, r0 ) = −

where r¯0 = (x0 , −y0 ) is the mirror point of r0 . The Burton–Miller formulation obtained using a linear combination of the conventional boundary integral equation (1) and NDBIE (2) can give a unique solution at all frequencies (45; 46), as follows:

3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where α is the coupling constant, when k > 1, α = i/k, and α = i for k < 1. Assuming that the boundary S is discretized into Ne boundary elements, the following system of linear algebraic equations can be obtained: [H − ikGB] p = pi

(5)

where the boundary condition for solving the above system of equation is matrix B is given as    β1    . .  B =  .   βn

∂p ∂n

= ikβp, and the

(6)

To simulate the pressure attenuation for different distributions of acoustic materials, a density and surface admittance should be assigned to each boundary element. Thus, βi denotes the normalized surface admittance of the i-th boundary element. S x represents the discretized element containing the source point r, so the boundary integral on S x is singular and requires special treatment. Herein, the Cauchy principal value and Hadamard finite part integral method (47) are used to solve the singular boundary integrals in Eqs. (1) and (2), as follows: Z Z Z   G(r, r s )p(r s ) dS (r s ) = G(ξ)J(ξ) − G0 (ξ)J0 (ξ) p(ξ) dξ + G0 (ξ)J0 (ξ)p(ξ) dξ, (7) Sx

Z

Sx

Sx

∂2G(r, r s ) p(r s ) dS (r s ) = ∂n(r)n(r s )

Z h Sx

Sx

i

F 1 (ξ)J(ξ) − F01 (ξ)J0 (ξ) p(ξ) dξ +

Z

Sx

F01 (ξ)J0 (ξ)p(ξ) dξ,

(8)

where ξ denotes local coordinate in every element, J denotes the Jacobian, and p(ξ) =

n X

φi pi ,

i=1

G0 (ξ) = lim G(r, r s ) = − r s →r

1 ln (kl0 ) , 2π

∂2G(r, r s ) , ∂n(r)n(r s ) " # k2 1 ik (1) F01 (r, r s ) = lim H1 (kl)n j (r)n j (r s ) = − ln (kl0 ) , 2 r s →r 4l 2kl0 4π F 1 (r, r s ) =

J0 (ξ) = lim J(ξ) , r s →r

l2 = (ξ − b)2 J ,

(9)

l02 = (ξ − b)2 J0 ,

where b denotes the local coordinate of the source point r located on the boundary element S x . φi represents the interpolation functions, pi is the nodal pressure, and Hn(1) is the n-th order Hankel function of the first kind. Ultimately, the singular parts of Eqs. (7) and (8) can be rewritten as follows: Z m X G0 (ξ)J0 (ξ)p(ξ) dξ = Bi pi , Sx

Z

Sx

i=1

F01 (ξ)J0 (ξ)p(ξ) dξ 4

=

m X i=1

(10)

Di pi ,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where m denotes the number of interpolation nodes inside the boundary element. When Lagrangian function interpolation is used to obtain the numerical solution, pi denotes the sound pressure at the corresponding collocation point, which lies on the boundary curve. However, for NURBS basis functions, it denotes the pressure coefficient at the corresponding control point, which does not necessarily lie on the boundary curve. By substituting Eq. (9) into Eq. (10), we obtain the following formulation: Z Z 1 J0 1 J0 φi ln |ξ − b| dξ , ln(J0 k) φi dξ − 2π 2π −1 −1 Z Z 1 Z 1 k2 J0 1 φi k2 J0 1 Di = − φi ln |ξ − b| dξ + dξ . ln(J0 k) φi dξ − 4π 4π 2πJ (ξ − b)2 0 −1 −1 −1 Bi = −

(11)

Then, we obtain the exact expressions of Eq. (11) by using the Cauchy principal value and Hadamard finite part integral method for different boundary element discretizations. Using different types of boundary elements based on Lagrangian function interpolation to discretize the boundary, we obtain different expressions for the coefficients Bi and Di ; see (25) for details. Herein, the boundary element with NURBS basis function interpolation is used to discretize the boundary, and the procedure is given in this section. The NURBS basis function can be written as follows: ¯ i N (ξ)w ¯ = Ri,p (ξ) ¯ = i,p φi (ξ) (12) ¯ W(ξ) where wi represents a weight associated with the control point, which may or may not lie on the ¯ is defined as boundary curve, and W(ξ) ¯ = W(ξ)

n X

¯ wi Ni,p (ξ)

(13)

i=0

where ξ¯ is a parameter, that is, ξ¯ ∈ [0, 1] along the boundary curve generated using NURBS basis function interpolation, and can be expressed as a function of the local coordinate parameter ξ, as follows 2ξ¯ − ξ¯e − ξ¯e+1 ξ= (14) ξ¯e+1 − ξ¯e where ξ¯ ∈ [ξ¯e , ξ¯e+1 ], ξ¯e and ξ¯e+1 refer to the parameter values of two endpoints on boundary ¯ in Eq. (12) is the B-spline basis function, which defined as element S x . Ni,p (ξ)

then, for p = 1, 2, 3, . . . ¯ = Ni,p (ξ)

  1 ¯ Ni,0 (ξ) =  0

if

ξ¯i ≤ ξ < ξ¯i+1

otherwise,

ξ¯ − ξ¯ ξ¯ − ξ¯i ¯ + i+p+1 ¯ N (ξ) N (ξ). ¯ξi+p − ξ¯i i,p−1 ¯ξi+p+1 − ξ¯i+1 i+1,p−1

(15)

(16)

where p denotes the order of the curve. More content about B-spline and NURBS basis functions can be found in (1; 14). By substituting Eqs. (12) and (14) into Eq. (11), Bi and Di based on 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

NURBS basis function expressed as: ! 2ξ¯ − ξ¯e − ξ¯e+1  ln(J0 k) Bi = −  Ri,p dξ¯ − ξ¯e+1 − ξ¯e ξ¯e π ξ¯e+1 − ξ¯e ! Z ξ¯e+1 2ξ¯ − ξ¯e − ξ¯e+1 J0 2ξ¯ − ξ¯e − ξ¯e+1   − b dξ¯ , Ri,p ln ξ¯e+1 − ξ¯e ξ¯e+1 − ξ¯e ξ¯e π ξ¯ − ξ¯ Z

J0

e+1

ξ¯e+1

(17)

e

! Z ξ¯e+1 k 2 J0 2ξ¯ − ξ¯e − ξ¯e+1   ln(J0 k) Di = − Ri,p dξ¯ − ξ¯e+1 − ξ¯e ξ¯e 2π ξ¯e+1 − ξ¯e ! Z ξ¯e+1 2ξ¯ − ξ¯e − ξ¯e+1 2ξ¯ − ξ¯e − ξ¯e+1 k 2 J0   Ri,p ln − b dξ¯ + ¯ ¯ ¯ ¯ ¯ ¯ ¯ ξ − ξ ξ − ξ ξ e+1 e e+1 e 2π ξe+1 − ξe e ! ! Z ξ¯e+1 −2 1 2ξ¯ − ξ¯e − ξ¯e+1 2ξ¯ − ξ¯e − ξ¯e+1   Ri,p −b dξ¯ . ξ¯e+1 − ξ¯e ξ¯e+1 − ξ¯e πJ0 ξ¯e+1 − ξ¯e ξ¯e

(18)

2.2. 2D acoustic state analysis based on high-frequency FMBEM The FMM has been widely applied for accelerated solution of the conventional BEM (23; 24; 25; 26). Herein, the key point is the use of a high-frequency FMM to accelerate the solution of the isogeometric BEM for 2D half-space acoustic problems with an absorbing material; the technique is called the IGA-FMBEM. The first step in implementing the IGA-FMBEM is to expand the Green function by plane wave expansion, as follows: ! !# I " −−→ −−−→ −−→ −− → −−1−→ i b 1 ikb k·r1 r 1 1 ikb k·rm1 rm 1 G(r, r s ) = e T θ, r s r + e T θ, r s rm e−ikk·rs rs dθ, (19) 8π where r1 is near from point r, r1s is near from r s , and rm stands for the mirror point of r. b k(θ) is expressed as b k(θ) = (cos θ, sin θ),

and

! X ! ∞ −−→ −−→ T θ, r1s r1 = e−inθ On r1s r1 .

(20)

(21)

n=−∞

The boundary S is divided into the far-field boundary S 1 and the near-field boundary S 2 . By substituting Eq. (19) into Eq. (1), the boundary integral on S 1 in Eq. (1) can be denoted by A1 : !# ! I " −−→ −−−→ −−−→ −−→ i b 1 b 1 (22) eikk·r r T θ, r1s r1 + eikk·rm rm T θ, r1s rm1 B(θ, r1s )dθ, A1 = 8π where B(θ, r1s ) is the high-frequency moments defined by # Z " −−→ −−→ b 1 ikβe−ikk·rs rs − E(r1s r s ) p(r s , r0 )dS (r s ), B(θ, r1s ) = S1

6

(23)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and −−→

b 1 −− → ∂e−ikk·rs rs 1 . E(r s r s ) = ∂n(r s )

(24)

The B2B, B2H and H2H translation formulas in the implementation process of the high-frequency FMBEM are defined as: b

−−2−→1

B(θ, r2s ) = e−ikk·rs rs B(θ, r1s ),

(25)

−−→ H(θ, r1 ) = T (θ, r2s r1 )B(θ, r2s ),

(26)

and b

−−1−→2

H(θ, r2 ) = eikk·r r H(θ, r1 ).

(27)

where r2s is a point near r s , and r2 is a point near r. To the end, the boundary integrals on S 1 in Eq. (1) can be expressed as "I # I −−→ −−−→ i ikb k·r2 r 2 ikb k·rm2 rm 2 e H(θ, r )dθ + e H(θ, rm )dθ , (28) A1 = 8π where rm2 is a point near the mirror point r¯. By differentiating Eq. (28) with respect to the normal vector n(r), we can obtain the new boundary integral formulas on S 1 in Eq. (2), as follows: I  −−→ −−−→ I b 2 b 2  i  ∂eikk·r r ∂eikk·rm rm 2 2 B1 = H(θ, r )dθ + H(θ, rm )dθ . (29)  8π ∂n(r) ∂n(r) When the FMM is applied to the solution of the isogeometric BEM, a similar operation is implemented for the B2B, B2H and H2H translation formulas. Because the NURBS basis function is used, the high-frequency moments in Eq. (23) are rewritten as B(θ, r1s ) =

nf Z N0 X X e=1 j=0

1

−1

" # −−→ −−→ b 1 ¯ j J(ξ)dξ ikβe−ikk·rs rs − E(r1s r s ) R j,p f (ξ)p

(30)

where the subset S 1 contains N0 NURBS elements. The local coordinates ξ¯ and ξ in Eq. (30) are not independent; see Eq. (14). Finally, we obtain the solution of the boundary integrals with the NURBS basis function in Eqs. (1) and (2) by using Eqs. (28) and (29). 3. 2D acoustic shape sensitivity analysis 3.1. 2D acoustic shape sensitivity analysis based on conventional BEM By differentiating Eqs. (1) and (2) with respect to the design variable determining the structural shape, we can obtain the following sensitivity equations: c(r) p(r,˙ r0 ) +

Z " S

# ˙ ∂G(r s , r) − ikβ(r s )G(r s , r) p(r s , r0 )ds(r s ) = G(r,˙ r0 ), ∂n(r s ) 7

(31)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and ikβc(r) p(r,˙ r0 ) +

# ˙ ˙ ∂2G(r s , r) ∂G(r, r0 ) ∂G(r s , r) − ikβ(r s ) p(r s , r0 )ds(r s ) = , ∂n(r s )∂n(r) ∂n(r) ∂n(r)

Z " S

(32)

where ()˙ is a derivative operator that derives the function for the design variable. The sensitivity boundary integrals on S x in Eqs. (31) and (32) are singular, and can be rewritten as Z

Sx

Z ˙ ˙   G(ξ)J(ξ) − G0 (ξ)J0 (ξ) p(ξ) dξ G(r, r s )p(r s ) dS (r s ) = Sx

Z

Sx

+

˙ G0 (ξ)J0 (ξ)p(ξ) dξ

(33)

and Z

Sx

Z h ˙ ˙ i ∂2G(r, r s ) p(r s ) dS (r s ) = F 1 (ξ)J(ξ) − F01 (ξ)J0 (ξ) p(ξ) dξ ∂n(r)n(r s ) Sx Z ˙ F01 (ξ)J0 (ξ)p(ξ) dξ,

+

Sx

(34)

The first sensitivity expression on the right-hand side of Eqs. (33) and (34) is non-singular and can be solved using the Gauss integral method. However, the second one on the right-hand side is singular. By differentiating Eq. (10), we obtain the following expression: Z

Sx

Z

Sx

n n ˙ X X B˙i pi , G0 (ξ)J0 (ξ)p(ξ) dξ = Bi p˙i + i=1

i=1

n n ˙ X X ˙ 1 Di p˙i , F0 (ξ)J0 (ξ)p(ξ) dξ = Di pi + i=1

(35)

i=1

where B˙i and D˙ i can be derived as: Z Z 1 J˙0 J˙0 1 ˙ Bi = − [ln(J0 k) + 1] φi ln |ξ − b| dξ , φi dξ − 2π 2π −1 −1 Z 1 Z Z 1 k2 J˙0 k2 J˙0 1 J˙0 φi [ln(J0 k) + 1] D˙ i = − φi dξ − φi ln |ξ − b| dξ − dξ . 2 4π 4π (ξ − b)2 2πJ −1 −1 0 −1

(36)

By substituting Eqs. (12) and (14) into Eq. (36), we obtain a new expression for B˙i and D˙ i based on the NURBS basis function. The sensitivity expressions of some functions used to solve

8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

the sensitivity integral equation (31) are derived as ik G(r˙s , r) = − H1(1) (kr)r˙ 4   ! ˙ ˙ !  r˙ ∂r ik (1) ∂r ∂G(r s , r) ∂r  ik2 (1)  = − H1 (kr)  + H (kr)r˙ + ∂n(r s ) 4 r ∂n(r s ) ∂n(r s )  4 2 ∂n(r s )   ! ˙ ˙ !  r˙ ∂r ∂G(r s , r) ik (1) ∂r ∂r  ik2 (1) = − H1 (kr)  + H (kr)r˙ + ∂n(r) 4 r ∂n(r) ∂n(r)  4 2 ∂n(r)

! ˙   ∂2G(r s , r) ik  ˙ ) ˙ l (r)nl (r s ) + H (1) (kr) nl ˙(r)nl (r s ) + nl (r)nl (r = −kH2(1) (kr)rn s 1 ∂n(r s )∂n(r) 4r 2 h i ∂r ik ∂r 2H2(1) (kr) − krH3(1) (kr) r˙ + 4r ∂n(r) ∂n(r s )   ! ˙ ˙ !  2  ∂r ik (1) ∂r ∂r ∂r  + H2 (kr)  +  4 ∂n(r) ∂n(r s ) ∂n(r) ∂n(r s ) 

(37) (38) (39)

(40)

When the NURBS basis function is used to represent the structural boundary, , the coordinates of points r and r s are obtained by using NURBS basis function interpolation, as follows: x(ξ) =

n X

Ri,p (ξ)Pi

(41)

i=0

where x(ξ) denotes the coordinates of boundary points, and Pi represents the coordinates of control points. By differentiating Eq. (41) with respect to the shape design variable, we obtain the following expression for the coordinate sensitivity: ˙ = x(ξ)

n X

Ri,p (ξ)P˙i

(42)

i=0

Using a similar procedure, we obtain the sensitivity expression of the coordinate derivative used for the boundary integral solution: "

˙ # X n dRi,p (ξ) ˙ dx(ξ) = Pi . dξ dξ i=0

(43)

After successfully solving the singular integral, combining Eqs. (31) and (32) to generate the Burton–Miller formulation, and discretizing the boundary S , we obtain the following sensitivity formulation:   ˙ − ikGB˙ p + [H − ikGB] p˙ = p˙ H (44) i Using the above equation, we obtain the sound pressure sensitivity values at the special control points used for the physical field approximation. Then, using Eq. (31), we obtain the sound pressure sensitivity p f in the acoustic domain. 9

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3.2. Diagonal formulations for half-space acoustic shape sensitivity analysis The sensitivity formulas of the boundary integrals in Eq. (31) with respect to the design variable can be rewritten as # Z " ˙ ∂G(r s , r) − ikβ(r s )G(r s , r) p(r s , r0 )ds(r s ) = D1 + D2 + D3 (45) ∂n(r s ) S where D1 =

Z "

# ˙ ∂G(r s , r) − ikβ(r s )G(r s , r) p(r s , r0 )ds(r s ) ∂n(r s )

(46)

D2 =

Z "

# ∂G(r s , r) − ikβ(r s )G(r s , r) p(r s˙, r0 )ds(r s ) ∂n(r s )

(47)

D3 =

Z "

# ∂G(r s , r) ˙ − ikβ(r s )G(r s , r) p(r s , r0 )ds(r s ) ∂n(r s )

(48)

and

S

S

S

By differentiating Eq. (19) with respect to the design variable, one can obtain the following expression: ! ! I  ˙ −−→ ˙−−1−→  b −−1→ −− → −−1−→ i  ikbk·r1 r ˙ b 1 1 1 ik k· r r G(r, r s ) = T θ, r s r + e m m T θ, r s rm  e−ikk·rs rs dθ+ e 8π ! !# I " −−→ ˙ −−1→ −−−→ −−→ −−−→ i b 1 b 1 eikk·r r T θ, r1s r1 + eikk·rm rm T θ, r1s rm1 e−ikbk·rs rs dθ (49) 8π Then, by substituting Eqs. (19) and (49) into the above three equations, we obtain the following formulations: ! ! I  ˙ −−→ ˙−−1−→  −− → −−1−→ i  ikbk·r1 r b 1 1 1  ik k· r r m m T θ, r s r + e T θ, r s rm  B(θ, r1s )dθ D1 = e  8π ! !# I " −−→ −−−→ −−→ −−−→ i b 1 b 1 + eikk·r r T θ, r1s r1 + eikk·rm rm T θ, r1s rm1 B1 (θ, r1s )dθ, (50) 8π i D2 = 8π

I "

i D3 = 8π

I "

and

−−→ ikb k·r1 r

! !# −−−→ −− → −−1−→ 1 1 ikb k·rm1 rm 1 T θ, r s r + e T θ, r s rm B2 (θ, r1s )dθ,

(51)

−−→ ikb k·r1 r

!# ! −−−→ −−1−→ −− → 1 1 1 ikb k·rm1 rm T θ, r s rm B3 (θ, r1s )dθ, T θ, r s r + e

(52)

e

e

where B1 (θ, r1s ) =

Z

S1

−−→ −− →˙ 1 ikβe−ikbk·rs rs − E(r1s r s )p(r s , r0 )dS (r s ), 10

(53)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

B2 (θ, r1s ) =

Z "

# −−→ −−→ b 1 ikβe−ikk·rs rs − E(r1s r s ) p(r s˙, r0 )dS (r s ),

(54)

B3 (θ, r1s ) =

Z "

# −−→ −−→ b 1 ˙ ikβe−ikk·rs rs − E(r1s r s ) p(r s , r0 )dS (r s ).

(55)

S1

and S1

The B2B, B2H, and H2H translation formulas for the above three equations are the same as Eqs. (25)– (27). Finally, D1 , D2 , and D3 can be expressed in terms of the local expansion coefficients as I  I ˙ −−→ ˙−−2−→  i  2 b b 2 2 ik k· r r ik k· r r m m D1 =  e H(θ, r )dθ + e H(θ, rm )dθ + 8π "I # I −−→ −−−→ i b 2 b 2 eikk·r r H 1 (θ, r2 )dθ + eikk·rm rm H 1 (θ, rm2 )dθ , (56) 8π i 8π

"I

i D3 = 8π

"I

D2 =

b

−−→ 2

eikk·r r H 2 (θ, r2 )dθ +

I

# −−−→ b 2 eikk·rm rm H 2 (θ, rm2 )dθ ,

(57)

and −−→ ikb k·r2 r

e

3

2

H (θ, r )dθ +

I

e

−−−→ ikb k·rm2 rm

H

3

(θ, rm2 )dθ

#

.

(58)

H 1 ,H 2 , and H 3 generated by B1 , B2 and B3 can be obtained by using Eqs. (25)- (27). 4. Optimization analysis of noise barriers with absorbing material In the above section, acoustic sensitivity analysis with respect to the design variables was presented in detail. For the shape optimization problem, the design variables are the coordinates of the control points, which determine the geometry. Further, the design variables are the acoustic admittance of the elements if we optimize the material distribution. Two types of optimization analysis are presented in this paper: optimization analysis of the shape and of the topological distribution of the absorbing material. In both types of optimization analysis, the OC method is used to update the design variables in order to obtain the final optimal solution. Note that a new material interpolation scheme for acoustic problems based on SIMP needs to be created because the interpolation variable is not the real structural density, but a fictitious material density that determines the normalized surface admittance. 4.1. Acoustic shape optimization analysis of noise barriers The optimal shape problem of minimizing the sound pressure on a specified reference plane is considered in this paper, and the optimization model can be formulated as:   min Πρ = p f p f      s.t. A(ϑ) − A0 6 0 (59)      l u  ϑi 6 ϑi 6 ϑi , i = 1, 2, . . . , Nϑ , 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where p f denotes the sound pressure vector at several field points located in a protected reference zone, p f denotes the complex conjugate vector of p f , and Πρ denotes the objective function. A0 is the area of the initial closed structure, and A(ϑ) is the area of the closed structure in the optimization procedure, where ϑ denotes the shape design variable. Nϑ is the number of shape design variables, and ϑli and ϑui are the lower and upper values of the shape design variables ϑi , respectively. The control points also determine the shape and size of the considered structure. Furthermore, directly setting the control points as the shape design variables can result in a simpler procedure for optimization analysis, and the associated optimization model Eq. (59) can be rewritten as   min Πρ = p f p f      s.t. A(x) − A0 6 0 (60)      i,l i i,u  x 6 x 6 x , i = 1, 2, . . . , N x , where x is the vector of the coordinates of the control points. N x denotes the number of control points that change with the shape variable ϑ, and xi,l and xi,u denote the lower and upper values, respectively, of the coordinate of the i-th control point xi . The sensitivity of the objective function Πρ to the shape design variable ϑ can be expressed as   ! ∂Πρ ∂ p f p f ∂p f = = 2< p f (61) ∂ϑ ∂ϑ ∂ϑ where Nx ∂p f ∂p f (x1 , x2 , . . . , xNx ) X ∂p f ∂xi = = ∂ϑ ∂ϑ ∂xi ∂ϑ i=1

(62)

By using Eq. (1) with c(r) = 1, we obtain the expression for the sound pressure vector p f in the acoustic domain: p f = − [H1 − ikG1 B] p + pi

(63)

where the matrices H1 and G1 are similar to the matrices in Eq. (5), but the source nodes are different from the special control points. By differentiating p f with respect to the design variables xi , the following formulation can be obtained: " # ∂p f ∂H1 ∂G1 ∂p ∂pi =− − ik i B p − [H1 − ikG1 B] i + i (64) ∂xi ∂xi ∂x ∂x ∂x ∂p f ∂p where i refers to p˙ in Eq. (44). By using Eqs. (44) and (64), we can obtain . The matrix∂x ∂xi vector products in Eqs. (44) and (64) can be easily handled using the FMM to save computational time. Finally, we obtain the sensitivity of the objective function using Eq. (61). After obtaining the value of the objective function and its sensitivity to the shape design variable, we use the OC method to update the design variables and search for the final optimal solution. When the OC method based on gradient analysis is used for optimization analysis, it is important to solve the constraint function A (the area of the closed structure) and its sensitivity 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

value efficiently as N

e 1X A= 2 e=1

and N

e ∂A 1 X = ∂ϑ 2 e=1

Z

ξe+1

ξe

"

Z

ξe+1

ξe

x(ξ) · nJ(ξ) dξ

# ∂n ∂J(ξ) ∂x(ξ) dξ · nJ(ξ) + x(ξ) · J(ξ) + x(ξ) · n ∂ϑ ∂ϑ ∂ϑ

(65)

(66)

4.2. Optimization analysis of absorbing material distribution The discretized optimization model of the distribution of the acoustic absorbing material on the structure is formulated as follows:   min Πρ = p f p f       Ne  X   (67) s.t. ρe ve − V0 6 0      e=1      06ρ 6ρ 61 min

e

Sometimes, it is more convenient to use the average sound pressure level (ASPL) to replace the conventional objective function presented in Eq. (67). The ASPL is defined as m pi 20 X log10 (68) AS PL = m i=1 p0

where m is the number of computing points located in the protected zone. pi denotes the sound pressure magnitude at the i-th computing point, and p0 is the reference sound pressure. ρi in Eq. (67) is not the real density of the structure, but the acoustic density governing the acoustic admittance βi . When the SIMP interpolation scheme is used for the topology optimization analysis, the relationship between ρi and βi can be expressed as βi = β0 ρni + β1 (1 − ρni ),

0 6 ρi 6 1

(69)

where n is the penalization factor. The use of larger n causes the intermediate density to rapidly approach 0 (void) or 1 (solid). The value of n may influence the efficiency and final solution of the optimization procedure and will be investigated in numerical examples. β0 and β1 in Eq. (69) denote the acoustic admittance of two different types of absorbing material. The aim of topology optimization analysis is to improve the performance of noise barriers by finding an optimal combinatorial distribution of the two types of absorbing material. To simplify the analysis, we set β0 to 1 and β1 to 0. Thus, the above equation can be rewritten as βi = β0 ρni

(70)

By differentiating the objective function with respect to the element’s acoustic density ρi , we obtain the sensitivities of the objective function in (67) as   ! ∂Πρ ∂ p f p f ∂p f = = 2< p f (71) ∂ρi ∂ρi ∂ρi 13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

By differentiating p f with respect to the element acoustic density ρi , the following formulation is obtained: " #  #  −1 " ∂p f ∂B ∂B = ikG1 p − H1 − ikG1 B H − ikGB p (72) ikG ∂ρi ∂ρi ∂ρi where operations such as A−1 b represent solving Ax = b, which can also be solved efficiently by the FMM and the iterative solver GMRES. The OC method, based on gradient approaches is applied to update the design variables and obtain the final optimization solution. First, the Lagrangian function of Eq. (67) based on the OC method can be expressed as  N N e e  X X λe2 (ρmin − ρe ) L =Πρ + λ1  ρe ve − V0  + e=1 e=1 (73) Ne X e + λ3 (ρe − 1) e=1

where λ1 , λe2 , and λe3 are the Lagrangian multipliers. By differentiating the above equation with respect to the design variables, we can obtain the following formulation: ∂Πρ ∂L = + λ1 ve − λe2 + λe3 = 0 ∂ρe ∂ρe

(74)

To simplify the analysis, the corresponding multipliers, λe2 and λe3 , are equal to zero in this paper. Thus, Eq. (74) can be rewritten as ∂Πρ ∂L = + λ1 ve = 0 ∂ρe ∂ρe

(75)

! ∂Πρ 1 1 Be = − =1 λ1 ∂ρe ve

(76)

We adjust the above equation as

An optimal solution is obtained when Be = 1. In this paper, the above equation is used directly to update the design variables, and the corresponding iteration updating scheme is defined as   η (k)   Le(k) if B(k) ρ(k)  e e 6 Le ,     η    (k) η (k) (k) (77) = ρ(k+1) ρ(k) if Le(k) < B(k) Be ρe e e e < Ue ,       η    Ue(k) if Ue(k) 6 B(k) ρ(k) e e ,

where η denotes a given tuning parameter, and ρ(k) e denotes the design variable in the k-th iteration. Le(k) and Ue(k) are the lower and upper bounds of the design variable ρ(k) e , respectively, and can be defined as   (k) Le(k) = max ρmin , ρ(k) (78) e − ml 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

  (k) Ue(k) = min 1, ρ(k) e + ml

where ml is the move limit of the design variable, which can be defined as   ml(k) = min α(k) ml(0) , mmin

(79)

(80)

where ml(0) is the initial move limit of the design variable, which is typically set to 0.15. The scale coefficient α is set to 0.96, and the minimal move limit mmin is 0.001. The iterative convergence criterion is whether the relative error in the objective function values is less than a prescribed value τ and is written as k+1 Π − Πk (81) err = < τ Πk

where Πk denotes the objective function at the k-th iteration step. The detailed optimization algorithm is as follows: 1. 2. 3. 4. 5. 6. 7.

Initialize the optimisation parameters: ρ = ρ0 , iter = 1, err = 1; Response analysis based on IGA-FMBEM by solving Eq. (5); Sensitivity analysis using direct method by solving Eqs. (71) and (72); Update the design variables ρi by using OC and SIMP method; Modify ρ using smoothed Heaviside function; If err > τ then goto step 2, else goto step 7; End.

5. Numerical examples 5.1. Scattering from an infinite rigid cylinder In this example, we consider acoustic scattering of a plane incident wave with unit amplitude on a rigid cylindrical shell with radius a = 1 m centered at point (0, 0), as shown in Fig. 1. By using this example with an analytical solution, we can validate the algorithm proposed in this paper. This model can be simplified to a 2D acoustic problem with a circular boundary. The analytical solution of the acoustic pressure at any point (r, θ) is given as (48) p(r, θ) = −

∞ X n=0

εn in

Jn0 (ka) Hn(1)0 (ka)

Hn(1) (kr) cos(nθ),

(82)

where εn denotes the Neumann symbols, i.e., ε0 = 1 ; εn = 2 when n is greater than 0. ( )0 stands for the differentiation with respect to ka. When a is chosen as the design variable, one can obtain the analytical solution of the sound pressure sensitivity by differentiating Eq. (82) with respect to the design variable, as follows: " 0 # ∞ X Jn (ka) 0 (1) ∂p(r, θ) =− εn in (1) Hn (kr) cos(nθ). 0 ∂a Hn (ka) n=0

(83)

Figure 2 shows the numerical and analytical solutions for different frequencies. According to this figure, the sound pressure obtained by the IGA-FMBEM does not agree with the analytical 15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

P 2

P 3

P 1

1 .0

0 .5

0 .0

N U R B S c u rv e C o n tro l p o in t

P 4

P 8

P 0

-0 .5

-1 .0

P 5

P 7

P 6 -1 .0

-0 .5

0 .0

0 .5

1 .0

Fig. 1: Example of a circle represented by a closed NURBS curve of degree 2 generated by 9 control points, p0 = (1, 0), p1 = (1, 1), p2 = (0, 1), p3 = (−1, 1), p4 = (−1, 0), p5 = (−1, −1), p6 = (0, −1), p7 = (1, −1), p8 = (1, 0).

solutions at fictitious eigenfrequencies, where the conventional single BIE is used. However, the IGA-FMBEM-BM (IGA-FMBEM based on the Burton–Miller method) yields excellent solutions at all frequencies. Sample internal points are evenly distributed on a circle with a radius of 2 m. Figure 3 shows that the numerical results are in perfect agreement with the analytical solutions, demonstrating the validity of the algorithm. Figure 4 shows that the solution converges when the boundary mesh is refined, indicating the accuracy of the presented algorithm compared with that of the constant BEM; the relative error is the surface error, which is solved according to (46). The CPU time used to calculate the sensitivity values at the test point is plotted in Fig. 5, which demonstrates the efficiency of the IGA-FMBEM for 2D acoustic design sensitivity analysis. The control points of the IGA determine the shape of the represented structure. Changing the positions of the control points can change the shape of the represented structure. Thus, it is necessary to investigate the effect of the coordinates of the control points on the scattering pressure from the represented structure. Figure 6 shows the sensitivity contours on a square plane with respect to different design variables, such as P0 , P2 , P4 , and P6 . We find that every control point has a large effect on the sensitivity value in the domain around it because each control point determines the boundary shape around it. 5.2. Shape optimization analysis of vertical noise Noise from expressways and railways is an important source of noise pollution. It can be rectified by erecting barriers between the noise source and the protected zone. For a long straight barrier and a line source, a 2D model can be used to predict the acoustic field. The barrier is assumed to have infinite length and a uniform cross section along the length. Figure 7 shows the parameters of the structure’s size and location. NURBS curves are used to represent the boundary. The red line on the left side of the vertical noise barrier will be optimized by inserting 16

1 .6 1 .4

S o u n d p re s s u re (P a )

1 .2 1 .0 0 .8 0 .6 0 .4

A n a ly tic a l IG A -F M B E M -C B IE IG A -F M B E M -B M

0 .2 0 .0 0

5 0 0

1 0 0 0 F re q u e n c y (H z )

1 5 0 0

2 0 0 0

Fig. 2: Sound pressure at point (2,0) in terms of frequency.

A m p litu d e o f s o u n d p re s s u re s e n s itiv ity

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

A A IG IG

2 .0

1 .6

n a ly n a ly A -F A -F

tic tic M M

a l f a l f B E B E

o r k = o r k = M fo M fo

π 2π r k = π r k = 2π

1 .2

0 .8

0 .4

0 .0

0 .0

0 .5

1 .0

1 .5

2 .0

θ/π

Fig. 3: Sound pressure at sample points with two different wave number k.

17

R e la tiv e e rro r

1 0

-2

1 0

-3

1 0

-4

1 0

-5

1 0

-6

1 0

-7

1 0

-8

IG A -F M B E M C B E M

1 0

2

1 0

3

D e g re e o f fre e d o m

Fig. 4: Relative error of the pressure sensitivity at point (2, 0) with k = 2π.

4

1 0

C B E M -G M R E S IG A -F M B E M

C P U tim e [s e c ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1 0

3

1 0

2

1 0

1

1 0

0

1 0

3

1 0

4

1 0

5

1 0

6

D e g re e o f fre e d o m

Fig. 5: CPU time used to calculate the pressure sensitivity values at point (2, 0) with k = 2π.

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a) P0 is the design variable

(b) P2 is the design variable

(c) P4 is the design variable

(d) P6 is the design variable

Fig. 6: Sensitivity contour on a square plane 10 × 10 m2 with respect to different shape design variables

19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

more control points and changing their locations. These inserted control points are distributed uniformly on the left side of the vertical noise barrier in the initial step; their coordinate on the Y axis is fixed, but that on the X axis changes. That is, the coordinate on the X axis represents the shape design variable. The initial value for the design variable is 5.0, and the lower and upper bounds are set to 4.9 and 5.1, respectively. The number of inserted control points is N, and using a larger N will generate a more complex curve.

P3 =P4

P1 =P2 (6,3)

(16,3)

reference plane 3m

10 × 10 points

the line to be optimized

source (0,1)

P5

P0

(6,0.1)

(16,0.1)

5m

Fig. 7: Example of a vertical noise barrier represented by NURBS of degree 2 generated by 6 control points, that is, p0 = (5.2, 0), p1 = (5.2, 3), p2 = (5.2, 3), p3 = (5, 3), p4 = (5, 3), p5 = (5, 0); the knot vector Ξg = {0.0, 0.0, 0.0, 0.45, 0.45, 0.55, 0.55, 1.0, 1.0, 1.0}.

Figure 8 shows the objective function and constraint function in terms of the iteration step for different numbers of inserted control points. Obviously, all the curves for the objective function converge very rapidly in terms of iteration step. A larger N produces a smaller optimal value of the objective function. In fact, a NURBS curve with more control points can more accurately represent the boundary of a complex structure. However, the computational efficiency of numerical solution based on the IGA-FMBEM will decrease as the number of control points increases. Note that choosing a suitable number of control points is important for representing the structure’s boundary by a NURBS curve, especially for shape optimization analysis, because it typically takes much more CPU time. Figure 9 shows the final optimal solution for different numbers of control points. We find that a NURBS curve with more control points can represent a more complex structural boundary. Further, this figure demonstrates the validity and efficiency of the algorithm proposed in this paper for shape optimization analysis of complex structural boundaries. All the results in the above subfigures are solved at a single frequency. It is worth considering the problem of whether the optimization solution changes with the computing frequency. Figure 10 reveals the answer. From this figure, we find that the optimization solutions for different computing frequencies are different and the optimized line becomes more complex and curved as the frequency increases. Figure 11 compares the sound performance of the initial and optimized shapes of the noise barrier. The plotted curve represents the values of the SPL at 10 × 10 points on the reference plane, where the values on the X axis denote the number of computing points on the reference plane from bottom to top. We find that the improvement in the performance of the optimized sound barrier at low frequencies, for examples, at 100 Hz, is not very clear, but it becomes more pronounced at higher computing frequencies. This phenomenon indicates that the optimal results are strongly frequency-dependent. In fact, the performance of the sound barrier is 20

-4

x 1 0 4 0

0 .6 0

N = 5 N = 1 0 N = 1 5 N = 2 0

3 5

0 .5 8 0 .5 6 2

)

3 0 2 5

A re a (m

O b je c tiv e fu n c tio n

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2 0 1 5

0 .5 4 0 .5 2 N = 5 N = 1 0 N = 1 5 N = 2 0

0 .5 0

1 0

0 .4 8 5 0

1 0

2 0

3 0 4 0 Ite ra tio n s te p

5 0

6 0

0 .4 6

7 0

(a) Objective function with iteration step

0

1 0

2 0

3 0 4 0 Ite ra tio n s te p

5 0

6 0

7 0

(b) Constraint function in terms of iteration step

Fig. 8: Objective function and Area function in terms of iteration step for different numbers of inserting control points; The computing frequency is 400 Hz.

3 .0

3 .0

N U R B S c u rv e C o n tro l p o in ts

2 .5 2 .0

2 .0

1 .5

1 .5

1 .0

1 .0

0 .5

0 .5

0 .0

4 .6

4 .8

5 .0

5 .2

5 .4

5 .6

5 .8

N U R B S c u rv e C o n tro l p o in ts

2 .5

N = 5

0 .0

6 .0

N = 1 0

4 .6

4 .8

5 .0

(a) N=5

3 .0

3 .0

N = 1 5

2 .0

1 .5

1 .5

1 .0

1 .0

0 .5

0 .5 4 .6

4 .8

5 .0

5 .2

5 .4

5 .6

5 .8

6 .0

5 .6

5 .8

N U R B S c u rv e C o n tro l p o in ts

2 .5

2 .0

0 .0

5 .4

(b) N=10

N U R B S c u rv e C o n tro l p o in ts

2 .5

5 .2

0 .0

6 .0

(c) N=15

N = 2 0

4 .6

4 .8

5 .0

5 .2

5 .4

5 .6

(d) N=20

Fig. 9: Final optimal solution for different number of control points with frequency 400 Hz

21

5 .8

6 .0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

more sensitive to incident sound waves with shorter wavelengths (at higher frequencies) because less diffraction occurs. On the other hand, this phenomenon suggests the difficulty of obtaining a suitable shape optimization solution that is useful at all frequencies. To overcome this difficulty, we formulate an average objective function in a frequency range, as follows Z ω2 1 Πρ dω (84) Π1 = ω2 − ω1 ω1 By using this operation, we replace the optimal solution at a single frequency with the average optimization solution for multipole frequencies. For example, a frequency range of 500–1000 Hz is given. The corresponding optimization solution and mean SPL for frequencies in this range are shown in Fig. 12.

3 .0

3 .0

N U R B S c u rv e C o n tro l p o in ts

2 .5

f= 1 0 0 H z

2 .0

1 .5

1 .0

1 .0

0 .5

0 .5 4 .6

4 .8

5 .0

5 .2

5 .4

5 .6

f= 4 0 0 H z

2 .0

1 .5

0 .0

N U R B S c u rv e C o n tro l p o in ts

2 .5

5 .8

6 .0

0 .0

4 .6

4 .8

(a) f=100 Hz

3 .0

1 .0

1 .0

0 .5

0 .5 4 .8

5 .0

5 .2

5 .4

(c) f=700 Hz

5 .8

6 .0

5 .6

f= 1 0 0 0 H z

2 .0 1 .5

4 .6

5 .6

N U R B S c u rv e C o n tro l p o in ts

2 .5

1 .5

0 .0

5 .4

3 .0

f= 7 0 0 H z

2 .0

5 .2

(b) f=400 Hz

N U R B S c u rv e C o n tro l p o in ts

2 .5

5 .0

5 .8

6 .0

0 .0

4 .6

4 .8

5 .0

5 .2

5 .4

5 .6

5 .8

6 .0

(d) f=1000 Hz

Fig. 10: Final optimal solution based on IGA-FMBEM with N = 15 and different frequencies.

5.3. Optimization analysis of absorbing material distribution In this example, the problem of acoustic scattering from a vertical sound barrier with an approximately circular upper structure is investigated to demonstrate the validity and applicability 22

7 0

6 0 In itia l v a lu e O p tim a l v a lu e

6 5

5 0 S P L (d B )

S P L (d B )

In itia l v a lu e O p tim a l v a lu e

5 5

6 0 5 5

4 5

5 0

4 0

4 5

3 5

4 0 0

2 0

4 0

6 0

8 0

3 0

1 0 0

0

2 0

(a) Frequency = 100 Hz

4 0

6 0

8 0

1 0 0

(b) Frequency = 400 Hz

5 5 5 0

In itia l v a lu e O p tim a l v a lu e

5 0

In itia l v a lu e O p tim a l v a lu e

4 5 S P L (d B )

4 5 S P L (d B )

4 0

4 0 3 5

3 5 3 0 3 0 2 5 2 5 0

2 0

4 0

6 0

8 0

1 0 0

0

2 0

(c) Frequency = 700 Hz

4 0

6 0

8 0

1 0 0

(d) Frequency = 1000 Hz

Fig. 11: Comparison of performance of optimized sound barrier on decreasing noise with initial shape with several computing frequencies.

3 .0

4 8

N U R B S c u rv e C o n tro l p o in ts

2 .5

4 4

f= 5 0 0 -1 0 0 0 H z

2 .0

M e a n S P L (d B )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1 .5 1 .0

4 0 3 6 3 2 In itia l v a lu e O p tim a l v a lu e

0 .5 2 8 0 .0

4 .6

4 .8

5 .0

5 .2

5 .4

5 .6

5 .8

6 .0

2 0 0

2 5 0

3 0 0 3 5 0 4 0 0 F re q u e n c y (H z )

Fig. 12: Final optimal shape and Mean SPL in a frequency range.

23

4 5 0

5 0 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2.8m 3.0m

Protected zone

d

Source 1.0m

11m

20m

10m

Fig. 13: Cross sections of the vertical barrier with a circle top structure with a radius r = 1m.

of the algorithm combining the FMBEM and OC method. The vertical barrier is erected on a rigid, flat, and sufficiently large surface. Figure 13 shows its cross-sectional contour over a planar surface 10 m from a coherent homogeneous monofrequency line source situated 1 m above the ground. The height of the vertical structure on the bottom of the sound barrier is set to 3.0 m, and the width is 0.2 m. The radius of the upper structure is 1 m. We choose 200 computing points that are distributed uniformly in the protected zone and use the sound pressure or SPL at these points as the objective function. The surface of the sound barrier is covered with an absorbing material. The relative parameters for the sound barrier structure and the initial values used for the optimization analysis are listed in Table 1. In the following examples, ρmin is set to 10−3 , and the iterative convergence criterion τ is set to 10−4 . Table 1: Relative parameters for sound barriers and optimization analysis.

Parameter

Symbol

Value

Unit

Density

ρ

1.2

kg m−3

Speed of sound

c

343

m s−1

Thickness of sound barrier

t

0.2

m

Acoustic density for admittance

ρi

0.5

-

Penalizaation factor

n Ne P

5

-

0.5

-

Constraint condition

V0 /

e=1

ve

Figure 14 shows the values of the objective function (the ASPL) in terms of the iteration number for several different frequencies. The ASPL decreases rapidly in the initial iteration steps for most of the data curves and then converges rapidly. The optimized ASPL value after convergence is smaller than the initial value, indicating the validity and efficiency of the optimization algorithm proposed in this paper. By observing the change in the absorbing material distribution in terms of the iteration number (Fig. 15), we can clearly see the performance of the proposed optimization algorithm. The 24

1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 1 0 0

5 2

O b je c tiv e fu n c tio n A S P L (d B )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4 8 4 4 4 0 3 6 3 2

H z H z H z H z H z H z H z H z H z 0 H z

2 8 2 4 2 0 1 6 0

5

1 0

1 5

2 0

Ite ra tio n n u m b e r Fig. 14: Objective function values in terms of iteration number with several different frequencies

parameters used to obtain the optimization solution are given in Table 1, and different colors in the figure represent the values of the acoustic density ρi . Black represents boundary elements having a material density of one. In contrast, yellow represents boundary elements having a material density of zero, and it indicates that no material should be located in these elements. Colors between black and yellow indicate the presence of gray elements, which have acoustic densities between 0 and 1. We want to avoid gray elements because they introduce some difficulty for industrial manufacture, although a Heaviside-like function method can be applied to overcome this difficulty; see (41) for detailed information. From Figs. 15a and 15b, we find that gray elements appear in the initial optimization procedure. However, these gray elements disappear rapidly as the iteration number increases, as shown in Figs. 15e and 15f. This phenomenon shows that the density operator method based on a smoothed Heaviside-like function used in this paper can be successfully applied to exclude gray elements. Figures 15d, 15e, and 15f are very similar, indicating that the optimization algorithm proposed in this paper has good convergence and stability for solving this type of problem. It is essential to study the effect of the initial values of design variables on the optimization solution. If the optimization solution changes greatly when a different initial value of the design variable is chosen, the proposed optimization algorithm does not exhibit good stability. We want to avoid this phenomenon. Figure 16 illustrates the stability of the proposed algorithm, where the parameter for the constraint condition is set to one. From this figure, we find that the final optimization solutions are very similar when the initial value is between 0.3 and 0.9, which demonstrates that the proposed algorithm has good adaptability and stability for different initial values. Of course, too small initial values, for example, 0.1 and 0.2, will produce a poor optimization solution, so we should avoid the use of too small initial values because they may cause the optimization analysis to converge to a local minimum. 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Absorbing material layout for iteration number iter=3 at frequency 500Hz 5

1

Absorbing material layout for iteration number iter=5 at frequency 500Hz 5

1

4.5

0.9

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

0

12

0

8

(a) Iteration number =3

9

10

11

12

(b) Iteration number =5

Absorbing material layout for iteration number iter=7 at frequency 500Hz 5

1

4.5

0.9

Absorbing material layout for iteration number iter=9 at frequency 500Hz 5

1

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

0

12

0

8

(c) Iteration number =7

9

10

11

12

(d) Iteration number =9

Absorbing material layout for iteration number iter=11 at frequency 500Hz 5

1

Absorbing material layout for iteration number iter=15 at frequency 500Hz 5

1

4.5

0.9

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

(e) Iteration number =11

12

0

0

8

9

10

11

12

(f) Iteration number =15

Fig. 15: optimisation procedure of material distribution for vertical sound barrier with absorbing material: Constraint condition is 0.5; Computing frequency is 500 Hz. 26

3 6 In itia In itia In itia In itia In itia In itia In itia In itia In itia

3 4

O b je c tiv e fu n c tio n A S P L (d B )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3 2 3 0 2 8 2 6 2 4

l v a l v a l v a l v a l v a l v a l v a l v a l v a

lu e lu e lu e lu e lu e lu e lu e lu e lu e

= 0 .1 = 0 .2 = 0 .3 = 0 .4 = 0 .5 = 0 .6 = 0 .7 = 0 .8 = 0 .9

2 2 2 0 1 8 1 6 0

2

4

6

8

1 0

1 2

1 4

1 6

1 8

Ite ra tio n n u m b e r Fig. 16: ASPL values in terms of iteration number for a series of different initial values of design variables

As mentioned in regard to Eq. 69, the penalization factor n, which determines the normalized surface admittance, may influence the optimization solution and computational efficiency. From Fig. 17, we find that different objective functions, ASPL values, and convergence iteration numbers are obtained for different n values. When n is 1 or 3, the ASPL values after optimization are larger than those for larger n. In fact, a too small penalization factor will not produce a good optimization solution because many gray elements will exist. Thus, the use of a very small n should be avoided. The ASPL values after convergence for n = 5 to n = 11 are very similar. However, as n increases further, we obtain the wrong solution with iteration number 3 for n = 15. In sum, we obtained the following results: a larger n will generate more efficient optimization solution, but a too large n will break the iteration solution; a smaller n may cause the optimization procedure to converge to a local minimum, and a solution may not be obtained efficiently when n is very small, for example, 1 in this case. Finally, a reasonable empirical value of n = 5 is given. It is still necessary to know whether the optimization solution changes with the computing frequency. If changing the computing frequency does not produce different optimization solutions, the algorithm proposed in this paper would have good adaptability and stability with respect to the computing frequency. Figure 18 shows that the optimization solutions for different computing frequencies are different. In fact, it indicates the difficulty of creating a suitable optimal material distribution that can be used for several frequencies. To overcome this difficulty, we make an average objective function in a frequency range, as shown in Eq. 84. The procedure for optimizing the absorbing material distribution in a frequency range is similar to that for shape optimization analysis in subsection 5.2. Thus, in the interest of brevity, the procedure for optimization in a frequency range is not given here. Figure 19 shows the optimal absorbing material distribution at 500 Hz for two constraint values. A larger material constraint value denotes that more absorbing material can be used, 27

3 6 n = 1 n = 3 n = 5 n = 7 n = 1 n = 1

3 4

O b je c tiv e fu n c tio n A S P L (d B )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3 2 3 0 2 8

, ite r= , ite r= , ite r= , ite r= 1 , ite r 5 , ite r

1 4 1 5 1 7 1 4 = 1 6 = 4

2 6 2 4 2 2 2 0 1 8 0

2

4

6

8

1 0

1 2

1 4

1 6

1 8

2 0

Ite ra tio n n u m b e r Fig. 17: ASPL values in terms of iteration number for different penalization factor

which yields a lower ASPL. The study of material topology optimization analysis with a suitable constraint function value is advantageous for guiding the industrial manufacture of structures. More solutions can be found in Table 2. We find that the optimized ASPL is smaller than the initial value at all frequencies, which demonstrates the validity and efficiency of the algorithm proposed in this paper. On the other hand, increasing the constraint value, that is, using more absorbing material, can efficiently improve the performance of the noise barrier. To further demonstrate the capability and adaptability of the algorithm proposed in this paper, examples of other types of noise barrier are given in Fig. 20. Herein, the parameter for the constraint condition is set to 0.7, and the computing frequency is 100 Hz. 6. Conclusion A novel algorithm based on the IGA-FMBEM is presented for the simulation of 2D halfspace acoustic sensitivity problems with absorbing boundary conditions. The objective function for sensitivity analysis of a design variable is derived. The Burton–Miller method is used to obtain correct solutions at all frequencies, and the hyper-singular boundary equation is solved directly using the Hadamard finite part integration method. When IGA is applied to the FMBEM, the accuracy of the numerical method can be efficiently improved by using a simple example with an analytical solution. After the sensitivity analysis, the OC method is used to update the design variable during the optimization procedure. Two types of optimization are studied, optimization of the shape and of the topological distribution of the material. For the material distribution optimization, a new concept called the acoustic density is presented, and the associated interpolation scheme based on SIMP is produced. Several examples are presented to validate the proposed algorithm for practical optimization problems. 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Absorbing material layout for constraint condition 0.5 at frequency 100Hz 5

1

Absorbing material layout for constraint condition 0.5 at frequency 300Hz 5

1

4.5

0.9

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

0

12

0

8

9

(a) 100 Hz

10

11

12

(b) 300 Hz

Absorbing material layout for constraint condition 0.5 at frequency 500Hz 5

1

Absorbing material layout for constraint condition 0.5 at frequency 700Hz 5

1

4.5

0.9

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

12

0

0

8

(c) 500 Hz

9

10

11

(d) 700 Hz

Fig. 18: The final material distribution of vertical sound barrier in terms of frequencies

29

12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Absorbing material layout for constraint condition 0.3 at frequency range 500Hz Absorbing material layout for constraint condition 0.7 at frequency range 500Hz 5

1

5

1

4.5

0.9

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

12

0

0

8

(a) Constraint condition 0.3

9

10

11

12

(b) Constraint condition 0.7

Fig. 19: The final material distribution for two different constraint condition with computing frequency 500 Hz

Table 2: Comparison of objective function value between initial state and optimized state for absorbing material distribution with two different constraint condition in terms of frequencies

Frequency [Hz] 100 200 300 400 500 600 700 800 900 1000

Constraint condition 0.4 Initial Optimized Decreasing 48.7 41.6 7.1 34.4 25.9 8.5 42.6 32.3 10.3 41.4 29.0 12.4 32.3 20.3 12.0 37.9 33.0 4.9 35.8 20.1 15.7 32.4 21.7 10.7 36.3 28.7 7.6 23.6 18.9 4.7

30

Constraint condition 0.7 Initial Optimized Decreasing 48.6 40.9 7.7 31.8 24.6 7.2 38.7 30.1 8.6 36.1 28.6 7.5 26.7 17.7 9 33.4 22.7 10.7 30.4 19.3 11.1 25.9 17.7 8.2 32.0 19.1 12.9 19.6 14.6 5.0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Absorbing material layout for vertical barrier with constraint condition 0.7 5

1

Absorbing material layout for half-Y barrier with constraint condition 0.7 5

1

4.5

0.9

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

12

0

13

0

8

(a) vertical noise barrier

9

10

11

12

(b) Half-Y-shaped noise barrier

Absorbing material layout for T-shape barrier with two wells

5

1

Absorbing material layout for circle top structure with constraint condition 0.7 5

1

4.5

0.9

4.5

0.9

4

0.8

4

0.8

3.5

0.7

3.5

0.7

3

0.6

3

0.6

2.5

0.5

2.5

0.5

2

0.4

2

0.4

1.5

0.3

1.5

0.3

1

0.2

1

0.2

0.5

0.1

0.5

0.1

0

0

8

9

10

11

12

13

0

0

8

(c) T-shaped noise barrier with two wells

9

10

11

12

(d) veritical noise barrier with a circle top structure

Fig. 20: The final material distribution of several types of noise barriers at 100 Hz.

31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Future work will further extend the proposed algorithm into 3D acoustic sensitivity analysis and optimization analysis of practical engineering problems, and discuss the interplay between wave length and the final shape or material distribution. Acknowledgements This study was funded by the National Natural Science Foundation of China (NSFC) under Grant No. 11702238, Henan Provincial Department of Science and Technology Research under Grant no. 172102210453, Key Scientific Research Project of Henan University under Grant no. 17A560009, and Nanhu Scholars Program for Young Scholars of XYNU. [1] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39) (2005) 4135–4195. [2] T. Takahashi, T. Matsumoto, An application of fast multipole method to isogeometric boundary element method for laplace equation in two dimensions, Engineering Analysis with Boundary Elements 36 (2012) 1766–1775. [3] J. Gu, J. Zhang, G. Li, Isogeometric analysis in bie for 3-d potential problem, Engineering Analysis with Boundary Elements 36 (2012) 858–865. [4] Y. Gong, C. Dong, X. Qin, An isogeometric boundary element method for three dimensional potential problems, Journal of Computational and Applied Mathematics 313 (2017) 454–468. [5] F. Auricchio, L. Beirao da Veiga, A. Buffa, C. Lovadina, A. Reali, G. Sangalli, A fully “locking-free” isogeometric approach for plane linear elasticity problems: a stream function formulation, Computer Methods in Applied Mechanics and Engineering 197 (1) (2007) 160–172. [6] D. Schillinger, L. Ded´e, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, T. J. R. Hughes, An isogeometric designthrough-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Computer Methods in Applied Mechanics and Engineering 249-252 (2012) 116–150. [7] R. N. Simpson, S. P. A. Bordas, J. Trevelyan, T. Rabczuk, A two-dimensional isogeometric boundary element method for elastostatic analysis, Computer Methods in Applied Mechanics and Engineering 209-212 (2012) 87– 100. [8] M. Yoon, S. Cho, Isogeometric shape design sensitivity analysis of elasticity problems using boundary integral equations, Engineering Analysis with Boundary Elements 66 (2016) 119–128. [9] Y. Bai, C. Dong, Z. Liu, Effective elastic properties and stress states of doubly periodic array of inclusions with complex shapes by isogeometric boundary element method, Composite Structures 128 (2015) 54–69. [10] X. Peng, E. Atroshchenko, P. Kerfriden, S. Bordas, Linear elastic fracture simulation directly from cad: 2d nurbsbased implementation and role of tip enrichment, International Journal of Fracture 204 (2017) 55–78. [11] A. Buffa, G. Sangalli, R. Vazquez, Isogeometric methods for computational electromagnetics: B-spline and Tspline discretizations, arXiv preprint arXiv:1209.2048. [12] R. V´azquez, A. Buffa, Isogeometric analysis for electromagnetic problems, IEEE Transactions on Magnetics 46 (8) (2010) 3305–3308. [13] Y. Bazilevs, V. M. Calo, T. J. R. Hughes, Y. Zhang, Isogeometric fluid-structure interaction: theory, algorithms, and computations, Computational mechanics 43 (1) (2008) 3–37. [14] J. A. Evans, T. J. R. Hughes, Isogeometric divergence-conforming B-splines for the steady Navier–Stokes equations, Mathematical Models and Methods in Applied Sciences 23 (08) (2013) 1421–1478. [15] L. Heltai, M. Arroyo, A. DeSimone, Nonsingular isogeometric boundary element method for stokes flows in 3D, Computer Methods in Applied Mechanics and Engineering 268 (2014) 514–539. [16] J. A. Cottrell, A. Reali, Y. Bazilevs, T. J. R. Hughes, Isogeometric analysis of structural vibrations, Computer Methods in Applied Mechanics and Engineering 195 (41) (2006) 5257–5296. [17] M. J. Peake, J. Trevelyan, G. Coates, Extended isogeometric boundary element method (XIBEM) for twodimensional helmholtz problems, Computer Methods in Applied Mechanics and Engineering 259 (2013) 93–102. [18] R. N. Simpson, M. A. Scott, M. Taus, D. C. Thomas, H. Lian, Acoustic isogeometric boundary element analysis, Computer Methods in Applied Mechanics and Engineering 269 (2014) 265–290. [19] R. Simpson, Z. Liu, Acceleration of isogeometric boundary element analysis through a black-box fast multipole method, Engineering Analysis with Boundary Elements 66 (2016) 168–182. [20] M. Peake, J. Trevelyan, G. Coates, Extended isogeometric boundary element method (xibem) for three-dimensional medium-wave acoustic scattering problems, Computer Methods in Applied Mechanics and Engineering 284 (2015) 762–780.

32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[21] S. Marburg, Developments in structural-acoustic optimisation for passive noise control, Archives of Computational Methods in Engineering 27 (2002) 291–370. [22] L. Chen, C. Zheng, H. Chen, A wideband fmbem for 2d acoustic design sensitivity analysis based on direct differentiation method, Computational Mechanics 52 (6) (2013) 631–648. [23] L. Chen, C. Zheng, H. Chen, Fem/wideband fmbem coupling for structural-acoustic design sensitivity analysis, Computer Methods in Applied Mechanics and Engineering 276 (12) (2014) 1–19. [24] L. Chen, H. Chen, C. Zheng, S. Marburg, Structural-acoustic sensitivity analysis of radiated sound power using a finite element/discontinuous fast multipole boundary element scheme, International Journal for Numerical Methods in Fluids 82 (12) (2016) 858–878. [25] L. Chen, L. Liu, W. Zhao, H. Chen, 2d acoustic design sensitivity analysis based on adjoint variable method using different types of boundary elements, Acoustics Australia 44 (12) (2016) 343–357. [26] 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, Journal of Computational Acoustics 25 (2017) 1750003. [27] C. Zheng, C. Zhang, C. Bi, H. Gao, L. Du, H. Chen, Coupled fecbe method for eigenvalue analysis of elastic 359 structures submerged in an infinite fluid domain, International Journal for Numerical Methods in Engineering 110 (2017) 163–185. [28] K. Li, X. Qian, Isogeometric analysis and shape optimisation via boundary integral, Computer-Aided Design 43 (11) (2011) 1427–1437. [29] H. Lian, P. Kerfriden, S. Bordas, Implementation of regularized isogeometric boundary element methods for gradient-based shape optimisation in two-dimensional linear elasticity, International Journal for Numerical Methods in Engineering 106 (12) (2016) 972–1017. [30] H. Lian, P. Kerfriden, S. Bordas, Shape optimization directly from cad: An isogeometric boundary element approach using t-splines, Computer Methods in Applied Mechanics and Engineering 317 (2017) 1–41. [31] K. Kostasa, A. Ginnis, C. Politis, P. Kaklis, Shape-optimization of 2d hydrofoils using an isogeometric bem solver, Computer-Aided Design 82 (2017) 79–87. [32] S. Lee, M. Yoon, S. Cho, Isogeometric topological shape optimization using dual evolution with boundary integral equation and level sets, Computer-Aided Design 82 (2017) 88–99. [33] M. Bendsøe, Optimal shape design as a material distribution problem, Structural optimisation 1 (4) (1989) 193–202. [34] R. Toledo, J. Azn¢rez, O. Maeso, D. Greiner, Optimization of thin noise barrier designs using evolutionary algorithms and a dual bem formulation, Journal of Sound and Vibration 334 (2015) 219–238. [35] R. Toledo, J. Azn´arez, D. Greiner, O. Maeso, Shape design optimization of road acoustic barriers featuring topedge devices by using genetic algorithms and boundary elements, Engineering Analysis with Boundary Elements 63 (2016) 49–60. [36] W. Zuo, K. Saitou, Multi-material topology optimization using ordered simp interpolation, Structural and Multidisciplinary Optimization 55 (2017) 1–15. [37] H. Long, Y. Hu, X. Jin, H. Yu, H. Zhu, An optimization procedure for spot-welded structures based on simp method, Computational Materials Science 117 (2016) 602–607. [38] O. Sigmund, Morphology-based black and white filters for topology optimisation. structural and multidisciplinary optimisation, Structural and Multidisciplinary optimisation 33 (4) (2007) 401–424. [39] S. Liu, Q. Li, W. Chen, R. Hu, L. Tong, H-dgtp—a heaviside-function based directional growth topology parameterization for design optimization of stiffener layout and height of thin-walled structures, Structural and Multidisciplinary Optimization 52 (5) (2015) 903–913. [40] S. Xu, Y. Cai, G. Cheng, Volume preserving nonlinear density filter based on heaviside functions, Structural and Multidisciplinary optimisation 41 (4) (2010) 495–505. [41] L. Kai, Z. Hongwei, A modified optimality criterion method for gray elements suppression, Journal of ComputerAided Design & Computer Graphics 22 (12) (2010) 2197–2201. [42] L. Shang, G. Zhao, Optimality criteria-based topology optimisation of a bi-material model for acoustic-structural coupled systems, Engineering optimisation 48 (6) (2016) 1060–107. [43] B. Hassani, M. Khanzadi, An isogeometrical approach to structural topology optimization by optimality criteria, Structural?and?Multidisciplinary Optimization 45 (2012) 223–233. [44] T. Ishizuka, K. Fujiwara, Performance of noise barriers with various edge shapes and acoustical conditions, Applied Acoustics 65 (2) (2004) 125–141. [45] A.J. Burton and G.F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary-value problems, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 323 (1971) 201–210. [46] S. Marburg, Six boundary elements per wavelength. is that enough ?, Journal of Computational Acoustics 10 (2002) 25–51. [47] G. Monegato, Numerical evaluation of hypersingular integrals, Journal of Computational Applied Mathematics 50 (1994) 9–31.

33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[48] M. Junger, D. Feit, Sound, structures, and their Interaction, Vol. 225, MIT press Cambridge, 1986.

34