2 Mathematical model

We consider an elastic flexible fibre immersed in a four-roll mill background flow (Young & Shelley Reference Young and Shelley2007; Wandersman et al. Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010). In the absence of a fibre or viscoelasticity, the velocity field of this cellular, doubly periodic flow on the domain $\tilde{\unicode[STIX]{x1D6FA}}=[-W,W]\times [-W,W]$ is

This background velocity may be generated by adding a time-independent external force density to the Navier–Stokes equations. Figure 1(a) shows this vortical background flow, where vertical compression and horizontal stretching occur near the central stagnation point. Within this domain, we will introduce a fibre of length $L$ . We choose the physical parameters as shown in table 1. These parameters largely follow the experiments of Wandersman et al. (Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010), except that here we choose a cellular length of $W=\unicode[STIX]{x03C0}~\text{cm}$ while their experiments were performed with $W=3~\text{cm}$ .

Figure 1. Fibre configurations at different times at $Re=2.5$ , $Wi=1$ and $\unicode[STIX]{x1D702}=370$ . This fairly stiff fibre approaches the stagnation point at the centre, bends, regains its straight configuration and then moves away from the stagnation point along the outgoing streamline to the right.

We choose our dimensional length scale to be the vortical cell diameter $W$ , and the time scale to be the inverse of the maximal extension rate $\dot{\unicode[STIX]{x1D6FE}}=\unicode[STIX]{x03C0}u_{0}/W$ . We assume that the viscoelastic fluid dynamics is described by the coupled non-dimensional Navier–Stokes/Fene-P equations in the domain $\unicode[STIX]{x1D6FA}=[-1,1]\times [-1,1]$ with periodic boundary conditions:

The non-dimensional velocity and pressure are $\boldsymbol{u}=(u,v)$ and $p$ , respectively, and the Reynolds number is $Re=\unicode[STIX]{x1D70C}\unicode[STIX]{x03C0}Wu_{0}/\unicode[STIX]{x1D707}$ . The Weissenberg number, $Wi=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{f}$ , is the ratio of the polymer relaxation time $\unicode[STIX]{x1D70F}_{p}$ and the time scale of the flow $\unicode[STIX]{x1D70F}_{f}=(\dot{\unicode[STIX]{x1D6FE}})^{-1}$ . The dimensionless stress tensor $\unicode[STIX]{x1D63E}$ has been scaled by $G$ , the isotropic stress arising from thermodynamic fluctuations of the polymer chains in the absence of flow. Note that $G=\unicode[STIX]{x1D705}TV$ , where $\unicode[STIX]{x1D705}$ is the Boltzmann constant, $T$ is temperature and $V$ is the number of polymer chains per unit volume (Larson Reference Larson1998). The Fene-P stress tensor $\unicode[STIX]{x1D64E}$ penalizes the extensibility of the infinitesimal polymer chains, and ${\mathcal{L}}$ represents the ratio of the maximal polymer length when stretched to the polymer length when coiled. Here $\unicode[STIX]{x1D64E}$ represents the extra stress generated by the transport and distension of the polymer field, and is coupled to the Navier–Stokes equation with the parameter $\unicode[STIX]{x1D6FD}=G\unicode[STIX]{x1D70F}_{f}/\unicode[STIX]{x1D707}$ . Note that $\unicode[STIX]{x1D6FD}\cdot Wi$ is a physical quantity that measures the ratio of polymer viscosity to solvent viscosity. As in Thomases & Shelley (Reference Thomases and Shelley2007), here we choose ${\mathcal{L}}^{2}=50$ and $\unicode[STIX]{x1D6FD}Wi=0.5$ .

In the momentum equation (2.2), $\boldsymbol{g}(\boldsymbol{x})$ is the external force density chosen to achieve the background velocity field,

and $\boldsymbol{f}$ is the force density imparted to the fluid by the fibre,

Here the position of the fibre $D$ is given by $\boldsymbol{X}(s,t)$ , with $s$ a material coordinate and $\unicode[STIX]{x1D6FF}$ the two-dimensional Dirac delta function. The fibre is taken to be a generalized Euler elastica whose equilibrium configuration is straight. The fibre supports forces $\boldsymbol{F}=-(\unicode[STIX]{x2202}E/\unicode[STIX]{x2202}\boldsymbol{X})$ that are derived from tensile ( $E_{t}$ ) and bending ( $E_{b}$ ) energies,

where $\unicode[STIX]{x1D70E}_{s}$ is a non-dimensional tensile stiffness and $\unicode[STIX]{x1D70E}_{b}$ is a non-dimensional bending stiffness. In the simulations presented below, we choose $\unicode[STIX]{x1D70E}_{s}$ large enough so that the fibre remains nearly inextensible, with deviations from its equilibrium arclength of no more than 3 %. The non-dimensional bending stiffness $\unicode[STIX]{x1D70E}_{b}$ may be varied to change fibre flexibility, and is related to the dimensional bending rigidity $EI$ of the fibre by

The relative importance of flow forces to elastic forces is measured by

where $c=-\!\log (\unicode[STIX]{x1D716}^{2}e)$ . As in Quennouz et al. (Reference Quennouz, Shelley, du Roure and Lindner2015), we refer to $\unicode[STIX]{x1D702}$ as the elasto-viscous number. This non-dimensional control parameter has also been referred to as ‘effective viscosity’ in Tornberg & Shelley (Reference Tornberg and Shelley2007), ‘effective flow forcing’ in Nguyen & Fauci (Reference Nguyen and Fauci2014) and ‘sperm number’ in Wandersman et al. (Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010). In the simulations presented below, we choose the bending stiffness $\unicode[STIX]{x1D70E}_{b}$ so that the range of the elasto-viscous number $\unicode[STIX]{x1D702}$ is as stated in table 1.

The system is closed by enforcing the no-slip condition that requires the material points of the immersed fibre to move at the fluid velocity defined at those points:

The formulation of the coupled fluid–fibre–viscoelastic system described above is readily discretized using an immersed boundary formulation, as in Teran, Fauci & Shelley (Reference Teran, Fauci and Shelley2010), Chrispell, Fauci & Shelley (Reference Chrispell, Fauci and Shelley2013) and Thomases & Guy (Reference Thomases and Guy2014). The fluid velocity, pressure and polymer stress tensors are defined on an $M\times M$ finite difference grid, and the centreline of the macroscopic fibre is discretized using $N$ points. The tensile and bending forces supported on the Lagrangian fibre are spread to the grid using discretized Dirac delta functions (Peskin Reference Peskin2002). Within this immersed boundary framework, we associate an effective radius of the fibre to be one grid spacing, which accounts for most of the support of the forces along the one-dimensional fibre centreline. A time-splitting method is used to evolve the momentum equation and the polymer stress equation. The momentum equation is updated using the formally second-order method of Lai & Peskin (Reference Lai and Peskin2000), and the polymer stress equation is updated using the square-root method of Balci et al. (Reference Balci, Thomases, Renardy and Doering2011). Both velocity and polymer stress are evolved using a second-order Runge–Kutta method. In the simulations presented below, we discretize the fluid domain using a lattice of $256\times 256$ grid points and the fibre, whose length is not varied, is discretized as a one-dimensional curve with $N=103$ points. Because of the properties of the discretized Dirac delta function, which the spreads the fibre forces to nearby grid points, the effective radius of the fibre is taken to be a single grid spacing. This interpretation is used to identify the elasto-viscous number of the coupled system.

3 A case study

Here we examine the motion of an elastic macroscopic fibre immersed in a polymeric fluid that is driven by external forcing. The Reynolds number $Re=2.5$ is small, but non-zero, the elasto-viscous number is $\unicode[STIX]{x1D702}=370$ and the Weissenberg number $Wi=1$ . We interpret the effective radius of the one-dimensional fibre to be the grid spacing, giving an aspect ratio $\unicode[STIX]{x1D716}=r/L\approx 0.02$ . This aspect ratio is at the high end of those reported in Quennouz et al. (Reference Quennouz, Shelley, du Roure and Lindner2015) and approximately three times the aspect ratio of the fibres used in the experiments of Wandersman et al. (Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010). At time $t=0$ the fibre is placed vertically in the domain, offset slightly from the midline (along $x=1.2\times 10^{-3}$ with its lowest point at $y=1.2\times 10^{-6}$ ). Initially, the fibre is in its equilibrium, straight and unstretched configuration (figure 1 a). The velocity field was initialized to be the cellular vortical flow given by (2.1), and the polymer stress tensor $\unicode[STIX]{x1D64E}$ was initialized to be the identity. As the fibre nears the stagnation point along the stable manifold, it buckles slightly (figure 1 b). At time $t=7$ we see that the fibre is nearly straight again and follows a path away from the stagnation point along the horizontal unstable manifold (figure 1 c,d). In related computational studies that examined motion of an inextensible fibre in a Newtonian–Stokes cellular flow, Young & Shelley (Reference Young and Shelley2007) found that the onset of a buckling instability occurred for elasto-viscous number $\unicode[STIX]{x1D702}\approx 328$ , and later experimental studies by Wandersman et al. (Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010) found this threshold value to be $\unicode[STIX]{x1D702}\approx 400$ . The elasto-viscous number of the fairly rigid fibre in figure 1 is near these buckling thresholds, and, while demonstrating some bending, does not exhibit dramatic shape changes.

While the Reynolds number chosen here $Re=2.5$ is in the range of the experiments in Wandersman et al. (Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010), the experiments of Quennouz et al. (Reference Quennouz, Shelley, du Roure and Lindner2015) were performed at lower values. We did run simulations as in figure 1 at the lower Reynolds numbers $Re=0.5$ and $1.0$ and found little difference in the fibre dynamics.

Figure 2. (a) Contour plot of the difference in vertical velocity components at $t=7$ between the coupled fibre–viscoelastic system shown in figure 1(c) and a simulation with no immersed fibre, normalized by the vertical velocity of the fibre-free case. (b) The relative L2 norm of the difference between the velocities (both horizontal and vertical components) of the coupled fibre–viscoelastic system and the fibre-free case for time $0\leqslant t\leqslant 15$ .

The presence of the fibre influences the local flow structure. In order to measure this influence, we ran a simulation using the same parameters as in the present case study, but without the immersed fibre. Figure 2(a) shows a contour plot of the difference in vertical velocity components at $t=7$ between the coupled fibre–viscoelastic system shown in figure 1(c) and the simulation with no immersed fibre, normalized by the vertical velocity of the fibre-free case. We see that the differences are localized near the fibre, and the vertical velocities near the fibre endpoints are altered nearly 25 %. Figure 2(b) shows the L2 norm of the difference between the two components of velocity, normalized by the L2 norm of velocity in the fibre-free case as a function of time. Interestingly, the fibre has the least effect on the velocity field at times that correspond to buckling events (near time $t=6$ , for instance).

Figure 3. Snapshots of contours of the $\log (\text{tr}\,\unicode[STIX]{x1D64E})$ and the macroscopic fibre. Note that $\text{tr}\,\unicode[STIX]{x1D64E}$ measures the mean squared distension of polymer coils. Here $Re=2.5$ , $Wi=1$ and $\unicode[STIX]{x1D702}=370$ .

Figure 3 shows contour plots of the log of the trace of the stress field $\unicode[STIX]{x1D64E}$ . Note that $\text{tr}(\unicode[STIX]{x1D64E})$ measures the mean squared distension of polymer coils. As demonstrated by Thomases & Shelley (Reference Thomases and Shelley2007) for cellular viscoelastic flows without an immersed macroscopic fibre, stress concentrates on a region emanating from the stagnation point in the extensional direction. For the small Weissenberg number $Wi=1$ , the polymers are slowly stretched and this stretch is limited by the Fene-P penalization. In fact, in all of the simulations presented here, the Weissenberg numbers are small enough and the runtime is small enough that the threshold of instability of the cellular flow is not reached. Figure 3(a) shows that when the fibre is almost straight and is being compressed, large polymer stresses are concentrated near its ends. At this instance, the fibre’s tensile energy is large as it tries to maintain its length (figure 7 a). The contours in figure 3(b) show that as the fibre starts to bend the polymer stresses decrease near its ends, and when it straightens in figure 3(c) the polymer stress is distributed along the fibre. As the fibre moves along the unstable horizontal direction, we again see large polymer stretching near its ends.

Figure 4. (a) At time $t=7$ , contours of the difference between $\text{tr}(\unicode[STIX]{x1D64E})$ in the fibre–viscoelastic system and that in the corresponding fibre-free viscoelastic system. (b) Values of $\text{tr}(\unicode[STIX]{x1D64E})$ at the stagnation point in both the fibre–viscoelastic system and the fibre-free viscoelastic system as a function of time. Here $Re=2.5$ , $Wi=1$ and $\unicode[STIX]{x1D702}=370$ .

In order to determine how the presence of the macroscopic fibre affects the polymer stress field, we examine the difference between $\text{tr}(\unicode[STIX]{x1D64E})$ with a fibre and $\text{tr}(\unicode[STIX]{x1D64E})$ in a simulation with no fibre. Figure 4(a) shows contours of this difference at time $t=7$ . We see positive values around the ends of the fibre, which demonstrate that there is extra distension of the polymers caused by the fibre locally at its ends. However, along the extensional horizontal direction around the stagnation point, negative contours indicate that the presence of the fibre limits the polymer distension there. Figure 4(b) plots the evolution of $\text{tr}(\unicode[STIX]{x1D64E})$ evaluated at the stagnation point for both the coupled fibre–viscoelastic system and the fibre-free system. The dotted blue line, indicating $\text{tr}(\unicode[STIX]{x1D64E})$ in the absence of a fibre, shows that this value stabilizes to a constant, as the Fene-P penalization limits the polymer stretching. The solid red curve indicates $\text{tr}(\unicode[STIX]{x1D64E})$ at the stagnation point for the fibre–viscoelastic system. Note that up until approximately $t=12$ , distension of the polymers at the stagnation point falls considerably below that in the fibre-free system. As the fibre moves away from the stagnation point, it carries along stretched polymers at its end to the fixed location of the stagnation point, and $\text{tr}(\unicode[STIX]{x1D64E})$ spikes. As the fibre exits the cell, the distension of polymers at the stagnation point relaxes to the value in the fibre-free case.

Figure 5. (a) Fibre and surrounding stress ellipse field at time $t=12$ . (b) Horizontal component of fluid velocity along $x$ -axis at time $t=12$ .

We may also visualize the symmetric, positive definite stress field $\unicode[STIX]{x1D64E}$ using ellipses whose axes align with eigenvectors of $\unicode[STIX]{x1D64E}$ and whose lengths are the two positive real eigenvalues. The major axis tells us how much and in what direction the polymers are stretched locally. Initially the stress fields are isotropic, and the added polymer stress ellipses are unit circles. Figure 5(a) shows the fibre at $t=12$ and the corresponding stress ellipses. Figure 5(b) shows the horizontal component of the fluid velocity along the $x$ -axis of the domain at the corresponding time $t=12$ . We see that this velocity is approximately constant in the interval of the axis that is occupied by the fibre at that time, indicating that it is moving uniformly to the right.

When this coupled fibre–viscoelastic system is run until time $t=80$ , the fibre meanders through the periodic array of vortical cells, and meets seven stagnation points. Figure 8(b) shows the trajectory of the centre of mass of the fibre. As it approaches a stagnation point along a stable manifold, it bends, samples the local velocity field near this hyperbolic point, and exits along an unstable manifold. The direction along which the fibre exits is not predictable, as slight differences in the bending of the fibre can cause it to move in either direction. We quantify the buckling behaviour of the fibre as it moves through the domain in two ways. One measure, that is easiest to capture in physical experiments, is the distance between the fibre’s endpoints, normalized by its equilibrium length. For a straight fibre, this scaled end-to-end length is one. In the context of our computational model, we can also measure fibre buckling by calculating its bending energy $E_{b}$ in (2.8). Figure 6 shows these two measures of the fibre’s bending over time. We can see that the fibre behaves similarly when it gets close to different stagnation points. As it moves towards the stagnation point, the scaled end-to-end length (figure 6 a) abruptly decreases but moves back to one as the fibre straightens. Figure 6(b) shows similar behaviour with bending energy abruptly increasing as the fibre meets the stagnation points. Both measures of buckling show that there is a nearly constant transit time for the fibre to get from one stagnation point to another, but that the evolving shape deformation of the fibre is not periodic.

The time evolution of the bending energy in figure 6(b) reveals a curious, repeatable, ‘blip’ during the straightening phase. The bending energy spikes as the fibre reaches its maximum bent shape near the stagnation point, quickly decreases as the fibre straightens, but then increases a small amount until decreasing again. We now focus on the fibre’s bending energy during its first trip in figure 7(a). Figure 7(b) shows the fibre configurations corresponding to the first peak of bending energy ( $t=6.1$ ), then at the smallest value of bending energy right before the blip ( $t=7.2$ ), and then at the second, smaller peak ( $t=7.9$ ). We see the fibre transition from its most bent configuration, straighten out, but then continue to bend in the recovery direction. One possibility for this overshoot in bending energy is the memory in the system due to both the inertia in the Navier–Stokes momentum conservation equation as well as the evolution of the stresses carried by the embedded polymers. Another possibility is that this blip is picking up the second peak in an exponentially damped sinusoidal solution of the classic elastic beam equation.

We note that in this immersed boundary framework, the fibre is not precisely inextensible, but imparts tensile forces to maintain its length. Although we choose a tensile coefficient large enough so that the fibre’s arclength does not deviate more than 3 % of its equilibrium length, its realized arclength does vary. Figure 7(a) also shows the evolution of the tensile energy $E_{t}$ (2.8) during its first trip. When the fibre approaches the stagnation point it is almost straight, but the compression due to the flow generates significant tensile energy. We see that, during the time interval when the fibre is most bent, the tensile energy is at its lowest – the fibre struggles less to maintain its length during buckling.

Figure 6. (a) The scaled end-to-end length of the fibre over time. (b) The bending energy of the fibre over time. Note that bending occurs periodically, but the details of the shape deformation are not periodic.

Figure 7. (a) Evolution of fibre bending energy and tensile energy during its first trip for time $0\leqslant t\leqslant 15$ . (b) Snapshots of fibre configurations during first trip that correspond to the first peak in bending energy, an intermediate time point, and then the second small peak in bending energy. Note that there is some overshoot of the straight equilibrium at the latest time point.

Figure 8. Trajectories of the centre of mass of the fibres for various $Wi$ : (a) Newtonian, (b)  $Wi=1$ , (c)  $Wi=2$ , (d)  $Wi=3$ , (e)  $Wi=4$ and (f)  $Wi=5$ . Here $Re=2.5$ and $\unicode[STIX]{x1D702}=370$ . Open circles indicate where they start and where they stop.

4 Effect of Weissenberg number $Wi$

Here we examine the effects of varying the Weissenberg number of the coupled fibre–viscoelastic system while keeping the same Reynolds number ( $Re=2.5$ ) and elasto-viscous number of the fibre ( $\unicode[STIX]{x1D702}=370$ ) as in the previous example. In the absence of a macroscopic fibre, Thomases & Shelley (Reference Thomases and Shelley2007) showed that the strain rate at the stagnation point at steady state in the viscoelastic cellular flow decreased monotonically with $Wi$ for the same background forcing. Since this strain rate governs the velocity time scales achieved in the cellular flow, we expect that fibre transport will slow down as viscoelasticity increases. We compare six simulations with Weissenberg number $Wi=0$ , 1, 2, 3, 4 and 5, noting that $Wi=0$ corresponds to a Newtonian fluid. Initially each fibre was placed vertically, close to the centre stagnation point, but slightly off the incoming streamline as in figure 1(a). At each $Wi$ , the qualitative motion of the fibres as they approached and moved away from the stagnation point was the same as that for the $Wi=1$ example above. The fibre in the Newtonian fluid bends and moves away from the stagnation point first, followed monotonically by the fibres at increasing Weissenberg numbers. This monotonic arrival at stagnation points is also consistent with the observation by Thomases & Shelley (Reference Thomases and Shelley2007) that steady shear rates decrease with $Wi$ . Figure 8 shows the trajectories of the centre of mass of each fibre at different $Wi$ up to time $t=80$ , expanded to include many cells in the periodic domain. Open circles indicate both initial and terminal positions. While each fibre moves away from the first stagnation point it meets along the same outgoing streamline, their trajectories diverge after meeting the second stagnation point. As discussed by Young & Shelley (Reference Young and Shelley2007), the outgoing streamline along which the fibre leaves the area of the stagnation point is chosen with apparent randomness. Small fluctuations in the fluid environment and small changes in fibre bending easily alter the direction (Quennouz et al. Reference Quennouz, Shelley, du Roure and Lindner2015). Figure 8 indicates that, in the given time interval, the fibres in the Newtonian and smaller $Wi$ flows are able to sample more of the spatial domain than the fibres in the higher $Wi$ flows. In particular, the fibres at $Wi=4$ and $5$ are trapped in the cells where they were initially launched. We will revisit this feature of trapped fibres below.

The trajectories in figure 8 do not give information about the time scale of fibre transport. Figure 9(a) shows the moment when the fibre arrives at successive stagnation points. This is measured by monitoring the distance between the fibre’s head and the stagnation points of the cellular domain. We see that for the first trip, the Newtonian one arrives at a stagnation point first, with the fibres at $Wi=1$ , 2, 3, 4 and 5 monotonically following. However, the evolving fibre geometries and deviation in chosen paths in the different $Wi$ cases cause this monotonicity to be destroyed at later stagnation points. Figure 9(b) shows the total distance that the fibre’s head travels as a function of time. Up until approximately $t=40$ , the distance travelled is monotonically decreasing in Weissenberg number. Deviations in the fibre paths and shape deformations again destroy this monotonicity at later times.

Figure 9. (a) Moment when the fibre arrives at successive stagnation points for different $Wi$ . (b) Distance that the fibre’s head travels up through time $t=80$ for different $Wi$ .

Because evolving fibre geometries and paths can vary greatly depending upon their initial placement and orientation (Quennouz et al. Reference Quennouz, Shelley, du Roure and Lindner2015), we restrict the comparison of the bending behaviours of fibres at different Weissenberg numbers to their first trip, as each is initialized identically. Figure 10(a) shows the bending energy of the fibres up to time $t=15$ . We see that the fibre in the Newtonian fluid bends earlier than the fibres at $Wi\neq 0$ , and also meets its second stagnation point during this time interval at approximately $t=11$ . The Newtonian fibre has a larger bending energy than the fibres at $Wi\neq 0$ as it meets the stagnation point. While the time of buckling is monotonic in Weissenberg number, the bending energy is not. One feature of note, however, is that the secondary peak of bending energy noted in the $Wi=1$ case above is much more pronounced in the $Wi\neq 0$ cases, with only a slight blip appearing in the Newtonian case (figure 10 a). This overshooting of curvature as the fibre straightens out is enhanced by the memory of polymer stresses distributed near the fibre.

Figure 10. (a) Bending energy for the fibre at different $Wi$ during the first trip. Here $Re=2.5$ and $\unicode[STIX]{x1D702}=370$ .

5 Effect of elasto-viscous number $\unicode[STIX]{x1D702}$

Here we examine the effects of flexibility on fibre transport. We vary the non-dimensional bending coefficient $\unicode[STIX]{x1D70E}_{b}$ of the fibre while keeping all the other parameters fixed for $Wi=0$ , 1 and 5. For bending coefficients $\unicode[STIX]{x1D70E}_{b}=0.1$ , 0.05, 0.025 and 0.0125, the elasto-viscous numbers are, respectively, $\unicode[STIX]{x1D702}=370$ , 740, 1480 and 2960. Early experiments by Forgacs & Mason (Reference Forgacs and Mason1959) classified the orbits of fibres in shear flow by tracking pulp, Dacron and elastomer fibres. Rigid rotations, springy rotations, snake turns and coil formations were described. We observe all these types of buckling behaviours when fibres with different elasto-viscous numbers are tracked at each $Wi=0$ , 1 and 5. Figure 11 shows snapshots of fibres of four stiffnesses as they traverse the cellular domain for $Wi=0$  (a) and $Wi=1$  (b). Each fibre was initially released at the same position. In each panel we superimpose successive positions of fibres of varying flexibility at equally spaced time intervals up through $t=17$ . The velocity field shown is the initial velocity field. As in the previous section, we are not able to discern a clear pattern of differences in fibre shape deformations with the same elasto-viscous numbers across $Wi$ by comparing bending energy or scaled end-to-end lengths. However, fibre trajectories at all stiffnesses show significant qualitative differences with respect to the viscoelastic properties of the fluid.

Figure 11. Fibre configurations in (a) Newtonian fluid and (b) $Wi=1$ at different times. For bending coefficients $\unicode[STIX]{x1D70E}_{b}=0.1$ , 0.05, 0.025 and 0.0125, the elasto-viscous numbers are, respectively, $\unicode[STIX]{x1D702}=370$ , 740, 1480 and 2960. Note that each panel represents the results of four simulations, with each of the fibre positions at the times $t=1$ , 3, 5, 7, 9, 11, 13, 15 and 17.

In a four-roll-mill Newtonian flow with no immersed fibres, streamlines are closed in each cellular region. A point particle will simply follow a vortical streamline and never escape the cell in which it was launched. However, a finite-length fibre generates forces to alter the surrounding velocity field and can cross streamlines, allowing for the possibility of escape. Young & Shelley (Reference Young and Shelley2007) studied the transport of fibres in Newtonian cellular flows and identified a region within the cell such that, when the fibre is launched in that region, it remains trapped. The experiments of Wandersman et al. (Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010) varied the distance between the fibre placement and the vertical centreline and measured the number of trajectories required for the fibre to escape. They found that flexible fibres escaped the original cell more frequently than rigid fibres. Quennouz et al. (Reference Quennouz, Shelley, du Roure and Lindner2015) expanded these experiments and concluded that escape is determined not only by the elasto-viscous number, but also by the fibre’s entry position and orientation.

Figure 12. Trajectories of centre of mass of fibres for different $Wi$ through time $t=80$ . (a) Stiffest fibre $\unicode[STIX]{x1D70E}_{b}=0.1$ ( $\unicode[STIX]{x1D702}=370$ ). (b) Most flexible fibre $\unicode[STIX]{x1D70E}_{b}=0.0125$ ( $\unicode[STIX]{x1D702}=2960$ ). In all cases the trajectories begin at the spatial position ( $3\times 10^{-3}$ , $3\times 10^{-6}$ ) (open circle).

For fibres launched close to the incoming streamline, as in figure 1(a), figure 12(a) shows the trajectories of the most rigid fibre ( $\unicode[STIX]{x1D702}=370$ ) in simulations of different $Wi$ and figure 12(b) shows the corresponding trajectories of the most flexible fibre ( $\unicode[STIX]{x1D702}=2960$ ). There are some important features to note. The most rigid and most flexible fibres at $Wi=4$ and $5$ remain trapped in a cellular region. All the flexible fibres at each $Wi$ are able to escape the original cell in which they were launched. This is consistent with the observations in Wandersman et al. (Reference Wandersman, Quennouz, Fermigier, Lindner and du Roure2010). However, for $Wi\neq 0$ , all of these flexible fibres remain trapped in the second cell they encounter. The most rigid fibres at $Wi=0,1$ do not remain trapped in any cell up through the simulation time $t=80$ . In contrast, the most rigid and most flexible fibres at $Wi=4$ and $5$ remain trapped in a cellular region as they leave the first stagnation point.

Figure 13 tracks the moment when the most rigid and most flexible fibres arrive at successive stagnation points (which could be repeated) in Newtonian, $Wi=1$ and $Wi=5$ flows. The rigid fibre always lags behind the flexible fibre in arrival times, and this is more pronounced at $Wi\neq 0$ . Note that, for the case $Wi=1$ , the rigid fibre is traversing multiple cells, while the flexible fibre remains trapped.

Figure 13. Times at which the most rigid fibre and the most flexible fibre arrive at successive stagnation points at different Weissenberg number: (a)  $Wi=0$ (Newtonian), (b)  $Wi=1$ and (c)  $Wi=5$ .

Figure 14. Contours of polymer force densities when the most rigid fibre (a,c,e, $\unicode[STIX]{x1D702}=370$ ) and most flexible fibre (b,d,f, $\unicode[STIX]{x1D702}=2960$ ) move in the cellular flow at $Re=2.5$ and $Wi=3$ : (a)  $y$ -component,  $\unicode[STIX]{x1D702}=370$ ,  $t=4$ ; (b)  $y$ -component,  $\unicode[STIX]{x1D702}=2960$ , $t=2.5$ ; (c)  $x$ -component,  $\unicode[STIX]{x1D702}=370$ , $t=15$ ; (d)  $x$ -component,  $\unicode[STIX]{x1D702}=2960$ , $t=14$ ; (e)  $y$ -component, $\unicode[STIX]{x1D702}=370$ , $t=25$ and (f)  $y$ -component,  $\unicode[STIX]{x1D702}=2960$ , $t=25$ .

We performed systematic studies that varied the initial placement of both the rigid fibre $\unicode[STIX]{x1D70E}_{b}=0.1$ and the flexible fibre $\unicode[STIX]{x1D70E}_{b}=0.0125$ within the vortical cell, at Weissenberg numbers $Wi=0$ , 1, 2, 3, 4 and 5. The fibres were always launched vertically in their straight equilibrium shape, along $x=b$ , where $b$ varied from $0.01$ to $0.09$ , and the lower point of the fibre was at $y=0.5$ . For all initial placements of the most rigid fibres, the fibres in the Newtonian fluid eventually escaped the original cell. When placed well within the vortical cell, even the fibre in the $Wi=1$ case remains trapped when $b\geqslant 0.07$ . For all values of $b$ , the rigid fibres in the $Wi=4$ and $5$ cases never escape. For the most flexible fibres, when initialized well within the cell, no fibre escaped – even in the Newtonian case. However, when placed further at the cell periphery, the flexible fibre does escape in the Newtonian case, but those in the non-zero $Wi$ cases do not.

To understand why the presence of the polymers in viscoelastic fluids hinders the fibre’s ability to escape, we examine the force densities $\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}$ that are generated by the polymers. Figure 14 shows the contour plot for these polymer force densities in the $x$ - and $y$ -directions for the most rigid fibre (left column) and the most flexible fibre (right column) at instances where the fibre is near the edge of the cellular region $(Wi=3)$ . In figure 14(a,b), when the fibres are near the bottom edge, large positive $y$ -components of polymer force density (outward $y$ -direction) are generated near the fibre’s head and push the fibre back into the cell. When the fibres are near the right edge (figure 14 c,d), large negative $x$ -components of polymer force density (inward $x$ -direction) are generated, again impeding escape. At time $t=25$ , figure 14(e) shows that large force densities in the inward $y$ -direction push back the most rigid fibre, which does get trapped in the cell as time progresses. Figure 14(f) shows that at the same time $t=25$ the most flexible fibre is already deep into the cellular region, and large force densities are not required to keep it confined to the cell.