MATHEMATICAL MODELING OF INTERACTION OF WET PARTICLES AND APPLICATION TO FLUIDIZED BEDS by Pirooz Darabi B.Sc., Sharif University of Technology, 2004 M.A.Sc., The University of British Columbia, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering) The University of British Columbia (Vancouver) May 2011 © Pirooz Darabi, 2011 - ii - ABSTRACT In many industrial operations, such as fluidized bed granulators, coaters, and fluid cokers, a binding or reacting liquid is introduced into the system. Due to the effects of liquid, the multi- phase transport phenomena of these systems are more complicated compared to conventional gas-solid fluidization systems. In this thesis, mathematical modeling is used to study the interaction of wet particles. First, a coalescence model is developed to describe the binary collision of wet particles. The model is in the form of a wet coefficient of restitution and is used to determine the critical velocitythe boundary between coalescence or rebound outcomesfor a range of capillary numbers. Model predictions are compared with the available experimental data and good agreement is found. The model accounts for both liquid viscosity and surface tension effects and is used to investigate the boundary between collisions with dominant capillary and respective viscous effects. Then, by incorporating time- and temperature-dependant variations of the viscosity and thickness of the liquid coating, the model is used to determine the agglomeration tendency of bitumen-coated coke particles in fluid cokers. A simplified mathematical model and numerical solution of the Navier-Stokes equations are used to study the rupture of stretching liquid bridges between two solid spherical particles. The simplified model considers the geometry of the problem in which the gas-liquid interface is represented with a parabola. The numerical simulations of the Navier Stokes equations are performed with FLUENT and are used to investigate the viscous, surface tension, inertial, gravitational, and contact angle effects on the rupture distance and liquid distribution. Finally, the interaction of multiple wet particles is addressed by implementing the wet coefficient of restitution proposed in this thesis, using MFIX, an open-source Discrete Element Method (DEM) tool. DEM simulations of a fluidized bed consisting of mono-sized solid spherical particles pre-coated with identical liquid coatings are performed, and the effect of coating viscosity and thickness on the fluidization behaviour is investigated. Snapshots of the instantaneous particle positions are presented, and time-averaged values of the bed centroid in the y-direction, wet coefficient of restitution, and relative normal collision velocity are analyzed. - iii - PREFACE This thesis entitled “Mathematical modeling of interaction of wet particles and application to fluidized beds” presents the research led and performed by Pirooz Darabi. The research conducted in this thesis was supervised, identified, and initiated by Martha Salcudean and co- supervised by Dana Grecov. • A version of Chapter 2 has been published. Darabi et al. (2009) A Novel Coalescence Model for Binary Collision of Identical Wet Particles. Chemical Engineering Science. 64: 1868-1876. The research, model development and coding of the equations, analyses of results, and drafting of the paper were performed by Pirooz Darabi. Konstantin Pougatch, Martha Salcudean, and Dana Grecov provided feedback and suggestions. • A version of Chapter 3 has been published. Darabi et al. (2010) Agglomeration of Bitumen- coated Coke Particles in Fluid Cokers. International Journal of Chemical Reactor Engineering. 8, A122. The research, model development and coding of the equations, analyses of results, and drafting of the paper were performed by Pirooz Darabi. Konstantin Pougatch, Martha Salcudean, and Dana Grecov provided feedback and suggestions. • A version of Chapter 4 (containing sections 4.2 up to 4.7) has been published. Darabi, P., Li, T., Pougatch, K., Salcudean, M., and Grecov, D. (2010) Modeling the Evolution and Rupture of Stretching Pendular Liquid Bridges. Chemical Engineering Science. 65: 4472-4483. The research, model development and coding of the equations, analyses of results, and drafting of the paper were performed by Pirooz Darabi. Also, some help with the initial setup of the numerical simulations performed with FLUENT and writing of the related methodology was provided by Tingwen Li. Konstantin Pougatch, Martha Salcudean, and Dana Grecov provided feedback and suggestions. • Section 4.8 of Chapter 4 has also been published. Darabi et al. (2010) Numerical Studies of Stretching Liquid Bridges between Two Solid Spherical Particles with Different Contact Angles. ASME International Mechanical Engineering Congress and Exposition (IMECE). IMECE2010- - iv - 37564. The research, model development, analyses of results, and drafting of the paper were performed by Pirooz Darabi. Konstantin Pougatch, Martha Salcudean, and Dana Grecov provided feedback and suggestions. • Section 4.9 of this chapter has been submitted for publication. Darabi et al. Numerical Investigations of Liquid Bridge Rupture in the Presence of Gravity. Submitted on November 10, 2010. The research, model development, analyses of results, and drafting of the paper were performed by Pirooz Darabi. Konstantin Pougatch, Martha Salcudean, and Dana Grecov provided feedback and suggestions. • A version of Chapter 5 has been submitted for publication. Darabi et al. DEM Investigations of Fluidized Beds in the Presence of Liquid Coating. Submitted on December 13, 2010. The research, model development and coding of the equations, analyses of results, and drafting of the paper were performed by Pirooz Darabi. Konstantin Pougatch, Martha Salcudean, and Dana Grecov provided feedback and suggestions. - v - TABLE OF CONTENTS Abstract.......................................................................................................................................... ii Preface........................................................................................................................................... iii Table of Contents .......................................................................................................................... v List of Tables ................................................................................................................................ ix List of Figures................................................................................................................................ x Nomenclature ............................................................................................................................. xiv Acknowledgements .................................................................................................................... xix Dedication .................................................................................................................................... xx Chapter 1. Introduction........................................................................................................... 1 1.1. Motivation....................................................................................................................... 1 1.2. Particle interactions......................................................................................................... 3 1.2.1. Problem statement................................................................................................... 3 1.2.2. Coalescence modeling ............................................................................................ 4 1.2.3. Rupture phenomenon of liquid bridges................................................................... 4 1.3. Numerical methods for simulating multi-phase flows.................................................... 5 1.3.1. Eulerian-Eulerian approach .................................................................................... 5 1.3.2. Eulerian-Lagrangian approach................................................................................ 6 1.4. Literature review............................................................................................................. 7 1.4.1. Inter-particle forces and coalescence modeling...................................................... 7 1.4.2. Size enlargement of particles in fluidized beds ...................................................... 9 1.4.3. Rupture and liquid distribution of liquid bridges.................................................. 11 1.4.4. Discrete Element Method simulations of wet fluidized beds ............................... 14 1.5. Thesis objectives........................................................................................................... 17 1.5.1. Approach............................................................................................................... 17 1.6. Thesis layout ................................................................................................................. 18 Chapter 2. A Coalescence Model for Binary Collision of Wet Particles........................... 20 2.1. Introduction................................................................................................................... 20 2.2. Model development ...................................................................................................... 21 2.2.1. Problem statement................................................................................................. 21 - vi - 2.2.2. Assumptions.......................................................................................................... 23 2.2.3. General equation of motion .................................................................................. 23 2.2.4. Methodology......................................................................................................... 24 2.2.5. Proposed coalescence model................................................................................. 25 2.2.6. Viscous limiting case ............................................................................................ 30 2.2.7. Capillary limiting case .......................................................................................... 32 2.3. Results and discussion .................................................................................................. 32 2.3.1. Comparison of analytical and numerical results ................................................... 32 2.3.2. Coalescence with dominant viscous effects.......................................................... 35 2.3.3. Coalescence with dominant capillary effects........................................................ 36 2.3.4. Comparison with experimental data ..................................................................... 37 2.4. Summary and conclusions ............................................................................................ 41 Chapter 3. Agglomeration of Bitumen-coated Coke Particles in Fluid Cokers ............... 43 3.1. Introduction................................................................................................................... 43 3.2. Fluid coker .................................................................................................................... 43 3.2.1. Flow regions within a fluid coker ......................................................................... 44 3.3. Mathematical model...................................................................................................... 45 3.3.1. Problem statement................................................................................................. 45 3.3.2. Proposed model..................................................................................................... 46 3.3.3. Solution methodology........................................................................................... 50 3.4. Results and discussion .................................................................................................. 51 3.4.1. Variations of the critical velocity with the reaction time...................................... 51 3.4.2. Estimating the maximum possible agglomerate size ............................................ 55 3.5. Summary and conclusions ............................................................................................ 56 Chapter 4. Modeling the Rupture of Stretching Liquid Bridges....................................... 57 4.1. Introduction................................................................................................................... 57 4.2. Problem statement......................................................................................................... 58 4.3. Non-dimensional analysis............................................................................................. 58 4.4. Simplified model........................................................................................................... 59 4.4.1. Model development .............................................................................................. 60 4.4.2. Solution methodology........................................................................................... 63 - vii - 4.4.3. Results and discussion .......................................................................................... 63 4.5. Numerical simulation of the Navier–Stokes equations................................................. 67 4.5.1. Model development .............................................................................................. 67 4.5.2. Solution methodology........................................................................................... 70 4.5.3. Grid independence study....................................................................................... 70 4.5.4. Results and discussion .......................................................................................... 71 4.6. Comparison between simplified model and numerical simulation............................... 76 4.6.1. Comparison between rupture distances predicted with simplified model and numerical simulation............................................................................................................. 76 4.6.2. Variations of rupture distance with Weber, capillary and Bond numbers ............ 77 4.6.3. Comparison between liquid bridge profiles of simplified model and numerical simulation.............................................................................................................................. 80 4.7. Limitations of simplified model.................................................................................... 83 4.8. Influence of contact angle............................................................................................. 84 4.8.1. Summary of cases ................................................................................................. 84 4.8.2. Grid independence study....................................................................................... 85 4.8.3. Results and discussion .......................................................................................... 86 4.9. Influence of gravity....................................................................................................... 92 4.9.1. Summary of cases ................................................................................................. 92 4.9.2. Results and discussion .......................................................................................... 92 4.10. Summary and conclusions ........................................................................................ 97 Chapter 5. DEM Investigations of Fluidized Beds in the Presence of Liquid Coating.... 99 5.1. Introduction................................................................................................................... 99 5.2. Model description ....................................................................................................... 100 5.2.1. Gas-phase............................................................................................................ 100 5.2.2. Solid-phase.......................................................................................................... 101 5.2.3. Model for interaction of wet particles................................................................. 103 5.2.4. Simulation setup.................................................................................................. 105 5.3. Results and discussion ................................................................................................ 109 5.3.1. Bed behaviour ..................................................................................................... 109 5.3.2. Average particle height ....................................................................................... 113 - viii - 5.3.3. Wet coefficient of restitution .............................................................................. 117 5.4. Summary and conclusions .......................................................................................... 119 Chapter 6. Summary, Conclusions, and Recommendations ............................................ 120 6.1. Summary and conclusions .......................................................................................... 120 6.1.1. Binary interaction of particles............................................................................. 120 6.1.2. DEM investigations of wet fluidized beds.......................................................... 122 6.2. Contributions............................................................................................................... 123 6.3. Recommendations....................................................................................................... 124 Bibliography .............................................................................................................................. 126 Appendix A. Viscous Force Expression between Two Spheres....................................... 138 Appendix B. Capillary Force Expression between Two Spheres.................................... 141 Appendix C. Derivation of the Overall Coefficient of Restitution.................................. 143 Appendix D. Derivation of the Velocity Expressions ....................................................... 145 Appendix E. Derivation of the Viscous Stokes Model...................................................... 148 - ix - LIST OF TABLES Table 2-1. Summary of the values of parameters. ........................................................................ 33 Table 3-1. Characteristic inter-particle collision velocity and temperature of the flow regions within fluid cokers. ....................................................................................................................... 45 Table 3-2. Properties of the wet coke particles............................................................................. 50 Table 4-1. Summary of the simulation conditions........................................................................ 69 Table 4-2. Values for the convergence-related parameters. ......................................................... 70 Table 4-3. Summary of grid independence study. ........................................................................ 71 Table 4-4. Summary of the cases considered in this section. ....................................................... 84 Table 4-5. Summary of the grid independence study for case 5................................................... 85 Table 4-6. Summary of the grid independence study for case 6................................................... 85 Table 4-7. Summary of the cases considered in this section. ....................................................... 92 Table 5-1. Summary of the constitutive equations of the wet coefficient of restitution model.. 104 Table 5-2. Summary of the simulation settings. ......................................................................... 107 Table 5-3. Summary of the coating cases. .................................................................................. 109 Table 5-4. Solid particle and liquid coating properties............................................................... 109 - x - LIST OF FIGURES Figure 1-1. A schematic diagram of a fluid coking system. ........................................................... 2 Figure 1-2. Binary interaction of wet particles: a) approach, b) coalescence, c) liquid bridge about to rupture, and d) liquid distribution (immediately after rupture)......................................... 4 Figure 2-1. Approach of two moving spherical particles and the formation of the liquid bridge. 21 Figure 2-2. Rebound of two moving spherical particles (a) and rupture of liquid bridge (b). ..... 22 Figure 2-3. Variations of the analytical and numerical critical velocities with respect to the capillary number. .......................................................................................................................... 34 Figure 2-4. Comparison between values of the critical velocity obtained with the general coalescence model, viscous limiting model, and capillary limited model.................................... 35 Figure 2-5. Comparison between the values of the critical velocity obtained with the general coalescence model, the viscous limiting model, and different versions of the Stokes model. ..... 36 Figure 2-6. Variations of the capillary work, rupture energy, and contact loss with respect to the dry coefficient of restitution.......................................................................................................... 37 Figure 2-7. Variations of the critical Stokes number with the respect to the initial height of liquid layer or three different particles sizes: a) R = 0.8 mm, b) R = 1.0 mm, and c) R = 2.0 mm. ....... 40 Figure 2-8. Comparison between experimental and analytical values of rupture energy............. 41 Figure 3-1. Collision of two moving spherical particles and the formation of the liquid bridge. 45 Figure 3-2. Variations of the thickness of bitumen coating with the normalized reaction time... 49 Figure 3-3. Variations of the critical velocity with the normalized reaction time for T = 400 oC.53 Figure 3-4. Variations of the critical velocity with the normalized reaction time for T = 503 oC.54 Figure 3-5. Variations of the critical velocity with the normalized reaction time for T = 530 oC.54 Figure 4-1. A stretching liquid bridge between two identical solid spherical particles................ 58 Figure 4-2. Liquid bridge between two spherical particles, shown at the instant of rupture........ 60 Figure 4-3. The liquid droplet on particle i immediately after rupture. ........................................ 62 Figure 4-4. Comparison of rupture distances predicted using the simplified model with the experimental data of Mason and Clark (using water as the suspending medium and di-n-butyl phthalate and liquid paraffin mixture for the liquid bridge) [114]................................................ 64 - xi - Figure 4-5. Comparison of rupture distances predicted using the simplified model and the experimental data of Mazzone et al. [85] (with measured liquid-solid contact angles of 10 o and using air as the suspending medium and dibutyl phthalate for the bridging liquid)..................... 65 Figure 4-6. Variations of the dimensionless rupture distance with respect to the contact angle and dimensionless bridge volume........................................................................................................ 66 Figure 4-7. Comparison of rupture distances predicted using the simplified model with the expressions of Lian et al. [81] and Mikami et al. [82].................................................................. 66 Figure 4-8. The computational domain used in this study. The blue and red regions represent the gas and liquid phases, respectively, and the boundary in-between represents the gas-liquid interface......................................................................................................................................... 68 Figure 4-9. Numerical mesh and boundary conditions. This is for an intermediate grid that initially has 10800 control volumes.............................................................................................. 69 Figure 4-10. Simulated bridge evolution in time in the absence of gravitational effects: a. t=0, b. t=0.025, c. t=0.05, d. t=0.075, e. t=0.1, f. t=0.125, g. t=0.1335, h. t=0.134, i. t=0.1341, j. t=0.1342, k. t=0.1343, and l. t=0.1344 ms. ................................................................................... 72 Figure 4-11. Simulated bridge evolution in time in the presence of gravitational effects: a. t=0.015, b. t=0.065, c. t=0.09, d. t=0.115, e. t=0.117, and f. t=0.1175 ms. .................................. 73 Figure 4-12. Comparison between our numerical simulations and the experimental data of Mason and Clark [114]. ............................................................................................................................ 74 Figure 4-13. Comparison between the liquid transfer fraction predictions of our numerical simulations and the numerical simulations of Darhuber et al. [118]. ........................................... 75 Figure 4-14. Comparison of the rupture distance predictions using both the numerical simulations and the simplified model with the experimental data of Mazzone et al. (1986). ...... 76 Figure 4-15. Variations of the simulated rupture distance with respect to the Weber number. ... 78 Figure 4-16. Variations of the simulated rupture distance with respect to the capillary number. 79 Figure 4-17. Variations of the simulated rupture distance with respect to the Bond number. ..... 80 Figure 4-18. Comparison between the bridge profiles predicted using the simplified model and numerical simulations at: (a) rupture instant, and (b) a few milliseconds before rupture, for θ = 30o. ................................................................................................................................................ 81 Figure 4-19. Contour plots of the velocity magnitude (in m/s) of the numerical simulations at: (a) rupture instant, and (b) a few milliseconds before rupture, for θ = 30o........................................ 82 - xii - Figure 4-20. Comparison of the bridge profiles predicted using the simplified model and numerical simulation at: (a) rupture instant, and (b) a few milliseconds before rupture, with θi = 30o and θj = 40o. ............................................................................................................................ 83 Figure 4-21. Variations of the rupture distance with respect to contact angle for cases 1 and 2.. 86 Figure 4-22. Variations of the rupture distance with respect to the capillary number for cases 3 and 4.............................................................................................................................................. 87 Figure 4-23. Variations of the rupture distance with respect to the capillary number for cases 5 and 6.............................................................................................................................................. 88 Figure 4-24. Variations of the liquid transfer fraction with respect to contact angle for cases 1 and 2.............................................................................................................................................. 89 Figure 4-25. Variations of the liquid transfer fraction with respect to the capillary number for cases 3 & 4. ................................................................................................................................... 91 Figure 4-26. Variations of the liquid transfer fraction with respect to the capillary number for cases 5 & 6. ................................................................................................................................... 91 Figure 4-27. Variations of the rupture distance with respect to the capillary number: (a) θl = 30 & θr = 40o, (b) θl = 30 & θr = 50o....................................................................................................... 94 Figure 4-28. Variations of the liquid transfer fraction with respect to the Bond number............. 95 Figure 4-29. Variations of the transfer fraction with respect to the capillary number: (a) θl = 30 & θr = 40o, (b) θl = 30 & θr = 50o....................................................................................................... 96 Figure 5-1: Schematic representation of the 2D simulation domain. ......................................... 106 Figure 5-2: Snapshots of the simulated transient bed evolution for case C3 at superficial velocity 2.8 m/s with: a) inviscid liquid coating, and b) coating with viscosity 100 mPa.s..................... 111 Figure 5-3: Snapshots of the simulated instantaneous particle positions for case C3 at superficial velocity 3.36 m/s with: a) inviscid liquid coating, and b) coating with viscosity 100 mPa.s..... 112 Figure 5-4: Time-averaged average particle height for different superficial velocities and different coating viscosities for: a) case C1, and b) case C3. ..................................................... 114 Figure 5-5: Time-averaged wet coefficient of restitution for different superficial velocities and different coating viscosities for: a) case C1, and b) case C3. ..................................................... 115 Figure 5-6: Time-averaged relative normal collision velocity for different superficial velocities and different coating viscosities for: a) case C1, and b) case C3. .............................................. 116 - xiii - Figure 5-7: Summary of the variations of the time-averaged wet coefficient of restitution with respect to coating viscosity and coating thickness at superficial velocity of 2.8 m/s................. 118 Figure 5-8: Summary of the variations of the time-averaged normal relative collision velocity with respect to coating viscosity and coating thickness at superficial velocity of 2.8 m/s......... 118 - xiv - NOMENCLATURE a constant, mm A constant, 2 0 3 2 R upiµ= Bo Bond number ( 2 /gRρ σ= ), - Ca capillary number ( /vµ σ= ), - d liquid bridge length, - D minimum half-separation distance, mm (i)D diameter of ith particle, mm gD strain rate tensor of the gas phase, s -1 Dm diameter of solid particle belonging to mth solid group, mm e dry coefficient of restitution, - 'e overall coefficient of restitution ne normal coefficient of restitution, - ne r unit vector normal to particle surface along axis of collision ue r unit vector in the velocity direction wete wet coefficient of restitution, - ( ) cF i net contact force acting on ith particle, N Fcap capillary force, Nµ (i) dF drag force acting on i th particle, N ( ) TF i total force acting on ith particle, N Fvis viscous force, Nµ g gravitational acceleration, m/s2 h spherical cap height, - h, h0 thickness of liquid coating, mm or mµ h0 initial thickness of liquid coating, mµ h1 final half-separation distance, mm - xv - ha average height of particle asperities (or surface roughness), mµ or mm hrup rupture distance, mm ruph% dimensionless rupture distance, - H average particle height, mm H(r) separation distance, mµ or mm ( )iI moment of inertia of the ith particle, kg.m2 Igm momentum transfer term between the gas phase and the mth solid phase, kg/(m2.s2) nk normal spring stiffness coefficient, N/m Lc contact loss, nJ Lt actual energy dissipation, nJ Lt,max maximum possible energy dissipation, nJ m particle mass, gµ ( )i m mass of ith particle, kg lm mass of liquid coating, kg sm mass of solid particle, kg wm mass of wet particle, kg M number of solid phases, - Np total number of solid particles, - Nm number of solid particle belonging to mth solid group, - Nx number of cells in x-direction, - Ny number of cells in y-direction, - P pressure, mPa Pg pressure of the gas phase, kg/(m.s2) r radial distance, mµ or mm rd radius of the dry particle, mm rw radius of the wet particle, mm rw wetting distance, mµ or mm R particle radius, mµ or mm Ragg agglomerate radius, mµ - xvi - Re Reynolds number ( / /Rv We Caρ µ= = ), - S separation distance, - S* dimensional rupture distance, m cS dimensionless rupture distance ( * /cS S R= ), - gS stress-tensor of the gas phase, kg/(ms2) St Stokes number, - St* critical Stokes number, - t normalized reaction time, - t time, s ti initial time, s tf final time, s ts simulation time, s T temperature, oC T droplet maximum cap height, - (i)T total torque acting on ith particle, N.m TF liquid transfer fraction, - u droplet radius, - u velocity, m/s *u critical velocity, m/s u0 initial velocity, m/s ua approach velocity, m/s uc contact velocity, m/s uc characteristic inter-particle collision velocity, m/s uf final velocity (after impact), m/s * pu peak critical velocity, m/s ur rebound velocity, m/s v relative separation velocity, m/s gv volume-averaged velocity of the gas-phase, m/s V, bV% dimensionless liquid bridge volume ( 3/bV R= ), - - xvii - (i)V linear velocity of ith particle, m/s bV dimensional liquid bridge volume, m 3 capV spherical cap volume, - lV volume of liquid coating, m 3 sV volume of solid particle, m 3 wV volume of wet particle, m 3 Wcap capillary work, pJ Wvis viscous work, pJ We Weber number ( 2 /v Rρ σ= ), - x x-coordinate (abscissa) xliq moisture content, - xmin liquid bridge x-coordinate at its thinnest point X(i) position of ith particle, m Xv correction term y y-coordinate (ordinate) z axial distance, mµ or mm Greek Letters t∆ time step, s gε volume fraction of the gas phase, - nη normal dashpot damping coefficient, N.s/m θ contact angle, radian θ , Θ granular temperature, m2/s gλ second coefficients of viscosity of the gas phase, kg/(m.s) µ viscosity of the liquid coating, mPa.s gµ dynamic viscosity of the gas phase, kg/(m.s) ρ liquid density, kg/m3 (i)ρ density of ith particle, kg/m3 gρ density of the gas phase, kg/m3 - xviii - lρ density of the liquid coating, kg/m3 pρ particle density, kg/m3 sρ density of the dry particle, kg/m3 smρ density of solid particle belonging to mth solid group, kg/m3 wρ density of the wet particle, kg/m3 σ surface tension of coating liquid, mN/m gτ shear stress tensor of the gas phase, kg/(m.s2) (i)ω angular velocity of ith particle, s-1 - xix - ACKNOWLEDGEMENTS I would not have completed this PhD thesis without the significant help of many great people who have made my time at UBC in the past few years a pleasant and successful experience. First and foremost, I would like to express my sincere gratitude to my outstanding supervisor, Dr. Martha Salcudean, for her continuous support, warm encouragement, and never-ending patience throughout the course of my graduate studies. I consider myself fortunate and honoured to have had the opportunity to work under her supervision. I would also like to thank my co- supervisor, Dr. Dana Grecov, for her encouraging guidance and support. For his technical support and beneficial criticism, I gratefully acknowledge Dr. Konstantin Pougatch, without whom I could not have achieved this work. I would like to thank my research committee members, Dr. John Grace and Dr. Mark Martinez, for their helpful advice and suggestions. I would also like to thank Maureen Phillips very much for reviewing my thesis and improving its writing. To all my wonderful friends and colleagues, I would like to express my gratitude. In particular and in alphabetical order, thanks to AmirAbbas Aliabadi, Ali Asadkarami, Bijan Azadi, Ehsan Azadi, Hessam Bavafa, Meysam Bavafa, Ehsan Bayaki, Mahdi Beheshtian, Edward Chan, Massimo Diciano, Saeed Fathololoumi, Ahmad Ghahremaninezhad, Asef Javan, Ali Kashani, Tingwen Li, Hamed Mahmudi, MohammadSadegh Masnadi, Nasim Massah, Farzad Moghimi, Amir Nejat, Nina Rajabinasab, Hassan Rivaz, Shirin Sabouri, Ramin Sahebjavaher, Mohammad Sepasi, Nazanin Shabani, Ali Vakil, Soheyl Vakili, Hamid Yamini, and Nima Ziraknejad. Finally, I would like to express my everlasting gratitude to my wonderful parents, Parvin and Behrouz, and my dear sisters, Parissa and Parinaz, for their continuous encouragement, support, and love during my education. I have always felt them beside me in Canada, even though they are living far away. - xx - DEDICATION To my parents - 1 - Chapter 1. INTRODUCTION 1.1. MOTIVATION Western Canada has an estimated 4.6 x 1011 cubic meters of bitumen in its oilsands deposits [1]. At the same time, the conventional oil reserves are depleting and worldwide demands are increasing drastically. This offers Canada an excellent opportunity as a major oil supplier in the world. However, bitumen differs from conventional crude oils that are produced by the majority of oil suppliers. It is the heaviest and thickest form of petroleum. It has a high density and its viscosity can be orders of magnitude higher than that of conventional crude oils. Bitumen contains heavy hydrocarbons and because of that it is considered as heavy oil. The existing conventional refineries have difficulty handling bitumen in its original state; therefore, some additional treatment is required. Any treatment of heavy oils with the aim of increasing its value is referred to as “upgrading.” The ultimate goal of an upgrading process is the conversion of the heavy oil into lighter components, and eventually production of synthetic crude oil (SCO). The synthetic crude oil has chemo-physical characteristics comparable to conventional light crude oils. The process of converting heavy hydrocarbon into lighter fractions is known as cracking, one of the main stages of the current upgrading processes of bitumen. Coking is one of the cracking processes used in Canada and has two types: “delayed coking” used by Suncor Energy Inc. and “fluid coking” used by Syncrude Canada Ltd. In delayed coking, bitumen cracks in a coke drum. The process occurs at 450-510 oC and low pressures, i.e. 0-300 kPa [1]. Fluid coking takes place in large fluidized beds and is further explained in the following paragraph. Fluid coking is an established commercial cracking process and is used for converting bitumen into lighter hydrocarbons. Fluid coker is one of the two key components of a fluid coking system, schematically shown in Figure 1-1. It is a large steam-fluidized bed containing hot coke particles in which the bitumen is injected through a large number of nozzles using steam-assisted atomization. The atomized bitumen droplets contact and coat the coke particles. Due to the high - 2 - temperature of the coke particles, the heavy hydrocarbon molecules of bituminous liquid thermally crack and are converted into lighter products. As a result of cracking reactions, the coke particles lose energy, cool down, and grow in size (as bitumen remaining on the particle surface is converted into coke). In order to be reheated, the large cold coke particles are directed into a second large fluidized bed known as the burner, and once hot again, the coke particles are circulated back into the fluid coker. Figure 1-1. A schematic diagram of a fluid coking system. A fluid coker is a very complex multi-phase flow reactor. In particular, numerous coupled processes at various scales occur in them at high temperatures and elevated pressures. These result in limitations on real-scale measurements and observations of the processes occurring inside fluid cokers. In addition, as reported in the literature [1-3], mainly due to the imperfect liquid–solid mixing, fluid cokers are prone to agglomeration. As a result of the poor liquid-solid mixing, “thick” bitumen layers coat some of the coke particles. The coating layers act as a binder - 3 - that leads to the coalescence of coke particles. Eventually, very large granules consisting of several up to hundreds of coke particles are formed; this is known as agglomeration. Agglomeration in fluid cokers is undesirable and causes reduced process efficiency as: 1) agglomerates can be heavy and difficult to fluidize, leading to partial or complete defluidization of the reactor; and 2) agglomeration results in reduced heat transfer to the reacting media and consequently lowers product yields. Learning about conditions leading to agglomeration can be beneficial to improving fluid coker operation. Agglomeration in fluid cokers is similar to wet granulation in fluidized beds [1]. Hence, this thesis is not limited to the agglomeration problem and is formulated such that it can be applied to wet granulation processes as well. In order to obtain a fundamental understanding of wet granulation processes, this thesis starts by investigating the problem at a micro-scale, i.e. the binary interaction of wet particles; after that, the problem is extended to a larger scale and the interaction of numerous wet particles in fluidized beds is studied. 1.2. PARTICLE INTERACTIONS 1.2.1. Problem statement Figure 1-2.a schematically shows two solid spherical particles coated with a liquid layer approaching each other. Once the two particles collide and the liquid bridge is formed, there are two possible outcomes. If the dissipating effects exceed the kinetic energy of the system of two particles, the particles stick together and coalesce (Figure 1-2.b); otherwise, they rebound. As a result of their rebound, particles separate with a non-zero velocity and the liquid bridge between them then ruptures. Figure 1-2.c shows the liquid bridge at the instant of rupture. Consequently, the liquid is distributed between the two solid particles (Figure 1-2.d). The key terms of the above explanation, i.e. coalescence, rupture of liquid bridge, and liquid distribution are very briefly described as follows. - 4 - Figure 1-2. Binary interaction of wet particles: a) approach, b) coalescence, c) liquid bridge about to rupture, and d) liquid distribution (immediately after rupture). 1.2.2. Coalescence modeling In order to determine the collision outcome of the particles, numerous models, usually known as coalescence models, have been proposed in the literature. These models account for the total energy dissipation and the total available kinetic energy of the system consisting of (two) particles and determine the collision outcome accordingly. 1.2.3. Rupture phenomenon of liquid bridges Upon the rebound of wet particles, the liquid bridge between them breaks and its liquid volume is partitioned between the particles. The main two parameters to quantify the rupture and post- rupture evolution of liquid bridges are the “rupture distance” and “liquid transfer fraction”. The - 5 - rupture distance indicates the distance at which the liquid bridge breaks. The liquid transfer fraction is a measure of how the liquid is partitioned between the solid bodies after rupture. Both of these parameters are important and can be used for modeling and simulating multi-phase systems in the presence of a liquid binder. 1.3. NUMERICAL METHODS FOR SIMULATING MULTI-PHASE FLOWS In the past decades, the combination of efficient solution algorithms and high-speed computing units has provided great opportunities for researchers and scientists to numerically study complicated problems. The popularity of the numerical simulation as a tool for studying and understanding multi-phase flows is ever growing. There are two main numerical approaches used for simulating multi-phase flows: Eulerian-Lagrangian and Eulerian-Eulerian. These models are briefly introduced as follows. 1.3.1. Eulerian-Eulerian approach In this approach, each phase is treated as a continuum; therefore, conservation equations for mass, momentum and energy are mathematically derived and solved. The Eulerian-Eulerian approach can be formulated in different forms depending on the treatment of the interface between the involved phases. One of the fairly simple versions of these models is known as the Volume-of-Fluid (VOF) model. In this model, the interaction of immiscible fluids is modeled by solving a single set of momentum equations and an additional volume fraction equation. The interface between fluids is reconstructed based on the volume fraction distribution. This model is suitable for cases where the position of the interface of the fluids is of interest; therefore, the model tracks the volume fraction of all of the fluids throughout the domain. The VOF model is most appropriate for free-surface flows, filling, sloshing, the motion of large bubbles in a liquid, the prediction of jet breakup, and the steady or transient tracking of liquid- gas interfaces. The more complex version of the Eulerian-Eulerian models is known as the Euler or multi-fluid approach. In this approach, different phases are assumed to be interpenetrating continua; - 6 - therefore, it involves the solution of a number of sets of conservation equations, equal to the number of phases involved. The coupling between the phases is achieved with pressure and other interfacial forces. For multi-phase problems involving a granular flow (e.g. solid particles), the system of equations is closed, for example, by applying the kinetic theory for granular flows (KTGF) [4]. The Euler KTGF approach is most commonly used for simulating dense fluidized beds. 1.3.2. Eulerian-Lagrangian approach Two of the variants of the Eulerian-Lagrangian approach used for simulating multi-phase flows are known as the discrete particle model (DPM) and discrete element method (DEM). In both of these approaches, the fluid phase is treated as a continuum for which the Navier-Stokes equations are solved, while the dispersed phase is solved by tracking the trajectories of dispersed entities (such as particles, bubbles, or droplets) through the calculated flow field. In the DEM approach individual entities are considered, while representative clusters are used in the DPM. The dispersed phase can exchange momentum, mass, and energy with the fluid phase. Similar to the Eulerian-Eulerian approach, the Eulerian-Lagrangian approach can be divided into subcategories. Considering the framework of the collision model used, there are two well-known approaches: Hard-sphere and soft-sphere. In the hard-sphere approach, the trajectories of particles are determined by conserving the momentum inside the system and assuming binary and instantaneous interactions between the particles. Therefore, the occurrence of multiple particle collisions cannot be considered in this framework. One of the first hard-sphere discrete particle simulations reported in the literature, used for studying granular systems, was by Campbell and Brennen [5]. Unlike the hard-sphere approach, the soft-sphere approach allows for multiple particle interactions where the net force is obtained by adding the results of all the binary interactions. One of the first publications using the soft-sphere approach was by Cundall and Strack [6]. Both of the above approaches, combined with the solution of the gas-phase, have been successfully used for simulating fluidized beds. Among the pioneers are Hoomans et al. [7] using the hard-sphere approach and Tsuji et al. [8] using the soft-sphere approach. - 7 - 1.4. LITERATURE REVIEW 1.4.1. Inter-particle forces and coalescence modeling Size enlargement of wet particles (granulation) occurs by the collision and successful coalescence of particles coated with a liquid binder. It is observed in a large number of industries, e.g. food, catalysts, ceramics, fertilizers, pharmaceuticals, and detergents. Traditionally, these industries relied on the knowledge of practitioners; however, a significant amount of research has recently focused on the fundamental understating of granulation processes [9]. In spite of the wealth of such granulation studies, there is still a significant gap between microscale studies and macroscopic models of granulation processes. The first step in modeling the granulation processes at a micro-scale is being able to determine the collision outcome of wet particles. Improved understanding of the interactions of wet particles can lead to an improved understanding of the granulation process [9, 10]. A good understanding of the governing inter-particle forces due to liquid bridging is key to the proper description of the interaction of wet particles at a micro-scale. Seville et al. [11] reviewed all the main inter-particle forces related to fluidization and estimated that the van der Walls and electrostatic forces are orders of magnitude smaller than the capillary forces. As a result, these two forces (i.e. van der Walls and electrostatic forces) are usually ignored in the studies related to the multi-phase flow systems, especially fluidized beds. Hence the focus is often on the viscous and capillary contributions associated with liquid bridges. There are two types of forces attributed to the liquid bridges connecting wet particles: static (capillary) and dynamic (viscous). Under static conditions, surface tension and the pressure deficiency due to the curvature of the interface induce an attractive capillary force. The strength of this capillary force has been the subject of extensive experimental and theoretical studies [10, 12-18]. For instance, Pitois et al. [15] proposed a new capillary force expression that depends on both the liquid bridge volume and separation distance; in addition, they compared the results of their new expression with their experimental measurements and found excellent agreement. When particles move relative to each other (under dynamic conditions), the viscous effects also contribute. The viscous contributions have also been studied extensively [15, 17, 19-24]. Chan and Horn [19] studied the drainage of the thin liquid film between solid surfaces experimentally, and found excellent - 8 - agreement with Reynolds lubrication theory. Ennis et al. [20] studied the liquid bridge between two oscillating particles and developed force expressions that compared well with their experimental data. Matthewson [22] modified the viscous force expression by considering a finite wetting distance for the liquid bridge between a sphere and a flat surface. Pitois et al. [15] showed that viscous forces can be more accurately predicted once the liquid volume effects are included in the viscous forces expression. Moreover, the relative importance of the capillary and viscous forces associated with a liquid bridge has been investigated by a number of authors [11, 15, 20]. Despite all the progress, a comprehensive work that considers both the capillary and viscous effects and compares their relative importance in related granulation processes is still lacking. Knowing the key forces attributed to liquid bridges, the collision outcome of particles (either coalescence or rebound) can be determined. Iveson [25] provides a summary of the available coalescence models in the literature. Basically, there are two modeling approaches for describing the coalescence of wet particles in the literature: approaches with dominant capillary effects (that ignore viscous effects) and approaches that assume viscous effects are dominant (and neglect capillary effects). In models based on the first approach [26-28], the coalescence of particles depends on the balance between the rupture energy of the liquid bridges and kinetic energy of particles. Consequently, the rupture energy has been the subject of a number of studies. For instance, Simons et al. [28] derived a simple expression for the rupture energy of liquid bridges that was in reasonable agreement with the limited data available at that time. Willett et al. [29] modified the rupture energy expression by incorporating the contact angle. Pitois et al. [26] proposed a simple analytical expression for the rupture energy of the liquid bridge between two spherical particles that compared very well with experimental results for a range of bridge volumes. Rossetti et al. [27] compared their experimental values for spherical agglomeration processes with the rupture energy expressions proposed by Simons et al. [28] and Pitois et al. [26] and indicated that the latter shows better agreement. However, the assumption of dominant capillary effects is not always correct. As argued by some researchers [20, 30, 31], in some cases the viscous forces can be dominant. Hence, this approach has been the subject of numerous publications. Barnocky and Davis [32] proposed a Stokes theory, assuming dominant viscous effects and showed good agreement with their experimental findings. Ennis et al. [21] used the - 9 - Stokes theory to describe the collision behaviour of two identical wet spherical particles. The viscous Stokes theory has been used by others as well [33, 34]. According to the above review, the current available collision models in the literature are only appropriate for limiting cases when either the viscous or capillary force is clearly dominant. There are no coalescence models in the literature that account for both the capillary and viscous contributions of the liquid bridge. In addition, there is no clear indication of when viscous effects are dominant and when capillary effects are dominant. Moreover, there are some shortcomings with the current coalescence models in the literature; for instance, the thin film analysis of Ennis et al. [20], used for the viscous force expression for liquid bridges between two spheres, fails when the separation distances between particles is large or the liquid bridge volume is insufficient. Pitois et al. [15] modified the capillary and viscous forces expressions by incorporating the liquid bridge volume effects, and found excellent agreement with experimental data. 1.4.2. Size enlargement of particles in fluidized beds Wet granulation is an important process in chemical, food, mining, and pharmaceutical industries and has been the subject of numerous experimental and modeling studies. The related experimental studies mostly focus on the effect of the process variables and chemo-physical properties of the primary particles (or agglomerates) and the liquid binder on the process [35-47]. In addition, various analytical models and numerical simulations have been used to study the granulation processes. The analytical models usually attempt to predict the collision outcome of particles (i.e. whether they coalesce or rebound); these models are usually based on the size of the primary particles and the chemo-physical properties of the binder [21, 26-28, 33, 48]. The numerical simulations usually have a broader objective and attempt to analyze the entire reactor, often using population balance modeling (PBM) [35, 39, 43] or discrete element modeling (DEM) [49-53]. As stated earlier, this thesis is motivated by the agglomeration problem in fluid cokers [1-3]. Fluid cokers have been the subject of numerous studies that particularly improve our understanding of the hydrodynamics and jet-bed interactions of three-phase flow inside the - 10 - reactor [54-69]. For instance, Tafreshi et al. [54] investigated the impact of liquid feed conditions on the coker efficiency and indicated that spray characteristics strongly depend on the flow regime inside the pipe. As highlighted by Gray [1], due to the imperfect liquid-solid mixing conditions inside fluid cokers, the liquid feed (bitumen) covers only a small portion of the coke particles with a relatively thick layer. Knapper et al. [3] conducted tracer studies to understand the fundamentals of jet penetration and to quantify the efficiency of liquid feed distribution among solid particles and concluded that 1) proper atomization and using high performance sprays improves the solid-liquid contacting, and 2) the overall contacting efficiency is low and the sprayed liquid droplets wet a portion of solid particles; these two observations eventually lead to agglomeration of coke particles in fluid cokers. In order to enhance jet-bed interactions and mixing inside fluid cokers, some have proposed to use the draft tube or an enhanced solid entrainment (ESE) device, which is an open-ended cylindrical tube placed shortly downstream from the nozzle and axially aligned with it [55, 59, 61, 69-72]. Some of these studies indicate the effectiveness of jet and draft tube configuration due to improved liquid-solid mixing; this could result in less agglomeration. For instance, House et al. [59] demonstrated that an ESE-nozzle configuration increases the number of contacted particles, provides a more uniform liquid-solid mixture, and improves liquid-solid contacting. In addition, McMillan et al. [61] showed that the ESE device provides better liquid-solid contacting due to the rapid and uniform distribution of sprayed droplets among solid particles. Bitumen consists of heavy hydrocarbons that are converted into lighter ones due to the cracking and coking reactions occurring inside the reactor. Physical properties and cracking reactions of bitumen at coking conditions have been the topic of a number of studies [73-77]. Gray et al. [73] experimentally investigated the effect of initial liquid film thickness on the reactor performance and concluded that limiting the film thickness to around 20 µm leads to optimum operation. According to Aminu et al. [75], due to the chemical reactions, bitumen viscosity could increase to an order of 104 mPa.s at reaction temperatures of 400-530 o C . Also, due to evaporation of lighter hydrocarbons as a result of the same chemical reactions, bitumen thickness continuously reduces with the reaction time [74]. The agglomeration of particles depends on both bitumen viscosity and thickness; the increase in bitumen viscosity aids agglomeration, while the reduction - 11 - in bitumen thickness counteracts. Hence, as proposed by Gray [1], in the intermediate time during coking reactions, a maximum agglomeration tendency should be observed. Despite the increased number of fluid coker studies within the past few years, there is practically no theoretical model that directly addresses the agglomeration of bitumen-coated coke particles in fluid cokers. In addition, the analogy between agglomeration in fluid cokers and wet granulation processes noted by Gray [1] allows established methods to be used for the analysis of fluid cokers. An example of this could be the application of coalescence models to study the collision of coke particles coated with bitumen layers; such models, indeed, require consideration of the time-dependant variations of bitumen viscosity and thickness due to chemical reactions. 1.4.3. Rupture and liquid distribution of liquid bridges The statics and dynamics of liquid bridges between solid surfaces, either flat or curved, is related to various multi-phase industrial processes such as wet granulation, coating, and printing; therefore, they have been studied for a long time. These investigations started with the pioneering works of Haines [78] and Fisher [12], studying the liquid bridges between two spheres on contact. The research on the liquid bridges was gradually developed in time and it was around the late 70s when the dynamics of liquid bridges generated more interest [79]. Mehrotra and Sastry [80] conducted a comprehensive chronological review of the related studies up to 1980 and showed that significant improvement could still be made in the area. Following this review, a substantial amount of research has been conducted; these studies have one of the following major goals in common: (1) improving the force and energy expressions associated with liquid bridges, e.g. [10, 15, 19-21, 21, 26-28, 31, 81-83], and (2) obtaining a better understanding of the stability and modeling the evolution of liquid bridges [81, 84-91]. The former was discussed in sub-section 1.4.1 and the latter is reviewed in this section. Due to the evolution of a liquid bridge, it can rupture and distribute liquid between the solid surfaces. Rupture of liquid bridges is an interesting phenomenon that has been studied for many years. These studies were initially on static liquid bridges [92-94]. The dynamics and rupture of liquid bridges was not extensively studied until the work of Fowle et al. [79]. Following their work, similar studies were performed on oscillating and stretching liquid bridges [95-99]. There - 12 - are two main parameters used to quantify the rupture and post-rupture evolution of liquid bridges: “rupture distance” and “liquid transfer fraction.” The rupture distance is the distance at which the liquid bridge breaks, and the liquid transfer fraction is a measure of how the liquid is distributed between the two solid surfaces after rupture. Several modeling approaches have been used to determine rupture; the related models are divided into the following three categories in this work: (1) Models using the Young-Laplace equation, (2) simplified models (basically relying on the geometry of the problem), and (3) models based on the governing Navier-Stokes equations. The Young-Laplace equation relates the pressure deficiency and liquid-solid surface tension to the mean curvature of the liquid bridge. In order to determine the rupture distance, the Young- Laplace equation is solved while minimizing the system free energy. However, the equation does not lend itself to an analytical solution except for a few special cases [14]. Therefore, efforts have been made to solve it numerically, e.g. Erle [100] and De Bisschop and Rigole [84]. Mazzone et al. [85] compared the theory developed by De Bisschop and Rigole with experimental data and argued that the theory underestimates the rupture distance; consequently, they developed a modified theory that considered gravitational effects and showed better agreement. Lian et al. [81] studied the liquid bridge between spherical bodies and determined the rupture distance by minimizing the system free energy. They showed that for identical particles with contact angles up to 40 degrees, the rupture distance can be written as a simple function of the liquid bridge volume and contact angle. Later, similar expressions were proposed by Mikami et al. [82] for the case of two identical particles, and for a particle and a wall. Using the Young- Laplace equation for determining the rupture distance usually requires constant and zero liquid- solid contact angles [87], which is one of the main limitations of this model. An alternative approach is approximating the liquid bridge shape with an a priori assumed profile and then solving the resulting system of equations; we refer to this approach as “the simplified model.” This model, contrary to the previous approach, is easily extendable to particles with varying sizes and/or with different contact angles. The simplified model has variations and has been used often in the past few years. Dai and Lu [86] used a circular liquid bridge profile and - 13 - employed the theory of liquid cylinder stability to predict the rupture distance. They used a correction factor to match their predictions with experimental data, which somewhat reduces the generality and reliability of their model. Pepin et al. [87] employed both the toroidal and parabolic approximations to model liquid bridges in the presence of wetting hysteresis. Pepin et al. [87] showed that even though the two profiles have different mathematical expressions, their approximations agree well with their own experimental data. They also pointed out that the toroidal approximation is slightly less accurate and more difficult to use; that is because one of the parameters of the toroidal approximation approaches infinity when the bridge shape converts from convex to concave, causing numerical difficulty. Shi and McCarthy [89] also used the parabolic approximation and determined the rupture distance and liquid distribution of liquid bridges between solid spherical particles with different sizes and different contact angles; finally, they implemented their findings in a numerical tool to study rotary drum spray-coating systems. Shi and McCarthy [89] also argued that the toroidal approximation is restricted to particles that have the same contact angles. In conclusion, the parabolic approximation is more suitable than the toroidal approximation for modeling the rupture and liquid distribution of liquid bridges. The simplified models have limitations. First, they assume that the bridge is stable throughout its evolution. Second, the capillary, viscous, inertial and gravitational effects are not explicitly accounted for in these models. Third, they use a pre-assumed profile for the gas-liquid interface. The above limitations of the simplified model need further investigation. Another model is the numerical solution of the governing Navier-Stokes equations. This model has been used much more often than the previous two models in the past two decades. This is mainly because of the significant improvement in computational power and numerical algorithms that have made numerical simulations popular. In order to simulate the dynamic evolution of stretching liquid bridges between two circular disks, Zhang et al. [99] used the slender jet approximation, similar to Eggers & Dupont [101] and Papageorgiou [102]. Zhang et al. [99] simplified the governing Navier-Stokes equations and derived two one-dimensional time- dependant differential equations for the axial velocity and bridge profile. They then performed extensive numerical investigations of continuous uniaxial stretching of liquid bridges between two circular disks. They concluded that the liquid bridge evolution is strongly influenced by the relative magnitude of the viscous, inertial, and gravitational effects compared to the surface - 14 - tension effects. Zhang [103] solved the full governing Navier-Stokes equations by using a finite- difference formulation and volume-of-fluid (VOF) method in an Eulerian mesh and predicted the evolution of a drop exiting a fine capillary tube. He accounted for the inertial, viscous, gravitational, and surface tension effects and showed excellent agreement with their experimental measurements. In order to model the micro-gravure-offset printing, Huang et al. [104] performed Computational Fluid Dynamics (CFD) studies of the liquid transfer between two parallel separating flat plates and between a trapezoidal cavity and a moving flat plate. They found that the liquid transfer fraction deviates from the theoretical value (obtained for the quasi- static conditions) as the capillary number increases. Dodds et al. [105] used a Galerkin finite element method to study the stretching liquid bridges between two flat plates with slipping contact lines. They showed that at fixed capillary numbers, the amount of liquid transferred to the plate with the lower contact angle increases as the contact angle of the other plate is increased, and for fixed contact angles, the liquid transfer improves (approaching uniform distribution) as the capillary number is increased. According to the above literature review, the simplified model has limitations that are yet to be quantified. In addition, the available numerical simulations of the governing Navier-Stokes equation for stretching liquid bridges in the literature are mainly focused on flat plates, and there is no evidence of any work for liquid bridges between solid particles. 1.4.4. Discrete Element Method simulations of wet fluidized beds Fluidized beds are extensively used in numerous industrial applications; they are complex multi- phase systems that are not fully understood yet. The basic step toward improved understanding of their operation could be obtaining more knowledge regarding the transport phenomena inside them. Traditionally, a large number of experimental studies have been conducted to investigate the transport phenomena in fluidized beds. However, there are practical limitations with experimentation of fluidized beds; for instance, the three-dimensional transport phenomena in fluidized beds are not readily visible. In addition, probe measurements are very difficult to achieve, and if accomplished, they locally disturb the flow inside the system and affect the final outcome of the experimental measurements. Because of the above, numerical simulations of fluidized beds have gained considerable attention in the past few years. In particular, numerical - 15 - simulations allow researchers to gain insight into the (three-dimensional) transport phenomena inside fluidized beds without disturbing them. A good numerical tool for obtaining detailed information of the transport phenomena inside fluidized beds is Discrete Element Method (DEM) simulation. With the aid of these simulations, valuable information of the entire system including the local velocities, volume fraction of solid phase, etc. could be obtained. Such detailed information, if not impossible, is extremely hard to obtain through experimentation. Furthermore, the results of the detailed small-scale DEM simulations (accounting for fluid- particle and particle-particle interactions) can be used to develop closure models (such as collision models) in order to improve the prediction capabilities of the current numerical tools used for simulating fluidized beds on a larger scale. On the other hand, DEM simulations are very expensive computationally; hence, due to the current limitations of computational systems, they are often used to simulate fluidized beds only on a small scale. The operation of conventional gas-solid fluidized beds becomes more complicated when a liquid phase is added to the system. This is primarily due to the formation of liquid bridges and their resulting cohesive forces between particles. Such cohesive forces can influence the transport phenomena inside the system; for instance, they could lead to agglomeration and inefficient operation of the fluidized beds. There are a number of publications that experimentally study wet fluidized beds [106-108]. For instance, McLaughlin and Rhodes [108] showed that the incremental addition of liquid to a fluidized bed of Geldart group B particles results in the transition of bed behaviour resembling a bed of group A particles, or even group C particles. Recently, DEM simulations have been used to study the effect of liquid on different industrial applications such as fluidized bed coaters, fluidized bed granulators, and tumbling granulators. For instance, granulation in fluidized beds has been the subject of a number of DEM investigations [49, 51, 52]; however, these studies do not effectively account for the cohesive effects of the added liquid in their simulations. However, there are DEM simulations in which the cohesive liquid effects are considered in further details. For example, Mikami et al. [82] performed 2-D DEM simulations of wet fluidized beds by including the capillary effects of liquid bridges. With the aid of their numerical tool, they were able to study the fluidization behaviour of both wet and dry particles; for example, they showed that due to the cohesive forces, the minimum fluidization velocity for wet particles is greater than that of dry particles. - 16 - Muguruma et al. [109] numerically studied the flow of wet mono-sized spherical particles in a centrifugal tumbling granulator by including the capillary effects of liquid in their DEM simulations. They showed that due to the liquid effects, the motion of particles in the systems is greatly affected; they also compared the calculated components of particle velocities with experimental measurements and found good agreement. The work of Muguruma et al. [109], similar to Mikami et al., ignores the viscous effects of the liquid bridge and only considers the capillary effects. Yang and Hsiau [110] performed DEM simulations to investigate the effect of liquid bridge volume and liquid viscosity on the motion of particles in a 2-D vibrated bed. According to them, the energy is mainly dissipated due to the viscous effects, interparticle friction, and inelasticity of particles rather than the capillary effect of liquid bridges. Kunal et al. [111] employed DEM simulations to investigate the effects of liquid-induced cohesion on gas- solid flows and showed a clear transition from free-flowing to cohesive behaviour of the system at a distinct value of Granular Capillary number, i.e. the ratio of the capillary force to the drag force. Recently, van Buijtenen et al. [53] conducted 3-D DEM simulations in order to study wet spout fluidized beds. For this purpose, they assumed that the coefficient of restitution varies depending on the moisture content of particles only. According to their work, due to the liquid effects, the coefficient of restitution decreases; this leads to the formation of a dense region inside the bed and consequently, the passage of gas through the system in the form of bubbles. Considering a coefficient of restitution only dependant on the moisture content is not realistic; in practice, parameters such as the liquid viscosity, liquid surface tension, and particle collision velocities matter as well. Despite the increased number of DEM studies on the operation of fluidized beds in the presence of liquid, our understanding of these systems is still incomplete and there is room for improvement. For instance, with the aid of DEM simulations, the effect of the liquid coating properties such as viscosity and surface tension on the bed hydrodynamics could be numerically investigated. - 17 - 1.5. THESIS OBJECTIVES As highlighted earlier, this thesis aims at obtaining an improved understanding of the interaction of wet particles and wet granulation phenomena. The main objectives of this work can be summarized as follows: - to model the collision of wet particles (coated with thin liquid films); - to obtain a better understanding of the agglomeration phenomenon in fluid cokers; - to model the rupture of stretching liquid bridges between solid particles; - to study the effect of liquid coating on the operation of wet fluidized beds. 1.5.1. Approach In order to achieve the above objectives, the problem is first explored at a micro-scale level for which only two wet particles are considered: • In order to describe the collision of wet particles, a coalescence model is developed. The model is written in the form of a wet coefficient of restitution and accounts for the collision velocity, liquid viscosity, liquid surface tension, liquid bridge volume, etc. • Then, in order to study the collision outcome of bitumen-coated particles in fluid cokers, the above collision model is further modified. The time- and temperature-dependant viscosity and thickness variations of the liquid coating are incorporated. The above modifications are incorporated in order to account for the actual operation conditions of fluid cokers in the model. • In order to study the rupture of stretching liquid bridges, two different approaches are used; i) a simplified mathematical model and ii) numerical simulations of the governing Navier-Stokes equations are used. The first has a number of simplifying assumptions, and thus has limitations. The second uses numerical simulations of the governing equations to investigate the rupture phenomenon of stretching liquid bridges. The numerical simulations are based on the VOF approach and account for the viscous, surface tension, inertial, gravitational, and contact angle effects on the rupture distance and liquid distribution. Comparison of the two methods allows assessing the simplified model and its limitations. - 18 - Finally, the interaction of multiple wet particles in a fluidized bed is considered: • An-open source numerical tool (MFIX) capable of simulating the gas-solid flow in fluidized beds is modified. MFIX is based on a soft-sphere DEM approach and is modified by implementing the wet coefficient of restitution proposed in this thesis. Finally, the modified code is used to perform DEM simulations in order to study the effect of liquid coating on the hydrodynamics of fluidized beds. 1.6. THESIS LAYOUT In Chapter 2, a coalescence model that incorporates both capillary and viscous contributions of the liquid binder, along with the liquid bridge volume effects, is developed. The model is based on an overall coefficient of restitution that is determined with the aid of the approximate analytical values of the maximum possible energy dissipation and a critical value that depends on the total initial kinetic energy of particles. With this model, the collision outcome of two spherical particles coated with a thin liquid film is predicted. In Chapter 3, a mathematical model is proposed to determine the agglomeration tendency of bitumen-coated coke particles in fluid cokers. The model calculates a theoretical critical velocity that depends on key parameters such as particle size, bitumen viscosity, and bitumen thickness; it also accounts for the temperature- and reaction-dependant variations of the bitumen thickness and viscosity. In Chapter 4, a simplified mathematical model is proposed and numerical simulations of the governing Navier-Stokes equations are used to predict the shape evolution, rupture distance, and liquid distribution of stretching pendular liquid bridges between two equal-sized spherical solid particles. In the simplified mathematical model, the bridge shape is approximated with a parabola, and it is assumed that the surface tension effects dominate the viscous, inertial, and gravitational effects. For the numerical simulations, a commercial Computational Fluid Dynamics (CFD) software package – FLUENT – is used. In the first part of this chapter, the - 19 - applicability range of the above two models is explored; the variations of the rupture distance with respect to the capillary, Weber and Bond number are investigated. Also, the shape of the liquid bridge interface between two equal-sized particles with identical and non-identical contact angles is studied. In the second part of this chapter, numerical simulations of the governing Navier-Stokes equations are used to study the influence of the contact angle on the rupture and liquid distribution. In the last part of this chapter, numerical simulations are used to investigate the influence of gravity on the rupture of stretching liquid bridges. In Chapter 5, the proposed wet coefficient of restitution, proposed in Chapter 2, is implemented in an open-source numerical tool - MFIX. The modified numerical tool is used to perform DEM simulations of a fluidized bed consisting of mono-sized solid spherical particles pre-coated with identical liquid coatings. Then, DEM simulations are performed to study the effect of the coating thickness and viscosity on the operation of a small-scale fluidized bed. Finally, the studies conducted in this thesis are summarized in Chapter 6, and the key conclusions drawn from them are presented. In addition, the main contributions of this thesis are highlighted and some recommendations for future work are provided. - 20 - Chapter 2. A COALESCENCE MODEL FOR BINARY COLLISION OF WET PARTICLES1 2.1. INTRODUCTION Wet granulation is a well-known process employed in various industries, e.g. chemical, food, catalysts, ceramics, fertilizers, pharmaceuticals, and detergents. It occurs when wet primary particles attach together and form larger entities, known as granules. Traditionally, the operation of wet granulating systems relied on the knowledge of practitioners; however, a significant amount of research has recently focused on the fundamental understating of granulation processes [9]. In order to obtain a good understanding of wet granulation processes, it is important to explain the problem at a micro-scale. The first step is to determine the governing forces between two wet particles and predict their collision outcome. According to the literature review in Chapter 2, the related forces (i.e. capillary and viscous) and collision models (also known as coalescence models) have been investigated in numerous publications. However, the available models in the literature are only appropriate for the limiting cases when one of the above forces is dominant. There is no coalescence model in the literature that accounts for both the capillary and viscous contributions of the liquid bridge. In this chapter, a coalescence model that incorporates both capillary and viscous contributions is developed. The idea is to describe the collision outcome with an overall coefficient of restitution. This coefficient of restitution depends on the maximum possible energy dissipation and a critical value. The maximum possible energy loss accounts for the capillary and viscous effects of the liquid bridge along with the inelasticity of dry particles, and the critical values represent the total kinetic energy of a system consisting of two particles. 1 A version of this chapter has been published. Darabi, P., Pougatch, K., Salcudean, M., and Grecov, D. (2009) A Novel Coalescence Model for Binary Collision of Identical Wet Particles. Chemical Engineering Science. 64: 1868- 1876. - 21 - 2.2. MODEL DEVELOPMENT 2.2.1. Problem statement Consider two identical solid spherical particles of radius R approaching from opposing directions, as schematically shown in Figure 2-1.a. For our problem, particles are approaching with initial normal transitional velocities of u0 and the tangential and angular velocities are neglected. The particles are assumed to have an average height of asperities with size ha and are coated with a well-dispersed Newtonian liquid with an initial height of h0. Once the liquid layers touch each other, a liquid bridge connecting the particles is formed instantaneously, as shown in Figure 2-1.b. A larger scale of the bridge is shown in Figure 2-1.c in which the two-dimensional system of cylindrical coordinates (r,z), the separation distance H(r) = 2D + r2/R, the minimum separation distance (2D), and the wetting distance (rw) are shown. Note that due to the geometric considerations, the wetting distance cannot exceed the particle radius. Figure 2-1. Approach of two moving spherical particles and the formation of the liquid bridge. - 22 - The collision is divided into three stages: approach, contact, and rebound. Each stage covers a certain range of distance during the entire collision. In the approach stage (see Figure 2-1), particles move from an initial half-separation distance equal to the initial height of the liquid layer (D = h0) until they contact. At contact, the half-separation distance is equal to the average height of the surface asperities of the particle (D = ha). The contact stage is self-explanatory. In the rebound stage, particles move away from an initial half-separation distance equal to the average height of the surface asperities of the particle (D = ha) to a certain final half-separation distance. Two choices are explored for the final half-separation distance: i) equal to the initial liquid height, or ii) equal to half of the rupture distance (0.5hrup). The rupture distance is the critical distance at which a liquid bridge ruptures. As depicted on the left-hand side of Figure 2-2, before liquid bridge rupture, the distance between particles is smaller than the rupture distance and liquid is still connecting the particles; however, as illustrated on the right-hand side of this figure, the separation distance is equal to the rupture distance, and the bridge is about to rupture. Further investigations related to the liquid bridge rupture are presented in Chapter 4. The liquid bridge formation is assumed to be an instantaneous phenomenon; therefore, the treatment of the liquid bridge effects in the approach stage is similar to that of the rebound stage. Figure 2-2. Rebound of two moving spherical particles (a) and rupture of liquid bridge (b). - 23 - 2.2.2. Assumptions In order to model the above collision problem, the following assumptions are considered. Many of these assumptions apply to the rest of this thesis as well. • Solid particles are spherical, non-porous, and of uniform size • Coating liquid is Newtonian and perfectly wets the particles • Gravitational forces, buoyancy forces, and gas viscous effects are neglected • The collision is normal; hence, the rotational and tangential components of the velocity profiles of the particles are neglected • There is no turbulence. 2.2.3. General equation of motion Applying Newton’s second law of motion and considering the viscous and capillary forces of the liquid bridge, the general equation of motion for one particle can be written as follows: vis cap duF F m dt + = r r r ( 2-1) where visF r is the viscous force, capF r is the capillary force, m is the particle mass, and ur is vector of the particle velocity. Using the thin liquid film (or lubrication) approximation, i.e., assuming the local radii of curvature are large compared to the minimum distance between the surfaces and the Reynolds number is low, the viscous force vector can be written as follows [19, 24]: 2 23 2vis v uF R X D piµ= − r r ( 2-2) where µ is the liquid viscosity, R is the particle radius, D is the half-separation distance, ur is the particle velocity vector ( / uu dD dt e= r r ), vX is the correction factor due to the liquid bridge volume effects ( 2 2( ) 1 /vX D D D a= − + ), 2a is a constant ( 2 / 2ba V Rpi= ), and bV is the liquid - 24 - bridge volume. The derivation of the above viscous force expression is provided in Appendix A. Writing the viscous force in the form of Eq. ( 2-2), its direction will be correctly determined in the opposite direction of the particle motion. Neglecting body forces (such as gravity and buoyancy), the resulting capillary force can be assumed to be the summation of the surface tension effects acting on the liquid-solid boundary and the reduced hydrostatic pressure effects, written as follows [12, 15]: 2 coscap v nF R X epi σ θ= r r ( 2-3) where σ is the liquid surface tension, θ is the solid-liquid contact angle, vX is the correction factor (as above), and ne r is the unit vector normal to the particle surface along the axis of collision. Further detail on the above capillary force expression is provided in Appendix B. Writing the capillary force vector in the form of Eq. ( 2-3), its direction will be correctly determined, indicating the attraction between the two particles. Approximating the liquid bridge connecting the two spheres with a cylindrical shape, its volume can be calculated as follows: 2 2 0 2 ( ) ( ) 4 2 wr b wV rH r dr R H r D pi pi = = − ∫ ( 2-4) 2.2.4. Methodology As mentioned earlier in this work, the collision outcome of particles is determined based on an overall coefficient of restitution. This coefficient depends on the maximum possible energy dissipation due to the viscous and capillary effects and inelasticity of dry particles. The losses are determined using an approximate analytical approach. An approximate approach is used here because the governing equation of motion as shown in Eq. ( 2-1) cannot be solved analytically. In order to determine the approximate analytical expressions for losses, first the capillary contribution is neglected and analytical velocity expressions are derived. These expressions lead to analytical expressions for the viscous works and contact loss. Independently, analytical - 25 - expressions for the capillary works are derived. Finally, with the aid of the superposition principle, the above analytical values are added, resulting in an approximate analytical expression for maximum possible energy dissipation. The superposition principle can be used as the problem is being solved at very small Reynolds numbers (equal to one or smaller); hence, the liquid acceleration term of the Navier-Stokes equation, leading to non-linearity, can be neglected. The general equation of motion is solved with a numerical method (coded in Matlab with a Runge-Kutta solver) and the losses due to each contribution are determined accordingly. These values are referred to as numerical results. 2.2.5. Proposed coalescence model The overall coefficient of restitution of two spherical particles colliding axially, approaching with equal initial velocities and without rotation can be written in the following form: 0 ' fu e u = − ( 2-5) where 'e is the overall coefficient of restitution, fu is the projection of final velocity of the particle after impact, and 0u is the projection of the initial velocity of the object before impact. The coefficient of restitution has values between 0-1, where 1 represents a perfectly elastic collision, and 0 represents a perfectly inelastic collision. The law of the conservation of energy for the system consisting of the above objects can be written as: 2 2 0 f tmu mu L= + ( 2-6) where tL represents the total energy losses. Substituting the after collision velocity expressions in Eq. ( 2-6) and solving the resulting equation for the overall coefficient of restitution, the following relation can be obtained: - 26 - 2 0 ' 1 tLe mu = − ( 2-7) The detailed derivation of this equation is provided in Appendix C. The above formulation is exact, provided the losses are correctly determined. However, to determine the actual losses, it is necessary to resolve the differential equation in order to find the exact location where the particle stops. This work though is not focussed on establishing the actual stopping point of the particles and the actual losses. Therefore, a simplified approach is proposed to evaluate the total energy dissipation. It is assumed that particles move up to a final separation distance. Then, based on this assumption, the maximum possible dissipation losses are calculated; as the name implies, the value of these losses is always greater than or equal to the actual losses. The proposed formulation results in a model with a simpler calculation and evaluation procedure. Finally, using the maximum possible energy dissipation, the overall coefficient of restitution can be written as follows: ,max 2 ,max 02 0 ' 1 ,t t L e L mu mu = − < 2 ,max 0' 0, te L mu= ≥ ( 2-8) where m represents the mass of particles, 0u represents the scalar value of the initial particle velocity, and ,maxtL represents the maximum possible energy dissipation. According to the above model, the maximum possible energy dissipation is compared to a critical value, which is equal to the summation of the kinetic energy of the two particles. In addition, although one can determine the overall coefficient of restitution for the rebound outcome, this is not addressed in this chapter. The application of the proposed coefficient of restitution including the rebound outcome is addressed in Chapter 5. Also, it is important to highlight that the above coefficient of restitution differs from the commonly used dry coefficient of restitution; the proposed expression addresses a more general problem, i.e. the collision of wet particles, and incorporates both the effects of the liquid coating and the inelasticity of dry particles. - 27 - The maximum possible energy dissipation for a system consisting of two particles can be written as follows: ,max , ,2( )t vis t cap t cL W W L= + + ( 2-9) where ,vis tW is the maximum possible viscous work, ,cap tW is the maximum possible capillary work, and cL is the contact loss. Maximum possible work expressions are used for the viscous and capillary contributions, as the maximum possible energy dissipation is being determined. 2.2.5.1. Viscous work In order to obtain the contact loss and viscous work for the approach and rebound stages, the velocity expression for each stage is required. The required velocity expressions are obtained by integrating Eq. ( 2-1) while only accounting for the viscous forces and ignoring the capillary force. The velocity cannot change direction under the action of the viscous force; thus, if the velocity reaches zero, the motion stops. This point is included in the derived expressions of velocity for the approach, contact, and rebound stages, as follows: 0 0 0 0 ( ) ( )11 ln , ln( ) ( ) ( )0, ln ( ) v v a v f h f h u St St f D f D u f hSt f D − > = ≤ r r ( 2-10) 0 0 0 0 ( ) ( )11 ln , ln( ) ( ) ( )0, ln ( ) v v a a c v a f h f h u St St f h f h u f hSt f h − > = ≤ r r ( 2-11) 0 0 0 0 ( ) ( )ln , ln( ) ( ) ( )0, ln ( ) c v v v a c a r v c a u eu St uf D f DSt St f h u eu f h u u f DSt eu f h − > = ≤ r r ( 2-12) - 28 - where au , cu and ru present the velocity expression for the approach, contact (right before contact) and rebound stage, respectively, 0u r is the initial velocity vector, vSt is the viscous Stokes number ( 202 / 3vSt mu Rpiµ= ), 0h is the initial liquid height, ah represents the average height of surface asperities of the particle (or surface roughness), ( )f D can be written as ( )22 2 2 2( ) /f D D D a D D a= + + + , 2 / 2ba V Rpi= , and bV is the liquid bridge volume. The above expressions are written for the bottom particle in Figure 2-1 and result in velocity directions consistent with the direction of motion of that particle in each stage. Similar velocity expressions with a negative sign can be written for the other particle. The derivation of the above velocity expressions is provided in Appendix D. The viscous work for the approach stage can be written as follows: 0 ,1 .( ) a h vis vis u h W F e dD= ∫ r r ( 2-13) Using the approach velocity expression, provided in Eq. ( 2-10), the viscous work can be written as follows: 20 0 ,1 ( ) ( )ln ln( ) 2 ( )vis a v a f h f hAW A f h St f h = − ( 2-14) where A is a constant ( 2 03 / 2A R upiµ= ). Similarly, the maximum possible viscous work for the rebound stage can be written as follows: 1 2 1 1 ,2 0 ( ) ( ) .( ) ln ln 2 ( ) ( ) a h c vis vis u v a ah Aeuf h f hAW F e dD St f h u f h − = = − ∫ r r ( 2-15) where 1h is the final half-separation distance that can either be the initial height of the liquid layer, as chosen by Ennis et al. ( 1 0h h= ), or half of the rupture distance ( 1 0.5 ruph h= ). Note that when calculating the second term on the right-hand side of Eq. ( 2-15), the contact velocity (right after contact) is in the opposite direction of the initial velocity; therefore, an additional negative sign will be automatically considered for this term. It is also important to highlight that the - 29 - viscous work expressions derived from our integral treatment, as shown in Eqs. ( 2-14) and ( 2-15), cannot be combined into one expression by modifying the integral bounds of either of them. This is because of the significant difference between the velocity expressions for the approach and rebound stage, see equations ( 2-10) and ( 2-12), respectively. The rupture distance of the liquid bridge connecting two particles can be determined from the following equation, as proposed by Lian [81]: ( ) 1/31 0.5rup bh Vθ= + ( 2-16) Summing the above work expressions, the maximum possible viscous work for one particle can be expressed as follows: , ,1 ,2vis t vis visW W W= + ( 2-17) 2.2.5.2. Capillary work The capillary work for the approach stage can be written as follows: ( )0 02 2,1 .( ) 2 cos a a h h cap cap u hh W F e dD R D D api σ θ= = − − +∫ r r ( 2-18) Similarly, the capillary work for the rebound stage can be written as: ( )0.5 0.52 2,2 .( ) 2 cosrup rup a a h h cap cap u hh W F e dD R D D api σ θ= = − +∫ r r ( 2-19) The maximum possible capillary work for one particle is the summation of the capillary works for the approach and rebound stages and can be determined as follows: ( ) 0 0.5 2 2 , ,1 ,2 2 cos ruph cap t cap cap h W W W R D D api σ θ= + = − + ( 2-20) Equation ( 2-20) can also be derived by integrating the capillary force, with respect to the half- separation distance, between the limits of the initial liquid height and the half-rupture distance. Using the derived capillary work expression as in Eq. ( 2-20), the rupture energy for two particles - 30 - connected via a liquid bridge, when separating from contact until the liquid bridge ruptures, can be written as follows: ( ) 0.52 2 0 4 cos ruph rupW R D D api σ θ= − + ( 2-21) The above equation can also be written in the following non-dimensional form [26, 27]: 2 2 22 cos b brup rup rup V VW h hpi θ pi pi = − + + % % % %% ( 2-22) where rupW% is the dimensionless rupture energy ( 2/rup rupW W Rσ=% ), ruph% is the dimensionless rupture distance ( ( ) 1/ 31 / 2rup bh Vθ= +% % ), and bV% is the dimensionless bridge volume ( 3/b bV V R=% ). Note that Eq. ( 2-20) is written for one particle only, while Eqs. ( 2-21) and ( 2-22) are written for a system consisting of two particles. 2.2.5.3. Contact Loss The contact loss for a dry inelastic collision can be calculated with the following equation: 2 21 (1 ) 2c c L m e u= − ( 2-23) where m is the particle mass, e is the dry coefficient of restitution and cu is the particle contact velocity. 2.2.6. Viscous limiting case For collision with dominant viscous effects, Ennis et al. [21] proposed the following model in order to determine the collision outcome of two identical particles: * 011 lnv a hSt e h = + ( 2-24) 08 9 p v u R St ρ µ = ( 2-25) - 31 - where *vSt is the critical Stokes number and vSt is the particle Stokes number. According to the model, whenever the viscous Stokes number exceeds the critical value, particles rebound and coalescence occurs whenever the reverse is true. The derivation of the critical Stokes number is provided in Appendix E. Assuming a critical condition ( *v vSt St= ), the following can be written: 0 08 11 ln 9 p a u R h e h ρ µ = + ( 2-26) Multiplying the above by a factor of 2 03 / 2R upi µ , the following is obtained: 2 2 0 0 0 3 11 ln 2 a h mu R u e h pi µ = + ( 2-27) where the left-hand side term represents the summation of the initial kinetic energy of the two particles and the right-hand side term represents the maximum possible energy dissipation when ignoring the capillary effects. Ignoring capillary effects and considering the same conditions as Ennis et al. [21], i.e. i) the Stokes number equal to the critical Stokes number, and ii) final half-separation distance equal to the initial liquid height ( 1 0h h= ), the maximum possible energy dissipation, as in Eq. ( 2-9), is reduced to the following expression: 2 0 ,max 0 ( )3 11 ln 2 ( )t a f hL R u e f hpi µ = + ( 2-28) Eq. ( 2-28) resembles the maximum possible energy dissipation expression directly derived from the viscous Stokes theory, see Eq. ( 2-27). The difference is that the maximum possible energy dissipation expression as in Eq. ( 2-28) contains an additional term that incorporates the liquid bridge volume effects. In addition, by comparing Eq. ( 2-28) with Eq. ( 2-24), two other critical Stokes numbers can be derived from our proposed model: * 0( )11 ln ( )v a f hSt e f h = + ( 2-29) - 32 - * 0 (0.5 )( ) 1ln ln( ) ( ) rup v a a f hf hSt f h e f h = + ( 2-30) In order to derive equations ( 2-29) and ( 2-30), the final half-stopping distance is equal to the initial liquid height ( 1 0h h= ) and half-rupture distance ( 1 0.5 ruph h= ), respectively. 2.2.7. Capillary limiting case Neglecting the viscous effects, the contact loss and expression for the maximum possible capillary work can be written as follows: 2 2 0 1 (1 ) 2c L m e u= − ( 2-31) ( ) 0 0.5 2 2 , 2 cos ruph cap t h W R D D api σ θ= − + ( 2-32) Using the above expressions, the maximum possible energy dissipation can be determined with the aid of Eq. ( 2-9). Working within the proposed framework, the coalescence of particles with dominant capillary effects can be determined with the aid of equations ( 2-31) and ( 2-32). Note that in Eq. ( 2-32), the capillary work for both the approach and rebound stages is considered; however, the available coalescence models available in the literature are based on the rupture energy [26, 27, 112], meaning the work of the approach stage is excluded. 2.3. RESULTS AND DISCUSSION 2.3.1. Comparison of analytical and numerical results It is stressed that in order to derive the velocity expressions for different collision stages (i.e. approach, contact, and rebound), the capillary contributions were neglected. This is a convenient simplification used for formulating the analytical velocity expressions and needs to be further examined. This simplification is evaluated by comparing the analytical and critical values of the critical velocity. The critical velocity is the approach velocity for which the maximum possible energy dissipation is equal to the total initial kinetic energy of the two particles i.e. 2 ,max 0tL mu= . - 33 - The analytical critical velocity is determined with the aid of the proposed model. In order to obtain the numerical critical velocity, the equation of motion, as in Eq. ( 2-1), is solved in Matlab using the “ode15s” solver. Since the approach velocity is an initial condition required for performing the calculations and is missing, the equation of motion is solved iteratively. This is achieved by repeatedly guessing the approach velocity until the critical condition is satisfied, i.e. the final velocity is equal to zero at the final half-separation distance (h1). The parameters used for the calculations performed in this chapter are summarized in Table 2-1 and remain constant unless otherwise specified. Particles are considered to be of glass, and all the liquid properties except for viscosity are assumed to be that of water at 25 oC. In order to study the problem for a range of capillary numbers, liquid viscosity is assumed to vary within the range of 0.1-104 mPa.s. In order to analyze the results, the capillary number is used. This number represents the ratio of the viscous effects to the capillary effects, and can be written as follows: 0uCa µ σ = ( 2-33) Table 2-1. Summary of the values of parameters. Parameter Value Unit Particle radius (R) 1 mm Particle density ( ρ ) 2500 Kg/m3 Liquid viscosity ( µ ) 0.1-104 mPa.s Liquid surface tension (σ ) 72 /mN m Contact angle (θ ) 0 radian Average height of surface asperities (ha) 0.005 mm Dry coefficient of restitution (e) 0.9 - Initial height of liquid layer (h0) 0.01 mm Figure 2-3 shows the variations of the analytical and numerical critical velocity with respect to the capillary number. The figure indicates a very good agreement for the range of capillary numbers explored. Likewise, an excellent agreement was found between the analytical and - 34 - numerical maximum possible energy dissipation and total initial kinetic energy. In conclusion, the proposed approximate model in this work is appropriate for determining the collision outcome of wet particles. In other words, the numerical solution to the equation of motion is no longer required in order to determine the collision outcome of wet particles. This is the major contribution of this work: the numerical solution of the equation of motion (for the binary collision of two wet particles) is replaced with an approximate analytical model. Figure 2-3. Variations of the analytical and numerical critical velocities with respect to the capillary number. Figure 2-4 demonstrates the variations of the critical velocity determined from three different coalescence models: the general coalescence model (including all the effects), capillary limiting case (ignoring viscous dissipation), and viscous limiting case (ignoring capillary effects). According to this figure, at large capillary numbers (equal to or greater than 100), the critical velocities obtained with the viscous limiting model match with those of the general model. On the other hand, a good agreement is found between the results of the capillary limiting model and the general model at very small capillary numbers (equal to or smaller than 0.01). In conclusion, at either of the above limiting conditions, the limiting versions of the proposed model, instead of its general form, can be used to determine the collision outcome of wet particles. However, at - 35 - intermediate capillary numbers, the general model (accounting for both viscous and capillary effects) needs to be used. Figure 2-4. Comparison between values of the critical velocity obtained with the general coalescence model, viscous limiting model, and capillary limited model. 2.3.2. Coalescence with dominant viscous effects Figure 2-5 depicts the comparison between the critical velocities obtained with the general coalescence model, the viscous limiting model, and three versions of the Stokes model for a range of capillary numbers. According to this figure, the general model, the viscous limiting model, and Eq. ( 2-30) match very well at very large capillary numbers. This is expected because: 1) at large capillary numbers, the viscous effects are dominant and the general model and the viscous limiting case are going to be the same, and 2) Eq. ( 2-30), similar to the general model and the viscous limiting model, incorporates the liquid bridge volume effects and assumes that the half-rupture distance is equal to the final separation distance. On the other hand, at small capillary numbers, the Stokes models and viscous limiting case underestimate the critical velocity. This is because all these models exclude the capillary contributions which are not negligible at the low range of capillary numbers. - 36 - Figure 2-5. Comparison between the values of the critical velocity obtained with the general coalescence model, the viscous limiting model, and different versions of the Stokes model. 2.3.3. Coalescence with dominant capillary effects Figure 2-6 shows the variations of the non-dimensional values of the capillary work, rupture energy, and contact loss with respect to the dry coefficient of restitution. All the values shown in the y-axis of this figure vary between 0 and 1. This is because they have been non- dimensionalized by the total energy dissipation with dominant capillary effects, i.e. the summation of the capillary work and the contact loss. According to this figure, the contact loss is negligible for near to elastic collisions, here collisions with a dry coefficient of restitution equal to or greater than 0.95. In other words, according to the figure, for collisions with a dry coefficient of restitution less than 0.95, the contact loss is not negligible with respect to the capillary work; hence, for such cases, the contact loss is required to be considered. In addition, the figure shows that the rupture energy is only slightly greater than the capillary work (for a system of two identical particles); hence, whenever required, these two terms can be used interchangeably with reasonable accuracy. In summary, for capillary-dominated collisions with negligible contact losses, consideration of the rupture energy will suffice. This confirms the - 37 - suitability of the coalescence models based only on the rupture energy for capillary-dominated collisions. Figure 2-6. Variations of the capillary work, rupture energy, and contact loss with respect to the dry coefficient of restitution. 2.3.4. Comparison with experimental data 2.3.4.1. Viscous dominated Barnocky and Davis [32] set up a ball-dropping apparatus in order to study the collision of particles on a wet flat surface. Their objective was to determine the minimum dropping height such that the ball would no longer stick to the plate; this dropping height corresponds to the critical (approach) velocity used in this chapter. In order to analyze their results, they reduced the liquid film height on the plate to a characteristic height equal to two-thirds of the actual height (due to a geometric consideration). Using this modified liquid thickness and the observed dropping heights (later used in calculating the critical velocity), they determined the experimental values of the critical Stokes number. Finally, they compared these experimental values with theory and found good agreement for Stokes numbers less than four. - 38 - In this section, the same experimental data are used to validate the three critical Stokes models that are proposed in this chapter. Note that in all of these models the actual liquid film thickness ( 0h ), rather than the reduced value is used for calculating the critical Stokes number. Eq. ( 2-24) represents the critical Stokes number in its original form, Eq. ( 2-29) represents the critical Stokes number considering the liquid bridge volume effects when the final seperation distance is equal to the initial height of liquid coating, i.e. hf = h0, and Eq. ( 2-30) represents the critical Stokes number considering the liquid bridge volume effects when the final seperation distance is equal to the half rupture distance, i.e. hf = 0.5hrup. The comparisons between these models and experimental data for three different ball sizes are shown in Figure 2-7(a-c). Note that in the x- axis of this figure, h0 and ha represent the liquid film thickness and plate surface roughness, respectively. In all of these figures, the experimental values are marked with squares, the calculated critical Stokes numbers of Barnocky and Davis (using the modified liquid film thickness) are marked with circles, and the calculated critical values using the modified versions of the Stokes model, i.e. equations. ( 2-24), ( 2-29), and ( 2-30) are marked with plus, star, and cross signs, respectively. According to these figures, the original Stokes theory, as in Eq. ( 2-24) and marked with plus signs, tends to overestimate the experimental measurements. However, the modified versions of the Stokes model that incorporate the liquid bridge volume effects, i.e. equations ( 2-29) and ( 2-30), show better agreement with the experimental data. In conclusion, including liquid bridge volume effects improves the prediction capability of the Stokes model. However, due to the limited comparisons made in this figure, it is not easy to determine which of these two modified versions provides a better prediction when compared to experimental findings. In order to be able to find the better model between equations ( 2-29) and ( 2-30), more experimental data and comparison is required. - 39 - (a) (b) - 40 - (c) Figure 2-7. Variations of the critical Stokes number with the respect to the initial height of liquid layer or three different particles sizes: a) R = 0.8 mm, b) R = 1.0 mm, and c) R = 2.0 mm. 2.3.4.2. Capillary dominated Pitois et al. [26] experimentally determined the rupture energy between two moving spheres connected with a liquid bridge. They compared the experimental values with their proposed expression for the rupture energy, Eq. ( 2-22), and found excellent agreement, as depicted in Figure 2-8. As mentioned earlier, for collisions with a coefficient of restitution close to unity, the contact loss is negligible in comparison to the rupture energy. Considering this explanation and the similarity between our capillary work and the rupture energy expression of Pitois et al. [26], Figure 2-8 can be regarded as a validation for the suitability of our proposed coalesence model for near to perfect elastic collisions with dominant capillary effects. - 41 - Figure 2-8. Comparison between experimental and analytical values of rupture energy. 2.4. SUMMARY AND CONCLUSIONS In this chapter, the collision behaviour of two particles coated with a thin liquid film was studied. A coalescence model was developed that considers the capillary, viscous, and liquid bridge volume effects. The model is based on an overall coefficient of restitution that was derived from energy conservation before and after collision. According to the proposed model, particles coalesce whenever the maximum possible energy dissipation is equal to or greater than the total initial kinetic energy of two particles. In order to determine the maximum possible energy dissipation, an approximate analytical treatment was chosen. The maximum possible energy dissipation was assumed to be the summation of the analytical expressions for the contact loss, capillary, and viscous works. Each loss/work term was derived separately. In order to derive the viscous and contact losses, the capillary effects were neglected; this makes the proposed analytical approach approximate. In order to verify the approximate analytical approach, the general equation of motion was solved with a numerical method. The numerical results were compared with the approximate analytical results, and excellent agreement was found. This means that the proposed model replaces the numerical solution of the governing equation of motion. - 42 - Two different limiting versions of the proposed model were formulated and compared with its general form. The results indicated the capillary contributions dominate at very small capillary numbers, while the viscous contributions are dominant at very large capillary numbers. In between, both of these contributions could be important, hence the general model incorporating both effects needs to be used. Moreover, for viscous-dominated cases, different modified versions of the Stokes theory were derived (from the proposed coalescence model) and were compared with the experimental data. The comparison indicated a good agreement between the available experimental data in the literature and the modified Stokes models that incorporate the liquid bridge volume effects. Also, the proposed model at capillary-dominated conditions was compared with experimental data and great agreement was found. For capillary limiting cases, the contact loss was compared with the capillary work in a range of dry coefficient of restitution. It was indicated that the contact loss is negligible for dry coefficients of restitution equal to or greater than 0.95, while it becomes significant when this coefficient is smaller than 0.95. It was also shown that the rupture energy is only slightly larger than the maximum possible capillary work. In summary, for capillary-dominated collisions and dry coefficients of restitution close to unity, the consideration of the rupture energy will suffice. In conclusion, the proposed model successfully predicts the collision outcome of two spherical particles coated with a thin liquid film. This work can be considered as one step towards narrowing the gap between microscopic and macroscopic modeling of wet granulation processes. Ultimately, the proposed model can be implemented in particulate phase motion simulations to describe the operation of granulation systems in the presence of a liquid. In Chapter 3, a version of the proposed model is used to determine the collision outcome of bitumen-coated coke particles. Also, in Chapter 5, the proposed coefficient of restitution is implemented in an open- source code to study the operation of fluidized beds in the presence of liquid coating. - 43 - Chapter 3. AGGLOMERATION OF BITUMEN-COATED COKE PARTICLES IN FLUID COKERS2 3.1. INTRODUCTION Fluid cokers are prone to agglomeration [1-3]. This occurs when bitumen-coated coke particles collide and stick together, and finally very large granules consisting of up to several hundred coke particles are formed. Agglomeration in fluid cokers is undesirable and causes reduced process efficiency. In order to improve the operation of fluid cokers, it is useful to learn more about conditions leading to agglomeration. As highlighted in the introduction chapter, fluid cokers have been the subject of numerous studies that: 1) aim at improving our understanding of the hydrodynamics and jet-bed interactions of the three-phase flow inside them, 2) study the effect of the draft tube on the jet- bed interactions and mixing inside them, and 3) investigate the physical properties and cracking reactions of bitumen at coking conditions. Despite the increased number of studies on the fluid cokers, there is no theoretical model that directly addresses the agglomeration of bitumen-coated coke particles. In this chapter, a simplified theoretical model is proposed to determine the collision outcome of bitumen-coated coke particles; the model incorporates the temperature- and time-dependant variations of the viscosity and thickness of bitumen that occur due to chemical reactions. The proposed model is then used to improve our understanding of the agglomeration process inside fluid cokers. 3.2. FLUID COKER The operation of fluid cokers has been described in the introduction chapter, and is presented once again here. The fluid coker is one of two key components of a fluid coking system. It is a large steam-fluidized bed containing hot coke particles inside which the bitumen is injected 2 A version of this chapter has been published. Darabi, P., Pougatch, K., Salcudean, M., and Grecov, D. (2010) Agglomeration of Bitumen-coated Coke Particles in Fluid Cokers. International Journal of Chemical Reactor Engineering. 8, A122. - 44 - through a large number of nozzles using steam-assisted atomization. The atomized bitumen droplets contact and coat the coke particles. Due to the high temperature of the coke particles, the heavy hydrocarbon molecules of bituminous liquid thermally crack and are converted into lighter products. As a result of cracking reactions, the coke particles lose energy, cool down, and grow in size (as bitumen remaining on the particle surface is converted into coke). In order to be reheated, the large, relatively cold coke particles are directed into a large fluidized bed, known as the burner, and once sufficiently hot again, the coke particles are circulated back into the fluid coker. 3.2.1. Flow regions within a fluid coker In this work, the flow inside the fluid cokers is classified into three regions: jet, intermediate, and well-mixed. These flow regions are suggested based on their difference in temperature and inter- particle collision velocities. Each jet introduces the bitumen at relatively low temperatures (usually around 350-400 oC) and very high velocities. In the well-mixed region, due to the contact with the hot coke particles, the bitumen temperature increases to around 530 oC, while due to the substantial solid volume fraction, the inter-particle collision velocities are significantly lower than in the jet region. In order to draw boundaries between the described regions, efforts are made to quantify the bed conditions within each region. Unfortunately, there are no available velocity or temperature measurements of the flow inside the fluid cokers, which is not surprising considering the process conditions. In such a situation, quantification is based on the common understanding of the process. As the fluidized bed is very well agitated, it seems that the range of inter-particle collision velocities from 0.2 to 0.3 m/s is reasonable; these collision velocities correspond to granular temperatures from 0.0057 to 0.0127 m2/s2, according to 24 9cuθ pi= [113]. This range is assigned to the well-mixed region. For the jet region, collision velocities greater than 1.5 m/s (corresponding to granular temperatures greater than 0.318 m2/s2) are assumed. Anything in- between falls into the intermediate region. We realize that the transition within the reactor is continuous; however, for the ease of analysis, the above regions are demarcated considering the proposed inter-particle collision velocities and reaction temperatures, as summarized in Table 3-1. - 45 - Table 3-1. Characteristic inter-particle collision velocity and temperature of the flow regions within fluid cokers. Flow region Characteristic inter-particle collision velocity, m/s Temperature, oC Jet >1.5 350-400 Intermediate 0.3-1.5 400-530 Well-mixed 0.2-0.3 530 3.3. MATHEMATICAL MODEL 3.3.1. Problem statement Consider two identical spherical particles of radius (R) approaching with equal initial transitional velocities (u0), as schematically shown in Figure 3-1.a. Each particle is assumed to be coated with a well-dispersed Newtonian liquid with a thickness of h . Once the liquid layers touch each other, a liquid bridge that connects the particles is instantaneously formed, as schematically shown in Figure 3-1.b. A larger scale of the bridge is shown in Figure 3-1.c, where the two- dimensional system of cylindrical coordinates (r,z), the separation distance changing with radial distance ( )H r , the minimum separation distance ( 2D ), and the wetting distance ( wr ) are shown. Figure 3-1. Collision of two moving spherical particles and the formation of the liquid bridge. - 46 - 3.3.2. Proposed model Previously, the above collision problem was studied in Chapter 2 and it was shown that the collision outcome depends on the initial kinetic energy of particles and the total losses due to the capillary and viscous effects of the liquid bridge along with the inelastic contact losses. The results were further analyzed and it was concluded that for dynamic collisions with large capillary numbers ( 0 /Ca u µ σ= ), the viscous effects dominate the capillary effects. It was also shown that by neglecting the capillary effects, the model reduces to a modified version of the Stokes criterion, introduced earlier by Ennis et al. [21]. In this work, it is assumed that due to the large bitumen viscosities, the viscous effects are dominant. Therefore, the modified Stokes model is used; see Chapter 2 for more details. The model can be mathematically written as follows: * 1 ( )1 ln ( )a f hSt e f h = + ( 3-1) 08 9 pu RSt ρ µ = ( 3-2) where *St is the critical Stokes number, e is the dry coefficient of restitution, ( )22 2 2 2( ) /f D D D a D D a= + + + , 2 / 2ba V Rpi= , bV is the liquid bridge volume, h is the liquid film thickness, ah represents the average height of surface asperities (or surface roughness) of the particle, St is the particle Stokes number, pρ is the particle density, 0u is the initial particle velocity, R is the particle radius, and µ is the liquid viscosity. The Stokes model is used as follows: Both Stokes numbers are calculated and compared. The particle velocity resulting by equating the above Stokes numbers is referred to as the “critical velocity”; this velocity is the demarcation between the rebound and coalescence outcomes. The critical velocity ( *u ) can be expressed as follows: *9 * 8 p u St R µ ρ = ( 3-3) - 47 - In order to determine the variations of the bitumen viscosity with the reaction time, Aminu et al. [75] performed experimental measurements at three reaction temperatures, i.e. T = 400, 503, & 530 oC. In this work, the linear interpolation of the experimental results of Aminu et al. [75] is used and the following expressions for the time-dependant bitumen viscosity at the above reaction temperatures are proposed: 26542 1 at T=400 ( ) 30131 1 at T=503 21012 1 at T=530 t C t t C t C µ + = + + o o o ( 3-4) where ( )tµ is the viscosity in mPa.s and t is the normalized reaction time. The normalized reaction time is defined as the ratio of the actual reaction time over the dry-out time. According to Aminu et al. [75], the dry-out time is the reaction time at which the bitumen film exhibits no further adhesive tendency; therefore, at reaction times greater than the dry-out time, the liquid effects are neglected. The dry-out time of a bitumen coating with a thickness of 20µm is 240, 24, and 14.4 seconds at reaction temperatures of T = 400, 503, and 530 oC, respectively [75]. According to the literature [1], bitumen is not uniformly distributed among the coke particles inside fluid cokers, hence the thickness of bitumen coating could vary from one particle to another. Considering the described mechanism of bitumen distribution among coke particles by Gray [1], the average bitumen coating on coke particles is around 22.5µm . In this work, this value is used as the average and it is assumed that the thickness of the initial liquid layer on the coke particles varies between 15 and 30µm . Due to the cracking and devolatilisation reactions, the thickness of bitumen coating continuously reduces with the reaction time. Gray et al. [74] performed experiments on thin films of bitumen (around 20µm ) at coking temperatures of 457, 503, and 530 °C. They also showed that the weight fraction of the feed decreases with the reaction time, meaning that the bitumen thickness decreases in time. Using the data from the above article, Aminu et al. [75] proposed an empirical expression to calculate the bitumen thickness as a function of reaction time. The expression shows the proper trend, i.e. the thickness of bitumen coating decreases in time. However, it does not fully account for the influence of the initial thickness on the coke yield (the remaining - 48 - material on the particle surface). According to Radmanesh et al. [77], the coke yield is dependant on the initial thickness of bitumen coating, while the proposed expression by Aminu et al. [75] does not exhibit such dependency. Therefore, a more general equation consistent with the experimental observations is required. Such an expression should comply with the following conditions: 1. The thickness of bitumen coating decreases with the reaction time. 2. The final coke yield is dependant on the initial thickness of bitumen coating. According to Radmanesh et al. [77], the final coke yield increases with the film thickness. Here, it is assumed that the coke yields for bitumen coatings of 10, 20, and 30 µm are 8, 10, and 12%, respectively. These values are in fair agreement with the findings in the literature [73, 74, 77]. 3. The final thickness of bitumen coating (at the dry-out time) is proportional to the initial thickness of the bitumen coating and the final coke yield. 4. According to Radmanesh et al. [77], for thin films (20µm and smaller), the coke yield is independent of the reaction temperature. Although film coatings up to 30µm are considered in this work, the above assumption is extended to all thicknesses as a first approximation. Therefore, in the proposed expression, the final coke yield is independent of the reaction temperature. Considering the above assumptions and maintaining the decaying trend of the original expression introduced by Aminu et al. [75], the following expression for the variations of the thickness of bitumen coating as a function of the normalized reaction time is proposed: 1/ 0.15340.1534 0.3 0( ) 0.4708 0.04h t h t t = − + ( 3-5) where h(t) is the time-dependant thickness of bitumen coating, h0 is the initial thickness of bitumen coating at time zero, and t is the normalized reaction time. Figure 3-2 shows the variations of the thickness of bitumen coating with respect to the normalized reaction time for initial bitumen coatings of 10, 20, and 30µm . The thickness of bitumen coating obtained with - 49 - the above equation is used to determine the liquid bridge volume and the critical velocity at each normalized reaction time. Figure 3-2. Variations of the thickness of bitumen coating with the normalized reaction time. Given the geometry in Figure 3-1.c, the separation distance between the two particles can be written as follows: 2 ( ) 2 rH r D R = + ( 3-6) Assuming a cylindrical shape for the liquid bridge, its volume can be approximated as follows: 2 2 0 2 ( ) ( ) 4 2 wr b wV rH r dr R H r D pi pi = = − ∫ ( 3-7) In order to obtain the liquid bridge volume, H(rw) and D must be determined. The liquid height at the wetting distance is assumed to be the summation of the thickness of the two liquid layers, i.e. ( ) 2 ( )wH r h t= , where h(t) is the time-dependant thickness of bitumen coating. Also, assuming that particles are in contact when the liquid bridge is formed, the minimum separation distance - 50 - will be equal to zero, i.e. D = 0. Given the above, Eq. ( 3-7) can be simplified to the following expression: 22 ( ( ))bV R h tpi= ( 3-8) It is assumed that: 1) due to the short collision times between particle pairs, the liquid bridge volume remains constant during each collision. This is because the collision time of particles is in the order of milliseconds while the required time for a significant change of bitumen properties is in the order of seconds; 2) the variations of the bitumen properties (e.g. viscosity and thickness) with the reaction time during each collision are negligible; 3) the dry coefficient of restitution of the coke particles is equal to 0.9; and 4) particle surface roughness increases with reaction time, and this increase is proportional to half of the thickness of the bitumen remaining on the particles after cracking. One should note that the latter assumption is used to account for the effect of the initial thickness of the bitumen coating on the final roughness of particles. This is an approximation that is proposed for this first time in this work and requires further investigation in future. Therefore, assuming an initial surface roughness of 6.5µm for coke particles (before reactions), the surface roughness of the particles will be 7.2, 7.5, and 8.3µm after cracking reactions. These surface roughness values correspond to the initial bitumen coatings of 15, 20, and 30µm , respectively. The values of the key parameters are summarized in Table 3-2. Table 3-2. Properties of the wet coke particles. Parameter Value Unit Radius of coke particles (R) 75 µm Density of coke particles ( pρ ) 1500 kg/m3 Initial thickness of bitumen coatings (h0) 15-30 µm Surface roughness of coke particles (ha) 6.5-8.3 µm Dry coefficient of restitution (e) 0.9 - 3.3.3. Solution methodology First, using equations ( 3-4) and ( 3-5), the thickness of bitumen coating and viscosity at the given normalized reaction time and reaction temperature are determined, respectively. Given the - 51 - thickness of bitumen coating and the size of coke particles, the liquid bridge volume is calculated with the aid of Eq. ( 3-8). Then, using equations ( 3-1) and ( 3-3), the critical Stokes number and the critical theoretical collision velocity (u*) are determined, respectively. Next, plots of the variations of the critical velocity for different initial bitumen coatings at different reaction temperatures are obtained. Knowing the reaction temperature, the corresponding flow region at the desired temperature is determined, for which a characteristic inter-particle collision velocity is provided in Table 3-1. Comparing this inter-particle collision velocity with the critical velocity calculated, obtained with Eq. ( 3-3), the collision outcome of the particles at that flow region can be determined: particles coalesce (and form an agglomerate) if the theoretical critical velocity is greater than the characteristic inter-particles velocity; otherwise, they rebound. The growth of the agglomerates through the process is not explicitly modeled in this work. Instead, using the peak critical velocity (for bitumen coatings with an initial thickness of 30 µm ) and considering the characteristic inter-particle collision velocity and the size of primary particles, the agglomeration tendency inside fluid cokers is quantified. For this matter, the maximum possible agglomerate size and the number of particles consisting agglomerates are estimated. This is performed for the well-mixed and intermediate regions that are prone to agglomeration. 3.4. RESULTS AND DISCUSSION 3.4.1. Variations of the critical velocity with the reaction time Knowing the variations of the viscosity and thickness of bitumen coating with the reaction time at each coking temperature, the variations of the critical velocity for the binary collision of bitumen-coated coke particles are determined. The results for the three coking temperatures of T = 400, 503, and 530 oC are shown in Figure 3-3, Figure 3-4, and Figure 3-5, respectively. These figures have the following in common: • As each reaction proceeds, the value of the critical velocity initially increases, reaches a maximum at a certain reaction time, and decreases afterwards. The existence of a peak is in - 52 - agreement with Gray [1] who predicted the maximum agglomeration tendency at an intermediate reaction time. • The peak value increases as the thickness of bitumen coating increases. This can be explained by considering the relationship between the critical velocity and the thickness of liquid film. As the thickness of bitumen film increases, the critical Stokes number increases, leading to a larger critical velocity. • The normalized reaction times corresponding to the peak critical velocities increase as the initial thickness of bitumen coatings increases. The peak critical velocities occur at the normalized reaction times of 0.01, 0.03, and 0.07 corresponding to the initial bitumen coating thicknesses of 15, 20, and 30µm , respectively. • For the normalized reaction times beyond a certain value, no critical velocity is calculated. This is because at such reaction times, the thickness of bitumen coatings is smaller than the particle surface roughness; therefore, the liquid coating is not effective. These normalized reaction times increase as the initial thickness of bitumen coating increases and are 0.02, 0.08, and 0.22 for the initial bitumen coatings of 15, 20, and 30µm , respectively. Comparing the values of the above theoretical critical velocities with the characteristic inter- particle collision velocities of each flow region, provided in Table 3-1, the following conclusions can be reached: - Theoretically, agglomeration of coke particles in the jet region (with T = 400 oC) could be expected. According to Figure 3-3, the largest critical velocity for this temperature is equal to around 2 m/s, corresponding to an initial bitumen coating thickness of 30µm and can be greater than the inter-particle collision velocity in the jet region (1.5 m/s); hence agglomeration is likely to occur. On the other hand, due to the very high velocities and short penetration lengths of jets, the actual residence time of coke particles in the jet regions is very low (around a second or less). In addition, the highest critical velocity at T = 400 oC occurs at a normalized reaction time of 0.07; considering the dry-out time of the thin bitumen coatings, i.e. around hundreds of seconds, - 53 - the above normalized time is equivalent to ten seconds or more. Finally, considering the occurrence time of the highest peak critical velocity and the residence time of coke particles in the jet region, it can be concluded that practically agglomeration of coke particles in the jet regions is not of any concern. - Agglomeration of coke particles in the intermediate region (with T = 503 oC) is very likely to occur. According to Figure 3-4, given the initial thickness of bitumen coatings of 20 and 30µm , the peak critical velocity is substantially greater than the average inter-particle collision velocities in the intermediate region. This conclusion is valid provided that the residence time of coke particles in the intermediate region is sufficient such that the peak critical velocities occur. This, in fact, is very likely to happen as fluid cokers could be modeled as continuously stirred- tank reactors (CSTR) for which the residence times are practically very high (and in theory, up to infinity). - Similar comments can be made for the well-mixed region (T = 530 oC). Figure 3-3. Variations of the critical velocity with the normalized reaction time for T = 400 oC. - 54 - Figure 3-4. Variations of the critical velocity with the normalized reaction time for T = 503 oC. Figure 3-5. Variations of the critical velocity with the normalized reaction time for T = 530 oC. - 55 - The above figures indicate that for identical initial bitumen coatings, the peak value for T = 503 oC is greater than that of T = 530 oC. For instance, for initial bitumen coating of 0 30µmh = , the peak critical velocity for T = 503 oC is * 2.50 m/su = , while for T = 530 oC it is around u* = 1.75 m/s. This can be translated as less agglomeration tendency at high coking temperatures. The difference between the critical velocities comes from the difference between the bitumen viscosities at these two temperatures at t = 0.07; the bitumen viscosity at T = 503 oC and t = 0.07 is equal to 2110 mPa.sµ = , while for T = 530 oC it is 1472mPa.sµ = . 3.4.2. Estimating the maximum possible agglomerate size Previously, it was shown that for the initial bitumen coatings of 30µm in the well-mixed region (T = 530 oC) and the intermediate region (T = 503 oC), the theoretical critical velocity could be significantly greater than the actual characteristic inter-particle collision velocities within the corresponding region, as provided in Table 3-1. This means that under such conditions agglomerates that are substantially larger than the primary coke particles could be formed. Accordingly to Eq. ( 3-3), the critical collision velocity and particle size are inversely correlated, meaning that if the particle size decreases, the critical velocity increases. With the aid of this equation, the maximum possible agglomerate size ( aggR ) can be determined to be: * p agg c u R R u = ( 3-9) where aggR is the radius of agglomerate, * pu is the peak theoretical critical velocity, cu is the characteristic inter-particle velocity of the flow region, and R is the radius of the primary coke particles. For instance, according to Figure 3-4, the peak critical velocity in the intermediate region with 0 30µmh = is u * = 2.50 m/s. Comparing this critical velocity with an average characteristic velocity of 0.5m/scu = in the intermediate region, agglomerates with a maximum size of 375µmaggR = could be obtained. Considering a porosity of 0.5 for the agglomerates, such agglomerates contain approximately 62 coke particles. Similarly, according to Figure 3-5, the peak critical velocity in the well-mixed region is around u* = 1.75 m/s. Comparing this critical velocity with a characteristic inter-particle velocity of 0.2 m/scu = for this region, agglomerates - 56 - of maximum size 656aggR = µm, approximately containing 335 coke particles, are obtained. These values are in agreement with the practical plant operations, for which agglomerates consisting of up to a few hundred coke particles are observed. The above characteristic inter- particle collision velocities are average values that are used to obtain an estimate of the size of agglomerates and shed some light on the agglomeration process inside the fluid coker reactor. In practice, the inter-particle collision velocities vary within a range inside the entire reactor; this eventually results in the formation of agglomerates with a size distribution. 3.5. SUMMARY AND CONCLUSIONS A simplified mathematical model that incorporates the time- and temperature-dependant variations of the bitumen viscosity and coating thickness was proposed to determine the collision outcome of the bitumen-coated coke particles inside the fluid cokers. The model calculates a critical velocity that is a function of a number of parameters such as the bitumen viscosity, the thickness of bitumen coating, the reaction time, and the size and surface roughness of the coke particles. The results indicated a peak for the variations of the critical velocity with respect to the normalized reaction time. It was found that agglomeration cannot occur in the jet regions, while for the intermediate and well-mixed regions agglomeration is very likely to occur. Also, the results show that the agglomeration tendency at T = 503 oC is greater than that at T = 530 oC. Further, the analysis of results suggests that large agglomerates containing up to a few hundred coke particles could be formed inside fluid cokers, which is in agreement with the actual plant observation. - 57 - Chapter 4. MODELING THE RUPTURE OF STRETCHING LIQUID BRIDGES3 4.1. INTRODUCTION The evolution of liquid bridges between solid surfaces is related to several industrial processes such as wet granulation, coating, and printing; hence, they have been the subject of numerous investigations. As highlighted in the introduction chapter, one of the main objectives of these investigations is to obtain a better understanding of the stability and evolution of liquid bridges [81, 84-89]. This has led researchers to develop mathematical models for the rupture of liquid bridges. There are two types of models that have been commonly used in the past few years: 1) a mathematical model based on a geometrical representation of liquid bridges; we refer to this group of models as the “simplified model,” and 2) numerical solution of the governing equations. As explained earlier in the introduction chapter, the simplified model has limitations; for instance, the liquid bridge shape is approximated with an a priori assumed profile and it does not allow for studying the effect of liquid properties (such as viscosity, surface tension, and density). On the other hand, the numerical simulations of the governing equations propose a more general modeling framework as they account for the surface tension, viscosity, and density of the bridging liquid; in addition, numerical simulations do not approximate the gas-liquid interface with a profile. In this chapter, a simplified model and numerical simulations of the Navier-Stokes equations are used to study the rupture of liquid bridges between two equal-sized solid spherical particles. This chapter can be divided into three parts. First, the simplified model and numerical simulations are developed and used to study a number of problems; both models are compared and the limitations of the simplified model are identified. The above part is mainly focused on the liquid bridges between particles with identical liquid-solid contact angles. In the second part, the numerical simulations of the governing equations are used to study stretching liquid bridges between particles with different liquid-solid contact angles. In the third part, the numerical simulations of the Navier-Stokes equations are used to study the effect of gravity on the rupture of stretching liquid bridges. 3 A version of this chapter has been published. Darabi, P., Li, T., Pougatch, K., Salcudean, M., and Grecov, D. (2010) Modeling the Evolution and Rupture of Stretching Pendular Liquid Bridges. Chemical Engineering Science. 65: 4472-4483. - 58 - 4.2. PROBLEM STATEMENT Figure 4-1 shows schematically a liquid bridge connecting two solid spherical particles. As indicated in the figure, the particles are separating with a constant relative velocity, v, as a result of which the liquid bridge is continuously stretching. There are two main parameters used for the analysis of the results in this chapter: the “rupture distance” and the “liquid transfer fraction.” The rupture distance (S*) is the distance at which the liquid bridge breaks, and the liquid transfer fraction (TF) is a measure of how the liquid is distributed between the two solid surfaces after rupture. Both of these quantities are quantified and discussed for a number of rupture problems in this chapter. Figure 4-1. A stretching liquid bridge between two identical solid spherical particles. 4.3. NON-DIMENSIONAL ANALYSIS Suppose that the rupture distance is a function of eight significant parameters, according to the following equation: * ( , , , , , , , )bS f R v V gρ θ σ µ= ( 4-1) where S* is the rupture distance (i.e. critical separation distance), R is the particle radius, v is the separation velocity, ρ is the liquid density, Vb is the liquid bridge volume,θ is the liquid-solid contact angle, σ is the liquid surface tension, µ is the liquid viscosity, and g is the gravitational - 59 - acceleration. Selecting{ , , }R v ρ as the scaling parameters, it can be shown with the Buckingham pi theorem that the original functional relationship, as in Eq. ( 4-1), can be written in the following equivalent form: * 3( , , , , )b VS f We Ca Bo R R θ= % ( 4-2) According to Eq. ( 4-2), the dimensionless rupture distance ( * /cS S R= ) is a function of the dimensionless liquid bridge volume ( 3/bV V R= ), contact angle ( θ ), the Weber number ( 2v /We R ρ σ= ), the capillary number ( /Ca vµ σ= ), and the Bond number ( 2 /Bo gR ρ σ= ). The latter three dimensionless numbers, i.e. We, Ca, and Bo, respectively, represent the ratio of the inertial, viscous, and gravitational effects over the surface tension effects. 4.4. SIMPLIFIED MODEL Figure 4-2 shows the liquid bridge representation as used in the simplified model, at rupture instance. For the simplified model, it is assumed that: a) the bridge is axisymmetric with respect to the x-axis; thus, only its upper section is considered in the model, and b) particles are smooth and non-porous. The key parameters of the model (as marked in the figure) are as follows: the spheres’ radii (R), the length of the liquid bridge (d), the critical separation distance, i.e. the rupture distance (Sc), the contact angles (θ ), and the spherical cap (h). In addition, parameters related to the left-hand side particle are marked with index i, and those for the right-hand side particle are marked with index j. - 60 - Figure 4-2. Liquid bridge between two spherical particles, shown at the instant of rupture. 4.4.1. Model development The proposed simplified model in this chapter is a combination of the system of equations developed by Shi and McCarthy [89] and the closure criterion by Pepin et al. [87]. First, the equations representing the parabolic shape of the bridge are presented; then, the closure criterion for determining the rupture distance is explained. Assuming a parabolic shape for the gas-liquid interface of the bridge, its profile can be approximated with the following second-order polynomial: 2y ax bx c= + + ( 4-3) where a, b and c are unknown. As depicted in Figure 4-2, the three phases (i.e. gas, liquid, and solid) intersect at two points with the coordinates: {0, (0)}y and { , ( )}d y d . Considering the problem geometry as given in Figure 4-2, the y-coordinate of these two contact points can be related to the cap heights (h) and particle radii (R) as follows: 2 2(0) ( )i i iy R R h= − − ( 4-4) - 61 - 2 2( ) ( )j j jy d R R h= − − ( 4-5) The cap heights can be related to the liquid bridge length and rupture distance with the following equation: c i jd S h h= + + ( 4-6) Considering the problem geometry as given in Figure 4-2, the contact angles (θ ) can be related to the liquid bridge profile as follows: (0) arctan( '(0)) arcsin 2i i yy R piθ = + − ( 4-7) ( ) arctan( '( )) arcsin 2j j y dy d R piθ = − − ( 4-8) The liquid bridge volume (V) is assumed to be constant throughout the liquid bridge evolution and is determined by: 2 , ,0 ( ). ( )d cap i cap jV y x dx V Vpi= − +∫ ( 4-9) where the volume of spherical caps can be written as follows: ( )2 2, 3 (0)6 icap i i hV y hpi= + ( 4-10) ( )2 2, 3 ( )6 jcap j j h V y d h pi = + ( 4-11) In summary, in order to determine the liquid profile and its rupture distance, the system of equations consisting of ( 4-4)-( 4-9) is solved. However, there are seven unknown variables ( , , , , , ,i j ca b c d h h S ) and only six equations available, meaning that one additional equation is required. The following section provides the closure of the above system of equations. - 62 - 4.4.1.1. Rupture point Figure 4-3 schematically shows the liquid droplet, immediately after rupture, on particle i. In this figure, notations S , L , and G stand for the solid, liquid, and gas phases, respectively. In this chapter, the assumption of Pepin et al. [87] is adopted and it is assumed that the rupture of the liquid bridge occurs when the area of the liquid bridge (before rupture) is equal to the area of the two droplets formed on the two particles (after rupture). This assumption can be written as follows: 2 0 2 ( ) 1 ( ). 2 ( )d i i j jy x y x dx u T u Tpi pi+ = +∫ & ( 4-12) Figure 4-3. The liquid droplet on particle i immediately after rupture. Eq. ( 4-12) is the final equation required to close the system of equations. In order to calculate the area of the two droplets, their radii and maximum cap heights are required. These two quantities for the left-hand side particle are determined with the following equations: 2 2 2 i i i i T y u T + = ( 4-13) ,2 2 2 26( 3 ) ( 3 )L ii i i i i i V T T y h h y pi + = + + ( 4-14) Similar equations are written for the liquid droplet on the right-hand side particle. - 63 - 4.4.2. Solution methodology In order to solve the problem with the simplified model, the system of equations consisting of Eqs. ( 4-4) - ( 4-14) is solved using the Matlab non-linear solver. The system of equations is solved given the radii of the particles, liquid-solid contact angles, and the liquid bridge volume. The final outcomes of the model are the rupture distance and the liquid transfer fraction of each particle. 4.4.3. Results and discussion All the results in this chapter represent the dimensionless rupture distance and dimensionless liquid bridge volumes; however, in order to avoid repetition, the word “dimensionless” is often omitted. For the problems solved in this section, i.e. the liquid bridge rupture between two equal- sized particles with identical liquid-solid contacts and no gravity, the liquid is uniformly distributed (i.e. TF = 0.5); hence, only the results concerning the rupture distance are presented. 4.4.3.1. Comparison with experimental data Figure 4-4 shows the comparison between the rupture distance predictions using the simplified model and the experimental data of Mason and Clark [114]. According to Mason and Clark, the measured contact angle was slightly greater than zero but the exact value was not provided; therefore, three series of predictions using three different contact angles of 0, 10 and 20 o are obtained and compared with the available experimental data. According to the figure, the predictions obtained with a zero contact angle underestimate the experimental measurements, while better agreement is found for contact angles of 10 and 20 degrees. - 64 - Figure 4-4. Comparison of rupture distances predicted using the simplified model with the experimental data of Mason and Clark (using water as the suspending medium and di-n-butyl phthalate and liquid paraffin mixture for the liquid bridge) [114]. Figure 4-5 displays the comparison between the rupture distance predictions using the simplified model and the experimental data of Mazzone et al. [85], with a measured contact angle of 10 o . According to this figure, the simplified model with 10θ = o over-predicts the measurements of Mazzone et al. This difference can be related to the gravitational effects. According to Mazzone et al. [85], the gravitational effects reduce the rupture distance; however, the gravitational effects are neglected in the simplified model and consequently, the rupture distance is over-predicted by this model. Further discussion regarding the impact of the gravitational effects on the rupture distance will be presented later in this chapter. - 65 - Figure 4-5. Comparison of rupture distances predicted using the simplified model and the experimental data of Mazzone et al. [85] (with measured liquid-solid contact angles of 10 o and using air as the suspending medium and dibutyl phthalate for the bridging liquid). 4.4.3.2. Effect of liquid bridge volume and contact angle on rupture distance Figure 4-6 shows the variations of the rupture distance with respect to the contact angle and liquid bridge volume. According to this figure, the rupture distance increases as the bridge volume and contact angle increase. These trends are consistent with the expression of Lian et al. [81], Eq. ( 4-15), and Mikami et al. [82], Eq. ( 4-16). Figure 4-7 shows the comparison between the results of Eqs. ( 4-15) and ( 4-16) and the rupture distance predictions using the simplified model. According to this figure, good agreement is found. In conclusion, the simplified model agrees well with other rupture distance expressions in the explored range of dimensionless liquid bridge volumes. 1/3(1 0.5 )cS Vθ= + ( 4-15) 0.34(0.99 0.62 )cS Vθ= + ( 4-16) - 66 - Figure 4-6. Variations of the dimensionless rupture distance with respect to the contact angle and dimensionless bridge volume. Figure 4-7. Comparison of rupture distances predicted using the simplified model with the expressions of Lian et al. [81] and Mikami et al. [82]. - 67 - 4.5. NUMERICAL SIMULATION OF THE NAVIER–STOKES EQUATIONS 4.5.1. Model development For the numerical simulations performed in this chapter, the full Navier-Stokes equations are solved with a CFD software package – FLUENT. In the current study, the volume-of-fluid (VOF) method combined with the interface reconstruction scheme is used for tracking the time- dependant evolution of the gas-liquid interface. To account for the surface tension effects, the continuum surface force (CSF) model, as explained in [115], is used. To consider the wall adhesion, the liquid-solid contact angle is used to obtain the normal vectors and surface curvature in the cells near the solid walls. The contact angle is assumed to be static and constant during each simulation. However, the value of this constant contact angle (either for the left-hand or the right-hand side particle) could vary from one simulation to another. Due to the movement of the two solid particles, the liquid bridge continuously elongates. To model the resulting stretching behaviour of the domain, the dynamic mesh model, as implemented in FLUENT, is used. With the dynamic mesh, layers of cells are continuously added to the mesh near the moving surfaces; this results in the continuous increase in the number of control volumes during each simulation. Further description of the above models and methods along with the governing equations can be found in the FLUENT documentation [116]. 4.5.1.1. Simulation set-up Figure 4-8 depicts the two-dimensional computational domain used in this study. The domain shows the gap between two equal-sized solid spherical particles that is filled with both gas and liquid phases. In this figure, the blue and red regions represent the gas and liquid phases, respectively; the left- and right-hand side curves represent the portion of each solid particle that forms a boundary of the domain. Since both particles are spherical, the problem is assumed to be axisymmetric. Therefore, only the 2D axial cross-section is used as the computational domain. The 2D axisymmetric domain represents a 3D domain, for which the gradients in the third direction (here, theta) are assumed to be zero. Each simulation is performed from a certain initial distance between the particles with a flat liquid bridge. The initial distance varies from one simulation to another, but in general it is chosen close to the rupture distance such that large computational times are avoided. As each simulation progresses, the shape of the liquid bridge - 68 - (which is initially flat) evolves. Its shape depends on the interaction of the liquid phase with the gas and solid phases. The flow is considered to be isothermal and incompressible. In addition, since the density and viscosity of the gas phase are orders of magnitude smaller than those of the liquid phase, the effect of the gas phase on the liquid phase is not significant in our numerical simulations. The solid particles are assumed to be smooth and non-porous. Figure 4-8. The computational domain used in this study. The blue and red regions represent the gas and liquid phases, respectively, and the boundary in-between represents the gas-liquid interface. Figure 4-9 shows one of the grids used in this work. The figure shows the gap between the two particles that is meshed with a structured quadrilateral grid. The computational domain is built and meshed with the geometry and mesh generation software package – GAMBIT. Figure 4-9 also includes the boundary conditions. For the walls, a no-slip boundary condition with a constant liquid-solid contact angle at the liquid-solid interfaces is used. A constant pressure boundary condition with zero gage pressure and an axisymmetric boundary condition are considered for the top and centerline (at the bottom), respectively. The generated mesh is then imported into FLUENT for further setup and simulation purposes. In this work, a segregated time-dependent, unsteady solver and the pressure-staggering option (PRESTO) scheme for pressure interpolation are used. The pressure implicit with splitting of operators (PISO) algorithm, by Issa [117], with skewness-neighbor correction is used to for the pressure-velocity coupling. For the momentum equation, the QUICK differencing scheme is used, and time integration is performed using a first-order implicit scheme. The Geo-Reconstruction scheme is used for the interface interpolation; this scheme represents the interface between the fluids using a piecewise-linear approach. - 69 - Figure 4-9. Numerical mesh and boundary conditions. This is for an intermediate grid that initially has 10800 control volumes. Table 4-1 presents a summary of the simulation conditions. In all the simulations, in order to impose symmetrical conditions, both particles are considered to move with the same velocity. The moving velocity of each particle is equal to half of the relative separation velocity, as given in Table 4-1. All these simulation conditions remain unchanged unless otherwise specified. Gravity is often neglected unless otherwise specified; whenever it is present, it acts along the line connecting the centers of the particles. Table 4-1. Summary of the simulation conditions. Parameter Value Unit Particle radius ( )R 1 mm Liquid bridge volume ( )bV -101 10× m3 Dimensionless liquid bridge volume 3( / )bV V R= 0.1 - Relative separation velocity ( )v 0.001 m/s Gravitational acceleration ( )g 9.81 m/s2 Liquid viscosity ( µ ) 0.001 mPa.s Liquid-solid contact angle ( )θ 30 degree - 70 - 4.5.2. Solution methodology The time step, the maximum number of iterations per time step, and the relaxation factors are carefully adjusted to ensure convergence of each simulation. The common values for the convergence-related parameters are summarized in Table 4-2. Table 4-2. Values for the convergence-related parameters. Parameter Value Time step (seconds) 65 10−× Maximum number of iterations per time step 100 Under-relaxation factor for pressure 0.3 Under-relaxation factor for density 1 Under-relaxation factor for body forces 1 Under-relaxation factor for momentum 0.3 Scaled convergence criteria for continuity and x- and y- momentum equations 10-3 4.5.3. Grid independence study In order to ensure grid convergence, three grids with different resolutions are considered. The initial numbers of control volumes (CV) for these grids are 5040, 10800, and 24300, respectively. Since a dynamic mesh is used, the number of CVs during each simulation gradually increases. The number of CVs at the beginning of each simulation is provided here. The grid information and predicted dimensionless rupture distances are summarized in Table 4-3. According to this table, as the grid becomes finer, the difference between the predicted rupture distances decreases, indicating grid convergence. Any of these grids can be used for our study; however, for most of the simulations in this chapter, the medium grid is used. Our choice represents a compromise between computational time and accuracy. - 71 - Table 4-3. Summary of grid independence study. Grid Average cell area (mm2) Initial number of CVs Predicted rupture distance Coarse 56.29 10−× 5040 0.638 Medium 52.94 10−× 10800 0.634 Fine 51.30 10−× 24300 0.633 4.5.4. Results and discussion 4.5.4.1. Bridge deformation and rupture Figure 4-10 illustrates the time-dependant evolution of the liquid bridge (ignoring gravitational effects). According to this figure, as the particles continuously move farther apart, the bridge continuously deforms and the neck becomes narrower. This continues until the liquid bridges no longer withstands any additional elongation and breaks. Figure 4-11 shows a similar evolution with the inclusion of the gravitational acceleration acting towards the right-hand side particle. This figure shows that due to the presence of the gravitational effects, more liquid accumulates on the right-hand side particle. - 72 - Figure 4-10. Simulated bridge evolution in time in the absence of gravitational effects: a. t=0, b. t=0.025, c. t=0.05, d. t=0.075, e. t=0.1, f. t=0.125, g. t=0.1335, h. t=0.134, i. t=0.1341, j. t=0.1342, k. t=0.1343, and l. t=0.1344 ms. - 73 - Figure 4-11. Simulated bridge evolution in time in the presence of gravitational effects: a. t=0.015, b. t=0.065, c. t=0.09, d. t=0.115, e. t=0.117, and f. t=0.1175 ms. 4.5.4.2. Comparison with experimental data Figure 4-14 compares the rupture distance prediction using the numerical simulations and the experimental data of Mason and Clark [114]. The simulation results presented in this figure are for very low capillary and Weber numbers, i.e. -4Ca = We = 10 , to make our simulation results comparable with the experimental results of Mason and Clark that are for quasi-static conditions. Mason and Clark [114] do not provide an exact contact angle value for their experiments and only mention that it is very close to zero. For our numerical simulations shown in Figure 4-14, the contact angle is assumed to be constant and equal to zero degrees. The figure indicates very good agreement (with the maximum relative difference of 5%) between our simulation results and the experimental measurements of Mason and Clark. Our numerical simulations are also compared with another set of experimental measurements presented later in this chapter. - 74 - Figure 4-12. Comparison between our numerical simulations and the experimental data of Mason and Clark [114]. 4.5.4.3. Comparison with another numerical tool Darhuber et al. [118] simulated the liquid transfer fraction between two plates at quasi-static conditions. They found that the quasi-static condition corresponds to low Reynolds and very low capillary numbers. The numerical simulations of Darhuber et al. were performed with the Surface Evolver software [119]. The software evolves the surface towards the minimal energy by a gradient descent method, and is used for studying surfaces shaped by energies such as surface tension, gravitational energy, etc. For our simulations, the full Navier-Stokes equations of the two-fluid interactions are solved with the VOF method; this means that the numerical tool used in this work differs from that of Darhuber et al. For the numerical simulations, particles are assumed to be equal-sized, and the contact angle on one particle is fixed at 30iθ = o while that of the other particle varies in a range of 5 60jθ = − o . Also, since the liquid bridge volume was not exactly specified by Darhuber et al., two different liquid bridge volumes are used for our simulations, i.e., V = 0.03 and 0.1. To comply with the quasi-static assumption of Darhuber et - 75 - al., the capillary and Weber numbers in our simulations are in the order of 310− and 410− , respectively. Also, the gravitational effects are neglected in our simulations. Figure 4-13 compares the liquid transfer fraction predictions of our numerical simulations with those of the numerical simulations of Darhuber et al. [118]. This figure indicates good agreement for the simulations with a liquid bridge volume of 0.03, and fair agreement for the simulations with a liquid bridge volume of 0.1. The better agreement found for the lower bridge volume corresponds with the assumption of Darhuber et al. of a very low amount of liquid in the bridge. Figure 4-13. Comparison between the liquid transfer fraction predictions of our numerical simulations and the numerical simulations of Darhuber et al. [118]. - 76 - 4.6. COMPARISON BETWEEN SIMPLIFIED MODEL AND NUMERICAL SIMULATION 4.6.1. Comparison between rupture distances predicted with simplified model and numerical simulation Figure 4-14 graphically compares the rupture distances predicted using both the numerical simulations and the simplified model with the experimental data of Mazzone et al. [85]. According to this figure, the rupture distance predictions of the numerical simulation are very close to the experimental data, while the simplified model tends to over-predict. The better agreement obtained with the numerical simulations can be attributed to the inclusion of the inertial, gravitational, and viscous effects in the simulations, while the simplified model does not explicitly account for them. For the reported numerical simulations, the Bond, capillary, and Weber numbers are equal to 1.254Bo = , 54.4 10Ca −= × , and 73.2 10We −= × , respectively. Due to the lack of experimental data, the results of the numerical simulations are only partially validated. A systematic experiential study is necessary for further validation of the numerical results in future. Figure 4-144. Comparison of the rupture distance predictions using both the numerical simulations and the simplified model with the experimental data of Mazzone et al. (1986). 4 There is no information on the errors of the experimental data. - 77 - 4.6.2. Variations of rupture distance with Weber, capillary and Bond numbers According to the non-dimensional analysis presented earlier, the dimensionless rupture distance depends on the dimensionless liquid bridge volume, contact angle, Weber number ( We ), capillary number ( Ca ), and Bond number ( Bo ). The relationship between the rupture distance and liquid bridge volume and contact angle was analyzed using the simplified model earlier in this chapter. Similar trends were observed with the numerical simulations. The two main objectives of this sub-section are as follows: (1) to investigate the relationship between the rupture distance and the Weber, Bond, and capillary numbers, and (2) to quantify the range of these dimensionless numbers in which the simplified model applies. The latter is achieved by performing three series of numerical simulations with different capillary, Weber, and Bond numbers. To keep the analyses simple, the gravitational effects are neglected in the first two series of simulations and are considered only in the last series (when investigating the effect of the Bond number). In addition, in all the simulations, three cases with different liquid bridge volumes and contact angles are studied: (1) 0.1V = and 30θ = o , (2) 0.2V = and 20θ = o , and (3) 0.3V = and 10θ = o . For the comparison purposes, the rupture distances for these cases are obtained with the simplified model and are 0.65, 0.76, and 0.78, respectively. Figure 4-15 shows the variations of the rupture distance with respect to the Weber number. In all these simulations, the capillary number is fixed at 310Ca −= . According to this figure, the rupture distance increases as the Weber number increases. Also, as the Weber number decreases, the rupture distance approaches an asymptotic value. The asymptotic values are 0.64, 0.83, and 0.90 for cases 1-3, respectively, and correspond to the rupture distance for Weber numbers equal to or less than 410− . Comparing these values with the predictions of the simplified model, it is clear that the agreement between the numerical simulations and the simplified model becomes worse when the liquid bridge volume is increased. Nevertheless, the asymptotic values are fairly close to the predictions of the simplified model. The interpretation of this observation is that at low Weber numbers, the surface tension effects dominate the inertial effects, hence the results of the numerical simulations agree with those of the simplified model. Further analysis shows that for Weber numbers up to 210We −= , the predicted rupture distances are within the 5% range of the - 78 - asymptotic values. In conclusion, for Weber numbers less than 10-2, the surface tension effects dominate the inertial effects, thus the simplified model applies. Figure 4-15. Variations of the simulated rupture distance with respect to the Weber number. Figure 4-16 shows the variations of the rupture distance with respect to the capillary number. In all these simulations, the Weber number is fixed at We = 10-4. According to this figure, the rupture distance increases as the capillary number increases. Also, as the capillary number decreases, the rupture distance approaches an asymptotic value. The asymptotic values are 0.62, 0.80, and 0.87 for cases 1-3, respectively, and correspond to the rupture distances for capillary numbers equal to or less than 10-4. Comparing these values with the predictions of the simplified model, it is again observed that the asymptotic values are fairly close to the predictions of the simplified model. The interpretation is that at low capillary numbers, the surface tension effects dominate the viscous effects and the results of the numerical simulation agree with those of the simplified model. Further analysis shows that for capillary numbers up to Ca = 10-3, the predicted rupture distances are within the 5% range of the asymptotic values. In conclusion, for capillary numbers less than 10-3, the surface tension effects dominate the viscous effects; hence the simplified model applies. - 79 - Figure 4-16. Variations of the simulated rupture distance with respect to the capillary number. Figure 4-17 shows the variations of the rupture distance with respect to the Bond number. In all these simulations, the capillary and Weber numbers are fixed at 310Ca −= and 410We −= . According to this figure, the rupture distance decreases as the Bond number increases. Also, as the Bond number decreases, the rupture distance approaches an asymptotic value. These asymptotic values are 0.64, 0.83, and 0.9 for cases 1-3, respectively, and correspond to the rupture distances for Bond numbers equal to or less than 10-1. It is stressed that comparing these values with the predictions of the simplified model, the agreement between the numerical simulations and the simplified model becomes worse as the liquid bridge volume increases. Nevertheless, the asymptotic values are fairly close to the predictions of the simplified model. The interpretation is that at low Bond numbers, the surface tension effects dominate the gravitational effects so that the results of the numerical simulation agree with those of the simplified model. Further analysis shows that for Bond numbers up to 1Bo = , the predicted rupture distances are within the 5% range of the asymptotic values. In conclusion, for Bond numbers less than 1, the surface tension effects dominate the gravitational effects and the simplified model applies. - 80 - Figure 4-17. Variations of the simulated rupture distance with respect to the Bond number. 4.6.3. Comparison between liquid bridge profiles of simplified model and numerical simulation One of the key assumptions of the simplified model is that the liquid bridge shape in the vicinity of the rupture instant can be approximated with a second-order polynomial. In addition, it is assumed that at any instant, the liquid bridge is stable, i.e., small perturbations do not grow. In this section, these two assumptions are analyzed in further detail. In Figure 4-18, the bridge profiles predicted using the simplified model (presented in green squares) and the numerical simulation at the rupture instant and a few milliseconds before rupture are shown. According to this figure, at the rupture instant, the simplified model does not represent the simulated bridge shape appropriately; however, a few milliseconds before rupture (here, 9 ms), the liquid bridge profile predicted using the simplified model resembles the simulated bridge profile. In order to further analyze this observation, the contour plots of the velocity magnitude of the numerical simulations are compared, as shown in Figure 4-19. According to this figure, at the rupture instant, the bridge is evolving with a very high velocity, as the maximum velocity is orders of magnitude greater than the prescribed separation velocity, i.e. 0.0005 m/s for each particle. However, a few milliseconds before rupture, the maximum - 81 - velocity magnitude (observed at the interface of liquid and solid) is comparable to the prescribed separation velocity. Comparing these two velocity contour plots, it can be concluded that at the rupture instant, the liquid bridge profile is unstable, whereas a few milliseconds before rupture, a stable liquid profile is observed. The liquid bridge profiles obtained with the numerical simulations were also fitted with second-, fourth-, and sixth-order polynomials. The confidence levels of the curve fittings for the bridge profile at the rupture instant were 91, 99, 99.9% for the second-, fourth-, and sixth-order polynomials, respectively. However, for the liquid bridge a few milliseconds before rupture, these confidence levels, in the same order, were 98, 99.95 and ~100%. It can be concluded that while a second-order polynomial accurately predicts the bridge shape a few milliseconds before rupture, it is not an appropriate choice for representing the bridge shape at the rupture distance. In other words, the second-order polynomial approximation is not an appropriate choice for attempting to match the more involved and complex computational models such as numerical solution of the governing Navier-Stokes equations. Figure 4-18. Comparison between the bridge profiles predicted using the simplified model and numerical simulations at: (a) rupture instant, and (b) a few milliseconds before rupture, for θ = 30o. - 82 - Figure 4-19. Contour plots of the velocity magnitude (in m/s) of the numerical simulations at: (a) rupture instant, and (b) a few milliseconds before rupture, for θ = 30o. In the above figures, the liquid-solid contact angle on both particles was constant and equal to θ = 30o. The impact of the contact angle on the liquid bridge shape was also studied for contact angles of 10, 20 and 40 degrees on both particles, and similar trends were found. In conclusion, the second-order polynomial, as used in the simplified model, is an appropriate representation of the gas-liquid interface (between particles with identical contact angles) only at stable conditions. The applicability of the parabolic profile was also investigated for particles with different contact angles. Figure 4-20 shows the comparison between the bridge profiles predicted using the simplified model (indicated in green squares) and the numerical simulation, for particles with θi = 30o and θj = 40o. According to this figure, even a few milliseconds before rupture, the bridge profile no longer resembles a second-order polynomial. It is concluded that a second-order polynomial is not an appropriate representation of the simulated liquid bridge profile for particles with different contact angles; therefore, the simplified model is not an appropriate model for predicting the bridge evolution between particles with different contact angles. - 83 - Figure 4-20. Comparison of the bridge profiles predicted using the simplified model and numerical simulation at: (a) rupture instant, and (b) a few milliseconds before rupture, with θi = 30o and θj = 40o. 4.7. LIMITATIONS OF SIMPLIFIED MODEL The liquid bridge shapes predicted using the numerical simulations and the simplified model were compared. It was concluded that a second-order polynomial is not an appropriate representation of the simulated bridge shape for two conditions: (1) unstable liquid profiles, and (2) the liquid bridges between particles with different contact angles. Therefore, for such conditions, another approach such as numerical simulation should be used. In addition, with the aid of the numerical simulations, the importance of the inertial, viscous, and gravitational effects in comparison to the surface tension effects was studied. It was concluded that the surface tension effects dominate the viscous, inertial, and gravitational effects for 310Ca −< , 210We −< , and 1Bo < , respectively. In other words, the simplified model is only applicable to the liquid bridges problems with dominant surface tension effect. If any of the above dimensionless numbers exceeds the limiting value (i.e. surface tension effects are not dominant), an alternative approach such as the numerical solution of the governing equations should be employed. - 84 - 4.8. INFLUENCE OF CONTACT ANGLE5 As shown in the previous section, the simplified model has limitations; for instance, it is not suitable for studying the rupture of stretching liquid bridges between particles with different contact angles and capillary numbers greater than 10-3. In this section, the numerical solutions are used for investigating the influence of the contact angle on the rupture distance and liquid distribution of stretching liquid bridges. The model and solution methodology have been explained earlier. Hence, they are not repeated in this section. 4.8.1. Summary of cases Six different simulation cases are considered in this section and are summarized in Table 4-4. In these simulations, the dimensionless liquid bridge volumes vary between 0.01 - 0.1. As before, the contact angles are static and constant. For the simulations associated with these cases, the contact angle of the left-hand side (LHS) particle is always equal to 30 degrees, while the contact angle of the right-hand side (RHS) particle is either fixed at one value for all the simulations (e.g. cases 3-6) or varies from one simulation to another (e.g. cases 1-2). Table 4-4. Summary of the cases considered in this section. Case # Dimensionless liquid bridge volume Contact angle of LHS particle, degrees Contact angle of RHS particle, degrees 1 0.01 30 5-60 2 0.03 30 5-60 3 0.05 30 10 4 0.05 30 20 5 0.1 30 40 6 0.1 30 50 5 A version of this section has been published. Darabi et al. (2010) Numerical Studies of Stretching Liquid Bridges between Two Solid Spherical Particles with Different Contact Angles. ASME International Mechanical Engineering Congress and Exposition (IMECE). IMECE2010-37564. - 85 - 4.8.2. Grid independence study Similar to before, the grid independence study is conducted on three different grids: coarse, intermediate, and fine, with 5040, 10800, and 24300 initial number of control volumes, respectively. This time both the rupture distance and liquid transfer fraction are addressed. Table 4-5 and Table 4-6 contain a summary of the simulation results for cases 5 and 6 with the three grids. In each table, the initial number of control volumes, the predicted rupture distance, and the predicted liquid transfer fraction for the right-hand side particle (that has a larger contact angle) is provided. Two series of simulations with different capillary numbers, i.e. Ca = 0.001 and 0.01 are considered for each case. According to these tables, as the number of control volumes increases (i.e. grid becomes finer), the difference between both the predicted rupture distances and the predicated liquid transfer fractions decreases; this indicates grid convergence. The relative difference between the predicted values for the fine and intermediate grids is less than 1%; therefore, in this entire section, the intermediate grid is used. Table 4-5. Summary of the grid independence study for case 5. Grid Rupture distance (Ca = 0.001) Transfer fraction (Ca = 0.001) Rupture distance (Ca = 0.01) Transfer fraction (Ca = 0.01) Coarse 0.593 0.162 0.708 0.392 Intermediate 0.596 0.160 0.713 0.396 Fine 0.598 0.159 0.718 0.398 Table 4-6. Summary of the grid independence study for case 6. Grid Rupture distance (Ca = 0.001) Transfer fraction (Ca = 0.001) Rupture distance (Ca = 0.01) Transfer fraction (Ca = 0.01) Coarse 0.557 0.050 0.698 0.196 Intermediate 0.560 0.049 0.703 0.202 Fine 0.562 0.049 0.707 0.208 - 86 - 4.8.3. Results and discussion 4.8.3.1. Effect of contact angle and capillary number on rupture distance Figure 4-21 shows the variations of the rupture distance with respect to the contact angle for cases 1-2. According to this figure, the maximum rupture distance is obtained when both particles have the same contact angles. In addition, this figure shows that the rupture distance decreases as the difference between the contact angles increases. The latter observation could be only valid for the quasi-static conditions, as cases 1-2 are performed for very low capillary and Weber numbers. More simulations with varying capillary numbers covering a wider range of conditions, rather than only for the quasi-static cases, are considered for cases 3-6 and are discussed as follows. Figure 4-21. Variations of the rupture distance with respect to contact angle for cases 1 and 2. Figure 4-22 and Figure 4-23 show the variations of the rupture distance with respect to the capillary number for cases 3-6. According to these figures: 1) similar to the previous findings (discussed earlier in this chapter), the rupture distance increases as the capillary number increases, 2) at low capillary numbers (i.e. quasi-static conditions), the case with the greater contact angle difference has the lower rupture distance, as also shown in Figure 4-21, and 3) the - 87 - impact of the contact angle difference on the rupture distance gradually vanishes as the capillary number increases. In other words, at high capillary numbers, the rupture distance does not depend on the contact angles. This is an interesting finding that could be explored further in future. Figure 4-22. Variations of the rupture distance with respect to the capillary number for cases 3 and 4. - 88 - Figure 4-23. Variations of the rupture distance with respect to the capillary number for cases 5 and 6. 4.8.3.2. Effect of contact angle and capillary number on liquid distribution Figure 4-24 shows the variations of the liquid transfer fraction with respect to the contact angle for cases 1-2. Note that these simulations are performed for the quasi-static conditions, i.e. for very low capillary and Weber numbers. According to this figure, under quasi-static conditions the liquid transfer fraction is strongly dependant on the contact angles of the two solid particles. The liquid is uniformly distributed between the two solid particles whenever they have identical contact angles. However, for particles with different contact angles, the liquid is non-uniformly distributed and more liquid is transferred to the particle with the smaller contact angle; the non- uniformity increases as the difference between the contact angles increases. - 89 - Figure 4-24. Variations of the liquid transfer fraction with respect to contact angle for cases 1 and 2. Figure 4-25 and Figure 4-26 show the variations of the liquid transfer fraction with respect to the capillary number for cases 3-6. These figures show that at low capillary numbers, the liquid distribution is strongly affected by the difference between the contact angles of the two particles. This is consistent with the previous finding for the quasi-static conditions (discussed earlier in this chapter). However, as the capillary number increases, the liquid distribution becomes less dependent on the difference between the contact angles, until it becomes almost independent at large capillary numbers. Nevertheless, according to the above figures, more liquid is always transferred to the particle with the smaller contact angle. Yakhnin and Chadov [120, 121] experimentally studied the liquid transfer between two flat plates under a variety of conditions. For their experiments, they considered surfaces made from different materials, different combinations of water/glycerol and ethanol/glycerol to obtain bridging liquids with different surface tension and viscous properties; in addition, they considered different separation velocities for their experiment. In order to analyze the results, they plotted the variations of the liquid transfer fraction with respect to the separation velocity - 90 - for their experimental measurements. Consistent with their trends, three different liquid transfer zones are found for our simulations, as shown in Figure 4-25 and Figure 4-26: (1) at high capillary numbers, the liquid transfer fraction is approaching an asymptotic value of 0.5. This indicates the insensitivity of the liquid distribution to the difference between the contact angles of the two solid surfaces (either flat plates or spherical particles) at high capillary numbers; we call this the “dynamic” zone; (2) at low capillary numbers, the liquid transfer fraction is also approaching an asymptotic value, for which the liquid transfer fraction depends on the difference between the contact angles of the two particles (a measure of how the contacting liquid-solid media interact). Similar to Yakhnin and Chadov, we refer to this zone as the “quasi-static” or “quasi-equilibrium” zone. For this zone, the higher the difference between the contact angles, more liquid is transferred to the surface with the smaller contact angle. A similar finding was earlier discussed for Figure 4-24; and (3) at intermediate capillary numbers, the curve connecting the above limiting conditions resembles an S-shape. We call this the “transition” zone, for which the liquid distribution depends on both the system dynamics (here, capillary number) and liquid- solid contact angles. In short, our numerical findings are in qualitative agreement with the experimental findings of Yakhnin and Chadov [120, 121]. In order to further quantify the matter, we analyzed our plots and concluded that: the quasi-equilibrium zone covers low capillary numbers up to 10-3 (which is consistent with our earlier numerical findings for particles with identical contact angles). The dynamic zone covers capillary numbers around 10-1 and higher and the transition zone lies between the above limits, i.e. capillary numbers from ~ 0.001 up to 0.1. Note that in these simulations the gravitational effects are neglected; the effect of gravity on the rupture of liquid bridges is discussed in the following section. - 91 - Figure 4-25. Variations of the liquid transfer fraction with respect to the capillary number for cases 3 & 4. Figure 4-26. Variations of the liquid transfer fraction with respect to the capillary number for cases 5 & 6. - 92 - 4.9. INFLUENCE OF GRAVITY6 In this section, numerical simulations of the governing Navier-Stokes equations are used to investigate the effect of gravity on the rupture of stretching liquid bridges between two equal- sized solid spherical particles. The gravity acts along the line connecting the centers of the particles. In addition, the model, its solution methodology, and the grid independence study have been discussed earlier in the chapter; therefore, they are not repeated. 4.9.1. Summary of cases For the simulations presented in this section, two different series of calculations are performed: in Series 1, the contact angle of both particles is identical; in the other series, particles have different contact angles. Both of the series are summarized in Table 4-7. Table 4-7. Summary of the cases considered in this section. Series # Dimensionless liquid bridge volume Contact angle of LHS particle [degrees] Contact angle of RHS particle [degrees] 1 0.1 – 0.3 10 – 30 = θLHS 2 0.1 30 40 or 50 4.9.2. Results and discussion 4.9.2.1. Effect of gravity on rupture distance The variations of the rupture distance with respect to the Bond number were discussed earlier (see Figure 4-17). According to this figure, the rupture distance decreases as the Bond number increases and the rupture distance approaches an asymptotic value for small Bond numbers. From this figure it was concluded that the gravitational effects on the rupture distance are negligible for Bond numbers smaller than 10-1. Figure 4-27 shows the variations of the rupture distance with respect to the capillary number for particles with different contact angles in the presence of gravity. The effect of the capillary 6 A version of this section has been submitted for publication. Darabi et al. Numerical Investigations of Liquid Bridge Rupture in the Presence of Gravity. Submitted on November 10, 2010. - 93 - number and contact angle on the rupture distance has been discussed earlier; therefore, only the effect of gravity is discussed here. For the simulations shown in this figure, the Bond number is equal to one. According to the figure, at small capillary numbers, the gravitational effects impact the rupture distance. The rupture distance decreases for cases in which the gravitational acceleration is directed to the particle with the smaller contact angle. Conversely, the rupture distance increases if the gravitational acceleration is directed to the particle with the greater contact angle. This can be explained as follows: for the first case, having a smaller contact and aiding gravitational effects both lead the liquid accumulation on the particle with the smaller contact angle; hence the lower rupture distance. On the other hand, as the aiding gravitational effects aim toward the particle with the greater contact angle, the rupture distance decreases. The above observation means that the rupture distance increases as the liquid distribution becomes more uniform. Further discussion regarding the liquid distribution is provided later. Figure 4-27 also shows that the effect of gravity on the rupture distance becomes less significant as the capillary number increases. (a) - 94 - (b) Figure 4-27. Variations of the rupture distance with respect to the capillary number: (a) θl = 30 & θr = 40o, (b) θl = 30 & θr = 50o. 4.9.2.2. Effect of gravity on liquid distribution Figure 4-28 shows the variations of the liquid transfer fraction with respect to the Bond number. The figure provides the transfer fraction of the left-hand side particle while the gravity is directed toward the other particle. Again, note the results obtained in this figure are for simulations at a very low capillary number between particles with identical contact angles. According to this figure, under such conditions the liquid transfer fraction is greatly affected by the value of the Bond number; uniform distribution is obtained at small Bond numbers and the liquid transfer fraction reduces as the Bond number increases. - 95 - Figure 4-28. Variations of the liquid transfer fraction with respect to the Bond number. Figure 4-29 shows the variations of the transfer fraction (of the particle with the greater contact angle) with respect to the capillary number for particles with different contact angles. Since the effect of the capillary number and contact angle on the liquid distribution has been discussed earlier, the effect of gravity only is discussed here. According to this figure, the gravitational effects shift the curve depending on the direction of gravity. If the gravity is directed to the particle with the greater contact angle, the curve is shifted to the top, meaning more liquid is transferred to the particle with the aiding gravitational effects. On the other hand, the curve is shifted to the bottom, as gravity is aimed toward the particle with the smaller contact angle, meaning further accumulation of liquid on that particle. The figure also indicates that the liquid transfer fraction curve with no gravitational effects lies approximately between the two above curves. - 96 - (a) (b) Figure 4-29. Variations of the transfer fraction with respect to the capillary number: (a) θl = 30 & θr = 40o, (b) θl = 30 & θr = 50o. - 97 - 4.10. SUMMARY AND CONCLUSIONS A simplified mathematical model and the numerical simulations of the equations were used to study the rupture of stretching liquid bridges between two smooth equal-sized solid spherical particles. For the numerical simulations, a commercial Computational Fluid Dynamics (CFD) tool, FLUENT, was used. The code solves the full Navier-Stokes equations, hence, it explicitly accounts for the viscous, inertial, gravitational, and surface tension effects. The simplified model was used to investigate the impact of liquid bridge volume and contact angle on the rupture distance. It was shown that the predictions of the simplified model were in good agreement with the related expressions available in the literature. The numerical simulations were used to simulate the time-dependant evolution of the liquid bridge. It was shown that as particles move farther apart, the liquid bridge stretches and its neck becomes narrower. This continuous elongation persists until the bridge can no longer withstand any additional elongation and ruptures. Several numerical simulations were used to investigate the effect of the capillary, Weber, and Bond numbers on the rupture distance. The results showed that for simulations with dominant capillary effects, i.e. 310Ca −< , 210We −< , and 1Bo < , the rupture distance is within 5% of an asymptotic value. For all these simulations, the asymptotic value was close to the rupture distances predicted using the simplified model. In conclusion, the simplified model is applicable to problems in which the surface tension effects dominate the viscous, inertial, and gravitational effects. Also, the bridge profiles predicted using the numerical simulations and simplified model were compared. It was shown that for liquid bridges between particles with identical contact angles and under stable conditions, the bridge profile can be represented with a second-order polynomial with enough accuracy. Therefore, the simplified model is appropriate for studying such rupture problems. However, for liquid bridges between particles with different contact angles, the parabolic assumption is no longer a good representation of the bridge profile; hence, the simplified model is not suitable. To summarize, the simplified model is a useful tool for studying the evolution of stretching liquid bridges, yet it has limitations. It is only applicable to bridges with dominant surface tension effects, under quasi-static conditions, and between particles with identical contact angles. However, the numerical simulations of the Navier-Stokes equations are more general. As a result, they can be - 98 - used for studying a wide range of problems, e.g. rupture problems with non-negligible viscous, inertial, and gravitational effects and between particles with different contact angles. Numerical simulations were performed to study the effect of the contact angle and capillary number on the rupture distance and liquid transfer fraction of liquid bridges between particles with different contact angles. The results showed that the rupture distance increases as the capillary number is increased and that the maximum rupture distance under quasi-static conditions (i.e. small capillary numbers) is obtained when both particles have identical contact angles. Also, it was observed that for the quasi-static conditions, the liquid distribution is highly dependant on the contact angles of the two particles. The liquid is evenly distributed when particles have identical contact angles. Otherwise, the liquid is non-uniformly distributed and more liquid is transferred to the particle with the smaller contact angle; the amount of the transferred liquid increases as the difference between the contact angles increases. It was shown that for ruptures with intermediate capillary numbers, the liquid distribution depends on both the difference between the contact angles and the capillary number. Finally, it was found that at high capillary numbers, the liquid distribution is almost uniform and is insensitive to the difference between the contact angles. Also, numerical simulations were used to study the effect of gravity on the rupture distance and liquid transfer fraction of stretching liquid bridges. The results showed that for particles with identical contact angles, the inclusion of gravity affects both the rupture distance and liquid transfer fraction. The rupture distance decreases and more liquid is accumulated on the particle with aiding gravitational effects as the Bond number increases. For particles with different contact angles, the inclusion of gravity either decreases or increases the rupture distance. This change depends on how liquid is distributed considering both the contact angles and gravitational effects. Both of these effects become less significant as the capillary number is increased. Also, due to the inclusion of gravity, the liquid transfer curves for particles with different contact angles are shifted (either up or down) toward further accumulation of liquid on the particle with aiding gravitational effects. - 99 - Chapter 5. DEM INVESTIGATIONS OF FLUIDIZED BEDS IN THE PRESENCE OF LIQUID COATING7 5.1. INTRODUCTION With the addition of liquid and its binding effects, the operation of fluidized beds becomes more complex. This is primarily because of the liquid effects that induce cohesion between particles, making the experimentation on the fluidization system more complicated than before. A good alternative for obtaining detailed information on the transport phenomena inside fluidization systems is performing DEM simulations. This is because they account for the motion of individual particles, particle-particle, and particle-wall interactions. DEM simulations have many advantages, but are very expensive computationally; hence, they are often used to investigate systems only on a small scale. Recently, DEM simulations have been employed to study the effect of cohesion between particles on a number of systems, e.g. fluidized beds [82], tumbling granulators [109], and vibrated beds [110]. However, as explained in the introductory chapter, the available DEM studies in the literature have limitations, leaving room for substantial improvement in this area. The main objective of this chapter is to perform DEM simulations in order to study the impact of liquid coating on the fluidized bed behaviour. It is assumed that the liquid coating affects the system behaviour by influencing the coefficient of restitution, used for determining the outcome of particle-particle interactions. The proposed wet coefficient of restitution, shown in Chapter 2 of this thesis, is implemented in an open-source numerical tool - MFIX. The modified numerical tool is used to perform DEM simulations of a fluidized bed consisting of mono-sized solid spherical particles pre-coated with identical liquid coatings. Numerous simulations are performed to investigate the effect of coating viscosity and coating thickness on the fluidization behaviour under a variety of conditions. In order to analyze the results, snapshots of the instantaneous particle positions are observed and time-averaged values of the average particle height, wet coefficient of restitution, and relative normal collision velocity are analyzed. 7 A version of this section has been submitted for publication. Darabi et al. DEM Investigation of Fluidized Beds in the Presence of Liquid Coating. Submitted on December 13, 2010. - 100 - 5.2. MODEL DESCRIPTION For the simulations performed in this chapter, an open-source code known as MFIX-DEM is used. MFIX stands for Multiphase Flow with Interphase eXchange and is a computer code developed in the National Energy Technology Laboratory (NETL); the code is used to study the complex flow in fluid-solid systems. In the MFIX-DEM code, the gas phase is resolved by using the governing conservation equations of mass and momentum. For the solid phase, the Discrete Element Method (DEM) is used where each individual particle of the dispersed phase is simulated. Due to the high computational cost associated with the particle neighbour search algorithm and very small time steps required for resolving the interaction of particles with the soft-sphere approach, MFIX-DEM simulations are limited to small-scale problems. Further description regarding the soft-sphere approach is provided later. In the following sub-sections, the governing equations for the gas and solid phases are described. 5.2.1. Gas-phase The governing conservation equations of mass and momentum for the gas phase, as implemented in MFIX [122], neglecting turbulence, phase change, chemical reactions, growth, aggregation, and breakage phenomena, are: ( ) ( )g. v 0g g g gt ε ρ ε ρ ∂ + ∇ = ∂ ( 5-1) ( )g 1 v . g IMg g g g g gmm D S Dt ε ρ ε ρ = = ∇ + −∑ ( 5-2) where gε is the volume fraction of the gas phase, gρ is the density of the gas phase, gv is the volume-averaged velocity of the gas-phase, gS is the stress-tensor of the gas phase, g is the gravitational acceleration, and gmI is the momentum transfer term between the gas phase and the m th solid phase. The stress-tensor of the gas phase can be written as: - 101 - g g gS P I τ= − + ( 5-3) where gP is the pressure of the gas phase and gτ is the shear stress tensor of the gas phase, given by: 2 .tr( )g g g g gD D Iτ µ λ= + ( 5-4) where ( )g g1/ 2 v v TgD = ∇ + ∇ is the strain rate tensor of the gas phase, and gµ and gλ are the dynamic and second coefficients of viscosity of the gas phase, respectively. Further description of the above governing equation is provided in MFIX documentation [122]. 5.2.2. Solid-phase As mentioned earlier, in MFIX-DEM individual particles are used to represent the dispersed solid phase. In the DEM approach, the solid phase can consist of M different groups of solid particles where the mth solid phase has Nm solid spherical particles with equal diameters Dm and densities smρ . Therefore, each solid group is differentiable from others, knowing the radius and density of its particles. Knowing that there are M group of particles and each group has Nm particles, the total number of particles will be 1 M p mm N N = =∑ . These Np particles at time t are represented by }{ (i) (i) (i) (i) (i)X ( ),V ( ), ( ), ( ), ( ), 1,..., pt t t D t t i Nω ρ = where (i)X ( )t designates the position of ith particle, (i)V ( )t and (i)( )tω denote its linear and angular velocity, respectively, and (i) ( )D t and (i)( )tρ represent its diameter and density, respectively. In the DEM approach, the trajectory of each particle is resolved using the Newtonian equations of motion. Such equations for determining the position, linear, and angular velocities of the ith particle can be written as follows: (i) (i)dX ( ) V ( )t t dt = ( 5-5) ( ) ( ) ( ) ( ) (i k,m) ( ) T d c V ( ) F g + F ( ) F ( ) i i i i itm m t t dt ∈ = = + ( 5-6) - 102 - (i) ( ) (i)( ) Ti tI dt ω = ( 5-7) where ( )im is the mass of the ith particle, ( )TF ( )i t is the total force (the summation of the drag, gravitational, and contact forces) acting on the ith particle, (i k,m)dF ∈ is the total drag force (pressure and viscous) on the ith particle that resides in the kth computational cell and belongs to the mth solid group, ( )cF ( )i t is the net contact force as a result of contact with other neighboring particles, ( )iI is the moment of inertia of the ith particle, and (i)T is the total torque acting on the ith particle. A detailed description regarding the calculation of the above forces and torques is provided in MFIX-DEM documentation [123]. In order to solve the above equations of motion, the mass and moment of inertia of the particle is required. As highlighted earlier, if a particle belongs to the mth solid phase, its diameter and density will be known and they are equal to Dm and smρ , respectively. Knowing the diameter and density of a spherical particle, its mass and moment of inertia can be calculated with the following equations: 3(i) ( ) (i) 6 i D m piρ= ( 5-8) 2( ) (i) ( ) 10 i i m DI = ( 5-9) One of the key features of the DEM approach that makes it advantageous over the continuous treatment of the solid phase is the explicit modeling of the contact force, accounting for the particle-particle interactions. Each contact force on a target particle is a net force and is the summation of all the contact forces that are determined based on the binary collision of the target particle with its neighbours. In MFIX-DEM, the soft-sphere approach is used for calculating the contact forces. The soft-sphere approach allows for multiple particle interactions where the net force is obtained by adding the results of all the binary interactions. In this work, the contact forces are calculated from the deformation history using a spring-dashpot mode model. This model is based on the formulation proposed by Cundall and Strack [6], often regarded as the first - 103 - soft-sphere publication in the literature. A complete description of the soft-sphere model of DEM-MFIX is provided in MFIX-DEM documentation [123]. The DEM solutions have been coupled with the numerical solutions of the governing Navier-Stokes equation in order to simulate the complex flow in different systems; the work of Tsuji et al. [8] is usually regarded as the first publication that uses the above approach to perform 2D simulations of fluidized beds. Some similar works have been referenced earlier in Chapter 1. 5.2.3. Model for interaction of wet particles In order to account for the effect of the liquid coating on solid particles, the particle interactions are modified by implementing a wet coefficient of restitution. This is achieved by including the coefficient of restitution, as proposed earlier in Chapter 2 of this thesis, in MFIX. Similar treatment has been proposed before in the literature [53] by implementing a wet coefficient of restitution expression that depends only on the moisture content. However, the proposed wet coefficient of restitution in this thesis is more general as it is based on a number of key properties such as the liquid viscosity, liquid surface tension, relative collision velocity, etc., rather than only the moisture content. This allows performing more comprehensive studies of fluidized beds in the presence of liquid. The model calculates the wet coefficient of restitution for the binary collision of two wet particles approaching with a non-zero relative normal collision velocity and is summarized in Table 5-1. The model determines the wet coefficient of restitution based on the total kinetic energy of the system of two particles (critical energy value) and approximate analytical values of maximum possible energy dissipation. The maximum possible energy dissipation accounts for the energy loss due to the viscous and capillary effects and inelastic collision. Further detail of this model is available in Chapter 2 of the thesis. In the MFIX-DEM, the wall-particle interactions are dealt with similar to the particle-particle interactions and can be modified with the same idea as above. However, in the current work, it is assumed that the liquid only affects the normal collision of particles and its impact on the coefficient of friction and wall-particle interactions is neglected. - 104 - Table 5-1. Summary of the constitutive equations of the wet coefficient of restitution model. Proposed Coalescence Model: ,max 2 .max 02 0 2 ,max 0 1 , 0, t wet t wet t L e L mu mu e L mu = − < = ≥ Maximum Possible Energy Dissipation: ,max ,1 ,1 ,2 ,22( )t vis cap c vis capL W W L W W= + + + + Viscous Work during Approach Stage: 20 0 ,1 ( ) ( )ln ln( ) 2 ( )vis a v a f h f hAW A f h St f h = − Viscous Work during Rebound Stage: 2 1 1 ,2 0 ( ) ( )ln ln 2 ( ) ( ) c vis v a a Aeuf h f hAW St f h u f h − = − Capillary Work during Approach Stage: ( ) 02 2,1 2 cos a h cap h W R D D api σ θ= − − + Capillary Work during Rebound Stage: ( ) 0.52 2,2 2 cos rup a h cap h W R D D api σ θ= − + Contact Loss: 2 21 (1 ) 2c c L m e u= − Contact velocity (right after contact): 0 0 ( )11 ln ( )c v a f h u u St f h = − − Constants and Functions: 08 9 pu RSt ρ µ = ( ) 2 2 2 2 2 ( ) D D af D D D a + = + + 2 bVa Rpi = 2 2( ) 4 2b w V R H r Dpi = − 2 ( ) 2 rH r D R = + 2 0 3 2 A R upiµ= In order to implement the model for the wet coefficient of restitution, some information is extracted from the MFIX-DEM code; in our implementation the relative normal collision velocity of the two wet particles and their sizes is taken from MFIX-DEM. The rest of the information required for calculating the wet coefficient of restitution is provided as input for the model; this includes the size of dry particles, viscosity and surface tension of the liquid coating, etc. Given the calculated wet coefficient of restitution for the normal collision of particles, the normal dashpot coefficient for the DEM simulations is calculated; this is performed using the following equation: - 105 - , , , 2 2 , 2 ln ln eff n ml n ml n ml n ml m k e e η pi = + ( 5-10) where ,n mlη is the normal dashpot damping coefficient of the mth and lth particles, /( )eff m l m lm m m m m= + is their effective mass, .n mlk is the spring stiffness coefficient, and ,n mle is the normal coefficient of restitution (and here is the calculated wet coefficient of restitution). The above is repeated to calculate the wet coefficient of restitution and the normal dashpot damping coefficient for all the binary collisions. Given the calculated normal dashpot damping coefficients, the contact forces due to the collision of particles are calculated and the rest of the calculations required for the DEM simulations are performed. Our proposed wet coefficient of restitution (accounting for the binary collision of wet particles) is compatible with the way the solid phase calculations are performed in MFIX-DEM; this is because in the implemented soft- sphere approach in MFIX-DEM, first the binary collisions are resolved and then the net contact force for each particle is calculated. Although the soft-sphere approach accounts for the particle deformations, the current implementation does not fully take advantage of this capability. This is because the proposed wet coefficient of restitution already incorporates the inelasticity of particles by assuming a (less than unity) dry coefficient of restitution. 5.2.4. Simulation setup For the simulations performed in this chapter, a fluidized bed with a central jet was considered; the simulations were performed in 2D and the schematic representation of the simulation domain is shown in Figure 5-1. The simulations were performed for three different central vertical jet velocities: 42 (base value), 46.2 (10% increase from base), and 50.4 (20% increase from base) m/s; considering the ratio of the size of the jet inlet with the bed width (i.e. 1/15), the above jet velocities correspond to superficial velocities of 2.8, 3.08, and 3.36 m/s, respectively. A summary of the simulation settings is provided in Table 5-2. According to this table, each simulation is performed for 20 s. In order to assure that this time is adequate, similar computations with a longer simulation time of 60 s were also performed. Both results were compared and no significant difference between the time-averaged values of the average particle height, wet coefficient of restitution, and normal relative collision velocity was found. For - 106 - example, by increasing the simulation time from 20 to 60 s, for case C3 with a superficial velocity of 2.8 m/s and coating viscosity of 10 mPa.s, the time-averaged average particle height remained unchanged, the time-averaged wet coefficient of restitution decreased from 0.39 to 0.37 (i.e. ~5% reduction), and the time-average normal relative collision velocity decreased from 0.101 to 0.100 m/s (i.e. ~1% reduction). In conclusion, a simulation time of 20 s is sufficient for the quantitative analyses performed in this work. In order to exclude the start-up effects of the fluidized bed (and avoiding large fluctuations of the above quantities), all the time-averaged values reported in this work are obtained between an initial time of 5 seconds (ti= 5 s) and a final time equal to the simulation time (tf = ts). This is because at the start-up of our simulations (here, from zero to 5 seconds) a substantial amount (around 50%) of the reported average particle heights lie beyond the “acceptable” range, with the acceptable range having upper and lower bounds equal to the average value plus the standard deviation and the average value minus the standard deviation, respectively. The above is less frequent for simulations times greater than 5 s (with about 25% of the data lying beyond the range); therefore; the first five seconds of each simulation are ignored in the reported time-averaged values in this work. Figure 5-1: Schematic representation of the 2D simulation domain. - 107 - Table 5-2. Summary of the simulation settings. Parameter Value Unit Number of solid particles (Np) 2400 - Number of cells in x-direction (Nx) 15 - Number of cells in y-direction (Ny) 45 - Simulation time (ts) 20 s Time step ( t∆ ) 5x10-4 s For the simulations performed in this work, 2400 particles are considered. For each simulation case, these particles are assumed to be pre-coated with identical coating thicknesses (h). In order to cover a reasonable range of coating cases, five different coating thicknesses are considered and are summarized in Table 5-3. Also for simplicity, it is assumed that the liquid coating is non- reactive and non-evaporative throughout each simulation and there is no liquid transfer due to particle-particle and particle-wall interactions. As a result of the above, the thickness of the liquid coating remains constant throughout each simulation. Table 5-3 also contains the moisture content of each particle which for our simulations is also the moisture content of the entire system, as the system consists of mono-sized spherical particles coated with identical liquid thicknesses. The moisture content ( liqx ) provided in Table 5-3 can be calculated from the following equation: l liq l s m x m m = + ( 5-11) where l l lm V ρ= and s s sm V ρ= are the masses of the liquid coating and solid particle, respectively. The densities of the liquid coating ( lρ ) and dry particle ( sρ ) are provided in Table 5-4 and the volume of the solid particle ( sV ) and liquid coating ( lV ) are determined as follows: 34 3s d V rpi= ( 5-12) 34 3l w s V r Vpi= − ( 5-13) - 108 - where dr and wr is the radii of the dry and wet particles, respectively, and are related with the following equation: w dr r h= + ( 5-14) Finally, the density of the wet particles ( wρ ), as given in Table 5-3, can be calculated as follows: w w w m V ρ = ( 5-15) where w l sm m m= + is the mass of the wet particle and 34 / 3w wV rpi= represents its volume. The radius and density of the particles, as given in Table 5-3, are used as the representative radius and density of the individual particles of each case in our DEM calculations. According to this table, the moisture contents of the systems considered in this work vary between 0.6 to 16.2 %. Note that the upper limit of the moisture content is close to the upper limit established by Flemmer [124], i.e. 13.6 % for randomly packed equal-sized spherical particles. According to Flemmer, for moisture contents greater than their indicated limit, the granular system transits from pendular to funicular regime. Table 5-4 contains the properties of the solid particle and liquid coating considered in this work that are used in the model for the wet coefficient of restitution. This model has been fully described in Chapter 2 of this thesis, and as explained in Section 5.2.3 of this chapter is used to model the interaction of wet particles. The model also requires the value of the half-separation distance (h1); in this work, it is assumed that it is equal to the thickness of the liquid coating (i.e. h1 = h). In order to investigate the effect of the liquid viscosity on the fluidized bed behaviour, the viscosity of the liquid coating is assumed to vary between 0 – 100 mPa.s (see Table 5-4). According to Table 5-4, the dry coefficient of restitution of particles is 0.95, which decreases due to the liquid effects and could be as low as 0 (meaning coalescence of particles); however, since the actual coalescence of particles is not considered in the current work, it is assumed that the minimum value of wet coefficient of restitution is equal to 0.001. - 109 - Table 5-3. Summary of the coating cases. Case No. 1 2 3 4 5 Coating thickness, mm 0.01 0.05 0.1 0.2 0.3 Moisture content, % 0.6 2.8 5.5 10.9 16.2 Radius of wet particle, mm 2.01 2.05 2.1 2.2 2.3 Density of wet particle, kg/m 2675 2580 2470 2280 2120 Table 5-4. Solid particle and liquid coating properties. Parameter Value Unit Radius of dry particle (rd ) 2 mm Density of dry particle ( sρ ) 2700 kg/m3 Density of coating liquid ( lρ ) 1000 kg/m3 Viscosity of coating liquid ( µ ) 0 - 100 mPa.s Surface tension of coating liquid (σ ) 72 /mN m Liquid-solid contact angle (θ ) 0 radian Average height of surface asperities (ha) 0.005 mm Dry coefficient of restitution (e) 0.95 - 5.3. RESULTS AND DISCUSSION 5.3.1. Bed behaviour Figure 5-2 shows the bed evolution of case C3 at fluidization velocity 2.8 m/s for an inviscid coating and a coating with viscosity 100 mPa.s. This figure illustrates the simulated instantaneous particle positions. According to this figure, there is an evident transition of the fluidization behaviour; for an inviscid liquid coating, due to the passage of gas, the bed expands and the solid particles are scattered everywhere in the bed, while as the viscosity of the liquid coating is increased (here to 100 mPa.s), the gas passes through the bed in the form of slugs. Similar transition of bed behaviour was observed for other simulation cases and other fluidization velocities considered in this work. For instance, Figure 5-3 shows the continuous bed evolution for case C3 at an increased fluidization velocity of 3.36 m/s; according to this figure, - 110 - consistent with our observation for a lower fluidization velocity, for an inviscid liquid coating the bed expands only while slugs are formed as the viscosity of liquid coating is increased to 100 mPa.s. The only major difference for the observations when having different fluidization velocities is that for the case with the higher superficial velocity, due to the additional amount of gas in the bed, more slugs are observed in the bed. The transition of the bed behaviour due to increased liquid coating viscosity can be explained by comparing the particle-particle collision related parameters; for instance, according to Figure 5-5 and Figure 5-6, the wet coefficient of restitution and the relative normal collision velocities are substantially reduced as the coating viscosity is increased from 0 to 100 mPa.s. In other words, as the coating viscosity is increased, particles tend to stay closer to each other; as a result, dense regions of particles are formed and gas passes through in the form of slugs. Further quantitative discussion regarding the effect of the liquid viscosity on the bed behaviour is presented later. The change in the fluidization behaviour due to the liquid effects, without immobilization of the system, has previously been confirmed both experimentally [125, 126] and numerically (via DEM simulations) [53] in the literature. According to van Buijtenen et al. [53], the addition of liquid results in particle collisions with low coefficient of restitution; this promotes the formation of a dense region, thus the passage of gas in the form of bubbles through the system. It is important to highlight that the proposed wet coefficient of restitution used in this work accounts for a number of key parameters such as liquid viscosity, relative normal collision velocity, etc. while the expression of wet coefficient of restitution used by van Buijtenen et al. depends only on the moisture content of particles. The experimental observations of Nagahashi et al. [125, 126] indicate fluidization enhancement with the addition of a small amount of liquid. This is explained by the reduction of the measured minimum spouting and fluidization velocities after the addition of the liquid to the system. The reported enhancement is significant for systems with relatively large (> 3 mm) and poorly wetted particles (i.e. with large liquid-solid contact angles) and when the added liquid has a density close to that of the solid phase. The observed enhancement of Nagahashi et al. cannot be predicted with our numerical tool. This is because for their case enhancement occurs when particles are poorly wetted while in all of our simulations, a well-dispersed liquid coating that perfectly wets the particles is considered. - 111 - a) b) Figure 5-2: Snapshots of the simulated transient bed evolution for case C3 at superficial velocity 2.8 m/s with: a) inviscid liquid coating, and b) coating with viscosity 100 mPa.s. - 112 - a) b) Figure 5-3: Snapshots of the simulated instantaneous particle positions for case C3 at superficial velocity 3.36 m/s with: a) inviscid liquid coating, and b) coating with viscosity 100 mPa.s. - 113 - 5.3.2. Average particle height In order to obtain the average particle height at time t ( ( )H t ), the y-coordinate of all the particles at that time were averaged, as follows: (i) 1 1( ) y ( )pN i p H t t N = = ∑ ( 5-16) Then, in order to obtain the time-averaged average particle height ( H ), the calculated average particle heights were averaged for the desired time duration. Figure 5-4 shows the time-averaged average particle height for cases C1 (left) and C3 (right) for different coating viscosities and different fluidization velocities. According to this figure, the average particle height decreases as the liquid viscosity increases; however, the average particle height increases as the superficial velocity increases. In order to explain the variations of the average particle height with respect to the liquid viscosity and fluidization velocity, the effect of these two parameters on the time- averaged wet coefficient of restitution and relative normal collision velocity is investigated, as shown in Figure 5-5 and Figure 5-6, respectively. According to these figures, at a constant superficial velocity, as the coating viscosity increases, the wet coefficient of restitution decreases and the relative normal collision velocity decreases; consequently, the average particle height decreases. As highlighted earlier, as a result, with increasing coating viscosity particles stay close, rather than being scattered, and the gas passes through in the form of slugs. Also, considering a constant viscosity, the wet coefficient of restitution and relative normal collision velocity increase as the superficial velocity increases; therefore, the average particle height increases. Similar observations were found for the rest of the coating cases considered in this work. - 114 - a) b) Figure 5-4: Time-averaged average particle height for different superficial velocities and different coating viscosities for: a) case C1, and b) case C3. - 115 - a) b) Figure 5-5: Time-averaged wet coefficient of restitution for different superficial velocities and different coating viscosities for: a) case C1, and b) case C3. - 116 - a) b) Figure 5-6: Time-averaged relative normal collision velocity for different superficial velocities and different coating viscosities for: a) case C1, and b) case C3. - 117 - 5.3.3. Wet coefficient of restitution Figure 5-7 summarizes the variations of the time-averaged wet coefficient of restitution with respect to the coating viscosity and coating thickness for a superficial velocity of 2.8 m/s. It is important to highlight that the time-averaged values of the wet coefficient of restitution are being reported in this figure; this is because the proposed wet coefficient of restitution in this work is also a function of the (time-dependant) normal collision velocity between the particles; hence, although particles have a constant coating during each simulation, the wet coefficient of restitution will be a time-dependant quantity. According to this figure, the wet coefficient of restitution decreases as the coating viscosity and thickness increase. This figure shows that at low coating viscosities (here, 1µ = mPa.s), the wet coefficient of restitution is very close to the prescribed dry coefficient of restitution (e = 0.95); however, it decreases substantially as the coating viscosity increases to 10 mPa.s and finally approaches zero as the liquid viscosity is further increased to 100 mPa.s. The figure also indicates the relationship between the coating thickness and wet coefficient of restitution; according to Figure 5-7, there is an obvious difference between case 1 (with the thinnest coating thickness) and other cases considered in this work. That is in fact expected as the thickness of case 1 is at least 5 times less than the coating thickness of the other cases and consequently, the liquid coating effects for case 1 are less significant than for the other cases. It also important to highlight that the observed wet coefficient of restitution also varies spatially within the system; this is because in the proposed model the wet coefficient of restitution is a function of interparticle collision velocity between the particles, which varies spatially within the fluidized bed. Figure 5-8 summarizes the variations of the time-averaged relative normal collision velocity with respect to the coating viscosity and coating thickness at fluidization velocity 2.8 m/s. As expected, the trend is closely related to that of Figure 5-7, i.e. the relative normal collision velocity decreases as the coating viscosity and thickness increase; also, there is an obvious difference between case 1 and other cases due to its substantially thinner coating. Similar trends were found for higher superficial velocities. - 118 - Figure 5-7: Summary of the variations of the time-averaged wet coefficient of restitution with respect to coating viscosity and coating thickness at superficial velocity of 2.8 m/s. Figure 5-8: Summary of the variations of the time-averaged normal relative collision velocity with respect to coating viscosity and coating thickness at superficial velocity of 2.8 m/s. - 119 - 5.4. SUMMARY AND CONCLUSIONS In this chapter, the impact of the liquid coating on the fluidized bed behaviour was studied by performing DEM investigations. It was assumed that the liquid coating affects the system behaviour by influencing the coefficient of restitution used to determine the outcome of particle- particle interactions. The proposed wet coefficient of restitution model, as shown in Chapter 2 of this thesis, was implemented in an open-source numerical tool. The model accounts for the kinetic energy of particles, the viscous and capillary contribution of liquid, and the inelasticity of particles. The modified numerical tool was then used to perform DEM simulations of a bed of mono-sized solid spherical particles pre-coated with identical liquid coatings. For the simulations performed in this work, three different fluidization velocities and five different coating thicknesses were considered, and the coating viscosity varied between 0 – 100 mPa.s. The simulation results were analyzed by observing the snapshots of the instantaneous particle positions and the time-averaged values of the average particle height, wet coefficient of restitution, and relative normal collision velocity. It was shown that as the coating viscosity increases, the fluidization behaviour changes where the fluidization gas passes through the system in the form of slugs; more slugs were formed as the fluidization velocity increased. Moreover, it was discussed that as the coating viscosity increases, the wet coefficient of restitution and the relative normal collision velocities decrease and consequently, the average particle height decreases. Conversely, it was shown that as the fluidization velocity increases, the wet coefficient of restitution and the relative normal collision velocity increase, leading to increased average particle height. Finally, it was observed that the wet coefficient of restitution and relative coating velocity decrease as the coating thickness increases. This is expected, as an increase in the coating thickness means more liquid effects. - 120 - Chapter 6. SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS 6.1. SUMMARY AND CONCLUSIONS Mathematical modeling has been used to study the interaction of wet particles. First, the binary interaction of particles is studied; then, the interaction of multiple particles in a fluidized bed is investigated. 6.1.1. Binary interaction of particles 6.1.1.1. Collision model A model is developed to determine the collision outcome of equal-sized wet particles coated with thin liquid films. This is achieved by introducing a wet coefficient of restitution that accounts for the liquid effects along with the inelasticity of dry particles. The key conclusions related to this part are as follows: • The proposed model successfully determines the critical collision velocity, below which particles coalesce, and above which they rebound. The results of the developed model are compared with the numerical results of the equation of motion and excellent agreement is found. This means that the proposed model is appropriate for determining the collision outcome of wet particles and the numerical solution of the equation of motion is no longer required within the indicated range. Moreover, model predictions are compared with the available experimental data in the literature and good agreement is found. • With the aid of the proposed model, the boundary between collisions with dominant capillary or viscous effects is determined. It is shown that for collisions with small capillary numbers, the surface tension effects dominate, while the viscous effects are dominant at large capillary numbers. For intermediate capillary numbers, both of these effects are important and should be considered. - 121 - • The proposed model has several advantages; for instance, it can be implemented in numerical tools as a sub-model to determine the collision outcome of wet particles. Also, it is more general compared to the other available models in the literature, for example. However, currently it is formulated for equal-sized spherical particles with identical liquid coatings and needs further modifications in order to be applicable to particles with different sizes and/or different liquid coatings. 6.1.1.2. Application of collision model to fluid cokers The above collision model is applied to fluid cokers by considering the interaction of bitumen- coated coke particles under practical operating conditions. In order to consider the actual coking conditions, the temperature- and time-dependant variations of the bitumen thickness and viscosity are incorporated into the model. The main conclusions of this work are as follows: • For bitumen-coated coke particles, the collision outcome is strongly affected by the reaction time. Due to the chemical reactions, the bitumen viscosity increases while its thickness decreases. Because these two effects are opposite, a peak critical velocity is observed at an intermediate reaction time. • By comparing the peak value of the predicted critical velocities with the characteristic inter- particle collision velocity within real-scale fluid cokers, the maximum possible size of agglomerates is estimated. The estimations suggest the formation of large agglomerates consisting up to a few hundred particles within the reactor. Such a prediction is in fair agreement with plant observations. 6.1.1.3. Rupture of liquid bridges A simplified mathematical model and numerical solution of the governing equations are developed to study the evolution of stretching liquid bridges between equal-sized solid spherical - 122 - particles. The simplified model is a geometric representation of the liquid bridge in which the gas-liquid interface is represented with a parabola. For the numerical simulations, the governing Navier-Stokes equations are solved using the Volume of Fluid (VOF) model in a Computational Fluid dynamics (CFD) tool – FLUENT. The key conclusions regarding this part of the work are highlighted as follows: • The simplified model has limitations; it is only suitable for problems with dominant surface tension effects, i.e. small capillary, Weber and Bond numbers. Also, it is inappropriate for liquid bridges between particles with different contact angles. • Numerical solutions of the Navier-Stokes equations explicitly account for the key properties related to the problem. This includes the liquid surface tension, liquid viscosity, liquid density, and the gravitational acceleration. Hence, a wider range of problems can be explored with the numerical simulations of the Navier-Stokes equations. In addition, numerical simulations, unlike the simplified model, are not restricted to the liquid bridge problems between particles with identical contact angle and negligible gravitational effects. As a result, the numerical simulations are used to investigate the variations of the rupture distance and liquid transfer fraction with respect to the capillary, Weber, and Bond numbers between particles with identical or different contact angles. 6.1.2. DEM investigations of wet fluidized beds In order to investigate the interaction of multiple wet particles, an open-source numerical tool – MFIX - is modified. The code is originally based on a soft-sphere Discrete Element Approach (DEM), accounting for multiple particle interactions, and is capable of simulating fluidized beds. In order to include the liquid effects, the wet coefficient of restitution model proposed in this thesis is implemented into the code. Then, the developed DEM tool is used to investigate the impact of liquid coating on the fluidization behaviour. The main conclusions of this work are as follows: - 123 - • The simulation results show that as the coating viscosity increases, particles tend to stay close together. Hence the fluidization behaviour changes and the air passes through the system in the form of slugs. • It is also shown that as the coating viscosity increases the time-average wet coefficient of restitution, normal relative collision velocities, and average particle height (i.e. bed centroid in the y-direction) decrease. • The effect of increasing the coating thickness is similar to the above, i.e. it results in a reduction of the time-averaged wet coefficient of restitution and normal relative collision velocities. 6.2. CONTRIBUTIONS • A model for determining the binary collision outcome of wet particles has been developed that includes both capillary and viscous effects. The model is written in the form a wet coefficient of restitution. Therefore, unlike most of the models available in the literature, the proposed collision model allows for determining both the coalescence and rebound outcomes of wet particles. • Based on the developed collision model, a simple mathematical model is proposed that can be used to determine the agglomeration tendency of bitumen-coated coke particles inside fluid cokers. • Stretching liquid bridges between solid spherical particles is investigated by numerically solving the governing Navier-Stokes equations. The model results are used to assess the - 124 - limitations of the simplified model. The conclusions can guide researchers to choose the proper model depending on the conditions of the rupture problem. • A modified DEM numerical tool is developed by incorporating the wet coefficient of restitution introduced in this thesis. The implemented model accounts for the main liquid properties, e.g. its viscosity or surface tension; therefore, the modified tool can be used to perform more comprehensive studies of fluidized beds in the presence of liquid. The proposed modification in this thesis is an efficient solution for accounting for the liquid properties in DEM simulations. This is because the current treatment requires only modifying the coefficient of restitution between particles in order to consider the liquid effects, and there is no need to complicate the coding process by modifying the particle-particle contact forces by including the liquid forces in the code. 6.3. RECOMMENDATIONS A number of future research topics are suggested below: Interaction of multiple particles: The proposed collision model and presented liquid evolution simulation results in this thesis are currently limited to the binary collisions. Extending them for the collision of multiple particles can broaden their applicability to practical problems. Further validation: Numerical simulations have been used to investigate the variations of the rupture distance and liquid distribution of stretching liquid bridges under a variety conditions. Due to the lack of experimental data, the above results cannot be fully validated. A systematic experimental study is necessary for further validation of such numerical results. Also, the DEM investigations presented in this work are not compared with quantitative experimental data; obtaining detailed experimental data for validating DEM simulations is not an easy task, but could be very useful for both validation and further model development. - 125 - Improvement of the DEM simulations: In the simulations performed in this thesis, the liquid transfer due to the particle-particle collisions is ignored; hence, the coating thickness of each individual particle remains unchanged. Liquid transfer could play a significant role in some fluidization systems. In addition, the coalescence of particles and breakage of the resulting granules should be implemented in the current DEM simulations. That would result in a much more powerful tool for the investigation of agglomerate formation in fluidized beds. Extension of the current models: the models presented in this thesis have some limitations that can be removed. For instance, in this work particles are assumed to be equal-sized. Therefore, one may consider extending the proposed collision model, the numerical simulations of the rupture phenomenon, and DEM investigations of wet fluidized beds for systems consisting of particles of different sizes in future. - 126 - BIBLIOGRAPHY [1] Gray, M. R., 2002, "Fundamentals of Bitumen Coking Processes Analogous to Granulations: A Critical Review," The Canadian Journal of Chemical Engineering, 80(3) pp. 393-401. [2] Dunlop, D. D., Griffin, L. I., and Moser, J. F., 1958, "Particle Size Control in Fluid Coking," Chem Eng Progr, 54(8) pp. 39-43. [3] Knapper, B., Gray, M. R., Chan, E. W., and Mikula, R., 2003, "Measurement of Efficiency of Distribution of Liquid Feed in a Gas-Solid Fluidized Bed Reactor," International Journal of Chemical and Reactor Engineering, 1pp. A35. [4] Savage, S. B., and Jeffrey, D. J., 1981, "The Stress Tensor in a Granular Flow at High Shear Rates," Journal of Fluid Mechanics, 110pp. 255-272. [5] Campbell, C. S., and Brennen, C. E., 1985, "Computer Simulation of Granular Shear Flows," Journal of Fluid Mechanics, 151pp. 167-188. [6] Cundall, P. A., and Strack, O. D. L., 1979, "A Discrete Numerical Model for Granular Assemblies," Geotechnique, 29(1) pp. 47-65. [7] Hoomans, B. P. B., Kuipers, J. A. M., Briels, W. J., and vanSwaaij, W. P. M., 1996, "Discrete Particle Simulation of Bubble and Slug Formation in a Two-Dimensional Gas- Fluidised Bed: A Hard-Sphere Approach," Chemical Engineering Science, 51(1) pp. 99-118. [8] Tsuji, Y., Kawaguchi, T., and Tanaka, T., 1993, "Discrete Particle Simulation of Two- Dimensional Fluidized Bed," Powder Technology, 77(1) pp. 79-87. [9] Simons, S. J. R., and Fairbrother, R. J., 2000, "Direct Observations of Liquid Binder-Particle Interactions: The Role of Wetting Behaviour in Agglomerate Growth," Powder Technology, 110(1) pp. 44-58. [10] Rossetti, D., and Simons, S. J. R., 2003, "A Microscale Investigation of Liquid Bridges in the Spherical Agglomeration Process," Powder Technology, 130(1-3) pp. 49-55. [11] Seville, J. P. K., Willett, C. D., and Knight, P. C., 2000, "Interparticle Forces in Fluidisation: A Review," Powder Technology, 113(3) pp. 261-268. - 127 - [12] Fisher, R. A., 1926, "On the Capillary Forces in an Ideal Soil; Correction of Formulae Given by WB Haines," Journal of Agricultural Science, 16pp. 492–505. [13] Heady, R. B., and Cahn, J. W., 1970, "Analysis of the Capillary Forces in Liquid-Phase Sintering of Spherical Particles," Metallurgical Transactions, 1(1) pp. 185-189. [14] Orr, F. M., Scriven, L. E., and Rivas, A. P., 1975, "Pendular Rings between Solids: Meniscus Properties and Capillary Force," Journal of Fluid Mechanics, 67(04) pp. 723-742. [15] Pitois, O., Moucheront, P., and Chateau, X., 2000, "Liquid Bridge between Two Moving Spheres: An Experimental Study of Viscosity Effects," Journal of Colloid and Interface Science, 231(1) pp. 26-31. [16] Rumpf, H., 1962, "The Strength of Granules and Agglomerates," Agglomeration, Interscience, New York, pp. 379-414. [17] Rynhart, P., McKibbin, R., McLachlan, R., and Jones, J. R., 2002, "Mathematical Modelling of Granulation: Static and Dynamic Liquid Bridges," Res. Lett. Inf. Math. Sci., 3pp. 199-212. [18] Pierrat, P., and Caram, H. S., 1997, "Tensile Strength of Wet Granula Materials," Powder Technology, 91(2) pp. 83-93. [19] Chan, D. Y. C., and Horn, R. G., 1985, "The Drainage of Thin Liquid Films between Solid Surfaces," The Journal of Chemical Physics, 83pp. 5311. [20] Ennis, B. J., Li, J., Tardos, G. I., and Pfeffer, R., 1990, "The Influence of Viscosity on the Strength of an Axially Strained Pendular Liquid Bridge," Chemical Engineering Science, 45pp. 3071-3088. [21] Ennis, B. J., Tardos, G., and Pfeffer, R., 1991, "A Microlevel-Based Characterization of Granulation Phenomena," Powder Technology, 65(1-3) pp. 257-272. [22] Matthewson, M. J., 1988, "Adhesion of Spheres by Thin Liquid Films," Philosophical Magazine A, 57(2) pp. 207-216. [23] McFarlane, J. S., and Tabor, D., 1950, "Relation between Friction and Adhesion," Proceedings of the Royal Society of London.Series A, Mathematical and Physical Sciences (1934-1990), 202(1069) pp. 244-253. - 128 - [24] Reynolds, O., 1886, "On the Theory of Lubrication and its Application to Mr. Beauchamp Tower’s Experiments, Including an Experimental Determination of the Viscosity of Olive Oil," Philosophical Transactions of the Royal Society of London, 177pp. 157–234. [25] Iveson, S. M., 2002, "Limitations of One-Dimensional Population Balance Models of Wet Granulation Processes," Powder Technology, 124(3) pp. 219-229. [26] Pitois, O., Moucheront, P., and Chateau, X., 2001, "Rupture Energy of a Pendular Liquid Bridge," The European Physical Journal B-Condensed Matter, 23(1) pp. 79-86. [27] Rossetti, D., Pepin, X., and Simons, S. J. R., 2003, "Rupture Energy and Wetting Behavior of Pendular Liquid Bridges in Relation to the Spherical Agglomeration Process," Journal of Colloid and Interface Science, 261(1) pp. 161-169. [28] Simons, S. J. R., Seville, J. P. K., and Adams, M. J., 1994, "An Analysis of the Rupture Energy of Pendular Liquid Bridges," Chemical Engineering Science, 49(14) pp. 2331-2339. [29] Willett, C. D., Zhang, Z., and Seville, J. P. K., 1998, "An Examination of the Effects of Contact Angle on the Properties of Exact and Toroidal Liquid Bridges," Proc.World Congress on Particle Technology, 3. [30] Adams, M. J., and Perchard, V., 1985, "The Cohesive Forces between Particles with Interstitial Liquid," Institute of Chemical Engineering Symposium, 91pp. 147-160. [31] Mazzone, D. N., Tardos, G. I., and Pfeffer, R., 1987, "The Behavior of Liquid Bridges between Two Relatively Moving Particles," Powder Technology, 51(1) pp. 71-83. [32] Barnocky, G., and Davis, R. H., 1988, "Elastohydrodynamic Collision and Rebound of Spheres: Experimental Verification," Physics of Fluids, 31pp. 1324. [33] Braumann, A., Goodson, M. J., Kraft, M., and Mort, P. R., 2007, "Modelling and Validation of Granulation with Heterogeneous Binder Dispersion and Chemical Reaction," Chemical Engineering Science, 62(17) pp. 4717-4728. [34] Litster, J. D., and Sarwono, R., 1996, "Fluidized Drum Granulation: Studies of Agglomerate Formation," Powder Technology, 88(2) pp. 165-172. - 129 - [35] Abberger, T., 2001, "The Effect of Powder Type, Free Moisture and Deformation Behaviour of Granules on the Kinetics of Fluid-Bed Granulation," European Journal of Pharmaceutics and Biopharmaceutics, 52(3) pp. 327-336. [36] Abberger, T., Seo, A., and Schæfer, T., 2002, "The Effect of Droplet Size and Powder Particle Size on the Mechanisms of Nucleation and Growth in Fluid Bed Melt Agglomeration," International Journal of Pharmaceutics, 249(1-2) pp. 185-197. [37] Dewettinck, K., and Huyghebaert, A., 1998, "Top-Spray Fluidized Bed Coating: Effect of Process Variables on Coating Efficiency," LWT-Food Science and Technology, 31(6) pp. 568- 575. [38] Bouffard, J., Kaster, M., and Dumont, H., 2005, "Influence of Process Variable and Physicochemical Properties on the Granulation Mechanism of Mannitol in a Fluid Bed Top Spray Granulator," Drug Development and Industrial Pharmacy, 31(9) pp. 923-933. [39] Saleh, K., Steinmetz, D., and Hemati, M., 2003, "Experimental Study and Modeling of Fluidized Bed Coating and Agglomeration," Powder Technology, 130(1-3) pp. 116-123. [40] Saleh, K., and Guigon, P., 2007, "Influence of Wetting Parameters on Particle Growth in Fluidized-Bed Coating and Agglomeration Processes," Particle & Particle Systems Characterization, 24(2) pp. 136-143. [41] Saleh, K., Cherif, R., and Hemati, M., 1999, "An Experimental Study of Fluidized-Bed Coating: Influence of Operating Conditions on Growth Rate and Mechanism," Advanced Powder Technology, 10(3) pp. 255–277. [42] Rajniak, P., Mancinelli, C., Chern, R. T., Stepanek, F., Farber, L., and Hill, B. T., 2007, "Experimental Study of Wet Granulation in Fluidized Bed: Impact of the Binder Properties on the Granule Morphology," International Journal of Pharmaceutics, 334(1-2) pp. 92-102. [43] Rajniak, P., Stepanek, F., Dhanasekharan, K., Fan, R., Mancinelli, C., and Chern, R. T., 2009, "A Combined Experimental and Computational Study of Wet Granulation in a Wurster Fluid Bed Granulator," Powder Technology, 189(2) pp. 190-201. [44] McDougall, S., Saberian, M., Briens, C., Berruti, F., and Chan, E., 2005, "Effect of Liquid Properties on the Agglomerating Tendency of a Wet gas–solid Fluidized Bed," Powder Technology, 149(2-3) pp. 61-67. - 130 - [45] Weber, S., Briens, C., Berruti, F., Chan, E., and Gray, M., 2006, "Agglomerate Stability in Fluidized Beds of Glass Beads and Silica Sand," Powder Technology, 165(3) pp. 115-127. [46] Weber, S., Briens, C., Berruti, F., Chan, E., and Gray, M., 2008, "Effect of Agglomerate Properties on Agglomerate Stability in Fluidized Beds," Chemical Engineering Science, 63(17) pp. 4245-4256. [47] Fan, X., Yang, Z., Parker, D. J., Ng, B., Ding, Y., and Ghadiri, M., 2009, "Impact of Surface Tension and Viscosity on Solids Motion in a Conical High Shear Mixer Granulator," AIChE Journal, 55(12) pp. 3088-3098. [48] Liu, L. X., Litster, J. D., Iveson, S. M., and Ennis, B. J., 2000, "Coalescence of Deformable Granules in Wet Granulation Processes," AIChE Journal, 46(3) pp. 529-539. [49] Goldschmidt, M. J. V., Weijers, G. G. C., Boerefijn, R., and Kuipers, J. A. M., 2003, "Discrete Element Modelling of Fluidised Bed Spray Granulation," Powder Technology, 138(1) pp. 39-45. [50] Goldschmidt, M. J. V., Beetstra, R., and Kuipers, J. A. M., 2004, "Hydrodynamic Modelling of Dense Gas-Fluidised Beds: Comparison and Validation of 3D Discrete Particle and Continuum Models," Powder Technology, 142(1) pp. 23-47. [51] Link, J. M., Godlieb, W., Deen, N. G., and Kuipers, J. A. M., 2007, "Discrete Element Study of Granulation in a Spout-Fluidized Bed," Chemical Engineering Science, 62(1-2) pp. 195-207. [52] Link, J. M., Godlieb, W., Tripp, P., Deen, N. G., Heinrich, S., Kuipers, J. A. M., Schönherr, M., and Peglow, M., 2009, "Comparison of Fibre Optical Measurements and Discrete Element Simulations for the Study of Granulation in a Spout Fluidized Bed," Powder Technology, 189(2) pp. 202-217. [53] van Buijtenen, M. S., Deen, N. G., Heinrich, S., Antonyuk, S., and Kuipers, J., 2009, "A Discrete Element Study of Wet Particle-Particle Interaction during Granulation in a Spout Fluidized Bed," The Canadian Journal of Chemical Engineering, 87(2) pp. 308-317. [54] Tafreshi, Z. M., Kirpalani, D., Bennett, A., and McCracken, T. W., 2002, "Improving the Efficiency of Fluid Cokers by Altering Two-Phase Feed Characteristics," Powder Technology, 125(2) pp. 234-241. - 131 - [55] Hulet, C., Briens, C., Berruti, F., Chan, E. W., and Ariyapadi, S., 2003, "Entrainment and Stability of a Horizontal Gas-Liquid Jet in a Fluidized Bed," International Journal of Chemical and Reactor Engineering, 1(A60). [56] Briens, C., McDougall, S., and Chan, E., 2003, "On-Line Detection of Bed Fluidity in a Fluidized Bed Coker," Powder Technology, 138(2-3) pp. 160-168. [57] Ariyapadi, S., Holdsworth, D. W., Norley, C. J. D., Berruti, F., and Briens, C., 2003, "Digital X-Ray Imaging Technique to Study the Horizontal Injection of Gas-Liquid Jets into Fluidized Beds," International Journal of Chemical and Reactor Engineering, 1(A56). [58] Song, X., Bi, H., JimLim, C., Grace, J., Chan, E., Knapper, B., and McKnight, C., 2004, "Hydrodynamics of the Reactor Section in Fluid Cokers," Powder Technology, 147(1-3) pp. 126-136. [59] House, P. K., Saberian, M., Briens, C. L., Berruti, F., and Chan, E., 2004, "Injection of a Liquid Spray into a Fluidized Bed: Particle-Liquid Mixing and Impact on Fluid Coker Yields," Industrial & Engineering Chemistry Research, 43(18) pp. 5663-5669. [60] Ariyapadi, S., Berruti, F., Briens, C., McMillan, J., and Zhou, D., 2004, "Horizontal Penetration of Gas-Liquid Spray Jets in Gas-Solid Fluidized Beds," International Journal of Chemical Reactor Engineering, 2(A22). [61] McMillan, J., Zhou, D., Ariyapadi, S., Briens, C., Berruti, F., and Chan, E., 2005, "Characterization of the Contact between Liquid Spray Droplets and Particles in a Fluidized Bed," Industrial & Engineering Chemistry Research, 44(14) pp. 4931-4939. [62] Ariyapadi, S., Berruti, F., Briens, C., Knapper, B., Skwarok, R., and Chan, E., 2005, "Stability of Horizontal gas–liquid Sprays in Open-Air and in a gas–solid Fluidized Bed," Powder Technology, 155(3) pp. 161-174. [63] Bi, H. T., Grace, J. R., Lim, C. J., Rusnell, D., Bulbuc, D., and McKnight, C. A., 2005, "Hydrodynamics of the Stripper Section of Fluid Cokers," The Canadian Journal of Chemical Engineering, 83(2) pp. 161-168. [64] Leach, A., Portoghese, F., Briens, C., and Berruti, F., 2008, "A New and Rapid Method for the Evaluation of the liquid–solid Contact Resulting from Liquid Injection into a Fluidized Bed," Powder Technology, 184(1) pp. 44-51. - 132 - [65] House, P. K., Briens, C. L., Berruti, F., and Chan, E., 2008, "Effect of Spray Nozzle Design on liquid–solid Contact in Fluidized Beds," Powder Technology, 186(1) pp. 89-98. [66] Li, T., Pougatch, K., Salcudean, M., and Grecov, D., 2008, "Numerical Simulation of Horizontal Jet Penetration in a Three-Dimensional Fluidized Bed," Powder Technology, 184(1) pp. 89-99. [67] Li, T., Pougatch, K., Salcudean, M., and Grecov, D., 2009, "Numerical Simulation of Single and Multiple Gas Jets in Bubbling Fluidized Beds," Chemical Engineering Science, 64(23) pp. 4884-4898. [68] Nowak, P., Pougatch, K., and Salcudean, M., 2010, "Simulation of the Steam-Bitumen Jet in the Fluidized Bed of Coke Particles using the Eulerian-Lagrangian Splitting Method," The Journal of Computational Multiphase Flows, 2(1) pp. 59-72. [69] Pougatch, K., and Salcudean, M., 2010, "Numerical Simulation of Liquid Spray Dispersion in a Fluidized Bed," The Canadian Journal of Chemical Engineering, 88(4) pp. 648-654. [70] Cui, M., Straatman, A. G., and Zhang, C., 2005, "A Computational Study of Gas-Solid Flow in an Enhanced Solid Entrainment (ESE) Nozzle System," International Journal of Chemical and Reactor Engineering, 3(A7). [71] Ariyapadi, S., McMillan, J., Zhou, D., Berruti, F., Briens, C., and Chan, E., 2005, "Modeling the Mixing of a gas–liquid Spray Jet Injected in a gas–solid Fluidized Bed: The Effect of the Draft Tube," Chemical Engineering Science, 60(21) pp. 5738-5750. [72] Zhou, D., Saberian, M., Briens, C., Berruti, F., Chan, E. W., and McDougall, S. L., 2004, "Enhancement of the Distribution of a Liquid Sprayed into a Fluidized Bed," International Journal of Chemical Reactor Engineering, 2(S2). [73] Gray, M. R., Le, T., McCaffrey, W. C., Berruti, F., Soundararajan, S., Chan, E., Huq, I., and Thorne, C., 2001, "Coupling of Mass Transfer and Reaction in Coking of Thin Films of an Athabasca Vacuum Residue," Industrial & Engineering Chemistry Research, 40(15) pp. 3317- 3324. [74] Gray, M. R., McCaffrey, W. C., Huq, I., and Le, T., 2004, "Kinetics of Cracking and Devolatilization during Coking of Athabasca Residues," Industrial & Engineering Chemistry Research, 43(18) pp. 5438-5445. - 133 - [75] Aminu, M. O., Elliott, J. A. W., McCaffrey, W. C., and Gray, M. R., 2004, "Fluid Properties at Coking Process Conditions," Industrial & Engineering Chemistry Research, 43(12) pp. 2929- 2935. [76] Gray, M. R., Le, T., and Wu, X. A., 2007, "Role of Pressure in Coking of Thin Films of Bitumen," Canadian Journal of Chemical Engineering, 85(5) pp. 773-780. [77] Radmanesh, R., Chan, E., and Gray, M. R., 2008, "Modeling of Mass Transfer and Thermal Cracking during the Coking of Athabasca Residues," Chemical Engineering Science, 63(6) pp. 1683-1691. [78] Haines, W. B., 1925, "A Note on the Physical Properties of Soils," Journal of Agricultural Science, 17pp. 264-290. [79] Fowle, A. A., Wang, C. A., and Strong, P. F., 1979, "Experiments on the Stability of Conical and Cylindrical Liquid Columns at Low Bond Numbers," In Proc. 3rd European Symp. Mat. Sci. Space, pp. 317-325. [80] Mehrotra, V. P., and Sastry, K. V. S., 1980, "Pendular Bond Strength between Unequal- Sized Spherical Particles," Powder Technology, 25(2) pp. 203-214. [81] Lian, G., Thornton, C., and Adams, M. J., 1993, "A Theoretical Study of the Liquid Bridge Forces between Two Rigid Spherical Bodies," Journal of Colloid and Interface Science, 161(1) pp. 138-147. [82] Mikami, T., Kamiya, H., and Horio, M., 1998, "Numerical Simulation of Cohesive Powder Behavior in a Fluidized Bed," Chemical Engineering Science, 53(10) pp. 1927-1940. [83] Soulie, F., Cherblanc, F., El Youssoufi, M. S., and Saix, C., 2006, "Influence of Liquid Bridges on the Mechanical Behaviour of Polydisperse Granular Materials," International Journal for Numerical and Analytical Methods in Geomechanics, 30(3) pp. 213-228. [84] De Bisschop, F. R. E., and Rigole, W. J. L., 1982, "A Physical Model for Liquid Capillary Bridges between Adsorptive Solid Spheres: The Nodoid of Plateau," Journal of Colloid and Interface Science, 88(1) pp. 117-128. - 134 - [85] Mazzone, D. N., Tardos, G. I., and Pfeffer, R., 1986, "The Effect of Gravity on the Shape and Strength of a Liquid Bridge between Two Spheres," Journal of Colloid and Interface Science, 113(2) pp. 544-556. [86] Dai, Z., and Lu, S., 1998, "Liquid Bridge Rupture Distance Criterion between Spheres," International Journal of Mineral Processing, 53(3) pp. 171-181. [87] Pepin, X., Rossetti, D., Iveson, S. M., and Simons, S. J. R., 2000, "Modeling the Evolution and Rupture of Pendular Liquid Bridges in the Presence of Large Wetting Hysteresis," Journal of Colloid and Interface Science, 232(2) pp. 289-297. [88] Pepin, X., Rossetti, D., and Simons, S. J. R., 2000, "Modeling Pendular Liquid Bridges with a Reducing Solid–Liquid Interface," Journal of Colloid and Interface Science, 232(2) pp. 298- 302. [89] Shi, D., and McCarthy, J. J., 2008, "Numerical Simulation of Liquid Transfer between Particles," Powder Technology, 184(1) pp. 64-75. [90] Concus, P., and Finn, R., 1998, "Discontinuous Behavior of Liquids between Parallel and Tilted Plates," Physics of Fluids, 10pp. 39. [91] Meseguer, J., Slobozhanin, L. A., and Perales, J. M., 1995, "A Review on the Stability of Liquid Bridges," Advances in Space Research, 16(7) pp. 5-14. [92] Haynes, J. M., 1970, "Stability of a Fluid Cylinder," Journal of Colloid and Interface Science, 32(4) pp. 652-654. [93] Gillette, R. D., and Dyson, D. C., 1971, "Stability of Fluid Interfaces of Revolution between Equal Solid Circular Plates," Chemical Engineering Journal, 2(1) pp. 44-54. [94] Coriell, S. R., Hardy, S. C., and Cordes, M. R., 1977, "Stability of Liquid Zones," Journal of Colloid and Interface Science, 60(1) pp. 126-136. [95] Sanz, A., 1985, "The Influence of the Outer Bath in the Dynamics of Axisymmetric Liquid Bridges," Journal of Fluid Mechanics, 156pp. 101-140. [96] Zhang, Y., and Alexander, J. I. D., 1990, "Sensitivity of Liquid Bridges Subject to Axial Residual Acceleration," Physics of Fluids A: Fluid Dynamics, 2pp. 1966-1974. - 135 - [97] Tsamopoulos, J., Chen, T. Y., and Borkar, A., 1992, "Viscous Oscillations of Capillary Bridges," Journal of Fluid Mechanics, 235pp. 579-609. [98] Mollot, D. J., Tsamopoulos, J., Chen, T. Y., and Ashgriz, N., 1993, "Nonlinear Dynamics of Capillary Bridges: Experiments," Journal of Fluid Mechanics, 255pp. 411-435. [99] Zhang, X., Padgett, R. S., and Basaran, O. A., 1996, "Nonlinear Deformation and Breakup of Stretching Liquid Bridges," Journal of Fluid Mechanics, 329pp. 207-245. [100] Erle, M. A., Dyson, D. C., and Morrow, N. R., 1971, "Liquid Bridges between Cylinders, in a Torus, and between Spheres," AIChE Journal, 17(1) pp. 115-121. [101] Eggers, J., and Dupont, T. F., 1994, "Drop Formation in a One-Dimensional Approximation of the Navier–Stokes Equation," Journal of Fluid Mechanics, 262pp. 205-221. [102] Papageorgiou, D. T., 1995, "On the Breakup of Viscous Liquid Threads," Physics of Fluids, 7pp. 1529-1544. [103] Zhang, X., 1999, "Dynamics of Drop Formation in Viscous Flows," Chemical Engineering Science, 54(12) pp. 1759-1774. [104] Huang, W., Lee, S., Sung, H. J., Lee, T., and Kim, D., 2008, "Simulation of Liquid Transfer between Separating Walls for Modeling Micro-Gravure-Offset Printing," International Journal of Heat and Fluid Flow, 29(5) pp. 1436-1446. [105] Dodds, S., da Silveira Carvalho, M., and Kumar, S., 2009, "Stretching and Slipping of Liquid Bridges Near Plates and Cavities," Physics of Fluids, 21pp. 092103-1-092103-15. [106] Passos, M. L., and Mujumdar, A. S., 2000, "Effect of Cohesive Forces on Fluidized and Spouted Beds of Wet Particles," Powder Technology, 110(3) pp. 222-238. [107] Vieira, M. G. A., and Rocha, S. C. S., 2004, "Influence of the Liquid Saturation Degree on the Fluid Dynamics of a Spouted-Bed Coater," Chemical Engineering and Processing, 43(10) pp. 1275-1280. [108] McLaughlin, L. J., and Rhodes, M. J., 2001, "Prediction of Fluidized Bed Behaviour in the Presence of Liquid Bridges," Powder Technology, 114(1-3) pp. 213-223. - 136 - [109] Muguruma, Y., Tanaka, T., and Tsuji, Y., 2000, "Numerical Simulation of Particulate Flow with Liquid Bridge between Particles (Simulation of Centrifugal Tumbling Granulator)," Powder Technology, 109(1-3) pp. 49-57. [110] Yang, S. C., and Hsiau, S. S., 2001, "The Simulation of Powders with Liquid Bridges in a 2D Vibrated Bed," Chemical Engineering Science, 56(24) pp. 6837-6849. [111] Jain, K., Shi, D., and McCarthy, J. J., 2004, "Discrete Characterization of Cohesion in gas– solid Flows," Powder Technology, 146(1-2) pp. 160-167. [112] Fairbrother, R. J., and Simons, S. J. R., 1998, "Modelling of Binder-Induced Agglomeration," Particle & Particle Systems Characterization, 15(1) pp. 16-20. [113] Goldschmidt, M. J. V., Beetstra, R., and Kuipers, J. A. M., 2002, "Hydrodynamic Modelling of Dense Gas-Fluidised Beds: Comparison of the Kinetic Theory of Granular Flow with 3D Hard-Sphere Discrete Particle Simulations," Chemical Engineering Science, 57(11) pp. 2059-2075. [114] Mason, G., and Clark, W. C., 1965, "Liquid Bridges between Spheres," Chemical Engineering Science, 20(10) pp. 859-866. [115] Brackbill, J. U., Kothe, D. B., and Zemach, C., 1992, "A Continuum Method for Modeling Surface Tension," Journal of Computational Physics, 100(2) pp. 335-354. [116] FLUENT Inc., 2006, "Fluent 6.3 User's Guide," Lebanon, New Hampshire. [117] Issa, R. I., Ahmadi-Befrui, B., Beshay, K. R., and Gosman, A. D., 1991, "Solution of the Implicitly Discretised Reacting Flow Equations by Operator-Splitting," Journal of Computational Physics, 93(2) pp. 388-410. [118] Darhuber, A. A., Troian, S. M., and Wagner, S., 2001, "Physical Mechanisms Governing Pattern Fidelity in Microscale Offset Printing," Journal of Applied Physics, 90pp. 3602. [119] Brakke, K. A., 1992, "The Surface Evolver," Experimental Mathematics, 1(2) pp. 141-165. [120] Chadov, A. V., and Yakhnin, E. D., 1979, "Investigation of the Transfer of a Liquid from One Solid Surface to another. I. Slow Transfer, Method of Approximate Calculation," Kolloidnyi Zhurnal, 41, 4pp. 817-820. - 137 - [121] Yakhnin, E. D., and Chadov, A. V., 1983, "Investigation of the Transfer of a Liquid from One Solid Surface to another. 2. Dynamic Transfer," Kolloidn.Zh, 45pp. 1183. [122] Syamlal, M., Rogers, W., and O'Brien, T. J., 1993, "MFIX Documentation: Theory Guide," Technical Note, DOE/METC-95/1013, . [123] Garg, R., Galvin, J., Li, T., and Pannala, S., 2010, "Documentation of Open-Source MFIX–DEM Software for Gas-Solids Flows,". [124] Flemmer, C. L., 1991, "On the Regime Boundaries of Moisture in Granular Materials," Powder Technology, 66(2) pp. 191-194. [125] Nagahashi, Y., Lee, D. H., Grace, J. R., Epstein, N., Yokogawa, A., and Asako, Y., 2003, "Enhancement of Large-Particle Gas-Fluidization by Adding Liquid," AIChE Journal, 49(3) pp. 675-681. [126] Nagahashi, Y., Epstein, N., Grace, J. R., Asako, Y., and Yokogawa, A., 2006, "Spouting Enhancement by Addition of Small Quantities of Liquid to Gas-Spouted Beds," The Canadian Journal of Chemical Engineering, 84(5) pp. 527-531. 138 Appendix A. VISCOUS FORCE EXPRESSION BETWEEN TWO SPHERES According to the thin liquid film (or lubrication) approximation, whenever the local radii of curvature are large compared to the minimum distance between surfaces, and at low Reynolds numbers, the equation relating the pressure in the thin liquid film in the region between the two surfaces can be written as follows [19, 24]: 3 ( )( ) 12 mdhd dP rrH r r dr dr dt µ = ( A-1) where ( )H r is the separation distance between two surfaces, P is the pressure in the thin film, µ is the liquid viscosity, mh is the minimum separation distance between the two surfaces, occurring at 0r = . For the rest of this document, for convenience and consistency, the minimum separation distance is replaced with the separation distance (2D) i.e. 2mh D= . For simplicity, the liquid bridge is assumed to be cylindrical; therefore, a cylindrical system of coordinates is used for describing the problem geometry. In order to derive Eq. ( A-1), it is assumed that the liquid film is very thin; therefore, the pressure is assumed to be constant across the film thickness and only varies in one direction (here, in the radial direction). Integrating Eq. ( A-1) and knowing that pressure should be finite at 0r = leads to the following equation: 3 ( ) 12 ( ) dP r r dD dr H r dt µ= ( A-2) The separation distance can be approximated as follows: 2 ( ) 2 rH r D R = + ( A-3) 139 where R is the sphere radius. Differentiating the above equation leads to 2rdr RdH= . Assuming a finite wetted radius ( wr ) and therefore, a finite amount of liquid in the bridge, Eq. ( A-2) can be integrated, leading to the following pressure profile in the liquid film: 3 2 2 1 1( ) 12 3( ) ( ) ( ) wr wr r dD dDP r dr R H r dt H r H r dt µ µ = − = − ∫ ( A-4) Given the above pressure expression, the viscous force expression can be determined as follows: 2 2 0 3 2 12 ( ) 1 2 ( ) wr vis w D dDF rP r dr R H r D dt pi piµ = = − ∫ ( A-5) Considering a cylindrical liquid bridge, its volume can be determined as follows: 2 2 0 2 ( ) ( ) 4 2 wr b wV rH r dr R H r D pi pi = = − ∫ ( A-6) Rearranging the above equation, the relationship between the liquid height at wr , liquid bridge volume ( bV ), and minimum half-separation distance (D) can be determined: 1/ 2 22( ) 4bw VH r D Rpi = + ( A-7) Substituting the above in Eq. ( A-5) yields to the following viscous force expression: 2 23 1 2vis v dDF R X D dt piµ= ( A-8) where 2 2 ( ) 1v DX D D a = − + , 2 2 bVa Rpi = , and bV is the liquid bridge volume. The above equation represents the viscous force expression in its scalar form. The viscous force in its vector form can be written as follows: 2 23 2vis v uF R X D piµ= − r r ( A-9) 140 where ur represents the particle velocity ( u dD u e dt = r r ). The above viscous force expression correctly determines the direction of the viscous force with respect to the particle motion. Figure A-1 shows the variations of the correction term (Xv) with respect to the half- separation distance for different liquid bridge volumes. According to this figure, the correction factor increases as the amount of liquid increases and the separation distance decreases. It can also be concluded that the correction factor significantly reduces the viscous contribution in the case of small liquid volumes and large separation distances. Note that the correction factor approaches unity as the liquid volume approaches infinity. Figure A-1: Variations of the correction factor (Xv) with respect to the separation distance for four different dimensionless liquid bridge volumes. 141 Appendix B. CAPILLARY FORCE EXPRESSION BETWEEN TWO SPHERES Neglecting the gravity and buoyancy effects, the resulting capillary force between two surfaces connected with a pendular liquid bridge can be the summation of two effects: i) the surface tension acting on the liquid-solid boundary, and ii) the reduced hydrostatic pressure in the bridge, and can be written as follows [12]: 2 2 22capF r r Ppi σ pi= + ∆ ( B-1) Writing the reduced pressure with the aid of the Laplace-Young equation, the capillary force can be expressed as follows [13]: 2 2 2 2 1 2 1 1 12 1cap rF r r r r r pi σ pi σ piγσ = + − = + ( B-2) where 1r and 2r are the radii of the meridian profile and at the neck, respectively (see Figure B-1). The above is known as the “gorge method.” Alternatively, the adhesive force can be evaluated at the boundary between the particle and liquid bridge. This method is known as the “boundary method” [14]. Pitois et al. [15] approximated the bridge with a cylindrical profile and proposed the following simplified capillary force expression (in its scalar form): 2 coscap vF R Xpi σ θ≈ ( B-3) where R is the particle radius, σ is the surface tension, θ is the contact angle (assumed to be constant), 2 2 ( ) 1v DX D D a = − + , 2 2 bVa Rpi = , and bV is the liquid bridge volume. 142 r θ r1 r2 R R 2D z Figure B-1: Liquid bridge between two particles. Pitois et al. [15, 26] compared their capillary force expression, as in Eq. ( B-3), against a more complicated capillary force expression and experimental data, and found an excellent agreement. Eq. ( B-3) represents the capillary force in its scalar form. Knowing that the capillary forces always act towards attracting the surfaces, its vector form can be written as: 2 coscap v nF R X epi σ θ= r r ( B-4) where σ is the liquid surface tension, θ is the solid-liquid contact angle, vX is the correction factor (as given earlier), and ne r is the unit vector normal to the particle surface along the axis of collision. Writing the capillary force as above, regardless of the direction of movement of the particle, the attraction caused by the capillary effect will be properly determined. For instance, if one considers the left-hand side particle in Figure B-1, the direction of the capillary force vector will be predicted in the positive direction during all the collision stages, meaning that the capillary force is attracting the particles together. Similar verification can be carried out for the right-hand side particle during the approach stage and for the rebound stage of both particles. 143 Appendix C. DERIVATION OF THE OVERALL COEFFICIENT OF RESTITUTION The overall coefficient of restitution of two objects colliding axially (in 1-D) can be defined as follows: 2, 1, 1,0 2,0 ' f fu u e u u − = − ( C-1) where 1, fu and 2, fu are the scalar final (after impact) velocity of the first and second objects, respectively, and 1,0u and 2,0u are the scalar initial (before impact) velocity of the first and second objects, respectively. When using the above equation, it is required to choose one direction as “positive” and the opposite as “negative”; then, use the scalar value of the velocity with the proper sign accordingly. The coefficient of restitution has values between 0 to 1, where 1 represents a perfectly elastic collision, and 0 represents a perfectly inelastic collision. The law of the conservation of linear momentum for the above collision can be written as follows: 1 1,0 2 2,0 1 1, 2 2,f fm u m u m u m u+ = + ( C-2) where 1m and 2m are the mass of the first and second objects, respectively. Solving the system of equations consisting of Eq.s ( C-1) and ( C-2), the after collision velocities can be determined, as follows: 2 2,0 1,0 1 2 1, 1 2 ( ' 1) ( ' ) f e m u u m e m u m m + + − = + ( C-3) 1 1,0 2,0 2 1 2, 1 2 ( ' 1) ( ' ) f e m u u m e m u m m + + − = + ( C-4) 144 The law of the conservation of energy for the system consisting of the above objects can be written as: 2 2 2 2 1 1,0 2 2,0 1 1, 2 2, 1 1 1 1 2 2 2 2f f t m u m u m u m u L+ = + + ( C-5) where tL represents all the losses. Substituting the after collision velocity expressions in Eq. ( C-5) and solving the resulting equation for the overall coefficient of restitution, the following expression can be found: 1 2 2 1 2 0,1 0,2 2 ( ) ' 1 ( ) tL m me m m u u + = − − ( C-6) According to the above, the objects will stick together whenever the losses exceed a critical value that depends on the mass and initial velocity of two objects, i.e. ' 0e = for 2 1 2 0,1 0,2 1 2 ( ) 2( )t m m u u L m m − ≥ + . Solving the previous problem for two identical particles approaching with the same initial velocities, the overall coefficient of restitution can be simplified to: 2 0 ' 1 tLe mu = − ( C-7) 145 Appendix D. DERIVATION OF THE VELOCITY EXPRESSIONS Neglecting the capillary effects and considering only the viscous effects, the general equation of motion in the vector form for the approach stage of particles can be written as follows: 2 23 2 v u duR X m D dt piµ− = r r ( D-1) where the velocity vector on the left-hand side of the above equation can be written as u dD u e dt = r r ; on the other hand, the derivative term on the right-hand side can be written as u d udu e dt dt = rr r . Substituting these two expressions in Eq. ( D-1), the above equation can be written in the following scalar form: 2 23 1 2 v dD duR X m D dt dt piµ = ( D-2) Integrating the above leads to the following velocity expression: 0 ( )ln( ) v u f D u St C = ( D-3) where vSt is the viscous Stokes number ( 02 2 3v muSt Rpiµ = ), 0u is the initial approach velocity ( ) 2 2 2 2 2 ( ) D D af D D D a + = + + , 2 2 bVa Rpi = , bV is the liquid bridge volume, and C is the integration constant that is determined with the aid of the boundary condition. Let u= u0 at x= h0 where h0 is the initial height of the liquid layer; this yields to the following scalar form of the velocity expression for the approach stage: 146 0 0 0 0 ( ) ( )11 ln , ln( ) ( ) ( )0, ln ( ) v v a a v a f h f h u St St f D f h u f hSt f h − > = ≤ ( D-4) Writing the approach velocity in the above form, the wrong negative velocity values for the bottom particle during this stage will be automatically avoided. The contact velocity (right before contact) using the above expression can be written as follows: 0 0 0 0 ( ) ( )11 ln , ln( ) ( ) ( )0, ln ( ) v v a a c v a f h f h u St St f h f h u f hSt f h − > = ≤ ( D-5) where ha represents the average height of sphere surface asperities. Similarly, using Eq. ( D-2), the velocity expression for the rebound stage can be written: 0 ( )ln( ) 'v u f D u St C = ( D-6) Knowing that cu eu= − (where e is the dry coefficient of restitution) at aD h= , the integration constant ( 'C ) can be determined and it will be equal to: 0 ' ( ) exp c va eu StC f h u = ( D-7) Eventually, the velocity profile for the rebound stage can be written as follows: 0 0 0 0 0 0 ( )( )ln , ln( ) ( ) ( )0, ln ( ) c v v v a c a r v c a u eu St u f hf D St St f h u eu f h u u f hSt eu f h − > = ≤ ( D-8) 147 Each particle in the rebound stage moves in the opposite direction of its approach stage velocity, so if for any reason the sign of the rebound velocity, at an arbitrary half- separation distance, becomes the same as the approach velocity, it will mean that at some point the rebound stage velocity has changed its sign. Also, the kinetic energy of particles for the rebound stage may not be adequate to overcome the resisting viscous effects during the rebound stage; therefore, particles may stop before reaching a separation distance equal to the initial height of the liquid layer ( 0D h= ). The above velocity expression avoids the wrong positive velocity value for the rebound stage. Note that in order to derive the above velocity expressions, the capillary effects are ignored. The above velocity expressions, i.e. Eq.s ( D-4), ( D-5) and ( D-8), are used for deriving the analytical expressions for the viscous work and contact loss in Chapter 2. - 148 - Appendix E. DERIVATION OF THE VISCOUS STOKES MODEL As mentioned in Chapter 2 of the thesis, the collision of particles can be separated into three stages: approach, contact, and rebound. In the approach stage, the minimum half-separation distance (D) decreases from D = h0 to D = ha while in the rebound stage this distance increases from D = ha to a final separation distance D = h1. Ennis et al. [21] assumed that the critical condition occurs as the rebound velocity becomes zero at the original half-separation distance i.e. u = 0 at D = h0. In order to derive the Stokes criterion as proposed by Ennis et al. [21], it is necessary to find the velocity expression for all the above stages. Similar to Appendix D, by the neglecting capillary effects, the equation of motion for the particle in its scalar form can be written as follows: 23 1 2 dD duF ma R m D dt dt piµ= ⇒ = ( E-1) The general solution to the above equation can be written as follows: 0 ln v u D u St C = ( E-2) where 202 / 3vSt mu Rpiµ= is the viscous Stokes number, 0u is the initial approach velocity, and C is the integration constant that is determined according to the boundary condition of each stage. The viscous Stokes number is a measure of the particle kinetic energy to the viscous dissipation. Since 34 / 3 pm Rpiρ= , the Stokes number can also be written in the following form: 08 9 p v u R St ρ µ = ( E-3) where pρ is the particle density, 0u is the initial particle velocity, R is the particle radius, and µ is the liquid viscosity. In order to determine the integration constant for the approach stage, let u = u0 at D = h0 where h0 is the initial height of the liquid layer; this yields to the following velocity expression for the approach stage as follows: - 149 - 0 0 0 11 ln , lnv v a h h u u St St D h = − > ( E-4) Therefore, the contact velocity (right before contact) can be written as follows: 0 0 0 11 ln , lnc v v a a h h u u St St h h = − > ( E-5) where ha represents the average height of sphere surface asperities. With the aid of similar arguments as in Appendix D, Eq. ( E-2) can be written and solved for the approach stage. Knowing that u = -euc at D = ha (where e is the dry coefficient of restitution), the integration constant ( 'C ) can be determined, as follows: 0 ' exp c va eu StC h u = ( E-6) Therefore, the velocity profile for the rebound stage can be written as follows: 0 0 0 0 ln , lnc v v v a c a u eu St u hD u St St h u eu h = − > ( E-7) According to Ennis et al. [21], the critical Stokes number is obtained whenever the rebound velocity becomes zero at D = h0 i.e. *0 0 0 * * 1ln 1 ln 0v v a v a u h h u eSt St h St h = − − = ( E-8) The above equation yields to the following equation: *0 0ln ln 0v a a h h eSt e h h − + = ( E-9) Rearranging the above equation, the critical viscous Stokes expression can be written as follows: * 011 lnv a hSt e h = + ( E-10) - 150 - Whenever the viscous Stokes number exceeds the critical value, it indicates that the particle will rebound and coalescence occurs whenever the reverse is true. Eq. ( E-3) and ( E-10) constitute the main body of the viscous Stokes model as proposed by Ennis et al. [21] for the binary collision of two wet identical solid spherical particles. By equating the viscous Stokes number to the critical value, the critical velocity can be determined as follows: * * 0 9 8 v p St u R µ ρ = ( E-11) The alternative and more general approach is to assume that particles stop at a distance other than the original liquid height. Assuming that particles stop at D = h1, Eq. ( E-9) can be written as follows: * 01ln ln 0v a a hh eSt e h h − + = ( E-12) Therefore, the general critical Stokes number can be written in the following form: * 0 11ln lnv a a h hSt h e h = + ( E-13) where h1 is the final stopping distance of the particle after rebound.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Mathematical modeling of interaction of wet particles...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Mathematical modeling of interaction of wet particles and application to fluidized beds Darabi, Pirooz 2011
pdf
Page Metadata
Item Metadata
Title | Mathematical modeling of interaction of wet particles and application to fluidized beds |
Creator |
Darabi, Pirooz |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | In many industrial operations, such as fluidized bed granulators, coaters, and fluid cokers, a binding or reacting liquid is introduced into the system. Due to the effects of liquid, the multi-phase transport phenomena of these systems are more complicated compared to conventional gas-solid fluidization systems. In this thesis, mathematical modeling is used to study the interaction of wet particles. First, a coalescence model is developed to describe the binary collision of wet particles. The model is in the form of a wet coefficient of restitution and is used to determine the critical velocity—the boundary between coalescence or rebound outcomes—for a range of capillary numbers. Model predictions are compared with the available experimental data and good agreement is found. The model accounts for both liquid viscosity and surface tension effects and is used to investigate the boundary between collisions with dominant capillary and respective viscous effects. Then, by incorporating time- and temperature-dependant variations of the viscosity and thickness of the liquid coating, the model is used to determine the agglomeration tendency of bitumen-coated coke particles in fluid cokers. A simplified mathematical model and numerical solution of the Navier-Stokes equations are used to study the rupture of stretching liquid bridges between two solid spherical particles. The simplified model considers the geometry of the problem in which the gas-liquid interface is represented with a parabola. The numerical simulations of the Navier Stokes equations are performed with FLUENT and are used to investigate the viscous, surface tension, inertial, gravitational, and contact angle effects on the rupture distance and liquid distribution. Finally, the interaction of multiple wet particles is addressed by implementing the wet coefficient of restitution proposed in this thesis, using MFIX, an open-source Discrete Element Method (DEM) tool. DEM simulations of a fluidized bed consisting of mono-sized solid spherical particles pre-coated with identical liquid coatings are performed, and the effect of coating viscosity and thickness on the fluidization behaviour is investigated. Snapshots of the instantaneous particle positions are presented, and time-averaged values of the bed centroid in the y-direction, wet coefficient of restitution, and relative normal collision velocity are analyzed. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-05-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0080685 |
URI | http://hdl.handle.net/2429/34297 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_fall_darabi_pirooz.pdf [ 8.26MB ]
- Metadata
- JSON: 24-1.0080685.json
- JSON-LD: 24-1.0080685-ld.json
- RDF/XML (Pretty): 24-1.0080685-rdf.xml
- RDF/JSON: 24-1.0080685-rdf.json
- Turtle: 24-1.0080685-turtle.txt
- N-Triples: 24-1.0080685-rdf-ntriples.txt
- Original Record: 24-1.0080685-source.json
- Full Text
- 24-1.0080685-fulltext.txt
- Citation
- 24-1.0080685.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0080685/manifest