1. Introduction
The study of the dynamics of elastic fibres in flow is of interest in diverse industrial applications as well as fundamental biological processes. At the microscale, experimental and numerical studies have focused on the shape deformations of slender fibres in linear shear flows, compressive flows and simple two-dimensional vortical flows (e.g. Young & Shelley Reference Young and Shelley2007; Du Roure et al. Reference Du Roure, Lindner, Nazockdast and Shelley2019; Chakrabarti et al. Reference Chakrabarti, Liu, LaGrone, Cortez, Fauci, du Roure, Saintillan and Lindner2020; Słowicka et al. Reference Słowicka, Xue, Sznajder, Nunes, Stone and Ekiel-Jeżewska2022). Bending and buckling similar to an Euler–Bernoulli structural element were observed and the tug-of-war between the viscous and bending forces experienced by a fibre has been described by an elastoviscous number that characterises the morphological transitions (Tornberg & Shelley Reference Tornberg and Shelley2004).
A separate set of experimental and numerical works have explored deformation, motion, alignment and dispersion of flexible passive fibres in turbulent flows in the past two decades (Marchioli, Rosti & Verhille Reference Marchioli, Rosti and Verhille2025). An intriguing example of flexible fibres in flow is diatom chains in the ocean (Young et al. Reference Young, Karp-Boss, Jumars and Landis2012). The classical assumption is that these phytoplankton experience turbulence as simple linear shear (Lazier & Mann Reference Lazier and Mann1989). This is based on the observation that phytoplankton cells are generally smaller than the diameters of the smallest coherent vortices of dissipating turbulence, and that viscosity quickly induces an approximately linear velocity gradient, leading to constant shear and vorticity across a scale of about 1 mm. However, progress in the understanding of turbulence as ‘a writhing tangle of vortex worms’ with vorticity diffusion, stretching and variable spatial shear challenges this classical assumption (Davidson Reference Davidson2004). Jumars et al. (Reference Jumars, Trowbridge, Boss and Karp-Boss2009) proposed the Burgers vortex as a simplified model of vortex structure and dynamics that closely resembles many dissipative vortices observed in the laboratory, in the field and in direct numerical simulations of turbulence. Recent experiments with copepods placed in a laboratory-generated Burgers vortex have demonstrated species-dependent behavioural changes (Elmi, Webster & Fields Reference Elmi, Webster and Fields2022). While interest is growing in exploring microswimmer interaction and trapping in vortical flows across a range of Reynolds numbers (Webster & Young Reference Webster and Young2015; Tanasijević & Lauga Reference Tanasijević and Lauga2022), questions remain as to understanding how passive fibres behave in a vortical setting with stretching.
In this paper, we focus on two primary objectives. First, we develop a Stokes flow analogue of the Burgers vortex using the superposition of regularised singular solutions of the Stokes equations. The resulting background flow, which we call a spiralet, matches that of the inertial Burgers vortex at the horizontal plane of the vortex core. More importantly, it provides us with spatial variation in vorticity and stretching, which is needed for the next primary objective: to study the morphological dynamics of passive flexible filaments in this flow, and to capture robust orbits of the fibres as they spin and deform.
2. Mathematical model
2.1. Burgers-like vortex
We begin by setting up a background flow that is an analogue of the Burgers vortex at the zero-Reynolds-number limit. The Burgers vortex is an analytical solution of the Navier–Stokes equation (Burgers Reference Burgers, von Mises and von Kármán1948) with steady radial, azimuthal and axial velocities:
\begin{equation} u_{\tilde {r}}=-\gamma \tilde {r}, \quad u^{}_{\theta }=\frac {\varGamma }{2\pi \tilde {r}}\left (1-{\textrm e}^{-\gamma \tilde {r}^2/2\nu }\right ), \quad u_z=2\gamma z \end{equation}
in a cylindrical coordinate system with
$\tilde {r}$
as the radial distance and kinematic viscosity
$\nu$
. The vortex parameters
$\gamma$
and
$\varGamma$
are the axial strain rate and the vortex circulation, respectively. Two properties of the Burgers vortex are:
-
(i) On planes where
$z$
is constant, the flow spirals inward towards the
$z$
axis. -
(ii) Along the
$z$
direction, the flow stretches outward from the origin, with flow magnitude proportional to
$|z|$
.
We seek to design a flow with similar local properties, but one that decays at infinity. In the context of zero Reynolds number, the fluid motion is described by the Stokes equations:
\begin{align} \mu \Delta \boldsymbol{u} & = \boldsymbol{\nabla }p - \boldsymbol{\mathcal{F}}, \end{align}
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} & = 0. \end{align}
Here
$\mu$
is the fluid viscosity,
$\boldsymbol{u}$
is the fluid velocity,
$p$
is the pressure and
$\boldsymbol{\mathcal{F}}$
represents the force (per unit volume) exerted on the fluid by an external mechanism. We use two types of forcing
$\boldsymbol{\mathcal{F}}$
: one type to generate the background flow and another to introduce the effect of the fibre.
The Stokes equations have a well-known fundamental solution (the Stokeslet) and other exact solutions derived from it. The background flow is generated by the sum of a rotational flow about the
$z$
axis (a rotlet) and an extensional flow (a stresslet) aligned with the
$z$
axis. Since these elements are placed in the fluid domain, we use regularised forcing as in the method of regularised Stokeslets (Cortez, Fauci & Medovikov Reference Cortez, Fauci and Medovikov2005), although here the regularisation parameter
$\epsilon$
need not be small. Instead,
$\epsilon$
represents the radius of the vortex stirring the fluid. Specifically, the force is
$\boldsymbol{\mathcal{F}}=\boldsymbol{\nabla }\times \rho \hat {\boldsymbol{e}}_3\ \phi _\epsilon (r)-\sigma \hat {\boldsymbol{e}}_3 (\hat {\boldsymbol{e}}_3\boldsymbol{\cdot }\boldsymbol{\nabla }) \phi _\epsilon (r)$
for the rotlet and stresslet, respectively, where
$\rho$
and
$\sigma$
are positive numbers,
$r=|\boldsymbol{x}-\boldsymbol{x_0}|$
and
$\epsilon$
is the regularisation parameter or blob size, which we also refer to as the vortex radius. Here
$\boldsymbol{x}$
is an evaluation point anywhere on the fluid domain,
$\boldsymbol{x_0}$
is the location of the regularised singularity and
\begin{equation} \phi _\epsilon (r)=\frac {15\epsilon ^4}{8\pi (r^2+\epsilon ^2)^{7/2}}. \end{equation}
The resulting background velocity is
\begin{equation} \boldsymbol{u}^{b}({\boldsymbol x})= \rho \left (\frac {2R^2+3\epsilon ^2}{8\pi \mu R^5} {\hat {\boldsymbol{e}}}^{}_3\times \boldsymbol{\tilde x} \right ) + \sigma \left (\frac {3\tilde {x}_3^2 -R^2}{8\pi \mu R^5} \boldsymbol{\tilde x} + \frac {3\epsilon ^2 \tilde {x}_3}{8\pi \mu R^5} \hat {\boldsymbol{e}}^{}_3 \right ), \end{equation}
where we have used the shorthand notation
$\boldsymbol{\tilde x} = \boldsymbol{x}-\boldsymbol{x}^{}_{\boldsymbol{0}}$
and
$R^2 = r^2+\epsilon ^2$
. In the current study, we choose
$\boldsymbol{x}_{\boldsymbol{0}} = (0,0,0)$
. Expressions in the cylindrical coordinates for this background flow are given in Appendix A.1.
Since the stresslet and the rotlet are solutions to the Stokes equations, their sum, also a solution, creates a simple analogue of the Burgers vortex that we name a spiralet (figure 1
a). It turns out that, for both the Burgers vortex and the spiralet at the symmetry plane, only the azimuthal component of velocity contributes to the vorticity, with no contribution from the stresslet or from the axial strain. We briefly outline how the spiralet parameters can be tuned so that the vorticity and azimuthal velocity profiles closely approximate those of a Burgers vortex near the core. For a given Burgers vortex (equation (2.1)), we can choose
$\rho$
and
$\epsilon$
so that the maximum azimuthal speeds match in value and in location.
Figure 1. (a) Streamlines in a spiralet vortex comprised of rotation (from the rotlet) and radial inflow and axial stretching (due to the stresslet). (b) Comparison of azimuthal speed and vorticity on the
$xy$
plane of a Burgers vortex and spiralet. Burgers vortex parameters used are
$\gamma ={3.08\times 10^5}\,\textrm{s}^{-1}$
,
$\varGamma ={9.55\times 10^3}\,\mathrm{\unicode{x03BC} m^2\,s^{-1}}$
. The corresponding spiralet parameters are set to
$\epsilon ={5}\,\mathrm{\unicode{x03BC} m}$
,
$\rho = \sigma ={0.15}\,\mathrm{g\,\unicode{x03BC} m^2\,s^{-2}}$
to match the location and magnitude of the maximum azimuthal speeds of the two vortices. Here viscosity
$ \mu = 10^{-6} \,\mathrm{g\,\unicode{x03BC} m^{-1}\,s^{-1}}$
and kinematic viscosity
$\nu = 10^{6}\, \mathrm{\unicode{x03BC} m^2\, s^{-1}}$
. (c) Spatial variation of the axial component of vorticity in the spiralet vortex for the chosen set of parameters.
We first match the radial distance at which the azimuthal speeds,
$u_\theta$
, of the vortex and the rotlet reach their peak. The rotlet flow is given by the first term on the right-hand side of (2.4). It has a maximum magnitude proportional to
$\rho /\mu \epsilon ^2$
at
\begin{equation} r_m =\frac {\epsilon }{2}\sqrt {-7+\sqrt {69}} \approx 0.5715\epsilon \end{equation}
on the
$xy$
plane. On the other hand, the azimuthal speed of the Burgers vortex in (2.1) attains its maximum value when
$r_m\approx 1.5852\sqrt {\nu /\gamma }$
. Thus, the Burgers vortex and the spiralet attain the maximum azimuthal speed at the same distance from the axis when we choose
$ \epsilon = 2.77357\sqrt {\nu /\gamma }$
. The maximum azimuthal speeds are matched by evaluating
$u_\theta$
at
$r_m$
. This gives the choice
$\rho = 8.7118 \varGamma \mu \sqrt {\nu /\gamma }$
.
On the plane of symmetry, the vorticity is of the form
$(0,0,\omega )$
, aligned with the vortex axis:
\begin{equation} \omega = \frac {1}{\mu }\frac {10\epsilon ^4-7\epsilon ^2\tilde {r}^2-2\tilde {r}^4}{8\pi (\tilde {r}^2+\epsilon ^2)^{7/2}}\rho ,\quad \omega _{max }=\frac {10 \rho }{8\pi \mu \epsilon ^3}, \end{equation}
where
$\tilde r = \sqrt {x^2+y^2}$
is the distance from the
$z$
axis to the evaluation point, and the peak occurs at the vortex centre. The values of
$\epsilon$
and
$\rho$
completely determine the rotlet vorticity, whose maximum value differs by about
$2\,\%$
from the Burgers maximum vorticity.
Figure 1(b) shows the azimuthal speed (blue solid curve) and the vorticity (red solid curve) as a function of radial distance from the centre for representative parameters. They are in good agreement with those of the Burgers vortex (dashed curves) located at the same plane. Away from the core, in the direction of the vortex axis, the spiralet velocity decays to zero (equation (2.4)), while in the Burgers vortex, the axial component of velocity increases linearly.
The axial component of vorticity is maximum at the centre of the spiralet. When moving either radially or axially away from the core, as shown in figure 1(c), this vorticity decays rapidly. In fact, within one vortex radius
$\epsilon$
, the vorticity falls to less than 1 % of the maximum. These large gradients in vorticity will have direct effects on shape deformations of flexible filaments that are long in comparison with the vortex core diameter, as we discuss below.
2.2. Kirchhoff rod model with regularised Stokeslet segments
We next describe the slender, passive fibres that will interact with the vortex. Using an unconstrained Kirchhoff rod framework based on Olson, Lim & Cortez (Reference Olson, Lim and Cortez2013), our fibre of length
$L$
is represented by a space curve
$\boldsymbol{X}(s,t)$
with a set of orthonormal triads,
$\boldsymbol{D}^1(s,t), \boldsymbol{D}^2(s,t), \boldsymbol{D}^3(s,t)$
. These are initially defined as unit normal, binormal and tangent vectors to the rod, respectively. Here,
$t$
is time and
$s$
is the arclength parameter (
$0\leqslant s \leqslant L$
). As the shape of the fibre responds to the background flow, it develops forces and torques along its length. The fibre is initiated as a straight rod with its centroid at the centre of the vortex with a slight spatial perturbation (tilted
$1^\circ$
relative to the horizontal plane).
If
$\boldsymbol{F}(s,t)$
and
$\boldsymbol{N}(s,t)$
are the internal forces and moments transmitted across a section of the fibre, then the momentum and angular momentum balances let us relate them with applied force density,
$\boldsymbol{f}(s,t)$
, and torque density,
$\boldsymbol{n}(s,t)$
, as
\begin{equation} 0 = \boldsymbol{f}+\frac {\partial \boldsymbol{F}}{\partial s} \quad \text{and}\quad 0 =\boldsymbol{n}+\frac {\partial \boldsymbol{N}}{\partial s}+\left (\frac {\partial \boldsymbol{X}}{\partial s}\times \boldsymbol{F}\right ), \end{equation}
where the constitutive relations between force, torque and triads can be derived from an elastic energy (Olson et al. Reference Olson, Lim and Cortez2013), and are given in component form by
\begin{equation} N_i=a_i\left (\frac {\partial \boldsymbol{D}^j}{\partial s}\boldsymbol{\cdot }\boldsymbol{D}^k-\varOmega _i\right ), \quad F_i= b_i\left (\boldsymbol{D}^i \boldsymbol{\cdot }\frac {\partial \boldsymbol{X}}{\partial s}-\delta _{3i}\right ), \end{equation}
where
$a_1$
and
$a_2$
are bending constants and
$a_3$
is the twist modulus. The parameters
$b_1$
and
$ b_2$
are shear moduli and
$b_3$
is a stretch modulus. Here we take
$a_1=a_2=a_3=EI$
and
$b_1=b_2=b_3=b$
, where
$EI$
is the bending rigidity and
$b$
is chosen large enough to keep the rod inextensible (see table 1). In (2.8),
$(i,j,k)$
is any cyclic permutation of indices
$(1,2,3)$
and
$\delta _{3i}$
is the Kronecker delta. The intrinsic shape of the elastic fibre is specified by the strain-twist vector
$(\varOmega _1(s,t), \varOmega _2(s,t), \varOmega _3(s,t))=(0,0,0)$
, so that it is natively straight. The resultant local curvature and torsion can be expressed in terms of the directors as
\begin{equation} \kappa \left (s,t\right ) = \sqrt {\left (\frac {\partial \boldsymbol{D}^2}{\partial s}\boldsymbol{\cdot }\boldsymbol{D}^3\right )^2+\left (\frac {\partial \boldsymbol{D}^3}{\partial s}\boldsymbol{\cdot }\boldsymbol{D}^1\right )^2}, \quad \xi \left (s,t\right )=\left |\frac {\partial \boldsymbol{D}^1}{\partial s}\boldsymbol{\cdot }\boldsymbol{D}^2\right | .\end{equation}
The force density can be interpreted as the transmission of forces and torques generated due to fibre shape deformation (per unit length, integrated over the whole length). The velocity they produce is the solution of the Stokes equations (2.2) with
\begin{equation} \boldsymbol{\mathcal{F}}(\boldsymbol{x},t)= \int _{0}^{L} -\boldsymbol{f}(s,t)\phi _{\delta }(\boldsymbol{x}-\boldsymbol{X}(s,t))\textrm{d} s+\frac {1}{2}\boldsymbol{\nabla }\times \int _{0}^{L}-\boldsymbol{n}(s,t)\phi _{\delta }(\boldsymbol{x}-\boldsymbol{X}(s,t))\textrm{d} s, \end{equation}
where we use the blob function defined in (2.3) with the regularisation parameter
$\delta$
to spatially distribute the forces and torques. The fluid motion generated by this
$\boldsymbol{\mathcal{F}}$
is
\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \int _0^L \boldsymbol{u}^0(\boldsymbol{x};\boldsymbol{X}(s,t),-\boldsymbol{f}(s,t))\textrm{d} s + \frac {1}{2}\boldsymbol{\nabla }\times \int _0^L \boldsymbol{u}^0(\boldsymbol{x};\boldsymbol{X}(s,t),-\boldsymbol{n}(s,t))\textrm{d} s, \end{equation}
where the integrands are regularised Stokeslets. For instance,
\begin{equation} \boldsymbol{u}^0(\boldsymbol{x};\boldsymbol{X}(s,t),\boldsymbol{f}(s,t)) = \frac {\delta ^2+R^2}{8\pi \mu R^3} {\boldsymbol f} + \frac {({\boldsymbol f}\boldsymbol{\cdot }\boldsymbol{\tilde x}) \boldsymbol{\tilde x}}{8\pi \mu R^3}. \end{equation}
The integrals are approximated by discretising the filament centreline into
$N$
linear segments of length
$\Delta s$
. The force and torque densities,
$\boldsymbol{f}$
and
$\boldsymbol{n}$
, are assumed to be linear in each segment. Under these assumptions, the integrals in (2.11) can be evaluated exactly at an arbitrary point
$\boldsymbol{x}$
using the method of regularised Stokeslet segments (Cortez Reference Cortez2018). The orthonormal triad update requires the integration of a rotlet and a dipole. By the linearity of the Stokes equation, the velocity at a point
$\boldsymbol{x}$
due to the
$N$
segments can be evaluated as a sum of the contributions due to each segment.
Table 1. Dimensional parameter values.
In the cases that we investigate, a no-slip condition is imposed on the immersed fibre so that it moves with the local fluid velocity (equation (2.13)) due to the steady background spiralet
$ (\boldsymbol{u}^b )$
as well as fluid flow created by the body forces and torques exerted by the fibre itself
$ (\boldsymbol{u} )$
. The resulting local angular velocity from these two contributions is used to rotate the associated triads (equation (2.14)):
\begin{align} \frac {\partial \boldsymbol{X}(s,t)}{\partial t} &=\boldsymbol{u}\left (\boldsymbol{X}(s,t),t\right )+\boldsymbol{u}^b(\boldsymbol{X}(s,t),t), \end{align}
\begin{align} \frac {\partial \boldsymbol{D}^i(s,t)}{\partial t} &= \left (\frac {1}{2} \boldsymbol{\nabla}\times \left ( \boldsymbol{u}(\boldsymbol{x},t)+\boldsymbol{u}^b(\boldsymbol{x},t) \right )\Big |_{\boldsymbol{x}=\boldsymbol{X}(s,t)}\right )\times \boldsymbol{D}^i(s,t) . \end{align}
3. Results
To characterise the relative importance of viscous and elastic forces, we define the elastoviscous number as
\begin{equation} \eta = \frac {8 \pi \mu \omega _{\textit{max}}L^4}{cEI}. \end{equation}
Here
$\omega _{\textit{max}}$
is the maximum vorticity at the vortex centre (defined in (2.6)),
$EI$
is the fibre’s bending rigidity,
$L$
is fibre length and
$c = -\ln (e\delta ^2/L^2 )$
, where we interpret the regularisation parameter
$\delta$
as the radius of the fibre and
$e$
is Euler’s number. The fibre is natively straight and torsion-free, but the background spiralet flow will induce shape deformations, and, hence, forces and torques that try to restore the intrinsic shape. First, we study the special case where the fibre is placed in the horizontal plane of symmetry, with its centroid at the vortex centre. The interplay between elastic forces, viscous forces and their respective time scales results in a succession of fibre shapes as it completes its excursion.
3.1. Symmetrically placed fibres
Figure 2 shows the orbits and shape evolution of three fibres (
$L = {20}\,\mathrm{\unicode{x03BC} m}$
) with increasing elastoviscous numbers in the same background spiralet flow as in figure 1. The fibres are initialised with their centres of mass at the vortex centre, slightly tilted out of the horizontal plane so that one-half of the fibre is above the
$xy$
plane, while the other half lies below. As the whole filament rotates, each half moves in opposite directions towards the vortex axis, causing it to move out of the
$xy$
plane, before it eventually reaches vertical alignment. Note that since the fibre is initialised in a straight configuration with its centre at the vortex core, the flow symmetry dictates that the centre remains at the vortex core, resulting in no fibre transport. For the three cases shown, the variation in the elastoviscous number is achieved by lowering the bending rigidity from
$EI = {6.99\times 10^{-22}}\,{\textrm{J m}}$
for figure 2(a) to
$EI = {4.96\times 10^{-25}}\,{\textrm{J m}}$
for figure 2(c). For the smallest elastoviscous number, this stiff fibre rotates much like a solid body in this vortical flow. Decreasing bending rigidity renders the fibre more flexible, and we see ‘S’ shapes in both horizontal and vertical orientations for the intermediate elastoviscous number. Increasing the elastoviscous number further results in more pronounced deformations, turning the nearly planar ‘S’ shapes into fully three-dimensional shapes that, nevertheless, retain symmetry due to their initial position. Because of the Kirchhoff rod formulation, we are able to track the torsion that develops along the fibres as they move in the spiralet (equation (2.9)). The fibres are coloured with their evolving torsion in figure 2. We note that the softest fibre in figure 2(c) shows pronounced variation in lengthwise torsion.
Figure 2. Flexible fibres initially centred at the vortex centre with increasing elastoviscous number. (a) Rigid rotation is seen in the first column for
$\eta ={3.54\times 10^2}$
with eventual vertical alignment due to stretching. (b,c) Fibres in the second and third columns (
$\eta ={4.99\times 10^4}$
and
$\eta ={4.99\times 10^5}$
, respectively) show more and more pronounced deformations. The colours on the fibre surface indicate local torsion
$(\xi )$
induced by the background flow. All the pertinent flow and fibre parameters are listed in table 1. The full dynamics is shown in supplementary movie 1 (case 2a), movie 2 (case 2b) and movie 3 (case 2c) available at https://doi.org/10.1017/jfm.2026.11342.
The kymographs in figures 3(a) and 3(b) track the curvature evolution of the two most flexible fibres in figure 2. In both cases, positions of maximum curvature first appear near the fibre’s centroid, where it is closest to the vortex centre. The high-curvature regions propagate outwards along the fibre towards opposite ends, until the fibres regain their straight posture. Comparing figures 3(a) and 3(b), we see that for the most flexible fibre, the maximal curvature regions are more localised. Eventually, the fibres straighten out and align with the vertical axis. Although the background flow reorients fibres in all three cases in figure 2 from horizontal to vertical, the time taken in doing so varies significantly between these cases – an observation that we discuss below. We chose to characterise the deformation of the fibres by a ‘tortuosity’ measure, or total curvature, along the fibre,
$\displaystyle \tau (t)= \int _0^L \kappa (s,t)\textrm{d} s$
, where
$\kappa$
is the unsigned local curvature. This quantitative measure of the shape evolution for the three fibres in figure 2 is shown in figure 3(c). Since all cases begin as a horizontal straight fibre and end as a vertical straight fibre, somewhere in between are the curved, deformed states. The two flexible cases (higher
$\eta$
) appear as skewed Gaussian shapes, with the most flexible case reaching the highest total curvature. In both cases, there is a steep rise in total curvature from zero to its maximum, with a more gradual decrease back to zero as the fibre reaches its straight, vertical equilibrium. For the stiffest case, the shape remains straight throughout. Hence, total curvature remains negligibly small compared with the other two cases.
Figure 3. Time–curvature plots showing travelling waves of maximal curvature positions for (a)
$\eta ={4.99\times 10^4}$
and (b)
$\eta ={4.99\times 10^5}$
(the two cases of figure 2
b,c). Curvature values are higher and peak curvature regions are more localised for this larger value of
$\eta$
. (c) Total curvature of the achieved shapes for the three cases in figure 2 as a function of time. Zero total curvature at the left indicates the straight horizontal configuration, while zero at the right is for the straight vertical configuration.
For fibres in linear shear flows, it has been shown that the elastoviscous number serves as a bifurcation parameter that dictates the observed shape modes (e.g. Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018). For instance, as
$\eta$
increases, fibres transition from rigid rotations to S-buckling to snaking. Based on our definition in (3.1), although
$\eta$
captures the maximum vorticity of the spiralet, the spatial variation in vorticity is not captured. In fact, we find that the ratio of the fibre length (
$L$
) and the vortex core diameter (
$d = 2\epsilon$
) also has a significant effect on fibre shape dynamics. Figure 4 shows snapshots from two simulations where an identical fibre was placed in two spiralets of different vortex core diameters. We see much more pronounced deformations in the vortex with a smaller core on the left. On the symmetry plane, the vorticity decays to less than 1 % of its maximum at a radial distance of
$\epsilon$
from the centre (figure 1
b). Therefore, the vorticity gradients felt by the fibre are stronger in spiralets with tighter diameters. For instance, when
$d/L\geqslant 1$
, the whole length of the fibre is placed inside the vortex core, but for
$d/L\lt 1$
only the midsection of the fibre is contained inside the vortex core. This midsection, then, is exposed to much greater vorticity than the two ends.
Figure 4. Comparison of shape deformations for two identical fibres put in spiralets of different diameters
$d=2\epsilon$
, while keeping
$\eta$
fixed: (a)
$d/L = 0.5,\ \eta = {1.18\times 10^6}$
and (b)
$d/L = 1.5,\ \eta = {1.18\times 10^6}$
. The fibre put in the smaller vortex is exposed to a greater vorticity gradient and deforms more.
We now study the maximal total curvature,
$\tau _{\textit{max}} = \max _ {t} \tau (t)$
, achieved during the fibre’s orbit as a function of
$\eta$
and
$d/L$
. Figure 5 shows the phase plot of
$\tau ^{}_{\textit{max}}$
using data from 60 different realisations of the fibre–fluid system for
$d/L={0.25}\textrm { to }{1.5}$
and
$\eta ={7.16\times 10^2}\textrm { to }{2.36\times 10^6}$
. Here, the most rigid cases that remain straight throughout appear at the left end of the phase plot (fibre length
${15}\,\mathrm{\unicode{x03BC} m}$
with
$EI={10^{-22}}\,{\textrm{J m}}$
). For all other cases, the fibres had a length of
${20}\,\mathrm{\unicode{x03BC} m}$
with
$EI={10^{-24}}\,{\textrm{J}\,\textrm{m}}$
, and the phase space was explored by varying the strength of the vortex
$\omega ^{}_{\textit{max}}$
, as well as the vortex diameter
$d$
. We remark that, as a check, we also performed simulations where the elastoviscous number was changed by varying fibre lengths
$L$
rather than
$\omega ^{}_{\textit{max}}$
, and the achieved
$\tau ^{}_{\textit{max}}$
values were nearly identical to those found at the same
$(\eta , d/L)$
within the 60 data points. Also shown in blue in figure 5 are some representative shapes corresponding to the particular values of
$\tau ^{}_{\textit{max}}$
indicated. Moving vertically in the phase plot corresponds to holding
$\eta$
fixed while increasing the ratio
$ d/L$
. For a fibre of fixed length
$L$
and bending rigidity
$EI$
, moving vertically means increasing
$\epsilon$
, while at the same time varying the rotlet strength
$\rho$
(equation (2.6)) to keep
$\omega _ {max}$
(and
$\eta$
) fixed. We note that shape deformation behaviour is very different as we move vertically in this phase plot for different
$\eta$
values. For low elastoviscous numbers, the vortex diameter does not have much effect on total curvature, as is evident from near-vertical contour lines in the darker region. However, for higher values of
$\eta$
, there is a greater variation in total curvature as
$d/L$
varies. For the smallest values of
$d/L = 0.25$
, where only one-quarter of the fibre initially overlapped with the vortex core, we see considerable deformations for values of
$\eta$
even as small as
$1\times 10^5$
. This can be attributed to the higher difference in viscous stresses that a fibre experiences along its length in a smaller vortex.
Figure 5. Phase plot for maximum total curvature
$(\tau _{\textit{max}})$
as a function of
$\eta$
and
$d/L$
, with isocurves of constant total curvature. Also shown are representative fibre shapes in different regions. Results from 60 simulations are shown in the range of
$\eta = {7.16\times 10^2}\textrm { to }{2.36\times 10^6}$
and
$d/L={.25}\textrm { to }{1.5}$
. Representative shapes are shown in different regions; the corresponding rotated three-dimensional shapes are shown in supplementary movie 4.
Variation in vorticity inherent to the background vortex also leads to the development of spatially and temporally varying torsion along the fibre (see the colourmap on fibres in figure 2). For the same set of simulations examined in the curvature phase plot in figure 5, we present a phase plot of maximal total torsion in figure 6. The total torsion along the fibre at time
$t$
is defined as
$\varXi (t)=\displaystyle \int _0^L \xi (s,t)\textrm{d} s$
, where
$\xi (s,t)$
is given in (2.9). We identify the maximum of this value over the duration of the fibre’s orbit from its initial horizontal position to its vertical equilibrium position, symmetric about the
$z=0$
plane as:
$\varXi _{\textit{max}} = \max _ {t} \varXi (t)$
. Similar to the curvature cases, we find the lowest levels of maximal torsion in the stiffest cases (low
$\eta$
). For each fixed
$d/L$
, we see a monotonic increase in maximal torsion as
$\eta$
increases. However, for
$\eta \geqslant {1\times 10^6}$
, total torsion shows a dependence on the relative size of the vortex, with the largest levels of total torsion centred around fibres with
$d/L = 0.75$
.
Figure 6. Phase plot for maximum total torsion
$(\varXi _{\textit{max}})$
as a function of
$\eta$
and
$d/L$
, with isocurves of constant total torsion. Data from 60 simulations are used in the range of
$\eta ={7.16\times 10^2}\textrm { to }{2.36\times 10^6}$
and
$d/L={0.25}\textrm { to }{1.5}$
. Also shown are four representative fibre shapes: A (
$d/L = 1, \eta = {2.95\times 10^5}$
), B (
$d/L = 1, \eta = {1.18\times 10^6}$
), C (
$d/L = 1, \eta = {2.06\times 10^6}$
) and D (
$d/L = 0.5, \eta = {2.06\times 10^6}$
) with colourmaps on the fibre surface highlighting qualitative intensity of local torsion
$\xi (s,t^*)$
(darker indicates higher values) at the time
$t^*$
when
$\varXi$
is maximum.
Also shown in figure 6 are four representative shapes of fibres (denoted by A, B, C and D) at the instant where they achieve the maximal total torsion. While the maximal total torsion occurs when each fibre is near its vertical, equilibrium position, in two of the cases, the fibre is already straight (A and D). However, in the other two cases (B and C), the fibres have not yet straightened out (they eventually will), and the ends are hook-shaped. Interestingly, the shape at the time when
$\varXi _{\textit{max}}$
is attained does not appear to be a function of either
$\eta$
or
$d/L$
. For instance, while
$d/L$
is fixed for cases A, B and C, case A attains
$\varXi _{\textit{max}}$
in the straight conformation, whereas the others are hook-shaped. Conversely, both cases C and D have different
$d/L$
but the same
$\eta$
; yet one is hook-shaped while the other is straight.
In all simulations of this special case, where fibres start slightly perturbed off the horizontal plane (tilted by
$1^\circ$
), centred at the vortex core, they rotate and deform, and eventually reach equilibrium, regaining their straight shape along the vertical axis. Remarkably, we found that the time needed to complete this orbit varied significantly with fibre flexibility. We examine the orbits of four different fibres of the same length, with different bending rigidities, placed within the identical spiralet background flow, with
$\eta$
ranging over four orders of magnitude. Snapshots of the fibres, superimposed, are shown in figure 7(a) at four different times (here we choose the non-dimensional time
$\tilde t = t\,{\omega _{\textit{max}}}$
). At time
$\tilde {t}=0$
, the fibres have the same initial position. At
$\tilde {t}=10$
, the stiffest fibre remains straight while the others, still positioned mostly in the horizontal plane, have realised the signature S-shape. At time
$\tilde {t} = 75$
, the three most flexible fibres have left the horizontal plane, and at time
$\tilde {t}=110$
, they are almost vertical, while the stiffest fibre is just beginning to emerge from the plane. We quantify the temporal evolution of the fibre’s angle with the horizontal plane in figure 7(b), and designate the alignment time
$\tilde {t}_{\textit{st}}$
as the time where the principal axis of the fibre is within
$1^{\circ}$
of the vortex axis. The orientation of the principal axis with respect to the symmetry plane,
$\theta$
, is computed from a gyration tensor-based measure (Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018) (see figure 11 in Appendix A.2). The alignment time is indicated for each of the four simulations in figure 7(b). Figure 7(c) shows the evolution of the total curvature
$\tau$
for the four different cases. As expected, we see distinctly higher total curvature levels for the most flexible (lower
$\eta$
) fibres. At each of the alignment times indicated in figure 7(b), the total curvature of the corresponding fibre has dropped to near zero as it regains its straight shape (see figure 7
c).
Figure 7. (a) Snapshots of four fibres of the same length
$L$
, but different bending rigidities
$EI$
. The fibres are placed within the same spiralet flow (
$\epsilon ={5}\,\mathrm{\unicode{x03BC} m},\ \rho = \sigma ={0.71}\,\mathrm{g}\,\mathrm{\unicode{x03BC} m^2\, s^{-2}},\ \omega _{\textit{max}}={2.28\times 10^3}\,\textrm{s}^{-1}$
) with
$d/L=0.5$
. (b) Time evolution of the fibre orientation
$\theta$
in each of the four simulations. Times needed for the fibres to evolve from a straight horizontal to a straight vertical orientation (alignment times)
$\tilde {t}_{\textit{st}}$
are indicated. (c) Time evolution of the total curvature of each of the four fibres as they reorient from the horizontal to vertical position inside the vortex. (d) Scaled alignment time
$(T=\log _{10}{\tilde {t}_{\textit{st}}})$
as a function of scaled elastoviscous number
$(\mathcal{X} = \log _{10}\eta )$
for three spiralets (
$\epsilon ={5},\ {7.5},\ {10}$
μm and
$\rho = \sigma ={0.71},\ {2.41},\ {5.72}\,\mathrm{g\,\unicode{x03BC} m^2\, s^{-2}}$
, respectively) with different vortex-diameter-to-fibre-length ratios. The inset shows a data collapse when the vortex-diameter-to-fibre-length ratio
$(d/L)$
is included in a shifted vertical axis as
$f(t_{\textit{st}}, d/L)=(\log _{10}\tilde {t}_{\textit{st}} - 1.25)(d/L)$
, where the purple dashed curve is based on (3.2) with
$\mathcal{A} = -0.014$
,
$\mathcal{B} = 0.593$
,
$\mathcal{C} = 0.095$
,
$\mathcal{D} = -1.6$
and
$\mathcal{E} = 6.5142$
.
Choosing the non-dimensional time,
$\tilde t = t \omega _{\textit{max}}$
, allows us to compare the alignment times of fibres between different spiralet flows. As such, figure 7(d) shows a log–log plot of alignment time
$\tilde {t}_{\textit{st}}$
versus
$ \eta$
for fibres in three spiralets with different vortex diameters. In each background flow, the ratio of vortex diameter to fibre length was held fixed (
$d/L = 0.5,\ 0.75,\ 1.0$
), while the bending rigidity was varied. We observe that the greatest variation in alignment time for different
$\eta$
occurs when only half of the fibre is within the vortex core (i.e.
$d/L=0.5$
). For the case
$d/L=1$
, where the whole fibre was initially within the core, the fibre experienced a lesser vorticity gradient. Here, the variation in alignment time is small, but still decreases with
$\eta$
. We also note that the alignment times decrease as
$d/L$
increases.
In each of the three
$d/L$
cases, 19 simulations were performed. The data in figure 7(d) indicate three regions in the functional relationship between the scaled time
$T = \log _{10}{\tilde {t}_{\textit{st}}}$
and the elastoviscous number
$\mathcal{X} = \log _{10}\eta$
. For the stiffest fibres (leftmost region) and the softest fibres (rightmost region), the differences in
$T$
are small. However, there is a transition region where alignment times vary significantly. As such, we assume that the relationship can be approximated by
\begin{equation} {T}\left (\mathcal{X}\right ) = \mathcal{A} \mathcal{X} + \mathcal{B} + \mathcal{C} \tanh \left ({\mathcal{D} \mathcal{X} + \mathcal{E}} \right ). \end{equation}
The values of the parameters in (3.2) were computed separately for each of the three
$d/L$
cases using a least squares fit. Those curves, also included in figure 7(d), show a satisfactory fit to the data. Furthermore, after shifting the data on the vertical axis to include
$d/L$
dependency, we collapse the data into a single curve (see inset in figure 7
d).
To understand why softer fibres take less time to reach equilibrium compared with stiffer ones, we examine the flow features that cause vertical alignment. On the symmetry plane of the spiralet, the background velocity in the vertical direction is zero. However, slightly above the plane, the vertical velocity points upwards, while slightly below the plane, the vertical velocity points downwards. This vertical (axial) velocity
$u_z^b$
is greatest near the vortex centre (see Appendix A.1). Initially, these fibres are placed approximately horizontally in the plane of symmetry, with their centres of mass at the origin. This position is perturbed so that one half of the fibre is slightly above the plane, while the other half lies below. For short times, while the fibres are largely in the plane, they experience shape deformations as in linear shear – stiffer fibres remain straight, while softer fibres exhibit the signature S-buckling about the origin. In cases with larger
$\eta$
values, the curvature of this initial S-buckling is larger, and the fibre is more closely concentrated near the origin. For this buckled fibre, even while mostly in the plane, there is a tug-of-war due to the axial velocity in opposite directions on each half of the fibre. This ‘pulling’ in opposite directions is more pronounced for fibres concentrated near the origin. We can quantify this geometric feature by using the gyration tensor to describe the evolving fibre shapes in our simulations, as in Liu et al. (Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018). Similar metrics are commonly used to characterise polymer shapes and conformation (Vymětal & Vondrášek Reference Vymětal and Vondrášek2011; Arkın & Janke Reference Arkın and Janke2013). We compute the squared radius of gyration
$ (R_{\textit{g}}^2 )$
by summing the eigenvalues of the gyration tensor (defined in Appendix A.2). The radius of gyration has the largest values when a fibre is straight and gets smaller as the fibre bends and more of its material points get closer to its centroid. For instance, in figure 7(a), at
$\tilde {t}=10$
the fibres have radii of gyration of
${5.84}$
,
${5.31}$
,
${4.01}$
and
${3.59}\,\mathrm{\unicode{x03BC} m}$
, respectively, with decreasing
$EI$
order.
Figure 8 shows the time evolution of the radius of gyration of the fibres in each of the 19 simulations for
$d/L=0.5$
whose alignment time is depicted in figure 7(d) where all 19 cases were subject to the same background vortex with
$\epsilon ={5}\,\mathrm{\unicode{x03BC} m}$
,
$\rho = \sigma ={0.71}\,\mathrm{g\,\unicode{x03BC} m^2\, s^{-2}}$
and
$\omega _{\textit{max}}={2.28\times 10^3}\,\textrm{s}^{-1}$
. Comparing the minimum radii of gyration
$ (R_{\textit{g,min}} )$
attained during the orbit from a horizontal straight configuration to a vertically aligned state, we see
$R_{\textit{g,min}}$
decreasing as
$EI$
is reduced. Within a given simulation, we see that the drop in
$R_g$
is also faster for floppier fibres. This ‘pliability’ of flexible fibres allows them to readily align with background vortex streamlines, leading to small alignment times. We find very little difference in the minimum radius of gyration for the stiffest fibres (darkest curves); hence the alignment times do not significantly vary. The inset of figure 8, that plots
$R_{\textit{g,min}}$
versus the bending rigidity (log scale), exhibits a logistic-like behaviour. The transition region (in figure 7
d) is where
$R_g$
and
$R_{\textit{g,min}}$
vary significantly as
$EI$
is changed. In general, inside a spiralet, the floppier fibres attain a smaller
$R_{\textit{g,min}}$
, reach a minimum faster than stiffer fibres and, hence, take a shorter time to complete their orbits.
Figure 8. Shape change of fibres reflected by the evolution of radius of gyration
$R_g$
. Note that the fibres are in straight configuration at the start and end of each simulation. This is indicated by the leftmost and the rightmost regions in the figure where the values of
$R_g$
, regardless of the bending stiffness, stay near the theoretical value of
$L/\sqrt {12}$
for a slender straight rod (see table 19-1 in Feynman, Leighton & Sands (Reference Feynman and Leighton1989)). At their most deformed state, softer fibres fit in a more compact spheroidal envelope, where the radius of gyration drops to as low as 30 % of the straight configuration. The inset shows the minimum attained radius of gyration
$R_{\textit{g,min}}$
as a function of bending stiffness. Representative shapes corresponding to four different stiffness cases are shown on the right.
3.2. Asymmetrically placed fibres
The results presented in the previous section were for fibres initialised with their centroid coinciding with the spiralet centre. This symmetry allowed the fibres to reach a vertical equilibrium position. We now demonstrate dynamics in a simulation where this symmetry is broken. Figure 9(a) shows an
$L= {20}\,\mathrm{\unicode{x03BC} m}$
fibre placed such that its centroid is shifted radially outward from the origin along its length by
${2}\,\mathrm{\unicode{x03BC} m}$
. In the ensuing rotation and deformation stages, although the fibre quite easily reaches the vertical alignment with the vortex axis, more than half of the fibre is now above the vortex plane, and starts moving upward with the background flow (third snapshot of figure 9
a at
$t={0.045}\,\mathrm{ s}$
). Some instances later, it suddenly buckles into a helical shape while still moving vertically upward as a whole. We use the curvature and torsion kymographs to piece together the deformation dynamics. Compared with the curvature kymograph for the symmetric case (figure 3
b), the asymmetry in lengthwise bending is immediately visible in figure 9(b). The bending phase, indicated by the diagonally moving high-curvature zones in the curvature kymograph, ends around
$t={0.04}\,\mathrm{ s}$
. From
$t={0.04}\,\mathrm{ s}$
to
$t={0.065}\,\mathrm{ s}$
, the fibre remains straight (zero-curvature region in the midsection of figure 9
b). However, we find that during this same time period, high levels of torsion build up in the midsection of the filament (figure 9
c). At
$t={0.065}\,\mathrm{ s}$
, a helicoidal buckling appears, which shows up as a sudden release in torsion and multiple simultaneous appearances of high-curvature regions along the length of the filament. We hypothesise that this buckling behaviour is induced by two different mechanisms. Build-up of local torsion along the length is known to cause helical buckling and plectoneme formation with accompanying release in torsion (Charles, Gazzola & Mahadevan Reference Charles, Gazzola and Mahadevan2019). At the same time, in our vortex set-up, the fibre is experiencing axial compression due to the upper part of the fibre being further away from the vortex centre, where the axial speed (in the positive
$z$
direction) is much smaller compared with the lower end, which is much closer to the vortex centre. Similar helical buckling was also reported in Chakrabarti et al. (Reference Chakrabarti, Liu, LaGrone, Cortez, Fauci, du Roure, Saintillan and Lindner2020) for fibres in compressional flows.
Figure 9. (a) Snapshots of an asymmetrically placed fibre with its centroid radially shifted from the origin by a small amount
$({2}\,\mathrm{\unicode{x03BC} m})$
(see supplementary movie 5). (b) Temporal evolution of local curvature and (c) local torsion along the fibre. The colour on the fibre surface in (a) indicates local torsion. Here
$d/L$
= 1 and
$\eta = {1.18\times 10^6}$
.
We next examined the effect of bending stiffness on helix formation in the cases when fibres of
$L = {20}\,\mathrm{\unicode{x03BC} m}$
were initialised in the vertical orientation with their centroid slightly offset at
$(x=0.06,y=0.14,z=3.72) \,\mathrm{\unicode{x03BC} m}$
(instead of at the origin). In each case, the background spiralet was fixed so that
$d/L= 1$
, but the fibre bending rigidity, and hence
$\eta$
, varied. Figure 10 shows snapshots from five different simulations at the same instant in time. While all fibres move upwards in space, we see that the stiffest fibre remains straight and the next stiffest forms a C-shape. The softest three cases buckle into helices whose characteristic wavelengths decrease as
$\eta$
increases, as seen in Manikantan & Saintillan (Reference Manikantan and Saintillan2015), Chakrabarti et al. (Reference Chakrabarti, Liu, LaGrone, Cortez, Fauci, du Roure, Saintillan and Lindner2020) and Sznajder et al. (Reference Sznajder, Zdybel, Liu and Ekiel-Jeżewska2024).
Figure 10. Effect of decreasing bending rigidity in influencing the buckling modes for fibres placed asymmetrically in the vortex. The stiffest case in (a) has
$EI = {5\times 10^{-23}}\,{\textrm{J m}}$
with the bending rigidity decreasing to
$EI = {5\times 10^{-24}}\,{\textrm{J m}}$
in (b),
$EI = {1\times 10^{-24}}\,{\textrm{J m}}$
in (c),
$EI = {4\times 10^{-25}}\,{\textrm{J m}}$
in (d) and
$EI = {1\times 10^{-25}}\,{\textrm{J m}}$
in (e). All background flow parameters and filament properties except for the bending rigidity are held constant.
Figure 11. Mean fibre orientation angle (
$\theta$
) for an arbitrarily deformed fibre. The cyan circle has the same radius as the fibre’s radius of gyration
$(R_g)$
.
4. Conclusion
We have introduced a simple superposition of two regularised solutions to the Stokes equations to produce a zero-Reynolds-number analogue of a Burgers vortex. This vortical flow has full three-dimensional features such as spatially varying vorticity, compression and stretching. Here, we performed computational studies of the shape deformations and orbits of flexible fibres immersed in this background flow. We mainly focused our attention on the special case where fibres were placed in the horizontal plane of symmetry of this spiralet vortex, and still exhibited a rich library of shape deformations, some of which are not attainable by applying a simple shear flow. These shapes were determined by two non-dimensional parameters: the elastoviscous number and the ratio of vortex diameter to fibre length. Moreover, in this special case of fibre placement, we demonstrated a robust spinning-to-standing orbit, whose completion time decreased with increasing fibre flexibility. Simulations with fibres initialised with small asymmetries result in an even wider range of shape deformations, including helical buckling. We also remark that this Kirchhoff rod framework allows us to track torsion along the fibre, and we can measure how the accumulation of torsion on a straight fibre gives way to buckling. Also, motivated by Yu & Graham (Reference Yu and Graham2021), we expect that shape deformations of two-dimensional sheets and membranes within these spiralet flows can be investigated. In addition, because of the simplicity of the spiralet background flow, we will investigate fibre dynamics in an ensemble of spiralets with random orientations and strengths. The hope is that this simple model could be used to assess the effects of turbulence on the transport and nutrient uptake of planktonic organisms at the microscale.
Promising progress in laboratory realisations of Burgers vortices (Webster & Young Reference Webster and Young2015; Aulnette et al. Reference Aulnette, Burshtein, Banaei, Brandt, Haward, Shen, Delmotte and Lindner2025) and flexible fibre deformation in homogeneous and wall-bounded turbulence (Marchioli et al. Reference Marchioli, Rosti and Verhille2025) has recently been made. We hope that the simulations presented here motivate laboratory experiments of flexible fibres in these fully three-dimensional vortical flows.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11342
Funding
This research was supported by the International Human Frontier Science Program Organization (grant no. RGP0017/2022; https://doi.org/10.52044/HFSP.RGP00172022.pc.gr.153610), the US National Science Foundation (grant no. DMS 2054333) and the Simons Foundation (grant no. SFI-MPS-SFM-00006482). L.F. also acknowledges the support of the Institut Henri Poincaré (UAR 839 CNRS-Sorbonne Université) and LabEx CARMIN (ANR-10-LABX-59-01).
Declaration of interests
The authors report no conflict of interest.
Appendix A
A.1. Velocity components of the spiralet in cylindrical coordinate system
For the background spiralet velocity field defined in (2.4), the azimuthal, axial and radial components can be expressed as
\begin{align} u_{\theta }^b(\tilde {r},\theta ,z,\rho ) &=\frac {\tilde {r}(5\epsilon ^2+2(\tilde {r}^2+z^2))}{8\pi \mu (\tilde {r}^2+z^2+\epsilon ^2)^{5/2}}\rho , \end{align}
\begin{align} u_z^{b}({\tilde {r}},\theta ,z,\sigma ) &= \frac {z(2\epsilon ^2-\tilde {r}^2+2z^2)}{8\pi \mu (\tilde {r}^2+z^2+\epsilon ^2)^{5/2}}\sigma , \end{align}
\begin{align} u_{\tilde {r}}^{b}({\tilde {r}},\theta ,z,\sigma ) & = \frac {\tilde {r}(-\epsilon ^2-\tilde {r}^2+2z^2)}{8\pi \mu (\tilde {r}^2+z^2+\epsilon ^2)^{5/2}}\sigma , \end{align}
where
$\tilde {r}$
is the radial distance from the vertical axis and
$\rho$
and
$\sigma$
are the strengths of the rotlet and stresslet, respectively.
A.2. Computing radius of gyration and fibre orientation
For a fibre with
$N$
segments and
$N+1$
discretised points (
$\boldsymbol{X}_i(t), i=0,\ldots ,N$
), the gyration tensor is defined as
\begin{align} \boldsymbol{G}(t) = \frac {1}{N+1} \!\begin{bmatrix} \sum _i (x_i(t) - x_{\textit{cm}})^2 & \sum _i (x_i(t) - x_{\textit{cm}})(y_i(t) - y_{\textit{cm}}) & \sum _i (x_i(t) - x_{\textit{cm}})(z_i(t) - z_{\textit{cm}}) \\ \sum _i (x_i(t) - x_{\textit{cm}})(y_i(t) - y_{\textit{cm}}) & \sum _i (y_i(t) - y_{\textit{cm}})^2 & \sum _i (y_i(t) - y_{\textit{cm}})(z_i(t) - z_{\textit{cm}}) \\ \sum _i (x_i(t) - x_{\textit{cm}})(z_i(t) - z_{\textit{cm}}) & \sum _i (y_i(t) - y_{\textit{cm}})(z_i(t) - z_{\textit{cm}}) & \sum _i (z_i(t) - z_{\textit{cm}})^2 \end{bmatrix}\!, \end{align}
where the centre of mass,
$(x_{\textit{cm}},y_{\textit{cm}},z_{\textit{cm}})$
, always remains at the spiralet core when symmetry is imposed by initial fibre placement. The squared radius of gyration is computed from the eigenvalues,
$\lambda _1, \lambda _2, \lambda _3$
, of
$\boldsymbol{G}$
as
\begin{equation} R_g^2=\lambda _1+\lambda _2+\lambda _3 .\end{equation}
The principal axis of the fibre is taken to be in the direction given by the eigenvector (
$\boldsymbol{V}_1$
) corresponding to the largest eigenvalue of
$\boldsymbol{G}$
. At any instant, the mean orientation of the fibre
$(\theta )$
is defined by the acute angle between its principal axis and the projection of the principal axis on the spiralet vortex symmetry plane.
Open access
