On the motion of non-spherical particles at high Reynolds number

Powder Technology 202 (2010) 1–13 Contents lists available at ScienceDirect Powder Technology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / p ow t e c Review On the motion of non-spherical particles at high Reynolds number Matthias Mandø ⁎, Lasse Rosendahl Aalborg University, Institute of Energy Technology, Pontoppidanstraede 101, 9220, Aalborg East, Denmark a r t i c l e i n f o Article history: Received 13 August 2009 Received in revised form 29 March 2010 Accepted 3 May 2010 Available online 7 May 2010 Keywords: Non-spherical particles Particle equation of motion Gas–solid interaction Dispersed multiphase flow a b s t r a c t This paper contains a critical review of available methodology for dealing with the motion of non-spherical particles at higher Reynolds numbers in the Eulerian–Lagrangian methodology for dispersed flow. First, an account of the various attempts to classify the various shapes and the efforts towards finding a universal shape parameter is given and the details regarding the significant secondary motion associated with nonspherical particles are outlined. Most investigations concerning large non-spherical particles to date have been focused on finding appropriate correlations of the drag coefficient for specific shapes either by parameter variation or by using shape parameters. Particular emphasis is here placed on showing the incapability of one-dimensional shape parameters to predict the multifaceted secondary motion associated with non-spherical particles. To properly predict secondary motion it is necessary to account for the noncoincidence between the center of pressure and center of gravity which is a direct consequence of the inertial pressure forces associated with particles at high Reynolds number flow. Extensions for non-spherical particles at higher Reynolds numbers are far in between and usually based on semi-heuristic approaches utilizing concepts from airfoil theory such as profile lift. Even for regular particles there seems to be a long way before a complete theory can be formulated. For irregular particles with small aspect ratio, where the secondary motion is insignificant compared to the effect of turbulence, the drag correlations based on onedimensional shape parameters come to their right. The interactions between non-spherical particles and turbulence are not well understood and modeling attempts are limited to extending methods developed for spheres. © 2010 Elsevier B.V. All rights reserved. Contents 1. Introduction . . . . . . . . . . . . . . . . 2. Classification of shape . . . . . . . . . . . 3. Drag correlations for translational motion . . 4. Classification of regimes of secondary motion 5. Orientation dependent models . . . . . . . 6. Interaction with turbulence . . . . . . . . . 7. Summary/conclusions . . . . . . . . . . . Appendix A. Equations of motion for non-spherical References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. Introduction Irregular non-spherical particles are found in most industrial particulate flows and similarly most engineering flows are turbulent. However, the vast majority of scientific investigations of particulate flows assume particles to be perfectly spherical particles. The exact governing equations for turbulent flow have been known for over a century but the utilization of these is significantly impeded by the ⁎ Corresponding author. E-mail addresses: mma@iet.aau.dk (M. Mandø), lar@iet.aau.dk (L. Rosendahl). 0032-5910/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.powtec.2010.05.001 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 3 4 6 10 11 11 12 need to resolve the smallest flow structures and time scales. Consequently, for most practical uses, turbulence is modeled using Reynolds or Favre averaging and the interaction with particles is handled by random walk models. Large non-spherical particles present their own set of particular problems in the context of Computational Fluids Dynamics (CFD): How to define and quantify the shape? How to deal with secondary motion? How well will the methodology developed for spheres work for highly non-spherical shapes? How to handle turbulence? This paper attempts to give an account of the present state of modeling the motion of large nonspherical particles. The relevance of this paper also becomes evident 2 M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 considering the increasing efforts towards the replacement of pulverized coal with biomass in existing and new power plants. Whereas pulverized coal particles are small and the spherical ideal is a fair approximation, pulverized biomass particles are characterized as large and with high aspect ratios due to their fibrous nature. This investigation has been limited to the Eulerian–Lagrangian methodology and to solid non-deforming particles in Newtonian fluids. Please refer to Sommerfeld et al. [68] for an updated outline on existing knowledge concerning multiphase flow and Chhabra [13] for particle motion in non-Newtonian fluids. 2. Classification of shape Particles come in all sort of shapes and sizes, in fact, due to the arbitrary nature of naturally occurring particles there are an indefinite number of possible shapes. This necessitates the need for a set of parameters to aid in the description of different particle shapes for the implementation in the numerical models and the relay of relevant information to other scientists. This information is available in many books on the subject, e.g. Rhodes [61] or Clift et al. [18] contain much useful information, and here we limit ourselves to focus on the most pertinent issues involved in the classification of shapes. One such issue seems to be the terminology used. The word non-spherical most often, and somehow also in the title of this paper, refers to all shapes which are not perfect spheres. However, the true meaning of the word spherical is: “shape which is sphere-like” which thus implies a subjective distinction. In the context of CFD it is useful to use this to objectively distinguish between shapes which with reasonable accuracy can be approximated as spheres and shapes which require a more intricate handling. Another often used terminology to characterize shapes which are not spheres1 is the word “irregular” whose true definition is counter to that of regular shapes. Table 1 outlines the categorization using the above discussed terminology. Table 1 serves as a reference for the following discussion of different simulation strategies in this work. Spherical particles have no or only little secondary motion associated with their trajectories, assume no preferred orientation but rather tumbles and thus it would be justified to base the simulation methodology on that used for spherical particles i.e. not considering orientation or shape induced lift. Possible extensions to the simulation methodology thus revolve around modifications to the drag coefficient considering the shape or by specifying an equivalent diameter and using drag correlations based on spheres. Non-spherical particles on the other hand are associated with shape induced lift, orientation dependent lift and drag forces, significant secondary motion and may assume a preferred orientation depending on the regime of motion. Thus, it is necessary to revise the usual strategy to properly capture these phenomena. This involves keeping track of the orientation and rotation of the particle as well as the formulation of appropriate orientation dependent lift and drag force on a per shape basis. For irregular non-spherical particles common strategies involve the approximation of the shape to a regular counterpart e.g. cylinder for a wood splinter, disk for a flake. The distinction between spherical and non-spherical particles is principally subjective and thus open for interpretation. Here, it is suggested that the distinction is made on the basis of the aspect ratio, β. This simple criterion can be easily measured via microscopy techniques and is a good representative for when secondary effects become important. According to Christiansen and Barker [15] and Clift et al. [18] a suitable value for this criterion is β = 1.7 which also roughly corresponds to the aspect ratio for a cube (based on the diagonal to the side length). Thus, particles below this ratio are considered spherical and can be treated with reasonable accuracy Table 1 One possible categorization of shapes. Regular Irregular Spherical Non-spherical Polygons, spheroids with low aspect ratio Pulverized coal, sand, many powders, particulate matter Cubes, cylinders, disks, tetrahedron, spheroids with high aspect ratio Pulverized biomass, flakes, splinters, agglomerates using a single drag correlation of choice. Particles above this ratio should be classified according to which generic shape they resemble the most e.g. cylinder, disk, spheroid, super-ellipsoid of revolution and treated accordingly. For spherical particles it is only necessary to specify an equivalent diameter and optionally a shape factor to account for the departure in shape from a sphere. Table 2 gives an outline of commonly used diameter definitions after Allen [1]. Note that projected area, Ferets and Martins diameters are determined directly from image analysis while area and volume equivalent diameters often are based on image analysis by assuming a thickness. The other diameters listed correspond to a particular analysis method e.g. Stokes diameter which is found from sedimentation techniques. The main difficulties are thus reduced to a matter of measurement and it seems appropriate to offer a few comments about available measurement methods. The basis for all methods is that they provide the same result when applied to a perfect sphere while marked differences occur as the shape becomes non-spherical due to the differences in the diameter definitions. In some scientific or industrial fields specific methods are prevailing due to the individual strengths of particular methods e.g. sieve analysis is often preferred whenever a wide size distribution is encountered while it may be unsuitable for very fine powders. Aerodynamic separation and sedimentation techniques are used for fine powders and particulate matter which tends to be spherical in nature and due to the diameter definitions the size distributions can be used directly in Lagrangian trajectory calculations using the drag coefficient of a sphere. Due to practical and theoretical considerations these methods are not used for particles larger than 50 μm and thus any discussion concerning large non-spherical particles becomes somewhat irrelevant. Image analysis is regarded as a benchmark compared to other techniques as this involves direct determination of the diameter. However, as image analysis is based on a two-dimensional measurement this method becomes increasingly biased as the particles deviate from the spherical ideal. For example the volume equivalent diameter of flake-like particles will be systematically overestimated if their thickness is assumed to be proportional to their 2D extent while Table 2 Commonly used diameter definitions. Aerodynamic/drag diameter Stokes diameter Projected area diameter Ferets diameter Martins diameter Area equivalent diameter Volume equivalent diameter Sieve/mesh diameter Laser diffraction diameter 1 The appropriate term non-sphere is surprisingly hardly ever used. Diameter of a sphere of unity density with the same terminal velocity as the particle Diameter of a sphere of same density and the same terminal velocity as the particle Diameter of a circle having the same area as the projection of the particle The mean value of the distance between pairs of parallel tangents to the projected outline of the particle The mean chord length of the projected outline of the particle Diameter of a sphere having the same surface area as the particle Diameter of a sphere having the same volume as the particle The width of the minimum square aperture through which the particle will pass Diameter is calculated according to the Mie or Fraunhofer diffraction theory M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 laying flat on a plate [59]. This bias can somewhat be rectified by analyzing images of particles in free fall but the failure to resolve the third dimension ultimately means that this method is associated with significant measurement uncertainty. Full 3D analysis which would involve the use of two or more cameras set at an angle to each other has the potential to provide accurate measurement of the shape and volume. However, the complicity of such a method has so far hindered the implementation into commercially available equipment. Popular methods to determine the size distribution of particles are by sieves and laser diffraction. Although these methods are very different they are associated with similar types of uncertainty. Sieve analysis allows slender particles to pass through a fine mesh compared to the particle volume equivalent diameter while flake-like particles will be stopped by a coarse mesh relative to the particle volume equivalent diameter. Similarly, the laser diffraction technique relates the diameter to the orientation of the particle crossing the measurement volume [3]. As the particle is allowed to rotate freely, this means that a slender particle can be assessed on the basis of its smallest dimension while the flake-like particle might be assessed on the basis of its largest dimension. As such both methods are associated with similar, wider size distributions [55]. Furthermore, the cross-sectional area averaged over all orientations for a non-spherical particle is larger than for an equal volume sphere [23]. Depending on predominant shape in a given sample sieve analysis might predict a larger or a smaller mean value compared to the true mean while laser diffraction tends to overestimate the true size distribution [55,58]. As a final remark on the use of single dimension definitions for non-spherical particles, it may be said that they are often reported indiscriminately and used without any regard to the requirement of the drag correlation [6]. Significant biases are associated with the measurement of particle sizes for non-spherical particles for all measurement methods and the problem severely deteriorates for increasing aspect ratios. By using equivalent diameters all data about the shape of the particle is essentially lost and to be able to retain this information additional shape factors have been suggested to quantify the geometry/ irregularity of non-spherical particles. These can be seen as a parallel to the roughness factors which are commonly used for pipe flow. Note that only image analysis is capable of supplying the additional information regarding the shape of the particle and that this information is often based on the assumption that 2D images of particles can be directly related to the 3D shape of the particle. Table 3 gives an outline of the most commonly used shape factors. These shape factors can be used for both regular and irregular particles, but especially suited for the latter since the shape of irregular particles cannot be expressed in any other way [44]. Many alternative shape factors [18,21,48,72] have also been suggested, but none has won greater acceptance or use despite clamed superiority. Fractal dimensions and harmonics have also been used to characterize the shape/morphology of irregular particles [16,67,76]. However, these have not been used in conjunction with CFD simulations and are not addressed further. Automated algorithms in image processing software allow for quick determination of shape factors as well as dimensions but the accuracy is limited by the previously mentioned assumptions for 2D images of 3D particles. Presently, the most Table 3 Commonly used shape factors [34,75]. Corey shape factor Volumetric shape factor Roundness Sphericity Ratio of the smallest principal length axis of the particle to the square root of the intermediate and longest principle length axis Ratio of the volume of the particle to the diameter of a sphere with the same projected area as the particle cubed Ratio of the average radius of curvature of the corners to the radius of the largest inscribed circle Ratio of the surface of a sphere with the same volume as the particle and the surface area of the actual particle 3 commonly used shape factor is the sphericity, ψ. This does not seem to be due to superior performance when used in correlations of the drag coefficient or because it is easier to measure than other shape factors. A closer look at the formulation of sphericity shows that it represents the inverse of a surface enhancement factor for a sphere with equivalent volume and can thus be used in combusting flows to additionally account for the surface area available for reactions. However, the true reason for the greater popularity of the sphericity is most likely that it simply seems to be the most elegant way to quantify the shape of an arbitrary particle. In lack of significantly better shape factors, evaluated on their ability to correlate the drag coefficient, the more elegant formulation has won predominance. 3. Drag correlations for translational motion Shape factors form the basis for most attempts to describe the motion of spherical and non-spherical particles at higher Reynolds numbers. Most of these correlations employ the volume equivalent sphere diameter, dVeq, as the characteristic size and the sphericity, ψ, to quantify the shape and is thus expressed as: CD = f ðRe; ψÞ ð1Þ where the characteristic size is usually taken as the diameter of a sphere with the same volume as the particle. Five different correlations of the drag coefficient for non-spherical particles have been compared against a large database of independent experimental data in the study of Chhabra et al. [14]. The average error reported varies between 16% and 43% whereas the maximum reported error for all correlations is above 100%. The largest errors are encountered for hollow cylinders and agglomerates of spherical particles. These shapes represent extremes in terms of the sphericity and they have little resemblance with a sphere. The general rule which can be drawn is that the further away from the spherical ideal the shape of the particles is, the poorer the correlations perform. Depending on the flow regime and the shape, particles which have the same value of sphericity might take on very different motion patterns or preferred directions and are thus associated with very different drag coefficients when the projected area used, is that of a sphere with the same volume as the particle. The classical example to illustrate this is by considering particles shaped as cylinders of different length to diameter ratio. The sphericity of a cylindrical particle can be expressed as: ψ=
2 = 3 3 2 β 2 1 + 2β ; β= L D ð2Þ where β is the aspect ratio expressed for a cylinder as the length, L, to the diameter, D, of the cylinder. From this expression it can be realized that both a cylinder with an aspect ratio less than one, commonly referred to as a disk, and a cylinder with an aspect ratio above unity can have the same value of sphericity. From the experimental data from McKay et al. [53] for the drag coefficient of falling cylinders, it can be realized that the drag coefficient for disks is distinctively different from that of cylinders, even when only small aspect ratios are considered. In Fig. 1 the difference between the measured drag coefficient and that calculated on basis of the correlation by Ganser [28], the most accurate of the correlations investigated by Chhabra et al. [14], is indicated as percentage error. To provide a relevant reference, the error obtained from using a correlation developed strictly for spheres [17] and a correlation for freely falling cylinders with finite length in liquids [39] is also indicated. It can be seen that the correlation by Ganser [28] provides an acceptable fit for aspect ratios below unity whereas for aspect ratios above unity the correlation becomes exceedingly poorer. Using correlations developed for specific shapes gives the most 4 M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 Fig. 1. The error using the correlations by Clift and Gauvin [17], by Isaacs and Thodos [39] and by Ganser [28] evaluated by the experimental data for disk and cylinders obtained by McKay et al. [53] with 1000 b Re b 30,000. satisfying result whereas the practice of using the diameter of a volume equivalent sphere in correlations of the standard drag curve for a sphere gives poor results. An overview of references on the drag for regular non-spherical shapes is provided in [18] and this topic is therefore not dealt with further. For irregular shapes, it is recommended to only use correlations based on the sphericity for particle shapes with a sphericity approaching unity. This corresponds to shapes with small aspect ratios and which thus only deviates slightly from the spherical ideal. However, if an investigation centers on a specific shape, the best result is obtained by making an empirical fit of the drag coefficient as a function of the Reynolds number for that specific shape. For distinctive non-spherical irregular particles, where the ratio of the maximum length to the minimum length is above 1.7, the shape should be classified either as rod-like, which can be approximated by a cylinder, or flake-like, which can be approximated by a disk [15]. Again it should be pointed out that although it is possible to specify a drag coefficient for a free falling non-spherical particle, which correlates with the primary translational motion, nonspherical particles are associated with significant secondary motion which again may alter their main trajectory. Note also that, although correlations exist for many regular non-spherical shapes, these are often based on scant data sets, which are associated with considerable scatter, due to the secondary motion. 4. Classification of regimes of secondary motion For spheres expansions of the equation of motion to higher Reynolds numbers are usually achieved by empirical fits of the drag coefficient. For spherical particles this approach works well because the motion pattern is not associated with noticeable secondary motion nor does it assume a preferred direction. The drag coefficient can be expressed by properly accounting for the increase in particle surface area either by using an aerodynamic equivalent diameter or by using shape factors. However, the drag on a non-spherical particle is dependent on its orientation. Primarily, the projected area, on which the drag is based, may differ by several orders of magnitude from one orientation to another but also the drag coefficient varies significantly depending on the orientation. Also, rotational effects are important when considering orientable particles and the equations for conservation of rotational momentum must be taken into consideration as translational motion depends directly on them. Non-spherical particles are associated with characteristic secondary motion depending on the Reynolds number regime and their shape. Moreover, in some Reynolds number regimes particles will take on a preferred direction. Most investigations of the motion of non-spherical particles deal with the generic shapes of ellipsoids, cylinders and disks since these can be made, by parameter variation, to resemble a great number of different shapes. Particles with an oblong shape, such as a prolate ellipsoid or a cylinder, are often used to resemble fibers, while particles with a flat shape, such as an oblate ellipsoid or a disk, can be used to represent flakes. For very low Reynolds number flow, Rep b 0.1 (Stokes flow), both oblong and flat particles in a shear flow will move in slow orbits, known as Jeffery's orbits, after G.B. Jeffery [40] who was the first to describe the motion. One restriction in this analysis is that the particles have to obey certain symmetry conditions which strictly speaking would exclude all irregular particles. Characteristic for nonspherical particles in Stokes flow is that, although they move in orbits the majority of the time, they will be aligned or be at a small angle to the flow [4]. In practical terms it is thus more useful to state that the particles tend to align themselves with the flow. This effect has also been observed for fibers used in the manufacturing of paper and thus it seems sensible to also assume that irregular particles would exhibit this behavior providing that they have a large aspect ratio. The motion of particles in Stokes flow represents the only purely theoretical approach to the motion of non-spherical particles and consequently most investigations on the motion of non-spherical particles dwell on this topic. The motion of particles in creeping flow has been extensively reviewed by Leal [49] and more recently in the work by Carlsson [12] and this topic is therefore not dealt with further. At moderate Reynolds number flow, 0.1 b Rep b 100, inertial effects become important, and a steady recirculation zone starts to build up in the wake of the particles. The pressure distribution on the particle, due to the recirculation zone, forces the particles to align themselves with their maximum cross-section normal to the flow. Generally this effect is more pronounced at higher Reynolds numbers and for particles with a more pronounced non-spherical shape. Since the particles are steadily aligned perpendicular to the flow, empirical data for the generic shapes, such as an infinite long cylinder in cross-flow, may be used to model the motion. For disks expressions for fixed disks in cross-flow can be used directly while for cylinders appropriate corrections for end effects should be applied for cylinders with finite length. High Reynolds number flow, Rep N 100, is characterized by significant secondary motion which is superimposed on the particles' steady fall or rise. The secondary motion is initiated by the onset of wake instability and also signals the beginning of vortex shedding from the wake of the particles. The secondary motion may be in the form of large periodic swings around a mean vertical path or chaotic tumbling which can take place at an angle to the vertical fall or rise direction. The oscillatory motion is coupled with the wake instability and photographic evidence using dye injected in the wake of a falling disk show that the end of each swing is followed by the shedding of a vortex [78]. Besides the Reynolds number the motion patterns have been shown to correlate well with the non-dimensional moment of inertia; here shown for a disk:  I = πρp Idisk β: = 64ρ ρD5 ð3Þ For a disk this is obtained by dividing the moment of inertia with the fluid density and disk diameter times five. Note that the dimensionless moment of inertia for a disk is thus transformed into an expression involving the density ratio and the aspect ratio. These two parameters are often used in the description of the motion of non-spherical particles. As such these two parameters were previously used to correlate the drag coefficient of a freely falling cylinder [39]. The dimensionless moment of inertia can similarly be defined for other regular shapes using the associated moment of inertia and the appropriate characteristic length for that shape. However, it is only for M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 5 oscillating cylinders in free fall Marchildon et al. [52] provide the following empirical fit of the Strouhal number: sffiffiffiffiffiffi ρβ Str = c ρp Fig. 2. Flowmap showing the behavior of disks (L/D b 0.1) as function of the Reynolds number and the dimensionless moment of inertia. Data from Stringham et al. [70], Willmarth et al. [78] and Field et al. [26]. Contours of constant Strouhal number are shown for periodic oscillations. disks that the dimensionless moment of inertia has been used to construct intricate flow maps as seen in Fig. 2. The dimensionless moment of inertia is also known as a stability number [70] and basically indicates the inertial resistance of nonspherical particle to rotate. Thus thin disks, with very low I*, will not perform a full rotation and tends to fall in great arcs of periodic sideward motion for high Re whereas more bulky disks, with high I*, provide little rotational stability and undergo full rotation, tumbling motion, for Re N 150. An intermediate regime where disks switch between periodic oscillations and tumbling motion is referred to as the glide–tumble regime. The different motion patterns for a disk are shown in Fig. 3. The steady fall, periodic oscillation and tumbling regimes can be categorized as stable regimes where the frequency of rotation, in regime IV, approaches the oscillation frequency of regime II. Furthermore, the Strouhal number for disks has been shown to a linear function of the dimensionless moment of inertia for constant Reynolds numbers [78]. The glide–tumble motion is unstable and it can be interpreted as a transition regime for periodic oscillations with intermittent bursts of tumbling [26]. For freely falling cylinders only two distinct motion patterns can be identified. Depending on the Reynolds number cylinders assume either steady falling or periodic oscillations with their maximum cross-section normal to the direction of the flow. For periodic ð4Þ where c = 0.095 and the Strouhal number have been based on the length of the cylinder as the characteristic dimension. This fit has since been verified by Sørensen et al. [69] for a wider range of Reynolds numbers although with a constant of proportionality approximately half of that given by Marchildon et al. [52]. Notice that for cylinders the aspect ratio and the density ratio are still important parameters but the relationship with the Strouhal number is not linear as it was for disks. Following Marchildon et al. [52] analysis it can realized that a Strouhal number based on a characteristic length which is a combination of the length and the diameter can reduce Eq. (4) to a relation only depending on the density ratio. However, the physical significance of such an approach is uncertain. It is also reported that the steady oscillatory motion around the horizontal plane is accompanied with rotation around the axis parallel to the fall direction or even a mean sideward motion [69,70]. There does not seem to be any clear pattern of this and the best explanation is possible that it relates to either initial condition, the release mechanism, or has to do with the vicinity of the walls of the experimental setup. The analysis up to now has only looked at particles with uniform mass distribution. However, large naturally occurring particles are also often characterized by having a non-uniform mass distribution which can be related to cavities in the shape. A prime example is the case of shredded straw which can be described as being mainly hollow, but where the presence of solid nodes seriously disrupts the uniformity of the mass distribution. Generally, for an otherwise symmetric particle, the movement of the center of mass away from the center of geometry acts to turn the particle to fall with its heaviest side downward. Clearly, this significantly alters the motion characteristics and can considerably increase the terminal velocity of that particle [64]. A particle with a non-uniform mass distribution but with a coincident location of the center of mass with the center of geometry, i.e. a straw particle with a node in the middle, will have the same resistance characteristics as the uniformly distributed case [5]. However, since the moment of inertia can be different this has the potential to affect the secondary motion pattern of that particle. Regular and irregular particles with aspect ratio close to unity falls with no preferred orientation and with a motion pattern which best can be described as tumbling. Indeed, if a dimensionless moment of inertia was calculated based on Eq. (3) for these particles it would Fig. 3. Regimes of motion for a disk. (I.) Steady fall. (II.) Periodic oscillations. (III.) Glide–tumble. (IV.) Tumbling. Modified from Stringham et al. [70]. 6 M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 indicate the instability of their motion. To sum up, it can be stated that the stability of a path, or resistance towards tumbling, increases with increasing aspect ratio. The density ratio is known to be correlated with the steady fall or rise velocity. Moreover, as the relative velocity of the particle increases so do the inertial effects acting to destabilize the particle. As the Reynolds number increases the recirculation zone in the wake of the particle becomes unstable and the particle initiates its secondary motion characteristic for that shape. In the absence of turbulence, the wake instability is the only source acting to promote the secondary motion. 5. Orientation dependent models Whether spherical or non-spherical, regular or irregular the motion of particles is derived by considering the conservation of linear and angular momentum. In differential form the equations can be stated as: d→ x =→ up ; dt mp d→ up → = ∑ Fi dt i ð5Þ → dθ =→ ωp ; dt → → → d ωp = ∑ Ti Ip dt i ð6Þ where x is the position vector, up is the particle linear velocity, mp is the particle mass, F is the force acting on the particle, θ is the angle between the principle axis of the particle and the inertial coordinate system, ωp is the angular velocity, Ip is the moment of inertia and T is the torque acting on the particle. Where Eq. (5) deals with the location and linear velocity of the particle Eq. (6) is responsible for the orientation and the angular velocity of the particle. Eqs. (5) and (6) nicely demonstrate the similarity between translational and rotational motion. However, these equations are only strictly correct for a particle which is symmetric around the center of mass (a sphere). For a generic non-spherical particle it is necessary to include an additional cross-linked term which addresses the difference in the moment of inertia in the different directions, see Appendix A. Generally, the particle translation is evaluated in the inertial coordinate system whereas the particle rotation is evaluated in a coordinate system parallel to the principle axis of the particle. Fig. 4 shows the relationship between the different coordinate systems. Since surface forces acting in the inertial system are based on the orientation of the particle this necessitates the use of transformation between coordinates. Moreover, besides solving for the particle position as well as the particle linear and angular velocity, it is required to solve for the particle orientation represented by the angles between the co-rotational and the co-moving coordinate systems; the so-called Euler angles. The entire procedure relating to the translation and rotation of a non-spherical particle has been outlined in Appendix A. There it can be seen that it is also necessary to use an additional ordinary differential equation for particles which undergo full rotation to avoid the singularity which would otherwise occur when the Euler angles are used in relation to the co-rotational coordinate system. Despite this, it may be stated that the additional evaluation of particle rotation and orientation only require approximately twice the processing power and memory. Thus, there is no strong argument not to consider particle rotation due to computational requirements! With regard to the mathematical procedure previously published studies usually translate the particle for a sufficiently short time interval with fixed orientation after which the particle is rotated for an equal time interval [27,79]. Physically, this implies that the translational and rotational motion is decoupled. Furthermore, such a procedure also assumes that the change in linear velocity is smaller than the change in the angular velocity. Remaining issues with regard to optimization of the numerical procedure relate to the possible use of different time steps for Eqs. (5) and (6), preferably as multiples of each other, and the use of accuracy control for each time step. The use of different time steps would ensure against redundant evaluations if the change in one velocity is much smaller than the other. It should also be noted that it is possible to use a weak coupling between the Lagrangian and Eulerian phase meaning that trajectories and the continuous phase can be updated independently during the iterations. Previous investigations have not focused on optimization of the numerical procedure but rather on the formulation of the forces and torques which act on the particle. Similar to the assumption of a spherical shape in most studies involving particles, most studies involving non-spherical particles assume Stokes flow. For non-spherical particles in Stokes flow it is possible to derive the forces and torques which act on the particle on a theoretical basis similar to the procedure used for spheres to derive the BBO-equation. The usual procedure for spheres, which involves the formulation of appropriate empirical coefficients to account for the difference from Stokes flow, is also applied for non-spherical particles. However, it is also necessary to account for the offset of the center of pressure in relation to the center of geometry, see Fig. 5. The pressure distribution on the surface of a particle inclined to the flow direction no longer follows the symmetry of that particle. This gives rise to an additional lift force as well as an additional torque due to the displacement of the center of pressure. Besides this, the main complication when considering non-spherical particles is the endless variations of the shape of the particle. To combat this, most investigations include some sort of parameter variation in the formulation of forces and torques. The most popular being the ellipsoid of revolution which can be used to resemble a large array of different shapes including flake-like particles and rod-like particles. A distinctive advantage of the ellipsoid of revolution is that it has no sharp edges which in a mathematical analysis would be seen as discontinuities. The groundbreaking work on the motion of ellipsoids was made by Jeffery [40] for suspension in uniform shear flow under Stokes conditions where the formulation for the resistance force and torque is derived. This analysis has later been expanded by Brenner2 in the 1960s to arbitrary flow fields although still only under Stokes flow conditions. Following the formal notation by Gavze [29] the equation of motion for a non-spherical particle can be formulated compactly as: t F = −ℝu−ℙ⋅u̇−∫0 Tðt−τÞ⋅u̇ðτÞdτ; F=     F U ; u= M ω ð7Þ where R, P and T are respectively the steady, potential and Basset tensors. However it should be noted that the coupling of the unsteady terms with the orientation of the particle is still a remaining challenge. Fig. 4. Relationship between the inertial (x,y,z), the co-rotational (x′,y′,z′) and the comoving (x″,y″,z″) coordinate systems. 2 The work of Jeffery was extended in the 1960s by Professor Howard Brenner in a series of publications: [7–10,32] and [11]. M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 7 Fig. 5. The pressure distribution around an inclined cylinder, the location of the center of pressure, the inclination angle α and the resulting forces acting on the particle. Pressure distribution generated by steady 2D CFD simulation. Red indicates high pressure whereas blue indicates low pressure. One study tried to derive the full equated of motion for creeping flow by simplifying the problem. As such, Lawrence and Weinbaum [46,47] conducted a study on a slightly eccentric ellipsoid of revolution with major semi-axis b = a(1 + ε), in oscillatory cross-flow, where only translational motion was considered. In addition to relevant expansions of the steady state, virtual mass and Basset a new time dependent term emerged related to the eccentricity. This shows the magnitude of the awaiting challenge and suggests that BBO-equation perhaps only is an asymptotic solution for a more general formulation as the shape goes towards complete symmetry around the center of geometry. When considering non-spherical particles in Stokes flow especially the work by Fan and Ahmadi [24,25] should be accentuated. There a complete formulation of the resistance forces as well as shear induced lift can be found along with a discussion of the importance of the individual terms. Omitting the advection part of the Navier–Stokes equations, allows for the analytical formulation of the time dependent equation of motion for spheres and the steady state solution for axis-symmetric shapes. This formulation, utilizing the Stokes flow assumption, is especially useful in the paper and pulp industry to predict the flow of fibers or in the description of blood flow. However, for many practical engineering applications it is necessary to also consider the effects of higher Reynolds number flow. For example Jeffery's [40] solution is only strictly valid for zero Reynolds number and even at Re ∼ O(10− 3) it has been proved that the inertial effect is sufficient to force nonspherical particle in a different orbit than that predicted by Jeffery [42,43]. For higher Reynolds numbers, Re N 0.1, the effect of flow separation will tend to slow down and stop any rotation caused by a velocity gradient [20]. Empirical expansions of especially the steady state term have long been the backbone in investigations at higher Reynolds number flow for both spheres and non-spherical shapes alike. For non-spherical particle this is usually done by inclusion of equivalence factors, such as the sphericity. However, as previously discussed such an approach does not address the secondary motion associated with high aspect ratio shapes. In order to model the primary and secondary motions of a steady falling non-spherical particle the following forces and torques can be proposed as the minimum required to explain the observations: mp → d up → → → → = F Drag + F Lift + F Buoyancy + F Other dt → → → d ωp → → = T resist + T offset + cross terms + T Other : Ip dt ð8Þ ð9Þ For a particle falling at its terminal velocity the steady state drag force is equal in magnitude to the buoyancy force. The lift force accounts for the sideward motion and is present when the particles principle axis is inclined to the main flow direction. With a concept taken from aerodynamics this can be explained as ‘profile’ lift. The resistance torques is the angular parallel to the steady state drag. Note that these always act to dampen the rotational motion. The torques resulting from the offset of the center of pressure from the geometric center accounts for the periodic oscillations of the particle around an axis parallel to the flow direction. Other forces acting on the particle are in this case related to the unsteady behavior of the particle. These forces will act as additional resistance towards acceleration and can generally not be assumed to be negligible. The reported secondary motion of non-spherical particles in a uniform flow field at higher Reynolds number flow, as outlined previously, was suggested to be caused by the wake of the particles. The pressure distribution is not symmetric and the particle is forced away from its initial horizontal alignment. Consider a particle failing at its terminal velocity in a uniform flow field as illustrated in Fig. 6. The pressure distribution, indicated with + and −, causes the resulting forces to work at the center of pressure rather than at the center of geometry. This noncoincidence of the center of pressure and center of gravity causes the sustained oscillations. Additionally, the pressure distribution also results in a lift force, known as profile lift, which moves the particle away from its otherwise vertical path. With regard to the drag force the main advantage for an orientation dependent calculation method that this is calculated on basis of the projected area instead of that of a sphere with the same volume as the particle:
 1 → → → → → F Drag = CD ρAp j u− up j u− up : 2 ð10Þ The challenge with regard to the drag force is the proper formulation of the drag coefficient which is applicable for a large range of Reynolds numbers, shapes and orientations. It has become a common practice to procure empirical fits at a range of Reynolds number for a specific shape. Some fits also include a parametric variation of the shape e.g. the aspect ratio of a cylinder or of an ellipsoid of revolution. However, these expressions are usually based on either a fixed orientation or a freely falling particle. Thus correlations of the drag coefficient, which consider the inclination angle, are not widely available. Two approaches have been proposed to address this predicament: The work of Rosendahl [63] suggests 8 M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 Fig. 6. Illustration of resultant forces and the pressure distribution of a particle at higher Reynolds numbers (Re N 100) in uniform flow. CP is the Center of Pressure and CG is the Center of Gravity/Geometry. FB is the buoyancy force, FL is the lift force and FD is the drag force. using a ‘blending’ function between the drag coefficient for flow normal and parallel to the major axis of the particle:
 3 CD ðαÞ = CD;α = 0 + CD;α = 90 −CD;α = 0 sin α ð11Þ where α is the angle between the major axis of the particle and the flow direction. Here the projected area at the evaluated orientation is used in the calculation of the drag force. Secondly, the work by Yin et al. [79] suggests using available drag correlations expressed by the sphericity and thus solely accounting for the dependence of orientation by using the projected area in the calculation of the drag force. Recently, a third option has presented itself. Based on a plethora of empirical data for fixed and freely falling particles Hölzer and Sommerfeld [37] came up with an expression which uses a cross-wise, ψ⊥, and lengthwise sphericity, ψ||, to account for the drag coefficient of different shapes at different orientations: CD = 8 1 16 1 3 1 0:4ð− logψÞ0:2 1 pffiffiffiffi + pffiffiffiffiffi + qffiffiffiffiffiffiffiffiffiffiffi + 0:4210 Re ψ∥ Re ψ Re ψ⊥ 3=4 ψ area. The cross-wise sphericity should thus aid in the correlation of the form drag while the lengthwise sphericity is expressive of the friction drag. Note that here the Reynolds number and the drag coefficient are based on the volume equivalent sphere. Fig. 7 shows the drag force for a cylinder at different orientations, normalized with the drag force at zero incidence angle, calculated using the three suggested methods and compared to the benchmark (lattice-Boltzmann simulations) by Hölzer and Sommerfeld [36]. Overall, it may be noted that the drag force increases with increasing incidence angles due to the increase in projected area. However, this alone is not sufficient to properly account for the observed results. The method by Rosendahl [63] provides a pragmatic way to calculate the drag force at different incidence angles but also relies upon the availability of experimental data. For regular shapes these can typically be found for particles at 90° incidence angle whereas empirical fits for particle at zero incidence angle are not widely available. In this regard it might be useful to refer to the studies by Militzer et al. [54] which provide a parametric fit for spheroids as a function of the Reynolds number and the aspect ratio as well as Isaacs and Thodos [39] which provides the same for disks and cylinders at 90° incidence angle. For the present benchmark data it may be noted that a ‘blending’ function using sin(α) instead of sin3(α) provides a superior fit. Hölzer and Sommerfeld [37] constitute a good fit of the present benchmark data and attractively address all possible shapes at all Reynolds numbers in a single expression. However, this also indicates that for some specific shapes such a correlation, similar to the one by Yin et al. [79], might be associated with a relatively large error compared to correlations developed for that specific shape. The theoretical and empirical basis of predicting the profile lift relies on much more scant information compared to that available for drag. For symmetric particles the lift is zero at both α = 0° and α = 90° and it assumes a maximum somewhere in between dependent on the shape and Reynolds number. The usual assumption has been to assume that the lift is proportional to the drag and that the dependence with the orientation is given by the so-called ‘crossflow principle’ with reference to Hoerner [35]: CL 2 = sin α⋅ cos α: CD ð13Þ This relationship was developed for infinite cylinders at Reynolds number in the Newton law regime. Fig. 8 shows data for a spheroid with small aspect ratio together with the cross-flow principle from Eq. (13). It can be seen that the cross-flow principle provides a fair fit to the present data at Reynolds numbers in the Newton law regime whereas the maximum lift/drag ratio diminishes as the Reynolds number ð12Þ Here the cross-wise sphericity is the ratio between the crosssectional area of the volume equivalent sphere and the projected area of the actual particle. The lengthwise sphericity is the ratio between the cross-sectional area of the volume equivalent sphere and the difference between half of the surface area and the mean projected Fig. 7. Evaluation of the different approaches to correlate the drag coefficient with the incidence angle. 9 M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 Table 4 The different expressions used to correlate the location of the center of pressure. Rayleigh [60] Marchildon et al. [52] Rosendahl [63] Yin et al. [79] Fig. 8. Lift/drag ratio at different Reynolds numbers. Data by Hölzer and Sommerfeld [36] and Rosendahl [62]. decreases. This is related to the relative importance of the friction and pressure drag at these intermediate Reynolds numbers. Here we provide the following fit to the present data set (30 b Re b 1500) to correlate the influence of the Reynolds number for the cross-flow principle: CL sin2 α⋅ cosα = CD 0:65 + 40Re0:72 ð14Þ This expression gives the correct asymptotic values for large and small Reynolds numbers but is based on a narrow data set with following low accuracy. It should also be noted that data shown in Fig. 8 is for a spheroid with relatively low aspect ratio. It seems as if the better the shape approximates an infinitely long cylinder, the clearer the resemblance with the cross-flow principle becomes. Once the lift coefficient is specified the lift force can be found using an expression equivalent to Eq. (10). In order to correctly predict the incidence angle, to estimate the forces and torques, it is of prime importance to locate the center of pressure. As previously stated, a non-spherical particle tends to fall with its largest cross-sectional area normal to the flow direction i.e. α = 90°. Here the center of pressure is coincidental with the geometric center and the lift force and torque is zero. Hence this can be described as the state of stable equilibrium of the particle. A non-spherical particle inclined to the flow direction with α = 0° will also experience no lift or torque but this can instinctively be perceived as an unstable equilibrium. At this extreme the center of pressure must therefore be non-coincidental with the geometric center to match observed behavior. Using concepts from airfoil theory the center of pressure at this extreme inclination is placed at the “quarter chord point” which in this case is equivalent to half the distance from the geometric center to the end of the particle which is oriented towards the flow [63,79]. Please refer to Fig. 5 for visual illustration. Marchildon et al. [52] provide a linear approximation to the derivation3 by Rayleigh [60] for the pressure distribution on an infinite flat plate to predict the center of pressure of a cylinder. This is reported by Marchildon et al. [52] to be valid for inclinations above α = 15° due to the uniformity of the pressure distribution above this angle. Both Rosendahl [63] and Yin et al. [79] present expressions which close the gap with regard to 3 Derived by the application of discontinuous potential flow theory. xcp = L xcp = L xcp = L xcp = L = ð3 = 4Þðsinαi Þ = ð4 + π cosαi Þ = ð90−α
i Þ = 480  = 0:25 1− sin3 αi = 0:25 cos3 αi the location of the center of pressure between the two extremes (Table 4). Fig. 9 shows an illustration of the different expressions and it can be seen that there is some discrepancy in the prediction of the center of pressure. More unfortunately, there seem not to be any guidelines towards which expression is most appropriate to use. A freely falling non-spherical particle will spend most of the time close to α = 90° and effort should thus be directed towards finding the best fit close to this point. Assuming that Rayleigh's derivation is valid for general nonspherical particles at intermediate Reynolds numbers it seems attractive to use the simple linear fit by Marchildon et al. [52]. Once the lift and drag forces are found as well as the location of their point of attack, i.e. the center of pressure, it is a small matter of calculating the resulting torque which is due to the offset from the geometric center, Toffset.
→ → →  → T offset = xcp F Lift + F Drag + F Other : ð15Þ The torque due to resistance can be directly derived by integration of the friction, caused by rotation, over the length of the particle. For spheroids subject to the Stokes conditions solutions have been known since Jeffery [40] and have since been expanded to other shapes [19]. Relevant expansions for higher Reynolds number can be found by incorporating appropriate fits for the drag coefficient in the definition of the drag force before the integration is performed.
2 → 2 L = 2→ L=2 T resist = 2∫0 F resist dl = ∫0 CD;cyl ρ ωf −ωp l Ap dl: ð16Þ This integral can be evaluated with increasing degrees of sophistication. Note that if the particle aspect ratio is sufficiently large the angular velocity will tend to be low and an assumption of creeping flow may suffice. For the completeness of this investigation Fig. 9. Location of the center of pressure for a cylinder with length L. 10 M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 an evaluation of Eq. (16) for a cylinder rotating around its minor axis, see Fig. 10, may be performed as: 0 B
2 B 1 → 4B T resist = ρD ωf −ωp L B + B64 @ CD;cyl 10 : =1+ Re2 = 3 1 3:36 1 C C C
 !2 = 3 C C ρD ωf −ωp L A μ for ð17Þ Here the drag coefficient, suggested by White [77], is valid for the entire subcritical region for a cylinder in cross-flow. The rotation around the major axis is not considered here since it does not influence the oscillating motion which is the prime feature desired to be modeled. The unraveling of orientation dependent models up to now constitutes a description of the minimum number of forces and torques which are required for the modeling of non-spherical particles. For specific problems it may be necessary to address additional forces and torques. For general fluid flow, these would be the forces caused by pressure and velocity gradients as well as unsteady forces such as virtual mass and Basset history force. Some of these forces may be evaluated by simple expansions of the equivalent expressions derived for a sphere whereas others, such as the Basset force, are utterly hopeless to evaluate for non-spherical particles even in creeping flow. As a general guideline these forces may be accounted for by using the projected area or an equivalent diameter as is suggested in the approach by Rosendahl [62]. Clearly, order-ofmagnitude estimates may be performed for the forces acting on nonspherical particle similar to those which it is custom to perform for spheres and thus for most gas–solid flows it is justified to neglect the unsteady forces. For a freely falling cylinder in water it is not possible to neglect the unsteady forces since these particles are oscillating. As such Sørensen et al. (2007) found that the terminal velocity of a steady falling cylinder varied slightly in tune with the larger oscillations of the angular velocity. For that investigation an intricate expansion of the drag force, depending on the angular acceleration was developed to account for the unsteady forces. However, the general application of this expression in the calculation procedure presented here is not possible. For small non-spherical particles it might be necessary to model non-continuum effects. This is addressed in the study by Fan and Ahmadi [24] who introduce both an additional Brownian force and a Brownian torque in the equations of motion to supplement the fluid dynamic forces. At the same time the fluid dynamic forces are modified by introducing approximations of the translational and rotational slip factors. There, in an Eulerian– Lagrangian framework the nature of Brownian motion is modeled as a Gaussian random process. Considering the similarities between Brownian and turbulent motion such an approach also indicates possible approaches for non-spherical particles in turbulent flow. Also Fig. 10. Resistance towards rotation. note that the effect of velocity gradients has already been incorporated into the expression for rotational resistance, Eq. (16), through the vorticity of the flow field. The present methods do not account for the disturbance which initiates the periodic oscillatory motion for an initial horizontal aligned particle. However, if placed in a turbulent environment the turbulence would provide this initial disturbance. 6. Interaction with turbulence The presence of turbulence significantly affects the motion pattern of a particle. Large uncertainty exists concerning the interaction between non-spherical particles and turbulence. Suffice to say that the presence turbulence may severely alter the motion pattern of nonspherical particles and similarly, the motion of non-spherical particles may alter the properties of the turbulence. Consequently, the treatment of this subject will here rely more on a discussion of the underlying mechanics and suggestions for implementation strategies in the Eulerian–Lagrangian framework than on a critical evaluation of existing approaches which simply do not exist. Overall, we distinguish between methods which resolve the turbulent structures directly and methods which use an average description of turbulence. Similarly, it is a common procedure to distinguish between methods which only consider the effect of turbulence on the particles (one-way coupling) and methods which additionally consider the effect of the particles of on the turbulence (two-way coupling). Typically, the former approach can only be justified at sufficiently low concentrations [22]. If the turbulent structures are resolved and one-way coupling is assumed the previously described methodology can be utilized without further ado. However, the prohibitive requirements for fully resolved DNS make this option less attractive. The use of LES and LES/RANS-hybrids lessens the requirements somewhat but imposes additional uncertainties regarding influence of the sub-grid stresses on the particles. To show the flight of non-spherical particles in a turbulent flow field the most popular approach has been to imitate the turbulence by means of a predefined flow field. For isotropic turbulence Fan and Ahmadi [25] and Olson [56] used a Gaussian random field where the instantaneous velocity field is given as series of Fourier nodes with zero mean and specified standard deviation. Similarly, Fan and Ahmadi [24] modeled the turbulent boundary layer using periodic vortical flow structures at various distances from the wall while Shin and Maxey [65] used a flow field consisting of four counter rotating 2D vortices. For spheres, the application of the Eulerian–Lagrangian methodology in the context of DNS and LES has recently been demonstrated by Vreman et al. [74]. Here the interaction with turbulence formed coherent structures of particles as well as a flattening of the mean velocity profile and an increase of the streamwise turbulence intensity. Clearly, similar simulation strategy could be utilized to show the equivalent impact of/on non-spherical particles. For practical applications it is more attractive to base the description of the turbulent flow field on the Reynolds averaged equations. Here, the conventional approach for spheres has been to model the turbulence as stochastic Markow-sequences; so-called random walk models. The most popular among these is the eddylifetime model which has been adjusted using empirical constants to predict the turbulent dispersion observed in a wide variety of multiphase flows [66]. For non-spherical particles this approach has only been applied in conjunction with drag correlation for translational motion using the sphericity factor [71,2]. For orientation dependent models a pragmatic approach could be to apply the eddy-life time model only on the translational motion and neglect the effect of turbulence on the rotational motion considering the lack of empirical data available. More correct would be to apply similar assumptions for rotational motion as used for the translational motion to form an expansion for the eddy-lifetime model. The main assumption of the eddy-lifetime model can be stated as: eddy M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 Fig. 11. Fiber alignment in the inter-vortex region. properties is constant for the entire eddy-lifetime, particles are smaller than the smallest eddy and eddy properties are a Gaussian random function of the turbulent kinetic energy, k, and dissipation rate, ε. An appropriate expansion for rotational motion would be: constant vorticity during the lifetime of the eddy. This enters nicely into the current equation set through the fluid vorticity term in Eq. (16). The fluid vorticity is formed from the present constants in the eddy-lifetime model by considering the characteristic size, le, and velocity, u′, of an eddy: ωf = u′ ε = le Ak ð18Þ where, A is an empirical parameter which ranges between 0.135 and 0.56 [66]. However, it should be noted that such an approach is entirely untested and is merely a suggestion by the authors to include the effect on the rotational motion in the description of turbulence. Whereas the effect of turbulence on particles is well known the effect of the particle on the turbulence is much less so. For spheres, general observations seem to suggest that small particles attenuate the carrier phase turbulence while larger particles tend to augment it [31]. From studies of power-spectral measurements of the fluctuating velocity it has been observed that the addition of particles results in a decrease of the turbulence energy in the high wave number region [41]. This is interpreted as a result of the transfer of turbulent kinetic energy from the eddies to particles which are accelerated by the eddies [73]. The production of turbulence is most often thought of as being due to the wake of the particles and as such should be a function of the velocity difference between the particles and carrier fluid [33]. In the case of turbulence modulation by non-spherical particles the type of interaction is likely to be even more complicated than that of spheres. Presently, there is no consensus concerning the modeling of turbulence modulation for spheres [50] and mechanisms for nonspherical particles must be considered as a speculation. That said, the secondary motion which is associated with all non-spherical particles while falling at higher Reynolds number, Rep N 100, in an otherwise quiescent environment, suggests that they are capable of transferring 11 mechanical energy into turbulent kinetic energy in more modes than is the case of spherical particles. On the other hand Klett [45] showed that otherwise steady falling non-spherical particles exposed to turbulence would experience a wobbling or chaotic motion depending on their size and the magnitude of the turbulence. This suggests that the secondary motion acts to attenuate the carrier phase turbulence by extracting turbulent kinetic energy into secondary motion. As such it was also revealed that non-spherical particles were able to both enhance and attenuate turbulence depending on the shape as well as the ratio between the particle diameter and the length scale of the turbulence [51]. Similarly, by considering the momentum coupling only, the additional consideration of shape leads to the conclusion that non-spherical particles have a greater effect on the turbulence, than the volume equivalent spheres, due to the larger drag coefficient [71]. Finally, using DNS to resolve the turbulent structures in the near wall region Paschkewitz et al. [57] showed how rigid slender fibers would align in inter-vortex regions as seen in Fig. 11. The large stresses generated to oppose the vortex motion thereby acted to dissipate the eddies. Drag reductions of up 26% were calculated depending on the aspect ratio and the concentration showing that the shape alone can significantly alter the turbulence characteristics. Clearly, the interaction between particles and the turbulent structures must be affected by the alignment and shape of the particle. 7. Summary/conclusions This outline of the motion of large non-spherical particles is made not only to give an overview of the present status of this topic but also to serve as a blueprint for future implementations of orientation dependent models. The additional consideration of orientation and angular velocity gives a number of decisive advantages. Firstly, by modeling the orientation dependent forces and torques it is possible to predict the secondary motion caused by the non-spherical shape. Secondly, the modeling of non-spherical particles in the Lagrangian reference frame, without the severe restriction of creeping flow, allows for the possibility to use this methodology on a variety of engineering flows which contain large non-spherical particles. Thirdly, the solution procedure is only around twice as computational intensive compared to the present implementation in commercial codes. Finally, it is postulated that the influence of turbulence on nonspherical particles can be addressed by an appropriate expansion of the popular eddy-lifetime model. Appendix A. Equations of motion for non-spherical particles When the linear and angular motion of particles which are not symmetric around the center of mass is considered it is necessary to use both inertial and co-rotational coordinate systems and account Fig. 12. Relationship between the inertial (x,y,z), the co-rotational (x′,y′,z′) and the co-moving (x″,y″,z″) coordinate systems. N = plane(x′,y′) ∩ plane(x″,y″). 12 M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 for the relation between them by transformation of coordinates. The particle position and velocity determined from the following differential equations: → dx → = up dt mp ð19Þ → d up → → → → = F Drag + F Lift + F Buoyancy + F Other dt ð20Þ where mp and up are respectively the mass and velocity of the particle, → F is a force acting on the particle, and x is the position vector expressed in the inertial frame according to Fig. 12. Notice that the evaluation of lift and drag forces is dependent on the orientation of the particle. The resulting lift and drag forces act in the center of pressure whereas the buoyancy force acts in the center of mass which for a particle with uniform mass is coincidental with the center of geometry. However, the center of pressure is generally not coincidental with the center of geometry and thus gives rise to additional torques acting on the particle. The rotational motion uses the corotational particle frame → x ′ = ½x′; y′; z′Š with origin at the particles mass center and its axis aligned with the primary axis of the particle while the co-moving coordinate → x ″ = ½x″; y″; z″Š has its axis aligned with that of the inertial frame. The differential equations for calculating the angular velocity are given by: ð21Þ where , Ix′, Iy′, Iz′, Tx′, Ty′, Tz′, ωx′, ωy′, and ωz′ are respectively the moments of inertia, the torques acting on the particle and the particle angular velocities around their principle axes. The additional terms in the angular momentum equation vanish for particles which are symmetric around the center of mass (a sphere) but needs to be retained for non-spherical particles. The main components which make up the torque are the resistance towards rotation and the offset between the center of pressure and geometric center. Notice that it is not possible to present this set of equations in vector format due to the cross-coupling of the angular velocity. The transformation between the co-moving and the co-rotational coordinates is accomplished by means of a transformation matrix, A [30]: → x ′ = A→ x″ ð22Þ where the elements in A represent the directional cosines of the angles [θ, ϕ, ψ] between the principle axis of the co-rotational and the co-moving coordinate system. These angles are also known as the Euler angles. However, these angles are not suitable for particles which undergo full rotation due to a singularity which occurs when they are used in relation to the angular velocities of the particle. Instead Euler's four parameters [ε1, ε2, ε3, η], which are also known as quaternions, are used. The four Euler parameters represent an expansion of the three Euler angles to eliminate the singularity. The transformation matrix using the Euler parameters is given by Hughes [38]:
 2 2 1−2 ε2 + ε3 6 6 6 A = 6 2ðε2 ε1 −ε3 ηÞ 6 4 2ðε1 ε3 + ε2 ηÞ 2ðε2 ε1 + ε3 ηÞ
 1−2 ε23 + ε21 2ðε3 ε2 −ε1 ηÞ Here the Euler parameters have been related to the Euler angles by the following relations: ϕ−ψ θ ϕ−ψ θ ϕ−ψ θ sin ; ε2 = sin sin ; ε3 = sin cos ; 2 2 2 2 2 2 ϕ−ψ θ η = cos cos : ð24Þ 2 2 ε1 = cos The time rate of change of the Euler parameters, used to update the orientation of the particles, is calculated by:
 dωx′ = ∑Tx′;i + ωy′ ωz′ Iy′ −Iz′ Ix′ dt dωy′ = ∑Ty′;i + ωz′ ωx′ ðIz′ −Ix′ Þ Iy′ dt
 dωz′ = ∑Tz′;i + ωx′ ωy′ Ix′ −Iy′ Iz′ dt 2 Fig. 13. Typical algorithm to solve for the translation and rotation of a non-spherical particle. 2ðε1 ε3 −ε2 ηÞ 3 7 7 7 2ðε3 ε2 + ε1 ηÞ 7: 7
5 2 2 1−2 ε1 + ε2 ð23Þ 3 dε1 6 dt 7 7 6 3 2 6 dε 7 ηωx′ −ε1 ωy′ + ε2 ωz′ 6 27 7 6 6 1 ε3 ωx′ + ηωy′ −ε1 ωz′ 7 6 dt 7 7: 7= 6 6 6 dε3 7 2 4 −ε2 ωx′ + ε1 ωy′ + ηωz′ 5 7 6 −ε1 ωx′ −ε2 ωy′ −ε3 ωz′ 6 dt 7 7 6 5 4 dη dt 2 ð25Þ A typical procedure for solving could be stated as: Fig. 13 illustrates a conventional algorithm to solve the trajectory of a non-spherical particle where the translational and rotational motion is decoupled. Similarly, the same fixed time interval is used for both the translation and rotation of the particle. References [1] T. Allen, Particle size Measurement, Chapman & Hall, 1981. [2] R.I. Backreedy, L.M. Fletcher, J.M. Jones, L. Ma, M. Pourkashanian, A. Williams, Cofiring pulverized coal and biomass: a modeling approach, Proceedings of the Combustion Institute 30 (2005) 2955–2964. [3] F.M. Barreiros, P.J. Ferreira, M.M. Figueiredo, Calculating shape factors from particle sizing data, Particle and Particle Systems Characterization 13 (6) (1996) 368–373. [4] O. Bernstein, M. Shapiro, Direct determination of the orientation distribution function of cylindrical particles immersed in laminar and turbulent shear flows, Journal of Aerosol Science 25 (1) (1994) 113–136. [5] W.K. Bilanski, R. Lal, Behavior of threshed materials in a vertical wind tunnel, Transactions of the ASAE 8 (1965) 411–416. [6] P. Bowen, J. Sheng, N. Jongen, Particle size distribution measurement of anisotropic—particles cylinders and platelets—practical examples, Powder Technology 128 (2–3) (2002) 256–261. [7] H. Brenner, The Stokes resistance of a slightly deformed sphere, Chemical Engineering Science 19 (8) (1964) 519–539. [8] H. Brenner, The Stokes resistance of an arbitrary particle. 2. An extension, Chemical Engineering Science 19 (9) (1964) 599–629. [9] H. Brenner, The Stokes resistance of an arbitrary particle. 3. Shear fields, Chemical Engineering Science 19 (9) (1964) 631–651. [10] H. Brenner, The Stokes resistance of an arbitrary particle. 4. Arbitrary fields of flow, Chemical Engineering Science 19 (10) (1964) 703–727. [11] H. Brenner, D.W. Condiff, Transport mechanics in systems of orientable particles. 3. Arbitrary particles, Journal of Colloid and Interface Science 41 (2) (1972) 228-&. [12] A. Carlsson, Orientation of fibres in suspensions flowing over a solid surface, licentiate thesis, Royal Institute of Technology, Stockholm (2007). M. Mandø, L. Rosendahl / Powder Technology 202 (2010) 1–13 [13] R.P. Chhabra, Bubbles, Drops, and Particles in Non-Newtonian FluidsSecond edition, CNR Press, 2006. [14] R.P. Chhabra, L. Agarwal, N.K. Sinha, Drag on non-spherical particles: an evaluation of available methods, Powder Technology 101 (3) (1999) 288–295. [15] E.B. Christiansen, D.H. Barker, Effect of shape and density on free settling of particles at high Reynolds number, AICHE Journal 11 (1) (1965) 145. [16] N.N. Clark, A new scheme for particle shape characterization based on fractal harmonics and fractal dimensions, Powder Technology 51 (3) (1987) 243–249. [17] R. Clift, W.H. Gauvin, The motion of particles in turbulent gas streams, Proc. Chemeca '70 1 (1970) 14–28. [18] R. Clift, R. Grace, M.E. Weber, Bubbles, Drops, and Particles, Academic Press, New York, 1978. [19] R.G. Cox, Motion of long slender bodies in a viscous fluid. Part 2. Shear flow, Journal of Fluid Mechanics 45 (1971) 625–657. [20] E.J. Ding, C.K. Aidun, The dynamics and scaling law for particles suspended in shear flow with inertia, Journal of Fluid Mechanics 423 (2000) 317–344. [21] A. Elfasakhany, Modeling of Pulverized Wood Flames, Ph.D. thesis, Lund University 2005. [22] S. Elghobashi, On predicting particle-laden turbulent flows, Applied Scientific Research 52 (4) (1994) 309–329. [23] G. Eshel, G.J. Levy, U. Mingelgrin, M.J. Singer, Critical evaluation of the use of laser diffraction for particle-size distribution analysis, Soil Science Society of America Journal 68 (3) (2004) 736–743. [24] F.G. Fan, G. Ahmadi, Wall deposition of small ellipsoids from turbulent air flows—a Brownian dynamics simulation, Journal of Aerosol Science 31 (10) (2000) 1205–1229. [25] F.G. Fan, G. Ahmadi, Dispersion of ellipsoidal particles in an isotropic pseudoturbulent flow-field, Journal of Fluids Engineering—Transactions of the ASME 117 (1) (1995) 154–161. [26] S.B. Field, M. Klaus, M.G. Moore, F. Nori, Chaotic dynamics of falling disks, Nature 388 (6639) (1997) 252–254. [27] I. Gallily, A. Cohen, On the orderly nature of the motion of nonspherical aerosol particles. II. Inertial collision between a spherical large droplet and an axially symmetrical elongated particle, Journal of Colloid and Interface Science 68 (2) (1979) 338–356. [28] G.H. Ganser, A rational approach to drag prediction of spherical and nonspherical particles, Powder Technology 77 (2) (1993) 143–152. [29] E. Gavze, The accelerated motion of rigid bodies in non-steady stokes flow, International Journal of Multiphase Flow 16 (1) (1990) 153–166. [30] H. Goldstein, Classical Mechanics, Addison-Wesley Press, 1980. [31] R.A. Gore, C.T. Crowe, Effect of particle-size on modulating turbulent intensity, International Journal of Multiphase Flow 15 (2) (1989) 279–285. [32] J. Happel, H. Brenner, Low Reynold's Number Hydrodynamics, Prentice-Hall, Englewood Cliffs, 1965. [33] G. Hetsroni, Particles turbulence interaction, International Journal of Multiphase Flow 15 (5) (1989) 735–746. [34] H. Heywood, Uniform and nonuniform motion of particles in fluids, Proceedings of the Symposium on the Interaction between Fluids and Particles, 1962, pp. 1–8. [35] J.F. Hoerner, Fluid-dynamics drag, Hoerner Fluid Dynamics, 1965 (Published by the author). [36] A. Hölzer, M. Sommerfeld, Lattice Boltzmann simulations to determine drag, lift and torque acting on non-spherical particles, Computers & Fluids 38 (3) (2009) 572–589. [37] A. Hölzer, M. Sommerfeld, New simple correlation formula for the drag coefficient of non-spherical particles, Powder Technology 184 (3) (2008) 361–365. [38] P.C. Hughes, Spacecraft Attitude Dynamics, Wiley, 1986. [39] J.L. Isaacs, G. Thodos, Free-settling of solid cylindrical particles in turbulent regime, Canadian Journal of Chemical Engineering 45 (3) (1967) 150-&. [40] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proceedings of the Royal Society of London.Series A, Containing Papers of a Mathematical and Physical Character 102 (715) (1922) 161–179. [41] B.H. Jou, H.J. Sheen, Y.T. Lee, Particle mass loading effect on a 2-phase turbulent downward jet flow, Particle & Particle Systems Characterization 10 (4) (1993) 173–181. [42] A. Karnis, H.l. Goldsmit, S.G. Mason, Flow of suspensions through tubes. V. Inertial effects, Canadian Journal of Chemical Engineering 44 (4) (1966) 181-&. [43] A. Karnis, S.G. Mason, H.L. Goldsmith, Axial migration of particles in Poiseuille flow, Nature 200 (490) (1963) 159-&. [44] G. Kaspers, Dynamics and measurement of smokes. I Size characterization of nonspherical particles, Aerosol Science and Technology 1 (2) (1982) 187–199. [45] J.D. Klett, Orientation model for particles in turbulence, Journal of the Atmospheric Sciences 52 (12) (1995) 2276–2285. [46] C.J. Lawrence, S. Weinbaum, The unsteady force on a body at low Reynoldsnumber—the axisymmetric motion of a spheroid, Journal of Fluid Mechanics 189 (1988) 463–489. [47] C.J. Lawrence, S. Weinbaum, The force on an axisymmetrical body in linearized, time-dependent motion—a new memory term, Journal of Fluid Mechanics 171 (1986) 209–218. [48] J.P. Le Roux, Settling velocity of ellipsoidal grains as related to shape entropy, Sedimentary Geology 101 (1–2) (1996) 15–20. 13 [49] L.G. Leal, Particle motions in a viscous fluid, Annual Review of Fluid Mechanics 12 (1) (1980) 435–476. [50] M. Mando, M.F. Lightstone, L. Rosendahl, C. Yin, H. Sørensen, Turbulence modulation in dilute particle-laden flow, International Journal of Heat and Fluid Flow 30 (2009) 331–338. [51] M. Mandoe, L. Rosendahl, Measurement of turbulence modulation by nonspherical particles, Proceedings of the 7th International Conference on Multiphase Flow, ICMF, 2010. [52] E.K. Marchildon, A. Clamen, W.H. Gauvin, Drag + oscillatory motion of freely falling cylindrical particles, Canadian Journal of Chemical Engineering 42 (4) (1964) 178-&. [53] G. Mckay, W.R. Murphy, M. Hillis, Settling characteristics of disks and cylinders, Chemical Engineering Research & Design 66 (1) (1988) 107–112. [54] J. Militzer, J.M. Kan, F. Hamdullahpur, P.R. Amyotte, A.M. Al Taweel, Drag coefficient for axisymmetric flow around individual spheroidal particles, Powder Technology 57 (3) (1989) 193–195. [55] V. Olaisen, N. Nesse, H. Volden, Use of laser diffraction for particle size distribution measurements in duodenal digesta, Journal of Animal Science 79 (2001) 761–765. [56] J.A. Olson, The motion of fibres in turbulent flow, stochastic simulation of isotropic homogeneous turbulence, International Journal of Multiphase Flow 27 (12) (2001) 2083–2103. [57] J.S. Paschkewitz, Y. Dubief, C.D. Dimitropoulos, E.S.G. Shaqfeh, P. Moin, Numerical simulation of turbulent drag reduction using rigid fibres, Journal of Fluid Mechanics 518 (2004) 281–317. [58] S. Paulrud, J.E. Mattsson, C. Nilsson, Particle and handling characteristics of wood fuel powder: effects of different mills, Fuel Processing Technology 76 (1) (2002) 23–39. [59] M.N. Pons, H. Vivier, K. Belaroui, B. Bernard-Michel, F. Cordier, D. Oulhana, J.A. Dodds, Particle morphology: from visualisation to measurement, Powder Technology 103 (1) (1999) 44–57. [60] L. Rayleigh, On the resistance of fluids, Philosophical Magazine Series 5 2 (13) (1876) 430–441. [61] M. Rhodes, Introduction to Particle TechnologySecond edn, Wiley, 2008. [62] Rosendahl, L.A. 1998, Extending the modelling framework for gas-particle systems, PhD thesis, Aalborg University. [63] L. Rosendahl, Using a multi-parameter particle shape description to predict the motion of non-spherical particle shapes in swirling flow, Applied Mathematical Modelling 24 (1) (2000) 11–25. [64] J.E. Shellard, R.H. Macmillan, Aerodynamic properties of threshed wheat materials, Journal of Agricultural Engineering Research 23 (3) (1978) 273–281. [65] H. Shin, M.R. Maxey, Chaotic motion of nonspherical particles settling in a cellular flow field, Physical Review E 56 (5) (1997) 5431–5444. [66] J.S. Shirolkar, C.F.M. Coimbra, M. Queiroz McQuay, Fundamental aspects of modeling turbulent particle dispersion in dilute flows, Progress in Energy and Combustion Science 22 (4) (1996) 363–399. [67] S.J.R. Simons, Modelling of agglomerating systems: from spheres to fractals, Powder Technology 87 (1) (1996) 29–41. [68] M. Sommerfeld, B. van Wachem, R. Oliemans, Best practice guidelines—on computational multiphase dynamics for turbulent dispersed multiphase flows, Version 16-Oct-07 edn, ERCOFTAC Special Interest Group on “Dispersed Turbulent Multi-Phase Flow”, 2007. [69] H. Sørensen, L. Rosendahl, C. Yin, M. Mandoe, Settling of a cylindrical particle in a stagnant fluid, Proceedings of the 6th International Conference on Multiphase Flow, ICMF, 2007. [70] G.E. Stringham, D.B. Simons, H.P. Guy, The behaviour of large particles falling in quiescent liquids, United States Geological Survey professional paper 562-C (1969) C1–C36. [71] L. Sun, J. Lin, F. WU, Y. Chen, Effect of non-spherical particles on the fluid turbulence in a particulate pipe flow, Journal of Hydrodynamics, Ser B 16 (6) (2004) 721–729. [72] T.L. Thompson, N.N. Clark, A holistic approach to particle drag prediction, Powder Technology 67 (1) (1991) 57–66. [73] J.Y. Tu, C.A.J. Fletcher, An improved model for particulate turbulence modulation in confined 2-phase flows, International Communications in Heat and Mass Transfer 21 (6) (1994) 775–783. [74] B. Vreman, B. Geurts, N. Deen, J. Kuipers, J. Kuerten, Two- and four-way coupled Euler–Lagrangian large-eddy simulation of turbulent particle-laden channel flow, Flow Turbulence and Combustion 82 (1) (2009) 47–71. [75] H. Wadell, The coefficient of resistance as a function of Reynolds number for solids of various shapes, Journal of the Franklin Institute 217 (4) (1934) 459–490. [76] B.W. Whalley, J.D. Orford, The use of fractals and pseudofractals in the analysis of two-dimensional outlines: review and further exploration, Computers & Geosciences 15 (2) (1989) 185–197. [77] F.M. White, Viscous Fluid FlowSecond edition, McGraw-Hill, 1991. [78] W.W. Willmarth, N.E. Hawk, R.L. Harvey, Steady and unsteady motions and wakes of freely falling disks, Physics of Fluids 7 (2) (1964) 197–208. [79] C.G. Yin, L. Rosendahl, S.K. Kaer, H. Sørensen, Modelling the motion of cylindrical particles in a nonuniform flow, Chemical Engineering Science 58 (15) (2003) 3489–3498.