Kinetic theory for rough spheres : numerical and experimental
study of dense gas-solid fluidized beds
Citation for published version (APA):
Yang, L. (2017). Kinetic theory for rough spheres : numerical and experimental study of dense gas-solid fluidized beds. Technische Universiteit Eindhoven.
Document status and date: Published: 05/07/2017
Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)
Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:
www.tue.nl/taverne
Take down policy
If you believe that this document breaches copyright please contact us at:
openaccess@tue.nl
numerical and experimental study of
dense gas-solid fluidized beds
PROEFSCHRIFT
ter verkrijging van de graad van doctor aan de
Technische Universiteit Eindhoven, op gezag van de
rector magnificus, prof.dr.ir. F.P.T. Baaijens, voor een
commissie aangewezen door het College voor
Promoties, in het openbaar te verdedigen
op woensdag 5 juli 2017 om 14:00 uur
door
Lei Yang
voorzitter: prof.dr.ir. R. Tuinier promotor: prof.dr.ir. J.A.M. Kuipers co-promotor: dr.ir. J.T. Padding
leden: prof.dr.ir. M. Van Sint Annaland prof.dr. J.G.M. Kuerten
prof.dr. R.F. Mudde (Technische Universiteit Delft) prof.dr.ir. S. Pirker (Johannes Kepler University Linz) prof.dr.-ing. S. Heinrich (Hamburg University of Technology)
Het onderzoek of ontwerp dat in dit proefschrift wordt beschreven is uitgevoerd in overeenstemming met de TU/e Gedragscode Wetenschapsbeoefening.
and co-promotor: dr.ir. J.T. Padding
This research was fully funded by the Europe Research Council (ERC), under its Advanced Investigator Grant scheme, contract number 247298 (Multiscale Flow).
Copyright c 2017 by Lei Yang, Eindhoven, The Netherlands
All the copyright reserved. No part of this publication may be reproduced or utilized in any form or by any means, electronic or mechanical, photocopying, recording, without the prior permission of the author.
A catalogue record is available from the Eindhoven University of Technology Library ISBN:978-90-386-4301-4
Flows of dense assemblies of granular materials are encountered in many industrial de-vices, such as hoppers, fluidized beds, circulating fluidized beds (CFBs), spouted beds, etc. It is well-known that gas-fluidized beds have high efficiencies in particle mixing and particle-fluid heat and mass transfer. Understanding the dynamics of fluidized beds is a key issue in improving efficiency, reliability and scale-up from laboratory models.
The so-called Euler-Euler Two-Fluid Model (TFM), incorporating the kinetic the-ory of granular flow (KTGF), has been widely used. Usually the KTGF is based on collisions between smooth frictionless spheres. In reality, however, granular materi-als are often frictional. In dense gas-solid fluidization processes, the particles interact through frictional contacts between multiple neighbours. After a collision between rough particles, the particles can rotate due to surface friction. Thus, translational kinetic energy may exchange with rotational kinetic energy. Consequently, it is nec-essary to establish an accurate rheological description of the solid phase to consider particle surface friction. Besides, in rapid granular flows (which often occur in fluidized beds), the interactions between the solid particles and the wall lead to a rapid succes-sion of almost instantaneous collisucces-sions. These collisucces-sions cause random fluctuations of the particle velocities, which determine the momentum and fluctuation energy transfer rates to the walls. Thus, the physics of particle-wall collisions should be incorporated into granular flows. In this thesis, we developed a kinetic for rough spheres and corre-sponding boundary conditions at frictional walls to physically account for the effect of particle-wall collisions.
First of all, a kinetic theory of granular flow for frictional spheres in dense systems, including rotation, sliding and sticking collisions is derived. We use the Chapman-Enskog solution procedure of successive approximations, where the single-particle ve-locity distribution is assumed to be nearly Maxwellian both translationally and ro-tationally. An expression for the first-order particle velocity distribution function is derived, which includes the effects of particle rotation and friction. Using a simple mo-ment method, balance equations for mass, momo-mentum and energy are derived together with closure equations for viscosities, and thermal conductivities and collisional energy
of the friction coefficient. The model has been incorporated into our in-house two-fluid model code for the modeling of dense gas-solid fluidized beds. For verification, a com-parison of the present model in the limit of zero friction with the classical (frictionless) KTGF mode is carried out.
Then, the hydrodynamics of a dense solid-gas fluidized bed is studied, using a two fluid model based on our newly developed kinetic theory of granular flow for rotat-ing rough particles. The TFM simulations are validated by comparrotat-ing with PIV-DIA experimental data and results obtained from discrete particle model simulations of bubbling gas-fluidized beds. However, the predicted rotational granular temperature (rotational fluctuating energy) in our TFM simulations shows an almost uniform dis-tribution in the bed as a result of the assumptions that both the local mean rotational velocity and the gradient of the rotational granular temperature at the wall are zero, indicating directions for future improvement.
Also, the effect of collisional parameters, like the normal restitution coefficient and friction coefficient, on the hydrodynamics of a dense bubbling solid-gas fluidized bed is investigated using a two fluid model based on our KTGF for rotating frictional particles. We performed a comparison between TFM simulations using the present KTGF model, a simpler KTGF model for rapid flows of slightly frictional, nearly elas-tic spheres, and discrete parelas-ticle model simulation results. This comparison reveals that both the coefficient of normal restitution and the friction coefficient play an im-portant role with respect to the homogeneity of a bubbling bed. Besides, variation of the normal restitution coefficient slightly influences the time-averaged solids distribu-tion and the gas-solids flow field. However, the fricdistribu-tion coefficient has a stronger effect on the flow pattern and the distribution of the solid phase.
Next, in order to account for particle-wall collision, we propose a theory for the momentum and energy transfer between a flow of granular particles and a frictional wall. The theory describes the physics of collisions between frictional particles and flat walls, including both rotational and translational granular temperature. We compare our theory for steady state flows with the data from literature. Based on the theory, new BCs for particle slip velocity, and the flux of translational and rotational fluctu-ation energy in the framework of our kinetic theory are derived. Both simulfluctu-ations of a bubbling pseudo-2D fluidized bed using a Two Fluid Model with these new BCs, as well as simulations using a Discrete Particle Model have been performed to investigate the performance of the new BCs.
The present TFM simulations with our rough wall BCs are validated by compar-ison with experimental data obtained from Magnetic Particle Tracking and results obtained from Discrete Particle Simulations of bubbling fluidized beds. It has been reported that two-dimensional simulations may yield significantly different results com-pared to three-dimensional simulations. Thus, both pseudo-two dimensional and full
due to a decrease of the frictional effect from the front and back walls. Our current TFM model produces results which are in very good agreement with experimental data and results obtained from discrete particle simulations in terms of the time-averaged bed hydrodynamics.
Finally, it is known that the most widely used bed geometries in industry are three-dimensional cylindrical beds. So far, detailed three-three-dimensional studies are rarely found in the literature because of high computational cost. From the experimental point of view, the flow visualization and measurements are difficult to perform. Thus three-dimensional simulation is of significant scientific and technological interest. To this end, our KTGF was incorporated into our in-house three-dimensional cylindrical TFM code.
Dichte granulaire stromingen komen voor in veel industri¨ele apparaten, waaronder hop-pers, geflu¨ıdiseerde bedden, circulerende geflu¨ıdiseerde bedden, spuitbedden, etc. Het is bekend dat gas geflu¨ıdiseerde bedden hoge effici¨enties vertonen voor het mengen van deeltjes en voor warmte- en massaoverdracht tussen de vaste deeltjes en het gas. Het begrip van de dynamica van geflu¨ıdiseerde bedden is van essentieel belang voor het verbeteren van de effici¨entie, betrouwbaarheid en opschalen van laboratoriummodellen. Het zogenaamde Euler-Euler Twee-Fase Model (TFM), gebaseerd op de kinetische theorie voor granulaire stroming (KTGF), wordt veel toegepast. Doorgaans wordt de KTGF gebaseerd op botsingen tussen gladde wrijvingsloze bollen. In werkelijkheid zijn granulaire materialen echter vaak ruw. In dichte gas-vast flu¨ıdisatieprocessen interac-teren de deeltjes via wrijvingscontacten met meerdere buren. Na botsingen tussen ruwe deeltjes kunnen de deeltje roteren ten gevolge van de oppervlaktewrijving. Zodoende kan translationele kinetische energie worden uitgewisseld met rotationele kinetische en-ergie. Daarom is het nodig een nauwkeurige reologische beschrijving van de vaste fase te vinden die deze wrijvingseffecten tussen de deeltjesoppervlakken in ogenschouw neemt. Bovendien zullen in snelle granulaire stromingen (welke vaak optreden in geflu¨ıdiseerde bedden) de interacties tussen de vaste deeltjes en de wanden leiden tot een snelle opeenvolging van bijna instantane botsingen. Deze botsingen leiden tot willekeurige fluctuaties in de deeltjessnelheden, welke de overdrachtssnelheden bepalen van impuls en fluctuatie energie naar de wanden. De fysica van deeltjes-wand botsingen moet dus ook meegenomen worden. In dit proefschrift hebben we een kinetische theorie on-twikkeld voor ruwe bollen en bijbehorende randvoorwaarden voor ruwe wanden die de effecten van deeltjes-wand botsingen vangen.
Als eerste leiden we een kinetische theorie af voor ruwe bollen in dichte systemen, in-clusief rotatie en glijdende en plakkende botsingen. We gebruiken de Chapman-Enskog oplossingsprocedure van opeenvolgende benaderingen, waarbij wordt aangenomen dat de ´e´n-deeltjes snelheidsverdeling bijna Maxwelliaans is, zowel translationeel als rota-tioneel. Een uitdrukking voor de eerste-orde deeltjessnelheidsdistributiefunctie wordt afgeleid, inclusief de effecten van deeltjesrotatie en wrijving. Gebruik makend van een
en dissipatiesnelheden voor rotationele en translationele energie ten gevolge van de deeltjesbotsingen. In de resulterende sluitingsrelaties worden de reologische eigen-schappen van de deeltjes expliciet beschreven in termen van de wrijvingsco¨effici¨ent. Het model is ingebouwd in ons eigen twee-fase model code voor het modelleren van dichte gas-vast geflu¨ıdiseerde bedden. Ter verificatie wordt een vergelijking gemaakt tussen het nieuwe model in de limiet van geen wrijving en het klassieke (wrijvingsloze) KTGF model.
Daarna wordt de hydrodynamica van een dicht gas-vast geflu¨ıdiseerd bed bestudeerd, gebruik makend van een twee-fase model gebaseerd op onze nieuwe kinetische the-orie voor granulaire stroming van roterende ruwe deeltjes. De TFM simulaties zijn gevalideerd door ze te vergelijken met PIV-DIA experimentele data en resultaten van discrete deeltjesmodel simulaties van bubbelende gas-geflu¨ıdiseerde bedden. Echter, de voorspelde rotationele granulaire temperatuur (rotationele fluctuatie energie) in onze TFM simulaties laat een bijna uniforme verdeling zien in het bed als gevolg van de aannames dat zowel de locale gemiddelde rotatiesnelheid en de rotationele granulaire temperatuur aan de wand nul zijn, wat de richting voor verdere verbeteringen laat zien. Ook het effect van de botsingsparameters, zoals de normale restitutieco¨effici¨ent, op de hydrodynamica van een dicht vast-gas geflu¨ıdiseerd bed is onderzocht door middel van een twee-fase model model gebaseerd op onze KTGF voor roterende ruwe deeltjes. We hebben een vergelijking gemaakt tussen TFM simulaties op basis van het huidige KTGF model, een simpeler KTGF model voor snelle stromingen van licht ruwe, licht inelastische deeltjes, en simulaties van een discreet deeltjesmodel. Deze vergelijking laat zien dat zowel de normale restitutieco¨effici¨ent als de wrijvingsco¨effici¨ent een be-langrijke rol spelen wat betreft de homogeniteit van een bubbelend bed. Bovendien is er een lichte invloed van de normale restitutieco¨effici¨ent op de tijdsgemiddelde deelt-jesverdeling en gas-vast stromingsveld. De wrijvingsco¨effici¨ent heeft echter een sterker effect op het stromingspatroon en de deeltjesverdeling.
Vervolgens, om rekening te houden met de deeltjes-wand botsingen, stellen we een theorie voor die de impuls- en energieoverdracht tussen een granulaire stroming en een ruwe wand voorspelt. De theorie beschrijft de fysica van botsingen tussen ruwe deelt-jes en vlakke wanden, inclusief de rotationele en translationele granulaire temperatuur. We vergelijken onze theorie voor constante stromingen met data uit de literatuur. Gebaseerd op deze theorie, zijn nieuwe randvoorwaarden voor deeltjes slipsnelheid en de flux van translationele en rotationele fluctuatie theorie afgeleid. Zowel simulaties van een bubbelend pseudo-2D geflu¨ıdiseerd bed d.m.v. een twee-fase model met deze nieuwe randvoorwaarden als simulaties d.m.v. een discreet deeltjesmodel zijn uitgevo-erd om de prestaties van de nieuwe randvoorwaarden te onderzoeken.
De huidige TFM simulaties met onze ruwe wand randvoorwaarden zijn gevalideerd door een vergelijking te maken met experimentele data verkregen d.m.v. Magnetic
Par-sultaten kunnen laten zien dan drie-dimensionale simulaties. Daarom zijn pseudo-twee dimensionale en volledig drie-dimensionale bedden onderzocht in dit werk. We hebben gevonden dat de deeltjes meer homogeen verdeeld zijn in een volledig 3D bed, vergeleken met een pseudo-2D bed, ten gevolge van een reductie van wrijvingseffecten met de voor- en achterwanden. Ons huidige TFM model produceert resultaten welke in goede overeenstemming zijn met experimentele data en resultaten uit discrete deelt-jessimulaties betreffende tijdsgemiddelde bed hydrodynamica.
Tenslotte is het bekend dat de meest gebruikte bed geometrie¨en in de industrie drie-dimensionale cylindrische bedden zijn. Tot nu toe worden gedetailleerde drie-dimensionale studies nauwelijks gerapporteerd in de literatuur vanwege de hoge com-putationele kosten. Vanuit experimenteel oogpunt zijn de stromingsvisualisatie en metingen moeilijk uit te voeren. Vandaar dat er een significante wetenschappelijke en technologische interesse is in drie-dimensionale simulaties. Daarom is onze nieuwe KTGF ook ingebouwd in onze TFM code voor drie-dimensionale cylindrische bedden.
Roman symbols
c particle translational velocity, m/s C particle fluctuating velocity, m/s
v particle mean translational velocity, m/s G particle relative contact velocity, m/s J impulse force, kg m/s
k unit normal vector at contact point t unit tangential vector at contact point P pressure, Pa
F force, N T torque, N/m
I moment of inertia, kg m2
g gravitational acceleration, m/s2
r position of the particle, m
f(0) single particle distribution function, first approximation
f(1) single particle distribution function, second approximation
f(2) pair particle distribution function t time, s
V mean contact velocity, m/s H height from the distributor, m m particle mass, kg
n particle number density, m−3 e normal restitution coefficient Umf minimum fluidization velocity, m/s
ε volume fraction ρ density, kg/m3
τ stress tensor, Pa
βA interphase momentum exchange coefficient
Θ granular temperature, m2/s2
κ thermal conductivity, kg/(m·s) γ energy dissipation rate, J/(m3s)
µ friction coefficient
β tangential restitution coefficient λ granular temperature ratio θ contact angle
σ particle diameter, m
ω particle rotational velocity, rad/s
Ω particle fluctuating rotational velocity, rad/s δ overlap, m
η damping coefficient
χc collisional source of particle properties
θc collisional flux of particle properties
Abbreviations and subscripts
TFM Two Fluid Model DPM Discrete Particle Model KTGF kinetic theory of granular flow 2D two dimension
3D three dimension BC boundary condition MPT Magnetic Particle Tracking i in the i direction j in the j direction t translation r rotation p particle w wall s solid g gas
Summary i
Samenvatting iii
Nomenclature ix
1 General introduction 1
1.1 Background and motivation . . . 1
1.2 Research objective . . . 3
1.3 Multiscale modeling . . . 4
1.4 Magnetic particle tracking . . . 5
1.5 Thesis outline . . . 5
2 Derivation and implementation of kinetic theory for rough spheres 7 2.1 Introduction . . . 8
2.2 Kinetic theory of rough spheres . . . 10
2.2.1 Binary collisions . . . 10
2.2.2 Transport equations and collisional production . . . 12
2.2.3 Balance laws . . . 14
2.2.4 Second approximation and linear theory . . . 17
2.2.5 Constitutive equations . . . 21
2.3 Numerical solution method . . . 24
2.3.1 Discretization . . . 24
2.3.2 Boundary conditions . . . 28
2.4 Validation results and discussion . . . 29
2.4.1 Implementation validation . . . 30
2.4.2 Influence of particle friction . . . 31
3 Numerical study of a bubbling fluidized bed using kinetic theory for
rough spheres 35
3.1 Introduction . . . 36
3.2 Numerical simulations . . . 37
3.2.1 Initial conditions, boundary conditions, and parameter settings . 37 3.2.2 Momentum exchange coefficient . . . 38
3.2.3 Measurements of particle height, granular temperature, and en-ergy budget . . . 38
3.3 Results and discussion . . . 40
3.4 Conclusion . . . 48
4 Investigation of sphere collisional parameters in fluidized beds 49 4.1 Introduction . . . 49
4.2 Simulation settings . . . 50
4.3 Results and discussion . . . 51
4.3.1 Influence of normal restitution coefficient . . . 51
4.3.2 Influence of friction coefficient . . . 54
4.4 Conclusions . . . 57
5 Partial slip boundary conditions for collisional granular flows at flat frictional walls 59 5.1 Introduction . . . 60
5.2 Theory . . . 62
5.2.1 Particle-wall collisions . . . 62
5.2.2 Probability distribution functions . . . 64
5.2.3 Balances at the wall . . . 65
5.2.4 Collisional rate of change . . . 66
5.2.5 Fluctuation heat flux . . . 72
5.2.6 Numerical implementation of wall boundary conditions . . . 75
5.3 Results and discussion . . . 77
5.4 Conclusions . . . 85
6 Investigation of bubbling dense solid-gas fluidized beds using kinetic theory for rough spheres 87 6.1 Introduction . . . 88
6.2 Numerical Methodology . . . 89
6.2.1 Two Fluid Model . . . 89
6.2.2 Discrete Particle Model . . . 89
6.3 Experimental Setup . . . 91
6.3.1 Single particle impact . . . 91
6.3.2 Magnetic particle tracking setup . . . 92
6.4 Results and discussion . . . 94
6.4.1 Pseudo-2D system . . . 94
6.4.2 3D system . . . 102
7.1 Introduction . . . 110
7.2 Numerical method . . . 111
7.3 Results and discussion . . . 112
7.3.1 Perfectly elastic particles . . . 112
7.3.2 Alumina beads . . . 112
7.4 Conclusions . . . 117
8 Conclusions and Outlook 119 Appendices 122 A Derivation of Maxwell transport equation for rough spheres 123 A.1 Boltzmann’s equation expressed in terms of the peculiar velocity . . . . 123
A.2 Maxwell’s transport equation . . . 124
B Table of Integrals 127 C Calculation of Collisional Changes 131 C.1 Calculation of shear stress tensor . . . 131
C.2 Collisional rate of translation . . . 133
C.3 Collisional rate of rotation . . . 135
Bibliography 145
Publications 147
Acknowledgments 149 Curriculum Vitae 151
General introduction
1.1
Background and motivation
Gas-solid fluidized beds find a widespread application in processes involving combus-tion, separacombus-tion, classificacombus-tion, and catalytic cracking (Davidson et al. (1985); Kunii and Levenspiel (1991)). It is well known that fluidized beds have high efficiencies in par-ticle mixing and parpar-ticle-fluid heat and mass transfer (Gilliland and Mason (1949)). Understanding the dynamics of fluidized beds is a key issue in improving efficiency, reliability and scale-up. Owning to the enormous increase in computer power and al-gorithm development, fundamental modelling of multiphase reactors has become an effective tool.
Flows of dense assemblies of granular materials are encountered in these fluidization processes, which defines a complex two-phase flow problem. Thus, understanding and modeling gas-particle flows involving dense assemblies of particles is of significant scien-tific and technological interest. As reported by Forterre and Pouliquen (2008), granular flows are often categorized into three different regimes: a dense quasi-static regime in which the deformations are very slow and the particles interact by frictional contacts (Roux and Combe (2002)); a gaseous regime in which the flow is very rapid and dilute, where the particles interact by collision (Campbell (1990)); and an intermediate liquid regime in which the flow is dense but still behaves like a liquid, the particles interacting both by collision and friction (Pouliquen and Chevoir (2002)). Among these regimes, the last liquid regime is most often encountered in applications.
Gi-daspow (2003)) is considered to be the most effective for large scale simulations. The two-fluid model (TFM) incorporating the kinetic theory of granular flow (KTGF) has been widely used. In dense gas-solid fluidization processes, the particles interact with each other through enduring frictional contact between multiple neighbours and, to a lesser extent, through collisions. In many cases, as pointed by Srivastava and Sundare-san (2003), gravity is sufficient to ensure that frictional interactions have a significant contribution to the particulate stress. As the particle volume fraction decreases, the collisional stress becomes more dominant. In these systems, the contributions from the frictional stresses and kinetic stresses are comparable. Thus, it is of significant interest to develop novel advanced KTGF models that combine the frictional and kinetic con-tributions.
Until now, most KTGF models focus on frictionless particles. In reality, however, granular materials are frictional. After a binary collision between rough particles, the particles can rotate due to surface friction. Thus, translational kinetic energy may exchange with rotational kinetic energy. In addition to the mass and momentum equa-tions, there is the requirement of an additional equation describing the particle spin. The roughness of the granular materials has been shown to have a significant effect on stresses at least in the quasi-static regime (Sun and Sundaresan (2011)). Figure 1.1 shows the Discrete Particle Simulation results with and without particle surface friction. It is evident that particle surface friction leads to a more pronounced hetero-geneous particle flow structure. Besides, the particle clustering and the downward solid flow in the main gas stream may also cause some gas back mixing which may impact the process efficiency. Consequently, it is of significant importance to characterize this friction effect.
In the literature, many researchers (Lun and Savage (1987); Lun (1991); Goldshtein and Shapiro (1995); Chou and Richman (1998); Kumaran (2006)) have tried to quan-tify the friction effect. However, these attempts have been somewhat limited. Thus, an accurate rheological description of the solid phase in kinetic theory of granular flow including particle friction and rotation is still lacking.
In addition, in rapid granular flows (which often occur in fluidized beds), the in-teractions between the solid particles and confining or internal walls lead to a rapid succession of almost instantaneous collisions. These collisions cause random fluctua-tions of the particle velocities, which determine the momentum and fluctuation energy transfer rates to the walls. Therefore, the physics of particle-wall collisions should also be incorporated into this type of granular flow simulations.
Many experimental and computational studies have been reported in the literature using pseudo two-dimensional (2D) beds (Taghipour et al. (2005); Yang et al. (2016d)). However, on the one hand, restricting the computational domain to two dimensions assumes no spatial dependence in the depth direction. An assumption which is ques-tionable because most experimental data has been obtained from thin fluidized beds to enable visual access. On the other hand, wall effects (particularly the particle-wall
friction) play a critical role in the bubble dynamics (Li et al. (2010)) in pseudo-2D labo-ratory and small-scale experiments. A comparison by Li et al. (2010) clearly shows that there are marked differences between two-dimensional and three-dimensional simula-tion predicsimula-tions. Those from the three-dimensional simulasimula-tions match the experimental results significantly better than those from the two-dimensional simulations. So far, the flow characteristics obtained for three-dimensional (3D) fluidized beds were mostly compared with 2D or pseudo-2D systems. Cranfield and Geldart (1974) also reported that 2D simulations may yield significantly different results compared with 3D simula-tions. Clearly, the most widely used bed geometries in industry are 3D. So far, detailed 3D studies are rarely found in the literature because of the high computational cost (Bakshi et al. (2014)). From an experimental point of view, flow visualization and measurements are difficult to perform in 3D and non-optical techniques must be used to study the bed hydrodynamics.
Figure 1.1: Influence of particle surface friction from DPM simulations, e = 0.97, β = 0.33, t = 1.8s
1.2
Research objective
The objective of this study is 1) to establish an accurate rheological description of the solid phase (the so-called kinetic theory for rough spheres) and corresponding boundary conditions, to be used in the frame work of the Two Fluid Model, and 2) apply and vali-date this theory by comparing the predicted hydrodynamics in dense gas-solid fluidized beds with experimental results and state-of-the-art Discrete Particle Simulations.
1.3
Multiscale modeling
Due to the large range of relevant scales in industrial size fluidization equipment, a multiscale modeling approach is necessary, as described in detail in Van der Hoef et al. (2004). The usual approach involves three levels of modeling (see Figure 1.2): the Direct Numerical Simulation (DNS), the Discrete Particle Model and the continuum Two Fluid Model. At the most detailed level, DNS (the fluid and particulate phases are treated by considering the Navier-Stokes equation and the Lagrangian equations of motion respectively) resolves the flow around the particles and can be used to establish drag closures. Those drag closures are then employed in DPM and TFM simulations.
Figure 1.2: Multiscale modeling scheme for dense gas-solid fluidized beds The Discrete Particle Model (DPM) is based on a combined Eulerian-Lagrangian approach. The solid particles are treated explicitly and their motion is described by Newtons second law. The forces on the particles arise from gravity, the interaction with the fluid phase, and collisions with other particles. For efficiency the averaged continuous fluid phase equations are solved on a grid which is much larger than the particle size. As a result, both a drag-force closure and a collision model have to be specified for the DPM. Details of the DPM model and its application are given by Ristow (2000), Deen et al. (2007), and Zhu et al. (2007).
The Two-Fluid Model is based on an Eulerian-Eulerian approach and can in prin-ciple be used for large scale simulations. In TFM the solid phase is treated as a second continuum, interpenetrating with the continuous gas phase. Balance equations
for mass, momentum and energy are solved by using additional closure equations for stress, viscosities, thermal conductivities, and energy dissipation rates (Kuipers et al. (1992); Wang et al. (2012)). The Euler-Euler approach has emerged as a very promis-ing tool as a result of its compromise between computational cost and amount of detail provided.
1.4
Magnetic particle tracking
Even though numerical simulations can be used to generate a detailed understanding of flow structures in fluidized beds, validation of the models using advanced and de-tailed experiments are still crucial. Due to the opaqueness of dense 3D fluidized beds, non-optical techniques must be used such as electrical resistance tomography, electri-cal capacitance tomography, positron emission particle tracking, magnetic resonance imaging, and Magnetic Particle Tracking (MPT). In particular, MPT has emerged as a promising technique to investigate the hydrodynamics of fluidized beds due to its long-term stability, safety, and relatively low costs (Mohs et al. (2009)). The method uses a magnetic tracer particle, which follows the bulk particle flow and whose mag-netic field is continuously detected by multiple magmag-netic field sensors located outside the bed. This allows the reconstruction of the instantaneous position and orientation of the magnetic tracer. Based on statistical analysis of the tracer trajectory, charac-teristic measures of the bulk particle flow, such as the average particle velocity and particle circulation pattern, can be determined as a function of operating conditions. The application of MPT in fluidization was initiated by Mohs et al. (2009) for the study of a spouted bed. Recently, MPT was improved by Buist et al. (2015) for the use in dense granular flows occurring in bubbling fluidized beds.
1.5
Thesis outline
The research reported in this thesis consists of several elements. In Chapter 2, a ki-netic theory of granular flow for frictional spheres in dense systems, including rotation, sliding and sticking collisions is derived. Using a simple moment method, balance equations for mass, momentum and energy are derived with closure equations for vis-cosities, and thermal conductivities and collisional energy dissipation rates of angular and translational kinetic energy. The model has been incorporated into our in-house two-fluid model code for the modeling of dense gas-solid fluidized beds. For verifica-tion, a comparison of the present model in the limit of zero friction with the classical (frictionless) KTGF model is carried out.
In Chapter 3, the hydrodynamics of a dense gas-solid fluidized bed is studied, using a Two Fluid Model based on our newly developed kinetic theory of granular flow for rotating rough particles. The TFM simulations are validated by comparing with PIV-DIA experimental data and results obtained from discrete particle model simulations of a bubbling fluidized bed.
The effect of collisional parameters on the hydrodynamics of a dense bubbling gas-solid fluidized bed using a Two Fluid Model based on our kinetic theory for rough spheres has been studied as well and is reported in Chapter 4.
In Chapter 5, a theory is proposed for the momentum and energy transfer between a flow of granular particles and a frictional wall in order to better take into account particle-wall collisions. Based on this theory, new BCs for particle slip velocity, and the flux of translational and rotational fluctuation energy in the framework of our kinetic theory are derived. Besides, simulations of a bubbling pseudo-2D fluidized bed using a two fluid model with the new BCs are compared with a Discrete Particle Model to investigate the performance of the new BCs.
In Chapter 6, the TFM simulations embedding our rough wall BCs are compared with experimental data obtained from Magnetic Particle Tracking and results obtained from discrete particle simulations of bubbling fluidized beds.
Finally, in Chapter 7, the incorporation of the new kinetic theory for rough spheres into our in-house three-dimensional cylindrical TFM code is reported. Finally, the model is verified by detailed comparison with the model from Jenkins and Zhang for the simulation of rough particles in cylindrical beds.
Derivation and
implementation of kinetic
theory for rough spheres
We derive a kinetic theory of granular flow (KTGF) for frictional spheres in dense sys-tems, including rotation, sliding and sticking collisions. We use the Chapman-Enskog solution procedure of successive approximations, where the single-particle velocity dis-tribution is assumed to be nearly Maxwellian both translationally and rotationally, as assumed by McCoy et al. (1966). An expression for the first-order particle velocity dis-tribution function is derived, which includes the effects of particle rotation and friction. Using a simple moment method, balance equations for mass, momentum and energy are derived with closure equations for viscosities, and thermal conductivities and col-lisional energy dissipation rates of angular and translational kinetic energy. Because the internal angular momentum changes are coupled to the flow field, the stress tensor contains anti-symmetric components which are associated with a rotational viscosity. In the resulting closure equations, the rheological properties of the particles are explic-itly described in terms of the friction coefficient. The model has been incorporated into our in-house two-fluid model (TFM) code for the modeling of dense gas-solid fluidized beds. For verification, a comparison of the present model in the limit of zero friction with the original (frictionless) KTGF model is carried out. Simulation results of both models agree well.
2.1
Introduction
Gas-solid fluidized beds find a widespread application in processes involving combus-tion, separacombus-tion, classificacombus-tion, and catalytic cracking (Davidson et al. (1985); Kunii and Levenspiel (1991)). Understanding the dynamics of fluidized beds is a key issue in the design and scale-up. With the improvements in computational power, Compu-tational Fluid Dynamics (CFD) simulation has become an effective tool for obtaining depth knowledge on fundamentals in multiphase flows.
For larger scale applications, the Euler-Euler approach (Kuipers et al. (1992); Gi-daspow (1994); Lu and GiGi-daspow (2003)) is considered to be the most effective. The two-fluid model (TFM) has emerged as a very promising tool as a result of its com-promise between computational cost and amount of detail provided. The challenge of the model is to establish an accurate hydrodynamic and rheological description of the solid phase. State-of-the-art closures have been obtained from the kinetic theory of granular flow, initiated by Jenkins and Savage (1983); Lun et al. (1984); Jenkins and Richman (1986a); Lun and Savage (1987); Ding and Gidaspow (1990); Lun (1991) and Nieuwland et al. (1996). This theory is based on the classical kinetic theory of dense gases (Chapman and Cowling, 1970).
The original KTGF models of Jenkins and Savage (1983); Lun et al. (1984); Jenkins and Richman (1986a) and Gidaspow (1994) are derived for frictionless and nearly elas-tic parelas-ticles with translational motion only. In reality, however, granular materials are frictional. After a binary collision between rough particles, the particles can rotate due to surface friction. Thus, translational kinetic energy may exchange with rotational kinetic energy. Attempts to quantify the friction effect have been somewhat limited.
Lun and Savage (1987) developed a kinetic theory for a system of inelastic, rough spherical particles to study the effects of particle surface friction and rotational iner-tia. Two collisional parameters were used. One was the normal restitution coefficient and the other was the roughness coefficient to characterize the effects of surface fric-tion. In their theory only dense systems, where the major stress contribution arises from particle collisions, was considered. However, at low solids concentrations the ki-netic stress contribution becomes dominant. Lun (1991) extended the kiki-netic theory to be appropriate for both dilute and dense granular flows including kinetic as well as collisional contributions for stresses and energy fluxes. Translational and rotational granular temperatures were involved to characterize the kinetic fluctuation energies. An extra energy balance equation for the rotational granular temperature was derived. Goldshtein and Shapiro (1995) developed a kinetic theory for rough inelastic spheres combining the translational and rotational granular temperatures into a total granular temperature. Their theory has been improved by Wang et al. (2012) to develop kinetic theory for granular flow of rough spheres. The constitutive relation for the granular flow of smooth, inelastic, spherical particles was derived by Kumaran (2004) in the adi-abatic limit where the length scale for the conduction of energy is small compared to the macroscopic scale. Later, Kumaran (2006) developed constitutive relations for the
granular flow of rough spheres in the limit where the energy dissipation in a collision is small compared to the energy of a particle. He also pointed out the antisymmetric component in the stress tensor for rough and partially rough particles. However, these models did not account for the dependence on the friction coefficient explicitly.
Abu-Zaid and Ahmadi (1990) introduced the friction coefficient into the kinetic theory model to include the effects of frictional energy dissipation during collision ex-plicitly. But they made the assumption that the particles were sufficiently small such that the effects of particle-spin could be neglected. Jenkins and Zhang (2002) derived a simple kinetic theory for collisional flows of identical, slightly frictional, nearly elastic spheres. An effective restitution coefficient which was given in terms of the collision parameters, namely normal coefficient of restitution e, friction coefficient µ, and tan-gential coefficient of restitution β. However, a comparison between experiment and simulation by Goldschmidt et al. (2001) for an effective restitution coefficient of 0.73 showed that the simulated bubble dynamics was too vigorous, which indicated that the effects of particle friction could not be replaced by using an effective (smaller) restitu-tion coefficient.
Rao et al. (2008) derived a kinetic theory for rough inelastic particles. They used a simple hard sphere collision model with two restitution coefficients, namely the nor-mal restitution coefficient e (0 < e < 1) and the tangential restitution coefficient β (−1 < β < 1). In effect, therefore, their theory applies to infinitely high friction coeffi-cient. However, in fluidized beds the particles mostly undergo sliding contacts. Thus, for the study of hydrodynamics in fluidized beds, inclusion of the friction coefficient is important because distinguishing between sliding and sticking collisions is necessary.
Recently, Zhao et al. (2013) developed a kinetic theory model for rough spheres. Both particle friction and rotation were considered for energy fluxes. However, they did not derive the expression for the first-order particle velocity distribution function and the effect of rotation on the particle velocity distribution function was not included. In the Hermite polynomials proposed by Grad (1949) to write the velocity distribution function, the coefficients given by Enwald et al. (1996) were used without considering the particle friction. Sun and Sundaresan (2011) proposed a new expression for the ra-dial distribution function at contact as well as modifications to the Garz´o-Dufty (GD) KT model (Garz´o and Dufty (1999)) expressions for shear stress and energy dissipation rate through data obtained from discrete element method (DEM) simulations. Their model relied on empiricism, particularly in the corrections for the shear stress. Accord-ing to their simulation, the granular temperature in the dilute region (solids volume fraction < 0.2) was overestimated due to neglect of the rotational degrees of freedom. The present work is concerned with further development of a kinetic theory for rough spheres including particle friction and rotation. We focus on the inertial regime (Chialvo et al. (2012)) so that kinetic theory can be applied without corrections for the drag caused by the interstitial fluid. Following Kumaran (2006), an antisymmetric contribution due to the convective transport of rotational momentum is incorporated
in the solid stress tensor. The present model is first incorporated into our in-house Euler-Euler code. Simulations of gas-solid fluidized bed are carried out. Comparison between the present model and the original model is then made.
2.2
Kinetic theory of rough spheres
2.2.1
Binary collisions
Figure 2.1: Sketch of a binary collision between frictional spheres.
We start the analysis by considering an inelastic collision between two frictional spheres (Figure 2.1). It is assumed that the interaction forces are impulsive and there-fore all other finite forces are negligible during a collision. In this work, we will mainly adopt the notation used by Walton (1993) and Foerster et al. (1994). The spheres have mass m, diameter σ, and moment of inertia I around their respective centers. The balance of linear momentum requires that the velocities c1 and c2 of the two spheres
before collision and the velocities c01, c02 after collision have the relation c01= c1+ J, c
0
2= c2− J (2.1)
where mJ is the collisional impulsive change of particle 1 induced by its collision with particle 2. The relations between angular velocities can be expressed by
ω01= ω1− mσ 2I k × J, ω 0 2= ω2− mσ 2I k × J (2.2) k is the unit vector directed from the center of particle 2 to the center of particle 1. The relative velocities of the centers of the spheres before and after a collision are
g = c2− c1 (2.3)
The relative velocity G at the contact points is given by G = g +σ
2 (ω1+ ω2) × k (2.5) J can be decomposed into its normal and tangential components,
J = Jnk + Jtt (2.6)
where the unit tangential vector t is defined as, t = (G × k) × k
|G × k| (2.7) Note that t is perpendicular to k, lies in the plane spanned by G and k, and that G · t ≤ 0 (Figure 2.1). At this point, three parameters are used to close the above constitutive relations. The first parameter is the normal restitution coefficient (0 ≤ e ≤ 1), i.e. the change in the normal component of the relative surface velocity after and before the collision:
G0· k = −e (G · k) (2.8) Furthermore, the coefficient of tangential restitution (−1 ≤ β0 ≤ 1) and friction
coefficient (µ ≥ 0) are defined as,
G0· t = −β0(G · t) , sticking (2.9)
|k × J| = µk · J, sliding (2.10) The above two equations should be interpreted as follows. The tangential com-ponent of the relative surface velocity is changed by a factor −β0. If β0 = −1 no
change in tangential surface velocity is taking place. In the other extreme, for β0= 1
the tangential surface velocity is fully inverted (like two cogwheels hitting each other tangentially). This corresponds to a sticking type of collision. However, the tangential impulsive change is limited to µ times the normal impulsive change, corresponding to a sliding (Coulomb-type) friction.
The normal component of the impulse vector can be determined by combining equa-tions 2.3, 2.4 and 2.8:
Jn =
1
2(1 + e) (G · k) (2.11) Whether the collision is of the sliding or sticking type depends on the angle θ between G and k (see Figure 2.1). For collisions of the sliding type (θ ≥ θ∗), the tangential impulsive change is given by
Jt= −
1
2µ (1 + e) (G · k) (2.12) For sticking collisions (θ ≤ θ∗), the tangential impulsive change is given by
Jt=
1 + β0
When θ = θ∗, the above two expressions of tangential impulse must agree, so using G · k = G cos θ and G · t = −G sin θ, we find
β0= −1 + µ(1 + e) 1 +mσ 2 4I cot θ∗ (2.14) Thus the critical angle is determined by
tan θ∗= µ1 + e 1 + β0 1 + mσ 2 4I (2.15) For homogenous spheres I =mσ2
10 , and thus tan θ∗= µ0= 7 2µ 1 + e 1 + β0 (2.16)
2.2.2
Transport equations and collisional production
The goal of any kinetic theory of granular flow is to obtain the transport balance equations for relevant hydrodynamic moments (e.g. mass-, momentum- and energy-density), along with closures for the transport properties such as shear viscosity and (pseudo-)thermal conductivity. In the following we will derive the Boltzmann equation for the probability distribution function of finding a particle at a specific location with a specific translational and rotational velocity. By taking ensemble averages we can then derive the corresponding Maxwell transport equations. Following Grad (1949), we assume that the state of the system can be described accurately by an extended, but limited, set of moments of the velocities and angular velocities. The balance equations for N th order moments do not form a closed system, because they generally contain fluxes and production terms which may involve moments of order (N + 1). Therefore, also following Grad (1949), we will construct a first approximation to the probability distribution function by expanding the unperturbed probability distribution function (which is a simple Maxwellian) into Hermite polynomials. The Hermite expansion co-efficients are then determined self-consistently, such that the first approximation of the probability distribution functions reproduces the expected moments. The system is closed because we consider as many Hermite coefficients are there are moments. Once the first approximation to the probability distribution function is known, we can com-pute fluxes and production terms for each of the hydrodynamic moments.
We denote the single-particle probability distribution function to encounter a par-ticle at position r with instantaneous velocity c and angular velocity ω at time t by f = f (r, c, ω, t). With this definition, the number of particles per unit volume can be obtained from n =R R f dcdω, and the mean (ensemble average) value of a quantity ψ is defined as:
n < ψ >= Z Z
ψf dcdω (2.17) Here, as in McCoy et al. (1966), we assume the unperturbed (zeroth-order) probability distribution function is of Maxwellian form for both the translational and rotational
velocity fluctuations: f(0)(r, c, ω, t) = n 1 2πΘt 3/2 I 2πmΘr 3/2 exp −1 2 C2 Θt + IΩ 2 mΘr (2.18) Here C (= c − v) and Ω (= ω − ω) are the peculiar translational and rotational ve-locities, defined relative to the local average particle velocity v and average angular velocity ω. The translational and rotational granular temperatures are a measure for the mean square peculiar velocities. More precisely they are given by Θt≡< C2> /3
and Θr≡ I < Ω2> /3m.
In deriving the Boltzmann equation it is assumed that only binary encounters are important. Based on Chapman and Cowling (1970), we consider a particle with ex-ternal force mF and exex-ternal torque I T, which may be a function of r and t (note that F and T are therefore the acceleration and angular acceleration, respectively, due to external forces and torques). For particles in fluidized beds, the dominant external forces and torques include gravity forces and fluid-particle interactions such as buoy-ancy, drag and lift. Between times t and t +dt, the velocities c, ω of any particle that does not collide with another will change to c+Fdt, ω+Tdt and its position vector will change to r+cdt. There are f (c, ω, r, t)dcdrdω particles at time t. After time interval dt, there are f (c + Fdt, ω + Tdt, r + cdt, t + dt)dcdrdω particles. The number of parti-cles in the second set differs from the first one. The net gain of partiparti-cles to the second set must be proportional to dcdωdrdt , and will be denoted by ∂ef
∂tdcdωdrdt. On
di-viding by dcdωdrdt and making dt tend to zero, Boltzmanns equation for f is obtained
F ·∂f ∂c + T · ∂f ∂ω + c · ∂f ∂r + ∂f ∂t = ∂ef ∂t (2.19) The quantity ∂ef
∂t is equal to the rate of change, owing to collisional encounters, in the
velocity distribution function f at a fixed point.
Multiplying Boltzmanns equation with a local property ψ(c, ω), which may depend on the particle velocity, angular velocity and/or position, and changing independent variables to the peculiar velocities C, Ω, instead of absolute velocities c, ω, we find after several transformations (similar to Chapman and Cowling (1970), p.48) the Maxwell transport equations: nD < ψ > Dt − n Fi− Dvi Dt ∂ψ ∂Ci + ∂ ∂ri (θc,i(ψ)+ < nCiψ >) +∂vi ∂rj n∂ψ ∂Ci Cj + θc,j ∂ψ ∂Ci + nDωi Dt ∂ψ ∂Ωi + n∂ψ ∂Ωi Cj + θc,j ∂ψ ∂Ωi ∂ωi ∂rj = χc(ψ) (2.20)
where for clarity the Einstein summation convention is used, i.e. a sum over all com-ponents (x, y, and z) is implied whenever the same component index (i or j) appears
twice in a product. Detailed derivation of Maxwell transport equations can be found in Appendix A. The rate of change ∂ef
∂t due to collisions can be written as follows
according to Jenkins and Richman (1986a), Z Z ψ∂ef ∂t dCdΩ = χc(ψ) − ∇ · θc(ψ) − ∇v · θc ∂ψ ∂C − ∇ω · θc ∂ψ ∂Ω (2.21) where the collisional flux and collisional source terms are
θc(ψ) = − σ3 2 Z G·k>0 (ψ0− ψ)G · kkf(2)(r + σk, c, ω, r, c 1, ω1, t)dkdc1dω1dcdω (2.22) χc(ψ) = σ2 2 Z G·k>0 (ψ0+ ψ01− ψ − ψ1)G · kf(2)(r + σk, c, ω, r, c1, ω1, t)dkdc1dω1dcdω (2.23) Note that the tensorial rank of the collisional flux term is one higher than that of ψ and the tensorial rank of the collisional source term is equal to that of ψ. As usual, we will assume that the 2-particle distribution function f(2) factorizes as a product of two
single-particle distribution functions, with position correlations included through the radial distribution function. Following a Taylor series gives the following expression for the 2-particle distribution function f(2)(r + σk, c, ω, r, c1, ω1, t)
f(2)(r + σk, c, ω, r, c1, ω1, t) = g0(s)f (r + σk, c, ω, t)f1(r, c1, ω1, t) = g0(s) f × f1+ 1 2σf f1k · ∇ ln f f1 (2.24) where g0is the radial distribution function at contact, which depends on the local solid
volume fraction s.
2.2.3
Balance laws
For convenience, in writing the balance laws corresponding to the higher moments of velocity fluctuation we adopt the following definitions
Mi1,i2...in≡ hCi1Ci2...Cini (2.25a)
Ni1,i2...in≡ hΩi1Ωi2...Ωini (2.25b)
χi1,i2...in≡ χ (mCi1Ci2...Cin) (2.25c)
θji1,i2...in≡ θj(mCi1Ci2...Cin) (2.25d)
Different balance equations can be obtained by inserting different properties ψ in equa-tion (Eq. 2.20).
(1) Let ψ be m D(mn) Dt + mn ∂vi ∂ri = 0 (2.26) which is the familiar continuity equation (mass conservation).
(2) Let ψ be mCi ρDvi Dt + ∂ ∂rj (θji+ ρMij) = ρFi (2.27)
Here, Pij= θji+ ρMij is the total pressure tensor (negative of the stress tensor). This
is the Navier-Stokes equation for the average particle momentum. We note that for smooth spheres θji = θij. However, for rough spheres this is not the case, i.e. the
pressure tensor contains asymmetric components. We will return to this point later. (3) Let ψ be IΩi nIDωi Dt + ∂ ∂rj θj(IΩi) = χ(IΩi) (2.28)
A collisional source term appears because the sum of angular momenta (around the respective spheres center of mass) is not conserved in binary collisions. Similar to Lun (1991), the collisional angular momentum flux θj(IΩi) is found to be zero.
(4) Let ψ be mCiCi, using the definition of granular temperature, the energy
bal-ance equation can be written as 3 2ρ DΘt Dt + 1 2 ∂ ∂rj (θjii+ ρMjii) + ∂vj ∂rk Pkj= 1 2χ(mCiCi) (2.29) Here ρMjii is the transported (kinetic) part of the energy flux vector, and θjii is the
collisional part.
(5) Let ψ be IΩiΩi, the energy balance equation for rotational granular temperature
can be written as 3 2ρ DΘr Dt + 1 2 ∂ ∂rj (θj(IΩiΩi)+ < nICjΩiΩi>) +∂ωj ∂rk (nI < ΩjCk> +θk(IΩj)) = 1 2χ(IΩiΩi) (2.30) We note that according to our calculations terms like < nICjΩiΩi >, nI < ΩjCk >
and θk(IΩj) vanish. Besides, in this work we will assume that in the bulk of the
system the mean rotational velocity is so small that we can set it to zero, in which case the production term due to gradients in the mean rotational velocity (last term on the left hand) side vanishes. So we do not have a production term for rotational
granular temperature in the present theory. This is in contrast with the work of Rao et al. (2008), who kept the production term due to the existence of a mean rotational velocity. Note that we do still have a collisional source/sink term on the right hand side. (6) Let ψ be mCiCj, the balance law for the mean of the second moment of velocity
fluctuation is derived by deducting the energy equation 2.29 1 2ρ D ˆMij Dt + ∂ ∂rk Qkij− 1 3δijQkpp + Pn(jvi),n− 1 3Pkmvm,kδij = 1 2χ(mCˆ iCj) (2.31) where we defined Qkij= 12(θkij+ ρMkij). Here and in the following a comma denotes
differentiation with respect to a spatial coordinate, e.g. ai,j = ∂r∂
jai, round brackets
around paired indices denotes symmetric averaging over all permutations of the indices, e.g. a(ibj)= 12(aibj+ ajbi) and a(ibjk)= 16(aibjk+ aibkj+ ajbik+ ajbki+ akbij+ akbji),
and the traceless part of a second-rank tensor is denoted by a hat above its symbol, e.g. ˆ
χ(mCiCj) = χ(mCiCj) −13χ(mCkCk)δij. Note that Eq. 2.31 describes the evolution
of the anisotropic part of the second moment of velocity fluctuations; the evolution of the isotropic part, which is the granular temperature, is given by Eq. 2.29.
(7) Let ψ be IΩiΩj, similarly, after deducing the angular energy equation (Eq.
2.30), the second moment of angular velocity fluctuation is 1 2ρ D ˆNij Dt + 1 2 ∂ ∂rk θk(IΩiΩj)+ < nICkΩiΩj> − 1 3δij(θk(IΩpΩp)+ < nICkΩpΩp>) +Pk(jωi),k− 1 3δij(θk(IΩp)+ < nIΩpCk>) ωp,k= 1 2χ(IΩˆ iΩj) (2.32) (8) Let ψ be mCiCjCk, the expression for the third moment of velocity fluctuation
can be written as ρDMijk Dt + 2Qlijk− 3 ∂ ∂rn Pn(iMjk)+ 6Qm(ijvk),m= χ(mCiCjCk) (2.33)
where we defined Qlijk= 12(θlijk+ ρMlijk).
(9) Let ψ be IΩiΩjCk, by using the linear moment balance law, the mixed third
moment of angular and translational velocity fluctuations can be expressed as nID < ΩiΩjCk > Dt + ∂ ∂rn (θn(IΩiΩjCk) + nIδijδkn< ΩiΩjCkCn>) + +∂vn ∂rm
(< δknnIΩiΩjCm> +δknθm< IΩiΩj >) + 2δinθm(IΩjCk)
∂ωn ∂rm − < δknIΩiΩj> ∂ m∂rm (θmn+ ρMnm) = χ(IΩiΩjCk) (2.34)
2.2.4
Second approximation and linear theory
The balance laws apply to averaged quantities < ψ >, and therefore depend on our as-sumption about the single particle distribution function f . Although in homogeneous systems it is reasonable to assume, to first order, a Maxwellian form, the distribution function is expected to alter under the influence of (relatively small) gradients in the average velocity, average angular velocity, translational granular temperature and ro-tational granular temperature fields. We therefore need a method to estimate these deviations. Following Grad (1949), we write the second approximation to the single particle distribution function with rotation as a Hermite polynomial series
f(1)(r, c, ω, t) = 1 − ai∂c∂i− bi∂Ω∂i + aij 2! ∂2 ∂ci∂cj + bij 2! ∂2 ∂Ωi∂Ωj +b 0 ij 2! ∂2 ∂ci∂Ωj − aijk 3! ∂3 ∂ci∂cj∂ck − cijk 2! ∂3 ∂ci∂cj∂Ωk −dijk 3! ∂3 ∂Ωi∂Ωj∂Ωk − bijk 2! ∂3 ∂ci∂Ωj∂Ωk+ ... f(0)(r, c, ω, t) (2.35)
Where f(0)(r, c, ω, t) is the (zeroth order) Maxwellian distribution function, Eq. 2.18. The coefficients ai, aij, aijk, bij, b
0
ij, bijk..., to be determined next, may depend on r and
t but crucially do not depend on c and ω. If the mean values of velocity fluctuations < Ci >, < Ωi >, < CiCj >, < ΩiΩj >, < CiΩj >,< CiCjΩk >, < ΩiΩjΩk > based
on the second approximation f(1)are carried out, it can be found that a
i= bi= bij =
b0ij = cijk = dijk= 0, aii= bii = 0. Also
Mij= Θtδij+ aij (2.36a)
Mijk= aijk (2.36b)
Mijkp= 3Θ2tδ(ijδkp)+ 6Θta(ijδkp)+ aijkp (2.36c)
bijk=< CiΩjΩk > (2.36d) < ΩiΩjCkCn> = δijδkn mΘtΘr I + aknδij mΘr I + bijkn (2.36e) Based on Jenkins and Richman (1986a), we derive the following expressions for coeffi-cients aijk and bijk
aijk= 1 5(aimmδjk+ ajmmδik+ akmmδij) (2.37a) bijk= 1 3bimmδjk (2.37b) Because aijk, bijkhave no direct physical interpretation, we focus on aimm, bimm, which
are directly related to the energy flux vectors. As a result, the second approximation for the distribution function can be expressed as
f(1)(r, c, ω, t) = 1 + aij 2Θ2 t CiCj+a10Θimm2 t C2 Θt − 5 Ci + Ibimm 6mΘtΘr IΩ2 mΘr − 3 Ci f (0)(r, c, ω, t) (2.38)
Inserting equations (Eqs. 2.36a-e) into the second and third moment balances, respec-tively, the second moment (Eq. 2.31) may be rewritten in terms of aij, aimm
ρDaij Dt + 1 5 ∂ ∂rj (ρaimm) + ∂ ∂ri (ρajmm) − 2 3 ∂ ∂rk (ρakmmδij) + 2ρΘtDˆij + (θkj+ ρakj) ∂vi ∂rk + (θki+ ρaki) ∂vj ∂rk −2 3(θkm+ ρakm) ∂vm ∂rk δij + ∂ ∂rk θkij− 1 3δijθkpp = χ(mCiCj) − 1 3χ(mCkCk)δij (2.39) where 2 ˆDij= ∂v∂rji + ∂vj ∂ri − 1
3∇ · v is the traceless and symmetric rate of strain tensor.
For the third moment balance equation (Eq. 2.33), let j = k (aijkp= 0) and insert
the expressions for Mij, Mijk, Mijkn in equations (Eqs. 2.36a-c) to obtain
ρDaikk Dt + 5 ∂(ρΘ2 t) ∂ri + ∂ ∂rl θlikk+ 7 ∂ ∂rl (ρΘtali) +1 5ρ 7aknm ∂vi ∂rk + 2ainn ∂vm ∂rm + 2aknn ∂vk ∂ri + ∂vi ∂rm θmkk+ 2 ∂vk ∂rm θmik − (5Θtδni+ 2ani) ∂ ∂rn (θnk+ ρΘtδnk+ ρank) = χ(CiCkCk) (2.40)
For the angular third moment (Eq. 2.34), let i = j (bijkn = 0) and insert equations
(Eqs. 2.36a, d, e) nIDbjjk Dt + ∂ ∂rn (θn(IΩjΩjCk) + δknρΘtΘr+ ρaknΘr) −δknΘr ∂ ∂rm (θmn+ ρ(Θtδmn+ amn)) + 2δinθm(IΩjCk) ∂ωn ∂rm +δkn ∂vn ∂rm
(nIbjjm+ θm< IΩjΩj>) = χ(IΩjΩjCk) (2.41)
Equations (2.39, 2.40, 2.41) are used to determine ρ, ui, Θt, Θr, aij, aijk, bijk. Some
assumptions regarding the order of magnitude of the spatial derivatives of the mean fluid field and granular temperature must be made to obtain a solution. Let L, U , Θ0
be the characteristic values of length, velocity, and granular temperature, and then use these to form dimensionless quantities proportional to the spatial gradients of granular temperature and mean velocity. In order to rewrite the moment balance equations in dimensionless forms, ρ, aij, aijk, bijk are divided by typical values ρ0, a(2), a(3), b(3).
We also make the assumption that σ/L, a(2)/Θ0, a(3)/Θ 3/2
0 , b(3)/Θ 3/2
0 are small and
of the same order of magnitude. The non-dimensional granular temperature gradient, the non-dimensional gradient of mean velocity and U/√Θ0 are all of order one.
Following Jenkins and Richman (1986a), the expression for χ(ψ) and θi(ψ) are
still linear in the perturbation and the spatial gradients of velocity and granular tem-peratures. The terms quadratic in aij, aijk and bijk that appear in the product of
f(1)(r, c1, ω1) with f(1)(r, c2, ω2) can be neglected when compared with linear terms,
while the products of aij, aijk and bijk with the spatial gradients of the mean fields
will be also ignored. Thus the flux and source terms can be written as
θi(ψ) = Ai(ψ) + Bi(ψ) + ajkBijk(ψ) + ajkpBijkp(ψ) + bjkpCijkp(ψ) (2.42)
χ(ψ) = E(ψ) + F (ψ) + ajkFjk(ψ) + ajkpFjkp(ψ) + bjkpDjkp(ψ) (2.43)
Using definitions f01= f(0)(r, c1, ω1, t), f02= f(0)(r, c2, ω2, t), we find for the flux
term Ai(ψ) = − σ3g0 2 Z G·k>0 (ψ0− ψ)G · kf01f02kidkdc1dω1dcdω (2.44a) Bi(ψ) = − σ4g 0 4 Z G·k>0 (ψ0− ψ)G · kf01f02k · ∇ ln f02 f01 kidkdc1dω1dcdω (2.44b) Bijk(ψ) = − σ3g 0 4 Z G·k>0 (ψ0− ψ)G · k f01 ∂2f 02 ∂c2j∂c2k +f02 ∂ 2f 01 ∂c1j∂c1k ! kidkdc1dω1dcdω (2.44c) Bijkp(ψ) = σ3g0 12 Z G·k>0 (ψ0− ψ)G · k f01 ∂3f 02 ∂c2j∂c2k∂c2p +f02 ∂ 3f 01 ∂c1j∂c1k∂c1p ! kidkdc1dω1dcdω (2.44d) Cijkp(ψ) = σ3g0 4 Z G·k>0 (ψ0− ψ)G · k f01 ∂3f 02 ∂Ω2j∂Ω2k∂c2p +f02 ∂ 3f 01 ∂Ω1j∂Ω1k∂c1p ! kidkdc1dω1dcdω(2.44e)
For the source term ∆χ(ψ) = ψ0+ ψ10 − ψ − ψ1, we find
E(ψ) = σ 2g 0 2 Z G·k>0 ∆χ(ψ)G · kf01f02dkdc1dω1dcdω (2.45a) F (ψ) = σ 3g 0 4 Z G·k>0 ∆χ(ψ)G · kf01f02k · ∇ ln f01 f02 dkdc1dω1dcdω (2.45b) Fjk(ψ) = σ2g 0 4 Z G·k>0 ∆χ(ψ)G · k f01 ∂2f 02 ∂c2j∂c2k + f02 ∂2f 01 ∂c1j∂c1k dkdc1dω1dcdω (2.45c) Fjkp(ψ) = − σ2g 0 12 Z G·k>0 ∆χ(ψ)G · k f01 ∂ 3f 02 ∂c2j∂c2k∂c2p +f02 ∂ 3f 01 ∂c1j∂c1k∂c1p ! dkdc1dω1dcdω (2.45d) Djkp(ψ) = − σ2g 0 4 Z G·k>0 ∆χ(ψ)G · k f01 ∂ 3f 02 ∂Ω2j∂Ω2k∂c2p +f02 ∂ 3f 01 ∂Ω1j∂Ω1k∂c1p ! dkdc1dω1dcdω (2.45e)
In order to get χ(ψ) and θi(ψ), we neglect the products of perturbations, aij, aijk
The characteristic values of length, velocity, and temperature are used to form di-mensionless quantities to determine which terms can be neglected. Then some rather cumbersome integrations can be carried out.
Finally, much simplified equations can be obtained by a term-by-term order of magnitude analysis of balance laws Eqs. 2.39, 2.40, 2.41 using the above assump-tions. Together with the integrals in the linear equations Eqs. 2.42, 2.43, the resulting equation for Eq. 2.39 will be satisfied provided that
µ ˆDij = −
1
2aij (2.46) where µ is a transport coefficient (a viscosity) defined as
µ ≡ −1 M2 {2 (1 + 2(1 + e)sg0) + 6sg0 2η1 2 5(1 + η1) + 4λ2+ 6λ + 1 5(1 + λ) + (1 + 6η1)(2λ + 1)A3− 2A4 (2.47) Here and in the following λ is defined as the granular temperature ratio, which for spheres is λ = 5Θr/2Θt, and η1 = −1+e2 , η2 = −1+β70. For perfectly elastic and
frictionless spheres, the viscosity is exactly the same as the original one from KTGF of smooth spheres. The expressions for A3, A4 (as well as the next A5, A6, A7, A9,
A11, A12, and A13) can be found in Appendix B, while the expression for M2(as well
as the next N1, N2, L1, L2 and L3) can be found in Appendix A. Note that Eq. 2.46
implies that the anisotropic deviations in the second moments of the particle velocities are controlled by the traceless rate of deformation of the particle velocity field. The expression of aijk can be determined from the third moment balance in a similar way,
leading to κt ∂Θt ∂ri =1 2aipp (2.48) Here, a transport coefficient for translational energy diffusion, similar to the coefficient of thermal conductivity for dilute gases is defined
κt≡ ρΘt
N2
2N1
(2.49) Again we find that for perfectly elastic and frictionless spheres without rotation, the above expression has the same form as the original model.
For rotational energy diffusion we find κr
∂Θr
∂ri
=1
Here, κr is the transport coefficient connected to rotational energy transport, defined as κr≡ ρ L3 2L1 (2.51) In the case of perfectly elastic and frictionless spheres without rotation, the (pseudo) rotational thermal conductivity is 0, corresponding to the original model. With these equations we are finally ready to cast our balance equations into more familiar forms.
2.2.5
Constitutive equations
The full expressions for the constitutive equations, including both kinetic and collisional parts, are listed in Table 2.1. Because the particles have rotational degrees of freedom, during particle collisions not only linear momentum is transferred but also angular momentum. These degrees of freedom are not necessarily equilibrated; the rate with which equilibration is achieved will be determined by the rotational viscosity coefficient. By identifying terms connected to the divergence of the average particle velocity field, and terms connected to the symmetric and anti-symmetric parts of the velocity gradient field, the expression for the solid stress tensor can be written as,
τs= −{(λs− 2 3µts)(∇ · vs)I + µts[∇vs+ (∇vs) T] + µ rs[∇vs− (∇vs) T]} (2.52)
Here, µts, µrs are shear (translational) and rotational viscosities and λsis the bulk
vis-cosity. For perfectly elastic and frictionless spheres without rotation, µrsbecomes zero,
and the expression of µts is the same as the original model. For inelastic and frictional
spheres, the present model is not the same as the old one since we consider particle friction, i.e. µts depends on the granular temperature ratio, which is not zero when
particle rotation exists. The rotational viscosity µrs is a consequence of the internal
ro-tational degrees of freedom of the particles: when particles collide with each other, the momentum transfer due to tangential collisions will be biased when the (average) parti-cle velocity field is rotating, i.e. when there is a non-zero rate of rotation ∇vs−(∇vs)T.
The balance equations of mass, momentum, and energy are used to simulate the flow of granular materials. For the translational energy equation, the first term on the RHS is the production term due to transport of stress by velocity fluctuations and collisions. The second term is the diffusion of energy. The third term is the dissipation of energy due to inelastic collisions. The last term is the energy exchange between the fluid and the solids phase. For the rotational energy equation, we neglect the rotational interphase interaction (rotational hydrodynamic torque). Besides, we made the assumption that in the bulk the mean rotational velocity is zero. As a result, there is no production term, but there is diffusion of rotational energy (first term on the RHS) and a collisional source/sink of granular temperature (second term on the RHS). A summary of the kinetic theory used in this work is shown in Table 2.1, while the summary of the closures in classical kinetic theory is given in Table 2.2.
Table 2.1: Constitutive equations of kinetic theory. Balance equations for fluid (f ) and solid(s):
∂(fρf) ∂t + ∂ ∂r· (fρfvf) = 0 ∂(fρfvf) ∂t + ∂ ∂r· (fρfvfvf) = −f∇Pf− ∇ · (fτf) + fρfg − βA(vf− vs) ∂(sρsvs) ∂t + ∂ ∂r· (sρsvsvs) = −s∇Pf − ∇ · (PsI + sτs) + sρsg + βA(vf − vs) 3 2 h∂( sρsΘt) ∂t + ∇ · (sρsvsΘt) i = −∇vs: (PsI + sτs) + s∇ · (κt∇Θt) − γt− 3βAΘt 3 2 h∂( sρsΘr) ∂t + ∇ · (sρsvsΘr) i = s∇ · (κr∇Θr) − γr
Solid pressure tensor: Ps= sρsΘt(1 + 2(1 + e)sg0)
Solid bulk viscosity: λs=43sρsσg0(1 + e)
q
Θt
π
Solid stress tensor:
τs= −{(λs−23µts)(∇ · vs)I + µts[∇vs+ (∇vs)
T] − µ
rs[∇vs− (∇vs)
T]}
Translational energy dissipation rate: γt= g0ρsΘt2s h −192 σ q Θt π (η1(1 + η1) − (2λ + 1)A1+ (1 + λ)A2) +12∇ · vs(η1(1 + η1) + 5(λ + 1)A4− (1 + 2λ)A3) i Rotational energy dissipation:
γr= g0ρsΘt2s h −96 σ q Θt π 5 2A2− λA1 + 120∇ · vs 5 2A4− λA3 i Shear viscosity: µts = µ(1 + µts,c) + 3 5λs, µts,c = −0.8sg0[2η1− 6(1 + 2λ)A1] Rotational viscosity: µrs = −8(1 + 2λ)σg0ρssA1 q Θt π
Pseudo translational thermal conductivity: κt= κt(1 + κts,c) +
3
2λs, κts,c = −sg0[2η1− 16(1 + 2λ)A1]
Pseudo rotational thermal conductivity: κr= κr(1 + κrs,c), κrs,c = 16λ(1 + λ)sg0 Θt Θr q Θt π A1
Gas particle drag relation (Beetstra et al. (2007)): βA= 180µg2s σ2g + 18µgs2g(1+1.5 √ s) σ2 + 0.413µgRess (4/3)σ2g h−1 g +3sg+8.4Re−0.343s 1+103sRe−0.5(1+4s)s i Radial distribution function at contact (Ma and Ahmadi (1986)): g0= 1 + 4s
1+2.5s+4.59042s+4.5154393s
[1−(s/maxs )3] 0.67802
For the rotational energy change per collision, the rate of dissipation due to inelastic collisions becomes γr= g0ρsΘt2s " −96 σ r Θt π 5 2A2− λA1 + 120∇ · vs 5 2A4− λA3 # (2.53) This vanishes when e = 1. In general, in the limit of perfectly elastic and perfectly smooth spheres (e → 1, µ → 0), all the closures reduce to the classical results of the first-order approximation from the dense smooth hard-sphere kinetic theory (Chapman and Cowling (1970)).
Table 2.2: Summary of the classical kinetic theory closures (Gidaspow (1994)). Translational energy dissipation rate:
γt= 3(1 − e2)2sρsg0Θt 4 σ q Θt π − ∇ · vs Shear viscosity: µs= 1.016965πρsσ q Θt π (1+8 5 1+e 2 sg0)(1+85sg0) sg0 + 4 5sρsσg0(1 + e) q Θt π
Pseudo translational thermal conductivity: κs= 1.0251338475πρsσ q Θt π (1+12 5 1+e 2 sg0)(1+125sg0) sg0 + 2sρsσg0(1 + e) q Θt π (a) (b)
Figure 2.2: Comparison of shear viscosity as a function of friction coefficient for different (σ=3 mm, s= 0.4, Θt=0.003 m2/s2).
Figure 2.2 shows the shear viscosity as a function of friction coefficient for different values of the granular temperature ratio λ for a system containing 3 mm particles at a
