Ellipsoids evaluated by Zastawny et al. [20]. *a* and *b* are the major and minor semiaxis, respectively.

## Abstract

In this study, the responsiveness of nonspherical particles, specifically ellipsoids and cylinders, in homogeneous and isotropic turbulence is investigated through kinematic simulations of the fluid velocity field. Particle tracking in such flow field includes not only the translational and rotational components but also the orientation through the Euler angles and parameters. Correlations for the flow coefficients, forces and torques, of the nonspherical particles in the range of intermediate Reynolds number are obtained from the literature. The Lagrangian time autocorrelation function, the translational and rotational particle response, and preferential orientation of the nonspherical particles in the turbulent flow are studied as function of their shape and inertia. As a result, particle autocorrelation functions, translational and rotational, decrease with aspect ratio, and particle linear root mean square velocity increases with aspect ratio, while rotational root mean square velocity first increases, reaches a maximum around aspect ratio 2, and then decreases again. Finally, cylinders do not present any preferential orientation in homogeneous isotropic turbulence, but ellipsoids do, resulting in preferred orientations that maximize the cross section exposed to the flow.

### Keywords

- kinematic simulations
- Lagrangian tracking
- nonspherical particles
- response behavior
- preferential orientation

## 1. Introduction

Nowadays, the use of numerical simulation techniques to assist the development and optimization of industrial processes dealing with turbulent multiphase flow has been included as one more step in their layout. Examples of them include pneumatic conveying, fluidized bed reactors, cyclones, classifiers, or flow mixers. Industrial sectors where such processes are important are the chemical, food, or paper industries as well as electric energy production. Due to the complexity of the involved flow, a great majority of simulations are carried out under Reynolds-averaged Navier-Stokes (RANS) in connection with an appropriate turbulence model to describe the turbulent dynamics of the carrier phase.

Two main frames are employed for the description of complex multiphase flows: the two-fluid model or Euler-Euler and the discrete particle models or Euler-Lagrange. In both of them, particles are approximated as point masses being transported in the carrier phase flow field; the solution of the flow around individual elements is usually too expensive and cannot be afforded. In the two-fluid model, both phases are conceived as two interpenetrating continua [1] whose properties are described by sets of partial differential equations. In the Euler-Lagrange approach, the discrete elements are considered as individual objects whose dynamics is governed by a Lagrangian motion equation. Therefore, to obtain the discrete phase variables in the computational domain, a large enough number of discrete element trajectories must be computed. On each particle, appropriate forces act reflecting the various microprocesses taking place at the element scale such as fluid-particle turbulent interaction, particle-(rough) wall interactions, and interparticle collisions [2]. Such technique is especially appropriate for the description of disperse multiphase flow, where usually particles have a size distribution, in confined domains where particle-wall collisions play a predominant role as pneumatic conveying, separation, and classification processes. Both techniques, two-fluid model and Euler-Lagrange, have been applied mainly considering spherical particles. This means that the forces due to the flow (drag and lift) as well as the microprocesses modeling, wall-particle and inter-particle interactions, are assumed to be for spherical-shaped elements [3]. In practical situations, however, nonspherical particles are encountered, either of irregular shape, either with well-defined shapes (fibers or granulates). For example, the paper industry uses large amounts of turbulent liquid to handle and transport the fibers that compose the paper pulp. Besides, such particles in the flow experience particle Reynolds numbers larger than one, *Re* > 1. For such particles, the most relevant transport mechanisms such as aerodynamic transport, wall-particle interactions, and interparticle collisions are substantially different than those for spherical particles.

With the objective of performing the numerical simulation of turbulent flows laden with nonspherical particles, additional information about the forces and torques due to flow (drag and lift forces and pitching and rotational torques due to the shear flow and particle rotation) is needed. It is known that for regular nonspherical particles, that is, ellipsoids or fibers, such forces depend on particle orientation with respect to the flow. For instance, fiber orientation plays a major role in chemical processes as injection, compression molding, or extrusion in which the mechanical properties of the suspensions are determined by the orientation distribution.

For the Stokes regime, particle Reynolds number much lower than 1, the behavior of the nonspherical particles can be determined by analytical methods. Forces and torques acting on an ellipsoidal particle were analytically computed by Jeffery [4]. In a series of papers, Brenner determined the forces due to the flow acting on arbitrary-shaped nonspherical particles in the Stokes regime under different flow configuration by means of theoretical methods [5]. In the creeping flow regime, also with particle Reynolds number much lower than 1, Bläser [6] computed the forces acting of the surface on an ellipsoid in free motion for different flow situations, which allow him to suggest a simple criterion for particle breakup.

The drag coefficient for particle Reynolds number higher than 1 must be obtained by experiments, physical or numerical, as the analytical methods are not applicable any more.

The experimental studies to determine the drag coefficients for nonspherical particles employ wind tunnels or sedimentation vessels. For a moderately wide particle Reynolds number range, there exist results for thin discs [7], isometric irregular particles [8], cylinders and plates [9], discs [10], and discs and cylinders [11]. Drag coefficients were developed in all cases only for certain particle orientations. Compiling such results, different correlations have been developed in terms of particle shape [12, 13, 14]. As representative parameter of the particle shape, two options appear to be dominant: the spherical particle equivalent diameter

The first numerical computations of very small nonspherical particles in pseudoturbulent flow were performed by Fan and Ahmadi [22] and Olson [23]. The hydrodynamic forces and torques were computed by the theoretical coefficients of the Stokes regime. Olson [23] estimated the time step for the translation and rotation motions in function of the fiber length, obtaining the corresponding dispersion coefficients. Fan and Ahmadi [22] showed that the dispersion of both, translation and rotation, was reduced with the fiber length. However, Olson [23] found a different result in the case of ellipsoidal particles. Lin et al. [24] investigated numerically the distribution of the orientation of the fibers in a developing mixing layer, comparing the obtained results with experiments. The fiber length was smaller than the Kolmogorov scale, so they employed the forces due to the flow of the Stokes regime. Zhang et al. [25], Mortensen et al. [26], and Marchioli et al. [27] studied the transport and deposition of ellipsoidal particles in a turbulent channel flow using direct numerical simulation (DNS). Again the hydrodynamic forces and torques were computed with Stokes regime expressions. Beyond the Stokes regime, van Wachem et al. [28] and Ouchene et al. [29] studied a turbulent channel flow laden with ellipsoidal particles using LES and DNS, respectively, employing the flow coefficients developed by themselves in previous works.

Rosendahl group developed a model for the numerical computation of cylindrical and superellipsoidal particles in laminar and turbulent flows in the intermediate Reynolds numbers regime [30, 31, 32]. Particle angular velocity and orientation were computed by means of the Euler parameters. Using a linear relationship between the drag coefficient and the ellipsoid parameters, it was possible to establish a correlation valid up to Reynolds numbers of 1000. To estimate the influence of the orientation, a correlation between the maximum (90°) and minimum (0°) drag was employed. Drag force was calculated using the projected area perpendicularly to the flow. The lift force was expressed in function of the drag coefficient and particle orientation. Other lift forces, such as those due to the fluid velocity gradients or particle rotation, were not considered. The study case was a combustion chamber with straw particles, which were quite well approximated by cylinders. It was found that straw particles were better dispersed than spheres [30], a fact that properly illustrates the importance of the correct modeling of nonspherical particles motion. In a later work [32], other forces such as added mass and pressure force were also included. Drag coefficient was computed using the Ganser [15] correlation, making it possible to numerically compute the biomass combustion chamber.

This contribution aims to study the motion of nonspherical particles immersed in homogeneous isotropic turbulent (HIT) velocity fields built from kinematic simulation at moderate Reynolds numbers. Computations were performed in a tailored in-house code. Properties analyzed include the Lagrangian time autocorrelation function, the translational and rotational particle response, and preferential orientation of the nonspherical particles in the turbulent flow, all of them in terms of particle aspect ratio and inertia.

## 2. Governing equations

### 2.1 Coordinate systems

To build the trajectory of a regular nonspherical particle, it is necessary to solve for its translational as well as rotational motion. However, whereas translation is solved in an inertial frame, rotation is solved referred to the so-called particle frame. Thus, the relevant coordinate frames and the transformations between them have to be introduced.

Figure 1 illustrates, in the case of a cylindrical particle, the employed coordinate systems:

Goldstein [33] gives the transformation between the comoving and particle coordinate systems, which is frequently employed in regular nonspherical particle tracking [31].

The time evolution of such Euler angles depends on the particle angular velocity regarding the particle frame axes. However, there is a difficulty in the sense that such time evolution equations present an unavoidable singularity. Therefore, instead of the Euler angles, the Euler parameters

And the transformation matrix

In the present study, the initial particle orientations are assigned by means of the Euler angles. From them, the corresponding Euler parameters are computed by Eq. (3), and with them, the initial transformation matrix is evaluated using Eq. (4). The Euler parameters evolve in time following Eq. (5), where the particle angular velocities are expressed in the particle frame of reference

### 2.2 Particle motion equations

The nonspherical particle motion equations in a general fluid flow [34] are written as:

Translational motion:

Rotational motion:

Here,

In addition to a sphere, the four ellipsoids of Zawstawny et al. [20] have been chosen. They have different sphericities and aspect ratio (see Table 1). In Table 1, *a* denotes the major semiaxis and *b* the minor semiaxis.

Shape | Aspect ratio | Sphericity |
---|---|---|

Ellipsoid 1 (prolate) | 0.88 | |

Ellipsoid 2 (prolate) | 0.99 | |

Disc (oblate) | 0.62 | |

Fiber (prolate) | 0.73 |

Using DNS for ellipsoidal particles immersed in a uniform flow, Zastawny et al. [20] determined correlations for the flow coefficients (drag

Here,

The developed correlations depend not only on particle Reynolds number

*a* values are listed in [20] as also coefficients *b*, *c*, *r*:

Moreover, for the cylinders, Vakil and Green [19] developed correlations for drag and lift coefficients depending on orientation, Reynolds number based on its diameter *AR*. In this case, such coefficients are expressed in terms of cylinder length L and diameter D:

The correlations are expressed as:

Coefficients

can be found in Vakil and Green [19].

However, expressions for the pitching and rotational torque coefficients are not provided in [19]. Therefore, for the cylinders, the approach of [31] has been assumed. In [31], the distance between the center of mass and the center or pressure in a cylinder, *AR* and

Then, the pitching torque

Particle motion equations and correlations for cylinders and ellipsoids presented in this section have been implemented in an in-house code. The numerical integration of the ordinary differential equations that govern the motion of nonspherical particles has been performed by a fourth-order Runge-Kutta method, with small enough time steps to avoid numerical instabilities [35, 36]. The fluid velocity field in which particles are immersed has been built by the kinematic simulation technique described in the next section. It is known that Runge-Kutta methods do not satisfy the time-reversal property, a fact that makes such methods inappropriate for integrating energy-conserving systems, for instance. However, particle equations are dissipative systems (as they include viscous drag forces) and, for them, Runge-Kutta algorithms can be used [37] provided that the time step is small enough to keep the errors bounded.

## 3. Kinematic simulation

There exist different options to calculate the Lagrangian properties in a turbulent flow. The starting point is the trajectory equation in which the position *t* = 0 is calculated solving:

Here,

The three-dimensional Eulerian velocity field to be employed in Eq. (18) is built as a series of random Fourier modes. The velocity field is solenoidal at each realization by construction. Moreover, the energy spectrum of the Fourier modes is prescribed, for example, by a power law, so the effects of small flow scales on Lagrangian statistics are directly included. Such KS velocity field is written as [39]:

*E*(*k*) [39]. Here, the energy spectrum has been the Kolmogorov decay law of −5/3.

Here,

Because of the construction of the velocity field given by Eq. (19),

To validate the spherical particles tracking in the KS velocity field, the values of particle Reynolds stresses (RS) in HIT have been selected. In this configuration, Hyland et al. [41] demonstrated that, as the fluid turbulence is homogeneous, particle RS only depend on time and they can be written as *β* is the inverse of particle relaxation time (see Eq. (21) below) and

## 4. Numerical simulation

Computations were performed in a tailored in-house code. The turbulent velocity field generated with KS resembles one of the fields worked in [39]. Such velocity field is characterized by a fluctuating velocity *u’* = 1 m/s, a fluid Reynolds number of 10^{4} resulting in a Kolmogorov length scale

The regular nonspherical particles studied have been the ellipsoids in [20] and the cylinders in [19]. In all cases, particles have the same particle volume equivalent diameter

If the Kolmogorov time scale *Re* > 1 are considered, an effective particle relaxation time is introduced as

Simulations proceed in the following way: for each KS realization of HIT fluid velocity field, a particle is located in the center of the domain with zero initial velocity; particle translational and rotational motion is computed from Eqs. (6) and (7), its orientation is calculated from Eq. (1), and its trajectory is built; particle tracking lasts for around 10 fluid integral time scales; and particle properties are stored every second for evaluation. Such process is carried out a sufficient number of times to reach significant statistical results. In this study, statistics has been performed based on 10^{5} KS realizations.

In the following section, the results of the particle Lagrangian time autocorrelation function, the translational and rotational particle response, and preferential orientation of the nonspherical particles in the turbulent flow are analyzed as function of their shape and effective Stokes number.

## 5. Results and discussion

The Lagrangian autocorrelation function

As it can be readily seen from Figure 3, higher inertia particles are characterized by larger integral Lagrangian time scales (ILTSs) (defined as the integral up to infinity of

Moreover, for intermediate inertia, the curves for all particle shapes nearly collapse in a single curve. On the other hand, a shape effect is noticeable for the heaviest particles, where the various shapes present differences in their curves. It is interesting to realize that ILTSs of higher aspect ratio (*AR*) are below those of smaller *AR*, for both ellipsoids and cylinders. This effect is a Reynolds number effect due to the dependence of drag coefficient on shape and *AR*: an interaction between translation and rotation motions occurs that results in a spreading of the particle effective Stokes number. As a consequence, particles with higher Reynolds numbers also have larger effective inertia (reflected on an increased Stokes number) and, therefore, their LAF decreases slower, implying a higher ILTS. In the ellipsoids case, it happens that those of typ. 2 present a

In an analogous way to translational LAF, a rotational autocorrelation function (RAF) *τ* as:

where the particle angular velocity is denoted by

The obtained results for the Lagrangian rotational autocorrelation function of the ellipsoids of Zastawny et al. [20] are shown in the left part of Figure 4, whereas those of the cylinders of Vakil and Green [19] are in the right side of such figure. Again, horizontal axis is the nondimensional time delay, *AR*, has the higher RAF curve of all prolate ellipsoids, while Ellipsoid 1 and the fiber have very similar rotational correlation functions.

In the case of cylinders (Figure 4, right), the RAF curves for all *AR* and inertias are different. For the smallest inertia particles, RAF decreases with increasing *AR*, similar to what was found for LAFs in the translational motion. Moreover, the *AR*; however, as Figure 4 (right) suggests, an asymptotic value for *L*/*D* seems to exist because the curves for *AR* = 10 and 20 are very close to each other. The change of behavior of the RAF with increasing inertia could be due to the fact that for the small and intermediate inertia, the particle relaxation time for rotation reduces with growing *AR*, whereas for heavy particles, such relaxation time behaves in the opposite way. Nevertheless, this fact must be further investigated, possibly using fully resolved simulations.

Next, the response of the nonspherical particles to the fluid fluctuating velocities is analyzed for both translation and rotation motions. Figure 5(a) shows the behavior of the particle’s relative linear root mean square (rms) velocity, that is, *AR*, the same fluctuating velocities as the fluid, indicating that they behave as fluid tracers.

As inertia increases, *AR* = 1.25. The trend of increasing particle fluctuating velocities with aspect ratio is consistent with the aforementioned fact that effective Stokes number tends to reduce with growing *AR*; therefore, particles with lower *AR* respond less to the fluid fluctuations than the more elongated ones. This result has been obtained for particles much smaller than Kolmogorov length scale; therefore, such particles can be considered to be immersed in a uniform flow field, just in the same conditions as the flow coefficient correlations were developed. In the study of Hölzer and Sommerfeld [42], nevertheless, a different result was found. Using a DNS based on the lattice Boltzmann method (LBM), Hölzer and Sommerfeld [42] obtained that relative particle fluctuating velocity reduced with increasing *AR*. The main difference with the present work is that Hölzer and Sommerfeld [42] employed particles with size well above *AR*. Let us remark that both results are not conflicting as the size range of the employed particles in the two studies is very different. Further work combining DNS with nonspherical point particles smaller than

Figure 5(b) presents the behavior of the particle angular rms velocity, *AR* and reaches up to a maximum of *AR*, the moment of inertia along the major axis reduces and along the minor axis increases, which would lead to higher and lower

In the following, the correlation relative velocity direction-particle orientation is analyzed depending on inertia and aspect ratio. A well-known fact is that regular nonspherical particles falling through a still liquid at intermediate Reynolds numbers tend to be oriented in a determined direction. Cylinders and prolate ellipsoids are prone to keep their symmetry axis (*z’* in Figure 1) perpendicular to the flow, thus maximizing drag. Differently, discs and oblate ellipsoids tend to move with the symmetry axis aligned with the flow, also maximizing drag [43]. However, spheroidal Stokes particles only show a preferential orientation if a persistent velocity gradient exists [27]. Therefore, in HIT flow where there are no mean velocity gradients, a Stokes particle will not have any preferred orientation.

Newsom and Bruce [44] analyzed the influence of turbulence on preferential alignment of quite elongated fibers with *Re* > 1, when the fiber has an oblique orientation with respect the flow, such distribution of the force is not any more symmetric and it experiences a net pitching torque. Such torque will promote a rotation that drives the fiber to adopt an orientation where its symmetry axis is orthogonal to the relative flow. Interestingly, if the fiber is oriented orthogonal or parallel to the flow, the net experienced torque is zero, due to the symmetry of the force distribution; however, in the first case, this situation is stable, while in the second case, where the centers of gravity and pressure do not coincide, this situation is unstable.

Previous reasoning is valid too for another kind of nonspherical shapes as disc-like, cylindrical, or ellipsoidal [45]. For high Reynolds numbers, that is, *Re* > 100 appears a secondary motion overimposed to the particles predominant movement direction. Such secondary motion is promoted by a wake instability and vortex detachment from the rear surface of the particles. Two main kinds of secondary motion can be observed: large quasi-periodic swings along the main path, and a more or less chaotic tumbling forming a definite angle with the main motion direction. There is a coupling between this kind of oscillatory motion and the wake instability [10]: a vortex detachment follows at the end of a particle swing. Nevertheless, in the present study, such secondary motions do not appear as the considered particle Reynolds number is not large enough, that is, *Re* < 40.

Let now *z’*. Therefore, ^{5} particles. Computed values of

As it is observed in Figure 6(a), it is found that prolate ellipsoids do manifest preferential orientation with respect to the relative velocity. Of course, spherical particles do not have a preferred orientation and the corresponding Pdf is a horizontal line (black color). Prolate ellipsoids have a preference for orientating its symmetry axis orthogonal to the relative flow, tending to maximize the drag, similar to what occurs in particle sedimentation studies. The orientation preference increases with inertia, which is quite similar for all aspect ratios considered in this study.

On the other hand, as it is presented in Figure 6(b), cylinders seem not to have any preferred orientation in the HIT KS velocity field, being all the curves pretty flat. Only for the case of higher *AR* and lowest inertia, the curves show a trend to be slightly higher for values of

For the discs, Figure 7 shows that there is a clear trend of the particle symmetry axis to be aligned with the relative flow, again maximizing drag, similar to the results obtained for sedimenting particles in stagnant fluid. Such trend is more marked when particle inertia increases.

## 6. Conclusions

In this study, regular nonspherical particle responsiveness to HIT flows has been investigated in combination with KS of fluid velocity field. The main results obtained are the following: the particle LAF reduces when particle *AR* is augmented, because effective particle inertia decreases if aspect ratio increases; this is true for both translational and rotational time autocorrelation functions. In the case of cylinders, *AR* in the considered case of particles much smaller than Kolmogorov length scale; such behavior is contrary, although not conflicting, to the findings of [42] for fully resolved particles with sizes larger than the Kolmogorov length scale. For both ellipsoids and cylinders, the particle angular rms velocity first increases with aspect ratio, reaches a maximum of

## Acknowledgments

The financial support of Universidad Autónoma de Occidente is gratefully acknowledged.