Rodlike fillers in a flow field of a viscous fluid may form complex structures after passing a sudden contraction. The rods start with a dilute distribution with random positions and orientations. Behind the contraction, a large amount of rods tumble in a spatially correlated way, such that orientations perpendicular to the flow field occur at regular distances along the channel. The correlated tumbling results from an interplay of several effects, the tumbling inferred by the space dependent shear flow, the accumulation of rods at a certain distance from the wall, and the rod alignment at the contraction. The system is studied numerically for rodlike fillers in a shearthinning viscous fluid.
I. INTRODUCTION
The dynamics of rodlike particles in the flow field of a viscous solution is of great importance for the production of composite materials. One important application is the creation of short fiber reinforced thermoplastics by injection or compression molding.^{1,2} The use of fillers like nano or microfibers opens up new opportunities for adjusting the matrix properties, such as viscosity and toughness.^{2,3} It is a widespread procedure to use composite approaches to expand the parameter range of materials, for example, to improve the tensile strength of weaker matrix materials like contemporary carbon fiber reinforced polymers.^{4} Rodlike fillers are also used to change electric and thermal conductivity.^{5,6} The physical characteristics of composite materials depend on the spatial distribution and alignment of the fillers. Especially, fillers with a preferred orientation induce anisotropic material properties.
Another application field is the 3Dprinting of hydrogels with rodlike fillers, which can potentially be used in tissue engineering approaches and biofabrication. In the field of biofabrication, the ink made of a hydrogelbased matrix containing short fragments with rod shape morphology will be extruded in a predefined morphology.^{7} The ink usually has viscoelastic properties. Typically, the aim is that fiber fragments in the ink are uniformly distributed and get aligned with the flow during the 3D printing in order to create a unidirectional distribution that modifies and enhances the properties of the printed material. The rodshaped fillers can improve material properties like the stability during the printing process or in the created biomaterial.^{8,9} In the printed structure, rodshaped fillers may serve as a guidance for anisotropic cell growth.^{8,9} Various other biomaterial properties depend on the orientational and spatial order of the rodlike fillers, so that a control of the filler dynamics during the printing process is of high relevance.^{4,10–12}
The flow of a solvent through a cavity has an amplitude that decreases toward the channel wall so that rodlike fillers are exposed to shear forces. This has an influence on the orientation of the rods; in many cases, one finds an average alignment of the rod axes in the direction of the flow field. However, especially for dilute rod systems, rods may still rotate but spend more time in the flow direction. The behavior of rodlike fillers in a confined flow has been investigated in several experimental studies.^{10,13–17}
In many molding or 3Dprinting processes, the melt or hydrogel flows through a contraction as it enters a circular or slotshaped nozzle, a narrow mold, or the needle of a 3Dprinter. In this article, we study this rather common scenario with computer simulations and find that dilute rod suspensions form a remarkable structure in the narrow channel behind the contraction.
In the simulations, we study a system that excludes all effects that are not essential for the discussed structure formation. We investigate the individual dynamics of rodlike fillers in a viscous, nonNewtonian matrix, flowing from a broad to a narrow channel through a sudden contraction. At the beginning, rodlike fillers in the broad channel start with random spatial and orientational distribution. We assume that the rod density is low so that rods hardly interact with each other. The rods stream with the flow field and pass the contraction into the narrow channel, as sketched in Fig. 1. In the narrow channel, they tumble in a way, that the average alignment of rods is strongly correlated with the distance from the contraction. A rod orientation perpendicular to the flow occurs at predictable points with fixed distance to each other.
It is known that individual rods in a shear flow of a fluid tumble with time forming socalled Jeffery orbits.^{18} The remarkable observation in our system is that a large amount of rods tumble with the same periodicity and phase shift, in spite of (a) their random starting conditions and (b) the fact that the shear rate of the flow field varies significantly along the cross section of the narrow channel.
Rods in a flow field have been investigated with various different computational methods. One direct approach is the use of coarsegrained molecular dynamic simulations (CGMD).^{19,20} Other methods are based on Brownian motion,^{21} the reciprocal theorem,^{22} or a grand mobility matrix that considers hydrodynamic interactions of particle surfaces.^{23} Another attempt includes smoothed particle hydrodynamics (SPH),^{24} which is sometimes combined with the discrete element method (DEM).^{25} For dilute systems, in which rods rarely meet each other, the motion of rods can be studied separated from each other. The dynamics of an ellipsoidal rod in a shear flow has been determined by Jeffery.^{18} A comprehensive listing of particlebased methods is given in the fourth section of an article by Kugler et al.,^{26} which otherwise focuses on macroscopic models for the orientation field.
Orientation fields are usually tensor fields that characterize the local orientation distribution as a function of position. In methods based on orientation fields, rods are not considered individually. Frequently, the dynamics of the tensor field follows an extension of Jeffery's equations,^{18} where additional terms consider the average rod–rod interactions and the effect of Brownian motion. Many of these numerical methods are based on the Folgar–Tucker model.^{27}
The motion of rodlike fillers in a viscous matrix depends on the flow field, the interaction with other rods and the impact of walls enclosing the fluid. The flow field depends on the geometry of the system, wall interactions, the viscosity and viscoelasticity of the matrix material, and the interaction with the rods. As discussed by Abtahi and Elfring,^{28} the viscosity of a shearthinning matrix may also be altered by a shear flow, induced locally by the rod. An explicit consideration of the rod impact on the flow field requires a selfconsistent solution of the rod and the fluid dynamics. If the influence of the rods on the fluid is restricted to short distances from the rods, one can calculate the undisturbed flow field and use it to determine the rod dynamics afterwards. If the Péclet number $ P e = \gamma \u0307 / D r$ with shear rate $ \gamma \u0307$ and rotational diffusion constant D_{r} is low, Brownian dynamics and diffusion get relevant and the dynamics of the rods changes qualitatively. Recently, orientation dependent diffusion has been studied for anisotropic particles in a Poiseuille flow, which depends significantly on the Péclet number.^{13} Experiments with rods flowing through a channel in a nonNewtonian fluid have reproduced periodic orbits predicted by Jeffery. With lower Péclet number, rods may show various types of oscillatory motions like tumbling, kayaking, or logrolling.^{29,30} The tumbling gets aperiodic and rods can rotate out of the shear plane.^{29,30}
In this article, the dynamics of rods is studied explicitly, using the following assumptions:

The mass of the fillers is small enough so that inertia terms in the equations of motion are negligible.

The Péclet number is large so that Brownian motion can be ignored.

The influence of a rod on the flow field is restricted to small distances from the rod.

The system of rods is dilute, and rod interactions play no role.

For the sake of clarity, we investigate the contraction of a planar channel, though we expect similar effects in cylindrical geometries.

We consider a steadystate flow field obtained for a shearthinning fluid. Rod induced shearthinning is neglected as a higher order effect.
We start our discussion with the simple example of a rod in the shear plane of a flow field $ u x ( y )$ with constant shear rate $ \gamma \u0307 = \u2202 u x \u2202 y$. In this special case, the center of the rod drifts with a constant velocity v_{0} in the x direction, while the rod axis tumbles periodically with an angle $\varphi $. For $ \varphi = 0$, the rod axis is parallel to the x axis. The same orientation occurs at $ \varphi = \pi $ and then the tumbling repeats. If the rotation from $ \varphi = 0$ to $ \varphi = \pi $ takes a time T, then the rod moves a length $ \lambda = T v 0$ until the tumbling repeats.
It is tempting to utilize the spatial periodicity of the rod dynamics to create a composite material with embedded rodlike fillers forming a regular pattern. In practice, there are two problems. (a) The tumbling rods need a synchronized phase in space. If rods tumble with the same frequency but start at a point x_{0} with random orientations, the resulting average orientation distribution will be spatially homogenous. (b) In nearly all molding or 3D printing processes, the shear rate is not constant in space. All viscous fluids that stream through a pipe or a planar channel have a spatially dependent shear flow like the Poiseuille flow field in the case of a Newtonian fluid. In this article, we show with numerical methods that both problems, (a) and (b), are strongly suppressed in the case of dilute rod distributions flowing through a sudden contraction. At the point where the fluid flows from a broad channel into a narrow one, the rodlike fillers tend to align with the direction of the narrow channel. One of the reasons is the acceleration, leading to an extension of the flow field. The rod alignment in extensional flows has been studied experimentally and numerically for rod distributions in confined flow fields.^{31,32} Of special interest for this article is the fact that this rod alignment is found in the accelerating flow field at channel contractions and pore entrances.^{33–35} At inverse contractions, rods orient away from the channel direction.^{34}
In our system, the extensional flow at the contraction leads to similar orientation of fillers behind the contraction, which leads to a reasonable synchronization of the rods, tumbling in the narrow channel. The distance λ that a rod travels during a half rotation depends on the distance y_{c} to the channel center. Therefore, rods with different y_{c} will desynchronize soon. This problem is strongly reduced by a second effect of the contraction: In the broad channel, a large number of rods are far away from the center of the channel. Passing the contraction, they end up in the outer region of the narrow channel. Here, they tumble and hit the wall so that further rotation shifts their center of mass to a distance half the rod length away from the wall, where the rod can rotate freely. This means a large fraction of the rods ends up at nearly the same distance from the wall, where they are exposed to the same shear rate.
II. THEORETICAL BACKGROUND
A. Flow dynamics
B. Rod motion
Inertial effects shall be negligible, so that the total force $ F = m v \u0307 c$ and the total torque $ M 0 = \Theta \Omega \u0307 0$ must vanish. Dynamic equations for $ v c$ and $ \Omega 0$ result from $ F = 0$ and $ M 0 = 0$.
We assume that the fluid sticks to the rod surface. Thus, at each surface point r, the fluid velocity is equal to $ v ( r )$, which typically differs from the flow field $ u ( r )$ in the absence of the rod.
The quantity $ \kappa L / D$ depends sensitively on the shape of the rod. Jeffery found $ \kappa L / D = 1$ for ellipsoids of revolution.^{18} Using $ c \u2225 / c \u22a5 = 1 / 2$ (Ref. 41), Eq. (21) leads to $ \kappa L / D = 3 / 2$, which for $ L / D = 8 \xb1 2$ is in the same range as the values found by Cox for cylindrical rods.^{42} The quantity c_{r} matches with the rotational friction coefficient.^{37}
1. Approximations for the integrals
C. Rod flow between planeparallel walls
It is useful to characterize the rod geometry by κ which includes the axis ratio and the absolute cylinder length L. In the following, we always use $ \kappa 2 = 1 / 108$, even if we vary the rod length L. First, we study the dynamics of rods, tumbling in the x, y plane in a planar channel.
The rod length appears only in $ q ( \varphi )$, a positive term that oscillates with $ sin \u2009 ( \varphi )$ and is proportional to L^{2}. For rods in the Newtonian fluid, the velocity is reduced by a term proportional to $ q ( \varphi )$, while, interestingly, $ \Omega ( \varphi )$ does not depend on L at all, as long as the rod does not collide with the channel wall.
Also, for rods in the shearthinning fluid, the velocity is decreased by a term proportional to $ q ( \varphi )$. The angular velocity Ω is increased by a term proportional to $ q ( \varphi )$. The lowest velocities are always found for $ \varphi ( x ) = \pi / 2 , \u2009 3 \pi / 2 , \u2026$ Examples of $ v x ( \varphi )$ and $ \Omega ( \varphi )$ of a rod in a planar channel are shown in Fig. 4.
For a shearthinning matrix, the integral is more complicated. In the following, a case of special interest is a rod in a shearthinning matrix with n = 1/2 at height $ y c = w / 2$ with rod length $ L \u2243 w$.
In general, λ decreases strongly with y_{c}. The rotation dynamics in a channel is visualized in Fig. 5, which shows $ \u2009 cos 2 ( \varphi ) ( x c , y c )$ for rods with length $ L = 0.9 w$ and $ \kappa 2 = 1 / 108$ that start with $ \varphi = 0$ at x_{c} = 0.
D. Rod flow through a sudden contraction
Now we come to the main subject of this article. We study rods in a fluid matrix that flows through a planar channel with a sudden contraction at x_{co} where the width reduces spontaneously from $ 2 w 0$ to 2w. A steadystate flow field was created with the help of the program OpenFoam for a system with $ 2 w 0 = 2.4 \u2009 mm , \u2009 2 w = 0.8 \u2009 mm , \u2009 x c o = 4 \u2009 mm = 10 w$, and a kinematic pressure decreasing with $ \Delta P / \Delta x = \u2212 1125 \u2009 ms \u2212 2$. The fluid velocity at the walls is $ u = 0$. The kinematic viscosity of the shearthinning matrix obeys the powerlaw equation (1) with flow behavior index n = 1/2 and flow consistency index $ K = 5.67 \xd7 10 \u2212 3 \u2009 m 2 s n \u2212 1$. The x and y components of the stationary flow field are shown in Fig. 6.
For our system, we can estimate the Péclet number by using the approach used by Kobayashi and Yamamoto.^{43} For the estimation, we use conditions in our model system, a temperature $ T = 300 \u2009 K$, and a mass density $ \rho = 1 \u2009 g / cm 3$, typical for hydrogels. Considering a particle in the thin cylinder that is one rod diameter away from the cylinder axis, we find a very large Péclet number of $ P e \u2265 10 9$. This is caused by the large size of the fillers in our system. For fillers in the micrometer range under the same conditions, one has $ P e \u2243 400$.
With the given flow field, we calculate the dynamics of rods, represented by spherocylinders. The diameter is $ d = 0.1 w$, and the cylindrical part of the rod has a length $ L = 0.9 w$ so that the total rod length is $ L tot = w$. Furthermore, we use $ \kappa 2 = 1 / 108$. At the walls, we assume a perfect slip motion of the rods. This means, at the collision point, the velocity normal to the collision axis does not change.
We investigate the formation of ordered structure by rods starting with a random distribution. The proper choice of accurate random starting conditions is not selfevident. We consider that in a region x < 0, rods have a random, homogenous orientational and spatial distribution. Thus, at each height y_{c}, rod centers have the same average distance $ \Delta x$ in the x direction. In a time interval $ \Delta t$, rods appear at x_{c} = 0 with an average frequency $ f ( y c ) = u x ( y c ) / \Delta x$. The amount of rods $ n 0 ( y c ) \Delta t$ per time $ \Delta t$ is proportional to $ u x ( 0 , y c ) \u221d 1 \u2212 ( y c / w 0 ) 3$. Therefore, in our simulation, the number of rods starting at $ ( 0 , y c , 0 , 0 )$ is weighted with the factor $ 1 \u2212 ( y c / w 0 ) 3$.
Approximately 32% of the rods end up in the upper region $ y c > y \u0303 c$. This percentage can be increased by choosing a larger $ w 0 / w$ ratio. As the rods continue to drift in the narrow channel, they keep rotating in the x, y plane, until they hit the channel wall, see Fig. 8. Wall repulsion presses the rod center away from the wall as the rod keeps rotating. As the rod gets perpendicular to the wall plane, the rod center reaches $ y c = y \u0303 c$, where it remains during the subsequent rotations. In Fig. 7, kinks in the paths indicate points at which rods hit the wall and the rod center reaches $ y c = y \u0303 c$. In the end, contraction, rotation, and wall repulsion let 32% of rods end up a distance, $ L tot / 2$ from the walls. Here, all these rods are exposed to the same u_{x} and $ \gamma \u0307$.
Now we study rods with 3D orientation, which start with random directors n, distributed homogeneously on a unit sphere. Figure 9 shows a logarithmic contour plot of the spatial distribution $ n c ( x , y )$ of rod centers captured from the whole paths of all rods. One can see a pronounced maximum at $ y \u2243 w \u2212 L tot / 2$. Analyzing $ n c ( x , y )$ shows that for $ y > 50 w$ about 30% of the rod centers are in the range of $ 0.5 w \u2264  y c  \u2264 0.55 w$. The orientation of the rods is reflected in the spatial distribution $ n tip ( x , y )$ of cylinder ends $ r c \xb1 L 2 n$, shown in Fig. 10.
At constant distances λ, there are spikes of $ n tip ( x , y )$ indicating that rod ends get close to the wall. At these positions, rods are perpendicular to the flow direction. In between, rods are nearly parallel to the x axis, most of the time. If this is the case, the cylinder ends are roughly at the same y position as the rod center, $ y c \u2243 \xb1 y \u0303$, which leads to the maxima of n_{tip} at $ \xb1 y \u0303$.
At the starting point, $ x c \u2243 0$, the quantities A_{xx}, A_{yy}, and $ A z z = 1 \u2212 A x x \u2212 A y y$ are 1/3, which correspond to an isotropic distribution. At $ x c \u2243 x c o$, one has $ A x x \u2243 0.78$. This indicates a distinct alignment in the flow direction. At fixed distances λ, A_{yy} shows small peaks, which are much more pronounced for $ A y y m$. The peaks coincide with sharp minima in $ A y y m$ and indicate that a larger amount of rods with $ y c \u2243 \xb1 y \u0303 c$ point in the y direction at these x positions.
At the contraction, the flow field accelerates, and the foremost part of a rod is torn in the flow direction, resulting in an alignment of the rods almost parallel to the x axis. The effect holds for all rods, independent of their starting point or initial orientation. When they point in the x direction, rods are exposed to the shear rate $ \gamma \u0307 ( y ) = \u2202 u x \u2202 y ( y )$ and tumble with an angular velocity perpendicular to the x, y plane. If their center of mass is above $ y \u0303 c$, they hit the wall and end up at $ y c \u2243 y \u0303 c$. The large amount of rods with $ y c = y \u0303 c$ is exposed to the same shear flow $ \gamma \u0307 ( y \u0303 c )$ and so they all rotate with the same angular velocity sequence. The rod alignment at the contraction leads to $ \varphi \u2243 0$ at this point so that many rods follow the same path with the same orientation sequence $ \varphi ( x c )$. The rods at $ y c = y \u0303 c$ contact the channel walls at regular distances, when the rods are perpendicular to the flow direction $ \varphi ( x ) = \pi / 2$.
We have investigated the behavior of rods entering a contraction for the case that the fluid is shearthinning. This was chosen because in many applications shearthinning fluids are used. At this point, it must be emphasized that shearthinning is not required for the observed structure formation. As an example, the spatial distribution $ n tip ( x , y )$ of cylinder ends is presented in Fig. 12 for rods passing a contraction in a Newtonian fluid. All other parameters are used as in the nonNewtonian case. One can see that the synchronized rotation of the rods is visible just like in the nonNewtonian case.
III. SUMMARY
Tumbling of rods in a viscous fluid that flows through a channel depends on the distance from the channel center y_{c} and the orientation at an initial point x_{co}. The distribution of these two parameters can be strongly narrowed by a regular sudden contraction. Passing the contraction, rods are strongly aligned with the flow direction, which serves as a reference point for the tumbling phase. Remarkably, even rod axes that deviate distinctly from the x, y plane are dragged at the contraction into the x direction from where they start a tumbling motion in the x, y plane. Shortly after the contraction, a large amount of rods is localized at a distance $ L tot / 2$ from a channel wall. Here, the rods are exposed to the same shear rate and tumble with the same frequency. Our simulations show a spatially synchronized tumbling of rods that have started in the wider channel with a random, uniform distribution of starting position and orientation. Rod orientations perpendicular to the flow direction occur at equally spaced positions $ x n = x 0 + n \lambda $ with $ n = 0 , 1 , 2 , \u2026$.
Those who aim for a spatially homogenous distribution of rods pointing in the flow direction must be aware of two aspects, if the filler density is low. (a) Fillers are far from being distributed homogeneously, instead the density has distinct local maxima at $  y  = w \u2212 L / 2$. Here, more than 30% of rods may be localized, depending on the width ratio of broad and narrow channel and the ratio of rod length and narrow channel width. (b) The fraction of rods that end up in the high density layers is not permanently aligned but tumble.
The tumbling of rods forms a regular pattern in space. This way the contraction provides a rather complex structure of fillers without any extra effort. If the matrix solidifies in the narrow channel, the resulting material may have spatially structured material properties with a periodicity λ that may be several hundred times larger than the rod length. Mechanical strength may vary in space periodically. For systems with electrically conducting fillers, it might be especially interesting that rods touch the surface of the matrix at regular distances. Thus, one ends up with a compound layer with spatially controlled contact points, see the conceptual visualization in Fig. 1.
In this article, we studied a somewhat idealized system, in which the reasons of the observed effects are clarified by leaving out various aspects. However, extended simulations indicate that the effects still occur in altered systems. They are found for different κ and other ratios of rod length L_{tot} to narrow channel width w. The effect does not vanish for rod densities with occasional rod interactions, and a synchronized tumbling is also found in contractions of cylindrical pipes. All these observations are of great interest, but go beyond the scope of this article and must be analyzed in independent, extensive studies.
Many other aspects are worth studying in the future. The rodwall interactions may be varied as well as the local geometry of the contraction. One should explore a critical density at which the system switches from a tumbling to a constantly aligned state. Length distributions of rods are another important point.
The requirements for performing the respective simulations are expounded here, in this article. Especially, it is shown how to treat the dynamics of rods in a flow field that may not be assumed to be linear on the length scale of the rods.
ACKNOWLEDGMENTS
This project is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project number 326998133 – SFB/TRR225 (subproject B03PI Sahar Salehi).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Thomas Gruhn: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Project administration (equal); Software (lead); Visualization (lead); Writing – original draft (lead). Camilo Ortiz Monsalve: Conceptualization (equal); Investigation (supporting); Methodology (supporting); Software (supporting); Writing – review & editing (equal). Sahar Salehi: Conceptualization (lead); Data curation (lead); Funding acquisition (lead); Project administration (equal); Validation (supporting); Writing – original draft (supporting); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.