Numerical Study of the Orientation of Cylindrical Particles in a

Nov 28, 2016 - Cylindrical particles become more oriented along rising height position in the riser, increasing the slenderness ratio and decreasing t...
0 downloads 0 Views 3MB Size
Subscriber access provided by Warwick University Library

Article

A Numerical Study of the Orientation of Cylindrical Particles in a Circulating Fluidized Bed Jie Cai, Charley Wu, Behdad Moghtaderi, Elham Doroodchi, Zhengbiao Peng, Xiaobao Zhao, and Zhu Lin Yuan Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b04022 • Publication Date (Web): 28 Nov 2016 Downloaded from http://pubs.acs.org on November 29, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

A Numerical Study of the Orientation of Cylindrical Particles in a Circulating Fluidized Bed Jie Cai1,2,3*, Charley Wu3*, Behdad Moghtaderi4, Elham Doroodchi4, Zhengbiao Peng4, Xiaobao Zhao2, Zhulin Yuan1, 1 School of Energy and Environment, Southeast University, Nanjing 210096, China 2 School of Energy and Mechanical Engineering, Nanjing Normal University, Nanjing 210042, China 3 Department of Chemical and Process Engineering, University of Surrey, Guildford GU2 7XH, UK 4 Discipline of Chemical Engineering, School of Engineering, The University of Newcastle, University Drive, Callaghan, NSW 2308, Australia

ABSTRACT:

Orientation

is

pivotal

in

circulating

fluidized

beds

using

cylinder-shaped particles. This study developed a three-dimensional multi-way coupling model for predicting the orientation of cylindrical particles in circulating fluidized beds. The coupling algorithm between motorial cylindrical particles and the turbulence was established by incorporating the correlation between Lagrangian time scales and the k − ε model. Collisions among cylindrical particles were solved by combining rigid body impact dynamics with the hard sphere model. The results showed that the majority of cylindrical particles are prone to align with the streamline, which is in good agreement with the experimental observation. Cylindrical particles become more oriented along rising height position in the riser, increasing slenderness ratio and decreasing distance to wall. Under turbulence conditions, the effect of Reynolds number on the orientation of cylindrical particles was found to be marginal. KEYWORDS: orientation, cylindrical particle, gas-solid flow, multi-way coupling, collisions among cylindrical particles

1. INTRODUCTION Circulating fluidized beds (CFD) involving cylindrical particles are often 1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 37

encountered in industrial processes such as direct combustion of biomass straws, molding and drying of cylindrical pills, and molding and machining of composite materials of short fibres. Orientation of cylindrical particles in these processes is of utmost importance as it directly determines the distribution of solids and hence the hydrodynamics of the gas-solid flow which in turn dictates the overall performance of the system

1-2

. However, a good understanding of the theory underlying the

gas-cylindrical particle flow is largely lacking due to the complex geometry of cylindrical particles, in particular in the presence of turbulence. Experimental study of the gas-cylindrical particle flow behavior has been proven difficult due to technological limitations regarding collection of data associated with the orientation of individual particles

3-4

. These limitations have greater impact in

situations where the turbulence is present and significantly affects the flow of cylindrical particles. As yet, no method of reliable measurement of the orientation of cylindrical particles has been published. Thus, development of accurate and efficient theoretical models to predict the cylindrical particle orientation under turbulence conditions becomes critical. Two approaches were generally employed in the literature for modeling the motion of cylindrical particles in the gas-solid flow: the slender-body theory the Discrete Element Method (DEM)

7-8

5-6

and

. However, the former was reported to have

difficulties in coping with the two-way coupling between cylindrical particles and turbulence as it cannot incorporate the correlation between Lagrangian time scales and the k − ε model

9-10

. Despite the dilute gas-solid flow in CFB, it was identified in

previous studies that the presence of cylindrical particles could significantly affect the fluid flow11. Therefore, failure in accurately solving the coupling interactions between the solid phase and the turbulence greatly limits the application of the slender-body theory. In combination with computational fluid dynamics (CFD), DEM firstly attempted by Tsuji

12

for modeling gas-solid flow, becomes increasingly popular to

model the flow of cylindrical particles. In DEM, a cylindrical particle is assembled by a string of discrete elements along its axis. The force acting on the cylindrical particle 2

ACS Paragon Plus Environment

Page 3 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

is the vector sum of the forces acting on each of the discrete elements. Using DEM, coupling interactions between the discrete elements of a cylindrical particle and turbulence could be readily defined based on the correlation between Lagrangian time scales and the k − ε model. Delaney et al.

13

generalized a particle breakage model

based on DEM allowing non-round particles to be broken into non-round progeny so as to predict the particle flow and compression breakage of non-round rock when passing through an industrial scale cone crusher. Cleary et al.

14

provided an efficient

collision detection model for three dimensional super-ellipsoidal particles based on DEM. Kodam, et al.

15

proposed a model of cylindrical object contact detection for

use in discrete element method simulations. Wachs et al.

16

suggested a novel variant

of DEM to simulate the flow dynamics of granular material made of non-spherical particles. The contact detection strategy relies on the use of the Gilbert–Johnson– Keerthi algorithm to compute the distance between two convex bodies. Boon et al.

17

proposed a new contact detection algorithm between three-dimensional potential non-spherical particles in DEM using a second degree polynomial function. Euler dynamics or rigid body dynamics proposed by Euler was firstly applied to theoretically solve the translation and rotation of a three-dimensional celestial rigid body in aerospace.18 Two coordinate systems namely the body axes and the fixed-reference frame are used in Euler dynamics. The celestial body is fixed on the body axes, and the orientation of the celestial body is described by relative orientation between the body axes and the fixed-reference frame. Euler dynamics uniquely locates the orientation of a celestial body by 3 Euler angles, i.e., precession angle (ψ ), nutation angle ( θ ) and spin angle ( φ ). A cylindrical particle is a typical celestial body and has 6 degrees of freedom involving 3 orientation vectors in addition to 3 position vectors. If a cylindrical particle is fixed in the body axes with its center locating at the origin and its axis in parallel with the ζ -axis of the body axes, the angle between the ζ -axis of the body axes and the z -axis of the fixed-reference frame is defined as the 3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

nutation angle. Fig. 1 illustrates the relative orientation between the body axes and the fixed-reference frame. As an important branch of Euler dynamics, rigid body impact dynamics provides a quasi-analytical solution of the collisions between cylindrical particles and between a cylindrical particle and walls. Yin et al.19 proposed a one-way coupling model for tracking non-spherical particles in a non-uniform fluid flow field and the rotation of a cylindrical particle was solved by rigid body dynamics. Grof et al. 20

simulated the breakage of needle-shaped particles under uniaxial compaction using

a one-way coupling gas-cylindrical particle flow based on rigid body dynamics. Feng et al.

21

presented a one-way coupling model of non-spherical particle flow based on

rigid body dynamics to predict the transport of non-spherical particles in complex internal shear flows. Efforts were also made on the numerical modelling of the gas-cylindrical particle flow under turbulentce conditions. Mortensen et al. 22 used a one-way simplified rigid body dynamics model to study the translation and rotation of small prolate ellipsoidal particles in a turbulent channel flow. Njobuenwu et al.

23

proposed an one-way

coupling Eulerian-Lagrangian approach based on large eddy simulation with a dynamic sub-grid scale model to predict a fully developed gas-solid flow at a shear Reynolds number Re = 300. In this model, the translation and rotation of non-spherical ellipsoidal particles were analyzed using the rigid body dynamics. However, a theoretical model that accurately calculates the collision dynamics of cylindrical particles and rigorously accounts for the coupling between the motion of cylindrical particles and turbulence is largely missing.

4

ACS Paragon Plus Environment

Page 4 of 37

Page 5 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 1. Nutation angle θ between body axes and the fixed-reference frame In this study, a three-dimensional multi-way coupling model has been developed to predict the orientation of cylindrical particles in a circulating fluidized bed. The two-way coupling algorithm between the motion of cylindrical particles and the turbulence was established by incorporating the correlation between Lagrangian time scales and the k − ε model. Collisions of cylindrical particles were solved using the rigid body impact dynamics based on the framework of hard sphere model, i.e., the modified Nabu collision probability method. The established model was validated using the experimental data obtained in this study. Detailed information on the gas-cylindrical particle flow characteristics has been obtained. Subsequently, the orientations of cylindrical particles in a circulating fluidized bed were investigated. Specifically, the axial and radial orientations of cylindrical particles in the riser were demonstrated and discussed. Further, effects of particle slenderness ratio and Reynolds number on the orientation of cylindrical particles were analyzed and explored.

5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 37

2. MATHEMATICAL MODELS 2.1. Modelling of turbulent flow The model of turbulent flow consists of continuity equation, time-averaged N-S equations and k − ε turbulence enclosure model. The key of the two-way coupling between the gas phase and the solid phase is the coupling of solids motion into the gas flow, i.e. adding source terms into the time-averaged N-S equations and

k − ε equations. The continuity equation is written as: ∂ ( ϕ g ρ g ) + ∇ ⋅ (ϕ g ρ g v g ) = 0 ∂t

(1)

n

where ϕ g is the void fraction of fluid, ϕg = 1 − ∑Vi / ∆V , Vi is the volume of i =1

discrete element in the fluid mesh, ∆V is the volume of fluid mesh; ρ g is the air density; vg is the velocity of fluid. The time-averaged N-S equation is expressed as: −1 ∂ 2   ϕg ρg vg ) + ∇ ⋅ (ϕg ρg vg vg ) = −ϕg ∇p + ϕg ( µ + µt ) ( ∇vg ) + ( ∇vg ) − δ ij ∇ ⋅ vg  + ϕg ρg g + fsg ( ∂t 3  

(2) where p is the pressure of fluid; µ is the dynamic viscosity; µt is the turbulent viscosity of fluid, µt = ρg Cµ

k2

ε

, C µ = 0.09 is an empirical constant; δ ij

is the Kronecker function; f sg is the source term of reaction force applied by discrete elements. Subscripts g and s represent the gas phase and the solid phase, respectively.

k − ε equations are given as:

6

ACS Paragon Plus Environment

Page 7 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

∂  µt  ∇k  + ϕgGk − ϕg ρgε + ϕg ρ g I k  (ϕ g ρ g k ) + ∇ ⋅ (ϕ g ρ g v g k ) = ∇ ⋅  ϕ g  ∂t  σk    ∂ ϕ ρ ε + ∇ ⋅ ϕ ρ v ε = ∇ ⋅  ϕ µ t ∇ε  + ϕ ε C G − C ρ ε + ϕ ρ I ( g g g )  g σ  g k ( 1ε k 2ε g ) g g ε  ∂t ( g g ) ε   

(3)

In Eq. (3), I k and I ε denote the effects of solid phase on turbulent energy and dissipation rate, respectively; Gk is the generation term of turbulent energy; σ k =1 ,

σ ε =1.33 , C1ε =1.44 , C2ε =1.92 are the empirical constants. Then we need to establish models of f sg , I k and I ε . The model, proposed by Crowe24, is used to calculate f sg : f sg = K sg (vs − vg )

(4)

In the equation, K sg is the momentum exchange coefficient, v s is the average velocity of all discrete elements in a mesh. K sg is an empirical correlation, proposed by Wen and Yu et al 25, and is given as:

 (1-ϕg )ϕg ρ g vs − vg -2.65 3 ϕg >0.8 K sg = CDs ϕg 4 ds   2 1-ϕg ) µ ρ g (1-ϕg ) vs − vg (  +1.75 ϕg ≤ 0.8 K sg =150 2 ds ϕg d s 

(5)

Where ds is the average equivalent diameter of discrete elements in a fluid mesh. The model of C Ds used in this paper was proposed by Fan et al. 26

ρ 24 C Ds cos φ " = 0.006983 + 0.6224 Res-1.046 ) ×  s ( ρ Res  g

Res =

  

−1.537

Ar 0.8524

(6)

ρ g vs − v g d s

(7)

µ

In this model, the relative intensity between free flow and forced flow is taken into

account,

so

the

modified

Archimedes

number

7

ACS Paragon Plus Environment

Ar

is

used,

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 37

Ar = d s3 ( ρs − ρg ) g / µ 2 ; φ "(α , β , γ ) is the angle between cylindrical particles and 2

the plane normal to the direction of velocity. The relation between φ "(α , β , γ ) and Euler angles is:

π

π

π

cos( − α ) = sin φ sin θ , cos( − β ) = − cosψ sin θ , cos( − γ ) = cos θ 2 2 2

(8)

The correlation between Lagrangian time scales τ t ,sg and the k − ε model was employed to establish the coupling relation between solid phase and the k − ε equations.

τ t is the characteristic time of energy-loaded turbulent vortex and is defined as:27 3 2

τ t = Cµ

k

(9)

ε

And the length scale of turbulence vortex Lt is defined as: 3

3 k2 Lt = Cµ 2 ε

(10)

And then Lagrangian time scales τ t ,sg is defined as:

τ t ,sg =

τt  vs − v g τ t   1 + Cβ    Lt  

2

Cβ = 1.8-1.35cos 2ϑ

(11)

(12)

ϑ is the angle between average velocity of discrete elements and average relative velocity. At the same time, the time scale and the length scale, representing movement, are used to estimate the turbulent energy of solid phase. The relaxation time of characteristic particles linking to inertia effect acting on discrete elements τ F ,sg is defined as:

8

ACS Paragon Plus Environment

Page 9 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

 ρs + CV ρ  g

τ F ,sg = (1-ϕg ) ρs K sg−1 

  

(13)

where CV =0.5 is the additional mass coefficient.

η sg is defined as the ratio between τ t ,sg and τ F ,sg :

ηsg =

τ t ,sg τ F ,sg

(14)

The turbulence scale of solid phase can be indicated by ksg , the covariance of velocities of solid phase and gas phase. The model was presented by Simonin.27  b + ηsg ksg = 2k   1+η sg 

  

(15)

where, ρ b = ( 1 + CV )  s + CV ρ  g

  

−1

(16)

Finally, the relative source term in k equation I k is obtained by:

Ik =

K sg

ϕg ρ g

(k

sg

− 2k )

(17)

And the source term in ε equation I ε , modelled by Elgobashi 28, is calculated by: I ε = C3ε

ε k

Ik

(18)

where C3ε =1.2 is the empirical constant. By inputting ϕ g , v s , vg , ψ , θ and φ into Eqs. (4), (5), (6) and (7), f sg can be achieved; together with inputting k and ε into Eqs. (17) and (18), I k and I ε can be calculated; and then by calculating time-averaged N-S equations and

k − ε equations, the turbulent flow coupled by solid phase is obtained.

9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 37

2.2. Modeling of the discrete phase (cylindrical particles) Discrete phase analysis consists of two strands: forces and the motion model of cylindrical particles in the fluid, and collisions among cylindrical particles. 2.2.1 Forces and the motion model The force analysis of a cylindrical particle in the fluid is shown in Fig.2.

Figure 2. Schematic diagram of force of a cylindrical particle in the fluid Mathematically, the force acting on a cylindrical particle F , of which length is l , is the integral of the forces acting on any part along its length F ( x ) . l  2 F =  ∫− 2l F ( x ) dx   l  M = 2 r x × F x dx ∫− 2l ( ) ( )  

(19)

where x is the position along the length of the cylindrical particle. After mathematical transformation, Eq. (19) is changed into: l 1   F = 2 ∫−1 F ( s ) ds  2  M = l 1 r ( s ) × F ( s ) ds  4 ∫−1

(20)

where F ( s ) is the force acting on unit length in the interval [-1, 1]; M is the torque acting on the centroid of a cylindrical particle in the fixed reference frame. 10

ACS Paragon Plus Environment

Page 11 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Then the integral formula should be transformed into relative form of discrete mathematics to exert the force analysis. Eq. (20) is discretized using Legendre-Gaussian integral method. The amount of discrete elements is 10. Then the relevant discrete form is given as: l 10   F = 2 ∑ λi Fi ( si ) / li  i =1  2 10 l M = ∑ r ( si ) × λi Fi ( si ) / li  4 i =1

(21)

Here Fi ( si ) is the component of force acting on discrete element i , Fi ( si ) = K sg v g − v ( si )  V ( si ) , V ( si ) is the volume of discrete element i ; r ( si ) is

the position vector between discrete element i and the centroid of the cylindrical particle in the fixed-reference frame; λi is the weight coefficient. Then the discrete elements are cylinders with different li but the slenderness ratios of discrete elements are small. The biggest and the smallest sizes of discrete elements are 7.15 mm and 1.31 mm, and the size of gird is 100 mm. The sizes of discrete elements are smaller enough than the size of fluid mesh. v ( si ) is the velocity of discrete element i: v ( si ) = u0 + ( A ⋅ ω * ) × r ( si )

(22)

where u0 is the velocity of a cylindrical particle, i.e. the centroid velocity; A is the cosine matrix of conversion of coordinates between two reference frames; ω* is the angular velocity in the body axes. In this paper, if a symbol has *, it means the symbol is in the body axes. Because the force and translation analysis of a cylindrical particle is in the fixed-reference frame and the rotation analysis of a cylindrical particle is in the body axes, the conversion of coordinates between two reference frames is necessary with matrix A. The relationship between the matrix A and the Euler angles can be expressed as, 18 11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

cosψ cos φ − sinψ cos θ sin φ A = sinψ cos φ + cosψ cos θ sin φ  sin θ sin φ

− cosψ sin φ − sinψ cos θ cos φ − sinψ sin φ + cosψ cos θ cos φ sin θ cos φ

Page 12 of 37

sinψ sin θ  − cosψ sin θ   cos θ

(23) The translation is tracked in the fixed-reference frame by: d( mu0 ) = F + mg dt

(24)

where m is the mass of a cylindrical particle, t is the time, and g is the acceleration of gravity. The rotation of a cylindrical particle is tracked in the body axes using Euler dynamics. Firstly, the torque in the fixed-reference frame should be transformed into the torque in the body axes: M * = A −1 M

(25)

Then the rotation equation of a cylindrical particle in the body axes is expressed as:

J* ⋅

dω* + ω* × J * ⋅ ω* = M * dt

(26)

2.2.2 Collision dynamics of cylindrical particles For dense or semi-dense gas-solid two-phase flow, collision events between cylindrical particles must be considered. Considering the dilute pneumatic conveying flow regime in the CFB, the hard sphere model was employed in this work. As the hard model was established initially for the analysis of collisions between spherical particles, modifications need to be made to make it work for non-spherical particles. In this paper, each discrete element of cylindrical particles is treated as an individual particle and contributes to the collision dynamics between cylindrical particles. Specifically, the working principle is presented as follows: (i) all discrete elements in a grid are sorted; (ii) in the grid, all discrete elements will be judged as whether they will collide with other discrete elements except those belonging to its own cylindrical particle; and (iii) if two discrete elements will collide with each other, the two 12

ACS Paragon Plus Environment

Page 13 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

cylindrical particles that the two discrete elements belong to will collide with each other and the two discrete elements are regarded as the corresponding collision points of the two cylindrical particles. The modified Nanbu collision probability is given as: N

N

j −1

j −1

Pi = ∑ Pij = ∑

n π Ds2Gij dt N

(27)

where n is the real number of discrete elements; N is number of sampled discrete elements; Ds is the equivalent volume diameter of discrete element i; and Gij is relative velocity between discrete element i and j. A random number R is extracted from a generator which has a uniform distribution ranging from zero to unity. If

0 < Pi < 1 and R >

j − Prij N

(28)

Then the collision between discrete element i and j is considered to occur. Collision dynamics between two cylindrical particles was calculated according to the model proposed by Smith 29. There are two arbitrarily-shaped rigid-bodies. The masses are mi and m j . The rotary inertias in the body axes are J i* and J *j , and the rotary inertias in the fixed-reference frame are J i and J j . The pre-collision velocities are Vi 0 and V j 0 , and the post-collision velocities are Vi and V j . The pre-collision angular velocities in the body axes are ωi*0 and ω*j 0 , and the post-collision angular velocities in the body axes are ωi* and ω*j . The relative vector positions between collision points and the centroids are ri and r j . The pre-collision velocities of collision points are ui and u j . The cosine matrix of coordinates conversion are Ai and A j . The pre-collision velocities of collision points ui and u j are given as:

13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

 ui = Vi 0 + ( Ai ⋅ ωi*0 ) × ri = Vi 0 − ri × Ai ⋅ ωi*0  * *  u j = V j 0 + ( A j ⋅ ω j 0 ) × rj = V j 0 − r j × A j ⋅ ω j 0

Page 14 of 37

(29)

V is defined as the velocity of collision point when the collision point is in the biggest deformation. Then ui − V and u j − V are vector decomposed into two parts: (ui − V )|| and (u j − V )|| that are parallel to ri and r j ; (ui − V ) ⊥ and ( u j − V ) ⊥

that are perpendicular to ri and r j . Then (ui − V )|| and ( u j − V )|| are given as: ( ui − V )|| = ri 0 [ri 0 ⋅ ( ui − V )] = Pi ⋅ (ui − V )  ( u j − V )|| = r j 0 [ r j 0 ⋅ (u j − V )] = Pj ⋅ ( u j − V )

(30)

where ri 0 and r j 0 are the unit vector of ri and r j ; Pi and Pj are the dyadic tensor of Vi 0 and V j 0 , respectively. The impulsive moments of (ui − V )|| and ( u j − V )|| are zero because the impulses of compression pass through the centroids of

cylindrical particles. And (ui − V ) ⊥ and ( u j − V ) ⊥ are calculated by: ( ui − V ) ⊥ = ( ui − V ) − ( ui − V )P = (Ι − Pi ) ⋅ ( ui − V )  ( u j − V ) ⊥ = ( u j − V ) − ( u j − V )P = (Ι − Pj ) ⋅ ( u j − V )

(31)

where I is the unit tensor in the fixed reference frame. Assuming the velocity variations of centroids of cylindrical particles on the directions that are perpendicular to ri 0 and r j 0 are Vi ⊥ and V j ⊥ , and the relative impulsive moments are qi ⊥ and q j ⊥ ,

 qi ⊥ = miVi ⊥   q j ⊥ = m jV j ⊥

(32)

Relative angular momentum variations of the impulsive moments are: * −1 *  Ai ⋅ J i ⋅ Ai ⋅ ωi ⊥ = ri × miVi ⊥  * −1 *  Aj ⋅ J j ⋅ A j ⋅ ω j ⊥ = rj × m jV j ⊥

(33)

Both sides of Eq. (33) plus Ai ⋅ J i* ⋅ Ai−1 = J i−1 or A j ⋅ J *j ⋅ A−j 1 = J −j 1 to transform the body axes to the fixed reference frame, then Eq. (33) is converted to 14

ACS Paragon Plus Environment

Page 15 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

ωi*⊥ = J i−1 ⋅ (ri × miVi ⊥ )  * −1 ω j ⊥ = J j ⋅ (rj × m jV j ⊥ )

(34)

And the variations of collision point velocities resulted in by variations of angular velocities are given as:

ωi*⊥ × ri = [ J i−1 ⋅ (ri × miVi ⊥ )] × ri = −ri × J i−1 ⋅ (ri × miVi ⊥ )  * −1 −1 ω j ⊥ × rj = [ J j ⋅ (rj × m jV j ⊥ )] × r j = −r j × J j ⋅ (rj × m jV j ⊥ )

(35)

So the total velocity variations on the directions that are perpendicular to ri 0 or r j 0 are: ( ui − V ) ⊥ = Vi ⊥ + ωi*⊥ × ri = Vi ⊥ − ri × J i−1 ⋅ ( ri × miVi ⊥ )  * −1 ( u j − V ) ⊥ = V j ⊥ + ω j ⊥ × r j = V j ⊥ − r j × J j ⋅ ( r j × m jV j ⊥ )

(36)

The equations are the vector equations about Vi⊥ and V j ⊥ . Incorporating with Eq. (30), restitution coefficient e , e = 0.5 , momentum conservation and cancelling

V , then the variations of velocities and angular velocities will be obtained. After vector transfer, Eq. (36) is converted as: ( ui − V ) ⊥ = Vi ⊥ − ri × J i−1 ⋅ ( ri × miVi ⊥ ) = Vi ⊥ − ( ri × J i−1 × r ) ⋅ miVi ⊥  i  −1 −1 ( u − V ) = V − r × J ⋅ ( r × m V ) = V − r × J × ( j j rj ) ⋅ m jV j ⊥  j ⊥ j⊥ j j j j j⊥ j⊥

(37)

Because Ι ⋅Vi ⊥ = Vi ⊥ , Ι ⋅V j ⊥ = V j ⊥ , then compare Eq. (31) with (37),

(Ι − Pi ) ⋅ (ui − V ) = Ι ⋅Vi ⊥ − Ti ⋅Vi ⊥  (Ι − Pj ) ⋅ (u j − V ) = Ι ⋅V j ⊥ − T j ⋅V j ⊥

(38)

In Eq. (38), Ti = mi ri × J i−1 × ri , T j = m j r j × J −j 1 × rj

(39)

Then Eq. (38) is converted as:

(Ι − Pi ) ⋅ (ui − V ) = ( Ι − Ti ) ⋅Vi ⊥  (Ι − Pj ) ⋅ (u j − V ) = ( Ι − T j ) ⋅V j ⊥

(40)

Vi ⊥ = ( Ι − Ti )−1 ( Ι − Pi ) ⋅ ( ui − V )   −1 V j ⊥ = ( Ι − T j ) ( Ι − Pj ) ⋅ ( u j − V )

(41)

After vector transfer,

15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 37

The total impulsive moment is the vector sum of normal and tangential impulsive moments, then

qi = mi [Vi ⊥ + (ui − V )|| ] = mi [(Ι − Ti )−1 ⋅ (Ι − Pi ) + Pi ](ui − V )  −1 q j = m j [V j ⊥ + (u j − V )|| ] = m j [(Ι − T j ) ⋅ (Ι − Pj ) + Pj ](u j − V )

(42)

Based on momentum conservation, qi + q j = 0 , then

mi [(Ι − Ti ) −1 ⋅ (Ι − Pi ) + Pi ](ui − V ) + m j [(Ι − T j )−1 ⋅ (Ι − Pj ) + Pj ](u j − V ) = 0

(43)

In order to analyze more conveniently, we define  M = (Ι − Ti ) −1 ⋅ (Ι − Pi ) + Pi  −1  N = (Ι − T j ) ⋅ (Ι − Pj ) + Pj

(44)

Then, Eq. (43) is transformed as:

mi M ⋅ (ui − V ) + m j N ⋅ (u j − V ) = 0

(45)

Therefore get the expression of V :

V = (mi M + m j N )−1 ⋅ (mi M ⋅ ui + m j N ⋅ u j )

(46)

Substitutes Eqs. (44) and (46) into Eq. (42), then V is cancelled. qi = mi M ⋅ [ui − (mi M + m j N ) −1 ⋅ ( mi M ⋅ ui + m j N ⋅ u j )]  −1 q j = m j N ⋅ [u j − ( mi M + m j N ) ⋅ ( mi M ⋅ ui + m j N ⋅ u j )]

(47)

The increments of velocities of cylindrical particles are: 1   ∆Vi = − m (1 + e)qi  i  1  ∆V j = (1 + e)q j  mj

(48)

And the increments of angular velocities of cylindrical particles are:

 Ai ⋅ J i ⋅ ∆ωi* = ri × [−(1 + e)qi ]  *  A j ⋅ J j ⋅ ∆ω j = r j × [(1 + e)q j ]

(49)

 ∆ω* = − (1 + e ) ( J * ) −1 ⋅ A −1 ⋅ ( r × q ) i i i i  i  − 1 ∆ω*j = (1 + e ) ( J *j ) ⋅ A j −1 ⋅ ( rj × q j ) 

(50)

With vector transfer,

16

ACS Paragon Plus Environment

Page 17 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

The post-collision velocities Vi and V j , and the post-collision angular velocities ωi* and ω*j are then given as:

1  Vi = Vi 0 − m (1 + e ) qi i  ω* = ω* − (1 + e ) J * −1 ⋅ A −1 ⋅ ( r × q ) ( i) i i i i0  i  V j = V j 0 + 1 (1 + e ) q j  mj  −1 ω*j = ω*j 0 + (1 + e ) ( J *j ) ⋅ A j −1 ⋅ ( rj × q j ) 

(51)

Ti and T j are calculated by inputting mi and m j , ri and r j , J i* and J *j , Ai and A j into Eq. (39), and then M and

N are obtained by Eq. (44), and then

qi and q j are calculated with Eq. (47), finally Vi and V j , ωi* and ω*j are achieved by Eqs. (51). 2.2.3 Boundary treatment Boundary treatment is also based on the rigid impact dynamics. As the mass of wall is infinite, the velocity and angular velocity of wall is set to be zero. 2.2.4 Numerical strategy and methodologies The multi-way coupling model of gas-cylindrical particle flow was programmed with C++. The turbulence model was solved by SIMPLE, and Euler dynamics equations are solved with the 4th order Runge-Kutta algorithm. The interaction between the gas phase and cylindrical particle phase was modeled by: (i) first solving the gas-phase flow (at that time ϕ g = 1 , f sg = 0 , I k =0 and Iε =0 ); (ii) Calculating the movement and collisions of cylindrical particles using equations (21)-(51) for 10 time steps for the single-phase flow; and (iii) calculating ϕ g , f sg , I k and I ε using equations (1)-(18); (iv) Adding ϕ g , f sg , I k and I ε into equations (1)-(3), and calculating the turbulent flow models to update turbulence; (v) calculating the movement and collisions of cylindrical particles for another 10 time steps based on 17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 37

the updated turbulence; and (vi) repeating the previous steps (iii) – (v). As such, the turbulence and solid-phase fields are solved alternately in a two-way coupling manner. The program of gas-cylindrical particle flow was run on the Dell precision T5610 workstation (Intel® Xeon® Processor E5-2637 v2: 4core HT, 3.5GHz Turbo, 15 MB, RAM 128 G). To complete a typical case with a physical gas-solid flow time of 5 s, the simulation took more than 10 days.

2.3. Model parameters The shape of the riser is a cuboid with a size of 0.5m×0.5m×3m. The density and diameter of the cylindrical particles is 600 kg/m3 and 0.004 m, respectively. The sampled cylindrical particles were fed and uniformly distributed in the base of bed. The 3 Euler angles are generated randomly within the range of 0 to π, and the initial velocity of cylindrical particles is zero. Other simulation conditions and parameters are given in Table 1 and Table 2 respectively.

Table 1. Conditions of simulation

ρg

υg

Vin

uini

* ωini

ρs 3

∆t

H riser

Wriser

(kg/m3)

(m2/s)

(m/s)

(m/s)

(rad/s)

(kg/m )

(s)

(m)

(m)

1.205

1.502×10-5

5

0

0~π

600

0.001

3

0.5

In Table 2, 0-0.5m, 0.5-1.0m, 1.0-1.5m, 1.5-2.0m, 2.0-2.5m and 2.0-3.0m represent height regions of A1, A2, A3, A4, A5 and A6, and 0.25-0.3m, 0.3-0.35m, 0.35-0.4m, 0.4-0.45m and 0.45-0.5m represent radial regions of R1, R2, R3, R4 and R5 shown in Fig. 1. RT is the slenderness ratio of cylindrical particle.

Table 2. Parameters of simulation

θ (°)

Ai (m)

Ri (m)

RT

0-15

0~0.5

0.25-0.3

6

15-30

0.5~1.

0.3-0.35

8

30-45

1. ~1.5

0.35-0.4

10

18

ACS Paragon Plus Environment

Page 19 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

45-60

1.5~2.

0.4-0.45

60-75

2. ~2.5

0.45-0.5

75-90

2.5~3.

12

3. EXPERIMENT FOR MODEL VERIFICATION 3.1. Experimental apparatus The experiment riser is made of high tensile strength organic glass, of which shape is cuboid and size is 0.5m×0.5m×3m, as shown in Fig.3. Given the complexity of CFB systems, the top of the experimental riser was left open. A number of cylindrical particles were placed at the bottom of the bed. The cylindrical particles are mono-sized biomass matchsticks with the density of 600 kg/m3 and the diameter of 0.004 m. With the aid of a high-speed camera and the CAD software, the orientation of cylindrical particles in the gas-solid flow was quantified.

3.2. Analysis of experimental data Some researchers experimentally investigated the orientation of fiber suspensions and obtained valuable data for conditions of lower velocities and low particle concentrations using a high-speed camera. The translational velocities, angular velocities during fluidization and the amount of cylindrical particles are all very high. Therefore, it is difficult to obtain the reliable information from the experiment photos. In this study the number concentration and velocity distribution of cylindrical particles was symmetrical in that the cross section of the riser was quadrate. The nutation angle of a cylindrical particle of which the axis parallels the wall could be measured. The cylindrical particles of which the axes approximately parallel the wall were trapped so as to acquire more data. An experimental photo was imported into AutoCAD and all cylindrical particles in the photo were measured with “dimlinear” widget of the CAD software. If its measured length is 90% - 100% of its actual length, the cylindrical particle would be considered to be in parallel with the 19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

shooting wall, and its nutation angle would be measured with the “dimangular” widget. Then the amount of cylindrical particles of which nutation angles within a certain nutation angle range was calculated. For high veracities, a large number of experiment photos have been handled to obtain the statistically averaged data.

1- Rotary blower; 2-Stop Valve; 3-Manometer; 4-Flow equalizer; 5-Wind room; 6-Riser; 7- High-speed camera; 8-News lamp; 9-A/D transfer; 10-Differential pressure transducer; 11-Datd acquisition computer, 12-Flowmeter

Figure 3. Schematic diagram of visual experimental system

4. RESULTS AND DISCUSSION 4.1. Model validation The turbulent gas-cylindrical particle flow is different from single-phase turbulent flow. The reasonable and accurate turbulence is one of the prerequisites to precisely analyze the movement of cylindrical particles in the fluid. The two-way coupling model can give more precise information of gas-cylindrical particle flow. Here the simulated velocity vector and void fraction of fluid on one section at one moment was obtained and then handled by Tecplot software. Fig. 4 shows the velocity vector and void fraction of fluid on the section along coordinate center of y-axis of the riser at moment of 2.47s. Figs. 4 (a) and 4(b) are the full figure and the partial 20

ACS Paragon Plus Environment

Page 20 of 37

Page 21 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

enlarged drawing of the section respectively. It can be seen that, especially in the partial enlarged drawing, at that moment cylindrical particles distribute randomly and unequally in the riser. The magnitude and direction of local velocity of fluid both will change if some cylindrical particles occupy this region, and the magnitude of local velocity of fluid will decline. The total velocity of fluid profile in the riser still confirms to the basic velocity of fluid distribution principle that the velocities of fluid in mainstream tend to be equal in horizontal direction, and decline quickly up to zero in near boundary layer region.

Figure 4. (a) Full figure

Figure 4. (b) Partial enlarged drawing

Figure 4. Velocity vector and void fraction of fluid on one section at moment of 2.47s ( Vin = 5m/s ) Figs. 5 and 6 show the simulated results and experimental snapshots of the 21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 37

fluidized cylindrical particles in the riser, respectively. Due to the visual angle of the camera, only one section of the experimental riser has been focused. The number of cylindrical particles that were fluidized and present in the entire riser was around 5000. We can see that the predicted gas-cylindrical particle flow features in terms of overall orientation and spatial distribution of cylindrical particles qualitatively agree with the experimental measurements. Due to the inevitable numerical error and the inherent flaws of measuring methodologies and techniques in experiments, minor deviation is observed between simulation and experiment. Specifically, the velocities of cylindrical particles in the experiments appear slightly greater than those obtained in the simulations. It can also be seen from Figs. 5 and 6 that the cylindrical particles rotate frequently while moving upward. The local number concentration of cylindrical particles in the riser changes significantly over time and the number concentration in near-wall regions is higher than that in central regions.

0.05s

0.23s

0.45s

0.70s

0.88s

1.15s

1.31s

1.46s

Figure 5. Numerical results of fluidized cylindrical particle in the riser 22

ACS Paragon Plus Environment

Page 23 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

( Vin = 5m/s )

Figure 6. Experimental photos of fluidized cylindrical particles in the riser ( Vin = 5m/s ) The orientation of cylindrical particles at several time instants was quantified in both simulations and experiments. All time instants were chosen after the gas-cylindrical particle flow reached the statistically steady state, i.e., after some cylindrical particles had passed through the exit of the riser. Accordingly, the statistically averaged values of the time instants after the gas-cylindrical particle flow reached the fully-fluidized state were used for analysis. Fig. 7 shows the statistically averaged values of the simulated orientation of cylindrical particles in the riser at moments of 1.48s, 1.51s, 1.54s and 1.57s and those of the experimental orientation at moments of 1.2s, 1.25s, 1.3s and 1.35s, respectively. In the following numerical studies, all time instants were chosen after the gas-cylindrical particle flow reached the fully-fluidized state, and the statistically averaged values are referred to as the values at moments of 2.1s, 2.15s, 2.2s and 2.25s. It is found that the numerical results are in good agreement with the experimental data. The majority of the cylindrical particles move upwards in the riser with a small nutation angle range of 0-45°. In other words, the axes of most of the upward moving cylindrical particles tend to be in parallel with the z-axis of the riser. It is also worth noting that the nutation angle of the vast majority of the cylindrical particles falls in the range of 15°-30°. Cylindrical particles that move upwards with nearly-horizontal or upright axes are very sparse. Fig.8 shows the experimental and simulated particles time consumption ratios arriving 23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

at the exit at different radial regions, and the time consumption at near-wall regions is set as 1 and as the reference time. It is found that in consideration of experimental error the simulated results are in good agreement with experimental results in a whole and that the particles time consumption increases with the incremental distance to the walls.

Figure 7. Experimental and simulated orientations of cylindrical particles in the rise ( Vin = 5m/s )

24

ACS Paragon Plus Environment

Page 24 of 37

Page 25 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 8. Experimental and simulated particles time consumption ratios arriving at the exit at different radial regions ( Vin = 5m/s )

4.2. Axial orientation of cylindrical particles in a riser

Figure 9. Orientations of cylindrical particles in different radial regions of the riser ( Vin = 5m/s ) 25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In the presence of turbulence, the gas-cylindrical particle flow system is very complex due to the particle-particle collisional energy dissipation and the non-linearity of coupling interactions between turbulent flow and cylindrical particles. The orientations of cylindrical particles at different height of the riser are different. Hence, it is important to investigate the orientation change of cylindrical particles along the height. The orientation parameters of cylindrical particles at different time instants during the fluidization were calculated. Fig. 9 shows the orientation of cylindrical particles at different heights of the riser. It is found that the orientation difference of cylindrical particles becomes more distinct at a higher elevation. Accordingly, as the height increases, the number percentage of cylindrical particles that are moving upwards with a nutation angle range of 15°-30° increases, whilst that of cylindrical particles with a nutation angle range of 60°-90° decreases. By comparison, it can also be seen that the orientation of cylindrical particles in the middle lower region of the riser changes more significantly than those in other regions. Analysis suggests that the difference of velocity profile of fluid field from the bottom to the top of the riser results in the situation. From the bottom to the top of the riser, the velocity profile changes from a horizontal line to a peak-pattern curve of which middle is horizontal and two franks are cliffy step by step. Moreover, the potential energy of cylindrical particles tends to get smaller and smaller during the fluidization so as to keep the system stable.

26

ACS Paragon Plus Environment

Page 26 of 37

Page 27 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

4.3. Radial orientation of cylindrical particles in a riser

Figure 10. Orientations of cylindrical particles in different axial regions of the riser ( Vin = 5m/s ) Wall effect is one of the key factors that affect the gas-solid flow, in particular for cases where cylindrical particles are present. In this study, the effect of the walls on the orientation of cylindrical particles has been investigated. As some orientation parameters of cylindrical particles in different radial regions of the riser at different time instants were obtained, the statistically averaged values were used for analysis. Fig. 10 shows the radial distribution of the orientation of cylindrical particles in the riser. It is found that the orientation of cylindrical particles is unique in the central region. The number percentage of cylindrical particles moving upwards with nutation angle range of 25°-50° is the highest, and that of cylindrical particles moving up with nutation angle range of 0-15° is the lowest. This might be due to the centrosymmetric and radial flow velocity gradient distribution in the region. In near-wall regions, the number percentage of cylindrical particles moving upwards with nutation angle range of 0-30° is the highest, namely, more cylindrical particles tend to parallel the z-axis during the fluidization probably due to the wall effect. In other words, the wall 27

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

enforces the cylindrical particles near it to parallel the wall. With the increment of distance, the effect becomes weaker. In other radial regions, the orientation characteristics of cylindrical particles are similar.

4.4. Orientation of cylindrical particles at different slenderness ratios

Figure 11. Orientations of cylindrical particles at different slenderness ratios ( Vin = 5m/s ) The slenderness ratio is one of the key parameters characterizing the shape of cylindrical particles. If the slenderness ratio approaches 1, the cylindrical particles can be regarded as a spherical particle. For high slenderness ratios, typical fluidization behavior such as rotation and asymmetrical collision would occur. In this study, the effect of slenderness ratio on the orientation of cylindrical particles has been explored. Fig. 11 shows the statistically averaged values of the orientation of cylindrical particles with slenderness ratios of 6, 8, 10 and 12. We can see that the orientation of cylindrical particles changes significantly as the particle slenderness ratio varies. It is evident that the slenderness ratio plays an important role in determining the orientation of cylindrical particles. For small slenderness ratios, the orientation of 28

ACS Paragon Plus Environment

Page 28 of 37

Page 29 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

cylindrical particles exhibits the single-hump distribution. Moreover, the number percentage of cylindrical particles moving upwards with a nutation angle about 45° is the highest, whilst the number percentage of cylindrical particles with a nutation angle range of 0-15° or 75°-90° is the lowest. As the slenderness ratio increases, the number percentage of cylindrical particles moving upwards with a smaller nutation angle range of 0-30° increases, whereas that of cylindrical particles with a broader nutation angle range decreases. The results imply that cylindrical particles with a greater slenderness ratio tend to be in parallel with the direction of the gas flow during fluidization probably because the longer cylindrical particles affected by velocity profile of fluid seemingly more deeply and it is easier for longer cylindrical particles to keep the system stable.

4.5. Orientation of cylindrical particles at different Reynolds numbers

Figure 12. Orientations of cylindrical particles in the riser at different Re Fig. 12 shows the orientation of cylindrical particles at different Reynolds numbers i.e., 166444, 183089, 199733 and 216378. The formula calculating the Reynolds number is Re =

uin d ed

υg , d ed is the equivalent diameter of the riser. 29

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Remarkably, it is found that there seems no much difference in the orientation of cylindrical particles under these Reynolds numbers. Therefore, for the cases examined in this study (i.e., turbulent gas flow conditions), the Reynolds number (i.e. inlet wind velocity) has marginal effect on the orientation of cylindrical particles, though a different Reynolds number might lead to distinct fluidization behavior of cylindrical particles. Analysis suggests that the velocity profiles in all turbulent flows with higher Reynolds numbers are similar so the orientations of cylindrical particles in different Reynolds numbers are similar. Because previous research indicated that the orientations of cylindrical particles in laminar flow and in turbulent flow are different due to evidently different velocity profiles30. However, further work is warranted to investigate the variation of the orientation of cylindrical particles under a broader range of Reynolds numbers.

5. CONCLUSIONS In this study a 3D multi-way coupling model has been developed to investigate the orientation of cylindrical particles in circulating fluidized beds. The two-way coupling between cylindrical particles and turbulent flow was solved using the correlation between Lagrangian time scales and the k - ε model. Collisions of cylindrical particles were solved by the rigid impact dynamics and the modified Nabu collision probability method. The model of gas-cylindrical particle flow was verified by some experimental data, including the general orientation and spatial distribution of cylindrical particles in a riser. The numerical results indicated that the majority of cylindrical particles moved upwards with a small nutation angle. Cylindrical particles moving upwards with nearly-horizontal or upright axes were very sparse. The number of cylindrical particles with a small nutation angle increased along the height of the riser. In new-wall regions, the number percentage of cylindrical particles with a nutation angle range of 0-30° reached the peak. The slenderness ratio played a very important role in determining the orientation of cylindrical particles. Cylindrical particles with a high 30

ACS Paragon Plus Environment

Page 30 of 37

Page 31 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

slenderness ratio were prone to orientate in parallel with the gas flow direction during fluidization. Under the turbulent gas flow conditions, the effect of Reynolds number on the orientation of cylindrical particles was found to be marginal.

AUTHOR INFORMATION Corresponding author: Jie Cai. Tel.: 00862585481263. E-mail: [email protected]; Charley Wu. Tel.: 01483683506. E-mail: [email protected].

ACKNOWLEDGEMENT This worked was supported by National Natural Science Foundation of China (51306094), supported by China Postdoctoral Science Foundation (2014M551490), and supported by Postdoctoral Science Foundation of Jiangsu Province (1401003A).

NOTATION Symbols

∆t = time step, s H riser , Wriser = height and width of riser, m Ai = axial region distribution in the riser, m Ri = radial region distribution in the riser, m RT = slenderness of cylindrical particle

Vin = inlet wind velocity, m/s vg = velocity of fluid, m/s

vs = average velocity of all discrete elements, m/s v ( si ) = velocity of a discrete element, m/s

uini = initial velocity of cylindrical particle, m/s 31

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2

g = acceleration of gravity, m/s

m = mass of a cylindrical particle, kg p = the pressure of flow field, Pa

f sg = reacted force cylindrical particles to turbulence, N Ar = the modified Archimedes number

C Ds = drag coefficient I k , Iε = the effect of solid phase on turbulent energy and turbulent dissipation rate F = drag force acting on the cylindrical particle, N

Fi ( si ) = drag force acting on the discrete element of cylindrical particle, N

Gij = relative velocity between discrete element i and j, m/s. M , M * = the torque acting on the centroid of a cylindrical particle in the fixed

reference frame and in the body axes, N.m Greek letters ξ , η , ζ = 3 axis symbols of the body axes x , y , z = 3 axis symbols of the fixed reference

ρ g , ρs

= air

density and stalk density, kg/m3

υg = kinematic viscosity of gas, m2/ s ψ , θ , φ = precession angle, nutation angle and spin angle, rad ϕg = void fraction of fluid µ , µt = air viscosity and turbulent viscosity, Pa.s

ϑ = the angle between average velocity of discrete elements and average relative velocity, rad

k = turbulent kinetic energy, m2/s2 32

ACS Paragon Plus Environment

Page 32 of 37

Page 33 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

ε = turbulent dissipation rate, m2/s3 τ t ,sg = Lagrangian time scale Lt = the length scale of turbulence vortex * = angular velocity and initial angular velocity of cylindrical particle in ω* , ωini

the body axes, rad/s Subscripts g = air s = cylindrical particle (solid phase) ini = initial t = turbulence in = inlet * = in the body axes

REFERENCES (1)

Zgheib, N; Bonometti, T; Balachandar, T. Direct numerical simulation of cylindrical particle-laden gravity currents. Comput. Fluids. 2015, 123, 23-31.

(2) Lewandowski, E.P.; Jr, M. C.; Botto, L.; Bernate, J. C., et al. Orientation and self-assembly of cylindrical particles by anisotropic capillary interactions. AIChE J. 2010, 26, 15142-15154. (3) Mongruel, A.; Lecoq, N.; Wajnryb, E.; Cichocki, B.; Feuillebois, F. Motion of a sphero-cylindrical particle in a viscous fluid in confined geometry. Eur. J. Mech. B/Fluids. 2011, 30, 405-408. (4) Dalabaev, U. Numerical investigation of the character of the lift on a cylindrical particle in Poiseuille flow of a plane channel. J. Eng. Phys. Thermophys. 2011, 84, 1388-1392. (5) Batchelor, G. K. Slender-body theory for particles of arbitrary cross-section in Stokes flow. J. Fluid Mech. 1970, 44, 419-440. (6)

Saintillan, D.; Shaqfeh, E. S. G.; Darve, E. The growth of concentration 33

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

fluctuations in dilute dispersions of orientable and deformable particles under sedimentation. J. Fluid Mech. 2006, 53, 347-388. (7)

Fan, L.; Mao Z. S.; Wang, Y. D. Numerical simulation of turbulent solid-liquid two-phase flow and orientation of slender particles in a stirred tank. Chem. Eng. Sci. 2005, 60, 7045-7056.

(8) Li, S. Q.; Marshall, J. S. Discrete element simulation of micro-particle deposition on a cylindrical fiber in an array. Aerosol Sci. 2007, 38, 1031–1046. (9) Camassa, R.; Leiterman, T. J.; Mclaughlin, R. M. Trajectory and flow properties for a rod spinning in a viscous fluid. Part 1: An exact solution. J. Fluid Mech. 2008, 612, 153-200. (10) Lubansky, A. S.; Boger, D. V.; Cooper-White, J. J. Batchelor’s theory extended to elongated cylindrical or ellipsoidal particles. J. Non-Newtonian Fluid Mech. 2005, 130, 57–61. (11) Paschkewitz, J. S.; Dubief, Y.; Dimitropoulos, C. D. Numerical simulation of turbulent drag reduction using rigid fibres. J. Fluid Mech. 2004, 518, 281-317. (12) Tsuji, Y.; Tanaka, T.; Ishida, T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technol. 1992, 71, 239-250. (13) Delaney, G.W.; Morrison, R.D.; Sinnott, M.D.; Cummins, S.; Cleary, P.W. DEM modelling of non-spherical particle breakage and flow in an industrial scale cone crusher. Miner. Eng. 2015, 74, 112-122. (14) Cleary, P. W.; Stokes, N. Efficient collision detection for three dimensional super-ellipsoidal particles. Computational Techniques and Applications: Process of 8th International Computational Technique and Applications Conference: CTAC97, Singapore, 1997, 139-144. (15) Kodam, M.; Bharadwaj, R.; Curtis, J.; Hancock, B.; Wassgren, C. Cylindrical object contact detection for use in discrete element method simulations Part I Contact detection algorithms. Chem. Eng. Sci. 2010, 65, 5852-5862. (16) Wachs, A.; Girolami, L.; Vinay, G.; Ferrer, G. Grains3D, A flexible DEM approach for particles of arbitrary convex shape—Part I: Numerical model and validations. 34

ACS Paragon Plus Environment

Page 34 of 37

Page 35 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Powder Technol. 2012, 224: 374-389. (17) Boon, C.W.; Houlsby, G.T.; Utili, S. A new contact detection algorithm for three-dimensional non-spherical particles. Powder Technol. 2013, 248, 94-102. (18) Ardema, M. D. Rigid dynamics. California USA: Springer. 2005. (19) Yin, C.; Rosendahl, Lasse.; Kӕr, S.K.; Sфrensen H. Modelling the motion of cylindrical particles in a nonuniform flow. Chem. Eng. Sci. 2003, 58, 3489-3498. (20) Grof, Z.; Kohout, M.; Štěpanek, F. Multi-scale simulation of needle-shaped particle breakage under uniaxial compaction. Chem. Eng. Sci. 2007, 62, 1418-1429. (21) Feng, Y.; Kleinstreuer, C. Analysis of non-spherical particle transport in complex internal shear flows. Phys. Fluids. 2013, 25, 091904. (22) Mortensen,P.H.; Andersson,H.I.; Gillissen,J.J.J.; Boersma, B.J. On the orientation of ellipsoidal particles in a turbulent shear flow. Int. J. Multiphase Flow. 2008, 34, 678-683. (23) Njobuenwu, D. O.; Fairweather, M. Dynamics of single, non-spherical ellipsoidal particles in a turbulent channel flow. Chem. Eng. Sci. 2015, 123, 265-282. (24) Crowe, C. T. On models for turbulence modulation in fluid-particle flows. Int. J. Multiphase Flow. 2000, 26, 719-727. (25) Wen, C. Y.; Yu, Y. H. Mechanics of fluidization. Chem. Eng. Prog. Sump. Ser.

1966, 162, 100e13. (26) Fan, L.; Mao, Z. S.; Yang, C. Experiment on settling of slender particles with large aspect ratio and correlation of the drag coefficient. Ind. Eng. Chem. Res. 2004, 43, 7664-7670. (27) Simonin, O.; Deutsch, E.; Minier, J. P. Eulerian prediction of the fluid-particle correlated motion in turbulent two-phase flows. Appl. Sci. Res. 1993, 51, 275-283. (28) Elghobashi, S.; Truesdell, G. C. On the two-way interaction between homogeneous turbulence and dispersed solid particles. I: Turbulence modification. Phys. Fluids A. 1993, 5, 1790-1801. (29) Smith, C. E. Predicting rebounds using rigid-body dynamics. J. Appl Mech. 1991, 35

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

58, 754-758 (30) Zhang, W. F.; Lin, J. Z. Research of the motion of particles in the turbulent pipe flow of fiber suspensions. Appl. Math. Mech. 2004, 25, 741-750.

36

ACS Paragon Plus Environment

Page 36 of 37

Page 37 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

For Table of Contents Only

Fluidization of cylindrical particles in a circulating fluidized bed and collisions among cylindrical particles were handled with rigid body dynamics

37

ACS Paragon Plus Environment