Journal of the Mechanics and Physics of Solids Tunable phononic bandgap materials designed via topology optimization

Topology optimization is used to design phononic bandgap materials that are tunable by mechanical deformation. A periodic media is considered, which due to the assumption of length scale separation, allows the dispersion relations to be obtained by analyzing a single unit cell subjected to Floquet–Bloch boundary conditions. A finite macroscopic deformation is applied to the unit cell to affect its geometry and hence dispersion. We tune the dispersion–deformation relation to our liking by solving a topology optimization problem using nonlinear programming. The adjoint method is employed to compute the sensitivities, and the non-differentiability of degenerate eigenvalues is avoided using symmetric polynomials. Several tunable phononic crystal designs are presented. Also, a verification analysis is performed, wherein the optimized design is interpreted and analyzed using a conforming finite element mesh.


Introduction
Phononic crystals are periodic structures engineered to control the propagation of elastic waves. The ability to tailor the wave propagation properties is partly attributed to Bragg scattering, in which a periodic arrangement of scatters allows waves with certain frequencies to propagate, whereas others, with frequencies in the so-called bandgap, cannot (Brillouin, 1953). An incident wave might also give rise to local resonance within the crystal, such that the energy is captured in the resonators and no propagation occurs (Liu et al., 2000;Raghavan and Phani, 2013).
Various applications can benefit from the control of elastic wave propagation. These include waveguides (Vasseur et al., 2007), sensors (Lucklum and Li, 2009) and multiplexing devices (Pennec et al., 2004). The investigation of bandgap elastic structures was first considered by Sigalas and Economou (1992) and Kushwaha et al. (1993), who found that bandgaps could be obtained by embedding stiff cylinders in a compliant matrix. The frequency range of the bandgaps can be tailored by the choice of the constituent materials, the lattice topology (square, hexagonal, etc.) and by the morphology of the inclusions.
The systematic design of phononic crystals is an obvious candidate for optimization. This was exemplified in the seminal work of Sigmund and Jensen (2003), whereby the relative size of bandgaps was maximized using gradient based topology optimization. The analogous design for photonic bandgaps, i.e. of electromagnetic wave propagation materials, was shortly thereafter investigated by Borel et al. (2004). Other examples of bandgap maximization using topology optimization include e.g. Hussein et al. (2007), Rupp et al. (2007), Men et al. (2014), Wormser et al. (2017) and Liu et al. (2020) and Swartz et al. (2021). A detailed literature review appears in Li et al. (2019).
Recently, the tunability of phononic crystals has attracted attention, wherein bandgaps are controlled by applying external stimuli in the form of e.g. mechanical deformation (Bertoldi and Boyce, 2008a), temperature variations (Jim et al., 2009) and piezoelectric effects (Hou et al., 2004). We focus on the former in this work. The effect of prestress on bandgaps was first investigated by Bertoldi and Boyce (2008a,b), where microstructural elastic instabilities were shown to disrupt the lattice periodicity, effectively allowing the bandgaps to open or close. The analysis was based on the small-on-large deformation framework (Ogden, 2007), wherein an incremental deformation is superimposed on a prestressed configuration in quasi-static equilibrium. This work was extended by Wang et al. (2014), whereby local resonance was considered, and by Zhang and Parnell (2017), whereby deformation independent bandgaps were obtained.
Topology optimization is a promising tool for tailoring phononic crystals which are tunable through elastic deformation. The method was pioneered by Evgrafov et al. (2008), for designing tunable wave-guides using a gradient based topology optimization framework. Hedayatrasa et al. (2016) demonstrated the approach using a genetic optimization algorithm (GA). A GA was also used by Bortot et al. (2018), to design dielectric elastomers actuated by external electric fields. The GA studies are subject to coarse design resolutions for obvious reasons. To address this limitation, the use of gradient based topology optimization is a must. To this end, the framework developed by Bortot et al. (2018) was repeated with gradient base optimization in the recent work by Sharma et al. (2021).
In this study, we use density based topology optimization to design hyperelastic phononic crystals with tunable wave propagation properties. More specifically, we generate designs for which a maximum bandgap between two preselected modes exists in the deformed configuration, but not in the undeformed configuration. We assume separation of length scales in the micro-to-macro transition to limit the analysis to a single unit cell. The deformed configuration is obtained by applying periodic boundary conditions on the unit cell, and the dispersion relations are identified by applying Floquet-Bloch boundary conditions to a sequence of generalized eigenvalue problems. The numerical analysis is based on a standard finite element framework. To solve the optimization problem, we employ the gradient based method MMA (Method of Moving Asymptotes, cf. Svanberg, 1987). This is a numerical study, i.e. we do not validate these designs. However, we believe that existing fabrication methods would allow us to do so in future research, cf. e.g. Vasseur et al. (1998).

Governing equations
We consider a macroscopic body with a periodic microstructure which occupies the region ∈ R 2 in the undeformed, reference configuration, cf. Fig. 1. The material response at each macroscopic material point ∈ is represented by the response of a microstructural unit cell ∈ R 2 , for which it holds that | | ≪ | |, where | ⋅ | denotes the area. 1 This assumption permits us to treat the macroscopic body as an infinite periodic array of unit cells with identical microstructural properties. Without loss of generality, we consider rectangular unit cells with lattice vectors 1 = 1 1 and 2 = 2 2 , where 1 and 2 are the orthonormal basis vectors and ( 1 , 2 ) ∈ (R × R). The unit cell boundary with outward unit normal is decomposed into the disjoint sets + = + 1 ∪ + 2 and − = − 1 ∪ − 2 , containing material points on opposite boundaries, cf. Fig. 1. The macroscopic body is deformed into the current equilibrium configuration, ∈ R 2 , by a smooth mapping ∶ → , when loaded. The macroscopic deformation is locally described by the macroscopic deformation gradient, = , where is the gradient operator with respect to the macroscopic reference location . The macroscopic deformation drives the microscopic displacement field ∶ → R 2 , such that each microscopic material point in the unit cell identified by its position vector ∈ , is displaced to the position = ( , ) + ∈ , where and ( , ) is the heterogeneous displacement fluctuation field, cf. Saeb et al. (2016). From (1), we find that the microscopic deformation gradient is additively decomposed into a uniform part and a fluctuation part, i.e.
where is the gradient operator with respect to the microscopic reference coordinates .

Micro-macro relations
To maintain consistency between the micro-macro derivations we enforce the conservation of deformation, i.e. we require assuming the absence of voids. A common way to satisfy this condition is to restrict to the space of periodic functions. We take this approach. As such, the following kinematic constraint is enforced on the displacement field A. Dalklint et al. To remove rigid body modes, a single material point ∈ is restricted from motion, i.e. ( ) = , wherefore the kinematically admissible microscopic displacement field belongs to the space  = { ∈ 1 ∶ ( ) = , ( + ) = ( + − )+( − ) for + ∈ + , = 1, 2}. To maintain consistency between the virtual works in the micro-macro relations, we invoke the celebrated Hill-Mandel condition (Hill, 1972;Mandel, 1971), which stipulates that where = is the microscopic first Piola-Kirchhoff stress tensor, is the strain energy and = . By choosing = and = for the arbitrary constant 2nd order tensor , (5) reduces to the stress homogenization relation If we instead equate = in (5), we obtain which yields the quasi-static local microscopic equilibrium equation assuming no body forces. To obtain the above, we use the product rule, the divergence theorem and the anti-periodicity of the traction ⋅ , which follows since ( + ) = − ( − ) and ( + ) = ( − ). This implies ( + ) = ( − ) and ( + ) = ( − ).

Elastic wave propagation in a 2D periodic solid
To investigate the elastic wave propagation in an infinite prestressed 2D periodic solid, we follow Bertoldi and Boyce (2008b) and consider an infinitesimal time-dependent deformation̊( , ) superimposed on a given equilibrium configuration, defined by the displacement . The local microscopic equation of motion at time follows from (8) where is the mass density. The split results in Performing a Taylor series expansion of the left hand side in (11) and ignoring higher order terms provides where L = 2 | | | = is the fourth order incremental elasticity tensor. Combining (8) and (12), and inserting the result into (11), gives Using separation of variables, the solutions to (13) are constructed from the harmonic oscillations where ( ) are the mode shapes, = √ −1 and we assume non-attenuating waves, i.e. that the angular frequencies ∈ R. The assumption (14) allows us to pose (13) as which, since − is spatially independent and non-zero, reduces to The above is an eigenvalue problem which we solve for the eigenpairs (− 2 , ).
According to Bloch's theorem, we can solve (16) over a single unit cell, to obtain the response over the infinite periodic domain (Kittel, 1976;Joannopoulos et al., 2011). The solutions to (16) take on Floquet-Bloch form which can be seen as a plane wave ⋅ which propagates in the direction described by the wave vector ∈ R 2 , and is modulated by a spatially periodic function ( ).
Due to the microstructural periodicity, the material properties at any point ∈ must be invariant with respect to any translation along the lattice vector = 1 1 + 2 2 , for integers 1 , 2 ∈ Z. Accordingly, we enforce which, in combination with (17), leads to the classical Bloch condition which represents a kinematical constraint on . We note that for = , (19) constitutes regular periodic boundary conditions. A wave vector is expanded as where are the reciprocal lattice vectors which are defined such that In the reciprocal space spanned by 1 and 2 , we define the reciprocal lattice vectors = 1 1 + 2 2 , for integers 1 , 2 ∈ Z, so that ⋅ = 2 , ∈ Z. Because of this, we see that i.e. the modes are invariant with respect to the increment of . Therefore, we do not need to consider all wave vectors in the lattice, rather only those that lie in the so-called Brillouin zone BZ = { ∈ R 2 ∶ = 1 1 + 2 2 , − Fig. 2. Effectively, for each ∈ BZ, we solve (16) to evaluate the full spectrum of eigenpairs (− 2 , ). If there is an interval [ 1 , 2 ] for which no resides, then no waves with these frequencies can propagate the media, creating a ''bandgap''.
We assume that the unit cell exhibits orthotropic symmetry, so we can further limit the space of wave vectors to the irreducible Brillouin zone (IBZ). It is however crucial to remember that the applied deformation affects these symmetries, as depicted in Fig. 1. Indeed, the lattice vectors in the current configuration, , are related to their undeformed counterparts , as = , = 1, 2 (Zhang and Parnell, 2017). We limit ourselves to macroscopic prestrains which retain domain symmetry with respect to reflections about the 1 and 2 axes. As such our  Most researchers only consider the -vectors that lie on the boundary IBZ of the IBZ, i.e. they only interrogate a finite subset of IBZ, constructed as where , 1 ≤ ≤ are the wave vectors along IBZ and ∈ N is the number of sampling points. However, according to Maurin et al. (2018), the band extrema may not be located on the boundary, and therefore the full IBZ must be interrogated to confirm the omnidirectionality of a band-gap. As a compromise, we perform the optimization by considering only those ∈  , but confirm the findings by searching the entire IBZ in a post-processing step, i.e. for all ∈  = { ∈ IBZ, 1 ≤ ≤ } where ≫ . As seen in (17), the incremental displacement field is in general complex-valued. In our implementation, we follow Åberg and Gudmundson (1997) and decompose into its real and imaginary parts so that the eigenvalue problem (16) reads Separating the real and imaginary parts of (25) results in which are subject to the kinematic Bloch constraints (19) that couple the two real problems in (26). The weak forms corresponding to (26) are stated as finding the ℜ( ) ∈ ℜ() and

Equilibrium computation,
To determine the equilibrium configuration , we discretize (7) as where ∈ R is the residual vector, ∈ R is the global degrees-of-freedom vector and are the internal and external force vectors, respectively.
Following Ivarsson et al. (2020) and Wallin and Tortorelli (2020), we employ a standard master-slave elimination to impose (4), wherein the slave degrees-of-freedom are explicitly eliminated from the global degrees-of-freedom vector ∈ R by defining the mapping from the master degrees-of-freedom vectořcorresponding to nodes on the boundary − and in the interior ⧵ , to . The matrix = [ 11 , 12 , 21 , 22 ] has unit entries for the degrees-of-freedom on the boundary + where the macroscopic displacement To mimic uniaxial loading conditions, we prescribe a non-vanishing uniaxial displacement increment 11 ≠ 0, and determine 12 , 21 and 22 such that = 11 1 ⊗ 1 . Due to the assumption of orthotropy, we know that 12 = 0 and 21 = 0 for any biaxial stress state, cf. Wallin and Tortorelli (2020). As such, 22 is defined such that 22 = 0 2 . Equivalently, in our analysis we enforce To facilitate the solution process, we rearrange (30) as ] . Ultimately, we find * such that * for all * . Noting that = due to the PBC, and using the arbitrariness of * , the above reduces to * = = 0, where the partitioned residual * = [ , 22 ] is introduced. To find the displacement vector which fulfills (34), we employ Newton's method, wherein the linearization of (34) is required * ( * + * ) ≈ * + * where = * , * ∶= and = . To summarize, we prescribe the macroscopic deformation 11 through (32), and iterate by solving (35) for * and updating ← + * until convergence, i.e. until ‖ * ‖ 2 ≈ 0. In the equilibrium configuration, , we investigate the wave propagation properties of the material.
Lastly, we mention that biaxial loading conditions are obtained by prescribing the non-vanishing biaxial displacements 11 ≠ 0 and 22 ≠ 0.
The analysis proceeds by substituting (37) in (36) and using the arbitrariness of • , to obtain the unconstrained eigenvalue problem 2 Note that we require 22 = 0, i.e. we use the first Piola-Kirchhoff stress rather than the Cauchy stress . However, since 12 = 21 = 0 the 22 = 0 condition is equivalent to 22 = 0.
By utilizing the notation (⋅) for block quantities, (38) becomes where Following the usual convention, the eigenvalues of (39) are sorted in ascending order, i.e. 1 ≤ 2 ≤ ⋯ ≤ , and the eigenvectors are mass-orthonormalized such that The (24) decomposition of the eigenvectors provides a real generalized eigenvalue problem (39), which is readily solved using traditional solvers.

Design dependent material model
We assume isotropic hyperelastic material behavior, and we use a compressible neo-Hookean material model in which the strain energy is where and are design dependent bulk and shear modulii, and = det( ). In our topology optimization, we optimally distribute two material phases in the unit cell. The material distribution is described by the continuous design volume fraction field ∶ → [0, 1]. The elastic material parameters and , and the mass density , are subsequently interpolated as where the subscripts refer to the two material phases. For standard compliance minimization, the linear interpolation (43) does not promote binary designs. However, for phononic bandgap maximization, Sigmund and Jensen (2003) obtained nearly binary designs when using (43). Indeed, this follows the intuition that stark material contrasts are favored for large bandgaps, cf. Joannopoulos et al. (2011).
Traditional density based structural topology optimization is plagued by ''checkerboard-patterns'' and mesh-dependency, which stem from the inherent ill-posedness of the underlying optimization problem, cf. Sigmund and Petersson (1998). These issues are often mitigated through the use of filtration techniques (Bruns and Tortorelli, 2003;, by which a minimal length-scale is introduced in the design. Fortunately, in bandgap topology optimization, these issues are mute, as evidenced by numerical results (Sigmund and Jensen, 2003;Borel et al., 2004;Men et al., 2014;Swartz et al., 2021). In addition, Sigmund and Hougaard (2008) and Li et al. (2018) showed that optimized bandgap materials could be obtained from simple geometric considerations, which further justify this claim. Nonetheless, we follow Sigmund and Jensen (2003) and Swartz et al. (2021) and introduce a filter to limit the minimal feature size for fabrication considerations and to hasten the convergence of the optimization. However, unlike Sigmund and Jensen (2003) who adopted the sensitivity filter (cf. Sigmund, 1997), and Swartz et al. (2021) who adopted the cone filter (cf. Bruns and Tortorelli, 2001), we impose a length-scale restriction by penalizing the fine-scale oscillations of via the Helmholtz PDE-filter proposed in Lazarov and . In this way, in (43) is replaced by the smooth field ∶ → [0, 1], obtained from the solution of the boundary-value problem subject to PBC. In (44), is the Laplacian with respect to and the amount of smoothing of is controlled by ∈ R + . To solve (44), we parameterize via the piece-wise uniform design field ∶ → [0, 1] over the finite elements ∈ , such that ( ) = for ∈ . However, the same mesh and shape functions that are used to interpolate and are used to interpolate via the nodal filtered density degrees-of-freedom vector . The PBC are implemented by enforcing where and * are analogous to and * in (30). Following Lazarov and Sigmund (2011), the discretized version of (44) is where * = , * = and * = . We restrict ourselves to unit cells that exhibit orthotropic symmetry. As such, we reduce the design space by one quarter and write where ∈ R * contains the * = ∕4 design variables corresponding to the volume fractions of the elements in the reduced domain, and ∶ R * → R provides the 4-fold design symmetry illustrated in Fig. 3.

Optimization problem
We use topology optimization to design structures with tunable bandgaps. More specifically, we seek the design which maximizes the bandgap in the deformed configuration , whilst minimizing it in the undeformed configuration . To this end, we define the objective function is the gap/midgap ratio (Joannopoulos et al., 2011) between the eigenvalues and +1 , that we want to separate. By minimizing the cost function defined in (48), we promotē( , ) → 0 and̄( , ) → ∞.
Unfortunately, the objective defined in (48) leads to convergence issues since the max/min operators in (49) are nondifferentiable. Also, (48) is not differentiable when degenerate eigenvalues occur (see e.g. Seyranian et al., 1994;Dalklint et al., 2020). To mitigate both of these issues, we first follow Men et al. (2010), and separate the eigenpairs into two sets for each wave vector ∈  and displacement where min > 1 and max < dof truncate the full eigenspace. These sets let us introduce the differentiable conservative approximations As seen above, two -norms are used to approximate the max/min operations, such that ) .
This approach was also used by Swartz et al. (2021) and Quinteros et al. (2021).
To further improve the rate of convergence, and obtain a more versatile problem formulation, we next follow Qian and Sigmund (2011) and introduce the auxiliary design variables and append these to the vector of design variables, i.e. = , , , , ] , so that (49)   To accommodate the (53) design variables we enforce four constraints in our optimization problem cf. Fig. 4. Unfortunately, the block composition of (39) poses problems in the subsequent sensitivity analysis, as all eigenvalues appear at least twice. Therefore, each eigenvalue of (39) satisfies either; 1) the eigenvalue is twice-degenerate and originally simple, or 2) the eigenvalue is 2 -degenerate and originally -degenerate. For case 1), the two associated eigenvectors simply differ by the scalar factor , wherefore any of the two eigenvectors could be used in the subsequent sensitivity analysis. The second case is more involved, since it is not obvious which eigenvectors differ by the factor , and which span the -dimensional complex hyperplane. However, since we approximate the constraints using the symmetric polynomials (51), we overcome this issue by simply including all 2 -tuples of eigenvalues of interest in N .

Sensitivity analysis
The gradient based nonlinear programming method MMA, cf. Svanberg (1987), is used to solve (56). The gradients of the cost and constraint functions with respect to the design variables are derived below. The sensitivities with respect to the auxiliary variables , , and are trivially obtained. The symmetry (47) and filter (44) operations in combination with the chain rule, provide the sensitivities with respect to , i.e.
In the above, is obtained via the adjoint method described in Appendix A and * = , cf. (45). The product * * is obtained by a second adjoint sensitivity analysis described in Lazarov and . Finally, = is readily obtained from (47).

Numerical examples
The unit cell with side-length 1 = 2 = = 10 mm, is discretized using 2D bilinear quadrilateral plane strain fully integrated elements, unless stated otherwise. We assume an elastomeric matrix (material phase 1) and aluminum inclusions (material phase 2), described by material parameters stated in Table 1, cf. (43).
We solve the optimization problem (P) using the MMA-scheme with default parameters, cf. Svanberg (1987). Following Swartz et al. (2021), we solve (P) on a sequence of four increasingly refined finite element meshes to improve the convergence and reduce A. Dalklint et al.  (55) are inactive, i.e. such that 1 = 2 = 3 = 4 = 10 −4 . In this way, we ensured a feasible initial design, which improves convergence.
In all examples, we sample the boundary of the Brillouin zone, i.e. we solve (39), = 35 times, and equate min = 1, max = 12, cf. (50), and = = 16. Also, to limit convergence oscillations, we explicitly control the moving asymptotes, i.e. the MMA approximation, for the , , and bound variables, cf. Verbart et al. (2017). The bandgaps of the optimized designs have been verified by querying the entire IBZ, i.e.  on a 11 × 11 grid.

Uniaxial loading
In the first example, we impose uniaxial loading with 11 = 1 mm. We equate = 3 in the objective function (48), wherefore we aim to tune the 3 -4 bandgap. The initial design is random. The topology optimization problem is nonconvex and it is well known that the MMA finds relative minimum that are ''close'' to the initial design. In addition, the occurrence of mode-crossing might inhibit the optimization from generating designs with bandgaps. For these reasons, the problem is solved for 20 initial designs. In Fig. 5 we depict the three best performing designs, together with their corresponding dispersion diagrams.
To gain physical insight of the results presented in Fig. 5, in Fig. 6 we show the eigenmodes that bound the bandgaps of the three considered designs. We find that the modes 1 and 3 of the design reported in Fig. 5(a) as well as modes 5 and 7 of the design reported in Fig. 5(b) localize at the thin aluminum members. Upon loading, the thin aluminum members stiffen, causing an overall increase in the = 4 band frequency. Focusing on the design of Fig. 5(c), we note that the horizontal elastomeric inclusions embedded in the large aluminum plates reduce the body's shearing and bending resistance. This is reasonable since the upper bound of the bandgap in the undeformed configuration consists of shearing and bending modes (cf. 9 ), and the bandgap should be minimized. As before, the horizontal aluminum bars are subject to stiffening upon loading, wherefore the frequency is marginally increased (cf. 11 ).
Differently, the modes 2 , 4 , 6 , 8 , 10 and 12 involves oscillations of the large aluminum plates. Interestingly, upon application of uniaxial tension, the frequencies of these modes decrease. This phenomena is readily explained by examining Fig. 7, where the 1 and 2 modes of Fig. 5(b) design are shown. We note that in the undeformed configuration, the 1 mode localizes to the area in between the unit cells in the -direction, whereas in the deformed configuration, the 2 mode is evenly distributed over the entire aluminum body. This leads to an increase of the effective mass that is oscillating in the deformed configuration, and thereby a decrease of the frequency (cf. (39)). This redistribution of the mode is attributed to the stiffening effect which occurs when applying the load. In general, we conclude that there are mainly three effects that enables the tunability of the bandgaps, namely the changes in 1) the local material tangent, 2) the stress distribution and 3) the effective mass.
The Fig. 5(a) design outperforms the others, and will therefore be subject to further analysis. The evolution of Fig. 5(a) unit cell design and performance during the mesh refinement iterations is depicted in Fig. 8. Each mesh encompasses finer details than its predecessor, while preserving the main features. As expected, we also observe a downward shift of the eigenfrequencies when refining the mesh, with the greatest shift occurring during the first mesh refinement. Also as expected, the objective function worsens when activating the filter.
Since Fig. 8(d) design has a diffuse boundary caused by the interpolation and filtering, we conduct a post-processing analysis to further verify its performance. To this end, a Heaviside projection (cf. Guest et al. (2004), Kawamoto et al. (2011)) with a large threshold parameter is imposed on the filtered field to obtain a strictly binary design. The internal boundaries of the post-processed design are identified using the marching squares algorithm (cf. Lorensen and Cline (1987)), which facilitates the subsequent triangularization performed in GMSH (Geuzaine and Remacle, 2009). We emphasize that the interpretation of Fig. 5(a) design is non-unique.
The conforming mesh consists of fully integrated triangular elements. To ensure that the unit cell remains orthotropic, a quarter of the unit cell is discretized and subsequently reflected to obtain the unit cell mesh. The results from the post-processing analysis of Fig. 5(a) design are depicted in Fig. 9.
In Table 2, we denote the (48) objective function value and (49) gap/midgap ratios for the ersatz and conforming analyses. This performance loss shows the importance of accurate interface modeling. Indeed, the greatest gap/midgap ratios in the and configurations are computed over for the conforming mesh, which is expected as a phononic bandgap benefits from a stark A. Dalklint et al.  material contrast. Naturally, if the gap/midgap ratio in and both increase, it is not obvious that the value of our objective function should improve. We conclude this example by depicting the deformed/undeformed 3 × 3 unit cell arrays corresponding to the conforming mesh, cf. Fig. 10. We also illustrate the bandgaps between 3 and 4 as functions of the uniaxial stretch .
The next example considers the same problem but for higher frequencies, by equating = 5 (Figs. 11 and 12) and = 6 (Figs. 13 and 14), cf. (48). It was difficult to find bandgap designs when starting from a random initial design. Therefore, we used stiff cross and stiff circle initial designs to obtain Figs. 11 and 13 designs, respectively. In Figs. 12 and 14, we see the corresponding 3 × 3 arrays of unit cells and their bandgaps versus stretch plots.
Finally, we note that several of the identified optimal designs exhibit small structural features. To investigate the effect of such small features on the dispersion relation, we focus on the design presented in Fig. 11(b) and conduct two additional deformationdispersion analyses in which some of these features are manually removed. First, we remove the aluminum members highlighted by the red ellipses in Fig. 15. As shown in Fig. 16(a), we find that the band gaps are minimally affected by these modifications in the design. Second, we remove all features highlighted by the red and black ellipses in Fig. 15. The band diagrams shown in Fig. 16(b) indicate that this modification has a profound effect on the bandgaps. As such, we conclude that some of the small structural features play an important role, while others do not. This means that additional analyses are required before concluding A. Dalklint et al.  whether or not various features are necessary and whether or not the designs require strict manufacturing tolerances. A way of addressing the latter issue is to introduce a robust formulation in the optimization formulation, cf. e.g. Wang et al. (2011).

Biaxial loading
In this example, we consider the biaxial loading conditions with 11 = 22 = 1 mm. Using a random initial design and = 3 renders Fig. 17 design. Choosing = 6 and an initial stiff circle design produces Fig. 19 design. The results further confirm that our framework is able to produce designs with the desired dispersion-deformation behavior. When performing the post-processing analysis, we again note a decrease in performance, cf. Figs. 17(b) and 19(b). As before, we depict the corresponding post-processed 3 × 3 arrays of unit cells and their bandgap versus stretch plots in Figs. 18 and 20.

Conclusions
We have demonstrated the use of density based topology optimization to design tunable phononic crystals. Several designs, with tunable wave propagation properties are produced. To solve the optimization problems, we employed the gradient based MMA scheme. A mesh refinement procedure is used to reduce the computational cost of the design optimization, and a minimal length scale of the design is imposed via the PDE filter to facilitate their fabrication.  We investigate the effect of the diffuse interfaces, by comparing their performance to their corresponding binary design in a post-processing analysis. The comparison demonstrates that diffuse boundaries have great influence on the dispersive properties of the structures. To address this issue, a Heaviside projection continuation scheme could be included in our optimization framework.
However, this greatly increases the number of design iterations, and thereby the computational burden to an already costly procedure due to the great number of eigenvalue solves. In addition, many different initial designs must be used due to the nonconvexity of the topology optimization problem and the occurrence of mode-crossing. In the future, we will investigate means to reduce the cost of the eigenvalue solves.
We limit ourselves to uniaxial and biaxial extension. In states of compression, the effects of pattern transformations caused by instabilities (cf. Bertoldi

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A
The sensitivity of the objective function (54) is trivial = .
A. Dalklint et al.  The sensitivities of the 1 and 2 constraints functions, cf. (51), are obtained from the chain rule, i.e.