Arbitrary-Order Sensitivity Analysis in Phononic Metamaterials Using the Multicomplex Taylor Series Expansion Method Coupled With Bloch’s Theorem

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 [23–29]. 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 [31–34]. 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⁻¹⁰. 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., i₁, i₂, and i₃. 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 [46–48], 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.

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⁻¹⁰ times, the smallest problem parameter will produce a truncation error below the machine precision.

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.

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⁻¹⁰) times the nominal value of the variable of interest.

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.

The resulting system after applying Bloch’s theorem is nonself-adjoint. Therefore, the QR algorithm is used to solve the eigenvalue problem [57–59]. The QR algorithm is recognized as one of the most important algorithms for eigenvalue computations due to its robustness [60–62]. 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.

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.

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.

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 N³, 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.

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.

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 m₂, 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 i_B. With this in mind, the CR matrix representation for the input parameters is as follows:

We implemented the multicomplex-step perturbation method to compute the dispersion relation in Eq. (24) and its sensitivities using the following input values, m₁ = 1, m₂ = 2, c = 1, and a = 1. With perturbation magnitude, h = 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⁻¹⁴ for the dispersion relation, 6 × 10⁻¹², and 10⁻¹⁰ for the first- and second-order sensitivities with respect to m₂, respectively, and 9 × 10⁻¹¹ for the second-order mixed sensitivity with respect to m₂ 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⁻⁹.

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 m₂. 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.

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.

A rectangular cell was analyzed that consisted of four masses m connected with horizontal and vertical springs with stiffness c_x and c_y, respectively. Springs with effective stiffness c_d 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.