1 Introduction↩︎

Critical velocities are a central concept in high-speed transportation and, more generally, in moving-load and waveguide dynamics [1], [2]. In railway engineering, the phenomenon is commonly associated with a rapid amplification of track and ground vibration when the train speed approaches a characteristic wave speed of the coupled track–embankment–soil system. Early theoretical work already linked this amplification to the generation of strong surface-wave radiation by superfast trains [3]. More recent studies have established that the relevant threshold is governed by the dispersive wave-propagation characteristics of the supporting system and, in particular, by minima of the phase-velocity spectrum rather than by a classical resonance of a finite structure [2], [4], [5]. Specifically, critical velocities occur when the phase velocity \(c_p\) of a propagating mode equals its group velocity \(c_g\). The same physical idea reappears in overhead contact systems, where the so-called catenary barrier reflects the interaction between the contact-point speed and the wave-propagation properties of the tensioned cable system [6], [7]. It also extends naturally to high-speed magnetically levitated vehicles, where wave-induced instability becomes relevant once the operating speed enters the supercritical regime [8]. Beyond transportation, the broader mechanics community has studied how wave speeds can be tailored or even self-controlled in nonlinear media, which further underlines the relevance of robust tools for locating characteristic propagation thresholds [9].

From a modeling perspective, many dynamic systems of interest are invariant, or at least locally periodic, in the direction of motion. Such structures include beams on elastic or viscoelastic foundations, layered soil profiles, periodically supported rails, overhead contact lines, and general prismatic structures. Such systems are naturally described in terms of dispersion relations between frequency and wavenumber. In the railway context, this viewpoint has been employed for critical-speed prediction in coupled track–ground models [2], [4], [5] and was recently extended to periodically varying viscoelastic foundations, where Floquet-type arguments become essential for determining the critical train speed [10]. In parallel, semi-analytical and dynamic-stiffness-based waveguide models have matured considerably. These include approaches based on the thin layer method (TLM) [11], [12], the scaled boundary finite element method (SBFEM) [13], [14], and semi-analytical finite elements (SAFE) [15], [16]. These three methods all employ a numerical discretization of a waveguide’s cross-section, typically based on conventional finite elements, spectral elements, or splines [17].

In contrast, the dynamic-stiffness method employs exact dispersion relations in the wavenumber-frequency domain for periodic structures and plate-built waveguides. Here, it was shown that the Wittrick–Williams algorithm can be adapted and enhanced for exact vibration and wave-propagation analyses, including explicit treatment of the \(J_0\) count, which is crucial for robust and efficient computations over broad frequency ranges [7], [18]–[20]. A common denominator of all the aforementioned approaches is that the underlying physics can be reduced, after semi-discretization in the cross-section, to a parameter-dependent eigenvalue problem of the form \(\mathbf{W}(k,\omega)\,\mathbf{u} = \mathbf{0}\), where the wavenumber \(k\) is the eigenvalue and the frequency \(\omega\) plays the role of an independent parameter, or vice versa. It is worth noting that the dynamic stiffness method generally leads to very small transcendental eigenvalue problems, whereas the semi-analytical methods yield somewhat larger but polynomial eigenvalue problems. Hence, both classes of methods require very different solution procedures and entail distinct benefits and drawbacks. In this paper, we focus on semi-analytical methods and exploit their simpler polynomial structure.

From an algorithmic point of view, critical velocities are usually still obtained indirectly. The standard procedure relies on computing dispersion curves by solving for \(k\) at a sequence of frequencies (or vice versa) to trace the branches, and then determining the critical points by searching for locations that satisfy the condition \(c_p=c_g\). For multimodal systems, this approach is cumbersome and potentially unreliable, because dispersion diagrams may exhibit crossings, osculations, and veering. These issues are well-documented in semi-analytical waveguide models [21], [22]. A closely related problem is the computation of zero-group-velocity (ZGV) points in waveguides, i.e., points on the dispersion curves with \(c_g=0\). It has recently been shown that such ZGV points can be computed directly as the solution of a multiparameter eigenvalue problem [23]. In a nutshell, the parameter-dependent eigenvalue problem mentioned earlier, together with the condition imposed on the group velocity, yields two coupled polynomial eigenvalue problems for the wavenumbers and frequencies at ZGV points. Such multiparameter eigenvalue problems are well understood [24], and numerical methods for the solution are available [25]. Other related applications of this framework were found in the direct computation of leaky waves in layered structures coupled to unbounded media [26], [27] as well as in the computation of Hopf bifurcations in fluid-conveying pipes [28].

Against this background, we will show in this paper that the approach outlined in [23] for ZGV points can be adapted to the computation of critical velocities, which corresponds to incorporating the condition \(c_g=c_p\) (rather than \(c_g = 0\)) into a multiparameter eigenvalue problem. This method allows the direct computation of all critical points in a single step rather than as a posteriori information extracted from traced dispersion branches. The resulting formulation is sufficiently general to cover simple benchmark systems such as a beam on an elastic foundation, layered media representing track support, and more general three-dimensional waveguides, while at the same time opening the door to future extensions to periodic structures, leaky-wave settings, and large coupled systems.

2 Problem statement↩︎

We consider elastic waves propagating along infinitely long structures of uniform cross-section. A classic example is given by guided waves in layered plates of constant thickness, usually referred to as Lamb waves [29]. However, the cross-section can generally be of any shape [30], see Fig. 1 for an example. Wave motion in such structures is characterized by modes, each propagating with a specific frequency-dependent axial wavenumber \(k\) as well as a mode shape \(\mathbf{u}\) that describes the displacement amplitudes on the cross-section.

Generally, \(k\) and \(\mathbf{u}\) are obtained as solutions to a parameter-dependent eigenvalue problem \[\label{eq:evp95discrete} \mathbf{W}(k, \omega)\, \mathbf{u} = \mathbf{0},\tag{1}\] where the frequency \(\omega\) plays the role of an independent parameter. Solving this eigenvalue problem for the wavenumbers \(k\) at varying \(\omega\) results in continuous curves \(k(\omega)\) for each mode, referred to as dispersion curves. Here, \(\mathbf{W}\) is a parameterized \(n\negthinspace \times\negthinspace n\)-matrix whose properties depend on the employed modeling approach, typically involving numerical approximations. A particularly relevant class of methods discretizes the cross-section by finite elements or closely related numerical methods, resulting in a quadratic matrix function of the form [11], [13], [15], [16] \[\label{eq:operator95discrete} \mathbf{W}(k, \omega) = -k^2 \mathbf{L}_2 + \mathrm{i}\mkern 1muk \mathbf{L}_1 + \mathbf{L}_0 + \omega^2 \mathbf{M}\tag{2}\] with real matrices \(\mathbf{L}_i\), \(\mathbf{M}\). We will also briefly discuss the following variant of a fourth-order eigenvalue problem that emerges in modeling, e.g., beams on an elastic foundation [2]: \[\label{eq:operator95discrete95order4} \mathbf{W}(k, \omega) = k^4 \mathbf{L}_4 + \mathbf{L}_0 + \omega^2 \mathbf{M}.\tag{3}\] Note that this problem is linear in \(k^4\) and \(\omega^2\). The formulations leading to Eqs. 2 and 3 are well established in the literature and do not require repetition (see the aforementioned references for details). Here, we are interested in the critical velocities, which are isolated points on the dispersion curves where the phase velocity equals the group velocity, i.e., \[\label{eq:cp95equals95cg} c_p \coloneq \frac{\omega}{k} = \frac{\partial \omega}{\partial k} \eqcolon c_g, \qquad \omega\in\mathbb{R},\quad k \in \mathbb{R}\negthinspace\setminus\negthinspace\{0\}.\tag{4}\] Note that, for the phase and group velocities to be well-defined, we assume real-valued frequencies and wavenumbers, which is true for propagating modes in lossless media. For a given eigenvector/eigenvalue pair corresponding to Eq. 2 and assuming a Hermitian structure¹, the group velocity can be computed a posteriori as \[\label{eq:groupVel} c_g = \frac{\mathbf{u}^{\mathsf{H}}(2k\mathbf{L}_2 - \mathrm{i}\mkern 1mu\mathbf{L}_1)\mathbf{u}}{2\omega\,\mathbf{u}^{\mathsf{H}}\mathbf{M}\mathbf{u}}.\tag{5}\]

3 Computation of critical velocities↩︎

To derive a formulation for the direct computation of critical velocities, we will differentiate the eigenvalue problem 1 with respect to the wavenumber and substitute the condition 4 . The resulting equation, together with the original eigenvalue problem, poses two coupled eigenvalue problems (i.e., a two-parameter eigenvalue problem), whose solutions contain the sought critical points. The basic formulation is similar to that previously presented for the computation of ZGV points [23], [31], except for a slightly different structure arising from an additional term due to the nonzero group velocity. Note that, in the following, the dash symbol always refers to derivatives with respect to the wavenumber \(k\).

We begin with the matrix operator defined by Eq. 2 . Differentiating the eigenvalue problem with respect to the wavenumber yields \[\mathbf{W}'\, \mathbf{u} + \mathbf{W}\, \mathbf{u}' = \mathbf{0}\] with \[\mathbf{W}'(k, \omega) = -2k\,\mathbf{L}_2 + \mathrm{i}\mkern 1mu\mathbf{L}_1 + 2\omega\,\omega'\,\mathbf{M}.\] Note that \(\omega'\) is the group velocity. At a critical point, we substitute the condition given by Eq. 4 and define the resulting matrix function as \[\mathbf{W}_c(k, \omega) \mathrel{\vcenter{:}}= -2k\,\mathbf{L}_2 + \mathrm{i}\mkern 1mu\mathbf{L}_1 + 2\,\frac{\omega^2}{k}\,\mathbf{M}.\] Hence, a solution of the eigenvalue problem with the additional property \(c_g=c_p\) must satisfy \(\mathbf{W}_c\, \mathbf{u} + \mathbf{W}\, \mathbf{u}' = \mathbf{0}\) as well as \(\mathbf{W}\, \mathbf{u} = \mathbf{0}\), or in block matrix form: \[\label{eq:criticalPointSystem} \begin{bmatrix} \mathbf{W} & \mathbf{0} \\ \mathbf{W}_c & \mathbf{W} \end{bmatrix} \begin{bmatrix} \mathbf{u} \\ \mathbf{u}' \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix}.\tag{6}\] We multiply the second equation by \(k\) to separate the parameters: \[\label{eq:criticalPointSystemK} \begin{bmatrix} \mathbf{W} & \mathbf{0} \\ k\mathbf{W}_c & \mathbf{W} \end{bmatrix} \begin{bmatrix} \mathbf{u} \\ k\mathbf{u}' \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix}.\tag{7}\] Substituting the explicit expressions for \(\mathbf{W}\) and \(\mathbf{W}_c\) yields \[\label{eq:criticalPointSystemExplicit} \left( -k^2 \begin{bmatrix} \mathbf{L}_2 & \mathbf{0} \\ 2\mathbf{L}_2 & \mathbf{L}_2 \end{bmatrix} + \mathrm{i}\mkern 1muk \begin{bmatrix} \mathbf{L}_1 & \mathbf{0} \\ \mathbf{L}_1 & \mathbf{L}_1 \end{bmatrix} + \begin{bmatrix} \mathbf{L}_0 & \mathbf{0} \\ \mathbf{0} & \mathbf{L}_0 \end{bmatrix} + \omega^2 \begin{bmatrix} \mathbf{M} & \mathbf{0} \\ 2\mathbf{M} & \mathbf{M} \end{bmatrix} \right) \begin{bmatrix} \mathbf{u} \\ k\mathbf{u}' \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix},\tag{8}\] which we abbreviate as \[\label{eq:criticalPointSystemShort} \big(-k^2 \mathbfcal{L}_2 + \mathrm{i}\mkern 1muk \mathbfcal{L}_1 + \mathbfcal{L}_0 + \omega^2 \mathbfcal{M}\big) \mathbf{v} = \mathbf{0}.\tag{9}\] The above equation constitutes an additional condition for an eigenvalue solution to be a critical point. This equation, together with the original eigenvalue problem 1 form a system of two eigenvalue problems in \(\omega\) and \(k\), whose solutions correspond to critical points. In other words, we obtain a quadratic two-parameter eigenvalue problem: \[\tag{10} \begin{align} (-k^2 \mathbf{L}_2 + \mathrm{i}\mkern 1muk \mathbf{L}_1 + \mathbf{L}_0 + \omega^2 \mathbf{M})\mathbf{u} & = \mathbf{0},\tag{11} \\ (-k^2 \mathbfcal{L}_2 + \mathrm{i}\mkern 1muk \mathbfcal{L}_1 + \mathbfcal{L}_0 + \omega^2 \mathbfcal{M}) \mathbf{v} & = \mathbf{0}. \tag{12} \end{align}\] For simplification and conciseness, it is rewritten as a linear three-parameter eigenvalue problem, defining the parameters \(\eta = (\mathrm{i}\mkern 1muk)^2\), \(\lambda = \mathrm{i}\mkern 1muk\), \(\mu = \omega^2\) and including a third equation to define the relationship between \(\eta\) and \(\lambda\): \[\tag{13} \begin{align} (\eta \mathbf{L}_2 + \lambda \mathbf{L}_1 + \mathbf{L}_0 + \mu \mathbf{M})\mathbf{u} & = \mathbf{0},\tag{14} \\ (\eta \mathbfcal{L}_2 + \lambda \mathbfcal{L}_1 + \mathbfcal{L}_0 + \mu \mathbfcal{M})\mathbf{v} & = \mathbf{0},\tag{15} \\ (\eta \mathbf{C}_2 + \lambda \mathbf{C}_1 + \mathbf{C}_0) \mathbf{w} & = \mathbf{0}\tag{16} \end{align}\] with \[\mathbf{C}_2 = \begin{bmatrix} 1 & 0 \\ 0& 0 \end{bmatrix},\; \mathbf{C}_1 = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix},\; \mathbf{C}_0 = \begin{bmatrix} 0 & 0 \\ 0 & 1 \end{bmatrix}.\] It can be easily verified that 16 incorporates the equation \[\eta = \lambda^2,\] because a solution to this eigenvalue problem satisfies \[\det(\eta \mathbf{C}_2 + \lambda \mathbf{C}_1 + \mathbf{C}_0) = \eta - \lambda^2 = 0.\] The system 13 is similar to those discussed in more detail in [23], [28], and we use the MATLAB implementation available at [25] to solve it. To briefly summarize the essential steps, we find solutions to the three-parameter EVP by solving related generalized eigenvalue problems (GEPs), see, e.g., [32] and the previous applications [23], [28]. These GEPs are \[\label{eq:delta95system} \mathbf{\Delta}_1 \mathbf{z} = \lambda \mathbf{\Delta}_0 \mathbf{z},\quad \mathbf{\Delta}_M \mathbf{z} = \mu \mathbf{\Delta}_0 \mathbf{z},\quad \mathbf{\Delta}_2 \mathbf{z} = \eta \mathbf{\Delta}_0 \mathbf{z},\tag{17}\] with the eigenvector \(\mathbf{z}\) obtained by the Kronecker product \(\mathbf{z}=\mathbf{u}\otimes \mathbf{v}\otimes \mathbf{w}\). The matrices in the above equations are referred to as operator determinants and are computed using the Kronecker product as: \[\begin{align} {2}\label{eq:delta01} \mathbf{\Delta}_0 = &\left|\begin{matrix} \mathbf{L}_2 & \mathbf{L}_1 & \mathbf{M} \cr \mathbfcal{L}_2 & \mathbfcal{L}_1 & \mathbfcal{M}\cr \mathbf{C}_2 & \mathbf{C}_1 & \mathbf{0} \end{matrix}\right|_\otimes\negthickspace\negthickspace,\quad \mathbf{\Delta}_1 = -&\left|\begin{matrix} \mathbf{L}_2 & \mathbf{L}_0 & \mathbf{M}\cr \mathbfcal{L}_2 & \mathbfcal{L}_0 & \mathbfcal{M}\cr \mathbf{C}_2 & \mathbf{C}_0 & \mathbf{0} \end{matrix}\right|_\otimes \negthickspace\negthickspace, \nonumber \\ \mathbf{\Delta}_2 = -&\left|\begin{matrix} \mathbf{L}_1 & \mathbf{L}_0 & \mathbf{M}\cr \mathbfcal{L}_1 & \mathbfcal{L}_0 & \mathbfcal{M}\cr \mathbf{C}_1 & \mathbf{C}_0 & \mathbf{0} \end{matrix}\right|_\otimes \negthickspace\negthickspace,\quad \mathbf{\Delta}_M= -&\left|\begin{matrix} \mathbf{L}_2 & \mathbf{L}_1 & \mathbf{L}_0 \cr \mathbfcal{L}_2 & \mathbfcal{L}_1 & \mathbfcal{L}_0 \cr \mathbf{C}_2 & \mathbf{C}_1 & \mathbf{C}_0 \end{matrix}\right|_\otimes \negthickspace\negthickspace. \end{align}\tag{18}\] Solving the first equation in 17 yields candidates for \(\eta\), \(\lambda\), \(\mu\) corresponding to critical points. Recovering the wavenumber \(k = -\mathrm{i}\mkern 1mu\lambda\) and substituting into the original parameter-dependent eigenvalue problem 1 provides the corresponding frequencies \(\omega\). It must be noted that 17 typically admits many additional solutions that are not critical points. Hence, we employ a simple postprocessing step, in which we select solutions with finite frequency that satisfy \(c_p\approx c_g\) to a given tolerance. Specifically, if \(k_\mathrm{cand}\), \(\omega_\mathrm{cand}\) are candidate solutions, we consider them a valid critical point if \(\omega_\mathrm{cand}\) is finite and \[\left|\frac{|c_{p,\mathrm{cand}}|}{|c_{g,\mathrm{cand}}|}-1\right| < 10^{-4}\] with \(c_{p,\mathrm{cand}} = \frac{\omega_\mathrm{cand}}{k_\mathrm{cand}}\) and \(c_{g,\mathrm{cand}}\) computed using Eq. 5 .

Since the GEPs are singular, the MATLAB toolbox MultiParEig [25] uses a rank-projection algorithm as proposed in [33]. Due to their construction via Kronecker products, the size \(n_\mathbf{\Delta}\) of the operator determinants increases quadratically with \(n\) (the size of the original parameter-dependent eigenvalue problem). Specifically, it is given by the product of the sizes of the three coupled eigenvalue problems in 13 , in our case \(n_\mathbf{\Delta} = n\negthinspace\cdot\negthinspace 2n\negthinspace\cdot\negthinspace 2 = 4n^2\). Computing all eigenvalues of the GEPs 17 is feasible on a current laptop computer for matrix sizes \(n\) up to about 50, which already corresponds to \(n_\mathbf{\Delta} = 10000\). In Fig. [fig:flowChart], we present an overview of the essential steps in computing the critical points based on the proposed approach.

3.1 Special case of a quartic eigenvalue problem↩︎

In principle, the approach outlined above can be straightforwardly extended to polynomial eigenvalue problems of arbitrary order (though the computational costs increase rapidly with the number of parameters). We present a special case of a fourth-order eigenvalue problem given by Eq. 3 , as it is encountered in the well-known formulation of a beam on an elastic foundation. Differentiating the eigenvalue problem analogously and introducing \(\xi = k^4\), \(\mu = \omega^2\) leads to the two-parameter eigenvalue problem for critical points: \[\tag{19} \begin{align} (\xi \mathbf{L}_4 + \mathbf{L}_0 + \mu \mathbf{M})\mathbf{u} & = \mathbf{0},\tag{20} \\ (\xi \mathbfcal{L}_4 + \mathbfcal{L}_0 + \mu \mathbfcal{M})\mathbf{v} & = \mathbf{0}\tag{21} \end{align}\] with \[\mathbfcal{L}_4 = \begin{bmatrix} \mathbf{L}_4 & \mathbf{0} \\ 4\mathbf{L}_4 & \mathbf{L}_4 \end{bmatrix}\] and \(\mathbfcal{L}_0\), \(\mathbfcal{M}\) defined as in 8 . Note that this simpler case leads to a two-parameter eigenvalue problem, as the original problem is already linear in \(k^4\). Consequently, the operator determinants are of size \(n_\mathbf{\Delta} = 2n^2\).

4 Numerical experiments↩︎

4.1 Illustrative example↩︎

We begin with an academic example featuring a simple solution that is easy to verify. Consider the parameter-dependent eigenvalue problem \[\label{eq:evp95minimalExample} (k^2 \mathbf{L}_2 + k \mathbf{L}_1 + \mathbf{L}_0 + \omega^2 \mathbf{M})\, \mathbf{u} = \mathbf{0}\tag{22}\] with \[\mathbf{L}_2 = -\begin{bmatrix} 153 & 64 \\ 64 & 57 \end{bmatrix},\quad \mathbf{L}_1 = \begin{bmatrix} 24 & -8 \\ -8 & 36 \end{bmatrix},\quad \mathbf{L}_0 = -\begin{bmatrix} 20 & 0 \\ 0 & 20 \end{bmatrix},\quad \mathbf{M} = \begin{bmatrix} 17 & 6 \\ 6 & 8 \end{bmatrix}.\] The matrix operator is of the form 2 , except that we define it directly in \(k\) rather than \(\mathrm{i}\mkern 1muk\) for conciseness. Solutions are given by the four eigencurves (including the positive and negative signs): \[\label{eq:minimalEigencurves} \begin{align} \omega_1 & = \pm \sqrt{5\,k^2 - 8\,k + 4}, \\ \omega_2 & = \pm \sqrt{9.25\,k^2 - k + 1}. \end{align}\tag{23}\] The modes each exhibit one critical point at \(k_1 = 1\) and \(k_2 = 2\), respectively, i.e., \[\omega_1'(k_1) = \frac{\omega_1(k_1)}{k_1} = 1,\quad \omega_2'(k_2) = \frac{\omega_2(k_2)}{k_2} = 3.\] Figure 3 shows the analytical eigencurves in terms of \(\omega(k)\), \(c_p(k)\), and \(c_g(k)\), together with the critical points (listed in Tab. 1), computed using the proposed approach. In this simple problem, the critical values are obtained exactly to machine precision.

4.2 Beam on an elastic foundation↩︎

Let us consider a classical scalar benchmark example: a beam of infinite length on an elastic foundation. This classical problem is often used as a model for rails on soil layers, and the critical velocities of such simplified systems are well understood, see, e.g., [2] and the references therein. Dispersion curves in terms of \(k\) and \(\omega\) are obtained as solutions to the scalar problem \[\label{eq:beamDispersion} \omega^2 - \omega_0^2 - k^4c_r^2R^2 = 0\tag{24}\] with the definitions \[\omega_0^2 \coloneq \frac{K}{\rho A},\qquad c_r^2 \coloneq \frac{E}{\rho},\qquad R^2 \coloneq \frac{I}{A},\] where \(E\), \(A\), \(\rho\), \(I\) denote the beam’s properties (Young’s modulus, cross-sectional area, mass density, moment of inertia), and \(K\) is the foundation’s stiffness per length. While this scalar problem clearly does not require numerical methods to obtain the critical speed, it is instructive for demonstrating the proposed method. In our notation established in Section 3, we treat 24 as the scalar version of 3 and pose the corresponding multiparameter eigenvalue problem of the form 19 with \[L_4 = -c_r^2R^2,\quad L_0 = -\omega_0^2, \quad M = 1,\] \[\mathbfcal{L}_4 = -c_r^2R^2 \begin{bmatrix} 1 & 0 \\ 4 & 1 \end{bmatrix},\quad \mathbfcal{L}_0 = -\omega_0^2\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix},\quad \mathbfcal{M}= \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}.\] The operator determinants are obtained as \[\mathbf{\Delta}_0 = -c_r^2R^2\begin{bmatrix} 0 & 0 \\ 2 & 0 \end{bmatrix},\quad \mathbf{\Delta}_1 = -\omega_0^2\begin{bmatrix} 0 & 0 \\ 2 & 0 \end{bmatrix},\quad \mathbf{\Delta}_2 = -c_r^2R^2\omega_0^2\begin{bmatrix} 0 & 0 \\ 4 & 0 \end{bmatrix}.\] In this simple case, it can be immediately observed that the operator determinants are, in fact, singular. Wavenumbers and frequencies corresponding to critical points satisfy the eigenvalue problems \[\mathbf{\Delta}_1 \mathbf{x} = k_c^4 \mathbf{\Delta}_0 \mathbf{x},\quad \mathbf{\Delta}_2 \mathbf{x} = \omega_c^2 \mathbf{\Delta}_0 \mathbf{x}.\] The above eigenvalue problems have exactly one finite solution at \[k_c^4 = \frac{\omega_0^2}{c_r^2\,R^2},\qquad \omega_c^2=2\omega_0^2,\qquad c_{c}=\frac{\omega_c}{k_c} =\sqrt{2\,\omega_0\,c_r R}.\] The dispersion curves obtained by computing wavenumbers \(k\) at discrete values of the frequency \(\omega\), as well as the critical point obtained by solving the multiparameter eigenvalue problem by virtue of the proposed approach, are presented in Fig. 4. The results are consistent with those presented in [2].

4.3 Layered medium↩︎

Moving to a more realistic, practical case, this example is adapted from that presented by Kausel [2], who studied critical velocities in stratified infinite media representing the foundation beneath rail systems. Specifically, we consider a system consisting of two layers of thicknesses 2 m and 3 m, respectively, representing ballast and embankment with the material parameters listed in Table 2 (in terms of shear wave velocity \(c_s\), mass density \(\rho\), and Poisson’s ratio \(\nu\)). The two layers are fully coupled (i.e., assuming continuous displacements across the interface), and the outer boundaries are traction-free. Following the semi-analytical approach described in [13], each of the two layers is discretized by only one higher-order finite element. For the presented frequency range up to \(60\,\unit{Hz}\) (\(377\,\unitfrac{rad}{s}\)), an element order of five is found to be sufficient, leading to a total of 22 degrees of freedom in a plane strain model. The matrices and dispersion curves have been computed using the open-source MATLAB code SAMWISE [34]. Results are presented in Fig. 5, showing three critical points. Interestingly, two of the critical points belong to the same mode and lie on a branch with small curvature. The numerical values of the three critical points are listed in Table 3.

4.4 Three-dimensional structure↩︎

As a final test, we consider a thin-walled prismatic waveguide whose cross-section is depicted in Fig. 6. The structure is assumed to be of infinite extent in the \(x\)-direction. Thin-walled built-up members of this type are representative of lightweight transport structures, where both in-plane and out-of-plane deformation mechanisms coexist. This example is therefore more demanding than the preceding layered-medium benchmark: it involves symmetries, dispersion branches interact closely, and multiple stationary points of the phase velocity occur within a moderate frequency range. Hence, the structure provides a useful test case for assessing whether the proposed multiparameter formulation can identify critical points without relying on branch tracing. The total width and height of the cross-section are 140 mm and 50 mm, respectively, and the thickness of each element is 4.2 mm. The material is defined by a Young’s modulus of 71 GPa, a mass density of 2700 kg/m\(^3\), and a Poisson’s ratio of 0.332. As the current approach requires the matrix size of the discretized model to be small, we use a coarse discretization, which nevertheless gives sufficiently accurate results for frequencies up to 5 kHz (\(\require{upgreek} 10^4\,\uppi\,\unitfrac{rad}{s}\)). Since the cross-section is symmetric with respect to the \(y\)-axis, we discretize half of it and apply, in turn, symmetric and antisymmetric boundary conditions (see [30] for implementation details). Half of the cross-section is discretized by five finite elements. Note that the elements use a ‘mixed-order-approach’ with linear interpolation across the thickness direction and quadratic interpolation along each member of the cross-section. This discretization leads to 18 nodes and, consequently, 54 degrees of freedom (50 and 46 after applying symmetric and antisymmetric boundary conditions, respectively). Again, the matrices and dispersion curves were computed using the software SAMWISE [34]. Results are presented in Fig. 7. This structure exhibits eight critical points within the selected frequency range. Numerical values of the computed critical points are listed in Table 4.

5 Conclusion↩︎

We have presented a direct approach to computing critical velocities in elastic waveguides. By incorporating the condition of coinciding phase and group velocities into the parameter-dependent eigenvalue problem arising in typical semi-analytical approaches, the identification of critical points is reformulated as a multiparameter eigenvalue problem. This formulation avoids the need for dispersion-curve tracing or complicated postprocessing.

We employed previously established approaches and algorithms to solve the resulting singular polynomial multiparameter eigenvalue problems in a straightforward manner. The numerical examples demonstrate that the proposed approach reliably identifies critical points in different scenarios, including simple benchmark systems, layered media, and three-dimensional prismatic structures. In particular, the method remains effective in situations with multiple critical points, mode crossings, and nearly flat (low-curvature) minima in the dispersion curves. Currently, the main drawback of the proposed approach is that its computational cost increases rapidly with matrix size, making it de facto feasible only for small systems. Hence, future work will focus on improving computational efficiency by assessing strategies similar to those proposed for computing ZGV points [23], [31], as well as on extending the framework to more general settings, including periodically structured media and leaky-wave problems.

References↩︎