Fully-resolved prolate spheroids in turbulent channel flows: A lattice Boltzmann study

Particulate suspensions are common in a wide range of industrial applications and environmental processes. Examples include paper industry, icy clouds, biomass combustion, crystallization and pharmaceutical processes.^1–3 The importance and vast application of particle-laden flows necessitate proper simulation and prediction of particle-fluid interactions.

In most practical cases, the flow regime is turbulent, leading to a chaotic behavior encompassing a wide range of time and length scales. The effect of finite-size particles on turbulence has been the subject of many studies.^4–7 However, these studies are mostly concerned with spherical particles. For example, Pan and Banerjee⁸ used a pseudo-spectral method to model particle-laden turbulent flows in open channels with particle volume fractions of $ϕ∼10−4$⁠. An increase of turbulence intensities and Reynolds stress was reported. Kajishima et al.⁹ simulated an upward turbulent flow in a vertical channel with solid spherical particles. The particle diameter and concentration were d_p/H = 0.125, (with H being the half channel height) and $ϕ∼10−3$⁠, respectively. Reduction of mean velocity and enhancement of velocity fluctuations was observed. Uhlmann¹⁰ simulated the motion of fully-resolved spherical particles in a vertical channel flow (⁠$ϕ=4.2×10−3$⁠) with direct numerical simulation (DNS). Enhancement of streamwise velocity fluctuations and turbulence intensity with respect to the single-phase flow at the same bulk Reynolds number was reported. Picano, Breugem, and Brandt¹¹ used immersed boundary method (IBM) to compute the back-influence of finite-size spheres on turbulent flow. For particle size of d_p/H = 1/9, increase of Reynolds shear stress at $ϕ=0.05,0.1$⁠, but its reduction at $ϕ=0.2$ was observed.

Though it is easier to treat numerically particle motion for spherical particles, most particles do not show a perfectly spherical shape in practical cases. In this work, we study the effect of prolate spheroids on the turbulent flow field. A spheroid is an ellipsoid with two diameters being equal, and a prolate spheroid is obtained by rotating an ellipse around its major axis. For such a configuration, most publications deal with laminar flows and consider a single suspended particle.^12–15 

The motion of non-spherical particles is inherently complex, even in laminar flows. In turbulent regimes, the problem is even more intriguing. During the last years, several researchers have modeled turbulent flows with non-spherical particles.^16–20 However, in most of these studies, the particle size is smaller than the Kolmogorov length scale (leading to so-called point particles). Moreover, only one-way or two-way coupled simulations were performed, which only hold for low volume fractions of the solid phase. Thus, particle-particle interactions are mainly ignored in these studies. Zhang et al.²¹ used DNS to compute the transport and deposition of ellipsoidal particles in dilute turbulent channel flows. The simulation was carried out by one-way coupling assumption. The accumulation of particles in the viscous sublayer and their alignment with the flow direction was reported in this work. Lin et al.²² studied a turbulent channel flow with cylindrical fibers. They showed that the flow rate of the fiber suspension is larger than that of a Newtonian fluid. Turbulent intensity and Reynolds stress were found smaller in the fiber suspension. Mortensen et al.²³ studied the distribution and orientation of small spheroids by means of DNS. The simulation was performed in an Eulerian-Lagrangian framework. Forces and torques were computed under creeping flow condition (Stokes flow), which is only valid for low particle Reynolds number. A preferential orientation of elongated particles in the streamwise direction was reported. Later, Marchioli, Fantoni, and Soldati²⁴ studied a turbulent channel flow with long particles and concluded that the aspect ratio has a minor effect on clustering, preferential distribution, and segregation of particles. Andersson, Zhao, and Barri²⁵ showed drag reduction effect of point particles with an aspect ratio of 5. Zhao and van Wachem²⁶ investigated the behavior of elongated particles in a turbulent channel with frictional Reynolds number of $Reτ=150$⁠. Particle-particle and particle-wall collisions were modeled in this work. However, the flow was in the dilute regime and the particles were smaller than the Kolmogorov length scale. Drag enhancement by particles was reported here. These authors later showed that, in the dilute regime, point particles of high aspect ratio and moderate Stokes number have minor drag reduction effect.²⁷ 

Investigations of finite-size anisotropic particles in turbulent flow are scarce. The effect of finite-size cylindrical fibers on turbulent flow was studied by Do-Quang et al.²⁸ They employed an external boundary force with a lattice Boltzmann method (LBM) to model the effect of rigid rods with a cylindrical shape for volume fractions up to 0.004. Very recently, Ardekani et al.²⁹ performed DNS of suspensions of oblate spheroids using IBM.

DNS of fully-resolved prolate spheroids in turbulent wall-bounded flows has not been performed yet to our knowledge. This is the purpose of the present work. All particles considered are finite-size and their volume fraction corresponds to the dense regime (⁠$ϕ>10−3$⁠). Moreover, particle-particle and particle-fluid interactions will be taken into account (four-way coupling). The outline of this paper is as follows: Section II describes the governing equations of LBM, IBM, and particle dynamics. The simulation set-up, boundary, and initial conditions are discussed in Sec. III. A validation is presented in Sec. IV. Section V describes the results of turbulent channel flows seeded with prolate spheroids, and the comparison with single-phase flow and spherical particles. Finally, conclusions are given in Sec. VI.

Lattice Boltzmann method is a promising alternative to classical CFD approaches that are based on Navier-Stokes equations. LBM is a mesoscopic approach in which probability distribution functions represent the fluid behavior. Three considerable advantages of LBM, 1) locality of all calculations (leading to excellent parallelization), 2) pressure deduced from equation of state (no Poisson equation), and 3) a linear convective term, have made it quite popular in recent years. These features contribute to simpler implementation and higher parallel scalability of the algorithm, together with a straightforward implementation of complex geometrical boundaries. It can be shown that a single relaxation time (SRT)-LBM formally recovers the Navier-Stokes equations in the low Mach number limit. This form of LBM is written as

where f_i is the distribution function of particles moving with speed $ci$ between lattice nodes and the right-hand side accounts for the Bhatnagar, Gross and Krook (BGK) collision term.³⁰ Here, $τ$ is the dimensionless relaxation time and $Δt$ is the time step. The equilibrium distribution function $fieq$ reads:

where $u$ is the macroscopic velocity vector, $ρf$ is the fluid density, $ωi$ is the weight associated with the velocity $ci$⁠, and the sound speed $cs=1/3$⁠. The macroscopic velocity $u$ in Eq. (2) must satisfy the requirement for low Mach number regime, i.e. $u/cs=Ma≪1$⁠. In the present work, the nineteen-velocity model (D3Q19) with discrete velocity vectors $ci=[cix,ciy,ciz]$ is used. The relaxation time $τ$ is related to the fluid kinematic lattice viscosity $ν$ through:

Finally, pressure is calculated from density via an equation of state: $p=ρfcs2$⁠.

In order to consider body forces in LBM, the force scheme presented by Shan and Chen³¹ is adopted here, where the equilibrium velocity is shifted.

The immersed boundary method was originally developed to model blood flow in the heart.³² In this work, the direct-forcing IBM is used to couple the fluid and particle phases.³³ Particles are thus represented by Lagrangian points, while the flow field domain is described on a fixed Cartesian grid.

It can be shown that, the particle-fluid interaction force density is computed from Eq. (4) by subtracting the Navier-Stokes equations with and without a force term. Thus,

where $unoF$ is the interpolated fluid velocity without a force term at each Lagrangian point, and $Ud$ is the particle velocity at the same Lagrangian point at previous time step. Here, $unoF$ is obtained by interpolation using Eq. (5):

In this equation, $Δh$ is the lattice size (= 1) and $D(xi,j,k−Xl)$ is a 4-point discrete delta function.³⁴ Parameters $Xl$ and $xi,j,k$ denote the location of Lagrangian and Eulerian points, respectively. The term $Ud$ at each Lagrangian point reads

where $Up$ and $Ωp$ are the translational and angular velocities of the particle, respectively, and $Xc$ is the location of particle center. Finally, the force term that is exerted on each Eulerian point (fluid), is calculated by spreading the Lagrangian force term:

with $ΔVl$ being the area of the relevant Lagrangian boundary point segment. This force term is then applied to the relevant flow field equations.

A prolate spheroid surface (Fig. 1) is described by the following equation:

where (⁠$x′$⁠, $y′$⁠, $z′$⁠) denotes a body-fixed coordinate system and a and b are the radii of the particle. In this study, a prolate spheroid is considered, i.e., $a>b$⁠, with particle aspect ratio of $λ=a/b>1$⁠, a spherical particle corresponding to $λ=1$⁠.

The general equations of motion for particles reads:

Equations (9) and (10) describe the translational and rotational motion of the particle, respectively. Here, $T$ is the total torque exerted on particle center, M is the mass, and I_p is the particle mass moment of inertia. The second term on the right hand side of Eq. (9) accounts for the mass of the fluid phase that is inside the particle volume. Subscripts f and p represent the fluid and particle, respectively, and $Mf=ρfVp$⁠. Particles experience a repulsive force $FR$ when they are closer than a threshold to the wall or other particles. Translational motion of the particle is updated by using the discretized form of Eq. (9).

To calculate the repulsive force $FR$⁠, first the gap between particles should be estimated. To find the closest points between two spheroids, the iterative algorithm proposed by Lin and Han³⁵ and implemented by Ardekani et al.³⁶ is employed. In this approach, first, an imaginary sphere is constructed inside each spheroid with its center coinciding with the particle center. The radius of this sphere is assumed equal to the lowest radius of curvature of the spheroid. This ensures that the sphere fits in the spheroid volume. Then, the intersection point of the connecting line of the centers of imaginary spheres and each spheroid surface is found. In the next step, another imaginary sphere is constructed at each of the already found intersection points; this sphere is tangential to the spheroid surface. The intersection of the connecting line of two new spheres centers and each spheroid surface is found again. This procedure continues until the change in the slope of the line gets smaller than a threshold. This efficient algorithm finds the smallest gap between two particles by few iterations.

After finding the closest points of two spheroids, the repulsive force must be applied. To do so, a combination of lubrication and spring force models is used here. The lubrication force model for two spheres of radii R_a and R_b reads³⁷ 

Here, R_a and R_b are the radius of curvature of each spheroid at the closest points, while $ξ$ and s are the gap between two spheres and the threshold distance, respectively. The parameters $xa$ and $xb$ denote the location of the two spheres center. When $ξ<s$⁠, the lubrication force model is activated. When the gap between two particles gets smaller than the Eulerian grid size, the repulsive force model of Glowinski et al.³⁸ is used.

To describe now the rotation of a prolate spheroid, two coordinate systems are employed: the inertial frame and the body-fixed frame. The inertial frame, $x=⟨x,y,z⟩$⁠, is the frame that spans the computational domain (Eulerian framework). The particle body-fixed frame, $x′=⟨x′,y′,z′⟩$⁠, is attached to the particle with its origin at the particle mass center and its axes along the major and minor axes of the particle. Equation (10) is thus written in the body-fixed frame:

In these equations, all variables labeled with a prime are referenced with respect to the body-fixed coordinate. The principal moments of inertia are $Ixx′=2Mpb2/5$⁠, $Iyy′=Mp(a2+b2)/5$ and $Izz′=Mp(a2+b2)/5$ with M_p denoting the particle mass. To solve Eqs. (12)–(14), four quaternion variables are introduced:

where $(θ,ϕ,ψ)$ are the Euler angles (see Fig. 2). The transformation from the inertial coordinate to the body-fixed frame is done by: first, a rotation by an angle $ϕ$ about the z-axis; second, a rotation by an angle $θ$ about the new x-axis (shown by N), and third, a rotation by an angle $ψ$ about the $z′$-axis.

The transformation matrix M from the inertial to the body-fixed coordinate reads

and thus

These two equations are used to map hydrodynamic moment and angular velocity vectors from the inertial to the rotating (body-fixed) frame.

The evolution of the quaternion parameters is governed by

In our CFD code, torques are first transformed from the inertial to the body-fixed coordinate using Eq. (20). Then, Eqs. (12)–(14) together with Eq. (22) (7 equations in total) are solved using a 4^th-order Runge-Kutta integration algorithm to obtain each quaternion q_i and the angular velocity $Ω′$ in the body-fixed coordinate system. Next, the rotation velocity in the inertial coordinate $Ω$ is obtained by Eq. (21) using the inverse of matrix M. Finally, the location of each Lagrangian point can be updated by $X=M−1X′$⁠.

The simulation set-up includes a turbulent flow between two parallel plates located at a distance of 2H from each other. Domain size is L_x = 4.6H, L_y = 2H and L_z = 2H, corresponding to the streamwise x, wall-normal y, and spanwise z directions, respectively. The domain is uniformly discretized by $600×260×260$ grid nodes. Walls experience a no-slip boundary condition using the half-way bounce back approach. Periodic boundary conditions are imposed on streamwise and spanwise directions. A constant driving force is applied along the x direction, which is given by $G=ρfg=−dp/dx=ρfuτ2/H$⁠. The Reynolds number based on the friction velocity is defined by $Reτ=uτH/ν=180$⁠, where $uτ=(τw/ρf)0.5$⁠. Here, $τw$ is the shear stress at the wall. This leads to a bulk Reynolds number based on the bulk velocity, U_b, and full channel height, 2H, equal to $Reb=Ub2H/ν=5600$⁠.

In both particle-laden and particle-free cases, the data are collected when the turbulent flow is fully developed and statistically reaches the steady-state. The data are then gathered in steps of $0.04H/uτ$⁠. For the particle-laden flow, the particles are initially randomly positioned throughout the domain, but avoiding any overlapping. A channel with fully developed turbulent single-phase flow is used as initial condition for the flow field.

The results will be presented in wall (inner) units unless otherwise stated, and the velocity values are scaled by $uτ$⁠. The root mean square of the velocity fluctuations is defined by $urms=u′u′¯1/2$⁠, where $u′$ represents the velocity fluctuations. The other diagonal components, $v$_rms and $w$_rms are defined similarly. All fluid statistics are computed excluding the Eulerian nodes that are located inside the particles. All simulations are performed by our in-house LB solver ALBORZ,^39–41 which has been successfully validated for a variety of single-phase and multiphase flows, either laminar or turbulent.

In the case of particle-laden flows, three rigid particle types are tested: two spheroidal (P1, P2) and one spherical (P3). All particles are neutrally buoyant. Neglecting the gravity ensures that particle sedimentation does not occur, so that particles distribution are statistically identical on the upper and lower parts of the channel. Particles size, aspect ratio, and Stokes number are listed in Table I. Particle Stokes number is defined as the ratio of the particle response time $τp$ to the characteristic time scale of the turbulent flow. Characteristic flow time scale can be either the Kolmogorov (⁠$τK=ν/𝜖$⁠) or the near-wall scale (⁠$τν=ν/uτ2$⁠). In Table I, the latter approach is employed (i.e., $St+=τp/τν$⁠). Particle response time for spherical particles is given by

In the case of spheroidal particles, particle response time (using $λ=a/b$⁠) reads⁴² 

recalling that b is the length of particle semi-minor axis.

A single-phase turbulent channel flow at the same frictional Reynolds number of $Reτ=180$ has already been considered in detail for validation purposes as described in Eshghinejadfard et al.⁴¹ This will not be repeated here.

In order to check the behavior of the developed numerical tool for particulate flows, the movement of a single spheroid particle will now be investigated. The rotation of a spheroidal particle in a Couette flow at small Reynolds numbers, as obtained from the equations presented in Sec. II C, is compared with the available analytical solution for a single spheroid. Then, higher Reynolds numbers are tested. Both test cases consider a neutrally-buoyant spheroid suspended between two parallel plates that are moving in opposite directions.

For the pure rotation of the particle around the z-axis, the Reynolds number is defined as:

where G is the shear rate. Jeffery⁴³ showed that, the rotation angle $φ$ and the angular rate of rotation $φ̇$⁠, when one of the principal axes of the ellipsoid is kept parallel to the vorticity vector reads in a shear flow without inertia (Re_p = 0):

For the first simulation, the computational domain is discretized by $120×120×60$ lattice nodes. The Reynolds number is Re_p = 0.5 and the particle radii are a = 6 and b = 4.5 in lattice units. The particle is located at the domain center and between two parallel plates that are at y = 0 and y = N_y and move with velocity U in opposite directions. The resulting shear rate is G = 2U/N_y. Periodic boundary conditions are applied in x and z directions. There is no constraint on particle translation and rotation but the simulation proved that the particle only rotates around the z-axis without any translational motion. Results of the rotational speed of the particle are shown in Fig. 3. Our simulation results agree very closely with the analytical solution of Jeffery.⁴³ 

The accuracy of the prediction is then examined for higher Reynolds numbers. Two domain resolutions are used: (N_x, N_y, N_z) = (40, 40, 40) and (80, 80, 80). Particle semi-major and semi-minor lengths are 12 and 3 lattice cells, respectively. The particle Reynolds number varies between 30 and 100 and the shear rate is G = 1/600. Figure 4 illustrates $G⋅Tr$ values (T_r: period of a complete tumbling motion of the particle) versus Reynolds number. For both domain resolutions, our results are compared with the published data of Rosén, Lundell, and Aidun.⁴⁴ The results are quite close; minor differences can be attributed to the different numerical approaches employed, since there is no analytical solution in this case. Considering this very good agreement, the developed direct forcing IB-LBM approach can be now applied to a collection of spheroidal particles in turbulent channel flows.

Figure 5 shows a snapshot of P1 particles distribution together with the velocity field on three boundary planes for $ϕ=1.5%$⁠. Particles are distributed in the whole channel, impacting the turbulent velocity field. An accurate estimation of local particles concentration in each part of the domain requires time-averaging and will be presented later. The following sections describe velocity statistics, particles distribution, and their orientation.

To quantitatively investigate the effect of particle shape on turbulence, fluid mean velocity is considered along the channel height. Figure 6 compares the normalized mean streamwise fluid velocity (⁠$u+=u¯/uτ$⁠) for volume fraction of $ϕ=1.5%$⁠. It is seen that all particle types decrease the mean streamwise velocity at this volume fraction. In the region close to the wall all velocity profiles coincide due to using the same driving force. Velocity reduction by spheroids (P1, P2) is less than that for spherical particles (P3). P1 and P2 reduce the mean velocity by only 2% and 1.6%, respectively, compared to 4.6% for P3. Comparison of velocity profiles of P1 and P3 prove that particle Stokes number (almost identical for P1 and P3, see Table I) is not the only parameter that contributes to pressure drop. Particle shape is another important factor. It is also interesting to see that velocity reduction by P2 is slightly lower than for P1, although the former is smaller and has a lower Stokes number. For spheres, smaller particles are known to induce higher velocity reduction, as discussed for instance in Eshghinejadfard et al.,⁴¹ Shao, Wu, and Yu.⁷ 

Comparing the fluid velocity profiles for P1 at $ϕ=0.75%$ and 1.5% in Fig. 7 reveals that by increasing the volume fraction of the spheroids, mean fluid velocity is reduced in the region of $15<y+<120$ but is increased beyond this region toward the channel center (i.e., $y+>120$⁠); the green curve locates above the red one. For spherical particles (P3), reduction in both areas is seen from the graph.

Particle effect on rms of velocity fluctuations for three different orthogonal directions at $ϕ=1.5%$ is shown in Fig. 8, where $urms+=urms/uτ$⁠, $vrms+=vrms/uτ$ and $wrms+=wrms/uτ$⁠. Adding the particles reduces the maximum streamwise velocity fluctuation for all particle types. However, the turbulence attenuation effect is more obvious for spherical particles. Maximum rms of streamwise velocity fluctuation of P2 is slightly lower than P1 because of its lower Stokes number. With respect to wall-normal (⁠$vrms+$⁠) and spanwise (⁠$wrms+$⁠) directions, the profiles for spheroid particles (P1, P2) are closer to a single-phase flow along most of the channel height. The effect of particle shape on wall-normal and spanwise velocity fluctuations is stronger close to the wall. These observations confirm qualitatively findings concerning small cylinders in a turbulent flow.²⁸ Enhancement of velocity fluctuations in the vicinity of the wall is due to small-scale vortices generated in this region. Furthermore, for all directions, the peak point of rms velocity is shifted toward the wall by adding the particles. This arises from local concentration of particles near the wall, in particular for spherical particles, as discussed later.

It must be noted that, the effect of finite-size particles on velocity fluctuations can be different in dilute and dense regimes. Uhlmann¹⁰ reported enhancement of streamwise velocity fluctuations and reduction of normal and spanwise velocity fluctuations for finite-size spheres in the dilute regime. This phenomenon was attributed to vortex shedding behind the particles. In the dense regime, particles are able to change the structure of large vortices and generate small-scale vortices as well. This leads to a different behavior. Inertial point particles have been shown to augment streamwise intensity in numerical simulations of Zhao, George, and Van Wachem²⁷ and of Mortensen et al.²³ 

The profile of Reynolds shear stress is not significantly influenced by particles for this volume fraction and therefore is not shown here. It is just noted that spherical particles again show a larger difference from the single-phase flow concerning this quantity.

Turbulent kinetic energy (TKE) for single-phase and multiphase cases are depicted in Fig. 9. Based on the definition of TKE (⁠$k=12(ux′2+uy′2+uz′2)$⁠) and on Fig. 8, this quantity is mainly affected by streamwise velocity fluctuations. Therefore, its peak point is shifted downward for particle-laden cases. Again, the reduction for spheroidal particles is less than for spherical ones. Regions close to the wall and channel center remain almost unchanged.

Now, particles distribution, orientation and velocity statistics are considered. To do so, each Eulerian point that is located inside a particle, is treated with the relevant particle data (e.g., translational and rotational velocity).

Local volume fraction of particles for two overall volume fractions of $ϕ=0.75$ and 1.5% is shown in Fig. 10. It can be seen that spherical particles have a local peak of distribution close to the wall. The peak point corresponds to y⁺ = 26 for P3. For spheroids, there is no clear local peak. Moreover, for both volume fractions, the spherical particles P3 show a higher volume fraction near the wall.

At $ϕ=0.75%$⁠, particles P1 show a maximum volume fraction at $y/H≃0.2$ before decreasing toward channel center. However, at $ϕ=1.5%$⁠, this particle type shows a quite uniform volume fraction after reaching its maximum at $y/H≃0.2$⁠. This distance from the wall corresponds to about 0.8l_maj, with l_maj being particle major length (= 2a). A similar effect is observed for P2 particles, which show a uniform concentration beyond y/H = 0.2 at $ϕ=1.5%$⁠. This means that spheroidal particles apparently do not show preferential concentration near the wall, when their volume fraction is sufficiently high. Similar phenomenon was observed in the experimental study of fiber suspensions by Capone, Miozzi, and Romano.⁴⁵ On the other hand, spherical particles stay mainly trapped close to the wall. The main factor that helps spheroids leaving this region is their increased collision probability with the wall, which generates a torque on the particle leading to particle rotation and affecting near-wall vortices.

Normalizing each profile of Fig. 10 by the overall volume fraction, P1 at $ϕ=1.5%$ shows $ϕ(y)/ϕ≥1$ close to the channel center. For volume fraction of $ϕ=0.75%$⁠, values $ϕ(y)/ϕ<1$ are observed close to the channel center. Thus, more particles tend to move toward channel center when increasing volume fraction.

Figure 11 shows the mean velocity difference between fluid and particle at $ϕ=1.5%$⁠. In the viscous sublayer and part of the buffer layer (⁠$y+≲18$⁠), all particle types move faster than the fluid. In this area, fluid velocity is low. Since particles have finite size, part of the particles might be located in regions with a higher velocity. This results in higher average velocity of particles compared to local fluid layers. Moreover, P1 moves faster than P3 in this region due to its lower projected area. It is later shown that spheroids tend to align themselves in streamwise direction. Lower projected area combined with streamwise alignment lead to higher velocity of spheroids. It is also seen that P1 and P2 always move faster than the fluid before reaching the same speed in the center. In contrast, P3 tends to move slower than the fluid in $18<y+<33$⁠. This can be justified by the local maximum concentration of spheres around this limit, while spheroids have no preferential concentration at $ϕ=1.5%$ (see again Fig. 10). For point particles with low Stokes number, Mortensen et al.²³ showed that the particles velocity is usually less than the fluid velocity near the wall, while for point particles with high Stokes number, particle velocity exceeds that of the fluid.

Figure 12 displays the rms of particle velocity fluctuations. All particle types have generally lower velocity fluctuations than the relevant fluid one (F) due to their finite size, except for a small region very close to the wall. It is also observed in Figs. 12a, b and c that spherical particles reach a higher peak in streamwise fluctuations compared to spheroids. Normal and spanwise velocity fluctuations are less affected by particle shape. However, spheroids have higher wall-normal intensity near the wall compared to spheres. This is probably due to stronger effect of particle-wall collisions on spheroids.

Another important issue for spheroids is their orientation. Particle orientation has a major influence on particle-fluid interactions. This issue has been investigated in the literature for point-size prolate and oblate spheroids^{20,24,26,46,47} as well as cylindrical fibers.^17,48 It is generally found that, in one-way coupled simulations, the orientation of tracer and inertial point spheroids are completely different near the wall. Tracer prolate spheroids (fibers) orient with their symmetry axis in the streamwise direction, whereas the symmetry axis of oblate spheroids (disks) align in the wall-normal direction. However, inertial point fibers demonstrate tumbling motion with the symmetry axis in the wall-normal/streamwise plane, while disks orient with their symmetry axis perpendicular to the wall-normal direction.⁴⁷ Orientation of finite-size prolates in four-way coupled regime will be investigated here.

Figure 13 depicts the absolute mean direction cosines $|cosθi|$ of P1 and P2 versus distance from the wall for x, y and z directions based on the definitions of Fig. 14, where $θi$ is the angle between the symmetry axis of the spheroid and the relevant direction of the inertial reference frame. It can be immediately seen that finite-size spheroids tend to be aligned with their symmetry axis along the streamwise direction, in particular close to the wall, where particles exhibit very strong preferential orientation in the streamwise direction (x). This preferential orientation decreases gradually toward the channel center. The minor fluctuations appearing on the curves are due to the discrete nature of the solid particles and to the limited amount of information available for statistics. High velocity fluctuations might be responsible for particle alignment along x-direction.^23,49 The peak point corresponds to $θx=41°$ for P1 and 28° for P2. P2 shows higher streamwise alignment due to its higher aspect ratio and smaller size. The probability of particles being oriented in spanwise and normal directions is lower than in streamwise direction. Particle orientation along wall-normal direction (y) is the least probable one, except for a very narrow region in the vicinity of the wall due to particle-wall collisions. Low values of $|cosθy|$ mean that spheroids are mainly confined in x-z plane when they are near the walls. Orientation of P1 in the spanwise direction (z) shows a plateau after $y+≃50$⁠. In the core region, the orientations are more isotropic and random orientation is more probable. This effect is due to higher isotropy of the flow field fluctuations around channel center. However, the preferential orientation is still present, but at a much weaker level. Interestingly, P1 and P2 have almost similar orientations beyond $y+≃110$⁠. Alignment of finite-size spheroids in our study is qualitatively similar to the results obtained for inertial point particles as reported in Refs. 26 and 27. It is expected that increasing the flow Reynolds number would decrease preferential orientation, since the Stokes number of the particles would increase.

The mean values of the absolute particle angular velocities of P1 are presented in Fig. 15a. The rotation rate in the spanwise direction is higher than streamwise and wall-normal directions. This is mainly due to particle collisions with the walls and high shear rates near the walls. The rotation velocities are quite isotropic far from the wall. High rotation rate of spheroids around z-direction together with particle alignment along x-direction implies that tumbling is the most frequent rotation mode of spheroids near the wall. At the same time, particles show a kayaking rotation. Close to channel center, particle rotation is more complex and a unique rotation mode cannot be identified. P2 particles show a similar behavior but with a higher z-rotation rate due to higher aspect ratio and smaller size.

Rms of spanwise angular velocity is still higher than the two other directions (Fig. 15b). The general behavior is similar to the mean component and isotropy is observed around channel center.

The probability density functions of the direction cosine of the P1 particle orientation angle, pdf$(|cosθi|)$⁠, are shown in Fig. 16. Solid lines are a 3^rd-order curve fit. This figure confirms the tendency of particles to be aligned in the streamwise direction. It is interesting to see that P1 has the same probability of alignment in all directions at $|cosθ|≈0.6$ or $θ=53°$⁠. Before this point, pdf$(|cosθy|)$ is higher than the two other directions. Beyond this point, the probability of alignment along x-direction suddenly rises.

The structure of the observed vortices in particle-free and particle-laden flows is illustrated in Fig. 17. There, the Q-criterion is used to plot the iso-surfaces at $Q/(u4/ν2)=0.006$⁠. It is seen that, addition of particles introduces many supplementary small-scale vortices in the near-wall region.

In this paper, IB-LBM has been used to simulate turbulent channel flows seeded with finite-size prolate spheroids. Test cases with two spheroidal particles of aspect ratio $λ=2$ and 4 have been successfully investigated. Comparing the obtained results with those of spherical particles it is found that:

• Finite-size spheroids have less velocity reduction effect compared to spheres. At a volume fraction of $ϕ=1.5%$⁠, spheroids reduce the mean velocity by only 2 and 1.6% compared to 4.6% for spheres. Increase of particle aspect ratio and decrease of its size increase the mean fluid streamwise velocity.

• Streamwise turbulence attenuation occurs for all particle types. However, the effect is more pronounced for spherical ones. For all diagonal directions, turbulence intensity for spheroids is closer to single-phase flow.

• Although spheres show a local peak of volume fraction near the wall, this is not clearly the case for spheroids. Local volume fraction of spheroids increases gradually before reaching a plateau far from the wall.

• Spheroids are mainly aligned along the streamwise direction. This preferential orientation is stronger close to the walls and increases with particle aspect ratio. Around channel center, particles show much less preferential orientation.

• Increasing the volume fraction of spheroidal particles decreases the mean streamwise velocity in the buffer layer but increases it close to the channel center. It is expected that even higher particle volume fractions might induce velocities higher than the single-phase one. Investigating what happens in very dense regimes is the next step of this study.

The financial support of the German research foundation (DFG) within the graduate college for “Micro-macro-interactions in structured media and particle systems” (GRK 1554) is gratefully acknowledged.