1 Introduction

Fiber suspension flows are found in many engineering applications such as fluidized bed, mixing, and paper-making. Among the fiber suspension flows, the turbulent jet flow is a kind of special flows. It represents a class of important flows, and provides a building block for many practical applications, e.g., combustion, mixing, and exhausting devices.

The particle suspension flows in the turbulent round jet have been widely studied in recent twenty years. Fleckhaus et al.^[1] found that turbulence attenuation increased with the particle size. Mastorakos et al.^[2] showed that the flow axial turbulence was produced by the interaction between the radial turbulence fluctuations and the cross-stream spatial gradients in the mean velocity of the particles. Sargianos et al.^[3] found that the fan-spreading turbulence was greatly affected by the mean and standard deviations of the particle distribution. Ramanujachari and Natarajan^[4] showed that the turbulent kinetic energy, Reynolds stress, jet spread, and entrainment rates decreased with the increases in the particle mass load and size. Prevost et al.^[5] pointed out that the production by the particle velocity gradient and transport terms affected the streamwise particle velocity, which was responsible for the increase in the streamwise particle velocity with the increase in the particle relaxation time. Kennedy and Moody^[6] indicated that the particle dispersion increased with the time. Lin and Zhou^[7] presented the interface instability curves for a moving jet at different particle parameters. Young et al.^[8] showed that the particle removal efficiency was directly but not linearly related to the exerted shear stress.

Actually, the studies shown above are all for spherical particles, whereas many industrial particles are non-spherical. Fiber suspension is a typical example. Compared with spherical particles, the fiber motion is more complicated due to its shape anisotropy and strong coupling effects of rotation and translation. Spatial and orientational distributions will affect the suspension and final properties of a homogeneous fiber material, e.g., thermal conductivity, absorption, and reflection of light. Therefore, it is necessary to understand the property of the fiber suspension in turbulent round jets in the optimization of the production processes and the design of the novel technology for fiber manipulation.

There are limited studies on the properties of the fiber suspension in a turbulent round jet. Filipsson et al.^[9] assumed that the Kelvin-Helmholtz (K-H) instability mechanism was important due to the presence of an inflexion point in the velocity profile, and suggested that the suppression of small-scale turbulence by fibers may be attributed to a damping of high frequency disturbances in a flow susceptible to the K-H instability. Lin et al.^[10] showed that the effect of turbulent fluctuation increased when the distance from the jet exit increased. Capone et al.^[11] studied experimentally a water turbulent jet laden with nylon fibers at two different mass concentrations, and found that fibers had an effect on the flow turbulence intensity in a distinct region. The effect was enhanced when the fiber concentration increased, whereas the results were proven to be Reynolds number independent. Capone et al.^[12] studied experimentally the near field of a turbulent round jet laden with the rod-like particles with the particle image velocimetry (PIV) at two mass fraction loads and the Reynolds number of 9 000, and found that, for the Stokes number around unity, the inertial effect gave rise to the velocity lag among the particles and fluid.

From the above literature survey, we can see that previous theoretical and numerical studies have never addressed the Brownian motion of fibers in the turbulent round jets, whereas the fibers with or without Brownian motion will show different dynamic characteristics. Therefore, the present paper aims to explore the effects of the Reynolds number, fiber volume fraction, and aspect ratio on the turbulent properties of the fiber suspension in a turbulent round jet with the consideration of the Brownian motion of fibers.

2 Flow and mathematical model 2.1 Governing equation of turbulent fiber suspension flow

The fiber suspension in a turbulent round jet flow is shown in Fig. 1, where the fiber orientation is also given. For the incompressible turbulent flow, the instantaneous velocity, pressure, rate-of-strain tensor, and orientation tensor can be written as a sum of the mean and fluctuation quantities. We substitute the sum into the instantaneous continuity equation and the modified Navier-Stokes equation^[13], and take the average. Then, we have

(1)

(2)

The last term on the right-hand side of Eq. (2) is the contribution of the fibers to the fiber suspension through exerting an additional stress on the suspension. In the above equations, U_i and P are the mean velocity and pressure of the fiber suspension, respectively. ρ is the density of the fiber suspension. µ is the fluid viscosity. is the Reynolds stress. ε_kl is the mean rate-of-strain tensor. I_ij is a unit tensor. a_kl and a_ijkl are the mean second-and fourth-order tensors of the fiber orientation, respectively. µ_a, as a function of the fiber concentration, aspect ratio, and orientation distribution, can be given by^[14]

(3)

where Φ is the fiber volume fraction, and λ is the fiber aspect ratio.

The density of the fiber suspension in Eq. (2) is

(4)

where subscripts fl and fi are used to represent fluid and fiber, respectively.

The mean second-and fourth-order tensors of the fiber orientation in Eq. (2) can be expressed as follows^[15]:

(5)

where p_i is a unit vector parallel to the fiber principal axis, ψ(p) is the mean probability density function for the fiber orientation, and p is the orientation vector.

2.2 Probability density functions for fiber orientation

As shown in Eq. (5), ψ(p) should be derived in advance in order to obtain the mean second-order tensor a_ij and the fourth-order tensor a_ijkl of the fiber orientation. There are two kinds of factors affecting the fiber orientation. One is the orienting factors embodied in the fiber rotational velocity, such as the fluid dynamic gradient and the external torque and force, and the other is the randomizing factors embodied within the fiber rotational diffusion tensor, i.e., the rotational Brownian motion and turbulent fluctuation.

The equation of probability density function ψ(p) for the fiber orientation is

(6)

where is the gradient operator projected onto the surface of the unit sphere, D_rB is the Brownian rotary diffusion coefficient as follows^[16]:

(7)

where T is the temperature, k_B is the Boltzmann constant, L_fi is the fiber length, and δ_rL and δ_rS are^[16-17]

(8)

In Eq. (6), is the particle angular velocity defined by^[18]

(9)

where ω_ij is the vorticity tensor defined by

η=(λ²−1)/(λ²+1) is related to the fiber aspect ratio λ, the last term on the right-hand side represents the interaction between the fibers, and D_rI is the rotary diffusion coefficient^[19] defined by , which provides that D_rI is isotropic^[20]. D_rI is equal to zero if the fiber is dilute.

We express the instantaneous probability density function ψ(p) and the fiber angular velocity as a sum of the mean fluctuation quantities, substitute these variables into Eqs. (6) and (7), and average the equation. Then, we have

(10)

where α_ψx and α_ψp are the dispersion coefficients of the linear and angular displacements, respectively, which are related to the fluid kinetic viscosity ν, the turbulent kinetic energy k, and the turbulent dissipation rate ε^[21], i.e.,

2.3 Turbulence kinetic energy and dissipation rate

In Eq. (2), the Reynolds stress tensor can be given as follows:

(11)

where µ_T is the eddy viscosity defined by

and it is assumed to be isotropic.

In order to solve Eqs. (1) and (11), the k-equation and ε-equation should be given. The fibers inevitably affect the turbulent kinetic energy k and the dissipation rate ε. Therefore, the k-equation and ε-equation should be modified with the consideration of the presence of fibers, i.e.,

(12)

(13)

where C₁ = 1.44, C₂ = 1.92, σ_k = 1.0, and σ_ε = 1.3. S_k and S_ε are the source terms resulting from the fibers. Based on the momentum equation, Lin and Shen^[22] derived the expressions of S_k and S_ε through a complex process as follows:

(14)

(15)

2.4 Equation for spatial distribution of fibers

Equations (3) and (4) include the particle volume fraction Φ, which usually has a non-uniform distribution in the flow because of the flow convection and diffusion. Therefore, it is necessary to build up and solve the equation for the spatial distribution of the fibers to get Φ for laying the foundation to solve Eqs. (3) and (4).

We express the instantaneous flow velocity u_j and fiber number density n(v) as a sum of the mean fluctuation quantities, substitute them into the instantaneous equation of the fiber number density, and average the equation. Then, we have

(16)

where n(v) is the mean fiber number density, and v is the fiber volume. On the left-hand side of Eq. (16), the second term is the convection term, the third term is the Brownian diffusion term, and the last term is the change in n(v) which is usually given by^[22]

where ν_t = 0.09k²/ε is the eddy viscosity. D_tB is the Brownian translational diffusion coefficient shown as follows^[19]:

(17)

where L_fi is the fiber length, and δ_t|| and δ_t⊥ are^[19-20]

(18)

Substituting Eq. (17) and into Eq. (15), we get the mean equation of the fiber number density as follows:

(19)

Multiplying the mean fiber number density by v^k, integrating over the entire volume distribution, and taking k=0, 1, we have^[23]

(20)

where M₀ = N is the total fiber number, M₁ = V is the fiber volume and directly proportional to the fiber mass at a definite fiber density. Based on the particle volume V, the fiber volume fraction Φ can be calculated.

2.5 Components of the fourth-order orientation tensor

According to the definition of θ and ϕ shown in Fig. 1, the unit orientation vector p can be presented as follows:

(21)

The components of the fourth-order orientation tensor in Eq. (21) have the following forms based on Eq. (5):

(22)

(23)

(24)

3 Numerical method and parameters

The main steps of the numerical simulation are as follows:

(ⅰ) Solve Eqs. (1)-(3) and (11)-(13) with Φ = µ_a = S_k = S_ε = 0 (i.e., no fiber) to get U_j, P, k, ε, and .

(ⅱ) Solve Eqs. (17)-(20) to obtain n(v) and Φ.

(ⅲ) Substitute Φ into Eq. (3) to get µ_a.

(ⅳ) Substitute U_j, k, ε and Eqs. (7)-(8) into Eq. (10), and solve the obtained equation to get ψ.

(ⅴ) Substitute ψ into Eqs. (5) and (22)-(24) to get a_ij and a_ijkl.

(ⅵ) Substitute µ_a, a_ij, and a_ijkl into Eqs. (1)-(2) and (11)-(15) to get U_j, P, k, ε, and .

(ⅶ) Turn to Step (ⅱ) based on the new values of U_j, P, k, ε, and if necessary.

Equations (2), (11)-(13), and (19) are solved by the finite difference approach. The central finite difference and the second-order upwind finite difference are used to the diffusion term and convective term, respectively. Equation (5) is integrated with the Simpson formula. The grid system consists of 120(x)×60(y)×60(z) = 432 000 grid points. The mesh size is uniform in the streamwise (x) direction, while it is refined in the y-and z-directions. The computational grid independence and the suitability of the grid size for the convergence results have been tested.

The fiber suspension is a mixture composed of water and nylon fiber. The water density and viscosity are ρ_fl = 998.2kg/m³ and ν = 1.007×10^-6 m²/s, respectively. For the fiber, the density is ρ_fi = 1140kg/m³. The values of the volume fraction Φ are 0.002%, 0.01%, and 0.1%. The values of the fiber aspect ratio λ are 3, 5, 15, and 30. The flow Reynolds numbers based on Re = DU₀/ν, where U₀ is the velocity, and is 5 000, 9 000, or 16 000 at the nozzle exit. The Boltzmann constant k_B is 1.38 × 10^-23 J/K.

In the present study, some results are compared with the experimental results^[12]. The related parameters in the experiment are as follows: the fiber density is 1 140 kg/m³, the fiber aspect ratio is 13.3, the Reynolds number is 9 000, and the fiber volume fractions are 0.001 7% and 0.005 2%, respectively.

4 Results and discussion 4.1 Visualization of numerical results at definite time

Figure 2 shows the jet flow and distribution of the fibers at a definite time. We can see that the fibers with different orientations are randomly distributed in the flow.

4.2 Mean streamwise velocity at the centerline along the jet axis

The turbulent properties of the fiber suspension depend on the fiber translational and rotational diffusion and the momentum transfer which are directly related to the flow shear rate. It is important to understand the effects of the fiber volume fraction, the fiber aspect ratio, and the Reynolds number on the mean streamwise velocity distribution because the flow shear rate is determined by the mean velocity distribution.

The numerical mean streamwise velocities at the centerline along the jet axis for the single-phase and fiber-laden cases are shown in Figs. 3-5, where the experimental results^[12] are also given. From the figures, we can see the effects of the fiber volume fraction, Reynolds number, and fiber aspect ratio on the mean streamwise velocity. The values of the velocity are nearly the same for different fiber volume fractions, Reynolds numbers, and fiber aspect ratios in the region close to the nozzle exit (x/D < 1.5), i.e., the fibers have a small effect on the velocity because the fibers do not have enough time to affect the flow. When the distance from the nozzle exit increases, a larger velocity for the fiber suspensions in comparison with the single phase can be observed, which becomes more obvious when the fiber volume fraction and aspect ratio increase and the Reynolds number decreases. Therefore, the fiber injection in the flow has an effect of delaying the mean streamwise velocity decay along the jet axis. This can be attributed to the shift in the orientation undergone by the fibers due to their interaction with the flow.

4.3 Radial distribution of the mean streamwise velocity

The radial distributions of the mean streamwise velocity for the single-phase and fiber-laden cases with different volume fractions at x/D=5 are shown in Fig. 6. We can see that the maximum velocity appears at the centerline, and the velocity decreases from the center to the outside. The flow with fibers shows an increase in the mean streamwise velocity. The effect of fibers on the velocity profile is obvious, especially in the region near the center, and the degree of effect is proportional to the fiber volume fraction. The mean velocity gradient is directly related to the shear stress which determines the momentum transfer between the vortices with different scales. The fiber alignment in the fiber suspension will result in a large extensional viscosity in the inter-vortex extensional flow and annihilate the small vortical structure, which makes a reduction in the momentum transfer through suppressing the turbulence and causing the momentum transfer to be less than that of single-phase. This may be the reason why there is an increase in the velocity for the flow with fibers.

Figure 7 shows that the mean streamwise velocity increases when the Reynolds number decreases. This is different from the case of varying the fiber volume fraction. When the Reynolds number increases, the phenomenon of fiber alignment becomes weaker, and the resulting extensional viscosity becomes smaller, which strengthens the small vortical structure and increases the momentum transfer through enhancing the turbulence.

The radial profiles of the mean streamwise velocity for different fiber aspect ratios are shown in Fig. 8. The variation tendency of velocity with the fiber aspect ratio is the same as that with the fiber volume fraction. The mean velocity increases when the fiber aspect ratio increases. When the fiber aspect ratio increases, the phenomenon of fiber alignment becomes more obvious and the resulting extensional viscosity becomes larger, which annihilates the small vortical structure and makes a reduction in the momentum transfer through suppressing the turbulence.

4.4 Radial distribution of the turbulent kinetic energy

Figure 9 shows the radial distributions of the root-mean-square (RMS) of the streamwise velocity fluctuation which is the main part of the turbulent kinetic energy. The RMS of the streamwise velocity fluctuation is the highest around the position of y/D=3, where the velocity gradient is also the largest. The RMS of the streamwise velocity fluctuation in the presence of fibers increases at all positions along the radial direction, and the increasing magnitude is proportional to the fiber volume fraction. The maximum increasing rate in the RMS is at the centerline, and the increasing rate decreases along the radial direction.

The distributions of the turbulent kinetic energy along the radial direction for different fiber volume fractions, Reynolds numbers, and fiber aspect ratios are shown in Figs. 10-12, where y_1/2 is the half jet width. A larger turbulent kinetic energy for the fiber suspensions in comparison with the single phase can be seen. This is because that the fibers increase the momentum transfer among the fibers and fluid by providing a solid link between the adjacent fluid layers. For spherical particles, there exists a critical scale, i.e., the ratio of the particle diameter to the turbulent length scale, which can be used to determine whether the turbulent kinetic energy caused by the addition of particles increases or decreases. Usually, large particles can enhance the turbulence, while small particles are just the opposite. For fibers, there is no characteristic length such as the diameter of the spherical particle, the complex dynamics of fibers stemming from their shape anisotropy and the strong coupling effect of rotation and translation induce the velocity fluctuation of the suspension, which leads to a similar outcome on the turbulent kinetic energy when the spherical spheres and the turbulence due to the wakes increase.

As shown in the figures, the enhancement degrees of the turbulent kinetic energy increase when the fiber volume fraction, Reynolds number, and fiber aspect ratio increase. Higher fiber volume fraction brings about a stronger effect on the flow. Larger Reynolds number is usually responsible for higher turbulent kinetic energy. It is easy for the fibers with a large aspect ratio to increase the momentum transfer through providing a solid link between the adjacent fluid layers, and longer fibers require smaller volume fraction to generate the same stress distribution as shorter fibers do.

4.5 Radial distribution of the Reynolds shear stress

The variations of the Reynolds stress along the radial direction are shown in Figs. 13 and 14 for different fiber volume fractions, Reynolds numbers, and fiber aspect ratios. The maximum value of the Reynolds stress occurs around the position of y/D=3, where the velocity gradient is the largest because the production of the Reynolds stress is directly related to the velocity gradient. The Reynolds stress is enhanced with the increase in the fiber volume fraction, Reynolds number, and fiber aspect ratio. The reason is analyzed in Subsection 4.4. The effect of fibers on the Reynolds stress is relatively weak in the regions near the center and outside of the jet.

5 Conclusions

The equations including the turbulent fiber suspension flow, probability density function for fiber orientation, and the spatial distribution of fibers are numerically simulated by the finite difference approach. Some numerical results are compared with the corresponding experimental data. The following conclusions can be drawn.

Fibers have a small effect on the mean streamwise velocity in the region close to the nozzle exit. When the distance from the nozzle exit increases, the fiber injection in the flow has a delay on the mean streamwise velocity decay along the jet axis, and such an effect becomes more obvious when the fiber volume fraction and aspect ratio increase and the Reynolds number decreases. The flow with fibers shows an increase in the mean streamwise velocity along the radial direction, and the increasing magnitude is directly proportional to the fiber volume fraction and fiber aspect ratio, and is inversely proportional to the Reynolds number.

The maximum values of the RMS of the streamwise velocity fluctuation and the Reynolds stress occur around the position of y/D=3. The flow with fibers shows an increase in the RMS, and the increasing rate is the largest at the centerline and is proportional to the fiber volume fraction. The presence of the fibers leads to an enhancement in the turbulent kinetic energy and Reynolds stress, and the enhancement degree increases when the fiber volume fraction, Reynolds number, and fiber aspect ratio increase. A higher fiber volume fraction brings about a stronger effect on the flow. The fibers with large aspect ratio are easy to increase the momentum transfer through providing a solid link between the adjacent fluid layers, and longer fibers require smaller volume fractions to generate the same stress distributions as shorter fibers do.