Abstract

Phononic metamaterials (PMs) exhibit frequency ranges at which elastic waves are attenuated called band gaps. However, this phenomenon is highly sensitive to geometrical variations and the unit cell’s mechanical properties. It is useful to have accurate sensitivity information to identify the variables that produce the highest impact on band gaps and guide the design of PMs with a desired wave propagation behavior. Current methodologies for sensitivity analysis in PMs, such as the finite difference method (FDM), are computationally inefficient, subjected to subtraction cancelation errors, and their accuracy is highly dependent on the magnitude of the perturbation step size. In this study, we introduce a new computational methodology to perform parameter sensitivity in the dynamic behavior of PMs using the multicomplex Taylor series expansion (ZTSE) coupled with Bloch’s theorem. The methodology allows one to obtain arbitrary-order sensitivities with high accuracy. In contrast to FDM, this methodology is computationally more efficient, eliminates the step size selection issue, and is not subjected to subtractive cancelation errors. Also, we show how the method can be applied using real algebra solvers. We limit our analysis to linear undamped PMs. The methodology using ZTSE with Bloch’s theorem is presented in numerical examples for the diatomic lattice and a 2D square lattice, where we compute up to third-order sensitivities. The results show a maximum normalized root-mean-squared deviation in the order of 10−9 for the diatomic lattice and in the order of 10−8 for the 2D square lattice when compared to the analytical solutions.

1 Introduction

Phononic metamaterials (PMs) are a class of manmade periodic materials that exhibit unusual dynamic behavior. The mechanical properties of PMs are determined by the geometry, topology, and materials of its repeating unit cell, yielding to Bragg reflections of acoustic or elastic waves. Therefore, these materials exhibit ranges of frequencies at which waves are either allowed to propagate (pass bands) or blocked in one or more directions (stop bands). These particular dynamic behaviors observed in PMs have drawn the attention of researchers in recent years, allowing advances in applications such as waveguides and filters [1], imaging [2], focusing [3], sound collimation [4], energy harvesting [5], negative refraction index [6], negative elastic modulus [7,8], locally resonant sonic materials [9], negative bulk modulus [10], negative mass density [11], peculiar transmission peaks [12], structures for acoustic collimation [13], slow-wave effect [14], and room temperature THz modulators for ultrafast wireless communication [15].

The wave propagation characteristics of PMs are described by their dispersion relations. These relations can be used to identify the wave group velocity, the directionality of the propagating waves, and the bandgaps presented in the material [16,17]. It is widely reported that the dispersion relations for PMs are highly sensitive to mechanical and geometrical variations [18,19]. Therefore, performing sensitivity analyses for these relations is essential for understanding the wave propagation behavior of periodic materials and design goals. For instance, the calculation of sensitivities is a crucial step for designing PMs through topology optimization strategies. However, optimizing the dispersion diagrams of phononic metamaterials is considered a complicated task as it involves calculating sensitivities of multiple eigenvalues [20]. Issues related to the calculation of eigenvalue sensitivities with respect to system design variables have been reported as one of the core limiting factors in the design of PMs using topology optimization [21]. Inaccurate sensitivities result in slow convergence ratios in optimization algorithms or can even halt the optimization process [22].

Different numerical methods have been developed to compute sensitivities in eigenvalue problems. Most of the seminal works in this area have been dedicated to the development of semi-analytical methods for the calculation of first-order and second-order sensitivities for the eigenvalues of structural systems [2329]. In these works, the eigenvectors are calculated using numerical methods and then used with the sensitivities of the stiffness and mass matrices in algebraic equations to obtain the eigenvalues’ sensitivities. It is important to highlight that these methods do not specify a procedure to obtain the sensitivities of the stiffness and mass matrices, and there are no expressions available for calculating higher-order sensitivities of the eigenvalues. Alternatively, one can use the finite difference method (FDM) for the calculation of sensitivities in eigenvalue problems [30]. In FDM, sensitivities are obtained by solving at least two eigenvalue problems: one at the nominal value of the input parameters and the other after applying a real-valued perturbation to the nominal value of the input parameters(s) of interest. The sensitivity of the eigenvalues with respect to the perturbed parameter(s) is then calculated by subtracting these two solutions and dividing by the magnitude of the perturbation used. The implementation of FDM is straightforward; however, the application of this method may lead to subtraction cancelation errors, and its accuracy is highly dependent on the magnitude of the applied perturbation [3134]. Despite the existence of methods to compute sensitivities of the eigenvalues (i.e., the modal method), these methods require one to obtain the eigenvectors in a first instance and then use them in conjunction with the sensitivities of the mass and stiffness matrices, which traditionally are obtained with FDM [35]. However, this means that current methods to obtain sensitivities on eigenvalue problems suffer from the same disadvantages inherited from FDM.

The sensitivities of eigenvalue problems can be also calculated using the complex Taylor series expansion (CTSE) method [36,37], a first-order numerical differentiation technique similar in spirit and concept to finite differencing but with significant numerical advantages. CTSE uses the orthogonality of the real and imaginary axes of the complex plane to calculate sensitivities without difference operations; hence, the specification of the step size is not an issue as long as it is sufficiently small, e.g., less that 10−10. Unlike finite differencing, CTSE does not require the difference of two analyses. Instead, the sensitivities are computed using a small perturbation along the imaginary axis of the input parameters, solving the eigenvalue problem, and then extracting the imaginary part of the eigenvalues and dividing by the magnitude of the perturbation. In this method, the sensitivity is accurate because no differencing operations are needed. Therefore, the step size can be made arbitrarily small with no concern about round-off error and also that the computation of first-order modal sensitivities is not hampered by repeated eigenvalues [38]. However, CTSE has limited application in the analysis of PMs, as the dispersions relations for these materials are obtained via Bloch’s theorem [39,40], which requires the solution of a complex value eigenvalue problem. Therefore, the imaginary axis required for the perturbation is already occupied with the numerical considerations resulting from the application of Bloch’s theorem.

In this study, we extend the CTSE method to calculate sensitivities in eigenvalue problems with hypercomplex mathematics. The extension of CTSE to the hypercomplex mathematics method is known as the multicomplex Taylor series expansion (ZTSE) and allows one to compute high-order and mixed sensitivities up to any order [41]. In essence, ZTSE consists of multiple independent imaginary axes that are independent and orthogonal, e.g., i1, i2, and i3. Then, the calculation of high-order and mixed partial sensitivities is obtained by perturbing the parameter of interest along the various imaginary axes. Besides computing higher-order sensitivities, one of the advantages of using hypercomplex mathematics is that one of the complex variable axes can be reserved to apply Bloch’s theorem. Therefore, this approach provides a new methodology to enable the sensitivity analysis of the dynamic behavior of linear undamped PMs by integrating multicomplex algebra with Bloch’s theorem. Also, as most open-source eigenvalue solvers cannot handle complex algebra operations, we demonstrate this new methodology using real algebra by taking advantage of the Cauchy-Riemann (CR) matrix representation of a complex number [42]. Finally, it is worth highlighting that the ZTSE has been verified in a number of application areas such as elastic analysis [43], fatigue [38,44], thermoelasticity [45], fracture mechanics [4648], plasticity [49], residual stresses [50], creep [51], heat transfer [52], variance estimates [53], and structural dynamics [32], among others; however, to the authors’ best knowledge, this is the first demonstration of ZTSE for phenomena described by physics involving complex value mathematics.

The outline for the rest of the article is as follows: Sec. 2 provides a summary of the background information, including the analysis of wave propagation in PMs using Bloch’s theorem, the calculation of first-order sensitivities using CTSE, the extension of CTSE to multicomplex mathematics, and the Cauchy-Riemann matrix representation of multicomplex numbers that eliminates the need of using complex variable solvers to obtain sensitivities using ZTSE. Section 3 introduces the methodology to calculate high-order sensitivities in PMs by integrating ZTSE with Bloch’s theorem and includes application examples in 1D and 2D mass-spring lattices. Finally, in Sec. 4, we provide the conclusions and future work.

2 Background

2.1 Wave Propagation in Phononic Metamaterials.

The wave propagation in PMs is traditionally studied by applying Bloch’s theorem [39,40]. According to this theorem, the dynamic behavior of a PM can be described by analyzing the wave motion of its fundamental unit cell. For a wave propagating without attenuation through such material with spatial periodicity, Bloch’s theorem states that the local (i.e., over a unit cell) change in wave amplitude does not depend on the specific location of the cell in the periodic ensemble. Then, the displacement field over a unit cell satisfies [54]:
u(x+a)=u(x)e±iκa
(1)
where κ is the wave vector, a is the lattice translation vector, and u(x) and u(x + a) describe the response at the field points x and x + a. The dispersion relation, ω(κ), of a material under free wave motion can be obtained by discretizing its equation of motion in a finite element framework. Considering a linear undamped system, the discretized equation of motion has the form:
ω2[M]{U}=[K]{U}
(2)
where ω is the frequency of harmonic wave propagation, K and M are the stiffness and mass matrices, respectively, and U is the nodal displacement vector for the unit cell. Applying the Bloch-Floquet Theorem to this relation results in the generalized eigenvalue problem described by
[KRω(κ)2MR]uR=0
(3)
where KR, MR, and uR are the reduced stiffness, mass matrices, and nodal displacement vector, respectively, resulting after removing the redundant equations from Eq. (2). Each assignment of κ yields a particular instance of the ω(κ). The relations ω(κ) provide essential information about the wave propagation behavior in a material, such as wave group velocity, wave directionality, and frequency band gaps. Therefore, obtaining sensitivities for ω(κ) is essential for understanding the behavior of periodic materials and for design goals.

2.2 First-Order Sensitivity Using CTSE.

The CTSE is a first-order differentiation technique that uses the orthogonality of the real and imaginary axis of the complex plane to calculate sensitivities. Unlike FDM, CTSE does not require the difference between the two analyses. Instead, the sensitivities are computed using a small (i.e., h ≤ 10−10) perturbation along the imaginary axis [36]. That is, variable X = x0 is perturbed to X = x0 + hi1, where i denotes an imaginary number and h denotes the perturbation step size. For reference, schematics of the perturbation schemes for FDM and CTSE are shown in Figs. 1(a) and 1(b), respectively. By using the complex Taylor series expansion, a complex function f(x + hi1) can be expanded as follows:
f(x+hi1)=f(x)+hi11!df(x)dx+(hi1)22!d2f(x)dx2++(hi1)nn!dnf(x)dxn+O(hn+1)
(4)
where dnf(x)dxn is the nth sensitivity of f(x) with respect to x. Then, taking the imaginary part of both sides of Eq. (4) and ignoring the higher-order terms yields an estimate of the first sensitivity:
df(x)dxIm1[f(x+hi1)]h+O(h2)
(5)

As there is no subtraction involved, the round-off error is eliminated. On the other hand, the truncation error is reduced by reducing h, e.g., using h ≤ 10−10 times, the smallest problem parameter will produce a truncation error below the machine precision.

Fig. 1
Perturbation diagrams, the dashed circles represent the unperturbed variable and the continuous circles represent the perturbed variable. (a) Real-step perturbation corresponding to FDM. (b) Perturbation corresponding to CTSE. (c) Multicomplex perturbation corresponding to ZTSE to compute second-order sensitivities. (d) Multicomplex perturbation corresponding to ZTSE for computing second-order cross sensitivities.
Fig. 1
Perturbation diagrams, the dashed circles represent the unperturbed variable and the continuous circles represent the perturbed variable. (a) Real-step perturbation corresponding to FDM. (b) Perturbation corresponding to CTSE. (c) Multicomplex perturbation corresponding to ZTSE to compute second-order sensitivities. (d) Multicomplex perturbation corresponding to ZTSE for computing second-order cross sensitivities.
Close modal

2.3 Extension of CTSE to Multicomplex Mathematics: High-Order Sensitivities.

Extending CTSE to multicomplex mathematics allows computing higher-order and mixed sensitivities up to any order [41]. Multicomplex numbers expand the concept of complex numbers by adding multiple imaginary axes to the real quantity (e.g., i1, i2, i3, …, in). These imaginary axes are independent and orthogonal [55]. For more details about the properties of multicomplex mathematics, the reader is directed to Ref. [42]. The calculation of high-order and mixed partial sensitivities using the ZTSE is achieved by perturbing the parameter of interest along the imaginary axes, and the number of imaginary directions perturbed determines the maximum order of the sensitivity to be computed. For instance, a second-order sensitivity of f(x) will require a bicomplex number (e.g., two independent imaginary axes) with a small perturbation along the two imaginary axes as shown in Fig. 1(c). By performing the Taylor series expansion of the bicomplex function f(x + h(i1 + i2)), we have
f(x+h(i1+i2))=f(x)+h(i1+i2)df(x)dx+h2(i12+2i1i2+i22)2!d2f(x)dx2+h3(i13+3i12i2+3i1i22+i23)3!d3f(x)dx3+O(h4)
(6)
The error of the approximation O() depends on the number of terms preserved from the complex Taylor series expansion before the truncation. Considering that multicomplex multiplication is commutative [56], we have that
inim={1n=m,inimnm,n,mN
(7)
Therefore, Eq. (7) can be simplified to
f(x+h(i1+i2))=f(x)+h(i1+i2)df(x)dx+h22(i1i21)2!d2f(x)dx2+O(h3)=(f(x)h2d2f(x)dx2)+i1(hdf(x)dx)+i2(hdf(x)dx)+i1i2(h2d2f(x)dx2)+O(h3)
(8)
Second-order sensitivities can be obtained by taking the imaginary part corresponding to i1i2 and ignoring the higher-order terms:
d2f(x)dx2Im12[f(x+h(i1+i2))]h2+O(h2)
(9)
Similarly, first-order sensitivities are also obtained by taking the imaginary part corresponding to i1 and ignoring the higher-order terms or by taking the imaginary part corresponding to i2 and ignoring the higher-order terms:
df(x)dxIm1[f(x+h(i1+i2))]h+O(h)df(x)dxIm2[f(x+h(i1+i2))]h+O(h)
(10)
In the case of cross sensitivities, each independent variable is perturbed along different imaginary axes, as shown in Fig. 1(d). Then, the sensitivities of the function f(x, y) are given by (assuming x is perturbed along i1 and y along i2),
f(x,y)xIm1[f(x+h1i1,y+h2i2)]h1+O(h1)f(x,y)yIm2[f(x+h1i1,y+h2i2)]h2+O(h2)2f(x,y)xyIm12[f(x+h1i1,y+h2i2)]h1h2+O(h1h2)
(11)
Extending this procedure to obtain nth-order sensitivities of a function f(x):
nf(x)xnIm1n[f(x+h(jnij))]hn+O(hn)
(12)

It is important to note that the numerical considerations of Bloch’s theorem in PMs are already expressed in complex variables (see Eq. (1)). Therefore, it is necessary to reserve one imaginary axis for these periodic boundary conditions and additional imaginary axes for the sensitivities. Then, to compute first-order sensitivities, it is necessary to represent all the input variables at least as bicomplex numbers (e.g., two imaginary axes). For second-order sensitivities, all the variables should be expressed at least as tricomplex numbers (e.g., three imaginary axes). Generalizing, if n imaginary axes are available, it is possible to calculate the sensitivities from order 1 to n−1, e.g., tricomplex numbers allow calculating first- and second-order sensitivities in PMs.

2.4 Cauchy-Riemann Matrix Notation of Multicomplex Numbers.

As most numerical packages do not support multicomplex variable operations, it is possible to take advantage of the multicomplex algebra isomorphism to represent any multicomplex number as a matrix with all real numbers. This representation of multicomplex numbers is known as CR matrix notation and allows the use of matrix functions and arithmetic operations as homologous of the operations between multicomplex numbers [42]. For reference, for a multicomplex number with n imaginary axes zη=z1η1+z2η1in, where zηCη and z1η1,z2η1Cη1, the general CR representation is given by Eq. (13), where [z1η1] and [z2η1] are recursively CR representations of multicomplex numbers with η − 1 imaginary axes:
[zη]=[[z1η1][z2η1][z2η1][z1η1]]=[[[z11η2][z12η2][z12η2][z11η2]][[z21η2][z22η2][z22η2][z21η2]][[z21η2][z22η2][z22η2][z21η2]][[z11η2][z12η2][z12η2][z11η2]]]
(13)
For example, in the case when η = 1, the CR representation of the complex number zη=1 = x1 + x2i1, where zη=1 ∈ ℂ1 and x1, x2 ∈ ℝ are shown in Eq. (14).
[zη=1]=[x1x2x2x1]
(14)
In the case of a bicomplex number with η = 2 imaginary axes, zη=2=z1η=1+z2η=1i2, where zη=2 ∈ ℂ2, and each zjη=1 is a complex number (z1η=1,z2η=1C1) of the form z1η=1=x1+x2i1 and z2η=1=x3+x4i1, where x1, x2, x3, x4 ∈ ℝ. The CR representation is given by Eq. (15):
[zη=2]=[[z1η=1][z2η=1][z2η=1][z1η=1]]=[x1x2x3x4x2x1x4x3x3x4x1x2x4x3x2x1]
(15)
For a tricomplex number with η = 3 imaginary axes, zη=3 = x1 + x2i1 + x3i2 + x4i1i2 + x5i3 + x6i1i3 + x7i2i3 + x8i1i2i3, where zη=3 ∈ ℂ3 and x1, x2, x3, x4, x5, x6, x7, x8 ∈ ℝ. The CR is given as follows:
[zη=3]=[x1x2x3x4x5x6x7x8x2x1x4x3x6x5x8x7x3x4x1x2x7x8x5x6x4x3x2x1x8x7x6x5x5x6x7x8x1x2x3x4x6x5x8x7x2x1x4x3x7x8x5x6x3x4x1x2x8x7x6x5x4x3x2x1]
(16)

3 Sensitivity Analysis in Phononic Materials Using ZTSE

We propose a new methodology that combines ZTSE with Bloch’s theorem described in Eqs. (1) and (3) to obtain sensitivity information on the dispersion diagrams that describe the dynamic behavior of PMs. The methodology is summarized in Fig. 2(a), and it requires converting all the parameters involved in the problem (e.g., mass, stiffness, geometry, and nodal coordinates) to multicomplex variables represented as real-value matrices using the CR notation. This conversion facilitates the operations between multicomplex quantities by allowing the use of matrix operations, eliminating the necessity of supporting multicomplex algebra by the numerical package used. When converting all the parameters to the CR notation, it is important to recall that one imaginary axis should be reserved to apply periodic boundary conditions using Bloch’s theorem. Therefore, n+1 imaginary axes should be considered to obtain sensitivities of order n. Then, to get sensitivities, imaginary perturbations are added to the parameters of interest. These perturbations must be incorporated using their CR notation. To minimize truncation errors, the perturbation magnitudes must be a low number (e.g., h ≤ 10−10) times the nominal value of the variable of interest.

Fig. 2
(a) Flowchart for the sensitivity analysis of PMs using ZTSE and Bloch’s theorem. (b) Schematic representation of the unit cell for a diatomic system. Each variable (mass value, spring stiffness, and nodal displacement) is represented as bicomplex variables using different superscripts for each axis (Re: Real, is: Imag 1, iB: Imag 2, and i12: Imag 12). Each degree-of-freedom is represented with different geometric shapes (circles: DOF1 and squares: DOF2). (c) Generalized eigenvalue problem corresponding to the discrete system shown in (b). Cauchy-Riemann matrix notation is used to express the bicomplex variables in the system of equations.
Fig. 2
(a) Flowchart for the sensitivity analysis of PMs using ZTSE and Bloch’s theorem. (b) Schematic representation of the unit cell for a diatomic system. Each variable (mass value, spring stiffness, and nodal displacement) is represented as bicomplex variables using different superscripts for each axis (Re: Real, is: Imag 1, iB: Imag 2, and i12: Imag 12). Each degree-of-freedom is represented with different geometric shapes (circles: DOF1 and squares: DOF2). (c) Generalized eigenvalue problem corresponding to the discrete system shown in (b). Cauchy-Riemann matrix notation is used to express the bicomplex variables in the system of equations.
Close modal
To exemplify this procedure, consider the discrete system with two masses shown in Fig. 2(b), representing the unit cell of an infinitely periodic monoatomic lattice. Each point mass is represented with a different shape: circles for m1, squares for m2, and the stiffness of the springs, c, is represented with parallelograms. Considering first-order sensitivities for any input parameter, it is necessary to use bicomplex numbers. Each axis representing the bicomplex numbers is identified with different colors: orange for the real axis, blue for i1, red for i2, and green for the cross-terms i12. In this example, we are reserving i1 for sensitivities with respect to any input parameter denoted from now on with is, and i2 is reserved for applying Bloch’s theorem and denoted with iB. In Fig. 2(c), we present the assembled generalized eigenvalue problem corresponding to Eq. (2). As shown, each position of the mass matrix is replaced with a CR representation for each bicomplex mass. The displacement vectors contain all the real and imaginary degrees–of-freedom, and the mass and stiffness matrices are constructed using the finite element formulation of spring elements as follows:
[M]=[m100m2][K]=[cccc]
(17)

After applying boundary conditions using Bloch’s theorem, the dispersion diagram for the diatomic system can be obtained by solving the generalized eigenvalue problem described in Fig. 3, where a is the size of the unit cell and κ is the wavenumber.

Fig. 3
Reduced eigenvalue problem
Fig. 3
Reduced eigenvalue problem
Close modal

The resulting system after applying Bloch’s theorem is nonself-adjoint. Therefore, the QR algorithm is used to solve the eigenvalue problem [5759]. The QR algorithm is recognized as one of the most important algorithms for eigenvalue computations due to its robustness [6062]. The most traditional numerical algorithms to solve eigenvalue problems achieve the stability for self-adjoint, Hermitian, or tridiagonal Hermitian eigenvalue problems [63]. This includes the Davidson method [64], the implicit restarted Lanczos [65], and the locally optimal block preconditioned conjugate gradient [66]. In this study, we aim to solve a nonself-adjoint eigenvalue problem. Nonself-adjoint eigenvalue problems are more difficult and computationally demanding to solve than self-adjoint problems [67]. Generally, self-adjoint eigenvalue problems are well conditioned, while nonself-adjoint problems are not [68]. Schur decomposition-based methods to solve eigenvalue problems, such as the QR algorithm, have been standard methods for solving small- to medium-sized eigenvalue problems [69,70]. The reader is referred to Ref. [63] for a thorough comparison of the available Open-source libraries to solve nonself-adjoint eigenvalue problems of large-scale power systems.

The QR algorithm requires converting the generalized eigenvalue problem into a standard eigenvalue problem. To this end, both sides of Eq. (3) are premultiplied by the negative of the inverse of the reduced mass matrix to obtain
[Mr1Krω(κ)2I]ur=0
(18)
The eigenvalues of Eq. (18) are obtained by assigning a value of κ and calculating the product [S]0 = − [Mr]−1[Kr]. Then, the QR decomposition is performed over [S]0, and the factors are multiplied and iterated as follows:
[S]j=[Q]j[R]j[S]j+1=[R]j[Q]j
(19)
If the size of the matrix [S]0 is N2, the QR algorithm requires N3 iterations to find the eigenvalues. In a traditional analysis, the QR algorithm would produce an upper triangular matrix with the eigenvalues located along the diagonal of the resulting matrix [S]N3. However, in the present analysis, the multicomplex eigenvalues along the diagonal are replaced with CR matrices. The number of CR matrices inside [S]N3 is equal to the number of degrees-of-freedom in the reduced system (Eq. (19)). To illustrate this point, Eq. (20) presents a system with N degrees-of-freedom, and each CR matrices along the diagonal of [S]N3 represents bicomplex numbers. In this system, the eigenvalues and the corresponding imaginary terms with the sensitivities information are in the first column of each CR matrix. For this example, Im1 is reserved for the first-order sensitivity with respect to any input parameter and denoted as ImS, and Im2 is reserved for the imaginary axis required by the periodic boundary conditions from Bloch’s theorem and denoted as ImB.
[S]N3=[λ1Reλ1ImSλ1ImBλ1Im12λ1ImSλ1Reλ1Im12λ1ImBλ1ImBλ1Im12λ1Reλ1ImSλ1Im12λ1ImBλ1ImSλ1ReλNReλNImSλNImBλNIm12λNImSλNReλNIm12λNImBλNImBλNIm12λNReλNImSλNIm12λNImBλNImSλNRe]
(20)
Given the matrix [S]N3, one can obtain the dispersion diagrams, ω(κ), and the corresponding sensitivities by calculating the square root of the resulting CR matrices of the eigenvalues using the Denman and Beavers iteration method [71]. Other methods used for calculating the square root of matrices [72] are known to fail when CR is used to represent multicomplex numbers [56].
[ω1Reω1ImSω1ImBω1Im12ω1ImSω1Reω1Im12ω1ImBω1ImBω1Im12ω1Reω1ImSω1Im12ω1ImBω1ImSω1Re]=sqrt([λ1Reλ1ImSλ1ImBλ1Im12λ1ImSλ1Reλ1Im12λ1ImBλ1ImBλ1Im12λ1Reλ1ImSλ1Im12λ1ImBλ1ImSλ1Re])
(21)
By applying Eqs. (10) and (11) to the color-highlighted terms shown in Eqs. (20) and (21), the eigenvalues and their corresponding first-order sensitivities are obtained as follows:
{ω}={ω1ReωNRe}{ωx1}={ω1ImShωNImSh}
(22)

For reference, the implementation to compute second-order sensitivities of this example is shown in the Supplementary Material I on the ASME Digital Collection.

The CR matrix notation of multicomplex numbers allows to expand the complex standard eigenvalue problem presented in Ref. [73] to include general damping in the formulation of the sensitivity analysis. However, the QR algorithm produced inaccurate results when compared with analytical expressions of PMs with damping. Hence, further research is necessary to better understand the impact of the eigenvalue problem solver used in the analysis of PMs with damping. For lucidity, we relegate this procedure for future work and present the implementation of ZTSE for linear undamped discrete systems.

3.1 Extension to Continuum Finite Elements.

The use of ZTSE can be extended to continuum finite elements, and this has been called the multicomplex finite element method (ZFEM) [74]. They key aspect of ZFEM is that the nodal coordinates are multicomplex variables, having real and imaginary components. The real components define the physical location of the nodes (similarly as in traditional finite elements method), and the imaginary components define the “perturbations” of the real (physical) nodal coordinates, with the specific perturbation defining the derivative to be computed. To illustrate this concept, Fig. 4 shows an example of a 16-noded quadrilateral element with four real nodes, four nodes corresponding to the first imaginary axis, four nodes corresponding to the second imaginary axis, and four more nodes corresponding to the mixed imaginary axis 1,2. The other input parameters, such as material properties (i.e., density, Young’s modulus, Poisson’s ratio) and boundary conditions are also multicomplex variables. In the process of computing the mass and stiffness matrices, each multicomplex variable is represented with the CRM notation of multicomplex numbers.

Fig. 4
Sixteen-noded isoparametric element (four real, four imaginary 1, four imaginary 2, and four mixed imaginary 1,2)
Fig. 4
Sixteen-noded isoparametric element (four real, four imaginary 1, four imaginary 2, and four mixed imaginary 1,2)
Close modal

ZFEM can be implemented in the sensitivity analysis of structural systems with large number of degrees-of-freedom, while preserving high accuracy; however, it is important to reiterate that the convergence rate of the QR algorithm is N3, where N is the total number of degrees-of-freedom, and as ZTSE multiplies the number of real degrees-of-freedom by 2η, where the total number of imaginary axis η is given by the highest order of sensitivity to compute plus the imaginary axis required by Bloch’s theorem. This means that the computational runtime grows as function of the number of degrees-of-freedom used and the highest order of sensitivity to compute. An additional consideration is that the accuracy of the derivatives is dependent on the accuracy of the real part; hence, increasing the number of nodes will increase the accuracy of the real part and its sensitivities. In this study, we limit our analysis to discrete systems due to their simplicity allowing an easy understanding and demonstration of the proposed methodology. The implementation and verification using other continuum elements will be relegated to future works.

4 Case Studies

To illustrate this new methodology for the sensitivity analysis of PMs, we apply the method to the sensitivity analysis of the diatomic and a 2D squared lattice in the following sections. The numerical results from the ZTSE and Bloch’s theorem methodology are compared with analytical equations, and the errors are measured using the normalized root-mean-square deviation (NRMSD) as described in Eq. (23). Also, a comparison of the runtime between the proposed methodology and FDM is presented.

NRMSD=1max(yAnalytical)min(yAnalytical)jN(yjAnalyticalyjNumerical)2N
(23)

4.1 Diatomic Lattice.

The diatomic lattice is an infinite array of alternating masses m1 and m2 interconnected with springs with linear stiffness c. Also, the two consecutive masses of the same kind are separated a distance a from each other, and each mass m1 and m2 possesses one translational degree-of-freedom us and vs, respectively. Figure 5 shows a schematic of the diatomic lattice. The analysis of wave motion for this system can be easily derived by obtaining the equations of motion for two neighboring masses, assuming the absence of external forcing, and assuming harmonic motion. Then, after applying the periodic conditions from Bloch’s theorem, the dispersion relation for this system is given by
ω={ωO=cm1m2(m1+m2+m12+2m1m2cos(κa)+m22)ωA=cm1m2(m1+m2m12+2m1m2cos(κa)+m22)
(24)
Fig. 5
Schematic for diatomic lattice
Fig. 5
Schematic for diatomic lattice
Close modal

Since the diatomic lattice has two active degrees-of-freedom, the dispersion relation has two branches called optical and acoustic modes denoted in Eq. (24) as ωO and ωA, respectively.

Next, we present the detailed numerical implementation of our method as illustrated in the flowchart of Fig. 2(a). First, we perturb the variable of interest, and second, we use CR matrix representation for each input variable. Out of briefness, we only present the procedure for first-order sensitivities of m2, but the procedure can be extrapolated for the other variables and higher-order sensitivities as well. In this example, we use the first imaginary axis for the sensitivities and the second imaginary axis for imposing Bloch’s periodic boundary conditions denoted with iB. With this in mind, the CR matrix representation for the input parameters is as follows:

[m1]=m1I4=[m10000m10000m10000m1];[c]=cI4=[c0000c0000c0000c];[m2]=m2+hi1=[m2h00hm20000m2h00hm2];[a]=aI4=[a0000a0000a0000a];[iB]=0+1i2=[0010000110000100];[κ]=κI4=[κ0000κ0000κ0000κ]
(25)
The next step is to assemble the global mass and stiffness matrices of the unit cell. For this example, we use the spring finite element formulation. The resulting generalized eigenvalue problem is expressed as follows:
ω2[[m1][12][0][0][0][m2][0][0][0][m1][12]]{{us}{vs}{us+1}}=[[c][c][0][c][2][c][c][0][c][c]]{{us}{vs}{us+1}}
(26)
where the vectors {us}, {vs}, and {us+1} contain the displacement degrees-of-freedom of the corresponding masses presented in Fig. 5. For the diatomic lattice, the wave propagation is in one dimension only, and therefore, a single wavenumber is enough to describe the irreducible Brillouin zone. For this reason, to apply Bloch-Floquet periodic boundary conditions, we have:
{{us}{vs}{us+1}}=[[1][0][0][1]exp([iB][κ][a])[0]]{{us}{vs}}
(27)
The matrix exponential function algorithm used is further detailed in Refs. [75,76] and yields correct results for CR matrix representation of multicomplex numbers. Finally, we solve the generalized eigenvalue problem for the reduced system:
2ω2[[m1][0][0][m2]]{{us}{vs}}=[[2][c]c[exp([iB][κ][a])+[1]]c[exp([iB][κ][a])+[1]][2][c]]{{us}{vs}}
(28)
Equation (28) is solved following the procedure in Eq. (19). The information of the modes and its sensitivities are found in the system matrix [S] as previously stated in Eqs. (20)(22).
[S]=[ωO2ωO2m2h00ωO2m2hωO20000ωO2ωO2m2h00ωO2m2hωO2ωA2ωA2m2h00ωA2m2hωA20000ωA2ωA2m2h00ωA2m2hωA2][S]=[[ωO2][0][0][ωA2]][ωO]=sqrt([ωO2])=[ωOωOm2h00ωOm2hωO0000ωOωOm2h00ωOm2hωO][ωA]=sqrt([ωA2])=[ωAωAm2h00ωAm2hωA0000ωAωOm2h00ωAm2hωO]ω={Real([ωo])=ωOReal([ωA])=ωAωm2={Imag1([ωO])h=ωOm2Imag1([ωA])h=ωAm2
(29)

4.1.1 Sensitivity Analysis of Diatomic Lattice.

We implemented the multicomplex-step perturbation method to compute the dispersion relation in Eq. (24) and its sensitivities using the following input values, m1 = 1, m2 = 2, c = 1, and a = 1. With perturbation magnitude, h = 10−10 times the nominal value of the perturbed variable. The error quantification was performed with Eq. (23). Figure 6 shows the comparison of the dispersion relation and first- and second-order sensitivities obtained using the symbolic algebra package of matlab and the numerical results obtained with ZTSE. Excellent agreement was found with both methods, resulting in NRMSD values of 5 × 10−14 for the dispersion relation, 6 × 10−12, and 10−10 for the first- and second-order sensitivities with respect to m2, respectively, and 9 × 10−11 for the second-order mixed sensitivity with respect to m2 and a. The remaining comparisons for the first-, second-, and third-order sensitivities can be found in the Supplementary Material II on the ASME Digital Collection. The largest NRMSD value obtain for third-order derivatives was 9 × 10−9.

Fig. 6
(a) Dispersion relation for diatomic lattice. (b) First-order sensitivity of the mass m2. (c) Second-order sensitivity of the mass m2. (d) Second-order mixed sensitivity of the mass m2 and the unit cell length a. Continuous lines represent results of the dispersion relation and its sensitivities for the acoustic and optical modes obtained with matlab’s symbolic algebra, circles and squares represent the numerical results obtained with ZTSE.
Fig. 6
(a) Dispersion relation for diatomic lattice. (b) First-order sensitivity of the mass m2. (c) Second-order sensitivity of the mass m2. (d) Second-order mixed sensitivity of the mass m2 and the unit cell length a. Continuous lines represent results of the dispersion relation and its sensitivities for the acoustic and optical modes obtained with matlab’s symbolic algebra, circles and squares represent the numerical results obtained with ZTSE.
Close modal

4.1.2 Runtime Comparison for Diatomic Lattice.

A comparison was conducted between ZTSE and FDM comparing the runtime for computing the first-, second, and third-order sensitivity of the dispersion relation with respect to the mass m2. Hence, the mean runtime of 1000 runs were used, and the results were normalized with respect to the runtime required by ZTSE to compute the first-order sensitivity. The results are listed in Table 1. We see a lower runtime necessary for ZTSE than FDM to compute the first- and second-order sensitivities. However, for third-order sensitivities, FDM took less computational runtime than ZTSE. These runtimes are equivalent for sensitivities on other input parameters.

Table 1

Mean runtime comparison for ZTSE and FDM over 1000 runs normalized with respect to the mean runtime of the first-order sensitivity using ZTSE

Methodω2m22ω2m223ω2m23
ZTSE×1.00×1.61×4.05
FDM×1.62×2.40×3.22
Methodω2m22ω2m223ω2m23
ZTSE×1.00×1.61×4.05
FDM×1.62×2.40×3.22

4.1.3 Perturbation Step Convergence for Diatomic Lattice.

Figure 7 shows the behavior of the error of the sensitivities as a function of the perturbation step size for FDM and ZTSE. The truncation error decreases as the perturbation step decreases for both methods until a critical point, from which the subtractive error starts dominating the inaccuracy of FDM, increasing the error of FDM as h decreases. In contrast, for ZTSE, the truncation error falls below machine precision, and the sensitivity reaches its maximum precision, from which further decreases in h do not result in accuracy increases for ZTSE.

Fig. 7
NRMSD as a function of the perturbation step magnitude comparison for FDM and ZTSE: (a) first-order, (b) second-order, and (c) third-order sensitivity of the mass m2
Fig. 7
NRMSD as a function of the perturbation step magnitude comparison for FDM and ZTSE: (a) first-order, (b) second-order, and (c) third-order sensitivity of the mass m2
Close modal

4.2 2D Square Lattice.

A rectangular cell was analyzed that consisted of four masses m connected with horizontal and vertical springs with stiffness cx and cy, respectively. Springs with effective stiffness cd connect the diagonals of the lattice. The rectangular cell has horizontal length a and vertical length b. The rectangular cell shown in Fig. 8 represents the unit cell of the infinitely periodic cell.

Fig. 8
2D square lattice with diagonal springs
Fig. 8
2D square lattice with diagonal springs
Close modal
The following expression gives the analytical dispersion relation for the 2D square lattice with diagonal springs:
ωx2=8mcxsin2(κxa2)+cdm(sin2(κxa+κyb2)+sin2(κxaκyb2))ωy2=8mcysin2(κyb2)+cdm(sin2(κxa+κyb2)+sin2(κxaκyb2))ω=ωx2+ωy2
(30)
Next, we present the detailed numerical implementation of our method, as illustrated in the flowchart of Fig. 2(a). First, we perturb the variable of interest, and second, we use CR matrix representation for each input variable. Out of briefness, we only present the procedure for first-order sensitivities of m, but the procedure can be extrapolated for the other variables and higher-order sensitivities. In this example, we use the first imaginary axis for the sensitivities and the second imaginary axis for imposing Bloch’s periodic boundary conditions. With this in mind, the CR matrix representation for the input parameters is expressed as follows:
[m]=m+hi1=[mh00hm0000mh00hm];[cd]=cdI4=[cd0000cd0000cd0000cd];[cx]=cxI4=[cx0000cx0000cx0000cx];[cy]=cyI4=[cy0000cy0000cy0000cy];[a]=aI4=[a0000a0000a0000a];[b]=bI4=[b0000b0000b0000b][κx]=κxI4=[κx0000κx0000κx0000κx];[κy]=κyI4=[κy0000κy0000κy0000κy];[iB]=i2=[0010000110000100]
(31)
The next step is to assemble the global mass and stiffness matrices of the unit cell. For this example, we use the spring finite element formulation. The resulting generalized eigenvalue problem is expressed as follows:
ω2[M]{{u1}{v1}{u2}{v2}{u3}{v3}{u4}{v4}}=[K]{{u1}{v1}{u2}{v2}{u3}{v3}{u4}{v4}}
(32)
where [M] and [K] are the unit cell’s global mass and stiffness matrices, respectively, filled with the bicomplex terms expressed in Eq. (31).
[M]=[[m][14][0][0][0][0][0][0][0][0][m][14][0][0][0][0][0][0][0][0][m][14][0][0][0][0][0][0][0][0][m][14][0][0][0][0][0][0][0][0][m][14][0][0][0][0][0][0][0][0][m][14][0][0][0][0][0][0][0][0][m][14][0][0][0][0][0][0][0][0][m][14]]
(33)
[K]=[[cx][cd][0][cx][0][cd][0][0][0][0][cy][cd][0][0][0][cd][0][cy][cx][0][cx][cd][0][0][0][cd][0][0][0][0][cy][cd][0][cy][0][cd][cd][0][0][0][cx][cd][0][cx][0][0][cd][0][cy][0][cy][cd][0][0][0][0][cd][0][cx][0][cx][cd][0][0][cy][0][cd][0][0][0][cy][cd]]
(34)

For the 2D square lattice, it is necessary to use two wavenumbers to describe the first Brillouin zone. For this reason, to apply Bloch-Floquet periodic boundary conditions, we have:

{{u1}{v1}{u2}{v2}{u3}{v3}{u4}{v4}}=[[1][0][0][1]exp([iB][κx][a])[0][0]exp([iB][κx][a])exp([iB]([κx][a]+[κy][b]))[0][0]exp([iB]([κx][a]+[κy][b]))exp([iB][κy][b])[0][0]exp([iB][κy][b])]{{u1}{v1}}
(35)
Finally, we solve the generalized eigenvalue problem for the reduced system:
ω2[Mr]{{u1}{v1}}=[Kr]{{u1}{v1}}
(36)

Equation (36) is solved following the procedure in Eq. (19). For this problem, each degree-of-freedom carries the information of the wave propagating along the x and y directions. Since our interest is on the magnitude of the frequency, it is necessary to add both degrees-of-freedom before using the matrix square root. The information of the modes ωx2 and ωy2 and their sensitivities are found in the system matrix [S] as previously stated in Eqs. (20)(22).

[S]=[ωx2ωx2mh00ωx2mhωx20000ωx2ωx2mh00ωx2mhωx2ωy2ωy2mh00ωy2mhωy20000ωy2ωy2mh00ωy2mhωy2][S]=[[ωx2][0][0][ωy2]][ω]=sqrt([ωx2]+[ωy2])=[ωωm2h00ωm2hω0000ωωm2h00ωm2hω]ω=Real([ω])=ωωm=Imag1([ω])h=ωm2hh
(37)

4.2.1 Sensitivity Analysis of the 2D Square Lattice.

We implemented the multicomplex-step perturbation method to compute the dispersion relation in Eq. (30) and its sensitivities using the following input values, m = 1, cx = 1, cy = 1, cd = 1, a = 1, and b = 1, with perturbation magnitude h = 1 × 10−10. The error quantification was performed with Eq. (23). Figure 9 shows excellent agreement between the numerical solution obtained with ZTSE and matlab’s symbolic algebra for the dispersion relation and first- and second-order sensitivities. It is important to note here that when κx = κy, for the input parameter values selected, ωx = ωy, which means the system produces repeated eigenvalues that were successfully obtained with our method. Excellent agreement was found with both methods, resulting in NRMSD values of 10−16 for the dispersion relation, 10−15 for both, the first and second-order sensitivities with respect m, respectively, and 5 × 10−15 for the second-order mixed sensitivity with respect to m and a. The remaining first-, second-, and third-order sensitivities can be found in the Supplementary Material III on the ASME Digital Collection. The largest NRMSD value obtain for third-order derivatives was 3 × 10−8.

Fig. 9
(a) Dispersion relation for the 2D square lattice. (b) First-order sensitivity of the mass. (c) Second-order sensitivity of the mass. (d) Second-order mixed sensitivity of the mass and the horizontal separation length.
Fig. 9
(a) Dispersion relation for the 2D square lattice. (b) First-order sensitivity of the mass. (c) Second-order sensitivity of the mass. (d) Second-order mixed sensitivity of the mass and the horizontal separation length.
Close modal
4.2.2 Runtime Comparison for the 2D Square Lattice.

A runtime comparison between ZTSE and FDM for the first-, second-, and third-order sensitivities of the dispersion relation with respect to the mass m was conducted. For this purpose, the mean runtime of 1000 runs were used and normalized with respect to the runtime required by ZTSE to compute the first-order sensitivities. The results are listed in Table 2. Overall, a lower runtime is required by ZTSE than FDM to compute the sensitivities of the dispersion relation for first-, second-, and third-order sensitivities. In contrast to the diatomic lattice, in the 2D square lattice, ZTSE took less computational runtime to compute the third-order sensitivity than FDM.

Table 2

Mean runtime comparison for ZTSE and FDM over 1000 runs normalized with respect to the mean runtime of the first-order sensitivity using CTSE

Methodω2m2ω2m23ω2m3
ZTSE×1.00×1.43×2.28
FDM×2.16×2.99×3.63
Methodω2m2ω2m23ω2m3
ZTSE×1.00×1.43×2.28
FDM×2.16×2.99×3.63
4.2.3 Perturbation Step Convergence for the 2D Square Lattice.

Figure 10 shows the behavior of the sensitivities error as a function of the perturbation step size for FDM and ZTSE. The truncation error decreases as the perturbation step decreases for both methods until a critical point, from which the round-off error starts dominating the inaccuracy of FDM, increasing the error of FDM as h decreases. In contrast, for ZTSE, the truncation error falls below the machine precision number, and the sensitivity reaches its maximum precision, from which further decrease in h does not result in an accuracy increase for ZTSE.

Fig. 10
NRMSD as function of the perturbation step magnitude comparison for FDM and CTSE: (a) first-order, (b) second-order, and (c) third-order sensitivity of the mass
Fig. 10
NRMSD as function of the perturbation step magnitude comparison for FDM and CTSE: (a) first-order, (b) second-order, and (c) third-order sensitivity of the mass
Close modal

5 Conclusions

In this study, we presented a new method that integrates the multicomplex Taylor series expansion with Bloch’s theorem to enable the arbitrary-order sensitivity analysis of phononic metamaterials. The capability of obtaining high-accuracy higher-order sensitivities has potential applications in improving current optimization strategies in the design of PMs and enabling the design of parameter insensitive PMs. We verified this new methodology by analyzing the sensitivities of two discrete PMs: (i) the diatomic lattice and (ii) a 2D squared lattice with diagonal springs. Both systems’ dispersion relations and sensitivities were calculated using our new methodology and compared with the finite differences method (FDM) and analytical functions.

We found that numerical sensitivities obtained with ZTSE and Bloch are significantly more accurate than those obtained with FDM. ZTSE proved to be perturbation step size independent. In contrast, FDM proved to be prone to truncation and round-off errors, leading to grossly inaccurate sensitivities. Moreover, first-, second-, and third-order numerical sensitivities obtained with ZTSE took less computational runtime than those obtained with FDM.

Additional advantages of the proposed method are that it does not require the computation of eigenvectors to obtain the sensitivities of the eigenvalues, nor to know the sensitivities of the stiffness and mass matrices a priori. Also, ZTSE allows for simultaneous computation of eigenvalues and their respective sensitivities in a single procedure, and it is general enough to compute arbitrary-order sensitivities. Furthermore, we showed that the coupled ZTSE-Bloch methodology can take advantage of the Cauchy-Riemann matrix notation for multicomplex numbers, eliminating the need to use eigenvalue problem solvers that handle complex variables. Discrete systems with local resonators can also be studied with ZTSE by increasing the number of degrees-of-freedom.

We posit that the method presented in this study will help close the breach between theoretical predictions and practical developments of PMs by providing a handy yet highly accurate method to perform sensitivity analysis of the dynamical behavior of PMs. Moreover, this methodology will help to understand better the impact of defects and other unintended variations in PMs and will aid the design of new PMs with unique dynamic properties.

Future research efforts will focus on extending this method into available real-value finite element software, including sensitivities to perform uncertainty quantification in PMs, and integrating this strategy to guide the design of new PMs with novel properties. The inclusion of viscous damping promises the ability to compute accurate sensitivities in more realistic analysis. However, further research is necessary to better understand the impact of the eigenvalue problem solver used for the analysis of PMs with damping. Similarly, the presented methodology is limited to the analysis of linear PMs. Future research efforts will commit on the implementation of ZTSE for the sensitivity analysis of nonlinear PMs.

Acknowledgment

The authors would like to acknowledge the support from ONR grant N00014-21-1-2428. D. R. and J. D. N. would also like to recognize the support from the GREAT seed grant awarded by the UTSA Office of Vice President for Research, Economic Development, and Knowledge Enterprise.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. Data provided by a third party listed in Acknowledgment.

Nomenclature

     
  • a =

    unit cell length

  •  
  • c =

    Springs’ stiffness

  •  
  • h =

    perturbation step magnitude

  •  
  • n =

    order of sensitivities

  •  
  • ℝ =

    set of real numbers

  •  
  • I =

    identity matrix

  •  
  • iB =

    imaginary axis reserved for Bloch’s theorem

  •  
  • is =

    imaginary axis reserved to the computation of sensitivities

  •  
  • m1 =

    mass corresponding to node 1

  •  
  • uR =

    reduced displacement vector

  •  
  • xi1 =

    imaginary part of x corresponding to first imaginary axis

  •  
  • xRe =

    real component of x

  •  
  • zη =

    multicomplex number of order η

  •  
  • 1 =

    set of complex numbers

  •  
  • 2 =

    set of bicomplex numbers

  •  
  • KR =

    reduced stiffness matrix

  •  
  • MR =

    reduced mass matrix

  •  
  • f(x) =

    arbitrary function evaluated at x

  •  
  • i1, i2, i3 =

    imaginary axes 1, 2, 3

  •  
  • u(x) =

    displacement vector at position x

  •  
  • Im1[] =

    part of a multicomplex number corresponding to the first imaginary axis

  •  
  • [K] =

    stiffness matrix

  •  
  • [M] =

    mass matrix

  •  
  • O() =

    error

  •  
  • {U} =

    displacement vector

  •  
  • [zη] =

    Cauchy-Riemann matrix notation of the multicomplex number z of order η

  •  
  • []j =

    matrix evaluated at the jth iteration of the QR algorithm

  •  
  • η =

    order of the multicomplex number

  •  
  • κ =

    wavenumber

  •  
  • λ =

    Eigenvalue

  •  
  • λ1 =

    Eigenvalue of the first degree-of-freedom

  •  
  • ω =

    frequency

References

1.
Sun
,
J.-H.
, and
Wu
,
T.-T.
,
2006
, “
Propagation of Surface Acoustic Waves Through Sharply Bent Two-Dimensional Phononic Crystal Waveguides Using a Finite-Difference Time-Domain Method
,”
Phys. Rev. B
,
74
(
17
), p.
174305
.
2.
Li
,
J.
,
Liu
,
Z.
, and
Qiu
,
C.
,
2006
, “
Negative Refraction Imaging of Acoustic Waves by a Two-Dimensional Three-Component Phononic Crystal
,”
Phys. Rev. B
,
73
(
5
), p.
054302
.
3.
Yang
,
S.
,
Page
,
J. H.
,
Liu
,
Z.
,
Cowan
,
M. L.
,
Chan
,
C. T.
, and
Sheng
,
P.
,
2004
, “
Focusing of Sound in a 3D Phononic Crystal
,”
Phys. Rev. Lett.
,
93
(
2
), p.
024301
.
4.
Kaya
,
O. A.
,
Cicek
,
A.
, and
Ulug
,
B.
,
2012
, “
Self-Collimated Slow Sound in Sonic Crystals
,”
J. Phys. D: Appl. Phys.
,
45
(
36
), p.
365101
.
5.
Yang
,
A.
,
Li
,
P.
,
Wen
,
Y.
,
Yang
,
C.
,
Wang
,
D.
,
Zhang
,
F.
, and
Zhang
,
J.
,
2015
, “
High-Q Cross-Plate Phononic Crystal Resonator for Enhanced Acoustic Wave Localization and Energy Harvesting
,”
Appl. Phys. Express
,
8
(
5
), p.
057101
.
6.
Zhang
,
X.
, and
Liu
,
Z.
,
2004
, “
Negative Refraction of Acoustic Waves in Two-Dimensional Phononic Crystals
,”
Appl. Phys. Lett.
,
85
(
2
), pp.
341
343
.
7.
Li
,
J.
, and
Chan
,
C. T.
,
2004
, “
Double-Negative Acoustic Metamaterial
,”
Phys. Rev. E
,
70
(
5
), p.
055602
.
8.
Fang
,
N.
,
Xi
,
D.
,
Xu
,
J.
,
Ambati
,
M.
,
Srituravanich
,
W.
,
Sun
,
C.
, and
Zhang
,
X.
,
2006
, “
Ultrasonic Metamaterials With Negative Modulus
,”
Nat. Mater.
,
5
(
6
), pp.
452
456
.
9.
Liu
,
Z.
,
Zhang
,
X.
,
Mao
,
Y.
,
Zhu
,
Y. Y.
,
Yang
,
Z.
, and
Chan
,
C. T.
,
2000
, “
Locally Resonant Sonic Materials
,”
Science
,
289
(
5485
), pp.
1734
1736
.
10.
Ding
,
Y.
,
Liu
,
Z.
,
Qiu
,
C.
, and
Shi
,
J.
,
2007
, “
Metamaterial With Simultaneously Negative Bulk Modulus and Mass Density
,”
Phys. Rev. Lett.
,
99
(
9
), p.
093904
.
11.
Mei
,
J.
,
Liu
,
Z.
,
Wen
,
W.
, and
Sheng
,
P.
,
2006
, “
Effective Mass Density of Fluid-Solid Composites
,”
Phys. Rev. Lett.
,
96
(
2
), p.
024301
.
12.
Zhao
,
D.
,
Wang
,
W.
,
Liu
,
Z.
,
Shi
,
J.
, and
Wen
,
W.
,
2007
, “
Peculiar Transmission Property of Acoustic Waves in a One-Dimensional Layered Phononic Crystal
,”
Phys. B: Condens. Matter
,
390
(
1–2
), pp.
159
166
.
13.
Shi
,
J.
,
Lin
,
S.-C. S.
, and
Huang
,
T. J.
,
2008
, “
Wide-Band Acoustic Collimating by Phononic Crystal Composites
,”
Appl. Phys. Lett.
,
92
(
11
), p.
111901
.
14.
Page
,
J. H.
,
Sheng
,
P.
,
Schriemer
,
H. P.
,
Jones
,
I.
,
Jing
,
X.
, and
Weitz
,
D. A.
,
1996
, “
Group Velocity in Strongly Scattering Media
,”
Science
,
271
(
5249
), pp.
634
637
.
15.
Chen
,
H.-T.
,
Padilla
,
W. J.
,
Zide
,
J. M. O.
,
Gossard
,
A. C.
,
Taylor
,
A. J.
, and
Averitt
,
R. D.
,
2006
, “
Active Terahertz Metamaterial Devices
,”
Nature
,
444
(
7119
), pp.
597
600
.
16.
Valencia
,
C.
,
Restrepo
,
D.
,
Mankame
,
N. D.
,
Zavattieri
,
P. D.
, and
Gomez
,
J.
,
2019
, “
Computational Characterization of the Wave Propagation Behavior of Multi-Stable Periodic Cellular Materials
,”
Extrem. Mech. Lett.
,
33
, p.
100565
.
17.
Meaud
,
J.
,
2018
, “
Multistable Two-Dimensional Spring-Mass Lattices With Tunable Band Gaps and Wave Directionality
,”
J. Sound Vib.
,
434
, pp.
44
62
.
18.
Schumacher
,
A.
, and
Vietor
,
T.
,
2018
,
Advances in Structural and Multidisciplinary Optimization
,
Springer International Publishing
,
Cham
.
19.
Arretche
,
I.
, and
Matlack
,
K. H.
,
2019
, “
Experimental Testing of Vibration Mitigation in 3D-Printed Architected Metastructures
,”
ASME J. Appl. Mech.
,
86
(
11
), p.
111008
.
20.
Injeti
,
S. S.
,
Celli
,
P.
,
Bhattacharya
,
K.
, and
Daraio
,
C.
,
2021
, “
Tuning Acoustic Impedance in Load-Bearing Structures
,”
arXiv preprint
, pp.
1
29
.
21.
Sigmund
,
O.
, and
Søndergaard Jensen
,
J.
,
2003
, “
Systematic Design of Phononic Band–Gap Materials and Structures by Topology Optimization
,”
Philos. Trans. R. Soc. London, Ser. A
,
361
(
1806
), pp.
1001
1019
.
22.
Martins
,
J. R. R. A.
, and
Hwang
,
J. T.
,
2013
, “
Review and Unification of Methods for Computing Derivatives of Multidisciplinary Computational Models
,”
AIAA J.
,
51
(
11
), pp.
2582
2599
.
23.
Wittrick
,
W. H.
,
1962
, “
Rates of Change of Eigenvalues, With Reference to Buckling and Vibration Problems
,”
J. R. Aeronaut. Soc.
,
66
(
621
), pp.
590
591
.
24.
Lancaster
,
P.
,
1964
, “
On Eigenvalues of Matrices Dependent on a Parameter
,”
Numer. Math.
,
6
(
1
), pp.
377
387
.
25.
Fox
,
R. L.
, and
Kapoor
,
M. P.
,
1968
, “
Rates of Change of Eigenvalues and Eigenvectors
,”
AIAA J.
,
6
(
12
), pp.
2426
2429
.
26.
Rogers
,
L. C.
,
1970
, “
Derivatives of Eigenvalues and Eigenvectors
,”
AIAA J.
,
8
(
5
), pp.
943
944
.
27.
Plaut
,
R. H.
, and
Huseyin
,
K.
,
1973
, “
Derivatives of Eigenvalues and Eigenvectors in Non-Self-Adjoint Systems
,”
AIAA J.
,
11
(
2
), pp.
250
251
.
28.
Rudisill
,
C. S.
,
1974
, “
Derivatives of Eigenvalues and Eigenvectors for a General Matrix
,”
AIAA J.
,
12
(
5
), pp.
721
722
.
29.
Jankovic
,
M. S.
,
1994
, “
Exact Nth Derivatives of Eigenvalues and Eigenvectors
,”
J. Guid. Control Dyn.
,
17
(
1
), pp.
136
144
.
30.
Garg
,
S.
,
1973
, “
Derivatives of Eigensolutions for a General Matrix
,”
AIAA J.
,
11
(
8
), pp.
1191
1194
.
31.
Lin
,
R. M.
,
Mottershead
,
J. E.
, and
Ng
,
T. Y.
,
2020
, “
A State-of-the-Art Review on Theory and Engineering Applications of Eigenvalue and Eigenvector Derivatives
,”
Mech. Syst. Signal Process.
,
138
, p.
106536
.
32.
Garza
,
J.
, and
Millwater
,
H.
,
2015
, “
Multicomplex Newmark-Beta Time Integration Method for Sensitivity Analysis in Structural Dynamics
,”
AIAA J.
,
53
(
5
), pp.
1188
1198
.
33.
Lim
,
K. B.
,
Junkins
,
J. L.
, and
Wang
,
B. P.
,
1987
, “
Re-Examination of Eigenvector Derivatives
,”
J. Guid. Control Dyn.
,
10
(
6
), pp.
581
587
.
34.
Iott
,
J.
,
Haftka
,
R. T.
, and
Adelman
,
H. M.
,
1985
,
Selecting Step Sizes in Sensitivity Analysis by Finite Differences
,
NASA Scientific and Technical Information Branch
, NASA Technical Memorandum 86382.
35.
Seyranian
,
A. P.
,
Lund
,
E.
, and
Olhoff
,
N.
,
1994
, “
Multiple Eigenvalues in Structural Optimization Problems
,”
Struct. Optim.
,
8
(
4
), pp.
207
227
.
36.
Squire
,
W.
, and
Trapp
,
G.
,
1998
, “
Using Complex Variables to Estimate Derivatives of Real Functions
,”
SIAM Rev.
,
40
(
1
), pp.
110
112
.
37.
Wang
,
B. P.
, and
Apte
,
A. P.
,
2006
, “
Complex Variable Method for Eigensolution Sensitivity Analysis
,”
AIAA J.
,
44
(
12
), pp.
2958
2961
.
38.
Voorhees
,
A.
,
Bagley
,
R.
,
Millwater
,
H.
, and
Golden
,
P.
,
2009
, “
Application of Complex Variable Methods for Fatigue Sensitivity Analysis
,”
50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference
,
Palm Springs, CA
,
May 4–7
.
39.
Bloch
,
F.
,
1929
, “
Uber Die Quantenmechanik Der Elektronen in Kristallgittern
,”
Z. Phys.
,
52
(
7–8
), pp.
555
600
.
40.
Floquet
,
G.
,
1883
, “
Sur Les Équations Différentielles Linéaires à Coefficients Périodiques
,”
Ann. Sci. l’École Norm. Supér.
,
12
, pp.
47
88
.
41.
Lantoine
,
G.
,
Russell
,
R. P.
, and
Dargent
,
T.
,
2012
, “
Using Multicomplex Variables for Automatic Computation of High-Order Derivatives
,”
ACM Trans. Math. Softw.
,
38
(
3
), pp.
1
21
.
42.
Price
,
G. B.
,
1991
,
An Introduction to Multicomplex Spaces and Functions
,
Dekker
,
New York, NY
.
43.
Voorhees
,
A.
,
Millwater
,
H.
, and
Bagley
,
R.
,
2011
, “
Complex Variable Methods for Shape Sensitivity of Finite Element Models
,”
Fin. Elem. Anal. Des.
,
47
(
10
), pp.
1146
1156
.
44.
Voorhees
,
A.
,
Millwater
,
H.
,
Bagley
,
R.
, and
Golden
,
P.
,
2012
, “
Fatigue Sensitivity Analysis Using Complex Variable Methods
,”
Int. J. Fatigue
,
40
, pp.
61
73
.
45.
Montoya
,
A.
, and
Millwater
,
H.
,
2017
, “
Sensitivity Analysis in Thermoelastic Problems Using the Complex Finite Element Method
,”
J. Therm. Stresses
,
40
(
3
), pp.
302
321
.
46.
Ramirez Tamayo
,
D.
,
Montoya
,
A.
, and
Millwater
,
H.
,
2018
, “
A Virtual Crack Extension Method for Thermoelastic Fracture Using a Complex-Variable Finite Element Method
,”
Eng. Fract. Mech.
,
192
, pp.
328
342
.
47.
Wagner
,
D.
,
Garcia
,
M. J.
,
Montoya
,
A.
, and
Millwater
,
H.
,
2019
, “
A Finite Element-Based Adaptive Energy Response Function Method for 2D Curvilinear Progressive Fracture
,”
Int. J. Fatigue
,
127
, pp.
229
245
.
48.
Ramirez-Tamayo
,
D.
,
Soulami
,
A.
,
Gupta
,
V.
,
Restrepo
,
D.
,
Montoya
,
A.
, and
Millwater
,
H.
,
2021
, “
A Complex-Variable Cohesive Finite Element Subroutine to Enable Efficient Determination of Interfacial Cohesive Material Parameters
,”
Eng. Fract. Mech.
,
247
, p.
107638
.
49.
Montoya
,
A.
,
Fielder
,
R.
,
Gomez-Farias
,
A.
, and
Millwater
,
H.
,
2015
, “
Finite-Element Sensitivity for Plasticity Using Complex Variable Methods
,”
J. Eng. Mech.
,
141
(
2
), p.
04014118
.
50.
Fielder
,
R.
,
Montoya
,
A.
,
Millwater
,
H.
, and
Golden
,
P.
,
2017
, “
Residual Stress Sensitivity Analysis Using a Complex Variable Finite Element Method
,”
Int. J. Mech. Sci.
,
133
, pp.
112
120
.
51.
Gomez-Farias
,
A.
,
Montoya
,
A.
, and
Millwater
,
H.
,
2015
, “
Complex Finite Element Sensitivity Method for Creep Analysis
,”
Int. J. Press. Vessel. Pip.
,
132–133
, pp.
27
42
.
52.
Monsalvo
,
J. F.
,
García
,
M. J.
,
Millwater
,
H.
, and
Feng
,
Y.
,
2017
, “
Sensitivity Analysis for Radiofrequency Induced Thermal Therapies Using the Complex Finite Element Method
,”
Finite Elem. Anal. Des.
,
135
, pp.
11
21
.
53.
Fielder
,
R.
,
Millwater
,
H.
,
Montoya
,
A.
, and
Golden
,
P.
,
2019
, “
Efficient Estimate of Residual Stress Variance Using Complex Variable Finite Element Methods
,”
Int. J. Press. Vessel. Pip.
,
173
, pp.
101
113
.
54.
Brillouin
,
L.
,
1953
,
Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices
,
Dover
,
New York
.
55.
Kantor
,
I. L.
,
Kantor
,
I. L.
, and
Solodovnikov
,
A. S.
,
1989
,
Hypercomplex Numbers: An Elementary Introduction to Algebras
,
Springer
,
New York
.
56.
Aguirre-Mesa
,
A. M.
,
Garcia
,
M. J.
, and
Millwater
,
H.
,
2020
, “
MultiZ: A Library for Computation of High Order Derivatives Using Multicomplex or Multidual Numbers
,”
ACM Trans. Math. Softw.
,
46
(
3
), pp.
1
30
.
57.
Francis
,
J. G. F.
,
1961
, “
The QR Transformation A Unitary Analogue to the LR Transformation—Part 1
,”
Comput. J.
,
4
(
3
), pp.
265
271
.
58.
Francis
,
J. G. F.
,
1962
, “
The QR Transformation–Part 2
,”
Comput. J.
,
4
(
4
), pp.
332
345
.
59.
Kublanovskaya
,
V. N.
,
1962
, “
On Some Algorithms for the Solution of the Complete Eigenvalue Problem
,”
USSR Comput. Math. Math. Phys.
,
1
(
3
), pp.
637
657
.
60.
Watkins
,
D. S.
,
1982
, “
Understanding the $QR$ Algorithm
,”
SIAM Rev.
,
24
(
4
), pp.
427
440
.
61.
Parlett
,
B. N.
,
2000
, “
The QR Algorithm
,”
Comput. Sci. Eng.
,
2
(
1
), pp.
38
42
.
62.
Watkins
,
D. S.
,
2008
, “
The QR Algorithm Revisited
,”
SIAM Rev.
,
50
(
1
), pp.
133
145
.
63.
Tzounas
,
G.
,
Dassios
,
I.
,
Liu
,
M.
, and
Milano
,
F.
,
2020
, “
Comparison of Numerical Methods and Open-Source Libraries for Eigenvalue Analysis of Large-Scale Power Systems
,”
Appl. Sci.
,
10
(
21
), p.
7592
.
64.
Davidson
,
E. R.
,
1975
, “
The Iterative Calculation of a Few of the Lowest Eigenvalues and Corresponding Eigenvectors of Large Real-Symmetric Matrices
,”
J. Comput. Phys.
,
17
(
1
), pp.
87
94
.
65.
Sorensen
,
D. C.
,
1992
, “
Implicit Application of Polynomial Filters in a k-Step Arnoldi Method
,”
SIAM J. Matrix Anal. Appl.
,
13
(
1
), pp.
357
385
.
66.
Knyazev
,
A. V.
,
2001
, “
Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method
,”
SIAM J. Sci. Comput.
,
23
(
2
), pp.
517
541
.
67.
Golub
,
G. H.
, and
Van Loan
,
C. F.
,
2013
,
Matrix Computations
,
The Johns Hopkins University Press
,
Baltimore, MD
.
68.
Trefethen
,
L. N.
, and
Bau
,
D.
, III
,
1997
,
Numerical Linear Algebra
,
Society for Industrial and Applied Mathematics
,
Philadelphia, PA
.
69.
Kundur
,
P.
,
Rogers
,
G. J.
,
Wong
,
D. Y.
,
Wang
,
L.
, and
Lauby
,
M. G.
,
1990
, “
A Comprehensive Computer Program Package for Small Signal Stability Analysis of Power Systems
,”
IEEE Trans. Power Syst.
,
5
(
4
), pp.
1076
1083
.
70.
Milano
,
F.
, and
Dassios
,
I.
,
2017
, “
Primal and Dual Generalized Eigenvalue Problems for Power Systems Small-Signal Stability Analysis
,”
IEEE Trans. Power Syst.
,
32
(
6
), pp.
4626
4635
.
71.
Denman
,
E. D.
, and
Beavers
,
A. N.
,
1976
, “
The Matrix Sign Function and Computations in Systems
,”
Appl. Math. Comput.
,
2
(
1
), pp.
63
94
.
72.
Deadman
,
E.
,
Higham
,
N. J.
, and
Ralha
,
R.
,
2013
, “Blocked Schur Algorithms for Computing the Matrix Square Root,”
Applied Parallel and Scientific Computing
,
Springer
,
Berlin/Heidelberg
, pp.
171
182
.
73.
Hussein
,
M. I.
, and
Frazier
,
M. J.
,
2010
, “
Band Structure of Phononic Crystals with General Damping
,”
J. Appl. Phys.
,
108
(
9
), p.
093506
.
74.
Millwater
,
H. R.
,
Wagner
,
D.
,
Garza
,
J. E.
,
Baines
,
A.
, and
Lovelady
,
K.
,
2013
, “
Application of the Complex Variable Finite Element Method, ZFEM, to Structural Mechanics Sensitivity Analysis
,”
54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference
,
Boston, MA
,
Apr. 8–11
.
75.
Al-Mohy
,
A. H.
, and
Higham
,
N. J.
,
2010
, “
A New Scaling and Squaring Algorithm for the Matrix Exponential
,”
SIAM J. Matrix Anal. Appl.
,
31
(
3
), pp.
970
989
.
76.
Higham
,
N. J.
,
2005
, “
The Scaling and Squaring Method for the Matrix Exponential Revisited
,”
SIAM J. Matrix Anal. Appl.
,
26
(
4
), pp.
1179
1193
.

Supplementary data