AIAA 2016-3830 AIAA Aviation 13-17 June 2016, Washington, D.C. 17th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference

Shape Optimization of Acoustic Metamaterials and Phononic Crystals with a Time-Dependent Adjoint Formulation: Extension to Three-Dimensions

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

Weiyang Lin1, James C. Newman III2, W. Kyle Anderson 3 and Xueying Zhang4 SimCenter: Center of Excellence in Applied Computational Science and Engineering University of Tennessee, Chattanooga, TN 37403

A time-dependent adjoint approach for obtaining sensitivity derivatives for shape optimizations of three-dimensional acoustic metamaterials and phononic crystals is presented. The acoustic wave propagation problem is solved in the time-domain using a stabilized finite element formulation. The gradient-based design procedure is suitable for large numbers of design variables, and results are shown for achieving effective material properties with a unit cell, and the broadband noise reduction with periodic arrays of stainless-steel cylinders.

I. Introduction

A

COUSTIC metamaterials are engineered materials that have acoustic properties that cannot easily be found in nature. Such properties include negative effective bulk modulus and mass density, and they can be applied to soundproof, acoustic cloaking, seismic shields, and so on. Metamaterials are dedicated to sub-wavelength structures, and when the size of the design units increases, the target material behaves more like a phononic crystal where the macro arrangement becomes more important. Phononic crystals are periodic arrays of inclusions embedded in a host matrix. This kind of composite media is frequently designed to achieve the “phononic band-gap” which can be employed to prevent acoustic/elastic waves in certain frequencies from propagating. A considerable amount of work has been accomplished for the designs of phononic crystals, where topology optimization has been widely used. Topology optimization treats geometries as if they are in a discrete domain of a Cartesian grid, and performs selections and de-selections of points based on the objective functional [1]. This methodology is well suited for the overall design of materials and structures. While proven successful, topology optimization is usually limited only to a conceptual level and typically requires adjustments to the final designs. Shape optimization is often employed along with topology optimizations for this purpose. Similar in both cases, gradient-based optimizers are used to minimize a specified objective function. While several techniques exist for evaluating the sensitivities derivatives for multidisciplinary simulations, the adjoint methods are of particular interest in this paper. For practical design, controls over various parts of the geometries are typically required, leading to a large number of design variables. The computational costs associated with an adjoint formulation scales with the number of cost function sensitivity derivatives and not the number of design variables. Therefore, the adjoint formulation is better suited for gradient-based design optimizations procedures that have large number of independent variables. However, for time-dependent sensitivity analysis, the adjoint formulation incurs additionally costs associated with the time integration, though it can be relieved by approximations [2, 3]. In this paper, the shape optimization for phononic crystals with time-dependent sensitivity analysis is extended to three dimensional applications from previous work [2]. To achieve high-order accurate solutions and designs on continuous geometries, the Streamline Upwind/Petrov Galerkin time-domain method is utilized. Sensitivity analysis is conducted in order to perform the shape optimization, and an adjoint-based formulation is enabled when necessary. The optimization procedure is applied to cylindrical periodic structures for broadband noise reduction, and a unit cell to achieve effective material properties. 1

Ph.D. Candidate, Graduate School of Computational Engineering, 701 E. M. L. King Blvd. Student Member. Professor, Graduate School of Computational Engineering, 701 E. M. L. King Blvd. Senior Member. 3 Professor (retired), Graduate School of Computational Engineering, 701 E. M. L. King Blvd. Associate Fellow. 4 Ph.D. Candidate, Graduate School of Computational Engineering, 701 E. M. L. King Blvd. Student Member. 1 American Institute of Aeronautics and Astronautics 2

Copyright © 2016 by Weiyang Lin. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

II. Stabilized Finite Element Method for Acoustic Metamaterials and Phononic Crystals The simulation in this paper is performed using a stabilized finite element formulation: Streamlined Upwind/Petrov Galerkin (SUPG) method for time-domain applications. Finite-difference time-domain (FDTD) methods have been widely adopted for solving different kinds of acoustic problems. As the speed of computing resources increases and as the geometry conforming becomes more of a concern, finite element methods (FEM) have become an alternative. While traditional finite-element time-domain (FETD) methods are used to solve for the second-order hyperbolic partial differential equations (PDEs), they have the disadvantage of resolving secondary variables with one order less accuracy than the primary variables. This problem can be overcome by a formulation based on the first-order governing PDEs recast from their second-order counterparts. A technique referred to as discontinuous Galerkin (DG) has recently been extended to acoustic and elastic problems based on the first-order PDEs. Both the SUPG and DG are the combination of finite volume methods (FVM) and FEM. They both have an origin in fluid mechanics, and are readily extensible to nonlinear problems. It has been shown that SUPG has the advantage of reduced number of unknowns in comparison with DG. In three-dimensions, the general acoustic wave propagation in a heterogeneous medium may be described by the equations  u v w  p  K0    (1) 0 t  x y z 

u 1 p  0 t 0 x

(2)

v 1 p  0 t 0 y

(3)

w 1 p  0 (4) t 0 z where p is the pressure , u, v and w are the velocity components, K0 is the bulk modulus of compressibility, ρ0 is the density. The speed of sound may be evaluated as c0  K0 0 . The governing equations in both cases can be written in a compact divergence form as

Q Q Q Q A B C 0 (5) t x y z where Q is the primitive variables, and A, B and C are functions of the material properties. While the wave equations can be expressed in a frequency-domain formulation for time-harmonic sources, the time-domain solution is of interest in this paper for its capability of solving broadband problems. The Petrov-Galerkin methods are formulated with a weighted residual method, which can be cast in the form of  Q Q Q Q  (6)     t  A x  B y  C z    0 where  is a weighting function. In the present work, the Streamlined Upwind Petrov-Galerkin (SUPG) method is used in defining the weighting function [4]  N  N N (7)    N  I     A   B  C     y z  x  where   is the stabilization matrix, and it can be obtained using the following definitions [5]

 

1

n

 k 1

N k N N  A  k  B  k C  x y z

N k N N 1  A  k  B  k C   T     T  x y z

(8) (9)

Here, T  and    are the right eigenvectors and eigenvalues, respectively, of the matrix on the left side. The inversion of local stabilization matrices   can be calculated using Gauss eliminations. While the quadrature rules for integrating the volume integrals are used for both volume and surface integrals, the surface integrals in the interior of the domain do not need to be evaluated due to the use of a Galerkin 2 American Institute of Aeronautics and Astronautics

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

formulation. On the boundaries of the physical domain, the appropriate boundary conditions are strongly enforced by incorporating them into the surface integral. At the interfaces of different materials, a Godunov's flux formulation is used to solve for the Riemann problem [6]. The temporal discretization of the governing equation is advanced by implementing a backward differentiation formula (BDF), and the solution is obtained with an implicit time marching approach. At each time step, the generalized minimal residual (GMRES) linear system equation solver [7] is utilized with an incomplete lower-upper (ILU) preconditioner to solve for the non-diagonally dominant system as a result of the finite element formulation. The proposed finite element solver is implemented for a distributed memory space using message passing interface (MPI). The solution procedure starts by decomposing a given unstructured mesh with the METIS library [8], and proceeds by conducting computation on each of the subdomains with overlapping nodes for communications. The partitioning of computational meshes is efficiently carried out by METIS, such that the surface-to-volume ratio is minimized, leading to an optimized parallel computation. Figure 1 shows a sample of a decomposed computational mesh using METIS, where each subdomain is represented with a different color. For verification, Fig. 2 shows a case of acoustic wave propagation initiated by a pulse in the -x direction of the sphere. The instantaneous pressure fields on a sphere and the X-plane are plotted, and the record of solution in nondimensional units at a sensor in the +x direction away from the sphere is shown in Fig. 3. The orders of accuracy of the proposed finite element solver are given in Fig. 4, and the slopes of the log-log plot indicate that the second order of accuracy is achieved by the linear elements (slope = 2.09), and the third order of accuracy (slope = 3.15) is approached with quadratic elements.

Figure 1. An example of the domain decomposition using METIS.

Figure 2. An instantaneous pressure field of a pulse hitting a sphere by the SUPG solver.

Figure 3. The time-domain solution of pressure at a sensor behind the sphere.

Figure 4. Orders of accuracy of the SUPG solver with linear and quadratic elements.

3 American Institute of Aeronautics and Astronautics

III. Time-Dependent Sensitivity Analysis and Adjoint Formulation

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

In gradient-based optimization frameworks, a functional I to be minimized is typically described with a set of design variables β. Depending on the applications, the cost function can be the transmitted or reflected pressure at a specified location, pressure at a given frequency [9], or a combination of different physical quantities. The sensitivity derivatives of the functional (i.e. the cost function) with respect to the design variables are utilized to evaluate an appropriate search direction for improving the design. A. Time-Dependent Sensitivity Analysis The sensitivity derivatives, or the gradients of the cost function with respect to the design variables, may be calculated in many different ways. The most straightforward method may be the finite-difference method, where a central-difference is given by dI I       I        O   2  (10) d 2 However, this method suffers from inaccuracy because of step-size problems and is usually prohibitively expensive. Another way to approximate the derivatives is through the use of complex Taylor series expansion (CTSE). This method perturbs the design variables in the complex part, resulting in higher-accuracy than finite-difference methods because of the elimination of subtractive cancellation error [10, 11] dI Im I     i  (11)   O   2  d  However, similar to the finite-difference methods, CTSE is generally too expensive to have practical implementations because of the requirement of multiple function evaluations. Alternatively, the sensitivity derivative can be computed using a forward mode direct differentiation by examining the functional dependencies of the cost function. Since the cost function is defined by a functional with dependencies on the design variables β, computational mesh X, and the solution quantities Q, i.e., I  I   , X , Q (12) the total differential of I with respect to β can be expressed by dI I I X I Q    (13) d  X  Q  The residual of the governing equation for a steady problem may be expressed as R   , X , Q  0 (14) the total differential of R with respect to β is given by dR R R X R Q    0 (15) d  X  Q  In some applications, the residual is not an explicit function of β, thus Eq. (14) may be rewritten in the form of the solution sensitivity 1

 R   R X  Q        Q   X   Substituting back into Eq. (13), the sensitivity derivative becomes

(16)

1

dI I I X I  R   R X  (17)        d  X  Q  Q   X   In time-dependent problems such as transient acoustics, the solutions from the previous time-steps have to be considered in the calculation of the residual of the governing equation. For a backward differentiation formula, BDF2 in this case, the total differential of R with respect to β at time-step k is expanded to Rk   , X , Qk , Qk 1 , Qk  2   0 (18) and therefore the total differential of R k with respect to β becomes dR k R k R k X R k Q k R k Q k 1 R k Q k  2      0 d  X  Q k  Q k 1  Q k  2  Again if the residual is not an explicit function of β, Eq. (19) may be rewritten as 4 American Institute of Aeronautics and Astronautics

(19)

1

 R k   R k X Q k R k Q k 1 R k Q k  2  (20)   k    k 1  k 2      Q  Q   X  Q Substituting it back into Eq. (13) results in the sensitivity derivative for time-dependent problems using the forward mode 1

 R k   R k X R k Q k 1 R k Q k  2  (21)  k 1  k 2   k     Q  Q   X  Q As seen from Eq. (19), the computational costs for the forward mode scale with the number of design variables. dI I I X ncyc I    d  X  k 1 Q k

B. Discrete Adjoint Formulations As the number of design variables increases, the implementation of an adjoint method for computing sensitivity derivatives becomes the more efficient formulation. Examining Eq. (17), it is found that the additional computational overhead is due to the repetitive calculations of the solution sensitivities in Eq. (16). The adjoint methodology eliminates this overhead by transposing the inverse of the Jacobian matrix Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

T

T T dI I I X   R   I    R X        d  X    Q   Q    X   and by defining the adjoint variable T

(22)

T

 R   I  (23) Q        Q   Q  Utilizing Eq. (23), the resulting final form of the adjoint equation for steady-state problems becomes  R X  dI I I X    QT  (24)  d  X   X   The adjoint-based formulation is more complex for time-dependent problems. Similar to the steady-state case, it is derived by transposing the inverse of the Jacobian matrices. Examine the third term in Eq. (21),

I  k  k 1 Q ncyc

1

 R k   R k X R k Q k 1 R k Q k  2   k 1  k 2   k     Q  Q   X  Q

k ncyc T  R X R k Q k 1 R k Q k  2    Qk      k 1  Q k  2   k 1  X  Q k k 1 k 2 ncyc T R X T Q T Q      Qk    1k    2k   X     k 1 

(25)

where

 R k    k  Q  k Q

T

 I   k  Q 

 R k 

T

(26)

T

 1k   k 1  Qk  Q   R k 

(27)

T

 2k   k  2  Qk  Q 

(28)

It is observed from Eq. (25) that the evaluation of each term involves the solution sensitivities from the earlier two time-steps, which are unfortunately not readily available. To overcome this problem, the adjoint variables of “newer” time-steps can be regrouped with the ones of “older” time-steps; that is, the adjoint variable is reformulated as T T T T   R k    I  Qk   k    k    1k 1    2k  2   (29)   Q    Q   The basic algorithm can thus be written as [2]: 5 American Institute of Aeronautics and Astronautics

Algorithm.

A discrete adjoint formulation for time-dependent sensitivity derivatives

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

(1) Set  1k 1 ,  2k 1 and  2k  2 to be zero. Set k to be ncyc (2) Solve Eq. (29) for the adjoint variable. (3) Set the sensitivity derivatives by k T R X dI dI I X   Qk   d d X  X  (4) Set k  k  1 . (5) Set  2k  2   2k 1 , solve Eqs. (26-27) for  1k 1 and  2k 1 . (6) If k  1 , stop; otherwise go to step 2.

(30)

The adjoint formulation requires the storage of sensitivity matrices for all time steps. For large-scale problems, this becomes prohibitive. The storage problem can be mitigated by an approximate formulation of Algorithm 1. By dividing the global time into several intervals, local-in-time sensitivities can be calculated, and the sum of the local sensitivities is found to be an approximation to the global sensitivities [2, 3]. The mesh sensitivities generally require the sensitivity derivatives of the surface mesh points with respect to the design variables and may be expressed as X 1 X B  K  (31)   where XB is the deformation on the surfaces and [K] is the stiffness matrix for solving the resulting mesh movement X as the deformations are propagated into the interior using linear elasticity. When the mesh movement derivatives of Eq. (31) are substituted into the sensitivity derivatives Eq. (24), it is realized that the calculation may be further reduced by introducing another adjoint variable [12]. For example, the steady-state adjoint formulation given in Eq. (24) can now be re-written as dI I I X R 1 X B    QT (32) K  d  X  X  The last term in Eq. (38) may be written as X B R 1 X B QT  XT K  (33) X   where T  T   R  X   K     Q  (34)   X     is the adjoint variable for mesh sensitivity. This may be referred to as a double-adjoint method, which makes the sensitivity derivatives fully independent on the number of design variables (ndv). The double-adjoint formulation, however, may not significantly reduce the computational expense in the context of time-dependent problems. By introducing the adjoint variable of mesh sensitivity to Eq. (29), the equation now becomes T X dI dI I X B   Xk   (35) d d    X  While Eq. (35) leads to a formulation independent of ndv as well, it requires ncyc (number of time steps) linear system solutions in order to obtain the adjoint variables  Xk . Therefore, the double-adjoint formulation is only an improvement for time-dependent problems when ncyc is smaller than ndv.

IV. Results In periodic structures of units at the level of wavelength, the phononic band-gap problem is of primary interest. While the basic principles of wave propagation in acoustic/elastic media are well established [13], the industrial applications of phononic band-gap materials is still a demanding computer-aided design procedure. For units at subwavelength levels, the medium behaves approximately as if it is piecewise homogeneous. In this section, the analyses and designs for phononic crystals and acoustic metamaterials are presented. Baseline cases are examined to establish the appropriate cost functions. 6 American Institute of Aeronautics and Astronautics

A. Acoustic Metamaterials A slab of sub-wavelength scale is considered as the baseline case for acoustic metamaterials. The unit cell is a square with edge length a = 1 cm. The simulation is run in the frequency range from 0 kHz to 4 kHz, periodic boundary conditions are applied to the transverse directions of the domain. The inclusion is a cylinder with material properties as  0 in  0.7  0 out and  K0 in  2  K0 out . Although the elastic wave equations could have been used for the simulation, the acoustic wave solver is used because the overall behavior of the composite in this case can be considered as that of a fluid instead of a solid. The transmission and reflection coefficients are calculated by transforming the time-domain to frequency domain, which are in turn used to retrieve the effective refractive index and acoustic impedance [14]  1  R2  T 2  cos 1   (36) 2T   n kd

1  R   T 2  2 1  R   T 2 2

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

Z

2

(37)

where k is the wave number, and d is the effective slab thickness. The surface is parameterized using a modified Hicks-Henne bump function [15] which is extended to twodimensional surfaces







bi  sin 4  xmi sin 4  y ni with

mi  ni 



ln  0.5

 

ln xM i

ln  0.5

 

ln yNi

(38)

(39)

(40)

where x and y are the coordinates on the surface, xM i and y Ni are pre-selected values corresponding to the locations of the maxima. To ensure the clustering at the end points, the locations of xM i and y Ni are selected appropriately. For example, selecting

xMi  0.5  1  cos i 

(41)

where i    i  N  1 , and N is the number of control points achieves the desired clustering. However, xM i can also be chosen to be the same as  i such that the control points are equally spaced. Figure 5 shows the shapes of the modified Hicks-Henne functions for a given domain. The weights of these shape functions bi are considered the design variables, which yields the calculation of the displacement in the direction of wave propagation as ndv

zn   i  bi  xn , yn 

(42)

i 1

The deformation with randomly chosen design variables is shown in Fig. 6 for illustrative purposes. The derivatives of the parameterization are given by

zn  bi  xn , yn  i

(43)

An optimization procedure is conducted to achieve the desired material properties. The optimization seeks to maximize the refractive index as well as minimizing impedance as 2 2 2 min I    Z  Z *    n  n*    d   1  (44) k s.t. R  0  L    U 7 American Institute of Aeronautics and Astronautics

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

with 1 = 0.5 kHz and 2 = 2 kHz, which corresponds to the wavelengths approximately equal to 9 to 69 unit cell sizes. The target values are Z* = 1.0 and n* = 0.91. An optimal shape is achieved after several design iterations, and the surface deformation is shown in Figs. 7 and 8. It is observed in Figs. 9 and 10 that with the optimal shape, the values of both effective refractive index and impedance are closer to the desired values than the original flat slab.

Figure 5. Modified Hicks-Henne bump functions for surface deformations.

Figure 6. Sample illustration of a surface deformation with the modified Hicks-Henne bump functions and randomly generated design variables.

Figure 7. Optimized surface of the metamaterial to achieve desired effective material properties.

Figure 8. A view of the metamaterial to achieve desired effective material properties.

Figure 9. The initial, target and final effective refractive index of the metamaterial.

Figure 10. The initial, target and final effective impedance of the metamaterial.

8 American Institute of Aeronautics and Astronautics

B. Phononic Crystals The analysis of the response of composite materials typically starts with calculating the band-structures of a unit

 

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

cell. Based on the Floquet-Bloch wave theory, the solutions to the periodic eigen-problems are eigenvalues  k as continuous functions of Bloch wave vector k , forming discrete bands. Furthermore, due to the symmetry of the unit cell in phononic crystals, the band-diagram can be plotted by restricting the wave vectors to the first Brillouin zone [13]. In practice, the response of a structure to the external source of excitation may be calculated. In a previous paper [2], the applications of the acoustic wave solver on two dimensional square lattices are considered. The simulations are performed by simulating a Gaussian pulse propagating in the normal incident x-direction to the square lattices, or along the  plane in the sense of the first Brillouin zone. The target is to diminish the noise level downstream of the lattices. While proven successful, the performance of the cylinders remains research interests when the real 3D structures are considered. Since the three dimensional effects may alter the actual behavior of the phononic crystals, the deformations of each cylinder in the axial x and transverse z directions are selected to be the design variables, again with Hicks-Henne functions. Periodic boundary conditions are applied to the transverse direction, and solid wall is considered at the bases of the cylinders in the y directions. The inclusions are considered to be stainless-steel surrounded by air, with material properties  0 in  6131 0 out and  K0 in  8 105  K0 out . The unit length of each lattice is 10 cm, the radius of the cylinders is 2.9 cm, and the volume fraction is thus f = 0.066. The attenuation of the sound pressure at the downstream (+x direction) of the phononic crystals in a frequency range is chosen to be the design target, with cost function given by

I

2

1

pˆ  d

(45) s.t. Rk = 0  L    U with 1 = 2.5 kHz and 2 = 3.0 kHz, and pˆ is the frequency domain solution. The original cylindrical shapes are shown in Fig. 11, whereas the optimized distributions are given in Fig. 12. It can be seen that the cylinders tend to move towards both the upstream (-x) and downstream (+x) directions. Note that symmetries are enforced in the transverse z direction where the periodic boundary conditions are applied. The transmissions of the pressure are plotted against the frequency in Fig. 13 for both the original and the optimal shapes. The black solid line represents the pressure with the original shape, and the blue dashed line represents the one with the optimization. It is shown from this figure that the sound pressures at the target frequency range are effectively reduced by the shape optimization. However, it is also noticed that the noise transmission at lower frequencies are increased as a consequence.

Figure 11. The original shapes of phononic crystals 9 American Institute of Aeronautics and Astronautics

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

Figure 12. The optimal shapes of phononic crystals for noise reductions in the frequency range from 2.5 kHz to 3.0 kHz

Figure 13. The transmissions of sound pressure at different frequencies.

V. Conclusion A procedure using adjoint-based sensitivity analysis for shape optimization has been presented. The timedependent sensitivity analysis enables the design for arbitrary cost functions and applicable to large numbers of design variables, and is suitable for designs in complex linear and nonlinear engineering problems such as casing treatments for the performance of compressors [16, 17]. The technique has been applied to acoustic metamaterials and phononic crystals in three dimensions. The optimized shape for phononic crystals has shown the attenuation of sound pressure at frequencies from 2.5 kHz to 3.0 kHz. Consequently, the shape of a type of metamaterial has been optimized to achieve desired effective bulk modulus and impedance.

Acknowledgments This work was supported by the THEC Center of Excellence in Applied Computational Science and Engineering. The support is greatly appreciated. 10 American Institute of Aeronautics and Astronautics

References 1.

2.

3. 4. 5. 6. 7.

Downloaded by Kidambi Sreenivas on July 2, 2016 | http://arc.aiaa.org | DOI: 10.2514/6.2016-3830

8. 9.

10.

11. 12. 13. 14. 15. 16. 17.

Lin, W., Newman III, J. C., and Anderson, W. K. "Broadband Shape and Topology Optimization of Acoustic Metamaterials and Phononic Crystals," 17th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, AIAA Aviation and Aeronautics Forum and Exposition 2016. Washington, D.C., 2016, AIAA Paper 2016-3216. Lin, W., Anderson, W. K., Newman III, J. C., and Zhang, X. "Shape Optimization of Two-Dimensional Acoustic Metamaterials and Phononic Crystals with a Time-Dependent Adjoint Formulation," AIAA SciTech. San Diego, CA, 2016, AIAA Paper 2016-1908. Yamaleev, N. K., Diskin, B., and Nielsen, E. J. "Local-in-time adjoint-based method for design optimization of unsteady flows," Journal of computational physics Vol. 229, No. 14, 2010, pp. 5394-5407. Bonhaus, D. L. "A Higher Order Accurate Finite Element Method for Viscous Compressible Flows." Vol. Ph.D., Virginia Polytechnic Institute and State University, 1998. Barth, T. J. "Numerical methods for gasdynamic systems on unstructured meshes," An introduction to recent developments in theory and numerics for conservation laws. Springer, 1999, pp. 195-285. LeVeque, R. J. Finite volume methods for hyperbolic problems: Cambridge university press, 2002. Saad, Y., and Schultz, M. H. "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems," SIAM Journal on scientific and statistical computing Vol. 7, No. 3, 1986, pp. 856-869. Karypis, G., and Kumar, V. "METIS--unstructured graph partitioning and sparse matrix ordering system, version 2.0," 1995. Sigmund, O., and Jensen, J. S. "Systematic design of phononic band–gap materials and structures by topology optimization," Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences Vol. 361, No. 1806, 2003, pp. 1001-1019. Newman III, J. C., Taylor III, A. C., Barnwell, R. W., Newman, P. A., and Hou, G. J.-W. "Overview of sensitivity analysis and shape optimization for complex aerodynamic configurations," Journal of Aircraft Vol. 36, No. 1, 1999, pp. 87-96. Anderson, W. K., Newman, J. C., Whitfield, D. L., and Nielsen, E. J. "Sensitivity analysis for Navier-Stokes equations on unstructured meshes using complex variables," AIAA journal Vol. 39, No. 1, 2001, pp. 56-63. Nielsen, E. J., and Park, M. A. "Using an adjoint approach to eliminate mesh sensitivities in computational design," AIAA journal Vol. 44, No. 5, 2006, pp. 948-953. Brillouin, L., and Parodi, M. Propagation of waves in periodic structures: Dover, New York, 1953. Fokin, V., Ambati, M., Sun, C., and Zhang, X. "Method for retrieving effective properties of locally resonant acoustic metamaterials," Physical review B Vol. 76, No. 14, 2007, 144302. Hicks, R. M., and Henne, P. A. "Wing design by numerical optimization," Journal of Aircraft Vol. 15, No. 7, 1978, pp. 407-412. Lin, W., Sreenivas, K., Webster, R. S., and Hyams, D. G. "Effect of Casing and Tip Modifications on the Performance of an Axial Flow Stage," AIAA Paper 2012-0475. Lin, W., Sreenivas, K., Webster, R. S., and Hyams, D. G. "Effect of Casing Groove Locations on the Performance of an Axial Flow Stage," AIAA Paper 2013-3632.

11 American Institute of Aeronautics and Astronautics