Rod-like 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 rod-like fillers in a shear-thinning viscous fluid.

The dynamics of rod-like 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 micro-fibers 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 Rod-like 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 3D-printing of hydrogels with rod-like fillers, which can potentially be used in tissue engineering approaches and biofabrication. In the field of biofabrication, the ink made of a hydrogel-based 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 rod-shaped fillers can improve material properties like the stability during the printing process or in the created biomaterial.8,9 In the printed structure, rod-shaped 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 rod-like 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 rod-like 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 rod-like fillers in a confined flow has been investigated in several experimental studies.10,13–17

In many molding or 3D-printing processes, the melt or hydrogel flows through a contraction as it enters a circular or slot-shaped nozzle, a narrow mold, or the needle of a 3D-printer. 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 rod-like fillers in a viscous, non-Newtonian matrix, flowing from a broad to a narrow channel through a sudden contraction. At the beginning, rod-like 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.

FIG. 1.

Conceptual visualization of a viscous fluid with rod-like fillers flowing through a contraction. Rods are set at rather close distance to make the structure more perceptible.

FIG. 1.

Conceptual visualization of a viscous fluid with rod-like fillers flowing through a contraction. Rods are set at rather close distance to make the structure more perceptible.

Close modal

It is known that individual rods in a shear flow of a fluid tumble with time forming so-called 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 coarse-grained 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 particle-based 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 rod-like 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 shear-thinning 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 self-consistent 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 = γ ̇ / D r with shear rate γ ̇ and rotational diffusion constant Dr 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 non-Newtonian 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 log-rolling.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 steady-state flow field obtained for a shear-thinning fluid. Rod induced shear-thinning 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 γ ̇ = u x y. In this special case, the center of the rod drifts with a constant velocity v0 in the x direction, while the rod axis tumbles periodically with an angle ϕ. For ϕ = 0, the rod axis is parallel to the x axis. The same orientation occurs at ϕ = π and then the tumbling repeats. If the rotation from ϕ = 0 to ϕ = π takes a time T, then the rod moves a length λ = 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 rod-like 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 x0 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 rod-like 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 yc to the channel center. Therefore, rods with different yc 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.

We consider that the rod dynamics is controlled by the given flow field. The solvent is incompressible, and inertial effects of the rods are negligibly small. For the contraction geometry, we calculate the flow field with OpenFoam,36 a frequently used program for computational fluid dynamics. Both, polymers used in molding processes and hydrogels deployed in bioprinting, are typically shear-thinning fluids; we consider a fluid material with a power-law relation for the viscosity μ,
μ = K | γ ̇ | n 1
(1)
with the shear rate γ ̇, the flow consistency index K, and the flow behavior index n. Since the viscosity depends on γ ̇, which is space-dependent in our case, the Navier–Stokes equation of an incompressible fluid reads
ρ D u D t = · ( μ ( r ) u ( r ) ) P ( r ) .
(2)
In this section, we provide a simple approach for a spherocylindrical rod. Our derivation follows roughly the concept used by Dhont and Briels37 to derive the motion of a stiff chain of spheres. The rod motion results from the surface forces applied by the matrix on the rod. For simplicity, we consider only surface forces on the cylindrical part with diameter D and length L, see Fig. 2. At a time t, the center of mass r c ( t ) moves with a velocity v c ( t ), and the rod axis is parallel to a unit vector, the so-called director n ( t ). The angular velocity is given by Ω 0 ( t ). Then, the velocity at some point r on the rod is
v ( r ) = v c + Ω 0 × ( r r c ) .
(3)
FIG. 2.

Scheme for the calculation of force and torque on a rod with velocity v ( r ) at surface points r in a flow field, which would be u ( r ) in the absence of the rod.

FIG. 2.

Scheme for the calculation of force and torque on a rod with velocity v ( r ) at surface points r in a flow field, which would be u ( r ) in the absence of the rod.

Close modal

Inertial effects shall be negligible, so that the total force F = m v ̇ c and the total torque M 0 = Θ Ω ̇ 0 must vanish. Dynamic equations for v c and Ω 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.

A rod moving with director n and velocity v in a flow field u is exposed to a force field37–40 
F fric = Ξ · ( u v ) ,
(4)
where the friction tensor Ξ has the components
Ξ i j = ( Ξ n i n j + Ξ ( δ i j n i n j ) )
(5)
with the Kronecker symbol δij. This considers the fact that (u- v) can be separated into a part parallel and a part perpendicular to the rod axis which make different contributions to the force, weighted by Ξ and Ξ .
For the case of an arbitrary moving particle in an arbitrary flow field, we use the same approach to define the local force density f at a point r on the rod surface
f ( r ) = Q · ( u ( r ) v ( r ) )
(6)
with a friction density tensor Q defined via
Q i j = ( c n i n j + c ( δ i j n i n j ) ) .
(7)
We neglect the force densities at the spherical caps, which is reasonable for suitably small D/L. Integrating f over the cylinder surface S provides the force
F = S f ( r ) d 2 r
(8)
and the torque
M 0 = S ( r r c ) × f ( r ) d 2 r
(9)
with respect to the rod center r c.
We define a vector r s : = r c + s n, which lies on the rod axis, a distance | s | from r c. A ring shaped path on the rod surface around r s is parametrized by r ( s , ψ ) : = r s + ρ ( ψ ), where ρ ( ψ ) : = ρ cos ( ψ ) e 2 + ρ sin ( ψ ) e 3 with 0 ψ < 2 π. The unit vectors e 2 and e 3 are perpendicular to n and to each other. The integrals over the cylinder surface become
F = L 2 L 2 0 2 π f ( r ( s , ψ ) ) ρ d ψ ds .
(10)
If u ( r ) varies little over lengths ρ, one has
u ( r s + ρ ) u ( r s ) + ρ · ( u ) ( r s )
(11)
and Eq. (10) becomes
F = 2 π ρ Q · L 2 L 2 ( u ( r s ) v ( r s ) ) d s .
(12)
Equation (3) leads to L / 2 L / 2 v ( r s ) d s = L v c. With F = 0, one finds
v c = 1 L L / 2 L / 2 u ( r s ) d s .
(13)
The rod can rotate around its axis n, but in the following, we are just interested in the dynamics of n itself. Therefore, we just consider
Ω : = ( 1 n n ) · Ω 0
(14)
the part of the angular velocity vector Ω 0 that is perpendicular to n. For calculating the dynamics of Ω, we introduce analogously
M : = ( 1 n n ) · M 0 .
(15)
Just like M 0 also M must be zero, which leads to the equation
M = L 2 L 2 0 2 π m ( s , ψ ) ρ d ψ ds = 0
(16)
with
m ( s , ψ ) : = ( 1 n n ) · [ ( r ( s , ψ ) r c ) × f ( r ( s , ψ ) ) ] .
(17)
Using Eqs. (6), (11), and (17), the integral over ψ in Eq. (16) can be calculated, which leads to
M = c r ( ( 1 + κ 2 ) Ω n × ( I 1 I 2 ) )
(18)
including two integrals
I 1 = 12 L 3 L / 2 L / 2 ( s u ( r s ) ) d s ,
(19)
I 2 = κ 2 L n · L / 2 L / 2 u ( r s ) d s
(20)
and two parameters
κ 2 : = 3 2 D 2 L 2 c c ,
(21)
c r : = 2 π ρ L L 2 12 c .
(22)

The quantity κ L / D depends sensitively on the shape of the rod. Jeffery found κ L / D = 1 for ellipsoids of revolution.18 Using c / c = 1 / 2 (Ref. 41), Eq. (21) leads to κ L / D = 3 / 2, which for L / D = 8 ± 2 is in the same range as the values found by Cox for cylindrical rods.42 The quantity cr matches with the rotational friction coefficient.37 

With M = 0 one ends up with
Ω = 1 1 κ 2 n × ( I 1 I 2 ) .
(23)
If D/L is small enough, one can set 1 1 κ 2 1.
In this article, we use Eqs. (13) and (23) to determine the dynamics of the rods. We just mention at this point the frequently considered case of a linear flow field for which u is constant. Then, the flow field can be written as u ( r ) = G · r with the fixed gradient velocity tensor G : = ( u ) T. Inserting u ( r ) into Eqs. (19), (20), and (23), one gets
Ω = 1 1 κ 2 n × ( G · n κ 2 ( G T · n ) ) .
(24)
This is the equation found by Jeffery,18 which is valid for all linear flow fields.

1. Approximations for the integrals

For numerical calculations, the integrals in Eqs. (13), (19), and (20) are impractical. In many cases, the undisturbed flow field u ( r c + s s ) with | s | < L / 2 can be approximated very well by a fourth order polynomial. For all fourth order polynomials p ( s ) = k = 0 4 a k s k, the following two equations are valid:
1 L L / 2 L / 2 p ( s ) d s = 4 9 p ( 0 ) + 5 18 [ p ( s 1 ) + p ( s 1 ) ] ,
(25)
1 L L / 2 L / 2 s p ( s ) d s = 5 18 s 1 [ p ( s 1 ) p ( s 1 ) ]
(26)
with s 1 : = 3 20 L. With a suitably smooth function u ( r ), the expressions Eqs. (13) and (23) can be approximated by
v c = 4 9 u ( r c ) + 5 18 [ u ( r + ) + u ( r ) ] ,
(27)
Ω = n × [ 10 3 L 2 s 1 [ u ( r + ) u ( r ) ] 4 9 κ 2 u ( r c ) · n 5 18 κ 2 [ u ( r + ) + u ( r ) ] · n ]
(28)
with
r ± = r c ± s 1 n .
(29)
First, we study a power-law fluid with flow behavior index n, flowing in the x direction with a steady-state flow field u = u x ( y ) e x in a channel with fixed, plane-parallel walls at y = ± w. With Eq. (1) and shear rate γ ̇ = u x y, Eq. (2) becomes
y ( K | u x y | n 1 u x y ) = P x ,
(30)
where the term on the right side is negative and constant. The solution of Eq. (30) with no-slip boundaries at the walls u x ( ± w ) = 0 is given by
u x ( y ) = u max ( 1 | y w | m )
(31)
with
m = n + 1 n ,
(32)
u max = n n + 1 w ( n w K | P x | ) 1 / n .
(33)
For a Newtonian fluid (n = 1), one gets the Poiseuille flow field
u lin ( y ) : = u max ( 1 y 2 w 2 ) .
(34)
A fluid with a viscosity μ | γ ̇ | 1 / 2 forms the steady-state flow field
u flow ( y ) : = u max ( 1 | y w | 3 ) .
(35)
Flow fields are shown in Fig. 3 in the range w / 2 y w / 2.
FIG. 3.

Flow velocities in a planar channel of width 2w in the range w / 2 y w / 2. (A) Flow field u x ( y ) for a power-law fluid with flow behavior index n = 1/2. This corresponds to the maximum velocity of rod centers with yc = y. The minimum velocities of rods with κ 2 = 1 / 108 and centers at y are shown for rod lengths (C) L = 0.5 w, (D) L = 0.7 w, and (E) L = 0.9 w. (B) The velocity profile of a Newtonian fluid.

FIG. 3.

Flow velocities in a planar channel of width 2w in the range w / 2 y w / 2. (A) Flow field u x ( y ) for a power-law fluid with flow behavior index n = 1/2. This corresponds to the maximum velocity of rod centers with yc = y. The minimum velocities of rods with κ 2 = 1 / 108 and centers at y are shown for rod lengths (C) L = 0.5 w, (D) L = 0.7 w, and (E) L = 0.9 w. (B) The velocity profile of a Newtonian fluid.

Close modal
We consider a rod in the x, y plane. The velocity v c and the angular velocity Ω of a rod in the planar channel result from Eqs. (13) and (23). The rod center flows in the x direction v c = v c e x, while Ω = Ω e z. Let ϕ be the angle between the rod axis and the x axis. In general, Ω and vc have contributions that oscillate with ϕ. For ϕ = 0, vc is maximum, while Ω is a small but finite value. For the Newtonian fluid, one has
Ω ( ϕ = 0 ) = 2 u max κ 2 w 2 y c ,
for the shear-thinning fluid with n = 1/2, one finds
Ω ( ϕ = 0 ) = 3 u max κ 2 w 3 | y c | y c .

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 κ 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.

Inserting Eqs. (35) and (34) into Eqs. (13) and (23), one reveals the velocity of the rod center v c , x and the angular velocity Ω as a function of ϕ. Setting q ( ϕ ) : = L 2 12 y c 2 sin 2 ( ϕ ), we find
v c , x ( ϕ ) = u x ( y ) y c 2 w 2 q ( ϕ ) ,
(36)
Ω ( ϕ ) = 2 u max y c w 2 ( κ 2 cos 2 ( ϕ ) + sin 2 ( ϕ ) )
(37)
for the Newtonian fluid. For the shear-thinning fluid, simple expressions are found if y c > L / 2. Then, one has
v c , x ( ϕ ) = u x ( y ) 3 y c 3 w 3 q ( ϕ ) ,
(38)
Ω ( ϕ ) = 3 u max y c 2 w 3 [ [ 1 + q ( ϕ ) ] κ 2 cos 2 ( ϕ ) + [ 1 + 3 5 q ( ϕ ) ] sin 2 ( ϕ ) ] .
(39)

The rod length appears only in q ( ϕ ), a positive term that oscillates with sin ( ϕ ) and is proportional to L2. For rods in the Newtonian fluid, the velocity is reduced by a term proportional to q ( ϕ ), while, interestingly, Ω ( ϕ ) does not depend on L at all, as long as the rod does not collide with the channel wall.

Also, for rods in the shear-thinning fluid, the velocity is decreased by a term proportional to q ( ϕ ). The angular velocity Ω is increased by a term proportional to q ( ϕ ). The lowest velocities are always found for ϕ ( x ) = π / 2 , 3 π / 2 , Examples of v x ( ϕ ) and Ω ( ϕ ) of a rod in a planar channel are shown in Fig. 4.

FIG. 4.

Motion of rods in a planar channel: (i) Angular velocities Ω and (ii) velocities v c , x of rod centers. The results are shown for rods with κ 2 = 1 / 108 for different lengths and different positions. The matrix is either shear-thinning with flow behavior index n = 1/2 or a viscous Newtonian fluid. Values are shown in length units of w and velocity units umax.

FIG. 4.

Motion of rods in a planar channel: (i) Angular velocities Ω and (ii) velocities v c , x of rod centers. The results are shown for rods with κ 2 = 1 / 108 for different lengths and different positions. The matrix is either shear-thinning with flow behavior index n = 1/2 or a viscous Newtonian fluid. Values are shown in length units of w and velocity units umax.

Close modal
One important characteristic of the rod dynamics is the distance λ : = x ( ϕ = π ) x ( ϕ = 0 ) that the rod travels as it fulfills a half rotation. Values of x ( ϕ ) are given by
x = 0 ϕ v ( ϕ ̃ ) Ω ( ϕ ̃ ) d ϕ ̃ .
For Newtonian fluids, one finds
λ = π 2 y 0 ( w 2 y 0 2 κ L 2 12 ( 1 + κ ) ) .

For a shear-thinning matrix, the integral is more complicated. In the following, a case of special interest is a rod in a shear-thinning matrix with n = 1/2 at height y c = w / 2 with rod length L w.

For this case, λ is approximately
λ 11 π 4 6 κ L .
(40)

In general, λ decreases strongly with yc. The rotation dynamics in a channel is visualized in Fig. 5, which shows cos 2 ( ϕ ) ( x c , y c ) for rods with length L = 0.9 w and κ 2 = 1 / 108 that start with ϕ = 0 at xc = 0.

FIG. 5.

Values of cos 2 ( ϕ ( x , y ) ) for angle ϕ between the rod axis and the x direction. Rods start at ( 0 , y ) with ϕ = 0 and drift in a shear-thinning fluid with n = 1/2 flowing through a planar channel.

FIG. 5.

Values of cos 2 ( ϕ ( x , y ) ) for angle ϕ between the rod axis and the x direction. Rods start at ( 0 , y ) with ϕ = 0 and drift in a shear-thinning fluid with n = 1/2 flowing through a planar channel.

Close modal

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 xco where the width reduces spontaneously from 2 w 0 to 2w. A steady-state flow field was created with the help of the program OpenFoam for a system with 2 w 0 = 2.4 mm , 2 w = 0.8 mm , x c o = 4 mm = 10 w, and a kinematic pressure decreasing with Δ P / Δ x = 1125 ms 2. The fluid velocity at the walls is u = 0. The kinematic viscosity of the shear-thinning matrix obeys the power-law equation (1) with flow behavior index n = 1/2 and flow consistency index K = 5.67 × 10 3 m 2 s n 1. The x and y components of the stationary flow field are shown in Fig. 6.

FIG. 6.

Flow field around the compression. The plots show (a) the x component and (b) the y component of the flow field. In the plots, u0 is the maximum absolute flow velocity in the whole region, and w is half the diameter of the small channel.

FIG. 6.

Flow field around the compression. The plots show (a) the x component and (b) the y component of the flow field. In the plots, u0 is the maximum absolute flow velocity in the whole region, and w is half the diameter of the small channel.

Close modal

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 K, and a mass density ρ = 1 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 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 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 κ 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 self-evident. We consider that in a region x < 0, rods have a random, homogenous orientational and spatial distribution. Thus, at each height yc, rod centers have the same average distance Δ x in the x direction. In a time interval Δ t, rods appear at xc = 0 with an average frequency f ( y c ) = u x ( y c ) / Δ x. The amount of rods n 0 ( y c ) Δ t per time Δ t is proportional to u x ( 0 , y c ) 1 ( 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 ( y c / w 0 ) 3.

Before we look at a system with random initial orientation in 3D, we discuss some aspects of the rod dynamics with initial rod axes lying in the x, y plane. In Fig. 7, paths of rods through the contraction are shown. One can see that rods that start above a height y ̃ c , 0 w 0 / 2 end up above a height y c ̃ : = w L tot / 2. With the given rate of rods n 0 ( y c ) that pass xc = 0 at height yc, the fraction of rods that end up above y ̃ c is about
y ̃ c , 0 y c , max n 0 ( y c ) d y c 0 y c , max n 0 ( y c ) d y c 0.32 .
FIG. 7.

Paths of centers of rods, passing a contraction. Rods starting at y c > y ̃ c , 0 end up above y c ̃ in the narrow channel and collide with the upper wall. Thick brown lines represent the upper wall. Dotted brown lines indicate the second half of the contraction.

FIG. 7.

Paths of centers of rods, passing a contraction. Rods starting at y c > y ̃ c , 0 end up above y c ̃ in the narrow channel and collide with the upper wall. Thick brown lines represent the upper wall. Dotted brown lines indicate the second half of the contraction.

Close modal

Approximately 32% of the rods end up in the upper region y c > y ̃ 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 ̃ 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 ̃ 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 ux and γ ̇.

FIG. 8.

Three rods passing the contraction. The length of the colored lines indicates the full length Ltot of the rods. The thick brown line denotes the wall. The motion of the rods continues in the lower plot.

FIG. 8.

Three rods passing the contraction. The length of the colored lines indicates the full length Ltot of the rods. The thick brown line denotes the wall. The motion of the rods continues in the lower plot.

Close modal

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 w 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 | y c | 0.55 w. The orientation of the rods is reflected in the spatial distribution n tip ( x , y ) of cylinder ends r c ± L 2 n, shown in Fig. 10.

FIG. 9.

Logarithmic contour plot of the density n c ( x , y ) of rod centers near the contraction. The fat brown line denotes the wall.

FIG. 9.

Logarithmic contour plot of the density n c ( x , y ) of rod centers near the contraction. The fat brown line denotes the wall.

Close modal
FIG. 10.

Contour plot of the probability density n tip ( x , y ) of the cylinder ends r c ± L 2 n of rods passing a contraction in a shear-thinning matrix. Thick brown lines indicate the walls.

FIG. 10.

Contour plot of the probability density n tip ( x , y ) of the cylinder ends r c ± L 2 n of rods passing a contraction in a shear-thinning matrix. Thick brown lines indicate the walls.

Close modal

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 ± y ̃, which leads to the maxima of ntip at ± y ̃.

Figure 11 shows the diagonal elements Axx and Ayy of the orientation tensor
A = n n ,
which is frequently used in the field of rod distributions.26,44 The component A x x ( x c ) = n x n x x c provides information on the average alignment in x-direction for rods with rod centers in the range ( x c , x c + Δ x ) with small Δ x. The same applies for A y y ( x c ) and A z z ( x c ). The plot in Fig. 11 shows values of Axx and Ayy (dashed lines), averaged over the whole range w < y c < w and A x x m and A y y m (continuous lines), averaged over rods with | y c y ̃ c | < 0.01 w, which belong to the maxima in Fig. 9.
FIG. 11.

Diagonal elements Axx (blue) and Ayy (red) of the orientation tensor. Orientations are sampled over the whole planar channel (dashed lines) or in the range | y c y ̃ c | < 0.01 w ( A x x m , A y y m, continuous lines).

FIG. 11.

Diagonal elements Axx (blue) and Ayy (red) of the orientation tensor. Orientations are sampled over the whole planar channel (dashed lines) or in the range | y c y ̃ c | < 0.01 w ( A x x m , A y y m, continuous lines).

Close modal

At the starting point, x c 0, the quantities Axx, Ayy, and A z z = 1 A x x A y y are 1/3, which correspond to an isotropic distribution. At x c x c o, one has A x x 0.78. This indicates a distinct alignment in the flow direction. At fixed distances λ, Ayy 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 ± y ̃ 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 γ ̇ ( y ) = u x y ( y ) and tumble with an angular velocity perpendicular to the x, y plane. If their center of mass is above y ̃ c, they hit the wall and end up at y c y ̃ c. The large amount of rods with y c = y ̃ c is exposed to the same shear flow γ ̇ ( y ̃ c ) and so they all rotate with the same angular velocity sequence. The rod alignment at the contraction leads to ϕ 0 at this point so that many rods follow the same path with the same orientation sequence ϕ ( x c ). The rods at y c = y ̃ c contact the channel walls at regular distances, when the rods are perpendicular to the flow direction ϕ ( x ) = π / 2.

We have investigated the behavior of rods entering a contraction for the case that the fluid is shear-thinning. This was chosen because in many applications shear-thinning fluids are used. At this point, it must be emphasized that shear-thinning 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 non-Newtonian case. One can see that the synchronized rotation of the rods is visible just like in the non-Newtonian case.

FIG. 12.

Contour plot of the probability density n tip ( x , y ) of the cylinder ends r c ± L 2 n of rods passing a contraction in a Newtonian fluid. Thick brown lines indicate the walls.

FIG. 12.

Contour plot of the probability density n tip ( x , y ) of the cylinder ends r c ± L 2 n of rods passing a contraction in a Newtonian fluid. Thick brown lines indicate the walls.

Close modal

Tumbling of rods in a viscous fluid that flows through a channel depends on the distance from the channel center yc and the orientation at an initial point xco. 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 λ with n = 0 , 1 , 2 , .

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 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 Ltot 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 rod-wall 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.

This project is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project number 326998133 – SFB/TRR225 (subproject B03-PI Sahar Salehi).

The authors have no conflicts to disclose.

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).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
L.
Jiang
,
Y.
Zhou
, and
F.
Jin
, “
Design of short fiber‐reinforced thermoplastic composites: A review
,”
Polym. Compos.
43
,
4835
4847
(
2022
).
2.
A.
Bernasconi
and
F.
Cosmi
, in
11th International Conference on the Mechanical Behavior of Materials (ICM11)
,
2011
.
3.
Y.
Wang
,
Z.
Wang
, and
L.
Zhu
, “
A short review of recent progress in improving the fracture toughness of FRP composites using short fibers
,”
Sustainability
14
,
6215
(
2022
).
4.
D.
Sonnleitner
,
S.
Schrufer
,
L.
Berglund
,
D. W.
Schubert
, and
G.
Lang
, “
Correlating rheology and printing performance of fiber-reinforced bioinks to assess predictive modelling for biofabrication
,”
J. Mater. Res.
36
,
3821
3832
(
2021
).
5.
D.
Kada
,
A.
Koubaa
,
G.
Tabak
,
S.
Migneault
,
B.
Garnier
, and
A.
Boudenne
, “
Tensile properties, thermal conductivity, and thermal stability of short carbon fiber reinforced polypropylene composites
,”
Polym. Compos.
39
,
E664
E670
(
2018
).
6.
R.
Ram
,
M.
Rahaman
,
A.
Aldalbahi
, and
D.
Khastgir
, “
Determination of percolation threshold and electrical conductivity of polyvinylidene fluoride (PVDF)/short carbon fiber (SCF) composites: Effect of SCF aspect ratio
,”
Polym. Int.
66
,
573
582
(
2017
).
7.
S.
Utech
and
A. R.
Boccaccini
, “
A review of hydrogel-based composites for biomedical applications: Enhancement of hydrogel properties by addition of rigid inorganic fillers
,”
J. Mater. Sci.
51
,
271
310
(
2016
).
8.
S.
Salehi
,
S.
Ostrovidov
,
M.
Ebrahimi
,
R. B.
Sadeghian
,
X.
Liang
,
K.
Nakajima
,
H.
Bae
,
T.
Fujie
, and
A.
Khademhosseini
, “
Development of flexible cell-loaded ultrathin ribbons for minimally invasive delivery of skeletal muscle cells
,”
ACS Biomater. Sci. Eng.
3
,
579
589
(
2017
).
9.
M. E.
Prendergast
,
M. D.
Davidson
, and
J. A.
Burdick
, “
A biofabrication method to align cells within bioprinted photocrosslinkable and cell-degradable hydrogel constructs via embedded fibers
,”
Biofabrication
13
,
044108
(
2021
).
10.
A.
Kosik-Koziol
,
M.
Costantini
,
T.
Bolek
,
K.
Szoke
,
A.
Barbetta
,
J.
Brinchmann
, and
W.
Swieszkowski
, “
PLA short sub-micron fiber reinforcement of 3D bioprinted alginate constructs for cartilage regeneration
,”
Biofabrication
9
,
044105
(
2017
).
11.
S.
Das
and
B.
Basu
, “
Extrusion‐based 3D printing of gelatin methacryloyl with nanocrystalline hydroxyapatite
,”
Int. J. Appl. Ceram. Technol.
19
,
924
938
(
2022
).
12.
L.
Ren
,
X.
Zhou
,
Q.
Liu
,
Y.
Liang
,
Z.
Song
,
B.
Zhang
, and
B.
Li
, “
3D magnetic printing of bio-inspired composites with tunable mechanical properties
,”
J. Mater. Sci.
53
,
14274
14286
(
2018
).
13.
A. S.
Khair
, “
Taylor dispersion of elongated rods at small and large rotational Péclet numbers
,”
Phys. Rev. Fluids
7
,
014502
(
2022
).
14.
J. C.
Rose
,
M.
Camara-Torres
,
K.
Rahimi
,
J.
Koehler
,
M.
Moeller
, and
L.
De Laporte
, “
Nerve cells decide to orient inside an injectable hydrogel with minimal structural guidance
,”
Nano Lett.
17
,
3782
3791
(
2017
).
15.
X.
Wang
,
Q.
Wang
, and
C.
Xu
, “
Nanocellulose-based inks for 3D bioprinting: Key aspects in research development and challenging perspectives in applications—A mini review
,”
Bioengineering
7
,
40
(
2020
).
16.
M. D.
Davidson
,
M. E.
Prendergast
,
E.
Ban
,
K. L.
Xu
,
G.
Mickel
,
P.
Mensah
,
A.
Dhand
,
P. A.
Janmey
,
V. B.
Shenoy
, and
J. A.
Burdick
, “
Programmable and contractile materials through cell encapsulation in fibrous hydrogel assemblies
,”
Sci. Adv.
7
,
eabi8157
(
2021
).
17.
T.
Kim
,
R.
Trangkanukulkij
, and
W. S.
Kim
, “
Nozzle shape guided filler orientation in 3D printed photo-curable nanocomposites
,”
Sci. Rep.
8
,
3805
(
2018
).
18.
G. B.
Jeffery
, “The motion of ellipsoidal particles immersed in a viscous fluid,”
Proc. R. Soc. A
102
,
161
179
(
1922
).
19.
J. H.
Park
and
Y. L.
Joo
, “
Tailoring nanorod alignment in a polymer matrix by elongational flow under confinement: Simulation, experiments, and surface enhanced Raman scattering application
,”
Soft Matter
10
,
3494
3505
(
2014
).
20.
G. N.
Toepperwein
,
N. C.
Karayiannis
,
R. A.
Riggleman
,
M.
Kroeger
, and
J. J.
de Pablo
, “
Influence of nanorod inclusions on structure and primitive path network of polymer nanocomposites at equilibrium and under deformation
,”
Macromolecules
44
,
1034
1045
(
2011
).
21.
Y.
Tao
,
W.
den Otter
, and
W.
Briels
, “
Kayaking and wagging of rods in shear flow
,”
Phys. Rev. Lett.
95
,
237802
(
2005
).
22.
J.
Zhang
,
X.
Xu
, and
T.
Qian
, “
Anisotropic particle in viscous shear flow: Navier slip, reciprocal symmetry, and Jeffery orbit
,”
Phys. Rev. E
91
,
033016
(
2015
).
23.
D.
Palanisamy
and
W. K.
den Otter
, “
Efficient Brownian dynamics of rigid colloids in linear flow fields based on the grand mobility matrix
,”
J. Chem. Phys.
148
,
194112
(
2018
).
24.
N.
Meyer
,
O.
Saburow
,
M.
Hohberg
,
A. N.
Hrymak
,
F.
Henning
, and
L.
Kaerger
, “
Parameter identification of fiber orientation models based on direct fiber simulation with smoothed particle hydrodynamics
,”
J. Compos. Sci.
4
,
77
(
2020
).
25.
K.
Wu
,
L.
Wan
,
H.
Zhang
, and
D.
Yang
, “
Numerical simulation of the injection molding process of short fiber composites by an integrated particle approach
,”
Int. J. Adv. Manuf. Technol.
97
,
3479
3491
(
2018
).
26.
S. K.
Kugler
,
A.
Kech
,
C.
Cruz
, and
T.
Osswald
, “
Fiber orientation predictions—A review of existing models
,”
J. Compos. Sci.
4
,
69
(
2020
).
27.
F.
Folgar
and
C. L.
Tucker
, “Orientation behavior of rigid fibers in concentrated suspensions,”
J. Rheol.
26
,
604
(
1982
).
28.
S. A.
Abtahi
and
G. J.
Elfring
, “
Jeffery orbits in shear-thinning fluids
,”
Phys. Fluids
31
,
103106
(
2019
).
29.
A.
Zoettl
,
K. E.
Klop
,
A. K.
Balin
,
Y.
Gao
,
J. M.
Yeomans
, and
D. G. A. L.
Aarts
, “
Dynamics of individual Brownian rods in a microchannel flow
,”
Soft Matter
15
,
5810
5814
(
2019
).
30.
J.
Einarsson
,
F.
Candelier
,
F.
Lundell
,
J. R.
Angilella
, and
B.
Mehlig
, “
Rotation of a spheroid in a simple shear at small Reynolds number
,”
Phys. Fluids
27
,
063301
(
2015
).
31.
V.
Calabrese
,
S. J.
Haward
, and
A. Q.
Shen
, “
Effects of shearing and extensional flows on the alignment of colloidal rods
,”
Macromolecules
54
,
4176
4185
(
2021
).
32.
Z.
Cui
,
A.
Dubey
,
L.
Zhao
, and
B.
Mehlig
, “
Alignment statistics of rods with the Lagrangian stretching direction in a channel flow
,”
J. Fluid Mech.
901
,
A16
(
2020
).
33.
A.
Ahmed
,
M.
Mansouri
,
I. M.
Joshi
,
A. M.
Byerley
,
S. W.
Day
,
T. R.
Gaborski
, and
V. V.
Abhyankar
, “
Local extensional flows promote long-range fiber alignment in 3D collagen hydrogels
,”
Biofabrication
14
,
035019
(
2022
).
34.
M.
Trebbin
,
D.
Steinhauser
,
J.
Perlich
,
A.
Buffet
,
S. V.
Roth
,
W.
Zimmermann
,
J.
Thiele
, and
S.
Foerster
, “
Anisotropic particles align perpendicular to the flow direction in narrow microchannels
,”
Proc. Natl. Acad. Sci. U. S. A.
110
,
6706
6711
(
2013
).
35.
J. P.
Lewicki
,
J. N.
Rodriguez
,
C.
Zhu
,
M. A.
Worsley
,
A. S.
Wu
,
Y.
Kanarska
,
J. D.
Horn
,
E. B.
Duoss
,
J. M.
Ortega
,
W.
Elmer
,
R.
Hensleigh
,
R. A.
Fellini
, and
M. J.
King
, “
3D-printing of meso-structurally ordered carbon fiber/polymer composites with unprecedented orthotropic physical properties
,”
Sci. Rep.
7
,
43401
(
2017
).
36.
H.
Weller
,
G.
Tabor
,
H.
Jasak
, and
C.
Fureby
, “
A tensorial approach to computational continuum mechanics using object-oriented techniques
,”
Comput. Phys.
12
,
620
631
(
1998
).
37.
J.
Dhont
and
W.
Briels
, in
Soft Matter Volume 2: Complex Colloidal Suspensions
, 1st ed., edited by
G.
Gompper
and
M.
Schnick
(
Wiley-VCH Verlag GmbH & Co KGaA
,
Weinheim
,
2005
), Part 3.1–3.5, pp.
147
180
.
38.
E.
Fischermeier
,
Simulations of Colloidal Liquid Crystals
, 1st ed. (
FAU University Press
,
Erlangen
,
2016
).
39.
M. M.
Tirado
and
J. G.
Garcia de la Torre
, “
Translational friction coefficients of rigid, symmetric top macromolecules. Application to circular cylinders
,”
J. Chem. Phys.
71
,
2581
2587
(
1979
).
40.
H.
Yamakawa
and
G.
Tanaka
, “
Translational diffusion coefficients of rodlike polymers: Application of the modified Oseen tensor
,”
J. Chem. Phys.
57
,
1537
(
1972
).
41.
A.
Meunier
, “
Friction coefficient of rod-like chains of spheres at very low Reynolds numbers—II: Numerical simulations
,”
J. Phys. II France
4
,
561
566
(
1994
).
42.
R.
Cox
, “
The motion of long slender bodies in a viscous fluid—Part 2: Shear flow
,”
J. Fluid Mech.
45
,
625
(
1971
).
43.
H.
Kobayashi
and
R.
Yamamoto
, “
Reentrant transition in the shear viscosity of dilute rigid-rod dispersions
,”
Phys. Rev. E
84
,
051404
(
2011
).
44.
C.
Tucker
and
F.
Folgar
, “
A model of compression mold filling
,”
Polym. Eng. Sci.
23
,
69
73
(
1983
).
Published open access through an agreement with Universität Bayreuth