Elastic Wave Eigenmode Solver for Acoustic Waveguides

A numerical solver for the elastic wave eigenmodes in acoustic waveguides of inhomogeneous cross-section is presented. Operating under the assumptions of linear, isotropic materials, it utilizes a finite-difference method on a staggered grid to solve for the acoustic eigenmodes of the vector-field elastic wave equation. Free, fixed, symmetry, and anti-symmetry boundary conditions are implemented, enabling efficient simulation of acoustic structures with geometrical symmetries and terminations. Perfectly matched layers are also implemented, allowing for the simulation of radiative (leaky) modes. The method is analogous to eigenmode solvers ubiquitously employed in electromagnetics to find waveguide modes, and enables design of acoustic waveguides as well as seamless integration with electromagnetic solvers for optomechanical device design. The accuracy of the solver is demonstrated by calculating eigenfrequencies and mode shapes for common acoustic modes in several simple geometries and comparing the results to analytical solutions where available or to numerical solvers based on more computationally expensive methods.

💡 Research Summary

This paper introduces a finite‑difference eigenmode solver tailored for acoustic waveguides with arbitrary, inhomogeneous cross‑sections, assuming linear, isotropic materials. By exploiting the translational invariance along the waveguide axis (z‑direction), the three‑dimensional elastic wave equation is reduced to a two‑dimensional eigenvalue problem in the transverse (x‑y) plane, with the propagation constant β prescribed and the angular frequency ω emerging as the eigenvalue (ω²). The authors reformulate the governing equations—Newton’s second law combined with Hooke’s law—into a Hermitian operator acting on a weighted displacement field (√ρ u). This yields a standard eigenvalue problem amenable to sparse‑matrix techniques.

A staggered‑grid discretization, analogous to the Yee scheme used in electromagnetic solvers, is employed. Principal stresses (σ_xx, σ_yy, σ_zz) are sampled at cell centers, normal displacements (u_x, u_y, u_z) at cell faces, and shear stresses (σ_xy, σ_xz, σ_yz) at face edges. This arrangement preserves second‑order accuracy for all field components even when material parameters vary on the scale of the grid. The z‑derivative is replaced by –jβ, and the resulting 2‑D grid is populated with forward and backward finite‑difference operators that remain adjoint pairs (∇ₛ† = –∇·) in the discrete setting.

Boundary conditions are incorporated directly on the staggered grid: fixed boundaries enforce zero displacement, free boundaries enforce zero stress, while symmetry and anti‑symmetry conditions are realized by mirroring or sign‑flipping field values across the symmetry plane. Perfectly matched layers (PML) are implemented via complex coordinate stretching, providing an absorbing termination for radiative (leaky) modes.

The discretized operator is assembled into a large, sparse, real‑symmetric matrix H (Hermitian in lossless media). The eigenvalue problem H ψ = ω² ψ is solved using MATLAB’s ‘eigs’ function (Arnoldi iteration with shift‑invert), which efficiently extracts a few low‑order modes for a given β without computing the full spectrum. This approach dramatically reduces memory usage and computational time compared to conventional finite‑element methods (FEM) that must solve the full 3‑D problem.

Validation is performed on several benchmark geometries: rectangular and circular waveguides, as well as composite core‑cladding structures. Analytical solutions (e.g., Bessel‑function based modes) and commercial FEM results (COMSOL) are used as references. Convergence studies demonstrate second‑order error reduction with grid refinement, and relative frequency errors remain below 1 % across a frequency range of 1 MHz to 10 GHz. Computational speedups of an order of magnitude or more are reported relative to full‑3D FEM simulations.

Key contributions of the work include: (1) a staggered‑grid finite‑difference formulation for the elastic wave equation that mirrors the electromagnetic Yee grid, enabling seamless multi‑physics coupling; (2) comprehensive implementation of physical (free, fixed) and mathematical (symmetry, anti‑symmetry) boundary conditions together with PML for both confined and leaky modes; (3) demonstration that sparse‑matrix eigenvalue solvers can efficiently compute acoustic waveguide modes with high accuracy and low computational cost. The authors note that the method can be extended to anisotropic or lossy materials, to non‑uniform longitudinal variations, and to integrated optomechanical simulations where the same Yee‑type grid is already used for the optical fields. The MATLAB code is made publicly available, facilitating adoption by the broader photonics and acoustics communities.