Attracting dynamical modes of highly elastic fibres settling under gravity in a viscous fluid

Abstract

The dynamics of a single highly elastic fibre settling under gravity in a very viscous fluid is studied numerically. We employ the bead model and multipole expansion of the Stokes equations, corrected for lubrication that is implemented in the precise Hydromultipole numerical codes. Four attracting regular dynamical modes of highly elastic fibres are found: two stationary shapes (one translating and the other rotating and translating), and two periodic oscillations around such shapes. The phase diagram of these modes is presented. It illustrates that the existence of each mode depends not only on the elasto-gravitation number but also on the fibre aspect ratio. Characteristic time scales, fibre deformation patterns and motion in the different modes are determined.

1. Introduction

Dynamics of flexible elongated objects in viscous fluid flows has been researched extensively for some time because of the importance of such objects in many biological systems and technological processes, especially at the scale of nanometres and micrometres. Examples could be diatom chains in the ocean (Nguyen & Fauci 2014), cellulose fibres in paper manufacturing (Hubbe & Heitmann 2007) or fibre hydrogels for wound healing and drug delivery (Perazzo et al. 2017). Deformation, orientation and motion of flexible fibres and flexible planar bodies have been studied in various fluid flows at small (Du Roure et al. 2019; Yu & Graham 2024), moderate (Shelley & Zhang 2011) and large Reynolds numbers, including turbulence (Dotto & Marchioli 2019). Instabilities have been investigated for sedimenting suspensions of flexible filaments (Manikantan et al. 2014; Manikantan & Saintillan 2016; Banaei et al. 2020).

It is known that the dynamics of not only flexible but also rigid objects of various shapes settling under gravity in the laminar regime are usually very complex (Witten & Diamant 2020), often owing to non-negligible, small inertia effects. For example, it was shown numerically (Ravichandran & Wettlaufer 2023) that the dynamics of a concavo–convex rigid body settling under gravity in a viscous fluid changes significantly with the increase of the Reynolds number $Re$. When consecutive critical threshold values of $Re$ are exceeded, in the range of O(1)–O(10$^2$), one stationary orientation changes into two, then back to one and next to an attracting limit cycle. Even more pronounced dependence on the Reynolds number and shape (different types of periodic motions and a chaotic behaviour) has been observed for the dynamics of sedimenting: thin curved rigid objects (experiments at $Re>100$) (Chan et al. 2021) and rigid discs (numerical simulation at $25 \le Re \le 300$) (Auguste, Magnaudet & Fabre 2013). Experiments performed at $31 < Re < 259$ (Lorite-Díez et al. 2022) showed that a settling elongated, finite-length, flexible but rather stiff cylinder exhibits periodic rigid-body oscillations and periodic bending oscillations, dependent on $Re$.

Systems of elastic microfilaments settling under gravity at the Reynolds number much smaller than unity are of special interest. To study them, a variety of simulation methods have been applied, including the slender-body theory (Xu & Nadim 1994; Li et al. 2013; Waszkiewicz, Szymczak & Lisicki 2021) and the bead-chain model (Cosentino Lagomarsino, Pagonabarraga & Lowe 2005; Llopis et al. 2007; Saggiorato et al. 2015; Bukowicki & Ekiel-Jeżewska 2018, 2019; Gruziel-Słomka et al. 2019; Shashank, Melikhov & Ekiel-Jeżewska 2023). The lack of inertia seems to simplify the dynamics of elastic fibres with large and moderate bending stiffness. Numerical simulations (Cosentino Lagomarsino et al. 2005; Schlagberger & Netz 2005; Li et al. 2013; Saggiorato et al. 2015; Shojaei & Dehghani 2015; Bukowicki & Ekiel-Jeżewska 2018) and experiments together with simulations (Marchetti et al. 2018; Cunha et al. 2022), without and with Brownian motion, for elastic filaments of large and moderate bending stiffness indicate the existence of stable, stationary, planar, U-shaped configuration oriented vertically.

However, numerical simulations (Cosentino Lagomarsino et al. 2005; Saggiorato et al. 2015; Shashank et al. 2023) and experiments (Shashank et al. 2023) with highly elastic fibres indicate that even for $Re \ll 1$, the U-shaped configurations are no longer stationary. A variety of different dynamical modes, with complex out-of-plane shapes, both stationary and with time-dependent deformations, have been reported for a highly elastic sedimenting loop (Gruziel-Słomka et al. 2019). Periodic deformations of a highly elastic fibre sedimenting in two-dimensional fluid have been reported (Shojaei & Dehghani 2015).

Therefore, in this work, we focus on the long-time three-dimensional dynamics of a highly elastic fibre settling under gravity in a very viscous fluid at the Reynolds number $Re=0$, i.e. within the Stokes regime. We look into the existence of different attracting dynamical modes. In § 2, the theoretical and numerical approach is outlined. Section 3 contains a phase diagram of different attracting dynamical modes, depending on the fibre aspect ratio and elasto-gravitation number, which estimates the ratio of gravitational and bending forces. The characteristic properties of these modes are discussed in §§ 4 and 5. The conclusions are presented in § 6.

2. Theoretical description and parameters of the simulation

The fibre is modelled as a chain of $N$ identical spherical beads of a diameter $d$, with the centres of the consecutive beads connected by springs. We employ a finitely extensible nonlinear elastic (FENE) model (Warner 1972; Bird, Armstrong & Hassager 1977) as the stretching potential energy,

(2.1)\begin{equation} U^{FENE} ={-} \frac{1}{2} k (l_{0}-d)^{2} \sum_{i=1}^{N-1}\ln \left[ {1 - \left( \frac{l_{0}-l_{i}}{l_{0}-d} \right)^{2}} \right]. \end{equation}

Here, $k$ is the spring constant, $l_{i}=|\boldsymbol {r}_{i+1}-\boldsymbol {r}_{i}|$ is the time-dependent distance between centres of two consecutive beads at the positions $\boldsymbol {r}_{i+1}$ and $\boldsymbol {r}_{i}$ and

(2.2)\begin{equation} l_0=1.01d \end{equation}

is the length of each spring at elastic equilibrium. In our simulations, the distance $l_i(t)$ between the centres of the consecutive beads stays relatively close to the equilibrium value $l_0$ for all the times $t$: typically, $\max _{t,i}l_i(t) \lesssim 1.017$ and $\min _{t,i}l_i(t) \gtrsim 1.003$.

Each triplet of the consecutive beads $j-1, j, j+1$ is straight at the elastic equilibrium, and it resists bending with the potential energy,

(2.3)\begin{equation} U^{b}_j = \frac{A}{2 l_{0}} \beta_{j}^2, \end{equation}

where $A$ is the bending stiffness and $\beta _{j}$ is the bending angle between $\boldsymbol {r}_{i}-\boldsymbol {r}_{i-1}$ and $\boldsymbol {r}_{i+1}-\boldsymbol {r}_{i}$.

The harmonic potential is commonly used in the literature to represent the bending potential energy in molecular modelling, in particular of proteins and DNA, as described, e.g. by MacKerell et al. (1998), Storm & Nelson (2003) and Frenkel & Smit (2023). For sedimenting elastic fibres, the Kratky–Porod potential is commonly used by most authors, including Schlagberger & Netz (2005), Manghi et al. (2006), Marchetti et al. (2018), Llopis et al. (2007) and Saggiorato et al. (2015). However, Bukowicki & Ekiel-Jeżewska (2018) showed that sometimes for larger bending angles, the Kratky–Porod potential may lead to spurious dynamics. For highly elastic fibres, large bending angles are expected, and therefore in this work, we use the harmonic bending potential that at large bending angles binds the fibre segments much stronger than the Kratky–Porod potential. It was checked by Bukowicki & Ekiel-Jeżewska (2018) and Bukowicki (2019) that the harmonic bending potential for highly elastic trumbbells does not lead to spurious dynamics. It was done by comparing with a logarithmic potential that blows up for increasing bending angles.

The bending stiffness

(2.4)\begin{equation} A = \frac{E_Y {\rm \pi}d^4}{64}=\frac{k d^{2} l_{0}}{ 16} \end{equation}

is related to the spring constant $k$ by the model of an elastic cylinder of diameter $d$ with the Young's modulus $E_{Y}$ (Bukowicki & Ekiel-Jeżewska 2018, 2019).

In addition to the elastic forces on a bead $i$,

(2.5)\begin{equation} \boldsymbol{F}_i^e ={-} \frac{\partial \left( U^{FENE} + \displaystyle\sum_{j=2}^{N-1} U^{b}_j \right)}{\partial \boldsymbol{r}_i}, \end{equation}

there is also a constant gravitational force that acts on each bead along the $z$-axis:

(2.6)\begin{equation} \boldsymbol{F}^g ={-} m g\, \hat{\boldsymbol{e}}_z. \end{equation}

Here, $m$ is the bead mass corrected for buoyancy, $g$ is the acceleration due to gravity and $\hat {\boldsymbol {e}}_z$ is the unit vector along the $z$-axis.

The fluid flow is assumed to obey the Stokes equations. The time-dependent positions $\boldsymbol {r}_i$ of the bead centres, $i=1,\ldots,N$, satisfy the set of the first-order ordinary differential equations,

(2.7)\begin{equation} \dot{\boldsymbol{r}}_{i} = \sum_{j=1}^{N} \boldsymbol{\mu}_{ij} \boldsymbol{\cdot} (\boldsymbol{F}_j^e + \boldsymbol{F}^g ). \end{equation}

Here, the 3-× Cartesian mobility matrix $\boldsymbol {\mu }_{ij}$ accounts for hydrodynamic interactions between all the beads and is evaluated from the multipole expansion by the precise numerical codes Hydromultipole (Cichocki, Ekiel-Jeżewska & Wajnryb 1999; Ekiel-Jeżewska & Wajnryb 2009). The multipole truncation order of 3 (Cichocki et al. 1994) was used in these computations, as previous modelling on similar systems showed that the improvement in accuracy was practically negligible when the multipole truncation order of 4 was used (Shashank et al. 2023).

A dimensionless elasto-gravitation number $B$ is introduced to relate gravitational and bending forces acting on the fibre (Cosentino Lagomarsino et al. 2005; Llopis et al. 2007; Saggiorato et al. 2015; Bukowicki & Ekiel-Jeżewska 2018, 2019; Marchetti et al. 2018),

(2.8)\begin{equation} B= L_0^2 G / A, \end{equation}

where $G=N m g$ is the gravity corrected for buoyancy acting on the whole fibre and $L_0 = (N-1) l_{0} + d$ is the equilibrium length of the fibre. With (2.2), the fibre aspect ratio is well-approximated as $N$.

In this study, the total number of beads varies in the range $12 \le N \le 40$, and the elasto-gravitation number is in the range $1000 \le B \le 11\,000$. Our choice is related to the results of Llopis et al. (2007) and Saggiorato et al. (2015), who assumed a moderate aspect ratio $N=30$. We extended for a range of $N$ around this value. In addition, the choice of $B$ is motivated by the above references. We include the upper limits for the range of values $B\le 5000$ from Llopis et al. (2007), and $B\le 3000$ from Saggiorato et al. (2015) and we also go above them to reach highly flexible fibres with much larger values of $B$.

Our choice of the initial conditions is also motivated by previous results. Cosentino Lagomarsino et al. (2005) and Llopis et al. (2007) investigated an elastic fibre initially straight and horizontal, and observed the evolution restricted to a vertical plane. Saggiorato et al. (2015) pointed out that highly elastic fibres (contrary to moderately elastic ones) are unstable to a perturbation out of this vertical plane, and demonstrated that the dynamics of highly elastic fibres with and without such perturbation are different (in-plane and out-of-plane, respectively).

Therefore, in this work, we follow Saggiorato et al. (2015) and apply a weak noise on the initially straight and horizontal fibre in the elastic equilibrium, with $l_{i}(0)=l_0$. The amplitude of the random perturbation is equal to $0.0001 d$. The evolution of the fibre is monitored until time $10^6\tau _b/N$, with $\tau _b = {\rm \pi}\eta d^2 / ( mg )$ chosen as the time unit, with $\eta$ being the dynamic viscosity of the fluid. For almost all cases, it takes nearly half of this time before the attracting mode is reached.

3. Phase diagram of different dynamical modes

For moderately elastic fibres, a vertical mode, with a stationary, symmetric V-like (or U-like) shape located in a vertical plane, is the only mode observed (Cosentino Lagomarsino et al. 2005; Schlagberger & Netz 2005; Saggiorato et al. 2015; Marchetti et al. 2018). For highly elastic fibres, i.e. for large values of $B$, the tilted and rotation modes, with stationary out-of-plane shapes, were reported by Saggiorato et al. (2015) (see figure 1 therein). Here, we confirm this finding and find new families of attracting dynamical modes (we call them crawling and rotation-crawling), using a different theoretical and numerical method. We also show that for very flexible fibres, the dynamics often become irregular.

Physically, there are two different reasons for large fibre flexibility: small ratio of the local elastic to the local gravitational forces, or large aspect ratio. In general, the dynamics are expected to depend on both these parameters. For moderately elastic fibres, the dynamics are often reported in the literature to depend only on a single parameter: the elasto-gravitation number that depends on both physical parameters. Indeed, in the phase diagram in figure 1, we observe that the transition between vertical and tilted mode corresponds to the same value of $B$ in the wide range of the aspect ratios $N$.

Figure 1. Attracting dynamical modes achieved by an elastic fibre with aspect ratio $N$ and elasto-gravitation number $B$.

However, the important finding of our manuscript is that the dynamics of highly elastic fibres depend on both physical parameters. Figure 1 is the map of the attracting dynamical modes of the highly elastic fibre, based on a systematic study of the dynamics in a wide range of the aspect ratio $N$ and the elasto-gravitation number $B$. The borders between the modes in the phase diagram in figure 1 indicate that there is no single parameter to determine the dynamics: two are needed.

It is interesting to study the transition from the tilted to the crawling mode, shown in figure 1, and extract the transition value $B_T$ of the elasto-gravitation number $B$ as a function of the transition value $N_T$ of the fibre aspect ratio $N$. Figure 2 depicts the transition points $B_T(N_T)$, which were identified employing the procedure described in Appendix A.

Figure 2. Transition elasto-gravitation number $B_{T}$ and corresponding aspect ratio $N_T$ for the transition from the tilted to the crawling mode. Data points represent values obtained from numerical simulations. The solid line depicts the fit based on (3.1). Error bars are comparable to the size of the symbol.

We expect a power law (Ekiel-Jeżewska 2014) dependence of the form

(3.1)\begin{equation} B_T - B_c = \alpha (N_T - N_c)^p. \end{equation}

We predict that there might exist a critical value of $N_c$ with qualitatively different dynamics for $N< N_c$ and $N>N_c$. This hypothesis is based on the dynamics of elastic loops, investigated numerically (Gruziel-Słomka et al. 2019). It was shown there that for a wide range of the bending stiffness $A$, sedimenting elastic loops made of $N=15$ overlapping beads (with the aspect ratio 9.4) orient horizontally and practically do not deform, in contrast to the loops with $N\ge 16$ and the same $A$, which tend to orient vertically, stay inclined or deform significantly.

The nonlinear least-square fit to the data presented in figure 2 was performed based on (3.1). The identified parameters $B_{c} = 3250 \pm 140$, $N_{c} = 11.8 \pm 0.4$, $\alpha = 44\,000 \pm 14\,000$ and $p = -1.32 \pm 0.14$ yield a function that adequately describes the data, as illustrated in figure 2. To determine the parameters in the power law (3.1) with greater precision, an estimation of $B_{T}$ with significantly higher accuracy is needed, which would require a significantly larger number of simulations near the transition values, and also simulations with the fibres having smaller numbers of beads.

4. Basic properties of the attracting dynamical modes

4.1. Comparison of typical shapes

Typical shapes of highly elastic fibres in different sedimentation modes are shown in figure 3. In the tilted and rotation modes, the shapes do not depend on time. In the crawling and rotation-crawling modes, the shapes change in time with a period $T$.

Figure 3. Snapshots of the fibre made of $N=36$ beads in different dynamical modes, reached after a long time. Positions of the fibre beads are normalised by their diameter. (a) The notation used. (b,c) Fixed shapes in the tilted ($B=3500$) and rotation ($B=6500$) modes, respectively. (d,e) Periodic evolution of shapes in the crawling ($B=7500$) and rotation-crawling ($B=8500$) modes, respectively. The reference frames are given in (4.2a–c) for (b), (4.3) for (c), (4.4a–c) for (d) and (4.5) for (e). See also the supplementary movies available at https://doi.org/10.1017/jfm.2024.729.

To quantify properties of the sedimentation modes, it is useful to compute the moment of inertia tensor (alternatively, the gyration tensor $({1}/{N})\sum _{j=1}^{N} r'_{j\alpha } r'_{j\beta }$ is often evaluated),

(4.1)\begin{equation} {I}_{\alpha \beta}=\sum_{j=1}^{N} (\delta_{\alpha \beta} {{r}'_j}^{2} - r'_{j\alpha} r'_{j\beta}), \end{equation}

where $\alpha = x, y, z,\ \beta = x, y, z$ and $\boldsymbol {r}'_j=(r'_{jx},r'_{jy},r'_{jz})$ is the position of $j$th bead centre in the centre-of-mass reference frame, $\boldsymbol {r}'_j=\boldsymbol {r}_j-\boldsymbol {r}_{CM}$, with the centre-of-mass position $\boldsymbol {r}_{CM}=(x_{CM},y_{CM},z_{CM})$ and $r_j'=|\boldsymbol {r}'_j|$. Knowledge of the moment of inertia tensor allows evaluation of its eigenvectors and eigenvalues. A unit vector $\boldsymbol {n}$, which is parallel to the eigenvector associated with the largest eigenvalue, can be introduced. For a planar shape, $\boldsymbol {n}$ would be perpendicular to this plane. For shapes slightly out of plane, $\boldsymbol {n}$ is perpendicular to a certain ‘average plane’ and determines its orientation. The unit vector $\boldsymbol {n}$ can be expressed in spherical coordinates with the polar angle $\theta$ that it makes with the vertical axis $z$ (parallel to gravity) and the azimuthal angle $\phi$ that its projection onto $xy$-plane makes with the $x$-axis (see figure 3a). The $x$ and $y$ axes are chosen arbitrarily and are not related to the initial fibre orientation that usually changes during the transient time of the evolution. In the following, the angles $\theta$ and $\phi$ are determined as functions of time and are used to characterise the dynamical modes.

The time-dependent end-to-end vector $\boldsymbol {r}_{1N} = \boldsymbol {r}_{N} - \boldsymbol {r}_{1}$ provides more information about the fibre shape and motion in different modes. Here we are interested in the radius $\rho$ and the altitude $\Delta z$, i.e. the cylindrical coordinates of $\boldsymbol {r}_{1N}$, with the cylindrical axis parallel to gravity and the origin at the point $\boldsymbol {r}_1$. The projection of the trajectory of the tip of the end-to-end vector $\boldsymbol {r}_{1N}$ on the plane of the cylindrical coordinates $\rho$ and $\Delta z$ will be analysed in § 5 to show the difference between the modes and bifurcations they undergo.

4.2. Tilted mode

The fibre shape in the tilted mode does not change with time, and it is not planar, as shown in figure 3(b) (compare with the similar numerical results in figure 1 of Saggiorato et al. 2015). The mode name, tilted, indicates that the ‘averaged plane’ of the fibre shape is inclined with respect to gravity. In figure 3(b), the horizontal axes of the side and front planes and vertical axis correspond respectively to

(4.2a–c)\begin{equation} (x-x_{CM})/d,\quad (y-y_{CM})/d,\quad (z-z_{CM})/d. \end{equation}

More information about the fibre shape, orientation and motion is provided in figure 4 (blue dashed lines). Its orientation is fixed: the polar angle $\theta$ and the azimuthal angle $\phi$ are constant in time, as shown in figure 4(a,d). The shape is symmetric with respect to reflection in the vertical plane parallel to the side plane of the snapshot shown in figure 3(b). In particular, the fibre ends are at the same horizontal level, with $\Delta z =0$, as visible in figure 4(b). Due to the inclined shape, i.e. the non-zero inclination angle $\theta$, the fibre moves horizontally along $x$. Therefore, the top view of the centre-of-mass trajectory is a straight line, shown in figure 4(c).

Figure 4. Basic features of the regular (a–d) and irregular (e,f) modes of fibre with $N=36$ (one example of each mode). Time dependence of (a) polar angle $\theta$, (b) vertical coordinate $\Delta z$ of the fibre end-to-end vector $\boldsymbol {r}_{1N}$, and (d) azimuthal angle $\phi$. Top view (c) of the centre-of-mass trajectory (the ranges of times are $1150 \tau _{b}$ for tilted, $1545 \tau _{b}$ for rotation, $1800 \tau _{b}$ for crawling and $3700 \tau _{b}$ for rotation-crawling modes). Evolution of an irregular mode: (e) azimuthal angle $\phi$ (with the results for crawling and rotation-crawling modes shown for comparison) and ( f) top view of the centre-of-mass trajectory.

4.3. Rotation mode

The rotation mode is another mode in which the shape of the fibre does not change with time. The shape is chiral and, therefore, it rotates around the gravity direction $z$ with a certain period $T_R$. In the reference frame with the origin on the rotation axis of the fibre centre of mass, $x_{CM}=R \cos ({2 {\rm \pi}t}/{T_{R}})$ and $y_{CM}=R \sin ({2 {\rm \pi}t}/{T_{R}})$, where $R$ is the radius of rotation. Projections of the fibre shape in the reference frame rotating and translating with the fibre and the origin on the rotation axis are shown in figure 3(c). Here, the horizontal axes of the side and front planes and the vertical axis correspond to

(4.3) \begin{equation} \left(\begin{matrix} \cos\left(\tfrac{2 {\rm \pi}t}{T_{R}}\right) & \sin\left(\tfrac{2 {\rm \pi}t}{T_{R}}\right) & 0 \\ -\sin\left(\tfrac{2 {\rm \pi} t}{T_{R}}\right) & \cos\left(\tfrac{2 {\rm \pi}t}{T_{R}}\right) & 0\\ 0 & 0 & 1 \end{matrix}\right)\boldsymbol{\cdot} \left(\begin{matrix} \left(x-R\cos\left(\dfrac{2 {\rm \pi} t}{T_{R}}\right)\right)/d \\ \left(y-R \sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right)\right)/d \\ (z-z_{CM})/d \end{matrix}\right), \end{equation}

respectively.

More information about the fibre shape, orientation and motion is provided in figure 4 (orange dotted lines). The polar angle $\theta$ remains constant in time, as visible in figure 4(a). The azimuthal angle $\phi$ increases with a constant slope that defines the characteristic period of rotation, with $T_{R}=1545 \tau _{b}$ for the orange dotted line in figure 4(d). A non-zero, time-independent value of the difference $\Delta z$ between vertical positions of the fibre ends can be seen in figure 4(b). The fibre centre of mass follows a helical trajectory with a circular cross-section, shown in figure 4(c). Another example of the fibre motion in the rotation mode was evaluated numerically by a different method (Saggiorato et al. 2015).

4.4. Crawling mode

The crawling mode has not been reported in the literature so far. The shape of the fibre in the crawling mode does not remain fixed but changes periodically with a period $T$. The snapshots of the fibre shape, taken at time intervals equal to $T/8$, are shown in figure 3(d), where the horizontal axes of the side and front planes and the vertical axis correspond respectively to

(4.4a–c)\begin{equation} (x-x_{CM})/d,\quad y/d,\quad (z-z_{CM})/d. \end{equation}

The fibre appears to be making crawling motions with its arms, hence the name of the mode. First, for half of the period $T$, the left arm is positioned above the right one and, then, for the next half of the period $T$, the right arm becomes above the left one. Moreover, the fibre shape at time $t +T/2$ is the mirror reflection of the shape at time $t$, with the reflection plane parallel to the side plane, as seen in figure 3(d).

More information about the fibre shape, orientation, and motion is provided in figure 4 (green solid lines). It is clear already in figure 3(d) that the shape is inclined, with the inclination angle $\theta > 0$, displayed in figure 4(a). This inclination causes a lateral drift of the fibre centre of mass along $x$, with small periodic side movements along $y$, as seen in figure 4(c). The azimuthal angle $\phi$ and the difference $\Delta z$ between the vertical positions of the fibre ends are periodic functions with the same period $T= 440 \tau _{b}$ (see the solid green lines in figure 4b,d). The polar angle $\theta$ is a periodic function with the period $T/2$ (see figure 4a).

4.5. Rotation-crawling mode

The rotation-crawling mode has not been reported in the literature so far. It is a more complex mode than those discussed previously. The rotation-crawling mode seems to be a quasi-periodic motion with two incommensurable characteristic time scales, $T_R$ and $T$. The shape is chiral, and the fibre rotates around $z$ with a period $T_R$, which is the time it takes for the time-averaged azimuthal angle $\phi$ to change by $2 {\rm \pi}$. In addition, similar to the crawling mode, the shape of the fibre changes with time, allowing to introduce a characteristic time $T$ required for the fibre to return to its initial conformation after completing a full cycle of deformation. To identify the period $T$, the time dependence of either $\theta$ or $\Delta z$ is used.

The period $T_R$ of the rotation-crawling mode is typically larger than the rotation period of the rotation mode with similar values of $B$ and $N$. Therefore, the shape oscillations influence the rotational dynamics.

The sequence of snapshots of the fibre shape, taken at the time intervals $T/8$, are shown in figure 3(e). The moving frame of reference is chosen in such a way that the horizontal axes of the side and front planes and the vertical axis are specified using the following equation:

(4.5)\begin{equation} \left(\begin{matrix} \cos\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & \sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & 0 \\ -\sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & \cos\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right) & 0 \\ 0 & 0 & 1 \end{matrix}\right)\boldsymbol{\cdot} \left(\begin{matrix} (x-x_{CM})/d \\ \left(y-R \sin\left(\dfrac{2 {\rm \pi}t}{T_{R}}\right)\right)/d \\ (z-z_{CM})/d \end{matrix}\right), \end{equation}

where $R$ is the radius of rotation of the time-averaged centre-of-mass trajectory $y_{CM}(x_{CM})$. Therefore, $y_{CM}$ oscillates around $R \sin ({2 {\rm \pi}t}/{T_{R}})$.

The snapshots in figure 3(e) are strikingly reminiscent of the sequence of snapshots from the crawling mode, shown in figure 3(d). In contrast to the crawling mode, for the rotation-crawling mode, the shapes at $t+T/2$ differ from the mirror reflection of shapes at time $t$.

The centre-of-mass trajectory in this mode shows the superposition of three motions: settling down, rotation around gravity and periodic sideways motion (see figure 4c). In the reference frame settling down with velocity $\dot {z}_{CM}$, two independent characteristic time scales describe this motion: a longer period $T_{R}$ that is associated with the rotation and a shorter period $T$ that characterises the smaller variations. For the example shown by red solid lines in figure 4, $T_R\approx 3690 \tau _{b}$ and $T\approx 390 \tau _{b}$. Both time scales are also seen in figure 4(d) as the azimuthal angle $\phi$ increases with the superimposed small periodic oscillations. The difference $\Delta z$ between the vertical positions of the fibre ends is also a periodic function with the period $T$, although its positive and negative amplitudes have different magnitudes, as opposed to the crawling mode (see figure 4b).

4.6. Irregular mode

Up to this point, we have discussed regular modes, which are attained after a sufficiently long time. However, for certain values of the parameters, the attracting modes have not been reached, even for a very long time. Such a behaviour will be called an irregular mode. An example of such an irregular evolution is shown in figure 4(e,f), from $t=0$ until $t=27\,800\tau _{b}$. The azimuthal angle $\phi (t)$ is plotted in figure 4(e). In the same plot, for comparison, we added two examples of the numerical trials, which, after a certain time, reach an attracting regular mode and later stay in this mode for a long time.

For all the numerical trials, the initial almost straight and horizontal configuration is followed by a short-time U-shape deformation of the fibre almost in a vertical plane (as for more stiff fibres Schlagberger & Netz 2005; Llopis et al. 2007). However, a vertical U-shape configuration is not stable for highly flexible fibres and, then, transient or irregular dynamics of the out-of-plane shapes are observed. For a range of large (but not too large) values of $B$, the dynamics after a relatively short time is transformed into a regular mode. However, for some of the larger values of $B$, the fibre moves non-systematically all the time. An example of the irregularity is visible in figure 4(e,f). The irregularity of the centre-of-mass trajectory shown in figure 4( f) contains certain parts with an irregular rotation and irregular oscillations, resembling the regular rotation-crawling mode on a very short time scale. The irregularity is also discussed in the next section and in Appendix C.

5. Time scales, asymmetry and bifurcations of the modes

In this section, we study in detail the properties of the modes. We first focus on the regular modes and determine the dependence on $B$ and $N$ of two characteristic periods: the rotation period, $T_R$, and the oscillation period, $T$.

The rotation period $T_{R}$ is the longest characteristic time that is present for rotation and rotation-crawling modes. $T_{R}$ as a function of $B$ and $N$ is shown in colours in figure 5(a). For both modes, $T_{R}$ varies in the comparable range of values. There is a clear dependence on $N$ and $B$ observed in the case of the rotation mode: $T_{R}$ decreases with an increase of $B$ for a fixed $N$, and it decreases also with growing $N$ for a fixed $B$. More flexible fibres form stationary shapes that rotate faster.

Figure 5. Characteristic features of the modes. (a) The rotation period $T_{R}$ for the rotation and rotation-crawling modes. (b) The oscillation period $T$ for the crawling and rotation-crawling modes. (c) Maximum difference between vertical positions of the fibre ends ${| \Delta z |}_{max}$ during the mode. (d) Average time $T_{100d}$ for a fibre to sediment by $z_{CM}=100d$. In (a) and (b), the white colour corresponds to the absence of rotation or oscillation, and the hatched area to the irregular mode; in (c), the white colour corresponds to $\Delta z = 0$ (exactly) for the tilted mode. In (a) $T_{R}=93\,280 \tau _{b}$ for $N=34$ and $B=5500$, and $T_{R}=20\,800 \tau _{b}$ for $N=22$ and $B=9000$.

For the rotation-crawling mode, the dependence of $T_{R}$ on $N$ and $B$ is not so simple. In most cases, a general tendency is observed, but there are exceptions. When $N=26$ is considered, $T_{R}$ for the rotation-crawling mode decreases with an increase of $B$, similar to the behaviour of $T_{R}$ for the rotation mode. This is also observed for $N=38$, though on two data points only. However, exceptionally, for $N=34$, the value of $T_{R}$ for $B=10\,000$ is much larger than for $B=8500$. In general, for the rotation-crawling mode, the dependence of $T_{R}$ on $N$ for a fixed $B$ is as for $B=8500$: a gradual increase of $T_{R}$ with increasing $N$ is observed, in contrast to the rotation mode. This difference is related to the time-dependent shape and time-dependent angular velocity of the fibre in the rotation-crawling mode. However, for $B=9000$, there is an exception from this general tendency: the value $T_{R} = 20\,800 \tau _{b}$ for $N=22$ is much larger than $T_{R}$ for a slightly larger $N=24$.

The characteristic oscillation period $T$ for the crawling behaviour as a function of $N$ and $B$ is shown in colours in figure 5(b). It is seen that, in general, $T$ increases with increasing $N$ for a fixed $B$, and $T$ decreases with increasing $B$ for a fixed $N$. But again, there are exceptions. Three regions stand out from the general tendency described previously. Two of them correspond to the exceptional behaviour of $T_R$. Namely, the value of $T$ for $N=22$ and $B=9000$ is much larger than the values of its neighbours. The values of $T$ for $N=32$, $34$ and $B=9500$, $10\,000$ are also larger than $T$ for smaller values of $B$. The third region is spanned by $30 \leq N \leq 34$ and $4500 \leq B \leq 6000$. In this range, extremely large values of $T$ are observed for the largest value $N=34$, close to the border with the rotation mode. It seems that the exceptions from the general tendency occur for values of $B$ and $N$ close to transitions between different modes.

It is useful to compare $T_R$ and $T$ to a characteristic time of the gravitational settling, chosen as the time $T_{100d}$ needed for a fibre to move along gravity the distance equal to $100d$. The time $T_{100d}$ is computed separately for each pair $B$ and $N$, based on the mean vertical centre-of-mass velocity in the respective attracting mode. In figure 5(d), we show $T_{100d}$ as a function of $N$ and $B$ with the use of colours.

It is clear that $T_{100d}$ is much smaller than the rotation period $T_{R}$ and also much smaller than the oscillation period $T$. The settling time $T_{100d}$ only weakly depends on $B$, but it systematically decreases with the increase of $N$. The increase of $N$ results in the increase of the fibre mass, but also in the increase of its drag. The competing nature of these two effects leads to a small increase of the fibre settling velocity and, therefore, a small decrease of $T_{100d}$ (for larger values of $N$, the settling velocity scales as $\ln N +C$ Adamczyk et al. 2010). Therefore, even though the fibre's mass, as well as the fibre's flexibility, change by up to a factor of three or more, figure 5(d) shows that the sedimentation time $T_{100d}$ changes by no more than 1.4 times. Recall that for a single bead, $T_{100d}=300 \tau _{b}$.

For the tilted mode, there is a clear dependence of $T_{100d}$ on $N$ and $B$: the sedimentation time $T_{100d}$ increases with increasing $B$ for a fixed $N$, but for a fixed $B$ it decreases with increasing $N$ as the fibre becomes heavier and heavier. In the case of the rotation mode, a similar decreasing tendency is observed with increasing $N$ for a fixed $B$. For this mode, on the other hand, $T_{100d}$ practically does not change with increasing $B$ for a fixed $N$. For the other sedimentation modes, the fibre shape changes over time, which significantly affects the mean sedimentation velocity and also $T_{100d}$.

More information about the time-dependent sedimentation velocity is provided in Appendix B. In figure 9, we plot the maximum sedimentation velocity and the difference between the maximum and the minimum sedimentation velocities for each mode with given values of $B$ and $N$.

The regular and irregular modes can be compared with each other by ‘an asymmetry parameter’: the maximum difference ${| \Delta z |}_{max}$ between the vertical positions of the fibre ends achieved during the observation of the mode with given values of $B$ and $N$. In figure 5(c), ${| \Delta z |}_{max}$ is normalised by $L_0$: the fibre length in the elastic equilibrium. In general, the increase in the fibre flexibility, i.e. the increase of $B$ or $N$, leads to an increase in ${| \Delta z |}_{max}/L_0$. However, due to the symmetry of the tilted mode, $\Delta z = 0$ for any values of $N$ and $B$. The rotation mode exhibits an increase in ${| \Delta z |}_{max}/L_0$ with increasing $N$ (or $B$) for a fixed $B$ (or $N$). For the crawling and rotation-crawling modes ${| \Delta z |}_{max}/L_0$ behaves similarly to the rotation mode, except the three regions of $N$ and $B$ which exhibit exceptional behaviour of $T_R$ or $T$, as described earlier.

To emphasise the complexity of the motion and shape deformation of different modes, and to highlight their differences, figure 6 shows the projection of the trajectories of the tip of the end-to-end vector $\boldsymbol {r}_{1N}$ on the plane of the cylindrical coordinates $\rho$ and $\Delta z$, for $N=36$ and $3000 \le B \le 8500$, changed by the step of 500. For the tilted mode, $| \Delta z | = 0$ and, therefore, the radial coordinate $\rho$ is constant in time, and the same is true for the length of the end-to-end vector $|\boldsymbol {r}_{1N} |$. As the shape of the fibre does not change in the rotation mode, the radial coordinate $\rho$ of the end-to-end vector does not change in time. Symmetry with respect to the reflection $\Delta z \rightarrow -\Delta z$ is present in the crawling mode when the $\rho -\Delta z$ projection of the tip of the vector follows a figure of eight: $\Delta z$ is a symmetric periodic function with a period $T$ (see also figure 4b). Finally, a distorted figure of eight is drawn on the $\rho -\Delta z$ plane by the tip of the end-to-end vector in the rotation-crawling mode, as $\Delta z$ is now a non-symmetric periodic function with a period $T$ (see also figure 4b).

Figure 6. Projection of the trajectories of the tip of the end-to-end vector $\boldsymbol {r}_{1N}$ on the plane of the cylindrical coordinates $\rho$ and $\Delta z$ for $N=36$ and $3000 \le B \le 8500$, changed by the step of $B_{step} = 500$. The colours of the curves in the online version correspond to the same colours of the modes as in figures 1 and 4. The values of $B$ are stated for the tilted ($B=3000,\ 3500$), crawling ($B=7500,\ 8000$) and rotation-crawling ($B=8500$) modes. The values $4000\le B \le 7000$ for the rotation mode correspond to the dots (orange online) under an arrow, which marks the increase of $B$.

Overall, the data analysis suggests that sedimentation modes with a non-stationary fibre shape cannot be easily described. There are several indications that this is the case. First, there exists the region $30 \leq N \leq 34$ and $4500 \leq B \leq 6000$, for which the behaviour of $T$ and ${| \Delta z |}_{max}/L_0$ differs from the behaviour for the other values of $N$ and $B$. Second, there are two other regions (with $N=22$ and $B=8500$, and $N=32$, $34$ and $B=9500$, $10\,000$), which break the monotonic dependence on $N$ and/or $B$ for $T_{R}$, $T$ and ${|\Delta z |}_{max}/L_0$. Finally, there exists the irregular mode for different, separated from each other regions of $N$ and $B$.

Therefore, in figure 7 we analyse the differences between the dynamics for $B=5500$ and different values of $N$. We show the projection of the trajectory of the tip of the end-to-end vector on the plane of $\rho$ and $\Delta z$ for $N=28,30$ (crawling mode), $N=32$ (irregular mode) and $N=34$ (rotation-crawling mode), close to the transition to the rotation mode, observed at $N=36$. Between $N=28$ and $N=30$, a period-doubling bifurcation takes place. A sequence of bifurcations between $N=30$ and $N=32$ might lead to a strange attractor and chaotic behaviour of the irregular mode. A sequence of period-halving bifurcations between $N=32$ and $N=36$ may lead to the formation of the rotation mode.

Figure 7. Bifurcations of the dynamics with the increasing aspect ratio $N$ (top-down: $N=28, 30, 32, 34$) for a fixed value $B=5500$. In each frame, points or curves for $N=36$ (rotation mode) are drawn in orange for comparison. (Left column) Cylindrical coordinate $\rho$ as a function of $\Delta z$, (centre) vertical coordinate $\Delta z$ as a function of time $t$, and (right) polar angle $\theta$ as a function of time $t$. For $N=32$, the parts of the curves corresponding to $20\,000 \le t/\tau _{b}\le 23\,000$ are highlighted in red on all three plots.

Figure 7 illustrates that the existence of the irregular mode is an inherent feature of the dynamics of highly elastic fibres. Transitions between the crawling, rotation-crawling and irregular modes, seem to be sensitive to even small changes in values of the parameters $B$ and $N$, as illustrated in Appendix C. To determine the nature of the bifurcations in the phase diagram, the Floquet stability analysis would be useful, but it goes beyond the scope of this paper.

6. Conclusions

In this work, we have numerically determined the dynamics of a single highly elastic fibre settling under gravity in a very viscous fluid, at the Reynolds number much smaller than unity. We focused on long-time behaviour and showed that there exist different attracting modes. We provided in figure 1 the phase diagram of these modes, obtained by systematically varying the elasto-gravitation number $B$ and the fibre aspect ratio $N$ in certain ranges. In this way, we demonstrated that for highly elastic fibres, not only $B$ but also $N$ significantly affects the fibre shapes and their dynamics. We confirmed the existence of two families of non-planar stationary configurations, shown by Saggiorato et al. (2015), one translating and one rotating, and we investigated their properties. We found two new families of the dynamical modes with periodic deformations of the fibre shape, one translating and one rotating. Finally, we showed that very flexible fibres often deform and move irregularly. This irregularity seems to be related to bifurcations of the dynamics relatively close to transitions between the crawling, rotation-crawling and irregular modes, with a possibility of chaotic dynamics. The transition between the vertical, tilted and crawling modes is more simple and regular.

The regular modes presented in this work are achieved after a very long time, as illustrated in figures 4(e) and 5(d). However, in our simulations, we often observe transient, short-time effects similar to the long-time attracting modes described previously. Examples of this similarity are visible by comparing our long-time modes with some features of the short-time dynamics observed in the numerical and experimental data reported by Shashank et al. (2023).

It is worth emphasising that the complex time-dependent shapes of highly elastic fibres and their motion found in this work are three-dimensional. Initially, we impose a small perturbation out of a vertical plane. Within our theoretical description of a highly elastic fibre in a three-dimensional unbounded fluid, we have not found any stationary or periodic solutions restricted to a vertical plane. Therefore, we have not recovered the non-symmetric periodic motions of a highly flexible filament within a vertical plane, found for a two-dimensional fluid and reported by Shojaei & Dehghani (2015). The absence of such solutions in our simulations may be related to a much smaller aspect ratio and, therefore, a much smaller value of $B$. Alternatively, it may be caused by the use of a more precise theoretical approach. Namely, we take into account the thickness of the fibre and the hydrodynamic interactions between all its segments.

The complex dynamics of highly flexible fibres found in this work differ significantly from the dynamics of moderately or weakly flexible fibres. It suggests the possibility of a qualitative change in the structure and transport properties of flexible fibre suspensions when their bending stiffness is decreased.