Orientation Dynamics of Rigid Fibers in a Microfluidic Burgers-like Vortex

Introduction

Flows laden with suspended particles are ubiquitous in nature and industrial processes, particularly when the particles are non-spherical or elongated [33]. Such anisotropic particle suspensions arise in applications ranging from paper manufacturing and composite material processing to drag-reduction strategies, and have gained renewed attention in the context of environmental microplastic pollution, where a large fraction of debris consists of fiber-like particles ([25, 9, 6, 36, 40]).

The dynamics of elongated particles are enriched by an additional orientational degree of freedom. In canonical flows, such as simple shear, extensional or steady flows, rigid fibers are advected along streamlines while continuously rotating and reorienting in response to the local velocity gradients ([18]). In more complex configurations, including mixed, unsteady and turbulent flows, the competition between strain and rotation gives rise to richer orientational dynamics ([43, 27]). However, the inherently unsteady nature of turbulence continually modifies the local velocity-gradient tensor, making the underlying alignment mechanisms difficult to isolate from instantaneous observations. Statistical studies have shown that elongated particles preferentially align with the direction of fluid stretching and, consequently, with the local vorticity ([30, 29, 31, 32, 34]). For sufficiently small neutrally buoyant particles in viscous flows, these dynamics are accurately described by Jeffery equations, which relates particle rotation to the local strain-rate and rotation-rate tensors ([18]). Jeffery’s framework has since been generalized to account for particle activity, flexibility, asymmetry and inertia ([20, 24, 12, 8, 10]).

Despite extensive investigations of fiber dynamics in canonical flows such as simple shear and elongation, as well as in turbulence, comparatively little attention has been devoted to stationary flow configurations combining strong rotation and strain, such as isolated vortices. Previous studies have examined the interaction of microswimmers with Burgers vortices, which serve as simplified models of turbulent vortex tubes ([19, 46, 39]). Vortical flows have also been exploited in biological applications, for example to enhance the production of extracellular vesicles from cell spheroids ([41]). However, the transport and orientation dynamics of passive rigid fibers in stationary stretched vortices remain largely unexplored.

In this paper, we investigate experimentally and numerically the orientation dynamics of rigid neutrally buoyant fibers in a stationary three-dimensional vortex at intermediate Reynolds number. This flow provides a controlled analogue of a Burgers-like vortex, relevant to vortical structures encountered in turbulence, enabling a precise characterization of the fiber orientation dynamics. Initial observations show that fibers tend to rotate around the vortex while aligning with a preferential direction (see experimental stacks in Figure 1), reminiscent of observations in turbulent flows. After introducing the experimental set-up and numerical methods in section 2, we adopt a textbook approach and derive the fiber orientation dynamics in a Burgers vortex from Jeffery equations and then compare this theoretical prediction with a complete experimental and numerical characterization in section 3.2 and section 3.3. Finally, we discuss deviations from Jeffery’s framework arising from finite particle size and fluid inertia.

Experimental and numerical methods

Experimental set-up and base flow

In the present study, a steady three-dimensional vortex is generated using a microfluidic cross-slot geometry (Figure 1 (a)). Two opposing inlet streams directed along the $y$-axis meet at the center of the junction, producing a stagnation point. Above a critical Reynolds number, $Re_c$, the base flow becomes unstable and develops into a Burgers-like vortex extending downstream into the outlet channels ([14, 4]).

The flow Reynolds number is defined as $Re = \frac{\rho_f Uw}{\eta}$, where $U = Q/(wd)$ is the mean flow velocity in each inlet and outlet channel, $Q$ is the imposed flow rate, and $\rho_f$ and $\eta$ are the fluid density and dynamic viscosity, respectively. Experiments are conducted in two microchannels of square cross-section ($w = d = 1.1$ mm, aspect ratio $AR = 1$): a glass device used for flow visualizations in the plane $x = 0$, referred as the Core Vortex Plane (CVP) [28] [3], and a PDMS device used for visualizations in the plane $z = 0$, referred as the Stretched Vortex Plane (SVP). Experimental chronophotographs are shown in Figure 1 (b): fibers rotate around the vortex core in the CVP and exhibit ”ribbon-like” dynamics in the orthogonal SVP. For this geometry, the instability onset occurs at $Re_c \simeq42$, and all experiments are performed in the range $40 < Re < 80$, for which the vortex remains stationary [5]. Further details regarding the experimental set-ups are provided in [2].

Quantitative measurements of the base flow are carried out in the CVP using microparticle image velocimetry ($\mu$-PIV) in a 25 wt% glycerol-DI water solution, and are complemented by direct numerical simulations of the single-phase flow [2]. The vorticity field displayed in Figure 1 (c) reveals the presence of separatrix streamlines in the CVP (highlighted in red), which delimit the regions that interact strongly with both the vortex core and the stagnation point as opposed to those that remain essentially unaffected by the vortical structure. Within this confined central region, vorticity is highly concentrated around the stagnation point, forming a narrow column aligned with the rotation axis $x$. In particular, the transverse vorticity profile in the vortex core closely follows a Gaussian distribution. These features are consistent with the analytical Burgers vortex solution, for which vortex stretching balances viscous diffusion to maintain a localized vorticity distribution. This agreement is further illustrated in Figure 1 (d).

Using the numerical simulations, the transverse vorticity profiles at different downstream positions $x/w$ are fitted with the Burgers vorticity profile $\omega_x = \omega_0 \exp\left(-z^2/r_\gamma^2\right)$ with $r_\gamma= \sqrt{2\nu/\gamma}$ the characteristic size of the Burgers vortex and $\nu= \eta/\rho_f$ the fluid kinematic viscosity, in order to estimate the local strain rate $\gamma= \partial u_x/\partial x$. The present flow exhibits a moderate downstream decrease of $\gamma$, from approximately 150 to 115 s$^{-1}$ over a distance of $0.5w$. Note also a small discrepancy between simulations and PIV estimation of $\gamma$, most likely due to slight errors on the position $x$ where PIV is performed due the finite depth of field. Although the strain rate is spatially uniform in an ideal Burgers vortex, this analytical solution remains the most relevant canonical framework for describing the structure and dynamics of the present flow.

Rigid fibers fabrication and characterization

Rigid elongated fibers are fabricated in bulk from an emulsion of SU8 droplets dispersed in a glycerol-ethanol mixture [1, 23]. During emulsification, shear elongates the SU8 droplets, which are subsequently photo-crosslinked under UV illumination to produce stable colloidal fibers. The particles are predominantly straight, with residual curvature below 5%, and possess aspect ratios ranging from 10 to 100. Their dimensions are controlled through the solvent composition and mixing rate. In the present study, a 70/30 wt% glycerol-ethanol mixture stirred at 300 rpm yields fibers of mean diameter $a \simeq2-4\ \mu\mathrm{m}$ and length $L \in[40,500]\ \mu\mathrm{m}$. Owing to the large Young’s modulus of crosslinked SU8 (0.9-7.4 GPa), the fibers can be considered rigid under the present flow conditions. The fibers are suspended in a 25 wt% glycerol-DI water solution ($\rho_f = 1059\ \mathrm{kg}\ \mathrm{m}^{-3}$, $\eta= 1.79\ \mathrm{mPa}\ \mathrm{s}$), resulting in nearly neutrally buoyant particles. The suspension is sufficiently dilute that fiber-fiber interactions are negligible, allowing the fibers to be treated as isolated particles advected by the flow.

It is relevant to this study to define the particle Reynolds number, $Re_p = (L/w)^2 Re$, using the largest particle dimension $L$ and the single phase flow Reynolds number $Re$. In our experiments, the particle Reynolds number is in the range $0.05 \leq Re_p \leq12$.

Numerical method

Simulations of isolated fibers of length $L$ and diameter $a$ are performed in the simulated cross-slot flow used in [2]. To explore similar initial conditions than the experiments, the fibers center of mass is initialized at $(y_0,z_0,x_0)/w = (-3/8,0,\epsilon)$, with $\epsilon\simeq10^{-4}$, $\epsilon\ne0$ so that particles are not trapped in the CVP for an extended period of time, and their initial orientation is varied in the ranges $0 < \beta_0 < \pi/2$ and $\pi/8 < \phi_0 < \pi/2$. In most experimental cases, $\beta_0 \simeq\pi/2$ (fibers enter parallel to the CVP) and $\phi\ne0$ due to their alignment with the slightly inclined flow streamlines caused by elongation along the separatrix. The fiber diameter is set to $a = 4$ microns and the fiber length spans from 40 to 500 microns to match the geometrical properties of the fibers used in the experiments. In the simulations, the particle Reynolds number is assumed to be sufficiently small for the flow around the fiber to remain in the viscous regime ($Re_p \ll1$). Although this assumption may appear restrictive, particularly for the longest fibers considered, the results presented in the following section demonstrate that viscous effects dominate the orientation dynamics and reproduce the experimental observations accurately. The effect of fluid inertia is discussed in greater detail in the Discussion section.

The elastohydrodynamic interactions between the fiber and flow are modeled with the bead-spring model, which has been widely used, and validated, in the literature to study the motion of rigid and flexible fibers in viscous flows [7, 26, 37, 38, 23, 22]. Full details of the method can be found in [26, 23], a brief overview is provided below.

The fiber is modeled as a chain of $n$ rigid spherical beads of radius $a/2$, connected by internal stretching and bending elastic forces $\mathbf{F}^{E}$ that enforce fiber inextensibility and rigidity over time. These forces are derived from an elastic potential $H$, whose expression is given in [26, 23]. This model can be used to simulate both flexible and rigid fibers within the same framework by tuning the Young modulus $E$ according to the situation. In our case, we chose $E = 0.16$ GPa, which is enough to prevent fiber bending in the present flow conditions.

The velocity of the fiber beads is then computed from the linear superposition principle [13]: $\mathbf{U}_i = \mathbf{u}_i^\infty+ \sum_{j=1}^{n} \mathbf{M}_{ij}\mathbf{F}_j^{E}$, where $\mathbf{U}_i$ is the total velocity of bead $i$, $\mathbf{u}_i^\infty$ is the background (cross-slot) flow velocity interpolated at the bead center, and $\sum_{j=1}^{n} \mathbf{M}_{ij}\mathbf{F}_j^{E}$ is the velocity induced by the nonhydrodynamics forces $\mathbf{F}^{E}$ acting on the $n$ fiber beads. $\mathbf{M}$ is the so-called mobility matrix that contains all hydrodynamic interactions between the fiber beads, i.e. their velocity induced by the flow disturbances generated by the nonhydrodynamic forces. In this work, we use the Rotne-Prager-Yamakawa mobility matrix [44], which is commonly used in the bead-model literature. The new position of the fiber beads $\mathbf{x}$ is then computed by integrating $d\mathbf{x}/dt = \mathbf{U}$.

Orientation dynamics of rigid fibers in a vortex

Initial observations

In this study, we investigate the orientation dynamics of rigid elongated fibers in a steady vortex flow. Experimentally, the fiber orientation is reconstructed from complementary measurements performed in the CVP and SVP, providing access to the azimuthal angle $\phi$ and to $\theta$, the projection of the inclination angle $\beta$ onto the SVP, as defined in Figure 1(c). Figure 2(a,c) show experimental chronophotographs in both planes together with reconstructed trajectories colour-coded by $\phi$ and $\theta$, respectively. Owing to the dilute nature of the suspension, only a few fibers are visible in each frame, such that particles may be considered isolated and individually tracked. These measurements are complemented by numerical simulations of confined fibers advected in the vortex. Representative chronophotographs are shown in Figure 2(b,d). Both experiments and simulations reveal that fibers spiral around the vortex axis while being advected downstream, undergoing simultaneous reorientation along their trajectories.

Fibers enter the cross-slot nearly perpendicular to the vortex axis and are initially confined to the CVP ($\theta= \beta\simeq\pi/2$). As they rotate around the vortex core and migrate towards it, the angle $\phi$ oscillates between $-\pi/2$ and $\pi/2$ at each revolution with a period $\tau_{\phi}$ (see Figure 2 (e,f)). Simultaneously, the fibers progressively align with the vortex axis, explaining the ”ribbon-like” structure observed in panel (c): $\beta$ and $|\theta|$ relax from $\pm\pi/2$ towards 0, corresponding to alignment, as shown in Figure 2 (g,h), where both angles are plotted as a function of time for experiments and simulations, respectively. The two angles are related by $\tan\theta= \tan\beta\cos\phi$ and therefore coincide at the maxima of the oscillation ($|\cos\phi| = 1$), indicated by $\star$ in panel (h). Note that, in the experiments, $\beta$ was reconstructed from the apparent fiber length and is therefore significantly noisier than the direct measurement of $\theta$, especially for small angles. Since $\theta$ follows the same alignment dynamics as $\beta$, differing only through the periodic modulation by $\phi$, we use $\theta$ throughout the experimental analysis.

Theory

To explain those dynamics, we use a simple approach combining Jeffery equations for the orientation of a fiber in a viscous flow and the analytical Burgers vortex velocity field.

In Jeffery’s work, anisotropic particles translate with the fluid, just as tracer particles do ($\dot{\mathbf{x}} = \mathbf{u}$), but their rotation depends now on their orientation with respect to the local fluid velocity gradients. We consider here prolate axisymmetric particles. Let us define $\mathbf{p}$ as the unit vector aligned with the symmetry axis of the elongated fiber (see inset in Figure 1 (a)). The overall rotation dynamics of such a particle can be split into a spinning part (rotation around the symmetry axis) and a tumbling part (rotation of the symmetry axis) $\dot{\mathbf{p}}$. In the absence of inertial and Brownian forces, Jeffery [18] found that the symmetry axis of an axisymmetric rigid particle suspended in a simple shear flow rotates along one of an infinite one-parameter family of closed periodic orbits, known as Jeffery orbits [18]. The tumbling $\dot{\mathbf{p}}$, i.e. the time evolution of the particle orientation, is given by:

with $\mathbf{\Omega}$ the fluid rotation rate tensor, $\mathbf{E}$ the fluid strain rate tensor and $\kappa= \frac{(L/a)^2-1}{(L/a)^2+1}$ a shape parameter depending on the aspect ratio of the length of the fiber along the symmetry axis $L$ and its diameter perpendicular to the axis $a$. Note that the second term between parentheses on the right hand side keeps $|\mathbf{p}| = 1$. In the case of very elongated fibers, $L/a \gg1$, the parameter $\kappa$ is close to 1.

Following this definition of $\mathbf{p}$, the unit vector aligned with the symmetry axis of the fiber, we apply equation (1) to the transport of an elongated fiber in a Burgers vortex velocity field defined in cylindrical coordinates by:

where $r = \sqrt{y^2+z^2}$ is the radial distance to the vortex center, $\gamma$ is the strain rate, $\Gamma$ the circulation, and $\nu$ the kinematic viscosity. The vorticity profile is given by:

The Burgers vortex is characterized by the concentration of vorticity around the vortex axis in an area of characteristic size $r_\gamma= \sqrt{2\nu/\gamma}$ (i.e. the standard deviation of the Gaussian vorticity profile).

Near the vortex core, i.e. for $r \ll r_\gamma$, the azimuthal velocity and the vorticity can be approximated by:

Each component of the particle orientation $\mathbf{p}$ provides the azimuthal and polar angles, $\phi$ and $\beta$ respectively (see Figure 1 (a)), defined as:

Finally, the Jeffery equations yields:

Substituting the angular parametrization (equation (5)) into the Jeffery equations and projecting onto the local polar direction gives:

Note that this equation is independent of the flow vorticity and of the radial position of the fiber. This is integrable and we can deduce the temporal evolution of the angle $\beta$:

In a Burgers vortex, the classical closed Jeffery orbits observed in a linear shear flow, are suppressed by the axisymmetric extensional strain. Instead, the fiber orientation exponentially relaxes towards two stable fixed points $\beta= 0$ and $\beta= \pi$ corresponding to alignment with the vortex axis. Therefore, an elongated fiber will align with the vortex axis with a characteristic timescale $(3\kappa\gamma)^{-1}$ that depends on the fiber aspect ratio and the vortex strain rate.

Substituting the angular parametrization ((5)) into the Jeffery equations and projecting onto one of the azimuthal directions, we show that the azimuthal angle evolves with the flow vorticity according to :

Once the fiber enters the vortex core, the fluid vorticity tends towards $\omega$, so we have $\phi(t) = \phi_0 + \omega t$. The fiber therefore undergoes uniform precession about the vortex axis while simultaneously aligning with it and both processes are decoupled. This is in qualitative agreement with experimental observation described in Figure 2. Note that, without the approximation made in (4), the vorticity increases along the particle trajectory so that it spins faster as it comes close to the core.

Comparing dynamics in the cross-slot vortex to the idealized Burgers vortex

We now compare the dynamics predicted for a Burgers vortex with the experimental and numerical results. Figure 3 summarizes the orientation dynamics for a range of fiber lengths and initial conditions by showing the temporal evolution of $|\tan\beta|$ or $|\tan\theta|$ (panels (a,b,c)) and the fluid and fiber rotation rates (panels (d,e,f)).

Following a short transient associated with the fiber entering the vortex core, the envelopes of $|\tan\beta|$ in theory (panel (a)) and simulations (panel (b)), and of $|\tan\theta|$ in the experiments (panel (c)), decay exponentially, as predicted by Jeffery equations (10). According to this prediction, the strain-induced alignment towards the vortex axis occurs on a characteristic timescale $(3\kappa\gamma)^{-1}$, providing a direct estimate of the local strain rate $\gamma$. For the highly slender fibers considered here, $\kappa\approx1$, such that the relaxation rate reduces to $3\gamma$. As expected, the Burgers vortex exhibits a unique relaxation rate, independent of fiber length (red line in Figure 3 (a)). Fitting the simulation results yields an average $\gamma\simeq120\ \mathrm{s}^{-1}$ (magenta dashed line in Figure 3 (b)), consistent with the range of strain rates obtained from the simulated base-flow analysis (Figure 1 (d)). Experimentally, we obtain an average $\gamma\simeq100\ \mathrm{s}^{-1}$ (magenta dashed line in Figure 3 (c)), again in good agreement with the strain rate inferred from the PIV measurements (Figure 1 (d)). The simulations reveal a slight dependence of the relaxation rate on the fiber length, a trend that is less clearly resolved in the experiments because of measurement noise, and is discussed further in section 4. We also note the presence of weak oscillations in $|\tan\beta|$ obtained from the simulations. These arise from slight deviations of the simulated base flow from an ideal Burgers vortex ([35]): the vortex core is not perfectly axisymmetric (Figure 1 (b)), causing the strain rate to vary azimuthally and producing weak modulations along the fiber trajectories. These deviations are discussed further in the Supplementary Material.

Finally, we compare the local fluid rotation rate, $\omega_x/2$, obtained from the Burgers solution ((2)) and from the measured and simulated velocity fields, with the fiber rotation rate (Figure 3 (d,e,f)). For the shortest fibers, the azimuthal rotation closely follows the local fluid rotation, accelerating until the fibers reach the vortex core, in excellent agreement with (11). As the fiber length increases, however, systematic departures from the local fluid rotation emerge, revealing finite-size effects that are discussed in section 4.

Discussion and concluding remarks

We have examined the orientation dynamics of elongated fibers in a three-dimensional stationary vortex through experiments and simulations. Overall, the experimental and numerical results demonstrate that orientation dynamics of elongated fibers in the present vortex are governed by the combined action of fluid rotation and axial strain. While the azimuthal motion is controlled by the local vorticity, the strain field progressively aligns fibers with the vortex axis. Despite the finite Reynolds number, the non-ideal nature of the flow and possible finite size effects, the observed dynamics remain remarkably well captured by the canonical Burgers-vortex framework and by Jeffery’s description.

We now discuss possible deviations from this simple picture. The Burgers-vortex description is a priori valid for tracer particles transported along streamlines. For such particles, the trajectories satisfy the power law $r = x^{-1/2}$, indicating that particles are progressively drawn towards the vortex axis as they are transported. Although the simulated fiber trajectories broadly follow this trend, larger fibers exhibit a weaker radial decay (Figure 4 (a)), revealing finite-size effects even in the absence of particle inertia. In the experiments, finite-size effects may be further enhanced by inertial lift forces, which are known to promote cross-streamline migration of finite-size particles ([2]; [21]; [42]; [15]). However, the particle Reynolds number remains small in the present experiments, suggesting that these inertial effects are weak. Together with the smaller range of fiber lengths investigated experimentally and the measurement noise, this prevents a clear identification of these deviations.

Furthermore, Jeffery equations relies on local velocity gradients estimated at the particle center of mass and assumes the flow is locally uniform over the particle length. Beyond a critical length, elongated objects may no longer satisfy this local description, as they sample varying velocity gradients along their symmetry axis. Our simulation results, which span a wider range of fiber lengths and provide higher resolution than the experiments, reveal a weak dependence of the azimuthal rotation period and alignment rate on the ratio of fiber length to vortex size $L/r_\gamma$ (Figure 4 (b,c)), with longer fibers rotating slightly slower and aligning slightly faster than predicted by the local Jeffery description.

In simple shear flows, the existence of an infinite number of marginally stable Jeffery orbits makes the orientation dynamics highly sensitive to weak perturbations arising from particle shape, inertia or other effects, which may drive particles across Jeffery orbits ([47]; [11]; [8]). In contrast, the stretched vortex considered here possesses a stable orientational attractor. The excellent agreement between theory, experiments and simulations indicates that the orientation dynamics remain remarkably robust despite finite particle size, finite Reynolds number, and deviations of the flow from an ideal Burgers vortex. Jeffery equations therefore provides an accurate description of the orientation dynamics over the range of particle Reynolds numbers and fiber lengths investigated here, suggesting that vortex stretching and strain stabilize the orientational dynamics. More generally, this robustness suggests that Jeffery’s description may remain applicable in more complex vortical flows such as turbulence ([45]).

Finally, a natural extension of this study would be to investigate the behavior of flexible fibers, particularly how their shape responds to the varying vorticity and strain encountered in stretched vortices ([27]; [17]; [16]). Such work would extend the present framework to natural and synthetic filaments, where flexibility may substantially alter both alignment and migration dynamics.

Acknowledgments

We thank Arash Belizad Banaei and Luca Brandt for their simulations data of the single phase base flow. We thank Noa Burshtein for fruitful discussions, as well as Simon J. Haward and Amy Q. Shen. AL and MA acknowledge support by the Institut Carnot IPPG (Institut Pierre-Gilles de Gennes), through the Carnot program funded by the Agence Nationale de la Recherche (ANR) and the Ministère de l’Enseignement Supérieur et de la Recherche (MESRI) and the Major Research Program of PSL Research University ”IPPG” launched by PSL Research University and implemented by ANR with the references ANR-10-IDEX-0001. This work benefited from the technical contribution of the joint service unit IPGG Technological Platform CNRS UAR 3750. BD acknowledges support from the French National Research Agency (ANR), under awards ANR-20-CE30-0006 and ANR-25-CE30-4710.

Declaration of Interests

The authors report no conflict of interest.