Shape evolution of long flexible fibers in viscous flows

The present work studies numerically the dynamics and shape evolution of long flexible fibers suspended in a Newtonian viscous cellular flow using a particle-level fiber simulation technique. The fiber is modeled as a chain of massless rigid cylindrical segments connected by ball and socket joints; one-way coupling between the fibers and the flow is considered while Brownian motion is neglected. The effect of stiffness, equilibrium shape, and aspect ratio of the fibers on the shape evolution of the fibers are analyzed. Moreover, the influence of fiber stiffness and their initial positions and orientations on fiber transport is investigated. For the conditions considered, the results show that the fiber curvature field resembles that of the flow streamline. It is found that the stiffer fibers experience not only a quicker relaxation phase, in which they transient from their initial shape to their “steady-state shape,” but they also regain their equilibrium shape to a larger extent. The findings also demonstrate that even a small deviation of fiber shape from perfectly straight impacts significantly the early-stage evolution of the fiber shape and their bending behavior. Increasing the fiber aspect ratio, when other parameters are kept fixed, leads the fiber to behave more flexible, and it consequently deforms to a larger extent to adjust to the shape of the flow streamlines. In agreement with the available experimental results, the fiber transport studies show that either the fiber becomes trapped within the vortices of the cellular array or it moves across the vortical arrays while exhibiting various complex shapes.


Introduction
The motion of long flexible fibers in fluid flows is central to many industrial applications, including pulp and paper making, textile, insulation industry, as well as microfluidic engineering and biotechnologies [1][2][3][4]. In addition to these applications, suspensions of long filaments are also a matter of interest in nature and some biological systems, for example in cell division [5] and swimming of microorganisms [6]. The flow-induced deformation of fibers is essential for fiber transport [7], affecting the isotropy and spatial dispersion of fibers. These impacts on the microstructure of the fibers suspension strongly affect the macroscopic properties (e.g., mechanical, thermal, electrical) of the final product [8][9][10]. Hence, unveiling the fundamental mechanisms and governing parameters behind these phenomena remains crucial to improving the design and optimization of products and processes.
There has been numerous experimental, theoretical, and numerical research to investigate the behavior of rigid spherical or non-spherical particles in various flow conditions. When objects have multiple degrees of freedom, which is the case for flexible fibers, the problem becomes increasingly challenging due to the formation of complex fiber structures created by flow-induced bending, buckling, twisting, and flocculation. Early experiments by Forgacs and Mason [11,12] with elastic fibers in simple shear flow demonstrated different shapes and motion behaviors depending on the characteristics of the suspension system, including fiber geometry, fiber stiffness, and flow properties. They identified various regimes of fiber motion, including the rigid motion in which the fiber follows Jeffery orbits [13], periodic orbits with springy and snake-like deformations, complex coiled orbits, and coiled motion with entanglement. With technological advances in microfluidic technology, the periodic shape deformations were studied experimentally for DNA strands and actin filaments [14,15], and complex coiling and entanglement of long actin fibers in shear flows [16]. Wandersman et al. [17] studied the deformation and transport of flexible fibers in a viscous cellular flow, created by a planar array of electromagnetically driven counter-rotating vortices. They observed both fiber trapping in vortices, and fiber transport and buckling through hyperbolic stagnation points for the sufficiently high fluid strain-rate. They also reported on the effect of fiber deformation on its translational motion, and identified the buckling threshold in a reasonable agreement with the predictions of linear analysis in a purely hyperbolic flow [17].
The Slender Body Theory (SBT) [18], Immersed Boundary Methods (IBM) [19,20], and bead and rod models [2,21,22] have been the three main methods used to study numerically the dynamics of suspended flexible fibers in a fluid flow [23]. Reviews [24][25][26][27] have presented the basics and algorithmic advances of these computational modeling approaches. Much progress has been made in employing SBT to describe the orbits and buckling of fibers in fluid flows [7,16,18,28]. In particular, Quennouz et al. [7] studied the dynamics, transport, and buckling of elastic fibers in a two-dimensional cellular flow experimentally and numerically using SBT. They investigated the probability of fiber buckling by varying the parameter so-called elastoviscous number, which represents the relative intensity of viscous forces to fiber elastic forces. The typical shape evolution of fibers in shear flow was studied by Nguyen and Fauci [29], Słowicka et al. [30], and LaGrone et al. [31]. They observed that an elastic fiber often tends to the flow-gradient plane, and its shape development in the shear flow depends on the fiber slenderness and the elasto-viscous number. The fiber rotates rigidly and nearly-rigidly in the shear for small elasto-viscous number; it buckles once the elasto-viscous number reaches the threshold value; and a transition to a "snaking behavior" occurs for larger elasto-viscous numbers [31]. By modeling a fiber as a chain of spherical particles, Kuei et al. [32] studied the complex coiling and formation of knots in shear flow. Kondora and Asendrych [33] employed a one-way coupled particle-level fiber method, where the fibers were modeled as a chain of spheres connected by ball and socket joints, to obtain the probability distribution function of fiber orientation in a converging channel of a papermachine headbox. The influence of fiber stiffness and the compressional flow field on the shape of a deformed fiber has been studied by Chakrabarti et al. [34]. In particular, they reported that flexible filament could buckle into a helical shape in a purely compressional flow when end moments and compressive forces are applied, standing out from classical helical buckling where the filament may obtain a chiral helicoidal morphology in the absence of any intrinsic twist or external moments.
Understanding how non-spherical structures interact with a turbulent flow is quite complicated. Beyond the complexity of turbulence of flows with spherical particles, the forces and torques that depend on the particle shape and orientation should also be considered [35,36]. The majority of the studies dealing with non-spherical particles in turbulence focus on rigid spheroids [3,37,38], reporting the tendency of the particles to align with the direction of the strongest Lagrangian fluid stretching [37,39]. It has been shown that the particle shape and size affect the inertia-driven phenomena such as preferential concentration and near-wall accumulation [36,38,40]. The flexibility is another parameter that can affect the behavior of high aspect ratio particles suspended in a turbulent flow since the particle may experience severe flow-induced velocity gradients over its length [41].
Few studies are available on the suspension of flexible particles in turbulent flows due to the complexities addressed earlier. Brouzet et al. [42] and Gay et al. [43] examined the spatial conformation of a flexible fiber in a homogeneous isotropic turbulence flow experimentally. They proposed a model for the transition from rigid to flexible regimes as the turbulence intensity increases or the fiber elastic energy decreases. They found that the straightening of the fibers, which becomes more evident as the length of the fiber increases, is due to the spatial and temporal correlations of the turbulent forcing. Niskanen et al. [44] performed both experiments and numerical simulations of a dilute suspension of flexible wood fibers in a contracting turbulent channel flow to identify the mechanisms affecting the development of fiber orientation. They observed that the orientation anisotropy obeys approximately a linear relation to the contraction ratio, meaning that the streamwise flow acceleration is the main mechanism affecting the development of the fiber orientation. From a numerical point of view, Andric et al. [45] studied the translation and re-orientation of the fibers in a turbulent channel flow. They found that the fiber motion is primarily governed by velocity correlations of the flow fluctuations. They also highlighted the tendency of the fibers to evolve into complex geometrical configurations. Dotto et al. [41] performed a statistical characterization of the bending dynamics of inertial flexible fibers dispersed in a turbulent channel flow, finding that fibers retain a more deformed configuration in the center of the channel while they stretch near the wall region by the effect of the mean shear. Moreover, they found that the fiber end-to-end distance reaches a steady-state regardless of fiber location with respect to the wall. Rosti et al. [46,47] studied the interaction between fiber elasticity and turbulent flow. They identified the flapping regime in which the fiber, despite its elasticity, is effectively slaved to the turbulent fluctuations, and therefore, they can be used as a proxy to measure various two-point statistics of turbulence.
To the best of the authors' knowledge, there is no study available analyzing the behavior of flexible fibers in the cellular flow that employs a particle-level fiber model where a fiber is represented as a chain of cylindrical rods. To contribute to fulfilling this limitation, the present study is built on the model proposed by Schmid et al. [48] and Switzer III and Klingenberg [49] to investigate the dynamics of long flexible fibers suspended in a Newtonian viscous cellular fluid flow. Such a flow setup has been proposed to capture some of the dynamics of biological polymers in motility assays [50]. Additionally, understanding the fiber configurational transitions has been proven to provide useful information for biological and fluid dynamics communities, for instance, in emerging nanotechnology for directed transport, assembling, sorting, and sensing of molecules [51][52][53][54], and nanometer-sized cargo in cellular-sized synthetic devices, as well as in sorting pulp fibers in the papermaking process [20]. In the model used, each fiber is treated as a chain of massless rigid cylindrical segments connected by ball and socket joints. The model includes realistic features of fibers, such as flexibility and irregularly deformed equilibrium shapes, and, in comparison with the previous studies, it can particularly provide additional information on the behavior of the fibers with respect to their equilibrium shape. Moreover, the mode considers one-way interaction between the fibers and the fluid flow while fiber-fiber interactions and Brownian motion are neglected. On each segment, the hydrodynamic forces and torques, the bending and twisting torques between each neighboring pair of segments, as well as the internal constraint forces and moments are applied. The purpose of this work is to study the dynamics and the deformation of fibers in a steady cellular flow, depending on the fiber stiffness, aspect ratio, and equilibrium shape. To this end, the shape evolution of fibers initially aligned with the flow direction is evaluated using the curvature of the fiber and the curvature of fluid flow at different distances from the center of vortices. Moreover, to gain a general understanding of the effects of fiber flexibility and initial positions and orientations on fiber transport, the analysis of fiber trajectories is performed.

Simulation method
The fiber model used in this work is based on the model proposed by Schmid et al. [48] and Switzer III & Klingenberg [49]. Thus, for the purpose of completeness, this Section describes the main features of the model and highlights the modifications applied in the present work.
As illustrated in Fig. 1, a flexible fiber is modeled as a chain of N rigid inextensible cylindrical segments connected by N − 1 hinges. All fiber segments have the same diameter D and length l, representing the fibers with aspect ratio r f = L/D = Nl/D. The position vector r i points to the center of mass of the segment i, and the unit vectorẑ i describes the segment orientation with respect to the non-inertial coordinate system. Fiber and fluid inertia are neglected because of the small Reynolds numbers under creeping flow conditions (here Re = ργ L D/μ 1 where ρ,γ , and μ are the density of the fluid, shear rate, and fluid viscosity, respectively) [48]. Moreover, given that the behavior of the fibers is studied for dilute flow condition, the hydrodynamic interactions between the fibers are neglected. The motion of each fiber segment i is governed by Euler's first (Eq. (1)) and second law (Eq. (2)) together with the connectivity constraint (Eq. (3)) as follows: In Eq. (1), F h i is the hydrodynamic drag, F b i is the body force, X i and X i+1 are the connectivity forces in the hinges at the two ends of the segment i. In Eq.
i ) and twisting torque ( Y t i ) in the joint i, and Y i+1 is the analogous for the joint i + 1. For the end segments X 1 = X N +1 = 0 as well as Y 1 = Y N +1 = 0.
The hydrodynamic force and torque on the segment are , respectively, where A i , C i , and H i are the resistance tensors that are approximated by the corresponding resistance tensors of a prolate spheroid with an equivalent aspect ratio proposed by Kim and Karrila [55]; i and Y t i ) represent the fiber resistance to the flow acting to restore its equilibrium shape. These torques are assumed to be proportional to the difference between the angle at the joint i formed by connected segments: here k b and k t are the bending and twisting stiffness constants,ê b i andê t i are the unit vectors in the bending and twisting torque directions. θ and φ are the bending and twisting angles of the segment whereas θ eq i and φ eq i are the equilibrium angles describing fiber equilibrium shapes when no forces or torques are acting on it (see Fig. 2). For an intrinsically straight fiber, θ eq i and φ eq i are zero, for a permanently deformed U-shaped fiber with an intrinsic radius of curvature R u , θ eq = 2tan −1 (L/(2R u )) is constant at every joint and φ eq i = 0, and in the case both θ eq i and φ eq i are non-zero constants at every joint, the fiber has an intrinsic helical shape [48].
The reader is referred to Ref. [48] for details of the discretization of the motion equations and the numerical algorithm to solve them.
Assuming the fluid flow velocity vector U (t) = u x , u y , u z , the curvature of the flow streamlines can be calculated [56]: Here × is the cross product, and a (t) is the acceleration field that for a steady flow can be computed as: To analyze the evolution of the fiber shape, the curvature of each fiber segment κ i is evaluated as where R i is the radius of the circle that is determined by the center of mass of the three consecutive segments.

Validation
The fiber model validation is performed by comparing the simulation results for both rigid and flexible fibers to well-known theoretical and experimental results. Jeffery [13] studied the motion of an isolated rigid ellipsoidal particle in a simple shear flow U = (γ y, 0, 0), whereγ is the shear rate of the flow (see Fig. 3). The periodic trajectories of the particle ending points, known as Jeffery orbits, can be described in spherical coordinates by where C is the orbit constant, k is the phase angle, T = 2π r s + r −1 s /γ is the orbit period, θ is the angle between the unit vector along the cylinder (ẑ) and the vorticity axis (z-axis), and φ is the angle between the projection of the unit vectorẑ on the flow-gradient plane and the x-axis (see Fig. 3). Bretherton [57] expanded Jeffery's solution to any axisymmetric particle when r s is replaced with an effective aspect ratio r e . Cox [58] found the semi-empirical correlation r e = 1.24r f / ln r f for a circular cylinder with an aspect ratio of r f . Figure 4 illustrates the trajectories of one end point of an isolated rigid fiber centered at the origin for different orbit constants C, obtained by varying the initial Euler angles (θ and φ). The solid lines are the path predicted from Jeffery's solution (Eqs. (8) and (9)), and the markers are from the simulations of a rigid fiber with N = 11, r f = 110, and B R = 2 (the parameter B R describes the fiber flexibility and can be calculated f where E Y is Young's modulus and r f is the fiber aspect ratio). It is seen that there is an excellent agreement between the numerical results and Jeffery's theory for the rigid fiber. Figure 5 shows the non-dimensional Jeffery orbit period Tγ as a function of fiber aspect ratio r f for both predictions from Cox's semi-empirical correlation and the one computed from the numerical model employed. The simulation results reported are from an isolated straight fiber (i.e., θ eq = φ eq = 0) with B R = 2 (corresponding to an essentially rigid fiber) and various fiber aspect ratios. The fiber was initially placed aligned with the direction of a simple shear flow, and, thus, it rotates in the gradient plane of the flow. As it can be seen in Fig. 5, a good agreement between the computed periods and the theoretical prediction is found for the range of fiber aspect ratios analyzed.
The implemented method was also used to reproduce the motion of an isolated flexible fiber in shear flow for different fiber stiffnesses. Forgacs and Mason [12] experimentally studied the motion of threadlike particles in shear flows. By varying shear rate, fiber stiffness and fiber length, they observed different rotation regimes such as rigid, springy, snake-like, helix, and coiled motion with self-entanglement. Figure 6 shows the orbiting    . 6 The time series of images from simulations showing the shape evolution of fibers in simple shear flow. Three types of orbits can be seen for a fiber with r f = 250, N = 25: a "Rigid," for a fiber with intrinsic radius of curvature R u = 100L (i.e., θ eq ≈ 0.01) with B R = 2; b "Springy," for a fiber with intrinsic radius of curvature R u = 100L (i.e., θ eq ≈ 0.01) with B R = 1; c "Snake-like S-shaped," for an intrinsically straight fiber with B R = 0.1; d "Snake-like C-shaped," for a fiber with intrinsic radius of curvature R u = 100L (i.e., θ eq ≈ 0.01) with B R = 0.1. For the sake of visualization, the fiber diameters are exaggerated Fig. 7 Surface plots of the (absolute value of) curvature as a function of time and fiber arclength; panels a-d correspond to the cases depicted in Fig. 6a-d, respectively behavior of a fiber with r f = 250, N = 25 (to ensure sufficient discretization for capturing fiber deformation) and three different values of the bending ratio B R = 2, 1, 0.1. It is predicted that a cylindrical fiber bends when B R < 1 [48], however, in agreement with Schmid et al. [48], the modeled intrinsically straight fiber behaves stiffer than actual fibers, and it bends at B R < 0.1. Additionally, an intrinsically straight fiber exhibits a snake-like S-shaped deformation (Fig. 6d). However, implementing a small permanent deformation to mimic the geometrical imperfection of physical fibers (here, R u = 100L [48]), leads to a Snake-like C-shaped fiber deformation (Fig. 6c), as reported in the experimental studies [12,59]. In conclusion, the numerical results are qualitatively in agreement with experimental observations of Forgacs and Mason [12]. Figure 7 depicts the surface plots of the evolving (absolute value of) curvature along the arclength of the fiber as a function of time, corresponding to the simulation runs depicted in Fig. 6. A rigid fiber (Fig. 7a) maintains its initially straight shape, and therefore there is no deformation in its curvature surface plot. For the S-shaped motion (Fig. 7c), the fiber buckles, and two bends occur during rotation. In this case, the positions of the maximal curvature do not travel along the arclength of the fiber while they occur as standing waves which appear periodically in each orbit. In the case of snake-like C-shaped deformation (Fig. 7d), in any orbit, the fiber bends from one endpoint and the point of maximal curvature travels along the fiber arclength toward another fiber endpoint until it regains its straight shape. In all simulations shown in Fig. 7, the fiber exhibits periodic motion; it remains most of the time in the straight configuration and aligned with the flow direction.

Results and discussion
In this work, the behavior of flexible fibers flowing in one cell of a lattice of counter-rotating vortices [17] is studied. Since the flow streamlines are closed, inertialess point-like particles became trapped in the vortices. Moreover, this vortical flow contains stagnation points generating compression-elongation rates, and the fiber may deform, which makes this flow an interesting choice to analyze the shape evolution and transport of flexible fibers.
The two-dimensional steady velocity field of the cellular flows is described by [17] Here x (resp. y) is the dimensionless coordinate π x/2W (resp. π y/2W ) in which 2W is the cell size. In this study, U 0 = 1 and W = 600π, to have the simulation box with side lengths of at least four times larger than the longest fiber considered in this study (2W ≥ 4L max ) [48]. Figure 8a depicts the velocity field, while Fig. 8b shows the vorticity field of one cell of the background cellular fluid flow in a square box computational domain of size 1200π. Figure 9a illustrates the streamline curvature field and Fig. 9b shows, as an example, the time-averaged fiber curvature field computed from a  In order to compute the fiber curvature field, the computational domain was divided into 80 × 80 identical square cells, and the fiber curvature was calculated by averaging the curvatures of the segments in which their centers of masses lie within the cell. By comparing Fig. 9a, b, it is seen that, for the conditions considered N , r f , DS, θ eq , φ eq = (11, 220, 0.1, 0.01, 0), the fiber curvature field is is in qualitative agreement with the streamline curvature field. Several numerical experiments were performed to gain a general understanding of the effects of initial fiber positions and orientations, as well as fibers stiffness on the transport and shape evolution of flexible fibers in a viscous cellular flow. The motion of fibers initially placed at various distances from the center of vortices b = 1000, 1200, 1500, 1700, 1800, 1850, 1870, 1880, and oriented at two different angles measured from the x-axis of α = 0 • (i.e., along x-axis) and α = 90 • (i.e., along y-axis) was tracked (see Fig. 8b for the definitions of α and b). Figure 10 shows the trajectories of the center of mass of the fibers. In addition, the observations of the selected cases with (b, α) = (1700, 0 • ), (1800, 0 • ), and (1850, 0 • ) are provided as the Supplemental Material Movie S1. In agreement with the study of Wandersman et al. [17], two main types of fiber dynamics were observed. In one dynamical state, the fiber is trapped within the vortex (Fig. 10a, b) and follows roughly the flow streamlines. The trapped fibers demonstrate oscillations about the mean trajectory, which can span the trajectory of the fiber toward the cell compression regions. Another dynamical state occurs for the fibers that are initially located either very close to or inside the compression region (Fig. 10c, d). In this state, the fibers meander across the vortical arrays and show transitioning from orbiting about one vortex to another one [7]. Similar to the observation of Wandersman et al. [17], fibers also show a staircase-like motion across the lattice of vortices (see the case in cyan in Fig. 10d). Figure 11 depicts successive snapshots taken from the instantaneous shape of two isolated fibers with different dimensionless stiffnesses DS = π B R/[32(ln 2r e − 1.5)] [48,58]. Both fibers are placed at b = 1850 and aligned with the y-axis (i.e., α = 90 • ). The rigid fiber with DS = 1 (on the left) does not buckle and maintains its straight shape, while the flexible fiber with DS = 1 × 10 −2 (on the right) buckles and shows different modes of deformation. Since the flow streamlines are closed and the fibers are without inertia, they remain trapped in the vortices of the cell.
As illustrated in Fig. 11 and in Supplemental Material Movie S1, depending on the fiber initial conditions and stiffness, the fiber undergoes different deformation modes, including straight, nearly straight, "J", "C", and "U" shapes.
The effects of fiber stiffness, fiber equilibrium shape, and fiber aspect ratio on the shape deformation of flexible fibers in cellular flow are analyzed as well. To perform such a parametric study, the fibers were considered to be made from vinylpolysiloxane, as reported in [7], with a density of 1.1 × 10 3 kg m −3 , a fixed radius of 41 μm, and different lengths ranging from 5.8 mm to 23 mm. The fibers considered in Ref. [7] have Young's modulus ranging from 20 kPa to 570 kPa ; however, this study analyses a wider range [0.  kPa to elucidate the effect of the fiber stiffness using a more inclusive series of stiffnesses. The fibers equilibrium shape can vary from straight fiber (θ eq = φ eq = 0) to a relatively U-shaped one by taking equilibrium bending angles of θ eq = 0.1 and 0.2 (see Fig. 2). The fluid medium was considered to be a mixture of polyethylene glycol 1000 (Fluka) and purified water, exhibiting a viscosity of 82 mPa s and density 1.25×10 3 kg m −3 [7]. In all the following simulations, 64 fibers each consisting of 11 segments are tracked for different combinations  Figure  12 shows the time evolution of the ratio between fiber curvature to the curvature of flow streamline ξ = κ f iber /κ f low (Eq. (5) and (7)) as a function of time and the distance from the center of vortices (b). The standard deviations of the curvature data are indicated by error bars. As shown at the top of the left column of Fig. 12, at dimensionless time of 1000, the value of ξ is very close to one for various distances from the center of vortices. The fibers are initially placed aligned with the flow direction and, therefore, at the early stages of the simulations, the fibers still have almost the same curvature as the flow streamline curvature (i.e., ξ ≈ 1). As time evolves, the restoring bending torques between each pair of neighboring segments (Eq. (4)) push the fiber to regain its equilibrium shape. In this regard, the fibers having higher dimensionless stiffness behave more rigid and restore their equilibrium shapes (θ eq = 0.01) to a larger extent, therefore, they exhibit lower κ f iber and ξ accordingly. In contrast, in the case of very flexible fibers (e.g., DS = 4.2 × 10 −6 ), the elastic forces cannot overcome viscous forces, and the fibers do not regain their equilibrium shape. Therefore, they remain almost aligned with the flow direction, with ξ ranging from 0.95 to 1 at any time.
As time progresses, the fiber relaxation phase occurs gradually, and the fiber deforms to regain its equilibrium shape. Figure 12 illustrates that, as time advances, the fiber straightens and eventually reaches its "steady-state shape" that is a slightly U-shaped fiber with θ eq = 0.01. The fibers, therefore, form smaller curvatures, and the value of ξ = κ f iber /κ f low decreases over time at any specific position (b). The simulation results presented in Fig. 12 suggest that the stiffer fibers exhibit a quicker relaxation phase. For instance, after 4 × 10 3 dimensionless time, the fibers with DS = 4.2 × 10 −2 reach their steady-state shape, while the fibers with DS = 4.2 × 10 −3 , reach their steady-state shapes just after 24 × 10 3 dimensionless time. Here, the steady-state shape is determined as the time of which the values of ξ became stable over time and their corresponding standard deviations are very small.
As it can be seen in Fig. 13, for the case with DS = 4.2 × 10 −3 at long times (e.g., 34 × 10 3 dimensionless time), the ξ is less than 0.1 for the fibers close to the center of vortices (b < 200), and it increases up to ξ = 0.7 for the fibers further from the center. This outcome is as to be expected since all the fibers have reached their  Fig. 13).

Effects of fiber shape
Simulations were performed for the systems containing U-shaped fibers for several different equilibrium bending angles. The results are presented for numerical experiments of flexible fibers with θ eq = 0, 0.1, 0.2 (see Fig. 2), while other parameters are kept fixed at (DS, r f , N ) = (4.2 × 10 −3 , 143, 11), corresponding to fibers with E y = 20 kPa, L = 11.72 mm, and D = 8.2 µm. The effects of the fiber equilibrium shape on the value of ξ as a function of time and the distance from the center of vortices b are illustrated in the middle column of Fig. 12. Although the fibers were initially placed aligned with the flow direction, even at the very beginning of the simulation (e.g., t = 1000), the value of ξ for all cases is between 0.5 and 0.55. This is due to the strong impact of the fiber equilibrium shape on the bending behavior of fibers. This effect can be observed in the Supplemental Material Movie S2, where the motion of the fibers with θ eq = 0.1 is displayed for the total dimensionless time of 50 × 10 3 . For the sake of visualization, the fiber diameters are exaggerated in the movie. There is a quick transition from the initial condition, where the fibers inherently represent ξ = 1, to ξ ≈ 0.5 after 2.5 × 10 3 dimensionless time.
In the system containing fibers with θ eq = 0, as time evolves, the fibers deform to regain their steady-state shapes, resulting in a gradual decrease in ξ . Figure 12 demonstrates that, after 24 × 10 3 dimensionless time the fibers are relaxed and almost perfectly straightened (i.e., κ f iber ≈ 0, hence ξ ≈ 0). However, as shown Fig. 12, the fibers with θ eq = 0.1 experience values of ξ ranging from 0.5 (at b ≈ 200) to 0.8 (at b ≈ 1200). Thus, a small amount of intrinsic curvature can have a significant effect on the bending behavior of the fibers.

Effects of fiber aspect ratio
The plots in the right column of Fig. 12  , φ eq ) = (11, 4.2 × 10 −3 , 0.01, 0) (here, Young's modulus of the fibers was considered to be E y = 20 kPa). As it can be seen in Fig. 12, at a given distance from the center of vortices, the fibers with a lower aspect ratio exhibit lower ξ . For the studied cases, this is to be expected; as the aspect ratio increases, the fiber behaves more flexible. Therefore, the fibers adapt more to the flow streamlines and represent a larger curvature compared to rigid fibers with a smaller aspect ratio. A similar discussion as the one made for the effect of fiber shape on the early-stage shape evolution of the fibers is valid here, too, for the early-stage values of ξ obtained (see Sect. 4.2).
After 24 × 10 3 dimensionless time, the value of ξ increases almost linearly with the distance from the center of vortices. As discussed in Sect. 4.1, once the fibers reach their steady-state shapes, they exhibit a fixed curvature across the computational domain, and since the streamline curvature decreases as b increases, the values of ξ increase accordingly.

Conclusions
In this work, simulations were performed to study the dynamics and shape evolution of flexible fibers in cellular flows. The particle-level rod-chain fiber model of Schmid et al. [48] and Switzer III & Klingenberg [49] was used to realistically capture the behavior of flexible fibers in a steady flow, while incorporating the realistic fiber features, including the fiber flexibility and non-straight fiber equilibrium shape. In all simulations, a total of 64 fibers, distributed uniformly in the computational domain, was tracked. The fibers were initially aligned with flow direction with the velocity of the fluid, and their behavior was analyzed for different combinations of fibers stiffness, aspect ratio, and equilibrium shape. In addition, the motion of isolated fibers was carried out to study the fiber transport for several initial conditions (position and orientation) and different fiber stiffnesses.
It was observed that depending on the fiber stiffness, and initial condition, the fibers either became trapped in a vortex of a cell or were freely transported through an array of counter-rotating vortices while exhibiting different modes of fiber deformations, including straight, nearly straight, J-shaped, U-shaped, and C-shaped geometry. These findings indicate the significant impact of initial conditions and mechanical properties of the fiber on the fiber buckling and transportation. The results are in agreement with the numerical and experimental results of Wandersman et al. [17].
The typical fiber relaxation times were long, for example, about 34 × 10 3 dimensionless time for the fiber with (number of segments N , dimensionless stiffness DS, aspect ratio r f , equilibrium bending angle θ eq , equilibrium twisting angle φ eq ) = (11, 4.2 × 10 −3 , 143, 0.01, 0). As expected, increasing the dimensionless stiffness led to a higher restoring bending torque and a more "rigid" behavior, which resulted in a quicker relaxation phase accordingly. Additionally, the findings showed that fibers with higher DS represent a lower ratio of fiber curvature to the flow streamline curvature (ξ = κ f iber /κ f low ) than the more flexible fibers, meaning that they have regained their nearly straight equilibrium shapes to a larger extent.
The results obtained showed a strong influence of the fiber equilibrium shape on the bending behavior of the fibers. For relatively flexible fibers with θ eq = 0 demonstrated a quick early-stage straightening from their initially flow-aligned configuration toward their straight equilibrium shape, and exhibit κ f iber = 0 across the computational domain after 24 × 10 3 dimensionless time. However, even a small deviation in fiber equilibrium shape from perfectly straight shape causes a significant change in the values of ξ . For example, a system containing fibers with the intrinsic radius of curvature R u = 10L (i.e., θ eq ≈ 0.1) exhibits ξ ranging from approximately 0.5 (at the distance from the center of vortices b ≈ 200) to ξ ≈ 0.8 (at b ≈ 1200) when the fibers reach their steady-state shapes.
The simulations performed for fibers with (N , DS, θ eq , φ eq ) = (11, 4.2 × 10 −3 , 0.01, 0) and various aspect ratios showed that fibers with higher r f experience higher ξ . This is to be expected, since for a fixed diameter longer fibers are more flexible than the shorter ones, and consequently as they move along the flow direction, they deform to a larger extent to assume the shape of streamlines.
The work presented herein showed that the shape development of the fiber leads to complex fiber transport dynamics. In further studies, it would be worthwhile to perform a detailed study on the origins of these dynamics, particularly the coupling between the fiber shape deformation and transport. It is recommended to perform further simulations with a wider range of fiber stiffnesses, equilibrium shape angles, and aspect ratios to gain a more detailed understanding of the fiber shape evolution in the studied flow field.