# A computationally efficient $k(\omega)$-spectral form for partial dispersion analyses within the wave finite element framework

Alvaro Gavilán Rojas, Qinghua Zhang, Christophe Droz

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alvaro Gavilán Rojas, Qinghua Zhang, Christophe Droz. A computationally efficient $k(\omega)$-spectral form for partial dispersion analyses within the wave finite element framework. Journal of Sound and Vibration, 2024, 593, pp.118652. 10.1016/j.jsv.2024.118652. hal-04750240

HAL Id: hal-04750240

https://hal.science/hal-04750240v1

Submitted on 23 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A computationally efficient $k(\omega)$-spectral form for partial dispersion analyses within the wave finite element framework

Alvaro Gavilán Rojas<sup>a,b</sup>, Qinghua Zhang<sup>a</sup>, Christophe Droz<sup>a</sup>

<sup>a</sup>Univ. Gustave Eiffel, Inria, COSYS-SII, I4S, Rennes, France

<sup>b</sup>Centre de Recherche Acoustique-Signal-Humain, Université de Sherbrooke, Sherbrooke, Canada

## Abstract

This paper addresses the computation of frequency-dependent dispersion curves (i.e., $k(\omega)$) and wave modes within the framework of the Wave Finite Element Method (WFEM) and in the context of high-dimensional periodic unit cell models. Numerous applications, ranging from phononics to vibroacoustics, now rely on dispersion analyses or wave expansion over a subset of eigensolutions –complex wavenumbers and Bloch waves– resulting from the resolution of an eigenvalue problem with a T-palindromic quadratic structure (T-PQEP). To exploit the structure of finite element models, various structure-preserving linearizations such as the Zhong-Williams and the $(S+S^{-1})$-transform have already been developed to achieve partial wave resolution of large T-PQEP, primarily targeting the dominating (least decaying) waves. In this paper we derive an alternative linearization of the T-PQEP for the $k(\omega)$ problem, which leads to enhanced targeting of the eigenvalues around the unit circle and reduces the inaccuracies induced by root multiplicity. A specific form of the problem is then proposed as an optimal compromise between ease of implementation, numerical stability, convergence and accuracy enhancement. The performance of our proposed linearization is compared against existing ones across various iterative eigensolvers, since the generalized eigenvalue problems involve complex non-hermitian matrices, which are not extensively included in eigensolvers. Results indicate that the proposed linearization should be favored for the WFEM, as it provides numerical enhancements in dispersion and wave vectors computation for large eigenvalue problems, as well as for further wave expansion applications.

**Keywords:** Wave propagation, Wave Finite Element Method, Periodic structures, Palindromic quadratic eigenvalue problem

## 1. Introduction

Periodic media are widely studied in physics, whether in the field of electromagnetism and acoustics –photonic and phononic crystals [1, 2]– or structural mechanics –composite materials, aircraft rib-skin structures [3], architected materials–. Their applications, at various scales, are diverse in engineering, since the advent of additive manufacturing promotes their industrial and commercial deployment. Some examples are photonic crystal fibers (PCF) [4], gradient index (GRIN) phononic crystals [5] (acoustic metamaterials) and locally resonant materials (LRM) [6]. Their study is realized mostly through numerical simulations, and it is of vital importance to develop a mature comprehension of the singular properties that arise in these complex structures which lead ultimately to the design of lightweight materials and structures as well as controlling wave propagation or vibroacoustic properties efficiently.

The numerical analysis of periodic media is mainly based on the examination of dispersion curves or band diagrams. These curves reveal band-gaps, where elastic waves are strongly attenuated, or do not propagate at all at certain frequencies. This analysis is performed for a single unit cell, representative of the entire periodic structure. Dispersion curves are often constructed by defining a wave vector – in the case of a 1D waveguide, the wavenumber $k \in \mathbb{C}$ – and solving an eigenvalue problem to find the $\omega$ frequencies. This approach is called $\omega(k)$ and it is widespread in the guided wave literature in forms such as the *plane wave expansion* (PWE) [7], or the *reduced Bloch mode expansion* (RBME) [8, 9, 10], both in the photonic and
phononic crystals field. However, the $\omega(k)$ problem has a limitation, since it does not take into account the wave spatial attenuation associated to the imaginary part of the wavenumber, which is essential to analyze dynamic motion in damped waveguides, both in a wave approach (complex band diagrams and Bloch modes [1], wave scattering through a coupling element [11]) and in a vibration approach (frequency response functions of finite waveguides [12, 13] and waveguide assemblies with joints, defects [14, 15] or piezoelectric actuators [1, 16], transient analysis [17] and vibroacoustic simulations [18]).

To deal with this wave spatial attenuation, a common strategy is to find the complex wavenumbers $k$ by fixing the frequency $\omega$, and so, generate complex band diagrams. This approach is called $k(\omega)$. Among the existing methods for wave propagation problems in the literature, the Wave Finite Element Method (WFEM) is capable of finding – following its initial formulation [19] –, from a finite element description of a unit cell, these complex wavenumbers $k$. The main equation of the WFEM is a T-palindromic quadratic eigenvalue problem (T-PQEP) [20] with the eigenpair $(\lambda, \phi)$. This eigenvalue problem has a symplectic spectrum, i.e., the eigenvalues come in reciprocal pairs $(\lambda, 1/\lambda)$ that represent positive- and negative-going free wave modes. For the many applications of the WFEM mentioned before, the eigenpairs for which $|\lambda| \approx 1$ – i.e., the propagative waves – are of interest. In the complex plane, these solutions lie near the unit circle. But the spectrum sparsity that comes from large differences between the eigenvalues of the propagative and evanescent waves (concentrated around $|\lambda| \ll 1$ and $|\lambda| \gg 1$) poses a challenge to the numerical solution of the eigenvalue problem [21, 22]. Nonetheless, the eigensolutions or free wave modes can be used to describe the dynamic motion of the periodic waveguide, and – in the same way as modal decomposition for modal analysis – a wave basis truncation can be done [11, 23] to approximate this motion with a selected set of waves – referred to as the wave approach (WA) or wave expansion – keeping the propagative modes and neglecting the majority of evanescent waves that do not contribute to the dynamic motion of the waveguide. Moreover, for large numerical models, the abundance of evanescent modes could significantly alter the numerical precision of the eigensolutions. Thus, it is preferable to perform a partial eigensolving instead of a full eigensolving that implies unnecessary higher computational burden.

The T-PQEP has received much attention by researchers on vibration of fast trains [20, 24, 25, 26] that worked on special linearizations of the eigenvalue problem for better accuracy and convergence [27, 28]. Although – to the authors’ knowledge – there is no specific method to find the solutions of this eigenvalue problem around the unit circle in the complex plane, a workaround is to obtain the partial spectrum with either the linearization proposed by Zhong and Williams [29] or the $(S + S^{-1})$-transform proposed by Huang et al. [30]. These formulations are based on the same linearization of the eigenvalue problem, i.e., $\mu = \lambda + 1/\lambda$, which preserves a symplectic structure. The stated purpose in the respective papers and preceding work [28, 27] were to improve the numerical condition of the eigenvalue problem (precision, convergence and computation time), developing structure preserving algorithms (SPA) that take benefit of the symplectic properties of the problem. These two linearizations are used constantly in the WFEM literature [22, 12, 31, 32] to build the mentioned set of selected waves, but apart from the GTSHIRA algorithm [30, 33] based upon the Huang et al. linearization, they were developed for a full eigenvalue problem resolution. Although very similar in their derivation, these linearizations are slightly different in their final form and the original eigenpair $(\lambda, \phi)$ recovery procedure, because Zhong-Williams uses a singular value decomposition (SVD) [21] to perform a linear combination of the $\mu$-related eigenvectors while Huang et al. propose a more direct linear combination using the $\lambda$ eigenvalues. However, these two linearizations provide repeated $\mu$ eigenvalues, and recently this was found to be prone to numerical errors for waveguides with cyclic symmetry in their cross-section [34] ($\lambda$ multiplicities higher than one).

These two linearizations in $\mu = \lambda + 1/\lambda$ are commonly used with partial iterative eigensolving techniques, but their structure (generalized eigenvalue problem (GEVP) with both complex non-hermitian matrices) is not commonly suitable for existing commercial eigensolvers. For instance, some eigensolvers such as `scipy sparse linalg eigs()` in Python and `geneig()` from the KrylovKit.jl package in Julia require a positive definite complex-hermitian matrix. Specialized eigensolvers for sparse matrices like SLEPc (Scalable Library for Eigenvalue Problem Computations) can address this particular GEVP, yet the matrices of the GEVP discussed in these article are not sparse but dense. General use eigensolvers capable to address the problem will be used in these article – namely, `eigs()` in MATLAB, `eigs()` from the Arpack.jl package in Julia and `partiaschur()/partialeigen()` from the ArnoldiMethod.jl package in Julia – in order to
evaluate their performance and robustness across the existent and the proposed T-PQEP linearizations. To summarize, the computation of complex wavenumbers for large-scale periodic structures still faces some unresolved numerical issues to this day.

In this paper, a novel $\mu(\phi)$ linearization is established for the partial resolution of the T-PQEP. The parameter $\phi$ represents a rotation in the complex plane and aims, firstly, to avoid the condition of repeated $\mu$ eigenvalues, and secondly, to increase convergence by closing the gap between the $\mu(\phi)$ partial spectrum and the desired $\lambda$ spectrum. This rotation parameter $\phi$ is introduced here in the context of $k(\omega)$-spectral forms of the wave propagation problem, for which the eigenvalues with $|\lambda| \approx 1$ represent the propagative waves of interest. Moreover, the $\eta$ linearization – obtained for the rotation of $\phi = \pi/2$ – is developed, since it can be described and implemented in a concise way. The proposed $\mu(\phi)$ linearization ensures accuracy and robustness to different iterative eigensolvers. In Section 2, a state-of-the-art review of the WFEM quadratic eigenvalue problem and their linearizations are presented. In Section 3, the proposed $\mu(\phi)$ and $\eta$ linearizations, as well as their derivation are explained in detail. In Section 4, two numerical examples are used to show the advantages – both in convergence and precision – of the proposed linearization over the existent ones, comparing also the three mentioned partial eigensolvers from different scientific computing packages. Conclusions are presented in Section 5.

## 2. The T-palindromic quadratic eigenvalue problem (T-PQEP)

The Wave Finite Element Method (WFEM) [19, 13], begins with a numerical finite element model description of the unit cell of a periodic waveguide, involving its mass, damping and stiffness matrices $M$, $C$, and $K$, all being real symmetric and of size $n_d$, the number of degrees of freedom of the unit cell. The size of the unit cell is $d$, that represents the distance between its left and right boundaries, which have the same mesh in order to numerically assemble the unit cells forming a waveguide.

<figure><img src="image_3.png" /><figcaption>Figure 1: Infinite 1D periodic waveguide.</figcaption></figure>

In harmonic regime, the equation of motion is written $\mathbb{D}\mathbf{q} = \mathbf{f}$, where $\mathbf{q} \in \mathbb{C}^{n_d}$ are the nodal displacements, $\mathbf{f} \in \mathbb{C}^{n_d}$ the nodal forces, and $\mathbb{D} = -\omega^2\mathbf{M} + j\omega\mathbf{C} + \mathbf{K}$ is the complex symmetric dynamic stiffness matrix with $\omega \in \mathbb{R}$ the angular frequency. The degrees of freedom of a unit cell are separated into three parts, $L$ and $R$ for those at left and right boundaries and $I$ for the inner ones. Accordingly, $\mathbb{D}$, $\mathbf{q}$ and $\mathbf{f}$ are partitioned and indexed by $L$, $R$ and $I$. Considering that external loads will be applied only at the side boundaries of the unit cell, i.e., $\mathbf{f}_I = \mathbf{0}$ for the internal nodes, the condensed equation of motion can be written as:

$$ \begin{bmatrix} \mathbf{D}_{LL} & \mathbf{D}_{LR} \\ \mathbf{D}_{RL} & \mathbf{D}_{RR} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_L \\ \mathbf{q}_R \end{Bmatrix} = \begin{Bmatrix} \mathbf{f}_L \\ \mathbf{f}_R \end{Bmatrix} \qquad (1) $$

where $\mathbf{D}_{jk} = \mathbb{D}_{jk} - \mathbb{D}_{jI}\mathbb{D}_{II}^{-1}\mathbb{D}_{Ik}$ $\forall j, k \in \{L, R\}$, $\mathbf{D}_{LL}$ and $\mathbf{D}_{RR}$ being complex symmetric and $\mathbf{D}_{LR} = \mathbf{D}_{RL}^T$.

An infinite waveguide can be assembled from the unit cell, as depicted in Figure 1. Free wave propagation is found when $\mathbf{q}_R = \lambda\mathbf{q}_L$ and $\mathbf{f}_R = -\lambda\mathbf{f}_L$, with $\lambda = \exp(ikd)$ being the propagation constant and $k \in \mathbb{C}$ the complex wavenumber. The complex displacement and force vectors at the interfaces are of size $n$, the number of degrees of freedom at the left and right boundaries. The internal force equilibrium at each of the interfaces of an assembled infinite waveguide, written $\lambda\mathbf{f}_L + \mathbf{f}_R = \mathbf{0}$, gives us the following quadratic eigenvalue problem (QEP):

$$ \left( \lambda \mathbf{D}_{LR} + (\mathbf{D}_{LL} + \mathbf{D}_{RR}) + \frac{1}{\lambda} \mathbf{D}_{RL} \right) \mathbf{\Phi} = \mathbf{0} \qquad (2) $$
The 2n eigenvalues $\lambda$ are the propagation constants, and the 2n eigenvectors $\phi \in \mathbb{C}^n$ are the nodal displacements at a side boundary from which the Bloch waves (wave modes of free wave propagation in periodic media) can be recovered.

Notice that replacing $\lambda$ by its reciprocal $1/\lambda$ and transposing Eq. (2) yields the same eigenvalue problem. Thus, Eq. (2) is called the T-palindromic quadratic eigenvalue problem [30, 35, 29] and its eigensolutions provide the propagation constants and their corresponding Bloch modes. This also means that the 2n eigenvalues occur in n pairs $(\lambda, 1/\lambda)$, and that they form a symplectic spectrum [36]. Half of the spectrum contains eigenvalues with $|\lambda| < 1$ associated to positive-going waves and the other half contains eigenvalues with $|\lambda| > 1$ associated to negative-going waves. Here and hereafter, damped waveguides are considered for the sake of simplicity, i.e., there are eigenvalues with $|\lambda| \approx 1$, but no eigenvalue with $|\lambda| = 1$.

To solve the T-PQEP (2), several linearizations have been proposed in the literature. It is worth briefly mentioning the transfer matrix approach, which is a standard eigenvalue problem with a symplectic transfer matrix, and the NL linearization that comes from the structural dynamics field [37, 38] and was adapted for the WFEM in different forms [19, 29, 39, 16]. However, these linearized forms involve $\lambda$ eigenvalues, and cannot be used with iterative eigensolvers to obtain just a few propagative solutions $|\lambda| \approx 1$. Two well-known linearizations present this advantage of partial eigensolving, and they will be discussed in the following.

## 2.1. Zhong & Williams linearization

Following the advances in the WFEM field, Zhong and Williams [29] proposed a linearization of the T-PQEP that preserves its symplectic structure. Although there is no mention about a partial resolution of the eigenvalue problem, it was immediately known [40] that this linearization allowed to obtain the propagative waves via partial eigensolving techniques [13, 21, 32, 34].

The Zhong-Williams GEVP is expressed as:

$$ \begin{bmatrix} \mathbf{D}_{LR} - \mathbf{D}_{RL} & -(\mathbf{D}_{LL} + \mathbf{D}_{RR}) \\ \mathbf{D}_{LL} + \mathbf{D}_{RR} & \mathbf{D}_{LR} - \mathbf{D}_{RL} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_L \\ \mathbf{q}_R \end{Bmatrix} = \mu \begin{bmatrix} \mathbf{0} & \mathbf{D}_{LR} \\ -\mathbf{D}_{RL} & \mathbf{0} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_L \\ \mathbf{q}_R \end{Bmatrix} \quad (3) $$

where $\mu = \lambda + 1/\lambda$, the propagation constants are given by $\lambda = \frac{\mu \pm \sqrt{\mu^2 - 4}}{2}$, and both matrices of Eq. (3) are complex skew-symmetric. It is important to note that the 2n $\mu$ eigenvalues give us 4n $\lambda$ eigenvalues, but, since each of the n symplectic pairs of eigenvalues $(\lambda, 1/\lambda)$ produce a single $\mu$, the n values of $\mu$ are repeated. Hence, each pair of the 2n eigenvectors associated to the same $\mu$ in Eq. (3) will be used in a linear combination to retrieve the eigenvectors $\phi$ of the original T-PQEP (2). These linear combinations are mentioned in the original paper [29] but are explained in detail by Duhamel et al. [13] and Waki et al. [21]. The latter proposes a singular value decomposition (SVD) to obtain the weights of the linear combination and correct some low-frequency issues related to nearly parallel eigenvectors.

Finally, it is important to note that this linearization $\mu = \lambda + 1/\lambda$ permits not only to preserve the symplectic structure of the T-PQEP, but to obtain the propagative waves ($|\lambda| \approx 1$) efficiently by using a partial eigensolver (e.g., ARPACK [41]) and extracting only a few solutions of interest, which are the $\mu$ eigenvalues of smallest absolute value.

## 2.2. Huang et al. linearization

Another linearization that is used in the WFEM community [22, 31, 42] is the one proposed by Huang et al. [30], which is founded on the previous works of Patel [28] and Lin [27]. It is based on the $(\mathbf{S} + \mathbf{S}^{-1})$-transform, where $\mathbf{S} = \mathbf{L}^{-1}\mathbf{M}$ is a symplectic matrix. To derive their proposed GEVP, the following linearization in $\lambda$ is needed:

$$ (\mathbf{M} - \lambda \mathbf{L}) \mathbf{z} = \left( \begin{bmatrix} \mathbf{D}_{RL} & \mathbf{0} \\ -(\mathbf{D}_{LL} + \mathbf{D}_{RR}) & -\mathbf{I} \end{bmatrix} - \lambda \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{D}_{LR} & \mathbf{0} \end{bmatrix} \right) \begin{Bmatrix} \mathbf{q}_L \\ \frac{1}{\lambda} \mathbf{D}_{RL} \mathbf{q}_L \end{Bmatrix} \quad (4) $$

With matrices $\mathbf{M}$ and $\mathbf{L}$ from Eq. (4) and the standard symplectic matrix $\mathbf{J}$, it can be verified that:
$$ \mathbf{M}\mathbf{J}\mathbf{M}^{\mathrm{T}} = \mathbf{L}\mathbf{J}\mathbf{L}^{\mathrm{T}} = \begin{bmatrix} \mathbf{0} & -\mathbf{D}_{RL} \\ \mathbf{D}_{LR} & \mathbf{0} \end{bmatrix} \quad \text{with} \quad \mathbf{J} = \begin{bmatrix} \mathbf{0} & \mathbb{I} \\ -\mathbb{I} & \mathbf{0} \end{bmatrix} \qquad (5) $$

which is the required condition for $\mathbf{M} - \lambda \mathbf{L}$ to be a symplectic pencil (see Patel definition 1.3 in Ref. [28]). Now, setting $\mathbf{z} = \mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s$ and rewriting the linearization in $\lambda$ and $1/\lambda$ as follows:

$$ \left\{ \begin{array}{l} \mathbf{M}\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s = \lambda\mathbf{J}\mathbf{L}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s \\ \mathbf{L}\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s = \frac{1}{\lambda}\mathbf{M}\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s \end{array} \right. \qquad (6) $$

using the identity of Eq. (5), left multiplying by $\mathbf{L}\mathbf{M}^{-1}$ the linearization in $1/\lambda$ and adding these two equations, the obtained matrix pencil is written as:

$$ (\mathbf{M}\mathbf{J}\mathbf{L}^{\mathrm{T}} + \mathbf{L}\mathbf{J}\mathbf{M}^{\mathrm{T}}) - \left(\lambda + \frac{1}{\lambda}\right) \mathbf{L}\mathbf{J}\mathbf{L}^{\mathrm{T}} \qquad (7) $$

Eq. (7) is multiplied by $-1$ to preserve the l.h.s. matrix of the Zhong-Williams GEVP. Finally, the Huang et al. GEVP is written as follows:

$$ \begin{bmatrix} \mathbf{D}_{LR} - \mathbf{D}_{RL} & -(\mathbf{D}_{LL} + \mathbf{D}_{RR}) \\ \mathbf{D}_{LL} + \mathbf{D}_{RR} & \mathbf{D}_{LR} - \mathbf{D}_{RL} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_R \\ -\mathbf{q}_L \end{Bmatrix} = \mu \begin{bmatrix} \mathbf{0} & \mathbf{D}_{RL} \\ -\mathbf{D}_{LR} & \mathbf{0} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_R \\ -\mathbf{q}_L \end{Bmatrix} \qquad (8) $$

where $\mu = \lambda + 1/\lambda$ and the propagation constants are given by $\lambda = \frac{\mu \pm \sqrt{\mu^2 - 4}}{2}$, as in the Zhong-Williams GEVP. The eigenvectors $\mathbf{z}_s = [\mathbf{z}_1^T, \mathbf{z}_2^T]^T$ were also derived from $\mathbf{z}_s = -\mathbf{L}^{-1/2} \mathbf{J} \mathbf{z}$, and both matrices of the GEVP are complex skew-symmetric.

There are only subtle differences between the Zhong-Williams and Huang et al. GEVPs of Eq. (3) and Eq. (8): The eigenvectors change their indexes and sign, and the r.h.s. matrix change the $L$ and $R$ indexes. Nonetheless, the advantage of this second linearization is that the original eigenvectors $\phi$ of the T-PQEP (2) are recovered directly from a linear combination (see Theorem 2.2 in Ref. [30]) between the two halves of $\mathbf{z}_s$, i.e., $\mathbf{z}_1$ and $\mathbf{z}_2$. The eigenvectors $\mathbf{z}_1 + \lambda \mathbf{z}_2$ and $\mathbf{z}_1 + \frac{1}{\lambda} \mathbf{z}_2$ are associated to the $1/\lambda$ and $\lambda$ eigenvalues respectively.

### 3. A novel $\mu(\phi)$ linearization for partial eigensolutions of the T-PQEP

The proposed linearization of the T-PQEP (2) derived here aims to enhance convergence and precision of the eigenpairs computation when a partial eigensolver is used. Since the eigenvalues of interest lie near the unit circle in the complex plane ($|\lambda| \approx 1$), the $\mu(\phi)$ linearization is based on the standard Huang et al. linearization that is used in the WFEM literature. In fact, the $\phi$ parameter of the $\mu(\phi)$ linearization represents a rotation of the $\mu = \lambda + 1/\lambda$ function in the complex plane, and if no rotation is applied, the $\mu(\phi = 0)$ linearization is equivalent to the standard Huang et al. one. Different rotation values ($\phi \neq 0$) permit to explore the influence of $\phi$ finding the solutions of the T-PQEP, and the particular rotation of $\phi = \pi/2$ was found to present advantageous conditions compared with the standard ($\phi = 0$) linearizations.

#### 3.1. Physical motivation for the $\phi$ rotation

The primary insight to develop the $\mu(\phi)$ linearization was to understand how propagative waves in elastic waveguides evolve in the frequency domain by looking at their propagation constants $\lambda \in \mathbb{C}$ located in the complex plane. To recall some physical meanings, the expression for the propagation constant is $\lambda = \exp(ikd)$, with $k \in \mathbb{C}$, so that the real part of the wavenumber is $\Re(k) = \arg(\lambda)/d$ and the imaginary part of the wavenumber is $\Im(k) = -\ln(\lambda)/d$. The imaginary part $\Im(k)$ is related to the spatial decay $\gamma = \exp(-\Im(k)d)$, whereas the real part $\Re(k)$ is related to the wavelength $\Lambda = 2\pi/\Re(k)$.

In non-dispersive media, it can be seen for propagative elastic waves that the smaller their frequency, the larger their wavelength, because the constant phase velocity $c_p$ is expressed as $c_p = \Lambda f$. This can be understood for a single wave mode at $f = 0$ Hz as if there was no motion, and so the wavelength $\Lambda$ is
infinite. The same applies for dispersive media at low or cut-on frequencies. Since the wavelength shortens as frequency increases, $\Re(k)$ begins to increase away from zero, and thus, $\arg(\lambda)$ is also near zero at cut-on frequencies. Hence, by inference, the first propagation constants $\lambda$ at low frequencies will appear at the point 1 in the complex plane. A typical evolution through the frequency domain at low – or cut-on – frequencies of a propagative wavemode can be seen in Figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2: Propagation constants of a propagative mode evolving in the frequency domain.</figcaption></figure>

It has been shown that both Zhong-Williams and Huang et al. use the $\mu = \lambda + 1/\lambda$ linearization. To form a set of selected waves with mainly propagative waves ($|\lambda| \approx 1$), it is common to use this linearization and search for the first few eigensolutions with the smallest $|\mu|$. Aiming to understand the reason why it works, the function $|\mu| = |\lambda + 1/\lambda|$ was plotted on the complex plane (see Figure 3a). It was noticed that, although the unit circle is entirely covered by this function before reaching $|\mu| = 2$, the global minima of the function are located at the points $-i$ and $i$.

<figure><img src="image_5.png" /><figcaption>Figure 3: Contour plots of the $|\mu|$ and $|\eta|$ functions of $\lambda$. The $|\eta|$ function can be seen as the $|\mu|$ function rotated by an angle of $\phi = \pi/2$ in the complex plane. The sought propagative solutions shown in Figure 2 are marked here as black diamonds. A global minimum of the $|\eta|$ function is coincident with these propagative solutions at the point 1 in the complex plane.</figcaption></figure>

Ideally, these global minima should have been at the points 1 and -1, because it has been shown in Figure 2 that propagative waves begin to evolve from their origin at the point 1. This was the main
motivation for the proposed $\mu(\phi)$ linearization. The $\phi$ parameter can be seen as a rotation of the $|\mu|$ function in the complex plane. For instance, if $\phi = \pi/2$ – which is called the $\eta$ linearization in this paper –, a rotation of $\pi/2$ in clockwise sense is applied to the $|\mu|$ function. The $|\eta|$ function is shown in Figure 3b. This rotated function $|\eta|$ has now its global minima at the points 1 and -1. By approaching the targeted $|\mu(\phi)|$ values to the partial $\lambda$ spectrum of our interest, improved convergence should be expected within the iterative eigensolver, and ultimately, enhanced accuracy of the eigenpairs.

## 3.2. Derivation of the $\mu(\phi)$ linearization

Based on the previous works of Lin [27], Patel [28] and Huang et al. [30], the proposed $\mu(\phi)$ linearization of the T-PQEP (2) is derived here. From the symplectic pencil $M - \lambda L$ of Eq. (4), the ML linearization $Mz = \lambda Lz$ can be rewritten as $L^{-1}Mz = \lambda \mathbb{I}z$ or $M^{-1}Lz = \frac{1}{\lambda} \mathbb{I}z$. With $S = L^{-1}M$, the so-called $(S + S^{-1})$-transform would consist on the matrix pencil $S + S^{-1} - \mu \mathbb{I}$ with $\mu = \lambda + 1/\lambda$. Here, a different linear combination is introduced, by stating $\lambda^* = e^{i\phi}\lambda$ (i.e., the $\lambda$ eigenvalues are rotated in the complex plane by an angle $\phi$ in the anti-clockwise sense and denoted as $\lambda^*$) and developing $\mu^* = \lambda^* + 1/\lambda^*$ (i.e., the $\mu$ eigenvalues are rotated in the complex plane by an angle $\phi$ in the clockwise sense and denoted as $\mu^*$). The $\mu(\phi)$ linearization can be written as:

$$ \mu(\phi) = \lambda e^{i\phi} + \frac{e^{-i\phi}}{\lambda} \qquad (9) $$

From the roots of this equation, the propagation constants $\lambda$ can be retrieved as:

$$ \lambda = \frac{\mu \pm \sqrt{\mu^2 - 4}}{2e^{i\phi}} \qquad (10) $$

Notice that, if a first root corresponds to $\lambda = \frac{\mu + \sqrt{\mu^2 - 4}}{2e^{i\phi}}$, then the reciprocal propagation constant from the symplectic pair $\frac{1}{\lambda} = \frac{\mu - \sqrt{\mu^2 - 4}}{2e^{i\phi}}e^{i2\phi}$ is the second root of the same equation multiplied by $e^{i2\phi}$, i.e, the second root is not anymore the reciprocal eigenvalue (obtained with a $\mu = \lambda + 1/\lambda$ linearization) but rather an out of phase reciprocal eigenvalue. This means that the $\mu(\phi)$ linearization will produce a $\lambda$ eigenvalue that is of our interest, and a fictitious, out of phase reciprocal eigenvalue $e^{i2\phi}/\lambda$ with no physical meaning. The fictitious eigenvalue will be easily filtered out because both $1/\lambda$ and $e^{i2\phi}/\lambda$ have the same norm, and for eigenvalues with the same norm, the one associated with a physically meaningful eigenvector will be the correct one. This condition will be advantageous to recover the $\lambda$ related eigenvectors $\phi$, as it will be seen in the next subsection.

To obtain the corresponding matrices of the GEVP, the obtained linearization of the eigenvalue problem is developed as:

$$ (e^{i\phi}\mathbf{S} + e^{-i\phi}\mathbf{S}^{-1})\mathbf{z} = \left(\lambda e^{i\phi} + \frac{e^{-i\phi}}{\lambda}\right)\mathbb{I}\mathbf{z} \qquad (11) $$

which can be further developed as:

$$ (e^{i\phi}\mathbf{L}^{-1}\mathbf{M} + e^{-i\phi}\mathbf{M}^{-1}\mathbf{L})\mathbf{z} = \mu(\phi)\mathbb{I}\mathbf{z} \qquad (12) $$

If $\mathbf{L}^{-1}$ is taken out as a factor and $\mathbf{z} = \mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s$ is set, the obtained linearization is:

$$ \mathbf{L}^{-1}(e^{i\phi}\mathbf{M} + e^{-i\phi}\mathbf{L}\mathbf{M}^{-1}\mathbf{L})(\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s) = \mu(\phi)\mathbb{I}(\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s) \qquad (13) $$

$$ (e^{i\phi}\mathbf{M}\mathbf{J}\mathbf{L}^{\mathrm{T}} + e^{-i\phi}\mathbf{L}(\mathbf{L}^{-1}\mathbf{M})^{-1}\mathbf{J}\mathbf{L}^{\mathrm{T}})\mathbf{z}_s = \mu(\phi)\mathbf{L}\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s \qquad (14) $$

With $\mathbf{S} = \mathbf{L}^{-1}\mathbf{M}$ being symplectic, the property $\mathbf{S}^{-1} = \mathbf{J}\mathbf{S}^{\mathrm{T}}\mathbf{J}^{-1}$ can be used to preserve the symplectic structure as follows:

$$ (e^{i\phi}\mathbf{M}\mathbf{J}\mathbf{L}^{\mathrm{T}} + e^{-i\phi}\mathbf{L}\mathbf{J}(\mathbf{M}^{\mathrm{T}}\mathbf{L}^{-T})\mathbf{J}^{-1}\mathbf{J}\mathbf{L}^{\mathrm{T}})\mathbf{z}_s = \mu(\phi)\mathbf{L}\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s \qquad (15) $$
and simplifying, the obtained linearization is:

$$ (e^{i\phi}\mathbf{M}\mathbf{J}\mathbf{L}^{\mathrm{T}} + e^{-i\phi}\mathbf{L}\mathbf{J}\mathbf{M}^{\mathrm{T}})\mathbf{z}_s = \mu(\phi)\mathbf{L}\mathbf{J}\mathbf{L}^{\mathrm{T}}\mathbf{z}_s \quad (16) $$

Finally, multiplying Eq. (16) by $-1$ to preserve the r.h.s matrix of the Huang et al. linearization, the obtained GEVP is written as:

$$ \begin{bmatrix} \mathbf{D}_{LR}e^{-i\phi} - \mathbf{D}_{RL}e^{i\phi} & -(\mathbf{D}_{LL} + \mathbf{D}_{RR})e^{-i\phi} \\ (\mathbf{D}_{LL} + \mathbf{D}_{RR})e^{i\phi} & \mathbf{D}_{LR}e^{-i\phi} - \mathbf{D}_{RL}e^{i\phi} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_R \\ -\mathbf{q}_L \end{Bmatrix} = \mu(\phi) \begin{bmatrix} \mathbf{0} & \mathbf{D}_{RL} \\ -\mathbf{D}_{LR} & \mathbf{0} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_R \\ -\mathbf{q}_L \end{Bmatrix} \quad (17) $$

Eq. (17) represents the $\mu(\phi)$ linearization proposed in this paper. The introduced original contribution is the $\phi$ parameter representing a clockwise rotation of the $|\mu|$ function in the complex plane, as it was shown in Section 3.1. In the previous section it has been mentioned that selecting $\phi = \pi/2$ will be the more convenient value to improve computation performance. This particular rotation value was the initial idea behind the proposed $\mu(\phi)$ linearization, and it yields a particularly easy form of the eigenvalue problem. Denoting $\eta = -i\mu(\pi/2)$, the resulting eigenvalue problem is written as:

$$ \begin{bmatrix} -(\mathbf{D}_{LR} + \mathbf{D}_{RL}) & \mathbf{D}_{LL} + \mathbf{D}_{RR} \\ \mathbf{D}_{LL} + \mathbf{D}_{RR} & -(\mathbf{D}_{LR} + \mathbf{D}_{RL}) \end{bmatrix} \begin{Bmatrix} \mathbf{q}_R \\ -\mathbf{q}_L \end{Bmatrix} = \eta \begin{bmatrix} \mathbf{0} & \mathbf{D}_{RL} \\ -\mathbf{D}_{LR} & \mathbf{0} \end{bmatrix} \begin{Bmatrix} \mathbf{q}_R \\ -\mathbf{q}_L \end{Bmatrix} \quad (18) $$

Now the l.h.s. matrix of the GEVP is complex symmetric and the r.h.s. one is complex skew-symmetric. The $\eta$ eigenvalues and the corresponding propagation constants $\lambda$ are written as:

$$ \begin{cases} \eta = \lambda - \frac{1}{\lambda} \\ \lambda = \frac{\eta \pm \sqrt{\eta^2 + 4}}{2} \end{cases} \quad (19) $$

The procedure to recover the $\lambda$-related eigenvectors is generalized in the next paragraphs for the $\mu(\phi)$ linearization, and can be readily applied for this particular $\eta$ linearization. Finally, it is worth mentioning that the proposed linearization $\mu(\phi)$ also preserves the symplectic structure of the T-PQEP (2) – because of the symplectic property used to derive Eq. (15) –, not loosing the numerical improvements achieved with the Zhong-Williams and Huang et al. linearizations.

### 3.3. Eigenvector recovery

It has been shown that the $2n$ $\mu$ eigenvalues lead to $4n$ $\lambda$ eigenvalues. When $\phi = 0$ – i.e., the Huang et al. linearization – $\mu$ gives us $\lambda$ and $1/\lambda$, and so there are $n$ values of $\mu$ that are repeated. But, from Eq. (10), a $\mu(\phi)$ value will give us $\lambda$ and $e^{i2\phi}/\lambda$, i.e., an out of phase reciprocal eigenvalue. Overall, this linearization will produce $2n$ $\lambda$ eigenvalues that are of our interest, and $2n$ $e^{i2\phi}/\lambda$ out of phase reciprocal eigenvalues. This condition gives us the benefit of choosing directly the correct $2n$ eigenvalues, because the symplectic pairs will give us a set of out of phase eigenvalues with the same norm as their reciprocals.

To compare the two sets of eigenpairs, an eigenvector error metric $\epsilon$ is needed, aiming to evaluate the accuracy of the obtained eigenpairs from whichever linearization of the T-PQEP (2). A particular solution of the T-PQEP would consist in the eigenvalue $\lambda$ – the propagation constant –, and its associated eigenvector $\phi$ – a nodal displacement vector $\mathbf{q}_L$ –. Nodal forces $\mathbf{f}_L$ and $\mathbf{f}_R$ can then be derived from Eq. (1) as:

$$ \begin{cases} \mathbf{f}_L = (\mathbf{D}_{LL} + \lambda \mathbf{D}_{LR}) \mathbf{q}_L \\ \mathbf{f}_R = (\mathbf{D}_{RL} + \lambda \mathbf{D}_{RR}) \mathbf{q}_L \end{cases} \quad (20) $$

With the left and right force vectors, the internal equilibrium at the interfaces is expressed as $\lambda \mathbf{f}_L + \mathbf{f}_R = \mathbf{0}$. A nodal force vector $\mathbf{f}$ can also be written as a resultant force vector with the three Cartesian components $f_x, f_y, f_z$, denoted $\mathbf{f}^*$. The relative error $\epsilon$ can then be computed dividing the residual norm $\lambda \mathbf{f}_L^* + \mathbf{f}_R^*$ by the norm of the force vector $\mathbf{f}_R^*$, which is affected by $\lambda$ just once. The eigenvector error metric $\epsilon$ is expressed more generally with the relative error between the two original force vectors:

$$ \epsilon = \frac{\|\lambda \mathbf{f}_L + \mathbf{f}_R\|}{\|\mathbf{f}_R\|} \quad (21) $$
Now that an error metric has been established, the eigenvector recovery procedure can be detailed in the following steps:

1. Compute the $2r \le 2n$ eigenpairs $(\mu, z_s)$ with an eigensolver on Eq. (17).

2. Compute the $4r$ eigenpairs $(\lambda, z_s)$ using Eq. (10).

3. Form two sets of $r$ eigenvalues $\lambda$ which have the same norm and verify $|\lambda| < 1$ – i.e., positive-going waves – and their corresponding $\phi^+ = z_1 + \frac{1}{\lambda}z_2$ eigenvectors.

4. For each of the $r$ pairs of $\phi^+$ corresponding to the same $|\lambda|$, compare their errors on the T-PQEP (2) computing $\epsilon$ from Eq. (21).

5. Normalize the $r$ chosen $\phi^+$ eigenvectors with the smallest $\epsilon$.

6. Repeat steps 3, 4 and 5 for negative-going waves ($|\lambda| > 1$) and collect the $r$ chosen $\phi^-$ eigenvectors.

# 4. Numerical examples

In this section, two numerical models are presented and their dispersion characteristics are evaluated by means of the three different partial eigensolving strategies discussed in the paper, namely the Zhong-Williams, the Huang et. al. and the $\mu(\phi)$ linearizations. Moreover, three different partial eigensolvers are employed to assess their performance on the three different linearizations of the eigenvalue problem, namely the `eigs()` function in MATLAB, `eigs()` from the Arpack.jl package and `partialschur()`/`partialeigen()` from the ArnoldiMethod.jl package in Julia. It is important to mention that, with the three GEVP having complex non-hermitian matrices at both sides, there is a limitation in the use of some other partial eigensolvers, e.g., `scipy=sparse.lINALG.eigs()` in Python or `geneig()` from the KrylovKit.jl package in Julia.

## 4.1. Finned tube heat exchanger

The selected waveguide aims to demonstrate the influence of cyclic symmetry – and the corresponding increase in the $\lambda$ multiplicities – on the dispersion characteristics, for which the Zhong-Williams linearization has been improved recently [34], in order to avoid the cyclic symmetry-related issues. It is an industrial finned tube heat exchanger (U-shaped fins welded to the tube). This type of structures is likely to be inspected by guided waves for damage detection [43, 44], and its numerical modeling helps in the selection of wave modes and frequency ranges before in situ non-destructive testing. A numerical finite element model built with Gmsh [45] is depicted in Figure 4a.

<figure><img src="image_6.png" /><figcaption>Figure 4: Unit cell of the finned tube heat exchanger. The unit cell is meshed with quadratic tetrahedra and its size is $d = 5$ mm. The number of degrees of freedom at the interfaces is $n = 1314$. The tube is made of aluminum, with Young's modulus $E = 69$ GPa, Poisson's ratio $\nu = 0.35$, density $\rho = 2700$ kg/m$^3$ and loss factor $\eta = 0.001$.</figcaption></figure>
### 4.1.1. Dispersion characteristics

Solving Eq. (2) gives us the propagation constants $\lambda$ and the nodal displacements $\phi$. As it has been shown, from the propagation constants $\lambda = \exp(ikd)$, the wavenumbers $k \in \mathbb{C}$ are obtained. The imaginary part of the wavenumber $\Im(k) = -\ln(|\lambda|)/d$ is a measure of the attenuation per unit length called *attenuation constant* [46], it is denoted $\alpha$ and measured in nepers per meter (1 Np $\approx$ 8,686 dB). To understand its physical meaning, it can be seen that a linear attenuation $\alpha = 10$ dB/m is equivalent to a reduction in nodal displacements to 31.6% of the displacements one meter forward. Figure 4b shows the phase velocities $c_p = \omega/\Re(k)$ of the first 40 Bloch modes as well as their degree of attenuation, colored by the attenuation constant. It can be appreciated that the modes corresponding to the longitudinal fins, which are triggered and get converted around 2.5 to 3.5 kHz, are slightly less propagative than the tube modes. The transition between $\alpha$ values is strong during these triggers or cut-ons, and the rest of the evanescent waves do not propagate beyond one or two unit cells.

From the $\phi$ nodal displacements, the Bloch modes in the waveguide can be recovered. Six propagative Bloch modes are illustrated in Figure 5, with colored stress fields (positive stresses in blue, negative stresses in red). For compressional and flexural waves, the stress field is $\sigma_{zz}$, and for torsional waves, the stress field is $\tau_{z\theta}$, being $z$ the axial direction of the waveguide. The modes are illustrated for an assembly of 100 unit cells ($L = 50$ cm).

<figure><img src="image_7.png" /><figcaption>Figure 5: Propagative Bloch modes of the heat exchanger model.</figcaption></figure>

### 4.1.2. Eigenvalue problem linearizations

The three linearizations suitable for partial eigensolving are: Zhong-Williams, presented in Eq. (3) and denoted $\mu_{ZW}$; Huang *et al.*, presented in Eq. (8) and denoted $\mu_H$; and the $\eta$ linearization presented in Eq. (18). These linearizations are tested with three eigensolvers: `eigs()` in MATLAB; `eigs()` from the Arpack.jl package and `partialschur()`/`partialeigen()` from the ArnoldiMethod.jl package in Julia.

Dispersion results are presented in Figure 6 in the form of phase velocities $c_p = \omega/\Re(k)$ calculated in the [0, 10] kHz frequency range with 200 iterations. The dispersion curves are colored by the eigenvector error
Eigenvalue problem linearization

<figure><img src="image_8.png" /><figcaption>Figure 6: Dispersion curves for different eigensolvers and linearizations of the T-PQEP.</figcaption></figure>

$\epsilon$ presented in Eq. (21). Overall computation time of the partial eigensolver altogether with the eigenvector recovery part for each one of the nine combinations are presented in the upper-right corner in a hh:mm:ss format.

Some important notes can be drawn from Figure 6. Firstly, that the $\mu_{ZW}$ linearization is more than three times slower than $\mu_H$ regardless of the eigensolver, and this is caused by the SVD computation to find the weights of the linear combination with the $\mu$ related eigenvectors (see Ref. [21] for more details). Secondly, that $\mu_H$ perform well – with lower eigenvector errors $\epsilon$ than $\mu_{ZW}$ across the whole frequency range – whereas $\mu_{ZW}$ present high errors in the low frequency range for the bending mode, which is not an
axisymmetric mode and, therefore, necessitates of two similar eigenvectors but differing in spatial phase in the cross-section (see Ref. [34], where this multiplicity-related issue was studied in detail). Thirdly –now introducing the proposed linearization– that η is around two times faster than μ<sub>H</sub>, it does not present any issues related to cyclic symmetric waveguides like those of μ<sub>ZW</sub>, its eigenvector errors ϵ are lower than those of the μ<sub>ZW</sub> – at low frequencies and for the local fin modes – and just slightly higher than the μ<sub>H</sub> errors. Lastly, that the eigensolver from the ArnoldiMethod.jl package in Julia gives correct results only with the η linearization (both only *partialschur()* and combined *partialschur()*/*partialeigen()* were tested), the other dispersion curves being polluted by non-convergent solutions across the whole frequency range. This will be explored further with the second numerical model.

Since it has been shown that $\mu_H$ is more efficient than $\mu_{ZW}$ because the linear combinations do not necessitate an SVD, a comparison between the standard $\mu_H$ and the $\mu(\phi)$ linearization is now treated. A parametric analysis of the rotation $\phi$ influence on the $\mu(\phi)$ linearization performance is shown in Figure 7, where the time reduction factor is the relationship between $\mu_H$ and $\mu(\phi)$ eigensolving times (e.g., the time reduction factor obtained from Figures 6 (b) and (c) is 00:21:46 / 00:09:21 = 2.33) and the number of required solutions to the eigensolver (eigs() in MATLAB, taken as a standard eigensolver) is varied from 10 to 40. It should be noted that $\mu(0)$ represents the $\mu_H$ linearization, and $\mu(\pi/2)$ represents the $\eta$ linearization.

<figure><img src="image_9.png" /><figcaption>Figure 7: Rotation $\phi$ impact on the numerical performance of the $\mu(\phi)$ linearization. The selected eigensolver is eigs() from MATLAB. The highest reduction factors (better performance) correspond to $\phi = \pi/2$, which is the $\eta$ linearization.</figcaption></figure>

It can be seen that there is a distinguishable trend in the rotation influence, with the numerical performance increasing as the $\phi$ rotation approaches to the $\eta$ linearization ($\phi = \pi/2$), for which the highest reduction factors are found. This confirms the motivation that was depicted in Figures 2 and 3, which was to close the gap between the partial $\lambda$ spectrum of interest and the targeted minima of the $|\mu(\phi)|$ linearizations to improve convergence. However, it is important to note that this gap is efficiently closed only when the partial $\lambda$ spectrum is concentrated around the point 1 in the complex plane. For instance, with a periodic unit cell, this means that the propagative waves evolving in the frequency domain do not approach to the Brillouin edge at $\Re(k) = \pi/d$ (point -1 in the complex plane). If this occurs, as $\mu(\pi/2)$ do not close a partial spectrum gap, it would perform just slightly faster than $\mu(0)$, and thus its advantage could be almost not noticeable.

## 4.2. Homogeneous beam

The selected waveguide for this second example, although a simple homogeneous beam, aims to demonstrate the benefit of the proposed $\mu(\phi)$ linearization, where the $\phi$ parameter can be chosen to simply avoid the condition of repeated $\mu$ eigenvalues, and consequently improve the eigensolutions. Also, this example serves to validate a forced response computation with a numerical finite element model of an assembled waveguide. The numerical model of the unit cell, built with Gmsh [45], is depicted in Figure 8a.
<figure><img src="image_10.png" /><figcaption>Figure 8: Unit cell of the homogeneous beam. The unit cell is meshed with quadratic triangles and its length is $d = 3$ mm. The beam thickness is $t = 10$ mm. The number of degrees of freedom at the interfaces is $n = 42$. The beam is made of steel, with Young's modulus $E = 210.0$ GPa, Poisson's ratio $\nu = 0.3$, density $\rho = 7850$ kg/m$^3$ and loss factor $\eta = 0.002$.</figcaption></figure>

### 4.2.1. Dispersion characteristics

Dispersion curves are presented in Figure 8b in the form of phase velocities $c_p = \omega/\Re(k)$ calculated in the [0, 400] kHz frequency range. These curves are colored by the eigenvector error $\epsilon$ presented in Eq. (21). The dispersion results were generated after solving the T-PQEP (2) with the proposed $\mu(\phi)$ linearization of Eq. (17) by utilizing the eigensolver from the ArnoldiMethod.jl package in Julia, requiring 16 eigenpairs – 8 in each direction of propagation –. This particular eigensolver can address a GEVP with both complex non-hermitian matrices, but repeated eigenvalues may be undesirable (as described in their documentation) for the eigendecomposition in their implementation. To the left, it can be appreciated that the $\mu(0) = \mu_H$ linearization has remarkably high eigenvector errors $\epsilon$ and even noticeable errors in the eigenvalues – which provide $c_p$. This has been circumvented by slightly shifting the $\phi$ parameter away from zero. The choice of a $\delta\phi = 10^{-6}$ has been proven to be useful and generate correct dispersion curves with low eigenvector errors $\epsilon$.

With the $\mu(\phi)$ linearization, we have that for $\phi \neq n\pi$ with $n \in \mathbb{N}$ the sought $\mu$ eigenvalues are not repeated. This condition of a simple spectrum enhances the numerical computation and convergence of the iterative eigensolver, although performance is not improved with $\phi \approx 0$. However, the $\phi$ impact on performance was also tested for this small model ($n = 42$) and it was found that $\eta$ performs around 44% faster than $\mu_H$, which, despite not being of great advantage, reinforces the motivation of the $\mu(\phi)$ linearization of approaching the targeted $|\mu|$ values to the partial $\lambda$ spectrum of interest.

### 4.2.2. Forced response

The forced response of structural waveguides is one of the main applications of the WFEM [13, 21, 12], since the wave expansion can be used to describe dynamic motion of assemblies, both subjected to harmonic or transient loads [17], and to extremal or distributed loads [47].

Aiming to validate the proposed $\mu(\phi)$ linearization and present its advantages over the existent ones, the vibration response of an homogeneous beam is presented. The waveguide is formed by the assembly of 100 unit cells, i.e., an homogeneous steel beam with length $L = 300$ mm and thickness $t = 10$ mm. The beam is subjected to Dirichlet-Dirichlet boundary conditions: an uniform vertical harmonic displacement of amplitude $10^{-5}$ mm is applied on the left side, and a clamped condition is applied on the right side. The frequency response functions (FRF) of the displacement at the middle of the beam calculated with three different linearizations and with the same eigensolver (from ArnoldiMethod.jl) are shown in Figure 9.

It can be seen that the errors on the dispersion results from Figure 8b have an impact on the FRF calculated with the $\mu(0) = \mu_H$ linearization, whereas the $\mu(\delta\phi)$ linearization does not present any issue, being perfectly overlapped with the FRF provided by the FE model.

In summary, the $\phi$ parameter has been advantageous to circumvent the issues of the eigensolver included in the ArnoldiMethod.jl package, improving dispersion curves and frequency response functions by slightly
<figure><img src="image_11.png" /><figcaption>Figure 9: Frequency response function (FRF) of the cantilever beam.</figcaption></figure>

avoiding the repeated $\mu$ eigenvalues with a $\delta\phi \approx 0$. In fact, since both Zhong-Williams and Huang et al. linearizations provide repeated $\mu$ eigenvalues, this particular eigensolver would not have provided any useful results unless the $\mu(\phi)$ linearization had been used – which has also been the case for the heat exchanger model, in Figure 6 –.

## 5. Conclusions

The two existent linearizations for partial eigensolving of the WFEM eigenvalue problem, namely the Zhong-Williams and Huang et al. ones, provide repeated $\mu$ eigenvalues. It has been shown that the Zhong-Williams linearization may present numerical issues in some cyclic symmetry cases. Moreover, these linearizations are generalized eigenvalue problems with both complex non-hermitian matrices, and their partial eigensolving is not widely implemented (e.g., `scipy.sparse.linalg.eigs()` in Python and `geneig()` from the KrylovKit.jl package in Julia do not consider it). This inconvenient has drawn attention to the performance of some remaining capable eigensolvers such as `eigs()` in MATLAB, `eigs()` from Arpack.jl and `partialschur()`/`partialdecomp()` from ArnoldiMethod.jl in Julia, which were tested in this paper on the mentioned linearizations.

In this work, a novel $\mu(\phi)$ linearization was established. It avoids the condition of repeated $\mu$ eigenvalues while preserving the symplectic structure of the eigenvalue problem in the same manner as the existent linearizations. This has proved to be advantageous for convergence and computation time in the heat exchanger numerical example, a homogeneous waveguide in which the sought eigenvalues were located around the point 1 in the complex plane. Firstly, the Huang et al. linearization provides faster and better results than Zhong-Williams. Secondly, the proposed $\eta$ linearization performed more than two times faster than the Huang et al. linearization and does not suffer the numerical instabilities of the Zhong-Williams’ one. Moreover, for a simple homogeneous beam model, it has been shown that even slightly avoiding the repeated eigenvalues with $\mu(\delta\phi)$ is useful to improve the band diagram and forced response computation precision when using the ArnoldiMethod.jl package in Julia. Hence, this $\mu(\delta\phi)$ permits us to keep eigensolvers performing well in WFEM problems, despite having issues related to repeated eigenvalues for a complex non-hermitian GEVP.

In summary, the proposed $\mu(\phi)$ linearization includes the standard Huang et al. linearization (with $\phi = 0$), whereas taking $\phi = \pi/2$ gives the proposed $\eta$ linearization, which has been proved to be advantageous. This rotation parameter $\phi$ could lead to further explorations and be an added value to the WFEM practice. In that sense, the authors recommend the $\eta$ linearization as a faster, numerically stable to iterative eigensolvers and general purpose scheme for their dispersion analyses involving high-dimensional unit cell models.
CRediT authorship contribution statement

Writing – original draft: A.G.R.; Software: A.G.R.; Visualization: A.G.R.; Methodology: A.G.R., Q.Z., C.D.; Supervision: Q.Z., C.D.; Writing – review & editing: A.G.R., Q.Z., C.D.; Conceptualization: A.G.R., C.D.; Project administration: C.D.; Funding acquisition: C.D.

Acknowledgements

The authors would like to thank Vincent Mahe (Inria Rennes) for his fruitful suggestions and Olivier Robin (Université de Sherbrooke) for his support during the writing of the article.

References

[1] E. J. P. Miranda, S. F. Rodrigues, J. M. C. D. Santos, Complex dispersion diagram and evanescent modes in piezomagnetic phononic structures, Solid State Communications 346 (2022) 114697. doi:10.1016/j.ssc.2022.114697.

[2] M. I. Hussein, M. J. Leamy, M. Ruzzene, Dynamics of Phononic Materials and Structures: Historical Origins, Recent Progress, and Future Outlook, Applied Mechanics Reviews 66 (2014) 040802. doi:10.1115/1.4026911.

[3] D. J. Mead, Vibration Response and Wave Propagation in Periodic Structures, Journal of Engineering for Industry 93 (3) (1971) 783-792. doi:10.1115/1.3428014.

[4] A. Bertoncini, C. Liberale, 3D printed waveguides based on photonic crystal fiber designs for complex fiber-end photonic devices, Optica 7 (11) (2020) 1487-1494, publisher: Optica Publishing Group. doi:10.1364/OPTICA.397281.

[5] G. Okudan, H. Danawe, D. Ozevin, S. Tol, Torsional wave focusing in cylindrical structures with the conformal gradient-index phononic crystal lens, Journal of Applied Physics 129 (17) (2021) 174902. doi:10.1063/5.0050295.

[6] M. Moscatelli, R. Ardito, L. Driemeier, C. Comi, Band-gap structure in two- and three-dimensional cellular locally resonant materials, Journal of Sound and Vibration 454 (2019) 73-84. doi:10.1016/j.jsv.2019.04.027.

[7] V. F. D. Poggetto, A. L. Serpa, Elastic wave band gaps in a three-dimensional periodic metamaterial using the plane wave expansion method, International Journal of Mechanical Sciences 184 (2020) 105841. doi:10.1016/j.ijmecsci.2020.105841.

[8] M. I. Hussein, Reduced Bloch mode expansion for periodic media band structure calculations, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 465 (2109) (2009) 2825-2848. doi:10.1098/rspa.2008.0471.

[9] D. Krattiger, M. I. Hussein, Generalized Bloch mode synthesis for accelerated calculation of elastic band structures, Journal of Computational Physics 357 (2018) 183-205. doi:10.1016/j.jcp.2017.12.016.

[10] A. Palermo, A. Marzani, A reduced Bloch operator finite element method for fast calculation of elastic complex band structures, International Journal of Solids and Structures 191-192 (2020) 601-613. doi:10.1016/j.ijsolstr.2019.12.011.

[11] J.-M. Mencik, M. Ichchou, Multi-mode propagation and diffusion in structures through finite elements, European Journal of Mechanics - A/Solids 24 (5) (2005) 877-898. doi:10.1016/j.euromechsol.2005.05.004.

[12] J.-M. Mencik, On the low- and mid-frequency forced response of elastic structures using wave finite elements with one-dimensional propagation, Computers & Structures 88 (11) (2010) 674-689. doi:10.1016/j.compstruc.2010.02.006.

[13] D. Duhamel, B. Mace, M. Brennan, Finite element analysis of the vibrations of waveguides and periodic structures, Journal of Sound and Vibration 294 (1) (2006) 205-220. doi:10.1016/j.jsv.2005.11.014.

[14] M. Ichchou, J.-M. Mencik, W. Zhou, Wave finite elements for low and mid-frequency description of coupled structures with damage, Computer Methods in Applied Mechanics and Engineering 198 (15) (2009) 1311-1326. doi:10.1016/j.cma.2008.11.024.

[15] D. Claro, V. Denis, J.-M. Mencik, Defect localization in waveguide assemblies with curved joints via wave finite elements and time of flight analysis, European Journal of Mechanics - A/Solids 97 (2023) 104814. doi:10.1016/j.euromechsol.2022.104814.

[16] Y. Fan, M. Collet, M. Ichchou, L. Li, O. Bareille, Z. Dimitrijevic, Enhanced wave and finite element method for wave propagation and forced response prediction in periodic piezoelectric structures, Chinese Journal of Aeronautics 30 (1) (2017) 75-87. doi:10.1016/j.cja.2016.12.011.

[17] C. Droz, R. Boukadia, W. Desmet, A multi-scale model order reduction scheme for transient modelling of periodic structures, Journal of Sound and Vibration 510 (2021) 116312. doi:10.1016/j.jsv.2021.116312.

[18] F. Errico, M. Ichchou, F. Franco, S. De Rosa, O. Bareille, C. Droz, Schemes for the sound transmission of flat, curved and axisymmetric structures excited by aerodynamic and acoustic sources, Journal of Sound and Vibration 456 (2019) 221-238. doi:10.1016/j.jsv.2019.05.041.

[19] D. J. Mead, A general theory of harmonic wave propagation in linear periodic systems with multiple coupling, Journal of Sound and Vibration 27 (2) (1973) 235-260. doi:10.1016/0022-460X(73)90064-3.

[20] A. Hilliges, C. Mehl, V. Mehrmann, On the solution of palindromic eigenvalue problems, ECCOMAS 2004 - European Congress on Computational Methods in Applied Sciences and Engineering (01 2004).

[21] Y. Waki, B. Mace, M. Brennan, Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides, Journal of Sound and Vibration 327 (1) (2009) 92-108. doi:10.1016/j.jsv.2009.06.005.

[22] J.-M. Mencik, D. Duhamel, A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models, Finite Elements in Analysis and Design 101 (2015) 1-14. doi:10.1016/j.finel.2015.03.003.
[23] D. J. Mead, The forced vibration of one-dimensional multi-coupled periodic structures: An application to finite element analysis, Journal of Sound and Vibration 319 (1) (2009) 282-304. doi:10.1016/j.jsv.2008.05.026.

[24] E. K.-W. Chu, T.-M. Hwang, W.-W. Lin, C.-T. Wu, Vibration of fast trains, palindromic eigenvalue problems and structure-preserving doubling algorithms, Journal of Computational and Applied Mathematics 219 (1) (2008) 237-252. doi:10.1016/j.cam.2007.07.016.

[25] I. C. Ipsen, Accurate eigenvalues for fast trains, SIAM News 37 (9) (2004) 1-2.

[26] C.-H. Guo, W.-W. Lin, Solving a structured quadratic eigenvalue problem by a structure-preserving doubling algorithm, SIAM Journal on Matrix Analysis and Applications 31 (5) (2010) 2784-2801. doi:10.1137/090763196.

[27] W.-W. Lin, A new method for computing the closed-loop eigenvalues of a discrete-time algebraic Riccati equation, Linear Algebra and its Applications 96 (1987) 157-180. doi:10.1016/0024-3795(87)90342-9.

[28] R. Patel, On computing the eigenvalues of a symplectic pencil, in: [1992 Proceedings of the 31st IEEE Conference on Decision and Control, 1992, pp. 1921-1926 vol.2. doi:10.1109/CDC.1992.371096.

[29] W. Zhong, F. Williams, On the direct solution of wave propagation for repetitive structures, Journal of Sound and Vibration 181 (3) (1995) 485-501. doi:10.1006/jsvi.1995.0153.

[30] T.-M. Huang, W.-W. Lin, J. Qian, Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration of fast trains, SIAM Journal on Matrix Analysis and Applications 30 (4) (2009) 1566-1592. doi:10.1137/080713550.

[31] T. Gras, Couplage de méthodes d'éléments finis standards (fem) et ondulatoires (wfem) pour le calcul de la réponse vibratoire d'une voie ferrée, Ph.D. thesis, Université de technologie de Compiègne (2017).

[32] A. Cicirello, B. R. Mace, M. J. Kingan, Y. Yang, Sensitivity analysis of generalised eigenproblems and application to wave and finite element models, Journal of Sound and Vibration 478 (2020) 115345. doi:10.1016/j.jsv.2020.115345.

[33] T.-M. Huang, W.-W. Lin, H. Tian, G.-H. Chen, Computing the full spectrum of large sparse palindromic quadratic eigenvalue problems arising from surface green's function calculations, Journal of Computational Physics 356 (2018) 340-355. doi:10.1016/j.jcp.2017.12.011.

[34] W. Wang, Y. Fan, L. Li, Extending Zhong-Williams scheme to solve repeated-root wave modes, Journal of Sound and Vibration 519 (2022) 116584. doi:10.1016/j.jsv.2021.116584.

[35] V. Mehrmann, H. Voss, Nonlinear eigenvalue problems: a challenge for modern eigenvalue methods, GAMM-Mitteilungen 27 (2) (2004) 121-152. doi:10.1002/gamm.201490007.

[36] M. A. Gosson, Symplectic Methods in Harmonic Analysis and in Mathematical Physics, Birkhäuser Basel, 2011. doi:10.1007/978-3-7643-9992-4.

[37] K. A. Foss, Co-Ordinates Which Uncouple the Equations of Motion of Damped Linear Dynamic Systems, Journal of Applied Mechanics 25 (3) (1958) 361-364. doi:10.1115/1.4011828.

[38] F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev. 43 (2001) 235-286.

[39] C. Droz, J.-P. Lainé, M. Ichchou, G. Inquiété, A reduced formulation for the free-wave propagation analysis in composite structures, Composite Structures 113 (2014) 134-144. doi:10.1016/j.compstruct.2014.03.017.

[40] D. J. Mead, Wave propagation in continuous periodic structures: Research contributions from Southampton, 1964-1995, Journal of Sound and Vibration 190 (3) (1996) 495-524. doi:10.1006/jsvi.1996.0076.

[41] R. B. Lehoucq, D. C. Sorensen, C. Yang, ARPACK Users' Guide, Society for Industrial and Applied Mathematics, 1998. doi:10.1137/1.9780898719628.

[42] R. Boukadia, Vibroacoustics of periodic structures: model order reduction, characterization, optimization, Ph.D. thesis, Université de Lyon (2021).

[43] O. M. Malinowski, M. S. Lindsey, J. K. V. Velsor, Ultrasonic guided wave testing of finned tubing, in: Proceedings of the ASME2015 Pressure Vessels and Piping Conference, 2015, p. V005T10A002.

[44] B. Zheng, J.-W. Su, Y. Xie, J. Miles, H. Wang, W. Gao, M. Xin, J. Lin, An autonomous robot for shell and tube heat exchanger inspection, Journal of Field Robotics 39 (8) (2022) 1165-1177. doi:10.1002/rob.22102.

[45] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309-1331. doi:10.1002/nme.2579.

[46] C. Droz, R. Boukadia, M. Ichchou, W. Desmet, Diffusion-based design of locally resonant sub-systems using a reduced wave finite element framework, in: Proceedings of ISMA2018 and USD2018, 2018, pp. 3071-3083.

[47] T. Hoang, D. Duhamel, G. Foret, Wave finite element method for waveguides and periodic structures subjected to arbitrary loads, Finite Elements in Analysis and Design 179 (2020) 103437. doi:10.1016/j.finel.2020.103437.