Locomotion Dynamics of a Vibration-Driven Robot with Spring–Damper Coupling via Non-Newtonian Fluids

1. Introduction

Vibration-driven robots are a class of locomotion mechanisms that utilize oscillatory motion to propel themselves across different surfaces [1,2]. Unlike conventional mobile robots that rely on wheels, legs, or tracks, vibration-driven robots exploit oscillation-induced frictional forces to generate net forward motion [3,4]. This principle eliminates the need for complex mechanical transmission systems and reduces structural redundancies, thereby making such robots simple to fabricate, lightweight, and cost-efficient [5,6]. Furthermore, vibration-driven robots can traverse diverse environments, including flat or inclined planes, rough terrains, and irregular surfaces, which highlights their versatility in both engineered and natural settings [8]. These attributes have positioned vibration-driven robots as attractive solutions for a range of tasks, such as post-disaster debris clearance, in-pipe inspection and maintenance, and even targeted medical applications at small scales [9,10,11].

A fundamental feature of vibration-driven robots is their reliance on internal oscillatory components, such as unbalanced rotating or reciprocating masses, piezoelectric actuators, pneumatic chambers, or electromagnetic drivers, to generate propulsion [12,13]. The locomotion efficiency of such robots strongly depends on how the oscillatory input energy is distributed and transformed into translational displacement through interactions with the ground [14]. An essential factor in this process is the coupling mechanism between different modules of the robot, which governs the exchange of energy and the resultant motion profile [15]. Most prior designs have relied on purely elastic coupling, typically springs, to connect modules. While springs can store and release oscillatory energy, they do not account for dissipative effects that naturally occur in real-world systems. This leads to simplified dynamics that may not capture the complex interplay between actuation, damping, and friction [16,17].

Previous research on vibration-driven robots has largely focused on simpler linear models. For instance, [11] analyzed a two-module vibration-driven system with elastic coupling and investigated the influence of friction asymmetry on net displacement. [7] derived governing dynamic equations for a two-body system connected by a spring and driven by unbalanced rotating masses, and validated their models with experiments. [3] explored dynamic modeling and open-loop control for two- and three-module vibration robots actuated by reciprocating internal masses, showing how excitation phase shifts affect steady-state velocity. [24] further analyzed locomotion due to stick–slip effects in two-module systems. While these studies provide valuable insight into the dynamics of vibration-driven systems, they all assume linear elastic coupling and neglect complex damping phenomena.

In order to more accurately represent and exploit such dynamics, damping mechanisms have been incorporated into various robotic systems. Conventional approaches often rely on linear viscous damping with constant coefficients [18,19,20]. However, in practice, damping is rarely linear, and materials or fluids often display nonlinear or rate-dependent behaviors. In particular, non-Newtonian fluids have emerged as a promising candidate for advanced damping applications [21]. Unlike Newtonian fluids, whose viscosity remains constant under varying shear rates, non-Newtonian fluids exhibit shear-dependent viscosity, which can be either shear-thinning (decreasing viscosity with increasing shear rate) or shear-thickening (increasing viscosity with shear rate) [22]. These unique rheological properties allow the design of damping elements whose resistance adapts dynamically to the excitation conditions [23].

Integrating non-Newtonian fluids into vibration-driven robots introduces a new paradigm in locomotion. When such a fluid is embedded in a damper connecting the robot’s modules, the effective damping becomes a nonlinear function of the excitation frequency and velocity. This provides a natural form of adaptive dissipation that can enhance stability, improve energy efficiency, or even generate novel locomotion regimes. Importantly, when the robot is actuated by an unbalanced rotating mass, the angular velocity of the rotor directly influences the shear rate in the non-Newtonian fluid, thereby coupling actuation and damping in a unique manner. As a result, the damping force is not only time-varying but also explicitly dependent on the actuator dynamics, creating a rich nonlinear interaction between oscillatory input, elastic restoring forces, nonlinear damping, and ground friction.

Some works have introduced damping explicitly but in linearized forms. [25] investigated single-module bristle-bots with unbalanced masses, analyzing velocity dependence on excitation and friction coefficients. [15] proposed a capsule-type vibration robot with an electromagnetic actuator and spring coupling, controlled via duty-cycle modulation. [26] presented a one-way wheeled vibration robot actuated by contra-rotating unbalanced masses, where locomotion efficiency was optimized near resonance frequencies. More recently, [5] conducted bi-objective optimization of a two-module vibration robot to simultaneously maximize forward velocity and minimize backward slip, using the NSGA-II algorithm [27]. Despite these advances, damping was always treated as either negligible or constant, leaving the potential of nonlinear, velocity-dependent damping largely unexplored.

The motivation of this work stems from these research gaps. Specifically, there has been limited exploration of how nonlinear damping, especially through non-Newtonian fluids, can influence the performance of vibration-driven robots. Given the strong velocity dependence of such fluids, incorporating them into spring–damper couplings offers the possibility of novel locomotion behaviors and enhanced energy efficiency. Moreover, when the damping characteristics are directly coupled to actuator velocity, the optimization and control strategies must be revisited, as linear assumptions no longer hold.

In this paper, we propose and analyze a two-module vibration-driven robot in which the modules are connected by a spring–damper system filled with a non-Newtonian fluid (as shown in Figure 1). The damper introduces nonlinear, velocity-dependent resistance that varies with the angular velocity of the unbalanced rotating mass. The governing dynamic equations are derived to capture this coupled nonlinear behavior, and the model is validated against numerical simulations. A parametric study is then conducted to assess the influence of key design variables, such as module mass, spring stiffness, and fluid properties, on the robot’s average velocity and energy consumption. Multi-objective optimization using NSGA-II is performed to determine parameter sets that maximize locomotion efficiency while minimizing energy use.

The main contributions of this work can be summarized as follows. First, we introduce a vibration-driven robot with a spring–damper coupling that incorporates a non-Newtonian fluid, resulting in nonlinear and velocity-dependent damping directly linked to the actuator dynamics. Second, we develop and validate governing dynamic equations that capture the interplay between elastic restoring forces, nonlinear damping, and frictional interactions with the ground. Third, we perform a multi-objective optimization of the robot parameters to simultaneously enhance locomotion velocity and reduce energy consumption.

This study advances the field of vibration-driven robotics by incorporating realistic nonlinear damping into the modeling and optimization framework, and by demonstrating that adaptive damping via non-Newtonian fluids can significantly influence locomotion performance and control.

2. Nonlinear Governing Equations of the Spring–Damper Coupled Robot

In this study, we investigate a vibration-driven modular robot composed of two rigid modules, denoted $M_1$ and $M_2$, connected via a spring–damper element. The coupling mechanism consists of a linear elastic spring of stiffness k and initial rest length d, in parallel with a nonlinear damper containing a non-Newtonian fluid. The damping coefficient is a function of the angular velocity $\omega$ of an unbalanced mass m mounted on a rod of length L pivoted at the center of the front module. The angular position, velocity, and acceleration of this unbalanced mass are $\varphi$, $\omega$, and $\alpha$, respectively. A schematic of the system is shown in Figure 2.

2.1. Assumptions and Modeling Considerations

To formulate the equations of motion, the following assumptions are made:

1.

The rod carrying the unbalanced mass is considered massless, so its rotational inertia is negligible.

2.

The modules are rigid, with equal masses $M_1=M_2=M$.

3.

The robot’s motion is constrained along a single horizontal axis, corresponding to translational locomotion.

4.

Static and kinetic coefficients of friction are identical, ${\mu }_s={\mu }_k=\mu$, to simplify the frictional model.

5.

The unbalanced mass is treated as a point mass, allowing simplification of the torque and centrifugal force expressions.

These assumptions allow the system to be modeled using a reduced-order description, while retaining the key nonlinear effects arising from the velocity-dependent damping. Except for the damping factor, the values of the parameters are based on [7] and presented in Table 1.

2.2. Horizontal Dynamics of the Front Module

Applying Newton’s second law in the horizontal direction to the front module yields the sum of all forces acting along the motion axis. These include the spring force, nonlinear damping, centrifugal and tangential contributions of the unbalanced mass, and the frictional resistance $R_2$. Denoting the displacements of the front and rear modules as $x_2$ and $x_1$, respectively, the governing equation reads:

$$\begin{array}{cc}\hfill (M+m){\ddot{x}}_2& +k(x_2-x_1-d)+c\left(\omega \right)({\dot{x}}_2-{\dot{x}}_1)\hfill \\ \hfill & =mL{\omega }^{2}cos\left(\varphi \right)+mL\alpha sin\left(\varphi \right)-R_2.\hfill \end{array}$$

(1)

The damping coefficient $c\left(\omega \right)$ accounts for the non-Newtonian nature of the fluid. For a shear-thickening fluid, $c\left(\omega \right)$ increases with $\omega$, whereas for a shear-thinning fluid, it decreases. A commonly used model is:

$$c\left(\omega \right)=c_0\left(1+{\beta \left|\omega \right|}^n\right),$$

(2)

where $c_0$ denotes the baseline viscosity-related parameter, $\beta >0$ is a scaling factor that modulates the degree of rate sensitivity, and $n>0$ defines the non-Newtonian regime, with $n<1$ corresponding to shear-thinning behavior and $n>1$ to shear-thickening behavior.

2.3. Vertical Equilibrium of the Front Module

The vertical forces acting on the front module consist of the module’s weight, the contribution of the unbalanced mass, and the normal reaction $N_2$ from the ground. To ensure continuous contact with the surface:

$$N_2=(M+m)g+mL\alpha cos\left(\varphi \right)-mL{\omega }^{2}sin\left(\varphi \right)\ge 0,$$

(3)

where g is the gravitational acceleration. This condition imposes bounds on the admissible angular velocity and acceleration of the rotor.

2.4. Dynamics of the Rear Module

The rear module experiences horizontal forces due to the spring, the nonlinear damper, and friction $R_1$:

$$M{\ddot{x}}_1+k(x_1-x_2+d)+c\left(\omega \right)({\dot{x}}_1-{\dot{x}}_2)=-R_1.$$

(4)

Its vertical equilibrium is trivial, as it does not carry the unbalanced mass:

$$N_1=Mg.$$

(5)

2.5. Frictional Forces

The frictional forces acting on each module depend on the instantaneous motion state $i=1,2$:

$$\begin{array}{c}\hfill R_i=\left\{\begin{array}{cc}\mathrm{applied}\mathrm{force},\hfill & \left|\mathrm{applied}\mathrm{force}\right|<\mu N_i,\left(\mathrm{static}\right)\hfill \\ \mu N_i\mathrm{sign}\left({\dot{x}}_i\right),\hfill & \mathrm{otherwise},\left(\mathrm{kinetic}\right)\hfill \end{array}\right..\end{array}$$

(6)

This piecewise definition ensures smooth transition from static to kinetic friction in simulations.

2.6. Matrix Form Representation

For compactness and suitability for numerical integration, the equations can be expressed in matrix form:

$$\mathbf{M}\left[\begin{array}{c}{\ddot{x}}_1\\ {\ddot{x}}_2\end{array}\right]+\mathbf{K}\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]+\mathbf{C}\left(\omega \right)\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\mathbf{F},$$

(7)

where

$$\mathbf{M}=\left[\begin{array}{cc}M& 0\\ 0& M+m\end{array}\right],\mathbf{K}=\left[\begin{array}{cc}k& -k\\ -k& k\end{array}\right],$$

(8)

and

$$\mathbf{C}\left(\omega \right)=\left[\begin{array}{cc}c\left(\omega \right)& -c\left(\omega \right)\\ -c\left(\omega \right)& c\left(\omega \right)\end{array}\right],$$

(9)

while the forcing vector is

$$\mathbf{F}=\left[\begin{array}{c}-R_1-kd\\ mL{\omega }^{2}cos\left(\varphi \right)+mL\alpha sin\left(\varphi \right)-R_2+kd\end{array}\right].$$

(10)

2.7. Actuator Power and Energy Consumption

The instantaneous power delivered by the unbalanced rotating mass is

$$P\left(t\right)=T\omega ,$$

(11)

where the driving torque is

$$T=m{L}^2\alpha -mL{\ddot{x}}_{2}sin\left(\varphi \right)+mgLcos\left(\varphi \right).$$

(12)

The total actuation energy consumed over a time interval $[0,t]$ is therefore

$$E={\int }_0^t\left|P\left(\tau \right)\right|d\tau .$$

(13)

2.8. Remarks on Nonlinear Coupling

The velocity-dependent damping term introduces strong nonlinear coupling between the rotational actuation and translational locomotion. Unlike purely elastic couplings, this system exhibits frequency-dependent locomotion effi- ciency, which must be carefully considered in control and design. In subsequent sections, the validity of these govern- ing equations will be verified through numerical simulations.

3. Optimization of the Robot Parameters with Non-Newtonian Damping

The optimization framework is adapted to a robotic platform in which the inter-module coupling consists of a linear elastic spring arranged in parallel with a nonlinear, non-Newtonian damper. The damping coefficient is not constant, but instead exhibits a strong dependence on the rotor angular velocity, reflecting the rheological properties of shear-thinning or shear-thickening fluids. The velocity-dependent damping is described by Equation (2). This formulation enables a systematic investigation of how ($c_0,\beta ,n$), together with the actuation frequency, influence the trade-offs between locomotion efficiency, stability, and energy consumption.

3.1. Performance Metrics and Evaluation Protocol

The evaluation of locomotion performance is based on a set of metrics that quantify both efficiency and motion quality. The average forward velocity and total actuation energy are retained as primary indicators, while smoothness and contact integrity are included to capture the additional effects introduced by velocity-dependent damping:

$$V_{\mathrm{avg}}={\displaystyle \frac{x_2\left(T_f\right)-x_2\left(T_i\right)}{T_f-T_i}},$$

(14)

$$E_{\mathrm{act}}={\int }_{T_i}^{T_f}\left|P\left(t\right)\right|dt,$$

(15)

$$J_{\mathrm{rms}}=\sqrt{{\displaystyle \frac{1}{T_f-T_i}}{\int }_{T_i}^{T_f}{\dot{x}}_2^2\left(t\right)dt}.$$

(16)

$$C=\underset{t\in [T_i,T_f]}{min}{N}_2\left(t\right).$$

(17)

where $N_2$ denotes the normal ground reaction force on the front module. The initial $N_{\mathrm{tr}}$ cycles of each simulation are discarded to eliminate transients, and metrics are computed over the subsequent $N_{\mathrm{ss}}$ cycles within the interval $[T_i,T_f]$.

3.2. Constraints

The feasible design space is further delimited by physical and mechanical constraints to ensure model validity and hardware safety:

$$C\ge 0,$$

(18)

$$|x_2\left(t\right)-x_1\left(t\right)-d|\le {\Delta }_{max},$$

(19)

$$0<{\omega }_{min}\le \omega \le {\omega }_{max},$$

(20)

$$0<m\le m_{max},0<k\le k_{max}.$$

(21)

3.3. Non-Dimensional Formulation

To facilitate comparisons across a wide range of operating points, the system dynamics are reformulated in terms of dimensionless groups that naturally incorporate the non- Newtonian damping law:

$$\mu ={\displaystyle \frac{m}{M}},\kappa ={\displaystyle \frac{k}{(M+m){\omega }^2}},$$

(22)

$${\Xi }_0={\displaystyle \frac{c_0\omega }{M+m}},\Xi \left(\omega \right)={\Xi }_0\left(1+{\beta \left|\omega \right|}^n\right).$$

(23)

$$\lambda ={\displaystyle \frac{L}{d}},\widehat{t}=\omega t,\widehat{x}={\displaystyle \frac{x}{d}}.$$

(24)

With $(\mu ,\lambda )$ fixed, performance can thus be primarily mapped as a function of $(\kappa ,{\Xi }_0,\beta ,n)$ and $\omega$, highlighting the interplay between stiffness, damping, and frequency.

3.4. Parametric Study

The influence of system parameters on locomotion performance is examined by considering their effects on $V_{\mathrm{avg}}$, $E_{\mathrm{act}}$, and $J_{\mathrm{rms}}$. The complete set of design variables is defined as $\theta =\left\{M,m,k,\omega ,c_0,\beta ,n\right\}.$

A Latin Hypercube Design (LHD) Viana, Venter and Balabanov (2010) of size $N_{\mathrm{samp}}$ is generated across prescribed parameter ranges:

$$\begin{array}{ccc}\hfill M& \in [M_{min},M_{max}],m\hfill & \hfill \in [m_{min},m_{max}],\\ \hfill k& \in [k_{min},k_{max}],\omega \hfill & \hfill \in [{\omega }_{min},{\omega }_{max}],\\ \hfill c_0& \in [c_{0,min},c_{0,max}],\beta \hfill & \hfill \in [{\beta }_{min},{\beta }_{max}],n& \in [n_{min},n_{max}].\hfill \end{array}$$

(25)

For each sample, the governing equations introduced in Sec. 2 are integrated with the velocity-dependent damping model of (2), transients are discarded, and the performance metrics (14)–(17) are evaluated. The results are analyzed through: (i) main-effects plots of each parameter against $V_{\mathrm{avg}}$ and $E_{\mathrm{act}}$, (ii) sensitivity indices based on Sobol or Morris screening to rank parameter influence, and (iii) two- dimensional response maps in the $(\omega ,\kappa )$ space for fixed $(\beta ,n)$ to reveal the frequency–damping interaction patterns.

3.5. Multi-Objective Optimization Problem

The optimization task seeks parameter configurations that achieve high forward speed with low energy consumption and smooth motion. This is formulated as a multi-objective problem:

$$\underset{\theta }{max}{V}_{\mathrm{avg}},\underset{\theta }{min}{E}_{\mathrm{act}},\underset{\theta }{min}{J}_{\mathrm{rms}}.$$

(26)

$$C\ge 0,|x_2-x_1-d|\le {\Delta }_{max},\theta \in \mathcal{B},$$

(27)

where $\mathcal{B}$ the box of admissible parameter bounds. For scalar selection of a single operating point, a weighted aggregate objective is considered:

$$J=-w_V{\displaystyle \frac{V_{\mathrm{avg}}}{V_{\mathrm{ref}}}}+w_E{\displaystyle \frac{E_{\mathrm{act}}}{E_{\mathrm{ref}}}}+w_S{\displaystyle \frac{J_{\mathrm{rms}}}{J_{\mathrm{ref}}}},$$

(28)

with ($w_V$, $w_E$, $w_S$) determined by the relative priority of speed, energy, and smoothness.

3.6. Optimization Algorithm

A population-based evolutionary strategy (NSGA-II) Deb et al. (2002) is employed to approximate the Pareto front, as described in Giagkiozis and Fleming (2014), in the objective space (Vavg,Eact,Jrms). Due to the enlarged design space arising from the non-Newtonian damping law, surrogate models such as radial basis functions, such as Buhmann (2000), or Gaussian processes, such as Seeger (2004), are trained on the LHD dataset. These surrogates are used for iterative infill based on Expected Improvement (EI) to accelerate convergence. At prescribed intervals, selected candidate solutions are re-evaluated with the high-fidelity dynamics model to ensure accuracy and correct for surrogate bias.

3.7. Reproducibility

Reproducibility is supported by providing explicit bounds for all parameters, initialization seeds for the LHD design, solver tolerances, and post-processing routines for calculating $V_{\mathrm{avg}}$, $E_{\mathrm{act}}$, and $J_{\mathrm{rms}}$. This ensures that all reported Pareto sets and optimization results can be independently validated and replicated.

4. Results and Discussion

The joint influence of spring stiffness k, baseline damping $c_0$, and rotor speed $\omega$ on the performance of the vibration–driven robot with non-Newtonian damping was evaluated numerically.

Representative time histories in Figure 3 show alternating stick–slip episodes and a persistent phase lead of the front module relative to the rear. The phase difference is central to net drift generation; it is amplified by appropriate stiffness and moderated by damping. No lift-off was observed at the operating point shown, consistent with the contact constraint satisfied in the simulations. The parameters are based on Table 1, where $c_0=1.4N/ms$, $n=1$, and $\beta =0.05$.

Given that the overall damping factor is modeled as a function of the angular velocity $\omega$, a comprehensive analysis was conducted to investigate how variations in $\omega$ influence key performance metrics of the robot. Specifically, the study examined the impact of changing angular velocity on drift velocity, total energy consumption, average power output, and the Cost of Transport (CoT). The CoT is a dimensionless parameter that quantifies the energy efficiency of locomotion, defined as the amount of energy required to move a unit mass over a unit distance.

To assess these effects, the angular velocity $\omega$ of the rotating mass was varied systematically from $20\mathrm{rad}/\mathrm{s}$ to $100\mathrm{rad}/\mathrm{s}$. The results, illustrated in Figure 4a–d, indicate a clear trend: as the rotational speed increases, there is a corresponding rise in drift velocity, energy consumption, and average power. With the robot almost not moving at angular velocities $\omega$ below 20 $\mathrm{rad}/\mathrm{s}$, this behavior suggests that while higher angular velocities enhance the translational motion of the robot (as reflected in increased drift velocity), they do so at the cost of greater energy expenditure and power demands, thereby influencing the overall energetic efficiency of the system.

However, the CoT surface in Figure 4d presents a valley that aligns with the drift ridge but is shifted toward slightly higher damping, indicating that minimal transport cost occurs where forward rectification remains effective while excessive oscillatory energy is curtailed by damping. This valley provides a practical design corridor for energy–limited operation.

Figure 5 demonstrates the relationship between total energy consumption and stiffness (k) for two level of damping coefficients, low: $c_0=0.30$ and high: $c_0=2.50$. For the lower damping coefficient ($c_0=0.30$), total energy increases with stiffness, peaking around $k=120\mathrm{N}/\mathrm{m}$, and then slightly decreasing as stiffness continues to increase. This suggests that with lower damping, the robot requires more energy as stiffness increases, but it reaches a maximum energy consumption at a certain stiffness level. On the other hand, for the higher damping coefficient ($c_0=2.50$), the total energy remains almost constant across the range of stiffness values, indicating that higher damping results in a more energy-efficient system, regardless of stiffness changes.

Figure 6 shows how drift velocity changes with stiffness for the same two damping values. For $c_0=0.30$, drift velocity increases significantly with stiffness, peaking around $k=100\mathrm{N}/\mathrm{m}$ before gradually decreasing with further increases in stiffness. This indicates that the robot’s drift velocity improves with higher stiffness up to a certain point, after which further stiffening has diminishing effects. However, for the higher damping coefficient ($c_0=2.50$), drift velocity remains nearly constant across the range of stiffness values, showing that higher damping drastically reduces the impact of stiffness on the robot’s motion, resulting in a less responsive system.

Figure 7 illustrates the average power required for the robot as a function of stiffness. For the lower damping coefficient ($c_0=0.30$), average power initially increases with stiffness, peaking and then decreasing at higher stiffness values, suggesting that maximum power consumption occurs at mid-range stiffness. In contrast, for the higher damping coefficient ($c_0=2.50$), the power consumption remains steady across stiffness levels, indicating a stable but less efficient performance. This behavior highlights that while high damping can provide stability, it leads to a less efficient power usage, particularly as stiffness increases.

The results demonstrate that locomotion efficiency arises from a three–way balance between elastic energy storage, velocity–dependent dissipation, and friction–mediated rectification. Quantitative ranges extracted from the simulations indicate average power between $\sim 0.052$–$0.067$ W and total actuation energy between $\sim 1.68$–$1.73$ J across the $(k,c_0)$ domain considered, while the cost of transport (CoT) collapses by orders of magnitude once $\omega$ exceeds a threshold near 50 rad/s.

The surface in Figure 8a shows that average power generally increases with both k and $c_0$, reflecting higher elastic forces and larger viscous losses. Non-monotonic regions appear where added damping stabilizes the stick–slip cycle and reduces wasted oscillation, producing shallow power minima for intermediate $(k,c_0)$.

Figure 8b exhibits a more uniform rise of total energy with stiffness, consistent with the larger work required to sustain oscillations against a stiffer coupling. Increasing $c_0$ produces a gradual energy increase due to enhanced dissipation. Together, Figure 8a,b indicate that tuning stiffness and damping independently is suboptimal; instead, an intermediate band in $(k,c_0)$ delivers acceptable power and energy budgets while maintaining viable oscillation amplitudes.

The corresponding drift–velocity map shown in Figure 8c reveals a ridge at moderate stiffness and low-to-moderate damping. Insufficient damping results in nearly symmetric oscillations with limited rectification, whereas excessive damping suppresses oscillation amplitude and reduces stride length. The observed ridge therefore reflects the parameter combination that maximizes phase-lag asymmetry between the modules under Coulomb friction.

The Pareto set in Figure 8d captures the trade–off between speed and energy. Knee–point candidates provide balanced operating points with significant velocity gains for modest energy penalties. Coloring by n indicates that shear–thinning cases ($n<1$) populate the upper–left portion of the Pareto front, corresponding to higher velocity for a given energy expenditure, whereas shear–thickening cases ($n>1$) cluster toward lower velocity and higher energy regions. This trend is consistent with stronger rate–dependent damping at elevated angular velocities $\omega$.

The results presented in Figure 10 illustrate the influence of the four key parameters $(k,\omega ,c_0,\beta )$ on the average velocity $V_{\mathrm{avg}}$ and actual energy expenditure $E_{\mathrm{act}}$. The plots demonstrate that $V_{\mathrm{avg}}$ is relatively insensitive to variations in stiffness k, with most data points concentrated near zero regardless of the stiffness values, except for a few isolated outliers. In contrast, $\omega$ exhibits a clearer trend, where increasing frequency generally leads to higher values of $V_{\mathrm{avg}}$, suggesting that oscillatory dynamics play a stronger role in enhancing system speed compared to stiffness adjustments.

Figure 9. Main effects of model parameters on average velocity ($V_{\mathrm{avg}}$) and actual energy expenditure ($E_{\mathrm{act}}$). The top row illustrates the variation of $V_{\mathrm{avg}}$ with stiffness (k) and angular frequency ($\omega$), whereas the bottom row shows the dependence of $E_{\mathrm{act}}$ on the damping coefficient ($c_0$) and nonlinearity parameter ($\beta$). The results indicate that $V_{\mathrm{avg}}$ is more sensitive to changes in $\omega$ than to variations in k, while $E_{\mathrm{act}}$ exhibits greater variability with respect to both $c_0$ and $\beta$, with occasional outliers influencing the observed trends. Overall, $\omega$ and $\beta$ appear to exert stronger influence on system performance than k and $c_0$.

Figure 9. Main effects of model parameters on average velocity ($V_{\mathrm{avg}}$) and actual energy expenditure ($E_{\mathrm{act}}$). The top row illustrates the variation of $V_{\mathrm{avg}}$ with stiffness (k) and angular frequency ($\omega$), whereas the bottom row shows the dependence of $E_{\mathrm{act}}$ on the damping coefficient ($c_0$) and nonlinearity parameter ($\beta$). The results indicate that $V_{\mathrm{avg}}$ is more sensitive to changes in $\omega$ than to variations in k, while $E_{\mathrm{act}}$ exhibits greater variability with respect to both $c_0$ and $\beta$, with occasional outliers influencing the observed trends. Overall, $\omega$ and $\beta$ appear to exert stronger influence on system performance than k and $c_0$.

Figure 10. Correlation-based sensitivity of $V_{avg}$ to ($(k,\omega ,c_0,\beta ,n)$).

Figure 10. Correlation-based sensitivity of $V_{avg}$ to ($(k,\omega ,c_0,\beta ,n)$).

On the other hand, the energy expenditure $E_{\mathrm{act}}$ does not exhibit a strong deterministic dependence on either $c_0$ or $\beta$, as indicated by the scattered distributions. While $c_0$ appears to influence $E_{\mathrm{act}}$ only weakly across its range, the effect of $\beta$ is slightly more pronounced, with occasional large increases in energy cost at higher $\beta$ values. These outliers suggest the possibility of nonlinear interactions or unstable operating regimes that may significantly elevate the energy requirements of the system.

These analyses highlight that while velocity performance is governed primarily by the oscillation frequency $\omega$, energy expenditure is more susceptible to nonlinear effects associated with $\beta$. This distinction implies that tuning $\omega$ is the most effective strategy for improving speed, whereas careful regulation of $\beta$ is crucial for maintaining energy efficiency.

Sensitivity analysis in Figure 10 ranks $\omega$ and k as the dominant contributors to drift velocity, with $c_0$ exerting a secondary influence and $(\beta ,n)$ becoming influential primarily in the high–frequency regime where $c\left(\omega \right)$ departs significantly from $c_0$. Although correlation–based indices provide only a first–order ranking, the resulting hierarchy aligns with the underlying physical interpretation: frequency determines the stride cadence, stiffness governs oscillation amplitude and phase relationships, and damping shapes rate–dependent energy dissipation and symmetry breaking.

5. Conclusions and Future Work

The findings indicate that the locomotion efficiency of the proposed vibration–driven robot, consisting of two modules linked by a spring and a damper whose resistance depends on angular velocity, arises from a subtle balance among actuation frequency, stiffness, and nonlinear damping. At low frequencies, inefficient translation results from limited rectification and small displacements per cycle. As the frequency increases, rectification becomes more effective, leading to rapid growth in average velocity and a corresponding reduction in the cost of transport. Beyond a certain threshold, however, further increases in frequency amplify the rate–dependent damping, which curtails oscillation amplitude and limits additional gains in performance. This establishes a transition regime in which efficiency improves sharply before saturating at higher frequencies.

The interaction between stiffness and damping further refines this balance. Intermediate stiffness values provide sufficient compliance for displacement while avoiding excessive elastic energy storage that would otherwise increase the power requirement. Similarly, nonlinear damping modulates system performance depending on whether the fluid operates in a shear–thinning or shear–thickening regime. In shear–thinning fluids, reduced effective damping at moderate frequencies permits larger stride lengths and higher velocities without proportional increases in energy expenditure. Conversely, in shear–thickening regimes, frequency–dependent damping suppresses oscillations, enhancing stability but reducing drift velocity unless compensated through simultaneous adjustments in stiffness and actuation frequency.

Furthermore, the resulting Pareto front was visualized using color maps corresponding to either n (representing shear–thinning versus shear–thickening behavior) or $\beta$ (representing damping sensitivity). Candidate solutions were highlighted through weighted aggregation of objectives and by identifying knee–points along the Pareto front. Constraint activity, such as the proportion of time during which $N_2\approx 0$, was also evaluated to exclude operating conditions approaching ground–contact loss.

Theoretical arguments and simulation results revealed several characteristic trends. Increasing either $\beta$ or n enhances $c\left(\omega \right)$ at elevated frequencies, which generally reduces oscillation amplitude and lowers $J_{\mathrm{rms}}$, but also suppresses forward velocity unless $\omega$ is adjusted concurrently. In shear–thinning regimes ($n<1$), moderate angular velocities result in lower effective damping, enabling larger stride lengths per actuation cycle and improved velocity at comparable energy expenditure. Consequently, the optimal combination of frequency and stiffness $(\omega ,k)$ is not fixed, but shifts according to $(\beta ,n)$, highlighting the need for integrated co–design of stiffness and damping characteristics.

These findings establish a clear design principle: optimal locomotion emerges when the excitation frequency is tuned slightly above the rectification threshold, where net displacement increases more rapidly than energy cost, while material parameters are selected to balance compliance with controlled damping. This emphasizes the importance of co–design strategies that integrate actuation frequency with structural and rheological parameters to achieve efficient transport. Moreover, the identification of parameter regimes that naturally align high velocity with low energy expenditure provides a rational basis for guiding future experimental validation and practical implementation of non–Newtonian vibration–driven locomotion systems.