Optimum design of phononic crystal perforated plate structures for widest bandgap of fundamental guided wave modes and maximized in-plane stiffness

Journal of the Mechanics and Physics of Solids 89 (2016) 31–58

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

Optimum design of phononic crystal perforated plate structures for widest bandgap of fundamental guided wave modes and maximized in-plane stiffness Saeid Hedayatrasa a,n, Kazem Abhary a, Mohammad Uddin a, Ching-Tai Ng b a b

School of Engineering, University of South Australia, Mawson Lakes, SA 5095, Australia School of Civil, Environmental and Mining Engineering, University of Adelaide, SA 5000, Australia

a r t i c l e i n f o

abstract

Article history: Received 10 November 2014 Received in revised form 7 November 2015 Accepted 19 January 2016 Available online 22 January 2016

This paper presents a topology optimization of single material phononic crystal plate (PhP) to be produced by perforation of a uniform background plate. The primary objective of this optimization study is to explore widest exclusive bandgaps of fundamental (first order) symmetric or asymmetric guided wave modes as well as widest complete bandgap of mixed wave modes (symmetric and asymmetric). However, in the case of single material porous phononic crystals the bandgap width essentially depends on the resultant structural integration introduced by achieved unitcell topology. Thinner connections of scattering segments (i.e. lower effective stiffness) generally lead to (i) wider bandgap due to enhanced interfacial reflections, and (ii) lower bandgap frequency range due to lower wave speed. In other words higher relative bandgap width (RBW) is produced by topology with lower effective stiffness. Hence in order to study the bandgap efficiency of PhP unitcell with respect to its structural worthiness, the inplane stiffness is incorporated in optimization algorithm as an opposing objective to be maximized. Thick and relatively thin Polysilicon PhP unitcells with square symmetry are studied. Non-dominated sorting genetic algorithm NSGA-II is employed for this multi-objective optimization problem and modal band analysis of individual topologies is performed through finite element method. Specialized topology initiation, evaluation and filtering are applied to achieve refined feasible topologies without penalizing the randomness of genetic algorithm (GA) and diversity of search space. Selected Pareto topologies are presented and gradient of RBW and elastic properties in between the two Pareto front extremes are investigated. Chosen intermediate Pareto topology, even not extreme topology with widest bandgap, show superior bandgap efficiency compared with the results reported in other works on widest bandgap topology of asymmetric guided waves, available in the literature. Finally, steady state and transient frequency response of finite thin PhP structures of selected Pareto topologies are studied and validity of obtained bandgaps is confirmed. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Phononic crystal Plate Topology optimization Guided wave Stiffness

1. Introduction Phononic crystals (PhCr) are acoustic meta-materials with promising manipulation capabilities on propagation of sound and elastodynamic waves (Deymier, 2011). PhCrs are indeed heterogeneous materials produced by periodic modulation of n

Corresponding author. E-mail address: [email protected] (S. Hedayatrasa).

http://dx.doi.org/10.1016/j.jmps.2016.01.010 0022-5096/& 2016 Elsevier Ltd. All rights reserved.

32

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

acoustic impedance in a lattice structure through either integration of two or more contrasting materials, or making void inclusions in a single material. The main feature of PhCrs making them known as acoustic bandgap materials is the existence of frequency ranges over which propagation of vibroacoustic waves is prohibited. This phenomenon is caused by constructive reflection and superposition of waves at the interface of periodic heterogeneities i.e. Bragg and Mie resonant scatterings (Olsson III and Elkady (2009)). Therefore the efficiency of PhCr is wave length dependent and the frequency range of bandgap is limited by its lattice periodicity constant or in other words its feature size. However, due to the strong wave length scale heterogeneity of PhCrs, they possess other astonishing wave manipulation properties. Introducing any kind of defect in phononic lattice e.g. by altering the features of a few adjacent cells, leads to advent of local resonance modes within phononic bandgap frequency. This capability is used to trap and guide waves inside defects specially tuned for desired frequencies (Zhu and Semperlotti, 2013, Olsson III and El-kady (2009)). Moreover, studying the equal frequency contours of PhCrs show flat and concave wave fronts below bandgap frequency applicable for self-collimation and steering of waves (Wang et al., 2013, Lin et al., 2009). Consequently, it is worthwhile to adjust the phononic properties so that the width and location of bandgap is optimized for application of interest. The efficiency of PhCr is principally defined by the features of its irreducible representative element (Unitcell) i.e. its geometry, topology and composition. The focus of this paper is to investigate the optimum topology of PhCr plate (PhP) unitcell with 2D periodicity for maximized bandgap width of guided waves. Guided waves are structural waves confined by traction free surfaces of thin walled structures, so called plate waves when guided by parallel faces of plates. In-plane symmetric (S) and anti-plane asymmetric (A) Lamb waves as well as symmetric shear horizontal (SHS) and asymmetric shear horizontal (SHA) in-plane waves are the well-known guided wave modes generated in a plate structure. Special characteristics of guided waves confined by finite thickness of such structures, make them ideal for non-destructive evaluation purposes (Su and Ye, 2009, Veidt et al., 2008) as well as production of low loss resonators, filters and waveguides (Mohammadi, 2010, Lin et al., 2014). The guided wave dispersion in PhP is governed by its transversal anisotropy i.e. plate’s thickness in addition to its planar anisotropy measured by lattice periodicity and unitcell constitution. The bandgaps of plate waves in 2D PhPs of prescribed topologies have been extensively studied (Pennec et al., 2010, Wu et al., 2011, Charles et al., 2006, Khelif et al., 2006). PhPs produced by either periodic insertion, attachment or perforation of cylindrical inclusions on a background plate have been generally considered in these studies. Phononic lattice with different patterns e.g. square, rectangular or polygonal have been considered. The phononic unitcell itself could contain a specific pattern of circular inclusions like the work by Li et al. (2009) who studied plate wave gaps in PhPs with Archimedean-like tilings. However, optimum topology of PhCrs and acoustic meta-materials have been widely investigated, e.g. (Olhoff et al., 2012, Hussein et al., 2006, Rupp et al., 2007, Manktelow et al., 2013, Krushynska et al., 2014, Hedayatrasa et al., 2015). Essentially topology with maximum relative bandgap width (RBW) between subsequent modes of interest is desired (Sigmund and Jensen, 2003). RBW is the ratio of bandgap width over midgap frequency, and so seeking maximum RBW results in a topology which can provide widest bandgap at lowest frequency range for specific unitcell size. Consequently, the relevant topology supports phononic wave manipulation over widest frequency range through miniature unitcells compared to wavelength. Most of topology optimization studies in relation to 2D periodic PhCrs have been concerned with bandgaps of in-plane and/or anti-plane bulk waves while only few works have been devoted to bandgaps of guided waves in 2D PhPs. Halkjær et al. (2006) studied the optimum topology of porous Polycarbonate PhP with rhombic unitcell for maximized RBW of first couple of flexural (asymmetric) plate waves. The Mindlin plate theory was implemented for definition of band structure of bending waves and gradient based optimization was performed through moving asymptotes method. Since the best topology for maximized RBW did not have acceptable stiffness, new objective was introduced to conversely search widest bandgap width at higher frequency ranges. Finally the discontinuities of optimized topology were locally modified for satisfactory stiffness and manufacturability. In another investigation by Bilal and Hussein (2012) the optimum topology of thin porous silicon PhPs for maximized RBW of basic flexural waves was studied. The Mindlin plate's theory was implemented and topology of square unitcell was optimized through genetic algorithm (GA) as an evolutionary based method. Single material porous unitcell with uniform through thickness constitution have been commonly considered, in which, the achieved topology could be simply produced by perforation of a uniform background plate. The high acoustic impedance contrast produced by the porosities, usually filled with air, leads to relatively wide acoustic bandgaps. Bandgaps of such porous materials are governed by wave reflection and scattering at the interface of inhomogeneities produced by perforation profile. Therefore the search for highest RBW naturally results in topologies with nearly isolated domains and in other words thin connectivity. The finer the topology resolution the thinner connectivity in the optimized topology for maximized bandgap. Moreover, porous topology with lower effective stiffness can produce bandgap at lower frequency range and so higher RBW due to its lower wave speed. Thus largest achievable RBW is extremely dependent on assumed unitcell’s resolution and relevant topology generally has low structural worthiness i.e. stiffness. Nevertheless, none of earlier works on topology optimization of porous PhCrs (Halkjær et al., 2006, Bilal and Hussein, 2012, Dong et al., 2014a) took into consideration the structural worthiness of achieved topologies. Although the mesh dependency of topology can be controlled by constraining the minimum length scale of the features, this technique cannot assure the optimality of stiffness. Furthermore, if the topology resolution is not fine enough, this approach may degrade the optimality of bandgap due to imposed feature size constraint to the entire topology area. In earlier studies the bandgaps of asymmetric guided wave modes have been solely explored (Halkjær et al., 2006, Bilal

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

33

and Hussein, 2012, Dong et al., 2014a) while it is also of great value to investigate the topology of PhPs for maximized RBW of symmetric modes as well as complete gap of mixed guided wave modes. Complete bandgap of guided waves can e.g. be used as vibration-less supporting structure of delicate dynamic device like gyroscope for immunity to noise sources. Guided waves can be also manipulated through exceptional flat and concave equal frequency contours produced near the bandgap frequency of corresponding mode, for self-collimation and wave steering purposes. If the bandgap is exclusive to either symmetric or asymmetric modes, then the modes can be handled individually. Allowed guided wave modes can be almost purely launched as the banned modes are filtered within the corresponding exclusive bandgaps. Bandgap mode on the other hand can be trapped or guided within appropriate defects in PhP lattice. So, detective wave beams of symmetric or asymmetric modes can be individually excited, steered and captured through PhP elements with promising application in structural health monitoring. Other potential applications are in acoustic energy harvesting and radio frequency (RF) communication (Olsson III and El-kady (2009)). Therefore, in this study the optimum topology of 2D porous PhP with square symmetry is explored for maximized RBW of fundamental (first order) guided wave modes for three different objectives: (i) exclusive bandgap of fundamental asymmetric wave modes A0 and SHA0, (ii) exclusive bandgap of fundamental symmetric wave modes S0 and SHS0, and (iii) complete bandgap of mixed wave modes. Besides, in-plane stiffness of the porous unitcell is incorporated as the second objective to be maximized in order to investigate the gradient of bandgap width versus stiffness and to obtain acoustic bandgaps through optimized structural stiffness. Multi-objective GA is employed to search for the optimum topology taking into account these two objectives. The modal band analysis and numerical homogenization of unitcell is performed through finite element modeling (FEM) for fitness calculation purpose. Specialized filtering is applied to the topologies in order to incline the search space towards feasible solutions with no discontinuities without compromising its diversity and randomness. Since the bandgap efficiency of PhP is substantially dependent on its transversal aspect ratio, thick and relatively thin unitcells are studied separately. The governing equations, objective functions and developed FEM-GA framework are firstly described in Section 2. Detailed explanation of FEM modeling of PhP unitcell as well as adopted GA algorithm and related challenges are discussed in Section 3. Then achieved optimization results are presented and discussed in Section 4. Selected Pareto topologies are given and gradient of RBW and elastic properties are investigated between the two Pareto front extremes (one with widest bandgap and the other one with highest stiffness). Finally, steady state and transient frequency response of finite thin PhP structures of selected Pareto topologies are studied in Section 5 and validity of obtained results is confirmed.

2. Analytical approach to the problem 2.1. Modal band analysis of unitcell 2D PhP is considered to be produced by two dimensional integration of the square basic unitcell of length a and transversal aspect ratio r = a/h, as shown schematically in Fig. 1 for an arbitrary perforation profile. The unitcell is homogeneous along thickness h and heterogeneous in x–y plane due to existence of transversally perforated region(s). The planar anisotropy can lead to constructive scattering of guided waves and consequently bandgaps of guided waves. The equation of motion in such heterogeneous plate structure with linear isotropic material properties can be expressed as:

¨ ∇(λ + μ)∇. U + ∇. μ∇U = ρU

(1)

where U = {u, v, w} is the displacement vector, λ and μ are isotropic elastic Lame constants, ρ is density and ∇ is the gradient operator. The Bloch-Floquet solution (Kittel, 1986) for 2D PhP with infinite periodicity can be developed by harmonic modulation of a periodic displacement filed as follows:

U(X, t ) = Up(X)ei(kx − ωt )

(2)

where X = {x, y, z} is the location vector, UP is a periodic function in x–y plane with periodic constant a according to the lattice periodicity, ω = 2πf is circular frequency,k = {k x , k y} is the wave vector and i = −1 . So after applying Bloch-Floquet periodic boundary condition on the plate unitcell with constraint free upper and lower surfaces, we will have:

Fig. 1. Illustration of a 2D PhP with arbitrary perforation profile, lattice constant a along x and y axes and thickness h along z axis.

34

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

UX + = UX −eik xa

(3)

Uy + = Uy −eik ya

(4)

where subscripts x þ and x  stand for corresponding points at two unitcell’s faces parallel to y–z plane at x = 0 and x = a , respectively. Similarly subscripts yþ and y stand for corresponding points at two unitcell's faces parallel to z–x plane at y = 0 and y = a . In this study, 3D FEM model of the unitcell is developed in ANSYS APDL (ANSYSs Academic Research, Release 14.5) for calculation of plate waves' modal band structure. Since the Bloch-Floquet boundary conditions (Eq. (3) and Eq. (4)) are complex, two superimposed meshes are modeled so that one represents the real term and the other one imaginary term of complex displacement domain (Aberg and Gudmundson, 1997). Consequently the relevant constraint equations connecting two coinciding meshes are expanded as follows for real (Re) and imaginary (Im) terms of Eq. (3), respectively:

(UX +)Re = (UX −)Im cos(k xa) − (UX −)Re sin(k xa)

(5)

(UX +)Im = (UX −)Im cos(k xa) + (UX −)Re sin(k xa)

(6)

and likewise for Eq. (4). So, by applying periodic boundary conditions on vertical boundaries and considering natural traction free boundary conditions for horizontal surfaces, the modal frequencies of the model are calculated versus wave vector defining the band structure of relevant plate waves. Theoretically wave vector k can take any value, but due to periodicity in Bloch-Floquet condition it can be limited to the first Brillouin zone k ∈ [ − π /a, π /a] and further to irreducible Brillouin zone based on the symmetry of topology (Fig. 2). According to the common practice, a discrete vector of wave numbers tracing only the border of irreducible Brillouin zone, shown by triangle ΓXM in Fig. 2, has proved to be adequate. In order to study the band structure of low order Symmetric (S0 þSHS0) and Asymmetric (A0 þ SHA0) plate waves individually, specific boundary conditions have been considered to stimulate the relevant modes exclusively. This approach relies on the nature of displacement filed through the thickness of plate. With regard to harmonic definition of guided waves in uniform plates (Kundu, 2004), the harmonic oscillations of a heterogeneous periodic plate in any planar direction in x–y plane, e.g. X axis, can be explained by symmetric and asymmetric terms with respect to the midplane as follows. Lamb wave oscillations in z–x plane are formulated as:

u = {ikP1(x)cos αz + βQ 1(x)cos βz}eikx + {kR1(x)sin αz + iβS1(x) sin βz}eikx

(7)

w = {iαP2(x)sin αz + kQ 2(x)sin βz}eikx + {αR2(x)cos αz + ikS2(x)cos βz}eikx

(8)

based on Stokes–Helmholtz decomposition, and anti-plane oscillations along y axis are defined by:

v = {P3(x)cos γz}eikx + {iQ 3(x)sin γz}eikx

(9)

where functionsPi, Q i, Ri andSi are a periodic amplitudes of relevant oscillations. Transversal wave numbers α , β and γ are indeed vertical components of in-plane normal, in-plane shear and anti-plane shear elastic wave numbers, respectively. The first and second terms of Eqs. (7–9) stand for symmetric and asymmetric plate modes. So, with regard to traction free surfaces of the plate and explained symmetry of displacement filed, these plate modes can be decoupled by applying appropriate boundary conditions to the mid-plane as follows (Karmazin, 2013):

Fig. 2. Wave number space and Irreducible Brillouin zone for unitcell with square symmetry.

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

{w , τxz , τyz} {u, v, σz}

z= 0

z= 0

= 0 ( Symmetricmode)

= 0 ( Asymmetricmode)

35

(10) (11)

where σz is the transversal normal stress and τxz , τyz are transversal shear stresses. Whereas, the displacement filed is directly or inversely symmetric with respect to the mid-plane for different polarizations, just half of unitcell from midplane z = 0 to top plane z = h/2 is modeled for mode decoupling purpose. So, on the plane of symmetry (z = 0), the constraint {u, v} = 0 is applied to enforce asymmetric modes, or constraint w = 0 is applied to enforce symmetric modes. Zero stress boundary conditions required as per Eqs. (10) and (11) are naturally satisfied by unconstrained degrees of freedom on the mid-plane in both cases.

2.2. Objective functions Maximum RBW and highest in-plane stiffness of PhP unitcell are two opposing objectives of this topology optimization study. Whereas bandgaps are governed by interfacial reflections of inhomogeneities produced by perforation profile, the search for highest RBW naturally results in topologies with nearly isolated domains and in other words thin connectivity. Finer topology resolution generally results in higher RBW by thinner connectivity, i.e. less stiffness. That is why the in-plane stiffness of the PhP unitcell is introduced as a second objective to study the variation of bandgap efficiency with respect to stiffness. As for first objective regarding RBW of PhP, it is fundamentally desired to manipulate the largest wave length possible through specific unitcell size. Hence the topology with largest bandgap width at lowest frequency range is sought during optimization, and the first objective function F1 is defined as:

F1 =

minin=k 1ωj2+ 1(ki ) − maxin=k 1ωj2(ki ) 0.5(minin=k 1ωj2+ 1(ki ) + maxin=k 1ωj2(ki ))

(12)

where F1 is the gap width over mean value of gap between Eigen values ω2 of two subsequent modal branches jth and (jþ 1)th of interest over the nk discrete search points k i on the border of irreducible Brillouin zone (Fig. 2). The second objective of optimization is to maximize the in-plane stiffness of PhP unitcell. Principally, the strain energy of the unitcell under specified loading condition represents its relevant stiffness. Higher strain energy implies lower stiffness or higher compliance. The strain energy W stored per volume V of a finite element subjected to in-plane stress in x–y plane is theoretically defined as (Kaw, 2005):

Wv =

dW 1 = (σxεx + σyεy + τxyγxy) dV 2

(13)

where {σx, σy, τxy} and {εx , εy , γxy} are orderly corresponding in-plane stress and strain tensors. Now by assumption of stress state σx = σy =

τxy b

= σ and with regard to linear elastic stress-strain relations, the Eq. (13) turns to:

⎛1 − ν b2 ⎞ e ⎟ Wv = σ 2⎜ + E 2 Ge ⎠ ⎝ e

(14)

where Ee , Ge and υe are the equivalent in-plane orthotropic elastic modulus, shear modulus and Poisson’s ratio of assumed square symmetric PhP unitcell, respectively, and b is stress state ratio. Finally, the second objective function F2 is defined as strain energy compliance of unitcell to be minimized:

⎛1 − ν ⎞ ⎛ 1⎞ e F2 = ⎜ ⎟ + 0.5b2⎜ ⎟ ⎝ Ee ⎠ ⎝ Ge ⎠

(15)

in order to attain maximum in-plane stiffness. The ratio of longitudinal compliance to shear compliance and their priority could be indeed controlled by the stress state ratio b . High order multiscale homogenization is commonly used in order to define the equivalent material properties of periodic lattice structures taking into account the unitcell scale variations (Gonella and Ruzzene, 2008, Oliveira et al., 2009). However the overall macromechanical stiffness of the PhP unitcell is desired to be maximized in this study. Hence, the classical first order homogenization is employed to define the equivalent elastic properties of PhP unitcell and calculate F2. For this purpose, periodic boundary condition is applied and overall test strain ε¯ is prescribed for the unitcell. Then the equivalent orthotropic stiffness De is derived based on average stress filed σ¯ of elements over the entire unitcell:

σ¯ =

1 V

∫V σdV

(16)

36

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

⎡C C C 0 0 0 ⎤ ⎥ ⎢ 11 12 13 ⎢ C12 C11 C13 0 0 0 ⎥ ⎢C C C 0 0 0 ⎥⎥ 13 13 33 De = σε ¯ ¯ −1 = ⎢ ⎢ 0 0 0 C44 0 0 ⎥ ⎥ ⎢ ⎢ 0 0 0 0 C55 0 ⎥ ⎢⎣ 0 0 0 0 0 C ⎥⎦ 55

(17)

Three different strain fields with pure longitudinal strain εx = 0.001, transversal strain εz = 0.001 and shear strain τxy = 0.001 are applied individually to define homogenized stiffness constants C11, C12, C13, C33&C44 of square symmetric topology and hence to calculate the required in-plane engineering elastic propertiesEe , Ge and υe . Detailed information regarding numerical homogenization of composite materials using FEM and periodic boundary condition can be found in Bendsoe and Sigmund (2003). The dispersion-less (low frequency) propagation speed of S0 and SH0 plate modes determined from initial slope of relevant dispersion curves of unitcell are governed by in-plane longitudinal and shear stiffness of unitcell (Kundu, 2004). One may alternatively consider them as representatives of unitcell's effective stiffness in definition of objective function. However, it has complications in orthotropic medium and in this work the numerical homogenization is employed for this purpose. 2.3. Multi-objective optimization strategy Evolutionary multiobjective genetic algorithm (GA) optimization is employed for this purpose based on non-dominated sorting GA (NSGA-II) (Pratap et al., 2002). In topology optimization of phononic bandgaps, it is well known that the highest contrast at the interface of heterogeneities and actually no intermediate area is desired for maximized interfacial reflections i.e. Bragg scattering. GA uses a coded string of topology to be mapped into the discretized design domain, searches over a population of potential solutions, reproduces better generations using crossover and mutation functions and iteratively produces better and better approximations to the solution (Weise, 2006, Hajela et al., 1993, Chapman et al., 1994). So it is a supreme tool for intelligent design of PhCrs with favored discrete design space free of intermediate pseudo domains. The discrete nature of GA also can easily handle definite topology constrains through image processing. Moreover, GA searches for optimum topologies within a population of designs in a single run and performs stochastic search over the entire search space and so unlikely to get stuck in a local optimum. Earlier works on topology optimization of PhCrs through GA have proved its competency in this area (Gazonas et al., 2006, Dong et al., 2014b, Dong et al., 2014c, Hussein et al., 2006, 2007, Bilal and Hussein, 2012, Liu et al., 2014, Manktelow et al., 2013). Multi-objective GA was implemented by the authors to study the optimum topology of 1D PhPs with respect to filling fraction of constituents and gradient 1D PhP structures were introduced (Hedayatrasa et al., 2015). Moreover, Dong et al. (2014a) used multi-objective GA for design of 2D PhCrs with maximum RBW of bulk waves and minimum mass. For more information regarding implementation of GA for structural topology optimization see e.g. (Wang and Tai, 2005, Denies et al., 2012, Gazonas et al., 2006, Kunakote and Bureerat, 2011). NSGA-II sorts GA population based on non-dominated fitness and crowding distance. In this way the designs with the most dominated fitness are placed at the first rank (Pareto front) and among them those with highest crowding distance are preferred. An individual dominates another one provided it has better fitness for both objectives. If exceptionally both individuals have the same fitness for just one objective, the one which has higher fitness for other objective is dominant. The crowding distance on the other hand is a measure to enforce even distribution and diversity of optimized solutions within Pareto front (Pratap et al., 2002).

3. Modeling parameters and adopted NSGA-II optimization algorithm Polysilicon with material properties Es = 169 GPa , υs = 0.22 and ρs = 2330 kg/m3 (Sharpe et al., 1997) is taken as the solid background material of PhP in this work which has been widely used as structural material for fabrication of micro-devices like micro PhCrs (Olsson III and El-kady (2009)). Unitcell with element resolution 32 × 32 in x–y plane is selected leading to a design domain with 136 independent design variables for square symmetric topology. Two transversal aspect ratios considered in this study are r = a/h = 2 and r = 10 in order to investigate the design and efficiency of thick and relatively thin plate unitcells, respectively. ANSYS linear solid element SOLID185 is used for FEM modeling of thick unitcells, and solidshell element SOLSH190 is used for modeling thin unitcells. The selected planar mesh resolution 32 × 32 provides adequate accuracy for FEM analysis of fundamental modal frequencies of interest, provided appropriate number of element layers is considered through the thickness of unitcell. Since GA optimization demands FEM analysis to be performed for a large number of iterations, it is important to develop a reliable model with minimized computational intensity. Proper mesh resolution along the thickness of unitcell h relies on its transversal aspect ratio r = a/h and the frequency range of highest plate wave mode to be calculated. SOLSH190 element usually gives more accurate predictions compared to classical shell elements (with Mindlin’s Plate formulation) when the

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

37

shell is thick (ANSYS 14.5 element reference). Thus, for modeling moderately thick PhP unitcells of aspect ratio r = 10 a single layer of solid-shell elements provides good accuracy for modal band analysis of low order modal branches. Two element layers are considered only for the unitcell of aspect ratio r = 2 when maximizing RBW of complete plate wave's gap. Of course, for calculation of the modal band structure and RBW of achieved topologies higher mesh resolutions of 10 and 4 layers are considered for aspect ratios 2 and 10, respectively. More details and relevant mesh refinement study are given in the Appendix A. For FEM evaluation of topologies, only the solid part of topology has to be modeled and appropriate boundary conditions should be applied to existing boundaries at any evaluation iteration. This approach is challenging to be included in an automated algorithm like GA and of course its complexity depends on the capabilities of implemented FEM solver. Also, as a general simplified approach a compliant material can be assigned to the void elements by a few orders less than the solid background. However, the degradation order of solid background properties to be considered as void and the constitution of randomly generated topologies are crucial factors in reliable calculation of RBW. The Appendix B explains the matter and scrutinizes associated complicacies. In the present study the density and elastic modulus of solid background are reduced by 9 orders of magnitude to be considered as void material. The following adopted NSGA-II algorithm is coded in MATLAB (MATLAB 8.0, The MathWorks Inc.) which iteratively executes the developed ANSYS macro file for FEM modeling and fitness evaluation purpose. Topology with square symmetry is prescribed for assumed square shape of PhP unitcell and a binary bit-string design space is considered for defining the material type of any pixel (column of elements along z axis). Specialized topology initializing, filtering and modifications are applied to eliminate interrupting modes and make optimized topologies converge towards feasible connected ones without disturbing the randomness and diversity of search domain. The adopted NSGA-II algorithm covering these techniques is as follows: 1. Create a random bit-string initial population for the solid-void problem defining the material type at each column of elements in FEM model as a design variable (gene). An identical fully solid population is randomly mutated with specified probability 0.2 to deliver an initial population of likely connected topologies. 2. Pass individuals for fitness evaluation provided that the relevant topology includes valid segment touching boundary and with face to face connected elements. Else, degrade the topology by assigning very low fitness for both objectives. 3. Determine the objective functions for passed individuals by considering the valid boundary connected segment only. The segments having single node (hinge) connection to the boundary connected segment are retained for further evolution. 4. If the bandgap objective F1 ≤ 0 (i.e. no bandgap) then totally compliant material (void) is considered to define stiffness objective F2. In this way the Pareto front delivers a spread of bandgap topologies merely. 5. Sort the population based on the rank and crowding distance of individuals (Pratap et al., 2002). 6. Create parent pool through tournament selection. Two randomly selected individuals are compared and the one with higher rank is selected as candidate parent. If both have the same rank, the one with highest crowding distance is selected. 7. Create offspring population through single point crossover with certain probability. Two Childs are produced per couples of randomly selected parents. 8. Mutate gens of offspring population with specified probability 9. Modify topology for checkerboard elimination. Randomly selected pixel is inversed (i.e. 0 →1 and vice versa) if surrounded by more than 2 dissimilar adjacent face to face pixels out of 4. This well-known technique alleviates the checker board problem which leads to overestimated stiffness when linear elements are used. 10. Modify topology in order to eliminate invalid disconnected segments of topology. Randomly selected disconnected segment (including any segment with hinge connection) is eliminated with specified probability. 11. Produce intermediate population by accumulating offspring and current population and sort them based on ranking and crowding distance. 12. Select individuals for next generation first based on their ranking and second, if the first rank is bigger than required population based on their crowding distance. 13. Go to Step-2 if termination condition (prescribed number of generations) is not satisfied, otherwise stop Due to rarity of perfectly connected topologies showing bandgap properties defined during population initialization and reproduction stages, the aforementioned algorithm is a robust approach for filtering invalid individuals while maintaining the diversity of search domain. Among other alternatives, generally, (i) fitness penalization disturbs optimality of Pareto front, (ii) accepting perfectly connected topologies restricts the diversity of search space, and (iii) filtering individuals to acceptable ones interferes with the randomness of GA search. In fact, in introduced approach existence of any disconnected segment in the topology does not interrupt its fitness evaluation and let it contribute to the evolution process. Single node hinge connection is also avoided. Actually this weak form of connection provides segments with degraded or virtually stiff bonding which opens wide bandgaps and therefore relevant topologies usually stand as elite solutions in the first rank. MATLAB Image Processing Toolbox (MATLAB 8.0, The MathWorks Inc.) is employed to perform required topology filtering and morphology assessments. The mutation of offspring as mentioned in Step-7 of the algorithm manipulates the next generation with specified probability. This probability proved to remarkably affect the convergence rate and efficiency of GA search by retaining

38

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Table 1 Employed GA settings for topology optimization studies. Stage

1 2 3 4

Probability

Generations

Crossover

Mutation

Disconnected segment elimination

0.9 0.9 0.9 0.9

0.3 0.1 0.05 0.0

0.0 0.2 0.5 1.0

5 40 75 20 50

diversity and preventing premature convergence to a local optimum. This is due to infrequency of usable topologies opening bandgap particularly between lowest modes of interest. Moreover, the effective mutation probability depends on the filling fraction of high ranked topologies. Lower filling fraction means more void segments with less effectiveness probability on the fitness when inversed via mutation and so requiring higher mutation to exploit new connected designs. High mutation probability on the other hand should be avoided as it disturbs the genome and prevents convergence to a unique optimized topology. Population size 200 is chosen for optimization studies and GA operation probabilities are varied during optimization in 4 steps as given in Table 1. In Stage-1 the design space is coarsely searched by high mutation probability 0.3 for only 5 generations to possibly have a good resource of initiating population. Then lower mutation probability 0.1 is used for another 40 generations which greatly mitigates the risk of getting stuck in a local optimum. Afterwards, low mutation probability 0.05 is kept for up to 75 new generations. The segment elimination probability is added and increased gradually to avoid disturbance of genomes in the early stages of optimization. In fact, any disconnected segment potentially can fill a considerable fraction of unitcell and so initial intense elimination can make the GA search turn biased. Finally, the optimization is continued for a number of generations with 100% segment elimination in the absence of mutation to get refined feasible topologies with checkerboard improvements.

4. Topology optimization results The topology optimization results obtained through aforementioned multi-objective GA are delivered in this section. The optimum design of PhP unitcell is achieved for 3 different bandgap objectives, as follows: I. Bandgap Objective-1: to maximize RBW of lowest gap of fundamental asymmetric guided wave modes including asymmetric lamb mode A0 and asymmetric shear horizontal mode SHA0 II. Bandgap Objective-2: to maximize RBW of lowest gap of fundamental symmetric guided wave modes including symmetric Lamb mode S0 and symmetric shear horizontal mode SHS0 III. Bandgap Objective-3: to maximize RBW of lowest complete gap of mixed fundamental guided wave modes (S0, A0, SHS0 and SHA0) These are separately considered as the first objective of optimization besides considering maximized in-plane stiffness as the second objective. Two transversal aspect ratios r = 2 and r = 10 are considered to study the design and efficiency of thick and relatively thin plate unitcells, respectively. In any optimization case, a set of 10–12 solutions are selected from Pareto front containing a variety of alternative optimized topology modes with different RBWs and in-plane stiffness. The two extreme topologies of Pareto front with widest bandgap and highest stiffness are included; however stiffest topology with very minor bandgap is ignored. The selected optimized Pareto topologies for 8 different objective settings are referred to as TA to TH, to be discussed later on. Since RBW is defined as the gap width over midgap value of either Eigen values (F1) or Eigen frequencies (Δω/ω¯ ) of desired mode branches in different works (Halkjær et al., 2006, Bilal and Hussein, 2012), both quantities are obtained and included in the results. In order to evaluate the relative effective stiffness of optimized topologies with respect to stiffness of constitutive solid material, relative elastic modulus Er = Ee/Es and relative shear modulus Gr = Ge/Gs are calculated. For nondimensional modal analysis of unitcell, dimensionless circular frequency is defined as ωd = ωa/Cs where Cs = Es/ρs = 8516.6(m /s ) is the elastodynamic wave speed in Polysilicon. 4.1. Thin PhP unitcell with aspect ratio r = 10 4.1.1. Bandgap Objective-1 and effect of stress state ratio b As for the PhP of aspect ratio r = 10, firstly the results concerning Bandgap Objective-1 are presented to maximize RBW of lowest possible gap between the first couple of asymmetric modal branches. Three different stress state ratios b = 2 ,

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. 3. Pareto front of optimization for Bandgap Objective-1 and stress state ratio b = its two extremes.

39

2 , and selected optimized topologies TA-1 to TA-12 (Fig. 4(a)) across

0 and 2 are examined in this case. With regard to Eq. (15) in relation to the stiffness objective of optimization F2, the stress state ratio b = 2 imposes even contribution of normal and shear stiffness. Other values b = 0&2 obviously lead to pure contribution of normal stiffness and weighted dominant contribution of shear stiffness, respectively. Although the dominancy of normal or shear stiffness of final solution depends on the nature of bandgap topologies, the Pareto front is inclined to either side in this way. The Pareto front of optimization for Bandgap Objective-1 and b = 2 , and relative location of selected optimized topologies TA-1 to TA-12 are shown in Fig. 3. These selected Pareto front topologies are shown in Fig. 4(a). The gradient of RBW (either F2 and Δω/ω¯ ), filling fractionVf

(

)

(

)

Fig. 4. Topology optimization results of PhP unitcell with aspect ratio r = 10 for maximized Bandgap Objective-1, stress state ratio b = topologies TA-1 to TA-12 are selected from relevant Pareto front (Fig. 3).

2 ; Optimized

40

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. 5. Modal band structure of plate wave modes for (a) topology TA-1, and (b) Topology TA-12 (Fig. 4); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and lowest bandgap of asymmetric modes is highlighted in gray.

and homogenized Poisson’s ratio νe across these Pareto front topologies are given in Fig. 4b. The variation of log(Er ) and log(Gr ) are also given in Fig. 4(c) to define the magnitude order of relative elastic modulus of selected topologies. Moreover, the modal band structure of two extreme topologies of Pareto front with maximum RBW (TA-1) and maximum in-plane stiffness (TA-12) are delivered in Fig. 5(a) and (b), respectively, to show the relevant bandgaps. Symmetric and asymmetric modes are displayed by solid and dash lines, respectively, and desired lowest bandgaps are highlighted in gray.

(

)

(

)

Fig. 6. Topology optimization results of PhP unitcell with aspect ratio r = 10 for maximized Bandgap Objective-1, stress state ratio b = 2; Optimized topologies TB-1 to TB-12 are selected from relevant Pareto front.

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

41

Fig. 7. Modal band structure of plate wave modes for (a) topology TB-1, and (b) Topology TB-12 (Fig. 6); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and lowest bandgap of asymmetric modes is highlighted in gray.

(

)

(

)

Fig. 8. Topology optimization results of PhP unitcell with aspect ratio r = 10 for maximized Bandgap Objective-1, stress state ratio b = 0 ; Optimized topologies TC-1 to TC-12 are selected from relevant Pareto front.

Likewise, the results of other stress state ratios b = 0&2 are delivered in following Figures Figs. 6–9. In any loading case a broad range of topologies is obtained with various bandgap and stiffness efficiencies. The stiffness of topologies is increased by around 2 orders of magnitude across the Pareto front with relatively slight rise of filling fraction specially for b = 2. It is easily concluded from the results that increasing the stiffness of PhP unitcell decreases its RBW. An almost uniform reduction of RBW is evident for b = 2 by simultaneous increase of both elastic and shear modulus, as expected. But for b = 2 (Fig. 6) the shear modulus rises dramatically compared to minor increase of elastic modulus. Similarly, for the loading case

42

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. 9. Modal band structure of plate wave modes for (a) topology TC-1, and (b) Topology TC-12 (Fig. 8); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and lowest bandgap of asymmetric modes is highlighted in gray.

b = 0 (Fig. 8) the elastic modulus is increased much more than the shear one. Although applying the stress state ratios b = 0&2 do not lead to remarkable rise in extreme values of intended stiffness, higher bandgaps are achievable for the same elastic modulus as compared to results of b = 2 . This is due to the fact that topologies obtained based on b = 2 have higher overall stiffness, both shear and normal, which degrades the bandgap efficiency. However, the optimality and supremacy of topologies depend on the desired functionality of PhP. Alternatively the stress state ratio b = 2 is considered for all remaining optimization cases of present work to get almost equally normal-shear stiffened topologies. According to the modal band structures of extreme topologies, all bandgaps are exclusive gaps of first two asymmetric modal branches with intruding symmetric modes. Ideally the topologies TA-1, TB-1 and TC-1 concerning widest RBW should be the same for all cases. However, affected by neighboring topologies, dissimilar topology TB-1 is obtained for b = 2 which shows thinner bandgap at lower frequency range with higher shear modulus compared to TA-1 or TC-1. 4.1.2. Comparison with other relevant works In order to compare the performance of obtained topologies with preceding works on bandgaps of asymmetric waves, the intermediate Pareto topology TA-6 is arbitrarily chosen and remodeled with Polycarbonate solid material of elastic modulus 2.3 GPa, density 1200 kg/m3 and Poisson’s ratio 0.35 and aspect ratio r = 1/0.0909 = 11.00 as per work by Halkjær et al. (2006). Accordingly the relevant frequency RBW of this topology is Δω/ω¯ = 0.328 which is more than double of that of optimized connected topology reported by Halkjær et al. (2006) with Δω/ω¯ = 0.16. Furthermore, the relevant Eigen value RBW of topology TA-6 (with 32 × 32 resolution) is F1 = 0.639 which is around 13% higher than that of optimized topology

Fig. 10. (a) Best Pareto topology for maximized RBW between 2nd and 3rd modal branches of symmetric guided wave modes and (b) relevant modal band structure; symmetric and asymmetric modes are shown by solid and dash lines, respectively, and bandgap of symmetric modes is highlighted in gray.

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

43

obtained by Bilal and Hussein (2012) with F1 = 0.5639 for refined square topology resolution of 64 × 64 . It is noteworthy that the bandgap performance of an intermediate Pareto topology (TA-6), and even not the one with widest RBW (TA-1), was compared with other works aimed at topology with absolutely maximized RBW. Moreover, the topologies of present study are optimized based on slightly different aspect ratio r = 10 < 11. However this comparison can only show the competency of implemented approach in digging out the desired optimum bandgap topology. Whereas the RBW of PhP is highly dependent on its in-plane stiffness, a fair judgment about relative bandgap-structural performance of referred topologies demands their in-plane stiffness to be evaluated and taken into account. 4.1.3. Bandgap Objectives 2 and 3 Alternatively, the stress state ratio b = 2 is considered for all remaining optimization cases of present work to get almost equally normal-shear stiffened topologies. As regard to Bandgap Objective-2 the search for lowest feasible bandgap of symmetric modes between 2nd and 3rd modal branches led to barely opened gap through limited spread of topologies. The best Pareto topology with maximum RBW of Δω/ω¯ = 0.166 and its modal band structure as given in Fig. 10, which show how the desired gap is opened by means of a locally resonant inclusion with weak (hinge) connection to the surrounding solid frame. Due to discreetness of the wave vectors on the border of Brillouin zone to be evaluated (30 points in this study) it is impossible to capture the sharp approaching points of the two modal branches located apart from the search points. This deficiency leads to invalid or overestimated bandgaps when the modes are tangential or gap width is tiny, like all Pareto topologies with face to face connection achieved for this case. Strengthened connectivity can also be obtained by modifying the bandgap objective to search for gaps at higher midgap frequencies (e.g. the work by Halkjær et al. (2006)), though the desired low frequency performance of produced gap is compromised in this way. Hence, preferably, the succeeding gap between 3rd and 4th symmetric modal branches is investigated to be maximized as the Bandgap Objective-2. The optimization results of this case for stress state ratio b = 2 are so given in Fig. 11. Obviously, significantly higher RBW is obtained through wide spread of well-connected topologies. The lowest frequency of bandgap is also less than what obtained for 2nd and 3rd modal branches. The elastic modulus is increased by almost two orders across the Pareto front within the filling fraction range of 0.66–0.92.

(

)

(

)

Fig. 11. Topology optimization results of PhP unitcell with aspect ratio r = 10 for maximized Bandgap Objective-2, stress state ratio b = topologies TD-1 to TD-12 are selected from relevant Pareto front.

2 ; Optimized

44

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. 12. Modal band structure of plate wave modes for (a) topology TD-1, and (b) Topology TD-12 (Fig. 11); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and lowest bandgap of symmetric modes is highlighted in gray.

(

)

(

)

Fig. 13. Topology optimization results of PhP unitcell with aspect ratio r = 10 for maximized Bandgap Objective-3, stress state ratio b = topologies TE-1 to TE-10 are selected from relevant Pareto front.

2 ; Optimized

Modal band structure of extreme topologies TD-1 and TD-12 presented in Fig. 12 show interruption of obtained symmetric bandgap by asymmetric modes and existence of several complete bandgaps of plate wave modes for topology TD-1. Finally the Bandgap Objective-3 is considered to get PhP topologies with maximized complete bandgap of mixed guided wave modes at lowest possible frequency range. The search for existence of any bandgaps within first 6 modal branches resulted in no noticeable gaps for assumed unitcell’s resolution and the first complete gap opened between 6th and 7th modal branches. The

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

45

Fig. 14. Modal band structure of plate wave modes for (a) topology TE-1, and (b) Topology TE-10 (Fig. 13); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and complete bandgap is highlighted in gray.

relevant optimization results for stress state ratio b = 2 are presented in Fig. 13 and modal band structure of extreme topologies TE-1 and TE-10 are given in Fig. 14. Compared to Bandgap Objectives 1 and 2, considerably lower spread of topologies is achieved in the Pareto front around the filling fraction 0.5. This could be due to rarity of topologies with bandgaps of both symmetric and asymmetric modes. Both normal and shear elastic modulus are increased by the order of around 0.5 over the entire range, regardless of the stiffest topology TE-10 which shows a relatively big ramp in both bandgap and stiffness. It is noteworthy that the

(

)

(

)

Fig. 15. Topology optimization results of PhP unitcell with aspect ratio r = 2 for maximized Bandgap Objective-1, stress state ratio b = topologies TF-1 to TF-12 are selected from relevant Pareto front.

2 ; Optimized

46

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. 16. Modal band structure of plate wave modes for (a) topology TF-1, and (b) Topology TF-12 (Fig. 15); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and lowest bandgap of asymmetric modes is highlighted in gray.

(

)

(

)

Fig. 17. Topology optimization results of PhP unitcell with aspect ratio r = 2 for maximized Bandgap Objective-2, stress state ratio b = topologies TG-1 to TG-10 are selected from relevant Pareto front.

2 ; Optimized

modal band structure of stiffest topology TE-10 includes a full bandgap of 2nd and 3rd symmetric modes. Comparing the modal band structures given in Figs. 5, 7 and 9 it is evident that the RBW of topologies for Bandgap Objective-1 are changed within the Pareto front mostly by the width of bandgap at almost the same frequency range. In contrary, the RBW of Bandgap Objectives 2 and 3 (Figs. 12 and 14) vary by both width and mid frequency of bandgaps.

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

47

Fig. 18. Modal band structure of plate wave modes for (a) topology TG-1, and (b) Topology TG-10 (Fig. 17); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and lowest bandgap of symmetric modes is highlighted in gray.

(

)

(

)

Fig. 19. Topology optimization results of PhP unitcell with aspect ratio r = 2 for maximized Bandgap Objective-3, stress state ratio b = topologies TH-1 to TH-10 are selected from relevant Pareto front.

2 ; Optimized

4.2. Thick PhP unitcell with aspect ratio r = 2 The optimum design of thick PhP unitcell with aspect ratio r = 2 is studied in this section. Higher thickness generally provides higher flexural stiffness for the PhP structure of the same lattice size. Moreover, thicker PhP principally shows more controllability on the lower frequency oscillations which have larger wavelength comparable to the plate's thickness. The

48

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. 20. Modal band structure of plate wave modes for (a) topology TH-1, and (b) Topology TH-10 (Fig. 19); symmetric and asymmetric modes are shown by solid and dash lines, respectively, and lowest complete bandgap is highlighted in gray.

topology optimization results for the three bandgap objectives based on stress state ratio b = 2 are presented in Figs. 15–20. For Bandgap Objective-1 the RBW of widest bandgap topology TF-1 (Fig. 15) is almost double of corresponding topology TA-1 for thin PhP (Fig. 4) around almost the same midgap frequency. Despite considering balanced stress state ratio b = 2 , the shear modulus shows considerably higher gradient compared to elastic modulus. The shear modulus is increased by around 2 orders across the Pareto front while the elastic modulus shows slight variation for a large portion of topologies and relatively lower rise at the end (Fig. 15(c)). It is worthy note that the normal stiffness is governed by elastic modulus and Poisson's ratio and so can vary by either of these elastic properties. As an example the topology TF-9 shows an exceptional mode of stiff topology within the Pareto front which possesses relatively higher elastic modulus but instead lower Poisson’s ratio (Fig. 15(b) and (c)). The modal band structures of extreme topologies TF-1 and TF-12 presented in Fig. 16 show the relevant exclusive bandgaps of first two asymmetric modal branches. Regarding the Bandgap Objective-2, dislike thin PhP of r = 10, the optimization for lowest possible bandgap between 2nd and 3rd symmetric modes led to acceptable topologies with adequate internal connections as demonstrated in Fig. 17(a). The shear modulus and elastic modulus increase by around 1.5 and 1 orders of magnitude, respectively, throughout the Pareto front around the filling fraction 0.6. Shear modulus and Poisson’s ratio both show a significant ramp when the topology changes from TG-4 to TG-5 with obviously different morphologies leading to higher in-plane stiffness. The modal band structure of widest bandgap topology TG-1 and stiffest topology TG-10 illustrated in Fig. 18 clearly show the relevant exclusive bandgap of desired symmetric modes. Of course, for the stiffest topology TG-10 the desired bandgap is almost closed and a wide exclusive bandgap is present between succeeding couple of symmetric modes. Finally the optimized topologies of thick PhP are studied for maximized Bandgap Objective-3 to get maximized RBW of mixed plate wave modes. Similar to the thin PhP, no noticeable bandgap was opened within first 6 modal branches during optimization studies and the first complete gap of plate waves was found to be between 6th and 7th modal branches. The relevant results of optimization are given in Figs. 19 and 20. Accordingly, the RBW of widest bandgap topology TH-1 (Fig. 19(a)) is slightly more than its corresponding thick PhP topology TE-1 (Fig. 13(a)) but with considerably higher gap width at higher frequency range. Based on the results of Fig. 19(b) the elastic modulus and shear modulus are increased by about 0.8 and 1.5 orders along the Pareto front over the filing fraction range 0.54–0.77. The topology mode is changed from TH-4 to TH-5 so that the Poisson's ratio is increased providing higher normal stiffness. The dip in elastic modulus corresponds to the topology TH-5 with maximal Poisson's ratio instead. The modal band structure of topologies TH-1 and TH-10 as shown in Fig. 20 demonstrates the complete gap of plate wave modes between 6th and 7th modal branches. The RBW is decreased by narrowing the gap width and around almost the same midgap frequency. An additional wide complete bandgap is present for widest band gap topology TH-1 between 10th and 11th modes.

5. Frequency response of finite PhP structure The optimum topology of PhP unitcell for bandgaps of fundamental guided waves was investigated in the preceding section. The modal band structure and RBW were calculated through Bloch-Floquet perfectly periodic boundary condition based on the assumption of infinite periodicity for the unitcell. However, in the real world, PhP structures with finite boundaries are to be used in different applications and it is essential to realize the validity and performance of achieved bandgap topologies in realistic conditions. Therefore representative finite thin PhP structures are modeled made by selected Pareto topologies of the 3 bandgap objectives and their steady state and transient frequency responses are studied.

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

49

Fig. 21. Modal band structure and mode shapes of plate wave modes at X for selected topologies (a) TA-6, (b) TD-6 and (c) TE-5; symmetric and asymmetric modes are shown by solid and dash lines, respectively, and bandgap highlighted in gray.

Intermediate topologies TA-6 (Fig. 4), TD-6 (Fig. 11) and TE-5 (Fig. 13) are chosen from the Pareto fronts of bandgap objectives 1, 2 and 3, respectively, with moderate stiffness and RBW. Lattice periodicity of a = 10 mm is assumed leading to plate’s thickness h = 1 mm for aspect ratio r = 10.

50

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. 22. (a) Finite square PhP structure with 20 cells of topology TA-6 (Fig. 4(a)) , r = 10 and h = 1mm on each side, and (b) its reduced symmetric model excited at the central point O, to be captured on straight sampling points A and B and diagonal ones I and J.

5.1. Modal band structure of selected topologies Initially the frequency band structure of selected thin PhP topologies with assumed structural thickness of 1 mm is calculated and the relevant bandgap frequencies are defined as shown in Fig. 21. The insets of Fig. 21 shows mode shapes of first few modes at arbitrary Brillouin zone point X representing wave modes corresponding to the associated modal branches. Accordingly, the lowest bandgap of asymmetric modes for topology TA-1 is opened between the fundamental asymmetric Lamb mode A0 and its resonated branch by Bragg reflection producing bandgap (Fig. 21(a)). Apparently this exclusive gap of asymmetric modes is interrupted by first symmetric shear horizontal mode SHS0. Fig. 21(b) shows the exclusive bandgap of symmetric modes for topology TD-6 encompassing several asymmetric modes. The lower symmetric modes of this gap are SHS0 and its folded branch as well as fundamental symmetric Lamb mode S0. The gap is confined on the upside by relevant resonated modes. The modal band structure of topology TE-5 is also given in Fig. 21(c) including first complete gap of mixed plate waves. The gap is surrounded by fundamental plate wave modes A0, S0, SHS0 and SHA0 from the downside and their resonated branches from the upside. 5.2. Steady state frequency response An original square plate lattice structure with 20 unitcells on each side is considered to be excited at the center (Fig. 22 (a)). Hence, the relevant reduced symmetric FEM model of assumed plate structure is developed like the one shown in Fig. 22(b) for selected topology TA-6. Unlike optimization procedure, no compliant material is considered and the solid part of topology is modeled only. The model is subjected to a probing load at the central point O. Apparently the geometrically central point O of lattice structure of selected topologies is void and so the diagonally nearest solid node is excited instead. Then the transmission of stimulated oscillations is captured at 4 different positions within the plate structure. Two straight points A and B as well as two diagonal points I and J are selected so that the transmission of plate waves to different distances in different angles is determined (Fig. 22). Harmonic load is applied to excitation point O of finite PhP structures over definite frequency range including the corresponding bandgap and then the resultant steady state frequency response is defined. Fig. 23 demonstrates the relative transmittance of excited load to introduced sampling points A, B, I and J (Fig. 22) over the selected frequency range in logarithmic scale dB. In-phase transversal loads (in z axis) and planar loads (in both x and y axes) are applied to both upper and lower sides of the point O for dominant excitation of asymmetric modes (Fig. 23(a)) and symmetric modes (Fig. 23(b)), respectively. Only upper side of excitation point O is stimulated transversally for generation of mixed plate wave modes (Fig. 23(c)). The transmittance of in-plane polarizations (with index xy) and transversal polarizations (with index z) are so defined based on the nature of excited wave. As for PhP structure of asymmetric bandgap topology TA-1, the transmission spectrum (Fig. 23(a)) of transversal polarization is highly attenuated in the same bandgap frequency predicted by the modal band structure of Fig. 21(a) in around 22–32 kHz. In order to model the finite structure of symmetric bandgap topology TD-6 a polished topology as shown in the inset of Fig. 21 (b) is employed. The transmission spectrum of Fig. 23(b) for in-plane polarization over 0–500 kHz presents the same wide bandgap as shown in the modal band structure of Fig. 21(b) in the frequency range 165–370 kHz. As expected, the lowest and the highest attenuations occur at sampling points A and J, respectively, due to their relative distance to the excitation point O. Lastly, the transmittance of both transversal and in-plane polarizations is determined for the PhP structure of topology TE-5 with complete gap of mixed plate waves. Despite the asymmetric nature of aforementioned loading for this case, the

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

51

Fig. 23. Frequency response of finite PhP structure with selected optimized topologies (a) TA-6, (b) refined TD-6 and (c) TE-5 in logarithmic scale dB; Left: transmission of in-plane (xy) and/or anti-plane (z) polarizations, Right: transmission contour of displacement vector at selected bandgap frequency.

related transmission spectrum given in Fig. 23(b) shows dominant symmetric excitation by a wide bandgap corresponding to the calculated gap of symmetric modes (Fig. 21(c)) over around 50–140 kHz. However, the spectrum of transversal

52

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

polarization which principally represents the oscillations of asymmetric modes is resonated at frequency 83 kHz corresponding to presence of resonated A0 mode at Γ(0, 0). For more clarity pure asymmetric modes are dominantly excited and the transmission spectrum to point J is calculated (index z2 in Fig. 23(c)). Expectedly the asymmetric bandgap of Lamb mode A0 in the frequency range 40–80 kHz is evident. Based on the relevant modal band structure of Fig. 21(c), two asymmetric modes A0 and SHA0 interrupt the symmetric bandgap. However, the mode SHA0 obviously was not dominant in either examined loading cases. The transmission contour of displacement vector at specified bandgap frequency is also given in any case to show the spacial gradient of frequency response. The contours clearly show the intense decay of oscillations’ amplitude all over the structure in all directions through successive number of PhP unitcells. 5.3. Transient frequency response The transient frequency response of finite PhP structure is also defined further for asymmetric bandgap topology TA-1 (Fig. 22) to inspect the attenuation of asymmetric modes as well as passage of symmetric modes within bandgap frequency. Guided wave modes are launched via a 3 cycles sinusoidal toneburst of bandgap frequency 27.5 kHz modified by Hanning window for high power concentration around this central frequency. The selected central frequency 27.5 kHz corresponds to the dip in Fig. 23(a). As for asymmetric modes, transversal displacement load wf with amplitude w0 = 0.01mm is applied to

Fig. 24. Transient response of the PhP structure of Fig. 22 under transversal Hanning window toneburst at central asymmetric bandgap frequency 27.5 kHz (Fig. 23a); time history and Gabor wavelet spectrum of excited and transmitted transversal oscillation w.

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

53

excitation point O as follows (Hedayatrasa et al., 2014, Veidt and Ng, 2011):

⎛ ⎛ ωt ⎞⎞ wf = ( w0 sin ωt )⎜ 0.5⎜ 1 − cos ⎟⎟ Nc ⎠⎠ ⎝ ⎝

(18)

over the time t = 109μs corresponding to Nc = 3 cycles of frequency 27.5 kHz. The second term of Eq. (18) is Hanning window function applied to the sinusoidal signal. Then the transient response of structure is defined up to t = 1200μs to let the asymmetric wave modes propagate throughout the plate. Displacement load is preferably applied to the structure to make sure the desired bandgap frequency is excited adequately despite natural resistance of structure to vibrate within gap frequency range. When force load was examined, the frequency spectrum naturally deviated from the intended bandgap frequency 27.5 kHz to around 22 kHz due to resonance of A0 mode confining lower gap frequency at 22 kHz. The time history of normalized amplitude of excited oscillations at point O and transmitted oscillations at points A and I are given in Fig. 24(a). The initial part of signal concerning loading point O shows the 3 cycles windowed exciting toneburst. The time– frequency spectrum of the signals is then obtained through Gabor wavelet transform. Wavelet transform is a great tool for time–frequency spectrum decomposition of transient response and among various mother wavelets the Gabor wavelet produces the best resolution in both time and frequency domains (Kishimoto et al., 1995, Quek et al., 2001, Hedayatrasa

Fig. 25. Transient response of the PhP structure of Fig. 22 under in-plane Hanning window toneburst at central asymmetric bandgap frequency 27.5 kHz (Fig. 23(a)); time history and Gabor wavelet spectrum of excited and transmitted in-plane oscillation in x–y plane.

54

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

et al., 2014, Ng et al., 2009). Gabor wavelet is actually a harmonically modulated Gaussian window in the complex form:

ψ (t ) =

4

π

2 −0.5( ωγ (t − tc ))2 iω(t − tc ) e e γ

(19)

where γ = π 2/ ln 2 and ω is circular modulation frequency and tc the temporal center of Gaussian function. Fig. 24 (b) represents the spectrum of excitation point which clearly shows power concentration within bandgap frequency and central frequency 27.5 kHz. However, in spite of intense excitation of bandgap frequency, the spectrum of transmitted oscillations to sampling points A (Fig. 24(c)) and I (Fig. 24(d)) confirms the strong attenuation of waves in the time domain within bandgap frequency range as compared to relevant band structure in Fig. 21(a). The logarithmic time–frequency spectrums of transmitted signals to points A and I are also given in Fig. 22(e) and (f), respectively, to discern the order of magnitude of spectrum. Wave attenuation in point A is obviously higher than point I, despite lower distance of point A to excitation point. This observation can be justified by wider bandgap width of asymmetric modes along Brillouin border Γ − X which corresponds to wave propagation through longitudinal axis x including point A. Symmetric guided wave modes are stimulated and analyzed further via the same procedure explained for asymmetric modes, but this time by exciting both in-plane polarizations u and v at point O. The time history of normalized amplitude of resultant excited oscillations at point O and transmitted oscillations at points A and I are given in Fig. 25(a). The bandgap frequency 27.5 kHz is properly excited according to the excitation spectrum given in Fig. 25(b) during excitation time of t = 109μs . The time–frequency spectrums of resultant in-plane oscillations transmitted to points A (Fig. 25(a) and (e)) and J (Fig. 25(a) and (f)) represent robust excitation and passage of symmetric modes throughout the PhP structure within bandgap frequency range 20–30 kHz of asymmetric modes. Despite power concentration of launched wavelet at bandgap frequency 27.5 kHz, the central frequency of transmitted waves is slightly changed. This is due to existence of partial bandgaps of symmetric modes acting at particular wave vectors, i.e. propagation directions, and resonance of relevant modal branches (Fig. 21(a)). The spectrum of diagonal point I is more noticeably raised to around 35 kHz. This observation can be explained by resonance of two SHS0 modal branches approaching at around frequency 35 kHz corresponding to symmetric Brillouin zone point M (Fig. 21(a)). However, transmitted symmetric waves to both points A and I show high spectrum power inside asymmetric bandgap frequency range 20–30 kHz.

6. Conclusion Evolutionary GA multi-objective optimization of PhP unitcell with square symmetry was performed successfully in order to maximize RBW of fundamental guided wave modes with maximal in-plane stiffness; where the obtained topology can be produced by perforation of a uniform solid plate. The optimum topology of such single material porous plate structure for widest bandgap naturally favors maximum isolation of scattering segments achievable through weak structural connections. Hence, the strain energy compliance of unitcell was considered as the second objective to be minimized to get bandgap topologies with optimized structural worthiness i.e. maximized in-plane stiffness. Complete bandgaps of guided waves as well as exclusive bandgaps of symmetric and asymmetric guided wave modes were explored. The resultant PhP structures have promising application in selective filtration, excitation, steering and capturing of relevant modes to be used in various areas e.g. structural health monitoring, acoustic energy harvesting or radio frequency communications. Polysilicon solid background was considered and very compliant material was assumed to model the void regions of topology for simplified FEM modeling of randomly generated topologies. The results of this study indicated that reducing the solid material properties, i.e. density and elastic modulus, by order of magnitude of more than 5 to be used as compliant material provides accurate results for modal band analysis of relevant porous unitcell. Bandgap topologies of fundamental guided wave modes were achieved successfully for both assumed aspect ratios. Thick PhP with aspect ratio 2 generally provides higher RBW compared to thin PhP with aspect ratio 10. Moreover, obtained Pareto solutions clearly showed degradation of bandgap efficiency by increasing the in-plane stiffness. Optimized topologies for bandgaps of asymmetric wave modes indicated superior bandgap efficiency compared to preceding works available in the literature. Different in-plane stress state conditions were also examined and proved to significantly influence relevant Pareto topologies. Finally, frequency response of thin finite PhP structures produced by selected intermediate Pareto topologies of symmetric, asymmetric and complete bandgaps were defined and validity of relevant topologies was confirmed. High attenuation of relative transmission to different sampling points throughout the plate structure was observed within relevant bandgap frequency under steady state harmonic excitation. Moreover, the transient response of finite PhP structure of asymmetric bandgap subjected to windowed toneburst was calculated. Time–frequency spectrum of transient response obtained via Gabor wavelet transform evidently indicated the attenuation of asymmetric wave modes and passage of symmetric wave modes in the time domain within relevant asymmetric bandgap frequency. Unitcell topology with square symmetry was studied in present work. However, it is worthwhile to explore bandgap topologies with reduced symmetric constrains e.g. rectangular symmetry which might have better bandgap-structural performance.

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

55

Appendix A. Transversal mesh resolution Proper mesh resolution along the thickness of unitcell h relies on its transversal aspect ratio r = a/h and the frequency range of highest plate wave mode to be calculated. SOLSH190 element usually gives more accurate predictions compared to classical shell elements (with Mindlin's Plate formulation) when the shell is thick . Thus, for modeling moderately thick PhP unitcells of aspect ratio r = 10 a single layer of solid-shell elements provides good accuracy for modal band analysis of low

Fig. A.1. Gradient of first few modal frequencies versus number of element layers through thickness of connected topology shown in Fig. B. 1(a) for transversal aspect ratio r = 2 and Ps/Pv = 1000 ; (a) whole unitcell modeled for mixed plate modes analysis, and (b) half of unitcell’s thickness modeled from the midplane for decoupling asymmetric modes.

order modal branches. As regard to thick unitcell with aspect ratio r = 2, Fig. A1 displays the modal frequencies at e.g. Brillouin zone point X (Fig. 2) versus the number of element layers along the thickness obtained for an arbitrary topology shown in Fig. A1(a) and material property ratio Ps/Pv = 1000 discussed later in Appendix B. Fig. A1(a) and (b) show the first 10 mixed plate wave modes and first 10 asymmetric modes, respectively (9 modes visible due to coincidence). Obviously, the steepest ramp occurs when increasing element layers from one to two, and the calculated asymmetric modes are less sensitive to transversal mesh resolution (Fig. A1(b)). This is due to the fact that just one half of plate’s thickness is modeled for mode decoupling which provides more accuracy compared to the model with one element layer for the entire thickness. In both cases the highest mode converges well after considering up to 10 layers. Although increasing the number of element layers leads to more accurate modal analysis, considering 1 element layer proved to be adequate for estimation and widening of lowest bandgap during optimization. Two element layers are considered only for the unitcell of aspect ratio r = 2 when maximizing RBW of complete plate wave's gap. Of course, for calculation of the modal band structure and RBW of achieved topologies higher mesh resolutions of 10 and 4 layers are considered for aspect ratios 2 and 10, respectively.

Appendix B. Modeling void segments and compliant material properties As a general approach to avoid complexity in modeling randomly generated topologies, void sections of porous PhP unitcell can be substituted by very compliant elements compared to solid elements. This approach has been frequently employed in the literature for modal band analysis and optimization of porous PhCrs but its challenges and factors influencing its efficiency have never been discussed. To give some sense, by considering air filled porosity with density ρv = 1.29(kg /m3) and sound velocity Cv = 343(m /s ) its equivalent elastic modulus would be Ev = 151.8(kPa) and thus ρs /ρv = 1.8 × 103 and Es/Ev = 1.11 × 106 . The appropriate order of degrading material properties of solid material to be used as compliant material is an essential feature of FEM modal band analysis. In order to investigate this matter, the sensitivity of modal frequencies of plate unitcell to the material property ratio Ps/Pv is studied as shown in Fig. B1. Unitcell with aspect ratio r = 2 and transversal mesh resolution 10 is modeled for this purpose. The variables Ps and Pv stand for material properties (density or elastic modulus) of solid and compliant void phases, respectively. Elastic modulus and density of solid material both are reduced with the same order to be considered as compliant material. So extraordinary mass to stiffness ratio leading to virtual low modal frequencies is avoided (Bendsoe and Sigmund, 2003).

56

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Fig. B.1. Topology and related modal response of unitcell at X (π /2, 0) versus ratio of material properties Ps/Pv , as well as the modal band structure of unitcell for Ps/Pv = 109 ; (a)–(c) connected topology; (d)–(f) disconnected topology; and (g)–(i) connected topology including disconnected segments.

With regard to the random nature of GA optimization and relevant individual topologies to be analyzed during fitness evaluation process, three distinct arbitrary topologies are studied as delivered in Fig. B1(a), (d) and (g). Then the gradient of dimensionless modal frequency ωd versus Ps/Pv at arbitrary Brillouin zone point X (π /2, 0) is presented in logarithmic scale in Fig. B1(b), (e) and (h) for these topologies, respectively. The first topology illustrated in Fig. B1(a) contains just one solid inclusion of filling fraction vf = 0.375 by which obviously a connected lattice plate could be produced. According to the relevant results in Fig. B1(b) the modal frequencies decrease dramatically until around Ps/Pv = 1000 and after slight variation up to around Ps/Pv = 105 show absolutely no variation any further. But, by insertion of an arbitrary disconnection in the topology as shown in Fig. B1(d), a bundle of low frequency modes keep reducing and tend to zero by increasing the ratio Ps/Pv (Fig. B1(e)). However, higher modal frequencies remain constant after an initial descend up to around Ps/Pv = 103 (Fig. B1 (e)). Fig. 3(g) is related to a connected topology similar to Fig. 3(a) but with a randomly added disconnected solid pixel. The gradient of modal frequency versus Ps/Pv for this case (Fig. B1(h)) is generally like the one for disconnected topology (Fig. B1 (e)). In fact, for material property ratios up to around 1000, the fundamental modal frequencies of unitcell are governed by both solid and compliant material phases. After this limit the compliant elements become extremely degraded causing the modal frequencies governed by these elements to tend to zero. Indeed, in this case the solid parts of topology are not supported by compliant elements anymore and their modal frequencies remain unchanged by further reduction of material property ratio (just like a model with no void elements). Consequently, by considering a high ratio Ps/Pv > 105 the modal response of connected unitcell topology can be accurately defined. However, with regard to modal band structures (excluding near zero modes) given in Fig. B1(c), (f) and (i) obtained

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

57

by high ratio Ps/Pv = 109 the challenge is how to distinguish modal branches of connected segment if any. Otherwise fictitious RBW will be calculated based on unvarying modes of disconnected segments. Principally disconnected segments of topology show boundary (wave number) independent modal response as confirmed by the modal band structure delivered in Fig. B1(f) and (i) containing virtual bandgaps. But if any connected segment touching the boundary is present in the topology, the relevant modal response is wave number variant (Fig. B1(c) and (i)) and uniquely has zero modal frequencies at Γ(0, 0) supporting the rigid body motion modes in three degrees of freedom. Considering relatively low material property ratio of around 1000 generally provides a good approximation of modal band structure for connected unitcell in the absence of instability and near zero modal frequencies. But, using this low ratio leads to impractical optimized topologies comprised of isolated segments predominantly affecting the produced bandgap (Halkjær et al., 2006). Actually the problem converges towards topologies producing bandgaps by locally resonant solid inclusions held by surrounding weak compliant elements. Hence, for accurate modal band analysis of the porous unitcell, sufficiently high degradation ratio Ps/Pv > 105 has to be considered, or void elements and relevant boundary conditions have to be removed from the model. However, in any case, a robust topology recognition and processing as implemented in present study is essential to assure the feasibility and validity of optimized topologies.

References Aberg, M., Gudmundson, P., 1997. The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure. J. Acoust. Soc. Am. 102, 2007. Bendsoe, M.P., Sigmund, O., 2003. Topology Optimization (Theory, Methods and Applications). Springer, Germany. Bilal, O.R., Hussein, M.I. 2012. Topologically evolved phononic material: Breaking the world record in band gap size, In: Photonic and Phononic Properties of Engineered Nanostructures. Chapman, C.D., Saitou, K., Jakiela, M.J., 1994. Genetic algorithms as an approach to configuration and topology design. J. Mech. Des. 116, 1005–1012. Charles, C., Bonello, B., Ganot, F., 2006. Propagation of guided elastic waves in 2D phononic crystals. Ultrasonics 44, e1209–e1213. Denies, J., Dehez, B., Glineur, F., Ahmed, H. B. 2012. Genetic algorithme-based topology optimization: performance improvement through dynamic evolution of the population size. International Symposium on Power Electronics, Electrical Drives, Automation and Motion. Deymier, P., 2011. Acoustic Metamaterials and Phononic Crystals. Springer, Germany. Dong, H.-W., Su, X.-X., Wang, Y.-S., 2014a. Multi-objective optimization of two-dimensional porous phononic crystals. J. Phys. D: Appl. Phys. 47, 155301. Dong, H.-W., Su, X.-X., Wang, Y.-S., Zhang, C., 2014b. Topological optimization of two-dimensional phononic crystals based on the finite element method and genetic algorithm. Struct. Multidiscip. Optim., 1–12. Dong, H.-W., Su, X.-X., Wang, Y.-S., Zhang, C., 2014c. Topology optimization of two-dimensional asymmetrical phononic crystals. Phys. Lett. A 378, 434–441. Gazonas, G.A., Weile, D.S., Wildman, R., Mohan, A., 2006. Genetic algorithm optimization of phononic bandgap structures. Int. J. Solids Struct. 43, 5851–5866. Gonella, S., Ruzzene, M., 2008. Homogenization and equivalent in-plane properties of two-dimensional periodic lattices. Int. J. Solids Struct. 45, 2897–2915. Hajela, P., Lee, E., Lin, C.Y., 1993. Genetic algorithms in structural topology optimization. In: Bendsøe, M., Soares, C.M. (Eds.), Topology Design of Structures. Springer, Netherlands. Halkjær, S., Sigmund, O., Jensen, J.S., 2006. Maximizing band gaps in plate structures. Struct. Multidiscip. Optim. 32, 263–275. Hedayatrasa, S., Bui, T.Q., Zhang, C., Lim, C.W., 2014. Numerical modeling of wave propagation in functionally graded materials using time-domain spectral Chebyshev elements. J. Comput. Phys. 258, 381–404. Hedayatrasa, S., Abhary, K., Uddin, M., 2015. Numerical study and topology optimization of 1D periodic bimaterial phononic crystal plates for bandgaps of low order Lamb waves. Ultrasonics 57, 104–124. Hussein, M., Hamza, K., Hulbert, G., Scott, R., Saitou, K., 2006. Multiobjective evolutionary optimization of periodic layered materials for desired wave dispersion characteristics. Struct. Multidiscip. Optim. 31, 60–75. Hussein, M.I., Hamza, K., Hulbert, G.M., SAITOU, K., 2007. Optimal synthesis of 2D phononic crystals for broadband frequency isolation. Waves Random Complex Media 17, 491–510. Karmazin, A., 2013. Time-Efficient Simulation of Surface-Excited Guided Lamb Wave Propagation in Composites. KIT Scientific Publishing, Germany. Kaw, A.K., 2005. Mechanics of Composite Materials. CRC press, USA. Khelif, A., Aoubiza, B., Mohammadi, S., Adibi, A., Laude, V., 2006. Complete band gaps in two-dimensional phononic crystal slabs. Phys. Rev. E 74, 046610. Kishimoto, K., Inoue, H., Hamada, M., Shibuya, T., 1995. Time frequency analysis of dispersive waves by means of wavelet transform. J. Appl. Mech. 62, 841–846. Kittel, C., 1986. Introduction to Solid State Physics. Wiley, New York. Krushynska, A.O., Kouznetsova, V.G., Geers, M.G.D., 2014. Towards optimal design of locally resonant acoustic metamaterials. J. Mech. Phys. Solids 71, 179–196. Kunakote, T., Bureerat, S., 2011. Multi-objective topology optimization using evolutionary algorithms. Eng. Optim. 43, 541–557. Kundu, T., 2004. Ultrasonic Nondestructive Evaluation – Engineering and Biological Material Characterization. CRC press. Li, J., Wang, Y.-S., Zhang, C., 2009. Finite element method for analysis of band structures of phononic crystal slabs with archimedean-like tilings. In: Ultrasonics Symposium (IUS), 2009 IEEE International. IEEE, Rome, pp. 1548–1551. Lin, C.-M., Hsu, J.-C., Senesky, D. G., Pisano, A. P. 2014. Anchor loss reduction in ALN Lamb wave resonators using phononic crystal strip tethers. Frequency Control Symposium (FCS), 2014 IEEE International. Lin, S.-C.S., Huang, T.J., Sun, J.-H., Wu, T.-T., 2009. Gradient-index phononic crystals. Phys. Rev. B 79, 094302. Liu, Z.-F., Wu, B., He, C.-F., 2014. Band-gap optimization of two-dimensional phononic crystals based on genetic algorithm and FPWE. Waves Random Complex Media 24, 1–20. Manktelow, K.L., Leamy, M.J., Ruzzene, M., 2013. Topology design and optimization of nonlinear periodic materials. J. Mech. Phys. Solids 61, 2433–2453. Mohammadi, S., 2010. Phononic Band Gap Micro/Nano-mechanical Structures for Communications and Sensing Applications. Georgia Institute of Technology, United States. Ng, C.-T., Veidt, M., Rajic, N. 2009. Integrated piezoceramic transducers for imaging damage in composite laminates. In: proceedings of the 2nd International Conference on Smart Materials and Nanotechnology in Engineering, International Society for Optics and Photonics, pp. 74932M–74932M-8. Olhoff, N., Niu, B., Cheng, G., 2012. Optimum design of band-gap beam structures. Int. J. Solids Struct. 49, 3158–3169. Oliveira, J., Pinho-Da-Cruz, J., Teixeira-Dias, F., 2009. Asymptotic homogenisation in linear elasticity. Part II: finite element procedures and multiscale applications. Comput. Mater. Sci. 45, 1081–1096. Olsson iii, R., El-kady, I., 2009. Microfabricated phononic crystal devices and applications. Meas. Sci. Technol. 20, 012002.

58

S. Hedayatrasa et al. / J. Mech. Phys. Solids 89 (2016) 31–58

Pennec, Y., Vasseur, J.O., Djafari-Rouhani, B., Dobrzyński, L., Deymier, P.A., 2010. Two-dimensional phononic crystals: examples and applications. Surf. Sci. Rep. 65, 229–291. Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evolut. Comput. 6, 182–197. Quek, S., Wang, Q., Zhang, L., Ong, K., 2001. Practical issues in the detection of damage in beams using wavelets. Smart Mater. Struct. 10, 1009. Rupp, C., Evgrafov, A., Maute, K., Dunn, M., 2007. Design of phononic materials/structures for surface wave devices using topology optimization. Struct. Multidiscip. Optim. 34, 111–121. Sharpe, W. N., Yuan, B., Vaidyanathan, R., Edwards, R. L. 1997. Measurements of Young modulus, Poisson ratio, and tensile strength of polysilicon. Micro Electro Mechanical Systems, MEMS'97, Proceedings, IEEE., Tenth Annual International Workshop on. IEEE. Sigmund, O., Jensen, J.S., 2003. Systematic design of phononic band-gap materials and structures by topology optimization. Phil. Trans. R. Soc. Lond. 361, 1001–1019. Su, Z., Ye, L., 2009. Identification of Damage Using Lamb Waves: From Fundamentals to Applications. Springer, Germany. Veidt, M., Ng, C.-T., Hames, S., Wattinger, T., 2008. Imaging laminar damage in plates using Lamb wave beamforming. Adv. Mater. Res. 47, 666–669. Veidt, M., Ng, C.-T., 2011. Influence of stacking sequence on scattering characteristics of the fundamental anti-symmetric Lamb wave at through holes in composite laminates. J. Acoust. Soc. Am. 129, 1280–1287. Wang, P., Shim, J., Bertoldi, K., 2013. Effects of geometric and material nonlinearities on tunable band gaps and low-frequency directionality of phononic crystals. Phys. Rev. B 88, 014304. Wang, S.Y., Tai, K., 2005. Structural topology design optimization using Genetic Algorithms with a bit-array representation. Comput. Methods Appl. Mech. Eng. 194, 3749–3770. Weise, T. 2006. Global Optimization Algorithms – Theory and Application, 〈http://www.it-weise.de/〉. Wu, T.-T., Hsu, J.-C., Sun, J.-H., 2011. Phononic plate waves. Ultrason. Ferroelectr. Freq. Control IEEE Trans. 58, 2146–2161. Zhu, H., Semperlotti, F., 2013. Metamaterial based embedded acoustic filters for structural applications. AIP Adv., 3.