1. Introduction
Bio-inspired underwater vehicles draw inspiration from the morphology and kinematics of fish to exploit natural fluid–structure interaction mechanisms for efficient propulsion (Guo, Fukuda & Asaka Reference Guo, Fukuda and Asaka2003; Zhu et al. Reference Zhu, White, Wainwright, Di Santo, Lauder and Bart-Smith2019). Unlike conventional marine propulsion systems based on rigid propellers, biomimetic approaches rely on oscillatory and deformable motions, offering advantages such as improved energy efficiency, enhanced manoeuvrability and reduced acoustic signatures (Sun et al. Reference Sun2022; Zhu et al. Reference Zhu, Kang, Tong, Ma, Tian and Fan2025a ). These characteristics make them particularly attractive for applications in environmental monitoring, marine exploration and search-and-rescue operations.
Among the diverse swimming strategies observed in nature, pitching-foil-based propulsion has emerged as a canonical reduced-order model for fish locomotion (Sfakiotakis, Lane & Davies Reference Sfakiotakis, Lane and Davies1999; Rosenberger Reference Rosenberger2001). Many batoid species, such as cownose rays, generate thrust through the oscillatory pitching motion of broad pectoral fins (Zhang et al. Reference Zhang, Huang, Pan, Yang and Huang2021; Xu et al. Reference Xu, Luo, Hou, Huang, Cao and Pan2023), a mechanism that has been widely adopted in the design and analysis of bio-inspired underwater vehicles (Niu et al. Reference Niu, Zhang, Bi and Cai2013; Liu et al. Reference Liu, Chen, Guo, Su, Li, Zeng, Fan and Cui2024). Extensive studies have shown that foil kinematics – particularly oscillation frequency and amplitude – play a central role in determining propulsion performance. Optimal thrust production is typically achieved within a narrow Strouhal number range, commonly reported as
$St = 0.25$
–
$0.35$
(Triantafyllou, Triantafyllou & Grosenbaugh Reference Triantafyllou, Triantafyllou and Grosenbaugh1993), a finding consistent across aquatic and aerial locomotion, including the pitching flight of birds, bats and insects (Taylor, Nudds & Thomas Reference Taylor, Nudds and Thomas2003).
As research has progressed, attention has increasingly shifted from kinematics alone to the associated wake dynamics, which provides the physical link between foil motion and force generation. Experimental measurements using particle image velocimetry (Godoy-Diana, Aider & Wesfreid Reference Godoy-Diana, Aider and Wesfreid2008), force transducers (Mackowski & Williamson Reference Mackowski and Williamson2015) and numerical simulations (Alam & Muhammad Reference Alam and Muhammad2020) have revealed that propulsion efficiency is closely tied to wake topology. A hallmark transition from a drag-producing Kármán vortex street to a thrust-generating reverse von Kármán wake signifies a fundamental change in momentum transfer. More detailed phase-space analyses further demonstrate that wake structure, force production and efficiency depend jointly on Strouhal number and oscillation amplitude, with multiple wake regimes coexisting across parameter space (Schnipper, Andersen & Bohr Reference Schnipper, Andersen and Bohr2009; Andersen et al. Reference Andersen, Bohr, Schnipper and Walther2017).
While these studies have substantially advanced our understanding of pitching-foil propulsion in free-stream conditions, bio-inspired swimmers in nature rarely operate far from boundaries. In near-bed or near-surface environments, ground effect plays a crucial role in shaping hydrodynamic performance. Biological systems such as birds (Blake Reference Blake1983), flying fish (Park & Choi Reference Park and Choi2010) and rays swimming close to the seabed (Hou et al. Reference Hou, Luo, Bao, Huang, Xing, Xu and Pan2024) are known to exploit boundary proximity, while engineered systems – including aircraft during takeoff and landing – exhibit analogous behaviour (Zheng & Ash Reference Zheng and Ash1996; Rozhdestvensky Reference Rozhdestvensky2006). In general, ground effect modifies pressure distributions, suppresses tip vortices and enhances lift-to-drag ratios (Luo & Chen Reference Luo and Chen2012; Babu et al. Reference Babu, Madhesh, Hasan and Mukesh2025).
Recent work has extended these concepts to unsteady, oscillatory propulsion near boundaries. Periodically pitching swimmers operating close to the ground have been shown to experience enhanced swimming speed (Fernandez-Prats et al. Reference Fernández-Prats, Raspa, Thiria, Huera-Huarte and Godoy-Diana2015), increased thrust without a proportional loss of efficiency (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014; Mivehchi, Dahl & Licht Reference Mivehchi, Dahl and Licht2016; Ryu et al. Reference Ryu, Park, Kim and Sung2016) and reduced power consumption at constant thrust (Sierra & Guo Reference Sierra and Guo2018). However, these benefits are not universal; under certain conditions, power requirements may increase (Blevins & Lauder Reference Blevins and Lauder2013). An equilibrium height has been identified at which propulsion performance is maximised, with deviations generating restorative lift forces that passively stabilise the swimmer (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014; Mivehchi et al. Reference Mivehchi, Dahl and Licht2016; Perkins et al. Reference Perkins, Elles, Badlissi, Mivehchi, Dahl and Licht2018; Kurt et al. Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019). Such performance variations are fundamentally governed by wake–boundary interactions. As wakes convect downstream, boundary-layer separation near the ground induces vortex shedding and deflection, forming single-row or distorted vortex structures (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014; Sierra & Guo Reference Sierra and Guo2018; Wakim et al. Reference Wakim, Brion, Dolfi-Bouteyre and Jacquin2020; Zhong et al. Reference Zhong, Han, Moored and Quinn2021). Control-volume-based scaling laws incorporating added-mass and circulatory effects have been proposed to rationalise these wake deflections (Mivehchi et al. Reference Mivehchi, Zhong, Kurt, Quinn and Moored2021; Han et al. Reference Han, Qiang, Amin, D.B. and Moored2023).
Despite this growing body of work, most existing studies idealise the boundary as rigid and impermeable, neglecting the complex nature of real seabeds. In many near-shore and aquaculture environments, the bed contains substantial amounts of synthetic and organic particles, including microplastics, fish faeces and uneaten feed pellets (Cromey et al. Reference Cromey, Nickell, Treasurer, Black and Inall2009; Kowalski, Reichardt & Waniek Reference Kowalski, Reichardt and Waniek2016; Auta, Emenike & Fauziah Reference Auta, Emenike and Fauziah2017; Hargrave et al. Reference Hargrave, Filgueira, Grant and Law2022). These particles are readily mobilised and resuspended by near-bed flows, influencing sediment transport, water quality and benthic ecosystems. While numerous studies have examined how boundaries affect propulsion performance, far fewer have considered the reverse influence – namely, how oscillating propulsors actively modify particle-laden and deformable beds.
In natural systems, this reverse interaction is evident. Rays disturb sediments through pectoral fin motion to facilitate camouflage, prey exposure and predator avoidance (Corn et al. Reference Corn, Farina, Summers and Gibb2018; Bo et al. Reference Bo, Wang, Xing, Qian, Song and Zheng2023; Ren et al. Reference Ren, Chen, Tian, Yu, Li and Tan2024), and may further fluidise the bed through suction generated by their mouths or gill slits (Gregory et al. Reference Gregory, Ballance, Gibson and Ayling1979; Wilga et al. Reference Wilga, Maia, Nauwelaerts and Lauder2012). Similar phenomena occur in engineered systems, where rotating propellers in sediment-laden environments significantly alter particle suspension and resettlement patterns (Tan & Yuksel Reference Tan and Yuksel2018; Cihan et al. Reference Cihan, Dog, Yilmaz, Ozan, Yildiz and Sahin2022). These examples highlight the inherently bidirectional nature of structure–environment interactions.
Such bidirectional coupling is a defining feature of fluid–structure–environment systems. Classical examples include scour induced by horseshoe vortices and wake shedding around marine piles (Sarpkaya Reference Sarpkaya1986; Obasaju, Bearman & Graham Reference Obasaju, Bearman and Graham1988; Sumer, Christiansen & Fredsoe Reference Sumer, Christiansen and Fredsoe1997), as well as feedback from evolving sediment beds that modify boundary layers and structural loading (Diaz-Carrasco et al. Reference Díaz-Carrasco, Croquer, Tamimi, Lacey and Poncet2021; Zhang et al. Reference Zhang, Zheng, Wang and Liang2022). Nevertheless, for oscillatory propulsion systems – particularly pitching foils operating above deformable, particle-laden beds – the dynamic feedback among fluid motion, structural kinematics and particle transport remains poorly understood. Although preliminary observations of sediment interaction effects on pitching propulsion have been reported (Aju et al. Reference Aju, Gong, Pham, Kaushik and Jin2022), a systematic, mechanistic understanding of this triphasic coupling is still lacking.
The present study aims to address this gap through numerical investigation of pitching-foil propulsion over particle-laden beds. We first characterise the propulsive performance and sediment transport response induced by the oscillating foil, and subsequently elucidate the underlying flow physics responsible for these behaviours. By explicitly linking foil kinematics, wake evolution and particle motion, this work provides a mechanistic framework for understanding coupled structure–fluid–particle interactions in near-bed, particle-laden environments (figure 1).
The remainder of the paper is organised as follows: § 2 describes the numerical methodology; § 3 characterises the mechanical response of the pitching foil; § 4 examines particle geometry and kinematics; § 5 analyses flow-field modifications and the wake dynamics; and § 6 summarises the key findings and discusses their broader implications.
The coupled structure–fluid–particle interactions.
2. Materials and methods
2.1. Model description
In this study, a NACA 0012 hydrofoil with a chord length of
$c$
= 0.1 m is used to investigate its dynamic behaviour. The hydrofoil moves at a free-stream velocity of
$U_\infty$
= 0.1 m s–1 in the negative
$x$
-direction, corresponding to a Reynolds number of
$\textit{Re} = ({\rho U_\infty c}/{\mu }) = 11\,000$
, where
$\rho$
is the fluid density and
$\mu$
is the dynamic viscosity. Simultaneously, the hydrofoil undergoes sinusoidal pitching motion about a pivot point located at 25
$\%$
of the chord length from the leading edge (figure 2
a). The pitching motion is described by the following equation:
\begin{equation} \theta (t) = \theta _0 \sin (2 \pi \textit{ft}), \end{equation}
where
$\theta _0$
is the pitching amplitude,
$f$
is the frequency and
$t$
is time.
Schematics of (a) the foil kinematics, (b) the computational domain and (c) the meshing strategy. To capture sufficient oscillation cycles of the foil at
$\textit{St}_c$
= 0.4, the length of the computational domain is extended to 40
$c$
.
The pitching motion of the hydrofoil is characterised by two key dimensionless parameters: the Strouhal number based on chord length,
$ \textit{St}_c = \textit{fc}/{U_\infty }$
, and the pitching amplitude,
$ \theta _0$
. In this study, measurements were conducted for five values of
$ \textit{St}_c$
, ranging from 0.4 to 1.2, and nine values of
$ \theta _0$
, spanning from
$ 8^\circ$
to
$ 24^\circ$
. The corresponding peak-to-peak trailing-edge displacement due to pitching is given by
$A = 2 \times 0.75c \sin \theta _0$
, where the factor 0.75 accounts for the distance from the pitching axis (located at 25 % chord) to the trailing edge. This geometric consideration leads to an amplitude-based Strouhal number,
$\textit{St}_\theta = \textit{fA}/{U_\infty }$
, which in this study ranges approximately from 0.0835 to 0.7321.This range reflects the natural kinematics of batfish pectoral-fin locomotion and encompasses regimes associated with efficient thrust production (Fish et al. Reference Fish, Dong, Zhu and Bart-Smith2017; Zhu et al. Reference Zhu, Liu, Han, Feng, Moored, Zhong and Quinn2025b
; Rosenberger Reference Rosenberger2001).
This study investigates the differences in propulsion performance and the wake dynamics of a pitching foil operating above a sediment bed compared with a rigid ground surface. The foil is positioned at a vertical distance
$ d$
from either the sediment bed or the solid boundary, with a non-dimensional spacing defined as
$ d^* = d / c$
. Three values of
$ d^*$
are considered: 0.4, 0.5 and 0.6. The sediment bed has a thickness of 20 mm and is composed of spherical particles with a diameter of
$ d_{\!p} = 2 \, \text{mm}$
and a density of
$ \rho _{\!p} = 1190 \, \rm{kg\,m}^{-3}$
. The selected particles represent synthetic and organic materials, such as microplastics and aquaculture-derived solids in near-shore environments.
The instantaneous thrust force
$ F_t(t) = -F_x(t)$
and lift force
$ F_l(t) = F_y(t)$
on the foil are measured during the simulations and normalised as thrust and lift coefficients (Zhong et al. Reference Zhong, Han, Moored and Quinn2021)
\begin{equation} C_t(t) = \frac {F_t(t)}{\frac {1}{2} \rho c U_\infty ^2}, \end{equation}
\begin{equation} C_l(t) = \frac {F_l(t)}{\frac {1}{2} \rho c U_\infty ^2}. \end{equation}
Thus, the time-averaged thrust and lift coefficients are define as
$ C_{\!T} = ({1}/{T})\int _0^T C_t \, {\rm d}t$
and
$ C_{\!L} = ({1}/{T})\int _0^T C_l \, {\rm d}t$
.
The instantaneous hydrodynamic power
$P(t)$
is calculated based on
$ P(t) = -F_t(t) U_\infty - M(t) \dot {\theta }(t)$
, where
$M(t)$
is the instantaneous torque about the pitching pivot point (Das, Shukla & Govardhan Reference Das, Shukla and Govardhan2021). The power coefficient is calculated as
\begin{equation} C_{\!p}(t) = \frac {P(t)}{\frac {1}{2} \rho c U_\infty ^3}. \end{equation}
Thus, the time-averaged power is given by
$ C_{\!P} = ({1}/{T})\int _0^T C_{\!p} \, {\rm d}t$
, and the propulsive efficiency
$\eta$
is obtained as
\begin{equation} \eta = \frac {C_{\!T}}{C_{\!P}}. \end{equation}
2.2. Numerical details
This study employs a two-way coupled computational fluid dynamics and discrete element method (CFD–DEM) framework to simulate the interactions between the fluid and solid phases (Balachandar & Eaton Reference Balachandar and Eaton2010; Zhang, Gutierrez & Li Reference Zhang, Gutierrez and Li2017). In this approach, the fluid phase is treated as a continuum using an Eulerian formulation within the CFD module, while the solid phase is represented as a discrete assembly of particles modelled via a Lagrangian approach in the DEM module. Coupling between the two phases is realised through interphase momentum exchange source terms that account for drag and other interaction forces.
Turbulence in the fluid phase is modelled using the realisable
$ k\text{-}\epsilon$
(
$k$
is the turbulent kinetic energy and
$\epsilon$
is the turbulent dissipation rate) model within the Reynolds-averaged Navier–Stokes (RANS) framework. This model has been shown to perform reliably for unsteady flows at moderate Reynolds numbers involving separation and rotation (Han, Sagaut & Lucor Reference Han, Sagaut and Lucor2012; Guo et al. Reference Guo, Hu, Guo, Zhang and Lin2020). The governing equations can be written as follows (Zhou et al. Reference Zhou, Kuang, Chu and Yu2010; Lu et al. Reference Lu, Liu, Xie, Sun, Yang and Guo2024):
\begin{equation} \frac {\partial (\varphi _{\!f}) }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\varphi _{\!f} \boldsymbol{u}) = 0, \end{equation}
\begin{equation} \frac {\partial (\varphi _{\!f}\boldsymbol{u}) }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\varphi _{\!f}\boldsymbol{u} \boldsymbol{u}) = -\frac {1}{\rho } \boldsymbol{\nabla }\!p + \nu \boldsymbol{\nabla} ^2 \boldsymbol{u} + \varphi _{\!f}\boldsymbol{g} + \boldsymbol{f}_{\!\textit{fp}}, \end{equation}
Here,
$\varphi _{\!f}$
is the volume fraction of the continuous phase,
$\boldsymbol{u}$
represents the fluid velocity vector and
$\boldsymbol{f}_{\!\textit{fp}}$
represents the force exerted by particles per unit volume on fluid.
For the solid phase, the motion of each individual particle is governed by Newton’s second law (Zhou et al. Reference Zhou, Kuang, Chu and Yu2010; Guo & Curtis Reference Guo and Curtis2015)
\begin{equation} m_{p,i} \frac {{\rm d} \boldsymbol{v}_{p,i}}{{\rm d}t} = \boldsymbol{F}_{\!pf,i} + \sum _{j=1}^{n} \boldsymbol{F}_{\!p,\textit{ij}} + \boldsymbol{F}_{\!g,i}, \end{equation}
\begin{equation} I_{p,i} \frac {{\rm d} \boldsymbol{\omega }_{p,i}}{{\rm d}t} = \boldsymbol{M}_{\!pf,i} + \sum _{j=1}^{n} \boldsymbol{M}_{\!p,\textit{ij}}, \end{equation}
where
$m_{p,i}$
and
$I_{p,i}$
are, respectively, the mass and rotational inertia of particle
$i$
, and
$\boldsymbol{v}_{p,i}$
and
$\boldsymbol{\omega }_{p,i}$
represent the translational and angular velocities.
The forces involved include the particle–fluid interaction force
$\boldsymbol{F}_{\!pf,i}$
, inter-particle forces
$\boldsymbol{F}_{\!p,\textit{ij}}$
and the gravitational force
$\boldsymbol{F}_{\!g,i}$
. Among these, the particle–fluid interaction force,
$\boldsymbol{F}_{\!pf,i}$
, is the cumulative effect of various forces exerted by the fluid on the particle. These include the drag force
$\boldsymbol{F}_d$
, the pressure gradient force
$\boldsymbol{F}_{\!\boldsymbol{\nabla }\!p}$
and the lift forces, such as the Sommerfeld force
$\boldsymbol{F}_{\!L,som}$
. And the buoyancy force is included in
$\boldsymbol{F}_{\!\boldsymbol{\nabla }\!p}$
. The inter-particle forces
$\boldsymbol{F}_{\!p,\textit{ij}}$
consist of the normal and tangential components,
$\boldsymbol{F}_{\!n,\textit{ij}}$
and
$\boldsymbol{F}_{\!t,\textit{ij}}$
, which are based on Hertz-–Mindlin no-slip contact theory (Zhao, Evans & Zhou Reference Zhao, Evans and Zhou2018; Harzanagh, Orhan & Ergun Reference Harzanagh, Orhan and Ergun2018). Specific properties are as follows: Young’s modulus
$E = 1.0 \times 10^{7}\,\mathrm{Pa}$
, Poisson’s ratio
$\nu = 0.35$
, static friction coefficient
$\mu _s = 0.61$
, normal restitution coefficient
$e_n = 0.5$
and tangential restitution coefficient
$e_t = 0.5$
.
The torque acting on particle
$i$
is primarily governed by inter-particle interactions
$\boldsymbol{M}_{\!p,\textit{ij}}$
, which can be decomposed into two dominant components: the rolling resistance torque
$\boldsymbol{M}_{\!r,ij}$
and other contact torques
$\sum \boldsymbol{M}_{\!c,ij}$
. The torque exerted by the fluid
$\boldsymbol{M}_{\!pf,i}$
is sufficiently small and can be neglected in comparison.
Figure 2(b) illustrates the configuration of the computational domain and the grid system employed in the simulations. The domain spans a total length and height of
$ 25c \times{} (8c + d)$
, where
$ c$
is the chord length of the hydrofoil and
$ d$
denotes the distance between the foil and the lower boundary. Symmetry boundary conditions are applied at the left, right and upper domain boundaries, while no-slip boundary conditions are imposed on both the hydrofoil surface and the bottom boundary.
To simulate the hydrofoil’s pitching motion, an overset mesh technique is adopted. As shown in figure 2(c), the domain is partitioned into a stationary background and a moving overset subdomain that encloses the hydrofoil. The overset region translates and rotates relative to the fixed background, with data exchange achieved through interpolation across overlapping mesh zones. This approach allows accurate resolution of moving boundaries without deforming the computational mesh.
The grids are locally refined in both the ground-adjacent region and the downstream wake region to capture fine-scale flow structures and the particle dynamics. Additionally, the surface meshes along the foil and the bottom boundary are refined to ensure wall-resolved simulations, with the first grid layer satisfying
$ y^+ \leqslant 1$
.
The base grid size is set to
$ 1.5 \times 10^{-2}\, \text{m}$
, and the uniform time step is
$ \Delta t = 5.0 \times 10^{-4}\, \text{s}$
. Details regarding grid independence and time-step convergence are provided in Appendix A. The CFD–DEM framework is further validated by comparing simulated particle patterns with experimental observations reported by Aju et al. (Reference Aju, Gong, Pham, Kaushik and Jin2022), as detailed in Appendix A. The resulting ‘Y’-shaped particle cloud is consistent with the experimental findings.
Simulations are performed for two-dimensional, incompressible and unsteady flows. The fluid phase is solved using a finite-volume-based segregated flow solver. Pressure–velocity coupling within this framework is handled using the semi-implicit method for pressure-linked equations algorithm. Spatial discretisation employs a second-order upwind scheme for convective fluxes and a second-order central differencing scheme for diffusive terms. Time advancement is performed using an implicit unsteady formulation with a second-order temporal discretisation (Guo et al. Reference Guo, Hu, Guo, Zhang and Lin2020). Particle positions and velocities are updated once per fluid time step using a first-order time integration scheme. To enhance numerical stability during strong particle–fluid interactions, particle tracking is sub-stepped, with a maximum of 100 sub-steps per fluid time step. The simulations are performed on a workstation equipped with 16 CPU cores. With a fixed time step of
$ \Delta t = 5.0 \times 10^{-4}\, \text{s}$
, simulating 1
$ \text{s}$
of physical foil motion requires approximately 0.5 hours of wall-clock computation time.
3. Propulsive performance near the sandy ground
3.1. Comparison of particle and particle-free cases
To assess the impact of sediment particles on propulsive performance, two representative cases are examined at Strouhal numbers
$ \textit{St}_c = 0.4$
and
$ \textit{St}_c = 1.0$
. For each case, three configurations are considered: (i) a pitching foil operating above a sediment bed (with particles), (ii) a pitching foil near a rigid, particle-free boundary (no particles, near ground) and (iii) a pitching foil far from any boundary (no particles, far ground). The time-averaged thrust coefficient
$ C_{\!T}$
, lift coefficient
$ C_{\!L}$
and propulsive efficiency
$ \eta$
are evaluated and compared across all configurations, as illustrated in figure 3.
At the lower Strouhal number (
$ \textit{St}_c = 0.4$
), the foil primarily behaves as a drag-producing body, as evidenced by negative values of the thrust coefficient
$ C_{\!T}$
across all configurations (figure 3
a). Among the three, the sediment case exhibits the greatest drag, indicating that particle–fluid interactions increase resistance to forward motion. In contrast, at the higher Strouhal number (
$ \textit{St}_c = 1.0$
), the foil transitions to a thrust-producing regime, with all configurations showing positive
$ C_{\!T}$
values. Nevertheless, the thrust generated in the presence of sediment remains consistently lower than that in the clean-ground and far-ground cases, suggesting that sediment presence suppresses thrust production through damping effects. This thrust suppression is attributed to a reduction in jet intensity, driven by interphase momentum exchange as particles extract energy from the flow (see § 5.2 for a detailed discussion based on streamline velocity analysis).
Comparison of propulsive performance across three configurations: (a) time-averaged thrust coefficient
$C_{\!T}$
(
$C_{\!T} \lt 0$
indicates net drag;
$C_{\!T} \gt 0$
indicates net thrust), (b) time-averaged lift coefficient
$C_{\!L}$
(
$C_{\!L} \lt 0$
indicates downward force toward the ground;
$C_{\!L} \gt 0$
indicates upward force) and (c) propulsive efficiency
$\eta$
. Red lines represent the sediment case, blue lines represent the clean-ground case and grey lines represent the far-ground case. Solid lines correspond to
$\textit{St}_c = 0.4$
, and dashed lines correspond to
$\textit{St}_c = 1.0$
.
The lift behaviour (figure 3
b) also reveals distinct modifications due to the presence of sediment. At
$ \textit{St}_c = 0.4$
, the foil produces a small positive lift in all configurations, resulting in a net force away from the boundary. However, the lift magnitude is reduced in the sediment case, suggesting a diminished upward momentum transfer due to particle–fluid interactions. At
$ \textit{St}_c = 1.0$
, a more complex trend emerges. At lower pitch amplitudes, both the sediment and clean-ground cases exhibit negative lift, indicating attraction toward the bed. With increasing amplitude, the clean-ground case undergoes a sharp reversal in
$ C_{\!L}$
, transitioning to strong upward lift. In contrast, the lift in the sediment case remains negative across all amplitudes, highlighting the suppressive influence of mobile particles on wake-induced lift generation. This reversal in
$ C_{\!L}$
stems from particles altering the surface pressure distribution (see § 5.3 for details).
The corresponding propulsive efficiency
$ \eta$
, shown in figure 3(c), exhibits similar trends. At
$ \textit{St}_c = 0.4$
, the system exhibits negative efficiency, consistent with net drag production. Notably, the sediment case yields a more negative
$ \eta$
, reflecting both enhanced drag and greater potential for energy dissipation or harvesting. At
$ \textit{St}_c = 1.0$
, the system transitions into a thrust-generating regime, and
$ \eta$
becomes positive. Nevertheless, efficiency remains substantially lower in the particle-laden configuration compared with the particle-free cases, indicating that sediment presence not only reduces thrust but also increases power consumption, thereby diminishing overall propulsive performance.
These findings collectively demonstrate that sediment particles alter both the hydrodynamic force and energy transfer mechanisms in pitching-based propulsion. Particle-induced damping and flow disruption suppress thrust and lift while degrading efficiency. The results underscore the necessity of accounting for the sediment dynamics in near-bed propulsion studies, particularly in biomimetic or environmentally relevant contexts.
Contour plots of (a)
$C_{\!T}$
, (b)
$C_{\!L}$
and (c)
$\eta$
at different distances
$d^*$
(0.4, 0.5 and 0.6). The black curves represent the iso-lines of the coefficients, and the solid curves indicate the equilibrium position.
3.2. Overall responses of particle effect
The overall propulsive behaviour of the pitching foil in the presence of sediment particles is illustrated in figure 4. Across all examined distances from the bed (
$ d^*$
), the thrust coefficient (
$ C_{\!T}$
), lift coefficient (
$ C_{\!L}$
) and propulsive efficiency (
$ \eta$
) exhibit consistent trends as functions of the Strouhal number (
$ \textit{St}_c$
) and pitch amplitude (
$ \theta _0$
). As
$ \textit{St}_c$
increases, the foil transitions from a drag-producing regime to thrust generation (figure 4
a), and from positive lift (pushing the foil away from the bed) to negative lift (drawing it toward the bed) (figure 4
b). The black contour line in each panel indicates the equilibrium point where net thrust or lift vanishes.
Effect of particles and foil kinematics on the time-resolved (a)
$C_t$
and (b)
$C_l$
. The solid curves represent the coefficient in the particle cases, while the dashed curves represent the particle-free cases. The red, blue and grey curves correspond to
$\theta _0$
=
$8^\circ$
,
$16^\circ$
and
$24^\circ$
, respectively. Here,
$d^*$
= 0.6.
Pitching amplitude
$ \theta _0$
plays a critical role in the propulsion dynamics. Larger amplitudes lead to greater thrust and lift magnitudes and cause earlier onset of thrust production – that is, the critical
$ \textit{St}_c$
at which net thrust becomes positive decreases with increasing
$ \theta _0$
. As thrust generally increases with both
$ \textit{St}_c$
and
$ \theta _0$
, the propulsive efficiency
$ \eta$
exhibits a similar trend. At higher
$ \textit{St}_c$
, where thrust is generated,
$ \eta$
also become positive – suggesting that energy is transferred to the streamwise flow and converted effectively into propulsion. This behaviour indicates an efficient locomotion regime (figure 4
c).
To further examine the influence of sediment particles, time-resolved force responses are analysed for two representative cases at
$ \textit{St}_c = 0.4$
and
$ \textit{St}_c = 1.0$
(figure 5). At
$ \textit{St}_c = 0.4$
, the foil experiences net drag and positive lift, while at
$ \textit{St}_c = 1.0$
, net thrust and typically negative lift are observed – although positive lift may still occur at higher amplitudes in the particle-free case.
For both
$ \textit{St}_c$
values, the time-resolved profiles of
$ C_t$
and
$ C_l$
follow predictable patterns:
$ C_t$
exhibits a primary peak during each half-oscillation, whereas
$ C_l$
oscillates once per full cycle. Increasing amplitude amplifies the peak magnitudes of both coefficients. However, the presence of sediment consistently suppresses
$ C_t$
throughout the oscillation cycle. At
$ \textit{St}_c = 1.0$
, a pronounced asymmetry emerges in the particle-free case – where a stronger positive
$ C_t$
peak (around
$ t^* = 0.37$
) and a weaker negative peak (around
$ t^* = 0.57$
) result in enhanced net thrust. This asymmetric thrust arises from the deflection of the wake away from the ground in the absence of particles, caused by complex, unstable interactions between the vortices and ground-induced secondary shear layer. In contrast, this asymmetry disappears in the particle case, where
$ C_t$
remains symmetric, suggesting that particle interactions damp the unsteady mechanisms responsible for peak thrust. This indicates that particle–flow interactions act as a momentum dissipator. The resulting attenuation of wake momentum not only reduces the thrust but also restores a symmetric temporal force distribution.
The effects on
$ C_l$
are more nuanced. At
$ \textit{St}_c = 0.4$
, both positive and negative lift peaks are slightly larger in the sediment case. At
$ \textit{St}_c = 1.0$
, however, this configuration shows a reduced positive lift peak and an enhanced negative peak, reinforcing the net downward force trend. This shift occurs because particles alter the surface pressure distribution, specifically by increasing pressure on the upper surface during the upstroke, and suppressing pressure recovery on the lower surface during the downstroke (detailed in § 5.3). Conversely, in the particle-free case with large
$ \theta _0$
, the positive
$ C_l$
peak becomes significantly amplified and occurs earlier (around
$ t^* = 0.46$
), leading to a time-averaged lift coefficient
$ C_{\!L}$
reversal from negative to positive. This phenomenon is likely driven by more pronounced flow separation associated with high-amplitude pitching, which is evidently suppressed by the particle-laden flow.
4. Particle motions and sediment transport patterns
The sedimentary response to foil pitching varies significantly with both Strouhal number (
$ \textit{St}_c$
) and pitch amplitude (
$ \theta _0$
), and can be categorised into three distinct geometric evolution patterns (figure 6
a): (i) no visible change in bed geometry, (ii) formation of localised peaks and valleys and (iii) development of a suspended particle cloud above the bed. In the third pattern, particles initially lift off in a clockwise trajectory relative to the foil motion, followed by the formation of a secondary, counterclockwise-rotating structure, ultimately forming a characteristic ‘Y’-shaped distribution.
To map the dependence of these patterns on
$ \textit{St}_c$
and
$ \theta _0$
, figure 6(b) presents a regime classification diagram based on observed sediment responses. Dashed contour lines of constant
$ \textit{St}_\theta$
are superimposed to delineate transition boundaries. Results show that a more complex sediment response is triggered as
$ \textit{St}_\theta$
increases. Furthermore, the critical
$ \textit{St}_\theta$
required increases with foil–bed separation distance
$ d^*$
. For example, the threshold for the onset of suspended particle clouds (pattern 3) rises from 0.33 to 0.48 as
$ d^*$
increases from 0.4 to 0.6.
Patterns for the sediment bed: (a) schematic; (b) pattern map, (
$\triangle$
) pattern 1, stationary beds, (
) pattern 2, local peaks/valleys, and (
$\bullet$
) pattern 3, suspended `Y’-shaped particle flows. The dashed black curves represent the pattern transition thresholds of
$\textit{St}_\theta$
.
Contour plots of the number of suspended particles: (a) leaving the sediment bed and (b) exceeding the equilibrium position. Particle counts are recorded from the foil leading edge to two chord lengths downstream of the trailing edge and then averaged over one cycle.
To elucidate the mechanisms underlying the observed sedimentary transformations, particle motion is further examined. Two indicators are introduced:
$ n_1$
, the number of particles detaching from the bed, and
$ n_2$
, the number of particles surpassing the equilibrium height of the foil. Both quantities increase with rising Strouhal number
$ \textit{St}_c$
and pitch amplitude
$ \theta _0$
, and decrease with increasing foil–bed distance
$ d^*$
(figure 7). The attenuation in particle activity at larger
$ d^*$
values is attributed to two primary factors: (i) reduced flow-induced disturbance at greater separations, and (ii) increased energetic requirements for particles to reach and exceed the foil’s equilibrium height.
Figure 8 presents the mean particle velocity
$v_{\!p}^*$
associated with each of the three identified sedimentary patterns. As expected, higher particle velocities are observed at elevated
$ \textit{St}_c$
, larger
$ \theta _0$
and smaller
$ d^*$
. In pattern 1, velocities remain near zero, reflecting negligible fluid–particle interaction. In patterns 2 and 3, velocity profiles exhibit dynamic evolution: a primary peak occurs around
$ t^* = 0.2$
, followed by a decline and a secondary peak near
$ t^* = 0.6$
. This bimodal signature is particularly pronounced in pattern 3 and indicates sustained interactions between the unsteady flow and suspended particles throughout the oscillation cycle.
Time-resolved mean particle velocity at various (a)
$\textit{St}_c$
, (b)
$\theta _0$
and (c)
$d^*$
. The mean particle velocity is calculated by averaging the velocity of particles from the foil leading edge to two chord lengths downstream of the trailing edge. Panel (a) shows
$\theta _0 = 16^\circ$
,
$d^*$
= 0.5. The red (
$\textit{St}_c$
= 0.8), blue (
$\textit{St}_c$
= 1.0) and grey (
$\textit{St}_c$
= 1.2) curves correspond to the sediment beds at pattern 2, pattern 3 and pattern 3, respectively. Panel (b) shows
$\textit{St}_c$
= 1.0,
$d^*$
= 0.6. The red (
$\theta _0 = 8^\circ$
), blue (
$\theta _0 = 16^\circ$
) and grey (
$\theta _0 = 24^\circ$
) curves correspond to the sediment beds at pattern 1, pattern 2 and pattern 3, respectively. Panel (c) shows
$\textit{St}_c$
= 1.0,
$\theta _0 = 16^\circ$
. The red (
$d^*$
= 0.4), blue (
$d^*$
= 0.5) and grey (
$d^*$
= 0.6) curves correspond to the sediment beds at pattern 3, pattern 3 and pattern 2, respectively.
To further explore particle kinematics, representative particle trajectories corresponding to the three sedimentary patterns are analysed. In pattern 1, particle motion is largely stochastic and remains confined near the bed, consistent with a low-disturbance regime. In pattern 2, particles primarily follow upward, arched trajectories, ascending and advecting upstream in a clockwise motion (denoted as path 1, the first column of figure 10). This directed transport leads to localised accumulation zones – corresponding to emergent bed peaks – while the original particle locations evolve into troughs due to sediment depletion.
Schematics of particle trajectories in the sediment-bed pattern 3 during a cycle. The blue and red arrows represent particle paths 2 and 3, respectively.
In pattern 3, particle motion becomes increasingly complex, exhibiting both clockwise and counterclockwise trajectories. A subset of particles continues along a clockwise loop (path 2, indicated by blue arrows in figure 9 and the second column of figure 10), ascending and subsequently advecting downstream. Another subset reverses direction after an initial clockwise ascent, forming a counterclockwise loop (path 3, red arrows in figure 9 and the third column of figure 10), resulting in an ‘8’-shaped trajectory. In path 3, particles initially move downstream (
$ v_x \gt 0, v_y \gt 0$
), then reverse their streamwise velocity (
$ v_x \lt 0$
) while producing a second, larger peak in vertical velocity (
$ v_y$
). This bidirectional motion and associated rotational structure contribute to the formation of secondary suspended particle clouds. Such complex, recirculating trajectories are characteristic of the suspension regime and distinctly differentiate it from the simpler sediment transport modes observed in patterns 1 and 2.
Particle trajectory (a) and velocity in (b) the
$x$
-axis and (c) the
$y$
-axis at various sediment-bed patterns. The particle positions are normalised by the foil chord
$c$
, and the average start position is set to (0,0). The particle velocities are normalised by
$U_\infty$
, the time is normalised by
$T$
and the start time is set to 0. The grey, blue and red curves represent particle paths 1, 2 and 3, corresponding to sediment-bed pattern 2 and two modes of pattern 3, respectively. The particles are selected from
$\theta _0 = 16^\circ$
,
$24^\circ$
and
$24^\circ$
at
$\textit{St}_c$
= 1.0,
$d^*$
= 0.6, respectively. For each type of particle path, 20 particles are selected, all starting at the same time.
5. Wake dynamics
Understanding the wake dynamics is essential for evaluating foil propulsion performance, particularly in particle-laden flows where suspended particles can substantially modify hydrodynamic forces by altering the wake structure. In the presence of particles, the wake becomes more intricate, with frequent particle–vortex and particle–flow interactions. These interactions can disrupt coherent vortex shedding, modify turbulence characteristics and ultimately influence thrust and drag generation. Investigating these effects is therefore critical for elucidating the mechanisms by which sediment particles impact propulsive efficiency and force production. Such insights are especially relevant to marine propulsion systems and other engineering applications involving complex multiphase environments.
5.1. Particle–vortex interactions in the wake
The wake dynamics of the pitching foil is primarily governed by the trailing-edge vortex (TEV). As shown in figures 11(a) and 11(b), a pair of counter-rotating vortices is shed during each oscillation cycle. During the upstroke (
$ t^* = -0.25$
to 0.25), a clockwise (negative) vortex forms, while the downstroke (
$ t^* = 0.25$
to 0.75) generates a counterclockwise (positive) vortex.
Effect of particles on the wake structure for
$\theta _0 = 20^\circ$
,
$d^* = 0.6$
and vorticity normalised as
$\omega ^* = \omega c / U_\infty$
; (a)
$\textit{St}_c$
= 0.4, von Kármán street,
$t^* = 0$
, (b)
$\textit{St}_c = 1.0$
, reverse von Kármán street,
$t^* = 0$
, (c) time-resolved vorticity at the monitoring points with
$\textit{St}_c$
= 0.4, (d) time-resolved vorticity at the monitoring points with
$\textit{St}_c$
= 1.0.
The alternating evolution of these TEVs plays a pivotal role in particle transport and in shaping the sediment bed. During the upstroke, the negative vortex induces a reverse flow near the bed, driving particles along a clockwise path. This mechanism leads to sediment accumulation peaks in pattern 2. As both
$ \textit{St}_c$
and
$ \theta _0$
increase, the size and intensity of these vortices also increase. When the disturbance generated by the negative vortex becomes sufficiently strong, particles detach from the bed and are advected downstream along the vortex path. During the subsequent downstroke, these suspended particles interact with the positive vortex, reversing their streamwise velocity (
$ v_x \gt 0 \rightarrow v_x \lt 0$
) and vertical velocity (
$ v_y \lt 0 \rightarrow v_y \gt 0$
), producing the characteristic ‘8’-shaped trajectories observed in pattern 3. The continual alternation of vortex shedding ultimately results in the formation of branched, ‘Y’-shaped sediment clouds (figure 6
a).
Importantly, the particles are not merely passive tracers; their presence exerts feedback on the flow field and modifies vortex evolution. In particle-laden cases, vortex dissipation downstream is accelerated, reducing both vortex strength and coherence compared with the particle-free case. Further analysis reveals that sediment dampens the sensitivity of vortex street evolution to changes in
$ \textit{St}_c$
. In the presence of particles, both positive and negative vortices align nearly linearly with respect to the bed (left columns of figures 11
a and 11
b), whereas in the absence of particles, the vortex cores are displaced vertically from the bed. At
$ \textit{St}_c = 0.4$
, a typical von Kármán vortex street forms (figure 11
a), with the negative vortex positioned slightly above the positive vortex. At
$ \textit{St}_c = 1.0$
, a reverse von Kármán street emerges (figure 11
b), with the positive vortex elevated above the negative.
To quantitatively assess the impact of particles on vortex evolution, four monitoring points are placed behind the trailing edge at
$ (0.05c, 0)$
,
$ (1c, 0)$
,
$ (1c, 0.3c)$
and
$ (1c, -0.3c)$
. These locations capture time-resolved vorticity across the wake key regions (figures 11
c and 11
d).
At point 1 (
$ 0.05c, 0$
), the peak vorticity in the particle-laden case is reduced by approximately 15 %, although the temporal structure remains consistent: negative vorticity dominates the first half-cycle, followed by positive vorticity in the second.
Further downstream, at one chord length behind the foil (points 2–4), sediment particles induce more pronounced differences. At
$ \textit{St}_c = 0.4$
(figure 11
c), the peak vorticity is reduced and delayed in the particle-laden case. At point 2, the peak vorticity decreases by nearly 50 %, with a delay of approximately
$ \Delta t^* = 0.08$
. At point 3, it decreases by approximately 70 % with a
$ \Delta t^* = 0.06$
delay. The temporal structure also changes: in the particle-laden case, vorticity alternates symmetrically between negative and positive peaks of similar magnitude across cycles. In contrast, the particle-free case shows asymmetric vorticity evolution. Here, successive vortices are connected by thin vorticity filaments (braids) and are influenced by ground-induced shear, producing two distinct negative peaks at
$ t^* = 0.38$
and
$ t^* = 0.47$
at point 2. At point 4 (
$ 1c, -0.3c$
), a sharp positive peak of
$ \omega _+^* = 9.5$
appears at
$ t^* = 0.92$
, approximately twice the magnitude observed in the particle-laden case.
At
$ \textit{St}_c = 1.0$
, the particle-induced effects on wake structure become even more distinct (figure 11
d). In the presence of particles, the monitoring points exhibit a symmetric pattern of alternating positive and negative vorticity. At point 2, the vorticity fluctuates between
$ \pm 7$
, while at points 3 and 4 – located on either side of the wake core – the peak values are reduced to approximately
$ \pm 4$
. In contrast, in the particle-free case, the reverse von Kármán street advects farther from the bed, significantly altering the vorticity distribution. At point 2, negative vorticity dominates, peaking at
$ \omega _-^* = -8.2$
at
$ t^* = 0.88$
, with only weak positive vorticity (
$ {\approx} 1.2$
) appearing between
$ t^* = 0.30$
and 0.60. At point 3, within the core of the negative vortex, vorticity is strongly asymmetric, with a negative peak of
$ \omega _-^* = -19.6$
at
$ t^* = 0.86$
, nearly four times the corresponding positive peak of
$ \omega _+^* = 5.4$
at
$ t^* = 0.12$
. At point 4, located beneath the vortex core, the vorticity remains weakly positive, likely influenced by boundary-layer detachment near the ground.
The observed reduction in vortex strength and coherence in the presence of particles can be attributed to interphase momentum exchange, through which particles continuously extract momentum and energy from the carrier flow (Balachandar & Eaton Reference Balachandar and Eaton2010; Chauchat & Lacaze Reference Chauchat and Lacaze2024). The associated dissipation mechanisms depend on the particle state and spatial distribution.
In the sediment-bed regime, the high local volume fraction and frequent particle contacts substantially restrict particle motion. Under these conditions, fluid flow through the bed induces enhanced pressure gradients across individual particles, leading to increased form drag (Powell Reference Powell2014; Lamb, Brun & Fuller Reference Lamb, Brun and Fuller2017). In addition, partial fluid penetration into the sediment layer weakens wall-induced secondary vorticity generation, promoting a transfer of energy from coherent wake structures toward smaller-scale dissipative motions (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006; Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015; Brandt & Coletti Reference Brandt and Coletti2022).
When particles are entrained into the flow and form a suspended cloud above the bed, momentum exchange is no longer localised near the boundary but distributed throughout the wake region (Balachandar & Eaton Reference Balachandar and Eaton2010). Increased slip velocities between phases, together with particle inertia, introduce pronounced local velocity gradients and phase decorrelation, which tend to disrupt the organisation of vortical structures (Mazzuoli, Balachandar & Chauchat Reference Mazzuoli, Balachandar and Chauchat2024).
Taken together, the combined effects of deposited and suspended particles act as an effective dissipative sink for wake energy, enhancing overall dissipation and leading to a progressive attenuation of turbulent kinetic energy and vortex coherence.
5.2. Streamwise velocity profile
The dynamics of the wake is intrinsically linked to the evolution of vortex structures, which directly shape the streamwise velocity field and thereby influence the foil’s propulsive performance. Figure 12 presents the time-averaged streamwise velocity profiles for three representative combinations of
$ \textit{St}_c$
and
$ \theta _0$
, highlighting the characteristic wake features and the detrimental impact of sediment particles on propulsion efficiency.
Effect of particles on the time-averaged streamwise velocity field, normalised to the relative velocity of the foil by
$U_\infty$
. The solid black curves show profiles along the cross-stream (
$y$
) direction at the streamwise location
$\varDelta$
$x^*$
= 1 (with the equilibrium trailing edge at
$\Delta x^*$
= 0). (a,b) The velocity fields at
$\textit{St}_c$
= 0.4 and
$\textit{St}_c$
= 1.0, showing the transition from velocity deficit to excess as
$\textit{St}_c$
increases. (b,c) The velocity fields at
$\theta _0$
=
$8^\circ$
,
$16^\circ$
, showing a more pronounced upward-inclined jet at higher
$\theta _0$
.
At
$ \textit{St}_c = 0.4$
(figure 12
a), the wake is dominated by a von Kármán vortex street, resulting in a pronounced velocity deficit indicative of drag generation. In the particle-laden case, this deficit intensifies, forming an extended region of reduced velocity downstream. By contrast, the particle-free case exhibits a milder, symmetric velocity deficit centred along the wake axis. Quantitatively, the presence of particles increases the peak velocity deficit by approximately 3.2 %, with the minimum deficit velocity
$\bar {u}_-^*$
decreasing from 0.94 to 0.91 at
$ \theta _0 = 8^\circ$
. This observation is consistent with previous results showing that particles attenuate vortex development and exacerbate wake momentum loss, thereby increasing drag.
At
$ \textit{St}_c = 1.0$
(figure 12
b), the emergence of a reverse von Kármán vortex street gives rise to a velocity excess, forming a jet-like profile characteristic of thrust production. However, the presence of particles significantly reduces the peak jet velocity
$\bar {u}_+^*$
from 1.37 to 1.02 at
$ \theta _0 = 8^\circ$
, corresponding to a 25.5 % decrease. Moreover, the spatial extent of the high-velocity region contracts, with the red jet core reducing to 1.28 chord lengths. At a higher pitch amplitude (
$ \theta _0 = 16^\circ$
), the peak velocity
$\bar {u}_+^*$
diminishes further by 35.6 % (from 1.91 to 1.23), and the jet length shortens to 2.35 chord lengths. These results emphasise the particle-induced suppression of thrust-related velocity structures.
In summary, interactions between particles and the surrounding flow alter the downstream momentum transfer. The resulting amplification of velocity deficit or attenuation of jet strength contributes to a marked reduction in propulsive efficiency. Specifically, particle presence introduces interphase drag acts as a momentum sink, dampening the vortex strength and suppressing asymmetric wake deflection. This ultimately leads to reduced thrust magnitudes and more symmetric propulsive signatures. This effect underscores the importance of accounting for particle–fluid coupling in the performance evaluation of bio-inspired propulsion systems operating in sediment-laden environments.
Effect of particles on the pressure distribution (
$C_{\!p} = {p}/({({1}/{2}) \rho U_\infty ^2})$
), with
$\textit{St}_c$
= 1.0,
$\theta _0 = 10^\circ$
and
$d^*$
= 0.6. (a–d) correspond to
$t^*$
= 0, 0.25, 0.50 and 0.75, respectively.
5.3. Pressure distribution
To elucidate the distinct differences in lift observed at high
$ \textit{St}_c$
and large
$ \theta _0$
– particularly the rapid increase in positive
$ C_{\!L}$
in the particle-free case – we examine the pressure distributions around the foil at
$ \textit{St}_c = 1.0$
for
$ \theta _0 = 10^\circ$
and
$ \theta _0 = 20^\circ$
.
Figure 13 compares the pressure fields for both particle-laden and particle-free conditions at
$ \theta _0 = 10^\circ$
. In the particle-laden case, the presence of particles increases pressure on the upper surface during the upstroke (figure 13
a), resulting in a slight increase in negative lift. The pressure on the lower surface is also enhanced during the downstroke, leading to an increase in positive lift. For instance, at
$ t^* = 0.5$
(figure 13
c), an intensified high-pressure region near the lower trailing edge, combined with a relatively weaker negative pressure near the upper leading edge, produces a larger vertical pressure differential, thereby enhancing positive lift. Consequently, the time-averaged lift coefficients for both particle-laden and particle-free cases remain similar (both yielding
$ C_{\!L} = -0.12$
). Moreover, pressure disturbances in the particle case are largely confined to within one chord length downstream, whereas the particle-free case exhibits a more extended and coherent vortex-induced low-pressure wake.
At a higher pitch amplitude
$ \theta _0 = 20^\circ$
, the pressure fields (figure 14) show markedly different behaviour. In the particle-laden case, the upstroke produces an intensified positive pressure on the upper surface (figure 14
a), leading to greater negative lift. However, during the downstroke, the positive pressure on the lower surface weakens (figure 14
c), thereby reducing the positive lift contribution. As a result, the time-averaged lift remains negative and exhibits a larger magnitude than at
$ \theta _0 = 10^\circ$
.
In contrast, the leading-edge vortex (LEV) affects surface pressure distribution in the particle-free case. During the upstroke, the separation of LEV on the upper surface reduces the surface pressure, evidenced by a reversed pressure gradient at
$ t^* = 0$
(figure 14
a) and a transition from negative to positive pressure near the leading edge at
$ t^* = 0.75$
(figure 13
d). This reduces the negative lift during the upstroke. Conversely, during the downstroke, an intensified positive pressure region near the lower trailing edge increases the positive lift (figure 14
c). As a result, the positive lift peak exceeds the negative lift peak, leading to a net increase in the time-averaged lift in the particle-free case at high
$ \theta _0$
.
Effect of particles on the pressure distribution (
$C_{\!p} = {p}/({({1}/{2}) \rho U_\infty ^2})$
), with
$\textit{St}_c$
= 1.0,
$\theta _0 = 20^\circ$
and
$d^*$
= 0.6. a–d correspond to
$t^*$
= 0, 0.25, 0.50 and 0.75, respectively.
6. Conclusions
This study reveals the critical influence of sediment bed on the propulsion performance of pitching hydrofoils. The presence of particles degrades thrust generation and propulsive efficiency by weakening vortex formation, suppressing momentum flux and modifying surface pressure distributions. In addition, the sediment bed undergoes dynamic deformation driven by flow-induced forces, with suspended particle motion playing a central role in sediment transport and morphological evolution. Notably, the foil–sediment interaction occurs through an indirect, fluid-mediated coupling rather than direct contact, whereby the hydrodynamic field generated by foil motion governs both particle response and bed evolution.
The results highlight a complex feedback mechanism linking foil kinematics, wake dynamics and sediment behaviour. This tripartite coupling between the moving structure, the surrounding flow and the particle phase governs both propulsion performance and sediment pattern formation. Each component influences the others, resulting in a dynamically evolving system that cannot be adequately captured using uncoupled or simplified models. These findings underscore the importance of explicitly accounting for fluid–structure–particle interactions when analysing bio-inspired propulsion in particle-laden environments.
Overall, this work provides new insights into the hydrodynamic and morphodynamic consequences of pitching hydrofoils operating near deformable beds. The results advance our understanding of propulsion in particle-laden flows and offer guidance for the design of biomimetic underwater vehicles intended for sediment-rich environments. In the present study, particles are modelled to represent synthetic or organic matter commonly encountered in near-shore and aquaculture settings, turbulence is treated within a RANS framework to focus on the dominant coherent wake structures governing the underlying flow physics, and idealised DEM parameters are employed to capture essential particle behaviour. These modelling choices are intended to provide a physically consistent and computationally tractable baseline. Future studies may extend this framework by systematically varying particle properties (e.g. size and density), concentration levels and foil kinematics (e.g. heaving or combined pitch–heave motions), thereby further elucidating interaction mechanisms and supporting the development of adaptive, high-efficiency propulsion systems for real-world marine applications.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Numerical validation and verification
To ensure the reliability of the numerical simulations, a systematic sensitivity analysis regarding mesh resolution and time step is performed.
A.1. Open water case
The accuracy of the numerical model is first evaluated using a two-dimensional NACA 0012 foil undergoing harmonic heaving and pitching in open water, replicating the experimental set-up of Read, Hover & Triantafyllou (Reference Read, Hover and Triantafyllou2003). The motion parameters are set to a heave amplitude of
$h_0/c = 0.75$
, a Strouhal number of
$St = 0.4$
, a maximum angle of attack of
$\alpha _{\max } = 20^\circ$
and a phase angle of
$\phi = 90^\circ$
. The numerical results for the thrust coefficient (
$C_{\!T}$
) and efficiency (
$\eta$
) are compared against the experimental values (
$C_{\!T} = 0.79$
and
$\eta = 55.6\%$
), with the results presented in terms of relative error percentages.
Figure 15(a) illustrates the five mesh configurations under evaluation, where the base size varies from
$ 1.1 \times 10^{-2}\, \text{}$
to
$ 1.9 \times 10^{-2}\, \text{m}$
in steps of
$ 2.0 \times 10^{-3}\, \text{m}$
. As the mesh size decreases, the relative errors of
$C_{\!T}$
(red line) and
$\eta$
(blue line) show a general convergence trend. Specifically, when the mesh base size is refined to
$ 1.5 \times 10^{-2}\, \text{m}$
or below, the errors remain within
$5\,\%$
, indicating that the numerical results are no longer sensitive to further mesh refinement. Furthermore, based on the selected mesh base size of
$ 1.5 \times 10^{-2}\, \text{m}$
, the independence of the time step
$\Delta t$
is validated using five different levels. The relative errors of
$C_{\!T}$
and
$\eta$
are consistently within the
$10\,\%$
range (figure 15
b).
Comparison of thrust coefficient
$C_{\!T}$
and efficiency
$\eta$
under different sets of (a) mesh and (b) time step.
A.2. Single-wall case
To further ensure the reliability of the simulation under ground effect conditions, the variations of the lift coefficient (
$C_l$
) and moment coefficient (
$C_m$
) over one cycle for a single-wall case are compared with experimental data from Su et al. (Reference Su, Miller, Mandre and Breuer2019) and numerical data from Yang et al. (Reference Yang, He, Mao, Mo and Ghassemi2023). Three mesh configurations with different base sizes are evaluated: mesh 1 (base size =
$ 1.8 \times 10^{-2}\, \text{m}$
), mesh 2 (base size =
$ 1.5 \times 10^{-2}\, \text{m}$
) and mesh 3 (base size =
$ 1.2 \times 10^{-2}\, \text{m}$
). The time step is set to
$ \Delta t = 5.0 \times 10^{-4}\, \text{s}$
. As shown in figure 16(a), the lift coefficient predicted by the three mesh configurations exhibits good consistency. However, the moment coefficient near
$t^* = 0.5$
predicted by the coarser mesh scheme, mesh 1, is higher than both the experimental result and the predictions from the finer meshes (figure 16
b). Additionally, the peak force in our results is slightly larger than that observed in the experiment but remains consistent with the numerical simulation result, likely due to the absence of energy loss associated with three-dimensional effect.
Comparison of (a) lift coefficient
$C_l$
and (b) moment coefficient
$C_m$
under different sets of mesh.
Comparison of (a) lift coefficient
$C_l$
and (b) moment coefficient
$C_m$
under different time steps.
Qualitative comparison between simulation and experiment visualisation.
The independence of the time step is also verified at
$ \Delta t_1 = 1.0 \times 10^{-3}\, \text{s}$
,
$ \Delta t_2 = 5.0 \times 10^{-4}\, \text{s}$
and
$ \Delta t_3 = 2.5 \times 10^{-4}\, \text{s}$
, using a mesh base size of
$ 1.5 \times 10^{-2}\, \text{m}$
(mesh 2). Similar trends to the mesh sensitivity analysis are observed, with our results showing good agreement with the literature data (figure 17). Considering both computational efficiency and accuracy, a mesh base size of
$ 1.5 \times 10^{-2}\, \text{m}$
and a time step of
$ \Delta t = 5.0 \times 10^{-4}\, \text{s}$
are adopted for the simulations.
A.3. Sediment-bed case
To further ensure the reliability of the simulation under sediment-bed conditions, a feasibility verification of the particle and sediment-bed modelling is performed by comparing the current numerical approach with the those presented by Aju et al. (Reference Aju, Gong, Pham, Kaushik and Jin2022). As shown in figure 18, under conditions of
$\textit{St}_D$
= 0.26 and
$A_D$
= 0.67, the suspended particles simulated using our approach closely match the experimental results. The hydrofoil advances and pitches over a sediment bed, causing the particles to lift off the bed and form a suspended flow. Notably, a ‘Y’-shaped sediment cloud behind the foil is accurately captured in our simulation.
In conclusion, the mesh and time step independence, as well as the simulation approach for the pitching foil over the sediment bed, are effectively verified, demonstrating the reliability and accuracy of the proposed numerical method.
