Tailoring a Three-Layer Track Model to Delay Instability and Minimize Critical Velocity Effects at Very High Velocities

1. Introduction

The railway enables the circulation of railway vehicles in a safe and comfortable manner. Prioritizing rail transport over road transport represents an effective approach to combating climate change and reducing the carbon footprint. In line with this, it is necessary to increase the capacity of the railway network by raising circulation velocity without compromising passenger safety and comfort. To achieve this, an adequate transfer of solicitations to the ground must be ensured. Given the high demands that rolling stock imposes on the rail, one of the most challenging problems is identifying the most suitable technical solutions capable of supporting an increase in traffic velocity under standardized conditions without compromising passenger comfort or vehicle ride quality.

The maximum attainable circulation velocity in railway systems is often constrained by the so-called critical velocity, which, in the context of a moving force, represents the lowest velocity at which waves can propagate through the supporting track structure. This velocity arises from the dynamic interplay between all structural components and corresponds to a resonance condition. Experimental data indicate that this threshold has already been reached under real operating conditions [1].

To deepen our understanding of such dynamic behavior, further investigations are conducted within the well-established framework of moving load problems—a subject that continues to attract significant interest in transportation engineering. To enhance computational tractability and facilitate optimization and sensitivity analyses, many studies adopt simplified models. These simplifications can pertain to either the structure, the vehicle, or both, depending on the specific aspect under study. Besides exploring the effect of key parameters on dynamic response, major areas of interest include resonance-related critical velocities and frequencies, system instability, and the interaction between moving neighboring inertial bodies.

Common structural simplifications have led to layered models comprising beams, concentrated masses, and spring-damper elements. These configurations are favored for their simplicity and computational efficiency, finding widespread application in both research and practice. A survey of foundation modeling approaches is available in the classical work [2], while the discrete–continuous modeling debate is discussed in [3]. Rail tracks are frequently represented as Euler–Bernoulli infinite beams supported by continuous elastic foundations with one [4,5,6,7], two [8,9], or three layers [10,11]. Alternatively, the Euler–Bernoulli formulation may be replaced with a Timoshenko beam model to account for shear deformation and rotary inertia [12,13,14]. To accurately reflect the influence of sleeper spacing, models incorporating discrete supports are also used [15]. Despite their simplicity, such idealized models have successfully captured various phenomena in railway dynamics, as demonstrated in [16]. Comprehensive studies on vehicle–track interactions can be found in [17,18,19], further affirming the practical value of simplified modeling approaches. A taxonomy of layered models is proposed in [20], while their experimental validation is discussed in [21]. Naturally, the inclusion of additional layers improves the model’s ability to capture real-world behavior.

The feasibility of three-layer models in railway applications is evaluated in [10], and comparative features of one-, two-, and three-layer systems are reviewed in [22]. Research aimed at determining model parameters often builds upon the ballast pyramid concept [23,24] and the Saller hypothesis from 1932 [25]. A detailed treatment of the pyramid model—referred to as the "stress cone model"—is presented in [26], where it is used to derive vertical ballast stiffness and identify dynamically activated ballast mass. Refinements include overlapping cones between adjacent sleepers in both longitudinal [27] and lateral directions [10].

A number of studies concerning constant forces acting on an infinite structure exploit semianalytical methods. Fully analytical formulations are derived in [28,29]. The objectives of these two works are very similar, with some differences in the characterization of the Winkler–Pasternak foundation and critical damping. Both papers are also focused on critical velocity.

Advanced methods such as wavelet transforms are used in [30] to investigate moving forces on a double-beam system. The Adomian decomposition method is applied to handle nonlinearity in the viscoelastic connection between beams, with the model later adapted into a two-layer version [31] and extended to include nonlinear foundation behavior [32,33]. These studies directly derive steady-state solutions and compare their results with experimental observations.

Other works address the limitations of classical beam theories [34,35,36]. Large deformations are handled in [37]; initial conditions are accounted for in [38]; non-uniform motion is analyzed in [39]. Related problems include changes in foundation stiffness [40,41] or analyses concerned with non-classical ballasted tracks [42].

Further studies related to this paper focus on sleeper spacing [43,44] or critical zones [45], where more accurate mechanical descriptions of the track components are necessary, thus requiring numerical methods. When dealing with foundation models, more general formulations should be applied [46].

Recent studies focus on the instability of moving inertial objects. In the book chapter [47], vibrational stability of mechanical oscillators on continuous beam–foundation systems is analyzed using the D-decomposition method, comparing two models: one with a stabilizer directly connected to the car body, and another with interconnected primary suspensions. Instability of a moving bogie, with the possibility of occurring in the subcritical velocity range, is addressed in [48]. Further subjects related to instability are discussed in the book chapters [49,50]. In [51], the importance of material characterization is demonstrated through laboratory (physical and mechanical) tests and finite element modeling using four different rock materials as ballast, aiming to reduce the use of stone materials and promote an optimized, sustainable concept for new railway projects. The effects of ballast thickness on subgrade stresses are analyzed in terms of subgrade bearing capacity and permanent deformation of geotechnical materials.

In [52], localized waves in an Euler–Bernoulli beam on a nonlinear elastic foundation with nonlinear structural damping under the action of a moving force are studied. New approximation formulas describing both linear and nonlinear damping characteristics are derived. The Ritz method and perturbation techniques are employed for this purpose. It is demonstrated that the nonlinear contribution reduces displacement amplitudes compared to the linear foundation; however, it also lowers the critical velocity. Additionally, there is never true resonance—displacements remain finite even in the absence of damping.

Machine learning and numerical approaches are compared in [53] with the aim of evaluating track performance. In [54], a novel perspective is presented for improving performance and reliability through material optimization based on wheel–rail interactions under variable load and track conditions in a Y25 bogie. As this involves more detailed numerical studies, finite element analysis is implemented.

In this paper, the geometry and material parameters of a three-layer railway track model are tailored to achieve favorable properties for the circulation of high-speed trains at very high velocities. Specifically, cases where the lowest critical velocity is replaced by a pseudo-critical value are optimized so that the vibration increase around this pseudo-critical velocity is negligible, and the next resonance as well as the onset of instability are as distant as possible. All results are presented (semi)analytically with very high precision and in dimensionless parameters to cover a wide range of scenarios. Then, real data are assigned to the optimal cases.

The paper is organized in the following way. In Section 2, the model, involved parameters, and the governing equations are presented. In Section 3, the terminology used is introduced, and issues related to the number of resonances and the role of critical and pseudo-critical velocities, including their similarities and differences, are explained, followed by the motivation for this work. In Section 4, first, one case study is presented, followed by several methods and optimization techniques such as design of experiments, Monte Carlo analysis, parametric search, and simulated annealing. Optimized data are connected to real parameters and compared. The present paper is concluded by a discussion of the results obtained, their utility, and future challenges. After that, conclusions are drawn. The paper includes one Appendix where all numerical data are summarized for ease of comparison.

2. Governing Equations and Parameters

The three-layer railway track model is quite popular due to its simplicity, easy manipulability, computational efficiency, and good agreement with reality [10]. It is schematically presented in Figure 1 [11].

To facilitate the analytical steps of the solution, all parameters are introduced into the governing equations in their distributed and dimensionless forms. Exploiting the symmetry, the model represents a half-track of a single-line straight section.

The governing equations in dimensionless parameters and moving coordinates read as (more details can be found in [11]):

$$\begin{array}{l}{\tilde{w}}_{,\xi \xi \xi \xi }+4\left({\tilde{w}}_{,\tau \tau }-2\alpha {\tilde{w}}_{,\xi \tau }+{\alpha }^2{\tilde{w}}_{,\xi \xi }\right)+4{\eta }_N{\tilde{w}}_{,\xi \xi }+8{\eta }_p\left({\tilde{w}}_{,\tau }-\alpha {\tilde{w}}_{,\xi }-{\tilde{u}}_{s,\tau }+\alpha {\tilde{u}}_{s,\xi }\right)\\ +4{\kappa }_p\left(\tilde{w}-{\tilde{u}}_s\right)=\left(8{\eta }_P-4{\eta }_M{\tilde{w}}_{,\tau \tau }\right)\delta \left(\xi \right)\end{array}$$

(1)

$$\begin{array}{l}{\mu }_s\left({\tilde{u}}_{s,\tau \tau }-2\alpha {\tilde{u}}_{s,\xi \tau }+{\alpha }^2{\tilde{u}}_{s,\xi \xi }\right)-2{\eta }_p\left({\tilde{w}}_{,\tau }-\alpha {\tilde{w}}_{,\xi }-{\tilde{u}}_{s,\tau }+\alpha {\tilde{u}}_{s,\xi }\right)\\ -{\kappa }_p\left(\tilde{w}-{\tilde{u}}_s\right)+2{\eta }_b\left({\tilde{u}}_{s,\tau }-\alpha {\tilde{u}}_{s,\xi }-{\tilde{u}}_{b,\tau }+\alpha {\tilde{u}}_{b,\xi }\right)+{\kappa }_b\left({\tilde{u}}_s-{\tilde{u}}_b\right)=0\end{array}$$

(2)

$$\begin{array}{l}{\mu }_b\left({\tilde{u}}_{b,\tau \tau }-2\alpha {\tilde{u}}_{b,\xi \tau }+{\alpha }^2{\tilde{u}}_{b,\xi \xi }\right)-{\eta }_s{\tilde{u}}_{b,\xi \xi }-2{\eta }_b\left({\tilde{u}}_{s,\tau }-\alpha {\tilde{u}}_{s,\xi }-{\tilde{u}}_{b,\tau }+\alpha {\tilde{u}}_{b,\xi }\right)\\ -{\kappa }_b\left({\tilde{u}}_s-{\tilde{u}}_b\right)+2{\eta }_f\left({\tilde{u}}_{b,\tau }-\alpha {\tilde{u}}_{b,\xi }\right)+{\tilde{u}}_b=0\end{array}$$

(3)

where

$$\chi =\sqrt[4]{{\displaystyle \frac{k_f}{4EI}}},v_{ref}=\sqrt[4]{{\displaystyle \frac{4k_{f}EI}{m^2}}}={\displaystyle \frac{1}{\chi }}\sqrt{{\displaystyle \frac{k_f}{m}}},\alpha ={\displaystyle \frac{v}{v_{ref}}}$$

(4)

$${\mu }_s={\displaystyle \frac{m_s}{l_{s}m}},{\mu }_b={\displaystyle \frac{m_b}{l_{s}m}},{\kappa }_p={\displaystyle \frac{k_p}{l_s{k}_f}},{\kappa }_b={\displaystyle \frac{k_b}{l_s{k}_f}},{\eta }_s={\displaystyle \frac{k_s{l}_s}{2\sqrt{EI{k}_f}}}$$

(5)

$${\eta }_N={\displaystyle \frac{N}{2\sqrt{EI{k}_f}}},{\eta }_p={\displaystyle \frac{c_p}{2l_s\sqrt{m{k}_f}}},{\eta }_b={\displaystyle \frac{c_b}{2l_s\sqrt{m{k}_f}}},{\eta }_f={\displaystyle \frac{c_f}{2\sqrt{m{k}_f}}}$$

(6)

where $EI$ and $m$ describe the Euler-Bernoulli beam on top of the model by its bending stiffness and mass per unit length. Further, $l_s$ is the sleeper spacing, $m_s$ and $m_b$ are masses of the sleeper (half) and ballast cone, $k_p$, $k_b$, $k_f$ and $k_s$ are the vertical stiffnesses of the rail pads, ballast cone, and foundation, and the shear stiffness of the ballast layer, respectively, where $k_f$ is usually defined as distributed, while the other values are discrete. The ballast cone theory according to [26], and refined in [27,10], is essential in the definition of the governing parameters. $N$ is the axial force and $P$ the vertical force acting on a moving mass $M$ that is traversing the model at constant velocity $v$. Unknown displacement fields of the beam $w$, the sleeper $u_s$, and the ballast $u_b$ are functions of spatial coordinate $x$ and time $t$. By switching to moving coordinates $x$ passes to $r=x-vt$ and then to their dimensionless counterparts:

$$\xi =\chi r,\tau =\chi v_{ref}t$$

(7)

further

$${\eta }_M={\displaystyle \frac{M\chi }{m}},{\eta }_P={\displaystyle \frac{P}{P}}$$

(8)

and other parameters in Equations (5) and (6) are masses (${\mu }_s$ and ${\mu }_b$), stiffnesses (${\kappa }_p$, ${\kappa }_b$ and ${\eta }_s$), damping ratios (${\eta }_p$, ${\eta }_b$ and ${\eta }_f$) and the axial force ratio (${\eta }_N$). To complete the list of symbols $\delta$ is the Dirac delta function, and derivatives are denoted by the symbol of the respective variable, preceded by a comma.

This means that a Winkler beam characterized by $EI$, $m$ and $k_f$ is selected as the reference beam. $\chi$ is inverse of its characteristic length, $v_{ref}$ is the critical velocity of a moving force, and $\alpha$ is therefore the velocity ratio. For identification of the moving mass instability, it is assumed that the mass is always in contact with the beam. Finally,

$$w_{ref}={\displaystyle \frac{P\chi }{2k_f}},\tilde{w}={\displaystyle \frac{w}{w_{ref}}},{\tilde{u}}_s={\displaystyle \frac{u_s}{w_{ref}}},{\tilde{u}}_b={\displaystyle \frac{u_b}{w_{ref}}}.$$

(9)

Considering the range of admissible parameters specified in detail in [11], it is possible to conclude that the admissible ranges of the leading parameters are approximately as follows:

$${\mu }_s\in 〈1;6〉,{\mu }_b\in 〈2;45〉,{\kappa }_p\in 〈0.03;4.2\times {10}^4〉,{\kappa }_b\in 〈0.04;6\times {10}^3〉$$

(10)

Nevertheless, it is also necessary to point out that there are general trends such as enlarging sleeper spacing, the introduction of novel lightweight materials, and other cost-saving measures, which can obviously affect these ranges.

Admissible intervals in Equation (10) were obtained in the following way: first, geometrical data needed for the stress cone theory were assumed as:

$$l_s=〈0.5;0.8〉,h_b=〈0.2;0.6〉,l_b=〈0.2;0.36〉,{\alpha }_b=〈{25}^{\mathrm{o}};{50}^{\mathrm{o}}〉$$

(11)

All values except for ${\alpha }_b$ are given in [m]. $h_b$, $l_b$ and ${\alpha }_b$ are ballast height, base sleeper width, and angle of the ballast cone, as exemplified in Figure 2, [10].

Further, $l_g$ and $l_e$ are the rail gauge and the sleeper effective length. For the calculations, $l_e=1\mathrm{m}$ is used as the most common value, and so $l_g$ does not have to be introduced. Then, by using the ranges from Equation (11), the dynamically activated ballast volume $m_{cn}$ is within the range (0.065–0.7) m³, and the parameter related to the vertical stiffness $k_{cn}$ ranges within (0.79–3.36) m (multiplication by the ballast Young´s modulus yields the correct base unit of N/m). According to [10,11], the range (40–160) kg is considered for $m_s$, and the ballast density ${\rho }_b$ is taken within the range (1200–2600) kg/m³. Further, values within the range (20–5000) MN/m are considered for $k_p$, and (50–400) MPa for the ballast Young’s modulus, $E_b$.

There is an obvious connection between ${\mu }_s$ and ${\mu }_b$, as well as between ${\kappa }_p$ and ${\kappa }_b$, as they are both connected to the reference beam. Their ratios must be bounded by the values determined as:

$${\left({\displaystyle \frac{{\mu }_s}{{\mu }_b}}\right)}^{+}={\displaystyle \frac{m_s^{+}}{m_{cn}^{-}{\rho }_b^{-}}}={\displaystyle \frac{160}{0.065\cdot 1200}}=2.05,{\left({\displaystyle \frac{{\mu }_s}{{\mu }_b}}\right)}^{-}={\displaystyle \frac{m_s^{-}}{m_{cn}^{+}{\rho }_b^{+}}}=0.022$$

(12)

$${\left({\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}\right)}^{+}={\displaystyle \frac{k_p^{+}}{k_{cn}^{-}{E}_b^{-}}}={\displaystyle \frac{5\times {10}^9}{0.79\cdot 5\times {10}^7}}=126.58,{\left({\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}\right)}^{-}={\displaystyle \frac{k_p^{-}}{k_{cn}^{+}{E}_b^{+}}}={\displaystyle \frac{2\times {10}^7}{3.36\cdot 4\times {10}^8}}=0.015$$

(13)

In what follows, “+” and “–“ will be used to identify the upper and lower bounds, respectively.

Additionally, it will be demonstrated that graphical representation in logarithmic scale for ${\kappa }_p$ and ${\kappa }_b$ is more advantageous, as already used in [22,11]. For the analysis in this paper, the full range does not have to be covered; thus, the following representation is proposed:

$${\kappa }_p=0.03\times {10}^{\left(in{d}_p-1\right)0.05},{\kappa }_b=0.04\times {10}^{\left(in{d}_b-1\right)0.05}$$

(14)

Which, for the ratios from Equation (13), defines a range for the difference between the indexes $in{d}_p$ and $in{d}_b$:

$${\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}=0.75\times {10}^{\left(in{d}_p-in{d}_b\right)0.05},-34.05\le \left(in{d}_p-in{d}_b\right)\le 44.55$$

(15)

3. Methods

3.1. Terminology

This subsection clarifies the adopted terminology. The moving force problem is defined as a longitudinally homogeneous supporting structure being traversed by a constant force at constant velocity. In such a case, the transient part of the solution is not very important, and if a steady-state solution exists, it is rapidly attained.

In this context: the critical velocity (CV) designates a velocity fulfilling four conditions:

(i) it marks a true resonance, meaning that when the force moves at this velocity, deflections tend to infinity in the undamped case (i.e., the steady-state solution is never reached);

(ii) this velocity marks a local minimum of the generalized mode number of an equivalent finite beam;

(iii) there is a sudden jump to zero deflection at the active point (AP) at a velocity infinitesimally higher; and

(iv) the CV value separates regions of instability in regular cases.

On the other hand, the false critical velocity (FCV), which is also a true resonance, marks a local maximum of the generalized mode number of an equivalent finite beam and has no influence on instability regions or the specific value at the AP.

In the three-layer model—as the three layers indicate—three CVs are expected, which means five resonances in the moving force problem: three corresponding to CVs and two to FCVs. In a regular case, there are five vertical asymptotes to the instability lines, marking two full branches in the first two regions delimited by the critical velocities; after that, the final branch tends to zero moving mass for infinite velocity [22,11].

Nevertheless, there are several exceptional cases. In fact, the number of resonances can be 1, 3, or 5. When it is less than 5, some FCVs do not exist, and CVs are replaced by pseudo-critical velocities (PCVs). Additionally, some cases at the boundary of separation regions may lead to a model exhibiting only three resonances, meaning two CVs and one FCV, or one CV and one PCV.

As already mentioned, when CVs are missing, they are replaced by PCVs. A PCV marks a peak in vibrations, but not an infinite value. Since no resonance occurs, even in the absence of damping, deflections never reach infinity. Unlike CVs, FCVs are never replaced, but they can be absent.

Determination of resonances can be done semianalytically as described in [22], but high precision must be implemented. It is extremely important to perform such determinations in the absence of damping, which can unequivocally distinguish a vibration peak from a true resonance. Many methods used by other authors—such as the Green’s function method—require damping for numerical stability and are therefore unsuitable.

PCVs can be determined by parametric analyses. However, when the peak is hardly visible, the peak at the AP, the maximum, and the minimum over the whole beam are usually attained for different velocities, and thus its value is ambiguous. Numerical precision is also essential here. Even if the calculation must be done without damping, MATLAB can be used because the roots required for calculating the steady-state shape can be separated into those defining propagating and evanescent waves. Evanescent waves can be analyzed in the vicinity of the AP. For propagating waves, it is possible to establish an analytical envelope function that can unequivocally identify the maximum and minimum over the entire beam without enlarging the tested region around the AP.

Namely, as our interest lies in the initial region of $\alpha$, then when ${\eta }_s=0$, there are two pairs of real roots and four complex roots symmetrically covering the four quadrants. Negligible damping can be introduced to correctly separate the real roots for the shape ahead of and behind AP. The four real roots define two harmonic waves with constant amplitude and zero phase. It is then possible to calculate the envelope function, which has clearly defined maximum and minimum values over the entire beam. It is then only necessary to check whether, in the vicinity of the AP, these values are exceeded due to evanescent waves. When ${\eta }_s\ne 0$, there is only one harmonic wave, so the amplitude is easily determined and an envelope function is unnecessary.

In the previous text, the term instability line was used. This relates to the instability of moving inertial objects. In such an analysis, the transient part of the solution is essential; without it, the problem of instability would remain hidden. Instability is not a resonance-type phenomenon; there is a velocity interval within which the object becomes unstable. At the extremities of this interval, the object remains stable, and instability is verified only within. This is characterized by an exponential increase in vibration amplitudes. Thus, although the steady-state solution usually exists and is finite, it is never reached.

Instability is determined by solving the roots of the characteristic equation. The roots always come in pairs, identifying the harmonic parts of the full solution. When this part grows exponentially, the roots are termed unstable. The instability line is traced in the ( $\alpha$– ${\eta }_M$) plane as a boundary between regions with a different number of unstable roots. In this paper, only one moving mass is considered. For this case, for a selected $\alpha$, it is possible to track real-valued frequencies until a real-valued equivalent stiffness (or flexibility) of the supporting structure is detected; then, ${\eta }_M$ is calculated from the characteristic equation. In this way, the analysis remains in the real domain, and the D-decomposition method used by other researchers can be avoided.

3.2. Determination of Number of Resonances and Pseudo-critical Velocities

Having the admissible range of dimensionless parameters, it is not easy to analytically indicate the regions corresponding to a particular number of resonances. Even though the base of the semianalytical solution is polynomial and all denominators are removed, the final degree is 16, and several roots are redundant. However, it is relatively easy to calculate the number of resonances semianalytically for a particular case using software with adaptable precision (e.g., Maple or Mathematica). MATLAB is not convenient in this context, as damping should not be included and high precision in root determination is required.

The number of resonances is always odd—1, 3, or 5—corresponding to 1, 2, or 3 critical velocities (CV) and 0, 1, or 2 false critical velocities (FCV), as explained in the previous section.

Four selected cases are shown as a demonstration in Figure 3. It was found convenient to select ${\mu }_s$, ${\mu }_b$ and ${\eta }_s$, or alternatively ${\eta }_N$ and to use the logarithmic scale for ${\kappa }_p$ and ${\kappa }_b$, as indicated in Equation (14). This means that both indexes $in{d}_p$ and $in{d}_b$ were selected for the semianalytical determination from discrete values ranging from 1 to 101 in increments of 1.

As already mentioned, these graphs can only be obtained semianalytically using Maple or similar software because MATLAB does not provide sufficient accuracy in root determination, which is essential.

Once the number of expected resonances is known, the pseudo-critical velocity (PCV) must be determined, as it is often lower than the first resonance. In such cases, the PCV plays the role of a critical velocity for track design. Specifically, in cases with fewer than 5 resonances, when the PCV is lower than the first resonance, it is identified as the velocity that should not be exceeded. Hence, it can be termed the design critical velocity (DCV).

PCVs, which can be determined by parametric analysis, may be well-pronounced or ambiguous in the sense that the peaks do not always occur at the same position. The following figures explain these issues graphically. The case from Figure 3(d) is selected first. Then, in addition to ${\mu }_s=4$, ${\mu }_b=12$ and ${\eta }_s=1$ (${\eta }_N=0$) and $in{d}_b=21$, $in{d}_p=31;41;51$ are used to identify cases with 1, 3, and 5 resonances..

In Figure 4, the case with 1 resonance is shown. This is an exceptional case characterized by a maximum of 3 resonances. Since only one resonance is present, the FCV is absent. The PCV is relatively well-pronounced at $\alpha =0.310$ for both the maximum and minimum over the whole beam, while the peak at AP occurs at $\alpha =0.303$. These values are determined with a precision of 0.001, which corresponds to the step size used in the parametric analysis. Resonance occurs at $\alpha =1.016877$.

In Figure 5, the case with 3 resonances is shown. Two CVs are clearly marked by a sudden jump to zero at the AP: $\alpha =0.317124$ and $\alpha =1.354864$. Then, $\alpha =0.327508$ is the FCV. There is also a barely visible PCV with peaks at $\alpha =0.3792$ for the maximum and minimum over the whole beam, while the peak at AP occurs at $\alpha =0.3706$. These values are determined with a precision of 0.0001, which is the step size used in the parametric analysis.

Finally, in Figure 6, the case with 5 resonances is shown. Three CVs are clearly marked by a sudden jump to zero at the AP: $\alpha =0.319163$, $\alpha =0.400814$ and $\alpha =1.806272$. Between them, two FCVs occur: $\alpha =0.340119$ and $\alpha =0.443505$. There are no indications of other vibration peaks.

Several other cases are included in [22,11].

Due to the challenges in determining the PCV—first, the necessity of parametric analysis instead of semianalytical methods, and second, the ambiguity in identifying clear peaks—it is worthwhile to explore other possibilities. In Figure 5, at least, the maxima and minima over the beam occurred at the same $\alpha$, but this is not always the case. Therefore, the PCV can be ambiguous when the peaks do not coincide.

Since the PCV replaces the CV in separating instability regions, it is natural to propose that it should be predicted by $\alpha$, which indicates a vertical asymptote in the instability line. This newly proposed approach exploits the regularity of the single moving mass problem, exemplified in the one-layer model in Figure 7.

It is demonstrated in Figure 7 that, as the damping ratio decreases, the instability lines approach the CV. Therefore, a very good estimate of the CV can be obtained by identifying the onset of instability for infinitely low damping. Such calculations can be performed semianalytically in Maple or similar software because a high number of digits of precision is required. Infinitely low damping and infinitesimally small real-valued frequency are used to find an $\alpha$-value leading to a real-valued equivalent stiffness (or flexibility). The corresponding ${\eta }_M$ tends to infinity.

Applying the same reasoning to the three-layer model, it should be possible to identify the lowest PCV. This was tested for the case shown in Figure 3(a) with ${\mu }_s=1$, ${\mu }_b=2$ and ${\eta }_s=0$. In Figure 8(a), $in{d}_b=26$, and in Figure 8(b), $in{d}_b=29$. The latter case is the only one where the region of 3 resonances extends fully, and 5 resonances occur only at the very beginning.

It can be concluded from Figure 8 that PCV identification works well: the line of green crosses is smooth, as expected. In regions with 5 resonances, the crosses coincide, while in other regions, green crosses indicate a much lower value. Nevertheless, there are cases—as already mentioned—where the maximum number of resonances is only 3. Then, the crosses match in the region of 3 resonances, as shown in Figure 9, where the results of the parametric analysis are also included to show that no additional peaks occur (higher velocities are not shown but were tested). The graph in part (a) is plotted only up to $in{d}_p=34$ for clarity.

In Figure 9(b), the resonances are: 0.339502; 0.377294 and 0.816655, corresponding to CV, FCV, and CV, respectively. The lowest asymptote is at $\alpha =0.339502$, which perfectly matches the lowest resonance. This shows that the proposed technique works well in most cases. However, there are exceptions where the lowest PCV is lower than the lowest asymptote. These cases are, in fact, the focus of this paper.

This section concludes by illustrating the effect of PCV, CV, and FCV on instability lines in typical cases. When there are 5 resonances, three of them correspond to CVs and two to FCVs. Five vertical asymptotes to the instability lines then appear, defining two full branches in the first two regions delimited by the CVs. Beyond that, the final branch tends to zero ${\eta }_M$ as $\alpha$ approaches infinity. When there are fewer resonances, CVs may be replaced by PCVs. In contrast, FCVs are never replaced and may be absent, as already mentioned.

The same cases as previously shown in [11] are selected for demonstration. These cases differ by ${\kappa }_p$, highlighting how this parameter affects the number of resonances. Specifically, both cases have ${\mu }_s=6$, ${\mu }_b=35$, ${\kappa }_b=7$ and ${\eta }_s=0$. For ${\kappa }_p=300$, there are 5 resonances (0.151–CV1; 0.152–FCV; 0.627–CV2; 0.851–FCV; 4.244–CV3). However, for ${\kappa }_p=30$, there is only one (~0.149–PCV1; 0.599–PCV2; 2.388–CV3). This is demonstrated in Figure 10.

It is clear from Figure 10 that the FCV does not influence the separation of instability lines into three distinct regions. It is further shown that PCVs act as CVs, and with decreasing damping, instability lines approach either the CV or PCV. In the middle region, larger deviations are observed from the right than from the left. To complete this analysis, the corresponding real-valued frequency $q_M$, which gives rise to Figure 10 as a root of the characteristic equation, is plotted in Figure 11.

Figure 11 again demonstrates that as damping decreases, CV and PCV are approached by vertical asymptotes for which $q_M$ tends to zero and ${\eta }_M$ tends to infinity. However, part (b) shows that in the first region, the damping level is not sufficiently low to clearly illustrate this behavior.

An interesting situation occurs when the above-described technique fails and such an asymptote does not exist. In these cases, only a slight increase in vibration is observed at the PCV. Alternatively, an asymptote may exist, but the instability line covers only unreasonably high ${\eta }_M$ values, so that instability cannot occur. In such cases, circulation velocities can exceed the expected critical velocity.

3.3. Motivation and Aim

As already mentioned, the motivation for this analysis lies in the observation that there are cases where the DCV is a PCV, which is not identified by a vertical asymptote. Thus, it appears possible to tailor the parameters in a way that reduces such vibrations, thereby extending the operational range into the supercritical region.

This suggests that it may be feasible to adjust the parameters of the model using optimization techniques so that the increase in vibrations remains minimal. Consequently, the DCV could be shifted to higher values without compromising safety—meaning, without the risk of excessive vibrations or potential instability. In more detail, the parameters must be adjusted such that the peak at the PCV is minimized, while simultaneously ensuring that the first resonance and the onset of instability for a reasonable value of ${\eta }_M$ and a moderate damping level occur at the highest possible velocity.

Mathematically, the first criterion is defined as the ratio between the highest upward displacement at the DCV = PCV, denoted $\left|{\tilde{w}}_{\min }^{PCV}\right|$, and the maximum downward displacement in the static case, ${\tilde{w}}_{st}$. The reasoning is twofold: first, because dimensionless displacements are scaled with respect to the static displacement of a reference beam, it would not be reasonable to use unity for comparison; second, because upward displacements are more critical, they were chosen as the primary criterion.

This criterion is combined with two others: the lowest resonance ${\alpha }_{res}$, and the velocity at the onset of instability ${\alpha }_{m50}$ for a representative mass ${\eta }_M=50$ under 1% damping ratio at all levels, ${\eta }_p={\eta }_b={\eta }_f=0.01$. Based on these, the performance index (PI) is defined as:

$$\mathrm{PI}={\displaystyle \frac{\left|{\tilde{w}}_{\min }^{PCV}\right|}{{\tilde{w}}_{st}}}+{\displaystyle \frac{1}{100\min \left(\frac{{\alpha }_{res}-{\alpha }_{PCV}}{{\alpha }_{PCV}},\frac{{\alpha }_{m50}-{\alpha }_{PCV}}{{\alpha }_{PCV}}\right)}}$$

(16)

This PI also serves as the objective function, with the aim of finding its minimum within the admissible domain.

Several optimization methods were tested, beginning with statistical analyses—in particular, design of experiments (DOE) and response surface methodology—followed by parametric search (PS), the Monte Carlo method (MC), and Simulated Annealing (SA). It was also examined whether a full analysis was necessary or if some pre-selection could be applied. Certain parts required computation in Maple due to the need for increased numerical precision, while other parts could be handled in MATLAB.

4. Particular Case Analysis

4.1. Preliminary Tests

In this section, a particular case is tested to check the viability of the PI definition given by Equation (16). It is obvious that such cases must have a sufficiently large region of $in{d}_p$, where the difference between PCV and ${\alpha }_{res}$ or ${\alpha }_{m50}$ is significant; otherwise, no improvement can be obtained. For this selection, the analysis of the number of resonances is useful, because these cases must have a reduced number of resonances. The selected case is characterized by ${\mu }_s=3$, ${\mu }_b=10$ and ${\eta }_s=0$. Two indexes are chosen to characterize ${\kappa }_b$, namely $in{d}_b=20$ and $in{d}_b=63$. $in{d}_p$ is tested for discrete values $in{d}_p=1:1:51$. Results are presented in Figure 12.

Several curves are plotted in Figure 12. In the legend: “peak” means $\alpha$ at the first peak in upward displacement, identifying the lowest PCV determined by parametric analysis in MATLAB.

It can be concluded from Figure 12(b) that the PCV line is not realistic, because it should be smooth. Very accurate calculation is used by implementing function envelopes for the propagating part of the solution; thus, this decrease is not a numerical error but an indication of a premature peak, which should be disregarded for PCV determination. With regularization, the graph in Figure 13 is obtained.

From Figure 12 and Figure 13, it can be concluded that there is a sufficiently large region of $in{d}_p$, where the difference between PCV and ${\alpha }_{res}$ or ${\alpha }_{m50}$ is quite large, confirming the feasibility of the adopted approach. Calculating PI for these cases, values of 1.3854 at $in{d}_p=19$ and 0.5412 at $in{d}_p=18$ are obtained for $in{d}_b=20$ and $in{d}_b=63$, respectively.

However, it is not sufficient to calculate PI alone, because the feasibility of the data used must be confirmed. First, Equations (12) and (15) must be fulfilled. Unfortunately, this is still not sufficient.

4.2. Real Data Identification

As already mentioned, even when the feasibility of the main dimensionless parameters is confirmed by Equations (12) and (15), this does not ensure the existence of a corresponding real data set within the intervals from Equation (11), nor within other ranges for $m_s$, $k_p$, ${\rho }_b$ and $E_b$, as indicated in Section 2 (text below Figure 2). It is also necessary to consider realistic values for $m$ ($m=54\mathrm{kg}/\mathrm{m}$ or $m=60\mathrm{kg}/\mathrm{m}$) and $k_f\in 〈0.22;1000〉{\mathrm{MN}/\mathrm{m}}^2$, [11], which complete the full set of parameters defining the intervals indicated in Equation (10).

The selection of real parameters can be conveniently carried out in the following way. First, the stress cone is parameterized according to:

$$l_s=0.5:0.1:0.8,h_b=0.2:0.05:0.6,l_b=0.2:0.02:0.36,{\alpha }_b={25}^{\mathrm{o}}:5^{\mathrm{o}}:{50}^{\mathrm{o}}$$

(17)

For each set of parameters, the dynamically activated ballast volume $m_{cn}$ is calculated in [m³], and the parameter related to the ballast layer vertical stiffness $k_{cn}$ is calculated in [m]. Knowing the optimal, i.e., target values ${\mu }_s$, ${\mu }_b$, ${\kappa }_p$ and ${\kappa }_b$, it must hold:

$$\begin{gathered} {\displaystyle \frac{m_s^{-}}{{\rho }_b^{+}}}={\displaystyle \frac{40}{2600}}={\displaystyle \frac{1}{65}}\le \left(m_{cn}{\displaystyle \frac{{\mu }_s}{{\mu }_b}}={\displaystyle \frac{m_s}{{\rho }_b}}\right)\le {\displaystyle \frac{m_s^{+}}{{\rho }_b^{-}}}={\displaystyle \frac{160}{1200}}={\displaystyle \frac{1}{7.5}}, \\ {\displaystyle \frac{k_p^{-}}{E_b^{+}}}={\displaystyle \frac{20}{400}}=0.05\le \left(k_{cn}{\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}={\displaystyle \frac{k_p}{E_b}}\right)\le {\displaystyle \frac{k_p^{+}}{E_b^{-}}}={\displaystyle \frac{5000}{50}}=100 \end{gathered}$$

(18)

Then, if these conditions are fulfilled, the selection may continue by defining admissible values:

$$\begin{gathered} {\widehat{m}}_s^{-}=\max \left(m_s^{-},{\rho }_b^{-}{m}_{cn}{\displaystyle \frac{{\mu }_s}{{\mu }_b}}\right)=\max \left(40,1200m_{cn}{\displaystyle \frac{{\mu }_s}{{\mu }_b}}\right), \\ {\widehat{m}}_s^{+}=\min \left(m_s^{+},{\rho }_b^{+}{m}_{cn}{\displaystyle \frac{{\mu }_s}{{\mu }_b}}\right)=\min \left(160,2600m_{cn}{\displaystyle \frac{{\mu }_s}{{\mu }_b}}\right), \end{gathered}$$

(19)

$$\begin{gathered} {\widehat{k}}_p^{-}=\max \left(k_p^{-},E_b^{-}{k}_{cn}{\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}\right)=\max \left(20,50k_{cn}{\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}\right), \\ {\widehat{k}}_p^{+}=\min \left(k_p^{+},E_b^{+}{k}_{cn}{\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}\right)=\min \left(5000,400k_{cn}{\displaystyle \frac{{\kappa }_p}{{\kappa }_b}}\right), \end{gathered}$$

(20)

Further:

$${\widehat{m}}^{-}={\displaystyle \frac{{\widehat{m}}_s^{-}}{l_s{\mu }_s}},{\widehat{m}}^{+}={\displaystyle \frac{{\widehat{m}}_s^{+}}{l_s{\mu }_s}},{\widehat{k}}_f^{-}={\displaystyle \frac{{\widehat{k}}_p^{-}}{l_s{\kappa }_p}},{\widehat{k}}_f^{+}={\displaystyle \frac{{\widehat{k}}_p^{+}}{l_s{\kappa }_p}},$$

(21)

Finally, if the interval $〈{\widehat{m}}^{-},{\widehat{m}}^{+}〉$ has some intersection with $〈m^{-},m^{+}〉=〈54,60〉\mathrm{kg}/\mathrm{m}$, and also the interval $〈{\widehat{k}}_f^{-},{\widehat{k}}_f^{+}〉$ has some intersection with $〈k_f^{-},k_f^{+}〉=〈0.22;1000〉{\mathrm{MN}/\mathrm{m}}^2$, then feasible design values can be found. As only two rail profiles are commonly used, either 54 or 60 kg/m is selected. For each of these profiles $EI$ is known. It was further confirmed during other optimization techniques that the optimal $k_f$ is usually very high; thus, in this procedure, it is more convenient to take the lowest feasible value. This completes the definition of the reference beam, and the missing values are simply calculated.

All results are summarized in tables, which are conveniently placed in the Appendix so that the results of all methods can be consulted together. In Table 1, the main optimized data are given in the form of optimized dimensionless parameters, PI, PCV, and an indication of possible extension, both by absolute $\alpha$-value and by percentage increase with respect to PCV. Table 2 summarizes the parameters of the stress cone; Table 3 lists the other values used to match the optimized data.

4.3. Design of Experiments and Response Surface

DOE and the response surface were the first optimization approaches used. Unfortunately, the response surface did not yield any improvement. When a larger design domain was applied, several cases with a set of fixed parameters ${\mu }_s$, ${\mu }_b$, $in{d}_b$ and ${\eta }_s$ produced undefined PI across the full range of $in{d}_p$. Consequently, it was necessary to use a penalizing value, which distorted the response surface. Even on a more restricted design domain encompassing the particular case from Section 4.1, a reasonable estimate of already known values was never obtained. A Monte Carlo search over the response surface indicated the best choice at one of the factorial design points. Residues were not well-behaved, and thus this type of analysis had to be disregarded from further steps.

The best value obtained from this analysis, after removing non-feasible cases, was PI=0.4956, achieved for ${\mu }_s=4$, ${\mu }_b=12$, $in{d}_b=40$, ${\eta }_s=0$ and $in{d}_p=16$. In this case, the full factorial design was performed for ${\mu }_s^{-}=2$, ${\mu }_s^{+}=4$, ${\mu }_b^{-}=8$, ${\mu }_b^{+}=12$, $in{d}_b^{-}=40$, $in{d}_b^{+}=80$, ${\eta }_s^{-}=0$, and ${\eta }_s^{+}=1$, and the best feasible choice corresponded to the combination ++--. Results of this analysis are presented in Figure 14, which shows that PCV=DCV can be increased by 44.05%. Numerical values are included in Table 1 (Appendix).

In Figure 14, the yellow dashed line marks the first peak in upward displacement indication in this way PCV=DCV. The second criterium is governed by ${\alpha }_{res}$, which is not marked by a vertical line because it is identified through the properties of CV—that is, a sudden jump to zero displacement at AP. The corresponding probability plot is shown in Figure 15.

The best design obtained during DOE (line 1 in Table 1) offers many possibilities for selecting a corresponding set of real data. As it is not feasible to present all of them, only two representative cases are selected, and the corresponding values are provided in Table 2 and Table 3.

5. Optimization and Parametric Search

A full analysis, as reported in Figure 12 and Figure 13 and performed in DOE, is quite time-consuming and therefore not suitable for use in metaheuristic optimization methods or parametric search. Consequently, these methods were applied only to a preselection of favorable cases, after which the best case was further analyzed over $in{d}_p$. Preselection was implemented in two ways: based on cases with the highest differences in ${\alpha }_{as}^{0.01}$ or in ${\alpha }_{res}$ over the range of $in{d}_p$. ${\alpha }_{as}^{0.01}$ can be easily extracted in MATLAB; however, lower damping levels would be prone to numerical errors due to MATLAB’s double precision. ${\alpha }_{res}$ had to be calculated in Maple because high precision in root determination is required. In what follows, $d{a}_{0.01}^{\%}$ will be used for percentage difference in ${\alpha }_{as}^{0.01}$ over the range of $in{d}_p$ and $d{r}^{\%}$ for the percentage difference in ${\alpha }_{res}$.

5.1. Monte Carlo Method

MC method was performed with preselection of the most favorable cases according to $d{a}_{0.01}^{\%}$. The design domain was defined as ${\mu }_s^{-}=1$, ${\mu }_s^{+}=6$, ${\mu }_b^{-}=2$, ${\mu }_b^{+}=45$, $in{d}_b^{-}=1$, $in{d}_b^{+}=101$, ${\eta }_s^{-}=0$ and ${\eta }_s^{+}=1$.

It was not useful to repeat the analysis, because the MATLAB function rand tends to reuse the same randomization, and several tested cases were duplicated. The best value obtained was 88.89%, achieved for ${\mu }_s=5.5244$, ${\mu }_b=44.1292$, $in{d}_b=44.8870$ and ${\eta }_s=0.1111$. This case was then analyzed for PI. Results of the analysis over $in{d}_p$ are presented in Figure 16. Asymptote and instability were only tested in the relevant region where PCV does not coincide with ${\alpha }_{res}$.

$in{d}_p$ at minimum value of PI was gradually increased util at least one feasible design appeared. This happened for $in{d}_p=11$ where PI=0.1905. Numerical data obtained are summarized in Table 1, Table 2 and Table 3. Final results are shown in Figure 17. ${\alpha }_{res}$ is again decisive for PI, and the possible increase in DCV is 205.76%.

5.2. Parametric Search

PS was also performed using preselection of the most favorable cases. Both $d{a}_{0.01}^{\%}$- and $d{r}^{\%}$-based options were used. Results of this preselection are shown in Figure 18 and Figure 19.

After preselection, the four best cases were analyzed for PI over $in{d}_p$. $in{d}_p$ at the minimum PI was gradually increased util at least one feasible design appeared. Final results are shown in Figure 20, Figure 21, Figure 22 and Figure 23. ${\alpha }_{res}$ is again decisive for PI.

PI values in these four cases are 0.4257; 0.4115; 0.2517and 0.2462, leading to a possible increase in DCV of 55.60%; 65.64%; 132.22% and 182.10%, respectively. Numerical data are summarized in Table 1, Table 2 and Table 3.

5.3. Simulated Annealing

SA was selected as another metaheuristic optimization method due to its ease of implementation in MATLAB. The predefined parameters governing the optimization process were kept unchanged. The design domain was defined as ${\mu }_s^{-}=1$, ${\mu }_s^{+}=6$, ${\mu }_b^{-}=2$, ${\mu }_b^{+}=45$, $in{d}_b^{-}=1$, $in{d}_b^{+}=50$, ${\eta }_s^{-}=0$, and ${\eta }_s^{+}=2$, and used in this order. Several initial data sets were used for preselection according to $1/d{a}_{0.01}^{\%}$, which was searched as a minimum.

Figure 24 shows the results of different runs.

The best value obtained is 1.46401, corresponding to 89.53%, which is better than that obtained by the MC method (88.89%) but worse than from the PS (89.90%). The preselected parameters ${\mu }_s=5.3333$, ${\mu }_b=42.8571$, $in{d}_b=37.3333$, and ${\eta }_s=0$ were then used for detailed analysis of PI. As before, $in{d}_p$ at the minimum value of PI was gradually increased until at least one feasible design appeared, which occurred for $in{d}_p=12$, for which PI=0.2055. The results of analysis over $in{d}_p$ are presented in Figure 25. Asymptote and instability were again tested only in the relevant region where PCV does not coincide with ${\alpha }_{res}$.

Final results are shown in Figure 26. ${\alpha }_{res}$ is again decisive. DCV=0.2437 could allow for an additional increase of 132.22%.

6. Discussion

In the previous sections, several feasible designs for which PI was evaluated as optimal have been presented. It is worthwhile to remark that ${\alpha }_{m50}$ was never decisive in these cases, and the extension of DCV was always guided by ${\alpha }_{res}$. From Table 1, it is evident that ${\mu }_s$ and ${\mu }_b$ exhibit higher dispersion, while $in{d}_p$ and $in{d}_b$ have very similar values in all cases. Finally, ${\eta }_s=0$ appears almost exclusively, except in the MC method—likely because a zero value was improbable to be selected randomly along with other optimal parameters. Furthermore, higher ${\mu }_s$ generally provides larger improvements, and $in{d}_p$ is always much lower than $in{d}_b$.

$in{d}_p$ at minimum value of PI was gradually increased util at least one feasible design appeared in all cases, except for DOE, where, in fact, no optimization was presented. Since these designs are the first to appear, this implies that some data will lie at the limit values of the considered ranges. As such, they might seem unrealistic, but it must be noted that further increase of $in{d}_p$ can lead to better selections without significantly compromising the PI value. Generally, for higher $in{d}_p$, the possible increase up to ${\alpha }_{res}$ is more favorable; however, the ratio $\left|{\tilde{w}}_{\min }^{PCV}\right|/{\tilde{w}}_{st}$ is increasing.

Parametric searches with a base value of ${\mu }_s=1$ can naturally only identify designs where this is possible—meaning designs using wooden sleepers. Nevertheless, this option was not disregarded, because it is essential to consider the possibility of using novel lightweight materials in the future.

Some of the real data from Table 2 and Table 3 may seem unreasonable, but this results from the initial ranges of possible values.

Due to the consistency of results obtained by different approaches, exploring other methods would not be useful.

The results confirm that preselection is useful and reduces calculation time, but the connections are not straightforward. It is not necessarily true that the best preselection yields the best PI and, at the same time, the highest possible increase in DCV.

When using the lowest feasible $in{d}_p$, combinations with relatively high $in{d}_b$ led to the appearance of the lowest $k_p$ found in the literature in all feasible data shown in Table 3. Similarly, for the same reason, $k_f$ is relatively high, resulting in quite high DCV—difficult to attain with current train designs. This implies that although the possible increase is large, it may be practically unimportant. Nevertheless, these results may inspire new approaches and provide guidance toward optimizing track design.

Further analyses are needed to verify whether PI was defined appropriately or if it should be modified. Since instability was never found to be decisive, it can likely be excluded from this analysis, allowing more detailed investigations.

Additional tests are also necessary because although several authors use this model without shear resistance, it must be verified whether such modeling is truly valid in practical applications.

7. Conclusions

In this paper, tailoring of input geometric and material data of the railway track was presented, with the objective of identifying cases where the first increase in vibrations is negligible and can therefore be bypassed.

Several optimization techniques and approaches were employed, leading to similar conclusions. However, the results require further examination in several aspects. In particular, the adequacy of the proposed PI must be evaluated. Since all results showed very high DCV with negligible increase in vibration, the significance of a further increase diminishes. It is thus necessary to investigate whether similar outcomes could be achieved for lower real velocities, where the potential for meaningful improvements is greater.

Results were obtained semianalytically with very high precision. This was crucial for accurately determining the number of resonances and identifying the first peak in upward vibration.

It was found that DCV can be significantly increased by optimizing the design parameters. However, such cases almost exclusively result in designs where the DCV is already very high. Therefore, the practical significance of further increase becomes questionable, as already mentioned above.

In all cases, very soft rail pads emerged in the final optimized design. It is essential to verify whether such configurations are feasible. Thus, the results suggest a tendency toward:

• Very high foundation stiffness,
• Low rail pad stiffness,
• Zero ballast shear stiffness.

Further research is particularly needed regarding ballast shear stiffness. In actual designs, eliminating it entirely may not be realistic. While it might be reduced via rail pretension, such an approach is not practical.

Due to the scope of this study, these aspects will be addressed in future research. It is believed that the theoretical findings presented here can contribute to the development of optimal track designs incorporating novel materials. Their implementation in real-world designs remains a subject for further investigation.