Flexible fibers in shear flow approach attracting periodic solutions

TL;DR

This study numerically investigates the complex three-dimensional behaviors of flexible fibers in shear flow, revealing multiple attracting states including steady configurations and periodic motions, with transitions influenced by bending to hydrodynamic force ratios.

Contribution

It identifies and characterizes new attracting periodic solutions and steady configurations of flexible fibers in shear flow, expanding understanding beyond extensional flow dynamics.

Findings

Existence of attracting steady and periodic motions.

Effective Jeffery orbits with exponential parameter drift.

Identification of a three-dimensional meandering periodic mode.

Abstract

The three-dimensional dynamics of a single non-Brownian flexible fiber in shear flow is evaluated numerically, in the absence of inertia. A wide range of ratios A of bending to hydrodynamic forces and hundreds of initial configurations are considered. We demonstrate that flexible fibers in shear flow exhibit much more complicated evolution patterns than in the case of extensional flow, where transitions to higher-order modes of characteristic shapes are observed when A exceeds consecutive threshold values. In shear flow, we identify the existence of an attracting steady configuration and different attracting periodic motions that are approached by long-lasting rolling, tumbling and meandering dynamical modes, respectively. We demonstrate that the final stages of the first and second modes are effective Jeffery orbits, with the constant parameter C replaced by an exponential function that either decays or increases in time, respectively, corresponding to a systematic drift of the trajectories. In the limit of C going to zero, the fiber aligns with the vorticity direction and in the limit of C going to infinity, the fiber periodically tumbles within the shear plane. For moderate values of A, a three-dimensional meandering periodic motion exists, which corresponds to intermediate values of C. Transient, close to periodic oscillations are also detected in the first stages of the modes.

Tables

Table 1. Table 1: Long-time dynamics of fibers for different l 0 subscript 𝑙 0 l_{0} & A 𝐴 A .

$l_0\backslash A$	9	10	11
1.005	rolling	close to periodic meandering	periodic   meandering
1.02	close to periodic meandering	periodic meandering	close to periodic meandering

Equations

$${\mathbf{v}}_0=\dot{\gamma }z{\mathbf{e}}_x,$$

$$A=E/(64\eta \dot{\gamma }),$$

$$\bar{\mu }=8(L/d{)}^4/A,$$

$$E_s=\frac{\hat{k}}{2}\sum _i({\ell }_i-{\ell }_0{)}^2.$$

$$E_b=\frac{\hat{A}}{2{\ell }_0}\sum _i({\hat{\mathbf{t}}}_{i+1}-{\hat{\mathbf{t}}}_i{)}^2,$$

$$\hat{A}=E\pi d^4/64.$$

$${\mathbf{F}}_i=-\frac{\partial }{\partial {\mathbf{r}}_i}(E_s+E_b).$$

$$\eta {\nabla }^2\mathbf{v}-\nabla p=0,$$

$$\displaystyle\left(\begin{array}[]{cc}\bm{U}\\ \bm{\Omega}\end{array}\right)$$

$$T=2\pi (r+1/r),$$

$$\tan {\Theta }_J(t)$$

$$\tan \Phi (t)$$

$${\cos }^2\Phi (t)=\frac{r^2\sin \left(\frac{rt}{1+r^2}\right)}{r^2+(1-r^2)\cos \left(\frac{rt}{1+r^2}\right)}.$$

$$&T\equiv T_{r}=459\hskip 27.03003pt\mbox{for the rolling mode,\hskip 14.22636pt}\\ T\equiv T_{t}=181\hskip 28.45274pt\mbox{for the tumbling mode,\hskip 14.22636pt}$$

$$&\tau\equiv\tau_{r}=-1253\hskip 14.22636pt\mbox{for the rolling mode},\hskip 14.22636pt\\ \tau\equiv\tau_{t}=170\hskip 28.45274pt\mbox{for the tumbling mode,\hskip 14.22636pt}$$

$$\tan \Theta (t)=\exp (t/\tau )\tan {\Theta }_J(t),$$

$$C(t)=C_0\exp (t/\tau ),$$

Full text

I Introduction

In nature and modern technologies, there are many systems containing elongated, flexible, micrometer- and nanometer-scale objects deforming and moving in a fluid flow duRoure . Examples are micro-swimmers such as bacteria and their flagella, actin Kantsler ; Harasim2013 ; Liu , large proteins or DNA molecules Chu , and micro- or nano-fibers Lindner_2010 ; Nunes_2012 ; lindner_2015 ; perazzo ; pawlowska . Flexibility leads to complex dynamics with buckling Hinch1976 , coil-stretch transitions Kantsler ; deGennes ; BeckerShelley ; YoungShelley , migration across the streamlines of a flow graham2004 ; graham2005 ; graham2006 ; Slowicka2012 ; Slowicka_2013 ; farutin ; Misbah2019 , knotting yeomans ; stone_2015 ; doyle ; gruziel and a variety of deformed shapes Arlov ; Skjetne ; Joung ; tornberg ; fauci1 ; fauci ; wang2013three ; zhang2019dynamics ; rost2019effective . Similar complexity of dynamics has been observed for flexible filaments in electrokinetic fields Doyle2016 or sedimenting under gravity Lagomarsino ; Netz ; Llopis ; Li ; Gompper ; BukowickiGruca ; Bukowicki ; Duprat2018 ; Bukowicki2 ; Gruziel2 .

A basic question in all of these configurations is how the dynamics and shapes of deformable elongated objects in flow depend on their flexibility. This problem has been investigated extensively at macro- and nano-scales for fibers in extensional, cellular and corner flows Kantsler ; YoungShelley ; lindner_2015 ; Autrusson . For example, it has been shown that the typical pattern of a fiber’s evolution in extensional flow is related to consecutive threshold values of the characteristic ratio of bending to hydrodynamic forces exerted by the fluid flow. When these values are exceeded, higher-order modes of the fiber shape are activated, with shorter characteristic length scales of elastic deformation Kantsler ; Harasim2013 ; Chu ; BeckerShelley ; YoungShelley ; Liu . In extensional flow, there exists a family of characteristic, well-defined shapes.

However, in general, the fiber deformation that occurs may depend on the type of flow Chu ; deGennes ; ICTAM . Shear flows typically occur owing to the presence of container walls, and so are of practical and fundamental interest. Therefore, as illustrated in Fig. 1, in this paper we study three-dimensional dynamics of a single flexible fiber immersed in steady shear flow with velocity

[TABLE]

where ${\bf e}_{x}$ is the unit vector along $x$ and $\dot{\gamma}$ is the shear rate. The Reynolds number of the system is assumed to be much smaller than unity and the fluid flow satisfies the quasi-steady Stokes equations.

In this article we provide a new perspective on the three-dimensional evolution of flexible fiber shapes in shear flow. We show that after a relaxation phase, a flexible fiber is attracted to one of several stationary, periodic or close to periodic solutions, with different typical sequences of shapes and orientations. Features of these characteristic solutions, their presence or absence, stability or instability, and basins of attraction depend on the ratio $A$ of local bending $E\pi d^{2}/64$ to hydrodynamic $\pi\eta\dot{\gamma}d^{2}$ forces, where $E$ is the Young’s modulus, $d$ is the fiber diameter and $\eta$ is the fluid’s dynamic viscosity. Therefore, by changing

[TABLE]

different structures and patterns of motion of dilute suspensions of flexible fibers can be obtained. For a fiber of length $L$ and diameter$d$ of the cross-section, the local bending-to-hydrodynamic ratio $A$ is related to the elasto-viscous number

[TABLE]

which depends on the fiber’s aspect ratio $L/d$ Liu .

A rigid straight fiber in a low-Reynolds-number shear flow follows a periodic Jeffery orbit Jeffery (th In this article we demonstrate that for a wide range of parameters and times, the motion and deformation of a flexible fiber can be interpreted as an effective Jeffery orbit that systematically drifts in time, owing to the exponential growth or decrease of the time-dependent parameter $C$. For moderate values of $A\approx$ 9-12, there exists a range of values of $C$ that correspond to periodic or close to periodic meandering motions.

A significant finding described below is that for a certain range of small values of $A\approx$ 4-5, fibers typically tend to align with the vorticity of the undisturbed flow. This result indicates a possibility to observe experimentally an ordered, dilute suspension of flexible fibers where all fibers are straight and aligned with the vorticity direction. We suggest that this ordered phase could be produced by adjusting the shear rate to reach the appropriate range of $A$. We are not aware of any such experimental observations.

As we document in this article, in shear flow (in contrast to extensional flow), the value of the bending-to-hydrodynamic ratio $A$ does not uniquely determine the type of the fiber shape. This conclusion is based on two main features of the dynamics. First, for the same value of $A$, depending on the initial configuration or orientation, an elastic fiber evolves to a different characteristic sequence of shapes. Moreover, long-lasting, chaotic transients are typical: we document close to periodic motions that later spontaneously change into periodic or effective Jeffery motions with different shape sequences.

II Theoretical model of a fiber in flow

II.1 Elastic fiber

The fiber is modeled Stark as a chain of $N$=40 spherical beads of diameter $d$. The time-dependent position of the center of a bead $i$ is denoted as $\bm{r}_{i}$. The centers of the consecutive beads are connected by springs of the equilibrium length $\ell_{0}\!=\!1.02d$ and the potential energy

[TABLE]

Here $\hat{k}$ is the spring constant and $\ell_{i}=|\bm{r}_{i}-\bm{r}_{i-1}|$ is the distance between centers of beads $i$ and $i\!-\!1$. In this paper we assume that the dimensionless elastic resistance ratio ${k}\!=\!{\hat{k}}/{(\pi\eta d\dot{\gamma})}$ is large, $k\!=\!1000$, which leads to an almost constant fiber length.

At equilibrium, the fiber is straight; its deformation costs energy

[TABLE]

dependent on the bending resistance

[TABLE]

Here $\hat{\bm{t}}_{i}$ is the unit vector parallel to the relative positions $\bm{r}_{i}-\bm{r}_{i-1}$ of the centers of beads $i$ and $i\!\!-\!\!1$.

The total external force $\bm{F}_{i}$ acting on bead $i$ is elastic,

[TABLE]

We assume that the dimensionless bending-to-hydrodynamic force ratio (relative bending stiffness) $A\!=\!{\hat{A}}/{(\pi\eta d^{4}\dot{\gamma})}$, given by Eq. (2), varies in the range of moderate values $4\!\leq\!A\!\leq\!40$, where most of the fibers subsequently deform and straighten while tumbling Slowicka_chaos_2015 . The length and time units in the simulations are, respectively, $d$ and $1/\dot{\gamma}$.

II.2 Initial fiber configurations

To study the fiber evolution, we analyze the time dependence of the end-to-end vector ${\bf n}(t)\!=\!(\delta x,\,\delta y,\,\delta z)$, shown in Fig. 1(b). We parameterize it by the standard spherical coordinates: the vector length $\Delta L(t)=|{\bf n}(t)|$, the angle $\Theta(t)$ between ${\bf n}(t)$ and the vorticity direction $y$, and the angle $\Phi(t)$ between the projection of ${\bf n}(t)$ on the $xz$ plane and the $x$ axis. Initially, ${\bf n}(0)\!=\!{\bf n}_{0}$, $\Theta(0)\!=\!\Theta_{0}$ and $\Phi(0)\!=\!\Phi_{0}$.

To investigate the characteristic features of the flexible fiber dynamics in shear flow, we consider the following family of initial conditions. We assume that at $t=0$ the fiber is straight and that all of the springs between the beads are at their equilibrium lengths ICTAM . The initial orientation of the fiber is given by the orientation of the end-to-end vector ${\bf n}_{0}$ that links the centers of the first and the last beads. The length of this vector is $|{\bf n}_{0}|\!=\!L_{0}\equiv(N-1)\ell_{0}$. The direction of ${\bf n}_{0}$ is parameterized by the spherical angles $\Theta_{0}$ and $\Phi_{0}$, as indicated in Fig. 1(a). We consider the whole range of the initial orientations.

In the following, we will systematically investigate how the fiber’s dynamics and shape evolution depend on the relative bending stiffness $A$ and the initial orientation $(\Theta_{0},\Phi_{0})$. In this way we will study the three-dimensional dynamics of a flexible fiber in shear flow, while most of the previous studies have focused on the two-dimensional dynamics in the shear plane.

II.3 Fiber dynamics in flow

The dynamics of an elastic fiber is determined by the external (i.e., elastic) forces exerted on each fiber bead, and the hydrodynamic interactions between them caused by the presence of the shear flow. We assume that the dynamic viscosity $\eta$ of the fluid that surrounds the fiber is large enough for the Reynolds number to be much less than unity. In this limit, the fluid velocity $\bm{v}$ and pressure p satisfy the Stokes equations batchelor1967 ; kim2013microhydrodynamics ,

[TABLE]

The no-slip boundary conditions at the beads surfaces are assumed. The fluid is unbounded, with the ambient flow velocity ${\bf v}_{0}$ given by Eq. (1); in the presence of the fiber, the fluid velocity $\bm{v}$ tends to ${\bf v}_{0}$ when the distance from the fiber goes to infinity.

To solve the Stokes equations in the presence of the fiber, we use the advanced theoretical algorithm cichocki1994friction ; ekiel2009precise , based on the multipole expansion corrected for lubrication cichocki , implemented in numerical codes Hydromultipole (cichocki ). The method, similar to Brady ; Ladd , is based on the boundary integral representation of the fluid velocity and the boundary integral equations for the surface density of the forces induced at the particle surfaces. These equations are projected on a complete set of elementary solutions of the Stokes equations (spherical multipole functions). The resulting set of linear algebraic equations is truncated at a controlled multipole order.The convergence of the multipole expansion is speed up by applying the lubrication correction cichocki ; Bossis ; Sangani .

In general, dynamics of $N$ spherical beads moving in an ambient Stokes flow have the form

[TABLE]

where $\bm{U}\!\!=\!\!(\bm{U}_{1},...,\bm{U}_{N})$ and $\bf\Omega\!\!=\!\!(\bf\Omega_{1},...,\bf\Omega_{N})$ are translational and rotational bead velocities, $\bm{F}\!\!=\!\!(\bm{F}_{1},...,\bm{F}_{N})$ and $\bm{{\cal T}}\!\!=\!\!(\bm{{\cal T}}_{1},...,\bm{{\cal T}}_{N})$ are external forces and torques (couples) exerted on the particles; $\bm{F}_{0}\!\!=\!\!(\bm{F}_{0,1}...,\bm{F}_{0,N})$ and $\bm{{\cal T}}_{0}\!\!=\!\!(\bm{{\cal T}}_{0,1},...,\bm{{\cal T}}_{0,N})$ are forces and torques exerted by the ambient flow on the motionless beads. Here, $\bm{\mu}_{ij}^{tt},\;\bm{\mu}_{ij}^{tr},\;\bm{\mu}_{ij}^{rt},\;\bm{\mu}_{ij}^{rr}$ are the 3$\times$3 translational-translational, translational-rotational, rotational-translational and rotational-rotational mobility matrices, respectively. They depend on positions $\bm{r}_{i}$ of the centers of all the beads $i=1,...,N$.

In this paper, Eq. (15) simplifies, because there are no external torques, $\bm{{\cal T}}\!\!=\!\!{\bf 0}$, and we are not interested in the bead rotations $\bm{\Omega}$. The Hydromultipole codes cichocki are used to evaluate with high precision $\bm{\mu}_{ij}^{tt},\;\bm{\mu}_{ij}^{tr},\;\bm{\mu}_{ij}^{rt},\;\bm{F}_{0},\;\bm{{\cal T}}_{0}$ for given positions $\bm{r}_{i}$ of the bead centers, and Eqs. (15) are solved numerically with the adaptive fourth-order Runge-Kutta method. More information about the numerical method and its accuracy is given in Appendix C.

III Attracting modes of the dynamics

One of our most significant findings is that, depending on the initial orientation and bending stiffness, the fiber is attracted to one of three distinct stationary or periodic solutions. The characteristic properties of these evolution patterns are illustrated in Fig. 2 and Movie 1, for $A\!\!=\!\!10$ and three different initial orientations of the end-to-end vector ${\bf n}_{0}$, $(\Theta_{0},\Phi_{0})\!\!=\!\!(5^{\circ},10^{\circ})$, $(45^{\circ},30^{\circ})$, $(10^{\circ},10^{\circ})$. The colors blue, red and green indicate time-dependent modes attracted to different stationary or periodic solutions. In the following these attractors will be called “rolling”, “tumbling” and “periodic meandering”, respectively.

The rolling mode (also called spin-rotation Skjetne , or log rolling) is shown at early times in Fig. 2(a) and Movie 1. In this mode, the fiber end-to-end vector tends at long times to the vorticity ($y$) direction, i.e. $\Theta\!\rightarrow\!0$. Moreover, all of the fiber beads rotate around $y$, and tend to align along $y$. The tumbling mode is shown at early times in Fig. 2(b) and Movie 1. In this mode, the fiber end-to-end vector tends at long times to the $xz$ plane ($\Theta\!\rightarrow\!90^{\circ}$). The convergence is confirmed by much longer simulations with $0\!\leq\!t\!\leq\!15000$.

The attracting character of the rolling and tumbling solutions have been identified previously Arlov ; Forgacs ; Skjetne ; Joung ; Wang_Yu ; ICTAM as the generic feature of the flexible fiber dynamics, and the corresponding dynamical modes have been called by the same names as their attractors Skjetne ; Joung ; stone_2015 .

To our surprise, we discovered that there exists a third dynamical evolution, illustrated in Fig. 2(c) and Movie 1 – an approach to an attracting three-dimensional “periodic meandering” solution, with $\Theta$ oscillating periodically for $1800\!\lesssim\!t\!\leq\!10^{5}$. To the best of our knowledge, this behavior has not been observed before.

The results in Fig. 2 are focused on the relaxation phase after which the attractors are reached. The relaxation time can be very long, of the order of thousands of dimensionless units, as can be seen in Fig. 2. Moreover, long-lasting, close to periodic, transient motions (different from tumbling and meandering), which we will call “squirming”, can appear in the relaxation phase, as shown in Fig. 2(c) and Movie 1. Squirming motions are typically present for a wide range of initial orientations and values of $A$, and in all of the modes - those that tend to meandering, rolling or tumbling solutions.

The approach to three different attracting solutions is also well illustrated by plotting the trajectories drawn by the tip of the end-to-end vector ${\bf n}(t)\!=\!(\delta x,\,\delta y,\,\delta z)$. This concept is explained in Movie 3. The results obtained for $A$=10, and the same initial orientations as in Fig. 2, are shown in Fig. 3. We use different scales for different modes to illustrate the detailed characteristic features of each evolving mode, which is attracted by a different solution: stationary rolling, periodic tumbling and periodic meandering, respectively. Approaching these solutions takes a long time. We highlight in Fig. 3(a)(ii) the regular, anisotropic spiraling towards the steady state. We demonstrate in Fig. 3(b) a fast convergence of ${\bf n}$ to the shear plane. We illustrate in Fig. 3(c) the existence of a transient, almost periodic squirming trajectory in the relaxation phase of the meandering mode.

The transition between the different modes is illustrated in Fig. 4. The columns (a)-(c) show trajectories of the fiber end-to-end vector for $A$=10 and three different initial orientations, which are very close to each other, but represent different modes; we keep the same scale in each row of the figure. It is clear that the dynamics can be sensitive to a small change of the initial conditions.

We next analyze the trajectories in Fig. 4(b,c) during the relaxation phase. By comparing them, we observe that at short and moderate times there appear a range of times with the characteristic squirming trajectory, which repeats almost periodically. This squirming end-to-end trajectory is shown separately in Fig. 5.

Summarizing, the transient periodic squirming motion is visible not only in the meandering mode (which ends up at the periodic meandering solution), but also in the tumbling mode (which ends up at the tumbling solution). For smaller values of the relative bending stiffness $A$, squirming motion is also often visible before the rolling solution is reached. The attractors – the rolling, tumbling and meandering solutions – will be shown in the next section.

IV The attractors: stationary and periodic solutions

The three modes presented in Figs. 2-4 tend to three different attracting solutions that are compared with each other in Figs. 6-7. As illustrated in Fig. 6, the ordering of the end-to-end vector is different for each of these attractors, with the corresponding mean values $\langle\delta y/L_{0}\rangle$ well-separated from each other.

In the rolling solution, the fiber is straight and it spins along the vorticity direction $y$. In the tumbling periodic solution (with the period $T_{t}\!=\!181$ for $A$=10), the fiber end-to-end vector is perpendicular to $y$. It tumbles in the shear plane around the vorticity direction while the fiber straightens along the flow and then becomes coiled, with large-amplitude oscillations of length of its end-to-end vector.

The periodic meandering orbit (with the period $T_{m}\!=\!735$ for $A$=10) has been observed even in very long simulations with more than 130 periods. Details about the periodic meandering motion are given in Appendix A. The evolution of shapes in the meandering and tumbling periodic solutions are compared with each other in Fig. 7(a,b) and Movie 2. In both periodic orbits, the $xz$ projections of the fibers have small amplitudes $\delta z$ along the flow gradient direction, but the corresponding shapes are essentially different from each other. Moreover, in the tumbling motion, all of the fiber beads stay in the shear plane (i.e., $\delta y$=0), while the meandering motion is three-dimensional, with large values of the end-to-end projection $\delta y$ along the vorticity direction. Therefore, the meandering and tumbling periodic motions differ significantly from each other.

V Dependence of the dynamics on the bending stiffness ratio $A$

In this paper, we have systematically investigated how the evolution pattern depends on the fiber’s initial orientation $(\Phi_{0},\Theta_{0})$ for different values of the bending stiffness ratio $A$. A summary of the most important findings is illustrated in Fig. 8(a-c) for $4\leq A\leq 40$. Flexible fibers with $A\!=\!4$, shown in Fig. 8(a), belong to the rolling mode (blue dots) for all the initial orientations other than in the shear plane $xz$ ($\Theta_{0}\!=\!90^{\circ}$). For rather stiff fibers with $A\!=\!40$, shown in Fig. 8(c), the tumbling mode (red dots) dominates for most of the initial orientations, except some of those that are very close to the vorticity direction and lead to the rolling mode (blue dots).

For fibers of a moderate stiffness $A\!=\!10$, shown in Fig. 8(b), there exists a range of the initial orientations that correspond to the periodic mode, with the same attracting periodic meandering orbit. This range, marked in Fig. 8 by green dots, is approximately contained between two black solid lines. Their meaning will be explained in the next section by comparing with the Jeffery orbits Jeffery . Close to the boundaries that separate different dynamical modes in Fig. 8, some irregularities appear and a relatively small change of initial orientations can result in a different dynamical mode. An example of such behavior for three close initial orientations was shown in Fig. 4. This sensitivity is related to chaotic properties of the transition between the modes.

We have analyzed the fiber dynamics for a wide range of values of the bending stiffness $A$. The periodic meandering solutions are observed for a certain range of $A$. In addition, for values of $A$ not far from this range, close-to-periodic meandering motions are observed. Examples of such motions are shown and discussed in Appendix B.

In Fig. 8(d), the dependence of the attracting solutions on the bending stiffness $A$ is shown for the initial orientations with $\Phi_{0}=10^{\circ}$ and different values of $\Theta_{0}$. The attracting periodic or close to periodic meandering motions (green dots) are visible for $9\leq A\leq 12$. Examples of close to periodic meandering solutions are shown in Appendix B. We observe that for $A=7.3$, some of the initial orientations are attracted to another three-dimensional periodic orbit, corresponding to the squirming solution (marked by green diamonds). This finding illustrates that the squirming periodic motion can be transient or attractive, depending on the value of the bending stiffness ratio $A$.

The results shown in Fig. 8 indicate that there exist well-defined ranges of the fiber bending stiffness in which different solutions dominate. This result is significant for dilute suspensions of flexible fibers where hydrodynamic interactions between fibers are weak. In particular, our results (Fig. 8a) predict the existence of an ordered phase with all (or almost all) the fibers straightened out and parallel to each other (and to the vorticity direction), for a narrow range of values of the ratio A of bending to hydrodynamic forces ($4\leq A<6$ for $N=40$). This range can be reached (or avoided) by adjusting the shear rate. Therefore, it should be possible to control the ordering of flexible fibers in dilute systems.

VI Comparison with Jeffery orbits

In this Section, we will demonstrate that the rolling and tumbling modes can be interpreted as effective Jeffery orbits of the flexible fiber end-to-end vector ${\bf n}$ (see Fig. 1(b) for the notation), with an effective amplitude that depends exponentially on time. We will also show that the meandering mode is essentially different from the Jeffery solution. For the sake of clarity, we will focus on a fixed value of the bending stiffness $A=10$. The results can be easily generalized for a wider range of values.

VI.1 Jeffery orbits

To compare, we will first remind a reader of the classical Jeffery equations of motion of a rigid prolate spheroid with an aspect ratio $r$, oscillating periodically in shear flow Jeffery ; Graham . We will keep the same parametrization of the spheroid orientation $(\Theta_{J}(t),\Phi(t))$ as for the end-to-end vector of the flexible fiber, and assume for simplicity that $\tan\Theta_{J}(0)\!=\!C_{0},\;\Phi(0)\!=\!90^{\circ}$. The period of the motion, in units of $\dot{\gamma}$, is

[TABLE]

and the evolution of the orientation angles satisfies

[TABLE]

We will use the relation, which follows from Eq. (18),

[TABLE]

VI.2 Rolling and tumbling modes

In this section, we will briefly recall basic features of Fig. 2(a-b) and then analyze them both qualitatively and quantitatively, in a more general context. The characteristic variables and parameters in the rolling and tumbling modes will be labeled by (r) and (t), respectively. For both modes, the time dependence of the fiber end-to-end orientation angle $\Theta(t)$, shown in Fig. 2(a-b), illustrates two important features. First, the characteristic oscillation time is almost constant in time. Second, the amplitude of the oscillations changes monotonically with time.

We highlight important features of the dynamics in Fig. 9. In the inset we confirm that the oscillation time is practically time independent, with

[TABLE]

respectively. In the main panel of Fig. 9, we plot $\ln|\tan\Theta|$ versus time and obtain linear dependence of the maxima and minima on time. In this way we demonstrate that for the rolling mode, the amplitude of $\tan\Theta$ decays exponentially while for the tumbling mode, the amplitude of $\tan\Theta$ grows exponentially with time. The characteristic relaxation times are

[TABLE]

as illustrated in Fig. 9 by straight lines of the corresponding inclinations. We have checked that the values of the characteristic oscillation times $T$ and the relaxation times $\tau$ are practically the same for other initial conditions leading to the same mode.

Moreover, we will now demonstrate that

[TABLE]

with the periodic function $\Theta_{J}(t)$, which will be matched to the Jeffery solution given by Eq. (18), and negative or positive values of $\tau$ given by Eq. (21).

Matching the period $T$ of $\Theta_{J}(t)$, given by Eqs. (20), with the Jeffery period of a spheroid with the aspect ratio $r$, given by Eq. (16), we obtain (for $A$=10) $r\!=\!r_{r}\!\!=\!73.0$ as ‘the effective aspect ratio’ of the flexible fiber in the rolling mode. Note that this value is significantly larger than the geometrical aspect ratio of the deformed fiber (the spring constant $k$ was chosen to be so large that the fiber is practically inextensible). Moreover, the shape of the fiber significantly changes in time, while the period of the damped oscillations and ‘the effective aspect ratio’ remain almost constant in time, decreasing only by a few percent. For the tumbling mode, the effective aspect ratio following from Eq. (16) is much shorter, $r\!\!=\!\!r_{t}\!\!=\!\!28.8$222The tumbling motion in the shear plane of the end-to-end vector of a flexible fiber was matched with an effective Jeffery orbit in Ref. Slowicka_chaos_2015 ..

In Fig. 10 we compare $\Theta_{J}(t)$=$\exp(-t/\tau)\tan\Theta(t)$ and $\cos^{2}\Phi(t)$ evaluated numerically (solid lines) with the effective Jeffery expressions (17) and (19) (dashed lines), in which time is shifted back (by 5 and 13 units) to match the initial conditions and to compensate for transients, including the small decrease of $T$ with time in the first stage of the evolution. Also, for the same reason, we use the fitted values of $C_{0}$ given in the caption of Fig. 10. They slightly differ from the corresponding initial values: $C(0)$= 0.0152 for the rolling mode and $C(0)$=0.501 and for the tumbling mode.

For the rolling mode, the numerical and analytical curves are superimposed. For the tumbling mode, there appear some differences, but still $\Theta_{J}(t)$ can by approximated remarkably well by the Jeffery solution.

In this way, we arrive to one of the main conclusions of this paper. The rolling and tumbling modes can be interpreted as ‘effective Jeffery orbits’. The word ‘effective’ means that the standard expressions (17)-(19) for the Jeffery trajectories still hold, but now with a constant value of the parameter $C_{0}$ replaced by the time dependent amplitude $C(t)$,

[TABLE]

with the negative and positive values of $\tau$ corresponding to two different attractors of the dynamics. The decay of $C(t)$ to zero for the rolling mode and the growth to infinity for the tumbling mode describe, respectively, the rates of convergence to the rolling stationary state and to the periodic two-dimensional tumbling restricted to the shear plane.

To analyze the two different patterns of the evolution, and inspect the transition between them, we investigate other initial conditions. In particular, in Fig. 11 we choose the initial values of the amplitude $C$ much closer to each other, i.e. $C_{0}$=0.0163 and $C_{0}$=0.268. The different time dependence of $C$ in both modes is clearly visible in Fig. 11

where we use the simulation results to plot $\Theta$ versus $\Phi$. In this plot, the change of $C$ can be traced along the trajectories as the consecutive values of $\arctan\Theta$ at $\Phi=90^{\circ}$. For different initial conditions with $C_{0}>$0.25 and $C_{0}<$0.0195, the evolution in time is qualitatively the same as shown in Fig. 11. In the next subsection, we will discuss the dynamics for 0.0195$<C_{0}<$0.25, which corresponds to the empty space that separates the rolling and tumbling modes in Fig. 11.

VI.3 Meandering modes

The essential property of the periodic meandering solution is that it differs significantly from the Jeffery orbits Jeffery . In Fig. 12, we plot (in green) $\Theta(t)$ versus $\Phi(t)$ for the meandering mode, evaluated numerically for $\Theta_{0}=10^{\circ}$, $\Phi_{0}=10^{\circ}$. The corresponding initial value of $C$ is $C_{0}$=0.0307, with $r\!=\!r_{m}$=117, as evaluated from Eq. (16) for $T_{m}$=735. It is clear that the meandering trajectory is much more complicated than the periodic Jeffery orbit. In this case there is no systematic drift of the trajectory, and therefore the concept of ‘an effective Jeffery orbit’ cannot be applied. In Fig. 12, we also plot two black solid lines that follow from the Jeffery relation (17) as in Fig. 11. The meandering trajectory is almost everywhere located inside the region determined by these lines.

We investigated numerically evolution of flexible fibers with different initial orientations. For most of the initial conditions with 0.0195$<C_{0}<$0.25, we have found the meandering mode, which fills the empty space that separates the rolling and tumbling modes in Fig. 11.

VI.4 Physical interpretation

The analysis of the effective Jeffery orbits and the meandering modes, performed in the previous subsections, provides a simple physical explanation for the dependence of the modes on the initial orientation, determined in Fig. 8(b) for a moderate bending stiffness ratio $A$=10. In this diagram, the phase space of the initial orientation angles separates into three distinct regions, each of them leading to a different dynamical mode, and finally to a different dynamical attractor. The borders between these regions are separated by the two special Jeffery solutions with $C_{0}$=0.0195 and $C_{0}$=0.25, plotted as the black lines in Figs. 8, 11 and 12. The effective Jeffery solutions, convergent to the tumbling and rolling modes exist for a sufficiently large and a sufficiently small value of $C_{0}$, respectively. However, they do not exist for a range of intermediate values of $C_{0}$. In this range, periodic or close to periodic motions dominate.

This generic classification of the modes is perturbed by the existence of transient, close to periodic motions for some of the orientations and bending stiffness ratios $A$, as illustrated in Figs. 3-4.

In future work we hope to provide analytical justification for these new observations and quantitive results.

VII Discussion and conclusions

One of the main results of this paper is that there exist periodic and close to periodic three-dimensional motions of a flexible fiber in shear flow. The meandering periodic solution is an attractor for a certain range of the relative bending stiffness $A$. The squirming motion is a transient in this range, and an attractor in a narrow range of smaller values close to $A$=7.3. Our results indicate that a change of the relative bending stiffness $A$ may trigger a transformation between different periodic (or close to periodic) solutions, and the corresponding change of the characteristic sequence of three-dimensional fiber shapes. We have demonstrated that the complexity of a flexible-fiber shape (and in particular the number of local maxima of the curvature) may change significantly with time. A more detailed study of periodic and close to periodic motions of flexible fibers in shear flow, for a wide range of the bending stiffness $A$, will be presented elsewhere.

We have shown also that the time-dependent rolling and tumbling modes of a flexible fiber in shear flow can be interpreted as effective Jeffery solutions, with the constant $C_{0}$ replaced by an exponential function of time $C(t)$ given by Eq. (23). Therefore, the effective Jeffery orbits drift towards one of two attractors: the fiber aligned with the vorticity direction and the fiber performing periodic motions entirely in the shear plane. There exist two thresholds for the initial value of $C(t)$: small $C(0)$ lead to the rolling mode, intermediate $C(0)$ to the meandering mode, and large $C(0)$ to the tumbling one. These thresholds are sensitive to the bending stiffness ratio $A$.

Unlike the slender body theory duRoure ; Liu , the Hydromultipole method used in this study takes into account the fiber thickness. This feature is important to study the three-dimensional dynamics in shear flow, and in particular, the effective Jeffery orbits and the tumbling motion, in contrast to the elastica models BeckerShelley ; YoungShelley that predict an infinite tumbling time of infinitely thin fibers.

Finally, we comment on possible comparisons with experiments. In Refs. Harasim2013 ; Liu , quasi-2D trajectories of actin were investigated, with the focal plane of the motion perpendicular to the vorticity direction. The elasto-viscous number $\bar{\mu}=8(L/d)^{4}/A\approx 2.2\times 10^{6}$, i.e., the largest value of $\bar{\mu}$ analyzed in Ref. Liu , corresponds to $A$=10 for our system. Indeed, for $A$=8 we recover a shape evolution similar to Movies 5 and 6 from Liu , with U turns and S turns that appear irregularly.

However, we observe that more stiff fibers, after a long time, approach the periodic tumbling motion, with the S-turns only, as illustrated in Fig. 13. These shapes are different from the C-turns observed in Liu . A more detailed study is needed to understand the reason for such a difference. The obvious differences are that the fibers

we study do not perform Brownian motion and are much shorter than actin. Moreover, in Ref. Liu , the simulations were performed in 2D and the experiments were carried out in a bounded geometry. Finally, we note that even a very small curvature of the external flow can have a drastic influence on flexible fiber shapes farutin .

As a last remark, we suggest searching for three-dimensional, periodic motions in experiments. Our estimates indicate that they appear in the range of the bending stiffness accessible in experiments. With this goal in mind, we highlight the need to perform video recordings in the $xy$ (flow-vorticity) plane where the characteristic repeatable evolution of shapes should be visible. We hope to report the results of such experiments in a future communication.

Acknowledgements.

We thank Olivia du Roure for helpful discussions. A.M.S. and M.L.E.J. were supported in part by Narodowe Centrum Nauki under grant No. 2014/15/B/ST8/04359. H.A.S. thanks the NSF CMMI-1661672. We benefited from the ITHACA project PPI/APM/2018/1/00045 financed by the Polish National Agency for Academic Exchange.

Appendix A Periodic meandering attractors

In this section, we provide more details about the periodic meandering solution. In Fig. 2(c) the maximal values of the relative length of the end-to-end vector are close to one. In the figure, by looking at the orientation $\Theta(t)$ we observe that during the periodic meandering motion, the fiber straightens while tumbling at an orientation that is not along the flow. This property is also visible in Fig. 7 where fiber shapes are shown. In particular, for the periodic meandering solution with $A$=10, the fiber is almost straight for $\Phi$=0 and $\Theta\!\approx\!78^{\circ}$. In contrast, for the tumbling solution, the fiber straighten along the flow, i.e. at $\Theta=90^{\circ}$.

In the meandering motion, every bead performs a periodic orbit superposed with translation along the shear flow with a velocity v${}_{\text{a}}$, equal to the mean velocity of the fiber center-of-mass, averaged over the period. The periodic meandering motion of the center-of-mass for $A\!\!=\!\!10$ is shown in Fig. 14 and Movie 4, in the frame of reference translating along the shear flow with velocity v${}_{\text{a}}$. Here, $A=10$, and the initial conditions are the same as in Figs. 2, 3, 6, 7 and Movie 2. The time range is $99188\!\!\leq\!\!t\!\!\leq\!\!99988$, with the period $T_{m}\!\!=\!\!735$. Owing to symmetry, trajectories of the end-to-end vector, displayed in Fig. 6, close after $T/2$. Our results show that the meandering period is longer than the transient squirming period and much longer than the tumbling period.

Appendix B Close-to-periodic meandering attractors

In Fig. 15 we illustrate that by changing values of the bending stiffness $A$ away from $A$=10, we observe close to periodic (rather than periodic) meandering motions, which are long-lasting and vary similarly to the periodic meandering solution for $A$=10. For $A$=9, departure from periodicity has a different pattern: the system oscillates back and forth between the twin meandering orbits, symmetric with respect to $\delta x\rightarrow-\delta x$ and $\delta z\rightarrow-\delta z$.

The end-to-end periodic meandering trajectory, shown in Fig. 6, is non-symmetric. Naturally, after rotation by

$\pi$ around the $y$ axis, one obtains another periodic orbit. Both appear in our simulations (compare Figs. 3-4 and 6). Moreover, we observed long time effects of flipping from one to the other trajectory, as shown by the example in Fig. 16 (again, for $\Theta_{0}=\Phi_{0}=10^{\circ}$).

The periodic meandering mode for $A$=10 seems to be stabilized already for times around 2000. However, as shown in Fig. 16, there appear some long-time effects: the trajectory later flips to another one, symmetric with respect to rotations along the vorticity direction (i.e., with respect to $\delta x\rightarrow-\delta x$ and $\delta z\rightarrow-\delta z$), and later remains in this shape.

Appendix C Accuracy of the simulations

The Stokes flow dynamics of an elastic fiber modeled as a chain of beads, evaluted using the Hydromultipole numerical algorithm, is more accurate for a larger multipole truncation order $LL$ and a smaller distance $l_{0}$ between the centers of consecutive beads. Numerical accuracy of the simulations of flexible fibers in shear flow, performed with the Hydromultipole numerical codes cichocki , was analyzed in detail in Ref. Slowicka_chaos_2015 for the periodic tumbling motion of fibers located in the $xz$ plane. The accuracy of the truncation at $LL\!=\!2$ was estimated there as 3%, and for the $l_{0}\!=\!1.02$ approximation as 10-15%. It is essential to keep the gap between consecutive bead surfaces small enough to prevent the beads from rotating under shear independently from each other. For small gaps, the lubrication effects serve this goal quite effectively.

Similar accuracy estimates are expected to hold for the three-dimensional motions considered in this article. However, in this case, the analysis of the accuracy is more subtle because we have already observed that the existence of periodic meandering solutions is limited to a narrow range of the bending stiffness. A small departure from this range gives rise to motions that are close to periodic (for example, see Fig. 15). Therefore, taking into account the chaotic nature of flexible fibers dynamics in shear flow Slowicka_chaos_2015 , the basic question is if three-dimensional periodic meandering motions are also observed for more precise (but computationally more demanding) bead models of a flexible fiber, with a smaller $l_{0}$.

Therefore, we choose $A\!=\!9,\,11$ as the values of the bending stiffness close to $A\!=\!10$, which for $l_{0}\!=\!1.02$ corresponds to the meandering attractor with the period $T\!=\!735$. In Fig. 17 and Table 1 we compare the long-time solutions obtained for $l_{0}\!=\!1.005$ and $l_{0}\!=\!1.02$. While changing $l_{0}$, we fix values of $N$ and $k$, and consider the identical initial orientations $\Theta_{0}\!=\!\Phi_{0}\!=\!10^{\circ}$.

As shown in Fig. 17 and Table 1, for $l_{0}\!=\!1.005$, the periodic meandering solution is not recovered at $A\!=\!10$; instead, for these values of the parameters, we find the close-to-periodic meandering motion, with a smaller maximal time $\tau\!=\!636$ between the consecutive tumblings. The decrease of the tumbling time for smaller values of $l_{0}$ was also observed for the tumbling periodic orbit Slowicka_chaos_2015 , and explained as the influence of the bead rotation. To minimize this effect, the gaps between the bead centers should be kept reasonably small. In practice, $l_{0}\!\lesssim\!1.02$ is sufficient.

The long-time periodic meandering solution, detected for $l_{0}\!=\!1.02$ and $A\!=\!10$ with the period $T\!=\!735$, is observed also for the smaller distance $l_{0}\!=\!1.005$, but at a larger value of the bending stiffness $A\!=\!11$, and with a smaller period $T\!=\!653$. Movie 5 illustrates that the evolution of shapes in these two periodic solutions is almost identical (here time is normalized by the two different periods, respectively). We also find that the long-time meandering quasi-periodic solutions for $l_{0}=1.02$ are similar to those obtained for $l_{0}=1.005$, but with larger values of the bending stiffness $A$, as shown in Fig. 17 and Table 1.

Appendix D Description of the movies

Movie 1.

Typical evolution of fiber shapes in the relaxation phase of the rolling (blue), tumbling (red) and meandering (green) modes. The bending stiffness $A=10$ and the initial orientations are the same as in the simulations shown in Fig. 2. The movie starts at the same time instants as the snapshots in Fig. 2 but it lasts longer; 500 time units are shown for each mode (time is synchronized). Close to periodic squirming motion in the relaxation phase of the meandering mode is clearly visible.

Movie 2.

Typical evolution of fiber shapes in the attracting periodic orbits: meandering (green) and tumbling (red) solutions are reached after a long time. The bending stiffness $A=10$ and the initial orientations are the same as in Figs. 2-3 and Movie 1. Time $t$ is synchronized, with $t_{0}\leq t\leq t_{0}+736$ where $t_{0}=99188$ for the meandering and $t_{0}=4653$ for the tumbling solutions.

Movie 3.

Illustration of the concept of the fiber end-to-end trajectories. A trajectory is drawn by the tip of the fiber end-to-end vector for the periodic meandering solution. Here $A=10$, and the initial conditions are the same as in Fig. 3 and Movie 2. The time range is $99188\leq t\leq 99988$.

Movie 4.

Periodic meandering motion of a flexible fiber center-of-mass, in the frame of reference translating along the shear flow with velocity v${}_{\text{a}}$, equal to the mean velocity of the fiber center-of-mass, averaged over the period. Here, $A$=10. The corresponding trajectory is shown in Fig. 14.

Movie 5.

Testing the accuracy of our bead model. The sequence of fiber shapes in the periodic meandering long-time solution for $l_{0}\!=\!1.02$ and $A\!=\!10$ is almost the same as for a smaller (and more accurate) value $l_{0}\!=\!1.005$ and $A\!=\!11$. The time is normalized by the corresponding periods: $T\!=\!735$ for $l_{0}\!=\!1.02$ and $T\!=\!653$ for $l_{0}\!=\!1.005$.