1. Introduction
Fluid–flexible-structure interaction problems are ubiquitous in nature (Dowell & Hall Reference Dowell and Hall2001), where deformable bodies interact with the surrounding fluid and thereby affect the dynamics of both. Depending on whether energy is actively supplied to the structure, these problems can be broadly divided into fluid–active-flexible-structure interactions and fluid–passive-flexible-structure interactions. In the active case, flexible bodies consume energy to generate prescribed motions for propulsion or force production, as observed in flying birds and insects (Wu Reference Wu2011; Liu et al. Reference Liu, Wang and Liu2024) and swimming fish (Triantafyllou et al. Reference Triantafyllou, Triantafyllou and Yue2000; Smits Reference Smits2019). In contrast, passive systems involve no active energy input, and the structure undergoes translation and deformation solely under fluid loading. These features naturally arise in a wide range of configurations, including plant reconfiguration (de Langre Reference de Langre2008), flapping flags (Shelley & Zhang Reference Shelley and Zhang2011), and fibre deformation and clustering in viscous flows (du Roure et al. Reference du Roure, Lindner, Nazockdast and Shelley2019; Marchioli et al. Reference Marchioli, Rosti and Verhille2026). Investigating the dynamics of such fluid–structure systems not only helps to elucidate the fundamental mechanisms underlying these natural phenomena, but also provides important inspiration for engineering applications (Nepf Reference Nepf2012; Lauder Reference Lauder2015; Liu et al. Reference Liu, Wang, Liu and He2023; Johansson Reference Johansson2024; Rong & Liu Reference Rong and Liu2025). As a result, passive fluid–flexible-structure interactions have been extensively studied in recent years (Wang et al. Reference Wang, He and Zhang2016; Griffith & Patankar Reference Griffith and Patankar2020).
The role of flexibility in enhancing the propulsive performance of actively flapping systems has been extensively investigated, together with the underlying physical mechanisms (Kang et al. Reference Kang, Aono, Cesnik and Shyy2011; Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012; Peng et al. Reference Peng, Sun, Yang, Xiong, Wang and Wang2022; Liu & Huang Reference Liu and Huang2024b
). It has been shown that a flexible plate with moderate bending stiffness can achieve optimal propulsion speed (Peng et al. Reference Peng, Huang and Lu2018), and this advantage persists even under intermittent flapping conditions (Liu et al. Reference Liu, Huang and Lu2020). In addition to improving swimming speed, flexibility can significantly enhance thrust production, with flexible foils generating up to three times the thrust of their rigid counterparts (Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012). From an energetic perspective, flexibility also plays a crucial role. Zhu et al. (Reference Zhu, He and Zhang2014) demonstrated that a flexible self-propelled plate achieves substantially higher propulsion efficiency than a rigid one when the mass ratio
$M$
is small. Similarly, Xu & Yu (Reference Xu and Yu2023) reported that bodies with intermediate stiffness attain the highest efficiency, whereas higher stiffness requires greater energy input to maintain the same flapping amplitude. In addition to constant structural properties, recent studies have further explored more complex forms of flexibility, including time-periodic stiffness variations (Shi et al. Reference Shi, Xiao and Zhu2020; Yudin et al. Reference Yudin, Floryan and Van Buren2023; Lu et al. Reference Lu, Zhang, Xu and Huang2025) and spatially distributed flexibility (Floryan & Rowley Reference Floryan and Rowley2020; Wang et al. Reference Wang, Huang and Lu2020; Xu et al. Reference Xu, Zhu, Hao and Yin2026), highlighting additional pathways for performance enhancement.
The passive motion of flexible bodies in various flow conditions has also attracted widespread attention (Eloy et al. Reference Eloy, Lagrange, Souilliez and Schouveiler2008; Jia & Yin Reference Jia and Yin2008; Ristroph & Zhang Reference Ristroph and Zhang2008). Note that in actively flapping systems, structural deformation is largely dictated by prescribed input, typically in the form of heaving, pitching or undulatory motions (Ramananarivo et al. Reference Ramananarivo, Godoy-Diana and Thiria2011; Floryan & Rowley Reference Floryan and Rowley2018; Xu et al. Reference Xu, Zhou, Yu and Yin2025). However, in passive systems, the motion and deformation of flexible bodies emerge entirely from fluid forcing, leading to much richer and more complex dynamics. The fluid can induce multiple response modes of the structure, while the structural deformation, in turn, modifies the flow and may further amplify flow instabilities. As a result, the effects of flexibility become more intricate, and a deeper understanding is still needed. A typical example is a conventional flag fixed at its leading edge and immersed in an incoming flow (Shelley et al. Reference Shelley, Vandenberghe and Zhang2005; Huang et al. Reference Huang, Shin and Sung2007; Alben & Shelley Reference Alben and Shelley2008). This system exhibits several distinct dynamical states, including stretched-straight, coherent-flapping and chaotic-flapping regimes (Zhang et al. Reference Zhang, Childress, Libchaber and Shelley2000; Connell & Yue Reference Connell and Yue2007; Shelley & Zhang Reference Shelley and Zhang2011). Another well-studied configuration is an inverted plate or filament clamped at the trailing edge (Kim et al. Reference Kim, Cossé, Huertas Cerdeira and Gharib2013; Deng et al. Reference Deng, Mao and Brandt2021). As the bending stiffness increases, the system transitions through a sequence of modes, including deflected–flapping, deflected, flapping, and straight states (Huang et al. Reference Huang, Wei and Lu2018). Despite the apparent complexity of these dynamical behaviours, underlying scaling laws can still be identified, revealing the governing physical mechanisms. For example, a compliant fibre tethered at its midpoint has been shown to exhibit a drag scaling proportional to
$U_\infty ^{4/3}$
at high Reynolds numbers (Alben et al. Reference Alben, Shelley and Zhang2002; Zhu Reference Zhu2008). Here
$U_\infty$
is the oncoming flow speed. This scaling differs markedly from the classical
$U_\infty ^2$
dependence for rigid bluff bodies (Batchelor Reference Batchelor1967), indicating that flexible bodies can achieve significant drag reduction through self-similar bending (Alben et al. Reference Alben, Shelley and Zhang2002).
In the aforementioned studies of passive systems, the presence of fixed points often constrains the overall motion of the structure. When such constraints are relaxed, freely moving flexible fibres can exhibit a much broader range of dynamical behaviours in flow (du Roure et al. Reference du Roure, Lindner, Nazockdast and Shelley2019). Forgacs & Mason (Reference Forgacs and Mason1959) experimentally investigated real fibres in Stokes shear flow, and identified five distinct modes, including rigid motion, springy motion, two types of snake turns, and coiled motion. Subsequent studies (Smith et al. Reference Smith, Babcock and Chu1999; Becker & Shelley Reference Becker and Shelley2001; Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018; Slowicka et al. Reference Slowicka, Stone and Ekiel-Jeżewska2020; Xu et al. Reference Xu, Liu, Liu, Xiong and Huang2022) have further confirmed and extended these classifications. Another representative problem is the sedimentation of flexible fibres in viscous fluids under external forcing such as gravity. In this case, fibres can develop a variety of complex shapes and settling trajectories (Cosentino Lagomarsino et al. Reference Cosentino Lagomarsino, Pagonabarraga and Lowe2005; Marchetti et al. Reference Marchetti, Raspa, Lindner, du Roure, Bergougnoux, Guazzelli and Duprat2018; Lau & Huang Reference Lau and Huang2021; Xu et al. Reference Xu, Li, Li and Liao2021).
Except these tethered and freely moving configurations, most recently, Liu et al. (Reference Liu, Gao and Lu2026) proposed a novel unconstrained configuration in which a flexible filament is free to move laterally. Five distinct motion modes were identified, namely static, quasi-static, translational, vibrational and chaotic modes (Liu et al. Reference Liu, Gao and Lu2026). The filament can undergo spontaneous translation or vibration without external mechanical actuation in most cases, with the required mechanical energy being extracted directly from the surrounding flow. This capability of converting fluid kinetic energy into the translational or oscillatory motion of a flexible structure may be relevant to the development of passive flow-energy-harvesting systems. For example, the sustained vibration of the filament could potentially be coupled with an appropriate electromechanical transducer to extract energy from a continuous fluid stream. In addition, the spontaneous and sustained translation observed in the translational mode suggests another possible application in the passive transport or manipulation of lightweight flexible materials in fluid environments, where their direction of motion could potentially be controlled through the flow conditions and structural properties.
Note that in the study of Liu et al. (Reference Liu, Gao and Lu2026), the incoming flow was normal to the filament. In practical fluid environments, however, the relative orientation between the incoming flow and a flexible structure is generally not fixed and may vary with the flow direction or the configuration of the structure (Bintein et al. Reference Bintein, Bense, Clanet and Quéré2019). Such variations can substantially influence the flow symmetry and lift generation of foils (Cleaver et al. Reference Cleaver, Wang and Gursul2012), the flapping modes of inverted cantilever plates (Huertas-Cerdeira et al. Reference Huertas-Cerdeira, Goza, Sader, Colonius and Gharib2021), and the energy-harvesting performance of travelling-wave plates (Liang et al. Reference Liang, Qian, Cui and Huang2025). Therefore, understanding how an unconstrained flexible filament behaves under an inclined flow represents a natural and practically relevant extension of the normal-flow configuration.
Inspired also by the upwind motion of droplets in a crosswind (Bintein et al. Reference Bintein, Bense, Clanet and Quéré2019), the present study investigates the dynamics of unconstrained filaments under more general inflow conditions, in which the incoming flow is inclined at an angle
$\delta$
relative to the horizontal direction, as illustrated in figure 1. The inclined flow breaks the geometric symmetry of the system and introduces a preferred lateral direction. Consequently, the filament may exhibit directional rather than randomly selected lateral motion, providing a simple framework for understanding and potentially controlling the passive transport of flexible structures in asymmetric flow environments. We systematically examine the effects of initial perturbations, inflow inclination angle, Reynolds number and bending stiffness on the dynamical characteristics of thin filaments. Furthermore, force decomposition and stability analyses are performed to elucidate the mechanisms underlying the bifurcation phenomena observed in the system.
The remainder of this paper is organised as follows. The physical problem and mathematical formulation are given in
$\S$
2. The numerical method and validation are presented in
$\S$
3. Detailed results are discussed in
$\S$
4, and concluding remarks are provided in
$\S$
5.
2. Physical problem and mathematical formulation
Figure 1 presents a schematic of the flexible filament considered in this study. The filament, of length
$L$
, is immersed in a uniform incoming flow with velocity
$\boldsymbol U_\infty$
, inclined at an angle
$\delta$
to the horizontal. Accordingly, the
$x$
- and
$y$
-components of
$\boldsymbol U_\infty$
are
$U_{\infty ,x} = U_\infty \cos \delta$
and
$U_{\infty ,y} = U_\infty \sin \delta$
, respectively, where
$U_\infty$
denotes the magnitude of
$\boldsymbol U_\infty$
. The two ends of the filament are constrained in the
$x$
-direction but remain free in the
$y$
-direction, whereas the rest of the filament moves freely under fluid–structure interactions. When
$\delta \neq 0$
, the system becomes asymmetric, and the filament is expected to undergo lateral translation. Intuitively, the filament may drift in the direction of the incoming flow, similar to a fallen leaf carried by the wind.
Schematic of a flexible filament in an inclined flow. A filament of length
$L$
is immersed in a uniform incoming flow with velocity
$\boldsymbol U_\infty$
, inclined at an angle
$\delta$
to the horizontal. Here,
$U_{\infty ,x} = U_\infty \cos \delta$
and
$U_{\infty ,y} = U_\infty \sin \delta$
denote the
$x$
- and
$y$
-components of
$\boldsymbol U_\infty$
, respectively, where
$U_\infty$
is the magnitude of
$\boldsymbol U_\infty$
. The two ends of the filament are fixed in the
$x$
-direction but are free to move in the
$y$
-direction. Here,
$s$
denotes the curvilinear coordinate along the filament, and
$V$
denotes the lateral velocity of the filament’s centre of mass, with
$V\gt 0$
and
$V\lt 0$
corresponding to downwind and upwind motion, respectively.
The incompressible Navier–Stokes equations are used to govern the fluid flow:
\begin{equation} \frac {\partial \boldsymbol{v}}{\partial t}+\boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{v}=-\frac {1}{\rho } \boldsymbol{\nabla } \!p+\frac {\mu }{\rho } {\nabla} ^{2} \boldsymbol{v}+\boldsymbol{f}_{\!b},\end{equation}
\begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{v}=0, \end{equation}
where
$\boldsymbol{v}$
is the velocity,
$p$
is the pressure,
$\rho$
is the density of the fluid,
$\mu$
is the dynamic viscosity, and
$\boldsymbol{f}_{\!b}$
denotes the body force acting on the surrounding fluid due to the immersed boundary.
The structural equation is adopted to describe the deformation and motion of the filament (Zhu & Peskin Reference Zhu and Peskin2002; Connell & Yue Reference Connell and Yue2007):
\begin{equation} \rho _{l} \frac {\partial ^{2} \boldsymbol{X}}{\partial t^{2}}-\frac {\partial }{\partial s}\left [\textit{Eh} \left (1-\left |\frac {\partial \boldsymbol{X}} {\partial s}\right |^{-1}\right ) \frac {\partial \boldsymbol{X}}{\partial s}\right ]+ \textit{EI} \frac {\partial ^{4} \boldsymbol{X}}{\partial s^{4}}=\boldsymbol{f}_{\!s}, \end{equation}
where
$s$
is the Lagrangian coordinate along the filament,
$\rho _{l}$
is the structural linear mass density,
$\boldsymbol{X}(s, t)=(X(s, t), Y(s, t))$
is the position vector of the filament, and
$\boldsymbol{f}_{\!s}$
is the Lagrangian force exerted on the filament by the surrounding fluid. Here,
$Eh$
and
$EI$
denote the structural stretching rigidity and bending rigidity, respectively.
At the two endpoints of the filament, the boundary conditions are
\begin{equation} -\textit{Eh}\left (1-\left |\frac {\partial \boldsymbol{X}}{\partial s}\right |^{-1}\right ) \frac {\partial Y}{\partial s}+\textit{EI}\,\frac {\partial ^{3} Y}{\partial s^{3}}=0, \quad X(t)=0, \quad \frac {\partial ^{2} \boldsymbol{X}}{\partial s^{2}}=0. \end{equation}
Note that the first boundary condition in (2.4) represents the vertically unconstrained condition. In addition,
$\boldsymbol{X}(s, 0)=(0,s),\ (\partial \boldsymbol{X} / \partial t)(s, 0)=(0,0)$
are the initial conditions for the filament.
In the present study, the fluid density
$\rho$
, dynamic viscosity
$\mu$
, and filament length
$L$
are held constant. To non-dimensionalise the governing equations,
$\rho$
,
$L$
and
$U_{\textit {ref}}=\lambda \mu / \rho L$
are selected as the characteristic density, length and velocity scales, respectively, where
$\lambda =200$
is a prescribed constant. Note that the oncoming flow velocity
$U_\infty$
is not adopted as the characteristic velocity because it varies among the cases considered. The particular choice of
$\lambda$
only sets the reference velocity scale, and other reasonable values could be used without altering the underlying physics or the conclusions of the present study (Liu & Huang Reference Liu and Huang2024a
). The corresponding characteristic time scale is therefore defined as
$T_{\textit {ref}}=L/U_{\textit {ref}}$
. Based on dimensional analysis, the dimensionless governing parameters are introduced as follows: the inflow inclination angle
$\delta$
, the Reynolds number based on the oncoming flow speed
$\textit{Re}=\rho U_\infty L / \mu$
, the filament-to-fluid mass ratio
$M=\rho _{l}/ \rho L$
, the stretching stiffness
$S =\textit{Eh} / \rho U_{\textit {ref}}^{2} L$
, and the bending stiffness
$K =\textit{EI} / \rho U_{\textit {ref}}^{2} L^{3}$
.
3. Numerical method and validation
3.1. Numerical method
The lattice Boltzmann method, which has been widely used to simulate complex flows (Chen & Doolen Reference Chen and Doolen1998), is employed to solve the Navier–Stokes equations numerically. The lattice Boltzmann equation with the Bhatnagar–Gross–Krook (BGK) collision model is
\begin{equation} f_i\left (\boldsymbol{x}+\boldsymbol{e}_i\, \Delta t, t+\Delta t\right )-f_i(\boldsymbol{x}, t)=-\frac {1}{\tau }\left [f_i(\boldsymbol{x}, t)-f_i^{e q}(\boldsymbol{x}, t)\right ]+\Delta t\, F_i, \end{equation}
where
$f_i(\boldsymbol{x}, t)$
is the distribution function associated with discrete particle velocity
$\boldsymbol e_i$
,
$\Delta t$
is the time step, and
$\tau$
is the non-dimensional relaxation time related to fluid viscosity, i.e.
$\nu =c_s^2(\tau -0.5)\, \Delta t$
. The equilibrium distribution function
$f_i^{eq}$
and the forcing term
$F_i$
are defined as
\begin{equation} f_i^{e q}=\omega _i \rho \left [1+\frac {\boldsymbol{e}_i \boldsymbol{\cdot }\boldsymbol{v}}{c_s^2}+\frac {\left (\boldsymbol{e}_i \boldsymbol{\cdot }\boldsymbol{v}\right )^2}{2 c_s^4}-\frac {\boldsymbol{v}^2}{2 c_s^2}\right ] \end{equation}
and
\begin{equation} F_i=\left (1-\frac {1}{2 \tau }\right ) \omega _i\left [\frac {\boldsymbol{e}_i-\boldsymbol{v}}{c_s^2}+\frac {\boldsymbol{e}_i \boldsymbol{\cdot }\boldsymbol{v}}{c_s^4} \boldsymbol{e}_i\right ] \boldsymbol{\cdot }\boldsymbol{f}_{\!b}. \end{equation}
Here,
$c_s=c/\sqrt {3}$
denotes the lattice speed of sound, where
$c=\Delta x/\Delta t$
is the lattice speed, and
$\Delta x$
is the grid spacing, while
$\omega _i$
denotes the corresponding weighting factor. In the present simulations, the D2Q9 lattice model is adopted, for which
\begin{equation} \boldsymbol{e}_i=c \boldsymbol{\cdot }\left [\begin{array}{ccccccccc} 0 & 1 & 0 & -1 & 0 & 1 & -1 & -1 & 1 \\ 0 & 0 & 1 & 0 & -1 & 1 & 1 & -1 & -1 \end{array}\right ] \end{equation}
and
\begin{equation} \omega _i= \begin{cases}4/9, & i=0, \\[3pt] 1/9, & i=1{-}4, \\[3pt] 1/36, & i=5{-}8 .\end{cases} \end{equation}
The macroscopic density, pressure and velocity are recovered from the distribution functions as
\begin{equation} \rho =\sum _i f_i, \quad p=c_s^2 \rho , \quad \rho \boldsymbol{v}=\sum _i \boldsymbol{e}_i f_i+\frac {1}{2} \boldsymbol{f}_{\!b}\, \Delta t. \end{equation}
The finite-element method (Doyle Reference Doyle2001) is used to simulate the motion and deformation of the flexible filament. Large filament displacements are resolved within a co-rotational framework (Doyle Reference Doyle2001), which is well suited for problems involving geometric nonlinearity. Fluid–structure coupling is realised through the immersed boundary method (IBM) (Zhu & Peskin Reference Zhu and Peskin2002; Mittal & Iaccarino Reference Mittal and Iaccarino2005), which enables bidirectional interaction between the fluid and structural domains.
In our IBM framework, the Lagrangian force
$\boldsymbol{f}_{\!s}$
is computed using a penalty scheme to enforce the no-slip boundary condition. Specifically,
$\boldsymbol{f}_{\!s}$
is evaluated as
\begin{equation} \boldsymbol{f}_{\!s}=\alpha _1\int _{0}^{t} \left [\boldsymbol{V}_{\!f}(s,t^{\prime })-\boldsymbol{V}_{\!s}(s,t^{\prime }) \right ]\mathrm{d}t^{\prime } + \alpha _2\left [\boldsymbol{V}_{\!f}(s,t)-\boldsymbol{V}_{\!s}(s,t)\right ]\!, \end{equation}
where
$\alpha _1$
and
$\alpha _2$
are penalty coefficients chosen to be sufficiently large to effectively minimise the velocity slip between the fluid and the solid structure (Hua et al. Reference Hua, Zhu and Lu2014; Zhang et al. Reference Zhang, Huang and Lu2020). Here,
$\boldsymbol{V}_{\!s}=\partial \boldsymbol{X}/\partial t$
represents the physical velocity of the filament at the Lagrangian points, and
$\boldsymbol{V}_{\!f}$
denotes the local fluid velocity obtained through interpolation from the Eulerian grid to the Lagrangian points, i.e.
\begin{equation} \boldsymbol{V}_{\!f}(s,t)=\int \boldsymbol{v}(\boldsymbol{x},t)\,\delta (\boldsymbol{x}-\boldsymbol{X}(s,t))\, \mathrm{d}\boldsymbol{x}. \end{equation}
The body force
$\boldsymbol f_b$
on the Eulerian points can be obtained from the Lagrangian force
$\boldsymbol{f}_{\!s}$
by applying the Dirac
$\delta$
function, i.e.
\begin{equation} \boldsymbol{f}_{\!b} (\boldsymbol{x}, t)=-\int \boldsymbol{f}_{\!s}(s, t)\, \delta (\boldsymbol{x}-\boldsymbol{X}(s, t))\, \mathrm{d} s. \end{equation}
In the present study, a four-point regularised
$\delta$
function is adopted (Lai & Peskin Reference Lai and Peskin2000; Peskin Reference Peskin2002), i.e.
\begin{equation} \begin{aligned} \delta &=\frac {1}{\Delta x\, \Delta y} \phi \left (\frac {x}{\Delta x}\right ) \phi \left (\frac {y}{\Delta y}\right )\!, \\[3pt] \phi (r) &= \begin{cases}\left (3-2|r|+\sqrt {1+4|r|-4 r^2}\right ) / 8, & |r|\lt 1, \\[7pt] \left (5-2|r|-\sqrt {-7+12|r|-4 r^2}\right ) / 8, & 1 \leqslant |r|\lt 2, \\[7pt] 0, & |r| \geqslant 2,\end{cases} \end{aligned} \end{equation}
where
$|r|$
is the distance between the Lagrangian point and the nearby Eulerian grid points.
3.2. Validation
In the simulations, the computational domain for the fluid flow is defined as
$[-15L,25L] \times [-15L,15L]$
in the
$x$
- and
$y$
-directions, respectively, which is sufficiently large to eliminate boundary-blockage effects. A uniform velocity condition
$\boldsymbol{v}=(U_{\infty ,x},U_{\infty ,y})$
is imposed at the inlet and side boundaries, while a convective boundary condition is specified at the outlet. Initially, if not otherwise specified, the fluid velocity is prescribed as
$\boldsymbol{v}=(U_{\infty ,x},U_{\infty ,y})$
throughout the domain. In § 4.1, the effects of the initial perturbation are considered. Without loss of generality, a sinusoidal perturbation is introduced into the lateral velocity component of the initial flow field (Liu et al. Reference Liu, Gao and Lu2026). Specifically, under the circumstances, the flow field at
$t=0$
is prescribed as
\begin{equation} \boldsymbol v(\boldsymbol x,0) = \left (U_{\infty ,x}, \ U_{\infty ,y} + \beta U_{\textit {ref}} \sin (2 \pi x/L) \right ), \end{equation}
where
$\beta$
is a dimensionless parameter representing the perturbation amplitude, and
$\beta = 0$
corresponds to the unperturbed initial condition. Note that the form of the initial perturbation is not unique, and other perturbation forms are equally admissible. Here, a sinusoidal perturbation is chosen for its simplicity and smoothness. Its role is to perturb the initially uniform flow field, and facilitate the emergence of different flow responses, thereby enabling the bifurcation behaviour and subsequent evolution of the system to be explored. The specific form of the perturbation is not expected to affect the intrinsic dynamics of the resulting stable state.
In addition, a finite moving computational domain is adopted in the
$y$
-direction to allow the filament to travel over a sufficiently long distance; specifically, whenever the filament moves by one lattice spacing in the
$y$
-direction, the computational domain is shifted accordingly (Liu et al. Reference Liu, Liu and Huang2022).
(a) Grid-independence and (b) time-step-independence studies for the case with
$\delta =30^\circ$
,
$\textit{Re}=100$
,
$K=0.01$
,
$M=0.5$
and
$S=10^4$
. The normalised streamwise force
$F_x/F_{\textit {ref}}$
is plotted as a function of time, where
$F_{\textit {ref}}=(1/2)\rho U_{\textit {ref}}^2L$
.
The numerical strategy employed in this work has been rigorously validated and successfully applied to a range of fluid–structure interaction problems, including the scaling laws governing the self-propelled performance of flexible plates (Liu et al. Reference Liu, Liu and Huang2022), the wake transition of an unconstrained self-propelled plate (Liu & Huang Reference Liu and Huang2024b ), and the flow-induced reconfiguration of flexible ribbons (Liu & Huang Reference Liu and Huang2024a ). Details of the numerical validation can be found in the corresponding references.
Figure 2 presents the results of the grid- and time-step-independence tests for the case with
$\delta = 30^\circ,\,\textit{Re} = 100$
,
$K=0.01$
,
$M=0.5$
and
$S=10^4$
. It is seen that
$\Delta x / L=0.01$
and
$\Delta t / T_{\textit {ref}}=0.00025$
are sufficient to ensure numerical accuracy. Therefore, these values are adopted in all subsequent simulations.
4. Results and discussion
The settings of all governing parameters used in the present simulations are summarised in table 1. The mass ratio
$M$
and stretching stiffness
$S$
are kept constant throughout the study. A sufficiently large value of
$S$
is selected to ensure that the filament remains nearly inextensible. The chosen value of
$M$
is consistent with those adopted in previous studies of flexible bodies in uniform flow (Connell & Yue Reference Connell and Yue2007; Zhu Reference Zhu2008; Huang & Sung Reference Huang and Sung2010), which typically lie within the range
$M \in [0.1, 1]$
. The primary parameters varied in this study are the inflow Reynolds number
$\textit{Re}$
, the bending stiffness
$K$
, the inflow inclination angle
$\delta$
, and the initial perturbation intensity
$\beta$
. Flexible filaments in inclined flow generally exhibit regular dynamics, but at sufficiently high Reynolds numbers, they may enter a chaotic regime characterised by highly irregular fluid forces and filament motions (Liu et al. Reference Liu, Gao and Lu2026). To focus on systematic bifurcation behaviour and avoid the additional complexity associated with chaos, we restrict
$\textit{Re}$
to the range
$[10, 200]$
.
The parameters in the simulations.
Table 1. Long description
The table lists parameters used in simulations, including mass ratio, stretching stiffness, inflow Reynolds number, bending stiffness, inflow inclination angle, and initial perturbation intensity. The mass ratio and stretching stiffness are kept constant, while the other parameters are varied. The inflow Reynolds number ranges from 10 to 200, bending stiffness values are 10^−4, 10^−2, and 1, inflow inclination angles vary from 0 to 50 degrees, and initial perturbation intensity ranges from −2 to 2.
In the following discussion, the time average of a variable
$q$
is defined as
\begin{equation} \langle q\rangle := \frac { 1 }{ t_1 - t_0 } \int _{t_0}^{t_1} q(t)\, \mathrm{d} t. \end{equation}
In our computations, the averaging interval
$t_1-t_0$
is chosen according to the nature of the quantity under consideration. For periodic quantities, it is taken as the fundamental period
$T=1/f_1$
, where
$f_1$
denotes the dominant frequency. For statistically steady quantities, the averaging interval may be chosen arbitrarily. For convenience, it is taken to be the reference time scale
$T_{\textit {ref}}$
in this study.
4.1. Emergent dynamics and bifurcations
The dynamic characteristics of filaments under an initially symmetric configuration, i.e.
$\delta = 0^\circ$
, were investigated thoroughly in a previous study (Liu et al. Reference Liu, Gao and Lu2026). When
$\delta \gt 0^\circ$
, the inherent symmetry of the system is broken, and the resulting dynamics therefore requires careful re-examination. In this subsection, we further examine the wake structures and motion states of the asymmetric system under different parameter settings.
The effects of the initial perturbation are investigated first. Figure 3(a) shows the time history of the filament lateral velocity
$V$
for various values of perturbation intensity
$\beta$
in the case
$\delta = 20^\circ$
,
$K=0.0001$
and
$\textit{Re} = 10$
. It can be seen that after sufficient time has elapsed, the influence of the initial perturbation gradually vanishes, and the values of
$V$
for different
$\beta$
eventually become identical and time-invariant. The corresponding instantaneous vorticity contours and streamlines are presented in figure 3(b). Shear layers of opposite-sign vorticity are generated at the two ends of the filament and extend downstream, forming a steady wake that is symmetric about the horizontal centreline (see supplementary movies available at https://doi.org/10.1017/jfm.2026.12215). This symmetry implies that the filament lateral velocity equals the
$y$
-component of the incoming flow velocity, namely
$V = U_{\infty ,y} = U_\infty \sin \delta = 0.017$
, as shown in figure 3(a). Apart from the non-zero lateral velocity, these wake characteristics are consistent with those of the static mode at
$\delta = 0^\circ$
. Therefore, this state is also referred to here as the static mode.
The evolution of
$V$
and the corresponding flow structures for the case with
$\delta = 20^\circ$
,
$K = 0.01$
and
$\textit{Re} = 50$
are shown in figures 3(c,d). As shown in figure 3(c), although a phase difference exists, the mean amplitude and period of
$V$
are identical for different
$\beta$
, indicating that
$\beta$
does not influence the equilibrium kinematics of the filament. Positive and negative vortices are shed alternately from the filament, forming a Kármán vortex street (see figure 3(d) or the supplementary movies). This periodic vortex shedding is accompanied by periodic oscillations of
$V$
. As in the static mode, the time-averaged lateral velocity of the filament,
$\bar V=\langle V\rangle$
, equals the transverse component of the incoming flow velocity, namely
$\bar V = U_{\infty ,y} = U_\infty \sin \delta = 0.086$
(see figure 3
c). Consequently, the wake remains symmetric about the horizontal axis behind the filament. Based on these characteristics, this case is classified as a vibration mode.
Time history of the lateral velocity of the filament under various initial disturbance intensities
$\beta$
for (a) the static mode
$\delta = 20^\circ$
,
$K=0.0001$
,
$\textit{Re} = 10$
, and (c) the vibrational mode
$\delta = 20^\circ$
,
$K=0.01$
,
$\textit{Re} = 50$
. The motion mode at the equilibrium state are identical for different
$\beta$
. (b,d) The corresponding instantaneous vorticity contours and streamlines. Animated visualisations are also provided in the supplementary movies.
Time history of the lateral velocity of the filament under various initial disturbance intensities
$\beta$
for (a) the quasi-static mode
$\delta = 40^\circ$
,
$K=0.0001$
,
$\textit{Re} = 50$
, and (d) the translational mode
$\delta = 30^\circ$
,
$K=0.01$
,
$\textit{Re} = 100$
. The motion mode at the equilibrium state splits into two branches for different
$\beta$
, i.e. the upwind branch (I) and the downwind branch (II). Panels (b,e) and (c, f) present the corresponding instantaneous vorticity contours and streamlines for the two branches. Animated visualisations are also provided in the supplementary movies.
Figure 4. Long description
Panel A: A line graph shows the time history of the lateral velocity of the filament under various initial disturbance intensities for the quasi-static mode. The x-axis represents time normalized by a reference period (t/Tref), and the y-axis represents the lateral velocity (V). Different lines represent different beta values (β = 0.3, β = 0, β = −0.5, β = −0.7). The graph splits into two branches for different beta values: the upwind branch (I) and the downwind branch (II). Panel B: A contour plot shows the instantaneous vorticity contours and streamlines for the upwind branch (I) of the quasi-static mode. The x-axis represents the horizontal position (x), and the y-axis represents the vertical position (y). Panel C: A contour plot shows the instantaneous vorticity contours and streamlines for the downwind branch (II) of the quasi-static mode. The x-axis represents the horizontal position (x), and the y-axis represents the vertical position (y). Panel D: A line graph shows the time history of the lateral velocity of the filament under various initial disturbance intensities for the translational mode. The x-axis represents time normalized by a reference period (t/Tref), and the y-axis represents the lateral velocity (V). Different lines represent different beta values (β = 1, β = 0, β = −1, β = −1.5). The graph splits into two branches for different beta values: the upwind branch (I) and the downwind branch (II). Panel E: A contour plot shows the instantaneous vorticity contours and streamlines for the upwind branch (I) of the translational mode. The x-axis represents the horizontal position (x), and the y-axis represents the vertical position (y). Panel F: A contour plot shows the instantaneous vorticity contours and streamlines for the downwind branch (II) of the translational mode. The x-axis represents the horizontal position (x), and the y-axis represents the vertical position (y).
Unlike the two wake modes mentioned above, for some cases, variations in
$\beta$
can substantially alter the filament dynamics. A representative example is the case with
$\delta = 40^\circ$
,
$K = 0.0001$
and
$\textit{Re} = 50$
, shown in figures 4(a–c). It is seen that the mean lateral velocity
$\bar V$
is 0.046 for
$\beta = 0$
and 0.3, but increases to 0.275 for
$\beta = -0.5$
and
$-0.7$
(see figure 4
a). In other words, depending on the initial disturbance level
$\beta$
, the filament may move either upstream or downstream relative to the inclined inflow. In these cases, however, the filament reaches a constant equilibrium velocity, and no vortex shedding occurs (see figures 4(b,c) or the supplementary movies). This regime is therefore classified as a quasi-static mode. For clarity, the states with lower and higher velocity magnitudes are denoted QS-I and QS-II (as marked in figure 4
a), respectively. Notably, in the QS-I mode, the mean lateral velocity of the filament
$\bar V$
is smaller than the lateral component of the incoming flow
$U_{\infty ,y}=0.16$
. The filament therefore moves upwind relative to the surrounding fluid (see figure 4
b), which may appear counterintuitive. In contrast, in the QS-II mode,
$\bar V$
exceeds
$U_{\infty ,y}$
, indicating that the filament moves downwind relative to the fluid. Consequently, the effective relative flows experienced by the filament in the two modes have opposite lateral components, resulting in oppositely inclined wakes, as shown in figures 4(b) and 4(c).
A similar phenomenon is observed for
$\delta = 30^\circ$
,
$K = 0.01$
and
$\textit{Re} = 100$
. Figures 4(d–f) show that the filament translates laterally while shedding vortices of opposite sign alternately, thereby forming an inclined Kármán vortex street (see figures 4(e, f) or the supplementary movies), which is characteristic of the translational mode. Varying
$\beta$
also gives rise to two distinct motion branches. Specifically, for
$\beta = 0$
and 1, the filament exhibits a negative mean lateral velocity
$\bar V = -0.4$
, corresponding to upstream motion against the incoming flow (see figure 4
e); this branch is termed mode T-I (as marked in figure 4
d). By contrast, for
$\beta = -1$
and
$-1.5$
, the filament undergoes the expected downstream motion, with
$\bar V = 0.9$
(see figures 4
d, f); this branch is termed mode T-II (as marked in figure 4
d). Notably, the vibration frequency is the same in both sub-modes (see figure 4
d).
It is worth noting that in the absence of an initial disturbance (i.e.
$\beta = 0$
), the filament is more likely to undergo upwind motion, as shown by the blue solid lines in figures 4(a,d). Downwind motion occurs only when the initial disturbance imparts a sufficiently large forward velocity to the filament, as in the cases with
$\beta = -0.7$
in figure 4(a) and
$\beta = -1.5$
in figure 4(d). This behaviour is rooted in the initial stability of the system, and will be explained in detail in
$\S$
4.4.
Time-averaged lateral velocity
$\bar V$
of the filament as functions of
$\textit{Re}$
for different
$\delta$
, with (a)
$K = 10^{-4}$
, (b)
$K = 10^{-2}$
and (c)
$K = 1$
.
To comprehensively investigate the bifurcation behaviour of filaments in inclined crossflow, figure 5 presents the time-averaged lateral velocity
$\bar V$
as a function of
$\textit{Re}$
for different
$\delta$
and
$K$
. In bifurcating cases, the solid and dashed lines represent the upwind (I) and downwind (II) branches, respectively. For
$K=10^{-4}$
in figure 5(a), only a single value of
$\bar V$
is observed at low Reynolds number (e.g.
$\textit{Re}=10$
). As
$\textit{Re}$
increases, the magnitude of
$\bar V$
gradually increases, and a dual-branch structure emerges, corresponding to the upwind and downwind motions of the filament. This is characteristic of a supercritical pitchfork bifurcation, in which bifurcation appears beyond a critical point (Strogatz Reference Strogatz2018). A similar bifurcation pattern is observed for
$K=0.01$
, as shown in figure 5(b). Moreover, except for the cases with
$\delta = 0$
, the two branches exhibit different velocity magnitudes:
$|\bar V|$
is larger for the downwind branch, and smaller for the upwind branch (see figures 5
a,b). When
$K$
increases to 1 (see figure 5
c), however, the reduced filament deformation suppresses the bifurcation, yielding a unique
$\bar V$
in each case.
Relative lateral velocity
$\bar V_r$
of the filament as functions of
$\textit{Re}$
for various
$\delta$
, with (a)
$K = 10^{-4}$
, (b)
$K = 10^{-2}$
and (c)
$K = 1$
.
The results in figure 3 also indicate that the mean filament velocity
$\bar V$
is closely related to the transverse inflow velocity
$U_{\infty ,y}$
. This observation motivates analysis of the filament kinematics in a reference frame translating laterally with the fluid. To this end, the relative lateral velocity
$\bar V_r$
is defined as
\begin{equation} \bar V_r = \bar V - U_{\infty ,y} = \bar V - U_{\infty }\sin \delta . \end{equation}
Figure 6 shows
$\bar V_r$
as a function of
$\textit{Re}$
for various
$\delta$
and
$K$
. Relative to the results in figure 5, the two branches recover symmetry in this moving frame. In particular, for small
$\textit{Re}$
(see figures 6
a,b) or large
$K$
(see figure 6
c),
$\bar V_r$
remains zero over a finite range of initial disturbances, indicating synchronous lateral motion of the filament and the surrounding fluid. As
$\textit{Re}$
increases beyond the critical value, bifurcation occurs, giving rise to two branches of equal magnitude in
$\bar V_r$
(see figures 6
a,b). Therefore, from the perspective of the moving frame, the upwind and downwind states are indistinguishable, implying that QS-I and QS-II, as well as T-I and T-II, are dynamically equivalent. Consequently, the filament in inclined crossflow can be classified into four dynamic modes: static (S), quasi-static (QS), translational (T), and vibrational (V), consistent with the results with
$\delta = 0$
(Liu et al. Reference Liu, Gao and Lu2026). Note that due to the relatively low
$\textit{Re}$
considered in the present study, no chaotic mode was observed within the parameter range investigated.
Figure 7 summarises the four motion states in the
$(\delta , Re, K)$
parameter space. It is observed that, the S mode occurs mainly at low Reynolds numbers, for example, at
$\textit{Re} = 10$
for
$K = 1$
. For
$K = 10^{-4}$
, increasing
$\textit{Re}$
leads to the sequential appearance of the S, QS and T modes. For
$K = 10^{-2}$
, the QS mode is replaced by the V mode, whereas for
$K = 1$
, only the V mode is observed at high Reynolds numbers. Although these modes can be distinguished directly from the flow field structures, quantitative indicators are introduced below to enable a more objective characterisation.
Three-dimensional phase diagram for the four motion states of the filament in the
$(\delta , Re, K)$
parameter space. Symbols
$\circ$
,
$\bigtriangleup$
,
$\bigtriangledown$
and
$\square$
represent the static, quasi-static, translational and vibrational modes, respectively.
Considering that the magnitudes of the relative velocity
$|\bar V_r|$
are identical for different branches in bifurcation cases, a relative translational Reynolds number
$\textit{Re}_r$
is introduced based on
$|\bar V_r|$
:
\begin{equation} \textit{Re}_r = \frac { \rho\, |\bar V_r|\, L }{ \mu }. \end{equation}
The Strouhal number
$St$
is defined as
\begin{equation} \textit{St} = \frac { f_1 L }{ U_{\infty , x} } , \end{equation}
where
$f_1$
denotes the filament vibration frequency, and
$U_{\infty ,x} = U_\infty \cos \delta$
is the streamwise component of the incoming flow velocity (see figure 1). To incorporate the effect of the inflow inclination angle
$\delta$
, the streamwise Reynolds number is defined as
\begin{equation} \textit{Re}_x = \frac { \rho U_{\infty , x} L }{ \mu } = \frac { \rho ( U_{\infty } \cos \delta ) L }{ \mu } = Re \cos \delta . \end{equation}
The wake inclination angle
$\theta$
is then calculated as
\begin{equation} \theta = \arctan \left( \frac { |\bar{V}_r| }{ U_{\infty , x} }\right)\!. \end{equation}
Figures 8, 9 and 10 present
$\textit{Re}_r$
,
$St$
and
$\theta$
as functions of
$\textit{Re}_x$
for various
$\delta$
and
$K$
, respectively. For each fixed
$K$
, the data corresponding to different
$\delta$
collapse onto a single curve, indicating that
$\textit{Re}_x$
serves as a more suitable control parameter and a more robust descriptor of the system than
$\textit{Re}$
.
The relative translational Reynolds number
$\textit{Re}_{r}$
as functions of the streamwise Reynolds number
$\textit{Re}_x$
for different
$\delta$
, with (a)
$K = 10^{-4}$
, (b)
$K = 10^{-2}$
and (c)
$K = 1$
.
Figure 8. Long description
Panel A: A scatter plot shows the relationship between the relative translational Reynolds number (Re_r) and the streamwise Reynolds number (Re_x) for a specific condition. The x-axis represents Re_x on a logarithmic scale ranging from 10^1 to 10^2, and the y-axis represents Re_r ranging from 0 to 250. The plot includes data points marked with green diamonds and a dashed line indicating a proportional relationship (Re_r ~ Re_x). The plot is divided into three regions labeled S, QS, and T. Panel B: Another scatter plot depicts the relationship between Re_r and Re_x for a different condition. The x-axis and y-axis are the same as in Panel A. The data points are marked with blue diamonds, and the plot also includes a dashed line indicating a proportional relationship (Re_r ~ Re_x). The plot is divided into three regions labeled S, V, and T. Panel C: The third scatter plot shows the relationship between Re_r and Re_x for various angles (δ) ranging from 0° to 50°. The x-axis and y-axis are the same as in the previous panels. The data points are marked with different symbols and colors representing different angles. The plot is divided into two regions labeled S and V.
The Strouhal number
$St$
as functions of the streamwise Reynolds number
$\textit{Re}_x$
for different
$\delta$
, with (a)
$K = 10^{-4}$
, (b)
$K = 10^{-2}$
and (c)
$K = 1$
.
The inclination angle
$\theta$
of the wake as functions of the streamwise Reynolds number
$\textit{Re}_x$
for different
$\delta$
, with (a)
$K = 10^{-4}$
, (b)
$K = 10^{-2}$
and (c)
$K = 1$
.
These statistics provide a quantitative basis for distinguishing the different modes, which are also highlighted by the background colours in figures 8–10. Specifically, for
$K = 10^{-4}$
, when
$\textit{Re}_x \leqslant 25$
, all three quantities, i.e.
$\textit{Re}_r$
,
$St$
and
$\theta$
, vanish; see figures 8(a), 9(a) and 10(a). This indicates that the filament undergoes neither relative lateral translation nor vibration, while the wake remains horizontal (see figure 3
b), characteristic of the S mode. Similar results are obtained for
$K = 10^{-2}$
and
$K = 1$
when
$\textit{Re}_x \leqslant 20$
; see figures 8(b,c), 9(b,c) and 10(b,c).
When
$\textit{Re}_x$
is moderate for
$K = 10^{-4}$
, i.e.
$25 \lt \textit{Re}_x \leqslant 40$
,
$\textit{Re}_r$
and
$\theta$
increase continuously from zero (see figures 8(a) and 10(a)), indicating the onset of relative lateral translation and the development of an inclined wake (see figures 4
b,c). Since
$St$
remains zero, the filament does not vibrate. These cases are therefore classified as the QS mode. In contrast, for
$K = 10^{-2}$
and
$K = 1$
at moderate
$\textit{Re}_x$
(i.e.
$20 \lt \textit{Re}_x \lt 55$
for
$K = 10^{-2}$
, and
$\textit{Re}_x \gt 20$
for
$K = 1$
), the QS mode is absent and replaced by the V mode, in which
$\textit{Re}_r$
and
$\theta$
remain zero, whereas
$St$
becomes clearly non-zero; figures 8(b,c), 9(b,c) and 10(b,c). This suggests that the filament undergoes vibration only, with no relative lateral translation, while the wake remains horizontal (see figure 3
d).
As
$\textit{Re}_x$
is further increased at
$K = 10^{-4}$
, i.e.
$\textit{Re}_x \gt 40$
in figures 8(a), 9(a) and 10(a), the filament enters the T mode. In this regime, the non-zero values of both
$\textit{Re}_r$
and
$St$
indicate that the filament undergoes lateral translation and vibration simultaneously. The wake also becomes inclined, as reflected by
$\theta \gt 0$
(see figure 10
a). Similar kinematic features are observed for
$K = 0.01$
; figures 8(b), 9(b) and 10(b). Moreover, for these translational cases,
$\textit{Re}_r$
obeys a simple scaling relation with
$\textit{Re}_x$
, namely,
$\textit{Re}_r \sim \textit{Re}_x$
, which is in agreement with our previous findings (Liu et al. Reference Liu, Gao and Lu2026).
4.2. Flow structure and the unsteady dynamics
To clarify the mechanism underlying the free motion of filaments in crosswind, in this section, we examine the vortical structures around the filaments and pressure distributions on the surfaces.
Figure 11 presents the time histories of the streamwise and lateral forces, i.e.
$F_x$
and
$F_y$
, for the S mode (
$\delta = 20^\circ$
,
$K = 0.0001$
,
$\textit{Re} = 10$
) and the QS mode (
$\delta = 40^\circ$
,
$K = 0.0001$
,
$\textit{Re} = 50$
), together with the corresponding instantaneous pressure contours. It is seen that in both modes, the forces remain time-invariant. In the S mode, the pressure distribution around the filament is symmetric about the horizontal axis (see the inset in figure 11
a), consistent with the vorticity contours shown in figure 3(b). This symmetry further implies that the filament undergoes no relative lateral motion (see figure 3
a), while for the QS mode shown in figure 11(b), the pressure symmetry is broken, resulting in a steady lateral velocity of the filament (see figure 4
a).
Time history of the streamwise force
$F_x$
and lateral force
$F_y$
for (a) the static mode
$\delta = 20^\circ$
,
$K=0.0001$
,
$\textit{Re} = 10$
, and (b) the quasi-static mode
$\delta = 40^\circ$
,
$K=0.0001$
,
$\textit{Re} = 50$
. Insets: the corresponding instantaneous pressure contours for the two modes.
In contrast to the two time-invariant modes discussed above, the V and T modes exhibit time-dependent forces on the filament due to periodic vortex shedding, as shown in figure 12. In these two modes,
$F_x$
and
$F_y$
fluctuate approximately sinusoidally. The filament experiences a net streamwise drag, while the cycle-averaged lateral force remains zero. In the V mode, the frequency of
$F_x$
is twice that of
$F_y$
, whereas in the translational mode, the two forces have the same frequency but are in antiphase (see figure 12). This distinction can be attributed to the different vortex structures associated with the two modes.
Time history of the streamwise force
$F_x$
and lateral force
$F_y$
for (a) the vibrational mode
$\delta = 20^\circ$
,
$K=0.01$
,
$\textit{Re} = 50$
, and (b) the translational mode
$\delta = 30^\circ$
,
$K=0.01$
,
$\textit{Re} = 100$
. Here,
$t_1{-}t_8$
are some typical instants, which will be discussed in figures 13 and 14.
(a–d) Instantaneous vorticity contours and (e–h) instantaneous pressure contours at
$t/T_{\textit {ref}} = t_1$
,
$t_2$
,
$t_3$
and
$t_4$
for the cases with
$\delta = 20^\circ$
,
$K=0.01$
,
$\textit{Re} = 50$
, corresponding to the vibrational mode. The four instants are labelled
$t_1{-}t_4$
in figure 12(a). Here, A and B mark the vortices.
Figure 13. Long description
Panel A: A contour plot showing instantaneous vorticity at a specific time instance. The x-axis ranges from −1 to 2, and the y-axis ranges from −1 to 1. The plot features two marked vortices labeled A and B. Panel B: Another contour plot of instantaneous vorticity at a different time instance, with the same axis ranges. Panel C: A contour plot of instantaneous vorticity at yet another time instance, maintaining the same axis ranges. Panel D: A contour plot of instantaneous vorticity at the final time instance, with the same axis ranges. Panel E: A contour plot showing instantaneous pressure at the first time instance, with the x-axis ranging from −1 to 2 and the y-axis ranging from −1 to 1. Panel F: Another contour plot of instantaneous pressure at the second time instance, with the same axis ranges. Panel G: A contour plot of instantaneous pressure at the third time instance, maintaining the same axis ranges. Panel H: A contour plot of instantaneous pressure at the final time instance, with the same axis ranges.
For the V mode, figure 13 shows the instantaneous vorticity and pressure contours at selected instants marked in figure 12(a). It is observed that at
$t/T_{\textit {ref}} = t_1$
, vortex A generates a low-pressure region behind the upper end of the filament and a high-pressure region immediately upstream (see figure 13
e), producing positive
$F_x$
and
$F_y$
(see figure 12
a). At
$t/T_{\textit {ref}} = t_2$
, both vortex A and the associated low-pressure region intensify significantly (see figures 13
b, f), resulting in a marked increase in
$F_x$
and
$F_y$
, with
$F_y$
attaining its maximum. At
$t/T_{\textit {ref}} = t_3$
, vortex A convects downwards and downstream relative to the filament (see figure 13
c), while the low-pressure region expands substantially (figure 13
g). This evolution strengthens
$F_x$
but weakens
$F_y$
, so that
$F_x$
reaches its maximum, whereas
$F_y$
decreases at this instant (see figure 12
a).
At
$t/T_{\textit {ref}} = t_4$
, vortex A convects farther downstream, and the associated low-pressure region is located relatively far from the filament (see figures 13
d,h). This weakens the pressure difference across the filament and hence reduces
$F_x$
(see figure 12
a). At the same time, vortex B near the lower end of the filament begins to strengthen, producing a low-pressure region there (see figures 13
d,h). The resulting downward suction force further decreases
$F_y$
(figure 12
a).
(a–d) Instantaneous vorticity contours and (e–h) instantaneous pressure contours at
$t/T_{\textit {ref}} = t_5$
,
$t_6$
,
$t_7$
and
$t_8$
for the cases with
$\delta = 30^\circ$
,
$K=0.01$
,
$\textit{Re} = 100$
, corresponding to the translational mode. The four instants are labelled
$t_5{-}t_8$
in figure 12(b). Here, A,
${\rm A}^\prime$
, B and
${\rm B}^\prime$
mark the vortices.
In the second half of the cycle, vortex B evolves in a manner analogous to vortex A, and affects the filament forces in the same way; the corresponding details are therefore omitted for brevity. These observations also clarify the frequency difference between
$F_x$
and
$F_y$
: the shedding of one vortex produces one full cycle in
$F_x$
, but only half a cycle in
$F_y$
.
Figure 14 presents the instantaneous vorticity and pressure contours for the T mode. It is seen that at
$t/T_{\textit {ref}} = t_5$
, vortex A induces a strong low-pressure region above the filament (see figure 14
e), generating a pronounced lateral attraction that causes
$F_y$
to reach its maximum value (see figure 12
b). By contrast, because the horizontal pressure difference across the filament is small at this instant (see figure 14
e),
$F_x$
attains its minimum value (see figure 12
b). At
$t/T_{\textit {ref}} = t_6$
, vortex A convects farther from the filament (see figure 14
b), reducing its influence and thereby decreasing
$F_y$
(see figure 12
b). At the same time, vortex B splits to form vortex
$\mathrm{B}^\prime$
, which creates a low-pressure region behind the filament (see figure 14
f). The resulting increase in pressure difference across the filament causes
$F_x$
to rise (figure 12
b).
At
$t/T_{\textit {ref}} = t_7$
, vortex A has convected sufficiently far from the filament that its influence becomes negligible (see figure 14
c). In contrast, vortex
$\mathrm{B}^\prime$
strengthens further, producing a large low-pressure region directly behind the filament, especially near the lower end, while the high-pressure region upstream is also intensified (see figure 14
g). Together, these pressure distributions generate strong backward and downward pressure differences, causing
$F_x$
and
$F_y$
to attain their maximum and minimum values, respectively (see figure 12
b).
At
$t/T_{\textit {ref}} = t_8$
, vortex
$\mathrm{A}^\prime$
, generated by the splitting of vortex A, induces a low-pressure region above the filament that offsets the low-pressure region below it (see figure 14
h), so that
$F_y \approx 0$
at this instant (see figure 12
b). At the same time, as vortex
$\mathrm{B}^\prime$
moves farther downstream, the low-pressure region behind the filament weakens (see figure 14
h), leading to a decrease in
$F_x$
(figure 12
b).
The above analysis suggests that in the T mode,
$F_x$
and
$F_y$
are controlled primarily by vortex A (
$\mathrm{A}^\prime$
) generated at the upper end of the filament, and vortex B (
$\mathrm{B}^\prime$
) generated at the lower end, respectively. The shedding frequencies of vortices A and B are identical; consequently, the corresponding frequencies of
$F_x$
and
$F_y$
coincide. In addition, the alternate shedding of vortices A and B gives rise to a phase difference
$\pi$
between
$F_x$
and
$F_y$
, i.e. the two forces are in antiphase (see figure 12
b).
4.3. Force decomposition and scaling laws
Deformation in flexible bodies redistributes aerodynamic or hydrodynamic forces and thus governs the filament’s mode of motion. Hence the forces acting on the filament are analysed below to elucidate the mechanism of its spontaneous motion.
Schematic diagram for force decomposition. Here,
$\boldsymbol{\tau }$
and
$\boldsymbol{n}$
are the local tangential and normal vectors, respectively, and
$\alpha$
is the angle between the local normal direction and the horizontal direction.
Figure 15. Long description
A diagram showing the decomposition of forces in a flapping system. The diagram includes vectors representing different force components. The local tangential vector is labeled as τ, and the local normal vector is labeled as n. The angle α is shown between the local normal direction and the horizontal direction. The force components are decomposed into normal and tangential directions, with fn representing the normal force, fτ representing the tangential force, and their respective components in the x and y directions. The diagram illustrates the relationships and directions of these forces in the context of a flapping system.
As shown in figure 15, the fluid force acting on the filament,
$\boldsymbol f_s$
, can be decomposed into a normal component
$\boldsymbol f_n$
, and a tangential component
$\boldsymbol f_{\tau }$
. The former is dominated by pressure, whereas the latter mainly originates from viscous effects. These forces are defined as follows (Peng et al. Reference Peng, Huang and Lu2018; Liu et al. Reference Liu, Huang and Lu2020, Reference Liu, Liu and Huang2022):
\begin{equation} \boldsymbol{f}_{\!s}=[-p \boldsymbol{I}+\boldsymbol{T}_{\!\mu} ] \boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{f}_{\!n}+\boldsymbol{f}_{\!\tau }, \end{equation}
\begin{equation} \boldsymbol{f}_{\!n}=\left (\boldsymbol{f}_{\!s} \boldsymbol{\cdot }\boldsymbol{n}\right ) \boldsymbol{n}= f_n \boldsymbol{n} = \left (f_{n, x}, f_{n, y} \right )\!, \end{equation}
\begin{equation} \boldsymbol{f}_{\!\tau }=\left (\boldsymbol{f}_{\!s} \boldsymbol{\cdot }\boldsymbol{\tau }\right ) \boldsymbol{\tau }= f_{\tau } \boldsymbol{\tau } = \left (f_{\tau , x}, f_{\tau , y} \right )\!, \end{equation}
where
$\boldsymbol{I}$
denotes the identity tensor,
$\boldsymbol{T}_{\!\mu}$
is the viscous stress tensor, and
$\boldsymbol{\tau }$
and
$\boldsymbol{n}$
are the unit tangential and normal vectors, respectively. The notation
$[{\cdot }]$
denotes the jump in a quantity across the immersed boundary. Furthermore,
$f_n$
and
$f_\tau$
represent the magnitudes of
$\boldsymbol f_n$
and
$\boldsymbol f_{\tau }$
, respectively, whereas
$f_{n,x}$
and
$f_{n,y}$
(
$f_{\tau ,x}$
and
$f_{\tau ,y}$
) are the corresponding components in the
$x$
- and
$y$
-directions. The time-averaged components are defined as
$\bar f_{n,x}=\langle f_{n,x}\rangle$
,
$\bar f_{n,y}=\langle f_{n,y}\rangle$
,
$\bar f_{\tau ,x}=\langle f_{\tau ,x}\rangle$
and
$\bar f_{\tau ,y}=\langle f_{\tau ,y}\rangle$
.
Distributions of (a) the
$x$
-component of the normal force
$\bar f_{n,x}$
, (b) the
$x$
-component of the tangential force
$\bar f_{\tau ,x}$
, (c) the
$y$
-component of the normal force
$\bar f_{n,y}$
, and (d) the
$y$
-component of the tangential force
$\bar f_{\tau ,y}$
along the filament for the cases with
$\textit{Re} = 50, 80, 100, 150$
,
$\delta = 30^\circ$
,
$K=0.01$
and
$\beta = 0$
.
Figure 16. Long description
Panel A: A line graph shows the distribution of the x-component of the normal force along the filament. The x-axis is labeled s/L and ranges from 0 to 1. The y-axis is labeled f-n,x and ranges from 0 to 0.03. Different symbols represent different Reynolds numbers: circles for Re = 50, upward triangles for Re = 80, diamonds for Re = 100, and downward triangles for Re = 150. The graph shows that the force increases to a peak and then decreases as s/L increases. Panel B: A line graph shows the distribution of the x-component of the tangential force along the filament. The x-axis is labeled s/L and ranges from 0 to 1. The y-axis is labeled f-τ,x and ranges from −6 x 10^−3 to 4 x 10^−3. The graph shows fluctuations around zero with different patterns for each Reynolds number. Panel C: A line graph shows the distribution of the y-component of the normal force along the filament. The x-axis is labeled s/L and ranges from 0 to 1. The y-axis is labeled f-n,y and ranges from −0.02 to 0.02. The graph shows that the force varies with s/L, with different patterns for each Reynolds number. Panel D: A line graph shows the distribution of the y-component of the tangential force along the filament. The x-axis is labeled s/L and ranges from 0 to 1. The y-axis is labeled f-τ,y and ranges from −5 x 10^−3 to 5 x 10^−3. The graph shows fluctuations around zero with different patterns for each Reynolds number.
To clarify the characteristics of these force components, figure 16 shows the distributions of
$\bar f_{n,x}$
,
$\bar f_{\tau ,x}$
,
$\bar f_{n,y}$
and
$\bar f_{\tau ,y}$
along the filament for four representative cases with
$\textit{Re} = 50$
, 80, 100 and 150, at
$\delta = 30^\circ$
,
$K = 0.01$
and
$\beta = 0$
. It is seen that both tangential force components, i.e.
$\bar f_{\tau ,x}$
and
$\bar f_{\tau ,y}$
, are negligibly small. In contrast, the normal force exhibits a distinct spatial variation. As
$s$
increases,
$\bar f_{n,x}$
first increases and then decreases, reaching a maximum near the midpoint
$s/L = 0.5$
(see figure 16
a). This behaviour results from the approximate up–down symmetry of the filament. Near the midpoint, the angle between the local normal direction and the horizontal direction,
$\alpha$
, is approximately zero, whereas its magnitude increases progressively towards both ends. Since
$f_{n,x}=f_n\cos \alpha$
and
$f_n$
is approximately uniform along the filament (Liu et al. Reference Liu, Gao and Lu2026),
$\bar f_{n,x}$
is largest near the centre and smaller near the endpoints. Conversely,
$|\bar f_{n,y}|$
is smallest at the midpoint and becomes larger towards both ends (see figure 16
c), owing to the relation
$f_{n,y} = f_n \sin \alpha$
.
Note that for the cases with
$\textit{Re}=80$
, 100 and 150 at
$\beta = 0$
, the filament ultimately undergoes upwind motion, namely
$\bar V_r \lt 0$
(see figure 4
d). In these cases,
$\bar f_{n,y}\lt 0$
over
$s/L \in [0,0.5]$
(see figure 16
c), i.e. the normal force contributes to the lateral thrust in this region because it is in the same direction as the relative velocity
$\bar V_r$
, whereas
$\bar f_{n,y}\gt 0$
over
$s/L \in [0.5,1]$
(figure 16
c), implying a contribution to the lateral drag. The lateral thrust and drag contributions of normal force can be quantified as
\begin{equation} F_{n,y}^{(t)} = \left | \int _0^{L/2} \bar f_{n,y}\, \mathrm{d} s \right |\!, \quad F_{n,y}^{(d)} = \left | \int _{L/2}^{L} \bar f_{n,y}\, \mathrm{d} s \right |\!, \end{equation}
where the superscripts
$(t)$
and
$(d)$
denote thrust and drag, respectively. By contrast, the tangential force contributes predominantly to the lateral drag, since
$\bar f_{\tau ,y}\gt 0$
over most of the filament except near the endpoints (see figure 16
d). The corresponding lateral drag contributions of tangential force can be calculated as
\begin{equation} F_{\tau ,y} = \int _{0}^{L} \bar f_{\tau ,y}\, \mathrm{d} s. \end{equation}
Likewise, in the horizontal direction, the total drag
$F_d$
and the drag contributions associated with the normal and tangential forces are expressed as
\begin{equation} F_{d} = \langle F_x \rangle , \quad F_{n,x} = \int _{0}^{L} \bar f_{n,x}\, \mathrm{d} s, \quad F_{\tau ,x} = \int _{0}^{L} \bar f_{\tau ,x}\, \mathrm{d} s. \end{equation}
Figure 17 presents these forces, namely
$F_{d}$
,
$F_{n,x}$
,
$F_{\tau ,x}$
,
$F_{n,y}^{(t)}$
,
$F_{n,y}^{(d)}$
and
$F_{\tau ,y}$
, for the four representative cases. It is seen that the horizontal drag is contributed almost entirely by the normal force, since
$F_d \approx F_{n,x}$
and
$F_{\tau ,x} \approx 0$
(see figure 17
a). In the lateral direction, for the translational cases (i.e.
$\textit{Re} \gt 50$
),
$F_{n,y}^{(t)} \gt F_{n,y}^{(d)}$
(see figure 17
b), meaning that the lateral thrust generated by the normal force exceeds the corresponding drag. This result reflects the up–down asymmetry of the normal force distribution induced by symmetry breaking, as shown in figure 16(c), where
$|\bar f_{n,y}|$
near the lower endpoint (
$s=0$
) is greater than that near the upper endpoint (
$s=L$
).
At equilibrium, the resultant lateral force on the filament is zero. It follows that
$F_{n, y}^{(t)} = F_{n, y}^{(d)} + F_{\tau , y}$
, as evidenced in figure 17(b). This formula can be further written as
\begin{equation} F_{\tau , y} = F_{n, y}^{(t)} - F_{n, y}^{(d)} = F_{n,y} = \left | \int _{0}^{L} \bar{f}_{n,y}\, \mathrm{d} s \right |\!, \end{equation}
where
$F_{n,y}$
denotes the net thrust arising from the normal force. Equation (4.13) therefore characterises the balance between the lateral thrust and drag of the filament.
(a) The drag
$F_d$
, the total
$x$
-component of the normal and tangential forces
$\boldsymbol f_n$
and
$\boldsymbol f_\tau$
(i.e.
$F_{n, x}$
and
$F_{\tau , x}$
), and (b) the lateral thrust and drag contributed by the
$y$
-component of the normal force
$\boldsymbol f_n$
(i.e.
$F_{n, y}^{(t)}$
and
$F_{n, y}^{(d)}$
), and the total
$y$
-component of the tangential force
$\boldsymbol f_\tau$
(i.e.
$F_{\tau , y}$
) for the cases with
$\textit{Re} = 50, 80, 100, 150$
,
$\delta = 30^\circ$
,
$K=0.01$
and
$\beta = 0$
.
These forces are intimately linked to filament motion and deformation, and their scaling behaviour therefore merits further examination. The drag in the horizontal direction is considered first. For a typical fixed rigid bluff body, the drag scales with the square of the oncoming flow velocity
$U_\infty$
(Batchelor Reference Batchelor1967). For a laterally unconstrained flexible filament, however, the relevant velocity is the effective velocity
$U_{\textit {eff}}$
(Liu et al. Reference Liu, Gao and Lu2026). Under the present inclined inflow,
$U_{\textit {eff}}$
is defined as
$U_{\textit {eff}}=\sqrt {U_{\infty ,x}^2+\bar V_r^2}$
. Moreover, owing to the pronounced bending deformation, the effective length for fluid–filament interaction is
$L_{\textit {eff}}=L_y$
rather than
$L$
, where
$L_y$
denotes the lateral projected length of the filament. Accordingly, the effective Reynolds number is defined as
\begin{equation} \textit{Re}_{\textit {eff}} = \frac { \rho U_{\textit {eff}} L_{\textit {eff}} }{\mu } = \frac {\rho \sqrt { U_{\infty ,x}^2 + \bar V_r^2 }\, L_y}{\mu }. \end{equation}
As shown in figure 18(a),
$F_{n,x}$
follows a clear
$\textit{Re}_{\textit {eff}}^2$
scaling over a broad range of Reynolds numbers, bending stiffnesses and inflow inclination angles for all four motion modes. This generalises the findings of Liu et al. (Reference Liu, Gao and Lu2026) from the symmetric inflow case to the more general inclined inflow configuration.
(a) The total
$x$
-component of the normal force
$F_{n, x}$
as a function of the effective Reynolds number
$\textit{Re}_{\textit{eff}}$
, and (b) the total
$y$
-component of the normal force
$F_{n, y}$
as a function of
$|\bar V_r|\, f_1 / \kappa _c^2$
for all cases.
Next, we focus on the lateral force. For the lateral thrust, the dominant contribution originates from the pressure difference across the filament. In a frame moving laterally with the oncoming flow, filament vibration and vortex shedding induce a small velocity
$u$
on the two sides of the filament, resulting in local flow velocities of
$\bar V_r + u$
and
$\bar V_r - u$
near the windward and leeward sides, respectively. By Bernoulli’s principle, the resulting unsteady thrust can be estimated as
\begin{equation} F_{n,y} \sim \rho \left [ (\bar V_r + u)^2 - (\bar V_r - u)^2 \right ] \mathcal{L} \sim u\, |\bar V_r|\, \mathcal{L}, \end{equation}
where
$\mathcal{L}$
is the characteristic horizontal length scale. Considering filament bending,
$\mathcal{L}$
is taken as
$1/\kappa _c$
, where
$\kappa _c=\left \langle ({1}/{L})\int _0^L \kappa\, \mathrm{d}s \right \rangle$
is the chord-averaged curvature, and
$\kappa$
is the local curvature. Since the characteristic time scale of filament vibration and vortex shedding is
$\mathcal{T}=1/f_1$
, the induced velocity scales as
\begin{equation} u \sim \mathcal{L}/\mathcal{T}=f_1/\kappa _c. \end{equation}
Hence (4.15) becomes
\begin{equation} F_{n,y} \sim u\, |\bar V_r|\, \mathcal{L} \sim \frac {|\bar V_r|\, f_1}{\kappa _c^2}. \end{equation}
Figure 18(b) shows
$F_{n,y}$
as a function of
$|\bar V_r| f_1/\kappa _c^2$
for all cases. The data, spanning four orders of magnitude in
$K$
, are in good agreement with the scaling (4.17). In particular, for vibrational cases such as
$K=1$
, both
$\bar V_r$
and
$F_{n,y}$
are zero (see the red symbols in figure 18
b), whereas in the translational mode,
$F_{n,y}$
is approximately proportional to
$|\bar V_r|\, f_1/\kappa _c^2$
(e.g. the blue symbols in figure 18
b).
4.4. Stability analysis
To clarify the mechanism underlying the system bifurcation, stability tests are performed for representative cases. In these tests, unlike the freely moving filaments considered previously, the filament is towed laterally at a prescribed relative velocity
$V_r$
, and the resulting lateral force
$F_y$
is measured. Figure 19(a) shows
$F_y$
as a function of
$V_r$
for
$\delta = 30^\circ$
,
$\textit{Re} = 80$
and
$K = 0.01$
or 1. The fitted curves based on the normal form of a supercritical pitchfork bifurcation (Strogatz Reference Strogatz2018), i.e.
$F_y = C V_r (V_r^2 - V_{r0}^2)$
, are also included, with
$C=-0.337$
and
$V_{r0}=0.55$
for
$K=0.01$
, and
$C=-0.046$
and
$V_{r0}=1.55{\rm i}$
for
$K=1$
(where
${\rm i}$
is the imaginary unit).
(a) Lateral forces
$F_y$
and (b) generalised hydrodynamic potentials
$\varPsi = - \int _{0}^{V_r} F_y (V_r^\prime )\, \mathrm{d} V_r^\prime$
of the towed filaments as functions of the relative lateral velocity
$V_r$
for the cases with
$\delta = 30^\circ$
,
$\textit{Re} = 80$
and
$K = 0.01, 1$
. The fitting curve in (a) is
$F_y = C V_r (V_r^2 - V_{r0}^2)$
, with
$C=-0.337$
,
$V_{r0} = 0.55$
and
$C=-0.046$
,
$V_{r0} = 1.55i$
for cases
$K=0.01$
and
$K=1$
, respectively.
For the large-stiffness cases (i.e.
$K=1$
), the lateral force vanishes when the filament is laterally stationary relative to the surrounding fluid, namely
$F_y=0$
at
$V_r=0$
(see figure 19
a). For
$V_r \neq 0$
, the lateral force always opposes the relative motion, e.g.
$F_y\gt 0$
when
$V_r\lt 0$
. In addition, the magnitude of
$F_y$
increases monotonically with
$|V_r|$
(see figure 19
a). This behaviour is consistent with physical expectation, as it reflects the resistance exerted by the fluid against the filament’s lateral motion. Consequently, the fitted
$F_y$
curve has only one zero, located at the origin (see the dashed line in figure 19
a). The negative local slope at this point indicates that a small positive or negative perturbation in
$V_r$
generates a restoring force of opposite sign. Thus the fluid force acts as an effective spring, stabilising the filament against small perturbations. This restoring mechanism explains why the filament remains stable under different initial perturbations in the S and V modes, as shown in figure 3.
For the cases with
$K=0.01$
, a slight increase in
$V_r$
from zero makes
$F_y$
positive (see figure 19
a). In other words, the lateral fluid force is aligned with the relative motion of the filament, indicating that the fluid promotes, rather than suppresses, the lateral motion. Once
$V_r$
exceeds a critical value
$V_{r0}$
, the direction of
$F_y$
returns to the conventional one, namely opposite to
$V_r$
. Beyond this threshold, the magnitude of
$F_y$
increases rapidly with increasing
$V_r$
. For
$V_r\lt 0$
, the response is exactly antisymmetric to that for
$V_r\gt 0$
(see figure 19
a). Under the circumstances, the fitted
$F_y$
curve has three zeros, i.e.
$V_r=\pm V_{r0}$
and 0 (see the solid line in figure 19
a). However, the positive slope at
$V_r=0$
indicates that small perturbations in
$V_r$
are amplified by a fluid force of the same sign, driving the filament toward a nearby equilibrium state. Hence the equilibrium at
$V_r=0$
is unstable, while the two non-zero roots (
$V_r=\pm V_{r0}$
) are stable fixed points, as their local slopes are negative.
This also provides a reasonable explanation for the tendency of the filament to move upwind when
$\beta =0$
, as shown in figure 4. In this case, the initial flow field is unperturbed, and the initial relative lateral velocity of the filament is
$V_r=-U_{\infty ,y}$
. Because this initial state is closer to the left equilibrium point, i.e.
$V_r=-V_{r0}$
, the filament evolves towards that state, and eventually undergoes upwind motion (see figure 4
d). By contrast, if the initial disturbance is sufficiently large to place the initial relative velocity closer to the right equilibrium point, the filament will evolve towards downwind motion (e.g.
$\beta =-1.5$
in figure 4
d). As a result, the filament response splits into two branches under different initial perturbations, which are symmetric in the frame moving laterally with the fluid (see figures 6
a,b).
Furthermore, the system stability can also be illustrated from an energy perspective by introducing the generalised hydrodynamic potential
$\varPsi$
, defined as
\begin{equation} \varPsi = - \int _{0}^{V_r} F_y (V_r^\prime )\, \mathrm{d} V_r^\prime . \end{equation}
Figure 19(b) presents
$\varPsi$
versus
$V_r$
for the towed filaments at
$\delta = 30^\circ$
,
$\textit{Re} = 80$
and
$K = 0.01, 1$
. Stable wells are observed in the resulting potential landscape. The well bottoms (where
$\varPsi$
reaches a local minimum) correspond to the stable equilibrium states, i.e.
$V_r = 0$
for
$K = 1$
, and
$V_r = \pm V_{r0}$
for
$K = 0.01$
, while the peak points (where
$\varPsi$
is locally maximal) represent unstable equilibrium states that are not observed in the simulations, e.g.
$V_r = 0$
for
$K = 0.01$
.
Figure 20 presents contours of
$\varPsi$
in the
$\textit{Re}$
–
$V_r$
plane for
$\delta = 30^\circ$
and
$K = 0.01$
or 1. The circular symbols represent
$V_r$
of the freely moving filaments (as shown in figure 6), whereas the solid and dashed lines denote the loci of the local minima and maxima of
$\varPsi$
, respectively. It is seen that for both
$K$
and a range of
$\textit{Re}$
, the equilibrium states of the freely moving filaments agree closely with the local minima of
$\varPsi$
(see figure 20). The slight discrepancy for
$K=0.01$
at large
$\textit{Re}$
(e.g.
$\textit{Re} = 200$
in figure 20
a) may result from the prescribed constant velocity of the towed filament, which differs from the motion of a freely moving filament.
For
$K=0.01$
(figure 20
a), the well bottom is located at
$V_r=0$
when
$\textit{Re}$
is small, similar to the case
$K=1$
(figure 20
b). As
$\textit{Re}$
increases beyond a critical value, the well bottom splits into two symmetric branches, corresponding to two equilibrium states with equal
$|V_r|$
but opposite directions. Meanwhile, a ridge appears at
$V_r=0$
, indicating an unstable equilibrium state (see the dashed line in figure 20
a). These features are characteristic of a supercritical pitchfork bifurcation (Strogatz Reference Strogatz2018), while no such bifurcation occurs for
$K=1$
(see figure 20
b).
Contours of the generalised hydrodynamic potentials
$\varPsi$
for (a)
$\delta = 30^\circ$
,
$K = 0.01$
, and (b)
$\delta = 30^\circ$
,
$K = 1$
. The solid and dashed lines denote the loci of the local minima and maxima of
$\varPsi$
, respectively. The circles represent the corresponding
$V_r$
of the lateral unconstrained filaments, as shown in figure 6.
5. Concluding remarks
In this study, the dynamics of laterally unconstrained flexible filaments in an inclined flow is investigated numerically. Due to the non-zero inclination angle of the oncoming flow, the system loses its natural symmetry, and the filament inevitably undergoes lateral translation in the laboratory frame. Remarkably, the filament can spontaneously develop and maintain upwind motion solely through hydrodynamic forcing. This counterintuitive behaviour is a manifestation of a supercritical bifurcation. In a frame moving laterally with the fluid, the upwind and downwind motions are equivalent, with the two branches having identical relative velocity magnitudes
$|\bar V_r|$
. Based on whether relative translation and vibration occur, the filament dynamics can be classified into four modes: static, quasi-static, translational and vibrational. These modes can be identified quantitatively using the relative translational Reynolds number
$\textit{Re}_r$
and the Strouhal number
$St$
, or qualitatively from the corresponding flow structures.
The force decomposition shows that the total normal force components, i.e.
$F_{n,x}$
and
$F_{n,y}$
, contribute mainly to the streamwise drag and lateral thrust, respectively. Corresponding scaling relations are derived as
$F_{n,x} \sim \textit{Re}_{\textit {eff}}^2$
and
$F_{n,y} \sim |\bar V_{r}| f_1 / \kappa _c^2$
. The stability analysis further clarifies the bifurcation behaviour. Before bifurcation, the system possesses a single stable equilibrium at
$V_r = 0$
. After bifurcation, two symmetric stable equilibria appear, and the origin becomes unstable. The generalised hydrodynamic potential
$\varPsi$
is introduced to characterise the system stability, with its local minima representing the stable equilibrium states. These results deepen the understanding of flexible-body dynamics in uniform flow.
Indeed, the present simulations are restricted to two dimensions. Three-dimensional effects, such as spanwise deformation and instabilities, may significantly alter the dynamics of flexible bodies. Moreover, a systematic three-dimensional investigation would introduce additional governing parameters, e.g. the aspect ratio. Exploring these effects therefore represents an important direction for future work.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.12215.
Acknowledgements
Numerical computations were performed at the Hefei Advanced Computing Centre.
Funding
This work was supported by the Natural Science Foundation of China (grant nos 12388101, 12502308, 12293002, 12293000, 12502284), the Postdoctoral Fellowship Programme of CPSF under grant no. GZC20241645, and the China Postdoctoral Science Foundation under grant no. 2025M781800.
Declaration of interests
The authors report no conflict of interest.
L
U∞
δ
U∞,x=U∞cosδ
U∞,y=U∞sinδ
x
y
U∞
U∞
U∞
x
y
s
V
V>0
V<0
δ=30∘
Re=100
K=0.01
M=0.5
S=104
Fx/Fref
Fref=(1/2)ρUref2L
β
δ=20∘
K=0.0001
Re=10
δ=20∘
K=0.01
Re=50
β
β
δ=40∘
K=0.0001
Re=50
δ=30∘
K=0.01
Re=100
β
V¯
Re
δ
K=10−4
K=10−2
K=1
V¯r
Re
δ
K=10−4
K=10−2
K=1
(δ,Re,K)
∘
△
▽
◻
Rer
Rex
δ
K=10−4
K=10−2
K=1
St
Rex
δ
K=10−4
K=10−2
K=1
θ
Rex
δ
K=10−4
K=10−2
K=1
Fx
Fy
δ=20∘
K=0.0001
Re=10
δ=40∘
K=0.0001
Re=50
Fx
Fy
δ=20∘
K=0.01
Re=50
δ=30∘
K=0.01
Re=100
t1−t8
t/Tref=t1
t2
t3
t4
δ=20∘
K=0.01
Re=50
t1−t4
t/Tref=t5
t6
t7
t8
δ=30∘
K=0.01
Re=100
t5−t8
A′
B′
τ
n
α
x
f¯n,x
x
f¯τ,x
y
f¯n,y
y
f¯τ,y
Re=50,80,100,150
δ=30∘
K=0.01
β=0
Fd
x
fn
fτ
Fn,x
Fτ,x
y
fn
Fn,y(t)
Fn,y(d)
y
fτ
Fτ,y
Re=50,80,100,150
δ=30∘
K=0.01
β=0
x
Fn,x
Reeff
y
Fn,y
|V¯r|f1/κc2
Fy
Ψ=−∫0VrFy(Vr′)dVr′
Vr
δ=30∘
Re=80
K=0.01,1
Fy=CVr(Vr2−Vr02)
C=−0.337
Vr0=0.55
C=−0.046
Vr0=1.55i
K=0.01
K=1
Ψ
δ=30∘
K=0.01
δ=30∘
K=1
Ψ
Vr
