Buckling of elastic fibers in a shear flow

Three-dimensional dynamics of flexible fibers in shear flow are studied numerically, with a qualitative comparison to experiments. Initially, the fibers are straight, with different orientations with respect to the flow. By changing the rotation speed of a shear rheometer, we change the ratio A of bending to shear forces. We observe fibers in the flow-vorticity plane, which gives insight into the motion out of the shear plane. The numerical simulations of moderately flexible fibers show that they rotate along effective Jeffery orbits, and therefore the fiber orientation rapidly becomes very close to the flow-vorticity plane, on average close to the flow direction, and the fiber remains in an almost straight configuration for a long time. This ‘ordering’ of fibers is temporary since they alternately bend and straighten while tumbling. We observe numerically and experimentally that if the fibers are initially in the compressional region of the shear flow, they can undergo compressional buckling, with a pronounced deformation of shape along their whole length during a short time, which is in contrast to the typical local bending that originates over a long time from the fiber ends. We identify differences between local and compressional bending and discuss their competition, which depends on the initial orientation of the fiber and the bending stiffness ratio A. There are two main finding. First, the compressional buckling is limited to a certain small range of the initial orientations, excluding those from the flow-vorticity plane. Second, since fibers straighten in the flow-vorticity plane while tumbling, the compressional buckling is transient—it does not appear for times longer than 1/4 of the Jeffery period. For larger times, bending of fibers is always driven by their ends.

Of particular interest is the compression of flexible fibers caused by the strain of an external fluid flow [10]. Such deformations of an elastic fiber have been investigated in simple shear [13] and compressional/extensional [18,26,35] flows. Recently these studies have been generalized for compression of elastic sheets [42,43] in the same flows; the elastic sheet geometry models, e.g. a graphene platelet [42,43]. The compressional buckling has been predicted theoretically using the Euler-Bernoulli beam model of an elastic slender filament (usually referred to as the elastica) [13], and also for an elastic slender sheet [42,43]. For simple shear flow, such an approximation is limited to very short times, smaller than half of the Jeffery period, because an infinitely thin object does not tumble in a simple shear flow.
In this paper we account for a non-zero thickness of an elastic fiber in a simple shear flow. We use a bead-chain model of a fiber. We apply a multipole expansion of the Stokes equations, corrected for lubrication [44,45], to determine hydrodynamic interactions between the segments of the fiber. The elasticity of the fiber is given by the harmonic stretching and bending potentials between pairs and triplets of consecutive beads [41]. With this theoretical and numerical tool [36,46], and with the use of experimental techniques [4-6, 24, 25] that allow fabrication of fibers of specified aspect ratios, and including flow visualization, we are able to trace the fiber dynamics. We focus on the time scales of the order of a few Jeffery periods. The outline of our paper begins in section 2, where we describe the experimental setup and the process of fiber formation. In section 3, we present the theoretical model. Then, in section 4, we demonstrate experimentally and numerically that fibers undergo a temporal ordering while performing effective Jeffery orbits: for a long time the fibers are almost straight and almost in the flow-vorticity plane; moreover, owing to the nature of the Jeffery orbits (see appendix), there are more fibers inclined at larger angles to the vorticity direction. Then, for a short time, the fibers bend significantly, and the process repeats. In section 5, we detect the compressional buckling in both experiments and numerical simulations, and we investigate its properties. We demonstrate that the compressional buckling happens only once, for time smaller than half of the Jeffery period, and for a very limited range of initial orientations. In section 6 we show that the most common bending pattern is different-bending originates locally at the fiber ends [41]. Finally, conclusions are presented in section 7.

Experimental setup 2.1. Fiber fabrication in a microfluidic channel
We used a standard microfluidic method to fabricate the fibers, where the properties of the fibers, such as the length, the diameter, and the modulus of the fibers, are highly reproducible and controllable [4-8, 24, 25]. In particular, we generated a uniform cylindrical jet of oligomer solution in a microfluidic channel, then exposed the jet to pulses of UV light to trigger the gelation of the jet and produced the fibers; see the schematic of the microfluidic device in figure 1(a).
The microfluidic device was fabricated using standard methods of soft lithography [1]. A polydimethylsiloxane (PDMS, Dow Sylgard 184; Ellsworth Adhesives) channel was plasma-bonded to a PDMS-coated glass slide using a Corona surface treater (Electro-Technic Products). The width and the height of the main channel were W = 200 μm and H = 135 μm, respectively. The oil solution (light blue in the schematic) and oligomer solution (violet in the schematic) were pumped into the microfluidic channel by syringe pumps (Harvard apparatus) [figure 1(a)]. The oil solution was, by weight, 11% Span 80 (Sigma-Aldrich) in light mineral oil (Sigma-Aldrich). The oligomer solution was composed of, by volume, 54% poly(ethylene glycol) diacrylate (PEG-DA, molecular weight = 575 g mol −1 ; Sigma-Aldrich), 38% de-ionized water, and 8% 2-hydroxy-2-methylpropiophenone (photoinitiator; Sigma-Aldrich). The volumetric flow rates of the oligomer and oil solutions were Q 1 = 0.1 ml h −1 and Q 2 = 0.6 ml h −1 , respectively. In the main microfluidic channel, the oligomer solution was focused and sheathed by the two symmetric streams of the oil solution to form a steady and uniform cylindrical jet. The size of the jet was approximately 40 μm, which is much smaller than the width (200 μm) and the height (135 μm) of the main channel. Therefore, the jet is cylindrical due to surface tension. Further, in the previous research [24], the same technique was used to produce the fibers, and their surface was imaged under SEM, and the fibers were indeed cylindrical. The diameter of the jet D ≈ 44.8 μm. Due to the effect of the surfactant (Span 80) in the oil solution, the jet of oligomer solution remained steady for several millimeters downstream in the main channel before it became unsteady and broke up into a series of small droplets (in the absence of UV light discussed in the next paragraph).
To trigger gelation in the jet of oligomer solution, the jet was exposed to pulses of UV light, which were generated by an ultraviolet light-emitting diode (LED) (UV LED, wavelength ≈365 nm; Thorlabs) and focused through a 20× objective. As a result, the UV light spot incident on the jet of oligomer solution was approximately a circle, 2.3 mm in diameter; see the schematic in figure 1(a) (the gray circle denotes the light spot). We set the current of the UV LED to 1.3 A, and the corresponding light intensity of the light spot was approximately 0.03 E/(m 2 s) [7]. The UV LED was pulsed periodically on for 20 ms and then off for 180 ms. Fibers with diameter D = 44.8 μm and length L = 1.61 mm were generated and were then collected downstream at the end of the microfluidic channel. The corresponding aspect ratio of the fibers r ≡ L/D ≈ 36 and was similar to the value used in the numerical simulations. As estimated in the previous research [24], the Young's modulus of the fibers E ≈ 10 5 Pa. The collected fibers were washed five times in 1 wt% Tween 80 aqueous solution (Sigma-Aldrich) and twice in 0.1 wt% Tween 80 aqueous solution, and were diluted and suspended in PEG-DA for further experiments (section 2.2). A more detailed discussion of the effect of different parameters in this setup (such as the period of the pulses of UV light) on the fabricated fibers can be found in reference [7].

Observing fiber dynamics in shear flow using a rheo-microscopy setup
A dilute fiber suspension in PEG-DA was placed in a rheometer (MCR 702; Anton Paar) to study the dynamics of fibers in a shear flow; see the schematic in figure 1(b). The concentration of the fibers in the suspension was kept at approximately 0.5 wt% in all the experiments, corresponding to a solid volume fraction about 0.005. The rheometer was equipped with two transparent glass parallel plates, and the radius of the plates was R = 21.5 mm. To image the fiber suspension between the parallel plates, a white LED was positioned above the glass plates. Light from the LED illuminated the fibers, was then reflected by a mirror below the glass plates, and was detected by a high-speed camera (v7.3; Phantom) mounted with an objective (5×; Mitutoyo) and a corresponding extension tube. The field of view of the images was near the outer edge of the parallel plates, with a size of 2.57 mm × 3.43 mm (600 pixels × 800 pixels). Since the composition of the fibers was similar to the ambient solution (PEG-DA), the refractive index of the fibers was very close to that of PEG-DA, and the fiber suspension was transparent under natural light. In our rheo-microscope setup, to enhance the contrast of the imaging, the LED was a point source, and the camera was set slightly out-of-focus, i.e. the focal plane of the imaging was slightly offset from the desired focal plane, which was in the middle of the fluid in the rheometer. An example of a typical image obtained from the rheo-microscope setup is displayed in figure 1(c). Initially, in the absence of flows, the fibers are straight in elastic equilibrium.
To image the dynamics of the flexible fibers, we applied a constant shear flow between the parallel plates. To track the dynamics of the fibers, we set the parallel plates to be counter-rotating [figure 1(b)], i.e. in the field of view, the top plate moved from the left to the right, and the bottom plate moved in the opposite direction, from the right to the left (the horizontal direction in the field of view was aligned with the velocity direction of the shear flow in the rheometer). The counter-rotating parallel plates created a stagnation plane in the flow, which was in the middle of the cell, and therefore we could track the fibers that were near the middle of the rheometer for a longer time. In the experiments, we set the gap size between the two parallel plates d to be approximately 0.5 or 1.5 mm, and the local average shear rate of the flow in the field of viewγ was, approximately, 375 or 750 s −1 . We note that it takes O(10 −3 ) seconds for the rheometer to accelerate to the desired shear rate. The viscosity of the solution η ≈ 55 mPa s. Therefore, in our experiments, the corresponding dimensionless stiffness ratio A ≡ E/(64ηγ) ≈ 80 or 40, respectively. These values of A in our experiments were similar to the values used in the numerical simulations below. Note that the length scale of the fiber (L = 1.61 mm) was comparable with the gap size between the parallel plates (d ≈ 0.5 or 1.5 mm), so there might be an effect of the confinement of the parallel plates on the dynamics of the fibers. However, in our experiments, we did not observe a significant difference in the dynamics of the fibers between the gap sizes d ≈ 0.5 and 1.5 mm (for a given shear rateγ). More discussion on the effect of confinement and the corresponding drifting of the fibers away from the walls can be found in reference [25].
Since d/R 1, the effective Reynolds number in this geometry is where ρ is the density of the solution. In particular, substituting the characteristic parameters in the experiments, ρ = 1.12 g cm −3 ,γ = 375 s −1 , d = 0.5 mm, η = 55 mPa s, and R = 21.5 mm, the resulting Re eff ≈ 0.02 1. Therefore the effect of inertia on flow in the gap is small. Note that equation (1) is a standard estimate of the effect of the inertia on the flow in a thin layer of the liquid between two parallel plates [47]. Alternatively, it is helpful to estimate the particle-scale Reynolds number by using the scale of the fiber (the diameter of the fiber, D = 44.8 μm in the experiments), and the corresponding Reynolds number which is also small. Brownian motion is irrelevant. Indeed, the Peclet number estimated based on fiber segments of length D is much larger than unity; Pe = 3πηD 3γ /(k B T) = O(10 9 ), where k B is the Boltzmann constant and T is the temperature in Kelvin.
In the recorded image sequence of the experiments, to enhance the contrast (and to remove the noise in the images), we applied a standard background subtraction process, i.e. each image was subtracted from the average intensity of the image sequence and was then re-scaled. Also, to qualitatively compare with the numerical results, we made the time in the experiments dimensionless by t =γτ , where τ is physical time in the experiments, and t is dimensionless time. The dimensionless time t was reported in all the experiments, and t = 0 denoted the time when the shear flow was applied.
Note that under shear, due to the density difference between the fibers and the ambient PEG-DA solution, the fibers migrate toward the center of the rheometer. However, the migration is very slow compared to our time window of interest. For example, for the shear rateγ = 750 s −1 and gap size d = 1.5 mm, the fibers drift approximately 1 mm (a third of the field of view in the y direction) in 2 s, and the corresponding dimensionless time t = 1500. Therefore, the long-time migration of the fibers toward the center of the rheometer does not affect the short-time observations of the fiber dynamics.
Since the fiber density was very close to, but slightly denser than the density of PEG-DA (ρ = 1.12 g cm −3 ), in the absence of shear the fibers sedimented to the bottom plate of the rheometer. The sedimentation of the fibers in the solution was slow. For example, it took tens of minutes for the fibers to sediment to the bottom plate; the corresponding sedimentation velocity v sed = O(10 −6 ) m s −1 and was much smaller than the speed of the shear flow [γd = O(10 −1 ) m s −1 ]. However, under the shear flow there was an opposite tendency: the fibers moved to the central plane of the rheometer, as discussed in [25]. We used this effect to minimize the effect of the confinement in the present study. Namely, prior to recording the fibers in the shear flow with the desired shear rate, we applied a pre-shear to the fiber suspension, i.e. we first applied a high local shear rateγ = 1500 s −1 for 10 s, then we waited for 30 s for the fibers to relax (and straighten), and then started the experiments with the desired shear rateγ and recorded a movie. It is essential that in this way we kept the fibers initially near the middle plane between the two walls due to the pre-shear, and mostly in the x-y plane; the z projection of the fiber length was relatively small. Moreover, the fibers were initially almost straight and their orientations were moderately random, which allowed to compare with the numerical simulations.

Theoretical model and its numerical implementation
The theoretical description of the fiber dynamics is based on the bead model, as in [22,36,41,46,48,49]. Each fiber consists of N identical spherical beads of diameter D. The centers of the consecutive beads are linked by springs of equilibrium length 0 , with the Hookean stretching potential energy, Herek is the spring constant, i = |r i − r i−1 | is the distance between the centers of beads i and i − 1 and the time-dependent position of the center of bead i is denoted as r i , as illustrated in figure 2.
The bending modulusÂ is determined by the Young modulus E, We assume the harmonic bending potential energy where θ i is the angle between the relative positions The total external force F i acting on bead i is the sum of the Hookean and bending forces, At the elastic equilibrium, the fiber is straight. The elastic fiber is immersed in the shear flow of a fluid with viscosity η. We neglect Brownian motion. Hydrodynamic interactions between the fiber beads are obtained from the Stokes equations by the multipole expansion [50,51]. The equations of motion for the positions r i of the beads are [36,41,46] with E ∞ = 1 2 (∇v 0 + (∇v 0 ) T ) denoting the rate-of-strain tensor of the external shear flow v 0 given by equation (7). The translational-translational and translational-dipolar mobility matrices, μ tt ij and μ td ij , respectively, are functions of positions of all the beads. Here, they are evaluated numerically by the Hydromultipole numerical program, based on the multipole expansion corrected for lubrication [44,45,52] to speed up the convergence and provide a high precision [44,45]. The equations of motion (8) are solved numerically by the fifth-order Runge-Kutta procedure.
We proceed with a dimensionless representation. Lengths are normalized by the bead diameter D and time by the inverse shear rate 1/γ. Values of the dimensionless parameters used in the simulations are: number of beads N = 40, distance 0 /D = 1.02 between the consecutive bead centers at the elastic equilibrium, aspect ratio r = (N − 1) 0 /D + 1 ≈ 40. We assume that the dimensionless elastic resistance ratio k =k/(πηdγ) is large, k = 1000, which leads to an almost constant fiber length. We perform computations for more than 300 values of the dimensionless bending stiffness ratio mostly for 7 A 200.

Initial conditions and early times
In the numerical simulations, the 3D dynamics of a single fiber are evaluated, in contrast to the previous 2D analysis [41,46]. Here we focus on short-time transient effects rather than studying the long-time approach to attracting modes of the dynamics as in [41]. The fiber is initially straight and at elastic equilibrium. The initial orientations are parameterized by the spherical angles Θ 0 and Φ 0 as shown in figure 2(a). The time evolution is determined for the whole range of the angles 0 Θ 0 π and 0 Φ 0 π with a grid of 5 • . Orientations (π + Φ 0 , π − Θ 0 ) are equivalent to (Φ 0 , Θ 0 ). Initially, in the experiments, fibers are straight with different, approximately random, orientations, as illustrated in figure 2(b). The evolution of a dilute suspension is observed in the xy (flow-vorticity) projection. When exposed to the shear flow, bending of the fiber ends slowly develops, with a characteristic bending time τ b that is typically long in comparison to the tumbling time (half-period of the fiber rotation) [41]. At early times t τ b , most of the fibers are almost straight and quickly orient close to the flow direction, as illustrated in figure 3, after which bending occurs. In the next section, we will trace the fiber evolution for longer times.  (7) is parallel to the x axis and has a gradient parallel to the z axis. In the numerical simulations the initial orientation of a single fiber is given by the spherical angles Θ 0 and Φ 0 . We determine the three-dimensional evolution of the fiber shape for the whole range of Θ 0 and Φ 0 .

Longer times
In the numerical simulations, we observe that for longer times, the fiber regularly bends and straightens out while rotating, with a certain characteristic tumbling time that can be interpreted as the dimensionless half-period T J /2 of an effective Jeffery orbit [36,46], with where s is a fiber 'effective' hydrodynamic aspect ratio (typically close to the geometrical value, N). For most initial orientations studied in this work, the period T J weakly depends on the bending stiffness ratio A, with T J = 180, 225, 252 for A = 10, 80, 160, respectively, as determined here from the periodic evolution of the fiber curvature versus time. In figure 3 we illustrate that a typical bending time τ b is a large fraction of T J ; τ b is larger than T J /4 and smaller than T J /2. Indeed, in figure 3 for time t 70 it is hard to see much bending at all-the bending time is larger than 70 while the Jeffery period is 225. In experiments and numerical simulations we observe that for most of the 3D initial orientations bending starts locally from the fiber ends, by analogy to previously reported experiments and numerical simulations of fibers initially aligned with the flow [27,29,41]. We first focus on our numerical simulations. To quantify how much the fibers are bent, the local curvature κ j+1 (t) is calculated as the inverse of the radius of a circle circumscribed on the centers of three consecutive beads located at r j (t), r j+1 (t) and r j+2 (t); the positions are normalized by the bead diameter 2a. The highest possible value of the local curvature in this model is √ 3, which corresponds to all three beads touching each other. The local curvature κ j+1 (t) at the (j + 1)th bead and average curvatureκ(t) averaged along the entire length of the To study evolution of the fiber orientation, we calculate the fiber inertia tensorÎ in its center-of-mass frame,Î where α = x, y, z, β = x, y, z and r j = (r jx , r jy , r jz ) is the position of jth bead in the center-of-mass reference frame. Then, we evaluate the eigenvalues I 1 , I 2 , I 3 ofÎ, and the corresponding normalized eigenvectors n 1 , n 2 , n 3 . Without loss of generality let I 1 > I 2 > I 3 . Then, the principal axis n 3 (t) is described by the spherical angles Θ(t) and Φ(t), in analogy to Θ 0 = Θ(0) and Φ 0 = Φ(0) shown in figure 2(a).
To analyze the main features of the evolution of the fiber shape, we evaluate numerically the dynamics for the ensemble of the initial orientations of straight fibers, within the full range of 0 Θ 0 π and 0 Φ 0 < 2π, with a grid π/36 = 5 • for both angles. For each case, we evaluate the ensemble-averaged, time-dependent mean curvature κ(t) and angle Φ(t) , defined by the following equations, whereκ(t|Θ p , Φ q ) is the average curvature at time t for a fiber that at time t = 0 was straight at the orientation Θ 0 = Θ p , Φ 0 = Φ q , and Θ p = p · 5 • , Φ q = q · 5 • with p = 0, . . .   In the experiments, it is observed in the flow-vorticity plane that at certain times, almost all of the fibers are straight and with orientations close to the flow direction. Then, almost all fibers are bent, and such a temporary configuration repeats in time, as illustrated in figure 5.

Compressional buckling
Some of the fibers that are initially straight and oriented so that they experience compressive stress from the flow, after a short time buckle globally along the whole length, with a large average curvature; this response can be predicted with the use of the elastica model in [13]. The compressional buckling happens only once, as illustrated in figure 6. Later, the fibers straighten out in the flow-vorticity plane, and slowly develop local bending at their ends, similarly as in the case of Θ 0 = 90 • described in [41]. Characteristic shapes of fibers in the compressional buckling and local bending modes differ from each other, and the maximum of the average curvature is much larger for the compressional buckling.
The compressional buckling can happen only for certain initial orientations Θ 0 and Φ 0 of the fiber, indicated in figures 7(a) and (b) as the regions of the light color (yellow online). For A = 10 and A = 160, the average curvatureκ(t) of the compressionally buckled shapes reaches a maximum exceeding 0.05 for short times t 20, and we used this criterion to determine the range of the initial angles Θ 0 and Φ 0 where the compressional buckling can take place. For Φ 0 175 • after times t > 20, the fibers are already out of the compressional regime of the flow. Therefore, for the dark (green online) regions of the initial orientation angles shown in figure 7, compressional buckling of a sufficiently large amplitude is not  observed. Also, for large Φ 0 , equal or very close to 180 • , the compressional buckling does not occur, which is indicated as a thin dark (green online) vertical stripe in figure 7. This finding explains that compressional buckling of moderately elastic fibers happens only once because such fibers straighten out again only while approaching the orientation with Φ = 0. Figure 7 illustrates that for a larger value of the bending stiffness A, the compressional buckling occurs for a smaller region of the initial orientations.
Compressional buckling is also observed in experiments. In the experiments, the fibers are initially near the middle plane between the two walls due to the pre-shear, and are mostly in the xy plane. The z projection of the fiber length is relatively small. Also, the buckling behavior is very rare in the experiments. We estimate that the ratio of the fibers that buckle is of the order of 1%.
Qualitative comparisons of experiments and simulations for similar parameter ranges are shown in figures 8 and 9. In figure 9, the influence of the compressional flow is stronger and leads to more curved shapes than in figure 8. In the numerical simulations shown in figures 9(b), there is a very small (of the order of 10 −6 ) initial random perturbation i of the position r i of every bead i in the straight fiber.
In figures 10 and 11 we perform more detailed analysis of the time-dependent xy projection of the fiber shapes for the experiments and simulations shown in figures 8 and 9, respectively. We measure length L ee of the xy projection of the end-to-end vector (as shown in figure 10(a)) of the buckling fiber as a function of time. In figures 10(b)-(e) and 11(a)-(d) this length L ee is normalized by the length of the fiber in a static fluid, L. We also measured the xy projection of the 'width' of the buckling fiber, δ, as illustrated in figure 10(a). The width δ is normalized by the diameter D of the fiber, and plotted in figures 10 and 11.  Also, in the xy projection we determined the angle α of the fiber end-to-end vector to the flow direction, as illustrated in figure 10(a). With time, the fiber rotates along an effective Jeffery orbit, and the angle α decreases. The rapid change of α takes place when the fiber orientation is in the compressional regime of the ambient flow. In this range, the buckling is observed as a sharp minimum of L ee and a sharp maximum of δ.
Our numerical simulations indicate that the shapes of the fibers during the compressional buckling significantly depend on the bending stiffness ratios A. The typical pattern of triggering higher modes of the deformed shapes for decreasing values of the bending stiffness ratio A is illustrated in figure 12. The fiber curvature grows while the fiber moves in the compressional region of the shear flow, and reaches a maximum for the angle Φ typically a little larger than 90 • , as shown in the left panels of figure 12. (In general, the beads are out of the plane of the figure and they can partially cover each other on the xy and xz projections shown there.)

Differences between compressional buckling and local bending
In the numerical simulations we observe two competing scenarios: compressional buckling and local bending. For a (moderate) fixed value of A, and a fixed, large enough value of the initial polar angle Θ 0 , a compressional mode is observed in a certain range of the initial azimuthal angles Φ 0 . For larger values of Φ 0 , very close to 180 • , the fiber has enough time to develop the local bending mode [41]. The main differences between the modes are illustrated in figures 13(a)-(c), (d)-(f). In the compressional buckling mode, the fiber deforms along the whole length; as the result, the principal axes of the fiber and of its central part have almost the same orientations during the entire bending process. Moreover, positions i of the maximum local curvature κ i do not change with time. In contrast, the local bending originates from the fiber ends and, with increasing time, the maximum local curvature moves toward the middle part of the fiber, as visible in figure 13(f) and discussed previously in [27,29]. During this process, orientations of the principal axes of the whole fiber and of its central part differ from each other significantly, as shown in figures 13(d) and (e). The above differences are visible also in figure 14, where typical evolutions of the fiber shapes showing the compressional buckling and local bending modes are compared with each other.
For larger values of Φ 0 , the compressional buckling and local bending modes coexist, as illustrated in figure 15, and their characteristic features can be identified. Even in the local bending evolution shown in figures 13(d)-(f) and 14(b) a certain (small) influence of the compressional bucking mode can be expected. The experimental and numerical results of this work indicate that for 3D dynamics of flexible fibers of a moderate length and a small bending stiffness ratio in shear flow, the dominant physical mechanism of bending is the local bending originating from the deformation of the fiber ends [41]. In contrast, the   compressional buckling [13] is a transient effect limited to short times and a narrow class of initial conditions.

Conclusions
Short-time dynamics of initially straight fibers in shear flow have been studied numerically and experimentally. We focused on fibers with the bending stiffness ratio A and the aspect ratio N large enough to allow fibers, after bending or buckling, to straighten out at later times. The experimental technique allowed the observation of fibers in the flow-vorticity plane. With this tool, some features of 3D dynamics of fibers that move or deform out of the shear plane have been measured. In particular, time-dependent ordering of the fiber orientations has been observed in experiments and simulations. The physical mechanism of this effect is that flexible fibers follow effective Jeffery orbits and, therefore most of the time remain almost straight and close to the flow-vorticity plane, with the angle Φ close to zero, as illustrated in figure 4 for an ensemble of initial orientations studied numerically. Moreover, owing to the nature of the Jeffery orbits (see appendix), for Φ ≈ 0 there are more fibers inclined at larger angles Θ to the vorticity direction.
Different fiber bending patterns have been observed in experiments and simulations, depending on the initial orientation of a straight fiber. For some initial orientations, fibers later undergo a strong compression along the whole length, followed by a decompression, as shown in figures 8 and 9. We have called this response the compressional buckling mode. For other initial orientations, this effect is not seen. The physical mechanism of the compressional buckling is the same as for the elastica, as described for motion in the shear plane by Becker and Shelley [13] for a small perturbation from a straight configuration. We have observed that this mode is very rapid (it lasts only a small fraction of the Jeffery period), but it is characterized by a significant deformation. These features are seen in figure 6 as a sharp peak of the numerically computed mean fiber curvature. We have shown numerically and experimentally that the compressional buckling happens only once at the very beginning of the evolution, after which it does not reoccur. In the numerical simulations, this transient mode can be seen only for a limited range of larger initial inclination angles Θ 0 and Φ 0 , excluding Φ 0 = 180 • (and Φ 0 slightly smaller than 180 • ), as shown in figure 7 as the region marked by the light color (yellow online). For larger values of the bending stiffness ratio A there are fewer initial orientations that allow for development of the short-time compressional buckling. It is worth noting that all the fibers that straighten out while tumbling reach the angle Φ = 0 • (the same as Φ 0 = 180 • ). Therefore the compressional buckling is a transient effect, limited to times smaller than one fourth of the effective Jeffery period.
It is important to note that the theoretical model takes into account the non-zero thickness of the fiber, hydrodynamic interactions between all the fiber beads, and the coupling of the local velocity of each bead with the rate-of-strain tensor of the ambient flow, as described in equation (8). As the result, the fiber reaches Φ = 0 in a finite time, and then tumbles, in contrast to the elastica model, for which time to approach Φ = 0 is infinite. Therefore our model allows us to determine that the buckling does not reappear after reaching Φ = 0.
For smaller initial orientation angles Θ 0 and Φ 0 , fibers for shorter times practically do not deform-they rapidly reach Φ close to 0 • while performing their Jeffery orbits. For the initial orientation angles Φ 0 180 • the fibers bend locally by a different physical mechanism [41]: the fiber bending slowly develops at the fiber ends and the maximum mean curvatureκ is smaller than in the compressional buckling mode, as visible in the numerical simulations used to construct figure 6, which shows that the local bending repeats in time periodically. These two mechanisms of bending compete with each other, and often coexist, depending on the ratio of the characteristic timescales, e.g. as illustrated in figures 13-15. Depending on the bending stiffness ratio A, the compressionally buckled mode has different but regular sinusoidal-like shapes. Smaller A yields a larger number of bends along the fiber (the smaller is the characteristic 'wavelength'), as depicted in the numerical simulations in figure 12. The wavelength also depends on the initial orientation angles and can increase with time.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix. Time-dependent orientation of the fiber principal axis
Time-dependent orientation (Θ, Φ) of the fiber principal axis n 3 (t) follows an effective Jeffery orbit, as illustrated by two examples in figure 16. Initially, the fiber is straight, with the angle Φ = 170 • . The angle Φ(t) of the fiber principal axis systematically decreases with time. The compression of the ambient shear flow acts until Φ(t) reaches 90 • , with the maximum at Φ = 135 • , It is interesting to observe numerically that deformation of the fiber still increases with time even when the fiber rotates to angles Φ smaller than 135 • , and therefore to smaller compression of the ambient flow. As shown in figure 16, the maximum of the fiber curvature takes place for an angle Φ smaller than 135 • and only a bit larger than 90 • where the compression of the ambient flow vanishes.
This pattern is observed for all initial angles Θ 0 > 0. Fibers out of the shear plane deform in both xz and xy planes. However, for smaller Θ 0 , the fiber feels a weaker compression of the ambient flow, and therefore the amplitude of its buckling decreases with the decrease of Θ 0 , as visible by comparing figures 8 and 9.