Direct Numerical Simulation of Concentration and Orientation Distribution of Fibers in a Mixing Layer

Abstract

The concentration and orientation of suspended fibers in a mixing layer are investigated numerically. Two cases (diffusive and nondiffusive) are investigated for the fiber concentration distribution. The fine structures of the instantaneous distributions under these two cases are very different due to molecular diffusion. Sharp front of concentration is observed in the nondiffusive case. However, there is no obvious difference in the mean concentration between the two cases. With regard to the orientation, a fiber may rotate periodically or approach an asymptotic orientation, which is determined by a determinant defined with the stain rate. The symmetric part of the strain rate tends to make a fiber align to an asymptotic orientation, while the antisymmetric part drives a fiber to rotate. When a fluid parcel passes through a region with relatively high shear rate, fibers carried by the fluid parcel are most likely to rotate incessantly. On the other hand, in the region of relatively high extension rate, fibers tend to align to some asymptotic orientation. Generally, fibers tend to align with the shear plane. This fact has significant implications in predicting the rheological properties of fiber suspension flows.

1. Introduction

Fiber suspension flows can be found in many processes, such as papermaking, polymer flows in melt-blowing extruders, and nanofibers in the respiratory system. Here, discussions are limited to rigid fibers, which are very slender bodies and can be treated as high aspect ratio cylinders or ellipsoids. The aspect ratio 𝑟_𝑐 is generally defined as the ratio of the maximum and the minimum characteristic sizes of a fiber. Because of the high aspect ratio, fiber suspensions generally exhibit anisotropic properties. From a macroscopic point of view, fiber suspension flows are generally nonNewtonian fluid, and the flows are described by the Navier-Stokes equation with appropriate constitutive equations for the strain rate and stress. In this scale, the effects of fiber additive are treated through ensemble average. In other words, continuous fields are used to describe the fibers states, most importantly fiber concentration and orientation distribution. From a microscopic point of view, the flow around a single fiber and the corresponding translational and rotational dynamics of a fiber are of concern. These two perspectives are integrated to provide a continuous constitutive equation for the strain rate and the stress through the statistical average on all the possible microscopic flow structures. This two-scale perspective leads the research in fiber suspension flow to focus on two main aspects.

One aspect is to investigate the movement of fibers in various flows. Since fibers are generally very small, the flow around a fiber can usually be seen as a creeping flow in the moving coordinates aligned with the fiber’s mass center. Other than the translational movement, the rotation of a fiber has received investigation long ago [1]. Since then, a lot of research works have been done on the orientation of fibers in various flows [2–15]. The other aspect is to investigate the new properties of flows caused by the addition of fibers. Intrinsically, it is the additional stress induced by the fiber suspensions that makes the fiber suspensions flows exhibit special rheological properties [16, 17], such as the drag reducing [18, 19], shear thinning [20]. In his serial works [21–23], Batchelor developed a well-accepted model about the additional stress, which closely relies on the fiber orientation distribution.

In this work, the translational and rotational movements of fibers in a canonical mixing layer are simulated. The mixing layer configuration is geometrically simple and is a very broadly used model to investigate the shear flow, which characterizes the most important fluid dynamics. A Lagrangian particles scheme is used to deal with the convection of fibers. The evolution of fiber orientation is tracked along Lagrangian trajectories. This new scheme provides insightful understanding of the fiber rotational dynamics, which helps understand and predict the rheological properties of fiber suspension flow. Meanwhile, it is computationally efficient and highly flexible in adjusting the discretization error on the orientation distribution. The paper is organized as follows. Models and methods are described in Section 2. The results on fiber concentration and orientation are presented in Section 3. Section 4 is the conclusion.

2. Models and Methods

2.1. Mixing Layer Configuration

This work is to simulate the concentration and orientation distribution of fibers in a canonical mixing layer. To save the computational time, the flow is assumed to be homogeneous in the spanwise direction. On the other hand, the feedback from the fiber additive is neglected; only one-way coupling is considered. Hence, a planar mixing is solved. The configuration is depicted in Figure 1. The spacial dimension is normalized by the momentum thickness of the mixing layer at the inlet. The computational methods for the fluid dynamics are the same as those in Zhou and He [24]. Finite difference scheme with structured grid is used. At the inlet, the streamwise velocity profile is imposed by combining two Blasius laminar boundary layers. The two streams have inlet velocities 𝑢₁=15m/s and 𝑢₂=5m/s, respectively. The Reynolds number at the inlet based on the velocity difference and the momentum thickness is 111. The normalized fiber concentrations at the two streams are 𝑐₁=1 and 𝑐₂=0, respectively.

Figure 1 

Schematic view of the mixing layer, with instantaneous fiber concentration 𝑐 in the mixing layer for the case Sc=1.

2.2. Fiber Concentration

It is assumed that the fiber concentration 𝑐 is a passive scalar which satisfies the convection diffusion equation 𝜕𝑐()𝜕𝑡+∇⋅𝐮𝑐=𝐷_𝑐∇²𝑐. (1) The assumption is valid when the fiber is very small and follows the fluid flow very well; additionally, the fiber concentration is very low and the additional stress due to the fiber additive is negligible. The main purpose is to investigate the Schmidt number effect on the fiber concentration in the mixing layer configuration here. The Schmidt number is Sc=𝜈/𝐷_𝑐, where 𝜈 is the kinematic viscosity. Two cases Sc=1 and Sc=∞ (i.e., 𝐷_𝑐=0) are investigated. The concentrations in the two streams of the mixing layer at the inlet are assumed to be 0 and 1 (nondimensionalized). The absolute magnitude of the concentration is immaterial, since the fiber concentration has no impact on the fluid flow (one-way coupling).

2.3. Fiber Rotational Dynamics

The rotation of a fiber is described by the Jeffery [1] equation 𝐩̇=𝝎⋅𝐩+𝜆(𝝐⋅𝐩−𝝐∶𝐩𝐩𝐩), (2) where 𝐩 is the unit vector aligned with the fiber axis (Figure 2), the dot over a variable denotes the time derivative, 𝝎=(∇𝐮^†−∇𝐮)/2 is the vorticity tensor, 𝝐=(∇𝐮^†+∇𝐮)/2 is the deformation rate tensor, and 𝜆=(𝑟²_𝑐−1)/(𝑟²_𝑐+1). The equation is originally derived for an ellipsoid in Stokes flows. It is found [1] that the force acting on the ellipsoid reduces to two couples, with one tending to make the ellipsoid adopt the same rotation as the surrounding fluid and the other tending to set the ellipsoid with its axes parallel to the principle axes of distortion of the surrounding fluid. This discovery was verified by the subsequent experiment of Taylor [25]. Furthermore, Bretherton [26] found that the same equation (2) also applies to a revolution body.

Figure 2 

Schematics of fiber orientation. The bold line segment denotes the fiber orientation vector.

As discussed in Section 2.1, the mixing layer is assumed to be homogeneous along the spanwise direction (𝑧 in Figure 2). Hence, only components 𝜔₁₂ in the vorticity tensor and 𝜖₁₁, 𝜖₁₂, and 𝜖₂₂ in the deformation rate tensor are present to determine the rotation of a fiber. Under this condition, the Jeffery equation (2) is simplified to𝑝̇₁=𝜔₁₂𝑝₂+𝜆𝜖₁₁𝑝₁+𝜖₁₂𝑝₂−𝜖₁₁𝑝³₁−2𝜖₁₂𝑝₁²𝑝₂−𝜖₂₂𝑝₂²𝑝₁,𝑝̇₂=−𝜔₁₂𝑝₁+𝜆𝜖₁₂𝑝₁+𝜖₂₂𝑝₂−𝜖₁₁𝑝²₁𝑝₂−2𝜖₁₂𝑝₁𝑝²₂−𝜖₂₂𝑝³₂,𝑝̇₃=−𝜆𝑝₃𝜖₁₁𝑝²₁+2𝜖₁₂𝑝₁𝑝₂+𝜖₂₂𝑝²₂. (3a)(3b)(3c)The orientation vector implicitly satisfies the condition 𝑝²₁+𝑝²₂+𝑝²₃=1, since the fiber is rigid, which cannot extend or shrink.

The rotational dynamic equations (3a), (3b), and (3c) are integrated along a Lagrangian trajectory (discussed in the next section) with the local vorticity and deformation rate by the solver DOPRI5 [27], which uses a fourth order Runge-Kutta method with adaptive time step based a fifth order error estimate.

2.4. Lagrangian Particles Scheme

For the translational movement of fibers, two cases Sc=1 and Sc=∞ are considered. For Sc=1, the WENO scheme [28] is used to solved the convection diffusion equation (1). For the nondiffusive case Sc=∞, the traditional schemes in Eulerian framework usually introduce numerical diffusion. To avoid the numerical diffusion, the Lagrangian particles scheme [29] is used for the case Sc=∞. Under the conditions of no diffusion and negligible inertial, the movement of the mass center of a fiber follows the local fluid parcel. Hence, the transport of fibers can be described by the Lagrangian particles method. In the Lagrangian particles method, the trajectories of a large number of Lagrangian particles are tracked simultaneously. Meanwhile, the Lagrangian particles carry the fiber concentration and orientation, which evolve along the Lagrangian trajectories independently. The fiber concentration is a scalar, which is easy to deal with. When calculating the orientation distribution at a spatial position, which is a function on the unit sphere (Figure 2), statistics of fiber orientation at the position are required. A Monte Carlo scheme is used to simulate the orientation evolution along a trajectory. A Lagrangian particle carries many fibers, and these fibers all obey the rotational dynamics equation (2). The initial orientations of these fibers are generated randomly on the unit sphere. The orientation distribution is approximated through statistical averaging on the orientations of all these fibers. Thousands of fibers are used to obtain statistically convergent orientation distribution. The main advantage of the Lagrangian particles scheme is that the convection of all the fibers along a Lagrangian particle is solved simultaneous in the same way. The number of fibers along a trajectory can be adjusted conveniently to balance the computational cost and the statistical error.

The dynamic equations for the Lagrangian particles and the concentration and orientation of fibers are as follows: d𝐱^𝑝d𝑡=𝐯^𝑝𝐱^𝑝,𝑝=1,2,…,𝑁_𝑝,d𝑐^𝑝d𝑡=0,𝑝=1,2,…,𝑁_𝑝, (4)(5) in addition to (3a), (3b), and (3c), where the strain rate is obtained along the trajectories determined by (4). Here 𝐱^𝑝 and 𝐯^𝑝 denote the location and the velocity of a Lagrangian particle 𝑝, and 𝑁_𝑝 is the number of Lagrangian particles.

In this study, the flow velocity is solved in the traditional Eulerian framework. Fiber concentration and orientation are evolved along Lagrangian trajectories. The Eulerian and Lagrangian frameworks are coupled together. On the one hand, trilinear interpolation is used to convert variables from the Eulerian framework to the Lagrangian framework. To calculate the particle trajectories, the particle velocity 𝐯^𝑝(𝐱^𝑝) is obtained by interpolating the Eulerian velocity field at position 𝐱^𝑝. The strain rate is also obtained from the Eulerian field through interpolating it in a grid cell. On the other hand, the reconstruction of Eulerian fields from the Lagrangian particles uses𝑐(𝐱)=1𝑁_Ω𝑝∈Ω𝑐^𝑝, (6) where Ω is the set of particles located in the cell 𝐱 and 𝑁_Ω is the number of such particles. Similar treatment for the orientation is carried out.

Lagrangian particles are injected at the inlet of the mixing layer at a rate proportional to the inlet velocity, to make the number of particle 100 in a cell on average. The number of particles be in a cell controls the discretization error in evaluating the orientation distribution, which can be conveniently adjusted. Increasing the number of particles does not affect the cost of solving the fluid dynamics but only increases the cost to track the Lagrangian trajectories and the integration of fiber rotational dynamics along these trajectories. It is found that 100 particles in a cell are enough to produce statistically convergent results.

3. Results and Discussion

3.1. Fiber Concentration

The instantaneous fiber concentrations for Sc=1 and Sc=∞ are given in Figure 3. The concentration distribution for the two cases has similar large structure. For Sc=1, the concentration changes smoothly from 0 to 1. However, for Sc=∞, the concentration is either 0 or 1; no intermediate value is in-between, which is typical for nondiffusive mixing process. The Lagrangian particles scheme captures the sharp front quite well. The front separating 𝑐=0 and 1 in Figure 3(b) is not smooth but exhibits staggered structure, which is related to the grid resolution. It is clear that high resolution is needed to resolve the fine scalar structure in high Schmidt numbers flows.

(a)

(b)

Figure 3 

Comparison of instantaneous fiber concentration distributions (clipped). (a) Sc=1, (b) Sc=∞. Same isocontours of 𝑐=0.9 and 𝑐=0.1 for Sc=1 are also added in (b) for comparison purpose. The fiber concentration has a sharp front for Sc=∞.

Although the instantaneous fine structures of the fiber concentration for Sc=1 and Sc=∞ are very different, the mean values are almost the same. Figure 4 shows the mean concentration cross profiles at various streamwise locations, where 𝑦⁺=𝑦/𝛿(𝑥^∗) is the normalized coordinate with the local momentum thickness (𝛿(𝑥^∗)) of the mixing layer. The profiles a, b, and c are for Sc=1 at 𝑥^∗=137, 𝑥^∗=228, and 𝑥^∗=319, respectively. The profiles A, B, and C are for Sc=∞ at the same locations corresponding to a, b, and c. There is no significant difference between the diffusive and nondiffusive cases. Only minor difference between a and A is observed. It is clear that the molecular diffusion has negligible effect in determining the mean concentration in turbulent flows. Generally, the Schmidt number for the translational movement of a fiber is much larger than unity. Hence, the molecular diffusion can be neglected in determining the mean concentration.

Figure 4 

Mean concentration a cross profiles at various crosswise locations (𝑥∗=137, 𝑥∗=228, and 𝑥∗=319) for Sc=1 (a, b, and c) and Sc=∞ (A, B, and C).

3.2. Fiber Orientation

The orientation vector has three components in the Cartesian coordinates. However, only two of them are independent. The third one can be determined through the normalization condition 𝑝²₁+𝑝²₂+𝑝²₃=1. Equivalently, the orientation vector can be described in the spherical coordinates with 𝜙=arctan(𝑝₂/𝑝₁) and 𝜃=arccos(𝑝₃) (see Figure 2). For a given constant strain rate, Zhou et al. [30] have derived the analytical solution for the rotational movement of a fiber in a planar flow. They found that a fiber either rotates periodically or approaches an asymptotic orientation according to the sign of a determinant Δ, which is defined with the strain rate. The asymptotic orientation (if Δ>0) and the period (if Δ<0) are also given there. In their derivation the shear stain is normalized by the component 𝑢_𝑦. Their formula cannot describe the case when 𝑢_𝑦=0. Here we extend their results slightly to present the results in a more general form (including the case for 𝑢_𝑦=0). The determinant is defined as Δ=𝜆²𝜖²₁₁+𝜖²₁₂−𝜔²₁₂. (7) If Δ>0, the fiber approaches an asymptotic orientation with 𝜙=−𝜆𝜖₁₁±2Δ𝜆𝜖₁₁−𝜔₁₂,𝜋𝜃=2if𝜃₀=𝜋2,0otherwise, (8a)(8b)where 𝜃₀ is the initial polar angle. The results state that except, for very special initial condition 𝜃₀=𝜋/2, a fiber will approach the shear plane (𝜃=𝜋/2) and the azimuthal angle given by (8a). If Δ<0, the fiber rotates periodically, with period 𝜋𝑇=|Δ|. (9) It is worth pointing out that the results on the asymptotic angle (if Δ>0) and the period (if Δ<0) differ from those of Zhou et al. [30] due to the difference in the definition of Δ (a constant factor 4). The results here are slightly more general in applicable conditions and more compact in form.

In this simulation, a large number of Lagrangian trajectories (over one million) are tracked. The evolution of fiber concentration and orientation along these trajectories is simulated. The trajectories reveal how fluid parcels move in the mixing layer. Figure 5 shows the Lagrangian trajectories that pass through the same spatial cell (at various time instants). The cell is located at 𝑥^∗=319 and corresponds to the median of the mean concentration, 𝑐=0.5. The median 𝑐=0.5 does not lie in the center line of the mixing layer but shifts to the slow stream side (here negative 𝑦). From a Lagrangian point of view, this can be explained by the ensemble average of the concentration along the trajectories passing through the cell. Since there is no diffusion, the concentration along a Lagrangian trajectory is constant (see (5)); the constant is determined by the initial value at the inlet of the mixing layer. From Figure 5 it is obvious that a fluid parcel originating from the fast stream (𝑦>0) is more likely to travel down to the slow stream region than a fluid parcel originating from the slow stream (𝑦<0) to travel up to the fast stream region. That is why the median 𝑐=0.5 is located in the slow stream region. On another aspect, the fluctuation of the trajectories also quantitatively reflects the strength of the mixing.

Figure 5 

Lagrangian trajectories passing through a point at various time instants.

The fiber orientation along a Lagrangian trajectory is determined by the vorticity and deformation rate that a fiber undergoes. Figure 6 shows the profiles of the vorticity and deformation rate along a sample Lagrangian trajectory. The profile of the determinant Δ from (7) is also given in the figure. Theoretical analysis demonstrates that a fiber will rotate periodically if Δ<0. From (7) it is clear that only 𝜔₁₂ contributes to the negative part of the determinant. In other words, the antisymmetric part of the strain rate (vorticity) drives a fiber to rotate, and the symmetric part (deformation rate) tends to push a fiber to a specific orientation. In fact, as Lipscomb et al. [31] pointed out, (2) may be interpreted physically as stating the 𝐩 rotates with the fluid according to term “𝝎⋅𝐩”, and simultaneously partially stains with the fluid according to term “𝜆𝝐⋅𝐩”. The term, −𝜆𝝐∶𝐩, compensates a change of length which resulted from the motion described by 𝝎⋅𝐩 and 𝜆𝝐⋅𝐩, because 𝐩 is of unit length. This physical interpretation is very informative. However, it is far from being straightforward to directly derive the conclusion from (2) per se. Here, the determinant Δ (see (7)) has very simple form, and it is straightforward to draw the conclusion from the form of Δ. By comparing the components of the vorticity and deformation rate in Figure 6, it is clear that the profile of 𝜔₁₂ is more irregular than those of 𝜖₁₁ and 𝜖₁₂ (notice the zig-zag structures). Along this sample trajectory, most of the Δ profile in Figure 6 is negative, indicating that a fiber is of periodic rotation most of the time along the sample Lagrangian trajectory but with changing period.

Figure 6 

Profiles of the vorticity and deformation rate along a sample Lagrangian trajectory. The continuous gray is the corresponding determinant Δ evaluated from (7).

Figure 7 shows the three components of the orientation vector for a fiber traveling along the same Lagrangian trajectory as discussed in Figure 6. Both 𝑝₁ and 𝑝₂ exhibit nearly periodic behavior. However, 𝑝₃ approaches 0. The initial values of 𝑝₁, 𝑝₂, and 𝑝₃ are not important; all will exhibit similar evolution behavior as shown in Figure 7; except that for a special initial value 𝑝₃=1 (and hence 𝑝₁=𝑝₂=0), the orientation vector will not change. The fact that 𝑝₃ approaches 0 shows that a fiber is very likely to align to the shear plane (𝑥-𝑦), insensitive to its initial orientation.

Figure 7 

Components of the fiber orientation vector along the same sample Lagrangian trajectory as used in Figure 6.

Figures 6 and 7 show the strain rate and orientation along the same sample Lagrangian trajectory. Along this selected trajectory, the magnitude of vorticity 𝜔₁₂ is relatively high compared to 𝜖₁₁ and 𝜖₁₂, most of the Δ value is negative, and hence a fiber mostly rotates periodically, with changing period predicted by (9). Investigation in a large number of Lagrangian trajectories shows that Δ is usually larger than zero, which means that fibers along these trajectories are most likely to approach an asymptotic orientation predicted by (8a) and (8b). However, since the strain rate is not constant along a trajectory, the asymptotic orientation also changes correspondingly. Hence, a fiber seems to swing incessantly when it tries to approach an evolving asymptotic orientation. Through averaging the orientation over a very large number of Lagrangian trajectories that pass through the same location (at various time instants), the orientation distribution at the location can be obtained. Figure 8 gives the orientation distribution probability density function (pdf(𝜙) and pdf(𝜃)) at various locations. It is worth pointing out that pdf(𝜙) is only shown in the range [−0.5𝜋,0.5𝜋], which has period 𝜋, since a fiber which orients at 𝜙 or 𝜙+𝜋 is indistinguishable in its distribution, and pdf(𝜃) is given in the range [0,0.5𝜋], which has period 0.5𝜋, since the flow field is symmetric with respect to the plane 𝜃=0. Distributions 1, 2, and 3 are at cross-section 𝑥^∗=228 and correspond to the mean concentration 𝑐=0.3, 𝑐=0.5, and 𝑐=0.7, respectively. Distributions 4, 5, and 6 are at cross section 𝑥^∗=319 and also correspond to 𝑐=0.3, 𝑐=0.5, and 𝑐=0.7, respectively. Through comparing the corresponding distributions at 𝑥^∗=228 and 𝑥^∗=319, no evident difference is found. It implies that the fiber orientation distribution (which is uniformly distributed at the inlet) has already reached a stable state, and it does not evolve along the streamwise direction any longer in the turbulent region. In Figure 8(a) the distributions pdf(𝜙) all have similar shape, a peak, and a valley. Under a simple shear flow, pdf(𝜙) will have a peak at 𝜙=0 [30, 32]. Here, due to the synthetic effect of vorticity and deformation rate the peak is in between 0 and 0.5𝜋, which is believed to be closely related to the steady asymptotic direction predicted by (8a) and (8b). On the other hand, the valley is believed to be related to the unsteady asymptotic direction in (8a) and (8b). For pdf(𝜃), it has a very high peak at 𝜃=0.5𝜋, which agrees with the analysis on the evolution of 𝑝₃ in Figure 7 that a fiber is very likely to turn to the shear plane (𝜃=0) and stays there. In Figure 8 the distributions 1 and 4 correspond to 𝑐=0.3. Both pdf(𝜙) and pdf(𝜃) are found to be flatter than the distributions corresponding to 𝑐=0.5 or 𝑐=0.7. This may be due to the position corresponding to 𝑐=0.3 is located in the slow stream of the mixing layer, and the magnitude of the strain rate is relatively smaller than those in the fast stream.

(a)

(b)

Figure 8 

Orientation distribution at six positions: 1, 2, and 3 are at 𝑥∗=228 and correspond to 𝑐=0.3, 𝑐=0.5, and 𝑐=0.7, respectively; 4, 5, and 6 are at 𝑥∗=319 and correspond to 𝑐=0.3, 𝑐=0.5, and 𝑐=0.7, respectively. (a) Azimuthal angle 𝜙; (b) polar angle 𝜃.

4. Conclusion

The concentration and orientation of suspended fibers in a mixing layer are investigated numerically. The flow is assumed to be homogeneous in the spanwise direction, and the effects of fiber additive on the flow are neglected (very dilute suspension). A Lagrangian particles scheme is used to deal with the convection of fibers. With this Lagrangian particles scheme, fiber concentration and orientation evolve along Lagrangian trajectories independently. Ensemble average over a large number of Lagrangian trajectories is used to obtain statistically steady values of concentration and orientation. This Lagrangian particles scheme is found to be very efficient to compute the fiber orientation, which is discretized by hundreds of points on the unit sphere to represent the fiber orientation distribution.

Two cases Sc=1 (diffusive) and Sc=∞ (nondiffusive) are investigated for the fiber concentration distribution. The fine structures of the instantaneous distributions under these two cases are very different due to the molecular diffusion. Sharp front of the concentration is observed in the nondiffusive case. However, there is no obvious difference in the mean concentration between the two cases. Molecular diffusion is negligible in determining the mean concentration of fibers in turbulent flows.

For the rotational dynamics of a fiber, the analytical solution of Zhou et al. [30] is slightly generalized, and the new solution is presented in a more compact and informative form. A fiber will rotate periodically if the determinant is negative, where the determinant is defined with the strain rate. A fiber will approach an asymptotic orientation if the determinant is positive. The sign of the determinant is determined by the relative magnitude of the deformation rate and the vorticity. The symmetric part of the strain rate tends to make a fiber align to an asymptotic orientation, while the antisymmetric part makes a fiber rotate. This general conclusion helps understand the evolution of the fiber orientation along Lagrangian trajectories. When a fluid parcel passes through a region with relatively high shear rate, fibers followed by the fluid parcel are most likely to rotate incessantly. On the other hand, in the region of relatively high extension rate, fibers tend to align to some asymptotic orientation. Analysis on the orientation distribution shows that the distribution is stable in the turbulent region, which does not change along the streamwise direction. It is also believed that the orientation is not sensitive to the initial distribution at the inlet of the mixing layer. The pdf of the azimuthal angle 𝜙 has the shape of a peak and a valley, which most likely correspond to the steady and unsteady asymptotic orientations predicted by the analytical solution. The pdf of the polar angle 𝜃 has a prominent peak at 𝜃=0.5𝜋, which shows that fibers are very likely to align with the shear plane. This fact has significant implications in predicting the rheology of fiber suspension flows.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

Kun Zhou would like to thank the help from A. Attili and F. Bisetti. Wei Yang is supported by the National Natural Science Foundation of China (Grant no. 11302110). Zhu He is supported by the National Key Technology R&D Program of Chian (Grant no. 2011BAK06B02).