Ellipsoids evaluated by Zastawny et al. .
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.
- kinematic simulations
- Lagrangian tracking
- nonspherical particles
- response behavior
- preferential orientation
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  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 . 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 . 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,
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 . 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 . In the creeping flow regime, also with particle Reynolds number much lower than 1, Bläser  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 , isometric irregular particles , cylinders and plates , discs , and discs and cylinders . 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 , and the sphericity, defined as the ratio between the surface of the spherical particle equivalent diameter and the actual surface of the nonspherical particle. None of such correlations takes into account the dependence of the drag coefficient with the particle orientation in the flow. There exist some correlations that consider such dependence [15, 16, 17] and, therefore, they are appropriated to be implemented in a Lagrangian computation scheme. Nevertheless, corresponding results for lift and pitching torque coefficients are still not equally available. One of the first results for the different coefficients in terms of the orientation of elliptic particles was obtained by Hölzer and Sommerfeld  using the lattice Boltzmann method (LBM). Vakil and Green  used DNS to study the drag and lift coefficients of cylindrical particles depending on their orientation and aspect ratio, for Reynolds number up to 40, providing a correlation for them. In such work, the underlying flow field was assumed to be uniform and the flow around the particle was completely resolved. More recently, Zastawny et al.  applied DNS in the frame of the implicit mirroring immersed boundary (MIB) method to obtain the flow coefficients for four different ellipsoids. The authors provide specific correlations for the drag, lift, pitching torque, and rotational torque coefficients depending on the particle Reynolds number and orientation but without including the effect of the aspect ratio. The covered Reynolds number range was up to 300. Ouchene et al.  determined with DNS the drag, lift, and pitching torque coefficients for prolate ellipsoids with aspect ratio up to 32 and adjusted their results to proper correlations that include the effects of particle orientation and aspect ratio up to a Reynolds number of 300.
The first numerical computations of very small nonspherical particles in pseudoturbulent flow were performed by Fan and Ahmadi  and Olson . The hydrodynamic forces and torques were computed by the theoretical coefficients of the Stokes regime. Olson  estimated the time step for the translation and rotation motions in function of the fiber length, obtaining the corresponding dispersion coefficients. Fan and Ahmadi  showed that the dispersion of both, translation and rotation, was reduced with the fiber length. However, Olson  found a different result in the case of ellipsoidal particles. Lin et al.  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. , Mortensen et al. , and Marchioli et al.  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.  and Ouchene et al.  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 , a fact that properly illustrates the importance of the correct modeling of nonspherical particles motion. In a later work , other forces such as added mass and pressure force were also included. Drag coefficient was computed using the Ganser  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: is the inertial frame; is the particle frame, whose origin is in the particle center of mass and its axes are the particle principal axes; and is the comoving frame, which has its origin at the same point than particle frame but its axes are parallel to the inertial frame axes. In the particle frame, the axis coincides with the particle symmetry axis and its position with respect to the comoving frame determines particle orientation.
is the orthogonal matrix that performs the transformation. Its components are the direction cosines of the particle axes in the comoving frame, written in function of Euler angles . Such Euler angles are defined according to the x-convention of :
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 are used instead:
And the transformation matrix is written as :
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  are written as:
Here, is the mass of the particle, is the translational velocity of the particle center of mass, referred to the inertial frame, and is the external forces acting on the particle. The moments of inertia with respect to the particle frame axes are , and are the torques experienced by the particle. It should be remarked that the equations for the translation motion are computed in the inertial frame but those of the rotation motion are expressed in the particle frame. In case of the torque experienced by the particle, it has two contributions: the pitching torque, due to the noncoincidence of the particle center of mass and center of pressure (same fact that happens in an airfoil), and the rotational torque, due to the viscous resistance experienced by a rotating body inside a fluid, generated by the differences between fluid and particle rotational velocities.
In addition to a sphere, the four ellipsoids of Zawstawny et al.  have been chosen. They have different sphericities and aspect ratio (see Table 1). In Table 1,
|Ellipsoid 1 (prolate)||0.88|
|Ellipsoid 2 (prolate)||0.99|
Using DNS for ellipsoidal particles immersed in a uniform flow, Zastawny et al.  determined correlations for the flow coefficients (drag , lift , pitching torque and rotational torque ). Such coefficients are written as :
Here, is the volume equivalent particle diameter or the diameter of a sphere with the same volume as the considered particle. The relative fluid velocity with respect to the particle is and is the fluid relative rotation with being the particle angular velocity.
The developed correlations depend not only on particle Reynolds number and particle rotation number , but also on orientation . They are written as :
Moreover, for the cylinders, Vakil and Green  developed correlations for drag and lift coefficients depending on orientation, Reynolds number based on its diameter and aspect ratio
The correlations are expressed as:
Coefficients and those in functions :
can be found in Vakil and Green .
However, expressions for the pitching and rotational torque coefficients are not provided in . Therefore, for the cylinders, the approach of  has been assumed. In , the distance between the center of mass and the center or pressure in a cylinder, , in terms of
Then, the pitching torque is just the cross-product between the particle orientation unitary vector and the resultant force acting on it times . Nevertheless, this torque is computed in the inertial frame of reference, so it should be transformed to the particle frame before being included in Eq. (7) to calculate the particle angular velocity. The approach to compute the viscous rotational torque is to integrate along the particle length the torque due to the drag force with respect the particle center of mass and it is described in . This torque is given directly in the particle frame.
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  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 of a particle released at point at time
Here, is the Eulerian velocity field. If it is known, it is possible to solve Eq. (18); however, finding is not an easy task. One possibility is to work with Lagrangian statistics but then it would be needed to close the relevant Lagrangian correlations. Another option to solve Eq. (18) is to use DNS to obtain ; however, this is computationally very expensive. A much more economical alternative is the use of kinematic simulation (KS) to compute the Lagrangian characteristics of turbulent flow fields. In this technique, stochastic fluid velocity fields are constructed in such a way that their statistical properties are in agreement with those extracted from experiments or reliable DNS. The main advantage of KS is that it employs an explicit continuous formula for computing , so it is not needed to perform interpolation of the fluid velocity field. Moreover, KS results of two particle statistics in HIT have been validated versus DNS showing good agreement .
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 :
represents the n-th wave number; coefficients , are random, uncorrelated vectors perpendicular to , whose amplitudes are chosen according to the prescribed energy spectrum
is the n-th frequency, which determines the unsteadiness of the corresponding mode; it is written proportional to the eddy-turnover time of the n-th mode:
Here, is a parameter of order 1 that governs the unsteadiness of the velocity field. In three-dimensional HIT flows, it has been demonstrated  that the statistical characteristics of two-particle diffusion are independent of . In particular, in this work, two values of the unsteadiness parameter of 0 and 0.5 have been tested, without significant differences in the computed statistical properties. Therefore, following the suggestions of , the value 0.5 has been adopted in the present computations.
Because of the construction of the velocity field given by Eq. (19), , it is solenoidal trajectory by trajectory. Moreover, as shown in , such field includes in each realization turbulent-like patterns as eddying, straining, and streaming regions.
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.  demonstrated that, as the fluid turbulence is homogeneous, particle RS only depend on time and they can be written as , that is, they are an isotropic tensor. Moreover, in the asymptotic limit, , where
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 . Such velocity field is characterized by a fluctuating velocity
The regular nonspherical particles studied have been the ellipsoids in  and the cylinders in . In all cases, particles have the same particle volume equivalent diameter , hence much smaller than . Therefore, such particles can be thought as immersed in a uniform flow field. The Stokes number has been modified by adjusting the material density of particles being the Stokesian particle relaxation time defined as:
If the Kolmogorov time scale is taken as the fluid time scale, particle Stokes number is defined as . According to this nondimensional number, three particle inertia classes are considered: light (), intermediate (), and heavy (). However, as cases with
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 105 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 for translational motion is expressed as:
represents the time delay. With the objective of making results independent of particle injection conditions, statistics are started to be collected after 2 s. The obtained results for the Lagrangian autocorrelation function (LAF) of the ellipsoids of Zastawny et al.  are presented in the left part of Figure 3, including the results for spherical particles, whereas those of the cylinders of Vakil and Green  are in the right side of such figure. Horizontal axis is the nondimensional time delay, . As the lighter particles have an autocorrelation function nearly equal to that of the fluid (tracer limit), the corresponding curves are not shown in Figure 3. Therefore, only the curves for intermediate and high inertia particles are presented. Moreover, in Figure 3 also the fluid Lagrangian (brown curve) and Eulerian (cyan curve) ’s are included for comparison.
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 ), as a result of their smaller responsiveness to the turbulent fluctuations. Same as in , all particle LAFs are mainly in between the fluid LAF and Eulerian autocorrelation function (EAF). As a consequence of inertia, for the smallest values of , the heavy overcomes the fluid EAF, differently from  who considered only noninertial particles.
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 (
In an analogous way to translational LAF, a rotational autocorrelation function (RAF) can be defined in terms of the time delay
where the particle angular velocity is denoted by .
The obtained results for the Lagrangian rotational autocorrelation function of the ellipsoids of Zastawny et al.  are shown in the left part of Figure 4, whereas those of the cylinders of Vakil and Green  are in the right side of such figure. Again, horizontal axis is the nondimensional time delay, . Similar to the case of translational motion, the angular velocities of heavy particles keep correlated for longer times than those of lighter particles. Such correlation time for ellipsoidal particles is much shorter than that of the translational motion. Also, for all inertia cases, the RAF of disc-like particles drops quicker than for the prolate ellipsoids. As mentioned for the translational correlations, the RAF curves for the prolate ellipsoids collapse for the lighter particles, but they show noticeable differences for the intermediate and large inertia particles demonstrating an effect of the aspect ratio on . Ellipsoid 2, with the smaller
In the case of cylinders (Figure 4, right), the RAF curves for all
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, , where is the fluid rms fluctuating translational velocity and denotes the same quantity but for the particles. The aspect ratio is in the horizontal axis, whereas the different curves correspond to the various inertia cases. On the one side, as it could be anticipated, reduces with increasing particle inertia because the more inertial particles are not able to follow all fluid velocity fluctuations. In fact, as it was found for the LAF, the less inertial particles present, for both ellipsoids and cylinders and for all values of
As inertia increases, decreases monotonically, as expected. However, it increases with growing aspect ratio for both cylinders and ellipsoids. There is one exception that spherical particles have values of slightly above those of the ellipsoid with
Figure 5(b) presents the behavior of the particle angular rms velocity, . In homogeneous and isotropic turbulence, the angular velocity of spherical particles is zero because of viscous damping and the absence of pitching torque. The situation is different in nonspherical particles because in them the geometrical and pressure centers do not coincide, so a net pitching torque is produced that promotes nonzero angular velocities. As illustrated in Figure 5(b), increases with
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 (
Newsom and Bruce  analyzed the influence of turbulence on preferential alignment of quite elongated fibers with . As explained by , preferential orientation of such fibers falling through a still fluid can only be clarified if fluid inertial effects are considered. In the Stokes regime, Khayat and Cox  demonstrate that the force distribution on the fiber is symmetrically distributed along its axis, independent of its orientation regarding the flow and, as a result, the fiber experiences a zero net torque. Beyond the Stokes regime,
Previous reasoning is valid too for another kind of nonspherical shapes as disc-like, cylindrical, or ellipsoidal . For high Reynolds numbers, that is,
Let now be the angle formed by the relative velocity, , and the particle symmetry axis,
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
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.
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
The financial support of Universidad Autónoma de Occidente is gratefully acknowledged.