UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 FLOW AND SEGREGATION BEHAVIOR OF FLUIDIZED SYSTEMS By AKHIL ARUN RAO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHI LOSOPHY UNIVERSITY OF FLORIDA 2011 PAGE 2 2 2011 Akhil Arun Rao PAGE 3 3 To my family and friends for their encouragement and support throughout graduate school PAGE 4 4 ACKNOWLEDGMENTS My sincere gratitude goes to my advisor, Dr. Jennifer Curtis for her constant guidance, supervision and encouragement. Her perspective has helped me gain a panoramic view on research. I also appreciate the intellectual advice provided by Dr. Carl Wassgren and Dr. Bruno Hancock throughout the course of my PhD. I would also like to thank Dr. Mark Pepple and Deepak Rangarajan without whom getting some of the experimental data would have been difficult. Deepak Rangarajan and I have worked together to obtained the data for the 1.5 m m glass particles (Chapter 3). I am also grateful for all my friends, roommates and office colleagues, and the pedantic lunch table discussion with them that kept me energized. Finally, I am indebted to my family and friends back home (in India). I shall cherish their love and support forever. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................. 4 LIST OF TABLES ............................................................................................................ 8 LIST OF FIGURES ........................................................................................................ 10 LIST OF ABBREVIATIONS ........................................................................................... 15 ABSTRACT ................................................................................................................... 21 CHAPTER 1 INTRODUCTION .................................................................................................... 23 Dilute Turbulent FluidSolid Fluidized Flows ........................................................... 23 Gas Solid Flows ............................................................................................... 23 Liquid Solid Flows ............................................................................................ 26 Fluidized Beds ........................................................................................................ 27 Binary Segregation ........................................................................................... 27 Wall Effects ...................................................................................................... 29 2 FLUIDIZED DILUTE TURBULENT GASSOLID FLOW ......................................... 30 Background ............................................................................................................. 30 Fluctuating Inter action Terms ........................................................................... 30 Drag Relation ................................................................................................... 32 Solid Phase Stress Description ........................................................................ 33 Mathematical Model ................................................................................................ 33 Fluid Stress (rz ) ............................................................................................... 34 Drag Force ( FD) ................................................................................................ 36 Granular Energy Dissipation ( ) ...................................................................... 38 Solid Phase Stress ( rz, rr, ) ....................................................................... 38 Granular Energy Flux ( qpTr) .............................................................................. 39 Fluctuating Interaction/Coupling Terms ( Ik, IT) .................................................. 41 Fluctuating Energy Transfer (FET) Model for the Fluctuating Interaction/Coupling Terms ........................................................................... 47 Boundary Conditions .............................................................................................. 50 Numerical Solution .................................................................................................. 51 Benchmark Data Sets ............................................................................................. 53 Results and Discussion ........................................................................................... 54 Comparing the Various Models for Interaction Terms and Cross Correlations ................................................................................................... 55 Comparing the Various Models for Drag Terms ............................................... 62 PAGE 6 6 Comparing the Various Models for Solid Stress Terms .................................... 62 Summary ................................................................................................................ 63 3 FLUIDIZED DILUTE TURBULENT LIQUIDSO LID FLOW ..................................... 88 Background ............................................................................................................. 88 Experimental Setup ................................................................................................ 90 Flow Loop ......................................................................................................... 90 PDPA System ................................................................................................... 92 Single Phase Validation ................................................................................... 92 Particles and Laser Setti ngs ............................................................................. 93 Two Phase Data Validation .............................................................................. 94 Eulerian TwoFluid Models for Dilute Turbulent LiquidSolid Flows ........................ 95 High ST, InertiaDominated Flow Regime ( ST > 40) ........................................ 95 Low ST, Viscous Dominated Flow Regime ( ST < 5) ......................................... 97 Intermediate ST, Transitional Flow Regime (5 < ST < 40) ................................ 98 Results and Discussion ......................................................................................... 100 High ST, InertiaDominated Flow Regime ...................................................... 101 Low ST, Viscous Dominated Flow Regime .................................................... 102 Intermediate ST, Transitional Flow Regime .................................................... 102 Summary .............................................................................................................. 105 4 BINARY SEGREGATION IN FLUIDIZED BED ..................................................... 120 Background ........................................................................................................... 120 Experiments .......................................................................................................... 124 Particles and Their Preparation ...................................................................... 124 Experimental Setup ........................................................................................ 125 Experimental Procedure ................................................................................. 126 Segregation Index ................................................................................................. 127 Results and Discussion ......................................................................................... 128 Mixture Types ................................................................................................. 129 Type A mixtures ....................................................................................... 130 Type B mixtures ....................................................................................... 132 Type C mixtures ....................................................................................... 132 Type D mixtures ....................................................................................... 134 Type E mixtures ....................................................................................... 134 Type F mixtures ....................................................................................... 134 Type G mixtures ....................................................................................... 135 Classification Diagram .................................................................................... 135 Experiments Validating Minimal Wall Effects ........................................................ 138 Effect of column H/D on Binary Segregation in Fluidized Bed .............................. 139 Summary .............................................................................................................. 140 PAGE 7 7 5 EFFECT OF COLUMN DIAMETER AND BED HEIGHT ON MINIMUM FLUIDIZATION VELOCITY ................................................................................... 156 Background ........................................................................................................... 156 Experimental Setup .............................................................................................. 157 Theory ................................................................................................................... 160 Resul ts and Discussion ......................................................................................... 167 Summary .............................................................................................................. 169 6 CONCLUSIONS AND RECOMMENDATIONS ..................................................... 175 Dilute Turbulent FluidSolid Fluidized Flows ......................................................... 175 Fluidized Beds Segregation of Binary Mixtures and Wall Effects ....................... 178 APPENDIX A DILUTETURBULENT FLUIDSOLID FLOW PROGRAM .................................... 181 B RAW FLUIDSOLID DATA .................................................................................... 227 C FIGURES FOR FLUIDIZED BED BINARY MIXTURES ........................................ 233 LIST OF REFERENCES ............................................................................................. 248 BIOGRAPHICAL SKETCH .......................................................................................... 254 PAGE 8 8 LIST OF TABLES Table page 2 1 Closure relations employed in the published gas solid flow models ................... 65 2 2 Summary of experimental data ........................................................................... 71 2 3 Time scale used in FET model for the different experimental data ..................... 72 2 4 Percentage Error in Predicted Gas Velocity Fluctuations Using the Various Combin ations of Different Interaction Terms and Velocity Cross Correlations ... 85 2 5 Percentage Error in Predicted Solid Velocity Fluctuations Using the Various Combinations of Different Interaction Terms and Velocity Cross Correlations ... 85 2 6 Percentage error in predicted gas velocity fluctuations using the various drag relations .............................................................................................................. 85 2 7 Percentage error in predicted gas velocity fluctuations using the various solid stress closures .................................................................................................... 87 3 1 PDPA settings .................................................................................................. 110 3 2 Summary of experimental data ......................................................................... 111 3 3 Summary of model equations in the high, low and intermediate ST flow regimes ............................................................................................................. 112 4 1 Summary of published work ............................................................................. 141 4 2 Experimental material and properties present study ...................................... 142 4 3 Experimental mixtures and thei r properties present study ............................. 142 4 4 Details of mixtures and their properties published work and present study ... 144 5 1 Experim ental materials and properties ............................................................. 170 B 1 Raw data for d = 1.5 mm, m = 0.0175 and Vfzcl = 3.07 m/s ............................... 227 B 2 Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 1 .................... 227 B 3 Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 2 .................... 228 B 4 Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 3 .................... 228 B 5 Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 4 .................... 229 PAGE 9 9 B 6 Raw data for d = 1.5 mm, m = 0.075 and Vfzcl = 3.04 m/s ................................. 229 B 7 Raw data for d = 1.5 mm, m = 0.0175 and Vfzcl = 5.03 m/s ............................... 230 B 8 Raw data for d = 1.5 mm, m = 0.0425 and Vfzcl = 4.90 m/s ............................... 230 B 9 Raw data for d = 1.5 mm, m = 0.075 and Vfzcl = 5.04 m/s ................................. 231 B 10 Raw data for d = 1.5 mm, m = 0.0175 and Vfzcl = 7 m/s .................................... 231 B 11 Raw data for d = 1.5 mm, m = 0.0425 and Vfzcl = 7.06 m/s ............................... 232 B 12 Raw data for d = 1.5 mm, m = 0.075 and Vfzcl = 7.02 m/s ................................. 232 PAGE 10 10 LIST OF FIGURES Figure page 2 1 Suspended particles enclosed in a box .............................................................. 65 2 2 Profile for the fluctuating energy transfer ............................................................ 66 2 3 Program flow chart ............................................................................................. 67 2 4 Comparison of solutions using various number of grids ..................................... 68 2 5 Comparison between the Bolio et al. code and the code used in the present study ................................................................................................................... 69 2 6 Comparison of prediction in the fluctuation velocity between (a) Benavides and van Wachem and (b) present study ............................................................. 69 2 7 Comparison between the MFIX code and t he present study .............................. 70 2 8 Present model predictions compared to Tsuji et al. ................. 73 2 9 Present model predictions compared to Tsuji et al. ................. 74 2 10 Present model predictions compared to Tsuji et al. ............... 74 2 11 Present model predictions compared to Tsuji et al. ............... 75 2 12 Pr esent model predictions compared to Jones et al. ................. 76 2 13 Present model predictions compared to Sheen et al. ......................................... 77 2 14 Louge et al. crosscorrelation predictions compared to Tsuji et al. ..................... 78 2 15 Koch and Sangani cross correlation predictions compared to Tsuji et al. .......... 78 2 16 Wylie et al. crosscorrelation predictions compared to Tsuji et al. ...................... 79 2 17 Simonin cross correlation predictions compared to Tsuji et al. ........................... 79 2 18 Sinclair and Mallo cross correlation predictions compared to Tsuji et al. ........... 80 2 19 FET model predictions compared to Tsuji 243 articles ............................... 80 2 20 Louge et al. cross .. 8 1 2 21 Koch and Sangani cross correlation p particles .............................................................................................................. 81 PAGE 11 11 2 22 Wylie et al. cross.... 82 2 23 Simonin cross ........ 82 2 24 Sinclair and Mallo cross particles .............................................................................................................. 83 2 25 Comparing the magnitudes of ksg (a), Ik (b) and IT(c) for the various m = 3.2 ............... 84 2 26 Comparing the magnitudes of FD (a) and (b) for the various drag models m = 3.2 ........................................................................ 86 2 27 Comparing the magnitudes of s/ f (a) and / f (b) for the various solid stress ........................................................... 86 3 1 Effect of ST on the solid path in the fluid .......................................................... 107 3 2 Flow loop .......................................................................................................... 108 3 3 Comparing singlephase fluctuating velocity to the singlephase kturbulence model ( fz vk ) .............................................................................. 109 3 4 Solid fluctuating velocity measurement for d = 1.5 mm, m = 0.0425 and Vfzcl ~ 3 m/s (averaged over four runs) ....................................................................... 111 3 5 Mean velocity ( Vfz and Vsz) measurement for d = 2.32 mm compared to flow predictions for high ST case ............................................................................. 113 3 6 k, T measurements for d = 2.32 mm compared to flow prediction for high ST case .................................................................................................................. 113 3 7 Mean velocity ( Vfz and Vsz) measurements for d = 0.5 mm compared to flow predictions for low ST case .............................................................................. 114 3 8 Fluctuating velocity ( vfz and vsz ) measurements for d = 0.5 mm compared to flow predictions for low ST case ....................................................................... 115 3 9 Mean velocity ( Vfz and Vsz) measurements for d = 1 mm compared to flow predictions for intermediate ST case ................................................................ 116 3 10 Fluctuating velocity ( vfz and vsz ) measurement for d = 1 mm compared to flow predictions for intermediate ST case ......................................................... 117 3 11 Mean velocity ( Vfz and Vsz) measurement for d = 1.5 mm compared to flow predictions for intermediate ST case ................................................................ 118 PAGE 12 12 3 12 Fluctuating velocity ( vfz and vsz ) measurement for d = 1.5 mm compared to flow predictions for intermediate ST case ......................................................... 119 4 1 Schematic of the experiment. ........................................................................... 143 4 2 Typical pressure drop profiles and segregation index behavior for a Type A mixture 50G38550G083. ................................................................................. 146 4 3 Process of fluidization for a homogeneous binary mixture which fluidizes at two points ......................................................................................................... 147 4 4 Typical pressure drop profiles and segregation index behavior for a Type B mixture 75S32825G328. ................................................................................. 148 4 5 Typical pressure drop profiles and segregation index behavior for a Type C mixture 50G19550G083. ................................................................................. 149 4 6 Pressure drop profiles for the mixture 75G23125G116 from an initiall y mixed state .................................................................................................................. 150 4 7 Typical pressure drop profiles and segregation index behavior for a Type D mixture 50G16550G083. ................................................................................. 151 4 8 Mixture type diagram. ....................................................................................... 152 4 9 Correlation for the minimum fluidization velocity ratio and the Archimedes number ratio. .................................................................................................... 153 4 10 Segregation profiles for the mixture 75G23125G116 ...................................... 153 4 11 Segregation profiles of 25G11675G231at V* =2.3 ........................................... 154 4 12 Segregation profiles of 50G13850G328 at V* =2.5 .......................................... 154 4 13 Segregation profiles of 50G13850G275 at V* =3 ............................................. 155 4 14 Segregation profiles of 70G11630P275 at V* =3.1 ........................................... 155 5 1 Example of a pressure drop profile (fluidization and defluidization) using G231 particles in the 1.6 cm diameter column ................................................. 170 5 2 Repmf as a function of H/D, for dif ferent glass particles in the column D =1.6 cm ..................................................................................................................... 171 5 3 Repmf as a function of H/D, for different glass particles in the column D =2.4 cm. .................................................................................................................... 171 5 4 Repmf as a function of H/D, for polystyrene particles in the column D =1.6 cm. 172 PAGE 13 13 5 5 Repmf as a function of H/D, for polystyrene particles in the column D =2.4 cm. 172 5 6 Sketch showing the fluidized bed coordinate system. ...................................... 173 5 7 A schematic showing how the vertical stress in the bed varies with bed depth at different fluid velocity V1, V2 and V3. ............................................................. 173 5 8 The minimum fluidization velocity as a function of the column diameter using the experimental data from Liu et al. ................................................................ 174 5 9 Comparison of the curves from equation 5 30 to the experimental data and existing correlations. ......................................................................................... 174 C 1 Pressure drop profiles and segregation index behavior for a Type A mixture 50G55050G083. .............................................................................................. 233 C 2 Pressure drop profiles and segregation index behavior for a Type A mixture 50G46250G083. .............................................................................................. 234 C 3 Pressure drop profiles and segregation index behavior for a Type A mixture 50G55050G116. .............................................................................................. 235 C 4 Pressure drop profiles and segregation index behavior for a Type A mixture 50G38550G083. .............................................................................................. 236 C 5 Pressure drop profiles and segregation index behavior for a Type A mixture 50G46250G116. .............................................................................................. 237 C 6 Pressure drop profiles and segregation index behavior for a Type B mixture 50G32850G083. .............................................................................................. 238 C 7 Pressure drop profiles and segregation index behavior for a Type B mixture 50G27550G083. .............................................................................................. 239 C 8 Pressure drop profiles and segregation index behavior for a Type C mixture 50G23150G083. .............................................................................................. 240 C 9 Pressure drop profiles and segregation index behavior for a Type C mixture 50G19550G083. .............................................................................................. 241 C 10 Pressure drop profiles and segregation index behavior for a Type D mixture 50G16550G083. .............................................................................................. 242 C 11 Pressure drop profiles and segregation index behavior for a Type D mixture 50G13850G083. .............................................................................................. 243 C 12 Pressure drop profiles and segregation index behavior for a Type D mixture 50G11650G083. .............................................................................................. 244 PAGE 14 14 C 13 Pressure drop profiles and segregati on index behavior for a Type B mixture 13S32887P328. ............................................................................................... 245 C 14 Pressure drop profiles and segregation index behavior for a Type B mixture 75S32825G328. .............................................................................................. 246 C 15 Pressure drop profiles and segregation index behavior for a Type E mixture 70G11630P275. .............................................................................................. 247 PAGE 15 15 LIST OF ABBREVIATIONS aE Ergun equation coefficient [kg m4] Aht area of heat transfer [m2] APL, BPL, CPL, DPL coefficients based on e for Peirano and Leckner 51 Aso, Bso parameter based on for Syamlal and OBrien model 49 bE Ergun equation coefficient [kg m3 s1] c1 coefficient of vertical stress equation [m1] c2 coefficient of vertical stress equation [kg m2 s2] CD friction factor Cmix composition of mixture cT1, cT2, cT3, c constants in k Cw coefficient for turbulence generation by vortex shedding c parameter based on angle between flow direction and slip velocities d particle diameter [m] dr ratio of particle diameter of jetsam to flotsam e particle particle coefficient of restitution ew particle wall coefficient of restitution Ew turbulence generation by vortex shedding [kg m1 s3] F F0, F1, F2, F3 dimensionless functions based on Rep and for Hill et al. 47, 48 model FB stress for ce on the bottom face [kg m s2] Fbu buoyancy force [kg m s2] FD drag force per unit volume [kg m2 s2] fD fluctuating drag force [kg m2 s2] Fd drag force [kg m s2] FET Fluctuation Energy Transfer PAGE 16 16 FF wall friction force [kg m s2] FT stress force on t he top face [kg m s2] fT1, fT2, fT3 f coefficients in kmodel FW weight of the element [kg m s2] g gravitational acceleration [m s2] G1k, G1c, dimensionless functions based on e and T G2k, G2c, dimensionless functions based on e and T G3k, G3c dim ensionless functions based on e and T go radial distribution coefficient H bed height [m] h instantaneous fixed bed height [m] ht heat transfer coefficient [kg s2 K1] Ik interaction term in the k equation [kg m1 s3] IT interaction term in the T equati on [kg m1 s3] k fluid turbulence kinetic energy [m2 s2] k1, k2, k1, k2 constants for the wall affected Ergun's equation kfb coefficient based on kJ Jassen's proportionality constant ksf fluid solid fluctuating velocity cros scorrelation [m2 s2] LDV LASER Doppler Velocitimetry m ratio of solid mass flux to fluid mass flux M Mixing Index MFBs Micro Fluidized Beds p pressure [kg m1 s2] PDPA phase doppler particle analysis PAGE 17 17 qpT granular energy flux [kg s3] R radius of pipe [m ] r radial component [m] Re Reynolds number based on mean flow Rep particle Reynolds number Repmf particle Reynolds number at minimum fluidization ReT Reynolds number based on T RT turbulent R eynolds number SI segregation index ST Stokes number T granular temperature [m2 s2] t Time and Volume Based Averaging [s] Tf real temperature of a fluid [k] TVBA Time and Volume Based Averaging Uc complete fluidization velocity [m s1] Umf minimum fluidization velocity [m s1] Umf lower lower minimum fluidization vel ocity of the two components [m s1] Ur ratio of minimum fluidization velocity of jetsam to flotsam UTO takeover velocity [m s1] U friction velocity [m s1] V average fluid velocity [m s1] V* ratio of operating velocity to compl ete fluidization velocity VBA Volume Based Averaging Vf mean gas velocity [m s1] v'f fluctuation fluid velocity [m s1] PAGE 18 18 Vfzcl fluid center line velocity [m s1] Vrm dimensionless velocity based on Rep and for Syamlal and OBrie n model 49 Vs mean solid velocity [m s1] v's fluctuation solid velocity [m s1] wHLK dimensionless parameter based on Rep and for Hill et al. 47, 48 model y+ dimensionless distance from the wall z axial component [m] Z paramet er based on velocity Greek drag coefficient [kg m3 s1] FET coefficient [kg m3 s1] dissipation of fluid turbulence [m2 s3] fluidized bed viodage parameter based on e t ratio of Lagrangian time scale to particle relaxation time scale a zimuthal component angle between flow direction and slip velocities sphericity conductivity of granular temperature modified by inelastic collisions [kg m1 s1] conductivity of granular temperature modified by inelastic collisions [kg m1 s1] mfp mean free path between particles [m] b fluid phase bulk viscosity [kg m1 s1] ef fluid phase viscosity modified by particles [kg m1 s1] PAGE 19 19 f fluid phase viscosity [kg m1 s1] s solid phase viscosity modified for inelastic collisions [kg m1 s1] s' solid phase viscosity [kg m1 s1] T fluid phase eddy viscosity [kg m1 s1] t turbulent fluidphase viscosity for vortex shedding [kg m1 s1] solid volume fraction 0 maximum packing solid volume fraction r parameter based on relative velocity and k avg average density of fluidsolid mixture [kg m3] f fluid density [kg m3] r ratio of density of jetsam to flotsam s solid density [kg m3] solid stress tensor [kg m1 s2] c collisional solid stress tensor [kg m1 s2] H horizontal stress due to friction [kg m1 s2] k kinetic solid stress tensor [kg m1 s2] k, constants in kmodel v vertical stress [kg m1 s2] fluid stress tensor [kg m1 s2] c collisional time scale [s] D particle relaxation time scale [s] e eddy time scale [s] f fluid response time scale [s] L Lagrangian time scale [s] PAGE 20 20 p solid response time scale [s] sf time scale for FET model [s] TUR turbulent component of fluidphase stress [kg m1 s2] dissipation rate of granular energy [kg m1 s3] specularity friction angle between particles and wall ratio of mass of solid in a unit volume to the mass of fluid i n the same volume function based on '' f unctions based on Rep, ReT dampening function PAGE 21 21 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy FLOW AND SEGREGATION BEHAVIOR OF FLUIDIZED SY STEMS By Akhil Arun Rao May 2011 Chair: Jennifer S. Curtis Major: Chemical Engineering A model for dilute turbulent gas solid flow with inelastic collisions in an Eulerian framework necessitates many different closures. The various closure models propos ed to date in the literature that are required for the fluctuating interaction terms, the solidfluid fluctuation velocity cross correlation, the solid stress and the drag formulation are compared. A new interpretation is provided for the fluctuating inter action terms, based on a fluctuating energy transfer mechanism. All the closures are evaluated against benchmark experiments for fully developed, dilute and turbulent gas solid flow in a vertical pipe. The Eulerian model is also extended to the case of dil ute turbulent liquidsolid flows in the viscous, transitional and inertial flow regimes. The model is validated with both published experimental data and with new experimental data obtained in this study. These experiments were performed in a three inch vertical pipe, with 1.5 mm diameter glass particles and Re from 2.0x105 to 5.0x105. The concentration of solids is varied from 0.7% to 3% by volume. Experimental investigations of binary mixtures in a gas fluidized bed have also been carried out in this study. Glass and polystyre are fluidized in a micro fluidized bed of diameter 1.6 cm. These investigations exhibit PAGE 22 22 varied pressure drop profiles and segregation patterns. D ifferent mixture types are mapped on a graph of density ratio versus size rat io. Categorizing binary mixtures in this manner gives a qualitative understanding of how the different mixtures behave upon fluidization. Further, fluidized bed experiments of particles with monosized distributions revealed that the minimum fluidization v elocity of particles increases as the diameter of the fluidization column is reduced, or if the height of the bed is increased. It is found that wall friction leads to an enhancement in the minimum fluidization velocity. A new, semi correlated model is proposed which incorporates Janssens wall effects in the calculation of the minimum fluidization velocity. The model predictions compare favorably to existing correlations and experimental data. Equation Chapter 1 Sectio n 1 PAGE 23 23 CHAPTER 1 INTRODUCTION Fluidized systems are encountered in an enormous number of industrial applications These applications vary from dilute and low mass loading conditions like pneumatic transport of powders, pulverized coal injection into entrained flow gasifiers, and cyclones, to dense fluidized bed applications like slurry transport, circulating fluidized beds, drying of pharmaceutical pow ders and fluidsolid catalytic reacti ons 1. Although fluidization technology has been around for decades, there are several aspects of fluidization about which not much is known. Dilute Turbulent FluidSolid Fluidized Flows GasSolid Flows Modeling of gas solid fluidized flows is essential for designing and optimizing, and for understanding the energy distribution of the various internal processes in equipment involving such flows. Unfortunately, even for the case of dilute fluidsolid flow there is a lack of fundamental understanding of the effect of the solid phase on the flow of fluid phase and vice versa. The experimental data considered in the present work for fully developed dilute turbulent gas solid flow in a vertical pipe 2, 3, 4, 5 is in the inert ia dominated flow regime, as most pneumatic transportation applications exhibit this flow regime. The f low regime of fluidsolid flows is det e rmined by the Stokes number ( ST ) of the particles. Louge et al. 6 develop ed an Eulerian two fluid model for dilute turbulent gas solid flow with particle particle interactions for the inertia dominated flow regime. They were the first to use a single k equation model to describe gas phase turbulence. This model was further developed by Bolio et al. 7 who employed a tw o equation kmodel to PAGE 24 24 describe the gas phase turbulence. Simonin 8 and Benavides and van Wachem 9 have suggested models along similar lines. All these models are able to describe the mean velocity profiles of the gas and solid phases well. The granular t emperature behavior for particles with low particle Reynolds numbers ( Rep ~10, Rep is based on particle diameter and slip velocity) is also captured fairly well by these models, with some deviations from the experimental data of Tsuji ** for 243 polysty rene beads. According to these data, as the mass loading (ratio of the solid mass flux to the fluid mass flux, m ) increases, the granular temperature decreases Furthermore, for smaller mass loadings ( m ~1) the granular temperature is virtually independent of the radial position, while for larger loadings ( m ~3 5) the granular temperature increases towards the wall. This behavior suggests that for larger loadings (but still dilute phase flow), particles tend to concentrate at the center of the pipe. These tre nds are correctly predicted by the models. Unfortunately, granular temperature data for large and intermediate Rep are not available in the literature. In contrast all of these models under predict the gas phase turbulence in the presence of the particles In fact, these models also predict that large particles (d>400 ) exhibit turbulence dampening, contrary to experimental observations as sum marized by Gore and Crowe 10. Further, Hestroni 11 showed that particles with small Rep (~10) exhibit turbulence dampening, while particles with large Rep (~1000) cause turbulenc e enhancement of the gas phase due to vortex shedding. Particles with intermediate Rep (~100) show inbetween behavior, exhibiting turbulence enhancement at the core of the pipe and turbulence dampening at the wall. The present study further ** Tsuji Y. Private Communication 1993. PAGE 25 25 advances the w ork of Bolio et al. 7 and corrects for the deficiencies in these other gas solid flow models. To generate reliable predictions for gas solid flow in an Eulerian framework, it is critically important to have accurate closure models. In particular, the gas p article turbulent interactions, the effect of the presence of particles on gas phase turbulence and vice versa, must be appropriately described. For stateof theart Eulerian, dilute gas solid flow models with particleparticle interactions, there are several important forces/interactions requiring constitutive modeling. The flow models vary from one another in the following key aspects The fluctuating interaction/coupling terms in the gas phase turbulent kinetic energy and granular temperature equations w hich includes the gas solid fluctuating velocity crosscorrelation The drag relation The solidphase stress description Although the closure models for the above terms have evolved over time, there is no single combination of closures for an Eulerian model for dilute turbulent gas solid flow with particle particle interactions that is generally accepted. In the present study, different models for the various closures are compared against the benchmark data sets of Tsuji et al. 2, Jones et al. 3, Sheen et al 4, and Lee and Durst 5. Further, a novel but simple interpretation for the fluctuating interaction term is suggested which leads to improved flow predictions and match the experimental data favorably as compared to the other models. Also, a new combination of various closures models is proposed such that there is minimum error between the predictions and the experimental data. PAGE 26 26 LiquidSolid Flows For the case of liquidsolid fluidized flows all of the three flow regimes namely viscousdominated flow r egime, inertiadominated flow regime and the transitional flow regime have industrial applications. However, very little experimental data are available and even few er modeling studies have been carried out. Reliable and detailed measurements of the of l iquid solid flow behavior are difficult as intrusive techniques disrupt flow behavior, and, due to the lack of data, model predictions cannot be validated. Alajbegovic et al. 12 used laser Doppler anemometer and a single beam ray densitometer nonintrusive techniques, to study a system of water laden with ceramic particles of size 2.32 mm Pepple 13 used the p haseDoppler particle analyzer (PDPA) to measure the mean and fluctuating velocities of glass particles of size 0.5 and 1 mm. The experiments performed herein, extend the work of Pepple 13 to obtain velocity data for 1.5 mm glass particles. The data for the 0.5 mm glass 13 show viscous dominated flow b ehavior, the 1 mm 13 and the 1.5 mm (present s tudy) glass exhibit transition flow regime behavior, and 2.32 mm ceramic particles 12 provide data for the inertiadominated flow regime. Modeling attempts of liquidsolid fluidized flows are rare. Hadinoto and Curtis 14 extended the work of Bolio et al. 7 to include fl uid lubric ation effects and employed their gas solid model to generate flow predictions corresponding to the operating conditions in the liquid solid data of Alajbegovic et al. 12. However, the Hadinoto and Curtis 14 model is prescribed for only particles with high inertia and is neither appropriate to describe the experimental data of Pepple 13 nor the experimental data obtained in the present study. The new combination of closure models for gas solid flows, proposed in PAGE 27 27 this work, is extended to describe the inertia dominated flow data of Alajbegovic et al. 12. The closures for solid stress and fluctuating interaction term are changed to that proposed by Chen and Wood 15 to describe the viscous dominated flow, and a bridge model between the above two models is used t o describe the intermediate flow regime. On increasing the loading of the solid phase (volume fraction of solids>10%) dense fluid so lid flow is observed, and the correlations of the fluctuating solid volume fraction and fluctuating velocity govern the flow behavior. Dense fluidsolid flows have not been studied in this work. Dense f luid solid flows in vertical pipes transition to fluidized beds at lower fluidization velocities. Fluidized Beds Binary Segregation Fluidized bed operations are gaining a lot of interest as they involve intimate contact between a fluid stream and a dense mixture of various solid materials leading to good heat and mass transfer coefficients An e normous amount of research on fluidized beds with monosized particles has been carried out and the flow patterns of such systems are well studied 1. In the industry however, monosized particles are rarely used and mixtures of particle with various sizes and densities are generally used. A common, observable phenomenon associated with multi solid beds is that they may segregate when subjected to fluidization. Segregation is primarily due to particle size and/or density differences. The segregation tendency of a powder mixture influences the overall process efficienc y. Much experimental data exist in the literature concerning fluidized bed segregation, yet the current understanding of the mechanisms controlling multi solid fluidization segregation is very poor, even for two component particle mixtures. In fact, PAGE 28 28 for a simple binary mixture fluidized by a gas, there are a variety of pressure drop profiles that have been observed and the different pressure drop profiles lead to different segregation patterns. For example, Formisani et al. 16 presented several binary mixtures having large density and/or size ratios that fluidize at two different velocity points, while Joseph et al. 17 has reported mixtures which fluidize at a single velocity point, just like a monocomponent powder. Marzocchella et al. 18 studied an extremely disparate mixture that m ixes well at low velocities, just beyond the point of fluidization, but segregates at larger velocities. Complicating matters further, the phenomenon of layer inversion has also been observed in fluidized mixtures containing smaller, denser particles and c oars er, less dense particles 19. In t he present study an attempt is made to categorize the various gas fluidized binary mixture types reported by various authors. In addition, a correlation between the minimum fluidization velocity ratio, m m fjetsam r fflotsamU U U and the density ratio, sjetsam r sflotsam and size ratio, jetsam r flotsamd d d is proposed that can be used to distinguish between these mixture types, the observed pressure drop profiles, and segregation patterns. Cases in which th e particle size dr and density ratio r (the causes for Ur) aid each other, r and dr are sufficient to give an idea of the level of disparity and, hence, make conclusions about the mixtures pressure drop or segregation behavior. But, in some cases where the particle size and density ratio oppose each other, it is difficult to gauge the level of disparity (nature of its pressure drop and segregation behavior) based only on r and dr. Thus, using Ur as a measure of disparity will be beneficial in such cases PAGE 29 29 Wall Effects As stated earlier, fluidized bed operations have become more popular in industry and, hence there i s a drive to study these beds in greater detail at laboratory scales. Reducing the size of fluidized beds is gaining interest because of th e better gas distribution and operability at smaller scales. Micro fluidized Beds (MFB), which were first proposed by Potic et al. 20, have even better mixing and heat transfer capabilities, making them a useful tool for studying various processes like dry ing, mixing or segregation 21, 22, 23, 24, as well as reaction kinetics 25. For example, MFBs are used to investigate fluidization segregation in the pharmaceutical industry where active ingredients are expensive and difficult to obtain during the early st ages of development. One parameter of particular interest when working with fluidized beds is the minimum fluidization velocity. Experimental data obtained in the present study indicate that the minimum fluidization velocity increases as the column diameter is decreased or if the fixed bed height is increased. Many correlations exist for calculating the minimum fluidization velocity 26 39 but n one of the correlations include the effect of the bed height or diameter of the column (wall effects) factors wh ich are likely to be of interest when MFBs are considered. Hence, a correlation is developed which encompasses the effect of aspect ratio and the ratio of column diameter to particle diameter to predict the elevation in the minimum fluidization velocity. Equation Chapter 2 Section 1 PAGE 30 30 CHAPTER 2 FLUIDIZED DILUTE TURBULENT GASSOLID FLOW Background Eulerian two fluid model for dilute turbulent gas solid flows, as discussed in Chapter 1, requires closure models for the following relations The fluctuating interaction term and the fluidsolid fluctuating velocity cross correlation within it The drag The solid stress In the literature many different closures have been suggested for each of these terms. Fluctuating Interaction Terms Louge et al. 6 suggested an approximation for closing the fluctuating Interaction term. The fluctuating interaction term they proposed was b ased on the gas turbulence, granular temperature, and the gas solid velocity crosscorrelation. Time and volume based averaging (TVBA) was used for developing expressions for the fluctuating interactions. Louge et al. 6 closed the gas solid velocity crosscorrelation, modifying Koch 40, in which it was assumed that a dilute gas solid suspension at very small particle Reynolds number in the limit where solid body collisions determine the particle velocity distribution function. TVBA and the modified Koch 40 gas solid velocity cross correlation were also used by Bolio et al. 7 and then tested by Benavides and van Wachem 9. Along similar lines, Igci et al. 41 used a crosscorrelation expression developed by Koch and Sangani 42, but Igci et al. 41 neglected gas phase turbulence in their analysis. More recently Hadinoto and Curtis 14 used TVBA for the interaction term and applied a closure model for the gas solid velocity crosscorrelation developed by Wylie et al. 43 for particles with high inertia and moderate fluid inertia. Both the Koch and PAGE 31 31 Sangani 42 and Wylie et al. 43 closure models are extensions of the original Koch 40 model, wherein the most important assumption is that the inertia of the fluid is negligible and the solid particle interactions occur for part icles in viscous fluids with small particle Reynolds number. Bolio et al. 6, Benavides and van Wachem 9 and Hadinoto and Curtis 14 observed that the abovementioned closure models could not predict the phenomenon of turbulence enhancement and under predicted the gas phase turbulence in the presence of particles A closure model for the gas solid velocity crosscorrelation based on time scales was developed by Simonin 8 which is discussed later (in the Mathematical Model section). He too used TVBA to dev elop closure models to describe these fluctuating interaction terms. The NETLMFIX code (Multiphase Flows with Interphase Exchanges) employs Simonin 8 closure model to simulate dilute, turbulent gas solid phase flow. Benavides and van Wachem 9 have also tested the Simonin 8 closure. Predictions from the MFIX code, as well as the work of Benavides and van Wachem 9, agai n showed that the Simonin 8 model cannot predict the phenomenon of turbulence enhancement and under predicts the gas phase turbulence in the presence of particles Benavides and van Wachem 9 and Hadinoto and Curtis 14 also test a simple closure model for the gas solid velocity cross correlation devel oped by Sinclair and Mallo 44. This model assumes the fluidsolid fluctuation velocity cross cor relation to be a geometric mean of the gas phase turbulence and the granular temperature. Finally, Crowe 45 developed a different form for the fluctuating interaction term. The Crowe 45 form includes additional generation based on the square of the relati ve velocity between the two phases. Volume based averaging (VBA) was used to develop PAGE 32 32 new relations for the fluctuating interaction terms in contrast to most of the other models which have used TVBA to develop expressions for the fluctuating interaction ter m s. Only Zhang and Reese 46 have used the Crowe 45 formulation, along with a modified Koch and Sangani 42 model for the fluidsolid velocity crosscorrelation, and observed a good match with available experimental data. A new closure model for the fluctuat ing interactions is introduced in this study based on a fluctuating energy transfer (FET) mechanism by analogy to a heat transfer mechanism. In this new closure, the Sinclair and Mallo 44 model i s used for the gas solid fluctuating velocity cross correlati on. Gas solid flow predictions based on this closure compare more favorably to experimental data for gas turbulence and granular temperature than all of the previously published models. This new model is capable of predicting the gas phase turbulence behav ior associated with a wide range of particle sizes such as turbulence dampening in the presence of particles, turbulence enhancement as well as in between turbulence behavior. Drag Relation Wen and Yu 39 is probably the most widely used and accepted partic le drag model. However, Hadinoto and Curtis 14 showed that the choice of the drag model may affect the predicted mean and fluctuating velocity profiles in dilute gas solid flow for small and low density particles at low velocities. In the present study, fl ow predictions employing the (1) Wen and Yu 39 drag model, based on experimental correlation, the (2) Hill et al. 47, 48 drag model, based on theory and Lattice Boltzman simulations, and the (3) Syamlal and OBrien drag model 49 also based on experimental correlation, are investigated and compared. PAGE 33 33 Solid Phase Stress Description Two different closures for the solidphase stress are considered in this study. Both of these originate from the kinetic theory for granular flow. The Lun et al. 50 granular flow theory was the first fundamental description for the solidphase stress based on analogy with molecular collisions in a dense gas. The second model, based on the work of Peirano and Leckner 51, includes the effect of gas phase turbulence and the effect of the interstitial fluid in the description for the solidphase stress. Table 21 summarizes the closures used by the various dilute, turbulent gas solid flow models in the literature. As outlined above, several model formulations obtained by combining the different interaction/coupling term formulations along with the various gas solid velocity cross correlations, three different drag force models, and two different solid phase stress models, are compared and tested against benchmark experimental data provid ed by Tsuji et al. 2, Jones et al. 3, Sheen et al. 4, and Lee and Durst 5. These data sets were obtained under the conditions of fully developed, dilute, turbulent gas solid flow in a vertical pipe using laser Doppler velocimetry. The mean and fluctuating velocity profiles of the gas and solid phases have been reported in these experiments. After evaluating all the combinations of the models, Table 21 puts forth the final model which yielded minimum error when compared to the various data sets Mathematical Model The continuum equations for an Eulerian description of dilute, turbulent fluidsolid flow have been well documented in many studies. The two fluid equations used in this work follow those given in Bolio et al. 7. The simplified equations are also consistent with the ones employed in MFIX and Benavides and van Wachem 9. Gas Phase Momentum Balance (z component) PAGE 34 34 1 0rzDz fdp rFg rr dz ( 2 1 ) Solid Phase Momentum Balance (z component) 1 0 ()rzDzsfrFg rr ( 2 2 ) Solid Phase Momentum Balance (r component) 1 0rrr rrr ( 2 3 ) Granular Energy Balance, 1 0sz PTrrz TV rq I rr r ( 2 4 ) Transport Equations for k, 21 01 11ef T Tfz k fkf fV k rI rr r r ( 2 5 ) 2 11 2 22 331 01 1 11. ef T Tfz TT ff f TT TTkV r cf rr r kr cfcfI kk ( 2 6 ) Equation 2 1 to equation 2 6 presents the Bolio et al. 7 governing equations for the simplified case of fully developed, dilute, turbulent gas solid flow in a vertical pipe of radius ( R) Fluid Stress (rz ) The total fluid stress ( rz) has two components. First, there is the viscous fluid component in which the intrinsic fluid viscosity ( f) is affected by the presence of the solids. The effect on f for very dilu te solids concentration cases i s given by Batchelor and Green 52. PAGE 35 35 fz rzefTV r ( 2 7 ) 2 012.57.61eff. ( 2 8 ) The second component in the total fluid stress is the t urbulent contribution (Reynold Stress) which is simplified using an eddy viscosity closure. For the twoequation kturbulence model, the eddy viscosity is described based on the turbulent kinetic energy ( k) and its dissipation rate ( ) 2 f Tcfk ( 2 9 ) cT1 =1.4, cT2 =1.8, cT3 =1.2, c=0.09k =1.4 =1.3,11Tf 7, 2 2 22 1exp 1exp 965T TR y f, ( 2 10 ) 3.45 1exp1 70Ty f R, ( 2 11 ) where, f fU yRr ( 2 12 ) and 2 f T efk R. ( 2 13 ) The turbulence model coefficients ( cT1, cT2, cT3, c, k, fT1, fT2 and f) are the same as given in Bolio et al. 7 and are summarized in equation 2 9 to equation 2 13. PAGE 36 36 Drag Force ( FD) The drag force is propor tional to the relative velocity and to the drag coefficient which is a function of the solids volume fraction ( ), and particle Reynolds number. Dz fzszFVV ( 2 14 ) Three drag models considered i n this study: Wen and Yu 39, Hill et al. 47, 48, and Syamlal and OBrien 49. 2.653 4 1g D fzszCVV d ( 2 15 ) 0.68724 10.15Re Re Dp pC for Rep < 1000, ( 2 16 ) 0.44DC for Rep > 1000, ( 2 17) 1 Reffzsz p fdVV ( 2 18 ) Equation 2 15 to equation 2 18 provide d etails for the Wen and Yu 39 model. 2181fF d ( 2 19 ) 3 1Re 8 pF for 0.01 and 2 31 Re 3/8pF F ( 2 20 ) 2 01RepFFF for 0.01 and 33102 14 Re 2pFFFFF F ( 2 21 ) 23RepFFF for 0.01 and 2 31 Re 3/8pF F ( 2 22) 23RepFFF for 0.01 and 33102 14 Re 2pFFFFF F ( 2 23) PAGE 37 37 0 3 2313/2135/64ln17.14 1 10 10.6818.488.16 1HLK HLKFw w ( 2 24) for 0.4 1 12 40 F for 0.1 ( 2 25) 2 3 2313/2135/64ln17.891 10 10.68111.0315.41 1HLK HLKFw w ( 2 26) for 0.4 30.95310.03667 F for 0.0953 ( 2 27) exp100.4/HKLw ( 2 28) 1 Re 2ffzsz p fdVV ( 2 29) Equation 2 19 to equation 2 29 provide the details for the Hill et al. 47, 48 model. 2 23 0.634.8/Re 4f rmpfzsz rmVVV Vd ( 2 30) 2 20.50.06Re0.06Re0.12Re2rm p p pSOSOVA BAA ( 2 31) 4.141SOA ( 2 32) and 2.651SOB, ( 2 33) Reffzsz p fdVV ( 2 34) Equation 2 302 34 provides the details for Syamlal and OBrien 49 model. PAGE 38 38 Granular Energy Dissipation ( ) Granular energy is dissipated due to inelastic collisions of particles with a restitution coefficient ( e ). This granular energy dissipation is given by Lun et al. 50 and is a function of the radial distribution function ( g0) the solid volume fraction ( ), and granular temperature ( T) 23/2 048 1sgT d ( 2 35) where 1 '' 3 sisiTvv ( 2 36) 1 2 e ( 2 37) 1/3 0 0 1/31/3 0g ( 2 38) Solid Phase Stress ( rz, rr, ) The solid phase stresses have two contributions: the kinetic (or translational) part ( k) and the collisional part ( c) Bolio et al. 7 modified the Lun et al. 50 stress expressions so that the kinetic contribution remains bounded in finitesized domains and d ampens according to a function proportional to the mean free path of the particles. Consequently, kc ( 2 39) where, 1 1/mfpR ( 2 40) and the mean free path, mfp is PAGE 39 39 62 mfpd ( 2 41) The descriptions for the normal solidphase stresses ( rr, ) are very similar in Lun et al. 50 and Peirano and Leckner 51, with the primar y difference being with the presence of the damping function in the Lun et al. 50 formulation. In the present study, the Bolio et al. 7 expressions have been followed for the normal solid stress when the Lun et al. 50 model is used. 11() rr skcGGT (for 50) ( 2 42) 11() rr skcGGT (for 51) ( 2 43) where, 1 kG, ( 2 44) 2 104cGg ( 2 45) The shear stress for the solid phase is expressed as a product of the solidphase viscosity ( s), and the solid velocity strain. Shear stress sz rzsV r ( 2 46) Granular Ener gy Flux ( qpTr) Similar to the solid phase shear stress, the granular energy flux is expressed as a product of the granular (or pseudothermal) conductivity ( ), and the gradient of granular temperature. pTr T q r ( 2 47) The expressions for the solids viscosity ( s) and conductivity ( ) given by Lun et al. 50 and Peirano and Leckner 51 are quite different. PAGE 40 40 22 sskcGG ( 2 48) '5 96ssdT ( 2 49) 20 018 132 25k Gg g ( 2 50) 2 0 20768 88 132 525 25cg Gg. ( 2 51) Equation 2 48 to equation 2 51 detail the Lun et al. 50 expressions for s. 2 0256 5 s bg. ( 2 52) Equation 2 52 gives the Lun et al. 50 expressions for b. 33 kcGG ( 2 53) 25 128sdT ( 2 54) 2 30 0812 143 41335k Gg g, ( 2 55) 2 3 009612 16 1 43 4133 541335 5cGgg. ( 2 56) Equation 2 53 to equation 2 56 detail the Lun et al. 50 expressions for 22 sskcGG ( 2 57) 1 2022 1 3PL k sgt PL Dc B GkTAg ( 2 58) 2028 5 ck T GgGd. ( 2 59) PAGE 41 41 Equation 2 57 to equation 2 59 detai l the Peirano and Leckner 51 expressions for s. 2 5 3 bscG ( 2 60) Equation 2 60 give s the Peirano and Leckner 51 expression for b. 33 skcGG ( 2 61) 1 30 93 9/5 1 102PL k sgt PL Dc D GkTCg ( 2 62) 30318 5 59ckT GgGd ( 2 63) where, 22/5131, 13/5, (1)21/100, 14933/100.PL PL PL PLAee Bee Cee Dee ( 2 64) Equation 2 61 to 2 64 detail the Peirano and Leckner 51 expressions for Here D is the drag time and c is the time between collisions, t is the ratio of Lagrangian time scale to the particle relaxation time scale and ks f is the fluid solid fluctuation velocity cross correlation (these variables are discussed in the next section). Fluctuating Interaction/Coupling Terms ( Ik, IT) Louge et al. 6 used TVBA and simplified the fluctuating interaction/coupling terms to the following form 12''k fisiI kvv, ( 2 65) 1''3T fisiIvvT, ( 2 66) PAGE 42 42 where Ik represents the interaction coupling between granular temperature and the gas phase turbulence equation and IT is the interaction coupling between the gas phase turbulence and the granular energy balance equation. For the sake of convenience, the crosscorrelation of the gas and solid fluctuating velocities '' fisivv will be denoted as ks f hereafter. The description for ks f is probably the most controversial point amongst the different closures required for gas solid flow modeling in an Eulerian framework. This controversy is due to the fact that there is a lack of fundamental understanding of the nature concerning the fluid solid interactions at the level of the velocity fluctuations. In the gas solid flow model of Louge et al. 6, the authors modified Koch 40 who rigorously derived an expression for ks f. 21 4fzsz sf sVV d k T. ( 2 67) However, the system under consideration by Koch 40 was assumed to have negligible fluid inertia with particles suspended in a viscous flui d. Consequently, the Koch 40 expression for ksf lacks the effect of gas phase turbulence on the cross correlation. Hence, the expressions for fluidsolid fluctuating interactions that have evolved from the Koch 40 model may not be appropriate in the case of turbulent gas solid flow. Igci et al. 41 employed the work of Koch and Sangani 42, an extension of the Koch 40 model, in a system where the fluid turbulence was again neglected. Koch and Sangani 42 model use a factor to summarize the effect of the solids concentration. PAGE 43 43 2 2 3 081 1ffzsz sf fVV k gdT ( 2 68) where, 2 2 23 13/2135/64ln17.14 13.55.910.6818.488.16. ( 2 69) The function replaces the solid volume fraction functionality obtained from using the drag term ( ) in the work of Louge et al. 6. Hadinoto and Curtis 14 used another form of the Koch 40 model based on the work of Wylie et al. 43. This model included the effect of large particle iner tia and moderate fluid inertia obtained by enhancing the factor with a functionality which depends on ReT, and Rep. 2 2 3 081 1'ffzsz sf fVV k gdT ( 2 70) where, 'Re''fbpK ( 2 71) 50.03360.1060.0116(1)fbK, ( 2 72) 2 24 2 24 2 2Re ReReReReRe 2 ''12 1exp ReRe ReRe2Re 2Red TTTTT pp pp p perf ( 2 73) Reffzsz p fdVV, ( 2 74) Ref T f dT ( 2 75) PAGE 44 44 Simon in 8 developed a model for the gas solid fluctuation velocity cross correlation based on the Lagrangian time scale and particle relaxation time scale. The Lagrangian time scale is the time for which the effect of a passing eddy can be felt on the rest of t he fluid and is defined as 21e L rc ( 2 76) where the eddy time scale ( e the turn over time of an eddy ) is 3 2 ek cf ( 2 77) and, 2 23 2fs rVV k ( 2 78) 21.81.35cos' c ( 2 79) The parameter is the angle between the mean fluid and solid velocity. The particle relaxation or the drag time ( D) is the time for which the effect of a passing particle can be felt on the neighboring fluid and is defined as 1s D. ( 2 80) Simonin 8 assumed that the ratio of the Lagrangian time scale (L) to particle relaxation time (D ) is equal to the ratio of the fluctuation energy of the cross correlation carried by the gas phase to the fluctuation energy carried by the entire mixture. That is, (1)2(1)3f sf L t D f savgsfk kTk ( 2 81) Equation 2 81 on rearrangement gives, PAGE 45 45 23 11t sf tk kT ( 2 82) where 1s f. ( 2 83) Crowe 45 argued that it was incorrect to take the averaged flow variables as being local flow variables and to treat the governing equations as if t hey represented singlephase flow with a local coupling term. Hence, Crowe 45 suggested a new approach for the development of the governing equations which is not based on temporal and volume averaging, but is based on volume averaging only (VBA). The equations he used for the fluctuation interaction terms Ik and IT, are different than those normally used by others, 21 13k fs sfIVV Tk ( 2 84) 13 T sfIkT. ( 2 85) This formul ation includes an additional generation term which is proportional to the square of the relative velocity between the two phases ( equation 2 84), arising due to the new averaging technique. This additional generation helps predict the enhanced fluid turbulence. The model of Zhang and Reese 46 employs equation 2 84 and equation 2 85 for the fluctuating interaction terms along with a modified Koch and Sangani 42 model to describe the cross correlation between the fluid and solid fluctuating velocity, 2 c sffs DkVV ( 2 86) where the collision time c is the time between collisions and is defined as 024cd gT ( 2 87) PAGE 46 46 Zhang and Reese 46 showed good predictions when compared with the Tsuji et al. 2 data. Sinclair and Mallo 44 developed a simple model for the cross correlation of fluidsolid velocity fluctuations using a simple geometric mean, 6sfkkT ( 2 8 8 ) Equation 2 88 is obtained by assuming that the correlation of the fluctuating gas velocity and the fluctuating drag force ( '.fiDivf ) divided by the gas fluctuation velocity is equal to the correlation of the fluctuating solid velocity and the fluctuating drag force ( '.siDi vf ) divided by the fluctuating solid velocity, i.e. 1/2 1/2'. '. '.''.'siDi fiDi fifi sisivf vf vvvv, ( 2 89) where fDi is the fluctuating drag force. The correlation between fluctuating gas velocity and the fluctuating drag force is generally approximated as Ik, the effect of solids on the gas turbule nce, and the correlation of fluctuating solid velocity and the fluctuating drag force is approximated as IT, the effect of gas turbulence on the granular temperature. Thus equation 2 89 can be interpreted as 23kTII kT ( 2 90) Substitutin g equation 2 65 and equation 2 66 in to equation 2 90 gives, 23 23 sf sfkkkT kT ( 2 91) Equation 2 91 on simplification gives equation 2 88. PAGE 47 47 Fluctuating Energy Transfer (FET) Model for the Fluctuating Interaction/Coupling Terms To develop the new fluctuating interaction model (FET model), we consider the hypothetical situation where: Particles are enclosed and suspended freely in a box (Figure 2 1) with no gravitational acceleration. The fluid has negligible viscosity; hence, there is no momentum diffusion and no turbulence dissipation. The collisions of the particles are elastic; consequently, there is no dissipat ion due to collisions. The diffusion of particle momentum is also negligible. There are no wall effects. Now, assume that at the start of an experiment the particles are stationary but the gas has an evenly distributed turbulent energy given by k. In this situation, the gas phase turbulence will initiate some velocity fluctuations in the particle phase that will continue to grow. In this hypothetical situation, the gas phase turbulent kinetic energy and granular energy equations will reduce to, 1 '12f k sfk I kk t ( 2 92) 3 '13 2 s T sf T IkT t ( 2 93) Eq uation 2 92 and equation 2 93 are analogous to the case of simple heat transfer between two fluids. The interaction terms are similar to convection heat transfer where ( 2k ks f) and ( ks f3T) are analogous to temperature differences ( Tf1Tf12) and ( Tf12Tf2), respectively ( Tf1 and Tf2 are the bulk temperatures of the two fluids and Tf12 is the temperature at the interface), and '1 is analogous to the heat transfer coefficient PAGE 48 48 multiplied by the surface area ( htAht). For the defined hypothetical situation, the energy lost by the gas phase is equal to the energy gained by the solid phase, 3 1 2 sf Tk tt ( 2 94) Figure 22 shows how the difference in the magnitude of the velocity fluctuations causes a transfer of energy from one phase to another, similar to the transfer of heat from one fluid phase to another. This picture of the energy transfer is consistent with the assumption of interpenetrating continua of the two phases and the use of the twofluid model. The energy exchange will continue until the fluctuati ons equalize, that is, ks f=2k=3T. Typical units of '1 are M1L3T1 and so '1f is a time scale ( sf) over which the transfer of fluctuation energy occurs. In the following sections it will be shown that for particles exhibiting turbulence dampening, the time scale for fluctuation energy transfer is D while for particles showing inbetween behavior and turbulence enhancement, the time scale for fluctuation ener gy transfer is c. This behavior occurs because, for particles exhibiting turbulence dampening the mechanism of the fluid being pulled along by the particles causes the transfer of fluctuation energy. On the other hand, particles showing inbetween behavior and turbulence enhancement, the change in direction of individual solid particles and the corresponding disturbance caused by the particle wake are responsible for the transfer of the fluctuation energy. Finally, to describe the vortex shedding phenomenon associated with gas phase turbulent enhancement, the fluctuating interaction/coupling term, Ik ( equation 2 92) is enhanced by the wake effect Ew as given by Lun 53. PAGE 49 49 212wt wCk E d ( 2 95) Reffzsz p fdVV ( 2 96) 0.017Ret pf for 150Re310p, ( 2 97) 21.20.0000Ret pf for 310Re610p ( 2 98) 0.029Ret pf for Re610p, ( 2 99) 10/3wC for 150Re310p ( 2 100) 24/3wC for Re310p ( 2 101) Equation 2 95 to equation 2 101 detail the vortex shedding term. The coefficient Cw is slightly modified from what was ori ginally prescribed in Lun 53 in order to produce a better fit to the experimental data. Thus, the following equations represent the complete fluctuating interaction/coupling terms for the new model 2f k sfw sfIkkE ( 2 102) 3f T sf sf IkT ( 2 103) where, sf is the time scale for the fluctuation energy transfer ( sfD, for particles which show turbulence damping and sfc, for particles which show inbetween behavior and turbulence enhancement). A couple of observations can be made about the cross correlation ks f: According to t he new model and from the heat transfer analogy, ksf must lie between 2k and 3T. PAGE 50 50 At the wall, as the instantaneous fluid velocity becomes zero due to the noslip condition, the cross correlation ksf must also become zero. Keeping these two conditions in m ind, the Sinclair and Mallo 44 model is used for the cross correlation ks f since it is the only closure that is consistent with both conditions. Boundary Conditions Appropriate boundary conditions are also of prime importance to generate reliable flow predictions. At the center of the pipe, the symmetry boundary condition is used for all the flow variables. For the solid phase, boundary conditions at the wall for shear stress ( rz) and the flux of granular energy ( qpT) follow Johnson and Jackson 55. For the fluid phase at the wall, two sets of conditions were employed in this study. The first set of conditions is based on a low Reynolds number kturbulence model as used by Bol io et al. 7 where turbulent transport equations for the fluid phase are integrated to the wall For low Reynolds number kturbulence model the mean and fluctuating fluid velocities are zero (noslip boundary condition). Dissipation of gas turbulence at the wall is given by, 2 2 fef kk I r ( 2 104) The second set of boundary conditions f or the fluid phase employ wall functions for the fluid phase and follow the MFIX code 54. These conditions assume some fluid velocity slip at the wall. Also the wall conditions for k and are modified to be consistent with this approximation. The wall functions for the fluidphase were used only when the code was being validated against the MFIX code and the results of Benavides and van Wachem 9. For all the figures and tables in the pres ent study, the low Reynolds PAGE 51 51 number kturbulence model is used ( i.e. no slip boundary condition along with equation 2 104). Numerical Solution To solve the governing equations numerically, one operating condition for the fluid phase and one for the solid phase are required. In the code developed for the present study, two options are available for the operati ng condition on the fluid phase: the centerline fluid velocity or the fluid flow rate. The choice of which to use is made based on the conditions specified in the experimental data. For the solid phase, solid mass loading m is used to specify the operating condition. The solid loading m is defined as 1ssz f fzVrdr m Vrdr. ( 2 105) The governing equations including the closure models, boundary and operating conditions are solved using a combination of marching and iterative schemes as described in the following paragraph. The continuum equations are discretized using the finite volume method approach described in Patankar and Spadling 56. Approxim ately one hundred discrete grid points were distributed in a nonuniform pattern throughout the domain cross section with higher grid resolution near the wall. Bolio et al. 7 showed that sixty grid points were sufficient to generate a gridindependent solution. Figure 23 depicts the flow chart for the code. First, a guess of single phase dimensionless fluid velocity Vf z, gas turbulence k, and gas phase dissipation at a specific Re (where 2 Reffzcl fVR and Vfzcl is the centerline fluid velocity ) are given as input to the code. The code then perturbs the Re towards the direction of the desired PAGE 52 52 operating velocity and generates the new singlephase profi les for the perturbed Re in an iterative scheme. The operating condition is solved along with the discretized version of the singlephase gas momentum equation. The operating condition provides an additional equation which allows the pressure drop ( p ) to be treated as an unknown parameter. Thus, a pressure drop guess is not needed. The perturbations in Re are made until the operating Re is achieved ( i.e. the operating velocity is achieved). At this point an assumption is made tha t the fluid profile at very low solids loading will match the single phase profiles. Further, the guess ed profiles for Vsz, and T are considered to be independent of r (the radial coordinate). The single phase profiles along w ith the guess for Vsz and T are taken as guesses for the very dilute case of twophase profiles. Like the fluid operating condition, the solid operating condition is solved along with the solid phase momentum balance (r component) equation for the solid volume fraction profiles. The simulated flow profiles are used as a guess for the very dilute case. The loading is increased in small discrete steps (profiles for each loading step were calculated in an iterative scheme) until the two phase flow operating conditions are obtained. Thus in this fashion, a stable numerical solution is obtained. The program is written in a modular fashion so that it is simple to toggle amongst the various closure models and boundary conditions (Appendix A details the Matlab program) Hence, the program developed herein has the ability to easily reproduce the predicted flow patterns (with slight differences arising due to the numerical technique used) given in Bolio et al. 7 ( Figure 2 5 ) Benavides and van Wachem 9 ( Figure 2 6 ), and MFIX (fol lowing the work of Simonin 8, Figure 2 7 ). Therefore, this code provides a useful tool to compare various dilute turbulent gas solid models against the same PAGE 53 53 benchmark experimental data sets on the same platform. Unfortunately, the code is unable to reproduce the r esults of Zhang and Reese 46 even after consultation with the authors. Figure 2 4 compares the solutions obtained from the model using Wen and Yu 39 darg relation, Peirano and Leckner 51 solid stress and FE T model (with Sinclair and Mallo 44 crosscorrelation) for a variety of number of grid s. Although it is determined that fifty grids are sufficient to get a grid independet solution, ninety nine grids are used for all simulation runs. Benchmark Data Sets As stated earlier, the data under consideration are for pneumatic transport of solids by gas (air) in a vertical pipe. These flows are dilute and turbulent, and the flow profiles are fully developed. The data sets of Tsuji et al. 2, Jones et al. 3, Lee and D urst 5, and Sheen et al. 4 span various sizes of polystyren e and glass bead particles (70 varying between 0 to 5. Since the experimental data by Tsuji et al. 2 is the most widely cited the data sets as a standard by which all the previously published models, as well as the present model, will be evaluated. However, all of the models have been evaluated against all the experimental data summarized in Table 22, and it is observed that the model utilizi ng the combination of Syamlal and O'Brien 49 drag force relati on, the Peirano and Leckner 51 solid stress model, and the FET closure compare more favorably to the data than the other models. The values for the coefficient of restitution of particleparticl e/particlewall collisions and the specularity PAGE 54 54 factor are the same as those given in Bolio et al. 7 and are approximately similar to the values employed in the other published gas solid flow models. Results and Discussion From the simulations, profiles of the gas and solid mean velocities Vf z and Vsz, along with the profiles of the gas turbulent kinetic energy k and the granular temperature T are obtained. Since experimental data report values for the axial solid velocity fluctuations vsz these values nee d to be extracted from the predicted granular temperature T To convert T to vsz isotropic solid velocity fluctuations are assumed. Measurements of axial and radial solid velocity fluctuations show that this assumption is a reasonable one for dil ute gas solid flow 57, 'sz vT ( 2 106) For the gas velocity fluctuations, Sheen et al 4 has shown that the gas turbulence is not isotropic in pipe flows C onsequently it is assumed that the radial and azimuthal fluctuations in velocity 'frv and 'fv respectively, are approximately half of the axial velocity fluctuations 4, '' 2 fz frfv vv, ( 2 107) which results in 'fz vk ( 2 108) Figures 2 8 to 2 11 compare the predicted profiles of the new model (whi ch in corporates the Wen and Yu 39 drag force, the Peirano and Leckner 51 solids stress, the new FET model for the interaction between the gas and particle velocity fluctuations a nd the Sinclair and Mallo 44 model for the solidgas velocity crosscorrelati on) to the PAGE 55 55 experimental data of Tsuji et al. 2. Figures 2 12 and 2 13 compare the predicted profiles from the same model to the Jones et al. 3 and Sheen et al. 4 data, respectively. In all of these predictions, when the Stokes number ( 2182fzcls fdV ST R ) is less than 100, the particle relaxation time scale D is used as the time scale for the fluctuation energy transfer. When ST is less than 100, the mechanism of the fluid being pulled along by the particles causes the transfer of the fluctuation energy. On the other hand, when ST is greater than 100, the particle collision time scale c is used because the change in direction of the individual solid particles and the corresponding disturbance in the particle wake are responsible for the transfer of the fluctuation energy. Finally, if Rep > 150 (large particles with large relative velocity), Ew is activated. When Rep is greater than 150, the turbulent boundary layer detaches from the particle surface and vortices are generated in the wake of the solid particles, enhancing gas phase turbulence. Table 2 3 summarizes the time scales used by the FET model for the various experimental data (detailed in Table 22). Overall, the simulation results for the mean gas and solid velocities, as well as the fluctuating gas and solid velocities, favorably match the experimental data with this proposed model. The predictions for gas solid flow shown in Figure 2 11 for 2780 m particles are not as good as the other cases. This discrepancy is likely due to an inadequate description for vortex shedding with such large particles. Comparing the Various Models for Interaction Terms and CrossCorrelations The different models for t he interaction terms are compared using the gas phase turbulence and granular temperature data for the PAGE 56 56 et al. 2. These experimental data are used as the standard for uniformly comparing model predictions since mos t (all but Zhang and Reese 46) of the other closure models do not predict gas phase turbulence enhancement in the presence of large particles (1.42 mm and 2.78 mm). In order to focus this comparison on the effect of the model predictions from the various descri ptions for interaction terms and the cross correlation, the same drag f orce relation 39 and the same solidphase stress 51 are used in generating the model predictions. Furthermore, only the gas fluctuation velocity ( vf z ) and the solid fluctuation velocit y ( vsz ) are compared, as variations in the predicted mean gas and solid velocity profiles ar e negligible. Figures 2 14, 2 15, 2 16 2 17 and 2 18 compare the predictions of vf z profile using Louge et al. 6, Koch and Sangani 42, Wylie et al. 43, Simonin 8, and Sinclair and Mallo 44 crosscorrelation (TVBA is used for interaction terms, Ik and IT) respectively, to the data of Tsuji et al. 2. Most models are able to capture the very dilute mass loading flow cases but under predict the cases with higher mass loading. Figures 2 19 compare the predictions of vsz profile using FET model (using Sinclair and Mallo, 44, crosscorrelati on). Further, Figures 220, 221 2 22, 2 23 and 2 24 compare the predictions of vsz profile using Louge et al. 6, Koch and Sanga ni 42, Wylie et al. 43, Simonin 8, and Sinclair and Mallo 44 crosscorrelation (TVBA is used for interaction terms, Ik and IT) res pectively, to the data of Tsuji **. In the case of vsz profile, most models again under predict the solid fluctuating velocity ** Tsuji Y. Private Communication 1993. PAGE 57 57 Figure 225 a compares the magnitude of the sfk (where is calcu lated from the Wen and Yu 39 drag correlation and ks f is calculated from the various cross correlations) for the various combinations of interaction terms and the cross correlation models for the Tsuji et al. 2 case of 243 particles with a solids loading of 3.2. Similarly, Ik ( equation 2 102) along with the generation 21fz TV r (calculated using FET model only) are compared in Figure 225b, and IT ( equation 2 103 ), along with the rate of granular energy dissipation ( equation 2 35 calculated again using FET model only) are compared in Figure 225c. From Figure 225b, it is observed that the generation is very small at the core of the pipe and sharply rises toward the wall. For the FET model, the Ik term dominates at the core of the pipe, and the remaining two terms, the generation and the dissipation (not shown in Figure 225b), play an important role near the wall of the pipe. Since Ik dominates at the core of the pipe, the gas phase turbulence k largely depends on the inter action term. While, at the wall the gas phase turbulence depends on the generation and the dissipation. Further, from Figure 225c, the FET model predicts that IT will be negative and all three terms, in the granular energy balance equation, granular temperature generation (not shown Figure 225c), dissipation, and interaction play equally important roles throughout the pipe. No one term dominates at any point. In cases in which the Louge et al. 6, Koch and Sangani 42, and Wylie et al. 43 cross correlation models (these models use TVBA for Ik and IT) are used, the error in predicating the gas phase turbulence profile (Figures 2 14, 2 15 and 2 16 respectively) PAGE 58 58 increases as the mass loading increases, especially at the core of the pipe. This result is because the magnitude sfk in equation 2 65 is underestimated, and consequently, Ik is small compared to the gas phase generation ( Figure 225b) even at the center of the pipe. Additionally, the IT values estimated form Louge et al. 6, Koch and Sangani 42, and Wylie et al. 43 models are more negative than what is predicted from the FE T model resulting in relatively poorer model predictions (Figures 220, 221 and 2 22). All the crosscorrelation models that contain a square of the relative velocity 2 fzszVV produce a profile for ksf that increase rapidly near the wall as the relative velocity between the gas and solid phases is very large However, sfk must be zero at the wall (as ks f = 0 no slip condition). Although the Simonin 8 crosscorrelation model yields values for Ik which dominate the kequation at the pipe center, the predictions of shape and magnitude of the fluctuation gas velocity (Figure 2 17 and 2 23) are not as good as those obtained from the FET model (which uses the Sinclair and Mallo 44 model for cross correlation), especially for higher mass loading. These predictions are because the Simonin 8 model predict s too much turbulence damping since the Ik values obtained from it are smalle r than those from the FET model ( Figure 225 b) On the other hand, the Sinclair and Mallo 8 closu re (using TVBA for Ik and IT) predicts the order of magnitude of k and T profiles well but fails to predict the shape of the k profile correctly (Figure 2 18 and 2 24). As mentioned earlier, the code was unable to reproduce the r esults of Zhang and Reese 46 even after consultation with the authors. Stable model solutions could not be obtained for the Tsuji et al. 2 PAGE 59 59 particles were poor. The additional generation term 21fzszVV in equation 2 84 produces values for Ik which are up to two orders of magnitude larger than the typical range of values for Ik and yields an unrealistically large increase in the gas phase turbulence resulting in unstable solutions in some cases Based on Figure 225b and 2 25c, all the models predict Ik to be positive, which is a net generation of gas turbulence at the core of the pipe, while IT is negative which i s a net dissipation of granular temperature. The interaction term IT in the granular temperature equation behaves as dissipation taking energy from the solid phase and transferring it to the fluid phase resulting in a generation of gas phase turbulence. This idea goes well with the concept of the fluctuation energy transfer mechanism. It is observed experimentally that for small particles ( d ST < 100) the solid velocity fluctuations are larger than the fluid velocity fluctuations at the core of the pipe. Furthermore, for such particles, if the solids loading increases, the magnitude of the solid velocity fluctuations reduces due to enhanced particle particle collisions and larger inelastic dissipation The reduction in the granular temperature T (solid velocity fluctuations) tends to further dampen the gas phase turbulence through the energy coupling terms ( Ik, IT). Hence, higher the mass loading, more is the gas phase turbulence dampening (Figures 2 8 b and 2 8 c). Flow predictions employing the FET interaction model are able to reproduce these experimental observations and trends very well. For the case of very large particl es ( d Rep > 150, Figure 2 10b and 2 11) with a large degree of slip ( VfzVsz), the vortex shedding Ew overpowers the effect of the interaction term Ik in the k equation which results in a high gas phase turbulence. PAGE 60 60 Als o, as the solids mass loading increases more vortex shedding takes place, this further enhances the gas phase turbulence. The predic ted fluid velocity fluctuations employing the FET model are greater than the solid velocity fluctuations ( 2k > 3T) resulting in a net dissipation of the fluid phase velocity fluctuations and a net generation for the solid phase. Thus, a counter effect is seen in these cases with larger particles in which increases in solid mass loading reduces T by increasing but this effect is offset by the increase in IT. Overall, the granular temperature T slowly increases as the solids mass loading increases Unfortunately, solid velocity fluctuation data are not available for these larger particles to validate these model predictions. The i ntermediate particle sizes (300 d ST > 100 and Rep < 150, Figure 2 9 b) display flow behavior between the two extremes, where the gas turbulence is enhanced at the core of the pipe but dampened at the wall as compared to singlephase flow. Also, the simulations incorporat ing the FET interaction term model predict that the solid velocity fluctuations are greater than the fluid fluctuations at low solid mass loadings but the fluid fluctuations become greater than those of the solid as the mass loading increases Again, ther e are no solid fluctuation data to validate these model predictions for the solidphase fluctuations, but the model favorably matches the experimental gas phase fluctuation data. All of these predicted qualitative trends obtained by using the FET interacti on term model, along wi th the Sinclair and Mallo 44 crosscorrelation, are in line with the well known wor k of Gore and Crowe 10. These authors suggest that for small particles, drag is responsible for the turbulence transfer between the phases and turbulence damping. In the FET model, the drag time scale D is used as the time scale for energy transfer PAGE 61 61 for small particl es. Furthermore, Hestroni 11 showed that particles with intermediate Rep display inbetween behavior (i.e., gas phase turbulence enhancement in the core of the pipe and dampening at the wall). In the FET model, this behavior occurs when c is the time scale but vortex shedding is not activated (i.e., ST > 100 but Rep < 150). Hestr oni 11 also showed that for very large Rep, vortex shedding is responsible for enhanced turbulence as is seen in the cases in which Rep > 150. The vortex shedding, in the present study is described using Ew. The FET model predictions for particle sizes ranging from 70 wide range of mass loadings compare very well with the experimental data. Hence, the FET model along wi th the Sinclair and Mallo 44 crosscorrelation is recommended over the other previously proposed closure models for the interaction ter ms and the solidgas fluctuating velocity cross correlation. Gas phase turbulence and granular temperature predictions using different fluctuation Interaction terms are compared in Table 2 4 and Table 2 5 respectively. These tables provide the percentage error values of the predictions from the experimental data. The percentage error for the gas phase is defined as 2'' % *100 'exptsim fzjfzj j expt fzj jvv error v ( 2 109) where 'expt fzjv is the jth experimental value of the fluctuating gas velocity and 'sim fzjv is the predicted value of the fluctuating gas velocity at the same radial pos ition r at its corresponding jth experimental value. Percentage error for the solid fluctuating velocity predictions is similarly defined for the solid phase. PAGE 62 62 Comparing the Various Models for Drag Terms Flow predictions using different drag force relations are compared in Table 2 6 which provides the percentage error values of the predictions from the experimental data. To compare the different drag relations, the FET interaction model, along wi th the Sinclair and Mallo 44 crosscorrelation a nd the Perino a nd Leckner 51 solid stress model were employed. Figure 226 shows that the magnitude of the drag forc e based on the Wen and Yu 39 and Syamlal and OBrien 49 m = 3.2 of Tsuji et al. 2 are very similar. The Hill et al. 47, 48 drag relation predicts a drag force magnitude slightly smaller than the other two models. T he percentage error values in Table 2 6 for the Hill et al. 47, 48 drag relation are slightly larger than the other two models, suggesting that the Hill et al. 47, 48 drag relation may under predict the drag coefficient and thereby under predict the drag force. Hadinoto and Curti s 14 studied the effect of the various drag relations for particles with low inertia on two phase flow predictions. They observed some influence in the predictions of the slip velocity d ue to changes in the drag model; however, no such effec t is seen in the present study, since the particles studied herein have large inertia. H ence, the effect of the choice of the drag model is not that significant. Since the predicted flow profiles found by changing the drag models ar e similar, the Wen and Yu 39 model, a r obust and widely used closure, is recommended for the case of dilute, turbulent, gas solid flow models. Comparing the Various Models for Solid Stress Terms Figure 22 7 compares the magnitude of the solid phase viscosity s and solid phase conductivity using the Lun et al. 50 and the Peirano and Leckner 51 solid stress m = 3.2 in the experiments of Tsuji et al. 2. PAGE 63 63 The FET interaction model and the Wen and Yu 39 drag relation were employed in this comparison. Both the solid viscosity s and conductivity values, as predicted by the Lun et al. 50 model, are insensitive to the radial position r Both the s and are primarily functions of the granular temperature, which is essentially constant over the pipe cross section. In contrast the predicted solid viscosity and conductivity from the Peirano and Leckner 51 solid stress closure exhibit some dependency on the radial position r as this solid stress model incorporates a direc t effect of the fluid on the s and values. Table 2 7 shows that the two solid stress closures produce similar error in the model predictions from the experimental data. Furthermore, the detailed shape of the mean and fluctuating solid and gas profiles obtained from the two solid st ress closures are hardly different. Since the Peirano and Leckner 51 solid stress closure directly incorporates fluid effects, it is chosen over the Lun et al. 50 solid stress closure. It is anticipated that this direct incorporation of fluid effects may y ield improved flow predictions in situations involving liquidsolid systems. Summary Experimental data considered in the present Chapter is for high ST gas solid flows in the inertia dominated regime. The flow behavior in this regime, for a wide range of p article sizes from 70 m to 3 mm, fluidized at Re ~ 3x104, is predicted using an Eulerian model. The Eulerian two fluid model for dilute turbulent gas solid flows with particle particle interaction (must be considered in inertia dominated flow regime) requires many closures. The fluid phase stresses which are enhance by gas phase turbulence are closed using the eddy viscosity assumption which in turn is computed with the help of a standard two equation kmodel. While granular kinetic theory is used PAGE 64 64 to com pute the solidphase stresses. Two different granular kinetic theory models 50, 51 are compared in the present study The drag force behaves as the interaction between the fluid and solid phases at the mean lev el. To compute this force innumerable models a re suggested in the literature. However, o nly three models ( Wen and Yu 39, Hill et al. 47, 48 and Syamlal and O'Brien 49) are compared in the present study as changing these models did not affect the predictions much. The interaction term at the fluctuating velocity level is not well understood and so a new model (FET model) is proposed. Using the FET model, along with the Sinclair and Mallo 44 gas solid fluctuating velocity cross correlation expression, in the Eulerian model improves the prediction capabil ities of the Eulerian model. Also, for the large particles the Lun 53 vortex shedding model is included in the Eulerian model to capture the enhanced turbulence. When the combination of the Syamlal and O'Br ien 49 drag relation, the Peirano and Leckner 51 solid stress closure, the FET model for fluctuating interaction terms and the Sinclair and Mallo 44 gas solid fluctuating velocity cross correlation expression are used, minimum error in prediction of the various velocity profiles is obtained. This model has the capability to describe and predict flow behavior of not only gas solid flows (at high ST ) but also of liquidsolid flows in the inertia dominated flow regime (Chapter 3) PAGE 65 65 Table 2 1. Closure relations employed in the published gas solid flow models Model Interaction Term ( Ik, IT) Cross Correlation ( ksf) Drag Force ( F D ) Solid Stress ( s, ) Bolio et al. 7 TVBA Louge et al. 6 Wen & Yu 39 Lun et al. 50 MFIX TVBA Simonin 8 Syamlal & OBrien model 49 Peirano & Leckner 51 Zhang & Reese 46 VBA, Crowe 45 Zhang & Reese 46 Zhang & Reese 46 Peirano & Leckner 51 Benavides & van Wachem 9 TVBA Louge et al. 6 Simonin 8, Sinclair & Mallo 44 Wen & Yu 39 Peirano & Leckner 51 Present study TVBA Louge et al. 6 Koch & Sangani 42, Wiely et al. 43, Simonin 8, Sinclair & Mallo 44 Wen & Yu 39 Syamlal & OBrien model 49, Hill et al. 47, 48 Lun et al. 50 Peirano & Leckner 51 VBA, Crowe 45 Zhang & Reese 46 FET Sinclair & Mallo 44 New proposed model FET Sinclair & Mallo 44 Syam lal & OBrien model 49 Peirano & Leckner 51 Figure 21. Suspended particles enclosed in a box PAGE 66 66 Figure 22. Profile for the fluctuating energy transfer PAGE 67 67 Figure 23. Program flow chart PAGE 68 68 Figure 24 Comparison of solutions using various number of grids, (a) fluid mean velocity, (b) gas phase turbulence (c) solid mean velocity (d) granular temperature, for the case of polystyrene particles d = 200 m, Vfzcl = 11.3 m/s and m = 3.2 PAGE 69 69 Figure 2 5 Comparison between the Bolio et al. 7 code and the code used in the present study, (a) mean velocity profiles, (b) gas phase turbulence and granular temperature profiles, for the case of polystyrene particles d = 200 m, Vfzcl = 11.3 m/s and m = 3.2 Figure 2 6 Comparison of prediction in the fluctuation velocity between (a) Benavides and van Wachem 9 and (b) present study for the case of polystyrene particles d = 200 m, Vfz cl = 10.8 m/s and m = 3.2 PAGE 70 70 F igure 2 7 Comparison between the MFIX code and the present study, (a) mean velocity profiles, (b) gas phase turbulence and granular temperature profiles for the case of polystyrene particles d = 243 m, Vfz cl = 12.4 m/s and m = 3.2 PAGE 71 71 Table 22. Summary of e xp erimental d ata Data Size Density Mass Loading V fzcl (m/s) R Profiles Provided kg/m 3 Cases Cases mm 1 2 3 1 2 3 Tsuji et al. 2 243 1020 0.5 1.3 3.2 13.1 12.8 10.8 15.25 0.9, 0.9, 0.002 V fz v fz ', v sz 500 10 20 0.7 1.3 3.4 12.2 13.3 10.7 V fz v fz 1420 1030 0.6 2 2 13.4 12.8 13.2 V fz v fz 2780 1020 0.6 2.3 3.4 14.5 13.8 14.2 V fz v fz Jones et al. 3 70 2529 1 2 4 18.1 17.6 16.5 7.112 0.94, 0.94, 0.002 V fz V sz k, v sz Lee & Durst 5 400 2500 1.5 5.77 20.9 0.94, 0.94, 0.002 V fz V sz 200 2500 1.31 5.84 Sheen et al. 4 275 1020 1.224 8.785 26 0.9, 0.9, 0.002 V fz V sz v fz 450 1020 0.885 9.22 V fz V sz v fz 800 1020 1.5 9.686 V fz V sz v fz For all cases, f = 1.8 105 Ns/m2, f = 1.2 kg/m3 and 0 = 0.65 PAGE 72 72 T able 2 3. Time scale used in FET model for the different experimental data Experimental Data Size Density ~Re p ~Stokes Time Scale Vortex Shedding kg/m 3 No. No. Tsuji et al. 2 243 1020 45 75 D no 500 1020 130 320 c no 1420 1030 730 2700 c yes 2780 1020 1800 11300 c yes Jones et al. 3 70 2529 12 50 D no Lee & Durst 5 400 2500 90 170 c no 200 2500 5 43 D no Sheen et al. 4 275 1020 30 40 D no 450 1020 85 115 c no 800 1020 220 380 c yes PAGE 73 73 Figure 2 8 Present model predictions compared to Tsuji et al. 2 PAGE 74 74 Figure 2 9 Present model predictions compared to Tsuji et al. 2 Figure 2 10 Present model predictions compared to Tsuji et al. 2 PAGE 75 75 Figure 2 11 Present model predictions compared to Tsuji et al. 2 PAGE 76 76 Figure 2 12 Present model predictions compared to Jones et al. 3 PAGE 77 77 Figure 2 13 Present model predictions compared to Sheen et al. 4 PAGE 78 78 Figure 2 14 Louge et al. 6 crosscorrelation predictions compared to Tsuji et al. 2 (a) Figure 2 15 Koch and Sangani 42 crosscorrelation predictions compared to Tsuji et al. 2 PAGE 79 79 Figure 2 16 Wylie et al. 43 crosscorrelation predictions compared to Tsuji et al. 2 (a) Figure 2 17 Simonin 8 crosscorrelati on predictions compared to Tsuji et al. 2 (a) 243 PAGE 80 80 Figure 2 18 Sinclair and Mallo 8 crosscorrelation predictions compared to Tsuji et al. 2 Figure 2 19 FET model predic tions compared to Tsuji ** articles ** Tsuji Y. Private Communication 1993. PAGE 81 81 Figure 220 Louge et al. 6 cross correlation predictions compared to Tsuji ** particles Figure 221 Koch and Sangani 42 crosscorrelation predictions compared to Tsuji ** 243 ** Tsuji Y. Private Communication 1993. PAGE 82 82 Fig ure 22 2 Wylie et al. 43 crosscorrelation predictions compared to Tsuji ** particles Figure 223 Simonin 8 crosscorrelation predictions compared to Tsuji ** particles ** Tsuji Y. Private Communication 1993. PAGE 83 83 Figure 2 24 Sinclair and Mallo 44 crosscorrelation predictions compared to Tsuji ** 243 ** Tsuji Y. Private Communication 1993. PAGE 84 84 Figure 225 Comparing the magnitudes of ksg (a), Ik (b) and IT(c) for the various m = 3.2 2 PAGE 85 85 Table 2 4 Percentage Error in Predicted Gas Velocity Fluctuations Using the Various Combinations of Different Interaction Terms and Velocity Cross Correlations Interaction term Cross Correlation Average error% m=0.5 m=1.3 m=3.2 m=0.7 m=1.3 m=3.4 FET Sin clair & Mallo 44 3.65 9.14 5.10 1.03 4.30 2.58 4.30 TVBA Louge et al. 6 5.49 8.97 17.57 4.11 6.84 16.70 9.95 TVBA Koch & Sangani 42 6.53 16.37 26.99 5.15 9.40 25.82 15.04 TVBA Wylie et al. 43 6.53 16.37 27.00 5.15 9.40 25.83 15.05 TVBA Simonin 8 6 .47 15.59 18.69 5.14 9.34 24.36 13.27 TVBA Sinclair & Mallo 44 3.03 11.22 11.42 1.44 3.02 7.86 6.33 Table 2 5 Percentage Error in Predicted Solid Velocity Fluctuations Using the Various Combinations of Different Interaction Terms and Velocity Cross Correlations Interaction term Cross Correlation Average error% m=0.5 m=1.3 m=3.2 FET Sinclair & Mallo 44 4.25 8.05 4.29 5.53 TVBA Louge et al. 6 9.41 14.29 8.41 10.70 TVBA Koch & Sangani 42 10.29 15.29 9.56 11.72 TVBA Wylie et al. 43 10.29 15.29 9.56 11.72 TVBA Simonin 8 10.22 15.24 9.32 11.59 TVBA Sinclair & Mallo 44 2.87 8.42 3.93 5.07 Table 2 6 Percentage error in predicted gas velocity fluctuations using the various drag relations Models Average erro r% m=0.5 m=1.3 m=3.2 m=0.7 m=1.3 m=3.4 Wen & Yu 39 3.65 9.14 5.10 1.03 4.30 2.60 4.30 Hill et al. 47, 48 3.76 10.25 8.41 0.88 4.00 4.12 5.24 Syamlal and OBrien 49 3.65 9.05 4.34 1.01 4.22 2.30 4.10 PAGE 86 86 Figure 226 Comparing the magnitudes of FD (a) and (b) for the various drag models m = 3.2 2 Figure 227 Comparing the magnitudes of s/ f (a) and / f (b) for the various solid 2 PAGE 87 87 Table 2 7 Percentage error in predicted gas v elocity fluctuations using the various solid stress closures Models Average error% m=0.5 m=1.3 m=3.2 m=0.7 m=1.3 m=3.4 Peirano & Leckner 51 3.65 9.14 5.10 1.03 4.30 2.58 4.30 Lun et al. 50 3.77 9.53 4.97 1.03 4.28 2.7 1 4.38 Equation Chapter 3 Section 1 PAGE 88 88 CHAPTER 3 FLUIDIZED DILUTE TURBULENT LIQUID SOLID FLOW Background For s lurry flows t heir dynamic behavior and transportation is a frequent issue in the industry Flow variables such as mean velocities, fluctuating velocities and volume fraction of both phases are necessary for understanding the physical processes involved an d the energy distribution. However, for macro sized particles in a liquid, very little experimental data are available and even fewer modeling studies have been carried out. Reliable and d etailed measurements of the two phase flow behavior are difficult as intrusive techniques disrupt the flow, and due to the lack of data, model predictions cannot be validated. For gas solid flows; Lee and Durst 5, Tsuji et al. 2, and Sheen et al. 4 used laser Doppler velocimeter (LDV) a nonintrusive technique, to success fully measure mean and fluctuating velocities of both phases in a vertical pipe. In liquid solid flow Alajbegovic et al. 12 also used laser Doppler anemometer and a singlebeam ray densitometer to respectively measure the velocities and volume fraction of both phases accurately. Alajbegovic et al. 12 studied a system of water (Re ~ 3x104) laden with ceramic particles of size 2.32 mm with density 2450 kg/m3 and solid concentration of ~23% by volume, i n a 30.6 mm diameter vertical pipe. Pepple 13 has carri ed out experiments with glass particles of size 0.5 and 1 mm, and density 2500 kg/m3. These particle wer e fluidized by water at Re ~ 2x1055x105 with a solid concentration ~ 0.73% by volume, in a vertical pipe of diameter 78 mm. LDV/ Phase Doppler particle analyzer (PDPA) wa s used to measure the axial mean and fluctuating velocities of the liquid and solid phases. The experiments performed herein, PAGE 89 89 extend the work of Pepple 13. Glass particles of size 1.5 mm and density 2500 kg/m3 were studied in the same ex perimental setup as that of Pepple 13 under similar operating conditions ( Re ~ 2x1055x105, solid concentration ~ 0.73% by volume and pipe diameter is 78 mm), although a different LASER setup was used. The data obtained in the present study along with the experiments performed by Alajbegovic et al. 12 and Pepple 13 provide benchmarks data sets for a wide range of particle sizes which exhibit either viscous, transitional or collisional flow regimes. The f low regime of fluidsolid flows depends upon the Stok es number ( ST ) of the particles. p fST, ( 3 1 ) where p and f are the particle and fluid response times, respectively. 218s p fd, ( 3 2 ) f fzclD V ( 3 3 ) here s is the solid density, d is the particle diameter, f is the intrinsic fluid viscosity, D is the pipe diameter and Vfzcl is the center line fluid velocity. At low ST (small and light particles), the solid response time is smaller than the fluid response time, the particles tend to follow the fluid streamlines (Figure 3 1) and the two phase flow is considered to be in a viscous dominated regime. Chen and Wood 15 have described an Eulerian two fluid model for dilute turbulent gas solid flows in which partic les tend to emulate the fluid motion. PAGE 90 90 At high ST the particle response time is much larger than the fluid response time. The particles have large inertia and they tend to travel in straight lines ( Figure 3 1) and collide with each other and the confining vessel In this case, the flow is considered to be i n inertiadominated flow regime. Hadinoto and Curtis 14 extended the work of Bolio et al. 7 to include fl uid lubrication effects and employed their gas solid model to generate flow predictions correspondi ng to operating conditions in the liquidsolid data of Alajbegovic et al. 12. Finally, at intermediate ST particles neither follow the fluid, nor do they travel in straight lines but instead their path is governed by both fluid eddies and particle inerti a. Hence, in this case, particles have curved paths (Figure 3 1). To date no Eulerian twofluid models have been proposed to describe this flow regime at intermediate ST In the present study, two fluid gas solid flow models are extended to predict liquid solid flow behavior The Chen and Wood 15 model describes liquid solid flows with low ST (viscousdominated regime). An Eulerian two fluid model using the combination of Syamlal and OBrien 49 drag force relation, Peirano and Leckner 51 solid stre ss description and the FET (Chapter 2) interaction term describes liquid solid flows with high ST (inertia dominated regime). For the intermediate ST (transitional regime) a combination of the two fluid models at the low and high ST limits describes these heterogeneous slurry flows. All the model solutions are validated against the benchmark data sets obtained in the present study and by Alajbegovic et al. 12. Experimental Setup Flow Loop The experiments are carried out in a versatile pilot scale flow loop system. The flow loop was designed inhouse with capabilities to handle a wide range of particle PAGE 91 91 sizes and solid concentrations at very high velocities of about 7 m/s. Enormous amount of time was spent on troubleshooting the flow loop and ensuring the accuracy of the measurements 13. Figure 3 2 presents the schematic of the flow loop. The loop is constructed with a nominal three inch (78 mm), schedule 40, and type 304, stainless steel pipe. Water from the storage tank can be circulated through the loop via a 50 hp centrifugal pump controlled by a variable frequency drive (ABB model number ACH550UH072A 4). The water is first pumped through a venturi eductor which when open entrains particles into the loop. An electromagnetic flow meter which can accurately measure the singlephase volumetric flow rate of the liquid is also installed in the loop. N ext, the flow passes through the 3.84 m vertical pipe section and then through the borosilicate glass test section, 0.91 m in length and an inner diameter of 78 mm. After t hat the flow is deposited into a particle separator through a by pass section (closed during normal operation) and a sampling tank (open during normal operation). A particle screen, mounted on top of the separator, allows only the water to overflow back i nto the storage tank, while the particles settle in the standpipe arrangement. If the eductor is open, the particles flow back into the loop, and if the eductor is closed the particles collect in the standpipe. The liquidsolid flow can also be sampled for the purpose of determining the solid mass loading ( m ) To sample, during normal operation, the drain at the bottom of the sampling tank is closed so that the sample starts collecting in the sampling tank. Next, after a measured time interval, the flow is diverted directly to the particle separator using the by pass. The volume of water in the tank is measured and the particles are filter out, dried and weighed. The ratio of the weight of the solids collected to the weight of the PAGE 92 92 water is equal to the mass loading. More information on the design and operation of the flow loop and details on the sampling technique are described in Pepple 13. PDPA System A LDV/PDPA system, laser based nonintrusive technique accurately measures the instantaneous and time averaged velocities of both the fluid and particle phases as well as the particle diameter with high spatial resolution. It is necessary for the particles to be perfect spheres to get accurate particle size data. The LDV/ PDPA system uses a n ArgonIon laser (Sp ectra Physics Stabilite 2017 6W laser ) and a Braggcell to avoid directional uncertainty. The power of each of the probes i s 4050 mV for the laser The transmitting and receiving optics and the real time signal analyzers are manufactured by TSI/Aerometics The transmitting and receiving optics are mounted on a traversing A 363 mm focal length transmitting lens and a 300 mm focal length receiver lens are used for the LASER to measure the seed particles in the f luid phase. SinglePhase Validation Before twophase measurements we re carried out, singlephase validation i s performed to establish a base case. R aw data obtained for the singlephase liquid fluctuating velocity 13 at three different centerli ne velocities (3, 5 and 7m/s) are compared to the standard singlephase kturbulence model in Figure 3 3a. There is a consistent difference between the raw data and the singlephase model prediction for the singlephase liquid flow that is not present in single phas e gas flow 58. Pepple et al. 59 have compared various benchmar k measurements of single phase flows published in the literature using a variety of experimental techniques (hot wire and LDV ). The authors observe that there is considerable variation in the values of PAGE 93 93 singlephase gas and singlephase liquid fluctuating v elocity amongst these commonly cited references These variations exceed the generic errors associated with the measurement techniques. Further more Pepple et al. 59 also show that the fluctuating velocities in water are consistently higher than the fluctuating velocities in air at the same Re. Additionally, the authors also conclude that, in air, the scaled turbulent velocity exhibits no dependency on Re far away from the wall, although the same i s not true for water. In fact, in water the scaled turbulent ve locities increase as the Re increases. In the present work effects similar to that of Pepple et al. 59 were observed, but instead of changing the parameters of a well established singlephase kturbulence model, all of the experimental data are multiplied by a constant factor of 0.85 to refine the data. Figure 33b compares the refined data to the singlephase kturbulence model predictions Thus, to be consistent, all of the twophase data are also refined through multip lication with the constant factor of 0.85. Particles and Laser Setting s Mov ing on to the two phase data, sodalime glass beads of size 1.5 mm (Jaygo Incorporated) in water are measured in this study. F or the se particles two traverses across the pipe radius are required. The first measures the solidphase and the second the fluidphase. The 1.5 mm particles also require larger fringe spacing, a transmitting lens of focal length 762 mm and a receiv er lens of 1000 mm focal lengt h for accurate size measurements The relevant optical settings are summarized in Table 3 1. Further information on the optical settings, the LDV/ PDPA system and the flow measurements are detailed in Pepple 13. PAGE 94 94 Table 3 2 details the various operating condi tions including Stokes numbers, system parameters, and properties of the glass particles used in the present study as well as the ceramic particles in the study by Alajbegovic et al. 12. and the glass particles used by Pepple 13. The values of the particl e particle coefficient of restitution ( e ) and particle wall coefficient of restitution ( ew), defined in a liquid 60 and the value of specularity coefficient ( ) are also given in Table 32. 27 experiments (0.5 and 1 mm, Pepple, 13 and 1.5 mm present study ) were conducted. For each of the particle size, the mass loading was held constant while the fluid centerline velocity varied. Three mass loading levels ( 0.0175, 0.0425 and 0.0750) were investigated for each particle size. In Alaj begovic et al. 12 study three experiments were conducted, each with a different mass loading and centerline fluid velocity. T he raw data for the 1.5 mm particles for the operating conditions described above is presented in Appendix B TwoPhase Data Validation Four replicates were performed t o gauge the reproducibility of the measurements of the solid velocity for the case of d = 1.5 mm with the operating conditions m = 0.0425 and Vfzcl = 3 m/s. Figure 3 4a shows that the deviation in the mean solid velocit y is extremely small (~1%). But Figure 3 4b shows that the deviation in the fluctuating solid velocity is ~5% in the core of the pipe and at the wall it i s ~10%. Near the wall, the laser light gets scattered and the resulting sign al interferes with the signal f r om the particles leading to distortion of the velocity peak at the wall. Hence, the data near the wall needs to be subranged appropriately. This subranging may lead to the larger errors near the wall versus the core of the pipe. PAGE 95 95 Eulerian Two Fluid Models for Dilute Turbulent LiquidS olid Flows T he fluid and solid fluctuating velocity data for the 0.5 mm particles 13 ( ST ~ 3) are approximately equal indicating that the particles follow the details of the fluid and hence exhibit viscous dominated flow behavior. While, for 2.32 mm particles 12, with ST > 40, the solid fluctuating data is lower than that of the fluid indicating an enery loss for the solid paricles due to collisions. This exhibits inertiadominated flow behavior. For the 1 mm 13 and 1.5 m m particles (5 < ST < 30) however, the solid fluctuating velocity is higher than that of the fluid exhibiting transitional behavior. Thus Stokes numbers 5 and 40 are assumed as the limits of the viscous and inertia regimes, respectively. The experimental data are described in detail in the results and discussion section of this C hapter High ST, InertiaDominated Flow Regime ( ST > 40) The continuum equations for the Eulerian model description of a dilute turbulent liquid solid flow model are being based on th e gas solid model discussed in Chapter 2. Equation 2 1 to equation 2 6 describe the continuum equations used for the liquid solid flow. The liquidphase stress rz is also described in the same fashion as that f or the gas solid flow (Chapter 2, Mathematical Model, Fluid Stress). The drag f orce ( FD) equation 2 14 is still valid for t he case of liquidsolid flow and for the drag coefficient is computed using t he Syamlal and OBrien 49 ( equation 2 30 to equation 2 34) drag relati on. For the case of high ST in liquid solid flow, granular kinetic theory is still applicable and hence, Peirano and Leckner 51 can be used for computing granular energy dissipation ( equation 2 35 to equation2 38), solid phase stresses ( rz, equation 2 46, rr and equation 2 43 to equation 2 45) and granular energy flux PAGE 96 96 ( qpTr, equation 2 47). Equatio n 2 57 to equation 2 51 are used for calculating s and equation 2 53 to equation 2 64 are used for calculati ng Finally, the fluctuating interaction terms ( Ik and IT) are bas ed on the FET model (Chapter 2 Mathematical Model, Fluctuating Interaction term) and the liquid solid fluctuating velocity cross correlation is based on the Sinclair and Mallo 44 model ( e quation 2 88). Liquid solid flows exhibit the phenomenon of vortex shedding at larger Rep than the gas solid flow. The only experimental data, of liquid solid flow, which exhibits vortex shedding has the operating conditions d = 2.32 mm Vfzcl = 1.81 m/s, m = 0.0556 12 and Rep ~ 540. The rest of the experimental conditions have Rep less than 500 and do not exhibit vortex shedding. The Lun 53 vortex shedding model i s modified for the liquidsolid flow as, 212wt wCk E d ( 3 4 ) Reffzsz p fdVV ( 3 5 ) 0.017Ret pf for Re500p, ( 3 6 ) 10/3wC for Re500p. ( 3 7 ) For the high ST cases, t he boundary conditions follow the gas solid model (Chapter 2 Boundary Conditions). The fluid center line velocity ( Vfz cl) is used as the operating condition for the fluidphase, and f or the solidphase, solid mass loading m is used as the operating condition. For the high ST model the solid volume fraction ( ) is calculated from r component of the solid phase momentum balance equation ( equation 2 3 ) as the equation of continuity for the solid phase yields a redundant equation. PAGE 97 97 Low ST, Viscous Dominated Flow Regime ( ST < 5) In two phase flows with low ST model particles are affected by the details of the fluid motion, and particles do not move in straight lines through turbulent eddies. Hence, granular kinetic theory is not appropriate to describe soli d phase fluctuations and the resulting solid phase stresses. Rather, sin ce the particles follow the details of the fluid motion, it is reasonable to assume that the fluid fluctuations and the solid fluctuations are similar Tk ( 3 8 ) The experimental data obtained in the Pepple 13 validate this assumption. In the case of low St okes number flow the Chen and Wood 15 expressions can be applied to describe the solidphase viscosity and the corresponding solid phase shear stress. (1)T sef p e, ( 3 9 ) where e is the time scale of the fluid eddy, For the fluctuation interaction term Ik, the Chen and Wood 15 equations are also followed. 12k sfIkk ( 3 10) where ksf is the fluidsolid velocity crosscorrelation, 2exp0.0825p sf ekk ( 3 11) In the case of low ST flow the solid continuity equation of Chen and W ood 15 reveal that the solids volume fraction i s not a function of the radial position. Hence, the PAGE 98 98 solid s mass loading operating condition ( equation 2 105) ca n be solved to yield the constant value for the solid volume fraction. (1)ssz ffzVrdr mVrdr, ( 3 12) which can be rewritten as, 1 1ssz ffzVrdr mVrdr ( 3 13) The mean fluid and solid v elocities as well as the fluid phase turbulence, follow the no slip condition at the wall Intermediate ST, T ransitional Flow Regime (5 < ST < 40) For intermediate ST flow, there are no published modeling studies For this regime, the modeling framework is similar to the high ST case. However, modifications are made to the granular dissipation ( ), the solidphase viscosity ( s), the solidphase conductivity ( ) and the time scale for fluctuation energy transfer ( sf ) since traditional granular kinetic theory as originally developed is not fully appropriate to describe the details of the solid motion in this flow regime. While the particles in this regime do have some degree of inertia, it is not large and hence the loss of energy on collisions is not that high. Thus, in the intermediate Stokes number regime, the granular dissipation ( ) is set to zero. 0 ( 3 14) PAGE 99 99 The solid viscosity and the solid conductivity are obtained by using a ST weighted combination of their respective expressions from the lo w ST model 15 and the high ST m odel 51. 22 40 5 405 405 (1)T s ef skc p eST ST GG, ( 3 15) 3340 5 405 405 (1)T ef skc p k eST ST GG ( 3 16) For the time scale for fluctuation energy transfer, a modified expression of the high ST model is developed based on the available experimental data in the intermediate ST flow regime (1 mm and 1.5 mm particles) 76.7*10Re370.62(1)f sfd D ( 3 17) here Re is the flow Reynolds number Refzclf f VD ( 3 18) The closure for the FET time scale (equation 3 17) is proposed such that, as the diameter of the particle increases, at a constant m and Vfzcl, the difference between the fluid and solid fluctuations increase. Also, upon increasing the operating Re, at a constant d and m, t he difference between the fluid and solid fluctuating velocity decreases. By i ncluding the drag coefficient in the FET time scale, a mass loading dependency is also obtained. This dependency is such that on increasing the mass PAGE 100 100 l oading, at a constant d and Vfzcl, the difference between the fluid and solid fluctuations decreases. All of these qualitative trends are consistent with experimental observations. Looking to the granular energy balance, the dissipation rate is set to zer o, the solid phase strain is small, and hence the generation term sz rzV r is small. Flow predictions in the intermediate ST regime strongly depend on IT and, consequently, they are dependent on f sf. Finally, the solid volume fraction is similarly based on a ST weighted combination of the resulting profiles from the low ST model ( equation 3 13) and the high ST model e xpress ions (equation 2 3 where solid phase normal stress are computed using 51). Intermediate High 40 1 5 405 405 1ST ST ssz ffzST ST Vrdr mVrdr ( 3 19) The boundary conditions applied for intermediate ST flow are similar to those for high ST flow. The only change is that the solidphase str ess and conductivity that appear in the solid velocity and granular temperature boundary condition follow equations 3 15 and 3 16 for the intermediate ST flow as mentioned above. Table 3 3 presents a summary of the model equations and constitutive relations for the hi gh, low and intermediate ST cases Results and D iscussion In Chapter 2 (Results and Discussion) for gas solid flows it is assumed that fluctuations in velocity are isotropic in all directions for the solid phase, and for the fluid phase the radial and azimuthal fluctuations in velocity 'frvand 'fv respectively, are PAGE 101 101 approximately half of the axial velocity fluctuations Similar assumptions are made for the liquidsolid flow cases too. Thus enabling the computation of the solid and fluid fluctuating velocities from T ( equation 2 106) and k ( equation 2 108 ) profiles, r espectively. Once the axial fluctuating velocity profiles are related to the corresponding k and T profiles, the flow predictions can be directly compared to the experimental data. High ST, InertiaDominated Flow Regime As indicated in Table 3 2, the data of Alajbegovic et al. 12 with ST > 40 falls in the inertia dominated flow regime. Figure 3 5 compares the predicted mean velocities of the fluid and solid phases from the high ST model to the Alajbegovic et al. 12 data, while Fig ure 3 6 compares the fluctuating velocities f or the liquid and solid phases. From the data for the high ST conditions a visible slip between the mean velocities of the liquid and solid phases (Fig ure 3 5) is clearly observed. This is indicative of the significant inertia the large particles ( d = 2.32mm) Th e Syamlal and OBrien 49 drag force relation captures the slip between the two phases for all the three operating conditions fairly well. In Figure 3 6 the experimental data further indicates that the granular temperature is lower than the liquidphase turbulence. I n the inertia dominated regime, particles with large inertia on collision lose energy, the granular dissipation increases and, as a result the granular temperature reduces. The application of the Peir ano and Leckner 51 model not only captures the reduction of the granular temperature, but also the shape of the T profile for all the three operating conditions. This significant slip velocity and reduction of granular temperature are not seen in the experimental data for t he viscous or transitional flow regimes. PAGE 102 102 Low ST, ViscousDominated Flow Regime F r om the experiment al data of 0.5 mm particles 13 ( ST < 5 ) the flow is anticipated to be in viscous dominated flow regime. Experimental data for the three centerline fluid oper ating conditions and the highest mass loading conditions in the core of the pipe, was difficult to obtain (Figure 3 7) since the particle number density ass ociated with these mass loadings was very large, leading to poor signal quality and extremely low d ata rate. From Figure 3 7 it can be observed that there is marginal slip between the mean liquid and solid phase velocities which is correctly captured by the low ST model. Further more the fluid velocity fluctuation essentially match the solid velocity f luctuation for all of the nine operating conditions (Figure 3 8), indicating that the particles are following the details of the fluid motion. The magnitude and the shape of mean and fluctuating velocity profiles for both phases at all nine operating condi tions are well predicted by the Chen and Wood 15 model. In addition, both the experimental data and the flow predictions reveal that there is very little effect on the magnitude and shape of the mean and fluctuating velocities as solid loading is increased from 0.0175 to 0.075 at a constant fluid centerline velocity Additionally, as the centerline fluid velocity is increased at a cons tant mass loading, and the magnitude of the liquid and solid velocity fluctuations increase (Figure 3 8). These trends are s imilar to those observed in singlephase fluid flow. Intermediate ST, Transitional Flow Regime Operating conditions associated with the 1 mm 13 and 1.5 mm particles have ST between 5 and 40, and exhibit transitional flow behavior. Figure 3 9 ( d = 1 mm) and Fig ure 3 11 ( d = 1.5 mm) present the experimental data and model predictions for the mean fluid and solid velocities of 1 mm and 1.5 mm particles, respectively. Figure 310 PAGE 103 103 and Figure 312 present the corresponding fluctuating fluid and solid velocities. Although there is a small slip in the mean velocity between the two phases (Figure 39 and Figure 310) large difference in magnitude in the solid and fluid velocity fluctuations exists (Fig ure 3 10 and 3 12). For the 1mm particles, as the solids loading i s increased from 0.0175 to 0.075 at a constant velocity, minimal change is seen in the mean velocity profiles (Figure 3 9) However, the magnitude of the difference in velocity fluctuations of both phases tends to decrease with the increasing mass loading (Fig ure 3 10). A s the centerline fluid velocity is increased at constant mass loading, there is no t r e nd in the fluid and solid velocity fluctuations In most cases the difference in magnitude between the fluctuation velocities of the two phases decreases while in some cases the difference increases. The mean velocity behavior of the 1.5 mm particles (Figure 311) is similar to that of the 1 mm particles. T he solid velocity fluctuations associated with the 1.5 mm particles are significantly larger than those associated with the 1 mm particles and greatly exceed the magnitude of the fluid velocity fluctuations As the solids mass loading is increased from 0.0175 to 0.075 at a constant velocity, the experimental data show that the magnitude and difference bet ween the velocity fl uctuations of both phases changes minimally (Fig ure 3 12). Also, as the centerline fluid velocity is increased at a constant mass loading, the magnitude of the velocity fluctuations f or both phases increases but the difference between the velocity fluctuations of both phases remains the same ( unlike the 1 mm particles ) As the particle size is increased from 0.5 mm to 1mm 13 and finally to 1.5 mm (at constant m and Vfzcl), minimal changes are observed in the mean velocity profiles PAGE 104 104 (Fig u re 3 7, 3 9 and 3 11). In fact, in all three cases marginal slip is observed because the Re is very high. Hadinoto et al. 61 experiments showed that as the Re is increased, the slip velocity reduced. This effect is also seen in Alajbegovic et al. 12 data (Fig ure 3 5) and the model also shows same behavior. On the other hand, the fluctuating velocities of both phases show different behavior for the different particle sizes. For the 0.5 mm particles there is no difference between the fluctuating velocities of both phases, while this difference increases dramatically a s the particle size increases from 0.5 mm to 1 mm 13 and then to 1.5 mm with the solid velocity fluctuations exceeding the fluid velocity fluctuations For the 2.32 mm particles 12, however, the solid velocity fluctuations are lower than that of the fluid. As mentioned earlier model predictions in the transitional flow regime are strongly dependent on IT. The closure proposed for IT (equation 3 17) is consistent with the following observed flow behaviors: As the particle diameter increases, at constant m and Vfzcl, the difference between the fluid and solid fluctuations increases As the flow Re of flo w increases, at constant d and m, the difference between the fluid and solid fluctuations decreases As t he mass loading increases, at constant d and Vfzcl, the difference between the fluid and solid fluctuations decreases In the fluctuating int eraction term, the dependency on the mass loading is through drag coefficient Although the model predict ions for the magnitude of the solid velocity fluctuations are not in perfect agreement with the experimental measurements, the intermediate ST model is able to capture the shape of the solid velocity fluctuation profile. This agreement indicates that the weighted expressions for the solidphase viscosity and PAGE 105 105 conductivity ( which govern the shape of this profile ) in equation 3 15 and 3 16 work well in this intermediate ST flow regime Summary There are many practical difficulties and li mitation due to which accurate measurement of data for liquidsolid turbulent flows is difficult. In this study LDV/ PDPA is used to measure the mean and fluctuating velocity of solid and liquid phases in a versatile, pilot scale flow loop. The flow loop setup can be used to for measuring flow properties of a wide range of particles (0.5 and 1 mm, Pepple 13, and 1.5 mm, present study) at a variety of flow rates ( Vfzcl = 3 m/s to 7m/s). The data obtained herein along wit h the data of 2.32 mm particles 12 and in Pepple 13 span all three flow regimes (v iscous, inertia and transitional flow regimes). The Eulerian two fluid model for dilute turbulent gas solid flows with particle particle int eraction described in Chapter 2 is used for predicting velocity profiles of the flows in the inertia dominated regime, Chen and Wood 15 model is used to describe flow in the viscous dominated flow regime, and a bridge model (using expressions from the high and low ST model) is proposed for the transitional flow regime. The predictions obtained for all three regimes compare well against the experimental data In the industry, generally the solid concentration is ~ 5% by volume for slurries, however the laser based setup can only measure a solid concentration up to 3% by volume, due to signal attenuation and low data rate. Efforts are been made to acquire data for higher mass loading conditions by matching the refractive indices of the liquid and solid phases. Even then, acquiring detailed dat a for dense liquidsolid flow (solid concentration > 10% by volume) using LDV based technique is difficult. Hence in the present work dense fluidsolid flows are not studied. Dense fluidsolid flows transition to PAGE 106 106 fluidized beds at low operating velocities Fluidized beds have many industrial applications ( Chapter 1 ) Some aspects of fluidized beds have also been studied in C hapter 4 and Chapter 5. PAGE 107 107 Fig ure 3 1 Effect of ST on the solid path in the fluid PAGE 108 108 Fig ure 3 2 Flow loop (A) w ater t ank ( B) p ump ( C) venturi e ductor ( D) e lectromagnetic f low m eter (E) test s ection (F) by p ass (G) s ampling t ank (H) p article screen (I) p article separator (modif ied from the original in Pepple 13) PAGE 109 109 Fig ure 3 3 Comparing singlep hase f luctuating velocity 13 to the sing le p hase kt urbulence m odel ( fz vk ) (a) r aw d ata (b) r efined d ata u sing the 0.85 f actor PAGE 110 110 Table 31 PDPA s ettings LASER setup Liquid Glass < 0.2 mm 1.5 mm Transmitting Optics Wavelength (nm) 514.5 514.5 Focal Length (mm) 36 3 762 3.7 7.8 LASER Beam Diameter (mm) 2.65 2.65 LASER Beam Intersection Angle 7.97 3.78 LASER Beam Separation (mm) 50 50 89.7 188 Probe Volume Length (mm) 1.29 5.71 Probe Volume Height (mm) 0.090 0.188 Pro be Volume Width (mm) 0.090 0.188 Receiving Optics 300 1000 150 150 Collection Angle 150 150 Software Setting High Pass Filter (MHz) 20 20 Frequency Shift (MHz) 36 36 Sampling Rate (MHz) 10, 20 10, 20 Low Pass Filter (MHz) 5, 10 5, 10 Burst Threshold (mV) 0.2 0.2 PAGE 111 111 Table 32 Summary of experimental data Data Density Size Mass Loading ( m ) V fzcl in m/s R ~ST Regime kg/m 3 mm Case mm 1 2 3 Pepple 13 2500 0.5 0.0175 3.12 4.93 7 .14 78 1.44 Viscous 0.0425 3.17 5.10 7.03 2.30 Viscous 0.0750 3.17 5.10 7.30 3.27 Viscous 2500 1 0.0175 3.11 5.04 7.06 5.66 Transitional 0.0425 3.02 5.03 7.04 8.98 Transitional 0.0750 3.16 4.68 6.83 12.8 Transitional Pres ent Study 2500 1.5 0.0175 3.12 5.04 7.00 78 12.5 Transitional 0.0425 2.97 4.91 7.05 20.5 Transitional 0.0750 3.04 5.04 7.02 28.9 Transitional Alajbegovic et al. 12 2450 2.32 0.0556 1.81 30.6 43.4 Inertial 0.0584 2.33 55.8 Inertia l 0.0670 2.83 67.8 Inertial For all cases: f=0.001 Ns/m2, f=1000 kg/m3, 0=0.65, e =0.54, ew=0.54, =0.002 Figure 34 Solid fluctuating velocity measurement for d = 1.5 mm, m = 0.0425 and Vfzcl ~ 3 m/s (averaged over four r uns) PAGE 112 112 Table 33 Summary of model equations in the high, low and i nt ermediate ST flow regimes V ariable High ST Low ST Intermediate ST Eq. 2 3 Solid phase continui ty equation (Chen & Wood 15) Eq. 3 19 Vfz Eq. 2 1 Eq. 2 1 Eq. 2 1 V sz Eq. 2 2 Eq. 2 2 Eq. 2 2 k Eq. 2 5 Eq. 2 6 Eq. 2 5 Eq. 2 6 Eq. 2 5 Eq. 2 6 T Eq. 2 4 Eq. 3 8 Eq. 2 4 Eq. 2 35 NA Eq. 3 14 s Peirano & Leckner 51 Eq. 3 9 NA Eq. 3 15 Eq. 3 16 I k I T Eq. 2 102 Eq. 2 103 Eq. 3 10 NA Eq. 2 102 Eq. 2 103 ssf sD for 40< ST <100 sc for ST >100 (FET model) Eq. 3 17 k sf Eq. 2 86 Eq. 3 11 Eq. 2 86 PAGE 113 113 Fig ure 3 5 Mean velocity ( Vfz and Vsz) m easurement for d = 2.32 mm 12 compared to flow predictions for high ST case Fi g ure 3 6 k, T measurements for d = 2.32 mm 12 compared to flow prediction for high ST case PAGE 114 114 Fig ure 3 7 Mean velocity ( Vfz and Vsz) measurements f or d = 0.5 mm 13 compared to flow predictions for low ST case PAGE 115 115 Fig ure 3 8 Fluctuating velocity ( vfz and vsz ) measurements for d = 0.5 mm 13 compared to flow predictions for low ST case PAGE 116 116 Fig ure 3 9 Mean velocity ( Vfz and Vsz) measurements f or d = 1 mm 13 compared to flow predictions for intermediate ST case PAGE 117 117 Fig ure 3 10 Fluctuating velocity ( vfz and vsz ) measurement for d = 1 mm 13 compared to flow predictions for intermediate ST case PAGE 118 118 Fig ure 3 11 Mean velocity ( Vfz and Vsz) measurement f or d = 1.5 mm compared t o flow predictions for intermediate ST case PAGE 119 119 Fig ure 3 1 2 Fluctuating velocity ( vfz and vsz ) measurement for d = 1.5 mm compared to flow predictions for intermediate ST case Equation Chapter 4 Section 1 PAGE 120 120 CHAPTER 4 BINARY SEGREGATION I N FLUIDIZED BED Background Fluidized bed with monosized particles have been studied in great detail 1 in the recent past. Nev ertheless, fluidized beds with multi solid mixtures are not well understood. This is due to the fact, multi solid fluidized beds on fluidization tend to segregate, thereby disturbing the flow patterns and leading to poor efficiencies. The phenomenon of par ticle segregation needs to be well understood in order to tackle the issue. It is widely accepted that segregation in a gas fluidized binary mixture occurs due to gas bubbles that are generated as soon as the bed leaves the packed state. Rowe et al. 62 showed that one of components of the binary mixture is selectively dragged upwards in the wake of the bubble to form the upper layer of the bed. This component is termed the flotsam. The second component, which sinks to the bottom of the column, is called t he jetsam. Despite the previous work, it is still not clear whether bubbling is responsible for the equilibrated segregation pattern, or if it is merely the mechanism that gives rise to segregation 16. Nienow et al. 63 experimentally studied the segregat ion behavior of nearly forty binary mixtures, although they reported on only a few of these. They concluded that mixtures with equal densities, but different sizes, mix easily. In addition, mixtures with equal sizes and with substantially different densiti es do not mix easily. They proposed a mixing index, M which is defined as the ratio of the average mass fraction of the jetsam in the uniformly mixed section to the average mass fraction of the jetsam through the PAGE 121 121 column (the exact definition of M is given in 63 64. M is a function of operating superficial gas velocity, U 1 [1exp()] M Z ( 4 1 ) 1 [1exp()] M Z ( 4 2 ) where, .TO mflowerUU Z UU ( 4 3 ) The parameter UTO is the takeover velocity, i.e. the superficial gas velocity at which the mixing index M is equal to 0.5, and Umf lower is the smaller of the two minimum fluidization velocities of each individual material. Equation 4 1 works well for mixtures in which the volume fraction of jetsam is less than 50% (volume concentration), and for particle size ratios less than three. Hence, equation 4 1 cannot describe the behavior of extremely disparate mixtures. Also, these studies did not consider the pressure drop profile of the mixtures. Nienow et al. 64 presented experimental data for not only binary mixtures, but also tertiary and quaternary systems. They studied the effect of different gas distributors on segregation patterns. Their conclusions for binary mixtures were similar to those of Nienow et al. 63. Garcia et al. 22 studied the effect of particle density ratio, r, for a binary mixture while maintaining the particl e size ratio ( dr) at unity. They observed that glass alumina and aluminapolyethylene systems ( r = 1.5 and 1.6, respectively) mixed better than the glass polyethylene system ( r = 2.39). PAGE 122 122 Hoffmann et al. 65 and Huilin et al. 23 performed experiments as wel l as numerical simulations of binary mixtures, mostly with the same particle densities but different particle sizes. The mixture for which dr = 1.5 65 mixed well whereas the other mixtures (1.9 < dr < 3.5) mixed completely only at relatively high gas veloc ities. The phenomenon of layer inversion in gas fluidized binary mixtures 19 has also been observed with some particle mixtures involving smaller, denser particles and coarser, less dense particles. In such mixtures, the smaller, denser particles behave as jetsam at low velocities just above the point of fluidization, but at higher velocities the coarser, less dense particles behave as jetsam. The authors associated this phenomenon with the change in mixture bulk density as a function of composition and gas velocity. Marzocchella et al. 18 studied transient fluidization behavior of a mixture consisting of particles with the same density, but very different sizes (dr = 4). The pressure drop profile from an initially mixed state for this system showed multiple peaks of fluidization (this phenomenon is described in more detail in the following sections) before the entire bed was completely fluidized and the pressure drop across the bed reached a steady value. The mixing index for this mixture initially increased as the gas velocity increased, but then reached a maximum value and decreased as the velocity of the gas further increased. In other words, the mixture started to segregate at relatively high velocities. Olivier et al. 66 investigated the transient fluidi zation behavior of less disparate mixtures. They observed that the pressure drop profile of these mixtures also had multiple peaks, but the degree of segregation reduced with an increase in velocity. PAGE 123 123 Formisani et al. 16, 67 studied mixtures with intermediate (1.8 < dr < 4) disparity and reported pressure drop profiles with multiple peaks when the mixtures were initially well mixed. The number of peaks reduced as the disparity between the two particle types diminished and, for two similar component systems, the pressure drop profiles showed a single peak before complete fluidization. The pressure drop profile from an initial segregated state of mixtures with intermediate disparity showed that the mixtures fluidized at two distinct points. In order to relate t he behavior of all the mixtures, the authors defined two velocities, Uif, which is the point of incipient fluidization, and Uff, the point at which the mixture completely fluidizes. Joseph et al. 17 performed experiments with mixtures having little dispari ty (1 < dr < 2). These mixtures fluidized at a single velocity point and showed segregation behavior similar to mixtures with a single peak in their pressure drop profiles 16, 67. The experiments performed in the present study involve both size segregating ( dr is varied from 1.4 to 6.6, in small increments) and density segregating ( r=3.1 and 7.4) mixtures. Pressure drop profiles, from initially mixed and segregated states, are reported along with axial segregation profiles at different velocities. These experiments encompass all of the previously reported types of pressure drop prof iles and segregation patterns. Table 41 summarizes the work of previous researchers as well as the experiments performed in the present study. The table presents the mixture type, the parameters studied, and whether or not the pressure drop profiles are r eported by the authors PAGE 124 124 The variety of possible pressure drop profiles, i.e ., two point fluidization, multi peak behavior, singlepeak behavior, and single point fluidization, have made it difficult to define the minimum fluidization velocity for a particl e mixture. In fact, the pressure drop profile is governed by the initial axial segregation profile. As a result, a variety of minimum fluidization velocities or incipient fluidization velocities can be defined based on the variety of pressure drop profiles The complete fluidization velocity, Uc, is perhaps the only parameter which remains sufficiently independent of the initial segregation profile. So, it is preferable to define the fluidization velocity of a mixture by its Uc the velocity at which the e ntire mixture, including both the jetsam and flotsam, is fluidized. Further, a parameter which would help categorize these different types of mixtures is also needed. Juxtaposing particle size ratio and particle density ratio is anticipated to be useful in predicting the fluidization behavior of these mixtures. As a result, parameters such as the ratio of the Archimedes numbers or the ratio of the minimum fluidization velocities of the individual components are expected to be of interest. Experiments Partic les and T heir Preparation Experiments were performed with glass (Mo Sci Corporation), polystyrene (Norstone, Inc.), and steel particles (Ervin Industries) with mean diameters ranging from f the particles are in the Geldart B class and on fluidization only exhibit either the bubbling or slugging flow regime. Thus, in all the cases when Umf is being compared, Umb (minimum bubbling velocity) is also being correspondingly compared. The particles were carefully sieved, dried in an oven for twelve hours and subject to an antistatic bar in order to eliminate PAGE 125 125 accumulated electrostatic charge. The particles were stored in a desiccator so that cohesive forces due to moisture and static effects were mi nimized. Table 42 summarizes the experimental materials and their properties. The notation used to represent the particle type has a letter indicating material type (G for glass, P for polystyrene, and S for steel) followed by the mean particle size in mi crons. The error in the measurement of the minimum fluidization velocities of all the particles was less than 10%. Table 43 presents the studied mixtures and their properties (composition, size ratio, density ratio, Archimedes number ratio, and the ratio of the minimum fluidization velocities). The notation followed for the particle mixtures is in two parts. The first part describes the jetsam percentage composition (by mass), material type (glass, polystyrene and steel), and particle size. The second part of the notation describes all of the same properties, but for the flotsam. Experimental Setup A fluidization segregation unit (Jenike and Johanson, Fluidization Material Sparing Segregation Tester, was used in the experiments. The tester has a column diameter of 1.6 cm and a height of 9.5 cm. There is a sliding disc assembly at the base of the column which can be used to divide the bed into multiple horizontal sections. Each of these sections may be transferred to a sampling container, one at a time, vi a a carousel arrangement. Details concerning the operation of the tester are given in Hedden et al 68 and ASTM D 6941 69. A schematic of the experimental set up is shown in Figure 41. A sintered metal columns. The air enters the column from the bottom, with its flow rate controlled by a PAGE 126 126 mass f low controller. The pressure drop across the entire setup was measured using a pressure transducer. The instantaneous pressure drop and velocity data were recorded on a computer. Experimental Procedure Prior to running an experiment, air was passed through the empty column to get a background pressure drop due to the column, diffuser, and the filter sections. In order to obtain the pressure drop profile for a segregated state, the material expected to be the jetsam was first weighed and a very small amount of antistatic powder (Larostat HTS 905 S, BASF Corporation, approximately 2 mg) was mixed with the particles and loaded into the column from the top. Next, the material expected to be the flotsam was weighed and antistatic powder was mixed into it and loaded into the column. The height, H to which the column was filled was recorded. The fixed bed height, H was approximately 4 cm for all experiments. The velocity was slowly increased at a rate of 0.0833 cm/s2 to a velocity much greater than that required to completely fluidize the mixture. The velocity was then decreased to zero at the same rate. For the pressure drop profile measurements from a mixed state, fresh amounts of the jetsam and flotsam, along with the antistatic powder (approximately 5 mg), were completely mixed either manually for disparate mixtures (mixtures with 2point fluidization), or for similar mixtures (mixtures with 1 point fluidization) complete mixing was obtained by maintaining the mixture at high air velocities (three times the complete fluidization velocity) for thirty minutes. In order to obtain the segregation profiles for the mixture at different velocities, the mixture was first completely mixed by following the same procedure used in the PAGE 127 127 pressure drop profile measurements from a mixed state. The fluidized bed was then maintained at the intended velocity for thirty minutes. At the end of the thirty minutes, the velocity was suddenly set to zero, and the bed collapsed to a fixed bed state. Next, the bed was sectioned axially by using the sampling disc assembly at the base of the column. Each section was collected in a sampling container via the carousel arrangement and its composition was analyzed. The composition for mixtures with different sizes was obtained by sieving, while m ixtures involving steel and glass or polystyrene were separated using magnets. Mixtures of glass and polystyrene were chosen such that they could be easily separated by sieves. Segregation Index Multiple definitions of an axial mixing index or segregation index for binary mixtures have been proposed 17, 63. The segregation index usually varies between zero and one, with zero indicating no segregation or uniform mixing and one indicating a completely segregated mixture. In order to define a segregation inde x, the feed composition of the jetsam by weight, xf, is first obtained as, fmassofjetsam x massofjetsammassofflotsam ( 4 4 ) In a similar manner, the final jetsam composition of the mixture, xm, is defined as the weight fraction of the jetsam at the end of the segregation experiment. The two weight fractions are defined in order to account for the material losses (losses in the present work are approximately 5% by weight of the feed composition). PAGE 128 128 Next, the weight fraction of the jetsam for each axial section of the column, xi, is determined. The segregation index, sii, for the ith axial section of the column is defined as, 21im i fxx si x ( 4 5 ) for sections in which xi > xf and, 2 im i fxx si x ( 4 6 ) for sections in which xi < xf. Finally, the overall segregation index, SI is computed as 1/2*.i imass of section SIsi mass of the entire bed ( 4 7 ) Since SI is a function of velocity, segregation profiles were obtained at various velocities and the corresponding segregation indices were calculated for each profile. It was observed that the reproducibility of the experiments increased as the disparity between the two components decreased, especially for size segr egated mixtures. Hence, it was only necessary to conduct three replicate experiments for two of mixtures with the greatest disparity (mixture numbers 1 and 2 in Table 43). Only a single set of experiments was carried out for the other mixtures. The error in the value of SI was found to be between 5 and 10%. Results and Discussion Table 44 provides a detailed summary of the mixture parameters reported in the literature, as outlined in Table 41, as well as the results from the experiments performed in the present study, as outlined in Table 43. The mixtures are arranged in descending order with respect to the jetsam to flotsam minimum fluidization velocity PAGE 129 129 ratio. This velocity ratio, Ur, is a good measure of particle mixture disparity as discussed later in this section. Throughout the published literature, as well as in this study, a wide variety of segregation behavior associated with different types of mixtures has been observed. Here, an attempt is made to qualitatively categorize these various mixtures based on the density ratio, size ratio, and the ratio of minimum fluidization velocities of the individual components. Although all of the published studies provide information on segregation profiles, not all include pressure drop information. Hence, in some cases, categorizing the mixtures involves hypothesizing some aspects of the mixture behavior based on other reported behavior for those and similar mixtures. By analyzing the pressure drop, flow, and segregation behavior of the various mixtures, seven different mixture types can be identified (and are listed in Table 44). In general, for binary mixture types A D, both the particle size and density ratios are equal to or greater than one ( dr r is not true for mixture types E G, whic h have one of the ratios less than one with the other greater than one. Mixture Types Each mixture type is described using a typical example of that particular type. The rest of the figures are presented in Appendix C where figures (a) and (b) are the pres sure drop profile of the fluidization and defluidization cycle for completely segregated and perfectly mixed states, respectively. Figure (c) represents the weight fraction of jetsam versus dimensionless height at different velocities and figure (d) shows how the SI varies as a function of velocity PAGE 130 130 Type A mixtures Mixtures with very large particle size ratio (dr > 4.5; Ur > 8) Type A mixtures fluidize at two distinct points when fluidized from a completely segregated state (Figure 42a). As the gas velocit y is increases, the entire segregated bed remains in a fixed state and the pressure drop linearly increases. Eventually, the flotsam becomes fluidized (point A in Figure 42a is the first point of fluidization), but the jetsam remains in a fixed bed state. At this point, the pressure drop curve follows a linear profile (as the velocity increases), but with a different slope. This change in slope is due to the partial fluidization of the bed. It is important to note that the velocity required to fluidize the flotsam is slightly greater than the minimum fluidization velocity of the flotsam alone. As the velocity is further increased, a point is reached at which both the jetsam and flotsam are fluidized (point B in Figure 42a, which is the second point of flui dization). The pressure drop across the entire bed remains constant for larger velocities. The velocity required to fluidize both the jetsam and flotsam is slightly larger than the minimum fluidization velocity of the jetsam alone. As the gas velocity is r educed, the bed height decreases and the jetsam settle at the bottom of the column (point C in Figure 42a is the first point of defluidization). The velocity at which point C occurs is generally smaller than both the velocity at point B and the minimum fl uidization velocity of the jetsam. Upon further reduction of the gas velocity, the entire bed eventually settles (point D in Figure 42a is the second point of defluidization). The velocity at this point is greater than the minimum fluidization velocity of the flotsam. Points A and D occur at similar velocities. When Type A mixtures are fluidized from an initial uniformly mixed state, there are multiple peaks observed in the pressure drop profile (oval E in Figure 42b highlights PAGE 131 131 this peaked behavior). Thes e peaks are a characteristic feature of mixtures with large particle size disparity and can be explained by the phenomenon of entrapment. As the gas velocity increases, the pressure drop across the bed increases linearly. At a certain velocity, only the fl otsam in the topmost layer of the mixture which is entrapped by the jetsam, gains sufficient momentum to fluidize. The jetsam falls back to the bottom of the column, entrapping the remainder of the flotsam in the lower layers (Figure 43). This phenomenon is associated with the first peak in the pressure drop. Beyond the first peak, the pressure drop again increases linearly with velocity until a second layer of flotsam escapes. In this manner, as additional layers of the flotsam escape, multiple peaks in t he pressure drop profile are observed. Eventually, a velocity is reached at which the entire amount jetsam fluidizes. At this point, the pressure drop across the bed stabilizes and the entire bed is completely fluidized as all of the flotsam is no longer e ntrapped. The number of peaks observed in the pressure drop profile corresponds to the number of escaping flotsam layers. The appearance of each peak corresponds very well with the visual observation of an escaping flotsam layer. Experiments were performed to determine that fewer peaks occur when the rate of increase of the gas velocity is higher. The magnitude of the pressure drop associated with each peak corresponds to the amount of flotsam escaping and the extent of each layer. The magnitude of the pres sure drop of the successive peaks reduces as the velocity increases, implying that each successive escaping layer of flotsam decreases. The nature of the defluidization curve for the mixture from an initially segregated or mixed state is the same (points F and G in Figure 42b are the first and second points of defluidization, respectively). In fact, the defluidization profile is independent of the initial PAGE 132 132 segregation profile (points C and D in Figure 42a coincide with F and G in Figure 42b). Hence, the defluidization curve is the most reproducible curve for a Type A mixture. Type A mixtures are not only characterized by 2point fluidization, but also by a minimum observed in the segregation index profile (Figure 42c). As the velocity increases beyond the complete fluidization velocity (points C and F in Figures 42a and 4 2b, respectively), the bed begins to mix and the SI decreases. However, as the gas velocity increases further, the bed expansion due to the flotsam is greater than that for the jetsam and the bed begins to segregate with a corresponding increase in SI Thus, for Type A mixtures, there exists an optimum velocity at which the SI is minimized. Type B mixtures Mixtures with significant level of disparity in particle size and density ( r > 3 o r 4.5 > dr > 3.3; 4.2 < Ur < 8) The pressure drop profiles for Type B mixtures are similar to those observed for Type A, where fluidization from a segregated state exhibits 2point fluidization behavior, and fluidization from an initially mixed state exhibits peaked behavior (points A D and E G in Figures 44a and 44.b, respectively, have the same definitions as the corresponding points in Figures 42a and 42b). The key difference between a Type A and Type B mixture is the behavior of the segregation index. For Type B mixtures, the segregation index decreases as the gas velocity increases (Figure 44c) rather than exhibiting a minimum. For Type B mixtures, the mixing quality improves as the gas velocity increases, although complete mixing is difficult to achieve and can be attained only at very large fluidization velocities. Type C mixtures Mixtures with intermediate level of disparity (2 < r < 3 or 2 < dr < 3.3; 2.5 < Ur < 4.2) PAGE 133 133 When Type C mixtures are fluidized from an initially segregated state, they exhibit 2 point fluidization (points A and B in Figure 45a are the first and second points of fluidization, respectively). However, when these mixtures are fluidized from a mixed state, they may either demonstrate single peak behavior (oval D in Figures 45b and 46) or single point fluidization behavior (shown in Figure 46). Single peak behavior is generally observed in either small diameter columns 1 or when the data acquisition system is sufficiently fast. It is this latter effect that is observed in the present study. When the mixtures were fluidized rapidly, fewer data points were obtained since the data acquisition rate remained the sam e and some of the key features (such as the peaks) of the pressure drop profiles were lost. Thus, it was necessary to fluidize the material slowly (the rate of velocity increase was 0.0833 cm/s2). The mixtures that have a single point or a single peak in t heir pressure drop profile show similar patterns. As an example, Figure 46 compares pressure drop profiles for the same mixture (75G23125G116) obtained in two different fluidized bed systems. One of the profiles is from the fluidized bed system used in t he present study (column diameter is 1.6 cm), which has a very high data acquisition rate, and the other profile was obtained using the Joseph et al. 17 fluidized bed system (column diameter is 12 cm) with a lower data acquisition rate. The pressure drop profile exhibits single pressure peak behavior, while the Joseph et al. 17 system shows single fluidization point behavior. Some combination of the effect of the data acquisition rate and column diameter is thought to have caused this difference in the behavior. PAGE 134 134 The segregation index profile (Figure 45c) shows that at low velocities the segregation index is large, but at higher velocities (approximately twice the complete fluidization velocity), the mixture mixes completely. Type D mixtures Mixtures with mi nimal disparity in particle size and density (1 < r < 2 or 1 < dr < 2; 1 < Ur < 2.5) These mixtures behave like single component particle beds and fluidize at a single point from both the initially segregated or well mixed state (Figures 47a and 47b). Type D mixtures also tend to mix easily at low vel ocities (Figure 47c). Furthermore, they do not exhibit segregation even at low velocities which are slightly above the complete fluidization velocity. Type E mixtures Mixtures with smaller, denser component as jetsam ( r > 2 and dr < 1) Type E mixtures c ontain smaller, denser particles and coarser, less dense particles such that the size difference opposes the density difference. Generally, the smaller, denser component behaves as jetsam when the density ratio is large. The pressure drop profile may show multiple peaks or a single peak depending upon the disparity level based on Ur 16, 17. Additionally, the segregation index decreases as the velocity increases. Type F mixtures Mixtures exhibiting layer inversion (1 < r < 1.5, 0.3 < dr < 1) Mixtures with smaller, denser particles and coarser, less dense particles having a low density ratio and a low to intermediate size ratio may exhibit the phenomenon of layer inversion. At lower velocities, the smaller, denser component behaves as jetsam. PAGE 135 135 However, at higher velocities, the coarser, lighter component behaves as jetsam. For these mixtures, pressure drop data and segregation index are not readily available in the literature and were not examined here. Type G mixtures M ixtures with Coarser, lighter component as jetsam (1 < r < 2, dr < 0.3) For mixtures with smaller, denser particles and coarser, less dense particles having a low density ratio and a large size ratio, the coarser, less dense components behave as jetsam. As an example, mixture number 15 from Table 44 exhibits this kind of behavior. When mixture number 15 was fluidized from an initially segregated state, it had two points of fluidization, and when it was fluidized from an initially mixed state, it showed multiple peaks. The segregation index of mixture number 15 reduced as the operating velocity was increased. Pressure drop and segregation index figures for all the binary mixture experiments performed in this study are provided in Appendix C Classification Diagram Figure 48 summarizes all of the data presented in Table 44 in a more concise manner. A loglog plot is used to due to the significant amount of available data for mixtures of smaller, denser particles and coarser, less dense particles ( r > 1 and 0.1 < dr < 1). In addition, there have been many exper iments performed for purely size segregating mixtures or density segregating mixtures. Hence, there are many data points along the x axis and y axis. The various mixture types can be classified via a plot of particle density ratio (y axis) versus particle size ratio (x axis). The boundary lines give an approximate range of the particle properties associated with each mixture type. The segregation index minimum phenomenon is not present for mixtures with large density difference ( r = PAGE 136 136 7.4). Since most practical cases involve density ratios less than this value, the boundary line for Type A mixtures is drawn as a straight line parallel to the y axis. As mentioned previously, the minimum fluidization velocity ratio is a useful measur e characterizing the level of particle disparity in the mixture. Figure 49 shows how the ratio of the Archimedes number for the jetsam to flotsam, 3 jetsamsjetsamg r flotsamsflotsamgd Ar d ( 4 8 ) which combines the effects of both the particle size and density ratios, is related to the minimum fluidization velocity ratio. For some mixtures, such as mixture numbers 15 and 34 in Table 44, Arr is less than one while Ur is greater one. In order to keep the data for all mixtures on a single plot, Arr is redefined as, 3 sflotsamg flotsam r jetsamsjetsamgd Ar d ( 4 9 ) in Figure 49 so that Arr is always greater than one. Also, for mixture numbers 39 to 43 in Table 44, both Arr and Ur are less than one. In these cases, both Arr and Ur are redefined so that both ratios are always greater than one. The Archimedes number is given by equation 4 9 while the minimum fluidization velocity ratio is min min flotsam r jetsamU U U ( 4 10) A simple correlation is developed between the two parameter s ( Arr and Ur), 0.491.02rrUAr ( 4 11) which is shown in Figure 49. If Ur > 8, then t he disparity level is extremely high (mixtures 1 to 8 in Table 44) and if 4.2 < Ur < 8, then there is a high level of disparity PAGE 137 137 (mixtures 9 18 in Table 44). Further, mixtures having Ur between 2.5 and 4.2 have intermediate disparity level (mixtures 19 to 26 in Table 44). And finally, if Ur varies from 1 to 2.5, there is a low level of disparity (mixtures 27 43 in Table 44). Hence the correlation given in equation 4 11 allows one to predict Ur, based on Arr, and thus give an indication of the level of mixture disparity. Figure 49 is a good fit for most mixtures except for numbers 5, 11, 14, and 18 in Table 44. Mixtures 5 and 14 are taken from the data of Nienow et al. 63 and Nienow et al. 64, re spectively, who reported values for minimum fluidization velocity for glass that were smaller than predictions from various correlations and previously reported and the present studys measurements. The data for m ixtures 11 and 18 are from the present study and have a higher minimum fluidization velocity for steel particles than expected from standard correlations like the Wen and Yu 39 correlation. Friction from the walls acting on the steel particles may be influencing these measurements. The friction coefficient between steel and acrylic (the column material) and the density of steel are both large which results in a la rge wall influence (Chapter 5). The minimum fluidization velocities of the other mixtures materials (glass and polystyrene) in the present study are negligibly affected by wall effects. For example, consider mixture number 22 (Table 44, i.e. 75G23125G116). This mixture has been studied in the present column ( D = 1.6 cm) as well as in a larger col umn, D = 12 cm 17. The Ur obtained for this mixture is 3.15 ( min min6 1.9jetsam r flotsamU U U ) for the present column and 3.11 ( min min5.6 1.8jetsam r flotsamU U U ) for the 12 cm diameter column. The difference in Ur, between both columns is ~1% and, assuming that Ur is the principal PAGE 138 138 parameter required to characterize different mixtures, wall effects on segregation patterns are negligible. Experiments were performed to compare segregation patterns from the smaller column ( D = 1.6 cm) to the segregation patterns from a wider column ( D = 12 cm 17) for mixture number 22. The segregation patterns for both columns were similar, validating that the wall has a minimal influence on mixture number 22. Experiments Validating Minimal Wall Effects It will be shown in Chapter 5 that it is possible for the minimum fluidization velocity of a single component system, Umf, in a smaller column to be larger than in a wider column. The increase in Umf depends upon the ratio of the particle size to column diameter, the ratio of the bed height t o column diameter, the friction coefficient between the particles and the confining wall, and the Archimedes number. For most of the particles studied, the increase in Umf was less than 14% with a bed height of 4 cm (Cha pter 5). Additionally, this small in crease in Umf for the individual components results in a small increase in Ur when the ratio of the minimum fluidization velocities of the individual components (the primary parameter used to characterize mixtures) is calculated. In addition, the segregati on profiles are only minimally altered due to the wall effect. As an example, Figure 410 compares the segregation profiles of the same mixture in two columns the column used in the present study ( D = 1.6 cm) and a larger column ( D = 12 cm 17). In this f igure, the operating velocity is scaled with the complete fluidization velocity, Uc, of the mixture and is denoted by V*. Generally, segregation profiles are given as weight fraction (or volume fraction) of either jetsam or flotsam versus dimensionless hei ght (instantaneous height, h divided by the overall bed height). However, in this case, to compare data from two different columns, weight PAGE 139 139 fraction is plotted against h/D and the zero of the dimensionless height, is defined to be at the center of the col umn. Although the segregation experiments were performed at different dimensionless operating velocities, the segregation profiles in the two columns match reasonably well. Similar results were also observed for mixture number 43 in Table 44 (30G11670P275). Effect of column H/D on Binary Segregation in Fluidized Bed Experiments investigating the effect of the ratio fixed bed height to the column diameter ( H/D) are also performed in this study. As the H/D increases, the complete fluidization velocity of the mixture also increases (especially in narrow columns). Hence to study the mixtures on the same platform the dimensionless velocity V* (as defined in the previous s ection) is maintained constant. Segregation behavior of four different mixtures namely 25G11675G231 at V* =2.3, 50G13850G328 at V* =2.5, 50G13850G275 at V* =3 and 70G11630P275 at V* =3.1, are presented in Figures 411, 412, 413 and 414 respectively. The four mixtures are from Type C category from a disparity point of view. Most of these exper iments were carried out in a column with D=1.6 cm except for the curves for mixtures 25G11675G231 at H/D=2 and 70G11630P275 at H/D=1.2 which are from Joseph et al. 17, D=12 cm). Again in these figures too the segregation profiles are plot against h/D ins tead of the conventional h/H and the zero is considered at the center of the fixed bed. This is done so that data from the two columns ( D= 1.6 cm and D=12 cm for Joseph et al. 2007) can be plotted together. From Figures 411, 412, 413 and 414 it is observed that even at very high fluidization velocities, as the H/D ratio increases, segregated tail develop at the two ends of the column. At the center however, the two components remain completely well mixed. Thus, for an infinitely tall column which is bei ng fluidized PAGE 140 140 by a very high velocity, there will be three zones, the center zone where it will be perfectly mixed and two end zones, one with a high jetsam concentration in the bottom end zone and the other with a high flotsam concentration in the top end zone. The length of the center perfectly mixed zone will depend on the level of disparity (i.e. dr, r and Ur) and the operating velocity. Summary Chapter 4 presents a new classification scheme for the pressure drop profiles, and segregation behavior of binary fluidized mixtures ( dr < 7 and r < 8) Seven mixture types are proposed. This classification s cheme is based on the particle size and density ratio of the two components and incorporates new data as well as previously published data exhibiting a wide range of fluidization behavior T he effect of H/D on segregation is also studied herein. It is obse rved that even mixtures with fairly low level of disparity being fluidized at high velocities if the H/D ratio of the column is indefinitely increased, segregated tails start appearing at the ends of the column. While conducting H/D experiments it was rea lized that as the height of the fixed bed was increased or the column diameter was reduced the minimum fluidization velocities of monosized particles were enhanced. Although the qualitative results presented in this C hapter are hardly influenced by the ele vation in minimum fluidization velocity due to the wall, there may be adverse effects on parameters like heat and mass transfer. In Chapter 5 the wall effects are quantified and a new semi correlation is proposed for the predicting the enhancement in minim um fluidization velocities, which include the wall effects. . PAGE 141 141 Table 41. Summary of published work (details in Table 44) Reference Type of Mixture Pressure Drop Profile Parameters studied Size difference Density difference Larger Denser Smaller Denser as jetsam Smaller Denser as flotsam Layer Inversion r=1 ,dr>1 r>1, dr=1 r>1, dr>1 r>1, dr<1 r<1, dr>1 r>1 ,dr<1 or Joseph et al. 17 yes yes no yes no no yes V, Cmix, H, SI Formisani et al. 16 yes yes no yes yes no yes V, Cmix, SI Formisani et al. 67 yes no no no no no yes V, Cmix, SI Marzocchella et al. 18 yes no no no no no yes V, Cmix, SI Hullin et al. 23 yes no no no no no yes V, Cmix, SI Hoffmann et al. 65 yes no no yes no no no V, SI Nienow et al. 63 yes no yes yes no no no V, Cmix, SI Nienow et al. 64 ye s no yes no no no no V, SI Olivieri et al. 66 no yes no yes no no no V, Cmix, SI Garcia et al. 22 no yes no no no no no V, Cmix, SI Rasul et al.19 no no no yes no yes no b Present study yes yes no yes no no yes V, H, SI PAGE 142 142 Table 42. Experimental m aterial and properties present study Material Diameter Density (kg/m3) Sphericity Umin (cm/s) Standard Deviation (cm/s) Notation Glass 75 89 2500 0.9 1.5 0.2 G083 104 125 2500 0.9 1.9 0.05 G116 125 152 2500 0.9 2.7 0.3 G138 152 178 2500 0.9 3 .6 0.3 G165 178 211 2500 0.9 4.6 0.3 G195 211 251 2500 0.9 6 0.6 G231 251 297 2500 0.9 8 0.2 G275 297 354 2500 0.9 11 0.3 G328 354 422 2500 0.9 13 0.2 G385 422 500 2500 0.9 19 0.3 G460 500 600 2500 0.9 25 0.8 G550 Polystyrene 251 297 1250 0. 9 4 0.3 P275 297 354 1250 0.9 7 0.5 P328 Steel 297 354 7800 0.85 46 2 S328 Table 43. Experimental mixtures and their properties present study Type Mixtures Size Ratio Density Ratio Ar no. ratio Umin Ratio Size Segregation 50G550 50G083 6.6 1 29 1 16.7 50G462 50G083 5.6 1 172 12.7 50G550 50G116 4.7 1 107 13.2 50G385 50G083 4.6 1 99.8 8.7 50G462 50G116 4.0 1 63.2 10.0 50G328 50G083 4.0 1 61.7 7.3 50G275 50G083 3.3 1 36.4 5.3 50G231 50G083 2.8 1 21.6 4.0 50G195 50G083 2.3 1 13.0 3.1 50G165 50G083 2.0 1 7.86 2.4 50G138 50G083 1.7 1 4.60 1.8 50G116 50G083 1.4 1 2.73 1.3 Density Segregation 13S328 87P328 1 7.40 7.40 6.6 75S328 25G328 1 3.10 3.10 4.2 Size and density against each other 70G11630P275 0.42 2.30 0.170 0.5 PAGE 143 143 Figure 41. Schematic of the experiment: (A) air inlet, (B) mass flow controller and pressure transducer, (C) computer, (D) column diameter, (E) air flowing into the setup, (F) pressure signal, (G) diffuser, (H) bed height, (I) exiting air, (J) air filter, (K) sliding disc assembly, (L) slice of the bed, (M) carousel arrangement, (N) sampling containers, (O) axis of rotation of the sliding disc assembly, and (P) side handle to help rotate the sliding disc assembly about axis O. PAGE 144 144 Table 44. Details of mixtures and their properties published work and present study Sr no. Type Ur Mixtures Size Ratio Density Ratio Ar # Ratio Reference J/F Mixtures Jetsam Flotsam J/F J/F J/F 1 A 16.67 50G55050G083 Glass Glass 6.63 1.00 291 Present Study 2 A 13.16 50G55050G116 Glass Glass 4.74 1.00 107 Present Study 3 A 12.94 50G50050Si125 Glass Silica Sand 4.00 0.98 63 Marzocchella et al. 18 4 A 12.67 50G46250G083 Glass Glass 5.57 1.00 172 Present Study 5 B 11.61 10S32490G165 Steel Glass 1.96 2.52 19.10 Nienow et al. 64 6 A 10.00 50G46250G116 Glass Glass 3.98 1.00 63.18 Present Study 7 A 8.67 50G38550G083 Glass Glass 4.64 1.00 99.80 Present Study 8 A 7.98 50G61250G154 Glass Glass 3.97 1.00 62.76 Formisani et al. 16 9 B 7.50 10G55090G165 Glass Glass 3.33 1.00 37.04 Nienow et al. 63, 64 10 B 7.33 50G32850G083 Glass Glass 3.95 1.00 61.71 Present Study 11 B 6.57 13S32887P328 Steel Polystyrene 1.00 7.43 7.43 Present Study 12 B 6.00 50G49950G172 Glass Glass 2.90 1.00 24.42 Formisani et al. 67 13 B 5.33 50G27550G083 Glass Glass 3.31 1.00 36.37 Present Study 14 B 5.30 10S27390G231 Steel Glass 1.18 2.52 4.16 Nienow et al. 64 15 G 5.00 50MS62450G154 Mol Sieves Glass 0.25 1.70 0.03 Formisani et al. 16 16 B 4.57 10Cp46190Q273 Copper Powder Quartz 1.69 3.34 16.10 Nienow et al. 63 17 B 4.35 10S39090S138 Steel Steel 2.83 1.00 22.57 Nienow et al. 64 18 B 4.18 75S32825G328 Steel Glass 1.00 3.12 3.12 Present Study 19 C 4.00 50G23150G083 Glass Glass 2.78 1.00 21.56 Present Study 20 C 3.56 50G56550G285 Glass Glass 1.98 1.00 7.79 Hoffmann et al. 65 21 C 3.33 55G549045G1590 Glass Glass 3.45 1.00 41.16 Huilin et al. 23 22 C 3.11 25G23175G116 Glass Glass 1.99 1.00 7.90 Joseph et al. 17 23 C 3.07 50G19550G083 Glass Glass 2.35 1.00 12.97 Present Study 24 C 3 .00 50G499 50G271 Glass Glass 1.84 1.00 6.24 Formisani et al. 16, 67 PAGE 145 145 Table 4 4. continued Sr no. Type Ur Mixtures Size Ratio Density Ratio Ar # Ratio Reference J/F Mixtures Jetsam Flotsam J/F J/F J/F 25 C 2.53 50S43950G428 Steel Glass 1.03 3.06 3.31 Formisani et al. 16 26 C 2.50 55G426045G2300 Glass Glass 1.85 1.00 6.35 Huilin et al. 23 27 D 2.40 50G16550G083 Glass Glass 1.99 1.00 7.86 Present Study 28 D 2.11 50G56550G365 Glass Glass 1.55 1.00 3.71 Hoffmann et al. 65 29 D 1.93 75G23125P231 Glass Polystyrene 1.00 2.33 2.33 Joseph et al. 17 30 D 1.86 46P32854P231 Polystyrene Polystyrene 1.42 1.00 2.86 Joseph et al. 17 31 D 1.80 50G13850G083 Glass Glass 1.66 1.00 4.60 Present Study 32 D 1.73 SS500 PP500 Silica sand Polypropylene 1.00 2.89 2.89 Olivieri et al. 66 33 D 1.67 50G375050PE3750 Glass Polyethylene 1.00 2.39 2.39 Garcia et al. 22 34 E 1.39 50S43950MS800 Steel Mol Sieves 0.55 5.21 0.86 Formisani et al. 16 35 E 1.35 50G59350MS624 Glass Mol Sieves 0.95 1.70 1.46 Formisani et al 16 36 D 1.33 50G375050A3750 Glass Alumina 1.00 1.57 1.57 Garcia et al. 22 37 D 1.27 50G11650G083 Glass Glass 1.40 1.00 2.73 Present Study 38 D 1.25 50A375050PE3750 Alumina Polyethylene 1.00 1.52 1.52 Garcia et al. 22 39 E 0.96 10B27390G461 Bronze Glass 0.59 2.89 0.601 Nienow et al. 63 40 E 0.69 80SG37520SS125 Silica Sand Silica Gel 0.33 4.33 0.160 Olivieri et al. 66 41 E 0.67 50B23550G565 Bronze Glass 0.42 3.49 0.251 Hoffmann et al. 65 42 F 0.52 FCCPumice FCC Pumice 0.65 1.12 0.311 Rasul et al. 19 43 E 0.45 70G11630P275 Glass Polystyrene 0.42 2.33 0.175 Joseph et al. 17 44 E FCCBagasse FCC Bagasse 0.32 2.89 0.095 Rasul et al. 19 45 E BagasseP2000 Bagasse Polystyrene 0.10 2.46 0.002 Rasul et al. 19 46 E PVC Bagasse PVC Bagasse 0.32 1 .93 0.063 Rasul et al. 19 47 F Cenolyte Bagasse Cenolyte Bagasse 0.32 1.40 0.046 Rasul et al. 19 PAGE 146 146 Figure 42. Typical pressure drop profiles and segregation index behavior for a Type A mixture 50G38550G083. (a) P ressure drop profile, initially segr egated state. (b) Pressure drop profile, initially mixed state. (c) Segregation index. PAGE 147 147 Figure 43. Process of fluidization for a homogeneous binary mixture which fluidizes at two points PAGE 148 148 Figure 44. Typical pressure drop profiles and segregation index behavior for a Type B mixture 75S32825G328. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation index. PAGE 149 149 Figure 45. Typical pressure drop profiles and segregation index behavior for a Type C mixture 50G19550G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation index. PAGE 150 150 Figure 46. Pressure drop profiles for the mixture 75G23125G116 from an initially mixed s tate. Comparison of the present work (column diameter is 1.6 cm) with the work of Joseph et al. 17 (column diameter is 12 cm). PAGE 151 151 Figure 47. Typical pressure drop profiles and segregation index behavior for a Type D mixture 50G16550G083. (a) Pressure dr op profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation index PAGE 152 152 Fig ure 48. M ixture type diagram. Type A: Very large particle size ratio. Type B: Significant level of disparity in particle size and densi ty. Type C: Intermediate level of disparity. Type D: Minimal disparity in particle size and density. Type E: Smaller, denser component as jetsam. Type F: Mixtures exhibiting layer inversion. Type G: Coarser, lighter component as jetsam. PAGE 153 153 Figure 49. Cor relation for the minimum fluidization velocity ratio and the Archimedes number ratio. Figure 410. Segregation profiles for the mixture 75G23125G116. (a) Present work ( D = 1.6 cm). (b) Joseph et al. 17 ( D = 12 cm) PAGE 154 154 Figure 411. Segregation profiles of 25G11675G231at V* =2.3 Figur e 412. Segregation profiles of 50G13850G328 at V* =2.5 PAGE 155 155 Figure 413. Segregation profiles of 50G13850G275 at V* =3 Figure 414. Segregation profiles of 70G11630P275 at V* =3.1 Equation Chapter 5 Section 1 PAGE 156 156 CHAPTER 5 EFFECT OF COLUMN DIAMETER AND BED HEIGHT ON MINIMUM FLUIDIZAT ION VELOCITY Background MFBs are gaining inte rest as gas distribution issues are minimized and due to the reduction in the quantity of solids required. Nonetheless there is a pitfall that comes along with the MFBs, i.e. wall effects may influence the observed results. As stated in Chapter 1 there ar e several correlations that exist for calculating the minimum fluidization velocity 2639. None of the correlations include the effect of the bed height or diameter of the column, factors which are likely to be of interest when MFBs are considered. The Wen and Yu 39 correlation, one of the more commonly used correlation for calculating the minimum fluidization velocity and is given as, 2 3Re Re 24.5 1650 0pmf pmfAr ( 5 1 ) where () Regmf pmf gdU, ( 5 2 ) 3 2()()sgg gdg Ar ( 5 3 ) where Rep mf is Reynolds number at minimum fluidiz ation velocity, Ar is the Archimedes number, is the sphericity, Umf is the minimum fluidization velocit y. Neither the Wen and Yu 39 nor any of the other correlations include the effect of the fixed bed height, H or diameter of the column, D factors which are likely to be of interest when MFBs are considered. PAGE 157 157 Two recent investigations, however, have examined the influence of column diameter and bed height. Di Felice and Gibilaro 70 described a method for predicting the pressure drop across a particle bed, taking into account the effect of column diameter. They considered the column to be comprised of two sections: an inner core where the voidage remains nearly constant, and an outer annular section where the voidage varies due to the presence of the wall. Since there is a difference in the voidage ov er the cross section of the column, the velocity also varies across the columns cross section. This fact was used to develop a modified Erguns equation, but this effect is only observed at very small column diameter ( D ) to particle diameter (d ) ratios, i .e. D / d < 15. Delebarre 71 assumed that the fluidizing gas density is a function of the instantaneous bed height which in turn influences the minimum fluidizing velocity. This effect only occurs if the column is very tall. For micro fluidized beds, where t he columns are short, this effect does not play much of a role. As will be shown in the following section, experiments in which the bed height and column diameter are varied show that both parameters influence the beds minimum fluidization velocity. Neith er Di Felice and Gibilaro 70, n or Delebarre 71 models are able to capture the experimental results. The remainder of this paper describes the theory and experiments in which a new mechanism is proposed for accounting for bed height and column wall effects. Experimental Setup Two fluidization segregation units (Jenike and Johanson, Fluidization Material Sparing Segregation Tester and Fluidization Segregation Tester), were used in the experiments. The first tester has a column diameter of 1.6 cm and a height of 9.5 cm PAGE 158 158 and the latter tester has a 2.4 cm diameter column and an 18.5 cm height. Details concerning the operation of the testers are given in Hedden et al. 68 and ASTM D 6941 69. A schematic of the experimental set up is shown in Figure 41. The experiments were carried out in both columns. A sintered metal plate with an average pore diameter from the bottom, and its flow rate is controlled by a mass flow controller The pressure drop across the entire setup was measured with a pressure transducer. The instantaneous pressure drop and velocity data were recorded on a computer. Experiments were performed with glass and polystyrene particles of different sizes, ranging particles were carefully sieved, then dried in an oven for twelve hours and finally subject to an antistatic bar (Takk industries) in order to eliminate accumulated electrostatic charge. The particles were stored in a desiccator so that cohesive forces due to moisture were minimized. To further reduce electrostatic charges that develop during fluidization, a very small amount (approximately 5 mg) of antistatic powder (Larostat HTS 905 S manufactured by BASF Corporation) was mixed with the particles, preceding the experiments. Table 51 summarizes the experimental materials and their properties. The notation used to represent the particle type is: G for glass and P for polystyrene, followed by the average particle size in microns. The values for the dynamic friction coefficients of glass and polystyrene on acrylic (the column wall material) are not readily available and so values of friction coefficients, for glass on PAGE 159 159 glass and polystyrene on polystyrene, obtained from Persson and Tosatti 72 were used in the model calculations. Prior to running an experiment, air was passed through the empty column to get the background pressure drop due to the column, diffuser, and the filter sections. N ext, the particulate material was weighed and loaded into the column from the top. The height, H to which the column was filled was recorded. The antistatic powder was mixed with the particles. The air velocity was slowly increased beyond the point of flu idization, and then decreased to zero to get the entire pressure drop profile. The point of intersection of the pressure drop line for the fixed bed and the horizontal line for the fully fluidized bed is typically defined as the minimum fluidization veloci ty. However, just before fluidization, the pressure drop across the bed overshoots the expected value and then decreases to a constant value. Such behavior is expected in columns with a small diameter. This overshoot was not observed during defluidization (Figure 51). Thus, the defluidization pressure drop curve was used to determine the point of fluidization 25, 73. A minimal of three cycles of fluidization and defluidization were carried out to determine an average minimum fluidization velocity. Each experiment was repeated twice and the procedure was repeated for all particle types. The maximum fluctuation observed in the value of Umf was ~ 7% (for G550 in D = 1.6 cm and H = 4.5 cm; Umf = 24.5 cm/s 0.9). The fairly consistent value of Umf observed over the three cycles and two repetitions show that the trends in Umf are not due to poor gas distribution but rather due to the wall effects. Additionally, gas channeling and maldistribution were not observed through the transparent columns during any of the experiments. PAGE 160 160 Figures 52 to 55 present the experimental data for the minimum fluidization velocity ( Umf), for the glass and polystyrene particles, obtained from both the columns. In these figures, the Reynolds number at Umf, is plotted against the aspect ratio of the bed ( H/D). In addition to the experimental data, the plots also include theoretical curves which will be described in the following section. The error bars represent the difference in the reading between two repetiti ons and not amongst different cycles. The experimental data (Figures 52 5 5) show that minimum fluidization velocity monotonically increases as the height of the bed is increased. Further, wall effects on the minimum fluidization velocity are prominent when column diameter, D = 1.6 cm (Figure 52 and Figure 54). While, wall effects on minimum fluidization velocity are much less pronounced for column diameter, D = 2.4 cm (Figure 5 3 and Figure 55). Although the model proposed by Di Felice and Gibilaro 70 predicts this trend, their model is limited to column to particle diameter ratios of less than 15. In the experiments performed here, this ratio can be as large as 150, well outside the range of validity in Di Felice and Gibilaro 70 model. As the height of the bed increases, the minimum fluidization velocity al so increases. Delebarre 71 model accounts for changes in the fluidizing gas density over the bed height and predicts at most a minimum fluidization velocity increase of approximately 0.5%, while the observed minimum fluidization velocity increase is greater than 20% in certain cases for the experiments examined here. Thus, this model does not predict the measured data. Theory Since existing correlations do not include the influence of column diameter or bed height, and because the recent models by Di Felice and Gibilaro 70 and Delebarre 71 do PAGE 161 161 not predict the measured data, a new model is developed in this section in order to account for the observed trends. As will be presented, a key component of the model is the stress between the column wall and the bed. The pressure drop, for a fixed bed of height, H is given by the semi empirical correlation of Er gun 1 as 2 EE P bVaV H ( 5 4 ) where 3(1) 1.75f Ea d and ( 5 5 ) 2 322(1) 150f Eb d. ( 5 6 ) The param eter V is the superficial fluid velocity and is the bed voidage. The bed void fraction and fluid density are assumed to remain constant throughout the bed. The pressure drop expression ( equation 5 4 ) can be used to predict the minimum fluidization velocity of a monosized material by balancing the pressure force acting on the bed by the effective bed weight. 2 2 33(1) (1) [1.75][150][(1)]pmf pmfRe Re Ar ( 5 7 ) Note that in the previous analysis, the stress between the wall and the bed has not been included. In the case of a fixed bed with no interstitial fluid, a fraction of the bed weight is supported by the walls of the column. This phenomenon is known as Janssens wall effect. To determine the wall force acting on the bed, consider a thin horizonta l slice of the particle bed of thickness dz (Figure 56). The vertical stress, v, acting on the top PAGE 162 162 face and the weight of the element are balanced by the vertical stress acting on the bottom face and the frictional forces applied by the wall. The resulti ng differential equation is, v vd cc dz12 ( 5 8 ) where 14Jtank c D and ( 5 9 ) 2(1)s cg ( 5 10) The parameter is the friction angle between the wall and particles. A horizontal frictional stress, H, arises from the vertical stress, v, produci ng a wall friction force of Htan 74. It is assumed here that this horizontal stress is directly proportional to the vertical stress with kJ as the proporti onality constant 74. HJvk ( 5 11) Solving equation 5 8 subject to the boundar y condition that the stress at the top of the bed is zero gives, cz vc e c12 1(1) ( 5 12) The stress increases nearly linearly near the free surface, but asymptotes to a constant value deeper into the bed. These wall induced vertical stresses will affect the minimum fluidization velocity in narrow columns. Now consider the situation when flui d is flowing through a fixed bed of solids with the same coordinate system as shown in Figure 56. At equilibrium, the force due to stress on the top face of the element, FT, and the weight of the element, FW, will be PAGE 163 163 balanced by the force due to the stress on the bottom face, FB, the wall friction, FF (Janssens effect), the drag force on the solids, Fd, and the buoyancy force, FBu. Thus, a force balance on the solids yields: ,BTFdBu WFFFFFF ( 5 13) which, when expanded is 22 22() 44 (1()). 44vvH v sgDdtanDdzDdp D gDdz ( 5 14) In the above expression Fw and FBu have been combined to give the last term on the right hand side of equation 5 14. Simplifying and rearranging the previous equation gives, 4 1()vH s fdtan dp g dzD dz ( 5 15) If HJvk is assumed, as in equation 5 11 4 1()f Jv v stank d dp g dzD dz ( 5 16) which may be written as 12 v vd cc dz ( 5 17) where 14Jtank c D and ( 5 18) 21()f sdp cg dz ( 5 19) PAGE 164 164 In equation 5 19 there are two additional terms as compared to equation 5 10 : the drag term, represented by dp/dz and the buoyancy term. The solution to equation 5 19 is equivalent to that of equation 5 8 cz vc e c12 1(1) ( 5 20) but with different c1 and c2. Figure 57 presents equation 5 20 in graphical form. When the velocity is zero ( i.e. there is no gas flowing through the column), then equation 5 20 simplifies to Janssens equation. As the velocity increases, the vertical contact stresses in the bed are smaller because the drag on the partic les supports a fraction of the beds weight. At and above the point of fluidization, the stresses in the column are zero (as contact does not exist amongst the particles), and thus v/dz is also zero. This situation is only possible if c2=0, 21()0f sdp cg dz ( 5 21) In equation 5 21, if Erguns pressure drop expression is substituted ( equation 5 4 ), the resultant expression will be the same as equation 5 7 and the wall effects will not be included in the analysis. It is now assum ed that the horizontal stress is not only a function of the downward vertical stress, but also a function the upward drag forces, caused by the gas flowing through the column. Since Erguns pressure drop equation is used to calculate the pressure drop, it is assumed that the structure of these new terms on which the horizontal stress depends is similar in form to those given by Erguns equation, but the PAGE 165 165 velocity terms are being scaled differently. These scales are chosen as they gave a favorable fit to the experimental data. 2 12 J f f HvkkV H D kV Dd H ( 5 22) where k1 and k2 are constants. Substituting equation 5 22 back into equation 5 15, yields: 1 2 24 4() 4()() 1().f f Jv v f stank d dzD H H tank tank dp dD D VVg D D dz ( 5 23) Comparing equation 5 23 to equation 5 16, there are two additional terms, proportional to V2 and V appearing on the right hand side of equation 5 23. These terms have come about due to the horizontal stress from the wall and the flow of the fluid through the bed. Substituting Erguns equation for pressure drop ( equation 5 4 ) into equation 5 23: 2 2 24 4() 4()() 1().Jv v f E Es f ftank d dzD H H tank tank dD D aV bV g DD ( 5 24) This equation is again of the form: v vd cc dz12 ( 5 25) where: 14Jtank c D and ( 5 26) PAGE 166 166 2 2 2 1 4()() 4()1().f E f E sH tank D ca V D H tank dD bVg D ( 5 27) Setting c2=0, the condition which is necessary for fluidization, and integrating over the length of the column: 1 2 24() 4()() 1()0f g EE f sH H tank tank dD D a Vb V g DD ( 5 28) where k1 and k2 are lumped parameters which also include integration coefficients. Substituting for aE and bE, which are Erguns constants, yields; 2 2 3 2 1 32(1) 1.754 () (1) 1504 ()1(). ()f s f fdH tank DDd V dH tank g V DDd ( 5 29) Multiplying throughout by the factor 3 2) (f fd in order to make th e equation dimensionless gives: '2 2 3 2 1 3(1) 1.754 ()Re (1) 1504 ()Re1.pmf pmfdH tank DD dH tank Ar DD ( 5 30) Equation 5 30 is a quadratic equation which can be easily solved for Repmf, with k1 and k2 remaining constant. Rewriting equation 5 30 results in PAGE 167 167 '' 2'',,,Re,,,Re10pmf pmfdH dH a b Ar DD DD ( 5 31) where a and b are functions not only of and but also of ( ) and ( H/D ). Hence, both the ratio of particle size to column diameter ( ) and the bed aspect ratio ( H/D) play a role in determining Umf. Results and Discussion The values of the universal constants k1 and k2 which give the minimum error for the experimental data set presented in this study were found out to be 610 and 30.1, respectively. In Figures 52 to 5 5 and Figures 58 and 59, the values of the universal constants k1 and k2 are kept at these values for all of the particle types, bed heights, and column diameters. Figures 52 to 55 show that including the wall effects in the prediction of Umf is successful in describing the effects of height versus column diameter ( H/D) and the effect of particle size to column diameter ( d/ D ). The model only slightly under predicts the increase in Umf for very small sized particles ( d D = 1.6 cm, which is likely due electro static cohesive and adhesive forces. For a fixed column diameter, as the H/D ratio increas es, Umf increases (Figures 52 to 55). Hence, as the column becomes taller in comparison to its diameter, the wall effects are more prominent, making it more difficult to fluidize the particles. Also, for a constant particle size, as d/D ratio increases ( due to changes in column diameter), Umf increases (compare Figures 52 and 53 or Figures 54 and 55). As the column diameter decreases in comparison to the particle diameter, the ratio of wall contact surface area to the bulk volume increases, leading to a more significant wall effect. This wall effect reduces as the column increases in diameter relative to the particles, and PAGE 168 168 eventually becomes negligible as d/D becomes very small. The effect of H/D and d/D on Umf are very similar; in fact, the increase i n Umf can be characterized by the product of H/D and d/D ( equation 5 30). Additionally, as particle size increases and approaches the same order of magnitude as the column diameter ( d/D ~ O(1)), then fluidization becomes difficult. Above a certain d/D ratio, the particle mixture is impossible to fluidize. This situation is predicted by the model, as the quadratic equation yields complex roots. Similarly, equati on 5 30 breaks down for large H/D ratio predicting that at large H/D the gas may not fluidize the bed. For example, the model predicts that if glass particles of D/d ~3) with H/D = 1, then the particles will not fluidize. Similarly, the model also predicts that if these particles are loaded into a 100 cm diameter column to a bed height of about 13 m ( H/D~13) the gas may not fluidize the bed. Figure 58 compares the experimental data from Liu et al. 25 to the values predicted by the model. The data are plotted in the same dimensional form as given in the original work of Liu et al. 25. The model does not predict the increase in Umf for the smaller particles ( d sizes. In this case also, the values of k1 and k2 are maintained at 610 and 30.1 respectively. For the largest particles ( d diameter column is not that good, but this may be due to experimental error, which is not reported in Liu et al. 25. Also, the values of voidage were not clearly stated; only values of bulk density were given. In the model, the values of bulk density were kept constant wi th respect to the column diameter, but it has been observed that this value changes depending on the column diameter 73, particularly for small column diameters. PAGE 169 169 Figure 59 compares the model results to the average predictions of the Reynolds number at Umf based on a large number of existing correlations 2639. The scatter bars on the correlation line represent the deviation in minimum fluidization velocity predictions using the different correlations. The experimental data obtained when the wall effects ar e minimum (at H =2cm and D =2.4cm) are consistent with the model prediction and with the averaged correlation curve. On the other hand, only the new model is able to predict experimental data obtained when the wall effects are significant (at H =10cm, D =1.6cm ), since the existing correlations for Umf do not take wall effects into account. Summary Experiments show that wall effects influence the minimum fluidization velocity. Existing correlations fail to incorporate these wall effects The model presented in t his study attempts to include the wall influence by introducing Janssens wall effect in the force balance during fluidization. This assumption leads to a modified Erguns equation. The new model fits experimental data reasonably well for the minimum fluidization velocity over a range of particles sizes (100 m 500 m), bed heights, and column diameters ( H/D > 6) Some deficiencies in the model predictions are noted for smaller e significance of cohesive and adhesive forces at this scale. PAGE 170 170 Table 51. Experimental materials and properties Material Diameter Density (kg/m 3 ) Sphericity Coefficient of friction Notation Glass 105 125 2500 0.9 0.4 G116 Glass 210 250 2500 0.9 0.4 G231 Glass 250 300 2500 0.9 0.4 G275 Glass 350 420 2500 0.9 0.4 G385 Glass 420 500 2500 0.9 0.4 G462 Glass 500 600 2500 0.9 0.4 G550 Polystyrene 250 300 1250 0.9 0.5 P275 Polystyrene 300 354 1250 0.9 0.5 P328 Figure 51. Example of a pressure drop profile (fluidization and defluidization) using G231 particles in the 1.6 cm diameter column. The minimum fluidization velocity is also shown in the figure PAGE 171 171 Figure 52. Repmf as a function of H/D, for different glass particles in the column D =1.6 cm Figure 53. Repmf as a function of H/D, for different glass particles in the column D =2.4 cm. PAGE 172 172 Figure 54. Repmf as a function of H/D, for polystyrene particles in the column D =1.6 cm. Figure 55. Repmf as a function of H/D, for polystyrene particles in the column D =2.4 cm. PAGE 173 173 Figure 56. Sketch showing the fluidized bed coordinate system. Figure 57. A schematic showing how the vertical stress in the bed varies with bed depth at different fluid velocity V1, V2 and V3. PAGE 174 174 Figure 58. The minimum fluidization velocity as a function of the column diameter using the experimental data from Liu et al. 25. The curves are from equation 5 30. Figure 59. Comparison of the curves from equation 5 30 to the experimental data and existing correlations. Equation Chapter 6 Section 1 PAGE 175 175 CHAPTER 6 CONCLUSIONS AND RECOMM ENDATIONS In this dissertation several aspects of fluidsolid fluidized systems are studied. Models and corelations are developed and many qualitative concl usions are drawn. The most important contributions of the present study are: Development of a dilute, turbulent, gas solid flow model which incorporates an improved description for i nteractions at the level of velocity fluctuations in both phases Acquisition of detailed and nonintrusive flow data for dilutephase liquidsolid flows in the transitional (intermediate ST ) regime Development of a dilute, turbulent, liquidsolid flow model for viscous, transitional and inertia flow regimes Acquisition of axial segregation data of binary mixtures for gas fluidized beds and its qualitative analysis Development of a semi correlation to predict enhancement in minimum fluidization velocity in fluidized bed due to wall effects Dilute Turbulent FluidSolid Fluidized Flows Experimental data for fully developed profiles of pneumatically conveyed solid part icles in a vertical pipe 2 have been available in the literature for more than twenty year s. Many authors have proposed Eulerian based, dilute turbulent gas solid flow models incorporating particleparticle interactions using a two equation kmodel to describe gas phase turbulence 7, 8, 46, 9 to simulate these data. These Eulerian models have used vari ous combinations of relations for drag, solidphase stress, and fluctuating interaction terms to successfully predict for the gas solid flows mean velocities. Unfortunately these models consistently under predict the gas turbulence and granular temperature. In the present study, the work of Bolio et al 7 is advanced to include a new closure relation for the fluctuating velocity interaction. PAGE 176 176 T he proposed new model (FET model along with Sinclair and Mallo, 44 cross correlation) for thefluctuating interaction term is formulated using an analogy with heat transfer. The time scales for the FET model depend on the Stokes number ( ST ) while activation of vortex shedding depends on Rep. If ST < 100, particle drag is responsible for the energy transfer and if ST > 100, then particle collisions are responsible for energy transfer. These observations are consistent with the f indings of Gore and Crowe 10 and Hestroni 11. The proposed new fluctuating interaction term model, along with the Wen and Yu 39 drag relation, the Peirano and Leckner 51 solid stress closure which includes fluid effects, is evaluated against several benchmark experimental data sets. The new model predicts the mean velocity profiles and also the fluctuations velocity profiles of gas and solid for both small and large particles. For particles with Rep > 150 vortex shedding is included in this model. It is also observed that the fluctuating interaction terms strongly influence the magnitude of gas turbulence away from the wall. Near the wall, turbulence generation and dissipation dominate over the fluctuating interaction term. In contrast, the predicted profiles are not sensitive to the choi ce of drag model or the solidstress closure. Nevertheless, for liquidsolid flows the Peirano and Leckner 51 model for solid stress shows an improvement in the flow predictions because it takes the fluid effects into account. Also, the correlation of Syam lal and O'Brien 49 predi cts the slip velocity more accurately in the case of the liquidsolid cases. T he mean and fluct uation velocity data for 0.5 mm and 1 mm 13 and 1.5 mm (present study) particles obtained in a vertical pipe under dilute and highly turbulent conditions are benchmark data sets. The 2.32 mm particles 12, with ST > 40, PAGE 177 177 demonstrate inertiadominated flow regime. In this flow regime, particleparticle interactions dominate the flow behavior The high ST model uses the Syamlal and OBrien 49 d rag relationship which predicts the velocity slip correctly the Peirano and Leckner 51 model for solidphase stress and the FET model which predict s the shape and magnitude of the fluctuating velocities well especially at the wall The 0.5 mm particles 13, with ST < 5, display viscous dominated flow behavior wherein the particles te nd to mimic the fluid motion. H ence, the flow behavior for the liquid and solid are similar. The Chen and Wood 15 model i s used to describe the behavior for the low ST flows Finally, the 1 mm 13 and 1.5 mm particles with 5 < ST < 40, exhibit transitional behavior. There is sudden increase in the solid velocity fluctuations with increasing ST which may be due to the curved paths followed by the solids. The data obtained for so me of the cases show nonmonotonic relationships with the mass loading and operating velocity. For the intermediate ST model, the closures described for solidp hase viscosity and conductivity predict the shape of the fluctuating velocity profiles appropriately In addition, the correlation developed for the time scale for fluctuating energy interaction term predicts the flow behavior of the transitional regime fairly well. To improve the predictions in the transitional regime, granular kinetic theory must be modified to include the effects of the deviations from the straight line paths followed by the solids. Also, it is shown in the present study that models developed for gas solid flows can be extended to liquidsolid flows. The models described herein g i ve insight i n to how fluid solid flows behave under various operating conditions and ST PAGE 178 178 In industry, fluidsolid flows are often operated at higher mass loadings than the mass loadings of the experimental conditions seen in the present study. In the future efforts should be made to measure the detailed velocity profile of fluidsolid flows at higher mass loading using a nonintrusive measurement technique. Although it is difficult to measure velocities under higher mass loading conditions (due to signal atte nuation and low data rate) matching refractive indices of the fluid and solid phases may help. The data acquired can be used to validate densephase flow models Furthermore, models for prediction of drag forces on nonspherical particles are also available 75. If these models are used to close the drag relation in the Eulerian two fluid model potentially particle shape can also be included in the analysis. Further the flow loop experiments can be extended to study radial segregation of binary mixtures under dilute and turbulent flow conditions. For the segregation experiments, a PDPA system must be used as it not only measures velocities but also measures the diameter distribution. This information can be used to develop an Eulerian model to predict flui d solid flows of binary mixtures. This model will also be usef ul because real situation s deal with polydispered particle systems Fluidized Beds Segregation of Binary Mixtures and Wall Effects In this study axial segregation of binary mixtures is studied in a fluidized bed. A new classification scheme for the pressure drop profiles and segregation behavior of binary fluidized mixtures is presented and seven mixture types are proposed. This classification scheme is based on the particle size and density r atio of the two components and incorporates new data as well as previously published data exhibiting a wide range of fluidization behavior PAGE 179 179 In addition, based on the Archimedes number ratio for the mixture, the ratio of minimum fluidization velocities of t he individual components can be estimated and the level of disparity can be identified. The knowledge of the mixture type and level of disparity in advance is a significant aid when one can select the size or density ratio in order to mitigate fluidization segregation and improve process efficiency. Further, identifying the jetsam and the flotsam in case of mixtures with a difference in size and density opposing each other may also help explain deviations from regular behavioral patterns due to layer invers ion. Additional experimentation will be necessary to further refine the boundari es for the seven mixture types. If the experiments are performed in a column with a large diameter (to minimize the influence of the wall), quantitative results can be obtained instead of qualitative observations, to precisely predict the segregation pattern. T he effect of H/D on segregation is also studied in the present work It is observed that even at high fluidization velocities and with mixtures with fairly low level of di sparity, if the H/D ratio of the column is indefinitely increased, segregated tails start appearing at the ends of the column. Experiments show that wall effects influence the minimum fluidization velocity of monosized particles If the H/D or the d/D rati o is increased the minimum fluidization velocity increases. Existing correlations fail to incorporate these wall effects and, hence, there is a need for a new model that can take these effects into consideration. The model presented in this study attempts to include the wall influence by introducing Janssens wall effect in the force balance during fluidization. It is assumed that the horizontal stresses acting at the wall are not only a function of the local vertical stress, PAGE 180 180 but also are a quadratic functi on of velocity. These new terms have the same structure as that of the drag term, i.e., the pressure drop, as given by the widely accepted Ergun equation. This assumption leads to a modified Erguns equation incorporating two universal constants. The new m odel fits experimental data reasonably well for the minimum fluidization velocity over a range of particles sizes, bed heights, and column diameters. This new model should greatly facilitate scaling results from micro fluidized beds to more traditional fluid bed sizes. This model can be further improved to include of cohesive and adhesive effects due to electrostatic forces. This may help predict minimum fluidization velocity for PAGE 181 181 APPENDIX A DILUTETURBULENT F LUIDSOLID FLOW PROGRAM PhD_Data.m %User Input File %Operating Conditions Re_center_line=yes; m=0.0000001; %Experimental set up R=0.01525; g=9.81; %Fluid muf=0.000018; rhof=1.2; Re=14.2*rhof*R*2/muf; %Solid rhos=1020; d=2780*10^ 6; sphere=1; nus0=0.65; e=0.9; ew=0.9; fi=0.002; %To specify a grid say yes %To use default grid say no User_Defined_Grid=no; %Is the flow upward or downward downward=no; PhD_Model_Options.m % Select Models % 1)The Interaction terms % 1.1)Drag Term description Ding_Gidaspow_1990=yes; % 1.1.1)Definition of CD Wen_Yu=no; Hill_Koch_Ladd=no; Syamlal_OBrien=yes; % 1.2)FKS & FKF description FET=yes; %FET + Sinclair and Mallo Louge_etal_1991=no; %TVBA + Louge Sinclair_and_Mallo_1998=no; %TVBA + Sinclair and Mallo Simonin_1996=no; %TVBA + Simonin Koch_Sangani_1999=no; %TVBA + Koch & Sangani Wylie_etal_2003=no; %TVBA + Wylie Zhang_Reese_2003=no; %VBA + Zhang and Reese Chen_Wood=no; %TVBA + Chen and Wood % 1.3)Wa ke enhancement in k E equation Lun_Wake_term =yes; % 2)Solid Stress Models GKT=yes; % 2.1)Description of diffusion coefficients for Vsz and T Lun=no; Pierano_Leckner=yes; % 3)Boundary Conditions PAGE 182 182 % 3.1)Vfzd Boundary Condition Vfz_no_s lip=yes; Vfz_wall_Mfix=no; % 3.2)Vszd Boundary Condition Vsz_no_slip=no; Johnson_Jackson_Vsz=yes; Vsz_wall_Mfix=no; fiw=0*pi/180; % 3.3)Td Boundary Condition Johnson_Jackson_T=yes; Td_wall_Mfix=no; fiw=0*pi/1 80; % 3.4)kd Boundary Condition k_no_slip=yes; k_wall_Mfix=no; % 3.5)Ed Boundary Condition Bolio_equation_2_62=yes; E_wall_Mfix=no; % 4)muef model Batchelor_Green=yes; PhD_Model_Parameters.m %Set Model Parameters and Weights %Drag term parameters wVfz(1)= 10; %weight of drag term on Vfz equation wVsz(1)=10; %weight of drag term on Vsz equation %Myoung and Kasagi Parameters sigmak=1.4; sigmaE=1.3; cT1=1.4; cT2=1.8; cT3=1.2; cmu=0.09; fT1=1; %Mfix Parameters % sigmak=1; % sigmaE=1.3; % cT1=1.44; % cT2=1.92; % cT3=1.22; % cmu=0.09; % fT1=1; wk(1)= 20; %weight of Generation term on k equation wE=[10 10]; %weight of 1)Generation; 2)Dissipation terms on E equation %PTE parameters wT=[ 10 10]; %weight of 1)Generation; 2)Dissipation terms o n T equation %Interacton terms parameters wk(2)= 10; %weight of Interaction term on k equation wT(3)= 10; %weight of Interaction term on T equation PhD_Marching_Parameters.m %Set Marching Parameters %Setting Marching Parameters %For Re Marching n_Re_march =10000; PAGE 183 183 Re_jump=0.05; n_Re_iterations=50000; %For m Marching n_m_march=10000; if rhof<200 %For Gas Solid flows m_start=0.05; m_jump=0.05; lower_jump_limit=0.001; else %For LiquidSolid flows m_start=0.2; m_jump=0.2; lower_jump_limit=0.00001; end n_m_iterations=20000; %output shows up after every n_m_working iterations n_m_working=500; %Tolerance Level nus_tolerance=0.0000000001; Vfzd_tolerance=0.0000000001; Vszd_tolerance=0.0000000001; Td_tolerance=0.0000000001; kd_tolerance= 0.0000000001; Ed_tolerance=0.0000000001; kTd_tolerance=0.0000000001; %Guess Vsz/Vfz at Center Line %keep 0 PAGE 184 184 if ST>5 && ST<=40 disp('*******************************************************') disp('Running Mid ST number case') disp('*******************************************************') Mid_ST=yes; end if ST>40 disp('*******************************************************') disp('Running High ST number case') disp('*******************************************************') High_ST=yes; end PhD_Guess_Profiles.m %Guess Profiles %Guess for Mfix wall con ditions if Vfz_wall_Mfix==yes  k_wall_Mfix==yes if Re<50001 load Re_50000_wall end if Re>50000 && Re<100001 load Re_100000_wall end if Re>100000 && Re<200001 load Re_200000_wall end if Re>200000 && Re<300001 load Re_300000_wall end if Re>300000 && Re<400001 load Re_400000_wall end if Re>400000 && Re<500001 load Re_500000_wall end if Re>500000 load Re_500000_wall end %Guess for Noslip condition else if Re<50001 load Re_50000 end if Re>50000 && Re<100001 load Re_100000 end if Re>100000 && Re<200001 load Re_200000 end if Re>200000 && Re<300001 load Re_300000 end if Re>300000 && Re< 400001 load Re_400000 end if Re>400000 && Re<500001 PAGE 185 185 load Re_500000 end if Re>500000 load Re_500000 end end PhD_Regriding.m %Adjusting guess grid to User spefified grid clear temp load User_Specified_Grid [n,dumm y]=size(r_User_Defined); %Finding number of new grid points rdim_UD=r_User_Defined/R; temp=zeros(3,n); %Interpolating values for the new grid for i=1:n temp(1,i)=interp1(rdim,Vfzd,rdim_UD(i)); temp(2,i)=interp1(rdim,kd,rdim_UD(i)); temp(3,i)=interp1(rdim,Ed,rdim_UD(i)); end clear Vfzd kd Ed dummy Vfzd=temp(1,:)'; kd=temp(2,:)'; Ed=temp(3,:)'; rdim=rdim_UD; clear temp rdim_UD %Clearing the unnecessary PhD_1_Strain.m %Calculating single phase strain for i=1:n if i==1  i==n if i==1 dVfzdr(i)=0; else h1=drb(n); h2=drb(n)+drb(n 1); dVfzdr(i)=h1*h2/(h2h1)*((Vfzd(i)Vfzd(i 1))/h1^2 (Vfzd(i)Vfzd(i 2))/h2^2); end else h1=drf(i); h2=drb(i); dVfzdr (i)=h1*h2/(h1+h2)*((Vfzd(i+1) Vfzd(i))/h1^2+(Vfzd(i) Vfzd(i 1))/h2^2); end end PhD_2_Strain.m %Calculating two phase strain for i=1:n if i==1  i==n if i==1 dVfzdr(i)=0; dVszdr(i)=0; else h1=drb(n); h2=drb(n)+drb(n 1); PAGE 186 186 dVfzdr(i)=h1*h2/(h2h1)*((Vfzd(i)Vfzd(i 1))/h1^2 (Vfzd(i)Vfzd(i 2))/h2^2); dVszdr(i)=h1*h2/(h2 h1)*((Vszd(i)Vszd(i 1))/h1^2 (Vszd(i)Vszd(i 2))/h2^2); end else h1=drf(i); h2=drb(i); dVfzdr(i)=h1*h2/(h1+h2)*((Vfzd(i+1) Vfzd(i))/h1^2+(Vfzd(i) Vfzd(i 1))/h2^2); dVszdr(i)=h1*h2/(h1+h2)*((Vszd(i+1) Vszd(i))/h1^2+(Vszd(i)Vszd(i 1))/h2^2); end end PhD_1_Viscosity.m %Calculating single phase viscosity % Calculating dimensionless friction velocity UTd=( dVfzdr(n))^0.5; %Calculating eddy viscosity %Using Wall condtions from MFIX if Vfz_wall_Mfix==yes  k_wall_Mfix==yes for i=1:n RT(i)=kd(i)^2/Ed(i); yplus(i)=UTd*(rd(n) rd(i)); f mu=1; muTd(i)=cmu*fmu*kd(i)^2/Ed(i); dmuTdkd(i)=(2*cmu*fmu*kd(i)/Ed(i)); dmuTdEd(i)=( cmu*fmu*kd(i)^2/Ed(i)^2); end %Using the noslip boundary conditions else for i=1:n RT(i)=kd(i)^2/Ed(i); yplus(i)=UTd*(rd( n) rd(i)); %For no slip at the wall if i==n fmu=0; muTd(i)=0; %for rest of the pipe else fmu=(1 exp( yplus(i)/70))*(1+3.45/RT(i)^0.5); muTd(i)=cmu*fmu*kd(i)^2/Ed(i); end dRTdkd=2*kd(i)/Ed(i); dfmudkd= (1 exp( yplus(i)/70))*(3.45/2*(RT(i)^ 1.5)*dRTdkd); dmuTdkd(i)=(2*cmu*fmu*kd(i)/Ed(i))+cmu*kd(i)^2/Ed(i)*dfmudkd; %Calculating dmuTd/dkd dRTdEd= kd(i)^2/Ed(i)^2; dfmudEd= (1 exp( yplu s(i)/70))*(3.45/2*(RT(i)^ 1.5)*dRTdEd); dmuTdEd(i)=( cmu*fmu*kd(i)^2/Ed(i)^2)+cmu*kd(i)^2/Ed(i)*dfmudEd; %Calculating dmuTd/dkd end PhD_2_Viscosity.m %Calculating two phase viscosities UTd=(abs(dVfzdr(n)))^0.5; %Calculating Friciton velocity %Calculating muefd if Batchelor_Green==yes for i=1:n muefd(i)=(1+2.5*nus(i)+7.6*nus(i)^2)*(1nus(i)/nus0); PAGE 187 187 end end %k E model for muT if Vfz_wall_Mfix==yes  k_wall_Mfix==yes for i=1:n RT(i)=kd(i)^2/Ed(i)/muefd(i); ypl us(i)=UTd*(rd(n)rd(i)); fmu=1; muTd(i)=cmu*fmu*kd(i)^2/Ed(i); dmuTdkd(i)=(2*cmu*fmu*kd(i)/Ed(i)); dmuTdEd(i)=( cmu*fmu*kd(i)^2/Ed(i)^2); end else for i=1:n RT(i)=kd(i)^2/Ed(i)/muefd(i); yplus(i)=UTd* (rd(n)rd(i)); if i==n fmu=0; muTd(i)=0; else fmu=(1 exp( yplus(i)/70))*(1+3.45/RT(i)^0.5); muTd(i)=cmu*fmu*kd(i)^2/Ed(i); %Calculating muTd end dRTdkd=2*kd(i)/Ed(i)/muefd(i); dfmudkd= (1 exp( yplus(i)/70))*(3.45/2*(RT(i)^ 1.5)*dRTdkd); dmuTdkd(i)=(2*cmu*fmu*kd(i)/Ed(i))+cmu*kd(i)^2/Ed(i)*dfmudkd; %Calculating dmuTd/dkd dRTdEd= kd(i)^2/Ed(i)^2/muefd(i); dfmudEd= (1 exp( yplus(i)/70))*(3.45/2*(RT(i )^ 1.5)*dRTdEd); dmuTdEd(i)=( cmu*fmu*kd(i)^2/Ed(i)^2)+cmu*kd(i)^2/Ed(i)*dfmudEd; %Calculating dmuTd/dEd end end %Granular Kinetic Thoery Parameters if GKT==yes PhD_mus_lambda_Models end PhD_mus_lambda_Models.m %lambda_mus_models %For Low Stokes number if Low_ST==yes for i=1:n if i==n musdG(i)=muefd(i)+0; lambdadG(i)=muefd(i)+0; else temp1=Rep^2/18*rhos/rhof*2/3/cmu*Ed(i)/kd(i); musdG(i)=muefd(i)+muTd(i)/(1+temp1); lambdadG(i)=muefd(i)+muTd(i)/sigmak/(1+temp1); end end end %For Mid and High Stokes number if Mid_ST==yes  High_ST==yes %High Stokes number relations if Lun==yes for i=1:n PAGE 188 188 musd(i)=5/96*pi^0.5*rhos/rhof*Rep*Td(i)^0.5; %Calculating musd lambdad(i)=25/128*pi^0.5*rhos/rhof*Rep*Td(i)^0.5; %Calculating lambdad w=(1+dR/6/2^0.5/nus(i))^ 1; g0=(nus0^(1/3))/(nus0^(1/3) nus(i)^(1/3)); G1k(i)=(1+8/5*eta*nus(i)*g0*(3*eta2))/ (eta*(2 eta)*g0); G1c(i)=(((1+8/5*eta*nus(i)*g0*(3*eta2))*8*nus(i))/(5*(2 eta)))+(768*nus(i)^2*g0*eta/25/pi); G2k(i)=nus(i); G2c(i)=4*eta*nus(i)^2*g0; G3k(i)=8*(1+12*eta^2*nus(i)*g0*(4*eta3)/5)/eta/(41 33*eta)/g0; G3c(i)=96*nus(i)*(1+12*eta^2*nus(i)*g0*(4*eta 3)/5+16*eta*nus(i)*g0*(41 33*eta)/15/pi)/5/(41 33*eta); alpha(i)=(w*G2k(i)+G2c(i)); %Calculating alpha dwdnus=((1+dR/6/2^0.5/nus(i))^ 2)*dR/6/2^0.5/nus(i)^2; dg0dnus=(nus0^(1/3))/3/((nus0^(1/3) nus(i)^(1/3))^2)/nus0^(2/3); dalphadnus(i)=w+nus(i)*dwdnus+8*eta*nus(i)*g0+4*eta*nus(i)^2*dg0dnus; %Calculating dalpha/dnus dmusdT(i)=5/96*pi^0.5*rhos/rhof*Rep*(0.5*Td(i)^ 0.5); d musdT(i)=dmusdT(i)*(w*G1k(i)+G1c(i)); %Calculating dmusd/dTd musdG(i)=musd(i)*(w*G1k(i)+G1c(i)); %Calculating musd*(wG1k+G1c) lambdadG(i)=lambdad(i)*(w*G3k(i)+G3c(i)); %Calculating lambdad*(wG1k+G1c) end end if Piera no_Leckner==yes cbeta=1.81.35; for i=1:n if i PAGE 189 189 temp3=muefd(i)+muTd(i)/sigmak/(1+temp1); end musdG(i)=(ST 5)/(40 5)*musdG(i)+(40 ST)/(40 5)*temp2; lambdadG(i)=(ST 5)/(40 5)*lambdadG(i)+(40 ST)/(40 5)*temp3; end end end PhD_Relative_Velocity.m %Relative Velocity Vrd=Vfzd Vszd; for i=1:n Vrd2(i)=(Vrd(i)^2+8*Td(i)/pi)^0.5; end PhD_Drag_Models.m %Drag models for Ding & Gidspow type formulation %Wen and Yu if Wen_Yu==yes for i=1:n Red=(1 nus(i)) *Rep*abs(Vrd(i)); if Red>1000 CD(i)=0.44; else CD(i)=24/Red*(1+0.15*Red^0.687); end beta(i)=3/4/Rep*CD(i)*nus(i)/(1 nus(i))^2.65*abs(Vrd(i)); end end %Hill Koch and Ladd if Hill_Koch_Ladd==yes for i=1:n Red=(1 nus(i))*Rep/2*abs(Vrd(i)); w_HKL=exp( 10*(0.4nus(i))/nus(i)); if nus(i)<0.4 F0=(1 w_HKL)*(1+3*(nus(i)/2)^0.5+135/64*nus(i)*log(nus(i))+17.14*nus(i))/(1+0.681*nus(i) 8.48*nus(i)^2+8.16*nus(i)^3)+w_HKL*1 0*nus(i)/(1 nus(i)^3); else F0=10*nus(i)/(1 nus(i)^3); end if nus(i)<0.1 F1=(2/nus(i))^0.5/40; else F1=0.11+0.00051*exp(11.6*nus(i)); end if nus(i)<0.4 F2=(1 w_ HKL)*(1+3*(nus(i)/2)^0.5+135/64*nus(i)*log(nus(i))+17.89*nus(i))/(1+0.681*nus(i) 11.03*nus(i)^2+15.41*nus(i)^3)+w_HKL*10*nus(i)/(1nus(i)^3); else F2=10*nus(i)/(1 nus(i)^3); end if nus(i)<0.0953 F3=0.09351*nu s(i)+0.03667; else F3=0.0673+0.212*nus(i)+0.0232/(1nus(i))^5; PAGE 190 190 end F4=(F3+(F3^2 4*F1*(F0 F2))^0.5)/2/F1; F5=(F2 1)/(3/8 F3); if nus>0.01 if Red>F5 F=F2+F3*Red; else F=F0+F1*Red^2; end else if Red>F4 F=F2+F3*Red; else F=1+3/8*Red; end end CD(i)=4/3*(17.3/2/Red+0.336); beta(i)=18*F/Rep^2*(1nu s(i))*nus(i); end end %Syamlal and O'Brien if Syamlal_OBrien==yes for i=1:n Red=Rep*abs(Vrd(i)); A_SO=(1 nus(i))^4.14; if (1 nus(i))>0.85 B_SO=(1 nus(i))^2.65; else B_SO=0.8*(1nus(i))^1.28; end Vrm=0.5*(A_SO 0.06*Red+((0.06*Red)^2+0.12*Red*(2*B_SO A_SO)+A_SO^2)^0.5); CD(i)=(0.63+4.8*(Vrm/Red)^0.5)^2; beta(i)=3/4*nus(i)/Vrm^2/Rep*CD(i)*abs(Vrd(i)); end end PhD_2_Generation.m %Calculating Generation and corresponding slope terms %Calculating Solid Stress for i=1:n if GKT==yes sigmarzd(i)= musdG(i)*dVszdr(i); dsigmarzdT(i)= dmusdT(i)*dVszdr(i); end end %Compiling the Generation terms for all the equations for i=1:n %The Pressure Drop / Gravity term is solved internaly G_Vfzd(i)=0; dG_Vfzd(i)=0; if downward==yes G_Vszd(i)=Ar/Rep^3*nus(i); else G_Vszd(i)= Ar/Rep^3*nus(i); end dG_Vszd(i)=0; PAGE 191 191 G_Td(i)= sigmarzd(i)*dVszdr(i); dG_Td(i)= wT(1)*d sigmarzdT(i)*dVszdr(i); G_kd(i)=(1 nus(i))*muTd(i)*dVfzdr(i)^2; dG_kd(i)=(1nus(i))*wk(1)*dmuTdkd(i)*dVfzdr(i)^2; if Vfz_wall_Mfix==yes  k_wall_Mfix==yes G_Ed(i)=(1 nus(i))*cT1*muTd(i)*dVfzdr(i)^2*Ed(i)/kd(i); dG_Ed(i)=(wE(1)* cT1*dVfzdr(i)^2/kd(i))*(1nus(i))*(muTd(i)+Ed(i)*dmuTdEd(i)); else if i==n G_Ed(i)=0; dG_Ed(i)=0; else G_Ed(i)=(1 nus(i))*cT1*fT1*muTd(i)*dVfzdr(i)^2*Ed(i)/kd(i); dG_Ed(i)=(wE(1)*cT1*fT1*dVfzdr(i)^2/kd(i))*(1 nus(i))*(muTd(i)+Ed(i)*dmuTdEd(i)); end end end PhD_2_Dissipation.m %Calculating Dissipation and corresponding slope terms %Calculating Granular Temperature disspation for i=1:n if GKT==yes g0=(nus0^(1/3))/(nus0^(1/3) nus(i)^(1/3)); if Kunn_mod_diff_coeff_table_4_1==yes gamma(i)=48/pi^0.5*etaf*(1 etaf)*g0*nus(i)^2/Rep*(rhos/rhof)*Td(i)^1.5; dgammadT(i)=48/pi^0.5*etaf*(1 etaf)*g0*nus(i)^2/Rep*(rhos/rhof)*1.5*Td(i)^0.5; els e gamma(i)=48/pi^0.5*eta*(1 eta)*g0*nus(i)^2/Rep*(rhos/rhof)*Td(i)^1.5; dgammadT(i)=48/pi^0.5*eta*(1eta)*g0*nus(i)^2/Rep*(rhos/rhof)*1.5*Td(i)^0.5; end end if i==n fT2(i)=0; else fT2(i)=(1 2/9*ex p( (RT(i)/6)^2))*(1 exp( yplus(i)/5))^2; end end %Compiling the dissipation terms for all the equations for i=1:n D_Vfzd(i)=0; dD_Vfzd(i)=0; D_Vszd(i)=0; dD_Vszd(i)=0; if Low_ST==yes D_Td(i)=0; dD_Td(i)=0; end if Mid_ST==yes %D_Td(i)=(40 ST)/(40 5)*0+(ST 5)/(40 5)*gamma(i); D_Td(i)=0; dD_Td(i)=0; end if High_ST==yes D_Td(i)=gamma(i); dD_Td(i)=wT(2)*dgammadT(i); PAGE 192 192 end D_kd(i)=(1 nus(i))*Ed(i); dD_kd(i)=0 ; if Vfz_wall_Mfix==yes  k_wall_Mfix==yes D_Ed(i)=(1 nus(i))*cT2*Ed(i)^2/kd(i); dD_Ed(i)= 2*(wE(2)*cT2*(1 nus(i))*Ed(i)/kd(i)); else if i==n D_Ed(i)=0; dD_Ed(i)=0; else D_Ed(i)=( 1 nus(i))*cT2*fT2(i)*Ed(i)^2/kd(i); dD_Ed(i)= 2*(wE(2)*cT2*fT2(i)*(1 nus(i))*Ed(i)/kd(i)); end end end PhD_2_Interaction.m %Calculating Interaction and corresponding slope terms %The Drag Term FD if Ding_Gidaspow_1990==yes for i=1:n Fd(i)=beta(i)*(1 nus(i))*(Vrd(i)); dbetadVfzd=3/4/Rep*CD(i)*nus(i)/(1nus(i))^2.65; dFdVfzd(i)=beta(i); dbetadVszd=3/4/Rep*CD(i)*nus(i)/(1nus(i))^2.65; dFdVszd(i)= beta(i); end end %FET model + Sinclair & Mallo if FET==yes for i=1:n kTd(i)=(6*kd(i)*Td(i))^0.5; dkTdT(i)=0.5*(6*kd(i))^0.5*Td(i)^ 0.5; dkTdk(i)=0.5*(6*Td(i))^0.5*kd(i)^ 0.5; g0=(nus0^(1/3))/(nus0^(1/3) nus(i)^(1/3)); if High_ST==yes if ST >100 beta2(i)=24*nus(i)*g0/Rep*(Td(i)/pi)^0.5; end if ST<=100 beta2(i)=(1 nus(i))*beta(i)*rhof/rhos/nus(i); end end if Mid_ST==yes temp1=(1 nus(i))*beta(i); temp2=(1 nus(i))*beta(i)*rhof/rhos/nus(i); beta2(i)=(6.7*10^ 7*Re18.44*dR+0.62)*temp1; end if Low_ST==yes beta2(i)=(1 nus(i))*beta(i); end bkd(i)=beta2(i)*2*kd(i); bTd(i)=beta2(i)*3* Td(i); bkTd(i)=beta2(i)*kTd(i); FKS(i)=bkTd(i) bTd(i); PAGE 193 193 dFKSdT(i)= beta2(i)*(dkTdT(i) 3); FKF(i)=bkd(i)bkTd(i); dFKFdk(i)=beta2(i)*(2 dkTdk(i)); if i==n dFKFdk(i)=0; end end end %TVBA + Louge if Louge_etal_1991==yes for i=1:n kTd(i)=4/pi^0.5*rhof/rhos*Rep*beta(i)*(1 nus(i))/nus(i)*(Vrd(i))^2/Td(i)^0.5; %dkTdT(i)=0.5*(4/pi^0.5*rhof/rhos*Rep*beta(i)*(1nus(i))/nus(i)*(Vrd(i))^2)*Td(i)^ 1.5; dkTdT(i)=0; dkTdk(i)=0; bkd(i)=beta(i)*(1 nus(i))*2*kd(i); bTd(i)=beta(i)*(1 nus(i))*3*Td(i); bkTd(i)=beta(i)*(1nus(i))*kTd(i); FKS(i)=bkTd(i) bTd(i); dFKSdT(i)= beta(i)*(1 nus(i))*(dkTdT(i) 3); FKF(i)=bkd(i)bkTd(i); dFKFdk(i)=beta(i)*(1nus(i))*(2 dkTdk(i)); if i==n %FKF(i)=0; %FKS(i)= bTd(i); dFKFdk(i)=0; end end end %TVBA + Sinclair & Mallo if Sinclair_and_Mallo_1998==yes; for i=1:n kTd(i )=(6*kd(i)*Td(i))^0.5; dkTdT(i)=0.5*(6*kd(i))^0.5*Td(i)^ 0.5; dkTdk(i)=0.5*(6*Td(i))^0.5*kd(i)^ 0.5; bkd(i)=beta(i)*(1 nus(i))*2*kd(i); bTd(i)=beta(i)*(1 nus(i))*3*Td(i); bkTd(i)=beta(i)*(1nus(i))*kTd(i); FK S(i)=bkTd(i) bTd(i); dFKSdT(i)= beta(i)*(1 nus(i))*(dkTdT(i) 3); FKF(i)=bkd(i)bkTd(i); dFKFdk(i)=beta(i)*(1nus(i))*(2 dkTdk(i)); if i==n %FKF(i)=0; %FKS(i)= bTd(i); dFKFdk(i)=0; end end end %TVBA + Simonin if Simonin_1996==yes cbeta=1.8 1.35; for i=1:n if i==n if Vfz_wall_Mfix==yes  k_wall_Mfix==yes zetar=3*(Vrd(i))^2/2/kd(i); PAGE 194 194 etat=3*cmu/2/(1+cbeta*zetar)^0.5*kd(i)/ Ed(i)*rhof/rhos*(1 nus(i))/nus(i)*beta(i); chi=nus(i)*rhos/rhof/(1nus(i)); kTd(i)=etat*(2*kd(i)+3*chi*Td(i))/(1+(1+chi)*etat); else kTd(i)=0; end else zetar=3*(Vrd(i))^2/2/kd(i); etat=3*cmu/2/(1+cbeta*zetar)^0.5*kd(i)/Ed(i)*rhof/rhos*(1 nus(i))/nus(i)*beta(i); chi=nus(i)*rhos/rhof/(1 nus(i)); kTd(i)=etat*(2*kd(i)+3*chi*Td(i))/(1+(1+chi)*etat); end bkd(i)=beta(i)*(1 nus(i))*2*kd(i); bTd(i)=beta(i)*(1 nus(i))*3*Td(i); bkTd(i)=beta(i)*(1nus(i))*kTd(i); FKS(i)=bkTd(i) bTd(i); dFKSdT(i)= beta(i)*(1 nus(i))*(dkTdT(i) 3); FKF(i)=bkd(i)bkTd(i); dFKFdk(i)=beta(i)*(1nus(i))*(2 dkTdk(i)); if i==n %FKF(i)=0; %FKS(i)= bTd(i); dFKFdk(i)=0; end end end %TVBA + Igci if Koch_Sangani_1999==yes for i=1:n g0=(nus0^(1/3))/(nus0^(1/3) nus(i)^(1/3)); Rdiss=1+(3* (nus(i)/2)^0.5)+(135/64*nus(i)*log(nus(i)))+(11.26*nus(i)*(1 5.1*nus(i)+16.57*nus(i)^221.77*nus(i)^3)) (nus(i)*g0*log(0.01)); if nus(i)<0.4 Rd=(1+(3*(nus(i)/2)^0.5)+(135/64*nus(i)*log(nus(i)))+(17.14*nus(i)))/(1+0.681*nus(i) 8.48*nus(i )^2+8.16*nus(i)^3); else Rd=10*nus(i)/(1 nus(i))^3+0.7; end phi=Rd^2/(1+3.5*nus(i)^0.5+5.9*nus(i)); bkd(i)=beta(i)*(1 nus(i))*2*kd(i); %bkd(i)=(36*nus(i)*muefd(i)*kd(i)*Rdiss/Rep^2); bTd(i)=beta(i )*(1 nus(i))*3*Td(i); %bTd(i)=(54*nus(i)*muefd(i)*Td(i)*Rdiss/Rep^2); bkTd(i)=(81*nus(i)*muefd(i)^2*(Vrd(i))^2/g0/(pi*Td(i))^0.5/Re^3*phi); FKS(i)=bkTd(i) bTd(i); dFKSdT(i)=0.5*(81*nus(i)*muefd(i)^2*(Vrd(i))^2/g0/(pi)^0.5/Re^3)*Td(i)^ 1.5+beta(i)*(1nus(i))*3; FKF(i)=bkd(i)bkTd(i); dFKFdk(i)=beta(i)*(1nus(i))*2; if i==n %FKF(i)=0; %FKS(i)= bTd(i); dFKFdk(i)=0; end end end %TVBA + Wylie PAGE 195 195 if Wylie_etal_2003==yes for i=1:n Red=Rep*abs(Vrd(i)); ReT=Rep*Td(i)^0.5; phi=(1+2*ReT^2/Red^2ReT^4/Red^4)*erf(Red/2^0.5/ReT)+(2/pi)^0.5*ReT/Red*(1+ReT^2/Red^2)*exp( Red^2/2/ReT^2); Kfb=0.0336+0.106*nus(i)+0.0116*(1nus(i))^ 5; Rd0=(1+3/(2*nus(i)^0.5)^0.5+135/64*nus(i)*log(nus(i))+17.14*nus(i))/(1+0.681*nus(i) 8.48*nus(i)^2+8.16*nus(i)^3); Rd=Rd0+Kfb*Red*phi; chi=(1+2.5*nus(i)+4.509*nus(i)^2+4.5154*nus(i)^3)/(1 (nus(i)/0.6436)^3)^0.6780; RS=(chi*(1+3.5*nus (i)^0.5+5.9*nus(i)))^ 1; bkd(i)=beta(i)*(1 nus(i))*2*kd(i); bTd(i)=beta(i)*(1 nus(i))*3*Td(i); bkTd(i)=(81*nus(i)*muefd(i)^2*(Vrd(i))^2*rhof/rhos/(pi*Td(i))^0.5/Re^3*RS*Rd^2); FKS(i)=bkTd(i) bTd(i); dFKSdT(i)=0.5*(81*nus(i)*muefd(i)^2*(Vrd(i))^2*rhof/rhos/(pi)^0.5/Re^3*RS*Rd^2)*Td(i)^ 1.5+beta(i)*(1 nus(i))*3; FKF(i)=bkd(i)bkTd(i); dFKFdk(i)=beta(i)*(1nus(i))*2; if i==n %FKF(i)=0; %FKS(i)= bTd(i); dFKFdk(i) =0; end end end %VBA + Zhang & Reese if Zhang_Reese_2003==yes; for i=1:n g0=(nus0^(1/3))/(nus0^(1/3) nus(i)^(1/3)); kTd(i)=(pi^0.5/24/g0*rhof/rhos*Rep*beta(i)*(1nus(i))/nus(i)^2*(Vrd2(i))^2/Td(i)^0.5); %dkTdT(i)=0.5*(pi^0.5/24/g0*rhof/rhos*Rep*beta(i)*(1 nus(i))/nus(i)^2*(Vrd2(i))^2)*Td(i)^ 1.5; dkTdT(i)=0; dkTdk(i)=0; bkd(i)=beta(i)*(1 nus(i))*2*kd(i); bTd(i)=beta(i)*(1 nus(i))*3*Td(i); bkTd(i)=beta(i)*(1nus(i))*kTd(i); FKS(i)=bkTd(i) bTd(i); dFKSdT(i)= beta(i)*(1 nus(i))*(dkTdT(i) 3); FKF(i)= beta(i)*(1nus(i))*Vrd(i)^2 bTd(i)+bkTd(i); dFKFdk(i)=beta(i)*(1nus(i))*(2 dkTdk(i)); if i==n %FKF(i)=0; %FKS(i)= bTd(i ); dFKFdk(i)=0; end end end %TVBA + Chen & Wood if Chen_Wood==yes for i=1:n if i==n kTd(i)=0; else temp1=Rep^2/18/rhos/rhof*2/3/cmu*Ed(i)/kd(i); PAGE 196 196 kTd(i)=2*kd(i)*exp( 0.0825*temp 1); end dkTdT(i)=0; dkTdk(i)=2; bkd(i)=beta(i)*(1 nus(i))*2*kd(i); bTd(i)=beta(i)*(1 nus(i))*3*Td(i); bkTd(i)=beta(i)*(1nus(i))*kTd(i); FKS(i)=bkTd(i) bTd(i); dFKSdT(i)= beta(i)*(1 nus(i))*(d kTdT(i) 3); FKF(i)=bkd(i)bkTd(i); dFKFdk(i)=beta(i)*(1nus(i))*(2 dkTdk(i)); if i==n %FKF(i)=0; %FKS(i)= bTd(i); dFKFdk(i)=0; end end end %Wake modeling for vortex shedding if Lun_Wake_term==yes for i=1:n Redcl=Rep*abs(Vrd(1)); Red=Rep*abs(Vrd(i)); if rhof<200 %For gas solid flows if Redcl<=150 temp1=1; temp2=0; end if Redcl>= 150 && Redcl<310 temp1=0.017*Red; end if Redcl>=310 && Redcl<610 temp1=1.2+0.00005*Red^2; end if Redcl>=610 && Redcl<1560 temp1=0.029*Red; end if Redcl>1560 temp1=0.029*Red; end if Redcl>150 && Redcl<310 temp2=10/3; end if Redcl>=310 temp2=24/3; end else %For liquid s olid flows if Redcl<500 temp1=1; temp2=0; end if Redcl<=500 temp1=0.017*Red; temp2=10/3; PAGE 197 197 end end Ew(i)=12*nus(i)*temp1*temp2*kd( i)/Rep^2; end end %Compiling the Interaction terms for all the equations for i=1:n I_Vfzd(i)= Fd(i); dI_Vfzd(i)=wVfz(1)*dFdVfzd(i); I_Vszd(i)=Fd(i); dI_Vszd(i)=wVsz(1)*dFdVszd(i); I_Td(i)=FKS(i); dI_Td(i)=wT(3)*dFKSdT(i); I_ kd(i)= FKF(i)+Ew(i); dI_kd(i)=wk(2)*dFKFdk(i); if Vfz_wall_Mfix==yes  k_wall_Mfix==yes I_Ed(i)=cT3*Ed(i)/kd(i)*( FKF(i)+Ew(i)); dI_Ed(i)=0; else if i==n I_Ed(i)=0; dI_Ed(i)=0; else I_Ed(i)= cT3*fT2(i)*Ed(i)/kd(i)*( FKF(i)+Ew(i)); dI_Ed(i)=0; end end end PhD_Thomas_Algorithims.m %Forward reduction of matrix for i=1:n if i==1 A(i,3)=A(i,3)/A(i,2); A(i,4)=A(i,4)/A(i,2); else A(i,3)=A(i,3)/(A(i,2) A(i 1,3)*A(i,1)); A(i,4)=(A(i,4) A(i 1,4)*A(i,1))/(A(i,2) A(i 1,3)*A(i,1)); end end %Back substitution to get new values for i=1:n if i==1 X(n+1 i)=A(n+1 i,4); else X(n+1 i)=A(n+1 i,4) A(n+1 i,3)*X (n+2 i); end end PhD_1_Vfz_Eqn.m %Solving the single phase Vfz equation %Constructing the matrix, using principals of Finite Volume %The matrix A is actually a tridiagonal matrix but is constructed in form %of 4 columnar arrys. The first column is the sub diagonal, the second %column is a daigonal and the third is the super diagonal, while the forth PAGE 198 198 %column is the B matrix in AX=B where A is the orignal tri diagonal variable=2; for i=1:n C(i)=1+muTd(i); %Constructing Diffusion term end for i=1:n if i==1  i==n if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)=0; A(i,2)=a2; A(i,3)= a2; A(i,4)=0; %Inserting operating condition %If the center line velocity is known if Re_center_line==yes B1(i)=1; %Average velocity is known else B1(i)=rd(i)*drf(i)/Re^2; end %Inserting Pressure Drop as a variable B2(i)=1; end if i==n PhD_1_BC %Wall Boundary Condition %For operating condition %If the center line velocity is known if Re_center_line==yes B1(i)=0; %Average velocity is known else B1(i)=rd(i)*drf(i 1)/Re^2; end end else rd2=((( rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; A(i,1)= a1; A(i,2)=a1+a2; A(i,3)= a2; A(i,4)=0; %For operating condition %If the center line velocity is known if Re_center_line==yes B1(i)=0; %Average velocity is known else B1(i)=rd(i)*(drf(i)+drf(i 1))/R e^2; end %Inserting Pressure Drop as a variable B2(i)=1; end end PAGE 199 199 %Solving the velocity equation temp3=0; %For operating condition %If the center line velocity is known if Re_center_line==yes temp4=1; %Average velocity is known else temp4=0.25; end %Modifying the system so that Thomas Algorithims can solve it for i=1:n if i==n temp1=B1(i)/A(i,2); B1(i)=B1(i)A(i,2)*temp1; temp4=temp4 A(i,4)*temp1; temp3=temp3 B2(i)*temp1; else temp1=A(i+1,1)/A(i,2); temp2=B1(i)/A(i,2); A(i+1,1)=A(i+1,1) A(i,2)*temp1; A(i+1,2)=A(i+1,2) A(i,3)*temp1; A(i+1,4)=A(i+1,4) A(i,4)*temp1; temp4=temp4 A(i,4)*temp2; B1(i)=B1(i)A(i,2)*temp2; B1(i+1 )=B1(i+1) A(i,3)*temp2; B2(i+1)=B2(i+1) B2(i)*temp1; temp3=temp3 B2(i)*temp2; end end temp1=temp3; temp3=temp3/temp1; temp4=temp4/temp1; for i=1:n temp1=B2(i); B2(i)=B2(i) temp3*temp1; A(i,4)=A(i,4) temp4*temp1; end %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values for i=1:n Vfzd(i)=X(i); end PhD_1_k_E_Eqn.m %Solving the k E equations for Single Phase flow %Solving the k equation for i=1:n C(i)=1+muTd(i)/sigmak; %Constructing the diffusion term D(i)=(muTd(i)*dVfzdr(i)^2)Ed(i); %Constructing the source/sink terms S(i)=wk(1)*dmuTdkd(i)*dVfzdr(i)^2; %Constructing the slope of source/sink terms end %Constructing the matrix, using principals of Finite Volume %The matrix A is actually a tridi agonal matrix but is constructed in form %of 4 columnar arrys. The first column is the sub diagonal, the second PAGE 200 200 %column is a daigonal and the third is the super diagonal, while the forth %column is the B matrix in AX=B where A is the orignal tri diagonal variable=5; for i=1:n if i==1  i==n %Inputing the Boundary Conditions if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)=0; A(i,2)=a2 S(i); A(i,3)= a2; A(i,4)=D(i) S(i)*kd(i); else %Wall Boundary condition PhD_1_BC end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; A(i,1)= a1; A(i,2)=a1+a2 S(i); A(i,3)= a2; A( i,4)=D(i) S(i)*kd(i); end end %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values kd=X; %Solving the E equation %if wall conditions are applied C, D & S are constructed differently at the wall as given in MFIX manual if Vfz_wall_Mf ix==yes  k_wall_Mfix==yes for i=1:n C(i)=1+muTd(i)/sigmaE; %Constructing the Diffusion term D(i)=(cT1*muTd(i)*dVfzdr(i)^2*Ed(i)/kd(i)) (cT2*Ed(i)^2/kd(i)); %Constructing the source/sink terms S(i)=(wE(1)*cT1*dVfzdr(i)^2/kd(i))*(muTd(i)+Ed(i)*dmuTdEd(i)) 2*(wE(2)*cT2*Ed(i)/kd(i)); %Constructing the slope of source/sink terms end %for all other conditions C, D & S are constructed as Bolio code else for i=1:n C(i)=1+muTd(i)/sigmaE; %Constructing the Diffusion term if i PAGE 201 201 %of 4 columnar arrys. The first column is the sub diagonal, the second %column is a daigonal and the third is the super diagonal, while the forth %column is the B matrix in AX=B where A is the orignal tri diagonal va riable=6; for i=1:n if i==1  i==n %Inputing the Boundary Conditions if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)=0; A(i,2)=a2 S(i); A(i,3)= a2; A(i,4)=D(i) S(i)*Ed(i); else %Wall Boundary condition Edwall=d^2(kd)/(d(rd))^2 PhD_1_BC end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; A(i,1)= a1; A(i,2)=a1+a2 S(i); A(i,3)= a2; A(i,4)=D(i)S(i)*Ed(i); end end %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values Ed=X; PhD_1_BC.m %Boundary Conditions %Boundary Condition for Vfz if variable==2 %No slip if Vfz_no_slip==yes A(i,1)=0; A(i,2)=1; A(i,3)=0; A(i,4)=0; B2(i)=0; end %Wall Function from Mfix if Vfz_wall_Mfix==yes h1=rd(i) rd(i 1); h2=cmu^0.25*kd(i)^0.5*h1/2; rd2=((rd(i)^2) ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1+rd(i)/rd2*0.42*cmu^0.25*kd(i)^0.5/log(9.81*h2); A(i,3)=0; A(i,4)=0; B2(i)=1; end end PAGE 202 202 %Boundary Condition for k if variable==5 %No slip if k_no_slip==yes A(i,1)=0; A(i,2)=1; A(i,3)=0; A(i,4)=0; end %Wall Function from Mfix if k_wall_Mfix==yes h1=rd(i) rd(i 1); h2=cmu^0.25*kd(i)^0.5*h1/2; t emp1=2*cmu^0.5*kd(i)*Vfzd(i)/h1/log(9.81*h2); temp2=2*cmu^0.75*kd(i)^1.5/0.42/h1; D(i)=temp1 temp2; S(i)=wk(1)*(2*cmu^0.5*Vfzd(i)/h1/log(9.81*h2)) 10*1.5*(2*cmu^0.75*kd(i)^0.5/0.42/h1); rd2=((rd(i)^2) ((rd(i)+rd(i 1))/2)^2)/ 2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1 S(i); A(i,3)=0; A(i,4)=D(i)S(i)*kd(i); end end %Boundary Condition for E if variable==6 %E at wall=muef/rhof *d^2k/dr^2 if Bolio_equation_2_62==yes A(i,1)=0; A(i,2)=1; A(i,3)=0; h1=drb(n); h2=drb(n)+drb(n 1); A(i,4)=((kd(i)kd(i 1))/h1(kd(i)kd(i 2))/h2)*(2/(h2 h1)); end %Wall Function from Mfix if E_wall_Mfix==yes h1=rd(i) rd(i 1); A(i,1)=0; A(i,2)=1; A(i,3)=0; A(i,4)=2*cmu^0.75*kd(i)^1.5/0.42/h1; end end PhD_2_nus_Eqn.m %Solving nues equation for the Two Phase flow %Constructing the matrix A %The matrix A is constructed in form of 4 columnar arrys. The first column %is unity representing a diagonal i.e. the variable itself. The second %column is the slope term and the third is valuve from the previous %iteration, while the forth column is the discretiz ed mass conservation equation. variable=1; %Filling in the operating condition for center line velocity is given if Re_center_line==yes PAGE 203 203 for i=1:n if i==1  i==n if i==1 A(i,4)=(Vszd(i)+M*Vfzd(i))*rd(i)*drf(i); end if i==n A(i,4)=(Vszd(i)+M*Vfzd(i))*rd(i)*drf(i 1); end else A(i,4)=(Vszd(i)+M*Vfzd(i))*rd(i)*(drf(i)+drf(i 1)); end end temp1=0; temp2=0; for i=1:n if i= =1  i==n if i==1 temp2=temp2+Vfzd(i)*rd(i)*drf(i); end if i==n temp2=temp2+Vfzd(i)*rd(i)*drf(i 1); end else temp2=temp2+Vfzd(i)*rd(i)*(drf(i 1)+drf(i)); end end temp2=M*temp2; %Filling in the operating condition for average velocity is given else for i=1:n if i==1  i==n if i==1 A(i,4)=(Vszd(i))*rd(i)*drf(i); end if i==n A(i,4)=(Vszd(i))*rd(i)*drf(i 1); end else A(i,4)=(Vszd(i))*rd(i)*(drf(i)+drf(i 1)); end end temp1=0; temp2=M*Re^2/4; end %For Low Stokes number if Low_ST==yes for i=1:n C(i)=(muefd(i )+muTd(i))/0.7; %Constructing the diffusion term end for i=1:n 1 if i==1 %Inputing the Boundary Conditions a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i)); A(i,1)=0; A(i,2)=a2; A(i,3)= a2; else PAGE 204 204 a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i) rd(i 1)); a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i)); A(i,1)= a1; A(i,2)=a1+a2; A(i,3)= a2; end end %Using matrix properties to get tri dagonal system for i=2:n 1 temp3=A(i,1); A(i,1)=A(i,1)/temp3; A(i,2)=A(i,2)/temp3; A(i,3)=A(i,3)/temp3; temp4=A(i 1,4); A(i 1,4)=A( i1,4) temp4*A(i,1); A(i,4)=A(i,4) temp4*A(i,2); A(i+1,4)=A(i+1,4) temp4*A(i,3); end A(n,1)=A(n 1,4); A(n,2)=A(n,4); A(n,3)=0; for i=1:n if i==n A(i,4)=temp2; else A(i,4)=0; end end %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values nus=X; end %For Mid and High Stokes number if Mid_ST==yes  High_ST==yes %Relations for High Stokes number for i=1:n D(i)=alpha(i); S(i)=dalphadnus(i); A(i,1)=1; A(i,2)= 1/S(i)/Td(i); A(i,3)=nus(i)D(i)/S(i); end %Using matrix properties to get updated nus values for i=1:n temp1=temp1 A(i,2)*A(i,4); temp2=temp2 A(i,3)*A(i,4); end constant=temp2/temp1; for i=1:n X(i)=A(i,3)constant*A(i,2); end %Updating the values of nus nus=X; %For Mid Stokes numbers if Mid_ST==yes PAGE 205 205 %Taking a weighted average between high and low Stokes number for i=1:n nus(i)=(ST 5)/(40 5)*nus(i)+(40 ST)/(40 5)*M/(M+1); end end end PhD_2_Vfz_Eqn.m %Solving the Vfz equation for the Two Phase flow %Constructing the matrix, using principals of Finite Volume %The matrix A is actually a t ridiagonal matrix but is constructed in form %of 4 columnar arrys. The first column is the sub diagonal, the second %column is a daigonal and the third is the super diagonal, while the forth %column is the B matrix in AX=B where A is the orignal tri diag onal variable=2; for i=1:n C(i)=muefd(i)+muTd(i); %Constructing Diffusion term D(i)=G_Vfzd(i) D_Vfzd(i)+I_Vfzd(i); %Constructing the source/sink terms S(i)=dG_Vfzd(i)+dD_Vfzd(i)+dI_Vfzd(i); %Constructing the slope of source/sink terms end for i =1:n if i==1  i==n if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)=0; A(i,2)=a2 S(i); A(i,3)= a2; A(i,4)=D(i) S(i)*Vfzd(i); %Inserting operating condition %If the center line velocity is known if Re_center_line==yes B1(i)=1; %Average velocity is known else B1(i)=(1 nus(i))*rd(i)*drf(i)/Re^2; end %Inserting Pressure Drop as a variable B2(i)=1; end if i==n PhD_2_BC %Wall Boundary Condition %For operating condition %If the center line velocity is known if Re_center_line==yes B1(i)=0; %Average velocity is known else B1(i)=(1 nus(i))*rd(i)*drf(i 1)/Re^2; end end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; A(i,1)= a1; PAGE 206 206 A(i,2)=a1+a2 S(i); A(i,3)= a2; A(i,4)=D(i)S(i)*Vfzd(i); %For operating condition %If the center line velocity is known if Re_center_line==yes B1(i)=0; %Average velocit y is known else B1(i)=(1 nus(i))*rd(i)*(drf(i)+drf(i 1))/Re^2; end %Inserting Pressure Drop as a variable B2(i)=1; end end %Solving the velocity equation temp3=0; %For operating condition %If the center line velocity is known if Re_center_line==yes temp4=1; %Average velocity is known else temp4=0.25; end %Modifying the system so that Thomas Algorithims can solve it for i=1:n if i==n temp1=B1(i)/A(i,2); B1(i)=B1(i)A(i,2)*temp1; temp4=temp4 A(i,4)*temp1; temp3=temp3 B2(i)*temp1; else temp1=A(i+1,1)/A(i,2); temp2=B1(i)/A(i,2); A(i+1,1)=A(i+1,1) A(i,2)*temp1; A(i+1,2)=A(i+1,2) A(i,3)*temp1; A(i+1,4)=A(i+1,4) A(i,4)*temp1; temp4=temp4 A(i,4)*temp2; B1(i)=B1(i)A(i,2)*temp2; B1(i+1)=B1(i+1) A(i,3)*temp2; B2(i+1)=B2(i+1) B2(i)*temp1; temp3=temp3 B2(i)*temp2; end end temp1=temp3; temp3=temp3/temp1; temp4=temp4/temp1; for i=1:n temp1=B2 (i); B2(i)=B2(i) temp3*temp1; A(i,4)=A(i,4) temp4*temp1; end %Solving the tri diagonal system PhD_Thomas_Algorithims for i=1:n PAGE 207 207 Vfzd(i)=X(i); end %Updating the Pressure Drop dPdz=temp4; PhD_2_Vsz_Eqn.m %Solving the Vsz equations for Two Phase f low variable=3; for i=1:n C(i)= musdG(i); %Constructing the diffusion term D(i)=G_Vszd(i)D_Vszd(i)+I_Vszd(i); %Constructing the source/sink terms S(i)=dG_Vszd(i)+dD_Vszd(i)+dI_Vszd(i); %Constructing the slope of source/sink terms end %Constru cting the matrix, using principals of Finite Volume %The matrix A is actually a tridiagonal matrix but is constructed in form %of 4 columnar arrys. The first column is the sub diagonal, the second %column is a daigonal and the third is the super diagonal, while the forth %column is the B matrix in AX=B where A is the orignal tri diagonal for i=1:n if i==1  i==n %Inputing the Boundary Conditions if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^ 2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)=0; A(i,2)=a2+S(i); A(i,3)= a2; A(i,4)= D(i)+S(i)*Vszd(i); else PhD_2_BC end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; A(i,1)= a1; A(i,2)=a1+a2+S(i); A(i,3)= a2; A(i,4)= D(i)+S(i)*Vszd(i); end end %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values Vszd=X; PhD_2_T_Eqn.m %Solving the T equations for Two Phase flow variable=4; %Constructing th e matrix, using principals of Finite Volume %The matrix A is actually a tridiagonal matrix but is constructed in form %of 4 columnar arrys. The first column is the sub diagonal, the second %column is a daigonal and the third is the super diagonal, while t he forth %column is the B matrix in AX=B where A is the orignal tri diagonal %For the case of Low Stokes number if Low_ST==yes PAGE 208 208 Td=kd; end %For the case of Mid and High Stokes number if Mid_ST==yes  High_ST==yes for i=1:n C(i)= lambdadG(i ); %Constructing the diffusion term D(i)=G_Td(i) D_Td(i)+I_Td(i); %Constructing the source/sink terms S(i)=dG_Td(i)+dD_Td(i)+dI_Td(i); %Constructing the slope of source/sink terms end for i=1:n if i==1  i==n %Inputing the Boundary Conditions if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2= ((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)= 0; A(i,2)=a2+S(i); A(i,3)= a2; A(i,4)= D(i)+S(i)*Td(i); else PhD_2_BC %Wall Boundary Condition end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/ 2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i) rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)= a1; A(i,2)=a1+a2+S(i); A(i,3)= a2; A(i,4)= D(i)+S(i)*Td(i); end end %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values Td=X; end PhD_2_k_E_Eqn.m %Solving the k E equations for Two Phase flow %Solving the k equation variable= 5; for i=1:n C(i)=(1 nus(i))*(muefd(i)+muTd(i)/sigmak); %Constructing the diffusion term D(i)=G_kd(i)D_kd(i)+I_kd(i); %Constructing the source/sink terms S(i)=dG_kd(i)+dD_kd(i)+dI_kd(i); %Constructing the slope of source/sink terms end %Constr ucting the matrix, using principals of Finite Volume %The matrix A is actually a tridiagonal matrix but is constructed in form %of 4 columnar arrys. The first column is the sub diagonal, the second %column is a daigonal and the third is the super diagonal while the forth %column is the B matrix in AX=B where A is the orignal tri diagonal for i=1:n if i==1  i==n %Inputing the Boundary Conditions if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; PAGE 209 209 a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)=0; A(i,2)=a2 S(i); A(i,3)= a2; A(i,4)=D(i) S(i)*kd(i); else %Wall Boundary condition Ph D_2_BC end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; A(i,1)= a1; A(i,2)=a1+a2 S(i); A(i,3)= a2; A(i,4)=D(i)S(i)*kd(i); end end %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values kd=X; %Solving the E equation variable=6; for i=1:n C(i)=(1 nus(i))*(mu efd(i)+muTd(i)/sigmaE); %Constructing the diffusion term D(i)=G_Ed(i)D_Ed(i)+I_Ed(i); %Constructing the source/sink terms S(i)=dG_Ed(i)+dD_Ed(i)+dI_Ed(i); %Constructing the slope of source/sink terms end %Constructing the matrix, using principals of Finite Volume %The matrix A is actually a tridiagonal matrix but is constructed in form %of 4 columnar arrys. The first column is the sub diagonal, the second %column is a daigonal and the third is the super diagonal, while the forth %column is the B matrix in AX=B where A is the orignal tri diagonal for i=1:n if i==1  i==n %Inputing the Boundary Conditions if i==1 %Center Line Boundary Condition of Symmetry rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+ 1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; A(i,1)=0; A(i,2)=a2 S(i); A(i,3)= a2; A(i,4)=D(i) S(i)*Ed(i); else PhD_2_BC %Wall Boundary condition end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; A(i,1)= a1; A(i,2)=a1+a 2 S(i); A(i,3)= a2; A(i,4)=D(i)S(i)*Ed(i); end end PAGE 210 210 %Solving the tri diagonal system PhD_Thomas_Algorithims %Updating values Ed =X; PhD_2_BC.m %Boundary Conditions %Boundary Condition for Vfz if variable==2 %No slip if Vfz_no_sl ip==yes A(i,1)=0; A(i,2)=1; A(i,3)=0; A(i,4)=0; B2(i)=0; end %Wall Function from Mfix if Vfz_wall_Mfix==yes h1=rd(i) rd(i 1); h2=cmu^0.25*kd(i)^0.5*h1/2; rd2=((rd(i)^2) ((rd(i)+rd( i1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1 S(i)+rd(i)/rd2*0.42*cmu^0.25*kd(i)^0.5/log(9.81*h2); A(i,3)=0; A(i,4)=D(i)S(i)*Vfzd(i); B2(i)=1; end end %Boundary Condition for Vsz if variable==3 %No slip if Vsz_no_slip==yes A(i,1)=0; A(i,2)=1; A(i,3)=0; A(i,4)=0; end if Johnson_Jackson_Vsz==yes rd2=((rd(i)^2) ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1+S(i) (rd(i)/rd2*pi*(rhos/rhof)*fi*Td(i)^0.5/(2*3^0.5)/(nus0/nus(i) (nus0/nus(i))^(2/3))); A(i,3)=0; A(i,4)= D(i)+S(i)*Vszd (i); ADiff(i,i 1)= a1; ADiff(i,i)=a1 (rd(i)/rd2*pi*(rhos/rhof)*fi*Td(i)^0.5/(2*3^0.5)/(nus0/nus(i) (nus0/nus(i))^(2/3))); end %Wall Boundary condition: from MFIX, Jenkins condition if Vsz_wall_Mfix==yes g0=(nus0^(1/3))/( nus0^(1/3) nus(i)^(1/3)); rd2=((rd(i)^2) ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1+S(i); PAGE 211 211 A(i,3)=0; A(i,4)= D(i)+S(i)*Vszd(i)+(rd (i)/rd2*nus(i)*rhos/rhof*Td(i)*(1+2*(1+ew)*nus(i)*g0)*tan(fiw)); ADiff(i,i 1)= a1; ADiff(i,i)=a1; end end %Boundary Condition for T if variable==4 if Johnson_Jackson_T==yes rd2=((rd(i)^2) ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1+S(i) (rd(i)/rd2*3^0.5*pi*(rhos/rhof)*(1 ew^2)*Td(i)^0.5/4/(nus0/nus(i) (nus0/nus(i))^(2/3))); A(i,3)=0; A(i,4)= D(i)+S(i)*Td(i)(rd(i)/rd2*pi*(rhos/rhof)*fi*(Vszd(i))^2*Td(i)^0.5/(2*3^0.5)/(nus0/nus(i) (nus0/nus(i))^(2/3))); ADiff(i,i 1)= a1; ADiff(i,i)=a1 (rd(i)/rd2*3^0.5*pi*(rhos/rhof)*(1ew^2)*Td(i)^0.5/4/(nus0/nus(i) (nus0/nus(i))^(2/3))); end %Wall Boundary condition: from MFIX, Jenkins condition if Td_wall_Mfix==yes rd2=((rd(i)^2) ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1+S(i) (rd(i)/rd2*nus (i)*rhos/rhof*(1+2*(1+ew)*nus(i)*g0)*(3*Td(i))^0.5*3/8*(7/2*(1+ew)*tan(fiw)^2 (1 ew))); A(i,3)=0; A(i,4)= D(i)+S(i)*Td(i); ADiff(i,i 1)= a1; ADiff(i,i)=a1; end end %Boundary Condition for k if variable==5 %No slip if k_no_slip==yes A(i,1)=0; A(i,2)=1; A(i,3)=0; A(i,4)=0; end %Wall Function from Mfix if k_wall_Mfix==yes h1=rd(i) rd(i 1); h2=cmu^0.25*kd(i)^0.5*h1/2; D(i)=(2*(1 nus(i))*cmu^0.5*kd(i)* Vfzd(i)/h1/log(9.81*h2)) (1 nus(i))*(2*cmu^0.75*kd(i)^1.5/0.42/h1)+I_kd(i); S(i)=wk(1)*(2*(1 nus(i))*cmu^0.5*Vfzd(i)/h1/log(9.81*h2)) 10*1.5*(1 nus(i))*(2*cmu^0.75*kd(i)^0.5/0.42/h1); rd2=((rd(i)^2) ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; A(i,1)= a1; A(i,2)=a1 S(i); A(i,3)=0; A(i,4)=D(i)S(i)*kd(i); end PAGE 212 212 end %Boundary Condition for E if variable==6 %E at wall=muef/rhof*d^2k/dr^2 if Bolio_equation_2_62==yes A(i,1)=0; A(i,2)=1; A(i,3)=0; h1=drb(n); h2=drb(n)+drb(n 1); A(i,4)=muefd(i)*((kd(i) kd(i 1))/h1 (kd(i)kd(i 2))/h2)*(2/(h2 h1))+I_kd(i)/(1 nus(i)); end %Wall Function from Mf ix if E_wall_Mfix==yes h1=rd(i) rd(i 1); A(i,1)=0; A(i,2)=1; A(i,3)=0; A(i,4)=2*cmu^0.75*kd(i)^1.5/0.42/h1; end end PhD_1_Convergence.m %Calculating convergence for single phase Vfzd_C(z)=norm(Vfzd Vfzd_o)/n orm(Vfzd); kd_C(z)=norm(kd kd_o)/norm(kd); Ed_C(z)=norm(EdEd_o)/norm(Ed); PhD_2_Convergence.m %Calculating convergence for two phase nus_C(z)=norm(nus nus_o)/norm(nus); Vfzd_C(z)=norm(Vfzd Vfzd_o)/norm(Vfzd); Vszd_C(z)=norm(Vszd Vszd_o)/norm(Vszd); Td_C( z)=norm(Td Td_o)/norm(Td); kd_C(z)=norm(kd kd_o)/norm(kd); Ed_C(z)=norm(EdEd_o)/norm(Ed); kTd_C(z)=norm(kTd kTd_o)/norm(kTd); PhD_Re_March.m %Iterating on Re for single phase n1=n_Re_iterations; for z=1:n1 Z(z)=z; %Saving old profiles Vfzd_o =Vfzd; kd_o=kd; Ed_o=Ed; PhD_1_Strain %Calculating strain PhD_1_Viscosity %Calculating eddy viscosity PhD_1_Vfz_Eqn %Solving Vfz equation PhD_1_k_E_Eqn %Solving k E equation PhD_1_Convergence %Calculating convergence %Chec king if anything is negative negative_profiles=no; for i=1:n PAGE 213 213 if Vfzd(i)<0 disp('Warning: negative Vfz') negative_profiles=yes; z_negative_profiles=z break end if kd(i)<0 disp('Warning: negative k') negative_profiles=yes; z_negative_profiles=z break end if Ed(i)<0 disp('Warning: negative E') negative_profiles=yes; z_negative_pr ofiles=z break end end if negative_profiles==yes break end %Testing if new solution has converged if Vfzd_C(z)< Vfzd_tolerance if kd_C(z) PAGE 214 214 kd_o=kd; Ed_o=Ed; kTd_o=kTd; PhD_Relative_Velocity PhD_2_Strain %Calculating Strain PhD_2_Viscosity %Calculating all the various viscosities PhD_Drag_Models %Calculating Drag Ph D_2_Generation %Calculating all the various Generation terms PhD_2_Dissipation %Calculating all the various Dissipation terms PhD_2_Interaction %Calculating all the various Interaction terms PhD_2_nus_Eqn %Solving the nus equation PhD_2_Vfz _Eqn %Solving the Vfz equation PhD_2_Vsz_Eqn %Solving the Vsz equation PhD_2_T_Eqn %Solving the T equation PhD_2_k_E_Eqn %Solving the k E equation PhD_2_Convergence %Calculating convergence %Checking if anything is negative negative_profiles=no; for i=1:n if nus(i)<0 disp('Warning: negative nus') negative_profiles=yes; z_negative_profiles=z break end if Vfzd(i)<0 disp('Warning: negative Vfz') negative_profiles=yes; z_negative_profiles=z break end if Vszd(i)<0 disp('Warning: negative Vsz') negative_profiles=yes; z_negative_profiles=z break end if Td(i)<0 disp('Warning: negative T') negative_profiles=yes; z_negative_profiles=z break end if kd(i)<0 disp('Warning: negative k') negative_profiles=yes; z_negative_profiles=z break end if Ed(i)<0 disp('Warning: negative E') negative_profiles=yes; z_negative_profiles=z break end end PAGE 215 215 if negative_profiles==yes break end %Saying Work in Progress if z2==z1 t2=cputime; disp('*******************************************************') disp('Marching on m') disp('Solving for...') Re_achieved=Re m_achi eved=m nu_C=nus_C(z) Vfz_C=Vfzd_C(z) Vsz_C=Vszd_C(z) T_C=Td_C(z) k_C=kd_C(z) E_C=Ed_C(z) kT_C=kTd_C(z) Time=t2 t1 Iteration=z end %Testing for convergence if nus_C(z) PAGE 216 216 else m_march_complete=no; end break end end end end end end end end PhD_Single_Phase.m %Single Phase Code %The Single Phase code is used to get Vfz, k & E at the exact operating Re. %These profiles will act as guesses for the Two Phase profiles. %Selecting a Guess for Vfz, k & E based on the operating Re Re_target=Re; PhD_Guess_Profiles % if User had specified grid, adjusting the guess grid to match the specified grid if User_Defined_Grid==yes PhD_Regriding end % Pre Defining matrices for Matlab optimization r=zeros(n,1); rd=zeros(n,1); drf=zeros(n,1); drb=zeros(n,1); Vfzd_o=zeros(n,1); kd_o=zeros(n,1); Ed_o=zeros(n,1); dVfzdr=zeros(n,1); fT2=zeros(n,1); muTd=zeros(n,1); dmuTdkd=zeros(n,1); dmuTdEd=zeros(n,1); C=ze ros(n,1); D=zeros(n,1); S=zeros(n,1); A=zeros(n,4); X=zeros(n,1); Solution_Single_Phase=zeros(n,1); r=R*rdim; %Defining r coordinate %Res is the guess Re %Begining of the marching from Res towards Re_target if Re_target PAGE 217 217 drf(i)=rd(i+1)rd(i); %Getting foward difference drb(j)=rd(j) rd(j 1); %Getting backward difference end PhD_Re_March %iterating on Re if negative_profiles==yes disp('Terminating Code') break end if Re==Re_target break end end else %Similar as above for q=1:n_Re_march Re=Res* exp(Re_jump*q); if Re>Re_target Re=Re_target; end r=R*rdim; V=Re*muf/rhof/2/R; rd=Re/2*rdim; for i=1:n 1 j=i+1; drf(i)=rd(i+1)rd(i); drb(j)=rd(j) r d(j 1); end PhD_Re_March if negative_profiles==yes disp('Terminating Code') break end if Re==Re_target break end end end if negative_profiles==yes break end for j=1 :4 for i=1:n if j==1 Solution_Single_Phase(i,j)=r(i); end if j==2 Solution_Single_Phase(i,j)=Vfzd(i)*V; end if j==3 Solution_Single_Phase(i,j)=kd(i)*V^2; end if j==4 Solution_Single_Phase(i,j)=Ed(i)*V^4*rhof/muf; end end end PAGE 218 218 PhD_Definitions.m %Pre defining orders of matrices for Matlab optimization %List is too long to be displayed PhD_Two_Phase %Main Code clear all clc format compact t1 =cputime; yes=1; no=0; PhD_Data %Taking in data from user PhD_Model_Options %Selecting the models PhD_Model_Parameters %Setting parameters and weights PhD_Marching_Parameters %Setting Marching parameters t2=cputime; %Running the Single Phase code disp('*** ****************************************************') disp('Starting to March on Re') Time=t2 t1 disp('*******************************************************') PhD_Single_Phase disp('*******************************************************') disp('Complet ed March on Re') disp('*******************************************************') PhD_Definitions %Pre Defining orders of matrices for Matlab optimization for i=1:n N(i)=i; %Constructing a grid point matrix end r=R*rdim; %Defining the radius of the pipe eta=(1+e)/2; etaf=(1+ef)/2; Ar=d^3*rhof*(rhos rhof)*g/muf^2; dR=d/R; Rep=Re/2*dR; ST=Rep/18*dR/2*rhos/rhof; PhD_ST_Classify m_target=m; M1=R*rhof*m/d/rhos; %Defining M1=R*rhof*m/d/rhos %M1 generally varies form 0 to 2 for most cases. %Initiating M1 at 0. 05 if necessary if M1>m_start M1=m_start; end M=M1*dR; m=M*rhos/rhof; for i=1:n %Guessing Vszd, Td & nus Vszd(i)=Ratio_V*Vfzd(1); nus(i)=M/(Ratio_V+M); Td(i)=kd(1); kTd(i)=(6*kd(i)*Td(i))^0.5; end nus_n=nus; Vfzd_n=Vfzd; PAGE 219 219 Vszd_n=Vszd; Td _n=Td; kd_n=kd; Ed_n=Ed; kTd_n=kTd; t2=cputime; disp('*******************************************************') disp('Starting to March on m') Time=t2 t1 disp('*******************************************************') %Begining of the marching from M1=0.05 towards m_target for q=1:n_m_march PhD_m_March if z==n_m_iterations %Give warning for overloading disp('*******************************************************') disp('Overloading of Solids') disp('*******************************************************') end if m_march_complete==yes %Breaking the loop break end if negative_profiles==yes  z==n_m_iterations disp('Reducing m_jump') m_jump=m_jump/2; if m_jump*d/R*rhos/rhof PAGE 220 220 if negative_profiles==yes  z==n_m_iterations else PhD_Post_Processing %Post Processing save('Output') end PhD_2_nus_Residual.m %Calculating nus Residual variable=1; %For Low Stokes number if Low_ST==yes for i=1:n C(i)=(muefd(i)+m uTd(i))/0.7; end %Constructing the diffusion matrix for i=1:n if i==1  i==n %Inputing Boundary Conditions if i==1 a2= ((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i)); ADiff (i,i)=a2; ADiff(i,i+1)= a2; else a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1)); ADiff(i,i 1)= a1; ADiff(i,i)=a1; end else a1= ((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i) rd(i 1)); a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i)); ADiff(i,i 1)= a1; ADiff(i,i)=a1+a2; ADiff(i,i+1)= a2; end e nd %Calculating the Diffusion term nus_Diff=ADiff*nus; nus_Residual=nus_Diff; end %For Mid and High Stokes number if Mid_ST==yes  High_ST==yes temp=constant; clear constant for i=1:n alphaTd(i)=alpha(i)*Td(i); constant(i)=temp; nus_Residual(i)=alpha(i)*Td(i) constant(i); end end PhD_2_Vfz_Residual.m %Calculating Vfz Residual variable=2; for i=1:n G_Vfzd(i)= dPdz; %letting Pressure Drop to be the generation term PAGE 221 221 C(i)=muefd(i)+muTd(i); %Calculating the Diffus ion co efficient D(i)=G_Vfzd(i) D_Vfzd(i)+I_Vfzd(i); %Compiling the source/sink terms end D(n)=0; %Constructing the Diffusion matrix for i=1:n if i==1 i==n %Inputing the Boundary Conditions if i==1 rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; ADiff(i,i)=a2; ADiff(i,i+1)= a2; end if i==n ADiff(i,i)=1; end else rd2=(((r d(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; ADiff(i,i 1)= a1; ADiff(i,i)=a1 +a2; ADiff(i,i+1)= a2; end end Vfz_Diff=ADiff*Vfzd; %Calculating the Diffusion term Vfz_Residual=Vfz_Diff D; %Calculating the Residual PhD_2_Vsz_Residual.m %Calculating Vfz Residual variable=3; for i=1:n C(i)= musdG(i); %Calculating the Di ffusion co efficient D(i)=G_Vszd(i)D_Vszd(i)+I_Vszd(i); %Compiling the source/sink terms end %Constructing the Diffusion matrix for i=1:n if i==1  i==n %Inputing Boundary Conditions if i==1 rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; ADiff(i,i)=a2; ADiff(i,i+1)= a2; else PhD_2_BC end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1) )/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; ADiff(i,i 1)= a1; ADiff(i,i)=a1+a2; ADiff(i,i+1)= a2; end end Vsz_Diff=ADiff*Vszd; %Calculating the Diffusion term PAGE 222 222 Vsz_Residual=Vsz_Diff+D; %Calculating the Residual PhD_2_T_Residual.m %Calculating T Residual varibale=4; if Bolio_equation__36==yes rd2=((rd(n)^2)((rd(n)+rd(n 1))/2)^2)/2; %Adding extra generation due to the wall G_Td(n)= G_Td(n)+(rd(n)/rd2*pi*(rhos/rhof)*fi*Vszd(n)^2*Td(n)^0.5/(2*3^0.5)/(nus0/nus(n) (nus0/nus(n))^(2/3))); end for i=1:n C(i)= lambdadG(i); %Calculating the Diffusion co efficient D(i)=G_Td(i)D_Td(i)+I_Td(i); %C ompiling the source/sink terms end %Constructing the Diffusion matrix for i=1:n if i==1  i==n %Inputing Boundary Conditions if i==1 rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i))/rd2; ADiff(i,i)=a2; ADiff(i,i+1)= a2; else PhD_2_BC end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; ADiff(i,i 1)= a1; ADiff(i,i)=a1+a2; ADiff(i,i+1)= a2; end end T_Diff=ADiff*Td; %Calculating the Diffusion term T_Residual=T_Diff+D; %Calculating the Residual PhD_2_k_E_Residual.m %Calculating k E Residual %k Residual variable=5; for i=1:n C(i)=(1 nus(i))*(muefd(i)+muTd(i)/sigmak); %Calculating the Diffusion coefficient D(i)=G_kd(i)D_kd(i)+I_kd(i); %Compiling the source/sink terms end D(n)=0; %Constructing the Diffusion matrix for i=1:n if i==1  i==n %Inputing Boundary Conditions if i==1 rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+ 1)))/(rd(i+1) rd(i))/rd2; ADiff(i,i)=a2; ADiff(i,i+1)= a2; else PAGE 223 223 ADiff(i,i 1)=0; ADiff(i,i)=1; end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; ADiff(i,i 1)= a1; ADiff(i,i)=a1+a2; ADiff(i,i+1)= a2; end end k_Diff=ADiff*kd; %Calculat ing the Diffusion term k_Residual=k_Diff D; %Calculating the Residual %E Residual variable=6; for i=1:n C(i)=(1 nus(i))*(muefd(i)+muTd(i)/sigmaE); %Calculating the Diffusion coefficient D(i)=G_Ed(i)D_Ed(i)+I_Ed(i); %Compiling the source/sink term s end D(n)=0; %Constructing the Diffusion matrix for i=1:n if i==1  i==n %Inputing Boundary Conditions if i==1 rd2=(((rd(i)+rd(i+1))/2)^2 (rd(i)^2))/2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1) rd(i ))/rd2; ADiff(i,i)=a2; ADiff(i,i+1)= a2; else ADiff(i,i 1)=0; ADiff(i,i)=1; end else rd2=(((rd(i)+rd(i+1))/2)^2 ((rd(i)+rd(i 1))/2)^2)/2; a1=((rd(i)+rd(i 1))/2)*(2*C(i)*C(i 1)/(C(i)+C(i 1)))/(rd(i)rd(i 1))/rd2; a2=((rd(i)+rd(i+1))/2)*(2*C(i)*C(i+1)/(C(i)+C(i+1)))/(rd(i+1)rd(i))/rd2; ADiff(i,i 1)= a1; ADiff(i,i)=a1+a2; ADiff(i,i+1)= a2; end end E_Diff=ADiff*Ed; %Calculating the Diffusion co efficient E_Residual=E_Diff D; %Calculating the Residual PhD_Dimensions.m %Giving dimensions oa all variables %Giving dimensions to solution variables Vfz=Vfzd*V; Vsz=Vszd*V; T=Td*V^2; k=kd*V^2; E=rhof*V^4*Ed/muf; UT=UTd*V; %Giving dimensions to all the terms in the equation if Low_ST==yes PAGE 224 224 nus_Residual=V*rhof*nus_Residual; Vfz_Diff=(V^2*rhof^2/muf)*V*Vfz_Diff; G_Vfz=(V^2*rhof^2/muf)*V*G_Vfzd; D_Vfz=(V^2*rhof^2/muf)*V*D_Vfzd; I_Vfz=(V^2*rhof^2/muf)*V*I_Vfzd; Vfz_Residual=(V^2*rhof ^2/muf)*V*Vfz_Residual; Vsz_Diff=(V^2*rhof^2/muf)*V*Vsz_Diff; G_Vsz=(V^2*rhof^2/muf)*V*G_Vszd; D_Vsz=(V^2*rhof^2/muf)*V*D_Vszd; I_Vsz=(V^2*rhof^2/muf)*V*I_Vszd; Vsz_Residual=(V^2*rhof^2/muf)*V*Vsz_Residual; k_Diff=(V^2*rhof^2/muf)*V^2*k_Diff; G_k=(V^2*rhof^2/muf)*V^2*G_kd; D_k=(V^2*rhof^2/muf)*V^2*D_kd; I_k=(V^2*rhof^2/muf)*V^2*I_kd; k_Residual=(V^2*rhof^2/muf)*V^2*k_Residual; E_Diff=(V^2*rhof^2/muf)*rhof/muf*V^4*E_Diff; G_E=(V^2*rhof^2/muf)*rhof/muf*V^4*G_Ed; D_E=(V^2*rhof^2/muf)*rhof/muf*V^4*D_Ed; I_E=(V^2*rhof^2/muf)*rhof/muf*V^4*I_Ed; E_Residual=(V^2*rhof^2/muf)*rhof/muf*V^4*E_Residual; end if Mid_ST==yes  High_ST==yes alphaT=V^2*alphaTd; constant=V^2*constant; nus_Residual=V^2*nus _Residual; Vfz_Diff=(V^2*rhof^2/muf)*V*Vfz_Diff; G_Vfz=(V^2*rhof^2/muf)*V*G_Vfzd; D_Vfz=(V^2*rhof^2/muf)*V*D_Vfzd; I_Vfz=(V^2*rhof^2/muf)*V*I_Vfzd; Vfz_Residual=(V^2*rhof^2/muf)*V*Vfz_Residual; Vsz_Diff=(V^2*rhof^2/muf)*V*Vsz_Diff; G_Vsz=(V^2*rhof^2/muf)*V*G_Vszd; D_Vsz=(V^2*rhof^2/muf)*V*D_Vszd; I_Vsz=(V^2*rhof^2/muf)*V*I_Vszd; Vsz_Residual=(V^2*rhof^2/muf)*V*Vsz_Residual; T_Diff=(V^2*rhof^2/muf)*V^2*T_Diff; G_T=(V^2*rhof^2/muf)*V^2*G_Td; D_T=(V^2*rhof^2/ muf)*V^2*D_Td; I_T=(V^2*rhof^2/muf)*V^2*I_Td; T_Residual=(V^2*rhof^2/muf)*V^2*T_Residual; k_Diff=(V^2*rhof^2/muf)*V^2*k_Diff; G_k=(V^2*rhof^2/muf)*V^2*G_kd; D_k=(V^2*rhof^2/muf)*V^2*D_kd; I_k=(V^2*rhof^2/muf)*V^2*I_kd; k_Residual=(V^2*rhof^2/muf)*V^2*k_Residual; E_Diff=(V^2*rhof^2/muf)*rhof/muf*V^4*E_Diff; G_E=(V^2*rhof^2/muf)*rhof/muf*V^4*G_Ed; D_E=(V^2*rhof^2/muf)*rhof/muf*V^4*D_Ed; I_E=(V^2*rhof^2/muf)*rhof/muf*V^4*I_Ed; E_Residual=(V^2*rhof^2/muf)*rhof/muf *V^4*E_Residual; end PhD_Post_Processing.m %Post Processing disp('Post Processing...') PAGE 225 225 disp('*******************************************************') ADiff=zeros(n,n); PhD_2_Strain %Calculating Strain PhD_2_Viscosity %Calculating various viscosities PhD_ 2_Generation %Calculating Generation terms PhD_2_Dissipation %Calculating Dissipation terms PhD_2_Interaction %Calculating Interaction terms %Calculating Residulas PhD_2_nus_Residual PhD_2_Vfz_Residual PhD_2_Vsz_Residual PhD_2_T_Residual PhD_2_k_E_Residual %Giving each term appropriate dimensions PhD_Dimensions %Calculating Pressure Drop and Friction velocity if downward==yes %Vfz= Vfz; %Vsz= Vsz; %Solution_Single_Phase(:,2)= Solution_Single_Phase(:,2); dPdz= (V^2*rhof^2/muf)*V*dPdz rhof*g; rhoavg=rhof*rhos*(m+1)/(rhos+m*rhof); Friction_Velocity=((dPdz+rhoavg*g)*R/2/rhoavg)^0.5; else dPdz=(V^2*rhof^2/muf)*V*dPdz rhof*g; rhoavg=rhof*rhos*(m+1)/(rhos+m*rhof); Friction_Velocity=(( dPdz rhoavg*g)*R/2/rhoavg)^0.5; end %Calculat ing average Velocity and average nus V_avg=0; nus_avg=0; for i=1:n 1 V_avg=V_avg+(r(i)*(1nus(i))*Vfz(i)+r(i+1)*(1 nus(i+1))*Vfz(i+1))*(r(i+1) r(i))/R^2; nus_avg=nus_avg+(r(i)*nus(i)+r(i+1)*nus(i+1))*(r(i+1) r(i))/R^2; end %Calculating various time scales cbeta=1.8 1.35; for i=1:n 1 g0=(nus0^(1/3))/(nus0^(1/3) nus(i)^(1/3)); zetar=3*(Vrd(i))^2/2/kd(i); time_collision(i)=d/24/nus(i)/g0*(pi/T(i))^0.5; time_turbulent_eddy(i)=3/2*cmu*k(i)/E(i); time_fluid_integral(i)=time_turbulent_e ddy(i)/(1+cbeta*zetar)^0.5; time_particle_relax(i)=nus(i)*rhos/(1 nus(i))/beta(i)/V^2/rhof^2*muf; ratio_c_r(i,1)=time_collision(i)/time_particle_relax(i); ratio_r_t(i,1)=time_particle_relax(i)/time_turbulent_eddy(i); ratio_c_t(i,1)=time_col lision(i)/time_turbulent_eddy(i); end for j=1:7 %Compiling the Solution matrix for i=1:n if j==1 Solution_Two_Phase(i,j)=r(i); end if j==2 Solution_Two_Phase(i,j)=nus(i); end PAGE 226 226 if j==3 Solution_Two_Phase(i,j)=Vfz(i); end if j==4 Solution_Two_Phase(i,j)=Vsz(i); end if j==5 Solution_Two_Phase(i,j)=T(i); end if j==6 Solution_Two_Phase(i,j)=k(i); end if j==7 Solution_Two_Phase(i,j)=E(i); end end end PAGE 227 227 APPENDIX B RAW FLUIDSOLID DATA Table B 1. Raw data for d = 1.5 mm, m = 0.0175 and Vfzcl = 3.07 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 1.86 0.357 2.21 0. 405 0.93 2.15 0.330 2.28 0.426 0.89 2.30 0.291 2.31 0.434 0.85 2.48 0.271 2.36 0.449 0.82 2.50 0.285 2.40 0.455 0.64 2.76 0.228 2.67 0.413 0.45 2.94 0.209 2.83 0.395 0.27 3.02 0.199 2.92 0.361 0.09 3.07 0.192 3.00 0.337 0.00 3.12 0.182 3.00 0.347 Table B 2. Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 1 r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 1.80 0.326 2.14 0.380 0.93 2.16 0.295 2.31 0.382 0.89 2.27 0.289 2.35 0.379 0.85 2.39 0.273 2.40 0.413 0.82 2.49 0.261 2.46 0.41 9 0.64 2.78 0.226 2.67 0.381 0.45 2.83 0.211 2.80 0.354 0.27 2.92 0.197 2.93 0.335 0.09 2.96 0.187 2.95 0.338 0.00 2.96 0.183 2.95 0.336 PAGE 228 228 Table B 3. Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 2 r/R V sz v sz m/s m/s 0.96 2.19 0.372 0.93 2.29 0.377 0.89 2.37 0.374 0.85 2.43 0.398 0.82 2.48 0.401 0.64 2.65 0.410 0.45 2.73 0.403 0.27 2.88 0.362 0.09 2.96 0.351 0.00 2.93 0.341 Table B 4. Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 3 r/R V sz v sz m/s m/s 0.96 2.19 0.470 0.93 2.24 0.455 0.89 2.29 0.459 0.85 2.33 0.430 0.82 2.39 0.409 0.64 2.67 0.392 0.45 2.80 0.357 0.27 2.90 0.322 0.09 2.95 0.331 0.00 2.94 0.324 PAGE 229 229 Table B 5. Raw data for d = 1.5 mm, m = 0.0425, Vfzcl = 2.96 m/s and run 4 r/R V sz v sz m/s m/s 0.96 2.16 0.445 0.93 2.20 0.435 0.89 2.26 0.449 0.85 2.31 0.434 0.82 2.29 0.437 0.64 2.70 0.380 0.45 2.80 0.372 0.27 2.91 0.332 0.09 2.94 0.328 0.00 2.97 0.326 Table B 6. Raw data for d = 1.5 mm, m = 0.075 and Vfzcl = 3.04 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 1.86 0.322 2.30 0.390 0.93 2.26 0.285 2.38 0.414 0.89 2.41 0.284 2.40 0.412 0.85 2.51 0.268 2.45 0.426 0.82 2.60 0.256 2.51 0.441 0.64 2.78 0.224 2.54 0.466 0.45 2.92 0.211 2.87 0.393 0.27 2.96 0.20 1 2.90 0.364 0.09 3.00 0.184 3.04 0.327 0.00 3.04 0.181 3.03 0.339 PAGE 230 230 Table B 7 Raw data for d = 1.5 mm, m = 0.0175 and Vfzcl = 5.03 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 2.99 0.574 3.79 0.512 0.93 3.52 0.455 3.95 0.565 0.89 3.86 0.424 4 .05 0.599 0.85 4.00 0.433 4.11 0.600 0.82 3.77 0.442 4.19 0.637 0.64 4.40 0.348 4.49 0.533 0.45 4.67 0.297 4.59 0.546 0.27 4.90 0.247 4.75 0.527 0.09 5.02 0.222 4.85 0.507 0.00 5.03 0.221 4.87 0.506 Table B 8. Raw data for d = 1.5 mm, m = 0.0425 and Vfzcl = 4.90 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 3.07 0.506 3.75 0.476 0.93 3.50 0.481 3.89 0.506 0.89 3.73 0.434 4.01 0.488 0.85 3.93 0.404 4.07 0.536 0.82 3.99 0.402 4.20 0.541 0.64 4.42 0.322 4.37 0.550 0.45 4.63 0.281 4.83 0.545 0.27 4.78 0.249 4.93 0.492 0.09 4.89 0.231 4.96 0.463 0.00 4.90 0.212 4.97 0.451 PAGE 231 231 Table B 9. Raw data for d = 1.5 mm, m = 0.075 and Vfzcl = 5.04 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 3.21 0.497 3.86 0.503 0.93 3.66 0.456 3.98 0.499 0. 89 3.83 0.419 4.09 0.525 0.85 4.03 0.417 4.17 0.535 0.82 4.16 0.416 4.24 0.557 0.64 4.50 0.321 4.48 0.557 0.45 4.76 0.281 4.82 0.486 0.27 4.92 0.244 4.91 0.513 0.09 5.02 0.234 5.07 0.453 0.00 5.04 0.218 5.09 0.485 Table B 10. Raw data for d = 1.5 mm, m = 0.0175 and Vfzcl = 7 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 4.19 0.680 5.16 0.732 0.93 4.90 0.590 5.42 0.757 0.89 5.40 0.585 5.63 0.758 0.85 5.61 0.549 5.62 0.745 0.82 5.70 0.524 5.75 0.729 0.64 6.21 0.427 6.39 0.621 0.45 6.57 0.3 76 6.71 0.572 0.27 6.79 0.308 0.09 6.97 0.268 0.00 7.00 0.272 PAGE 232 232 Table B 11. Raw data for d = 1.5 mm, m = 0.0425 and Vfzcl = 7.06 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 4.16 0.772 5.26 0.732 0.93 4.97 0.671 5.47 0.712 0.89 5.40 0. 570 5.60 0.687 0.85 5.61 0.553 5.73 0.709 0.82 5.69 0.530 5.87 0.746 0.64 6.33 0.438 6.43 0.672 0.45 6.62 0.382 6.74 0.601 0.27 6.85 0.325 6.98 0.521 0.09 7.04 0.283 6.95 0.513 0.00 7.06 0.263 6.94 0.474 Table B 12. Raw data for d = 1.5 mm, m = 0. 075 and Vfzcl = 7.02 m/s r/R V fz v fz V sz v sz m/s m/s m/s m/s 0.96 4.42 0.639 5.03 0.732 0.93 5.07 0.618 5.32 0.739 0.89 5.43 0.569 5.49 0.737 0.85 5.63 0.552 5.67 0.717 0.82 5.80 0.532 5.84 0.683 0.64 6.27 0.429 6.15 0.775 0.45 6.60 0.365 6.66 0.613 0.27 6.85 0.315 6.88 0.564 0.09 7.00 0.278 6.96 0.532 0.00 7.02 0.273 7.06 0.576 PAGE 233 233 APPENDIX C FIGURES FOR FLUIDIZED BED BINARY MIXTURE S Figure C 1. Pressure drop profiles and segregation index behavior for a Type A mixture 50G55050G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 234 234 Figure C 2. Pressure drop profiles and segregation index behavior for a Type A mixture 50G46250G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 235 235 Figure C 3. Pressure drop profiles and segregation index behavior for a Type A mixture 50G55050G116. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 236 236 Figure C 4. Pressure drop profiles and segregation index behavior for a Type A mixture 50G38550G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 237 237 Figure C 5. Pressure drop profiles and segregation index behavior for a Type A mixture 50G 46250G116. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 238 238 Figure C 6. Pressure drop profiles and segregation index behavior for a Type B mixture 50G32850G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 239 239 Figure C 7. Pressure drop profiles and segregation index behavior for a Type B mix ture 50G27550G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation index. PAGE 240 240 Figure C 8. Pressure drop profiles and segregation index behavior for a Type C mixture 50G23150G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 241 241 Figure C 9. Pressure drop profiles and segregation index behavior for a Type C mixture 50G19550G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation index. PAGE 242 242 Figure C 10. Pressure drop profiles and segregation index behavior for a Type D mixture 50G16550G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 243 243 Figure C 11. Pressure drop profiles and segregation index behavior for a Type D mixture 50G13850G083. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 244 244 Figure C 12. Pressure drop profiles and segregation index behavior for a Type D mixture 50G11650G083. (a) Press ure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 245 245 Figure C 13. Pressure drop profiles and segregation index behavior for a Type B mixture 13S32887P328. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 246 246 Figure C 14. Pressure drop profiles and segregation index behavior for a Type B mixture 75S32825G328. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. PAGE 247 247 Figure C 15. Pressure drop profiles and segregation index behavior for a Type E mixture 70G11630P275. (a) Pressure drop profile, initially segregated state. (b) Pressure drop profile, initially mixed state. (c) Segregation profiles. (d) Segregation index. Equation Chapter (Next) Section 1 PAGE 248 248 LIST OF REF E RENCES 1. Kunii D Levenspiel O Fluidization Engineering, 2nd ed. Newton MA: Butterworth Heinemann 1991. 2. Tsuji Y Morikawa Y Shiomi H LDV measurements of air solid twophase flow in a vertical pipe Journal of Fluid Mechanics 1984; 139: 417 434. 3. Jones E Yurteri C Sinclair J. Effect of mass loading on gas solids motion and particle segregation patterns Handbook of Powder Technology 2001; 10: 843 849. 4. Sheen H C haing Y Chaing Y Two dimensional measurements of flow structure in a two phase vertical pipe flow Proceedings of National Science Council Republic of China (A) 1993; 17: 200 213. 5. Lee S Durst F On the motions of particles in turbulent duct flow Inter national Journal of Multiphase Flow 1982 ; 8 : 125 146. 6. Louge M Mastorakos E Jenkins J. The role of particle collisions in pneumatic transport Journal of Fluid Mechanics 1991; 231: 345 359. 7. Bolio E Yasuna J, Sinclair J. Dilute turbulent gas solid flow in risers with particle particle interaction. AIChE Journal 1995; 41: 1375 1388. 8. Simonin O Continuum modeling of dispersed twophase flow in combustion and turbulence in two phase flow Von Karman Institute of Fluid Dynamics Lecture Series 1996. 9. Benviad es A van Wachen B Numerical simulation and validation of dilute turbulent gas particle flow with inelastic collisions and turbulence modulation Powder Technology 2008; 182: 294 306. 10. Gore R Crowe C Effect of particle size on modulating turbulent intensity International Jornal of Multiphase Flow 1989; 15: 279 285. 11. Hestroni G Particles turbulence interaction. International Jornal of Multiphase Flow 1989 ; 15: 735 746 12. Alajbegovic A Assad A Bonetto F Lahey R Phase distribution and turbulence struct ure for solid/fluid upflow in a pipe. International Journal of Multiphase Flow 1994; 20: 453 479. 13. Pepple M Benchmark data and analysis of dilute turbulent fluidparticle flow in viscous and transitional regimes PhD Thesis University of Florida. 2010. PAGE 249 249 14. H adinoto K Curtis J Numerical simulation of turbulent particleladen flowa with significant fluid to particle inertia ratio. Industrial Engineering Chemistry Research. 2009; 48: 5874 5884. 15. Chen C Wood P Turbulence closure model for dilute gas particle f lows The Canadian Journal of Chemical Engineering. 1985; 63: 349 360. 16. Formisani B Girimonte R Longo T Development of the approach based on the fluidization velocity interval Powder Technology 2008; 185: 97 108. 17. Joseph G Leboreiro J Hrenya C Steven A Experimental segregation profiles in bubbling gas fluidized bed. AIChE Journal 2007; 53: 2804 2813. 18. Marzocchella A Salatino P Di Pastena V Lirer L Transient fluidization and segregation of binary mixtures of particles AIChE Journal 2000; 53: 2175 2182. 19. Rasul M Rudolph V Carsky M Segregation potential in binary gas fluidized beds Powder Technology 1999; 103: 175 181. 20. Potic B Kerstn S Ye M Van der Hoef M van Swaaij W Fluidization of hot compressed water in micro reactors Chemical Engineering Science 2005; 60: 5982 5990 21. Chiba S Nienow A Chiba T Kobayashi H Fluidized binary mixtures in which denser component may be floatasam Powder Technology 1980; 26: 1 10. 22. Garcia F Romero A Villar J, Bello A A study of segregation in a gas soli d fluidized bed: Particles of different density Powder Technology 1989; 58: 169 174. 23. Huilin L Yurong D Gidaspow D Lidan Y Yukun Q Size segregation of binary mixture of solid in bubbling fluidized bed. Powder Technology 2003; 134: 86 97. 24. Geldart D Baeyens J, Pope D Van De Wijer P Segregation in beds of large particles at high velocity Powder Technology 1981; 30: 195 205. 25. Liu X Xu G Gao S Micro fluidized beds: Wall effects and operability Chemical Engineering Journal 2008; 137: 302 307. 26. Babu S Shah B Talwalker A Fluidization correlation for coal gasification materials minimum fluidization velocity and fluidized bed expansion. AIChE Symposium Series 1978; 76: 176 186. 27. Baerg A Klassen J, Gishler P Heat transfer in a fluidized solid bed. Can J Res Section F. 1950; 28: 287 306. PAGE 250 250 28. Broadhurst T Becker H Onset of fluidization and slugging in beds of uniform particles AIChE Journal 1975; 21: 238 247. 29. Davies L Richardson J. Gas interchange between bubbles and the continuous phase in a fluidi zed bed. Trans Inst Chem Eng. 1966; 44: 293 305. 30. Goroshko V Rozenbaum R Todes O Hydrodynamics and Heat Transfer in F luidized B eds Cambridge MA : M I T Press 1966. 31. Grewal N Saxena S Comparision of commonly used correlations for minimum fluidization velocityof small solid particles Powder technology 1980; 26: 229 234. 32. Johnson E Institute of Gas Engineers Report Publication No 318, London, 1949. 33. Kumar P Sen Gupta P Prediction of minimum fluidization velocity for multicomponent mixtures ( short commu nication) Indian Journal Technology 1974; 12 (5) : 225 251. 34. Leva M Fluidization. New York NY: McGraw Hill, 1959. 35. Miller C Longwinuk A Fluidization studies of solid particles Industrial and Engineering Chemistry 1951; 43: 1220 1226. 36. Pillai B Raja Rao M Pressure Drop and minimum fluidization velocities of air fluidized bed. Indian Journal Technology 1971; 9 : 77 89. 37. Saxena S Vogel G The measurement of incipient fluidization velocities in a bed of coarse dolomite at temperature and pressure. Trans I nst Chem Eng. 1977; 55: 184 189. 38. van Heerdeen C Nobel A Krevelen V Studies on fluidization ll heat transfer Chemical Engineering Science. 1951; 1 (2) : 51 66. 39. Wen C Yu Y Mechanics of fluidization. Chemical Engineering Progress Symposium Series 1966; 62: 100 111. 40. Koch D Kinetic theory for a monodispersed gas solid suspension. Physics of fluids 1990 ; 2 : 1711 1722. 41. Igci Y Andrew IV A Sundaresan S Pannala S O'Brien T Filtered twofluid models for fluidized gas particle suspensions AIChE Journal 2008 ; 54: 1431 1448. 42. Koch D Sangani A Particle pressure and marginal stability limits for homogeneous monodisperse gas fluidized bed: kinetic theory and numerical simulations Journal of Fluid Mechanics 1999; 400: 229 263. PAGE 251 251 43. Wylie J, Koch D Ladd A Rhel ogy of suspension with high particle inertia and moderate fluid inertia. Journal of Fluid Mechanics 2003 ; 480: 95 118 44. Sinclair J, Mallo T Describing particle turbulence interaction in a twofluid modeling framework ASME Fluids Engineering Division Summ er Meeting. New York NY: Washington DC: ASME press 1998; 4 : 7 14 45. Crowe C On models for turbulence modulation in fluidparticle flows International Journal of Multiphase Flow 2000 ; 26: 710 727. 46. Zhang Y Reese J. Gas turbulence modulation in a twofluid model for gas solid flows AIChE Journal 2003; 49: 3048 3065. 47. Hill R Koch D Ladd J. ModerateReynolds number flows in ordered and random arrays of spheres Journal of Fluid Mechanics 2001; 448 : 243 278. 48. Hill R Koch D Ladd J. The first effects of flu id inertia on flows in ordered and random arrays of spheres Journal of Fluid Mechanics 2001; 448: 213 241. 49. Benyahia S Syamlal M O'Brien T Summary of MFIX equations 2007, July Retrieved from www mfix org: https://www.mfix.org/documentation/MFIXEquations20054 4.pdf 50. Lun C Savage S Jeffery D Chepurniy N Kinetic theories for granular flow: inelastic particles in couette flow and slightly inelastic particles in general flow field Journal of Fluid Mechanics 1984 ; 140: 223 256. 51. Peirano E Leckner B Fundamentals of turbulent gas solid applied to circulating fluidized bed combustion. Progress in Energy and Combution Sciences 1998; 24: 259 296. 52. Batchelor G Green J The determ ination of bulk stress in a suspension of spherical particle to order of c2. Journal of Fluid Mechanics 1972; 56: 401 427. 53. Lun C Numerical simulation of dilute turbulent gas solid flows International Journal of Multiphase Flow 2000 ; 26: 1707 1736. 54. Beny ahia S Syamlal M O'Brien T Study of the ability of multiphase contimuum models to predict coreannulus flow AIChE Journal 2007; 53: 2549 2569. 55. Johnson P Jackson R Frictional collisional constitutive relations for granular materials with applications to plane shearing. Journal of Fluid Mechanics 1987 ; 176: 67 93. 56. Patankar S Spadling D Heat and Mass Transfer in Boundary Layers London: MorganGrampain Press 1967. PAGE 252 252 57. Jones E An experimental investigation of particle size distribution effect in dilute phase gas solid flow. PhD Thesis Perdue University 2001. 58. Hrenya C, Bolio E, Chakrabarti D, Sinclair J. Comparison of low Reynolds k turbulence model in predicting fully developed pipe flow. Chemical Engineering Science 1995; 50: 1923 1941. 59. Pepple M C urtis J, Yurteri C Variation in mean turbulence intensity in fully developed pipe flow International Review of Chemical Engineering. 2010; 2 : 337 342. 60. Gondert P Lance M Petit L Bouncing motion of spherical particles in fluids Physics of Fluids 2002; 14: 643 652. 61. Hadinoto K Jones E Yurteri C Curtis J. Reynolds number dependence of gas phase turbulence in gas particle f lows International Journal of Multiphase Flow 2005; 31: 416 434. 62. Rowe P Nienow A Agbim A Preliminary quantitative study of part icles segregation in gas fluidised beds binary systems of near spherical particles Trans Inst Chem Eng. 1972; 50: 325 333. 63. Nienow A Rowe P Cheung L A quantitative analysis of the mixing of two segregation powders of different densities in a gas fluidized bed. Powder Technology 1978; 20: 89 97. 64. Nienow A Naimer N Chiba T Studies of segregation/mixing in fluidised beds of different size particles Chem Eng Comm 1987; 62: 53 66. 65. Hoffmann A Janssen L Prins J. Particle segregation in fluidised binary mixtures Chemical Engineering Science. 1993; 48: 1583 1592. 66. Olivieri G Marzocchella A Salatino P Segregation of fluidized binar mixtures of granular solids AIChE Journal 2004; 50: 3095 3106. 67. Formisani B De Cristofaro G Girimonte R A fundamental approach to the phenomenology of fluidization of size segregating binary mixtures of solids Chemical Engineering Science. 2001; 56: 109 119. 68. Hedden D Brone D Clement S McCall M Olsofsy A Patel P Prescott J, Hancock B. Development of an improved fluidization segregation tester for use with pharmaceutical powders Pharmaceutical Technology 2006; 30: 54 64. 69. ASTM International Standard practice for measuring fluidization segregation tendencies of powders ASTM International D6941 03. PAGE 253 253 70. Di Felice R Gibilaro L Wall effects for the pressure drop in fixed beds Chemical Engineering Science. 2004, 59: 3037 3040. 71. Delebarre A Does the Minimum Fluidization Exist? Journal of Fluids Engineering 2002; 124: 595 600. 72. Persson B Tosatti E Physics of Sliding Frictio n (NATO Science Series E) New York: Springer Publishing Company 1996. 73. Loezos P Costamagna P Sundaresan S The role of contact stresses and wall friction on fluidization. Chemical Engineering Science. 2002; 57: 5123 5141. 74. Janssen H Versuche uber getrei dedruck in silozellen. Z Ver Deutsch Ing 1895; 39: 1045 1049. 75. Chhabra R Agarwal L Sinha N Drag on nonspherical particles: an evaluation of available methods Powder Technology 1999; 101: 288 295. PAGE 254 254 BIOGRAPHICAL SKETCH Akhil Rao was born in 1986 and br ought up in Mumbai India. He first studied at Bombay Scottish School (Indian Certificate of Secondary Education, Delhi board) and graduated from there in 2001. The same year, he enrolled in Ramnivas Ruia Junior College and two years later obtained the Hig her Secondary Certificate (Maharashtra board). In 2003, he began his undergraduate studies at the Institute of Chemical Technology (Mumbai, India) and successfully graduated from there with a b achelor s in chemical e ngineering in 2007. In August 2007 he joined University of Florida (Gainesville, Florida) to pursue his PhD in chemical engineering. 