Computing guided waves

We are interested in the dispersion relation of an homogeneous elastic waveguide, infinite in the direction of propagation, with a two-dimensional cross-section. We solve the three-dimensional elastodynamics equation,

\[ \boldsymbol{\nabla} \cdot (\boldsymbol{C} : \boldsymbol{\nabla} \boldsymbol{u}) = \rho \ddot{\boldsymbol{u}} \tag{1}\]

where $\boldsymbol{C}$ is the elasticity tensor, $\rho$ is the density, and $\boldsymbol{u}$ is the displacement field, subjected to the appropriate boundary conditions at the edges of the waveguide. After inserting a propagative displacement ansatz, the dispersion relation is given by the non-trivial solutions of this boundary value problem.

We detail below how to obtain these non-zero solutions using the spectral collocation method that relies on discretizing the cross-section of the waveguide, transforming the boundary value problem into a polynomial eigenvalue problem. This method has already been presented for plates in Kiefer, 2022 and for strips with a rectangular cross-section in Delory et al., 2024. We focus here on its implementation in cylindrical coordinates for tubes and curved strips, as well as on the implementation of fluid coupling.

Tube

We consider a tube with mean radius $R$ and thickness $h$. We work in cylindrical coordinates $(x,r,\theta)$ and the propagative displacement ansatz reads,

\[ \boldsymbol{u} = \boldsymbol{u}(k_x,r,m,\omega)e^{i(k_x x + m\theta - \omega t)}, \tag{2}\]

evidencing that a one-dimensional discretization, in the radial direction, allows us to compute axial guided waves in a tube.

We recall that in cylindrical coordinates $\boldsymbol{\nabla} = \boldsymbol{e}_x \partial_x + \boldsymbol{e}_r \partial_r + \boldsymbol{e}_\theta \frac{1}{r} \partial_\theta$, so that,

\[ \boldsymbol{\nabla u} = \boldsymbol{e}_x \partial_x \otimes \boldsymbol{u} + \boldsymbol{e}_r \partial_r \otimes \boldsymbol{u} + \boldsymbol{e}_\theta \frac{1}{r} \otimes \boldsymbol{A} \cdot \boldsymbol{u}, \tag{3}\]

where $\boldsymbol{A} = \boldsymbol{I}\partial_\theta + \boldsymbol{e}_\theta \otimes \boldsymbol{e}_r - \boldsymbol{e}_r \otimes \boldsymbol{e}_\theta$ with $\boldsymbol{I}$ the three-dimensional identity tensor.

Inserting $(2)$ in $(1)$, we get,

\[\begin{split} \left[ (ik_x)^2 \boldsymbol{c}_{xx} + ik_x \left( (\boldsymbol{c}_{xr} + \boldsymbol{c}_{rx}) \partial_r + \frac{1}{r} \left(\boldsymbol{c}_{rx} + \boldsymbol{c}_{x\theta} \cdot \boldsymbol{A} + \boldsymbol{A}\cdot\boldsymbol{c}_{\theta x}\right)\right)\right. + \\ \left. \boldsymbol{c}_{rr} \partial^2_r + \frac{1}{r} \left(\boldsymbol{c}_{rr} + \boldsymbol{c}_{r\theta} \cdot \boldsymbol{A} + \boldsymbol{A} \cdot \boldsymbol{c}_{\theta r} \right) \partial_r + \frac{1}{r^2} \boldsymbol{A} \cdot \boldsymbol{c}_{\theta\theta} \cdot \boldsymbol{A} + \rho \omega^2 \boldsymbol{I}\right] \cdot \boldsymbol{u} &= 0, \end{split} \tag{4}\]

where $\boldsymbol{c}_{ij} = \boldsymbol{e}_i \cdot \boldsymbol{C} \cdot \boldsymbol{e}_j$ and $\boldsymbol{A}$ now reads $\boldsymbol{A} = i m \boldsymbol{I} + \boldsymbol{e}_\theta \otimes \boldsymbol{e}_r - \boldsymbol{e}_r \otimes \boldsymbol{e}_\theta$.

The elastodynamics equation is complemented by boundary conditions at $r = R - h/2$ and $r = R + h/2$. Here, we describe the case of a tube in vacuum where we impose homogeneous Neumann boundary conditions,

\[ \left.\boldsymbol{e}_r \cdot \boldsymbol{C} : \boldsymbol{\nabla u}\right|_{r = R \pm h/2} = \left[ik_x \boldsymbol{c}_{rx} + \boldsymbol{c}_{rr} \partial_r + \frac{1}{r} \boldsymbol{c}_{r\theta} \cdot \boldsymbol{A}\right] \cdot \boldsymbol{u} = 0. \tag{5}\]

The boundary value problem defined by equations (4) and (5) is now discretized using the spectral collocation method. The displacement vector is approximated at the $N$ Chebyshev-Gauss-Lobatto collocation points $r_i, i \in \{1\ldots N\}$ located in the interval $[R-h/2, R+h/2]$, becoming a vector of size $3N$, and the first and second order differentiation matrices along the radial direction, $D_r$ and $D_{rr}$ of size $N \times N$, are computed using chebdif.jl.

The discrete form of equations (4) and (5) is,

\[\begin{aligned} \left[(ik_x)^2L_{kk} + ik_x L_k + L_0 + \omega^2 M\right]u &= 0, \quad \mathrm{for} \; r \in [R-h/2,R+h/2]\\ \left[ik_x B_k + B_0 \right]u &= 0, \quad \mathrm{for} \; r = R \pm h/2 \end{aligned} \tag{6}\]

with,

\[\begin{aligned} &L_{kk} = \boldsymbol{c}_{xx} \otimes I_N, \quad L_k = (\boldsymbol{c}_{xr} + \boldsymbol{c}_{rx}) \otimes D_r + \left(\boldsymbol{c}_{rx} + \boldsymbol{c}_{x\theta} \cdot \boldsymbol{A} + \boldsymbol{A}\cdot\boldsymbol{c}_{\theta x}\right) \otimes R_N^{-1}, \\ &L_0 = \boldsymbol{c}_{rr} \otimes D_{rr} + \left(\boldsymbol{c}_{rr} + \boldsymbol{c}_{r\theta} \cdot \boldsymbol{A} + \boldsymbol{A} \cdot \boldsymbol{c}_{\theta r} \right) \otimes R_N^{-1}D_r + \boldsymbol{A} \cdot \boldsymbol{c}_{\theta\theta} \cdot \boldsymbol{A} \otimes R_N^{-2}, \quad M = \rho \boldsymbol{I} \otimes I_N, \\ &B_k = \boldsymbol{c}_{rx} \otimes I_N, \quad B_0 = \boldsymbol{c}_{rr} \otimes D_r + \boldsymbol{c}_{r\theta} \cdot \boldsymbol{A} \otimes R_N^{-1}, \end{aligned}\]

where $\otimes$ denotes the Kronecker product between two matrices, $I_N$ is the identity matrix of size $N \times N$, $R_N^{-1}$ is the matrix $\mathrm{Diagonal}(1/r_i, i =1 \ldots N)$, and $u$ is the discretized displacement vector.

Finally, we obtain a polynomial eigenvalue problem of size $3N \times 3N$ by replacing the rows of the elastodynamics equation in (6) corresponding to boundary collocation points (i.e. the rows $1$, $N$, $N+1$, $2N$, $2N+1$, and $3N$) by the corresponding rows of the boundary condition. We solve this problem for the eigenpair $(k_x,u)$ parametrized in $\omega$ using the polyeig function.

Tube with inner fluid

We consider that the elastic tube is filled with fluid. As the domain $r \in [0, R - h/2]$ is finite, we choose to discretize it.

The only difference lies in the choice of the collocation points. Indeed, the Chebyshev-Gauss-Lobatto points include the domain boundaries, leading to a singularity at $r = 0$ which makes them unsuitable to discretize the fluid domain. Instead, we choose to use $P$ Legendre-Gauss-Radau collocation points which lie in the $]0, R - h/2]$ interval. These collocation points are computed using the function lgrpointsleft and the differentiation matrices $Df_r$ and $Df_{rr}$ (with size $P \times P$) using the function lgrdiffleft.

We couple the two domains at $r = R - h/2$ by expressing the continuity of the normal traction and of the normal displacement,

\[\begin{aligned} \boldsymbol{e}_r \cdot \boldsymbol{C}_s : \boldsymbol{\nabla} \boldsymbol{u}_s - \boldsymbol{e}_r \cdot \boldsymbol{C}_f : \boldsymbol{\nabla} \boldsymbol{u}_f = 0, \\ \boldsymbol{e}_r \cdot \boldsymbol{u}_s - \boldsymbol{e}_r \cdot \boldsymbol{u}_f = 0, \end{aligned} \tag{7}\]

where the subscripts $s$ and $f$ stand for solid and fluid.

Practically, this amounts to constructing two polynomial eigenvalue problems, one for the solid with size $3N \times 3N$ and one for the fluid with size $3P \times 3P$, as described in the section above. The matrices $L^s_{kk}$, $L^s_{k}$, $L^s_{0}$, $M^s$, and $L^f_{kk}$, $L^f_{k}$, $L^f_{0}$, $M^f$ (for the solid and fluid problem, respectively) are then assembled into four matrices of size $3(N+P)$ so that the elastodynamics equation reads,

\[\begin{aligned} \left[(ik_x)^2\begin{pmatrix}L^s_{kk} & 0_{3N\times 3P} \\ 0_{3P\times 3N} & L^f_{kk} \end{pmatrix} + (ik_x)\begin{pmatrix}L^s_{k} & 0 \\ 0 & L^f_{k}\end{pmatrix} + \begin{pmatrix}L^s_{0} & 0 \\ 0 & L^f_{0}\end{pmatrix} + \omega^2\begin{pmatrix}M^s & 0 \\ 0 & M^f\end{pmatrix}\right] \begin{pmatrix} u_s \\ u_f \end{pmatrix} = 0. \end{aligned} \tag{8}\]

The Neumann boundary condition at $r = R + h/2$ is then incorporated by replacing the rows $1$, $N+1$, and $2N+1$ of equation (8) by the corresponding rows of of the boundary condition in (6) taking care of adding $3P$ trailing zeros to accommodate the size increase of the coupled displacement vector.

The coupling boundary conditions given by the equations (7) are discretized similarly as above and injected in the rows $N$, $2N$, $3N$, $3N+1$, $3N+P+1$, and $3N+2P+1$ of the assembled matrices.

Finally, we obtain a polynomial eigenvalue problem of size $3(N+P)$ which can be solved for the eigenpair $(k_x,u)$ parametrized by $\omega$.

Tube with outer fluid

We consider an elastic tube surrounded by an infinite fluid medium. This time it is not convenient to discretize the fluid domain, and we choose to take into account its presence in the boundary condition at $r = R + h/2$ which now reads:

\[ \boldsymbol{e}_r\cdot(\boldsymbol{C} : \boldsymbol{\nabla} \boldsymbol{u})|_{r=R+h/2} = \frac{\rho_f\omega^2}{\gamma}\frac{H_m\left(\gamma r\right)}{H'_m\left(\gamma r\right)} u_r \boldsymbol{e}_r, \tag{9}\]

where $\rho_f$ is the fluid density, $H_m$ is the Hankel function of the first kind of order $m$, and $\gamma = \sqrt{k_f^2 - k_x^2}$ with $k_f$ the wavenumber in the fluid. Directly discretizing equation (9) and incorporating it in equation (6) transforms the polynomial eigenvalue problem into a non-linear eigenvalue problem.

As solving such non-linear eigenvalue problems is challenging, we turn to the iterative approach proposed by Gravenkamp et al., 2014. We first use the following approximate boundary condition,

\[ \boldsymbol{e}_r\cdot(\boldsymbol{C} : \boldsymbol{\nabla} \boldsymbol{u})|_{r=R+h/2} = -\frac{\rho_f\omega^2}{k_x} u_r \boldsymbol{e}_r, \tag{10}\]

which can be discretized as,

\[ \left[(ik_x)^2 \boldsymbol{c}_{rx} \otimes I_N + ik_x\left(\boldsymbol{c}_{rr} \otimes D_r + \boldsymbol{c}_{r\theta} \cdot \boldsymbol{A} \otimes R_N^{-1} \right)+ \rho_f \omega^2 \boldsymbol{I}e_r^T \otimes I_N \right] u = 0, \tag{11}\]

where $e_r$ is the vector $(0,1,0)$, and included in the relevant rows of equation (6).

We can solve the eigenvalue problem,

\[\begin{split} \left[(ik_x)^2L_{kk} + ik_x L_k + L_0 + \omega^2 M\right]u &= 0, \quad \mathrm{for} \; r \in [R-h/2,R+h/2]\\ \left[(ik_x)^2 B_k + ikB_0 + \omega^2 B_M\right]u &= 0, \quad \mathrm{for} \; r = R + h/2 \\ \left[ik_x B_k + B_0 \right]u &= 0, \quad \mathrm{for} \; r = R - h/2 \end{split} \tag{12}\]

with $B_M = i\rho_f \boldsymbol{I}e_r^T \otimes I_N$, and obtain approximate solutions that we denote $(k_x^+, u^+)$.

We now consider the exact eigenvalue problem,

\[ \left[(ik_x)^2L_{kk} + ik_x L_k + L_0 + \omega^2 M(k_x)\right]u = 0, \tag{13}\]

where the appropriate rows of the matrices have been replaced to accommodate the Neumann boundary condition at $r=R-h/2$ and the exact boundary condition given by equation (9) at $r=R+h/2$. We highlight the fact that the eigenvalue problem becomes non-linear by making the dependence of $M$ on $k_x$ explicit.

We perform a Taylor expansion of this eigenvalue problem in the vicinity of $k_x^+$,

\[\begin{aligned} (ik_x)^2L_{kk} + ik_x L_k + L_0 + \omega^2 M(k_x) &= ik_x\left(2ik_x^+L_{kk}+L_k\right) - \left((ik_x^+)^2L_{kk} - L_0 - \omega^2M(k_x^+)\right) \\ & = ik_x A_2(k_x^+) - A_1(k_x^+), \end{aligned} \tag{14}\]

which we inject in equation (10) to obtain the following generalized eigenvalue problem,

\[ ik_xA_2(k_x^+) u = A_1(k_x^+) u. \tag{15}\]

Inverse iteration is then used to improve the approximate eigenpairs $(k_x^+,u^+)$, starting from their initial value $k_x^{+(0)}$ and $u^{+(0)}$ by computing $u^{+(1)}$

\[ \left(A_1(k_x^{+(0)}) - ik_x^{+(0)}A_2(k_x^{+(0)}) \right)u^{+(1)} = A_2 u^{+(0)}, \tag{16}\]

and $k_x^{+(1)}$,

\[ ik_x^{+(1)}A_2(k_x^{+(0)})u^{+(1)} = A_1(k_x^{+(0)})u^{+(1)}, \tag{17}\]

using the backslash operator to invert the matrices in equations (16) and (17). This procedure is iterated until convergence is reached.

Curved strip

We now focus on describing the propagation of axial waves in curved strip with radius $R$, thickness $h$, and curvilinear width $w$. We work in cylindrical coordinates $(x,r,\theta)$ and the propagative displacement ansatz reads,

\[ \boldsymbol{u} = \boldsymbol{u}(k_x,r,\theta,\omega)e^{i(k_x x - \omega t)}, \tag{18}\]

indicating that we perform a two-dimensional discretization, unlike in the case of a tube.

Inserting this ansatz in equation (1) still yields equation (4),

\[ \begin{split} \left[ (ik_x)^2 \boldsymbol{c}_{xx} + ik_x \left( (\boldsymbol{c}_{xr} + \boldsymbol{c}_{rx}) \partial_r + \frac{1}{r} \left(\boldsymbol{c}_{rx} + \boldsymbol{c}_{x\theta} \cdot \boldsymbol{A} + \boldsymbol{A}\cdot\boldsymbol{c}_{\theta x}\right)\right)\right. + \\ \left. \boldsymbol{c}_{rr} \partial^2_r + \frac{1}{r} \left(\boldsymbol{c}_{rr} + \boldsymbol{c}_{r\theta} \cdot \boldsymbol{A} + \boldsymbol{A} \cdot \boldsymbol{c}_{\theta r} \right) \partial_r + \frac{1}{r^2} \boldsymbol{A} \cdot \boldsymbol{c}_{\theta\theta} \cdot \boldsymbol{A} + \rho \omega^2 \boldsymbol{I}\right] \cdot \boldsymbol{u} &= 0, \end{split} \]

but this time the matrix $\boldsymbol{A}$ still involves an azimuthal derivative as $\boldsymbol{A} = \boldsymbol{I}\partial_\theta + \boldsymbol{e}_\theta \otimes \boldsymbol{e}_r - \boldsymbol{e}_r \otimes \boldsymbol{e}_\theta$.

Boundary conditions are required at $r = R \pm h/2$ and at $\theta = \pm w/(2R)$. We show here the case of a strip clamped on its lateral edges in vacuum,

\[\begin{aligned} \left.\boldsymbol{e}_r \cdot \boldsymbol{C} : \boldsymbol{\nabla u}\right|_{r = R \pm h/2} = \left[ik_x \boldsymbol{c}_{rx} + \boldsymbol{c}_{rr} \partial_r + \frac{1}{r} \boldsymbol{c}_{r\theta} \cdot \boldsymbol{A}\right] \cdot \boldsymbol{u} &= 0, \\ \left.\boldsymbol{u}\right|_{\theta = \pm w/(2R)} = \boldsymbol{I} \cdot \boldsymbol{u} &= 0. \end{aligned} \tag{19}\]

We discretize this boundary value problem on a grid of $N \times P$ Chebyshev-Gauss-Lobatto collocation points and calculate the radial $N \times N$ differentiation matrices $D_r$ and $D_{rr}$ as well as the azimuthal $P \times P$ differentiation matrices $D_\theta$ and $D_{\theta\theta}$ using chebdif. The discrete form of the boundary value problem given by equations (4) and (19) reads,

\[\begin{split} \left[(ik_x)^2L_{kk} + ik_x L_k + L_0 + \omega^2 M\right]u &= 0, \quad \mathrm{for} \; r \in [R-h/2,R+h/2], \;\theta \in [-w/(2R), w/(2R)]\\ \left[(ik_x) B_k + B_0 \right]u &= 0, \quad \mathrm{for} \; r = R \pm h/2 \\ \left[\boldsymbol{I} \otimes I_{NP} \right]u &= 0, \quad \mathrm{for} \; \theta = \pm w/(2R), \end{split} \tag{20}\]

with,

\[\begin{aligned} &L_{kk} = \boldsymbol{c}_{xx} \otimes I_{NP}, \quad L_k = (\boldsymbol{c}_{xr} + \boldsymbol{c}_{rx}) \otimes I_P \otimes D_r + \left(\boldsymbol{c}_{rx} \otimes I_P + (\boldsymbol{c}_{x\theta} + \boldsymbol{c}_{\theta x}) \otimes D_\theta + (\boldsymbol{c}_{x\theta}A + A\boldsymbol{c}_{\theta x}) \otimes I_P \right) \otimes R_N^{-1}, \\ &L_0 = \boldsymbol{c}_{rr} \otimes I_P \otimes D_{rr} + \left(\boldsymbol{c}_{rr} \otimes I_P + \boldsymbol{c}_{\theta r} \otimes D_\theta + A\boldsymbol{c}_{\theta r} \otimes I_P \right) \otimes R_N^{-1}D_r + (\boldsymbol{c}_{r\theta} \otimes D_\theta + \boldsymbol{c}_{r\theta}A \otimes I_P ) \otimes (D_rR_N^{-1} + R_N^{-2}) + \\ &\quad \left(\boldsymbol{c}_{\theta\theta} \otimes D_{\theta\theta}+ (\boldsymbol{c}_{\theta\theta}A + A\boldsymbol{c}_{\theta\theta}) \otimes D_\theta + A\boldsymbol{c}_{\theta\theta}A \otimes I_P\right) \otimes R_N^{-2}, \quad M = \rho \boldsymbol{I} \otimes I_{NP}, \\ &B_k = \boldsymbol{c}_{rx} \otimes I_{NP}, \quad B_0 = \boldsymbol{c}_{rr} \otimes I_P \otimes D_r + (\boldsymbol{c}_{r\theta} \otimes D_\theta + \boldsymbol{c}_{r\theta}A \otimes I_P) \otimes R_N^{-1}, \end{aligned}\]

where $A = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix}$,

and all the matrices of the discretized polynomial eigenvalue problem are of size $3NP \times 3NP$.

We validate this implementation by comparing predictions for almost flat strips (i.e. with $R\gg w$) to Comsol simulations of the in-plane guided modes for free and clamped flat strips. We also check the influence of $R$ against Comsol simulations, focusing on the first out-of-plane mode, obtaining excellent agreement in both test cases.

Curved strip with approximate coupling to a semi-infinite fluid domain

We discuss the coupling of a curved strip to a semi-infinite fluid medium. We introduce an approximate boundary condition that reads,

\[ \left.\boldsymbol{e}_r \cdot \boldsymbol{C} : \boldsymbol{\nabla u}\right|_{r = R + h/2} = \rho_f \omega^2 \frac{\boldsymbol{u}\cdot \boldsymbol{e}_r}{i\boldsymbol{k}_f\cdot\boldsymbol{e}_r}, \tag{21}\]

where we used the continuity of the normal displacement at the interface, and where $\boldsymbol{k}_f$ is the wavevector of an inhomogeneous plane wave radiated into the fluid. We further simplify this boundary condition by considering that only evanescent waves propagate in the fluid, and by neglecting the transverse wavenumber in the $\boldsymbol{e}_\theta$ direction. Under these assumptions, the approximate boundary condition at $r=R+h/2$ can be discretized as,

\[ \left[(ik_x)^2 B_k + ik_x B_0 + i\rho_f \omega^2 \boldsymbol{I}e_r^T \otimes I_{NP} \right] \cdot \boldsymbol{u} = 0. \tag{22}\]

After inserting the rows of equation (22) in the rows corresponding to collocation points at the top boundary in the discretized eigenvalue problem defined in (20), we can determine the eigenpairs $(k_x,u)$ corresponding to the problem of a curved strip coupled to a semi-infinite fluid medium.