# Reduced Frequency-driven Bloch Wave Decomposition for Harmonic Analysis of Finite Periodic Structures

Christophe Droz¹, Alvaro Gavilan-Rojas¹

¹ I4S team, Inria - National Institute for Research in Digital Science and Technology, Rennes, France, christophe.droz@inria.fr, alvaro.gavilan-rojas@inria.fr

**Abstract:** A frequency-driven, reduced wave finite element method is proposed for periodic structures which are solely defined by their dynamic stiffness matrices. Using an S-S⁻¹ transform combined with an adaptive eigenvector sampling and a ω-driven Bloch wave decomposition, we achieve dispersion and harmonic analyses with enhanced accuracy and efficiency.

Since phononic structures become more intricate both in terms of physics and geometry, their numerical modeling demands increasingly high-dimensional unit-cell finite element (FE) descriptions. Moreover, certain models incorporate frequency-dependent effects that deviate from standard linear frequency polynomials. Specifically, the Dynamic Stiffness Matrix (DSM, noted D(ω)) in such cases no longer takes the form: D(ω) = K + jωC - ω²M, where M, C and K stand for the mass, damping, and stiffness FE matrices, respectively. This arises in various scenarios, including viscoelastic materials, fluid-structure interactions, nonlinear damping mechanisms, or when effective masses are derived from the homogenization of locally resonant systems. Additionally, such scenarios can be encountered when D(ω) is derived from experimental or numerical data, linearized formulations, meshless methods, or other black-box modeling approaches that preclude its decomposition into ω-independent system matrices. In such cases, Model-Order Reduction (MOR) methods leveraging modal analyses (i.e., ω(k)) to build a reduced Bloch mode description of the unit-cell (e.g., Hussein et al.¹) may become unavailable, since one has to solely rely on these large, sparse, D(ω) input matrices for band diagrams calculations, vibration and vibro-acoustic simulations. Several works have been targeting this problem, most often aimed at the computation of complex dispersion relations (e.g., EPWE², Floquet-Bloch decomposition³, or reduced D(ω) techniques⁴). In the framework of the Wave Finite Element Method (WFEM) in 1d-periodic structures, denoting λ=e⁻⁴kd the propagation constant, the dispersion relation writes:

$$[D_{RL} \lambda + (D_{RR}+D_{LL}) + D_{LR} \lambda^{-1}] \Psi(\omega) = 0 \quad (1)$$

This eigenvalue problem has in fact a T-palindromic quadratic structure (T-QEP), and provides complex-valued (Λ(ω), Ψ(ω)). These solutions can later be used to achieve a ω-driven Bloch Wave Decomposition (denoted ω-BWD) of the structure's displacements u_n, which takes the following form in a finite N-unit-cells assembly:

$$u_n(\omega) = [\Psi^+] [\Lambda]^n a^+ + [\Psi^-] [\Lambda]^{N-n} a^- \quad (2)$$

Although this method has the potential to enable fast complex dispersion and dynamic analyses in finite or infinite phononic structures subjected to harmonic excitations, Eq. (1) still raises multiple numerical challenges which currently jeopardize the method's accuracy, versatility and computational scalability.

We address these issues by proposing a modelling framework for ω-driven periodic structure models, incorporating MOR, well-conditioned T-QEP resolution schemes and Bloch wave decomposition of the finite or infinite waveguide subjected to Dirichlet or Neumann-type boundary conditions.

## Purely ω-driven MOR for the WFEM

Model-order reduction usually relies on a sampling of the Brillouin Zone in terms of given wave-numbers, where the propagating modes are collected up to a pre-determined frequency threshold. These strategies are akin to a proper orthogonal decomposition, where the linear or quadratic EVP relies on the separation between M, C and K matrices. The major limitation lies in the difficulty to determine the ωmax threshold to ensure that evanescent waves can be described within the bandwidth of interest. Instead, in the WFEM, frequency snapshots are used, which consists in the resolution of the ω-driven T-QEP on a few samples, filtered up to a desired evanescence |λ|min threshold to build a reduced unit-cell. Selected solutions are shown in the dispersion diagram Fig. 1.b).
## An alternative S-S⁻¹ transform for solving the ω-driven T-QEP

To compute these snapshots, Huang et al.⁵ introduced the S+S⁻¹ transform, a structure-preserving algorithm that enables the use of Arnoldi method to extract the few dominant Bloch solutions from the T-QEP. However, this approach artificially doubles the geometric multiplicity of the system, leading to longer resolution times and degraded eigenvector solution accuracy. An alternative approach, introduced by Gavilan-Rojas et al.⁶, involves rotating the T-QEP by the angle π while preserving its symplectic structure, thereby introducing an S-S⁻¹ transform. This modification halves the geometric multiplicity, ensuring solution accuracy through partial eigensolvers (i.e., *eigs*) and reducing resolution time per snapshot, making it even shorter than in the ω(k) scheme.

<figure><img src="image_1.png" /><figcaption>**Figure 1** Finite Element description (46k DOF per unit-cell, meshed with quadratic tetrahedra) (a). Complex dispersion diagram obtained through ω-driven MOR. Solutions retained for MOR on 11 ω-snapshots are marked with x and spatial decay on colormap (b). The mobility (i.e., FRF) is computed by ω-BWD for a finite 5-unit-cells structure under harmonic force excitation (c).</figcaption></figure>

## Application of the reduced ω-BWD for steady-state analysis of a finite chiral metamaterial

Once the reduced *wave* finite elements are derived, the amplitudes $a^+$, $a^-$ and dynamic response can be computed under various harmonic boundary conditions (finite or infinite) using the WFEM framework. We apply the proposed MOR strategy to the chiral, syndiotactic-type phononic crystal investigated by Bermagini et al.⁷ whose unit-cell model is depicted in Fig. 1.a). This syndiotactic model exhibits a full bandgap between 180 and 500 Hz. The reduced-order model takes 10 seconds to build, including the S-S⁻¹ transformed ω-snapshots and the SVD of the selected eigenvectors. Computing the complex dispersion diagrams (Fig. 1.b) and the harmonic response of the finite structure (Fig. 1.c) over 5000 frequencies can then be achieved in less than 1 second. Noteworthy, this ω-sampling strategy – which can be seen as a binary search for eigenvector matching – is in accordance with the premise that fewer snapshots are needed within bandgaps, since Bloch solutions remain evanescent, but cut-on frequencies or complex wave phenomena (e.g., veering, coupling or locking) need an enriched local sampling.

## Conclusion

Dynamic analyses of high-dimensional finite periodic structures can therefore remain computationally efficient without having to presuppose the frequency dependence of the dynamic stiffness matrix. These MOR strategies hold significant potential for multi-indicator simulation of intricate phononic crystals.

## References

¹ M. I. Hussein, *P. Roy. Soc. A.*, **465**(2109), 2825-2848 (2009).

² Y.C. Hsue, A.J. Freeman and B.Y. Gu, *Phys. Rev. B.*, **72**(19), 195118 (2005).

³ M. Collet, M. Ouisse, M. Ruzzene and M. N. Ichchou, *Int. J. Sol. Struct.*, **48**(20), 2837-2848 (2011).

⁴ C. Droz, J.P. Lainé, M.N. Ichchou, and G. Inquiété, *Compos. Struct.*, **113**, 134-144 (2014).

⁵ T.-M. Huang, W.-W. Lin, J. Qian, *SIAM J. Matrix Anal. A.*, **30**(4), 1566-1592 (2009).

⁶ A. Gavilan-Rojas, Q. Zhang, and C. Droz, *J. Sound Vib.*, **593**, 118652 (2024).

⁷ A. Bergamini, M. Miniaci, T. Delpero, D. Tallarico, B. Van Damme et int. and A. Zemp, *Nat. Commun.*, **10**(1), 4525 (2019).