1. Introduction
Underwater robots are essential for marine exploration and resource exploitation, with their performance largely governed by propulsion efficiency and agility (Yoerger et al. Reference Yoerger2021; Hartmann et al. Reference Hartmann, Baskaran, Raynaud, Benbedda, Mulleners and Shea2025). Conventional propeller-based systems often suffer from inefficiency in complex aquatic environments. Inspired by the exceptional swimmers of nature, bioinspired propulsion offers a promising solution by mimicking the locomotion of aquatic organisms (Lane, Sfakiotakis & Davies Reference Lane, Sfakiotakis and Davies1998; Li et al. Reference Li2021). Among these, body and/or caudal fin (BCF) propulsion emulates fish undulation and delivers exceptional thrust and agility, making it particularly well-suited for advanced underwater robotics (Webb Reference Webb1984).
Early biomimetic efforts focused on replicating the undulatory motion of fish. These systems typically employed sequences of servo motors to prescribe joint trajectories, generating undulatory movement (Triantafyllou & Triantafyllou Reference Triantafyllou and Triantafyllou1995; Hu et al. Reference Hu, Liu, Dukes and Francis2006; Liu & Hu Reference Liu and Hu2007). While these approaches succeeded in mimicking gross body waveforms, their propulsion efficiency remained limited due to the absence of dynamic coupling between the swimmer and the surrounding fluid. In contrast to strictly pre-programmed kinematics, biological swimmers continuously adapt their body movements through dynamic interactions with the surrounding fluid environment (Long Reference Long and John1998; Altringham & Ellerby Reference Altringham and Ellerby1999). Specifically, they utilised neuromuscular deformation to oscillate their bodies, disturbing the surrounding flow field and generating hydrodynamic forces. This, in turn, influences the body’s oscillations and forward motion (Pabst Reference Pabst1996; Dickinson et al. Reference Dickinson, Farley, Full, Koehl, Kram and Lehman2000). Inspired by this mechanism, researchers have introduced elastic elements at joint interfaces to confer mechanical adaptability. These elements, such as compliant hinges (Matthews et al. Reference Matthews, Zhu, Wang, Dong, Bart-Smith and Lauder2022) and soft actuators (El-Atab et al. Reference El-Atab, Mishra, Al-Modaf, Joharji, Alsharif, Alamoudi, Diaz, Qaiser and Hussain2020), enable dynamic responses to external hydrodynamic forces, driving body deformation and motion. This strategy enhances the coordination between the propulsors and the fluid, thereby improving propulsion efficiency (Li et al. Reference Li, Xiao, Liu, Hu, Li, Li, Liu, Hu and Wen2018b ; Renda et al. Reference Renda, Giorgio-Serchi, Boyer, Laschi, Dias and Seneviratne2018). Building on this foundation, subsequent studies have employed continuously flexible bodies to better capture the smooth, uninterrupted deformations characteristic of fish locomotion (Arora et al. Reference Arora, Kang, Shyy and Gupta2018; Zhu et al. Reference Zhu, White, Wainwright, Di Santo, Lauder and Bart-Smith2019; Wu et al. Reference Wu, Zhang, Tian, Li and Lu2020). These flexible designs are capable of adapting to the surrounding flow, closely mimicking the aquatic swimmers. Nevertheless, most biomimetic systems are built upon linear stiffness models (Wang et al. Reference Wang, Ren, Hu, Li and Sitti2021), which fail to capture the complex nonlinear time-varying stiffness observed in biological swimmers during locomotion.
Schematics of a marine swimmer (using yellowfin tuna as an example). (a) Biological model: caudal peduncle connects anterior body muscles to the tail. (b) Tendon variations within the caudal peduncle during the coupling heave/pitch motions exhibited by the tail. (c) Conceptual model: the body musculature provides energy, the caudal peduncle exhibiting spring-like linear/torsional behaviour, and the tail conducting heave/pitch motion for propulsion. (d) Mechanical model of the self-propulsor, where
$ y_1(t) = h_0 \sin (2\pi f t)$
represents the energy input,
$(k_y, c_y,k_\theta , c_\theta )$
model the mechanical properties of the caudal peduncle, the black flat plate represents the tail fin, and the solid red arrow indicates the stiffness tuning concept model
$ k_{\textit{Non}}$
along the heave direction.
In contrast, biological organisms dynamically modulate their body stiffness in real time, a strategy that significantly enhances their locomotor performance (Quinn & Lauder Reference Quinn and Lauder2021; Liu et al. Reference Liu, Liu, Liang, Ren and Ren2023). For instance, fish generate forward thrust through a highly coordinated interplay between body undulation and tail motion, as exemplified by tuna species (figure 1 a). Moreover, fish modulate their body stiffness through muscle–tendon interactions, as is evident in figure 1(b) (Zhong et al. Reference Zhong, Zhu, Fish, Kerr, Downs, Bart-Smith and Quinn2021). This real-time stiffness modulation represents a functional distinction between biological systems and most existing biomimetic designs. While several empirical studies suggest that stiffness modulation contributes to enhanced propulsive efficiency (Li et al. Reference Li, Jiang, Wang and Yu2018a ; White, Lauder & Bart-Smith Reference White, Lauder and Bart-Smith2021), the underlying physical mechanism remains elusive. Hence a key challenge lies in developing a model that captures the tunable stiffness characteristics of fish bodies and elucidating its role in enhancing propulsive performance. Revealing this mechanism requires solving a complex fluid–structure interaction (FSI) problem, which fundamentally governs thrust production. To address this challenge, traditional hydrodynamic models, such as quasi-steady lift theory (Chopra Reference Chopra1974), Lighthill’s slender body theory (Lighthill Reference Lighthill1969) and fluctuating plate theory (Wu Reference Wu1961), use empirical formulas to approximate hydrodynamic forces. However, they often fail to capture the interactions between body deformation and wake dynamics, which are essential for understanding unsteady hydrodynamic behaviour. Computational fluid dynamics methods (Wei, Mark & Hutchison Reference Wei, Mark and Hutchison2014; Liu et al. Reference Liu, Ren, Dong, Akanyeti, Liao and Lauder2017) offer high-fidelity solutions but often come with significant computational costs, particularly when applied to nonlinear problems arising from tunable stiffness. This limitation poses challenges for their use in exploring physical mechanisms. Thus balancing modelling accuracy with computational efficiency becomes essential for uncovering the underlying physical mechanisms. Among various approaches, the vortex sheet method emerges as a widely adopted and efficient technique for capturing boundary layer development and wake vorticity dynamics. It has been effectively utilised to analyse flow structures around swimming bodies, providing key insights into their unsteady hydrodynamic behaviour (Nitsche & Krasny Reference Nitsche and Krasny1994; Jones Reference Jones2003; Huang et al. Reference Huang, Ristroph, Luhar and Kanso2018). Building on this approach, a modified nonlinear vortex sheet method (NVSM) is developed to better address instability challenges in nonlinear FSI simulations.
In this study, a bioinspired conceptual model of a self-propulsor is developed to reflect the functional differentiation observed across body regions in swimming fish. Moreover, a parabolic stiffness-tuning strategy is proposed to emulate the stiffness modulation exhibited by the musculature in the caudal peduncle. The propulsion benefits enabled by tunable stiffness are examined through analyses of the propulsor’s dynamic response and the associated vortex dynamics. And the NVSM is developed to address the nonlinear FSI issues. This approach balances computational efficiency with the ability to reveal propulsive mechanism. Numerical results demonstrate that propulsors with nonlinear tunable stiffness can double their speed and efficiency across nearly the entire frequency range. The underlying mechanism is attributed to a snap-through instability, which promotes the formation of quasi-quadrilateral swimming kinematics and amplifies vortex generation, thereby enhancing propulsive performance. The paper is organised as follows. Section 2 presents the development of the bioinspired conceptual model along with the coupled fluid–structure numerical algorithm. In § 3, a detailed analysis of the impact of stiffness tuning on propulsive performance across the entire frequency range is provided. Additionally, the influence of stiffness modulation on the system’s dynamic response and wake evolution is discussed, revealing the underlying mechanism of high-efficiency biomimetic propulsion. Finally, § 4 summarises the main findings and offers concluding remarks.
2. Physics problem and mathematical formulation
This section presents a bioinspired conceptual model that captures the functional differentiation along the body of swimming fish, incorporating a parabolic stiffness-tuning strategy to emulate the time-varying stiffness of the caudal musculature. Based on this model, the governing equations for both the structure and the fluid are derived, and a fully coupled FSI model is established. An efficient numerical algorithm is further developed to resolve the nonlinear interactions between the propulsor and its wake, enabling accurate prediction of unsteady flow features while maintaining computational efficiency.
2.1. Bioinspiration and model development
Aquatic swimmers employing the BCF propulsion mode, such as the yellowfin tuna, typically consist of three anatomically and functionally distinct regions: the body musculature, the caudal peduncle, and the tail fin, as illustrated in figure 1(a). During the swimming process, energy is primarily generated by the body musculature, transmitted through the compliant caudal peduncle, and ultimately manifested as oscillatory motion of the tail fin (Ren et al. Reference Ren, Yang, Wang and Wen2016). This motion interacts with the surrounding fluid, shedding vorticity and forming a vortex wake (Taneda Reference Taneda1965; Triantafyllou & Triantafyllou Reference Triantafyllou and Triantafyllou1993). Variations in the flow field near the swimmer lead to changes in the pressure difference between the two sides of the structure, thereby propelling the fish forwards (Liao et al. Reference Liao, Beal, Lauder and Triantafyllou2003). Notably, fish modulate stiffness via muscle contractions to achieve efficient swimming performance. As conceptually depicted in figure 1(b), the tendons reach minimal stiffness at state 1 (the relaxed position), while progressively hardening with increasing distance from state 1, as exemplified in states 2 and 3 (Lauder et al. Reference Lauder, Madden, Tangorra, Anderson and Baker2011).
Drawing inspiration from the swimming process, we propose a conceptual model that captures the coordinated interaction of various anatomical components, as illustrated in figure 1(c). The axial musculature is represented as a periodic excitation source, while the caudal peduncle is abstracted as a compliant transmission unit. This unit modulates energy transfer and governs the tail fin’s oscillatory motion in the heave and pitch directions.
To investigate the role of stiffness modulation in propulsion performance, a simplified structural model is developed, as shown in figure 1(d). The system is analysed in a global inertial reference frame
$(x, y)$
. The input from the axial musculature is idealised as a sinusoidal heaving displacement given by
$ y_1(t) = h_0 \sin (2\pi f t)$
, where
$ h_0$
denotes the amplitude, and
$ f$
denotes the excitation frequency. The caudal peduncle is modelled using two independent spring–damper pairs:
$(k_y, c_y)$
along the heave direction, and
$(k_\theta , c_\theta )$
along the pitch direction. Additionally, a parabolic stiffness-tuning strategy is incorporated along the
$y$
-direction, as indicated by the solid red arrow in figure 1(d). The tail fin is idealised as a rigid flat plate of length
$c$
, thickness
$e \leqslant c$
, and uniform density
$\rho$
. Point
${A}$
marks the leading edge of the plate, while point
${C}$
denotes its centre of mass. The translational velocity at
${C}$
is given by
$(\dot {x}_{{C}}, \dot {y}_{{C}})$
, and the angular velocity about the joint is
$\dot {\theta }$
, where the overdot indicates time differentiation.
2.2. Governing equations and vortex sheet model
In the context of bionic swimming, the fluid is assumed to be incompressible with low dimensionless viscosity, corresponding to a high Reynolds number regime. The swimmer is immersed in this fluid environment and is assumed to propel itself in the negative
$x$
-direction. The governing equations of motion for the swimmer can be derived from Newton’s second law:
\begin{equation} \left \{ \begin{aligned} & m{{\ddot {x}}_{{C}}} = -F\sin \theta + D\cos \theta ,\\ & m{{\ddot {y}}_{{C}}} + {{k}_{y}}({{y}_{{A}}} - {{y}_{1}}) + {{c}_{y}}({{\dot {y}}_{{A}}} - {{\dot {y}}_{1}}) + {{F}_{\textit{Non}}} = F\cos \theta + D\sin \theta ,\\ & {{I}_{{A}}}\ddot {\theta } + {{c}_{\theta }}\dot {\theta } + {{k}_{\theta }}\theta = {{M}_{{A}}} + {{M}_{\textit{Inertia}},} \end{aligned} \right . \end{equation}
where
$m, {{I}_{{A}}}$
are the swimmer’s mass and the moment of inertia about the leading edge, respectively. Here,
${y}_{{A}},\dot {y}_{{A}}$
and
${y}_{{C}},\dot {y}_{{C}}$
denote the displacement and velocity at points
${A}$
and
${C}$
, respectively. They satisfy the relationships
${{y}_{{A}}}={{y}_{{C}}}-({c}/{2})\sin \theta$
and
${\dot {y}_{{A}}}={\dot {y}_{{C}}} - ({c}/{2})\cos \theta\, \dot {\theta }$
. The added moment of inertia
${M}_{\textit{Inertia}} = m \operatorname {Im} [-\boldsymbol{AC} \boldsymbol{\cdot} \overline {{a}_{{A}}} ]$
accounts for the inertial contribution when considering rotation about a point other than the centre of mass. Here,
$\operatorname {Im}$
represents the imaginary part in complex fields,
$\boldsymbol{AC}$
denotes the vector from point
${A}$
to point
${C}$
, and
$\overline {{a}_{A}}$
is the complex conjugate of the acceleration
${a}_{A}$
at point
${A}$
, which can be written as
\begin{equation} {{a}_{{A}}} = {{\ddot {x}}_{C}} + \cos \theta\, {{\dot {\theta }}^{2}} + \sin \theta\, \ddot {\theta } + \text{i} \big ({{\ddot {y}}_{C}} + \sin \theta\, {{\dot {\theta }}^{2}} - \cos \theta\, \ddot {\theta } \big ), \end{equation}
where
$\text{i} = \sqrt {-1}$
, the first three terms represent the acceleration along the
$x$
-direction, and the remaining terms represent the acceleration along the
$y$
-direction.
To simulate the muscle tension contraction behaviour observed in biological swimmers, a stiffness tuning concept model
$ k_{\textit{Non}} = k_y ( -a + 3b y_{{A}}^2)$
as depicted by the red arrow in figure 1(d) is introduced. Here,
$ F_{\textit{Non}}$
arises from the stiffness tuning model
\begin{equation} F_{\textit{Non}} = k_y(-ay_{{A}} + b y_{{A}}^3). \end{equation}
The total equivalent stiffness
$ k_r$
, combining both linear and nonlinear components, is given by
\begin{equation} k_r = k_y (1 - a) + 3b k_y y_{{A}}^2. \end{equation}
This equivalent stiffness yields a parabolic stiffness tuning model characterised by low stiffness at the centre and higher stiffness at larger displacements (Chen et al. Reference Chen, Liao, Xia, Sun and Zhou2023), accurately capturing the stiffness variation in muscle-tendon tissues as shown in figure 1(d). The model facilitates the adjustment of parameters
$ a$
and
$ b$
to achieve various stiffness profiles, such as a quasi-zero-stiffness model suitable for vibration control (Jin et al. Reference Jin, Zhang, Zheng, Lu, Xu and Greaves2024), and a bistable model designed for a high-efficiency actuator (Chi et al. Reference Chi, Li, Zhao, Hong, Tang and Yin2022b
).
From an energy perspective, stiffness modulation also alters the system’s potential energy (further illustrated in figure 8 below), expressed as
\begin{equation} U_p = \frac {1}{2} k_y (1 - a) y_{{A}}^2 + \frac {1}{4} k_y b y_{{A}}^4. \end{equation}
The flow is governed by the incompressible Euler equations
\begin{equation} \left \{ \begin{array}{l} \boldsymbol{\nabla} \boldsymbol{u} = 0, \\ \dfrac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{u} = -\dfrac {1}{\rho _f} \boldsymbol{\nabla} p. \end{array} \right . \end{equation}
Here,
$ \boldsymbol{u}$
is the velocity vector,
$ p$
is the static pressure field, and
$ \boldsymbol{\nabla}$
denotes the gradient operator. However, directly solving the governing equations in partial differential form presents significant computational challenges. To capture the essential flow dynamics in bioinspired propulsion, a vortex dynamics approach is adopted, where vorticity evolution is tracked to reduce computational complexity. Specifically, the vortex sheet model is employed, in which the velocity field is reconstructed as a superposition of bound vortices
$(\zeta (s,t), \bar {\gamma }_b(s,t))$
distributed along the body surface, and free vortices
$(\zeta (s,t), \bar{\gamma }_f(s,t))$
shed into the wake (Nitsche & Krasny Reference Nitsche and Krasny1994), as illustrated in figure 2. The induced velocity field
$\omega (s,t)$
is given by the Biot–Savart law
\begin{equation} \overline {\omega }(s,t) = U(t) + \frac {1}{2\pi \text{i}} \int _{-1}^{s_{{max}}} \frac {\overline {\gamma }(s',t) \, \text{d}s'}{\zeta (s,t) - \zeta (s',t)}, \end{equation}
where the bar denotes the complex conjugate,
$s'$
is the integrative variable,
$-1 \leqslant s \leqslant 1$
denotes the contours of the bound vortex sheets, and
$1 \leqslant s \leqslant s_{{max}}$
denotes the contours of the free vortex sheets. The first term
$U$
on the right-hand side of (2.7) is the prescribed background flow, which is set to zero in this work. This approach effectively transforms the problem of solving the full velocity field into tracking the evolution of discrete vortex elements, offering a computationally efficient method for unsteady flow analysis.
Schematic of the vortex sheet model for a self-propulsive swimmer.
2.3. Boundary conditions and fluid–structure coupling
To complete the computation of the above equations, boundary conditions are imposed to ensure the consistency of the structural and fluid dynamics equations.
1. No-penetration condition, which ensures that the fluid does not pass through the structure’s surface. Mathematically, the normal component of the fluid velocity at the surface must match the normal velocity of the structure:
\begin{equation} \text{Re}\left ( \boldsymbol{n} \left ( \overline {\frac {\partial \zeta (s,t)}{\partial t}} - \overline {\omega (s,t)} \right ) \right ) = 0,\quad -1 \leqslant s \leqslant 1, \end{equation}
where
$\text{Re}$
indicates the real part in complex fields, and
$\boldsymbol{n} = \text{i}\boldsymbol{s} =\text{i}\, {\rm e}^{\text{i}\theta }$
denotes the normal vector of the swimmer.
2. Kelvin’s circulation theorem, which imposes the conservation of total circulation in the flow. It governs the shedding of free vortices from the trailing edge:
\begin{equation} \Gamma (t) = \int _{-1}^{1} \overline {\gamma }_b(s,t) \, \text{d}s = -\int _{1}^{s_{{max}}} \overline {\gamma }_f(s',t) \, \text{d}s', \end{equation}
where
$\Gamma (t)$
denotes the total circulation around the swimmer.
3. Pressure continuity at the trailing edge, which ensures smooth vortex shedding at the trailing edge of the structure (Jones Reference Jones2003):
\begin{equation} [p](1,t) = 0, \end{equation}
These conditions link the swimmer’s motion to the surrounding fluid. Based on this foundation, the surface pressure distribution is determined by applying the unsteady Bernoulli equation within the potential flow framework (Alben Reference Alben2009):
\begin{equation} [p](s,t) = \int _{1}^{s} \partial _t \bar {\gamma }_{b}(s',t) \, \text{d}s' + \left ( \mu (s,t) - \xi (s,t) \right ) \bar {\gamma }_{b}(s,t) + \dot {\Gamma }(t), \quad -1 \leqslant s \leqslant 1, \end{equation}
where
$\mu (s,t) = \operatorname{Re}\,(\boldsymbol{s}\, \overline {\omega }(s,t))$
denotes the component of the average fluid velocity tangent to the body surface, and
$\xi (s,t) = \operatorname {Re}(\boldsymbol{s}\, \partial _t \overline {\zeta }(s,t))$
denotes the component of the body’s velocity in the same tangential direction. The specific derivation process can be found in Fang et al. (Reference Fang, Ho, Ristroph and Shelley2017). The total hydrodynamic force
$F$
and torque
$M_A$
acting on the body are obtained by integrating the pressure distribution
$p$
along the body surface:
\begin{equation} \left \{ \begin{aligned} & F = F_x + \text{i}F_y = \int _{-1}^{1} \boldsymbol{n} [p] \, \text{d}s, \\ & M_{{A}} =\int _{-1}^{1} \boldsymbol{n} [p] (s+1) \, \text{d}s, \end{aligned} \right . \end{equation}
where the component
$F_x = \int _{-1}^{1} [p](-\sin \theta ) \, \text{d}s$
denotes the thrust force along the
$x$
-direction, and
$F_y = \int _{-1}^{1} [p](\cos \theta ) \, \text{d}s$
is the lift force along the
$y$
-direction.
Additionally, the drag force
$D$
in (2.1), emulating the effect of skin friction due to fluid viscosity, is introduced by (Fang et al. Reference Fang, Ho, Ristroph and Shelley2017)
\begin{equation} D = C_d \big( \overline {U_+}^{3/2} + \overline {U_-}^{3/2} \big), \end{equation}
where
$C_d$
denotes the skin friction coefficient; its value is set at
$C_d = 0.04$
based on the experiments of Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016).
Here,
$\overline {U_\pm }$
denotes the spatially averaged tangential slip velocity of the fluid relative to the body surface. It is evaluated separately on the upper (
$\mu _+$
) and lower (
$\mu _-$
) sides of the swimmer:
\begin{equation} \overline {U_\pm } = \frac {1}{s} \int _{-s/2}^{s/2} \mu _\pm (s',t) \, \text{d}s' - \xi (s,t). \end{equation}
Notably, evaluating the hydrodynamic force
$F$
, torque
$M_{{A}}$
and drag
$D$
requires knowledge of the pressure distribution (2.12) and the velocity field (2.13). These flow characteristics are determined by the vortex sheet model as described in (2.7), which in turn depends on the instantaneous motion of the swimmer. However, the swimmer’s motion
$({{x}_{{C}}},{{y}_{{C}}},\theta )$
itself is governed by the hydrodynamic forces via (2.1). This mutual dependency establishes a tightly coupled feedback loop, underscoring the nature of FSIs. Moreover, the presence of nonlinear elements, such as nonlinear forces
${F}_{\textit{Non}}$
and trigonometric functions
$\sin \theta ,\cos \theta $
, exacerbates the complexity and the labour-intensive nature of the computations. The strong nonlinear coupling between fluid and structure poses a significant computational bottleneck. To address this challenge, the NVSM is employed, which avoids conventional time-stepping schemes such as Runge–Kutta methods (Huang et al. Reference Huang, Ristroph, Luhar and Kanso2018; Hang et al. Reference Hang, Heydari, Costello and Kanso2022). Instead, the NVSM employs a tailored two-stage strategy: (i) discretising the governing equations in differential form into a set of nonlinear algebraic equations; and (ii) performing a concurrent iterative solution. The specific procedures of the NVSM are described in Appendix A.
2.4. Simulation parameters and metrics of swimming performance
This subsection details the numerical parameters used in the study, as well as the performance metrics adopted to quantify swimming efficiency and locomotion characteristics.
All governing equations and boundary conditions are expressed in dimensionless form. Assuming a neutrally buoyant swimmer, the characteristic scales for non-dimensionalisation are selected as
${c}/{2}$
,
$\rho$
and
$f$
. A summary of the main simulation parameters used in this study is provided in table 1. To facilitate the analysis, the following dimensionless variables are introduced.
Main parameters of the proposed swimming model.
\begin{align} \begin{aligned} & \text{Structure:}\quad {{c}^{*}}=c/(c/2);{{m}^{*}}=m/(c\rho /2);h_{0}^{*}={{h}_{0}}/(c/2);{{f}^{*}}=f/(f)=1;\\ & I_{{A}}^{*}={{I}_{{A}}}/({{c}^{3}}\rho /8); k_{y}^{*}={{k}_{y}}/(\rho {{f}^{2}}c/2);k_{r}^{*}={{k}_{r}}/(\rho {{f}^{2}}c/2);k_{\theta }^{*}={{k}_{\theta }}/(\rho {{f}^{2}}{{c}^{3}}/8). \\[3pt] & \text{State:}\quad {{t}^{*}}=t/(1/f);x_{{A}}^{*}={{x}_{{A}}}/(c/2);y_{{C}}^{*}={{y}_{{C}}}/(c/2);y_{{A}}^{*}={{y}_{{A}}}/(c/2); \\ & \dot {x}_{{C}}^{*}={{{\dot {x}}}_{{C}}}/(cf/2);\dot {y}_{{C}}^{*}={{{\dot {y}}}_{{C}}}/(cf/2);{{{\dot {\theta }}}^{*}}=\dot {\theta }/f; \\ & {{F}^{*}}=F/(\rho {{f}^{2}}{{c}^{2}}/4);{{D}^{*}}=D/(\rho {{f}^{2}}{{c}^{2}}/4);F_{x}^{*}={{F}_{x}}/(\rho {{f}^{2}}{{c}^{2}}/4);F_{y}^{*}={{F}_{y}}/(\rho {{f}^{2}}{{c}^{2}}/4); \\ & M_{{A}}^{*}={{M}_{{A}}}/(\rho {{f}^{2}}{{c}^{3}}/8). \\[3pt] & \text{Energy:}\,\,\, P_{in}^{*}={{P}_{in}}/(\rho {{f}^{3}}{{c}^{3}}/8);P_{out}^{*}={{P}_{out}}/(\rho {{f}^{3}}{{c}^{3}}/8);U_{P}^{*}={{U}_{P}}/(\rho {{f}^{3}}{{c}^{3}}/8). \end{aligned} \end{align}
To assess swimming performance, this study introduces three key metrics (Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998; Heydari & Kanso Reference Heydari and Kanso2021). The swimming speed, expressed in body lengths per second (
$\text{BL}\ \text{s}^{-1}$
), is defined as
$\bar {V} = \left | \int _{t}^{t+T} \dot {x}_{{C}}(t) \, \text{d}t \right | / (cT)$
at the steady state (where
$T = 1/f$
denotes the motion period). The period-averaged Froude efficiency
$\eta$
is introduced because swimming speed alone cannot adequately reflect energy utilisation efficiency:
\begin{equation} \eta = \frac {\int _{t}^{t+T} P_{{out}}(t) \, \text{d}t}{\int _{t}^{t+T} P_{\textit{in}}(t) \, \text{d}t}, \end{equation}
where the output power is
${P}_{out} = {F}_{x} \dot {x}_{{C}}$
, and the input power is similarly defined as
${P}_{in}$
, which equals the product of the base excitation velocity
${\dot {y}}_{1}$
and the reaction force. For the swimmer with linear stiffness, the input power
${P}_{in}$
is given by
\begin{equation} P_{\textit{in}} = \left [ k_y (y_1 - y_{{A}}) + c_y (\dot {y}_1 - \dot {y}_{{A}}) \right ] \dot {y}_1. \end{equation}
For the swimmer with nonlinear tunable stiffness, the input power
${P}_{in}$
is
\begin{equation} P_{\textit{in}} = \left [ k_r (y_1 - y_{{A}}) + c_y (\dot {y}_1 - \dot {y}_{{A}}) \right ] \dot {y}_1. \end{equation}
Also, the angle of attack
$\alpha$
is defined as the angle between the swimmer’s body gesture and the velocity vector
$(\dot {x}_{{C}}, \dot {y}_{{C}})$
, and is introduced to investigate the relationship between swimming motion and performance:
\begin{equation} \alpha = \tan ^{-1} \left ( \frac {\dot {y}_{{C}}}{-\dot {x}_{{C}}} \right ) - \theta . \end{equation}
3. Results and discussion
In this section, the swimming behaviours and vortex field evolution are systematically analysed to elucidate the mechanism of high-efficiency propulsion enabled by nonlinear tunable stiffness.
3.1. How well the swimmer performs: high swimming performance induced by nonlinear tunable stiffness
To illustrate the role of nonlinear tunable stiffness, a representative case with parameters
$ a^{*} = 1.5$
and
$ b^{*} = 2$
is analysed. For comparison, the linear stiffness model is chosen as the baseline, where the component
$ k_{\textit{Non}}$
is disregarded. As shown in figure 3(a), the equivalent stiffness varies parabolically with the heave displacement. Specifically, the system exhibits lower stiffness than the linear case near
$ y_{{A}}^{*} = 0$
, but exceeds it at larger displacements. This tunable stiffness aligns with the adaptive properties of biological muscles illustrated in figure 1(b).
Figure 3(b) quantifies the impact of stiffness modulation on swimming performance. It compares the average speed and efficiency of swimmers with linear and nonlinear tunable stiffness across a range of frequencies. It is evident that in the linear case, both swimming speed and efficiency exhibit a non-monotonic dependence on actuation frequency, resembling a resonance-like behaviour along the
$y$
-direction. Similar frequency-dependent performance peaks have been observed in both numerical and robotic models, such as flexible plate simulations (Peng et al. Reference Peng, Sun, Yang, Xiong, Wang and Wang2022) and fish-inspired robots (Quinn & Lauder Reference Quinn and Lauder2021). In the nonlinear scenario, a significant improvement in swimming performance is observed, with both swimming speed and efficiency highlighted in purple consistently outperforming their linear counterparts across the entire frequency range. The peak swimming speed of the nonlinear tunable stiffness model (
$\bar {V} = 5.6 \ \textrm{BL}\ {\textrm s^-{^1}}$
) is nearly twice that of its linear counterpart (
$\bar {V} = 3 \ \textrm{BL}\ {\textrm s^-{^1}}$
). One key factor underlying this performance gain is the reduction in angle of attack, as shown in figure 3(c). Specifically, the incorporation of tunable stiffness reduces the angle of attack across the entire frequency range, with the reduction particularly pronounced at low excitation frequencies. This suppression is crucial for maintaining high performance, as excessively high angles of attack may lead to dynamic stall of the tail, thereby diminishing propulsion efficiency (Walker & Westneat Reference Walker and Westneat2002).
To further elucidate the influence of the angle of attack on swimming performance, figure 3(d) presents its time history at
$f = 0.3 \ \text{Hz}$
. In addition, representative swimming postures and velocity vector diagrams at three evenly spaced time intervals are provided. Notably, the angle of attack
$\alpha$
is defined as the angle between the swimmer’s body posture and the velocity vector
$(\dot {x}_{{C}}, \dot {y}_{{C}})$
. The results indicate that in the linear case, the swimmer remains predominantly horizontal throughout the swimming cycle. In this orientation, the fluid pressure acting normal to the body surface is directed primarily along the
$y$
-axis. Consequently, it cannot be effectively redirected in the forward direction, resulting in significantly reduced thrust, and overall performance as illustrated in figure 3(b). The introduction of nonlinear tunable stiffness influences the heave motion, which directly affects the angle of attack as described by (2.19). The snapshots presented in figure 3(d) illustrate that this effect results in a tuned phase relationship between the pitch and heave motions, reducing the angle of attack. This phase tuning enables a smaller angle of attack to be achieved while maintaining a large pitching amplitude. Consequently, the fluid pressure acting perpendicular to the swimmer is more effectively redirected forwards, enhancing propulsion efficiency.
(a) Equivalent stiffness as a function of heaving motion for linear and nonlinear tunable stiffness models. (b) Period-averaged swimming speed and efficiency as functions of frequency for both models. (c) Maximum angle of attack over the full frequency range. (d) Time history of the angle of attack at
$ f = 0.3$
Hz, with snapshots illustrating the swimmer’s body posture (solid lines) and velocity vectors (red dashed arrows) at three evenly spaced intervals. All cases use
$ a^{*} = 1.5$
,
$ b^{*} = 2$
.
As shown in figure 3(b), the two models exhibit distinct optimal frequencies for peak swimming speed. The tunable stiffness model achieves its maximum velocity at
$ f = 1\, \text{Hz}$
, whereas the linear model peaks at
$ f = 1.2\, \text{Hz}$
with
$ \bar {V} = 3 \, \text{BL}\,{\textrm s^-{^1}}$
. Notably, even at
$ f = 1.2\, \text{Hz}$
, the tunable stiffness model maintains high speed
$ \bar {V} = 5.4 \, \text{BL}\,{\textrm s^-{^1}}$
, significantly outperforming the linear counterpart. To clarify the characteristics of swimming performance, figure 4 presents column-wise comparisons of key hydrodynamic metrics
$\left | \dot {x}_{{C}}^* \right |$
,
$F_x^*$
,
$\Gamma$
,
$P_{\textit{in}}^*$
and
$P_{{out}}^*$
between swimmers with linear and nonlinear tunable stiffness. Figure 4(a,b) correspond to
$f = 1\, \text{Hz}$
, and figure 4(c,d) correspond to
$f = 1.2\, \text{Hz}$
.
A closer examination of the second column in figure 4 shows that the tunable stiffness model generates significantly greater thrust than its linear counterpart. This enhanced thrust directly contributes to improved propulsion performance, emphasising that the introduction of tunable stiffness leads to superior propulsion, independent of the chosen frequency. To understand these differences, the vortex structures are analysed, given their strong correlation with propulsion force as described by Kelvin’s circulation theorem (Drucker & Lauder Reference Drucker and Lauder2001). The swimmer equipped with nonlinear tunable stiffness shortens the vortex rolling duration and generates a wake characterised by minimal vortex rings and slight curvature. This distinct wake structure effectively reduces energy loss by minimising the capture of wake energy by large vortices, thereby enhancing momentum transfer into forward thrust (Sarpkaya Reference Sarpkaya1975; Hang et al. Reference Hang, Heydari, Costello and Kanso2022). Despite differences in flow fields, the fundamental mechanisms underlying the propulsive effects of tunable stiffness remain unclear. A deeper understanding of how tunable stiffness enhances propulsion will be pursued in the following subsections.
Instantaneous swimming performance for the swimmer with (a) linear stiffness at
$ f = 1 \, \text{Hz}$
, (b) nonlinear tunable stiffness at
$ f = 1 \, \text{Hz}$
, (c) linear stiffness at
$ f = 1.2 \, \text{Hz}$
, and (d) nonlinear tunable stiffness at
$ f = 1.2 \, \text{Hz}$
. First column shows time-dependent speed
$| \dot {x}_{{C}}^{*} |$
; second column shows thrust
$F_{x}^{*}$
; third column shows circulation
$\Gamma$
; fourth column shows input power
$P_{in}^{*}$
and output power
$P_{out}^{*}$
. Solid lines represent the instantaneous values, and dashed lines represent time-period averages, with all the results obtained from the steady state starting at
$t_s = 10T$
.
3.2. Why swimmers perform better: influence of nonlinear tunable stiffness on swimming behaviour and vortex evolution
To reveal the underlying high-efficienct propulsive mechanism, the dynamic behaviours and vorticity evolution in swimmers with both linear and nonlinear tunable stiffness are analysed at reference frequency
$ f = 1\, \text{Hz}$
.
Figure 5 (left-hand column) illustrates the phase diagrams of heave and pitch motions. Additionally, the corresponding hydrodynamic pressure distributions at four evenly spaced time instants within one cycle are illustrated. Meanwhile, the vorticity attached to the structure and the trailing-edge vortex at the corresponding instants are illustrated in the right-hand columns. Red represents an anticlockwise vortex (positive vorticity), while blue indicates a clockwise vortex (negative vorticity). These shed vortical structures undergo convection and rolling up, ultimately forming spiral-shaped vortex sheets (Zhu et al. Reference Zhu, Wang, Xu and Xiao2023).
In the linear case, the phase diagram in the left-hand column of figure 5(a) displays an elliptical shape, indicating a short residence time at large pitch angles. Correspondingly, the right-hand columns of figure 5(a) show that during the upward motion (
$ t = 0$
to
$t = 0.25T$
), clockwise vortex sheets are generated and shed from the upper surface at the trailing edge. These vortices are evenly spaced in the wake, leading to a region of lower velocity and higher pressure on the swimmer’s upper surface. The resulting downward pressure, coupled with a negative pitching angle, facilitates the generation of forward thrust (Liao et al. Reference Liao, Beal, Lauder and Triantafyllou2003) during this phase. As the leading edge moves downwards (from
$ t = 0.25T$
to
$ t = 0.5T$
), anticlockwise vortex sheets are shed from the bottom surface, altering the flow field, and increasing pressure on the swimmer’s lower surface. Consequently, the hydrodynamic force projected along the
$ x$
-axis contributes to forward propulsion. However, due to the consistently small pitching angle in this phase, the pressure-induced force is not efficiently converted into thrust, limiting the overall propulsive effectiveness. Furthermore, the second half of the cycle (from
$ t = 0.5T$
to
$ t = T$
) exhibits characteristics similar to those of the first half.
In the nonlinear case shown in the left-hand column of figure 5(b), the phase diagram exhibits a quasi-quadrilateral pattern that indicates that the swimmer sustains a high pitch angle while undergoing pronounced heave motion near
$ t = 0.25T$
and
$0.75T$
. This behaviour reflects effective phase tuning between pitch and heave motions, enabling the swimmer to maintain a thrust-conducive posture that promotes efficient forward projection of hydrodynamic pressure. Furthermore, the introduction of nonlinear force
$ F_{\textit{Non}}$
promotes larger heaving amplitudes (black curve) and pitch angles (green curve). These changes enhance vortex strength (Li et al. Reference Li, Feng, Wang and Liang2023) and increase hydrodynamic pressure. However, the swimmer with nonlinear tunable stiffness also faces a challenge: reduced propulsion due to the small pitching angle at
$ t = 0$
and
$0.5T$
.
Structural kinematics and flow evolution of the swimmer with (a) linear and (b) nonlinear tunable stiffness at
$f = 1 \, \text{Hz}$
. Left column: Phase diagrams of heave and pitch motion, with embedded pressure distribution snapshots at four representative time instants. Right column: Flow field visualizations showing velocity fields (grey arrows) at synchronized phases. All data collected after steady state (
${{t}_{s}}=10T$
).
The issue of thrust loss at small pitching angles is effectively answered through the following analysis. Figure 6 illustrates the time series of the pitch/heave motion, with time corresponding to the results in figure 5. In the linear case, the time interval between the maximum clockwise pitch (marked in green at
$ t = 0.05T$
) and the maximum anticlockwise pitch (at
$ t = 0.55T$
) spans half of the oscillation cycle. During this interval, the swimmer maintains a small pitching angle, which limits the effective conversion of hydrodynamic pressure into forward thrust, as discussed previously.
In contrast, the swimmer with nonlinear tunable stiffness exhibits a markedly faster pitch transition. The interval between maximum clockwise and anticlockwise pitch amplitudes is reduced to approximately one-quarter of a cycle (highlighted in blue from
$ t = 0.35T$
to
$ t = 0.6T$
). This rapid reversal minimises the duration spent at small pitch angles, thereby enhancing propulsion efficiency. A similar rapid response is observed in the heaving motion, as indicated by the steeper slope of the purple curve between
$ t = 0.4T$
and
$ t = 6T$
(highlighted in blue). Notably, the rapid downward heaving coincides with an extended period of large pitching angles (from
$ t = 0.6T$
to
$ t = 0.85T$
), representing the right flank of the quasi-quadrilateral swimming motion shown in figure 5(b). This synchronisation of heaving and pitching dynamics is closely associated with improved propulsive performance.
Time histories of the pitching/heaving motion of the swimmer with (a) linear stiffness and (b) nonlinear tunable stiffness at
$ f = 1 \, \text{Hz}$
.
The potential energy as a function of the heaving motion
$ y_{{A}}^{*}$
for a swimmer with nonlinear tunable stiffness. The cases
$ a^{*} = 1.5,\ b^{*} = 2$
,
$ a^{*} = 1.5,\ b^{*} = 1.5$
and
$ a^{*} = 1.8,\ b^{*} = 2$
correspond to the purple, gold and blue regions, respectively. The dashed lines represent the time-domain trajectories of
$ y_{{A}}^{*}$
and
$ U_{p}^{*}$
under the same external excitation conditions, given by
$ y_1(t) = h_0 \sin (2\pi ft)$
with
$ h_0 = 0.3$
and
$ f = 1 \, \text{Hz}$
. The snapshots within the purple region illustrate the hydrodynamic pressure distribution at representative instants, corresponding to the results presented in figures 5 and 6.
The exceptional dynamic behaviour fundamentally arises from the snap-through instability induced by nonlinear tunable stiffness (Overvelde et al. Reference Overvelde, Kloek, D’haen and Bertoldi2015), which can be understood from an energy perspective. Figure 7 illustrates the variation of potential energy
$ U_{p}^{*}$
throughout the heaving motion
$ y_{{A}}^{*}$
as given by (2.5), with dashed lines representing the time-domain trajectory of the leading edge. The figure delineates three distinct regions, each corresponding to different nonlinear stiffness parameters. It illustrates the existence of two stable potential wells (associated with equilibrium states, separated by an unstable equilibrium at
$ y_{{A}}^{*} = 0$
, forming a potential barrier defined as
$ \Delta U_{p}^{*}$
. The height of
$ \Delta U_{p}^{*}$
determines whether the potential energy barrier can be overcome.
Within the purple-highlighted region, the observed dynamic behaviour aligns with the results presented in figures 5 and 6. At
$ t = 0.35T$
, the swimmer is located near a potential well, where the potential energy reaches a minimum. As the swimmer progresses from
$ t = 0.35T$
to
$ t = 0.4T$
(highlighted by the purple dashed line), it approaches the maximum displacement
$ y_{{A}}^{*}$
, where potential energy peaks. Subsequently, from
$ t = 0.4T$
to
$ t = 0.6T$
, the swimmer rapidly descends due to the sudden release of stored potential energy, similar to a ball descending from a high point. This descent is further accelerated by hydrodynamic forces that generate additional downward motion. Consequently, the swimmer exhibits a rapid switching phenomenon in rotational direction between
$ t = 0.35T$
and
$ t = 0.6T$
, as illustrated in figure 6. This behaviour corresponds to the upper black curve in the quasi-quadrilateral phase diagram of figure 5, indicating that the swimmer avoids the limitation of small pitching angles that would fail to generate effective thrust. Further analysis reveals that between
$ t = 0.6T$
and
$ t = 0.85T$
, the swimmer experiences an upward hydrodynamic force that induces an anticlockwise moment on the leading edge, thereby prolonging the duration of large pitching angles. This process aligns with the right-hand-side green curve in the quasi-quadrilateral phase diagram of figure 5, demonstrating the swimmer’s ability to tune the phase relationship between heave and pitch motions. This phase tuning allows the swimmer to maintain a large pitching angle during the vertical motion, enhancing propulsion efficiency. Hence this reveals the high efficient propulsive mechanism based on nonlinear tunable stiffness.
3.3. Swimming characteristics altered by nonlinear parameter variations
Previous analysis has demonstrated that nonlinear tunable stiffness enhances swimming performance. However, its characteristics are governed by specific model parameters as given in (2.4), highlighting the need for investigating how parameter variations influence swimming behaviour. Hence the influence of parameter changes is further examined from an energetic perspective.
As illustrated by the gold and blue regions corresponding to
$ a^{*} = 1.5,\ b^{*} = 1.5$
and
$ a^{*} = 1.8,\ b^{*} = 2$
in figure 7, the potential barrier height increases with varying nonlinear parameters. As a result, the swimmer may only occasionally cross the potential barrier, or may be entirely confined within a single potential well.
For the parameters
$ a^{*} = 1.5$
and
$ b^{*} = 1.5$
(highlighted in gold), the potential barrier height increases compared to the previously discussed case (highlighted in purple). Consequently, the swimmer cannot traverse the barrier as freely as before. Instead, the system exhibits mixed-well oscillations characterised by inter-well oscillation from
$ t = 0$
to
$ t = 0.3T$
, and from
$ t = 0.4T$
to
$ t = 0.8T$
, along with intra-well oscillation from
$ t = 0.3T$
to
$ t = 0.4T$
, and from
$ t = 0.8T$
to
$ t = 0.9T$
. For the parameters
$ a^{*} = 1.8$
and
$ b^{*} = 2$
(highlighted in blue), the potential barrier further increases, preventing the swimmer from overcoming it. As a result, the swimmer exhibits periodic reciprocating motion within a single potential well, akin to a ball undergoing localised oscillations that cannot surpass the barrier.
These findings underscore the critical role of nonlinear parameters in determining a swimmer’s ability to traverse the potential barrier. The following analysis further explores the influence of varying nonlinear stiffness parameters on swimming behaviour and vortex field structures. Figure 8 illustrates the swimming behaviours alongside the corresponding vortex states corresponding to the gold and blue regions in figure 7, respectively.
In figure 8(a), the trajectory of
$ y_{{A}}^{*}$
exhibits mixed-well motion (both inter-well and intra-well mixing), consistent with figure 7. These oscillations reduce heaving amplitude, weakening vortex strength and hydrodynamic pressure. The wake structure resembles a superposition of vortex sheets from the linear and inter-well cases analysed in figure 5. In contrast, figure 8(b) depicts intra-well motion, where the swimmer oscillates within a potential well with weak amplitude and velocity. The resulting wake structure mirrors that of the linear case, featuring evenly spaced vortex sheets with alternating signs. However, the confined reciprocating motion produces weaker vortices, leading to significantly lower hydrodynamic pressure.
Time-domain kinematics and flow evolution o the swimmer with (a)
${a^{*}}=1.5,\ b^{*}=1.5$
and (b)
${a^{*}}=1.8,\ b^{*}=2$
at
$ f = 1 \, \text{Hz}$
. Left column: Time histories of pitch-heave motion with embedded pressure distribution snapshots at four representative time instants. Right column: Flow field visualizations showing velocity fields (grey arrows) at synchronized phases. All data collected after steady state (
${{t}_{s}}=10T$
).
Distinct vortex patterns significantly influence swimming performance, as shown in figure 9. For the mixed-well case (figure 9
a), the swimmer achieves speed
$3.8\ \text{BL}\ \text{s}^{-1}$
with efficiency 44 %, falling between the linear and inter-well cases analysed in figure 4. This aligns with the observed flow field, where vortex sheets from both regimes superimpose, contributing to intermediate propulsion efficiency. In contrast, the intra-well motion leads to a significant decline in performance, with the swimming speed and efficiency reducing to
$2.1\ \text{BL}\ \text{s}^{-1}$
and 11 %, respectively. This deterioration stems from the reduced oscillatory amplitude and the weakened vortex strength inherent to intra-well dynamics, which diminish propulsion effectiveness.
Instantaneous swimming performance for the swimmer with nonlinear parameters (a)
${a^{*}}=1.5, b^{*}=1.5$
and (b)
${a^{*}}=1.8,\ b^{*}=2$
at
$ f = 1 \, \text{Hz}$
. First column shows time-dependent speed
$|\dot {x}_{{C}}^{*}|$
; second column shows thrust
$F_{x}^{*}$
; third column shows circulation
$\Gamma$
; fourth column shows input power
$P_{in}^{*}$
and output power
$P_{out}^{*}$
. Solid lines represent the instantaneous values, and dashed lines represent time-period averages, with all the results obtained from the steady state starting at
$t_{s}=10T$
.
The analysis reveals that swimming performance is influenced by the characteristics of nonlinear tunable stiffness. Depending on whether the swimmer is able to overcome the potential energy barrier, three distinct vertical motion patterns – inter-well, mixed-well and intra-well motions – are observed. Each pattern leads to different swimming behaviours, with inter-well motion generally associated with superior propulsion performance. Notably, the nonlinear stiffness characteristics are governed by specific structural parameters in practical design. Therefore, selecting these parameters appropriately is essential for achieving optimal swimming efficiency.
4. Conclusions
This study presents a self-propelled physical model with parabolic tunable stiffness to emulate aquatic swimmers using the BCF propulsion mode. Additionally, the NVSM was developed to address nonlinear fluid–structure coupling. Numerical results show that the introduction of nonlinear tunable stiffness significantly improves swimming speed and propulsive efficiency, with the most notable enhancement observed at low actuation frequencies.
This improvement is attributed to the ability of stiffness modulation to tune the phase between heave and pitch motions, forming a quasi-quadrilateral pattern that enhances thrust. In addition, stiffness tuning intensifies vortex shedding and accelerates roll-up, further reducing energy expenditure. Energy based analysis reveals that the observed swimming behaviours are driven by snap-through instability, which enables the swimmer to traverse potential energy barriers and access high performance. Further analysis shows that variations in nonlinear stiffness reshape the potential energy landscape, giving rise to distinct dynamic regimes, including mixed-well oscillations and intra-well motions. As such, careful selection of these parameters is crucial for optimising swimming efficiency in practical design applications. These findings deepen the understanding of nonlinear mechanisms underlying propulsion enhancement, and offer theoretical guidance for the design and optimisation of high-performance, bioinspired robotic swimmers. However, the current model is based on a two-dimensional framework, which overlooks important three-dimensional effects. Additionally, although the model relies on numerical simulations and theoretical analysis, it has not yet been validated experimentally. These limitations highlight important areas for future work, including refining the model and enhancing its practical applicability.
Funding
This research work was supported by the National Natural Science Foundation of China (grant nos 12272128, 12425206 and 12272206), the Guangdong Basic and Applied Basic Research Foundation (grant no. 2024A1515010683) and the Natural Science Foundation of Jiangsu Province (grant no. BK20243019).
Declaration of interests
The authors declare no conflict of interest.
Appendix A. The NVSM
This appendix details the numerical procedure used to solve the nonlinear FSI model formulated in § 2. The algorithm involves time discretisation of the governing equations, and concurrent iterative solution for updating the vortex sheet and structural states from
$t_n$
to
$t_{n+1}$
.
The first step independently updates the free vortex sheet. At each time step, a new segment is shed from the swimmer’s trailing edge and appended to the existing wake, as shown in the orange-highlighted region in figure 10. At time
$t_n$
, the segment
$[\zeta _f^{j-1,n}, \zeta _f^{j,n}]$
corresponds to vorticity shed at
$t_j$
(no shedding at
$t_0$
). From
$t_n$
to
$t_{n+1}$
, the sheet is updated as follows: segments generated at or before
$t_{n-1}$
are advanced using a second-order Adams–Bashforth scheme, while the most recent endpoint
$\zeta _f^{n,n}$
is advanced using the Euler method:
\begin{equation} \begin{aligned} & \zeta _{f}^{j,n+1} = \zeta _{f}^{j,n} + 1.5\, \Delta t\, {D}_{t} \zeta _{f}^{j,n} - 0.5\, \Delta t\, {D}_{t} \zeta _{f}^{j,n-1}, \quad j = 0, \ldots , n-1, \\ & \zeta _{f}^{n,n+1} = \zeta _{f}^{n,n} + \Delta t\, {D}_{t} \zeta _{f}^{n,n}. \end{aligned} \end{equation}
The swimmer is discretised into
$m$
Chebyshev–Gauss–Lobatto points. Accordingly, there are
$m+6+1$
unknowns to solve from
${t}_{n}$
to
${t}_{n+1}$
, i.e.
$m$
bound vortex strengths
$({\bar {\gamma }_{b}^{(1,\ldots ,m)}})_{n+1}$
, six structural state variables
${({{x}_{{C}}},{{\dot {x}}_{{C}}},{{y}_{{C}}},{{\dot {y}}_{{C}}},\theta ,\dot {\theta })}_{n+1}$
, and vortex circulation
${(\Gamma )}_{n+1}$
. These unknowns are determined by enforcing the structural governing equations (2.1) and the boundary conditions described in § 2.3.
The no-penetration condition (2.8) at each Chebyshev point yields
$m$
constraints:
\begin{equation} \operatorname {Re}\big(\boldsymbol{n}\big(\overline {{{\partial }_{t}}\zeta ({{s}_{i}},t)}-\overline {\omega ({{s}_{i}},t)}\big)\big) = 0, \quad i = 1, \ldots , m. \end{equation}
The Kutta condition (2.9) provides an additional constraint at the trailing edge:
\begin{equation} \dot {\Gamma }(t) + (\mu (1,t) - \xi (1,t))\, \bar {\gamma }_{b}(1,t) = 0. \end{equation}
The remaining six constraints are obtained by discretising the structural dynamics (2.1) using the second-order Crank–Nicolson method (Crank & Nicolson Reference Crank and Nicolson1947):
\begin{align} & x_{{C}}^{n+1} - x_{{C}}^{n} - \frac {\Delta t}{2} (\dot {x}_{{C}}^{n+1} + \dot {x}_{{C}}^{n}) = 0, \notag\\[1.5pt] & \dot {x}_{{C}}^{n+1} - \dot {x}_{{C}}^{n} - \frac {\Delta t}{2m} ( - F^{n+1} \sin \theta ^{n+1} + D^{n+1} \cos \theta ^{n+1} - F^{n} \sin \theta ^{n} + D^{n} \cos \theta ^{n} ) = 0, \notag\\[1.5pt] & y_{{C}}^{n+1} - y_{{C}}^{n} - \frac {\Delta t}{2} (\dot {y}_{{C}}^{n+1} + \dot {y}_{{C}}^{n}) = 0, \notag\\[1.5pt] & \dot {y}_{{C}}^{n+1} - \dot {y}_{{C}}^{n} - \frac {\Delta t}{2m} ( F^{n+1} \cos \theta ^{n+1} + D^{n+1} \sin \theta ^{n+1} + F^{n} \cos \theta ^{n} + D^{n} \sin \theta ^{n} ) \notag\\ & \quad {}- F_{\textit{Non}}^{n+1} - F_{\textit{Non}}^{n} - k_{y} (y_{{C}}^{n+1} - \sin \theta ^{n+1} - y_{1}^{n+1}) - k_{y} (y_{{C}}^{n} - \sin \theta ^{n} - y_{1}^{n}) \notag\\ & \quad {}- c_{y} (\dot {y}_{{C}}^{n+1} - \cos \theta ^{n+1} \dot {\theta }^{n+1} - \dot {y}_{1}^{n+1}) - c_{y} (\dot {y}_{{C}}^{n} - \cos \theta ^{n} \dot {\theta }^{n} - \dot {y}_{1}^{n}) = 0, \notag\\ & \theta ^{n+1} - \theta ^{n} - \frac {\Delta t}{2} (\dot {\theta }^{n+1} + \dot {\theta }^{n}) = 0, \notag\\ & \dot {\theta }^{n+1} - \dot {\theta }^{n} - \frac {\Delta t}{2 I_{{A}}} \big ( M_{{A}}^{n+1} + M_{\textit{Inertia}}^{n+1} - c_{\theta } \dot {\theta }^{n+1} - k_{\theta } \theta ^{n+1} \notag\\ &\qquad {}+ M_{{A}}^{n} + M_{\textit{Inertia}}^{n} - c_{\theta } \dot {\theta }^{n} - k_{\theta } \theta ^{n} \big ) = 0. \end{align}
Consequently, as illustrated in the red dashed box of figure 10, the formulation defines a complete set of physical constraints governing the system’s evolution. Solving these nonlinear algebraic system via Broyden’s method (Broyden Reference Broyden1965) yields the target state at the next time step.
A flowchart of the NVSM solving the nonlinear fluid–structure interaction problems.
To ensure reliability, the NVSM is validated by comparing its results with those obtained from a conventional vortex sheet method using a simple swimming model. Specifically, a pure pitching model is included as a baseline (Heydari & Kanso Reference Heydari and Kanso2021), where the swimmer experiences active sinusoidal pitch motion
${{\theta }_{a}}(t)={{\theta }_{0}}\sin (2\pi ft)$
around the leading edge (where
${\theta }_{a}$
denotes the angle relative to the swimming direction,
${\theta }_{0}$
is the flapping amplitude, and
$f$
is the flapping frequency). As shown in figure 11, the result of the NVSM agrees well with the traditional vortex sheet method, demonstrating the accuracy and reliability of the NVSM for nonlinear fluid–structure interactions. This consistency across critical performance metrics reinforces the robustness of the NVSM in addressing nonlinear coupling effects, thereby enhancing confidence in its application to complex swimming models.
Comparison of the NVSM and conventional vortex sheet method results for a swimmer undergoing active sinusoidal pitching motion.
