On the Emergence of Negative Effective Density and Modulus in 2-phase Phononic Crystals

On the Emergence of Negative Effective Density and Modulus in 2-phase Phononic Crystals

Amir Ashkan Mokhtari    Yan Lu    Ankit Srivastava asriva13@iit.edu Department of Mechanical, Materials, and Aerospace Engineering Illinois Institute of Technology, Chicago, IL, 60616 USA
August 2, 2019

In this paper we report metamaterial properties including negative and singular effective properties for what would traditionally be considered non locally resonant 2-phase phononic unit cells. The negative effective material properties reported here occur well below the homogenization limit and are, therefore, acceptable descriptions of overall behavior. The material property combinations which make this possible were first revealed by a novel level set based topology optimization process which we describe. The optimization process revealed that a 2-phase unit cell in which one of the phases is simultaneously lighter and stiffer than the other results in dynamic behavior which has all the attendant characteristics of a locally resonant composite including negative effective properties far below the homogenization limit. We investigate this further using the Craig-Bampton decomposition and clarify that these properties emerge through an interplay between the fundamental internal modeshape of the unit cell and a rigid body mode. Through explicit numerical calculations on 1-D, 2-phase unit cells, we show that negative effective properties only appear for the specific material property combination mentioned above. Furthermore, we provide a proof which supports this conclusion. The concept is also shown to hold for 2-D unit cells where we show that an appropriately designed hexagonal unit cell made of 2 material phases exhibits negative effective shear modulus and density in an appropriate frequency regime in which it also exhibits negative refraction. As in the 1-D case, we show that the homogenized description in the 2-D case is rigorous as well. An important conclusion of this paper is that the class of unit cells expected to result in negative properties can be expanded beyond the classic unit cell (three-phase unit cells with an explicit locally resonant phase) to include topologically simpler 2-phase unit cells as well.

Metamaterial, Phononics, Homogenization, Negative effective Density, Craig-Bampton, Topology Optimization
thanks: Corresponding Author

I Introduction

Metamaterials have been an active area of research over roughly the last two decades Srivastava (2015a). They are architectured materials which promise material properties which are generally impossible to find in naturally occurring materialsFang et al. (2006); Li and Chan (2004). As far as acoustic and elastodynamic metamaterials are concerned, the broad objective for their design and existence is to control the flow of sound or stress waves in manners which were not possible before their advent. Some of the wave control objectives include the design of cloaksMilton et al. (2006); Cummer and Schurig (2007); Norris (2008), negative refraction and superlensesKe et al. (2007); Morvan et al. (2010); Hu et al. (2004).

The predominant mechanism for the design and realization of acoustic and elastodynamic metamaterials in literature has been through the introduction of locally resonant phases in the unit cell of a phononic crystals. Under appropriate conditions, the internal resonance of these phases corresponds to singularities and otherwise anomalous behavior in the effective mass, effective density, and/or effective compliance tensors of these composites Milton and Willis (2007); Liu et al. (2000); Huang et al. (2009). These locally resonant crystals are predominantly composed of three phases which includes an explicitly locally resonant phase. In some studies, 2-D, 2-phase unit cells composed of an inclusion inside a host material have been noted to display resonance like behavior Wang et al. (2004) but they have never been associated with negative effective properties. In fact, the presence of resonance is a necessary but not sufficient condition for the emergence of negative effective properties. One of the simplest examples is the case of a 1-D phononic crystal with 2 phases. To our knowledge, this simple unit cell has never been shown to exhibit negative effective properties.

One of the main contributions of this paper is the elucidation of conditions under which a 1-D, 2-phase phononic crystal exhibits negative effective density and modulus. These conditions were computationally discovered by a multi-phase level set based topology optimization process capable of searching very large spaces in the design domain. Topology optimization has been used in the area of both phononics for the optimization of bandgapsLu et al. (2017); Bilal and Hussein (2011), and in metamaterials for the optimization of metamaterial behaviorLu et al. (2013); Krushynska et al. (2014); Yang et al. (2016). We have implemented a comparatively crude form of the level set method geared towards optimizing effective properties. At this point there are different methods available which can be used to calculate the effective properties in metamaterials. Readers are directed towards relevant review papers Srivastava (2015a); Sébastien Guenneau (2013) . Here we have implemented a field/ensemble averaging based homogenization scheme for the purpose of optimization Nemat-Nasser and Srivastava (2011a, b); Nemat-Nasser et al. (2011); Srivastava and Nemat-Nasser (2012, 2014); Willis (1983, 2009, 2011); Sridhar et al. (2018); Muhlestein and Haberman (2016). The optimization method as implemented by us is crude but what the method lacks in sophistication, it makes up in decreasing the size of optimization problem and allowing us to deploy global search schemes such as the basin-hopping technique Wales and Doye (1997). Using off the shelf codes and algorithms Jones et al. (2014), our optimization algorithm results in interesting optimized unit cells. This includes both locally resonant multi-phase unit cells as well as non-locally resonant 2-phase unit cells exhibiting metamaterial properties.

We further explore the dynamics of the 2-phase unit cell suggested by the optimization algorithm through the use of the Craig-Bampton decomposition. This method was introduced to express the dynamic behavior of a structure as a function of its (or its substructures’) internal vibrational modes Craig and Bampton (1968). In the phononics area, it has been used recently to develop highly efficient algorithms for bandstructure and associated calculations Krattiger and Hussein (2014). In the metamaterials area, it has found natural usage in connecting the effective properties to the internal resonance modes of unit cells Sridhar et al. (2018). In the present application, the Craig-Brampton decomposition gives invaluable insights into the nature of the interplay between the internal vibration modes of a unit cell and its rigid body motion. For example, the method illuminates that on the second passband, the behavior of the unit cell is dominated by the first vibrational modeshape of the unit cell under either a fixed-fixed or guided boundary conditions – a conclusion which has been noted earlier using analytical techniquesMead (1975). We use the decomposition to highlight the close correspondence between the dynamics of a locally resonant unit cell and a 2-phase unit cell made of an appropriate combination of material properties (suggested by the optimization scheme.) The insights gained from the Craig-Bampton analysis allow us to formally consider the problem of the emergence of negative effective properties in 2-phase phononic crystals. Through explicit numerical simulation and analytical proof we show that appropriately designed 2-phase phononic crystals can exhibit negative effective properties far below the homogenization limit, thus constituting a valid description of overall behavior much like the negative effective properties for locally resonant unit cells.

In the subsequent sections we first introduce the effective property calculation scheme based on ensemble average of field variables along with two examples, one of which possesses negative effective properties and the other shows only positive values. We then define the level set topology optimization framework in Section III and present the multiphase optimization results as well as a result which shows a two phase unit cell exhibiting negative effective density and elastic modulus. Craig-Bampton decomposition is presented in Section IV followed by its application to the analysis of unit cell dynamics, revealing the central mechanisms which gives rise to negative properties – an interplay between the first internal vibration mode and rigid body motion. Finally, we apply the essential findings in layered composites to the design of two phase, two dimensional hexagonal unit cell in Section V, where we show negative effective properties on the optical branch as well as a simulation of negative refraction in the associated frequency region.

Ii Effective Properties

We consider harmonic waves traveling in a layered composite with a periodic unit cell . Waves are assumed to be normal to the interfaces and in this case the field variables (displacement, velocity, strain, stress, and momentum) are scalar and take the following Bloch form:


where functions , , , , , and are periodic with the periodicity of the unit cell. The dynamic equilibrium and the strain-rate/velocity relations give


One way of defining the effective properties is through the technique of field averaging Nemat-Nasser et al. (2011) which is a subset of ensemble averaging Srivastava (2015a) for periodic composites. Specifically we define:


It can be shown that the above effective properties satisfy both the average equations of motion and the dispersion relation of the composite:


The effective parameters defined by (3) are real-valued only if the unit cell is symmetric Nemat-Nasser and Srivastava (2011a); Amirkhizi (2017); Srivastava and Nemat-Nasser (2012); Willis (2012); Srivastava (2015a). For non-symmetric unit cells, the coupling among the field variables renders the effective stiffness, , and mass-density, , complex-valued. Two simple 1-D examples are given in Fig. (1). Fig. (1a) shows a locally-resonant 3-phase unit cell which is composed of a stiff and heavy central core coated with a compliant phase, all enclosed in a stiff matrix. The sub-figures below the unit cell show the bandstructure over the first two bands, the corresponding frequency dependent effective density () and effective modulus () respectively. The effective properties have been calculated according to Eq. (3). These figures show that there exists a frequency region on the second band of the band-structure where the effective properties become negative. This behavior does not occur in the case of a regular non-resonant 2-phase unit cell (Fig. 1b) where the effective properties are never less than zero in any region over the passbands. The effective properties are always real in the passbands, however, they are complex in the stopbands. Only the real parts are plotted in Fig. (1).

Figure 1: Dispersion relations and effective properties of 1-D layered structures. Material properties are excerpted from RefNemat-Nasser and Srivastava (2011b). Total thickness of each material phase is listed. a. Locally resonant structure. b. conventional 2-phase structure. Zero effective properties are marked by the red dot-dashed lines.

Our main motivation in this paper is to elucidate conditions under which 2-phased unit cells can exhibit negative properties. The system which makes this possible was discovered by a topology optimization framework which we describe in the next section.

Iii Topology Optimization Framework

Topology optimization has evolved rapidly in recent years as a form-finding methodology for structural and materials designDeaton and Grandhi (2014); Sigmund and Maute (2013); Cadman et al. (2013); Osanov and Guest (2016); Lu et al. (2017). It seeks to optimize the distribution of material resources across a design domain such that a defined objective function is minimized (or maximized) while the constraints are satisfied. Typically, finite element methods are used to discretize the design domain and a material relative density ranging continuously from to , is assigned to each element, with and indicating the presence of only material 1 or material 2 in the element, respectively. Intermediate values represent mixtures of the two material phases and are prevented by penalizing their existence, such as through the Solid Isotropic Penalization MethodBendsøe (1989a); Rozvany and Zhou (1991).

Another category of approaches for topology optimization is based on the level set method. Level set methods were introduced by Osher and Sethian Osher and Sethian (1988) and they express material boundaries as zero level sets of a function in a higher dimension Allaire et al. (2002, 2004); Osher and Santosa (2001). With level set methods, there is no need to explicitly parameterize the moving surfaces and boundaries. For example, considering the topology optimization of a 2-D 1-phase (solid-void) structure for minimum compliance, the level set method proceeds by defining a function whose level set defines the boundary between the sold-void phases. The function is changed through optimization techniques which implicitly changes the boundaries of the solid-void phase. Here, we use a version of the level set method which is primarily geared towards speed and a potential global search. This is possible to do in our case since we are dealing only with a 1-D case. A description of the technique is presented as follows.

Figure 2: (a) Level set functions (b) material distribution corresponding to the level set functions

Consider the topology optimization problem of distributing three different material phases in a 1-D unit cell for achieving some metamaterial objective. Specifically, we want to find a material distribution in this unit cell that leads to negative effective density at a desired point on the irreducible Brillouin zone. The material distribution is governed by the zeros of two polynomials and which are defined over the axis using the locations of their roots. These two polynomials divide the domain into 3 sub-domains: , , and (Fig 2). These three sub-domains are identified as the three different phases of the materials in the design domain. In this technique, the boundaries of these three phases are crisply defined as opposed to the SIMP method of topology optimization where the boundaries are blurry Bendsøe (1989b). The material distribution is changed by changing the root locations of and . In our implementation of the level set method, we have only a few design variables (root locations) for the optimization problem and we can, therefore, use global optimization methods to search very large design-domains. For the optimization algorithm, we use a stochastic technique called Basin-hopping Wales and Doye (1997) which attempts to find the global minimum of an objective function.

Figure 3: (a) Initial random material distribution (b) density distribution in the unit cell (c) effective Density(d) Young modulus distribution

As the first example, the three phases being optimized have material properties which correspond to a traditional locally resonant unit cell. One of the phases is stiff and dense, another phase is compliant and light, and the third phase has stiffness and density values between . For this problem, the objective function to be minimized is the effective density at a desired wavenumber and passband (3). Fig. (3-a) shows the random density distribution (under a relaxation function Belytschko et al. ()) over the unit cell taken as the initial condition for the optimization algorithm. Fig. (3b,d) are the material property distributions after optimization showing the clean separation of the three phases. When periodicity is considered, this configuration is a traditional local resonance type unit cell with a hard inclusion in the middle (), a soft coating layer around it (), and all surrounded by a matrix phase (). Compared with Fig.( 1-a), Fig. (3-c) shows the effective density calculated over the first two passbands. From Fig. (3-c), it can be seen that the effective density goes to a large negative value at a frequency value of around 1600 kHz. The level set optimization algorithm, therefore, naturally finds the characteristic metamaterial unit cell configuration when the optimization objective is the minimization of the effective density. The optimized structure exhibits a negative density region which is expected from an appropriate locally-resonant phononic crystal.

Figure 4: (a) Stiffness distribution in the unit cell (b) density distribution in the unit cell (c) band-structure (d) effective density

The above results are not surprising as the topology optimization framework merely reiterates the unit cell configuration which is traditionally understood to produce negative effective properties in metamaterials. As the next step, in addition to the root locations of the level set polynomials, we also took the material properties (density and Young modulus) in those phases as design variables in the optimization scheme. The optimization objective was to find the material property and geometric distribution combination which would minimize the effective density at a desired wavenumber in a passband. Although a three-phase structure was anticipated as the optimized result in this case, the optimizer, starting with a random distribution of material properties, instead produced a 2-phase structure which exhibited both negative and singular effective density values in some frequency ranges (Fig. 4). The material combination for this structure was different from a conventional 2-phase structure in the following way: in general, there exists a positive correlation between modulus and density of materials Ashby (2005). In the optimized 2-phased structure produced by the optimizer, one of the material phases was both lighter and stiffer than the other phase ( and ). The material distribution and effective density plots of the optimized structure are shown in Fig. 4. Figs. (4a,b) show the stiffness and density distributions in the optimized unit cell, showing that the optimized structure is a 2-phased composite and that it is made of the specific combination of materials that is mentioned above. Figs. (4c,d) show the bandstructure of this unit cell and the effective density, respectively, underlying the metamaterial nature of this unit cell with the negative and singular values of effective density in some frequency ranges. In the next section, we relate the emergence of negative and singular effective properties to the properties of the internal resonance characteristics of the 1-D unit cell using Craig Bampton reduction method.

Iv Craig-Bampton Analysis

Craig-Bampton decomposition Craig and Bampton (1968) is a technique which was first introduced in structural dynamics for the analysis of complex structures in terms of their sub-components. In this technique, the structure is typically divided in two regions - its boundary and its interior (called the substructure.) The dynamic behavior of the structure is described in terms of a set of generalized coordinates which include constraint coordinates and normal mode coordinates. The constraint modes relate the boundary nodal displacement of the structure to the interior displacements whereas the normal mode coordinates are the vibration modes of the interior with all boundary nodes fixed. This approach, when applied to the present problem, can elucidate the effect of internal resonance on the overall Bloch wave. We summarize the relevant equations for the 1D case below (See Bloch Mode Synthesis Krattiger and Hussein (2014) for more details). Consider the discretized form of the dynamics of a unit cell in a 1D periodic structure:


Where and are the nodal stiffness and mass matrices respectively and , and represent the nodal displacements, acceleration and applied forces respectively. The displacement vector can be divided into a set of interior nodes () and left and right boundary nodes (.) The partitioned form of the Eq 5 is:


Before applying the Bloch boundary condition on the end nodes and , Craig-Bampton model reduction is used to express the dynamics of the interior nodes in terms of the normal modes and constraint modes as described above. Using the constraint and normal modes, one can replace with a new vector through the following transformation


where is a matrix whose columns are the vibrational modes of the interior of the unit cell when the boundary nodes are fixed, is a column vector consisting of the amplitude of each normal mode, and is a matrix describing the constraint modes. Since are vibrational modes corresponding to Dirichlet boundary conditions on the unit cell boundary, we have


The governing equation (Eq. 5) can now be transformed into the new degrees of freedom which involve the new C-B transformed stiffness and mass matrices:


At this point in the current problem, the Bloch boundary condition is applied on the left and right boundary nodes and by the following transformation:


Where is the wavenumber and is the unit cell length. This further transforms the stiffness and mass matrices into:


The final form of the eigenvalue problem is as follows:


The eigenvalues resulting from the solution of the eigenvalue problem give rise to an ordered sequence of radian frequencies with corresponding eigenvectors . In any eigenvector , the magnitude represents the relative contribution of the vibrational mode in the eigenvalue-eigenvector pair. With this, the periodic part of the displacement for a given value and over the branch, , can be expressed as where


where , , and are concatenated to form the total mode shape and is the number of normal vibration modes considered in Eq. (7). is the contribution coming purely from the internal modes of vibration and is the contribution coming from constraint mode. As it can be seen from Eq. (16), is a function of unit cell free vibration normal modes and constraint modes with Bloch boundary condition. When defining effective properties from the solutions above, it is critical that the effective properties, defined in terms of unit cell averages, make physical sense. This essentially translates into the requirement that we should be able to identify a matrix phase which exhibits minimal deformation at the frequency where we are trying to define the effective properties. If this frequency is then this roughly translates into the requirement that should be much larger than the size of the matrix phase. In addition, it is required that the wavelength inside each inclusion should be larger than or on the order of the size of the inclusion Pham et al. (2013).

Locally resonant unit cells

Consider a 3-phase locally resonant unit cell in a symmetric configuration (Fig.1a.) The matrix, inclusion, and the inclusion-coating are indicated by the subscripts respectively. If the length of each phase is kept constant, we can follow a concrete set of strategies to induce negative effective properties on the second branch of the band-structure. It has already been shown that for 1-D symmetric unit cells, the frequency limits of the passbands are governed by the natural vibration frequencies of the unit cell in the fixed-fixed and free-free configurations Mead (1975). As we will show below, the second branch of locally resonant phononic crystals is primarily influenced by the first vibrational mode of its unit cell in the fixed-fixed configuration which also roughly governs the frequency range of the second branch. If the frequency of this mode is then . Generally would change with the unit cell configuration and material properties but it is possible to make it largely independent of large changes in certain properties of the unit cell. Consider the case when the matrix and the inclusion are much stiffer than the coating. In this case, the modeshape of the fundamental vibrational mode in the fixed-fixed configuration of the unit cell displays a behavior where almost all the deformation is in the coating layer. The inclusion moves almost rigidly and there is negligible deformation in the matrix due to its ends being fixed. If we now fix the material and geometric properties of the inclusion and the coating, it is clear that will be largely independent of the matrix properties due to the fact that there is almost no deformation in the matrix. In fact, will be roughly given by . The wavelength in the matrix at this frequency is and as long as , the unit cell can be homogenized. Ignoring the term, this simplifies into the relation for homogenization to be applicable where and . If the lengths of the material phases are such that then this further simplifies into or .

Figure 5: (a) Scaled modeshape , and such that , (b) modeshapes plotted for the different values of (c) and for

These ideas can be illustrated by taking an example. Consider the set of cases where and . Fig. (5-a) shows the displacement profile, , of the unit cell for and on the second branch when a specific combination of material properties is assumed (). Note that for this case so that homogenization can be performed on the second branch. It also shows the contributions from the internal vibration modes and the boundary constraint mode , and the first internal vibration mode of the unit cell . There are two points of interest here. First, is almost entirely made up of the first internal mode of vibration of the unit cell or , and the boundary constraint mode is almost a rigid body mode or . Since velocity is directly proportional to displacement, we have near where is a scalar and is the velocity profile corresponding to . The momentum profile is . Note that in Fig. (5-a), is a negative number which elucidates the central mechanism of locally resonant unit cells – inclusion moving out of phase with the matrix. If is the velocity of the inclusion in the modeshape then the average of the velocity, , which appears in effective density equation can be approximated as . This is possible since the matrix is nearly stationary and the coating layer is thin. Without any loss of generalization this can be written as where and since can always be normalized to set . Similarly, the average momentum is where is a positive real number. and will have opposite signs if . The requirement for this equation to hold is (if which is true in the present case.) is inversely proportional to which means that decreases (increases) with increasing (decreasing) . Changing has no appreciable effect over since the matrix is nearly stationary anyway. However, changing serves to change . This is illustrated in Fig. (5-b) where is plotted for several values of ( are kept constant.) Unsurprisingly, increasing (decreasing) serves to decrease (increase) the influence of the rigid body component of motion , in turn, decreasing (increasing) in the process. Fig. (5-c) shows as functions of indicating the condition derived above for negative effective properties is indeed satisfied for a range of . In summary, as long as , homogenization on the second branch is likely acceptable. Within this constraint the density of the matrix phase can be tuned to modulate the balance between the rigid body component and the first internal mode of vibration component of motion thus giving rise to negative effective density (at least near .) Since the effective properties satisfy the dispersion relation of the composite (Eq.4), it is automatically assured that in such regions the effective stiffness will also be negative. Of course, in the above, it is not necessary to assume that the coating layer is thin, however, thin coating layer serves to render the equations in a simple form.

Two phase unit cells

Figure 6: (a) Effective stiffness, (b) effective density (c) mode shape , normal mode shape and constraint mode for the unit cell made of materials with , , and
(d) effective stiffness, (e) effective density (f) mode shape , normal mode shape and constraint mode for the unit cell made of , , and

Consider a 2-phase unit cell in the symmetric configuration like Fig. (1). In this case, we cannot separate the system into matrix, inclusion and coating layer like the 3-phase case. Instead it has two parts: the soft (with subscript ) and hard phase (with subscript ). Like the 3-phase case, for a 2-phase unit cell with periodic boundary condition the behavior on the second passband is largely influenced by the first natural vibration modeshape of the unit cell with the fixed-fixed or guided configuration. If the stiffness of one of the phases is large enough compared to the other phase () then the first free vibration modeshape of the unit cell will comprise of deformation which is concentrated in the soft layer. Again, if the stiffness contrast is high enough, the first natural frequency is determined by the properties of the soft layer which is roughly given as . This is the result of an eigenvalue problem emerging from a second order homogeneous differential equation in the soft layer with Dirichlet boundary conditions set to zero displacement. The wavelength in the hard layer at this frequency is and the relaxed requirement for homogenization is that the wavelength in the matrix should be much larger than the length of the matrix but larger or comparable to the length of the soft phase (inclusion phase.) Note that under the relaxed homogenization requirement, this case is indeed homogenizable. In the vicinity of , the wavelength in the hard phase is much larger than the length of the phase and, given the modeshape of the first vibrational mode, it is almost twice the length of the soft phase. The nature of the deformation also allows us to consider the harder layer as the matrix. Defining and , homogenization is allowed if or equivalently .

Consider a specific 2-phase unit cell with made of two materials (equal lengths ) with , , and , which the second material is heavier and softer than the first one. For the wavenumbers near , the normal modeshape is almost entirely made up of the first mode of the unit cell natural vibration () and the harder layer is forced into a rigid body motion () as depicted in the Fig. (6-c). Like the 3-phase case, we have where is a scalar and is the velocity profile corresponding to . The momentum profile is as well. If is the velocity of the soft layer in the modeshape then the average of the velocity, , which appears in effective density equation can be approximated as where . Similarly the average momentum is where . Following the results from the previous section, effective density would be negative if (if ) or (if .) Of course, since this is a 2-phase system, the requirement for negative properties is if and if . The effective property and unit cell deformation results for the present case are graphically shown in Figs. (6a-c) showing the regions of negative effective properties on the second branch. Figs. (6d-f) show the corresponding results for another case where the only difference is in the density of the hard phase. Specifically, , , and . Comparing the mode shapes for both cases (Figs. 6c,f), we note that most of the displacement is in the softer layer in both cases and the harder layer is in a quasi-static state. However, in the first case – in which the softer phase has also a larger density – the mode shape is significantly divided into two parts compared to the second case. This is similar to the behavior of a relatively heavier mass vibrating inside a soft material which is characteristic of 3-phase locally resonant unit cells. In contrast, for the second case, in the softer layer most the displacement is in the region. In effect, modulating serves to modulate thus giving rise to the potential of negative effective properties.

Figure 7: (a) Effective density and (b) effective modulus distribution within the frequency range of the second dispersion curve of 2-phase unit cells. Arbitrary units are applied. (c) Effective properties corresponding to the points marked in (a).

Fig. (7) shows the effective property distribution () as a function of for a 2-phase unit cell equal material phase lengths. The material phases have properties respectively and have been defined as and . Figs. (7a,b) are log-log plots where the horizontal axis is and the vertical axis is . is the main diagonal. Points correspond to cases where and are, thus, situated in homogenizable regions. Points correspond to cases where which are also homogenizable when one interchanges the two material phases – something which is possible in a two phase layered composite. In Figs. (7a,b), effective properties are evaluated for the second dispersion bands and we only collect the positive and negative values with the smallest absolute values among the discretized frequencies over which the effective properties are calculated. The values are appropriately normalized so that the boundaries can be clearly shown. First it should be noted that Figs. (7a,b) seem to indicate that negative effective properties are achieved only in 2 of the 4 quadrants. The top right quadrant corresponds to and the bottom left quadrant corresponds to . In both these cases, one can conclude that negative effective properties are only achieved when one of the phases is simultaneously stiffer and lighter than the other phase. This is explicitly shown for the marked points in Fig. (7a) by calculating the effective density and modulus over the first two branches. These points marked are paired and separated by the horizontal and vertical boundaries of the quadrants and they are chosen to demonstrate that by only changing one of the two material property ratios from to , or vice versa, the effective properties can be made to flip signs. These points are also chosen to satisfy or so that homogenization of the second branch may be acceptable. The effective property results, shown in Fig. (7c), confirm that negative effective properties emerge on the second branch for points but do not for points .

Required Material Properties for Negative Effective Properties

Figure 8: (a) Unit cell with material properties: , , and (b) Modeshapes calculated using two different methods: Craig-Bampton and Eq. (18) (c) Bandstructure of the unit cell (d) plot and the corresponding location of the points A,B and C.

The result that negative effective properties only emerge for certain combination of can be shown analytically. Assume that the lengths of the phases is the same (Fig. 8a). Also assume that and . At and on the second branch, we have shown that the total modeshape is a combination of a rigid body mode (because ) and the first vibrational frequency mode (). Let’s assume that the displacement profile in just the soft layer is . The form of the modeshape suggest that can be determined by considering the governing equation in the soft layer:

with the boundary condition:


which imposes Dirichlet boundary conditions corresponding to a rigid body motion to the two ends of the soft layer. In the above, is the wave velocity in the soft layer. The solution to this equation is:


The full displacement field in the unit cell which emerges from the above is plotted in Fig. (8b) showing that it compares very well with the Craig-Bampton solution ( term has been ignored in the plot). Having the displacement fields in both of the soft and hard layers, we can find the expression for the effective density as follows:


where and which is the wavenumber in the soft layer. Fig. (8c) shows the bandstructure of a typical unit cell. corresponds to point B in the figure and we intend to show that in the vicinity of B, , given in the equation above, is negative if is greater than 1 and sufficiently large. To show this, we will assume that and show that in the vicinity of point B, the denominator is positive and . This would then automatically mean that will be negative if is sufficiently large (from Eq. 19).

The first resonance frequency of the unit cell with the ends fixed coincides with the endpoint of the second pass band (point A in the Fig. 8c) and the second resonance frequency coincides with the beginning of the third passband (point C in Fig. 8-c)Mead (1975). If we assume that the deformation in the hard layer is negligible, then the frequency at the point A is related to the wavenumber in the soft layer. At this wavenumber the term . This point is shown in the Fig. (8-c). Starting from point A on the second passband and moving toward the point C on the plot, we are moving from (for ) to zero (for ). Thus, the tangent term in Eq. (19)is negative on the second passband.

Now for point B where , we have to prove that . We can write:


The mark means that the proof is required for the statement. As explained above, point B is in the range where , so . Thus, we have:


Now we consider the RytovRytov (1956) dispersion solution:


Where and . The RHS is equal to 1 in the above at point B. If the hard layer is behaving quasistatically then we can show that at point B and . Since the RHS is equal to 1 in the above, we have to show that:


to complete the proof. From the above we must show:


The last part of Eq.(23) is true given the assumption that . We have, therefore, shown that if then the denominator of Eq. (19) is positive whereas the tangent term is negative. By making sufficiently large, we can make the numerator negative whereas the denominator is assured to be positive thus giving rise to negative effective density. In fact is the range in which the denominator of the Eq. (19) is positive because . Point B is in this range for the considered combination of and as proved. This mathematical approach elucidates how the material combination of and leads to negative effective properties.

V Negative effective properties in 2-D 2-phase unit cells

The essential lessons from the last section also carry on for the design of 2-phase, 2-D unit cells which exhibit negative effective properties. Determination of effective properties in higher dimensions is a much more complicated task than it is in 1-D (Srivastava and Nemat-Nasser, 2012). Therefore, we choose an example which is particularly simple but which serves to elucidate the concept. We consider the propagation of antiplane shear waves in a 2-D hexagonal lattice. The displacement field variable of interest is where is periodic with the unit cell. The displacement field gives rise to stresses with unit cell periodic parts respectively. In addition we have velocity (), momentum (), and strain fields () with periodic parts respectively. The equations of motion are:


We can define effective properties through field averaging:


It can be shown that the above effective properties satisfy the average equations of motion:

In general, and will not be the same. Such anisotropy leads to complications as far as determining when a unit cell can be said to exhibit negative effective properties is concerned Srivastava (2015b). However, if we can find cases where wave propagation is largely isotropic in a certain frequency range, then in that range we can expect . In such cases if then the composite can unambiguously be said to exhibit negative properties. In such cases, we can also show that the effective properties satisfy the dispersion relation:

Figure 9: (a) Hexagonal unit cell, (b) Band structure and effective properties, (c) Modeshape on the optical branch at 2548, (d) Equi-frequency contour, (e) Negative refraction at 2548kHz

For the hexagonal unit cell under consideration, the lattice includes a matrix phase with shear modulus and density . The unit cell contains a symmetric circular inclusion made of a material with shear modulus and density (Fig. 9a). The matrix is chosen to be stiff which ensures that it exhibits minimal deformation up to and around the first vibrational frequency of the circular inclusion. Conversely, the circular inclusion is made of a compliant material which ensures that it deforms far easier than the matrix. Furthermore, its density is chosen to be suitably high in order to ensure that is small and, thus, ensuring the applicability of homogenization around . This choice of material property directly follows from the last part of the last section where it was found that material inclusion phases which are simultaneously more compliant and heavier than matrix phases preferentially lead to negative effective properties. The unit cell for a specific material combination satisfying these properties () is shown in Fig. (9)a with lattice constant and filling fraction 0.5. In Fig. (9)d, we show the equi-frequency contours (EFCs) for both the acoustic and the optical branches for the unit cell. In the contour maps, circular contours represent frequency zones where wave propagation is isotropic. Unsurprisingly, there are circular contours on the acoustic branch signifying regions of quasi-static behavior. More interestingly though, there are circular contours on the optical branch as well. These appear at higher frequencies () on the optical branch and signify regions where wave propagation is isotropic. Fig. (9b), shows the effective properties calculated along the direction and over the first two branches. It should be noted that the effective properties are negative on the second branch. More importantly though, they are negative in the region around and in that region is ostensibly equal to . Similar to the 1-D cases in the previous section, in these frequency regions deformation occurs predominantly inside the circular inclusion while the matrix remains in a quasistatic state. This is shown, for example in Fig. (9c), where the modeshape of the unit cell at 2548 has been plotted. The circular equi-frequency contours on the second branch also predicts negative refraction. Negative refraction is explicitly shown by a frequency domain Finite Element simulation in Fig. (9e). Here the incident wave is excited by the lower left line source at 2548 and it propagates predominantly perpendicular to the shorter edge of the prism. The prism is made of repeating unit cells of the hexagonal crystal described above. The homogeneous media in which the prism is embedded is the same as the matrix material of the phononic crystal. The corresponding wave number in the homogeneous medium is 101.253rad/m and in the crystal is 67.992rad/m. When the wave reaches the upper interface between the prism and the homogeneous media at 60, the outward negative refraction angle is predicted to be 35.56 from Snell’s law. This matches the results from the FE simulation.

Vi Conclusions

A novel level set topology optimization method has revealed that a 2-phase phononic crystal unit cell in which one of the material phases is simultaneously lighter and stiffer than the other can give rise to negative and singular effective properties on the second branch. The proposed topology optimization method seeks to optimize the material distribution governed by the zeros of level set polynomials. Although in a comparatively crude form, this method allows us to search very large design spaces using global optimizers such as Basin-hopping. Craig-Bampton decomposition is used to gain insight into the emergence of negative effective properties in 2-phase composites. The decomposition elucidates that the second passband is dominated by the fundamental mode of vibration of the soft layer under Dirichlet boundary conditions and it further elucidates the central mechanism of local resonance in layered composites – the interplay between the deformation in the soft phase with the rigid body motion of the other phase. The method also shows that by tuning the material combinations, the balance between the fundamental mode of vibration and the rigid body component is modulated, thus giving rise to the possibility of negative effective properties. When considering 2-phase unit cells with equal layer thicknesses, we can understand the effect of material contrast ratios on effective properties through explicit calculations. These calculations show that simultaneous negative effective density and stiffness happens when material combination appears in two of the four quadrants where one material phase is softer and lighter than the other phase. By slightly changing the contrast ratio from one quadrant to another the effective properties can be made to flip signs. We have provided a mathematical proof which also supports this conclusion. The insights gained from the 1-D case also applies to the 2-D case. We design a 2-D, 2-phase hexagonal unit cell with simple circular inclusions and appropriate material properties and show that it results in negative effective density and shear modulus within an appropriate frequency range. In this range, the wave propagation is nearly isotropic and the matrix behaves quasistatically thus homogenization is acceptable. As in the 1-D case, in this case as well, the deformation occurs predominantly within the soft inclusion and a wave incident at an interface between the crystal and a homogeneous medium suffers negative refraction.

With the results presented in this paper, we can expand the class of phononic unit cells which exhibit negative effective properties to include simple 2-phase unit cells as well.

A.S. acknowledges support from the NSF CAREER grant #1554033 to the Illinois Institute of Technology and NSF grant #1825354 to the Illinois Institute of Technology


  • Srivastava (2015a) Ankit Srivastava, “Elastic metamaterials and dynamic homogenization: a review,” International Journal of Smart and Nano Materials 6, 41–60 (2015a).
  • Fang et al. (2006) Nicholas Fang, Dongjuan Xi, Jianyi Xu, Muralidhar Ambati, Werayut Srituravanich, Cheng Sun,  and Xiang Zhang, “Ultrasonic metamaterials with negative modulus,” Nature Materials 5, 452 EP – (2006).
  • Li and Chan (2004) Jensen Li and C. T. Chan, “Double-negative acoustic metamaterial,” Phys. Rev. E 70, 055602 (2004).
  • Milton et al. (2006) Graeme W Milton, Marc Briane,  and John R Willis, “On cloaking for elasticity and physical equations with a transformation invariant form,” New Journal of Physics 8, 248 (2006).
  • Cummer and Schurig (2007) Steven A Cummer and David Schurig, “One path to acoustic cloaking,” New Journal of Physics 9, 45 (2007).
  • Norris (2008) Andrew N Norris, “Acoustic cloaking theory,” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 464, 2411–2434 (2008)http://rspa.royalsocietypublishing.org/content/464/2097/2411.full.pdf .
  • Ke et al. (2007) Manzhu Ke, Zhengyou Liu, Zhigang Cheng, Jing Li, Pai Peng,  and Jing Shi, “Flat superlens by using negative refraction in two-dimensional phononic crystals,” Solid State Communications 142, 177–180 (2007).
  • Morvan et al. (2010) Bruno Morvan, Alain Tinel, Anne-Christine Hladky-Hennion, Jérôme Vasseur,  and Bertrand Dubus, “Experimental demonstration of the negative refraction of a transverse elastic wave in a two-dimensional solid phononic crystal,” Applied Physics Letters 96, 101905 (2010).
  • Hu et al. (2004) Xinhua Hu, Yifeng Shen, Xiaohan Liu, Rongtang Fu,  and Jian Zi, “Superlensing effect in liquid surface waves,” Physical Review E 69 (2004), 10.1103/physreve.69.030201.
  • Milton and Willis (2007) GW Milton and JR Willis, “On modifications of newton’s second law and linear continuum elastodynamics,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 463, 855 (2007).
  • Liu et al. (2000) Z. Liu, X. Zhang, Y. Mao, YY Zhu, Z. Yang, CT Chan,  and P. Sheng, ‘‘Locally resonant sonic materials,” Science 289, 1734 (2000).
  • Huang et al. (2009) HH Huang, CT Sun,  and GL Huang, “On the negative effective mass density in acoustic metamaterials,” International Journal of Engineering Science 47, 610–617 (2009).
  • Wang et al. (2004) Gang Wang, Xisen Wen, Jihong Wen, Lihui Shao,  and Yaozong Liu, “Two-dimensional locally resonant phononic crystals with binary structures,” Phys. Rev. Lett. 93, 154302 (2004).
  • Lu et al. (2017) Yan Lu, Yang Yang, James K Guest,  and Ankit Srivastava, “3-d phononic crystals with ultra-wide band gaps,” Scientific Reports 7 (2017).
  • Bilal and Hussein (2011) Osama R. Bilal and Mahmoud I. Hussein, “Ultrawide phononic band gap for combined in-plane and out-of-plane waves,” Phys. Rev. E 84, 065701 (2011).
  • Lu et al. (2013) Lirong Lu, Takashi Yamamoto, Masaki Otomori, Takayuki Yamada, Kazuhiro Izui,  and Shinji Nishiwaki, “Topology optimization of an acoustic metamaterial with negative bulk modulus using local resonance,” Finite Elements in Analysis and Design 72, 1–12 (2013).
  • Krushynska et al. (2014) A.O. Krushynska, V.G. Kouznetsova,  and M.G.D. Geers, “Towards optimal design of locally resonant acoustic metamaterials,” Journal of the Mechanics and Physics of Solids 71, 179 – 196 (2014).
  • Yang et al. (2016) Xiong Wei Yang, Joong Seok Lee,  and Yoon Young Kim, “Effective mass density based topology optimization of locally resonant acoustic metamaterials for bandgap maximization,” Journal of Sound and Vibration 383, 89 – 107 (2016).
  • Sébastien Guenneau (2013) Richard V. Craster Sébastien Guenneau, Acoustic Metamaterials: Negative Refraction, Imaging, Lensing and Cloaking, 1st ed., Springer Series in Materials Science 166 (Springer Netherlands, 2013).
  • Nemat-Nasser and Srivastava (2011a) Sia Nemat-Nasser and Ankit Srivastava, “Overall dynamic constitutive relations of layered elastic composites,” Journal of the Mechanics and Physics of Solids 59, 1953–1965 (2011a).
  • Nemat-Nasser and Srivastava (2011b) Sia Nemat-Nasser and Ankit Srivastava, “Negative effective dynamic mass-density and stiffness: Micro-architecture and phononic transport in periodic composites,” AIP Adv 1, 041502 (2011b).
  • Nemat-Nasser et al. (2011) S Nemat-Nasser, JR Willis, A Srivastava,  and AV Amirkhizi, “Homogenization of periodic elastic composites and locally resonant sonic materials,” Physical Review B 83, 104103 (2011).
  • Srivastava and Nemat-Nasser (2012) Ankit Srivastava and Sia Nemat-Nasser, “Overall dynamic properties of three-dimensional periodic elastic composites,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 468, 269–287 (2012).
  • Srivastava and Nemat-Nasser (2014) Ankit Srivastava and Sia Nemat-Nasser, “On the limit and applicability of dynamic homogenization,” Wave Motion 51, 1045–1054 (2014).
  • Willis (1983) JR Willis, “The overall elastic response of composite materials,” Journal of Applied Mechanics 50, 1202–1209 (1983).
  • Willis (2009) JR Willis, “Exact effective relations for dynamics of a laminated body,” Mechanics of Materials 41, 385–393 (2009).
  • Willis (2011) JR Willis, “Effective constitutive relations for waves in composites and metamaterials,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 467, 1865–1879 (2011).
  • Sridhar et al. (2018) A. Sridhar, L. Liu, V.G. Kouznetsova,  and M.G.D. Geers, “Homogenized enriched continuum analysis of acoustic metamaterials with negative stiffness and double negative effects,” Journal of the Mechanics and Physics of Solids 119, 104 – 117 (2018).
  • Muhlestein and Haberman (2016) Michael B. Muhlestein and Michael R. Haberman, “A micromechanical approach for homogenization of elastic metamaterials with dynamic microstructure,” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 472 (2016), 10.1098/rspa.2016.0438http://rspa.royalsocietypublishing.org/content/472/2192/20160438.full.pdf .
  • Wales and Doye (1997) David J Wales and Jonathan PK Doye, “Global optimization by basin-hopping and the lowest energy structures of lennard-jones clusters containing up to 110 atoms,” The Journal of Physical Chemistry A 101, 5111–5116 (1997).
  • Jones et al. (2014) Eric Jones, Travis Oliphant,  and Pearu Peterson, “SciPy: open source scientific tools for Python,”  (2014).
  • Craig and Bampton (1968) R-R. Craig and M C C Bampton, “Coupling of Substructures for Dynamics Analyses,” AIAA Journal  (1968), 10.2514/3.4741.
  • Krattiger and Hussein (2014) Dimitri Krattiger and Mahmoud I. Hussein, “Bloch mode synthesis: Ultrafast methodology for elastic band-structure calculations,” Physical Review E - Statistical, Nonlinear, and Soft Matter Physics  (2014), 10.1103/PhysRevE.90.063306.
  • Mead (1975) D.J. Mead, “Wave propagation and natural modes in periodic systems: I. mono-coupled systems,” Journal of Sound and Vibration 40, 1 – 18 (1975).
  • Amirkhizi (2017) Alireza V. Amirkhizi, “Homogenization of layered media based on scattering response and field integration,” Mechanics of Materials 114, 76–87 (2017).
  • Willis (2012) JR Willis, “The construction of effective relations for waves in a composite,” Comptes Rendus Mécanique 340, 181–192 (2012).
  • Deaton and Grandhi (2014) Joshua D Deaton and Ramana V Grandhi, “A survey of structural and multidisciplinary continuum topology optimization: post 2000,” Structural and Multidisciplinary Optimization 49, 1–38 (2014).
  • Sigmund and Maute (2013) Ole Sigmund and Kurt Maute, “Topology optimization approaches,” Structural and Multidisciplinary Optimization 48, 1031–1055 (2013).
  • Cadman et al. (2013) Joseph E Cadman, Shiwei Zhou, Yuhang Chen,  and Qing Li, “On design of multi-functional microstructural materials,” Journal of Materials Science 48, 51–66 (2013).
  • Osanov and Guest (2016) Mikhail Osanov and James Guest, “Topology Optimization for Architected Materials Design,” Annual Review of Materials Research 46, 211–233 (2016).
  • Bendsøe (1989a) Martin P Bendsøe, ‘‘Optimal shape design as a material distribution problem,” Structural Optimization 1, 193–202 (1989a).
  • Rozvany and Zhou (1991) G I N Rozvany and Ming Zhou, “The COC algorithm, part I: cross-section optimization or sizing,” Computer Methods in Applied Mechanics and Engineering 89, 281–308 (1991).
  • Osher and Sethian (1988) Stanley Osher and James A Sethian, “Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations,” Journal of Computational Physics 79, 12–49 (1988).
  • Allaire et al. (2002) Grégoire Allaire, François Jouve,  and Anca-Maria Toader, “A level-set method for shape optimization,” Comptes Rendus Mathematique 334, 1125 – 1130 (2002).
  • Allaire et al. (2004) Grégoire Allaire, François Jouve,  and Anca-Maria Toader, “Structural optimization using sensitivity analysis and a level-set method,” Journal of Computational Physics 194, 363 – 393 (2004).
  • Osher and Santosa (2001) Stanley J. Osher and Fadil Santosa, “Level set methods for optimization problems involving geometry and constraints: I. frequencies of a two-density inhomogeneous drum,” Journal of Computational Physics 171, 272 – 288 (2001).
  • Bendsøe (1989b) M. P. Bendsøe, “Optimal shape design as a material distribution problem,” Structural Optimization 1, 193–202 (1989b).
  • (48) T. Belytschko, S. P. Xiao,  and C. Parimi, “Topology optimization with implicit functions and regularization,” International Journal for Numerical Methods in Engineering 57, 1177–1196.
  • Ashby (2005) Michael F. Ashby, “Materials Selection in Mechanical Design,” Design  (2005), 10.1016/B978-1-85617-663-7.00011-4.
  • Pham et al. (2013) K. Pham, V.G. Kouznetsova,  and M.G.D. Geers, “Transient computational homogenization for heterogeneous materials under dynamic excitation,” Journal of the Mechanics and Physics of Solids 61, 2125 – 2146 (2013).
  • Rytov (1956) SM Rytov, “Acoustical properties of a thinly laminated medium,” Soviet Physics-Acoustics 2, 68 (1956).
  • Srivastava (2015b) Ankit Srivastava, “Causality and passivity in elastodynamics,” in Proc. R. Soc. A, Vol. 471 (The Royal Society, 2015) p. 20150256.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description