Topology design and optimization of nonlinear periodic materials

Author's Accepted Manuscript

Topology Design and Optimization of Nonlinear Periodic Materials Kevin L. Manktelow, Michael J. Leamy, Massimo Ruzzene

www.elsevier.com/locate/jmps

PII: DOI: Reference:

S0022-5096(13)00139-7 http://dx.doi.org/10.1016/j.jmps.2013.07.009 MPS2325

To appear in:

Journal of the Mechanics and Physics of Solids

Received date: 11 April 2013 Revised date: 9 July 2013 Accepted date: 18 July 2013 Cite this article as: Kevin L. Manktelow, Michael J. Leamy, Massimo Ruzzene, Topology Design and Optimization of Nonlinear Periodic Materials, Journal of the Mechanics and Physics of Solids, http://dx.doi.org/10.1016/j.jmps.2013.07.009 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Topology Design and Optimization of Nonlinear Periodic Materials Kevin L. Manktelow Michael J. Leamy Massimo Ruzzene George W. Woodruff School of Mechanical Engineering Georgia Institute of Technology, 771 Ferst Drive NW. Atlanta, GA 30332-0405, USA July 24, 2013 Abstract This article explores optimal topologies yielding large band gap shifts in one- and two-dimensional nonlinear periodic materials. The presence of a nonlinearity in a periodic material system results in amplitude-dependent dispersion behavior, leading to novel wave-based devices such as tunable filters, resonators, and waveguides. The performance of these devices over a broad frequency range requires large, tunable band gaps, motivating the present study. Consideration of a one-dimensional bilayer system composed of alternating linear and nonlinear layers shows that optimal designs consist of thin, compliant nonlinear layers. This is at first surprising considering the source of the shift originates from only the nonlinear layer; however, thin layers lead to localized stresses that activate the nonlinear character of the system. This trend persists in two-dimensional materials where optimization studies are performed on plane-stress models discretized using bilinear Lagrange elements. A fast algorithm is introduced for computing the dispersion shifts, enabling efficient parametric analyses of two-dimensional inclusion systems. Analogous to the one-dimensional system, it is shown that thin ligaments of nonlinear material lead to large dispersion shifts and group velocity variations. Optimal topologies of the two-dimensional system are also explored using genetic algorithms aimed at producing large increases in complete band gap width and shift, or group velocity variation, without presupposing the topology. The optimal topologies that result resemble the two-dimensional inclusion systems, but with small corner features that tend to enhance the production of dispersion shift further. Finally, the study concludes with a discussion on Bloch wave modes and their important role in the production of amplitude-dependent dispersion behavior. The results of the study provide insight and guidance on selecting topologies and materials which can yield large amplitude-dependent band gap shifts and group velocity variations, thus enabling sensitive nonlinear devices.

1

Introduction

Structures with periodically varying geometry and material properties are known to exhibit advantageous wave dispersion leading to negative refractive indices, frequency band gaps and spatial beaming [1, 2]. On small scales, periodic compositions of materials with contrasting mass and stiffness properties are termed phononic crystals [3, 4, 5]. Phononic crystals are typically considered to operate in regimes where a linear constitutive relationship (e.g., stressstrain) provides an adequate representation. For higher intensity wave propagation, however, weak nonlinearities may influence performance. It is now known that cubic nonlinearities give rise to amplitude-dependent frequency shifting and thus a shift in band gap location [6]. Amplitude-tunable band gaps and group velocity can lead to tunable waveguides, acoustic filters, acoustic mirrors, transducers, and diodes [2, 3, 4, 5, 7, 8, 9]. Elastic one-dimensional (1D) layered systems and two-dimensional (2D) plane stress systems are practically important from a design perspective in that various topologies may be readily fabricated. Layered systems often serve as fundamental building blocks in topology optimization studies [10, 11]. Compound layered systems (sometimes termed sequentially ranked laminates), for example, find use in topology optimization for microstructure distribution [12]. However, structures composed of sequentially ranked laminate microstructures derive their overall response from homogenized properties since typical wavelengths are much greater than the unit cell length. Optimization of

1

the one-dimensional layered system has been considered by others with regard to optimizing band gap size/location (e.g. Refs [13, 14]). Topology optimization of a linear binary inclusion system was considered in [11] with the goal of optimizing dispersion band structure and frequency response characteristics. Genetic algorithms have been used extensively by the phononic/photonic crystal community for identifying configurations which achieve complete band gaps [15, 16, 17, 18, 19]. Genetic algorithms identify optimal solutions through iterative population generation, mutation, crossover, and selection processes. A fitness function ranks individual members of a population based on a specified criteria. The genetic algorithm handles binary material systems directly, whereas gradient-based optimization routines usually require continuous design variables and explicit gradient expressions. This partially explains why genetic algorithms have gained such popularity in phononic and photonic crystal analysis where it is known a priori that stark material contrasts are essential to producing the desired results. For nonlinear problems with many degrees of freedom, evaluating the jacobian is computationally intensive and time consuming. Since genetic algorithms do not require a jacobian and implementations are readily available, a genetic algorithm has been employed as a tool to search design spaces for preferred designs. We note that there has been much work devoted to efficient gradient-based optimization for phononic and photonic crystals recently, as well, including algorithm development as well as practical application implementations (waveguide design, energy control, etc.) [11, 10, 20, 21, 8, 22]. These methods may be more appropriate for problems with more degrees of freedom and their use is recommended for additional study in future nonlinear dispersion analysis. In this paper we focus on binary material distributions where one material contains a cubic stress-strain nonlinearity. A previously-developed perturbation analysis for Bloch wave propagation provides an analytical expression for the frequency correction as a function of linear and nonlinear material parameters, Bloch wave modes, and frequencies [23]. Unit cells which contain only a few degrees of freedom may be analyzed in closed-form. However, finiteelement discretized domains with nonlinear constitutive laws lead to expressions requiring numerical evaluation. To make computations practical, an algorithm which efficiently evaluates nonlinear frequency shifts is provided. This is especially important in topology optimization routines where nonlinear dispersion relations are repeatedly evaluated. Using the aforementioned tools, the influence of constitutive nonlinearity on a two-layer unit cell with variable layer widths is first investigated parametrically. This approach offers a useful picture of general trends associated with nonlinear dispersion shifts which may be used in designing/optimizing more complex systems. An example amplitudetunable filter exemplifies the utility of topologies optimized for sensitive band gap shifts. The system considered resembles the system discussed in [11], but considers nonlinear stress-strain relationships. Two-dimensional designs informed by study of the one-dimensional layered system are considered next with regard to group velocity and first band gap width, followed by a genetic algorithm-based topology optimization study aimed at optimizing unit cell designs with regard to band gap width. The study concludes with a discussion on the role Bloch modes play in increasing dispersion shift sensitivity, and general guidelines for designing unit cells with large amplitude-dependent band gap shifts and group velocity variations.

2 2.1

Theoretical background System model

The developments that follow assume a periodic elastic structure whose dynamics are governed by the conservation of momentum in the absence of externally applied forces ∇ · σ (r) = ρ

∂ 2 u(r) , ∂t 2

∀r ∈ Ω,

(1)

where σ denotes the elastic stress tensor, ρ is the material density, and u is the displacement field vector. We consider one-dimensional (1D) and two-dimensional (2D) computational domains. In the following, we report in detail the equations governing the 2D case, with the 1D case being omitted as a particular case. In the 2D case, each material point with in-plane position r = x e1 + y e2 for (x , y) ∈ Ω undergoes in-plane motion. Then, the displacement vector u = [u(r), v(r)] describes motion along the e1 and e2 directions, respectively. With

2

these assumptions, Eqs. (1) are expressed as: ∂ 2u ∂ σxx ∂ σxy + =ρ 2, ∂x ∂y ∂t ∂ σxy ∂ σyy ∂ 2v + =ρ 2, ∂x ∂y ∂t

∀(x, y) ∈ Ω,

(2a)

∀(x, y) ∈ Ω,

(2b)

where σi j denotes the stress tensor, h denotes thickness in the e3 direction, and ρ denotes volume density. We adopt Voigt notation such that components of the stress tensor are given by σi = [σxx , σyy , σxy ]T . A nonlinear elastic constitutive law is considered      σxx c11 (εxx ) c12 0 εxx σyy  =  c12 c22 (εyy ) 0   εyy  = Cεε , (3) σxy 0 0 c66 (εxy ) 2εxy where C denotes the stiffness matrix and strain components ε are defined according to       εxx ∂ /∂ x 0 u  εyy  =  0 ∂ /∂ y = Tu, v 2εxy ∂ /∂ y ∂ /∂ x

(4)

where T denotes the derivative matrix. Strain-dependence in off-diagonal terms c12 is neglected. The coefficient functions c11 , c22 , c12 , and c66 describe an isotropic material with cubic nonlinearities and are given by [24] E1 + Γx εxx 2 1 − ν2 E1 c22 = + Γy εyy 2 1 − ν2 E1 + Γxy εxy 2 c66 = 2(1 + ν) νE1 c12 = , 1 − ν2

c11 =

(5a) (5b) (5c) (5d)

where E1 and ν denote the linear modulus and Poisson ratio, respectively. The cubic nonlinearity coefficients Γx , Γy , and Γxy may be derived from experimental data, or from well-known nonlinear constitutive laws such as the NeoHookean, Mooney-Rivlin [25], or Murnaghan Potential models [26, 27]. Small strains and weak nonlinearity are enforced in the constitutive equations by introducing Γˆ x , Γˆ y , and Γˆ xy according to Γx = ε Γˆ x , Γy = ε Γˆ y , and Γxy = ε Γˆ xy . These terms result in a weakly nonlinear contribution when the stresses due to the nonlinear terms are much less than the stresses due to linear terms. The updated plane stress equations are      ∂ ∂u ∂v ∂ ∂u ∂v ∂ 2u c11 + c12 + c66 + + ε fxNL = ρ 2 (x, y) ∈ Ω (6a) ∂x ∂x ∂y ∂y ∂y ∂x ∂t      ∂ ∂u ∂v ∂ ∂u ∂u ∂ 2u c66 + + c22 + c12 + ε fyNL = ρ 2 (x, y) ∈ Ω. (6b) ∂x ∂y ∂x ∂y ∂x ∂x ∂t and the nonlinear functions fxNL and fyNL are given by "   # "  3 # 3 ∂ ∂ u ∂ ∂ u ∂ v fxNL = Γˆ x + Γˆ xy + ∂x ∂x ∂y ∂y ∂x "  "   #  # ∂ ˆ ∂u ∂v 3 ∂ ˆ ∂u 3 NL fy = + + . Γxy Γy ∂x ∂y ∂x ∂y ∂y

3

(7)

2.2

Finite-element discretization

The nonlinear analysis of the continuous equations of motion is impractical for general periodic unit cells; therefore, finite element-based discretized equations of motion will be employed herein. As the objective of the study is the dispersion analysis for periodic systems, the computational domain Ω of interest is a unit cell, which is discretized into S a suitable number of elements Ω(e) , such that Ω = (e) Ω(e) . Consider for a moment the plane stress equation (2) for an element domain Ω(e) . The displacement fields u(x, y,t) and v(x, y,t) over the element domain Ω(e) are represented with n interpolation functions Ni , i = 1..n and nodal displacements ui (t) and vi (t) according to u(x, y,t) = ∑ni=1 Ni (x, y)ui (t) and v(x, y,t) = ∑ni=1 Ni (x, y)vi (t). We choose to employ quadrilateral elements primarily because manual assembly is straight-forward and produces a structure mesh. A structured mesh improves evaluation speed of the nonlinear force vector, which is particularly beneficial in optimization routines (to be discussed later). However, a weak formulation of the governing equations is next presented without bias towards any particular element type. Each of the governing equations is multipled by an interpolation function Ni (x, y) and integrated over the element domain to produce the weak form of the governing equations: (see [28] for complete details)       ZZ ∂u ∂ Ni ∂v ∂u ∂v ∂ Ni NL c11 + + ε fui + ρNi u¨ dΩ(e) + c12 c66 + h ∂x ∂y ∂y ∂y ∂x Ω(e) ∂ x (8a)      Z ∂v ∂u ∂v ∂u + c12 nx + c66 + ny ds = 0 −h Ni c11 ∂x ∂y ∂y ∂x ∂ Ω(e) and, similarly for the second equation:       ZZ ∂ Ni ∂u ∂v ∂ Ni ∂u ∂v h c66 + + c12 + c22 + ε fvNL + ρN v ¨ dΩ(e) i i ∂y ∂x ∂y ∂x ∂y Ω(e) ∂ x

(8b)

      ∂u ∂v ∂u ∂v + nx + c12 + c22 ny ds. = 0 −h Ni c66 ∂y ∂x ∂x ∂y ∂ Ω(e) Z

Equations (8) represent the weak form of the plane stress equations. The first term in either equation is the field equation, while the latter term represents boundary conditions. Integration has been carried out over the entire domain, including the e3 direction, which is responsible for producing the thickness h in the expressions. The nx and ny terms denote components of the boundary normal in e1 and e2 directions, respectively. However, as the system considered is infinite in extent the latter term is not explicitely considered beyond this point. Elemental mass, stiffness, and nonlinear force vectors are formed by introducing interpolation functions and nodal coordinates ui and vi into the weak form of the governing equations Eq. (8), carrying out the integration, and collecting terms as appropriate. The discretized matrix form of the governing equations for an element can be expressed Mu¨ + Ku + εfNL = 0

(9)

where M, K, fNL , and u and denote element mass and stiffness matrices, nonlinear force vector, and displacement vector. Exact expressions for the mass and stiffness matrix terms are given by standard finite element expressions M=h

Z Ω(e)

ρNT NdΩ(e)

and K = h

Z Ω(e)

BT CBdΩ(e) ,

(10)

where N is the vector of shape functions and B = TN is the matrix of shape function derivatives [28, 29]. Each of Eqs. (8) contains a nonlinear force term that arises from nonlinearities in the constitutive law. Their contribution to the elemental nonlinear force vector, for the special case of orthotropic nonlinearity, are given by       ZZ ∂ Ni ˆ ∂ u 3 ˆ ∂ u ∂ v 3 = h fuNL + Γ + dΩ(e) Γ x y i ∂x ∂y ∂x Ω(e) ∂ x (11)   3  3  ZZ ∂ N ∂ u ∂ v ∂ u i fvNL =h Γˆ x + + Γˆ y dΩ(e) , i ∂y ∂x ∂y Ω(e) ∂ y where fuNL and fvNL denote nonlinear restoring forces applied to the ith node in the e1 and e2 directions, respectively. i i An assembly of nonlinear elastic plane stress elements of the form given in Eq. (10) produces a system of weakly 4

nonlinear governing equations which fits the general form required for the perturbation analysis. In what follows, the domain Ω is defined as the union of 9 adjacent unit cells as depicted in Fig. A.1 and matrices M, K, u, and fNL are redefined such that they refer to assembled matrices (rather than elemental matrics). Then, a nonlinear Bloch wave analysis follows directly as decribed in Appendix A. The evaluation of the nonlinear force vector can be a computationally intensive task that prevents viability of the topology optimization. The cost of this operation is greatly reduced herein through the symbolic evaluation of the integrals in Eq. (11) for each interpolation function Ni . The use of uniform meshes with equal-size elements leads to symbolic expressions for the elemental nonlinear forces that are identical for all elements. This further streamlines the nonlinear term evaluation, whose symbolic expression is saved in the form of functions whose arguments are the elemental nodal displacements and the relevant constitutive parameters. Nonlinear band structure calculations require repeated assembly and evaluation of fNL for multiple wave vectors. Optimization codes that use band structure as part of an objective function may require evaluation of fNL many thousands of times. Further improvements are obtained by using compiled codes that utilize parallel assembly and efficient evaluation of c1 , the first Fourier coefficient of the nonlinear force vector fNL (as described in Appendix A). An efficient algorithm for evaluating the first Fourier coefficient c1 is described next.

2.3

Efficient algorithm for calculation of Fourier coefficients

The nonlinear frequency correction requires a Fourier representation of fNL , and specifically the evaluation of the first coefficient c1 . Previous work has evaluated the function at discrete time instants τi ∈ [0, 2π] and numerically computed the first Fourier coefficient using either numerical integration or FFT algorithms (see Refs. [6, 30]). Both of these methods require fine temporal discretization over a single fundamental period to obtain accurate estimates for the c1 coefficient. Individual band structure computations do not suffer significantly from this method; however, this becomes impractical for parametric studies and topology optimization. An improved evaluation method takes advantage of the polynomial nature of the constitutive law, and thus reduces the number of function evaluations to four. Individual elements of the (cubic) nonlinear force vector may generally be expressed in the form fNL i (u(τ)) =

∑ αlmn ul um un ,

(12)

l,m,n

where each nodal displacement ui is 2π periodic in dimensionless time τ. Thus, the functional form is known a priori to be fNL (τ) = A1 cos(τ) + A3 cos(3τ) + B1 sin(τ) + B3 sin(3τ) (13) as a result of trigonometric identities. This equation is linear in the four unknown amplitude vectors A1 , A3 , B1 , B3 . The first complex Fourier coefficient is determined from the expression c1 = (A1 − iB1 )/2.

(14)

Thus, it remains only to compute the vectors A1 and B1 . The unknown amplitude vectors are determined by evaluating the nonlinear force vector at four distinct times to obtain a consistent set of equations. Strategically chosen times τi = [0, π/2, π/3, π/6] reduce the complexity of the resulting expressions. The system of equations obtained using these τi is given by       NL f (0) I I 0 0 A1      NL 0 0 I −I    A3  = fNL (π/2) , (15) cos(π/3)I −I sin(π/3)I 0 B1  f (π/3) NL cos(π/6)I 0 sin(π/6)I −I B3 f (π/6) where I denotes the N × N identity matrix and 0 denotes the N × N zero matrix, and N denotes the total number of degrees of freedom in the system. The unique solution to this linear system of equations provides the desired c1 coefficient. The solution for A1 and B1 is 1 1 A1 = fa − √ fb , B1 = fb − √ fa , (16) 3 3 where the substitutions fa = fNL (0) + fNL (π/3) and fb = fNL (π/2) + fNL (π/6) have been made. This enhanced formulation dramatically decreases the computation time by reducing the number of fNL evaluations to an absolute minimum of four. 5

2.4

Shift sensitivity analysis

Recall that the systems considered all contain cubic nonlinearities or exhibit one following Taylor expansion. Recall also, as described in Appendix A, that O(ε 0 ) Bloch wave solutions take the form u(0) = Aφφ /2 exp(−iτ) + c.c. Since the nonlinear force term has been defined as the cubic contribution of a general nonlinear force function in Eq. (5d), it is possible to factor out A3 from each element of c1 without loss of generality (see Eq. (12)). In doing so, the nonlinear dispersion corrections may be computed using the expression µ )A2 , ω1, j = Kω (µ

Kω ≡

φ Hj c1

(17)

ω0, j φ Hj Mφφ j

where we term Kω the shift sensitivity, and redefine c1 such that it is no longer parametrized by A. The shift sensitivity describes the quantitative and qualitative nature of the dispersion shifts for a given nonlinearity, independent of wave mode amplitude. Large Kω are found for systems where large frequency shifts occur for relatively small amplitudes and are desirable for designing tunable systems which respond to signal gain. Softening nonlinearities shift the dispersion curves downward in frequency, implying Kω < 0, while hardening nonlinearities shift the dispersion curves upward, implying Kω > 0, for a given wave vector. Using Eq. (37) and the sensitivity Kω , the updated dispersion relation is given by µ ) = ω0, j (µ µ ) + εA2 Kω (µ µ ) + O(ε 2 ). ω j (µ (18) Note that a group velocity calculation based on Eq. (18) yields a correction to the group velocity, (0)

(1)

cg = ∇µ ω0, j + εA2 (∇µ Kω ) + O(ε 2 ) = cg + εcg + O(ε 2 ),

(19)

(0)

where ∇µ denotes a gradient with respect to the components of the wave vector, cg denotes the group velocity of the (1)

linear system, and cg = A2 ∇µ Kω denotes the group velocity correction.

3 3.1

Strain-induced dispersion shift in 1D layered systems Model description

The bi-layer structure of Fig. 3.1 provides a straightforward context for exploring the manner in which topological changes can enhance advantageous behavior present in nonlinear band structures. In this example, material A includes a cubic nonlinearity in its stress-strain relationship, while material B includes only a linear term. The topology of this 1D system is defined by the length fraction a/d for a given unit cell length d. We consider, for material A, a

Figure 3.1: Unit cell p for a bi-layer system hardening-type nonlinearity as previously described, while material B is governed by a linear stress-strain constitutive relationship. Analytical solutions for the dispersion relation of a linear bi-layer system are widely available and can be expressed as [31, 32, 33],           ωa ωb 1 ZA ZB ωa ωb cos(µ1 d) = cos cos − + sin sin , (20) cA cB 2 ZB ZA cA cB 6

p p where c j = E1 j /ρ j denotes the phase speeds in a homogeneous material, and Z j = E1 j ρ j the impedance of the jth layer, with j being either A or B. Equation (20) indicates that band gaps arising in the linear layered system derive from an impedance mismatch caused by stiffness and density contrast at the material interface. The influence of stiffness contrast on the shift sensitivity parameter Kω is considered next to identify configurations exhibiting optimized band gap width and frequency shifts.

3.2

Parametric analysis of dispersion shifts

Dispersion shifts for layered structures generally fall into two categories: individual branch shifts and collective band gap shifts (both center frequency and widths). Group velocity shifts are arguably an additional category; however, large dispersion branch shifts generally imply large group velocity shifts. Figure 3.2 schematically illustrates some of the effects that “small” and “large” shifts may have on dispersion. The middle figure, for example, illustrates that large acoustic branch shifts which maintain the band gap width, but alter the center frequency, may be possible. In contrast, large dispersion shifts to the acoustic branch combined with small optical branch shifts might result in band gap closure (and a shifted center frequency).

Figure 3.2: Optimizing 1D layered structures may achieve large shifts to individual dispersion branches (middle) and also to band gap widths (right). Dotted lines indicate a low-amplitude (linear) system, while solid lines indicate the nonlinearly shifted dispersion relation. As an initial investigation into topology optimization of nonlinear layered systems, the length fraction a/d and linear stiffness ratio E1B /E1A are varied to determine their dispersion effect. Consideration is given to metrics related to the first band gap. In particular, ∆ω1 and ∆ω2 denote the frequency shift for the first and second dispersion branches evaluated at the edge of the Brillouin zone (µ1 = π). The domain for the bi-layer rod is discretized into 200 linear finite elements to obtain fine resolution for small layer fractions a/d. The system parameters used in Ref. [34] (E1A = 1, E1B = 4, and ρA = ρB = 1) are chosen for comparison purposes, which also contains an account of the finite-element discretization procedure. Figure 3.3a depicts the computed frequency shift of the acoustic branch ∆ω1 ∝ Kω (π). The maximum acoustic branch shift for a given material stiffness ratio occurs in nearly all cases for nonlinear length fractions a/d  0.5. The white line indicating the locus of maximum frequency shift tends towards a length fraction of approximately 3% at high stiffness ratios. Figure 3.3b, on the other hand, reveals that maximum corrections to the second dispersion branch ∆ω2 appear for larger length fractions a/d ≈ 0.35. Therefore, the length fraction corresponding to a maximum first band gap width lies somewhere between the two optimum shift values for any given stiffness ratio. Shifts as depicted in Fig. 3.2, for example, may be obtained for locating layer fractions and stiffness ratios where shifts have similar magnitudes. A configuration with a stiffness ratio of 12 and layer fraction of 0.67 yields this type of result (there are many other possibilities). Figure 3.4 further details these results by considering a cross-section S1 of the frequency shift maps taken at E1B /E1A = 4. Figure 3.4a depicts the frequency shifts to the first and second dispersion branches. The maximum shift to the acoustic branch occurs at a layer thickness of approximately 12% in this case. The second dispersion branch, in contrast, does not exhibit a discernible optimum shift value. Instead, the correction ∆ω2 increases towards an asymptote representing the shift expected for a single-layer nonlinear system

7

(a) First branch shift, ∆ω1 (π)

(b) Second branch shift, ∆ω2 (π)

Figure 3.3: (color online) Frequency shift maps evaluated at the edge of the Brillouin zone for a bi-layer unit cell (ε = 0.1). White lines indicate the locus of maximum shift for a given linear modulus or layer thickness. (a/d = 1.0), reaching a slight maximum value near a/d = 0.9. A perturbation calculation based on a continuous, homogeneous system gives the value of this asymptote (details in Appendix B). Indeed, an examination of the frequency shift map in Fig. 3.3b shows that at stiffness ratios E1B /E1A . 7, a clear maxima is not discernible. A relatively sharp transition between small and large frequency shifts as the length fraction increases is noted for the second dispersion branch that may be advantageous in nonlinear acoustic switches, diodes, and superprism crystals [35]. Figure 3.4b depicts how shifts to the first and second dispersion branches combine to vary the band gap width. The maximum gap width for a bi-layered system using these materials occurs at a/d ≈ 0.30. These results indicate that, contrary to what might be expected, less nonlinear material (i.e., small a/d nonlinear layer fractions) corresponds to larger shifts for this example case, and highlights the importance of systematic topology optimization and design. We note that the combined effect of a large Kω for the acoustic branch and a small Kω for the optical branch (evaluated at the band edge) may result in a complete closure of the first band gap at large amplitudes for hardening nonlinearities, or a significant size increase for softening nonlinearities. Figure 3.5 depicts the resulting dispersion diagram for a system with a/d = 0.05 and E1B /E1A = 20. Several important features of this nonlinear dispersion diagram may influence the design of nonlinear phononic crystals: 1. The band gap width of the nonlinear system decreases with amplitude due to a hardening nonlinearity. The change in band gap width is dominated by ∆ω1 since ∆ω2 is negligible (small Kω ). Operation at frequencies above 10 rad s−1 , for low amplitudes, results in complete attenuation. A slight amplitude increase shifts the dispersion curve upward, allowing unattenuated wave propagation. 2. Frequency upshifts are greatest near the edge of the Brillouin zone where the group velocity is zero. Similar results have been demonstrated by the photonic crystal community for Kerr-type nonlinear electromagnetic materials [36, 35]. 3. The group (and phase) velocity for wave numbers µ1 < π are also significantly shifted; in fact, at greater amplitudes the material appears less dispersive because there exist a greater range of wave numbers with the same phase speed. This concept may be important in realizing soliton solutions and mitigating unwanted dispersion. These features are explored for use in an amplitude-tunable filter next.

3.3

Application: Amplitude-tunable filter

A time-domain simulation demonstrates how nonlinear layered systems behave as simple tunable filters. We consider a finite-element model of a bi-layered system composed of 80 elements per unit cell with material properties ρA = ρB = 1, E1A = 1, E1B = 20, and hardening nonlinearity Γx = 1. A length fraction a/d = 0.05 achieves a large acoustic branch shift (see Fig. 3.3a). A harmonic stress σxx (0,t) = σ0 cos(ωt) at the left boundary excites Bloch wave propagation through the system. The right boundary is located sufficiently far from the excitation such that no reflections are produced. 8

Nonlinear Shift

0.4 ∆ω1, Branch 1

0.3

∆ω2, Branch 2

0.2

a ≈ 12%

0.1 0

0

0.2

0.4 0.6 Layer fraction, a/d

0.8

1

0.8

1

1st Bandgap Width

(a) Nonlinear branch shifts

2 1.5 a ≈ 30%

1 0.5 0

0

0.2

0.4 0.6 Layer fraction, a/d (b) First band gap width

Figure 3.4: Cross-sections of the frequency shift maps shown in Fig. 3.3a and 3.3b evaluated at E1B /E1A = 4. The gap width for low- and high-amplitude systems are shown in Fig. 3.4b with dashed and solid lines, respectively, to illustrate the effect of the nonlinearity on the band gap width at this cross-section.

Figure 3.5: Dispersion diagram (a/d = 0.05 and E1B /E1A = 20) for the simulated system evaluated for low- and high-amplitude excitations. The marker (◦◦) shows verification results obtained via a 2D-FFT of the finite element simulation response. Figure 3.5 depicts dispersion diagrams for low (dashed) and high (solid) amplitude excitations using normalized p frequency ωd/cA (cA = E1A /ρA ) and normalized wavenumber µ. The combination of amplitude (σ0 ) and nonlinearity (Γx ) considered yields a shift to the acoustic branch that is on the order of the first band gap width. Indeed, increasing the propagation amplitude would close this band gap entirely. This effect allows the first band gap to oper-

9

0.5

Displacement

Displacement

ate as an amplitude-dependent filter or switch that permits propagation at certain frequencies only for high amplitudes. The excitation frequency ω = 11.0 rad s−1 falls within the band gap for an essentially linear system (low-amplitude excitation); in contrast, ω = 11.0 rad s−1 falls well within the pass band for a weakly nonlinear system (high-amplitude). We note that a softening nonlinearity (Γx = −1), for example, produces shifts opposite to those pictured in Fig. 3.5 and would increase the bandgap width. However, caution is necessary when considering softening nonlinearities due to inherent instabilities at large amplitudes without higher-order hardening terms. Two simulations are performed: the first at a low-amplitude excitation where the propagation is expected to decay exponentially, and a second at high-amplitude excitation where propagation is expected. Figure 3.6 depicts time traces, for both low- and high-amplitude excitations, at the input boundary (first unit cell) and at the 10th unit cell. The

0

−0.5

0

2

4

6

8

0.5

0

−0.5

10

Time

0

2

4

6

8

10

Time

(a) Low amplitude

(b) High amplitude

Figure 3.6: Time domain results for probes located at the input boundary (unit cell 1, solid gray) and at the 10th unit cell (solid black). The low-amplitude signals are immediately attenuated and consist primarily of transient broadband content. In contrast, the high-amplitude exictation signals clearly propagate, though with some waveform distortion as a result of the nonlinearity used to enable transmission. low-amplitude simulation (left subfigure) shows immediate exponential attenuation and small transient propagation broadband signal content with frequency content in the first few passband regions. The response for higher-amplitude excitation is markedly different: the input signal is propagated, but with some harmonic distortion as the signal progresses through the system. Fourier transforms1 of the displacement signal at the input and 10th unit cell show a 25% decrease in the fundamental frequency ω accompanied by growth of higher harmonics (e.g., 6x increase in 3ω).

Figure 3.7: Normalized average power distributions (unit in-plane thickness) for low- and high-amplitude excitations (left subfigures). Evanescent decay of power in space is observed, as expected, for the low-amplitude system. The amplitude of the Fourier transform for each signal at ω (black, left axis labels) and 3ω (red, right axis labels) is shown as a function of space in the adjacent subfigures. A two-dimensional fast Fourier transform provides frequency-domain information from the simulation data that corresponds directly to the dispersion diagram, and facilitates a one-to-one comparison between the predicted disper1A

flat top window was applied in order to extract accurate amplitudes

10

ω

0

Normalized Response

10

High−amplitude Low−amplitude

3ω

5ω

7ω

9ω

−2

10

−4

10

Input frequency 0

50

100 Frequency [rad/s]

150

200

(a) Displacement FFT

Figure 3.8: Displacement signal frequency response in the 10th unit cell. The low amplitude system exhibits a response nearly two orders of magnitude less than the high-amplitude system. The frequency content of the low-amplitude response is distributed throughout the pass bands of the linear system. Super-harmonics appear with larger energy in the high-amplitude response. sion relationship and the numerical simulation. The marker in Fig. 3.5 resulting from the 2D-FFT shows that the pair ω = 11.0 rad s−1 , µ1 ≈ 2.5 dominates the signal content of the forced response simulation, in agreement with the theoretical dispersion relation. Figure 3.7 further elucidates the qualitative disparity between the two simulations by considering the distribution of power Π throughout the domain (only the first 10 unit cells are shown). Power distributions are normalized to the highest power Πmax occurring to highlight the spatial distributions. The power distribution for the low-amplitude simulation (top) decays quickly in space, while for the high-amplitude excitation (bottom) power is evenly distributed through the system. Some modulation in the harmonic amplitudes is evident in the high-amplitude case that results from nonlinear energy transfer between the first and higher harmonics. The right subfigures of Fig. 3.7 depict the evolution of the signal amplitudes at ω and 3ω as the wave propagates throughout the system (note the different scales used for ω and 3ω). Thus, the nonlinearity provides an amplitude-tunable filtering mechanism. Normalized fast Fourier transforms at the end of the 10th unit cell provide further comparison of the two simulations (Fig. 3.8). Nearly all of the frequency content of the low-amplitude signal is due to wave transients with frequency content located within the first pass band (Fig. 3.5). In contrast, the high-amplitude signal contains a dominant frequency at ω = 11.0 rad s−1 , as well as small third harmonic and sum-and-difference frequencies due to the cubic nonlinearity.

4

Analysis and optimization of 2D plane stress systems

The results of Sec. 3 demonstrate that 1D multilayer periodic structures and materials exhibit optimum configurations with respect to shift sensitivity and band gap tunability. Next, we investigate optimal configurations of 2D phononic crystals constructed with a nonlinear matrix material. The arrangement of the nonlinear material is optimized in terms of band gap width and nonlinear shifts through both parametric analysis and topology optimization. The selected optimizer is a genetic algorithm that explores a design space defined by a unit cell Ω of the 2D periodic structure with fixed dimensions.

4.1

Model description

Figure 4.1 depicts the systems under consideration. The material is constructed of a periodic arrangement of arbitrarily distributed materials A and B. Figure 4.1 depicts inclusions ΩB inserted into a matrix of ΩA ; equally valid is the inverse arrangement, or arrangements other than inclusions in an underlying matrix (e.g., random phases, etc.). The domain Ω (i) (i) is defined as the union of multiple simply-connected domains ΩA and ΩB . The unit cell domain Ω is subdivided into S (e) a number of finite elements such that Ω = (e) Ω as shown in Fig. 4.2. Each element comprises a sub-domain which takes on a binary value of 0 or 1 corresponding to material A (ΩA , nonlinear) or material B (ΩB , linear), respectively. 11

Figure 4.1: Schematic depiction of a periodic material constructed of two phases labeled A and B. Each sub-domain is considered a design variable, referenced by gi j with (i, j) = 1..48 for a total of 2(48×48) = 2.8e+14 potential design combinations. Application of a 1/8 symmetry constraint eliminates ambiguity in the design space due to translational-invariance. Moreover, 1/8 symmetry eliminates designs for which partial band gaps are preferentially larger in either the e1 or e2 directions [15, 16]. This reduces the design space from 2(48×48) to 2300 = 2e+9. It is more convenient to index the domains in a linear fashion whereby gi j is instead referenced with gi , i = 1..300. The numbering scheme indexes elements from left to right, increasing from bottom to top as shown in Fig. 4.2. 300

48

Symmetric

24

47

25

1 1

24

48

1

24

Figure 4.2: Structured unit cell mesh and numbering scheme for design variables gi j

4.2

Analysis of strain-induced dispersion shifts

A specific subset of the possible designs inspired by the results of Sec. 3 are considered next. Figure 4.3 depicts a prototypical unit cell (enclosed with dashed lines) with dimensions dx and dy describing the width and height. ΓXM symmetry calls for a square unit cell with a = b and d = dx = dy such that the lattice vectors are equal, a1 = a2 . White material represents a nonlinear matrix (Material A) within which stiffer inclusions (Material B, black) are embedded. The white sections are termed nonlinear ligaments because their purpose is to nonlinearly couple the motion of adjacent inclusions. Recall that 1D layered systems exhibit sensitive tunability and large dispersion shifts when the nonlinear material is a) restricted to thin layers, and b) compliant relative to inclusions. These qualities eliminate some common phononic crystal structures (such as air/silicon) as potential design options. Indeed, nonlinear dispersion analysis of a similar two-dimensional structure confirms that the combination of thin, compliant2 layers produces Bloch wave modes with 2 Compliance

is described here as a property of linear systems, since for large strains a hardening material may not be compliant.

12

d

y

b

a dx

Figure 4.3: A specific subset of important unit cell designs are characterized by nonlinear ligament widths a and b. The subset considers unit cells with 1/8 symmetry as indicated by the red line. localized high-strain regions. An example unit cell is highlighted in Fig. 4.4. This figure graphically depicts the dilitational strain field ε (x, y) = εxx e1 + εyy e2 produced in the unit cell for a Bloch wave mode on the ΓX portion of the acoustic dispersion surface with an approximate wavelength three to four times greater than d. Strain-localizations

Figure 4.4: Small nonlinear ligaments result in high-strain areas localized in the nonlinear ligament sections. Dashed lines indicate boundaries of rigid inclusions. These high-strain Bloch wave modes are often found in the acoustic band. enter the nonlinear force vector fNL through a few large entries acting on nodes in the ligament, rather than the cumulative effect of many smaller contributions (as in a completely homogeneous material). The specific relationship between ligament geometry and associated dispersion shifts is quite complex. The following analysis interrogates the effects of ligament width and stiffness contrast on group velocity and band gap width through systematic parameter variations. Table 1 contains the material properties used for the following analysis. Shear stiffness nonlinearities Γxy are assumed negligible. The prototypical unit cell considered has dimensions d = 1.0 and is subdivided into 48x48=2304 linear elements Ω(e) . Dispersion computations use a strain-normalization scheme such that each Bloch wave mode is normalized by the maximum strain e0 = max(εxx , εyy , εxy ) = 0.1. Nonlinear ligament thickness is increased from a = 0 (100% Material B) to a = d (100% Material A) for various material stiffness contrasts. Figures 4.5 and 4.6 depict the variation in group velocity and band gap width as a function of layer fraction a/d and 13

Material A Material B

Table 1: Material properties E1 Γx Γy Γxy 1.0 4.0 4.0 0 {5.0, 10.0, 20.0, 40.0} 0 0 0

ρ 1.00 1.00

ν 0.34 0.34

h 0.01 0.01

material stiffness contrasts. Figures 4.5a and 4.5b display group velocity shifts evaluated on the acoustic dispersion surface at µ = [π/6, 0]. The transverse acoustic branch experiences no nonlinear shift because nonlinear shear moduli Γxy are zero. Low stiffness contrasts (curve A) result in a maximum group velocity shift for layer fractions of a/d = 1 (homogeneous nonlinear material), with a small local maximum appearing near a/d = 0.25. As the stiffness contrast (1) increases (curves B–D), a global maximum appears where thin nonlinear ligaments (a/d = 0.18) produce larger cg shifts than a purely homogeneous material. Increasing material contrasts tend to push the local maximum towards smaller layer fractions, as was seen in the bi-layer system.

Relative group velocity shift

Group velocity shift, [m/s]

0.4 D

0.3

C

0.2

B 0.1

0

A Mat’l B, (Linear) 0

0.2

Mat’l A, (Nonlinear)

0.4 0.6 Layer fraction, a/d

0.8

0.15 D 0.1

B

(a) Absolute group velocity shift

A

0.05

0 1

C

0

0.2

0.4 0.6 Layer fraction, a/d

0.8

1

(b) Relative group velocity shift

Figure 4.5: (a) Absolute group velocity shift and (b) relative group velocity shifts for the longitudinal acoustic mode along the ΓX-direction. Curves A–D show results for material stiffness ratios E1B /E1A of 5 (A), 10 (B), 20 (C), and 40 (D). (0)

We note that the linear group velocity cg varies with topology because the physical system is changing. Thus, the (1) group velocity shift cg quantified relative to the linear group velocity provides another metric for defining optimal (1) (0) (0) designs. Figure 4.5b depicts the relative group velocity shift defined by cgr = |cg |/|cg | such that cg = cg (1 + cgr ). All material combinations converge to a single value (0.15) as the layer fraction increases toward a homogeneous nonlinear material. This result can be validated by considering Eq. (49) of Appendix B, repeated here with Γx for convenience with the substitution ω (0) = cµ for phase speed c,   3Γx 2 ω = cµ 1 + ε . (21) 8E1 0 As the homogeneous system exhibits no dispersion, the group velocity and phase velocity are identical such that (0) cg = c. Therefore, the relative group velocity is easily obtained from the group velocity cg = dω/dµ. This ratio is precisely 3Γx 2 chomog = ε . (22) gr 8E1 0 The relative group velocity shift is chomog = 0.15 when evaluated with the system parameters provided in Table 1 and gr ε02 = 0.1. Thus, the relative group velocity shift calculated using a discrete wave-based perturbation approach in two dimensions is validated by a perturbation analysis for the corresponding continuous system. 14

Figures 4.6a and 4.6b reveal that near-optimal configurations that achieve complete band gaps follow similar trends. Complete band gaps do not exist for a/d > 0.5. The first complete band gap is located between the 3rd and 4th dispersion branches. The width of the band gap is given as ∆ωi ≡ min(ωi+1 ) − max(ωi ),

µ ∈ ΓXM,

(23)

while band gap shifts are defined as the difference between the linear and nonlinear band gap widths (1)

Band gap shift ≡ |∆ωi

(0)

− ∆ωi |.

(24)

The hardening nonlinearity tends to shift the third branch more than the 4th branch for optimal configurations, resulting in a net decrease in band gap width. Similarly, a softening nonlinearity increases the band gap by downshifting the 3rd dispersion branch. Black lines in Fig. 4.6a show results for a hardening nonlinearity, while the lowamplitude (linear) systems are shown with nearby gray lines for comparison.

0.7 Bandgap shift, [rad/s]

Bandgap width, [rad/s]

4 C

3

2 B 1 A 0

0

0.1

0.2 0.3 Layer fraction, a/d

0.4

0.6

0.4 B 0.3 0.2

0.5

(a) Linear and nonlinear band gap width

C

0.5

A 0

0.1

0.2 0.3 Layer fraction, a/d

0.4

0.5

(b) Nonlinear band gap shift

Figure 4.6: Complete band gap widths for material systems A–C are shown in (a) where black lines indicate the band gap width of the nonlinear system and nearby gray lines indicate the width of the corresponding linear system. The change in complete band gap width is shown in(b), corresponding to the difference between nonlinear and linear band gaps. Some material/geometry combinations may be effectively used to close the band gap at high amplitude. For instance, material system A corresponding to a geometry of a/d = 0.16 tends to exhibit a small band gap ∆ω3 ≈ 0.3 rad/s. For larger amplitudes, this band gap closes due to the frequency shift. The magnitude of this band gap shift is shown in Fig. 4.6b. Indeed, as the band gap shift is greater than the band gap width, the band gap closes. Both small layer fractions and high material contrasts tend to increase the band gap shift. Figure 4.7 depicts the dispersion band structure for the configuration corresponding to curves C with layer fraction a/d = 0.1. Constrained to the predefined geometry, this configuration represents the optimal bandgap width (∆3 = 2.87 rad s−1 ) and near-optimal bandgap shift (0.53 rad s−1 ). As in one-dimension, the maximum band gap width and maximum band gap shift for a given material system do not appear at the same layer fraction. However, both maxima appear for small layer fractions and thus this may be used to guide design of systems which are sensitive to the Bloch wave amplitude. While the heuristic design utilized for parametric analysis produced simple configurations, the design space was restricted to those systems corresponding to the predefined geometry defined by the ligament width a. This constraint is removed in the following section whereby the selected optimizer assigns material properties to each of the elements in the considered domain.

4.3

Genetic algorithm optimization for tuning dispersion shifts

In this section we seek to simultaneously optimize low-frequency band gap width and band gap shift (similar to Fig. 4.6). A genetic algorithm implementation offers a convenient and viable solution for strategically assigning material

15

Normalized Frequency, ωd/cA

15

10

5 L.A. T.A. 0

Γ

X

M Wavenumber

Γ

Figure 4.7: Dispersion diagram for unit cell with optimized bandgap width and near-optimal dispersion shift (a/d = 0.1, curve C) within the predefined configuration space. Solid lines indicate the nonlinear system while dotted lines indicate the low-amplitude (linear) system. to each of the element domains within the 1/8 symmetry constraint Ω(e) defined by g = single aggregate objective function (fitness function) defined according to F(g) = α1 F1 (g) + α2 F2 (g),

S

gi . We first introduce a (25)

where F1 (g) is a metric on the band gap width and F2 (g) is a metric on the band gap shift. The linear system band gap width appearing between branches 3 and 4 is calculated according to (0)

F1 (g) = ∆ω3 (µ1 , µ2 ),

(µ1 , µ2 ) ∈ ΓXM,

while the corresponding nonlinear band gap shift is obtained through the expression (0) F2 (g) = ∆ω3 (µ1 , µ2 ) − ∆ω3 (µ1 , µ2 ) , (µ1 , µ2 ) ∈ ΓXM,

(26)

(27)

where the dispersion relation is evaluated along the first Brillouin zone contour ΓXM. The absolute value signs for F2 (g) ensure that shifts which both open and close low-amplitude band gaps are equally awarded. Since almost any randomly generated design will not exhibit a band gap, the coefficients α1 and α2 are defined piecewise according to whether or not a band gap exists, similar to [16]. α1 = 10e4, α1 = 1,

α2 = 10e3 α2 = 0

if F1 (g) > 0 (band gap exists) . if F1 (g) ≤ 0 (no band gap)

(28)

This fitness function promotes individuals which exhibit large band gaps and, to a lesser extent, those which exhibit large band gap shifts. The piecewise jump in the α1 function ensures that individuals exhibiting band gaps are always selected over those individuals which exhibit no band gap. Until a configuration exhibiting a band gap is reached, the band gap shift function F2 has no meaning, and so α2 is set to zero for configurations with no band gap. The optimization problem for simultaneously maximizing band gap width and shift is expressed as ˜ max F(g) subject to: g ∈ Ω, g

16

(29)

˜ denotes the collection of sub-domains g within the 1/8 symmetry constraint. A priori knowledge suggests where Ω choosing a high contrast stiffness ratio E1B /E1A = 20 for the analysis. The algorithm begins by generating an initial population of 20 individuals. A uniform, random distribution of individual sequences (with all elements of g taking either 0 or 1) ensures that subsequent evolutionary iterations are not biased towards any particular design. Each individual of the population is evaluated according to the fitness function; individuals with greater fitness are selected to move on to the next generation (iteration), while those with lower fitness scores are discarded. New individuals enter into the population through both mutation and crossover. During mutation, small modifications (bit-flips from 0 to 1, or 1 to 0) to existing individuals are introduced through random selection. Crossover individuals spawn from combinations of two- or more parents. This implementation uses a single-point crossover scheme whereby two parents swap a randomly selected portion of their genes (sub-domains). The genetic algorithm terminates when the change in objective function is less than 1e-6 (typically about 1000 – 1500 iterations). Figure 4.8 illustrates the evolution of the fitness function F(g) as a function of generation number. The algorithm begins with a random distribution of material as seen in the subfigure depicting the initial stages of evolution; such designs contain no bandgaps. The transition from a no-bandgap-configuration to a bandgap-configuration is noted after approximately 40 generations where F(g) transitions from a modestly decreasing region to a more pronounced trend. 5

0

x 10

F(g)

−2

−4

−6

0

100

200

300

400

500 600 Generation

700

800

900

1000

Figure 4.8: Unit cell evolution as a function of generation number (the best value of F(g) is shown for each generation). The final design is reached after 1000 iterations. We note that the transmission from designs without bandgaps to those with bandgaps is clearly evident around generation 40 where a decreasing trend in the fitness function begins. Snapshots of the unit cell evolution included. Figure 4.9a depicts the final design surrounded by partial neighbor unit cells. The design consists of a square inclusion of rigid material surrounded by a thin nonlinear ligament of width a/d = 0.125. In contrast to the simple designs presented earlier, the optimized configurations contain rigid beams protruding from the corners of the unit cell. The appearance of these corner features may have been predicted, since previous (gradient-based) studies have shown that corner features tend to distinguish optimal designs from sub-optimal designs in linear band gap optimization problems [8]. A brief optimization study using the described genetic algorithm approach and linear material properties is provided in Appendix C as a validation study. The dispersion relation evaluated along the first Brillouin zone contour is shown in Fig. 4.10. Points P3 and P4 bound the complete band gap of the nonlinear system with width ∆ω3 = 3.191 rad s−1 . The band gap shift (0.43 rad s−1 ) results from nonlinear shifts at both of points P3 and P4 . However, the branch shift may be roughly approximated by the difference between frequencies at P2 and P3 . Practical implementation of the optimal design may be effected through slight modification of the final design produced by the genetic algorithm as a post-processing procedure. A similar post-processing procedure has been employed in bandgap-based optimization of Mindlin plates structures [20]. Figure 4.9b depicts one such design whereby a majority image filter has been applied. The resulting design has a slightly smaller band gap width at ∆3 = 3.046 rad s−1 , but a larger band gap shift (0.44 rad s−1 ). The modified design retains longitudinal and transverse stiffness at the inclusion corners and exhibits no qualitative difference in dispersion band structure. An alternative approach for generating designs with improved manufacturing ability utilizes a third term in the 17

(a) Unit cell before processing

(b) Unit cell after processing

(c) Unit with manufacturing constraints

Figure 4.9: Final unit cell design (enclosed by dashed line) surrounded by partial neighboring cells. (a) Depicts the raw output from the genetic algorithm, while (b) depicts a post-processed version amenable to fabrication. (c) Depicts the results from an optimization study which includes manufacturing constraints as part of the optimization process. White regions correspond to nonlinear material A, while black regions correspond to material B. fitness function which penalizes designs that contain elements with dissimilar von Neumann or Moore neighborhoods. The unit cell resulting from an optimization/design study using this approach are shown Fig. 4.9c. The band gap width ∆ω3 = 3.234 rad s−1 is an improvement over the results without any manufacturing constraints (Fig. 4.9a), while the band gap shift has decreased (0.396 rad s−1 ). Although successful in generating a reasonable unit cell, implementing these manufacturing constraints in a genetic algorithm optimization and design procedure has practical challenges since the unit cell evolves from a random material distribution. With a large number of elements, the likelihood of evolving finer structures (as in the corner features of Fig. 4.9c) into more manufacturable designs decreases. Additionally, more terms in the objective/fitness function increase the number of evolutions needed to reach a near-optimal design. For these reasons, it is preferable to leave manufacturing design as a post-processing operation.

4.4

Discussion of wave modes and their effect on dispersion shifts

An understanding of wave modes which exhibit maximum shifts may be used to speed up computation times by reducing the number of wave vectors evaluations, or may assist in design strategies for nonlinear phononic systems and metamaterials. An examination of the strain fields for Bloch wave modes at points P1 , P3 and P4 provides some intuition about where and why minimum and maximum shifts occur. The O(ε 0 ) strain fields εxx (x, y), εyy (x, y), and εxy (x, y) evaluated at these points are shown in Fig. 4.11. Points P1 and P4 exhibit negligible shifts relative to the shift at point P3 . The εxx strain field at point P3 is much stronger, relative to points P1 and P4 , while the εyy components are approximately the same order of magnitude. This localized tension and compression in the nonlinear ligaments is responsible for producing large stresses due to the Γx coefficient. This results in an effective stiffness increase, and therefore a shift in the dispersion curves. Similarly, points along the longitudinal acoustic branch, labeled L.A., receive shifts which produce effective group velocity variation. In contrast, points P1 and P4 exhibit greater shear strain εxy than the wave mode at P3 . As the nonlinear shear coefficient Γxy is zero, these modes exhibit no nonlinear shifts. In the same vein, those modes along the transverse acoustic branch (labeled T.A.) also exhibit no shift since they are in nearly pure shear. The reason for the rigid fingers protruding from the corner of the cell is less apparent; one possible explanation is that they may tend to compress/tension ligament material which would otherwise see less strain.

5

Conclusions

Nonlinear systems such as those examined herein exhibit amplitude-dependent dispersion which introduces new design opportunities for tunable devices such as filters and waveguides. A frequency shift sensitivity Kω quantifies the shift sensitivity for a particular configuration by isolating system-dependencies from amplitude. The nonlinear force

18

Normalized Frequency, ωd/cA

15

P4 10

P3

P1 P2 5 L.A. T.A. 0

Γ

X

M

Γ

Figure 4.10: Dispersion diagram for optimized unit cell with manufacturing improvements applied as post-processing operations (Fig. 4.9b). Solid lines indicate the nonlinear system while dotted lines indicate the low-amplitude (linear) system. vector presents the biggest technical hurdle because its functional time-dependence is required to obtain c1 . However, techniques such as the one presented herein offer simple solutions to this problem. We have shown that highly tunable nonlinear systems may be realized by arranging rigid inclusions in a compliant matrix of nonlinear material. In both 1D and 2D settings, arrangements with thin nonlinear layers resulted in highly sensitive group velocity and band gap width which may be practical in many applications (e.g., diodes, filters, and waveguides). A genetic algorithm topology optimization routine confirms that such designs are, indeed, optimal for achieving high sensitivity. The genetic algorithm implemented herein tended toward a square inclusion with nominal dimensions a/d = 0.125 and small fingers protruding from the corners. A heuristic analysis of critical wave modes provides some insight into why some dispersion branch points result in large shifts, while others result in no shift. It is shown that high strain regions in the linear Bloch wave modes are evidence of such sensitivity, and may be used to enhance design strategies for nonlinear phononic crystals and metamaterial devices.

Acknowledgement The authors would like to thank Dr. Eduardo A. Misawa and the National Science Foundation for supporting this research under Grant No. (CMMI 0926776).

References [1] J. Joannopoulos, P. Villeneuve, and S. Fan, “Photonic crystals,” Solid State Communications, vol. 102, no. 2-3, pp. 165–173, 1997. [2] R. H. Olsson and I. El-Kady, “Microfabricated phononic crystal devices and applications,” Measurement Science and Technology, vol. 20, p. 012002, November 2009.

19

0.2 0.1 0 −0.1 −0.2 εxx

εyy

εxy

(a) Strain fields at P1

0.2 0.1 0 −0.1 −0.2 εxx

εyy

εxy

(b) Strain fields at P3 (and P2)

0.2 0.1 0 −0.1 −0.2 εxx

εyy

εxy

(c) Strain fields at P4

Figure 4.11: Mode shape deformation and strain at P1, P3 (P2), and P4. [3] K. Bertoldi and M. C. Boyce, “Mechanically triggered transformations of phononic band gap elastomeric structures.,” Physical Review B 77, pp. 1–10, February 2008. [4] B. Liang, B. Yuan, and J.-c. Cheng, “Acoustic diode: Rectification of acoustic energy flux in one-dimensional systems,” Phys. Rev. Lett., vol. 103, p. 104301, Sep 2009. [5] R. H. Olsson III, I. F. El-Kady, M. F. Su, M. R. Tuck, and J. G. Fleming, “Microfabricated VHF acoustic crystals and waveguides,” Sensors and Actuators A, pp. 87–93, 2008. [6] K. Manktelow, R. K. Narisetti, M. J. Leamy, and M. Ruzzene, “Finite-element based perturbation analysis of wave propagation in nonlinear periodic structures,” Mechanical Systems and Signal Processing, pp. –, 2012. [7] J. Jang, C. K. U. Ullal, T. Gorishnyy, V. V. Tsukruk, and E. L. Thomas, “Mechanically tunable three-dimensional elastomeric network/air structures via interference lithography,” Nano Letters, vol. 6, no. 4, pp. 740–743, 2006. [8] O. Sigmund and J. S. Jensen, “Systematic design of phononic band-gap materials and structures by topology optimization,” Philosophical Transactions: Mathematical, Physical, and Engineering Sciences, vol. 361, pp. 1001– 1019, May 2003. 20

[9] Y. Yun, G. Miao, P. Zhang, K. Huang, and R. Wei, “Nonlinear acoustic wave propagating in one-dimensional layered system,” Physics Letters A, vol. 343, no. 5, pp. 351 – 358, 2005. [10] S. Halkjær, O. Sigmund, and J. S. Jensen, “Inverse design of phononic crystals by topology optimization,” Zeitschrift fur Kristallographie, vol. 220, no. 9–10, pp. 895–905, 2005. [11] J. S. Jensen, “Topology optimization problems for reflection and dissipation of elastic waves,” Journal of Sound and Vibration, vol. 301, no. 12, pp. 319 – 340, 2007. [12] C. Le, T. E. Bruns, and D. A. Tortorelli, “Material microstructure optimization for linear elastodynamic energy wave management,” Journal of the Mechanics and Physics of Solids, vol. 60, no. 2, pp. 351 – 378, 2012. [13] J. S. Jensen, “Space-time topology optimization for one-dimensional wave propagation,” Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 58, pp. 705 – 715, 2009. [14] M. Hussein, K. Hamza, G. Hulbert, R. Scott, and K. Saitou, “Multiobjective evolutionary optimization of periodic layered materials for desired wave dispersion characteristics,” Structural and Multidisciplinary Optimization, vol. 31, no. 1, pp. 60–75, 2006. [15] O. Bilal, M. El-Beltagy, M. Rasmy, and M. Hussein, “The effect of symmetry on the optimal design of twodimensional periodic materials,” in Informatics and Systems (INFOS), 2010 The 7th International Conference on, pp. 1 –7, march 2010. [16] O. Bilal and M. Hussein, “Ultrawide phononic band gap for combined in-plane and out-of-plane waves,” Physical Review E, vol. 84, no. 6, p. 065701, 2011. [17] M. I. Hussein, G. M. Hulbert, and R. A. Scott, “Dispersive elastodynamics of 1d banded materials and structures: Design,” Journal of Sound and Vibration, vol. 307, no. 35, pp. 865 – 893, 2007. [18] H. Meng, J. Wen, H. Zhao, and X. Wen, “Optimization of locally resonant acoustic metamaterials on underwater sound absorption characteristics,” Journal of Sound and Vibration, vol. 331, no. 20, pp. 4406 – 4416, 2012. [19] R. A. Wildman and G. A. Gazonas, “Genetic programming-based phononic bandgap structure design,” Army Research Laboratory, pp. 1 – 36, September 2011. [20] S. Halkjær, O. Sigmund, and J. S. Jensen, “Maximizing band gaps in plate structures,” Structural and Multidisciplinary Optimization, vol. 32, no. 4, pp. 263–275, 2006. [21] J. S. Jensen, O. Sigmund, L. H. Frandsen, P. I. Borel, A. Harpoth, and M. Kristensen, “Topology design and fabrication of an efficient double 90 photonic crystal waveguide bend,” Photonics Technology Letters, IEEE, vol. 17, no. 6, pp. 1202–1204, 2005. [22] A. A. Larsen, B. Laksafoss, J. S. Jensen, and O. Sigmund, “Topological material layout in plates for vibration suppression and wave propagation control,” Structural and Multidisciplinary Optimization, vol. 37, no. 6, pp. 585–594, 2009. [23] R. K. Narisetti, M. Ruzzene, and M. J. Leamy, “A perturbation approach for analyzing dispersion and group velocities in two-dimensional nonlinear periodic lattices,” Journal of Vibration and Acoustics, vol. 133, no. 6, p. 061020, 2011. [24] W. C. Elmore and M. A. Heald, Physics of waves. New-York: Dover, 1985. [25] M. J. Leamy and O. Gottlieb, “Internal resonances in whirling strings involving longitudinal dynamics and material non-linearities,” Journal of Sound and Vibration, vol. 236, no. 4, pp. 683 – 703, 2000. [26] J. J. Rushchitsky and C. Cattani, “Evolution equations for plane cubically nonlinear elastic waves,” International Applied Mechanics, vol. 40, pp. 70–76, January 2004. [27] J. J. Rushchitsky, “Interaction of waves in solid mixtures,” Applied Mechanics Reviews, vol. 52, no. 2, pp. 35–74, 1999.

21

[28] J. Reddy, An introduction to the finite element method, vol. 2. McGraw-Hill New York, 2nd ed., 1993. [29] T. Hughes, The finite element method: linear static and dynamic finite element analysis, vol. 682. Dover Publications New York, 2000. [30] R. K. Narisetti, M. Ruzzene, and M. J. Leamy, “Wave propagation in membrane-based nonlinear periodic structures,” Proceedings of ASME 2011 IDETC, August 28-31 2011. [31] D. S. Bethune, “Optical harmonic generation and mixing in multilayer media: analysis using optical transfer matrix techniques,” J. Opt. Soc. Am. B, vol. 6, pp. 910–916, May 1989. [32] P. Luan and Z. Ye, “Acoustic wave propagation in a one-dimensional layered system,” Phys. Rev. E, vol. 63, p. 066611, May 2001. [33] M. Shen and W. Cao, “Acoustic bandgap formation in a periodic structure with multilayer unit cells,” Journal of Physics D: Applied Physics, vol. 33, p. 1150, 2000. [34] K. Manktelow, M. J. Leamy, and M. Ruzzene, “Comparison of asymptotic and transfer matrix approaches for evaluating intensity-dependent dispersion in nonlinear photonic and phononic crystals,” Wave Motion, vol. 50, no. 3, pp. 494 – 508, 2013. [35] N. Panoiu, M. Bahl, R. Osgood Jr, et al., “Optically tunable superprism effect in nonlinear photonic crystals,” Optics letters, vol. 28, no. 24, pp. 2503–2505, 2003. [36] M. Soljaˇci´c, S. Johnson, S. Fan, M. Ibanescu, E. Ippen, and J. Joannopoulos, “Photonic-crystal slow-light enhancement of nonlinear phase sensitivity,” JOSA B, vol. 19, no. 9, pp. 2052–2059, 2002.

A

Nonlinear dispersion analysis review

We consider materials comprised of repeating unit cells as shown in Fig. A.1. The governing equations describing dynamics of a unit cell and its neighbors post-discretization are Mu¨ + Ku + εfNL (u) = 0,

(30)

where |ε|  1 is a small parameter, and M and K denote assembled mass and stiffness matrices. The nonlinear force vector fNL arises from a nonlinear force-displacement or stress-strain relationship, and is therefore a nonlinear function of the displacement vector u. As the nonlinear restoring force fNL contributes weakly to the overall system response, only the quadratic and cubic terms from its Taylor expansion are retained. Quadratic terms do not contribute frequency shifts according to first-order perturbation theory; thus, the remainder of the discussion will consider the effect of cubic terms on a system’s dispersion properties.

Figure A.1: Central unit cell (colored red) surrounded by neighboring unit cells for 1D, 2D, and 3D periodic systems

22

We note that Eq. (30) contains enough information to truly represent wave propagation in an infinite system since all internal interactions on the central unit cell are represented. This is not true for the surrounding unit cells where internal forces are unknown. The perturbation approach described next resolves this issue. Note that unlike linear systems where only a single unit cell is sufficient, neighbors must be included in order to properly identify nonlinear forces. The first step in the nonlinear dispersion analysis (NDA) is the introduction of dimensionless time τ = ωt and standard asymptotic expansions ω = ω (0) + εω (1) + O(ε 2 ) u = u(0) + εu(1) + O(ε 2 ).

(31)

We also expand the nonlinear force vector fNL in a Fourier series as +∞

fNL =

∑

cn exp (inτ).

(32)

n=−∞

The expansions (31)–(32) are substituted into the equation of motion (30) to yield the ordered perturbation equations O(ε 0 ) : O(ε 1 ) :

∂ 2 u(0) + Ku(0) = 0 ∂ τ2 +∞ ∂ 2 u(0) ∂ 2 u(1) (1) (0) (1) + Ku = −2ω ω M − (ω (0) )2 M ∑ cn exp (inτ). ∂ τ2 ∂ τ2 n=−∞ (ω (0) )2 M

(33) (34)

The Bloch theorem applied to the O(ε 0 ) equations produces a set of equations describing an infinite [linear] system that depends on the Bloch wavevector µ 2 (0) ˆ ∂ uˆ + K ˆ uˆ (0) = 0, ω02 M (35) ∂ τ2 ˆ and K ˆ describing the central unit where uˆ (0) holds only the displacements of the central unit cell and the matrices M cell may, in general, be functions of the Bloch wavevector µ . If we seek solutions with time-harmonic displacement uˆ (0) =

A φ exp(iτ) + c.c., 2 j

(36)

(0) µ ). The disthen Eq. (35) yields a standard eigenvalue problem for eigenvalues λ j = (ω j )2 and eigenvectors φ j (µ (0) µ ), where the jth eigenvalue corresponds to the persion relation for the continuous system is approximated by ω j (µ jth branch on the dispersion curve. Removal of secular terms from the O(ε 1 ) ordered equation leads to a frequency correction term for the jth branch. This salient result underpinning NDA describes the amplitude/intensity variation of the dispersion relationship due to nonlinearity, (1)

ωj =

φ Hj c1 (0) ˆ φj Aω j φ Hj Mφ

.

(37)

Equation (37) has shown good agreement with numerical simulations for |ε| < 0.25 [34]. However, as amplitudes or nonlinear coefficients increase, more harmonic generation and other nonlinear effects should be expected to compete with this perturbation analysis and the actual dispersion shift may vary.

B

Dispersion shift for a homogeneous rod

The dispersion shift arising from a one-dimensional rod with a constitutive nonlinearity may be obtained through perturbation of the continuous PDE governing wave propagation ∂σ ∂ 2u =ρ 2. ∂x ∂t

23

(38)

The constitutive law governing the strain displacement relationship is defined as  3 ∂u ∂u σ = E1 . + E3 ∂x ∂x

(39)

Substitution of Eq. (39) into Eq. (38) results in a nonlinear wave equation. The small parameter ε is introduced through the relationship E3 = ε Eˆ3 to enforce weak nonlinearity, where Eˆ3 is of the same order as the linear modulus E1 .  3 ∂ ∂ 2 u(x,t) ∂u ∂ 2 u(x,t) ˆ3 (x) + ε E . (40) E1 (x) = ρ(x) ∂ x2 ∂x ∂x ∂t 2 Nondimensional time τ = ωt and asymptotic expansions are introduced according to ω = ω (0) + εω (1) + O(ε 2 ) u = u(0) + εu(1) + O(ε 2 ).

(41)

Substitution of the expansions Eq. (41) into the weakly nonlinear wave equation (40) produces a set of equations ordered by the small parameter ε O(ε 0 ) : O(ε 0 ) :

2 (0) ∂ 2 u(0) (0) 2 ∂ u = ρ(ω ) ∂ x2 ∂ τ2  (0) 2 2 (0) 2 (1) ∂ u ∂u ∂ u ∂ 2 u(1) ∂ 2 u(0) E1 + 3Eˆ3 = ρ(ω (0) )2 + 2ρω (0) ω (1) . 2 2 2 ∂x ∂x ∂x ∂τ ∂ τ2

E1

(42) (43)

The solution to the O(ε 0 ) is well-known, and given by the plane wave A exp i(kx − τ) + c.c = A cos(kx − τ). (44) 2 p where the wavenumber k is defined through the relationship ω (0) = ck, c = E1 /ρ denotes phase speed, and c.c. denotes complex conjugate terms. Subsequent substitution of u(0) (x,t) into the O(ε) equation results in the equation u(0) (x,t) =

E1

2 (1) ∂ 2 u(1) (0) 2 ∂ u − ρ(ω ) = −2ω (0) ω (1) ρA cos(kx − τ)+ ∂ x2 ∂ τ2 3ˆ 3 4 E3 A k cos(kx − τ) + O(3kx − 3τ), 4

(45)

where terms higher frequency terms are indicated by O(3kx − 3τ). The linear kernel of the O(ε 1 ) equation is identical to the linear kernel of the O(ε 0 ) equation. Removing secular (those multiplying cos (kx − τ)) results in an equation which may be solved for a first-order frequency correction ω (1) 3 −2ω (0) ω (1) ρA + Eˆ3 A3 k4 = 0. 4

(46)

The solution for the first frequency correction is ω (1) =

3 Eˆ3 A2 k4 . 8 ρω (0)

(47)

The strain amplitude e0 = kA and the relationship ω (0) = ck may be used to obtain alternate forms of the frequency correction, namely 3 Eˆ3 2 (0) ω (1) = e ω . (48) 8 E1 0 Thus, the first-order corrected frequency is given by     3 Eˆ3 2 (0) 3E3 2 ω = ω (0) + ε e0 ω = ω (0) 1 + e0 . (49) 8 E1 8E1 This relationship is useful for validating results produced from nonlinear dispersion analysis routines. 24

C

Band gap optimization results for linear materials

The genetic algorithm implemented as a tool to search the design space was tested with linear material properties to validate its ability to generate near-optimal unit cells. The material properties of [8] were used in order to generate a comparison between the results presented therein and optimized unit cells produced from the genetic algorithm. This good qualitative and quantitative agreement suggests the genetic algorithm implementation can yield near-opti-

Figure C.1: Unit cell optimized for a complete first bandgap with linear material properties. mal configurations for the chosen material properties and finite-element discretization. Finer discretization levels may require the development of a gradient-based approach to handle the larger number of design variables. The unit cell dimension used was different; however this scales the frequency axis and relative values are comparable to results reported in literature. An optimized unit cell for linear material properties as well as its dispersion diagram is depicted in Fig. C.1. The relative bandgap width (total width divided by center frequency) is approximately 0.20. These results are in good agreement with those presented in [8] where the relative bandgap width as derived from an efficient gradient-based approach is reported as 0.21. The small difference may be atributed to finer finite-element discretization (48x48 elements has been used in Fig. C.1) as well as finer evaluation of the dispersion diagram.

25