A level-set-based topology optimisation for acoustic–elastic coupled problems with a fast BEM–FEM solver

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521 www.elsevier.com/locate/cma

A level-set-based topology optimisation for acoustic–elastic coupled problems with a fast BEM–FEM solver Hiroshi Isakari a,∗ , Toyohiro Kondo b , Toru Takahashi a , Toshiro Matsumoto a a Nagoya University, Japan b Toyota Industries Corporation, Japan

Received 8 June 2016; received in revised form 19 October 2016; accepted 4 November 2016 Available online 16 November 2016

Abstract This paper presents a structural optimisation method in three-dimensional acoustic–elastic coupled problems. The proposed optimisation method finds an optimal allocation of elastic materials which reduces the sound level on some fixed observation points. In the process of the optimisation, configuration of the elastic materials is expressed with a level set function, and the distribution of the level set function is iteratively updated with the help of the topological derivative. The topological derivative is associated with state and adjoint variables which are the solutions of the acoustic–elastic coupled problems. In this paper, the acoustic–elastic coupled problems are solved by a BEM–FEM coupled solver, in which the fast multipole method (FMM) and a multi-frontal solver for sparse matrices are efficiently combined. Along with the detailed formulations for the topological derivative and the BEM–FEM coupled solver, we present some numerical examples of optimal designs of elastic sound scatterer to manipulate sound waves, from which we confirm the effectiveness of the present method. c 2016 Elsevier B.V. All rights reserved. ⃝

Keywords: Topology optimisation; Topological derivative; Level set method; Acoustic–elastic coupled problem; BEM–FEM coupled solver; Wave scattering

1. Introduction Computer simulations play an important role in modern product manufacturing process in various engineering fields. The concept of computer aided engineering (CAE) is now widely accepted in some industries, which utilises numerical analysis to aid in tasks to evaluate the performance of engineering products. With the help of CAE, the total cost and period for product developments have considerably been reduced. In these days, use of computer simulations is not limited to performance evaluation, but is extended to design process. As such an attempt, we can mention a structural optimisation, which is classified into sizing, shape and topology optimisation [1]. The topology optimisation is considered as the most powerful design method in structural optimisations since it can design not only the shape ∗ Corresponding author.

E-mail address: [email protected] (H. Isakari). http://dx.doi.org/10.1016/j.cma.2016.11.006 c 2016 Elsevier B.V. All rights reserved. 0045-7825/⃝

502

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

but also the topology of devices, i.e., the topology optimisation allows a nucleation of a new material and/or hole in its process. Hence, the obtained optimal design by the topology optimisation is less affected by an initial guess than the other structural optimisations. After a pioneering work in [1], the topology optimisation is intensively studied mainly in the field of structural mechanics in order to design light but stiff structural members [1–3]. Recently, one of the main interests in the topology optimisation community is to widen the applicability of the topology optimisation to problems in various engineering fields other than structural mechanics, such as thermal problems [4–6], fluid problems [7–9], elastodynamic problems [10–13], electromagnetic wave problems [14–18], and so on. There are also many efforts to extend the topology optimisation into various design problems in acoustics such as an acoustic horn to maximise sound level [19], an acoustic metamaterial which realises a material with negative effective bulk modulus [20], and a poroelastic sound-proofing material [21–25]. We think, however, the applicability of these existing methods for industrial designs is still limited because these topology optimisations use the finite element method (FEM) to solve boundary value problems involved in sensitivity analysis. Since the acoustic problem is often defined in an unbounded domain, the unbounded domain is approximated with a large one in FEM, which leads to an unexpectedly large scale problem. Especially, this problem may be crucial when sound sources and/or observation points are located far away from the design object because, for FEM, both sources and observation points as well as the design objects must be covered by finite element mesh. Further, an artificial boundary condition such as perfect matched layer (PML) is required to make sure that the scattered wave does not reflect on the truncated boundary. On the other hand, when the boundary element method (BEM) is employed to solve wave scattering problems, only the boundary of the domain is needed to be discretised, which considerably reduces the number of elements. Even when sources and/or observation points are far away from the design objects, the boundary element mesh is required only on the boundary. Also, the scattered fields by the BEM automatically satisfy the radiation condition, i.e., no artificial boundary condition is required to deal with the unbounded domain with the BEM. Thus, the BEM is more suitable for topology optimisations in wave problems than the FEM. As pioneering works on BEM-based topology optimisations, we can mention Abe et al. [26] and Du and Olhoff [27], In the first one, they have solved a two-dimensional shape optimisation problem by the BEM to design a sound barrier. In the second one, they have solved a topology optimisation problem for a noise reduction device from vibrating structures, in which they use an approximated boundary integral formulation for high frequency problems. The applicability of these BEM-based method is limited to either two-dimensional problem [26] or high frequency problem [27] because a naive BEM for three-dimensional realistic scale problems is too expensive. It is inevitable to accelerate the BEM by, for example, fast multipole method (FMM) [28,29], H-matrix algebra [30] and fast direct solver [31] for topology optimisation in sound problems. In order to realise a topology optimisation for three-dimensional realistic design problems of wave devices, we have been investigating level-set-based topology optimisations with the BEM accelerated by the FMM. In our methodology, a candidate for optimal configuration is expressed with a level set function which is iteratively updated with the help of the topological derivative [32–35] to find an optimal distribution of scatterers. Since we have adopted the level set method and BEM, our methods are greyscale free, i.e. the boundary of the design object is clearly expressed. In our methods, the boundary elements are re-meshed from the distribution of the level set function at the every step of the optimisation. In [32], we have investigated a topology optimisation for rigid materials to minimise sound pressure on some observation points. We have extended the methodology to find an optimal allocation of sound absorbers in [36], in which a sound absorbing material is modelled with the impedance boundary condition. The impedance boundary condition is, however, not appropriate to model sound absorbing material in some applications. For example, in analysis with the impedance boundary condition, penetrated sound waves in the sound absorbing materials cannot explicitly be observed, and vibrations in the sound absorber itself are neglected. In this study, to further enhance the applicability of our methodology, we present a level-set-based topology optimisation in acoustic–elastic coupled problems, with which the vibrations of sound scatterers made of elastic materials are explicitly evaluated. In order to solve the acoustic–elastic coupled problem, we adopt a BEM–FEM solver which solves the acoustic and elastic field by BEM and FEM, respectively. This choice is reasonable since, with our settings, the elastic material is in a bounded domain while the acoustic host matrix is an unbounded domain. Although the acoustic–elastic coupled problem can appropriately be solved by the BEM [37,38], we use the BEM–FEM solver [39] since the solver can naturally be extended to deal with elastic material other than the isotropic one such as anisotropic material and Biot’s poroelastic material [40,41]. So far, some fast techniques [42]

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

503

for the BEM–FEM solver are proposed. We here propose another acceleration technique for the BEM–FEM coupled solver, in which FMM and a multi-frontal solver for sparse matrices are efficiently combined. It is worth mentioning that several attempts have already been made to solve the topology optimisation problem for acoustic–elastic problems with a multi-phase modelling [23–25]. In these methods, acoustic and elastic materials are expressed with “two-phase” materials. One of the most distinguished aspects of the present method from the multi-phase modelled method is that the boundary between elastic and acoustic material is explicitly expressed by boundary elements, with which we can evaluate in an accurate manner the boundary values of acoustic properties on the boundary such as sound pressure, sound flux, etc. The rest of the paper is organised as follows. After presenting the statement of the acoustic–elastic coupled problem and related optimisation problem in Section 2.1, we derive the relevant topological derivatives in Section 2.2. In Section 2.3, we propose a new fast method with the FMM and a multi-frontal solver to solve algebraic equations stemmed from the BEM–FEM coupling for the acoustic–elastic coupled problem. After briefly reviewing a formulation and algorithm of the present optimisation method in Sections 2.4 and 2.5, we present some numerical examples which verify the efficiency of the proposed methods in Section 3. Specifically, we check the computational cost for the present BEM–FEM coupled solver in Section 3.1, numerically verify the topological derivatives in Section 3.2, and present two optimal designs of elastic sound scatterers in Section 3.3. In Section 4, we conclude the paper, and discuss the remaining issues to be addressed in the future. In the appendix, we briefly show the meshing algorithm from the distribution of the level set function, which is involved in the proposed optimisation method. 2. Formulations 2.1. Statement of the optimisation problem in acoustic–elastic coupled problems We consider a three-dimensional acoustic–elastic coupled problem in which an incident sound wave from sources on xisrc (i = 1, . . . , M src ), where M src is the number of the sources, is scattered by elastic scatterers filled in a bounded domain Ω c . The sound field in Ω := R3 \ Ω c and the transmitted elastic field in Ω c are governed by the following boundary value problem: p, j j (x) + kf2 p (x) +

src M 

  src Asrc =0 m δ x − xm

x ∈ Ω,

(1)

m=1

σ ji, j (x) + ρs ω2 u i (x) = 0 ti (x) + p (x) n i (x) = 0 ∂ p(x) q (x) := = ρf ω2 u i (x) n i (x) ∂n(x) |x| (q (x) − ikf p (x)) → 0

x ∈ Ω c, x ∈ Γ := ∂Ω

(2) ∩ ∂Ω c ,

(3)

x ∈ Γ,

(4)

as |x| → ∞,

(5)

where is the intensity of mth sound source, and δ is the Dirac delta, p is the sound pressure, σi j is the stress, u i is the displacement, and ti is the traction defined as Asrc m

ti (x) = Ci jkℓ u k,ℓj (x)n j (x),

(6) Ωc

where n(x) is the exterior normal vector with respect to on x ∈ Γ , and Ci jkℓ is the elastic tensor which, for the case that an isotropic material is concerned, has the following representation: Ci jkℓ = λδi j δkℓ + µδik δ jℓ + µδiℓ δ jk ,

(7)

where λ and µ are the Lam´e constants, and δi j is the Kronecker delta. Also, ρf and ρs are the densities for acoustic and elastic materials, respectively, and kf is the wave number for the acoustic wave defined as  ρf , (8) kf = ω Λf where Λf is the bulk modulus of the acoustic material, and ω is the angular frequency with which the time dependency of physical quantities is assumed to be e−iωt .

504

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

Fig. 1. Settings for the topology optimisation in acoustic–elastic coupled problem.

Fig. 2. An infinitesimal elastic inclusion Ωε is introduced in the design domain D.

The boundary value problem in Eqs. (1)–(5) is not uniquely solvable for certain frequencies, called eigenfrequencies [38,43]. With the eigenfrequency ω, there exists a non-trivial solution of the homogeneous boundary src ) in (1). For a real eigenfrequency called the “Jones frequency”, value problem, i.e. for Asrc m = 0 (m = 1, . . . , M the corresponding non-trivial solution satisfies p = 0 in Ω and u · n = 0, t = 0 on Γ . It is also known that complex eigenfrequencies with negative imaginary part may exist, which may affect the accuracy of numerical methods when the relevant frequency is close (even when not identical) to one of the eigenfrequencies. We henceforth assume that the frequency ω is real and far away from any of the eigenfrequencies. Our optimisation problem is defined as to find an optimal distribution of elastic material(s) Ω c ⊂ D (Fig. 1) which minimises the objective function J defined as follows: J=

obs M 

f

   p xobs , m

(9)

m=1

where D, which is so called design domain, is bounded, and xobs m ̸∈ D is mth observation point on which the sound pressure is evaluated, and M obs is the number of the observation points. Also, f is a real-valued functional of the sound pressure p. 2.2. The topological derivative for acoustic–elastic coupled problems In this subsection, we present the topological derivative for the objective function in Eq. (9). We here derive the topological derivative TΩ which characterises the sensitivity of J to an appearance of an infinitesimal spherical elastic material Ωε in the acoustic matrix Ω . The topological derivative TΩ c with respect to appearance of an acoustic material Ωε in the elastic inclusion Ω c can similarly be obtained. Let us assume that an infinitesimal spherical elastic material Ωε , whose elastic properties are identical to those of Ω c , appears in D (Fig. 2). We henceforth denote the centre and the radius of the infinitesimal sphere Ωε as x0 and ε, respectively. Due to the appearance of the small elastic inclusion Ωε , the functions p, q, u i , ti , σi j in Eqs. (1)–(5) suffer from perturbations which are henceforth denoted as in Table 1. The perturbations are governed by the following boundary value problem: δp, j j (x) + kf2 δp (x) = 0 δσ ji, j (x) + ρs ω δu i (x) = 0 2

x ∈ Ω \ Ω ε,

(10)

c

(11)

x∈Ω ,

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

505

Table 1 Perturbed physical quantities due to the appearance of the elastic inclusion Ωε . In Ω \ Ω ε

In Ω c

In Ωε

Sound pressure: p + δp Sound flux: q + δq –

Displacement u i + δu i Stress σi j + δσi j Traction: ti + δti

Displacement uˆ i Stress σˆ i j Traction: tˆi

δti (x) + δp (x) n i (x) = 0

x ∈ Γ,

(12)

2

δq (x) = ρf ω δu i (x) n i (x)

x ∈ Γ,

(13)

σˆ ji, j (x) + ρs ω2 uˆ i (x) = 0

x ∈ Ωε ,

(14)

tˆi (x) + ( p + δp) (x) n i (x) = 0

x ∈ Γε ,

(15)

(q + δq) (x) = ρf ω uˆ i (x) n i (x) |x| (δq (x) − ikf δp (x)) → 0

x ∈ Γε ,

(16)

as |x| → ∞,

(17)

2

where the exterior normal vector on Γε := ∂Ωε ∩ ∂Ω is defined with respect to Ωε . The objective function in (9) also suffers from a perturbation δ J due to the appearance of Ωε as       obs M   ∂ f p xobs m . δJ = ℜ δp xobs (18) m ∂ p m=1 Note that the direct evaluation of δ J with Eq. (18) is impractical since it involves the perturbation of the sound pressure δp on all observation points xobs m which is the solution of the boundary value problem in Eqs. (10)–(17). In this paper, we use the adjoint variable method to evaluate δ J without going through δp(xobs m ). The adjoint problem is defined as follows:     obs M   ∂ f p xobs m 2 p˜ , j j (x) + kf p˜ (x) + δ x − xobs = 0 x ∈ Ω, (19) m ∂p m=1 σ˜ ji, j (x) + ρs ω2 u˜ i (x) = 0 t˜i (x) + p˜ (x) n i (x) = 0 ∂ p(x) ˜ = ρf ω2 u˜ i (x) n i (x) q˜ (x) := ∂n |x| (q(x) ˜ − ikf p(x)) ˜ →0

x ∈ Ω c,

(20)

x ∈ Γ,

(21)

x ∈ Γ,

(22)

as |x| → ∞,

(23)

where p, ˜ q, ˜ u˜ i , σ˜ i j and t˜i are the adjoint sound pressure, the adjoint sound flux, the adjoint displacement, the adjoint stress and the adjoint traction, respectively. According to the reciprocity of the state variable p and the adjoint variable p˜ in Ωε , we have the following identity:  ˜ − pq) ˜ dΓ = 0. (24) ( pq Γε

A similar procedure to uˆ i and p˜ in Ωε together with Eqs. (15), (16), (19) and (24) gives the following identity:     ˜ − pδq) ˜ dΓ = Λf uˆ i,i p˜ , j j − p˜ ,i j σˆ ji + ρs ω2 p˜ ,i uˆ i − ρf ω2 uˆ i p˜ ,i dΩ . (qδp Γε

(25)

Ωε

We also have the following reciprocal relation between δu i and p˜ in Ω c combined with the boundary condition in Eq. (12):    δu i t˜i + δp u˜ i n i dΓ = 0. (26) Γ

506

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

With the reciprocal theorem between p and p˜ in Ω \ Ωε , and Eqs. (13), (21), (22), (25) and (26), we can evaluate δ J as follows:     δJ = ℜ Λf uˆ i,i p˜ , j j − p˜ ,i j σˆ ji + ρs ω2 p˜ ,i uˆ i − ρf ω2 uˆ i p˜ ,i dΩ . (27) Ωε

The expression (27) can further be simplified with the help of the Gauss theorem as      2 δJ = ℜ − ρf ω uˆ r p˜ + tˆi p˜ ,i dΓ ,

(28)

Γε

where uˆ r := uˆ i n i is the radial component of the displacement on Γε . Note that the expression in Eq. (28) does not involve the perturbations of the state variables on the observation points. We just need to evaluate the perturbations on the boundary Γε of an infinitesimal elastic inclusion. In the following, we evaluate the asymptotic behaviour of δ J in Eq. (28) as ε → 0. To this end, p˜ and its gradient are respectively expanded as p(x) ˜ = p˜ 0 + ε p˜ ,0j n j (x) + o (ε)

x ∈ Γε ,

(29)

+ ε p˜ ,i0 j n j (x) + o (ε)

x ∈ Γε ,

(30)

p˜ ,i (x) =

p˜ ,i0

where f 0 denotes f 0 = f (x0 ). The state variables p and q can also be expanded as q(x) =

p(x) = p 0 + εpi0 n i (x) + o (ε)

x ∈ Γε ,

(31)

p,i0 n i (x) + εp,i0 j n i (x)n j (x) + o (ε)

x ∈ Γε ,

(32)

respectively. The asymptotic expansions of the radial displacement uˆ r and the traction tˆi on Ωε in (28) are then evaluated. Although uˆ r and tˆi are the solution of the boundary value problem in (10)–(17), it is sufficient to solve the following approximated one [33]: x ∈ Ωεc := R3 \ Ωε ,

(33)

σˆ j,i j (x) + ρs ω uˆ i (x) = 0

x ∈ Ωε ,

(34)

(q + δq) (x) = ρf ω2 uˆ i (x) n i (x) δti (x) + ( p + δp) (x) n i (x) = 0

x ∈ Γε ,

(35)

x ∈ Γε ,

(36)

as |x| → ∞,

(37)

δp, j j (x) + kf2 δp (x) = 0 2

|x| (δq(x) − ikf δp(x)) → 0

since the asymptotic behaviour as ε → 0 is now concerned. For the sake of reference, we rewrite the boundary conditions (35) and (36) as follows: q (x) + δq (x) = ρf ω2 uˆ r (x) p (x) + δp (x) + σˆ rr (x) = 0 σˆ r θ (x) = σˆ r φ (x) = 0

x ∈ Γε ,

(38)

x ∈ Γε ,

(39)

x ∈ Γε ,

(40)

where σˆ rr , σˆ r θ and σˆ r φ denote the polar representations of the stress σˆ i j . The solutions of the boundary value problem (33), (34), (37)–(40) can be written in terms of spherical functions [44] as follows: δp (x) =

∞  n  n=0 m=−n ∞  n 

dnm h n(1) (kfr ) Pnm (cos θ )eimφ

x ∈ Ωεc ,

(1)

∂h n (kfr ) m Pn (cos θ )eimφ x ∈ Ωεc , ∂r m=−n n=0   ∞  n  1 cnm n m n uˆ r (x) = an U1 (r ) − U3 (r ) Pnm (cos θ )eimφ x ∈ Ωε , r k T n=0 m=−n   ∞  n  2µ cnm n m n σˆ rr (x) = a T (r ) − T (r ) Pnm (cos θ )eimφ x ∈ Ωε , n 11 13 2 k r T n=0 m=−n δq (x) =

(41)

dnm

(42) (43) (44)

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

507

   ∞  n  n+m m im m 2µ m n m θ) − a T (r ) n cot θ P θ + bnm T42 (r ) P P (cos θ ) (cos (cos ) n 41 n n−1 2 sin θ sin θ n r n=0 m=−n   n+m m m n m (45) P + cn T43 (r ) n cot θ Pn (cos θ) − (cos θ ) eimφ x ∈ Ωε , sin θ n−1  ∞  n  2µ im m n σˆ r φ (x) = P (cos θ) anm T51 (r ) 2 sin θ n r n=0 m=−n σˆ r θ (x) =

 im  m n cos θ Pnm (cos θ ) − (n + m) Pn−1 (cos θ ) sin θ  cnm n im m + T (r ) P (cos θ ) eimφ x ∈ Ωε , kT 53 sin θ n n − bnm T52 (r )

(46)

(1)

where h n is the nth spherical Hankel function of the first kind, Pnm is the associated Legendre function defined as follows: dm Pn (x) (m ≥ 0), dxm (n − m)! m Pn−m (x) = (−1)m P (x) (m ≥ 0), (n + m)! n Pnm (x) = (1 − x 2 )m/2

(47) (48)

where Pn is the Legendre polynomial, and (r, θ, φ) represents the spherical coordinate of the point x. Also, U1n , U3n , n , T n , T n , T n , T n , T n , T n and T n are the functions defined as follows: T11 13 41 42 43 51 52 53 U1n (r ) = n jn (kLr ) − kLr jn+1 (kLr ) ,

(49)

U3n (r )

(50)

= n (n + 1) jn (kTr ) ,   1 n T11 (r ) = n 2 − n − kT2 r 2 jn (kLr ) + 2kLr jn+1 (kLr ) , 2 n T13 (r ) = n (n + 1) ((n − 1) jn (kTr ) − kTr jn+1 (kTr )) ,

(51)

n T41 (r )

= (n − 1) jn (kLr ) − kLr jn+1 (kLr ) , 1 n T42 (r ) = r ((n − 1) jn (kTr ) − kTr jn+1 (kTr )) , 2   1 n T43 (r ) = n 2 − 1 − kT2 r 2 jn (kTr ) + kTr jn+1 (kTr ) , 2 n n T51 (r ) = T41 (r ),

(53)

n n T52 (r ) = T42 (r ),

(57)

n (r ) T53

(58)

=

n T43 (r ),

(52)

(54) (55) (56)

where jn is the nth spherical Bessel function, and kL and kT are respectively the wave numbers for the longitudinal and the transverse wave which have the following expressions:  ρs kL = ω , (59) λ + 2µ  ρs kT = ω . (60) µ Also, anm , bnm , cnm and dnm ∈ C are coefficients of the spherical expansion. Substituting (45), (46) into (40) gives the following relation for anm , bnm and cnm : bnm = 0, n T41 (ε)anm +

n (ε) T43 cm kT n

= 0,

(61) (62)

508

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

by which Eqs. (43) and (44) are reduced as   ∞  n n (ε)  T41 anm n n U1 (r ) − n uˆ r (x) = U3 (r ) Ynm (θ, φ) , r T (ε) 43 n=0 m=−n   n ∞  n (ε)  T41 2µ m n n T (r ) − σˆ rr (x) = a T (r ) Ynm (θ, φ) , 11 n (ε) 13 2 n T r 43 n=0 m=−n

(63) (64)

respectively. By substituting (32), (42) and (63) into (38), substituting (31), (41) and (64) into (39), and exploiting the orthogonal properties of the spherical harmonics, we obtain a system of algebraic equations to determine dnm and anm , from which we obtain the asymptotic expansions of uˆ r and σˆ rr as kL2 3 0  p 0 ε + o (ε) ,  n + p (x) j ,j 2µkT2 + ρf ω2 µ 4kL2 − 3kT2   0 p, j j p,i0 j 3 0 0 σˆ rr (x) = − p + − p, j n j (x) − n i (x)n j (x) ε + o (ε) , 9 3 2µkT2 + ρf ω2 uˆ r (x) =

(65) (66)

respectively. The asymptotic expansion of the traction tˆi can easily be calculated from Eqs. (40) and (66). With the evaluated asymptotic expansions of the state variables, we can evaluate the asymptotic expansion of δ J as      4 3 3 (ρs − ρf ) 0 0 Λs − Λf 2 0 0 δ J = ℜ πε p, j p˜ , j − ρf ω p p˜ + o ε3 , (67) 3 2ρs + ρf Λs Λf where Λs is the bulk modulus for elastic material Ω c defined as 2 Λs = λ + µ. 3

(68)

The topological derivative is defined as the coefficient of the leading term of the asymptotic expansion (67) of the objective function J [32–35] as follows: δ J = TΩ (x) v (ε) + o (v (ε)) ,

(69)

where v (x) is a monotonically increasing function in x > 0. By comparing Eqs. (67) and (69), we obtain the topological derivative TΩ with respect to an appearance of an elastic material Ωε in the fluid matrix Ω as follows:   3 (ρs − ρf ) Λs − Λf 2 TΩ (x) = ℜ p, j (x) p˜ , j (x) − ρf ω p (x) p˜ (x) , (70) 2ρs + ρf Λs Λf where v(ε) in Eq. (69) is chosen as v (ε) = 34 π ε 3 . The topological derivative TΩ c related to appearance of an acoustic material Ωε in the elastic inclusion Ω c can similarly be obtained as follows:    TΩ c (x) = ℜ ρf ω2 (A − B) σ˜ ii (x)σii (x) + 3B σ˜ i j (x)σi j (x) − (ρs − ρf ) ω2 u˜ i (x)u i (x) , (71) where the coefficients A and B are defined as 3 (Λf − Λs ) (λ + 2µ)  , (3λ + 2µ) 12µ (λ + Λf ) + 9λ2 + 4µ2 5(λ + 2µ) B= , 2µ (9λ + 14µ) A=

respectively. Note that the result in (71) is consistent with the one in Guzina and Chikichev [45].

(72) (73)

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

509

2.3. A fast BEM–FEM solver for acoustic–elastic coupled problems In order to evaluate the topological derivatives in Eqs. (70) and (71), we need to calculate the sound pressure p and its gradient p, j in Ω , and the displacement u i and the stress σi j in Ω c , and their adjoint counterparts. Although these quantities can appropriately be calculated by the boundary element method (BEM), a BEM–FEM (finite element method) coupled solver is utilised in this paper. The proposed solver deals with the acoustic field and the elastic field by BEM and FEM, respectively. This is because, in our future publications, we plan to extend the present topology optimisation for elastic materials other than the isotropic one, e.g. anisotropic material and Biot’s poroelastic material for which the FEM is more suitable than the BEM. In this section, we also present a fast algorithm for the BEM–FEM coupled solver, in which the fast multipole method (FMM) and a multi-frontal solver for sparse matrices are efficiently combined. We present the formulation of the BEM–FEM coupled solver for the forward problem in Eqs. (1)–(5). The adjoint problem in Eqs. (19)–(23) can be solved in the same manner. From Eqs. (1) and (4), we have the following boundary integral equation: obs   M  ∂G (x − y) p (x) obs 2 G (x − y) u ℓ (y) n ℓ (y) dΓ (y) + = G(x − xm ) − ρf ω p (y) dΓ (y) , 2 ∂n(y) Γ Γ m=1

(74)

where G(x) = eikf |x| /4π |x| is the fundamental solution of three dimensional Helmholtz’ equation. The weak form in Ω c and the boundary condition in Eq. (3) give the following equation:    ∗ ∗ 2 u i n i pdΓ − u i, j σ ji dΩ + ρs ω u i∗ u i dΩ = 0, (75) Γ u i∗

Ωc

Ωc

where is a test function. The proposed method solves the system of integral equations (74) and (75). In the Nfe Nbe Γj, Ωe , where Ωe is a tetrahedron, and Γ is divided as Γ = ∪ j=1 discretisation, Ω c is divided as Ω c = ∪e=1 where Γ j is a triangular patch which coincides to a surface of a tetrahedron Ωe . In this study, the sound pressure p is approximated by locally constant functions on Γ j while the displacement u i is approximated by locally linear functions in Ωe . With these settings, a standard collocation for Eq. (74) and the Galerkin discretisation for Eq. (75) in which the test functions are also expanded by the locally linear functions in Ωe give the following system of algebraic equations:      src  I S 2 − D  p = p , (76) u 0 N K − ρs ω 2 M where S and D are the coefficient matrices for single and double layer potentials, respectively. K and M are the finite element stiffness and mass matrices, respectively. Also, the matrix N is stemmed from the first term of Eq. (75). The vectors p and psrc respectively contain the total and the incident sound pressures on the collocation points x j ∈ Γ j ( j = 1, . . . , Nbe ), and the vector u is composed of nodal displacements in Ω c . In the following, we propose a fast solver for the algebraic equation (76). The basic idea of the proposed solver is to combine the fast multipole method (FMM) [28,29] for calculations of the BEM matrices, and a multi-frontal solver to factorise the FEM matrix in by exploiting the block structure in Eq. (76). The second row in (76) can be written as u = (K − ρs ω2 M)−1 Np,

(77)

provided that ω2 is not an eigenvalue of a homogeneous Neumann problem in Ω c . Note that the matrix K − ρs ω2 M is a sparse matrix which can efficiently be factorised by a multi-frontal solver. By substituting Eq. (77) into the first row of Eq. (76), we obtain the following algebraic equation:   −1  I − D + S K − ρs ω2 M N p = psrc . (78) 2 From the Calderon identity, one observes that almost all the eigenvalues of the matrices S and D are close to 0 [46–48]. Thus, the condition number of the coefficient matrix in Eq. (78) is expected to be small, which leads fast convergence of an iterative solver such as GMRES. Matrix–vector products involved in an iterative solver can efficiently be

510

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

performed by either the FMM or a multi-frontal solver. The algorithm for solving the boundary value problem in Eqs. (1)–(5) is summarised as follows: 1. The matrix K − ρs ω2 M in (77) is factorised by a multi-frontal solver, and the factorised matrices are stored. 2. Eq. (78) is solved with an iterative solver to obtain the sound pressure p on the collocation point x j ∈ Γ j ( j = 1, . . . , Nbe ). The product of the coefficient matrix in Eq. (78) and a vector x, which is required in the algorithm of the iterative solver, is performed as follows: (a) The matrix–vector product z := Nx is calculated. In the calculation, the sparsity of the matrix N is exploited. (b) The vector y := (K − ρs ω2 M)−1 z is calculated by solving the algebraic equation (K − ρs ω2 M)y = z with a multi-frontal solver. Note that the coefficient matrix K − ρs ω2 M has already been factorised in 1.   (c) 2I − D x + Sy is calculated by the FMM. 3. The matrix–vector product z′ := Np is calculated. In the calculation, the sparsity of the matrix N is, again, exploited. 4. The displacement u on the nodal point of the finite elements Ωe is obtained by solving (K − ρs ω2 M)u = z′ . Again, note that the coefficient matrix K − ρs ω2 M has already been factorised in 1. 5. The sound pressure p and its gradient p,i are calculated by integral representations on arbitrary points x ∈ Ω . 6. The displacement u and the stress σi j are calculated by interpolation with shape functions on arbitrary points x ∈ Ω c. In the implementation of the FMM, a low-frequency FMM [32] is employed. 2.4. A level-set-based topology optimisation In this subsection, we briefly review a level-set-based methodology to solve the optimisation problem to find an optimal distribution of elastic material(s) Ω c ⊂ D which minimises the objective function in Eq. (9) subject to the constrain conditions in Eqs. (1)–(5) (see also Fig. 1). The reader is referred to the original paper [3] and our previous papers [17,32,49] for further details. In the level set method, domains Ω and Ω c , and its boundary Γ are recognised as Ω c = {ξ | 0 < φ(ξ ) ≤ 1},

(79)

Γ = {ξ | φ(ξ ) = 0},

(80)

Ω = {ξ | −1 ≤ φ(ξ ) < 0},

(81)

respectively. With the level set function φ, the topology optimisation problem is converted into the problem to find an optimal distribution of φ in the design domain D which minimises the objective function in Eq. (9) under the constrains in Eqs. (1)–(5). We explore the optimum distribution of φ using the topological derivative TΩ and TΩ c in Eqs. (70) and (71) from an initial distribution φ0 (ξ ) as follows: ∂φ(ξ , t) = −CTΩ (ξ , t) + τ L 2 ∇ 2 φ(ξ ) ∂t ∂φ(ξ , t) = CTΩ c (ξ , t) + τ L 2 ∇ 2 φ(ξ ) ∂t

for φ(ξ , t) < 0,

(82)

for φ(ξ , t) > 0,

(83)

φ(ξ , 0) = φ0 (ξ ),

(84)

where t represents a fictitious time, C > 0 is a constant, and L is a characteristic length of the design domain D. In the present method, the topological derivatives are used to modify the distribution of the level set function φ. TΩ (ξ , t) in Eq. (82), for example, works to allocate a small elastic scatterer on ξ when TΩ is negative by increasing φ(ξ , t). Also, τ > 0 is a parameter which prescribes the complexity of the geometry of Ω c [3]. The following boundary condition for the time evolution equations (82) and (83) is also defined: φ(ξ , t) = φ(ξ , t) < 0

for ξ ∈ ∂ D and t > 0,

(85)

where φ¯ is a known function. The boundary condition (85) is imposed so that Ω c ⊂ D holds. Thus, the optimisation problem is now converted to the initial–boundary value problem (82)–(85). The initial–boundary value problem is numerically solved with FEM in which the temporal derivative is discretised by the forward difference.

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

511

A rigorous mathematical justification of the optimality of the solution by (82)–(85) may be discussed in our future publications. A relevant discussion can be found, for example, in [50]. 2.5. Algorithm of the present topology optimisation Combining all the techniques presented above, the algorithm of the proposed topology optimisation is summarised as follows: N

D 1. The fixed design domain D is divided into finite elements (voxels) as D = ∪e=1 ΩeD , where ΩeD is a voxel, and ND is the number of the voxels. 2. An initial distribution of the level set function φ(ξ , 0) is given on nodes (lattice points) of the finite elements ND ∪e=1 ΩeD . 3. A set X of the points x such that φ(x) = 0 on the lattice edge is explored. 4. By the procedure presented in the Appendix, the elements of X are appropriately connected, and triangular patches that cover the iso-surface of zero value of the level set function, and the triangular patches are stored in STL format. Nfe Nbe 5. The STL data is converted to the boundary elements Γ = ∪ j=1 Γ j and the finite elements Ω c = ∪e=1 Ωe . For the meshing, NETGEN [51] is used in this study. Note that NETGEN can not only generate boundary/finite meshes but also improve the quality of the meshes. 6. The boundary value problem (1)–(5) (the forward problem) defined in Ω ∪ Ω c is solved by the BEM–FEM coupled solver presented in Section 2.3, and evaluate the objective function in Eq. (9). When the objective function converges, stop. 7. The boundary value problem (19)–(23) (the adjoint problem) defined in Ω ∪ Ω c is solved by the BEM–FEM coupled solver in Section 2.3. Note that the procedure 1 in the algorithm of the BEM–FEM coupled solver can be skipped since matrix K − ρs ω2 M has already been factorised in the forward analysis. Note also that some of the quantities in the FMM algorithms such as tree structures, direct interactions and M2L operators, etc. calculated in the forward analysis can be recycled in the adjoint procedure. 8. With the state and adjoint variables calculated in 6 and 7, the topological derivatives in Eqs. (70) and (71) on the all lattice points expanding D are evaluated. 9. The initial–boundary value problem (82)–(85) is solved by FEM. Go to 3.

3. Numerical examples In this section, we present some numerical examples to confirm the validity and the efficiency of the proposed method. We first state common issues to all the examples to follow: PARDISO routines by Intel MKL is used as the multi-frontal solver involved in the present BEM–FEM solver. GMRES is used as the iterative solver involved in the present BEM–FEM solver. The tolerance of GMRES in the present BEM–FEM solver is set to be 10−5 . Truncation numbers for infinite series in the FMM are numerically determined so that the truncation error is less than 10−5 . • All numerical experiments were run on a PC with Intel Xeon CPU E5-4650 with 32 cores. The code is OpenMP parallelised. • The host acoustic matrix is assumed to be water. This is partly because the interaction between air and metal or resin can appropriately be modelled by sound-hard approximation.

• • • •

3.1. Performance tests of the present BEM–FEM coupled solver In this subsection, we check the performance of the proposed BEM–FEM coupled solver presented in Section 2.3 with benchmark problems. To this end, the algebraic equations (76) corresponding to benchmark problems are solved by conventional direct and iterative solvers as well as by the present solver. We used a LAPACK routine ZGESV and GMRES as the conventional direct and iterative solvers to naively solve Eq. (76), respectively. In the conventional GMRES, the FMM is not employed to accelerate the matrix–vector products involved in the algorithm of GMRES, and the tolerance of the GMRES is set to be 10−12 , which is set by numerical experiments so that the total errors for

512

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

Fig. 3. Computational time for the present and conventional solvers against the number of finite element nodes Np for a sound scattering problem by an elastic sphere.

the sound pressures and the displacements are comparable to those by the present method and the conventional direct solver. As the first benchmark problem, we consider a sound scattering by an elastic sphere whose analytical solution is available in [44]. We set an elastic sphere Ω c := {x | |x| < 0.25} in an acoustic host matrix. We here assume that the elastic inclusion and host acoustic matrix are composed of a tungsten and water, respectively. The parameters for the tungsten are set as the density ρs = 64.85, Young’s modulus E = 174.57 and Poisson’s ratio ν = 0.3, which are normalised by the material parameters of water as the density ρf = 1.0 and the bulk modulus Λf = 1.0. As the incident wave, we used a plane wave propagating in x3 direction with the frequency ω = 1.0. The amplitude of the incident wave is set to be 1.0. With these settings, the algebraic equations (76) are solved by the present method, in which the surface of Ω c is divided into 2456 boundary elements, and Ω c into finite elements of 3007 nodal points. The average of the relative errors for the sound pressures on collocation points and for nodal values of the displacements were 0.0628% and 0.107%, respectively. The accuracy of the present method is comparable to that of a conventional direct solver and a conventional iterative solver. We then discuss the timing. Fig. 3 shows the computational time for the present and conventional solvers against the number of the nodes of finite elements Np . The computational time is measured using an OpenMP run-time library routine omp get wtime. One observes that the computational time for the present method scales as Np , and the present method is the fastest among the tested solvers even when the degrees of freedom is relatively small. The computational complexities of the conventional iterative and direct solvers are O(Np2 ) and O(Np3 ), respectively. We think that the bad condition of the coefficient matrix in Eq. (76) makes the convergence of the conventional iterative solver slow. Indeed, the iteration number for the conventional GMRES is 390 in the case of Np = 3007. On the other hand, the condition of the coefficient matrix in Eq. (78) in the present method is well as discussed in Section 2.3. The iteration number for Eq. (78) is 3 for this case. Table 2 shows computational time for each procedure in the present BEM–FEM solver in the case of Np = 194 536. One finds that calculating the FEM matrices N and K − ρs ω2 M, and factorising K − ρs ω2 M take almost half of the whole computational time. As indicated in Section 2.5, after solving the forward problem in (1)–(5) in the process of the optimisation, we need not repeat these procedures in solving the adjoint problem in Eq. (19)–(23) since the FEM matrices N and K − ρs ω2 M for the adjoint problem are common to those for the forward problem. With these observations, it is expected that the computational time for solving the adjoint problem is as approximately half as that for solving the forward problem. Thus, the present BEM–FEM solver can efficiently be applied to the topology optimisation. As the second benchmark problem, we consider elastic scatterers shown in Fig. 4 to check the applicability of the present BEM–FEM solver to a complex-shaped domain. The elastic scatterers in Fig. 4 are taken from [32], which are obtained in a process of a topology optimisation. The elastic scatterers are expressed with finite elements of

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

513

Table 2 Computational time for each procedure in the present BEM–FEM solver for sound scattering problem by an elastic sphere in the case of Np = 194 536. Procedure

Comp. time [s]

Calculation of matrices N and K − ρs ω2 M Factorisation of K − ρs ω2 M by PARDISO GMRES for solving Eq. (78) Others Sum

24.44 51.78 70.36 6.05 152.63

Table 3 Computational time of the present and a conventional iterative solver for sound scattering problems by elastic scatterers. Method

Sphere (s)

Complicated shape (s)

Present Conventional GMRES

6.04 421.44

52.87 2147659709.14

Fig. 4. A complex-shaped elastic scatterers taken from [32].

Np = 3845. The material parameters and the incident wave are set as in the previous test. In this benchmark problem, we compare the performance of the present solver with that of a conventional iterative solver; GMRES is naively employed to solve the algebraic equation (76). Table 3 shows the computational time for scattering by the complexshaped scatterers and by a sphere with Np = 3007 for comparison. The computational time of the present solver for the complicated shape is as only 9 times as that for the sphere, while the conventional iterative solver for the complicated shape is much slower than the case for a sphere. Thus, the performance of the present solver is less affected by shape of scatterers than that of the conventional iterative solver. This fact can also be seen in Table 4, which shows the iterative number of the GMRES. We confirm that, for complicated domain, the condition number of Eq. (76) can be quite large, while that of Eq. (78) is kept relatively small. This is because the present method exploits the spectral properties of the boundary element matrices S and D according to the Calderon formulae. Thus, the proposed BEM–FEM solver is efficient for complex-shaped domain which is often encountered in topology optimisation. 3.2. Verification of the topological derivative In this subsection, we numerically verify the topological derivative TΩ in Eq. (70) derived in Section 2.2. We here consider to put a spherical elastic material of infinitesimal radius in R3 filled with an acoustic host matrix. We assume that the elastic sphere and the acoustic host matrix are composed of acrylonitrile butadiene styrene (ABS) resin and water, respectively. The ABS resin is a typical elastic material and is now often used as “ink” for 3D printer. The material parameters for ABS resin is normalised as the density ρs = 1.1, the Young modulus E = 1.17 and the

514

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521 Table 4 The number of iterations for the present and a conventional iterative solver for sound scattering problems by elastic scatterers. Method

Sphere

Complicated shape

Present Conventional GMRES

3 390

74 16 092

Fig. 5. Definitions for (left:) the observation points in (86) and (right:) the evaluation points for the topological derivative for the verification of the topological derivative in Eq. (70).

Poisson ratio ν = 0.369 with the ones for water as the density ρf = 1.0 and the bulk modulus Λf = 1.0. The incident wave is assumed to be a plane wave propagating in x3 direction with the angular frequency ω = 1.0. The amplitude of the incident wave is set as 1.0. We define the following sum of the sound norm on the 9 observation points xobs m (defined in the left figure of Fig. 5) as the objective functional: J=

9 1 2 | p(xobs m )| . 2 m=1

(86)

Figs. 6 and 7 show the distribution of the topological derivatives for the objective functional in Eq. (86) on lines ℓ1 and ℓ2 in the right figure of Fig. 5, respectively. In the figures, the “topological differences” D are also plotted, which are defined as follows: JΩ \Ω ε − JΩ , (87) D= v(ε) where v(ε) = 4π ε 3 /3, and JΩ and JΩ \Ω ε represent the objective function before and after a small spherical elastic scatterer of radius ε is introduced, respectively. In the case that ε is small in Eq. (87), D is expected to agree with the topological derivative. We used the BEM for the calculation of JΩ and JΩ \Ω ε with ε = 0.03 whose surface is divided into 2000 boundary elements. One confirms that the topological derivatives derived in this paper agree well with the reference values. 3.3. Optimal designs In this subsection, we show two examples of optimal design of sound scatterers which reduce the sound norm on some preset observation points. In the first example, we consider a hard elastic material which may appropriately be modelled by a rigid one. We show the obtained configuration of the hard elastic material is similar to that of rigid one, with which the validity of the proposed topology optimisation is confirmed. In the second example, we consider a design problem of a soft elastic sound scatterer which cannot be solved by conventional topology optimisation methods. 3.3.1. Hard scatter We explore, with the present topology optimisation, an optimal distribution of hard elastic material in a design domain D := {x | 0 ≤ xi ≤ 2.5 (i = 1, 2, 3)}, which minimises the sound norm on observation points. The observation points are set as (1.25, 1.25, 5.0) and (1.25, 1.25, −2.5), and 20 points on the circles, whose centre are

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

515

Fig. 6. Topological derivatives on the line ℓ1 in Fig. 5.

Fig. 7. Topological derivatives on the line ℓ2 in Fig. 5.

these points and radius is 1.0. The circles are parallel to x1 x2 plane (Fig. 8). Sound sources are set on (5.0, 1.25, 1.25) and (−2.5, 1.25, 1.25) whose intensity and frequency are set as 70 and 2π , respectively. As the initial guess for the elastic material, we used a sphere whose centre and radius are (1.25, 1.25, 1.25) and 0.25, respectively. The elastic material and the host acoustic material are respectively assumed to be a tungsten and water whose material parameters are listed in Section 3.1. In the optimisation algorithm in Section 2.4, we divide the design domain D into 100 × 100 × 100 finite elements, and set τ in Eqs. (82) and (83) as 10−4 . We show the obtained configuration in Fig. 9. For comparison, we also show in Fig. 10 a result of the optimisation problem in which the elastic material is replaced by a rigid one. One finds that the obtained configuration of tungsten is similar to that of rigid material, which is reasonable since the tungsten is quite hard and heavy compared to water. Thus, the tungsten embedded in water can appropriately be approximated by rigid material. This fact can also be confirmed by the expression of the topological derivative in (70). By taking limits as ρf /ρs → 0 and Λf /Λs → 0 (sound-hard limit), the topological derivative TΩ becomes   3 TΩ (x) → ℜ p, j (x) p˜ , j (x) − ρf ω2 p(x) p(x) ˜ , (88) 2 which is identical to the topological derivative for the rigid material [32,33].

516

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

Fig. 8. Settings for a topological optimal design for sound scatterer with hard elastic material. In the figure, the initial configuration of the elastic material is also plotted.

Fig. 9. The obtained configuration of tungsten.

Fig. 10. The obtained configuration of rigid material.

Fig. 11 shows the distribution of the squared sound norm in x2 = 1.25 plane for the case that optimal elastic scatterers are allocated in the design domain D. One observes that the sound norm on the observation points are reduced. This can also be confirmed by the objective function, which was J = 42.19 for the initial configuration in Fig. 8, is J = 2.84 for the optimal configuration. Thus, the present method can reduce the sound norm on the observation points. 3.3.2. Soft scatter We next consider a topology optimisation problem for a soft elastic material to manipulate sound waves. We use the same design domain and the initial configuration of the elastic scatterer as the one in the previous example, and put a point sound source on (1.25, 1.25, −2.5) whose frequency and intensity are 2π and 70, respectively. The objective function is defined as the sum of the sound norm on 33 observation points set on a hemisphere 3 {x | i=1 (xi − 1.25)2 < 2.52 , 1.25 ≤ x3 } (Fig. 12). We here consider a silicone rubber immersed in water. The material parameters for the silicone rubber are normalised as the density ρs = 0.97, the Young modulus E = 0.018

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

517

Fig. 11. The squared sound norm | p|2 around the design domain D in x2 = 1.25 when optimal elastic scatterers as Fig. 9 are allocated.

Fig. 12. Settings for a topological optimal design for sound scatterer with soft elastic material. In the figure, the initial configuration of the elastic material is also plotted.

and the Poisson ratio ν = 0.49 with the ones for water as the density ρf = 1.0 and the bulk modulus Λf = 1.0. Note that the present calculations ignore viscoelastic effect related to silicone rubber. The viscoelastic effect can be considered by complexifying the Lam´e constants [43], which may be addressed in our future publications. In the optimisation algorithm in Section 2.4, we divide the design domain D into 100 × 100 × 100 finite elements, and set τ in Eqs. (82) and (83) as 10−4 as in the previous example. We show the obtained configuration in Fig. 13 and sound norm around the optimal scatterers in x2 = 1.25 in the left figure of Fig. 15. For comparison, we also show the figures in Figs. 14 and 15 (right) results of the optimisation problem in which the silicone rubber is replaced by a rigid material. One observes that the optimal configuration of silicone rubber is different from that of rigid material, and, for both cases, the sound norms on the observation points are small. The objective functions for silicone rubber (resp. rigid material) are reduced from J = 18.49 (resp. 17.82) to 13.76 (resp. 5.35). These figures show that, since the material properties of silicone rubber are quite different from those of rigid material, such soft materials cannot be designed with conventional topology optimisation with rigid approximation. With these examples, we conclude that the present methodology can efficiently design elastic materials to manipulate sound waves.

518

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

Fig. 13. The obtained configuration of silicone rubber.

Fig. 14. The obtained configuration of rigid material.

Fig. 15. The squared sound norm | p|2 around the design domain D in x2 = 1.25 when (left:) optimal elastic scatterers (right:) optimal rigid scatterers are allocated.

4. Conclusion We have developed a new topology optimisation for elastic material to reduce sound level, in which a fast BEM–FEM coupled solver is employed to evaluate the topological derivative. The derivation of the topological derivative for acoustic–elastic coupled problems is described. We have confirmed that the present topology optimisation method can efficiently design elastic sound scatterers. In this paper, we have tested the proposed method in pure elastic problem to reduce sound norm defined at observation points. Applications of the proposed method to viscoelastic design, objective functions defined on boundary such as energy flux on sound absorbing device are, however, still remain to be investigated. Also, we need to investigate the behaviour of the present BEM–FEM solver near resonance (including fictitious resonance). In our future publications, we plan to extend the present method to deal with anisotropic and/or the Biot sound absorbers.

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

(a) Fixed design domain D.

(c) Boundary elements are generated in a voxel.

519

(b) D is divided into voxels.

(d) Boundary elements.

Fig. A.16. Retrieving a boundary element mesh from the distribution of the level set function.

Acknowledgements This work was supported by JSPS Grant-in-Aid for Scientific Research (B) (Grant No. 16H04255) and JSPS Grantin-Aid for Challenging Exploratory Research (Grant No. 15K13856). Appendix. Meshing algorithm We briefly explain, in this appendix, a procedure to retrieve geometrical information on the elastic inclusion Ω c from the distribution of the level set function. As mentioned in Section 2.5, the distribution of the level set function φ in the fixed design domain D (Fig. A.16(a)) is stored on the lattice points expanding D (Fig. A.16(b)). In this study, we use the following procedure which is similar but not identical to marching cubes: for All voxels do if The level set function has the same sign on all 8 nodes then continue // no boundary is included in the voxel. // else for All edges in the voxel do if The level set function has the same sign on both ends then continue // no boundary element node is on the edge. // else By linearly interpolating the value of the level set function on the edge, we search an intersection of φ = 0 and the edge. On the intersection, we put a “node” on the edge (see Fig. A.16(c)). end if end for for All surface of the voxel do By connecting “nodes” and the lattice edge on which the level set function is positive, we create a polygon and divide the polygon into several triangles. end for By connecting “nodes”, we create a polygon and divide the polygon into several triangles. end if

520

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

end for for All triangle’ edges do if The edge is not sheared by two triangles then Delete the triangle. continue end if end for for All triangles do We calculate the normal vectors of the triangles. // required for STL format. // end for Note that, although the aspect ratio of the obtained triangles is not necessarily small (Fig. A.16(d)), the quality of the mesh can be improved by NETGEN. References [1] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (2) (1988) 197–224. [2] Katsuyuki Suzuki, Noboru Kikuchi, A homogenization method for shape and topology optimization, Comput. Methods Appl. Mech. Engrg. 93 (3) (1991) 291–318. [3] T. Yamada, K. Izui, S. Nishiwaki, A. Takezawa, A topology optimization method based on the level set method incorporating a fictitious interface energy, Comput. Methods Appl. Mech. Engrg. 199 (45) (2010) 2876–2891. [4] T. Yamada, K. Izui, S. Nishiwaki, A level set-based topology optimization method for maximizing thermal diffusivity in problems including design-dependent effects, J. Mech. Des. 133 (3) (2011) 031011. [5] G. Jing, H. Isakari, T. Matsumoto, T. Yamada, T. Takahashi, Topological sensitivity of the objective function defined on morphing boundaries of two-dimensional heat conduction problems, Bound. Elem. Other Mesh Reduct. Methods XXXVII 57 (2014) 3. [6] Baotong Li, Jun Hong, Xiangyang Tian, Generating optimal topologies for heat conduction by heat flow paths identification, Int. Commun. Heat Mass Transfer 75 (2016) 177–182. [7] K. Yaji, T. Yamada, M. Yoshino, T. Matsumoto, K. Izui, S. Nishiwaki, Topology optimization using the lattice Boltzmann method incorporating level set boundary expressions, J. Comput. Phys. 274 (2014) 158–181. [8] E.M. Papoutsis-Kiachagias, K.C. Giannakoglou, Continuous adjoint methods for turbulent flows, applied to shape and topology optimization: Industrial applications, Arch. Comput. Methods Eng. (2015) 1–45. [9] Kentaro Yaji, Takayuki Yamada, Masato Yoshino, Toshiro Matsumoto, Kazuhiro Izui, Shinji Nishiwaki, Topology optimization in thermalfluid flow using the lattice Boltzmann method, J. Comput. Phys. 307 (2016) 355–377. [10] J.S. Jensen, Topology optimization problems for reflection and dissipation of elastic waves, J. Sound Vib. 301 (1) (2007) 319–340. [11] O. Sigmund, J.S. Jensen, Systematic design of phononic band–gap materials and structures by topology optimization, Phil. Trans. R. Soc. A 361 (1806) (2003) 1001–1019. [12] Z.D. Ma, N. Kikuchi, H.C. Cheng, Topological design for vibrating structures, Comput. Methods Appl. Mech. Engrg. 121 (1) (1995) 259–280. [13] Y. Noguchi, T. Yamada, M. Otomori, K. Izui, S. Nishiwaki, An acoustic metasurface design for wave motion conversion of longitudinal waves to transverse waves using topology optimization, Appl. Phys. Lett. 107 (22) (2015) 221909. [14] Jeonghoon Yoo, Noboru Kikuchi, Topology optimization in magnetic fields using the homogenization design method, 2000. [15] F. Abe, H. Isakari, T. Takahashi, T. Matsumoto, A topology optimisation in two-dimensional electromagnetics with the level set method and the boundary element method, Trans. JASCOME 13 (2013) 37–42 (in Japanese). [16] Y. Kourogi, H. Isakari, T. Takahashi, T. Yamada, T. Matsumoto, On a topological sensitivity analysis of three-dimensional electromagnetic wave problems with the boundary element method and its application to a level set based structural optimization, Trans. JASCOME 13 (2013) 55–60 (in Japanese). [17] H. Isakari, K. Nakamoto, T. Kitabayashi, T. Takahashi, T. Matsumoto, A multi-objective topology optimization for 2D electro-magnetic wave problems with the level set method and BEM, Eur. J. Comput. Mech. 25 (1–2) (2016) 165–193. [18] Namjoon Heo, Jeonghoon Yoo, Dielectric structure design for microwave cloaking considering material properties, J. Appl. Phys. 119 (1) (2016) 014102. [19] Eddie Wadbro, Martin Berggren, Topology optimization of an acoustic horn, Comput. Methods Appl. Mech. Engrg. 196 (1) (2006) 420–436. [20] Lirong Lu, Takashi Yamamoto, Masaki Otomori, Takayuki Yamada, Kazuhiro Izui, Shinji Nishiwaki, Topology optimization of an acoustic metamaterial with negative bulk modulus using local resonance, Finite Elem. Anal. Des. 72 (2013) 1–12. [21] Joong Seok Lee, Eun Il Kim, Yoon Young Kim, Jung Soo Kim, Yeon June Kang, Optimal poroelastic layer sequencing for sound transmission loss maximization by topology optimization method, J. Acoust. Soc. Am. 122 (4) (2007) 2097–2106. [22] T. Yamamoto, S. Maruyama, S. Nishiwaki, M. Yoshimura, Topology design of multi-material soundproof structures including poroelastic media to minimize sound pressure levels, Comput. Methods Appl. Mech. Engrg. 198 (17) (2009) 1439–1455. [23] Joong Seok Lee, Yeon June Kang, Yoon Young Kim, Unified multiphase modeling for evolving, acoustically coupled systems consisting of acoustic, elastic, poroelastic media and septa, J. Sound Vib. 331 (25) (2012) 5518–5536. [24] Joong Seok Lee, Peter G¨oransson, Yoon Young Kim, Topology optimization for three-phase materials distribution in a dissipative expansion chamber by unified multiphase modeling approach, Comput. Methods Appl. Mech. Engrg. 287 (2015) 191–211.

H. Isakari et al. / Comput. Methods Appl. Mech. Engrg. 315 (2017) 501–521

521

[25] Yuki Noguchi, Takayuki Yamada, Takashi Yamamoto, Izui Kazuhiro, Shinji Nishiwaki, Topological derivative for an acoustic-elastic coupled system based on two-phase material model, Mech. Eng. Lett. 2 (2016) 16–00246. [26] K. Abe, S. Kazama, K. Koro, A boundary element approach for topology optimization problem using the level set method, Commun. Numer. Methods Eng. 23 (5) (2007) 405–416. [27] Jianbin Du, Niels Olhoff, Minimization of sound radiation from vibrating bi-material structures using topology optimization, Struct. Multidiscip. Optim. 33 (4–5) (2007) 305–321. [28] V. Rokhlin, Rapid solution of intergral equations of classical potential theory, J. Comput. Phys. 60 (1985) 187–207. [29] L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (2) (1987) 325–348. [30] M. Bebendorf, Hierarchical Matrices, Springer, 2008. [31] PG. Martinsson, V. Rokhlin, A fast direct solver for boundary integral equations in two dimensions, J. Comput. Phys. 205 (1) (2005) 1–23. [32] H. Isakari, K. Kuriyama, S. Harada, T. Yamada, T. Takahashi, T. Matsumoto, A topology optimisation for three-dimensional acoustics with the level set method and the fast multipole boundary element method, Mech. Eng. J. 1 (4) (2014) CM0039. [33] M. Bonnet, N. Nemitz, FM-BEM and topological derivative applied to acoustic inverse scattering, in: Boundary Element Analysis, Springer, 2007, pp. 187–212. [34] A.A. Novotny, R.A. Feij´oo, E. Taroco, C. Padra, Topological sensitivity analysis, Comput. Methods Appl. Mech. Engrg. 192 (7) (2003) 803–829. [35] A. Carpio, M.L. Rapn, Solving inhomogeneous inverse problems by topological derivative methods, Inverse Problems 24 (4) (2008) 045014. [36] T. Kondo, H. Isakari, T. Takahashi, T. Matsumoto, A topology optimisation in 3d-acoustics for scatterers with impedance boundaries with the level set method and the fast multipole boundary element method, Trans. JASCOME 14 (2014) 19–24 (in Japanese). [37] Hiroshi Isakari, Periodic fmms and Calderon’s preconditioning in acoustics and elastodynamics, 2012. [38] Arno Kimeswenger, Olaf Steinbach, Gerhard Unger, Coupled finite and boundary element methods for vibro-acoustic interface problems, in: Domain Decomposition Methods in Science and Engineering XXI, Springer, 2014, pp. 507–515. [39] Gordon C. Everstine, Francis M. Henderson, Coupled finite element/boundary element approach for fluid–structure interaction, J. Acoust. Soc. Am. 87 (5) (1990) 1938–1947. [40] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 168–178. [41] Maurice A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 179–191. [42] Matthias Fischer, Lothar Gaul, Fast bem–fem mortar coupling for acoustic–structure interaction, Internat. J. Numer. Methods Engrg. 62 (12) (2005) 1677–1690. [43] C.J. Luke, P.A. Martin, Fluid-solid interaction: acoustic scattering by a smooth elastic obstacle, SIAM J. Appl. Math. 55 (4) (1995) 904–922. [44] A. Cemal Eringen, Erdogan S Suhubi, C.C. Chao, Elastodynamics, Vol. II, Linear theory, J. Appl. Mech. 45 (1978) 229. [45] Bojan B. Guzina, Ivan Chikichev, From imaging to material identification: a generalized concept of topological sensitivity, J. Mech. Phys. Solids 55 (2) (2007) 245–279. [46] Snorre H. Christiansen, Jean-Claude N´ed´elec, A preconditioner for the electric field integral equation based on Calderon formulas, SIAM J. Numer. Anal. 40 (3) (2002) 1100–1135. [47] Kazuki Niino, Naoshi Nishimura, Preconditioning based on Calderon’s formulae for periodic fast multipole methods for Helmholtz’ equation, J. Comput. Phys. 231 (1) (2012) 66–81. [48] Hiroshi Isakari, Kazuki Niino, Hitoshi Yoshikawa, Naoshi Nishimura, Calderon’s preconditioning for periodic fast multipole method for elastodynamics in 3D, Internat. J. Numer. Methods Engrg. 90 (4) (2012) 484–505. [49] G. Jing, H. Isakari, T. Matsumoto, T. Yamada, T. Takahashi, Level set-based topology optimization for 2D heat conduction problems using bem with objective function defined on design-dependent boundary with heat transfer boundary condition, Eng. Anal. Bound. Elem. 61 (2015) 61–70. [50] Samuel Amstutz, Analysis of a level set method for topology optimization, Optim. Methods Softw. 26 (4–5) (2011) 555–573. [51] NETGEN - automatic mesh generator. http://www.hpfem.jku.at/netgen/.