UFDC Home  myUFDC Home  Help 
Material Information
Thesis/Dissertation Information
Subjects
Notes
Record Information

Material Information
Thesis/Dissertation Information
Subjects
Notes
Record Information

Full Text  
PAGE 1 1 COMPUTATIONAL STUDY OF HEAT TRANSFER IN FLUIDSOLID FLOWS By LIKE LI 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 2013 PAGE 2 2 2013 Like Li PAGE 3 3 To my parents and my wife PAGE 4 4 ACKNOWLEDGMENTS I would like to express my most sincere gratitude to my advisor, Professor Renwei Mei, and my coadvisor, Professor James F. Klausner, for their constant advice, support, inspiration and encouragement This dissertation would have been impossible without their expert guidance and enthusiasm on the subject. I would like to express my sincere gratitude to my com mittee member, Professor David W. Hahn, for his concern and help with Chapter 2 about the heat conduction computation between colliding particles. I would also like to express my sincere gratitude to my committee member, Professor Jrg Petrasch, for his g uidance and help with the material presented in Chapter 6 about the coupled radiationconduction model. M y sincere gratitude also goes to Professor Jason E. Butler, from the Department of Chemical Engineering at the Universi ty of Florida for serving as an external member in my committee and for his valuable advice and suggestions regarding the particlelevel simulations. I would like to thank Abhishek Singh in our department for developing the radiation model and providing m e the required heat flux boundary conditions for the coupled radiationconduction model in Chapter 6. We have discussed and collaborated a lot and learned from each other during the process. I would also like to take this opportunity to thank my labmates: Drs. Richard Parker, Fotouh Al Raqom, Ayyoub Mehdizadeh, Nicholas AuYeung Bradley Bon, Chang kang Guan, and Fadi Alnaimat, and those who are on the way of pursuing a Doctors degree, including Nima Rahmatian, Chen Chen, Kyle Allen, Benjamin Greek, Prasanna Venuvanalingam and Amey Barde. PAGE 5 5 Finally, I am deeply grateful for the enduring and endless support and love of my wife, Jing Hu, my parents, and my brother. PAGE 6 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................. 4 LIST OF TABLES .......................................................................................................... 10 LIST OF FIGURES ........................................................................................................ 11 LIST OF ABBREVIATIONS ........................................................................................... 17 ABSTRACT ................................................................................................................... 18 CHAPTER 1 INTRODUCTION TO HEAT TRANSFER IN FLUIDSOLID FLOWS ...................... 20 1.1 Background and Motivations ............................................................................. 20 1.2 Objectives ......................................................................................................... 21 1.2.1 Heat Transfer between Colliding Solid Particles ...................................... 22 1.2.2 Thermal Boundary Conditions on Curved Boundaries ............................. 22 1.2.3 Heat Transfer on Curved Boundaries ...................................................... 24 1.2.4 MultipleRelaxationTime Lattice Boltzmann Model for the Axisymmetri c ConvectionDiffusion Equation ................................................ 25 1.2.5 Thermal Analysis in a Porous Structure in a Solar Thermochemical System .......................................................................................................... 26 2 HEAT TRANSFER BETWEEN C OLLIDING SURFACES AND PARTICLES ......... 28 2.1 Introduction ....................................................................................................... 28 2.2 Formulation of the Collisional Heat Transfer Problem ....................................... 30 2.2.1 Elastic Collision and Deformation ............................................................ 30 2.2.2 Governing Equations for Heat Transfer ................................................... 32 2.2.3 Self Similar Solution ................................................................................ 35 2.2.4 Heat Transfer .......................................................................................... 36 2.3 Results and Discussion ..................................................................................... 37 2.3.1 Self Similar Solution at Small Time ......................................................... 37 2.3.2 Behavior of the Singularity at ( r z ) = ( R ( t ), 0) .......................................... 38 2.3.3 Asymptotic Limit at Fo 0 ...................................................................... 42 2.3.4 Numerical Results for Finite Fo ............................................................... 47 2.3.5 Particles of Different Materials ................................................................. 49 2.4 Summary and Conclusions ............................................................................... 52 3 BOUNDARY CONDITIONS FOR THERMAL LA TTICE BOLTZMANN EQUATION METHOD ............................................................................................. 64 3.1 Literature Review of Thermal Lattice Boltzmann Equation (TLBE) Models ....... 64 PAGE 7 7 3.2 Literature Review of Boundary Conditions in TLBE Simulations ....................... 66 3.3 Present Boundary Condition Treatments in TLBE Method ............................... 70 3.3.1 The ConvectionDiffusion Equation ......................................................... 71 3.3.2 The D3Q7 and D2Q5 Thermal LBE models ............................................ 72 3.3.3 Thermal Boundary Condition Treatments ................................................ 77 3.3.3.1 Dirichlet boundary condition ........................................................... 79 3.3.3.2 Neumann boundary condition ........................................................ 82 3.4 Numerical Validation and Discussion ................................................................ 89 3.4.1 TwoD Steady State Channel Flow ......................................................... 90 3.4.1.1 Dirichlet boundary condition ........................................................... 90 3.4.1.2 Neumann boundary condition ........................................................ 93 3.4.2 OneD Transient Conduction in an Inclined Semi Infinite Solid ............... 94 3 .4.3 TwoD Steady State Heat Conduction Inside a Circle ............................. 96 3.4.3.1 Dirichlet boundary condition ........................................................... 96 3.4.3.2 Neumann boundary condition ........................................................ 97 3.4.4 TwoD Transient Heat Conduction Inside a Circle ................................... 99 3.4.5 ThreeD Steady State Circular Pipe Flow .............................................. 101 3.4.5.1 Dirichlet boundary condition ......................................................... 101 3.4.5.2 Neumann boundary condition ...................................................... 1 02 3.5 Summary and Conclusions ............................................................................. 103 4 HEAT TRANSFER EVALUATION ON CURVED BOUNDARIES IN THERMAL LA TTICE BOLTZMANN EQUATION METHOD .................................................... 122 4.1 Introduction ..................................................................................................... 122 4.2 Heat Transfer Evaluation ................................................................................ 125 4.2.1 Boundary Heat Flux ............................................................................... 125 4.2.2 Heat Transfer on Curved Boundaries .................................................... 126 4.2.3 Temperature Gradient in the Interior Field ............................................. 128 4.3 Nu merical Results and Discussion .................................................................. 129 4.3.1 TwoD Steady State Channel Flow ....................................................... 130 4.3.1.1 Channel flow with a Dirichlet boundary condition ......................... 131 4.3.1.2 Channel flow with a Neumann boundary condition ...................... 133 4.3.2 OneD Transient Conduction in an Inclined Semi Infini te Solid ............. 133 4.3.2.1 Dirichlet Boundary condition f ( t ) = kt1/2 ......................................... 134 4.3.2.2 Dirichlet Boundary condition f ( t ) = 1 exp( ) ............................ 134 4.3.3 TwoD Transient Heat Conduction Inside a Circle ................................. 135 4.3.4 ThreeD Steady State Circular Pipe Flow with a Dirichlet Condition ..... 136 4.3.4 TwoD Steady State Natural Convection in a Square Enclosure with a Circular Cylinder in the Center .................................................................... 137 4.4 Summary and Conclusions ............................................................................. 141 5 MULTIPLE RELAXATION TIME LATTICE BOLTZMANN MODEL FOR THE AXISYMMETRI C CONVECTION DIFFUSION EQUATION .................................. 150 5.1 Introduction ..................................................................................................... 150 PAGE 8 8 5.2 Lattice Boltzmann Model for the Axisymmetric ConvectionDiffusion Equations ........................................................................................................... 156 5.3 Numerical Val idation and Discussion .............................................................. 159 5.3.1 Steady convectiondiffusion in an annulus ............................................ 159 5.3.2 Steady convectiondiffusion in a circular pipe ........................................ 160 5.3.2.1 Dirichlet boundary condition on the wall ....................................... 161 5.3.2.2 Neumann boundary condition on the wall .................................... 162 5.3.3 Steady heat conduction inside a sphere ................................................ 163 5.3.3.1 Dirichlet boundary condition on the wall ....................................... 164 5.3.3.2 Neumann boundary condition on the wall .................................... 165 5.3.4 Unsteady heat conduction in a circle ..................................................... 166 5.3.5 St eady natural convection in an annulus between two coaxial vertical cylinders ...................................................................................................... 167 5.3.6 Swirling thermal flows in a vertical cylinder ........................................... 169 5.4 Summary and Conclusions ............................................................................. 172 6 ENERGY TRANSPORT IN A SOLAR THERMOCHEMICAL REACTOR FOR HYDROGEN PRODUCTION ................................................................................ 185 6.1 Introduction ..................................................................................................... 185 6.2 Horizontal Cavity Receiver Solar Reactor ....................................................... 188 6.3 Heat Transfer Analysis .................................................................................... 189 6.3.1 Radiative Heat Transfer ......................................................................... 189 6.3.1.1 Radiation between the cavity wall and the absorbers .................. 189 6.3.1.2 Radiation is the porous structure inside the absorbers ................ 191 6.3.2 Convective Heat Transfer ...................................................................... 192 6.3.3 Conductive Heat Transfer ...................................................................... 192 6.3.4 Endothermic Chemical Reaction ........................................................... 192 6.4 TLBE Simulation of Heat Transfer in Absorbers ............................................. 193 6.4.1 Scaling from the physical units to the LBE units .................................... 194 6.4.2 Rescaling for highdiffusion coefficient CDE in TLBE ............................ 195 6.4.3 Computational Procedure for the Coupled Model .................................. 196 6.5 Simulation Results and Discussion ................................................................. 197 6.6 Summary and Conclusions ............................................................................. 200 7 SUMMARY AND FUTURE WORKS ..................................................................... 207 APPENDIX A VALIDITY OF THE SEMI INFINITE MEDIA ASSUMPTION ................................. 210 B ASYMPTOTIC ANALYSIS OF SECOND ORDER ACCURATE BOUNDARY SCHEMES ............................................................................................................ 215 C EXTENSION OF ZHOUS BGK AXISYMMETRIC LATTICE BOLTZMANN MO DEL (ZHOU 2011) TO AN MRT LATTICE BOLTZMANN MODEL .................. 222 LIST OF REFERENCES ............................................................................................. 226 PAGE 9 9 BIOGRAPHICAL SKETCH .......................................................................................... 236 PAGE 10 10 LIST OF TABLES Table page 2 1 Heat flux integral **()wqt based on the self similar solution at small time ( *1t ) ............................................................................................................... 53 3 1 Accuracy of the TLBE simulations with Dirichlet/Neumann conditions ............. 105 4 1 Surfaceaveraged Nusselt numbers at different grid resolution. ....................... 142 4 2 Comparison of the present surfaceaver aged Nusselt numbers with previous results. .............................................................................................................. 142 5 1 Surfaceaveraged Nusselt numbers at different grid resolution. ....................... 174 5 2 Compa rison of the present surfaceaveraged Nusselt numbers with published results. .............................................................................................. 174 6 1 Baseline simulation parameters for cavity aspect ratio Lcav, in / Dcav = 0.69, 1.01, and 1.16 for Cases I, II and III, respectively. ............................................ 201 6 2 Solar to fuel efficiency calculations for Cases I, II and III. ................................ 201 PAGE 11 11 LIST OF FIGURES Figure page 2 1 Variation of the contact area during collision. ..................................................... 54 2 2 Computational domain and boundary conditions at z = 0. .................................. 54 2 3 Temperature contours on the surface of particle 1 based on the self similar solutions at small time ( *1 t ) .......................................................................... 55 2 4 Dimensionless heat flux across the contact area at small time ( *1t ) based on self similar solutions. ..................................................................................... 56 2 5 (a) Local solution of ( 0.5)/( 1/2) around the singularity point at ( ) = (1, 0). (b) The value of w in the local solution of around the singularity point at ( ) = (1, 0). ...................................................................................................... 57 2 6 Local transient temperat ure profiles at ( > 1, = 0) when Fo = 1.0. .................. 58 2 7 Compasion of numerical results and analytical fit for fi( ) in Eq. (2 58). ............. 58 2 8 Variation of the dimensionless heat flux integral **()wqt over the contact area at Fo = 0. ............................................................................................................ 59 2 9 Variation of the multiplier **/At in Eqs. (2 30b) and (231c). .......................... 59 2 10 Variations of the dimensionless heat flow rate **()Qt at Fo = 0; the self similar solution is **()Qt = 0.2820948 **/ At ................................................... 60 2 11 Comparison of the dimensionless heat flux integral **()wqt at different Fourier numbers. ............................................................................................................ 60 2 12 Heat flux **(,) qt in the contact area at various times for Fo = 1.0. ................... 61 2 13 Temperature variations along direction ( = 0) at various times for Fo = 1.0. ..................................................................................................................... 61 2 14 Heat transfer as a function of the impact Fourier number (particles of the same material). ................................................................................................... 62 2 15 Correction factor comparison between computation and the approximate expression at ( a ) k2/ k1=1.0, ( b ) k2/ k1=0.1 and ( c) k2/ k1=10.0. .............................. 63 3 1 Discrete velocity set { e} for the D3Q7 TLBE model. ........................................ 106 PAGE 12 12 3 2 Schematic depiction of the lattice link intersected by a boundary wall. ............. 106 3 3 Layout of th e rectangular lattice and a curved boundary (solid circles: field nodes; solid squares: boundary nodes; open circles: exterior nodes). ............. 106 3 4 Schematic representation of the lattice in a 2D value. ................................................................................................................ 107 3 5 (a) Lower stability bound of the adjustable parameter cd 1 versus ( D 0.5), and (b) variations of cd 1 as a function the lattice fraction. ................................. 107 3 6 Temperature contours in a 2D channel with a Dirichlet boundary condition using D H = 64. ................................................................. 108 3 7 Relative L 2 error norm E2 versus the grid resolution, 1/ H for the channel flow Dirichlet problem. ...................................................................................... 109 3 8 Relative L 2 error norm E2 channel flow Dirich let problem at Ny = 34. ....................................................... 110 3 9 Comparison of E2 obtained using Scheme 2 with those using Ginzburgs schemes (Ginzburg 2005) for the channel flow Dirichlet problem. .................... 110 3 10 Wall heat flux error E2_qw versus the grid resolution, 1/ H for the channel flow D = 0.75, (b) and 0.99 at D = 0.75. ....................................................................................... 111 3 11 Interior gradient error E2_grad versus the grid resolution, 1/ H for the channel D = 0.75, (b) 0.01, and 0.99 at D = 0.75. ............................................................................... 111 3 12 Temperature contours in a 2D channel with a Neumann boundary condition using D H = 64 ............................................................... 112 3 13 Relative L 2 er ror norm E2 versus the grid resolution, 1/ H for the channel D = 0.51, (b) D = 0.55, and (c) D = 0.75. .............................................................................................................. 113 3 14 Wall temperature error E2_tw versus the grid resolution, 1/ H for the channel flow Neumann problem. (a) = 0.50, 0.25, and 0.75 at D = 0.75, (b) = 0.50, 0.01, and 0.99 at D = 0.75. ...................................................................... 114 3 15 Interior gradient error E2_grad v ersus the grid resolution, 1/ H for the channel flow Neumann problem. .................................................................................... 114 3 16 Layout of the lattice around the inclined semi infinite solid. .............................. 115 PAGE 13 13 3 17 (a) Temperature variations with time at point P ; (b) temperature profiles in the ydirection going through P at t = 500 ........................................................... 115 3 18 along the azimuthal angle at r0 = 30.5. .......................................................... 116 3 19 Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Dirichlet condition (results based on Ginzburgs schemes (2005) are included). .................................................. 116 3 20 Wall heat flux error E2_qw versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Dirichlet condition (results based on Ginzburgs schemes (2005) are included). ....................................................... 117 3 21 Interior gradient error E2_grad versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Dirichlet condition (results based on Ginzburgs schemes (2005) are included). .................................................. 117 3 22 Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Neumann condition. ...................... 118 3 23 Wall temperature error E2_tw versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Neumann condition. ...................... 118 3 24 Interior g radient error E2_grad versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Neumann condition. ...................... 119 3 25 Comparison of the error norm E ( t*) defined in Eq (3 71) for the transient heat conduction inside a circle using Schemes 1, 2 and 3 with r0 = 5.8. .................. 119 3 26 Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the transient heat conduction inside a circle. ......................................................................... 120 3 27 Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the pipe flow with a Dirichlet condition. .................................................................................. 120 3 28 Wall heat flux error E2_qw versus the grid resolution, 1/ r0, for the pipe flow with a Dirichlet condition. ......................................................................................... 121 3 29 Relative L 2 norm error E2 versus the grid resoluti on, 1/ r0, for the pipe flow with a Neumann condition. ............................................................................... 121 4 1 Illustration of the heat transfer evaluation on a curved boundary (closed circles: field nodes; closed squares: boundary nodes; open circles: exterior nodes). ............................................................................................................. 143 PAGE 14 14 4 2 Wall heat flux error E2_qw versus the grid resolution, 1/ H for the channel flow and 0.99. ................................................................................................................. 143 4 3 Wall heat flux error E2_qw Ny = 34 for the channel flow Dirichlet problem. ......................................................................... 144 4 4 Wall heat flux error E2_qw versus ( D 0.5) at Ny = 66 for the channel flow Dirichlet problem. .............................................................................................. 144 4 5 Interior gradient error E2_gx (in xdirection) versus the grid resolution, 1/ H for 0.50, 0.01, and 0.99. ......................................................................................... 145 4 6 Interior gradient error E2_gy (in ydirection) versus the grid resolution, 1/ H for 0.50, 0.01, and 0.99. ......................................................................................... 145 4 7 Interior gradient errors of (a) E2_gx, and (b) E2_gy, respectively, vers us the grid resolution, 1/ H for the channel flow Neumann problem. .................................. 146 4 8 Variations of the heat transfer rate with time for the Dirichlet condition T ( l =0)= t1/2 at the end of an inclined sem i infinite solid (symbols: LBE results; lines: exact solutions). ................................................................................................ 146 4 9 Variations of the heat transfer rate with time for the Dirichlet condition T ( l = 0) = 1 exp( 0.1 t ) at the end of an i nclined semi infinite solid (symbols: LBE results; lines: exact solutions). .......................................................................... 147 4 10 Variations of the error E ( t*) defined in Eq. (4 16) with time at different grid resolution for the transient heat conduction inside a circle. .............................. 147 4 11 Heat transfer error, E2, defined in Eq. (4 17) versus the grid resolution, 1/ r0, for the transient heat conduction inside a circle. ............................................... 148 4 12 Heat transfer error 2wQE_ versus the grid resolution, 1/ r0, for the circular pipe flow with a Dirichlet boundary condition. ........................................................... 148 4 13 Schematic depiction of the computational domain ( r0 = 0.2 L ) and the Dirichlet thermal and velocity boundary conditions on the inner and outer walls. ........... 149 4 14 Isotherms in t he field between the circular cylinder and the cavity walls at (a) Ra = 103, (b) Ra = 104, (c) Ra = 105, and (d) Ra = 106. .................................... 149 4 15 Streamlines in the field between the circular cylinder and t he cavity walls at (a) Ra = 103, (b) Ra = 104, (c) Ra = 105, and (d) Ra = 106. .............................. 149 PAGE 15 15 5 1 Discrete velocity set { e} for the D2Q5 LB model. e = (0, 0), (1, 0), and (0, 1). ................................................................................................................... 175 5 2 Relative L 2 error norm E2 for the interior field of versus the grid resolution, 1/ Ri, for the 1 D annulus convectiondif fusion problem. ................................... 175 5 3 Schematic of the square lattice in the r z plane of a 3D circular pipe ............. 176 5 4 Contours of in the r z plane at Nr = 42, Nz convectiondiffusion in a pipe with a Dirichlet boundary condition on the wall. 176 5 5 Relative L 2 error norms versus the grid resolution, 1/ R0, for the steady convectiondiffusion in a pipe with a Dirichlet boundary condition. ................... 177 5 6 Relative L 2 error norms versus the grid resolution, 1/ R0, for the steady convectiondiffusion in a pipe with a Neumann boundary condition. ................ 178 5 7 Schematic representation of the computational domain and the lattice layout for the heat conduction inside a sphere. ........................................................... 179 5 8 Relative L 2 error norms E2 and E2_qw for the interior field of and the boundary flux, respectively, versus the grid resolution, 1/ R0, for the steady heat conduction inside a sphere with a Dirichlet boundary condition. ............... 179 5 9 Relative L 2 error norms E2 and E2_tw for the interior field and boundary values of respectively, versus the grid resolution, 1/ R0, for the steady heat conduction inside a sphere with a Neumann boundary condition. .................... 180 5 10 Time averaged error norm E2 for the interior field of versus the grid resolution, 1/ R0, for the unsteady heat conduction in a circle with a temporal D irichlet boundary condition. ............................................................................ 180 5 11 Schematic depiction of the computational domain ( Ro/ Ri = 2.0, and H /( Ro Ri) = 2.0) and the thermal boundary conditions on the walls all placed in the ................................................................... 181 5 12 Streamlines in the annulus between two coaxial vertical cylinder s at (a) Ra = 103, (b) Ra = 104, and (c) Ra = 105. .................................................................. 182 5 13 Isotherms in the annulus between two coaxial vertical cylinders at (a) Ra = 103, (b) Ra = 104, and (c) Ra = 105. .................................................................. 182 5 14 Contours of the azimuthal velocity u (a) and isotherms (b) of the confined swirling flow at Re = 2000 and Ri = 0. .............................................................. 183 5 15 Contours of the azimuthal velocity u (a) and isotherms (b) of the confined swirling flow at Re = 2000 and Ri = 1.0. ........................................................... 183 PAGE 16 16 5 16 Velocity components and temperature profiles along the vertical line at r / R = 0.8 for Re = 2000: (a) radial velocity ur, (b) azimuthal velocity u, (c) axial velocity uz, and (d) temperature. ....................................................................... 184 6 1 Schematical depiction of the horizontal solar reactor configuration: (a) cross section view, and (b) front view. ....................................................................... 202 6 2 (a) Variations of the average bed temperature ( Tbed) with time, and (b) variations of Tbed and solar power input ( solarQ ) with time at 1200 o Tbed 1450 oC for Cases I, II and III. .......................................................................... 202 6 3 I ncident heat flux distribution (in W/m2) on the surface of Absorber 1 (see Fig. 6 1) for (a) Cases I, (b) Case II, and (c) Case III. ............................................. 203 6 4 Net heat flux distribution (in W/m2) on the surface of Absorber 1 for (a) Cases I, (b) Case II, and (c) Case III. .......................................................................... 204 6 5 Temperature distribution (in K) on the surface of Absorber 1 for (a) Cases I, (b) Case II, and (c) Case III. ............................................................................. 205 6 6 Comparison of t emperature profiles (in K) along ( = 0, z ) on the surface of Absorber 1 for the three cases simulated. ........................................................ 205 6 7 Interior temperature distribution on the central xz plane of Absorber 1 for (a) Cases I, (b) Case II, and (c) Cas e III. ............................................................... 206 A 1 Isothermals on the surface of particle 1 at Fo = 50 when (a) t = 0.50, (b) t = 0.75, and (c) t = 0.95. ...................................................................................... 214 PAGE 17 17 LIST OF ABBREVIATIONS CDE convectiondiffusion equation CFD computational fluid dynamics LB lattice Boltzmann LBE lattice Boltzmann equation MCRT Monte Carlo ray tr acing MRT multiple relaxationtime SRT single relaxationtime TDF temperature distribution function TLBE thermal lattice Boltzmann equation PAGE 18 18 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Parti al Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONAL STUDY OF HEAT TRANSFER IN FLUIDSOLID FLOWS By Like Li August 2013 Chair: Renwei Mei Cochair: James F. Klausner Major: Mechanical Engineering Efficient and ac curate heat transfer evaluation is essential to the study of thermal flows and energy exchange in fluidsolid systems. For fluidsolid flows with moving solid particles, heat transfer between the colliding particles is investigated computationally using a finite difference method and analytically using various asymptotic methods. An efficient heat transfer evaluation technique for any curved boundary has also been proposed and validated in this dissertation based on the thermal lattice Boltzmann equation (T LBE) method with an improved boundary condition treatment. Two specific topics are further studied with the proposed boundary schemes and heat transfer evaluation technique, including the development of a multiplerelaxationtime (MRT) lattice Boltzmann (LB) model for the axisymmetric convectiondiffusion equation (CDE), and a radiationconduction coupled model for the energy transport in a solar thermochemical reactor. From the computationally determined heat transfer during the collision, a closedform fo rmula is developed to predict the heat transfer as a function of the impact Fourier PAGE 19 19 number, the thermal diffusivity ratio and the thermal conductivity ratio of the impacting particles. To develop a general and efficient technique for evaluating the heat tr ansfer on a curved boundary, the TLBE method is applied for its simple algorithm and capability to deal with complex geometries. An improved thermal boundary condition treatment is proposed based on the bounceback and interpolation of the temperature di stribution functions (TDFs) of the TLBE model. The exact geometry on a boundary is preserved with the proposed boundary condition treatment so that the boundary heat fluxes in the discrete velocity directions of the TLBE model can be directly obtained from the TDFs at the neighboring lattice nodes. The total heat transfer on the boundary is simply a summation of all the discrete heat fluxes with a constant surface area. As applications of the proposed boundary schemes and heat transfer evaluation technique in the TLBE method, 1) a quasi two dimensional (2D) axisymmetric LB model with an MRT collision operator is proposed to simulate the general axisymmetric CDE in 3 D, and 2) a radiationconduction coupled model is developed to simulate the energy transport in a solar reactor. PAGE 20 20 CHAPTER 1 INTRODUCTION TO HEAT TRANSFER IN FLUIDSOLID FLOWS 1.1 Background and Motivations Transport phenomena including momentum, heat and mass exchange are routinely investigated in many engineering systems such as heat exchangers, heat pipes, electronic cooling devices, fluidized beds and chemical reactors. Efficient modeling tools are required for system design and optimization purposes (Faghri and Zhang 2006) There has also been growing interest in simulating the fluid flow s and heat and mass transfer on the porescale in porous media, micro channels, and micro reactors (Sullivan et al. 2005, Manjhi et al. 2006, Verma et al. 2007, Verma et al. 2008, Kao et al. 2007, Garimella and Sobhan 2003, Kang et al. 2006, Kang et al. 2007, W ang et al. 2007, Chen et al. 2012, Maier and Bernard 2010) due to the growth of supercomputing power and the need to gain a fundamental understanding of the local transport processes. The present dissertation is concentrated on the computational study of h eat transfer in fluidsolid systems. Heat transfer between two impacting solids, such as particleto particle and particle to wall, initially at different temperatures is the first topic addressed. For most engineering applications, the collision period is very short and the contact area is much smaller than the characteristic length of the solids (Sun and Chen 1988, Li et al. 2012). Under the assumptions of semi infinite media and no thermal resistance, the contact area is characterized by the Hertzian elastic theory. With a coordinate transformation employed the computational domain becomes fixed for this moving boundary problem and the heat transfer through the contact area is conveniently evaluated. PAGE 21 21 For heat transfer evaluation on an arbitrarily shaped c urved boundary, a general, efficient and accurate technique that preserves the local geometry is desirable. The thermal lattice Boltzmann equation (TLBE) method, which is a direct extension of the lattice Boltzmann equation (LBE) for fluid dynamics simulat ions, has proved to be a promising candidate considering the inherent benefits of the LBE method and the continuous development of effective TLBE models ( Alexander et al. 1993, Chen et al. 1994, Shan 1997, He et al. 1998, Lallemand and Luo 2003, Guo et al. 2007, Shi and Guo 2009, Mezrhad et al. 2010, Yoshida and Nagaoka 2010). The evolution of the microscopic temperature/energy distribution functions is simulated in the TLBE method and the macroscopic temperature is obtained from the moment of the temperature distribution functions (TDFs). Accurate and robust boundary conditions for the TDFs that conserve the macroscopic Dirichlet and Neumann conditions are essential to the integrity of the TLBE simulation and the ov erall heat transfer evaluation. 1.2 Object ives The main objective of the present research is to investigate the thermal field and heat transfer in fluidsolid systems. To gain a fundamental understanding of the local thermal behavior, a computational study on the scale of the impact solid particles is conducted to analyze the heat transfer between colliding particles and surfaces initially at different temperatures, and a kinetic based mesoscopic lattice Boltzmann equation (LBE) method is applied to investigate the heat transfer on the fluidsolid interface and the temperature distribution inside a porous structure in a solar thermochemical system Specifically, the following topics are studied in detail. PAGE 22 22 1.2.1 Heat Transfer between Colliding Solid Particles There have been continuous interest and e ffort in evaluating and predicting the transient heat transfer across the interface between two contacting solids initially at different surface temperatures (Heasley 1965, Barber 1989, Ammar et al. 1992, Li and Mason 2000 Kostoglou and Konstandopoulos 2002). This transient problem is encountered in frictional heating where surface temperatures are different (Heasley 1965), in thermal machine components whose surfaces are meeting and separating periodically, such as the exhaust valves of internal combustion engines (Howard and Sutton 1970), and in the heat exchange between impacting particles in fluidized beds, where the heat conducted during the deformation of the impacting bodies must be considered (Sun and Chen 1988, Kaviany 1988). In this work we aim at obtaining an accurate closedform expr ession to predict the heat transfer between two colliding particles as a function of the impact Fourier number, the thermal diffusivity ratio, and the thermal conductivity ratio of the particles (Li et al. 2012a) 1.2 .2 Thermal Boundary Conditions on Curved Boundaries The kinetic based lattice Boltzmann equation (LBE) method has become an alternative and attractive numerical approach for fluid dynamics simulations because it is explicit in formulation, easy to implement and compatible with parallelization compared with solving the Navier Stokes equations. It has been successfully applied in a variety of fluid flow problems (Qian et al. 1992, Chen and Doolen 1998, Succi 2001, Yu et al. 2003) and has been extended to complex flows, such as multi phase and multicomponent flows ( Luo and Girimaji 2003 Zheng et al. 2006, Premnath and Abraham 2007, Asinari 2008, Aidun and Clausen 2010), turbulent flows (Ltt et al. PAGE 23 23 2006, Yu et al. 2006, Aidun and Clausen 2010), micro flows (G uo et al. 2006 Aidun and Clausen 2010), and flows through porous media ( Guo and Zhao 2002a Pan et al. 2006). Because of the advantages of the LBE method, various thermal LBE (TLBE) models have been developed and successfully used for thermal transport pr oblems (Alexander et al. 1993, Chen et al. 1994, Shan 1997, He et al. 1998, Lallemand and Luo 2003, Guo et al. 2007, Shi and Guo 2009, Mezrhad et al. 2010, Yoshida and Nagaoka 2010). Accurate boundary condition treatment is essential for the overall integr ity of the LBE simulations. For the fluid velocity field, the bounceback scheme has been intensively studied (Ziegler 1993, Noble et al. 1995, Ginzbourg and dHumires 1996, Chen et al. 1996, He et al. 1997, Zou and He 1997, Ginzburg and dHumires 2003) and extended to curved boundaries ( Filippova and Hnel 1998 Mei et al. 1999, Mei et al. 2000, Bouzidi et al. 2001, Guo et al. 2002). The standard bounceback method is second order accurate for noslip wall boundaries when the distance from the wall to the nearest fluid node is half of the lattice unit spacing. For other hydrodynamic boundary conditions such as given density, pressure, slipping velocity or their derivatives, the bounceback scheme can be modified with the flow properties at on the boundary included in the treatments ( Ginzbourg and dHumires 1996, Chen et al. 1996, He et al. 1997, Zou and He 1997 Filippova and Hnel 1998, Mei et al. 1999, Mei et al. 2000, Bouzidi et al. 2001, Guo et al. 2002, Ginzburg and dHumires 2003). Spatial interpolation or extrapolation has also been incorporated to improve the degraded accuracy of the bounceback scheme when the boundary is located in arbitrary position in the lattice cell (Chen et al. 1996, Filippova and Hnel 1998, Mei et al. 1999, Mei et al. PAGE 24 24 2 000, Bouzidi et al. 2001, Guo et al. 2002). With appropriate modifications, some of the athermal boundary condition treatments have been extended to TLBE simulations as well. With different types of thermal boundary condition treatments available in the literature, we notice that discussions about the applicability and accuracy of these treatments in curvedboundary situations are rather limited. The objectives of the present work are twofold. The first is to develop secondorder accurate boundary condition treatments for both the Dirichlet and Neumann conditions on straight boundaries that can be placed at arbitrary locations in the lattice. The second is to extend the present boundary condition treatments to curvedboundary situations and to investigate t he effect of curvedgeometry on the accuracy of the present boundary condition treatments. 1.2.3 Heat Transfer on Curved Boundaries With conventional computational fluid dynamics (CFD) techniques, heat transfer on the complex fluidsolid interface is freq uently simulated with a multi scale analysis in order to capture the local thermal behavior on the interface. Accurate integration is often required to determine the total heat transfer across the boundary. The mesoscopic nature of the LBE method makes it very efficient in examining the local flow and momentum exchange at the fluidsolid interface. Mei et al. (2002) investigated two different methods for evaluating the fluid force on the fluidsolid boundary based on momentum exchange and stress integration, respectively, in the lattice Boltzmann equation. The momentum exchange approach was found to be reliable, accurate and much easier to implement than the stress integration approach for both twodimensional and threedimensional flows. PAGE 25 25 Following the cons ervation and exchange ideas for force evaluation in the LBE method (Mei et al. 2002), an energy exchange approach for heat transfer evaluation on the fluidsolid interface in the TLBE is proposed in this work. The boundary heat flux in the discrete velocity directions of the TLBE model is directly obtained from the temperature distribution functions at the lattice nodes, rather than using any finitedifference schemes with the obtained temperature field. The heat fluxes in the discrete velocity directions can be integrated over the lattice surface area to obtain the heat transfer on the interface, thus avoiding the calculations of the normal heat flux on the boundary and the surface area approximations 1.2.4 Multiple Relaxation Time Lattice Boltzmann Model for the Axisymmetric ConvectionDiffusion Equation Although 3D LB models can be directly applied to simulate the axisymmetric CDE when the curved boundary condition is properly handled, quasi 2 D LB models, which r equire far less computational effort and bypass the curved boundary treatment for straight pipes, are attractive alternatives for solving the axisymmetric CDE as demonstrated by the LB models developed for axisymmetric athermal flows (Halliday et al. 2001, Lee et al. 2006, Reis and Phillips 2007, Reis and Phillips 2008, Zhou 2008, Chen et al. 2008, Huang and Lu 2009, Guo et al. 2009, Wang et al. 2010, Li et al. 2010, Zhou 2011) and thermal problems (Peng et al. 2003, Huang et al. 2007, Chen et al. 2009, Li et al. 2009, Zheng et al. 2010a, Zheng et al. 2010b). It is noted that all the existing LB models for axisymmetric thermal flows in (Chen et al. 2009, Li et al. 2009, Zheng et al. 2010a, Zheng et al. 2010b) used the BGK collision operator. A robust and str aightforward MRT LB model that can simulate the general axisymmetric CDE is PAGE 26 26 highly desired. Furthermore, we notice that the numerical validation regarding the order of accuracy of the axisymmetric thermal LB models has not been well established. In this wo rk, we propose an MRT LB model for the general axisymmetric CDE for scalar transport such as temperature, mass concentration and azimuthal velocity component in rotational flows. When the CDE is coupled with the hydrodynamic equations, the radial axial vel ocity field is solved using Zhous axisymmetric LB model (Zhou et al. 2011) with the BGK collision model replaced by an MRT based model. 1.2. 5 Thermal Analysis in a Porous Structure in a Solar Thermochemical System The energy transport in a solar thermoche mical high temperature system is numerically studied and the temperature distribution in the porous reacting bed inside the absorbers sitting inside a cavity receiver is simulated. We apply the Monte Carlo ray tracing technique to analyze the radiative transfer between the cavity wall and the absorber surfaces. The temperature field inside the absorbers is modeled with the TLBE method and the proposed thermal boundary condition treatments for curved boundaries; endothermic chemical reaction is treated as a heat sink based on the reaction rate. These two modules are explicitly coupled with the reemission from the absorber surfaces included in the radiative transfer analysis. An optimal design of the cavity reactor can be pursued by investigating the solar en ergy absorption efficiency and the temperature distribution inside the absorbers with different configurations. The present dissertation is constructed as follows. Chapter 2 presents the heat transfer computation between colliding solid particles and surfaces with a closedform formula for heat transfer prediction obtained in the end. Chapter 3 reviews the boundary conditions in TLBE method in the literature and improved boundary condition treatments are proposed for the Dirichlet, Neumann and mixed thermal boundary conditions. The PAGE 27 27 stability and accuracy of the proposed boundary schemes are validated with detailed numerical tests. Based on the improved boundary condition treatments in Chapter 3, an efficient and accurate heat transfer evaluation technique is proposed and numerical validated in Chapter 4. Chapter 5 presents an MRT LB model for the axisymmetric CDE. Next in Chapter 6, a radiationconduction coupled numerical model is developed for simulating the energy transport in a hightemperature solar ther mochemical reactor based on a MonteCarlo ray tracing radiation model and the TLBE conduction model. Chapter 7 summarizes the dissertation and discusses the related future works. PAGE 28 * This work has been published in ASME Journal of Heat Transfer (Li et al. 2012). 28 CHAPTER 2 HEAT TRANSFER BETWEEN COLLIDING SURFACES AND PARTICLES* 2.1 Introduction Transient heat transfer across the interface between two contacting solids initially at different temperatures has been investigated in many previous works ( Heasley 1965, Barber 1989, Ammar et al. 1992, Li and Mason 2000, Kostoglou and Konstandopoul os 2002). Although some approximate solutions for the transient heat flow (Heasley 1965) and an asymptotic solution for the short time transient heat conduction (Barber 1989) have been presented, no analytical expression was provided in their work for predicting the amount of heat exchange between contacting bodies. Sun and Chen (1988) carried out a computational study using a finite difference method to solve the unsteady conduction equations in the region near the contact area of the colliding surfaces. B ecause both the contact area and the collision time are typically very small, the heat transfer media were assumed to be semi infinite which greatly simplified the formulation of the problem and the subsequent computation. They gave a closedform solution for the heat transfer in the limit of zero impact Fourier number ( Fo based on the maximum contact radius and the contact duration). Typical numerical results from the finite difference solution were presented graphically over a range of Fourier numbers, t hermal diffusivity ratios, and conductivity ratios. Zhou et al. (2008) employed the finite element method to solve the heat conduction problem between colliding particles. The temperature distributions on the entire body of the colliding spheres were computed and the total amount of heat transfer during the impact was determined. A closedform approximate expression for the heat transfer was developed based on their computed results for a special case when the colliding materials have the same thermal PAGE 29 29 cond uctivity. Close examination also indicates that the energy transfer computed by Zhou et al. (2008) is noticeably smaller than that by Sun and Chen (1988). The main objective in the present research is to develop an accurate closedform expression to predic t the amount of heat transfer between two colliding particles with arbitrary impact Fourier number, thermal diffusivity ratio, and thermal conductivity ratio. Hertzian contact theory of elasticity (Timoshenko and Goodier 1970, Johnson 1987) is used to char acterize the particle impact and associated deformations The semi infinite media assumption ( Sun and Chen 1988) is also employed and the governing conduction equation is then transformed to a moving coordinate so that the computational domain remains fixed. Grid stretching in both the radial and axial directions is applied to resolve the high temperature gradient and singular regions. A self similar solution of the temperature field during the initial contact period is obtained to provide an accurate initi al condition for the simulation of the entire collision period and to gain an insight into the asymptotic state of the heat transfer. In addition, the singular behavior of the discontinuity at the edge of the contact area is analyzed based on a local analy tical solution, which also offers an insight into the evolution of the temperature field in the whole spatial domain and the heat flux during the collision period. The asymptotic solutions for the thermal field and the transient heat flux over the contact area at small time and largetime are obtained in the zeroFourier number limit. The behavior of the transient heat flux and heat flow rate in the contact area is fully simulated and analyzed. Finally, a closedform expression for predicting the heat exchange between impacting particles of various thermal properties is provided based on various asymptotic limits and computational results over a wide range of Fourier PAGE 30 30 numbers, thermal diffusivity ratios, and thermal conductivity ratios of the colliding materi als. 2.2 Formulation of the Collisional Heat Transfer Problem 2.2.1 Elastic Collision and Deformation When two spherical particles of masses m1 and m2, radii R1 and R2, elastic moduli E1 and E2, and Poisson ratios 1 and 2, impact with a normal relative velocity v12, the contact area A is circular and its rate of change during the compression process can be expressed as below according to the Hertzian elastic theory (Timoshenko and Goodier 1970, Johnson 1987) 1/2 2 5/2 12 1212 124 () 5 E dA dt m (2 1) where 12 12 12RR R RR (2 2) 12 22 11224/3 (1)/(1)/ E EE (2 3) and 12 12 12mm m mm (2 4) If the second surface is flat, R2 should be set to infinity; and if the second surface is the interior wall of a pipe, R2 should take the negative value of the pipe radius. The decompression is an inverse process of the compression. The maximum contact area Ac is obtained by setting the derivative in Eq. (21) to zero. PAGE 31 31 2/522 4/51212121254ccmRA E (2.5) Equation (21) can also be integrated to yield a dimensionless relationship between the transient contact area and the collision time as *5/21/2 0(1)Adx x (2 6) where 2/5 1/5 12 1212 124 5 E m (2 7) and 22//ccAAARR (2 8) In the compression process Eq. (2 6) is numerically integrated from A* = 0 to A* = 1 to yield a dimensionless time of c = 1.4716376. The total contact duration tc is then found from Eq. (2 7) as 2/51/512121212524ccmtE (2 9) Defining t*= t / tc= /(2 c) (2 10) the relationship between A* and t* can be numer i cally obtained; the result is shown in Fig. 2 1. To carry ou t asymptotic analyses for the heat flux in the contact area, the following are noted for t* * 7/2*7/2 6*611(2)(2)()(2)()...7 112cc cA (2 11a) PAGE 32 32 ** *5/2*5 12 **1()()... tdA btbt Adt (2 11b) with 5/2 51255(2),(2)14 784ccb (2 11c) Equation (211a) also implies that **At and hence *1/2() Rt as t* (2 12) 2.2.2 Governing Equations for Heat Transfer The transient heat conduction equation for the two colliding particles in the axisymmetric cylindrical co ordinates (Fig. 22) is (Sun and Chen 1988, Zhou et al. 2008) 2211() (1,2)iiiiTTTri (2 13) The semi infinite media assumption was used by Sun and Chen (1988). Its validity is closely examined in this study and the details are given in Appendix A It is shown that this assumption is valid over a wide range of conditions and thus is used in the present work. Hence for a contact area with instantaneous contact radius R ( t ), the boundary conditions are summarized as follows 12 121 2 12, () 0,0, () TT TTkkrRt zz TT rRt zz at z = 0 (2 14) T1 = T01 as z T2 = T02 as z (2 15) T1 = T01 and T2 = T02 as r (2 16) 120 TT rr at r = 0 (2 17) where k1 and k2 are the thermal conductivities of the two particles, T01 and T02 are the surface temperatures of t he particles before contact ( T01 > T02). PAGE 33 33 The initial condition is assumed to be T1 = T01 and T2 = T02 at t = 0 (2 18) The above formulation is the same as described in Sun and Chens work (1988). It is important to note that the discontinuity of the boundary conditions in Eq. (214) occurs at r = R ( t ) which changes with time. Noting that the contact area varies with time during the entire collision period, a coordinate transformation is applied so that the computation can be more conveniently carried out i n a well defined, fixed coordinate system. To this end, the following dimensionless variables are introduced 02 0102 1, () TT zr RtTT (2 19) The governing equations become 2 ** 1111 1 2*1 () 22 tR Fo t (2 20a) 2***2222 2* 2*1[ ()]22 tRt (2 20b) where ** **2 tdRtdA R RdtAdt (2 21) = 2/ 1 (2 22) and the characteristic impact Fourier number is 1 2 c c Fo R (2 23a) Using Eqs. (2 5) and (29), the Fourier number can also be expressed in terms of the equivalent radius and the impacting velocity as PAGE 34 34 1 1 1212 12122 2.9432752c Fo Rv Rv (2 23b) The boundary conditions given by Eq. (214) now become 1222121, 1 at = 0 (2 24a) 120,0 > 1 at = 0 (2 24b) It is important to note that the discontinuity in the boundary conditions in E q. (2 14) now occurs at = 1 which does not vary with time. The conditions given by Eqs. (215 2 17) can be similarly modified to: 11 and 20 as (2 25) 11 and 20 as (2 26) 120 at = 0 (2 27) In order to maintain sufficient resolution near the contact area (especially near the singularity at = 1) and to obtain a computational domain large enough to avoid numerical interference of the far field boundaries, grid stretching in both directions is e m ployed. Equation (220) is solved using the alternative direction iteration (ADI) scheme. A second order upwind scheme is used for the first order derivative terms in Eq. (220). This avoids the grid level oscillation which normally occurs in the presence of large advection when central difference is used. This also avoids heavy numerical damping when a first order upwind is used. For finite Fo near t* = 1, it is found that the integral of the heat flux over the contact area increases rapidly due to the detachment of the PAGE 35 35 particles. This necessitates a gradual reduction of the time step to ensure accuracy of the solution. 2.2.3 Self S imilar Solution While the initial condition for the temperature distribution given by Eq. (218) is theoretically correct at t = 0, it results in persistent oscillation of the temperature during the initial period when Eq. (218) is applied directly. Thus a self similar solution for 1 t is sought. As t* 0, t* i t* 0, from Eqs. (2 11) and (221) we also have R* 1 and t*/ A* 1/(2 c), thus the transient conduction problem has a self similar steady solution in the transformed coordinate ( ) in the small time limit, then Eq. (220) becomes 2211[+()] (1,2)222iiii ic (2 28) The boundary conditions for the self similar solution are the same as given by Eqs. (2 24 2 27). The above equation can be further simplified when 1 For example, only half of the domain needs to be solved and 12(,0)(,0)0.5 can be enforced directly for 1 on the contact area. The converged solution obtained through an iterative procedure is used as an initial condition for solving Eq. (2 20) on the same grids. Of particular interest is the limiting case of Fo = 0 with 1 In this special case, the governing equation (228) for 1 becomes 2 111 222 The solution for 1 at = 0 is the standard error function: 1(0,)(1erf(/2))/2 The solution for 1 at = 1 is, however, singular. A transformation using Y = /1 gives PAGE 36 36 2 2(1)(1/2) 0 22 Y At = 1, this equ ation becomes 2204YY ; the solution is another error function of 1[1erf()]222Y Thus near = 1 the leading order behavior for the temperature is 1 ~ 11 [1erf()] 2 221 The derivative on the interface for Fo = 0 near = 1 is 1 0(~1) ~ 11 221 (2 29) 2.2.4 Heat Transfer Denoting q ( t r ) as the instantaneous heat flux at the interface and Q ( t ) as the instantaneous heat flow rate across the cont act area, we have 101021110 01()(,) zkTTTqtrkz (2 30a) () 1 10102 1 0 0 *0 12() ()2(,) Rt c ckTTA A Qt t (2 30b) From Eq. (2 30) we can further define the dimensionless heat flux **(,) qt the dimensionless heat flux integral **()wqt and the dimensionless heat flow rate **()Qt as **10(,)qt (2 31a) 1**100()wqt (2 31b) *** ***()()wAQtqtt (2 31c) The total energy transfer during the impact is PAGE 37 37 0()cteQtdt (2 32) 2.3 Results and Discussion First, a grid refinement study is carried out using the self similar solutions as a surrogate. The results for the dimensionless heat flux integral **()wqt defined in Eq. (231b) at t* ent mesh sizes are compared in Table 21. The iterative solutions are considered to be converged when the maximum difference in the temperature field between two continuous iterations is less than 1010. Typically, 320x128 inter vals in the and directi ons are chosen to obtain the temperature field, respectively, on each side of the co n tact. For the purpose of extracting the parameters associated with the singular behavior near = 1 throughout this paper, grids much finer than 320x128 have been used to ensure the reliability of the numerical sol u tions. For finite Fo and finite t*, a nominal =1/800 is used; at t* is used in order to capture the rapid rise in the heat flux near the detachment. For Fo = 0, = 1/800, 1/1600, and 1/3200 were used in order to obtain an accurate value of the corre ction factor C0 (see Eq. (2 61)). The infinity boundary in the direction is placed at = 15. In the direction, the boundary location varies depending on the Fourier number. For Fo = 0, = 1 is used and for Fo > 1, = 35 is used. 2.3.1 Self S imilar Solution at Small Time The self similar solution of the temperature field serves as an excellent initial condition to obtain accurate and smoothly varying temperature fields. M ore importantly, it is found that the dimensionless heat flow rate **()Qt obtained from the self similar solution remains very close to the time dependent solution for t* < 0.2 (see Fig. 210 later). Thus a great deal of insight can be gained by examining the much simpler self  PAGE 38 38 similar solution first. Since in most of the practical situations heat transfer occurs between particles of the same material, only the results corresponding to 1 will be presented here f or brevity. Figure 23 ( a d ) shows the temperature contours for Fo ranging from 0.02 to 10 based on the self similar solutions Fig 2 3( b ) also includes a schematic of the grid distribution near the discont i nuity. For Fo = 0.02 and Fo = 0.1, the contours are very similar, indicating the existence of an asymptotic limiting solution as Fo 0 which can be representative of the solution for Fo 0.1 to the leading order. Fig. 2 3 (c d) clearly shows that as Fo increases, the heat conduction in the r direction becomes more and more significant as part of the t o tal heat conduction. From the self similar solution, the dimensionless heat flux **(,) qt in the contact area from = 0 to = 1 can be obtained and the results at different Fourier nu mbers are presented in Fig. 24. As Fo increases, there is a strong increase in **(,) qt throughout the entire contact area. It is also clear that there exists a singularity for **(,) qt at = 1, as shown analytically by Eq. (2 29) for the case of Fo = 0. This singularity has a significant effect on the computation of the total heat flow rate across the contact area. The nature of this singularity for Fo > 0 is further investigated next. 2.3.2 Behavior of the Singulari ty at ( r z) = ( R ( t ), 0) Because of the discontinuity in the boundary conditions at ( ) = (1, 0), the solutions to are expected to be singular near ( ) = (1, 0) for Eq. (2 20) as well as for Eq. (2 28). In order to gain an insight, 1 is considered for simplicity and thus only the temperature on the surface of particle 1 is simulated. Thus, the boundary condition at = 0 becomes PAGE 39 39 (1)0.5 and (1)0 at = 0 (2 33) To investigate the behavior of the singularity at ( ) = (1, 0), the following variables are introduced 2211/2c1, ,tan(/)(/2) Fo (2 34) With the boundary conditions in Eq. (233), the local solution to is expected to be 1~0.5(), ~()nn (2 35) One can also deduce that the dominant terms in Eq. (228), in the limit of become 22 2221 ~0 (2 36) The boundary involved in Eq. (233) for the local singular solution dictates that n = 1/2 in Eq. (235). Thus the solution to the Laplaces Eq. (236) su b jected to Eq. (2 33) is = 0.5 + w 1/2cos( /2) (2 37) Hence the leading order solution to near = 0 is 1/2221/4 1(/2)(1) 10.5[ ]cos[tan ]/2 2 1ccFo (2 38) where w is a constant that depends on Fo and needs to be determined by matching the local solution to the finite difference solution in the o uter region for different values of Fo It is observed from the local solution, Eq. (238), that PAGE 40 40 1/4(/2)(,0) for 121cFo (2 39) thus there is an integrable squareroot type of singularity in the heat flux, **(,) qt near = 1 (or r = R ( t )) in the contact region for all values of Fo as shown in Fig. 24. Although the foregoing analysis is presented for particles of the same material ( 1 ), and for small time, the nature of the singularity is the same for part icles of different materials and for the entire duration of collision as well. Figure 25 (a) shows the variations of ( 0.5)/( 1/2) given by Eq. (238) from ) is another coordinate system related to ( ) with 22(1) and = tan1( /( 1)) at Fo = 0.02, 1.0, and 10.0 together with the finite difference results obtained from the self similar solution. The values of w in Eq. (2 38) are determined by fitting the analytical solution to the numerical solution near the singularity. Excellent agreement between the analytical and numerical solutions c an be observed for the variation in the azimuthal direction. Although not shown here, similar agreement is observed for the variation of in the direction at = 0 for a wide range of Fo values. This ensures the validity of the local solution given by Eq. (238). For very small Fo there is a rapid decrease of as consistent with the temperat ure contour in Fig. 23( a ) near = 1 and = 0. From > 1 to < 1, for a given local radius there is a large drop in in the tangential direction because of the discontinuity in the boundary conditions in Eq. (233). Comparing Eq. (2 29) with Eq. (2 39) near = 1, it can be deduced that w ~ 1/4(2/)/2c = 1/40.52254/Fo as Fo 2 5 (b) shows the variation of the w value with Fo in the local singularity from the self similar solution. As Fo PAGE 41 41 Fig. 2 5 (b) agrees very well with the prediction w ~ 1/40.52254/Fo As Fo becomes larger, the stronger diffusion together with the advection in the direction causes a much stronger heat flow in the direction toward the contact point. Hence w increases with increasing Fo and for Fo > 5, the numerical result can be fit with w ~ 1/20.620.05772Fo very well. For larger value of t* (0 < t* < 1), self similarity no longer exists. The dominant terms in the governing equation (220) for part icles with the same material near the discontinuity point becomes 2* 2*1 ()0 Fo (2 40) The transient local solution can be obtained through similar procedure as 2 **1/221/4 1**(1) 1(/)~0.5[ ]cos[tan ]/21 FotA (2 41) It should be noted that the constant w in Eq. (2 41) now varies not only with the Fourier number but also with time during the collision and should be determined by matching the local solution to the timevarying outer region finite difference solution. Figure 26 shows the temperature variations along the axis ( = 0) near = 1 when Fo = 1, at t* = 0.40, 0.60, 0.70, 0.80, 0.90, and 0.95 from the finite difference and local analytical solutions. The values of w used in the local solution are: 0.645, 0.545, 0.480, 0.430, 0.370 and 0.330 respectively, whic h are determined again from the numerical solutions for the temperature field near the singularity point. The transient behavior of the temperature field is largely the same for t* < 0.25 as given by the self similar solution. In the region close to the di scontinuity ( = 1), the computed results from finite difference solution agree well with the analytical solutions given by Eq. (2 PAGE 42 42 41). The decrease in w as t* increases suggests that while the nature of the singularity remains the same, the strength of the singularity continues to weaken during the decompression process. It is noted that the increasing diffusion in the direction (see the second terms on the right hand side of Eq. (220)) tends to smear the discontinuity and decrease effectively the strength of the si ngularity. The increasing larger advection in the positive direction given by the second terms on the left hand side of Eq. (220) for t* > 0.5 represents the sweeping from the cold region ( < 1) with = 0.5 to the warmer region ( > 1) and reducing ef fectively the temperature values after the contact point. This is in general true for all other values of Fo The insight gained from the local solution for the thermal field is essential to the understanding on the behavior of the global thermal field. 2. 3.3 Asymptotic Limit at Fo 0 In the present study, the numerical result at very small Fo indicates that when 1 *~1/(2)wq is obtained from the self similar solution (see Fig. 2 4). Thus the asymptotic limit (as Fo 0) for the heat transfer is pursued by solving Eq. (220) while setting Fo = 0. Because of the loss of the diffusion term in the direction in Eq. (220) when Fo = 0, the solution must be separated into two stages depending on the sign of R*. When R* > 0, inspection of Eqs. (2 20) and (2 24) quickly yields 121,0 for ( 2 42) 2***21 (1,2)22iiiiitR it for 1. (2 43) PAGE 43 43 It is clear that in the absence of the 22/i term, Eq. (2 43) is parabolic in the direction. The boundary condition by Eq. (227) must be abandoned for R* > 0, thus the only boundary condition in the direction should be 11 and 20 at = 1 (2 44) At each time step, the temperature *(,,)i for < 1 can be obtained numerically by marching from = 1 to = 0. On the other hand, when R* < 0 ( t* > 1/2), Eq. (2 27) should be applied as the only boundary condition in the direction and the solution is obtained by marc hing from = 0. While the foregoing procedure leads to correct numerical solution for the Fo = 0 case, further insight in the asymptotic limits can be gained by considering a perturbation solution for the temperature distribution. At small time t* 1 the temperature can be sought in the form of *2.5 *5,0 ,1 ,2(,,)(,)(,)()(,)()...iiii (2 45) The powers of t* are determined based on the expansion of A*( t*) given by (211b). Substitution of Eq. (245) into Eq. (243) leads to a set of steady state parabolic equations for 1,0, 1,1 and 1,2 which can be solved numerically using a finite difference method. These solutions in turn give the heat flux integral in the contact zone as 111,0 1,1 1,2** *2.5 *5 *1000000() ()(), 1wqt (2 46) The result is qw *( t*) ~ 0.2820948 0.430065( t*)2.5 0.236343( t*)5 *1t (2 47) PAGE 44 44 Similarly, as t* 1 at the largetime stage, one can deduce from Eq. (211) that *****tdARAdt *5/2 *5 12 *[1(1)(1)...] ~ 1 tbtbt t the solution for the largetime behavior of the temperature near the detachment point can thus be sought in the form of ***,0 1,1 2,2(,,)(,)()(,)()(,)...,iiii t* 1 (2 48) where the leading term i, 0 must be of O (1) due to the finite boundary conditions while 1 and 2 are the gauge functions for the next two higher order terms wit h 2( t*) 1( t*) 1. Substitution of the above expansion into Eq. (243) gives ,0 ,1 ,2 *'' *3/2 1.12.2 1 1 2 *1 [1(1)]( ...)[1(1)...][ ...] 12iii ii t t 222 ,0 ,1 ,2 ,0 ,1 ,2 12 12 222 1[ ...][++...] (1,2) 2iii iii i (2 49) the largest term in the above is of O (1/(1 t*)). The balance of the equation to O (1/(1 t*)) requires that ,00 (1,2)ii (2 50) Equation (250) simply states that i, 0 is a function of only. Since *0 R for *1/2t ,0()i can be de termined at = 0. Setting = 0 in Eq. (243) for ,0()i gives 2,0 ,021 (1,2)2iiii (2 51) The solution to the above can be expressed in terms of the error functions. For particles of the same material ( 1 ) we will have 1,0(,)(1erf(/2))/2 (2 52) PAGE 45 45 Th e second largest term in Eq. (249) is of O (1) and the balance of the equation to O (1) gives 2,1 ,0 ,0'11.1*21 (1,2)122iiiii t (2 53) which requires 1 = 1 t* (the continuation of the same procedure would yield 2 = (1 t*)2). Then Eq. (2 53) becomes 2,1 ,0 ,0.121 (1,2)22iiiii (2 54) Based on Eq. (2 51) the above can be further simplified to ,1 ,10 (1,2) 2i i (2 55) so that 2,1() (1,2)ii (2 56) With the numerical solutions to Eq. (243) and the expansion (248), fi( ) can be numerically estimated from *2 ,0()~[(,,)(,)]/[(1)]iiif (2 57) Figure 27 shows fi( ) for 1 at t* = 0.99 with = 0.2, 0.4, 0.6 a nd 0.8. The numerical results collapse very well into a single curve. Furthermore, the data can be fit almost exactly using 21()exp()44if (2 58) so that the second term in Eq. (248) can be approximated as 22 1,1(,)exp() 4 4 for 1 (2 59) PAGE 46 46 Thus the large time behavior of the heat flux integral for 1 can be approximated as 1 ** 1 0 05 (1) ~ 16w t qt (2 60) Figure 28 compares the computationally obtained **()wqt wit h the perturbation solutions given by Eqs. (247) and (260) for particles of the same material at Fo = 0. Excellent agreement in the limits can be observed. To compute the heat flow rate, the heat flux integral is multiplied by the timevarying term **/At as shown in Eq. (231c). Fig. 2 9 shows the variation of the multiplier of **/At and Fig. 2 10 shows the heat flow rate *Q ( t*) during the entire collision period for particles of the same material at Fo = 0. The rapid rise of *Q ( t*) during the initial contact is associated with the thin thermal diffusion layer on the order of ( t*)1/2 and the contact area on the order of t* as indicated by Eq. (212). For compari son, the result of *Q ( t*) obtained from the self similar solution is also included in Fig. 2 10. Integrating *Q ( t*) over the entire collisional period gives the energy transfer at Fo 0 1/20010201/2 1/2111 222()()()ccCTTAte (2 61) where the subscript 0 for e0 indicates that this is the energy transfer for the zeroFourier number case. For particles of the same material, which will be used to normalize the energy trans fer under general condit ions, extensive computation with high grid PAGE 47 47 resolution yields C0 = 0.87093. This is in exact agreement with the earlier result of C0 = 0.87 given by Sun and Chen (1988) based on their asymptotic analysis for Fo = 0. 2.3.4 Numerical Results for Finite Fo To understand the effect of the characteristic Fourier number on the heat conduction, the dimensionless heat flux integral **()wqt in the contact area is examined using the full numerical solution to Eq. (220). Figure 211 shows the variations of wq with time at different Fo values. When Fo to that of Fo = 0. However, for Fo at the late stage, and its peak v alue occurs at the end of the collision when the diffusion in the r direction becomes dominant. Thus the heat flux **(,) qt throughout the contact area as well as the temperature variations along the axis ( = 0) at various times are examined. Fig s 2 12 and 21 3 show **(,) qt and *(,1) respectively at Fo = 1. As shown in Fig. 21 2, as *t increases the overall heat flux **(,) qt increases. Towards the end of the collision, **(,) qt increase more rapidly in the central region of the contact. In the meantime, the strength of the singularity at = 1 becomes weaker as the distribution along becomes less steep as *t 2 13 for (130) at different times, the strength of the singularity continues to decrease as () spre ads out further which is caused by the increasing diffusion in the direction, 1 *1 () t Fo A and the increasing advection from the cold contact region, *2i R as discussed at the end of Section 2.3.2. In fact, Fig. 21 3 is a large scale view of Fig. 26; together they show that the large scale temperature distribution in Fig. 213 is actually a consequence of the PAGE 48 48 local singular temperature distribution near = 1. As *t **/tA so that the diffusion of the heat flow in the presence of the singularity in the direction becomes extremely effective; this causes a rapid rise of **(,)qt as *t 1 in the c entral portion of the contact area as shown in Fig. 212. In order to develop a general expression for the heat transfer during the collision, the thermal fields for a large range of Fourier numbers are computed and the heat transferred, e is normalized by e0 in Eq. (2 61) as C e / e0 = func(, )Fo (2 62) Asymptotic analysis of the results for ''1 indicates that there are two limits for the correction factor C as a function of Fo : (,1,1)CFo = 1 + 1. 07646 Fo when Fo 1 (2 63a) (,1,1)CFo = 0.605039 + 1.08748 Fo1/2, when Fo 1 (2 63b) Combin in g the above, the following formula for C is obtained (,1,1)()CFo 0.332298e2 Fo) Fo ]1/2 (2 64) Fi gure 214 compares the computed e/e0 = g ( Fo ) versus the approximation given by Eq. (2 64). Excellent agreement is observed. The computational results obtained by Sun and Chen (1988) and Zhou et al. (2008) are also included in Fig. 214. The result of Sun and Chen (1988) is very close to the present work for Fo 10, which indicates that the finite difference solutions by Sun and Chen (1988) are highly accurate despite the low resolution used in their work. At high Fourier numbers, the heat exchange based on a finite element method by Zhou et al (2008) is much smaller than the present result and that in Sun and Chens (1988) work. PAGE 49 49 2.3.5 Particles of Different Materials When particles of different materials collide, besides the Fourier number, the heat transf er will be affected by two additional dimensionless parameters the thermal diffusivity ratio and the thermal conductivity ratio k as observed in Eqs. (220b) and (2 24a), respectively. The asymptotic lim it for e0 given by Eq. (261) will still be used to normalize e under general conditions. It is first noted that when two particles have the same value of thermal diffusivity ( = 1), the governing equations for the two particles are identical: 2 ** 2*1 ()(1,2) 22iiii i tR Fo t (2 65) when = 1, it is convenient to introduce a conductivity weighted average temperature ( t*, ) that satisfies Eq. (2 65) ( t*, ) = ( k11+ k22)/(k1+ k2) ( 2 6 6 ) The relevant boundary conditions at = 0 and 00 (2 67) ( t*, ) = k1/( k1+ k2) (2 68) At t* = 0, the initial condition is (0, ) = k1/( k1+ k2) (2 69) Thus by inspection, the solution to must be ( t*, ) = k1/( k1+ k2) (2 70) That is k11( t*, ) + k22( t*, ) = k1 ( 2 7 1 ) Combing Eq. (271) with 1( t*, 0) = 2( t*, 0) in the contact area, it is revealed that PAGE 50 50 1( t*, 0) = 2( t*, 0) = k1/( k1+ k2) (2 72) Equation (272) shows that at = 1 the temperature in the contact area is constant when the conductivity ratio 21/kkk is fixed; a similar phenomenon was also pre sented by Heasley (1965), where they derived that the steady temperature at the edge of the contact spot is constant. Here we have obtained in Eq. (272) a more general conclusion that the temperature throughout the contact area at anytime is constant, whi ch was also numerically verified for all Fo and k cases in our study. Furthermore, the temperature difference between = 0 and for particle 1 is k2/( k1+ k2), which drives the wall heat flux at = 0. Thus the heat flux integral **()wqt and the heat transfer follows **()wqt k1k2/( k1+ k2) (2 73) e0 k1k2/( k1+ k2), and e k1k2/( k1+ k2) (2 74) Hence C = e / e0 is independent of k when = 1 and Eq. (264) is also valid for arbitrary values of k when = 1. This fact is directly verified in the present computational study for the entire range of Fourier numbers (from 0 to 50) and this observation is critical to deduce the form postulated in Eq. (277) la ter on for > 1. It is noticed that this fact was not represented in the closedform formula provided by Zhou et al. (2008), which is only applicable to the special case of k = 1. For 1, it is convenient to introduce a relative correction factor based on the results from 1 as (,,) ,,/, 1, 1,,/ (2 75) PAGE 51 51 where g ( Fo ) is given by Eq. (2 64). An accurate approximation for allows one to obtain C = e / e0 from Eq. (2 75). For > 1, in the limit of very large Fo the computational results show that has the following asymptotic behavior (,,)(,,)/(,1,1)~/(1)/(1) (2 76) It should be noticed that the Fo 76) is bounded by the upper limit as restricted by the semi infinite media assumption (see Appendix A ). In our simulations, the asymptotic behavior of Eq. (2 76) is obtained at Fo = 50. For finite Fo the relative correction factor can be approximated as 2(1)4(,,)11(1) kk ( > 1) (2 77) with 1/20.95 0.950.500.12 Fo FoFo (2 78) 210(log) 1/31/20.0777 [1] 0.7666 FoFo (2.79) Since switching the definition of the subscripts of the part icles does not affect the amount of heat transfer, there exists a symmetric relation for < 1 (,,) (, 1/, 1/)CFo (2 80) Thus the correction factor for < 1 can be simply written as (, 1,) (, 1/, 1/)() CFo (2 81) where f is given by Eq. (277). Figure 2 15 shows the correction factor C from the finite difference method and the above approximate expressions for k and ranging from 0.1 to 10. 0. Good PAGE 52 52 agreement suggests that the approximate formula developed in this study based on numerical computations and asymptotic analyses is accurate. It can be used as a practical tool to predict the heat transfer between two colliding particles or between particles and surfaces. 2.4 Summary and Conclusions Heat conduction between two colliding particles of various materials is computationally and analytically investigated over a large range of Fourier numbers with high quality resolution. The singular natur e of the thermal field and heat flux at the edge of the contact area is elucidated through analytical solutions. A self similar solution for the thermal field during the initial period of contact is developed. It serves as an accurate initial condition for the heat transfer simulation during the entire collisional period. A twodimensional asymptotic analysis is presented for the thermal field and heat flux in the contact area in the small Fourier number limit to facilitate the understanding of collisional heat transfer. A closedform expression is developed for the heat transfer prediction as a function of the Fourier number, the thermal diffusivity ratio and the thermal conductivity ratio of the impacting particles. PAGE 53 53 Table 21. Heat flux integral **()wqt based on the self similar solution at small time ( *1 t ) 320 128 640 256 1280 512 Fo = 0.1 0.300430 0.300085 0.299867 Fo = 1.0 0.396332 0.396192 0.395993 Fo = 10.0 0.772609 0.775793 0.777692 PAGE 54 54 Figure 21. Variation of the contact area during collision. Figure 22. Computational domain and boundary conditions at z = 0 PAGE 55 55 ( a ) Fo =0.02 ( b ) Fo =0.1 & schematic of grids ( c) Fo =1.0 ( d ) Fo =10.0 Figure 23. Temperature contours on the surface of particle 1 based on the self similar solutions at small time ( *1 t ) at ( a ) F o =0.02, ( b ) F o =0.1, ( c) F o =1.0 and ( d ) F o =10.0. PAGE 56 56 Figure 24 Dimensionless heat flux across the contact area at small time ( *1t ) based on self similar solutions PAGE 57 57 00.511.522.533.5 00.511.522.53 00.20.40.60.811.21.4 Fo=0.02 Fo=10.0 Fo=1.0 symbols: finite differencelines: analytical (a) 0.60.70.80.911.11.21.31.41.5 0.010.1110100 Numerical 0.5225Fo 0.62+0.05772Fo Fo 1/4 1/2 w (b) Figure 25 (a) Local solution of ( 0.5)/( 1/2) around the si n gularity point at ( ) = (1, 0) ( b ) The value of w in the local solution of around the singularity point at ( ) = (1, 0). PAGE 58 58 Figure 26 Local transient temperature profiles at ( > 1, = 0) when Fo = 1 .0. Figure 2 7 Compasion of numerical result s and analytical fit for fi( ) in Eq. (2 58) PAGE 59 59 Figure 28 Variation of the dimensionless heat flux integral **()wqt over the contact area at Fo = 0 Figure 29 Variation of the multiplier **/ At in Eqs. (2 30b) and (231c) PAGE 60 60 Figure 210. Variations of the dimensionless heat flow rate **()Qt at Fo = 0; the self similar solution is **()Qt = 0.2820948 **/ At Figure 211. Comparison of the dimensionless heat flux integral **()wqt at different Fourier numbers PAGE 61 61 Figure 212. Heat flux **(,) qt in the contact area at various times for Fo = 1 .0. Figure 213. Temperature variations along direction ( = 0) at various times for Fo = 1 .0. PAGE 62 62 Figure 214. Heat transfer as a function of the impact Fourier number (particles of the same material) PAGE 63 63 ( a ) ( b ) ( c) Figure 215. Correction factor comparison between computation and the approximate e xpression at ( a ) k2/ k1=1.0, ( b ) k2/ k1=0.1 and ( c) k2/ k1=10.0. PAGE 64 * This work has been published in Journal of Computational Physics (Li et al. 2013a) 64 CHAPTER 3 BOUNDARY CONDITIONS FOR THERMAL LATTICE BOLTZMANN EQUATION METHOD 3.1 Literature Review of Thermal Lattice Boltzmann Equation (TLBE) Models The LBE method simulates the particle interactions on a discrete phase space with a specific set of discrete velocities, and the particle dynamics is represented by the particle distribution functions (also called particle density distribution function or particl e velocity distribution function). The velocity distribution functions of the particles are governed by the lattice Boltzmann equation, from which the macroscopic Navier Stokes equations can be recovered. Particle collision in the LBE method is approximated by a relaxation process toward the equilibrium state of the particle distributions T here are basically two collision models for the LBE methods: the first and most commonly used single relaxationtime (SRT), or the Bhatnagar Gross Krook (BGK) model (Qia n et al. 1992, Chen et al. 1992) based on the work of Bhatnagar et al. ( 1954); and the multiplerelaxationtime (MRT), also known as the moment model (dHumires 1992, Lallemand and Luo 2000, dHumires et al. 2002). The MRT LBE model overcomes some obvious defects of the BGK model such as fixed Prandtl number and fixed ratio between the kinematic and bulk viscosities (dHumires et al. 2002) and shows better numerical stability ( Lallemand and Luo 2000) and reduced numerical dispersion (Yu et al. 2003). Bec ause of the advantages of the LBE method, there have been continuous interest and effort in developing thermal LBE (TLBE) models for thermal transport problems (Alexander et al. 1993, Chen et al. 1994, Shan 1997, He et al. 1998, Lallemand and Luo 2003, Guo et al. 2007, Shi and Guo 2009 Mezrhad et al. 2010, Yoshida and Nagaoka 2010). In general, the TLBE models fall into four categories: the PAGE 65 65 passive scalar approach, the multispeed approach, the hybrid approach, and the doublepopulation distribution functi on (DDF) approach. In the passivescalar based TLBE model, the macroscopic temperature is considered as a passive scalar that is advected by the velocity field and satisfies the convectiondiffusion equation (Shan 1997). The main shortcoming of this model is that the temperature field is driven by the flow velocity, while the temperature does not affect the flow field. The multispeed approach is a direct extension of the athermal LBE model so that only the velocity distribution function is used. In order to obtain the macroscopic temperature, additional discrete speeds and higher order velocity expansion for the Maxwell Boltzmann equilibrium distribution are required in this model (Chen et al. 1994). In the hybrid approach, the mass and momentum conservation equations are solved with the athermal LBE model, whereas the advectiondiffusion equation for the temperature is solved separately by using the finitedifference method or other conventional numerical techniques ( Lallemand and Luo 2003). The DDF approach utilizes two different distribution functions, one for the velocity field and the other one for the temperature or energy distribution. With this approach, the interaction between the flow fields (density and velocity) and the temperature field can be properly accounted for. Based on the methods of defining the temperature distribution functions, the DDF approach can be further classified into two groups. The first derives the temperature distribution function from the continuous Boltzmann equation as that in the athermal LBE model. Thus these two sets of distribution functions are explicitly coupled in the thermal LBE. For example, He et al. (1998) used an internal energy based distribution function and Guo et al. (2007) used a total energy based distribut ion function to PAGE 66 66 simulate the thermal flows. In the other group for the DDF approach, however, the temperature distribution function is defined with the objective of recovering the macroscopic convectiondiffusion equation (CDE), in which the other terms of the energy equation such as pressure work and heat dissipation are conveniently treated as source terms. With this approach the macroscopic velocities, rather than the velocity distribution functions, are included in the evolution equation of the temperat ure distribution functions. Thus no explicit coupling between the velocity and temperaturedistribution functions is needed. To this end, Shi and Guo (2009) proposed a thermal LBE model for nonlinear CDE simulations, in which the equilibria of the temperature distribution functions are defined to recover the correct CDE. Recognizing the advantages of the MRT model over the BGK model, Mezrhab et al. (2010) and Yoshida and Nagaoka (2010) proposed their MRT TLBE models for the CDE. To recover the CDE, the equilibria of the thermal moments are specified in Mezrhab et al.s model while the thermal equilibrium distribution functions are specified in Yoshida and Nagaokas model. A detailed derivation of the macroscopic CDE from the MRT TLBE model can be found in Y oshida and Nagaokas work (2010). The present work also uses the MRT LBE model by Yoshida and Nagaoka. 3.2 Literature Review of Boundary Conditions in TLBE Simulations Accurate boundary condition treatment is essential for the overall integrity of the LBE simulations. For the fluid velocity field, the bounceback scheme has been intensively studied (Ziegler 1993, Noble et al. 1995, Ginzbourg and dHumires 1996, Chen et al. 1996, He et al. 1997, Zou and He 1997, Ginzburg and dHumires 2003) and extended to handle curved boundaries ( Filippova and Hnel 1998, Mei et al. 1999, Mei et al. 2000, Bouzidi et al. 2001, Guo et al. 2002). The standard bounceback PAGE 67 67 method is secondorder accurate for noslip wall boundaries when the distance from the wall to the nearest fluid node is half of the lattice unit spacing. For other hydrodynamic boundary conditions such as given density, pressure, slipping velocity or their derivatives, the bounceback scheme can be modified with the flow properties at on the boundary incl uded in the treatments ( Ginzbourg and dHumires 1996, Chen et al. 1996, He et al. 1997, Zou and He 1997, Filippova and Hnel 1998, Mei et al. 1999, Mei et al. 2000, Bouzidi et al. 2001, Guo et al. 2002, Ginzburg and dHumires 2003). Spatial interpolation or extrapolation has also been incorporated to improve the degraded accuracy of the bounceback scheme when the boundary is located in arbitrary position in the lattice cell (Chen et al. 1996, Filippova and Hnel 1998, Mei et al. 1999, Mei et al. 2000, Bo uzidi et al. 2001, Guo et al. 2002). With appropriate modifications, some of the athermal boundary condition treatments have been extended to TLBE simulations as well. He et al. (1998) extended the bounceback rule of the nonequilibrium distributions prop osed by Zou and He (1997) to impose thermal boundary conditions. This nonequilibrium extrapolation approach was also applied by Tang et al. (2005) Huang et al. (2006) and Jeong et al. (2010). Liu et al. (2010), however, used correction functions for the unknown energy distributions based on local known energy distributions, which is an extension of their hydrodynamic boundary condition treatment (Ho et al. 2009). It is noticed that those thermal boundary condition treatments in Refs. (Tang et al. 2005, Hu ang et al. 2006, Jeong et al. 2010, Liu et al. 2010) were directly proposed for Dirichlet boundary condition and the Neumann (flux) condition was transformed into Dirichlet condition by various finitedifference schemes. This transformation will become inc onvenient when PAGE 68 68 mixed thermal boundary conditions are encountered. As a consequence, the requirement of a direct boundary condition treatment for the Neumann condition becomes natural. DOrazio et al. (2003) used a counter slip approach to handle thermal boundary conditions, where a counter slip thermal energy density is introduced to represent the given temperature or heat flux condition on the boundary. Kuo and Chen (2009) proposed a nonequilibrium mirror reflection scheme to implement both the isothermal and adiabatic thermal boundary conditions. Ginzburg ( 2005) extended the multi reflection approach ( Ginzburg and dHumires 2003) to model both Dirichlet and Neumann boundary conditions. The recent 2D MRT TLBE model (Mezrhab et al. 2010) implemented the bounceback scheme to handle Dirichlet boundary conditions where they provided an explicit formula to determine the unknown populations based on the distributions at the neighboring fluid nodes and the temperature on the boundary. And a simple mirror ing treatment was used for the adiabatic condition. In the MRT TLBE model by Yoshida and Nagaoka (2010) for both 2D and 3D cases, boundary condition treatments for both the Dirichlet and the Neumann conditions were provided for the halfway case where the physical boundary is placed in the center of the lattice cell. Although thermal boundary condition treatment is an important issue, the discussions about the applicability and accuracy of these treatments in curvedboundary situations in the literature are rather limited. Huang et al. (2006) extended the nonequilibrium approach to deal with curved boundary condition, where the equilibrium part of the unknown population is evaluated according to the Dirichlet condition and linear extrapolation is used f or the nonequilibrium component; for Neumann condition, they PAGE 69 69 also transferred it into Dirichlet condition with a finite difference scheme and a simple projection rule was used to convert the normal boundary flux into the heat fluxes in the discrete veloci ty directions of the TLBE model. Secondorder accuracy of the temperature field was observed for Dirichlet condition on the curved boundary while no examples were provided regarding their heat flux conversion and the accuracy of the transferred Neumann condition treatment in their work. For curved boundaries with Neumann conditions, Yoshida and Nagaoka (2010) realized that the halfway thermal bounceback scheme leads to considerable disagreement between their simulation and analytical solutions. Thus they presented a modified treatment with the intersection area accounted to alleviate the discrepancy. Noticeable improvement was achieved in their numerical example but the accuracy was not investigated in detail. In the multi reflection approach for arbitr arily shaped surfaces, Ginzburg (2005) decomposed the equilibrium populations into the symmetric and anti symmetric components. The symmetric part was tuned to build secondorder and thirdorder accurate Dirichlet boundary conditions and the anti symmetr ic part was used for Neumann conditions. It was also pointed out that the Neumann condition with zero tangential heat flux yielded second order accuracy while only linear convergence was achieved when nonzero tangential heat flux exists (Ginzburg 2005). I t is also noted that the heat flux in one of the discrete velocity directions is required in (Ginzburg 2005) in order to remove the tangential flux constraint and to obtain linear convergence. In principle, however, only the derivative in the normal direct ion is specified in Neumann problems. Ginzburg (2005) did not provide a means to evaluate the heat fluxes in the discrete lattice velocity directions based on the normal derivative on the curved boundary. PAGE 70 70 3.3 Present Boundary Condition Treatments in TLBE M ethod The objectives of the present work are twofold. The first one is to develop second order accurate boundary condition treatments for both the Dirichlet and Neumann conditions on straight boundaries that can be placed at arbitrary locations in the lat tice. The second is to extend the present boundary condition treatments into curvedboundary situations and to investigate the effect of curvedgeometry on the accuracy of the present boundary condition treatments. The present boundary condition treatment s are independent of the specific TLBE model and the collision operator; thus for illustration purposes, we formulate the problem using Yoshida and Nagaokas MRT TLBE model ( Yoshida and Nagaoka 2010). For the Dirichlet condition, there is an adjustable par ameter in the treatment with second order accuracy. Three particular choices are discussed and their performances are examined. For the Neumann problem, the secondorder accurate boundary condition is unique. When the present boundary treatments are extend ed to handle inclined or curved boundary walls, the Dirichlet condition can be directly implemented. For Neumann problems on curved boundaries the given heat flux condition in the normal direction alone cannot be used to determine the fluxes in the discret e velocity directions in the TLBE model. Only for the special case where the tangential temperature gradient does not contribute to the fluxes in the discrete directions, i.e., the tangential gradient is known to be zero or the boundary is perpendicular to one of the discrete velocity directions, can a simple formula for the distribution function be obtained with second order accuracy based on the normal derivative. For general cases the boundary heat fluxes in the discrete velocity directions must be coupl ed with the unknown temperature value on the boundary which can be independently evaluated on PAGE 71 71 each lattice velocity direction. An expression for the heat flux determination in the discrete velocity directions is derived from the solution of the coupled equations. The accuracy of the proposed boundary condition treatments is examined using a series of numerical tests for which analytical solutions are available. The test cases include: i) 2D steady state channel flow with a constant and uniform velocity and a sinusoidal wall temperature or heat flux boundary condition; ii) 1D transient heat conduction in an inclined semi infinite solid with a constant normal heat flux at the end; iii) 2D steady state heat conduction inside a circle with a sinusoidal wall t emperature or heat flux boundary condition; iv) 2D transient heat conduction inside a circle with a transient wall temperature boundary condition; and v) 3D steady state circular pipe flow with a constant and uniform velocity and a sinusoidal wall temper ature or heat flux boundary condition. Because of the irregularity of the lattice fractions cut by the curved boundary, second order accuracy for the temperature field is realized for Dirichlet problems, but only first order accuracy can be obtained for Neumann problems on curved boundaries. 3.3.1 The Convection Diffusion Equation Following Yoshida and Nagaokas formulation ( Yoshida and Nagaoka 2010), the standard convectiondiffusion equation (CDE) can be written as ()()j ij jijvDG txxx (3 1) where is a scalar variable that can be either temperature in thermal simulations or mass concentration in mass transfer problems, vj is the velocity component in the xj direction, Dij is the diffusion coefficient, and G represen ts the source term. The Dirichlet boundary condition is d (3 2) PAGE 72 72 and f or the Neumann boundary condition involving normal derivatives it is in general written as iij jjn jnDnv x (3 3) where d and n are given functions of time and space on the boundary, ni and nj denote the i and j components of the unit normal vector pointing inward to the domain. 3.3.2 The D3Q7 and D2Q5 Thermal LBE models The discrete phase space in the TLBE model is defined by a regular lattice in D dimensions together with a set of discrete velocities { e = 0, 1, m }. Analogous to the athermal MRT LBE model, Yoshida and Nagaoka (2010) proposed a D3Q7 thermal LBE model (Fig. 31) with an MRT collision operator for the CDE, in which they have carried out an asymptotic analysis based on diffusive scaling to show that the TLBE model is first and secondorder accurate in time and space, respectively. This MRT TLBE model is also used in the present work. The TLBE governs the set of temperature distribution functions { g( x t ) = 0, 1, m }, from which t he leading order solution of the macroscopic temperature is obtained with 0mg In the D3Q7 TLBE model, the discrete velocity set is defined as e = (0,0,0), (0), (1,0,0),(0,1,0),(0,0,1),(1,2,3,4,5,6), (3 4) The evolution of the TLBE can be written as g( x + e t t + t ) g( x t ) = L g( x t ) + G t (3 5) PAGE 73 73 where is the time step, is the weight coefficient defined to satisfy the CDE in Eq. ( 3 1), and L is the collision operator. In the SRT BGK model, the collision operator is represented by L g( x t ) = 1 ( g( x t ) eqg ( x t )), (3 6) where is a coefficient that represents the relaxation time relative to the time step. In the BGK thermal model by Yoshida and Nagaoka (2010), the following equilibrium distribution functions are proposed to recover the CDE: eq()j g (3 7) where is the grid spacing in the lattice. The weight coefficients are set to be = 1/4,(0),1/8,(1,2,3,4,5,6), (3 8) and the constant = 1/4 for the 3D model. In the MRT model ( Yoshida and Nagaoka 2010), the collision operator is given by L g ( x t ) = M1SMQ g ( x t ). (3 9) Thus the evolution of the TLBE becomes g( x + e t t + t ) g( x t ) = [ M1SMQ g ( x t ) ] + G t (3 10) and the matrix Q is explicitly given as below based on Eq. (37) Q = 016[(,,...,)T jjjjeT001166(,,...,)](1,1,1,1,1,1,1) I77, (3 11) with I77 being the 77 identity matrix. The transformation matrix M for the D3Q7 TLBE model was given as ( Yoshida and Nagaoka 2010) PAGE 74 74 M = 1111111 0110000 0001100 0000011 6111111 0221111 0001111 (3 12) and the relaxation matrix S is S1 = 0456000000000000000000000000000000000000xxxyxzxyyyyzxzyzzz (3 13) As pointed out by Yoshida and Nagaoka (2010), the off diagonal components in the relaxation matrix can take nonzero values to account for anisotropic diffusion. The relaxation coefficients are related to the diffusion coefficient matrix as 21 2()ijij ij (3 14) in order to recover the solution of the CDE to the leading order. In the above, ij is the Kroneckers delta. As in the athermal MRT LBE model, the equilibria of the conserved hydrodynamic moments and the nonconserved kinetic moments, are explicitly defined (dHumires et al. 2002). The equilibria of the conserved thermal mom ents and the nonconserved kinetic moments { eqm  = 0, 1, m } in the MRT TLBE model can also be explicitly specified, as demonstrated in the D2Q5 model (Mezrhab et PAGE 75 75 al. 2010). In the following we show that the D3Q7 TLBE model ( Yoshida and Nagaoka 2010) can be equivalently formulated with the equilibrium moments explicitly defined. The thermal LBE is rewritten as g( x + e t t + t ) g( x t ) = [ M1S ( m ( x t ) meq( x t ))] + G t (3 15) where the transformation matrix M a nd the relaxation matrix S are the same as that in Eqs. (3 12) and (313), respectively; the thermal moments m ( x t ) are related to the distribution functions g ( x t ) by m = M g (3 16) and the equilibria of the moments are defined as meq = eqeqeqeqeqeqeq0123456(,,,,,,)Tmmmmmmm = uTvTwTaTT(0,,,,,0,0) (3 17) where u v, w are the velocity components and a is a constant related to the coefficients by a = (7 0 1) = 3/4. (3 18) Now the conservation of the first moment results in the macroscopic temperature m0 = T = 60g (3 19) Comparing Eq. (3 15) with Eq. (310), we can relate these two equivalent forms of evolution by m ( x t ) meq( x t ) = MQ g ( x t ). (3 20) Similarly, keeping only the first five discrete velocity components from Fig. 3 1 for the D2Q5 TLBE model yields the following from Eq. (38) e = (0,0),(0),(1,0),(0,1),(1,2,3,4), (3 21) PAGE 76 76 and the transfor mation from the distribution functions g to the thermal moments m now becomes m = 01234mmmmm = 012341111101100000114111101111ggggg = M g (3 22) The corresponding equilibria of the moments are defined as meq = eqeqeqeqeq01234(,,,,)Tmmmmm = (0,,,,0)TuTvTaT (3 23) In the D2Q5 model, the weight coefficients are chosen to be = 1/3,(0),1/6,(1,2,3,4), (3 24) and the constant a in Eq. (3 23) is determined by a = (5 0 1) = 2/3. (3 25) The relaxation matrix S for the D2 Q5 model now becomes S1 = 034000000000000000000xxxyxyyy (3 26) and the conservation of the first moment results in the macroscopic temperature m0 = T = g40 (3 27) For efficient computation, the TLBE is solved in two steps: collision step g ( x t ) = g( x t ) [ M1S ( m ( x t ) meq( x t ))] + G t and (3 28) PAGE 77 77 streaming step g( x + e t + ) = g ( x t ), (3 29) where g represents the post collision state. 3.3.3 Thermal Boundary Condition Treatment s The lattice cut by the boundary wall with arbitrary intersection distance is illustrated in Fig. 3 2. In the present work, the lattice spacing is set to be unity ( = = 1) and the fraction of the intersected link in the field to be solved is denoted xf xw / xf xe. Boundary condition treatment for the TLBE simulation is required during the streaming step when unknown populations from the undefined exterior nodes ( xe in Fig. 3 2) to the field nodes ( xf) are encountered. It should be noted that the CDE can also be used to govern the energy equation in the solid phase by setting the velocity to zero; thus the TLBE method can be extended to simulate the thermal field in the solid phase. Hence we replace the conventional fluid node wit h the more general field node notation in the computational domain. The standard momentum bounceback scheme is derived for the case where the distance from the boundary wall to the nearest field lattice node in the computational domain is equal to hal f of the lattice spacing. The basic idea of the momentum bounceback can be applied to derive the thermal boundary conditions in TLBE. Considering the (2010) proposed the following treatments for the Dirichlet a nd Neumann boundary conditions for walls perpendicular to one of the lattice velocity directions in the D3Q7 model. An asymptotic analysis was also given in their work to show that secondorder accuracy is obtained for both conditions ( Yoshida and Nagaoka 2010) PAGE 78 78 Dirichlet condition: gt (,)(,) xx (3 30) Neumann condition: gt (,)(,)(/) xx (3 31) where n is parallel to e and e If n is not parallel to e one would need to replace n in Eq. ( 3 31) with n which in turn involves the unknown heat flux in the tangential direction normal to e It is important to have secondorder accurate boundary conditions for the Dirichlet and Neumann problems back idea on the boundary and spatial interpolation, one can develop various boundary schemes with different number of neighboring field nodes and their populations in different discrete directions. The schemes i n Eqs. (3 30) and (331) above use only one post collision value g ( xf, t ) at the node adjacent to the boundary in the e direction; thus it is considered as a localized scheme. The multi reflection schemes (Ginzburg 2005) utilize five post collision neighboring values involving three field nodes, and thus up to third order accurate boundary condition schemes were derived for the Dirichlet condition. In this work we show that secondorder accuracy for both the Dirichlet and Ne umann boundary conditions can be achieved by using only three post collision neighboring values, g ( xf, t ), g ( xff, t ), and g ( xf, t ) for each unknown population g ( xf, t + ). Furthermore, we compare the present secondorder Dirichlet boundary condition treatment with the thirdorder scheme (Ginzburg 2005) by examining the errors of the temperature fields, the boundary and interior heat fluxes with numeri cal examples and PAGE 79 79 observe that: i) both show quadratic convergence since the TLBE model is secondorder accurate in the interior space; ii) the thirdorder scheme yields close, quantitatively even higher for some cases, numerical errors comparing with the present treatment under the same conditions; iii) both give very similar performance for curvedboundary simulations such as quadratic convergence of the temperature field and linear convergence of the heat flux on the boundary; thus the degradation of th e accuracy due to the curvedboundary cannot be mitigated by increasing the order of accuracy on the boundary by using more lattice nodes and populations in the field. The construction of the present boundary condition treatment is given next in detail. 3. 3.3.1 Dirichlet boundary condition As shown in Fig. 32, the present treatment for the Dirichlet boundary condition is given by gt 1234 (,)(,)(,)(,) xxxx (3 32) where cd 1, cd 2, cd 3, and cd 4 second order accuracy, the following relations are required, leaving cd 1 adjustable (see Eq. ( B 19) in Appendix B ) 12131412222121dddddddccccccc (3 33) Furthermore, in order to minimize num erical instability one may require 1dic ( i = 1, 2, and 3), (3 34a) PAGE 80 80 so that the following constraint may be desired 111dc (3 34b) A stability study with a numerical example in Sec. 3.4.1 shows that ( 3 34b) would always guarantee numerical stability for reasonable parameters used. It is also shown that while the upper bound for the numerical stability region is indeed cd 1 = 1, the lower bound can be extended from 1 to certain values depending on the relaxation coefficient ij Note that further effort can be devoted to construct more constraints on the parameter cd 1 by using more distributions on extra lattice nodes in the interpolation, such as the local thirdorder Dirichlet schemes provided in (Ginzburg 2005). However, considering the LBE is secondorder accurate in the interior lattice nodes, we aim at only obtaining a secondorder boundary treatment scheme in this work. Furthermore, the choice of the relaxation coefficient s ij will affect the error terms in the TLBE solution. A comprehensive study on the effects of the relaxation coefficients and optimization of the parameters in Eq. (332) is thus beyond the scope of this work. To demonstrate the applicability and the second order accuracy of the present boundary condition treatment, we provide the following three particular choices for the Dirichlet condition treatment. Scheme 1: gtgtgtgtgt(2(,)111()(,)(1)(,)()2 xxxxx (3 35a) Scheme 2: PAGE 81 81 gt 2(2(,)2( 2xxx gt 2 2()(,)() 2 x and (3 35b) Scheme 3: gt2(,)(,)(,)(,)2xxxx (3 35c) The three schemes in Eq. ( 3 35) all have secondorder accuracy in space and they reduce to the ordinary bounceback scheme in Eq. ( 3 all computationally stable. One may also observe that S cheme 1 is very similar to the hydrodynamic boundary condition treatment (Bouzidi et al. 2001) which is based on the bounceback rule and linear interpolation. The present analysis and the derivation in Appendix B can also offer some insights into the understanding of the various interpolation schemes for velocity boundary conditions as those in ( Ginzburg and dHumires 2003, Chen et al. 1996, Filippova and Hnel 1998, Mei et al. 1999, Mei et al. 2000, and Bouzidi et al. 2001) We notice that, however, a rigorous secondorder accurate scheme for the Dirichlet velocity boundary condition cannot be similarly formulated, even though for some special cases such as pressuredriven c hannel flows second coefficients carefully chosen based on interpolation. The reason is that in athermal LBE simulations, in principle, the leading order and first order moments of the velocity distribution functions must be able to recover the mass and momentum equations, respectively, which necessitates more constraints on the velocity distribution functions. Some of the relationships for the temperature distribution functions i n TLBE do not have PAGE 82 82 their corresponding forms for the velocity distribution functions in athermal LBE simulations. Ginzburg showed that with five populations incorporated, a thirdorder Dirichlet condition treatment can be constructed with one adjustable coefficient. With the present notation, the oneparameter family of thirdorder conditions (Ginzburg 2005) becomes gt 222111 33(,)[(122 22111 1(222 211 1(22 2xxxxxx (3 36) where xfff = xff + e = xf + 2 e and is an adjus table coefficient. Two particular solutions were provided in Ginzburgs work by choosing = 0 and = 2(12 in Eq. (3 36). In addition, the particular secondorder Dirichlet condition applied in (Ginzburg 2005) used the same three populations as in the present treatment and can be obtained with cd 1 = (1/2 34), which yields gt 1 11 (,)( 2 22 xxxx (3 37) 3.3.3.2 Neumann boundary condition In order to develop an appropriate boundary condition for a curved wall with an aligned with e and e as in Eq. (3 31). Hence, the boundary condition for gt(,)x can be represented as gt1234(,)(,)(,)(,)(/)xxxx (3 38) where cn 1, cn 2, cn 3, and cn 4 n PAGE 83 83 Different from the adjustability of the treatment for the Dirichlet condition, the second order accurate scheme for the Neumann boundary condition is determined by (see Eq. ( B 28) in Appendix B ) 12341222222nnnncccc (3 39) Thus Eq. (3 38) can be explicitly expressed as gt 2 (,)(,)(,)(,) 2 xxxx (3 40) Equation (340) reduces to the ordinary bounceback scheme in Eq. (331 0.5 and its coefficients satisfy the stability requirement 1nic ( i = 1, 2, and 3). Remark 1: The present Dirichlet and Neumann boundary condition treatments given by Eqs. (332) and (338) with their coefficients determine d by Eqs. (333) and (3with their normal direction parallel to the discrete velocity direction e For curved boundaries, Eq. (332) can be directly extended for Dirichlet problems and Eq. (338) constitutes an excellent starting point for general Neumann problems. Remark 2: Eq. (3 38) is also valid for Dirichlet problems, and Eq. (332) is also applicable to Neumann problems. Thus the combin ation of Eqs. (3 32) and (338) gives a means to evaluate heat flux on the wall along e direction for Dirichelt problems and to evaluate the wall temperature at the cuts for Neumann problems. With the PAGE 84 84 coefficients given by Eqs. (333) and (339), the following expressions can be obtained on the wall intersected by the lattice vector e d n ddc gtgtgt 1 111 22 2 11 xxx (3 41a) or d d dd gtgtgt 1 111222 111 xxx (3 41b) where cd 1 is the aforementioned adjustable var iable in Eq. (332) and its values based on the three different choices in Eq. (335) will be examined subsequently. Remark 3: The temperature gradient within the interior of the field can be evaluated using the formula ( Yoshida and Nagaoka 2010) 2 01 ()m i ij jv x (3 42) This allows for an assessment of the error propagation in terms of heat flux into the field due to the boundary condition treatment. Remark 4: For Neumann problems with curved boundaries, the heat flux n in the discrete velocity direction must be obtained based on the given n in the absence of explicit knowledge on the tangential derivative, t along the wall for curved boundaries (see Fig. 3 3). A practical strategy for the numerical determination of n is discussed next. PAGE 85 85 Consider a curved boundary wall intersected by the lattice as shown in Fig. 33. The wall cuts the horizontal lattice link along the e direction between xe and xf at xw, and the heat flux along the link direction ( e ) from exterior to interior, n is to be derived. The angles between e and the wall tangent and normal are t and n respectively. Along e the two neighboring field lattice nodes are xf and xff. Along the direction perpendicular to e we designate x f = xw and x ff = xw + e in which e is the lattice vector from exterior to the interior along the direction perpendicular to e The projection of the boundary fluxes gives n In general, t is needed in the calculation of n unless cos0t (i.e. the wall is perpendicular to e ) or t is implicitly known in very special cases. However, well posedness of the mathematical problem for CDE does not even allow the specification of t in addition to n on the same boundary. Thus the Neumann boun dary condition, or the more general mixed boundary condition in the form of 123 n ( 10 ) (3 43) involving the normal derivative on the boundary requires a fundamentally different approach than the Diri chlet problem. By rearranging Eq. (343), the boundary heat flux in the normal direction can be written as ntt tdDDDn332211 11 (3 44) The standard Neumann condition is recovered by setting 20 in the above. PAGE 86 86 For brevity, only the D2Q5 TLBE model in 2D is considered here as the extension to 3D is straightforward. As in Fig. 3 3 the normal heat flux can be expressed as n n (3 45) where n is the angle between e and the normal n Additionally, we can eliminate gt(,)x by combining Eqs. (332) and (338) to obtain a relation between d and the heat flux n in the discrete velocit y direction e on the boundary (see Eq. 341b) ddnd c 411 22 33 4 xxx (3 46a) It is important to recognize that a similar relation can be derived between d and the heat flux n in the e direction as well. That is ddnd c ''''''''''' 411 22 33 4 xxx (3 46b) where dic nic are the coefficients corresponding to dic nic ( i = 1, 2, 3, 4) defined in Eqs. (3 33) and (339), respectively, with the fraction = 0 in the e direction. The distribution functions gt'(,)x gt'(,)x and fgt'(,)x at f'x and ff 'x can be obtained using a quadratic interpolation or extrapolation of the distributions at the closest interior lattice nodes. Further elimination of d from Eqs. (3 46 a b) results in the following relation between n and n nd dccgtccgtccgtc c11 22 33 4 41 [()(,)()(,)()(,)(/) xxx PAGE 87 87 nd dccgtccgtccgtc c'''''''''' 11 22 33 4 41 [()(,)()(,)()(,)(/) xxx (3 47) Finally, combining Eqs. (345) and (347) leads to the desired expression for n in terms of the normal heat flux and the post collision distributions in the neighboring lattice nodes, n dccgtccgtccgtc'''''''''11 22 33'41xxx nd dccgtccgtccgtc11 22 3341[()(,)()(,)()(,)]sinxxx nnn nn dddccc c '' 444 '' 444 (3 48) Equati on (348) can be directly applied to Neumann problems on curved boundaries. For the special case where t or cos0t a simpler formula can be used to relate n to n n (3 49) The above is equivalent to that used by Huang et al. (2006), but it should not be used in general. For mixed boundary conditions, n is not specified. Instead, n is related to d through Eqs. ( 3 44 ) and ( 3 46) Substitution of Eq. (3 46a) into Eq. (344) yields nttnd dDDccgtccgtccgtc 3211 22 33 4114 xxx (3 50) Combining Eqs. (348) and (350) gives PAGE 88 88 n dccgtccgtccgtc'''''''''11 22 33'41xxx nd dccgtccgtccgt c11 22 33 41 [()(,)()(,)()(,)]sin xxx n t nd ddc D ccgtccgtccgt c 42 11 22 33 414 [()(,)()(,)()(,)] xxx nnnnntn dddddccccc Dc' ''234 4 4 4 42' ''414 4 4 41}/[sin cos()] (3 51) Once n is obtained the boundary condition treatment in Eq. (340) can be used to complete the streaming step. It is noted that n is not needed in Eq. (3 40) as xw in Fig. 3 3 is located on the link along the e direction. The same procedure can be applied to derive n when xw is located on the link along the e d irection. For more general cases such as those lattice models that are more complicated than D3Q7 (for 3D) or D2Q5 (for 2D) and when a lattice velocity vector e is not in the axial x, yor z direction, the same proce dure can be used to compute the axial heat fluxes nx ny and nz in a coupled manner first. The nonaxial heat flux n can be subsequently obtained using nx ny and nz This method of determining n for Neumann or mixed boundary condition is hereinafter referred to as Cartesian decomposition method. The accuracy of the present boundary condition treatments will be examined through a series of numerical tests in the next section, and the asymptotic analysis is given in Appendix B When the node xff is outside the field, we recommend using the localized schemes by Yoshida and Nagaoka (2010). PAGE 89 89 3.4 Numerical Validation and Discussion To verify the applicability and accuracy of the present thermal boundary condition treatments, we have conducted several numerical tests involving both 2D and 3D geometries. Both conv ection and diffusion are considered in the first example, where the Dirichlet and Neumann boundary conditions are separately applied and investigated. The second example shows how to convert the Neumann condition in the normal direction of an inclined wall into the heat fluxes in the discrete velocity directions when there is no tangential flux. From the third example through the fifth example, curved geometries are considered. Steady state and transient boundary conditions are given in the third and fourth problems; while the fifth test is for the simulation of steady state convection and diffusion inside a circular (3D) pipe. In this study, only isotropic diffusion is considered, i.e., Dij = Dtij; thus the off diagonal components in the relaxation matri x are set to be zeros, and the constraint for the diagonal coefficients is maintained so that Eq. (314) can be rewritten as 21 () 2()ijDij tij (3 52) The positivity of the diffusion coefficient Dt requires that D > 1/2. A typical value of D = 0.75 is used in each of the numerical examples to demonstrate the accuracy of the present boundary condition treatments. Although no detailed stability analysis is performed for the proposed Dirichlet and Neumann boundary schemes, we have numerically c hecked their stability by using D values close to 1/2, for which numerical instability is most likely to occur. For brevity, the results for D = 0.55 and D = 0.51 are presented for some cases. The relaxation coefficient 0 for the conserved quantity has no e ffect on the numerical solution, and the remaining adjustable coefficients do not PAGE 90 90 contribute to the leading order solution of the CDE since they only affect the higher order error terms of the LBE solutions; thus 0 = 0, and P = 1 ( p = 4, 5, 6 in 3D and p = 3, 4 in 2D) are used in the rest of this work. 3.4.1 Two D Steady State Channel Flow First we consider a plug flow of velocity U in a 2D channel of constant height, H (Fig. 3 4). Each wall is placed at a distance of ( y) from the first interior latt ice in the y direction. The CDE for the temperature inside the channel is 2222()tTTTUDxxy (3 53) This thermal flow problem is characterized by the Pclet number Pe = HU/ Dt; and Pe = 20 is used in this study. Periodic boundary conditions are imposed in the xdirection so that f ( x+ L ) = f ( x) for both the temperature T and distribution functions g. Both Dirichlet and Neumann thermal boundary conditions on the wall are considered. 3.4.1.1 Dirichlet boundary condition On the wall at y = 0 and y = H T ( x, y = 0) = T ( x, y = H ) = cos( kx), with k = 2 / L (3 54) An analytical solution for the thermal field can be obtained. When expressed with complex variables, the exact solution is Tex( x, y) = Real ikx ee e e 1 (3 55) where Real means taking the real part of a complex variable and 1tiUDk Before evaluating the accuracy of the present Dirichlet schemes, the numerical stability is investigated first for this problem. The work of (Mei et al. 1999) indicates that PAGE 91 91 the numerical instability associated with boundary condition treatments in using the LBE method is m ore pronounced in channel flows than in curved geometries. Hence guidance on numerical instability from solving the thermal problem in the channel flow will be useful. Fig. 3 5 (a) shows the lower stability bound for the adjustable coefficient cd 1 at Pe = 20 and Ny = 34 ( Ny denotes the number of lattice nodes in y direction). Very close results for the lower bound are obtained also for ( Pe Ny ) = (20, 130), (100, 34) and (100, 130). The upper bound for cd 1 is 1.0 for all cases, noting that dc11 as required in Eq. (333). An asymptotic lower bound of cd 1 tested as D approaches 0.5 in Fig. 35 (a). On the other hand, when D increases, the lower bound drops, resulting in a broader stable region while the ove rall relative error of the temperature field increases (see Figs. 37 and 38 below). Fig. 35 (b) shows the dependence of cd 1 0 for the three schemes in relation to the lower stability bound D = 0.50001 and the upper stability bound. Clearly, all three schemes are well within the stable region. As an illustration, Fig. 3 6 shows the temperature contours obtained using D = H = 64. To assess the accuracy of the LBE solution of the temperature TLBE the following relative L2 error norm is defined: E2 = 1/2 2 ,,()/xy xyTTT2 LBEex ex (3 56) Figure 37 (a f) shows the log log plots of E2 versus the dimensionless grid resolution, 1/ H over a complete range of for the three schemes given in Eq (3 35 a c). Three relaxation coefficient values D = 0.51, 0.55 and 0.75 are used. Second order accuracy is observed for all cases. At D = 0.51, all three schemes give very close results as PAGE 92 92 shown in Fig. 37 (a, b). As D increases, the differences among the three schemes become more noticeable and Scheme 2 gives smallest relative error for each case as shown in Fig. 37 (c f). The variations of E2 in the range of from 0 to 1 are shown in Fig. 3 8 with Ny = 34. It is noticed that with D = 0.75, all the three schemes have back scheme is applied. When D increases to 1.0, the relative error also increases for each scheme when ecomes more complicated. The results in Figs. 3 7 and 38 clearly verify the existence of the adjustable parameter cd 1 in using the TLBE method for CDE subjected to Dirichlet boundary conditions; and the proposed Dirichlet schemes are stable and of secondorder accurate for D > 0.5. For comparison purposes, the results obtained using the thirdorder and second order schemes (Ginzburg 2005) for = 0.25, 0.5 0, and 0.75 are compared with that from Scheme 2 in Fig. 3 9 (a, b) where the Ginzburg 3rdorder # 1 and Ginzburg 3rd order # 2 represent the two particular schemes whose coefficients are constructed with = 0 and = ()/ 2122 in Eq. (3 36), respectively, and the coefficients for the secondorder scheme are given in Eq. (337). Fig. 3 9 demonstrates that the thirdorder Dirichlet schemes (Ginzburg 2005) also lead to secondorder convergence since the TLBE method is secondorder accurate in the interior In addition, Scheme 2 gives smaller error than the local thirdorder schemes in (Ginzburg 2005) for each value tested in Fig. 3 9 (a, b). compared with those obtained using Ginzburgs schemes. To gain further insight into the impact of the boundary condition treatment on the accuracy of the TLBE solution, we compute the boundary heat flux and interior temper ature gradient using the formulas given by Eqs. (341a) and (342), respectively. PAGE 93 93 The relative L2 norm errors are computed using the respective derivatives obtained from the analytical solutions given by Eq. (355), E2_qw = 1/2 22[()()]/() TTT yyyLBE ex ex boundarynodes boundarynodes (3 57) E2_grad = 1/2 22[()()]/() TTT yyyLBE ex ex interiornodes interiornodes (3 58) The results of E2_qw for the wall heat flux and E2_grad for the interior temperature gradient at D = 0.75 are shown in Figs. 310 and 311, respectively. It is clear that the numerical solutions for boundary heat flux and interior temperature gradient are of secondorder accuracy for Dirichlet problems for straight walls perpendicular to one of the l attice vectors. Scheme 2 gives slightly smaller error than the other two. 3.4.1.2 Neumann boundary condition For the Neumann problem, the normal heat flux on the wall is given as cos()/wtyHtT qDDkxH y (3 59) The exact solution for the temperatur e is Tex( x, y) = Real ikxeee1 (3 60) Again, the numerical solution of the temperature field at D H = 64 is shown in Fig. 3 12 as an example; in addition to the interior temperature and temperature gradient error norms defined in Eqs. ( 3 56, 3 58), the relative L2 error norm for the boundary temperature is defined as PAGE 94 94 E2_tw = 1/222()/TTTLBEex exboundarynodes boundarynodes (3 61) Figures 313 (a c) show E2 (Eq. (3 56)) versus 1/ H for ranging from 0.01 to 0.99 and D = 0.51, 0.55 and 0.75 at Pe = 20, and Figs. 3 1 4 and 31 5 show E2_tw and E2_grad, respectively, for D = 0.75. The results confir m that the present Neumann condition scheme is stable, and secondorder accuracy is achieved for the temperature field, boundary temperature, and interior heat flux when the straight wall is perpendicular to one of the lattice vectors. Extreme cases including D 0.001, 0.999, 0.0001, 0.9999 at Ny = 34 have also been examined. The computations using present Neumann boundary scheme are all stable and give good numerical results for all cases. 3.4.2 OneD Transient Conduction in an Inclined SemiInfinite Solid In order to demonstrate the applicability of the present Neumann boundary condition treatment to inclined walls, we consider the transient heat conduction of a constant heat flux into an inclined semi infinite solid as shown i n Fig. 3 1 6 For a given wall inclination angle with respect to the yx y are calculated for each field point adjacent to the boundary. Note there is still no tangential variation in this case; the Neumann boundary condition with nonzero tangential gradient will b e investigated in Sec. 3.4.3. The initial temperature is zero in the field (solid), and the constant heat flux condition is imposed as 02twtlDTqDn (3 62) Hence the transient temperature along the solid is solved to be PAGE 95 95 2/41/2(,)[ ()]22tlDtttlTtlteDtDtexerfc (3 63) where l is the distance in the normal direction and erfc is the complementary error function. The temperature on the wall is f ( t ) = T ( t l = 0) = t1/2. With no temperature variation along the tangential direction, the boundary heat f luxes in the discrete velocity directions can be directly calculated using Eq. (349), 2coscos 2t wD qq (3 64a) and 3sinsin 2t wD qq (3 64b) Figure 31 7 (a) shows the variations of the temperature with time at point P in Fig. 3 1 6 for inclination angles = tan1(1.0), tan1(1.2), tan1(1.5), and tan1(2.0). For each inclination angle, the temperature profile along the vertical line going through point P at t = 500 is compared with the exact solution in Fig. 31 7 (b), where the distances in the y coordi nate have been projected into that in the normal direction of the physical boundary. Very good agreement between the simulation and the analytical solution is noticed in Figs. 31 7 (a) and (b). These results verify that Eq. (349) is correct when the tangential temperature gradient on the boundary is zero. This test also serves as a good example to show that the present Neumann boundary condition treatment can be readily implemented for a boundary that is intersected by the rectangular lattice with arbitrar also the subject of the next example. PAGE 96 96 3.4.3 Two D Steady State Heat Conduction Inside a Circle Figure 31 8 (a) shows the schematic layout of the lattice covering the entire surface and the inter ior of a circle of radius r0. The specific x or y value of each including both x and y, as a function of the azimuthal angle of the boundary nodes at r0 = 30.5 is depicted in Fig. 31 8 (b). 3.4.3.1 Dirichlet boundary condition First we consider the steady state Dirichlet boundary condition T ( r = r0) = f ( ) = cos( n ) ; n = integer. (3 65) The exact solution for the temperature field is Tex( r ) = ( r / r0)ncos( n ). (3 66) Figure 31 9 shows the relative L2 norm error E2 defined in Eq. (356) for all the interior nodes versus the grid resolution, 1/ r0, on the log log sc ale with n = 4 used in Eq. (3 65) The present S chemes 1, 2 and 3 are based on Eq. (335 a c ) respectively, and the secondorder and thirdorder schemes from (Ginzburg 2005) are used for comparison (those schemes are rewritten in Eqs. (336, 337)) Sec ondorder accuracy for the temperature field is observed in Fig. 31 9 for each scheme examined. It is seen that Scheme 2 gives better accuracy comparing with the other two given by Eq. ( 3 35). The results for the error norm E2_qw of the boundary heat flux n according to Eq. ( 3 41a) are shown in Fig. 3 20, where on average only first order convergence is observed. It should be noted that in Fig. 3 20 the results of Ginzburgs thirdand second order schemes are obtained by combining Eqs. ( 3 36) and ( 3 37) with our PAGE 97 97 second order Neumann scheme in Eq. ( 3 40), respectively, and following the same procedure as that for Eq. ( 3 41a). Fig. 3 21 gives the results of the error norm of the interior temperature gradient E2_grad defined in Eq. (5 8). Realizing Scheme 2 has better accuracy than the other two given by Eq. (35) in the previous examples only the results using Scheme 2 are presented in Figs. 3 20 and 3 21. The finitedifference case represents a quadratic finitedifference scheme based on the numerical values of the temperature field at the lattice nodes and the given wall temperature. One valuable insight to be gained from Figs. 3 20 and 321 is that for the boundary heat flux and interior temperature gradient in Eqs. (341a) and (34 2), respectively, these various Dirichlet schemes all yield linear convergence for the boundary heat flux, and the error norm of the interior temperature gradient E2_grad is proportional to (1/ r0)1.5, both of which are lower than the secondorder convergence as shown in Figs. 310 and 3 11 respectively for straight boundaries. The degradation of the convergence rates for the derivative related physical quantities is due to the curved boundary at which the int (see Fig. 31 8 ( b ) ). On the whole, the present Dirichlet condition treatment on a curved boundary is able to produce secondorder accurate temperature field and first order boundary heat flux, and the interior temperature gradient converges to the exact solution with a superlinear order close to 1.5. This conclusion also offers an insight into the accuracy of the Neumann boundary condition treatment in curvedboundary situations, which will be discussed next. 3.4.3.2 Neumann boundary condition Now consider the following Neumann (flux) condition applied on the boundary PAGE 98 98 00cos()/wtrtT qDnDn r (3 67) The solution of the temperature field is the same as that in Eq. (366) and the case of n = 4 is used. To account for the contribution to n from the tangential temperature gradient, Eq. (3 48) is used. After n is obtained numerically, Eq. (340) is implemented to complete the streaming step for Neumann problems. The error norms E2, E2_tw and E2_grad of the interior temperature field, the boundary temperature, and the interior temperature gradient as defined in Eqs. (356, 3 61, and 358) versus the grid resolution, 1/ r0, are shown in Figs. 322 3 24, respectively, where the present coupled schemes 1, 2 and 3 represent the different coefficient sets cdi and 'dic in Eq. (3 46) determined by the three schemes in Eq s. (3 35 a c), respectively. From the analytical soluti on for the temperature field, Eq. (366), the exact boundary heat flux n in the discrete velocity direction can be obtained with projection of the normal and tangential heat fluxes into the xor ycoordinate. For comparison, the results of the three types of error norms E2, E2_tw and E2_grad based on this exact boundary flux are also presented in Figs. 322 3 24 respectively. First of all, we notice that with the exact boundary flux n applied, the us e of Eq. (3 40) for the evaluation of the distribution function leads to secondorder accuracy for all the three quantities, i.e., the interior and boundary temperatures and the interior temperature gradients. However, only first order accuracy is obtained for both the interior and boundary temperatures for all the three coupled schemes. The average decrease of the error norms in Figs. 322 and 323 is proportional to the grid resolution, PAGE 99 99 1/ r0, on the log log scale and considerable oscillations are observ ed even at high grid resolution. The interior temperature gradient as shown in Fig. 324 still converges with a superlienar order around 1.5, which is similar to the result of E2_grad in Fig. 3 21 with a Dirichlet condition. Compared with the secondorder convergence in the case of the exact boundary flux, the degradation of the order of accuracy in the coupled schemes is due to the heat flux approximation of n curved boundary. Recal ling in Fig. 3 20 only a linear convergence of the boundary heat flux is obtained for the Dirichlet problem, and the heat flux formula Eq. ( 3 41a, or 3 46a) is used in deriving the ultimate heat flux approximation in Eq. ( 3 48), i t is believed that some se cond order errors terms O ( 2x introduced in the heat flux approximation. Those secondorder errors cannot be always cancelled out due to the random distribution s (Fig. 3 18 (b)) resulting in averaged linear convergence and p ersistent oscillations even at high resolution as shown in Figs. 322 and 3 23 Based on the multi reflection approach, Ginzburg (2005) also found that only linear convergence can be achieved for the Neumann condition on the curved boundary when tangenti al temperature gradient exists. The orders of accuracy for the Dirichlet and Neumann conditions on straight and curved boundaries for the temperature and its derivative related quantities are summarized in Table 1 for 2D steady thermal problems. 3.4.4 Two D Transient Heat Conduction Inside a Circle In this case, a transient boundary condition is applied on the curved boundary to examine the temporal accuracy of the present boundary condition treatment for timedependent problems. The geometry is the same as that in Fig. 3 18 (a). PAGE 100 100 The temperature is uniform in the azimuthal direction and the transient Dirichlet boundary condition reads T ( t r = r0) = f ( t ) = sin( t ), (3 68) where is the frequency of the imposed thermal boundary condition. This transient problem is characterized by the Stokes number St defined as 2202ttrS (3 69) The analytical solution for this 2D problem can be fo und as ( Hahn and Ozisik 2012) 22()001010()2(,)() [(0) ()]()tn tntDnnnnJTrtft feedfrex 22202221010()cos()sin()2sin() [ ]()()tnDntn tnnnn tnr (3 70) where Ji( ) (i = 0, 1) is the solution to the Bessels function of the first kind of order i and n is determined by 00()0nJ For error assessment purposes for the this transient problem we define two error norms E ( t*) = 21/2 ,1 [()]xy xyTT NNLBEex (3 71) E2 = 2/ 121/2 2*1/200,,11[()][()]2xy xyxy xyTTdt TTdt LBEex LBEex ,(3 72) where t* = /2 and the product of ( NxNy) accounts for the total number of the lattice nodes inside the circle. To eliminate the effect of initial transient, the computation is performed for a sufficiently long period of time ( t > 1). The Stokes number St is set to be unity in all PAGE 101 101 simulations. The results of E ( t *) defined in Eq. ( 3 71) at D = 0.51 and 0.75 are presented in Fig. 3 25 using the three schemes in Eq s. ( 3 35 a c). The grid resolution here is relatively poor as r0 = 5.8. The average error E2 defined in Eq. ( 3 72) is evaluated over one single period after the dynamic steady state is reached. The results of E2 with different resolution are shown in Fig. 3 26 as a function of 1/ r0 for D = 0.75. The decrease of the error norm E2 is proportional to (1/ r0)2 on the loglog scale for each of the three schemes. It can also be observed in Fig. 3 26 that Schemes 2 and 3 have very close errors and they are more accurate than Scheme 1 in general. When D = 0.51 is used, the three schemes give almost identical results for E2 and the magnitude is very close to that in Fig. 3 26 for D = 0.75 thus they are not shown here for brevity. Overall, Scheme 2 appears to give smaller error than the other two schemes given by Eq. (35). It is thus confirmed that the present Dirichlet boundary condition treatment can be applied in transient problems with cur ved geometry, and secondorder accuracy is preserved for the temperature field. 3.4.5 ThreeD Steady State Circular Pipe Flow To test the TLBE model and the boundary condition treatments in 3D cases, the convection and diffusion inside a circular pipe are considered in this example. The temperature is assumed to be uniform in the azimuthal direction in order to obtain analytical solutions. 3.4.5.1 Dirichlet boundary condition The plug flow with constant velocity U in the axial z direction is assumed. The p eriodic boundary condition in the z direction is applied as T ( r z = L ) = T ( r z = 0), (3 73) and the Dirichlet boundary condition is imposed as PAGE 102 102 T ( r = r0, z ) = cos( kz ), with 2/k (3 74) Following the same procedure for 2D cha nnel flow with Dirichlet boundary conditions in Sec. 3.4.1.1, the analytical solution for the temperature field is solved to be Tex( r z ) = Real [ eikz 0 00() () I I ] (3 75) where I0( ) is the modified Bessels function of the first kind of order 0 and 1tiU Dk The relative error norms E2 and E2_qw are extended from 2D to 3 D cases E2 = 2 21/2 ,, ,,[()/]r TTTLBEex ex (3 76) E2_qw = 002 21/2,, ,,[(()())/()]rr TTTrrrLBE ex ex (3 77) Figures 32 7 and 32 8 show the results of E2 and E2_qw for this 3D case, respectively, versus the grid resolution, 1/ r0, with the three schemes given in Eq s. (3 35 a c) implemented. The characteristic Pclet number is Pe = (2 r0) U/Dt = 20 and the relaxation coefficient is D = 0.75. Similar to the 2 D example in Sec. 3.4.1.1, this 3D example with Dirichlet condition on the curved boundary gives quadratic convergence of the temperature field and only a linear convergence of the boundary heat flux. The secondorder accuracy of the present Dirichlet boundary condition treatment for the temperature field is therefore verified for 3D thermal flows with curved geometry. 3.4.5.2 Neumann boundary condition For the Neumann problem, the heat flux is given as PAGE 103 103 00cos()/wtrrtT qDDkzr r (3 78) and the other conditions are the same as in Sec. 3.4.5.1.The analytical solution for the temperature field is Tex( r z ) = Real [ eikz 0010()()I ], (3 79) where I1( ) is the modified Bessels function of the first kind of order 1. There is no azimuthal temperature gradient along the curved boundary; the nonzero tangential gradient is on the z direction, which is aligned with e5 and e6 of the lattice velocity vector s. Thus the exact flux conversion formula in Eq. (349) can be directly applied. The error norm E2 defined in Eq. (376) is again examined for the Neumann problem and the results using D = 0.51, 0.55 and 0.75 are shown in Fig. 32 9 For each of the D va lues tested, the present Neumann boundary condition treatment leads to a stable computation. The decrease of the relative error E2 is proportional to (1/ r0)2 on the log log scale. The temperature field obtained is thus second order accurate. Comparing wit h the linear convergence of the temperature field in Sec. 3.4.3.2, it is confirmed that with zero temperature gradient along the curved boundary, secondorder accuracy of the temperature field can be achieved with the present Neumann boundary condition tre atment for 2D and 3D problems. 3.5 Summary and Conclusions Thermal boundary condition treatments based on the bounceback scheme and spatial interpolation are presented for both the Dirichlet and Neumann conditions. When second order accuracy is pursued, there is an adjustable parameter for the Dirichlet problem, while the Neumann condition treatment is unique. Three particular schemes PAGE 104 104 for the Dirichlet condition treatment are presented and their stability and accuracy are examined. All three Dirichlet schemes are stable and Scheme 2 gives smaller error in most cases investigated. When extended to curved boundaries, the present Dirichlet condition treatment can be applied directly and it leads to secondorder accuracy of the temperature field and first order accuracy of the boundary heat flux. In contrast the proposed Neumann condition treatment requires the determination of heat fluxes in the discrete velocity directions of the TLBE model. A general formula to convert the normal heat flux into the fluxe s in the discrete velocity directions is derived. Different from the Dirichlet problem, only first order accuracy is obtained for the temperature field with Neumann conditions due to the combined effect of the curved boundary and the heat flux conversion. For the special case with no tangential variation of temperature along the curved boundary, the heat flux conversion becomes exact and secondorder accuracy of the temperature field can be obtained. PAGE 105 105 Table 31. Accuracy of the TLBE simulations with Dirichlet/Neumann conditions Dirichlet condition Neumann condition parallel straight boundary curved boundary parallel straight boundary curved boundary curved boundary Interior temperature field 2ndorder 2ndorder 2ndorder 2ndorder 1st order Boundary heat flux 2ndorder 1st order Boundary temperature 2ndorder 2ndorder 1st order Interior temperature gradient 2ndorder superlinear order ~ 1.5 2ndorder 2ndorder superlinear order ~ 1.5 the boundary normal is parallel with a lattice vector; exa ct boundary flux; coupled schemes PAGE 106 106 0 1 2 3 4 6 5 x y z Figure 3 1 Discrete velocity set { e} for the D3Q7 TLBE model. xexfxffxw field exterior boundary wall ^ ^ ^ g(xf, t + t ) g g(xf, t ) (xf, t ) g(xf f t ) Figure 32. Schematic depiction of the lattice link intersected by a boundary wall. xxexfxffxw field exterior boundary wall x'fx'ff 5 3 1 ey en t nt ntnn 2 4 8 6 7 e Figure 3 3 Layout of the rectangular lattice and a curved boundary (solid circles: field nodes; solid squares: boundary nodes; open circles: exterior nodes ) PAGE 107 107 x y L y y H eFigure 3 4 Schematic representation of the lattice in a 2v alue. (a) (b) Fig ure 35. (a) Lower stability bound of the adjustable parameter cd 1 versus ( D 0.5), and (b) variations of cd 1 Dirichlet schemes at D = 0.50001 for the channel flow Dirichlet problem (the lower and upper stability bounds are included). PAGE 108 108 Fig ure 36. Temperature contours in a 2D channel with a Dirichlet boundary condition using D H = 64. All three schemes reduce to the bo unceback condition. PAGE 109 109 (a) (b) (c) (d) (e) (f) Fig ure 37. Relative L 2 error norm E2 versus the grid resolution, 1/ H for the channel f low Dirichlet problem. (a) D = 0.51, (b) D = 0.51, D = 0.55, (d) D = 0.55, D = 0.75, (f) D = 0.75. PAGE 110 110 Fig ure 38. Relative L 2 error norm E2 channel flow Dirichlet problem at Ny = 34. (a) (b) Fig ure 39. Comparison of E2 obtained using Scheme 2 with those using Ginzburgs schemes (Ginzburg 2005) for the channel flow Dirichlet problem at D = 0.75. back scheme), (b) 0.25 and 0.75. PAGE 111 111 (a) (b) Fig ure 310. Wall heat flux error E2_qw versus the grid resolution, 1/ H for the channel flow Dirichlet problem. D = 0.75, (b) 0.01, and 0.99 at D = 0.75. (a) (b) Fig ure 311. Interior gradient error E2_grad versus the grid resolution, 1/ H for the ch annel D = 0.75, (b) 0.01, and 0.99 at D = 0.75. PAGE 112 112 Fig ure 312. Temperature contours in a 2D channel with a Neumann boundary condition using D H = 64 PAGE 113 113 (a) (b) (c) Fig ure 313. Relative L 2 error norm E2 versus the grid resolution, 1/ H for the channel flow Neumann problem with various (a) D = 0.51, (b) D = 0.55, and (c) D = 0.75. PAGE 114 114 (a) (b) Fig ure 314. Wall temperature error E2_tw versus the grid resolution, 1/ H for the channel flow Neumann problem. (a) = 0.50, 0.25, and 0.75 at D = 0.75, (b) = 0.50, 0.01, and 0.99 at D = 0.75. Fig ure 315. Interior gradient error E2_grad versus the grid resolution, 1/ H for the channel flow Neumann problem. PAGE 115 115 y x Pnx y l qnq3q2qn semiinfinite field exterior inclined wall at different angles Figure 3 16. Layout of th e lattice around the inclined semi infinite solid. (a) (b) Figure 3 17. (a) Temperature variations with time at point P ; (b) temperature profiles in the ydirection going through P at t = 500 PAGE 116 116 r0 x P y field (a) (b) Figure 3 18. (a) Schematic layout of the lattice around a circle; ( b values along the azimuthal angle at r0 = 30.5. Fi gure 3 19. Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Dirichlet condition (results based on Ginzburgs schemes (2005) are included) PAGE 117 117 Figure 3 20. Wall heat flux error E2_qw v ersus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Dirichlet condition (results based on Ginzburgs schemes (2005) are included) Figure 3 21. Interior gradient error E2_grad versus the grid resolution, 1/ r0, fo r the steady state heat conduction inside a circle with a Dirichlet condition (results based on Ginzburgs schemes (2005) are included) PAGE 118 118 Figure 3 22. Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the steady state heat conduction insi de a circle with a Neumann condition. Figure 3 23. Wall temperature error E2_tw versus the grid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Neumann condition. PAGE 119 119 Figure 3 24. Interior gradient error E2_grad versus the g rid resolution, 1/ r0, for the steady state heat conduction inside a circle with a Neumann condition. Figure 3 25. Comparison of the error norm E ( t*) defined in Eq. (3 71) for the transient heat conduction inside a circle using Schemes 1, 2 and 3 with r0 = 5.8. PAGE 120 120 Figure 3 26. Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the transient heat conduction inside a circle. Figure 3 27. Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the pipe flow with a Dirichlet condi tion. PAGE 121 121 Figure 3 28. Wall heat flux error E2_qw versus the grid resolution, 1/ r0, for the pipe flow with a Dirichlet condition. Figure 3 29. Relative L 2 norm error E2 versus the grid resolution, 1/ r0, for the pipe flow with a Neumann condition. PAGE 122 * This work has been accepted f or publication in ASME Journal of Heat Transfer (Li et al. 2013b). 122 CHAPTER 4 HEAT TRANSFER EVALUATION ON CURVED BOUNDARIES I N THERMAL LATTICE BOLTZMANN EQUATION METHOD* 4.1 Introduction Accurate and efficient evaluation of hydrodynamic forces and heat and mass transfer on boundaries and interfaces involving curved geometry is critical to the growth of the lattice Boltzmann equation (LBE) method for modeling and analyzing transport phenomena. For momentum transfer on curved boundaries, an accurate force evaluation is crucial, and there are basically two different approaches to determine the fluid dynamic forces in the LBE method: the momentum exchange approach (Ladd 1994, Behrend 1995, Mei et al. 2002, Yin et al. 2012) and the surface stress integration approach (Mei et al. 2002, He and Doolen 1997, Filippova and H nel 1998) Mei et al. (2002) compared these two different methods with numerical tests and concluded that the momentum exchange approach is reliable, accurate and much easier to implement than the stress integration approach for both twoand threedimensional fl ows. Following the momentum exchange idea for force evaluation in the LBE method in (Ladd 1994, Mei et al. 2002) an energy exchange approach for total heat transfer evaluation on curved boundaries in the TLBE method is proposed in this work. At the present time, most of the existing TLBE simulations focus on the temperature field, while the heat transfer on the boundaries and fluidsolid interfaces, especially those with curved geometries, has not been systematically studied. It is the main objective of the present study to establish an accurate and efficient numerical approach to determine the boundary heat transfer based on the TLBE temperature distribution functions PAGE 123 123 instead of approximating the heat transfer from the solved temperature field using finit e difference schemes. In this manner, the heat transfer simulation is conducted on the same level as that of the temperature distribution functions, and thus the mesoscopic nature of the LBE method is preserved in the heat transfer evaluation. To formulate the problem, the TLBE model with a multiplerelaxationtime (MRT) collision operator for the convectiondiffusion equation (CDE) proposed by Yoshida and Nagaoka (2010) is used in the present work. This approach has proved useful due to its second order accuracy in space and convenience for implementation compared with those TLBE models (Chen et al. 1994, He et al. 1998, Guo et al. 2007) that were derived from the continuous Boltzmann equation. The thermal boundary condition treatment in Chapter 3 is imple mented as it preserves the exact geometry on the boundary, and has secondorder accuracy for both the Dirichlet and Neumann conditions on straight walls and at least first order accuracy on curved boundaries (Li et al. 2013a) In the present heat transfer evaluation method, the boundary heat fluxes in the discrete lattice velocity directions of the TLBE model are directly obtained from the temperature distribution functions at the lattice nodes and the temperature values on the boundary, rather than using any finite difference schemes based on the solution for the temperature field. From the conservation of energy, integration of the discrete heat fluxes with effective surface areas gives the energy exchange on the boundary, thus bypassing the need for deter mining the boundary heat flux in the normal direction. In addition, due to the uniform lattice spacing in the TLBE model, the effective surface area for each of the discrete heat fluxes is uniform and constant. The surface area approximation in the heat fl ux integration, which is laborious in general and impractical PAGE 124 124 in complex geometries such as porous structures, is thus avoided. These two major benefits make the present heat transfer evaluation very efficient in curvedboundary simulations. It is worth noting that the present heat transfer evaluation method can be simply extended to mass transfer computations on curved boundaries as the general CDE also governs the species mass concentration. The Dirichlet, Neumann and mixed thermal boundary conditions can be simply replaced with corresponding concentration boundary conditions. A series of numerical test cases are presented to validate the proposed heat transfer evaluation scheme, including (i) 2dimensional (2D) steady state plug flow in a channel with a constant axial velocity and a sinusoidal wall temperature distribution, (ii) 1 D transient heat conduction in an inclined semi infinite solid, (iii) 2 D transient heat conduction inside a circle with a temporal sinusoidal wall temperature variation, (iv) 3D steady state thermal flow in a circular pipe with a constant axial velocity and a sinusoidal wall temperature variation along the flow direction, and (v) 2D steady state natural convection in a square enclosure with a circular cylinder located at the c enter. For cases (i) through (iv) analytical solutions are available and the convergence of numerical results shows that the proposed heat transfer evaluation scheme is secondorder accurate for straight boundaries perpendicular to one of the discrete latt ice velocity vectors, while for curved boundaries only first order accuracy of the heat transfer can be obtained due to the irregularly distributed lattice fractions cut by the curved boundary. For case (v), the surfaceaveraged Nusselt numbers computed us ing the proposed technique for heat transfer evaluation on curved boundaries agree well with published numerical results in the literature. PAGE 125 125 4.2 Heat Transfer Evaluation The D3Q7 and D2Q5 TLBE models by Yoshida and Nagaoka (2010) and the thermal boundary condition treatment by Li et al. (201 3a) are implemented in the present boundary heat transfer evaluation. The details for the TLBE models and thermal boundary condition treatment are presented in Sections 3.3.2 and 3.3.3 of Chapter 3. 4.2.1 Boundary Heat Fl ux The present heat transfer evaluation is based on the integration of the boundary heat flux n in the discrete velocity direction e with constant effective surface area in the TLBE model instead of integ rating the normal heat flux n As discussed in Remark 2 of Section 3.3.3, t he boundary heat flux n can be conveniently obtained by combining the Dirichlet and Neumann treatments in Eqs. (3 32) and (338) to yield dn ddcgxt gxt gxt 111122 2 11 (4 1) where cd1 is the adjustable v ariable discussed in Eq. (3 32) (see also Eq. (341a) for the notations) For mixed boundary conditions, n is not explicitly known; thus Eq. ( 4 1 ) for n cannot be used. Instead, the Cartesian decomposition method discussed in Section 3.3 should be employed to compute n The computed values of n at all the boundary intersections along all the relevant lattice vector directions then form the basis for heat transfer evaluation on the boundary surface. PAGE 126 126 4.2.2 Heat Transfer on Curved Boundaries Intuitively heat transfer on a boundary can be computed following a threestep procedure for g iven Dirichlet or mixed boundary conditions in the TLBE method. First, the heat flux n in the discrete lattice velocity direction is obtained along the boundary. Then the relationship between n and n can be applied to obtain the normal heat flux, n on the boundary. Integration of the normal heat flux along the boundary wall gives the heat transfer rate wnnQqdAdA (4 2 ) where dA is the local surface area. There are two main difficulties associated with the heat transfer rate integration in Eq. (15), especially for that on curved boundaries. First, the heat flux conversion from n to n is laborious in general. Second, the integration of the normal heat flux requires a local surface area approximation of the actual boundary. When the boundary nodes, i.e., the intersection locations on the boundary by the lattice, are not uniformly distr ibuted (see Fig. 4 1 ), the effective surface area approximation for the normal heat flux would require a great amount of effort for curved boundaries and might be impractical for irregular geometries such as thermal simulations in porous media. To avoid the boundary heat flux conversion and surface area approximation, we propose the following heat transfer evaluation technique based on the idea of energy exchange. Following the principle of energy conservation, the integration of the normal heat flux at a b oundary, which accounts for the energy transferred into the field, must be equal to the integration of the heat flux in the discrete velocity directions of the TLBE model with effective surface areas. For each of the boundary nodes such as xw in Fig. PAGE 127 127 4 1 t here is a discrete boundary heat flux, n or n that needs to be determined in TLBE simulations, according to whether the boundary is cut by the link along the discrete velocity vector, e or e respectively. When the boundary is intersected exactly at a lattice no n and n need to be determined at the same boundary node. Because the lattice nodes are distributed in the Cartesian coordinates with uniform lattice spacing in the standard TLBE models, the effective surface areas for n and n are ( ) and ( ), respectively, with the lattice spacing in the direction perpendicular t o the plane. When = = = 1, the effective surface area is unity for all the discrete heat fluxes, and thus the heat transfer on the boundary is simply a summation of all the discrete heat fluxes multiplied by the unit effective surface area. For ef ficient computation, two scalar arrays, w ( i j ) and wb( i j ) are introduced, where i and j are the indices of the lattice nodes in a 2D problem The extension to 3D problems is straightforward. We choose w ( i j ) = 1 for the all the field nodes xf, and w ( i j ) = 0 for all the exterior nodes xe. The array wb( i j ) is set to zero everywhere except for those field nodes that satisfy xf + e = xe (see Fig. 4 1 ) where a value of 1 is assigned. By summing the boundary heat fluxes from all discrete velocity directions into the field nodes with wb( i j ) = 1, the heat transfer across the boundary into the field over one time step is bwn w Q 2 10 xe (4 3 ) It is noticed that the fractio =  xf xw / xf xe is not explicitly included in Eq. (43 ), but it is implicitly taken into account in the determination of n in the boundary PAGE 128 128 condition treatment. Thus it is expected that the heat transfer evaluated from Eq. ( 4 3 ) will converge with the same order of accuracy as that of the corresponding boundary condition treatment. This will be verified in the Section 4.3 with numerical tests. 4.2. 3 Temperature Gradient in the Interior Field In athermal LBE simulations, Mei et al (2002) gave a formula for deviatoric stress evaluation as ij ()11 (1)(,)( ) 2neqxee (4 4 ) where ()fneq is the nonequilibrium component of the velocity distribution function f ( x , t ), e and e are the Cartesian components of the discrete velocity vector e and D is the spatial dimensions. Similarly, with the definition of the nonequilibrium part of the temperature distribution function () ()gggneq eq the following formula for the temperature gradient inside the field can be derived mmijij j ij D eg egx () ()112neq neq (4 5 ) The physical basis for the shear stress or temperature gradient in terms of respective ()fneq or ()gneq is that only the nonequilibrium component s of the distribution functions contribute to the velocity or temperature gradient evaluations. Yoshida and Nagaoka (2010) also provided a for mula for the temperature gradient approximation in terms of the total temperature distribution g as 2 01 ()m i ij j x (4 6 ) PAGE 129 129 The formula in Eq. (46 ) was also used in Chapter 3. Recalling the relationship between ij and Dij in Eq. (3 14), Eq. (4 6 ) can be rewritten as mmijij j ij D egv egvx 002 (4 7 ) With the definition of g()eq in Eq. (3 7), m egv(),,0eq can be derived and Eq. (4 7 ) becomes the same as Eq. (4 5 ). Thus the present interior temperature gradient evaluation (4 5) is equivalent to that proposed by Yoshida and Nagaoka (2010) in Eq. (4 6 ). 4.3 Numerical Results and Discussion A series of numerical tests are conducted to validate the applicability and accuracy of the pr oposed heat transfer evaluation method. The first test is for thermal flows in a 2D channel with straight walls. The boundary walls in this test are placed in tran sfer evaluation is fully investigated. Interior temperature gradient evaluation is also examined in this test. The second test is for the transient heat conduction in an inclined semi Curved boundaries are considered in the third to fifth tests. A transient Dirichlet condition is imposed on the surface of a circle in the third test, and thermal convection diffusion in a 3 D circular pipe is studied in the fourth test. Analytical solut ions are available for tests (i) to (iv), and thus the accuracy of the present heat transfer evaluation can be fully assessed. The fifth test involves the natural convection simulation in a square enclosure with a circular cylinder placed in the center and is intended to demonstrate the PAGE 130 130 applicability of the proposed heat transfer evaluation method when the TLBE model is coupled with an athermal LBE model for the velocity field. For simplicity, thermal diffusion is assumed to be isotropic in all cases, and t hus the off diagonal components in the relaxation matrix are zeros and the diagonal elements are related to the diffusion coefficient Dt as 21 2()Dt (4 8 ) A typical value of D = 0.75 is used in each case; t he other relaxation coefficients are P = 1 ( p = 0, 4, 5, 6 in 3D and p = 0, 3, 4 in 2D) since they do not contribute to the leading order solution of the CDE. More details about the numerical stability and accuracy of th e temperature field computation s for those examples have been presented in Chapter 3. 4.3.1 Two D Steady State Channel Flow First we consider the 2 D steady state plug flow as discussed in Section 3.4.1. T he boundary wall s are straight and perpendicular to the discrete velocity vector e (see Fig. 34), thus the boundary heat flux in the discrete velocity direction is equal to the normal heat flux, i.e., n and the surface area is exact Hence we only need t o examine the heat flux on the boundaries for heat transfer evaluation with Dirichlet conditions. In addition, the interior temperature gradients in both xand ydirections are separately examined when Dirichlet and Neumann boundary conditions are respec tively imposed. It is worth mentioning that the present heat transfer evaluation method in Eq. ( 4 3 ) can be directly applied to determine the heat transfer over any portion of a complete boundary wall. Hence even though the overall heat transfer on the boundary walls from x = 0 to x = L in this test is zero due to the sinusoidal boundary conditions PAGE 131 131 applied to obtain analytical solutions, the accuracy assessment of the boundary heat flux on the whole boundary is representative for the heat transfer evaluatio n on any portion of the boundary wall and for any other boundary conditions applied as well. The following results are presented at the Pclet number, Pe = HU/ Dt = 50. 4.3. 1 .1 Channel flow with a Dirichlet boundary condition The Dirichlet problem in Section 3.4.1.1 is further studied to evaluate the heat transfer evaluation. With the analytical solution of the temperature field available in Eq. ( 3 55), the exact boundary heat fluxes on the walls and the interior temperature gradients in both directions can be calculated and the following L2 error norms are computed to assess the accuracy of the LBE solutions E2_qw = n 1/222LBE ex exboundarynodes boundarynodes (4 9 ) E2_gx = TTT xxx1/2 2 2/interiornodes interiornodes LBE ex ex (4 10) E2_gy = TTTyyy1/222/interiornodes interiornodesLBE ex ex (4 1 1 ) The results of E2_qw versus the grid resolution, 1/H, at different values are shown in Figs. 4 2 (a, b). The Schemes 1, 2 and 3 throughout the simulation results in this work correspond to the three different choices of the adjustable coefficient cd 1 that are specified in Eq s. ( 3 35 a c) and used in Eqs. ( 3 32) and ( 4 1 ), respectively, for the Dirichlet condition implementation and boundary heat flux evaluation. Clearly, quadratic convergence is observed for all cases, demonstrating the secondorder accuracy for the boundary heat flux and thus the overall heat transfer evaluation. Fi g. 4 3 shows the PAGE 132 132 error norm E2_qw as a function of the link fraction for the three schemes at Ny = 34 ( Ny denotes the number of lattice nodes in ydirection in Fig. 34 ). For both D = 0.75 and 1.0, it is observed that the variations of the error norm E2_qw are very small for Scheme 3 however, Scheme 2 shows larger errors than the other two. It is thus more favorable to place the boundary at the locations of accuracy of the heat transfer. For curved boundaries, the values are irregularly distributed thus the three schemes are expected to have similar overall performance. The computational results of the enti re temperature field in Section 3.4.1.1 showed that Scheme 2 gave smaller errors than Schemes 1 and 3 in most cases studied thus Scheme 2 is also recommended in the proposed heat transfer evaluation schemes. Furthermore, the dependence of E2_qw on the relaxation parameter D is illustrated in Fig. 4 4 where Ny = 66 is used and the results are presented with respect to ( D 0.5) since the positivity of the diffusion coefficient and numerical stability require that D > 1/2 (Yoshida and Nagaoka 2010, Li et al. 2013a) All t hree schemes approach an asymptotic minimum E2_qw when D /2. It is noted that smaller D values correspond to smaller time steps; thus more time steps are required in general to reach steady state. On the other hand, the error norms increase almost lin early with D for D > 0.60 as observed in Fig. 4 4 For accurate computations it is thus preferable to choose D < 1.0. The interior temperature gradients in both xand ydirections are evaluated with the formula in Eq. ( 4 5 ) and their error norms E2_gx and E2_gy defined in Eqs. ( 4 10) and PAGE 133 133 ( 4 11), versus the grid resolution, 1/ H at different values are shown in Figs. 4 5 and 4 6 respectively, with the three Dirichlet schemes in Eqs. ( 3 35 a c) used. Secondorder accuracy is obtained for each case in Figs. 4 5 and 4 6 and Scheme 2 shows better performance than the other two for most cases studied. 4.3. 1 2 Channel flow with a Neumann boundary condition The Neumann problem in Section 3.4.1.2 is again used here to study the interior temperature gradient e rrors E2_gx and E2_gy defined in ( 4 10) and ( 4 11 ) respectively. Fig s. 4 7 (a, b) show the results of E2_gx and E2_gy versus 1/H for ranging from 0.01 to 0.99. Again, secondorder accuracy is obtained. It is thus concluded that the proposed boundary heat transfer and interior temperature gradient evaluation methods are secondorder accurate on straight b oundaries perpendicular to one of the lattice vectors. 4.3.2 OneD Transient Conduction in an Inclined Semi Infinite Solid The layout of the lattice around a semi infinite solid with different inclination angles is shown in Fig. 316. The width of the solid end is W and identical boundary conditions are given along the solid end for this oneD conduction problem. With zero initial temperature and Dirichlet condition T ( l = 0) = f ( t ) at the end, the transient temperature is solved to be 22222(,) ()4tltDtlTtlftedDex (4 12) where l is the distance in the normal direction of the solid end, and the heat transfer rate at the boundary is 0(/)wntlQWqWDTnex (4 13) PAGE 134 134 Two transient boundary conditions for f ( t ) are considered in this section. The three Dirichlet schemes in Eqs. ( 3 35 a c) give very similar numerical solutions; thus only the results from Scheme 2 are presented for brevity. 4.3.2.1 Dirichlet Boundary condition f ( t ) = kt1/2 This temperature boundary condition with k being a constant corresponds to a constant heat flux on the solid end (Hahn and Ozisik, 2012) 02ntltTkqDnex (4 14) In the T LBE simulations, the boundary heat fluxes n in the discrete velocity directions obtained from Eq. (41) are summed up along the boundary to yiel d the heat transfer rate wQ as in Eq. (43 ). The numerical results of wQ are compared with exact solutions in Fig. 4 8 for different inclination angles = tan1(1.0), tan1(1.2), tan1(1.5), and tan1(2.0) with k = 1. Because of the initial temperature discontinuity at t = 0 on the boundary, there exists oscillation in the numerical results of wQ at small time in Fig. 4 8 for eac h case tested. After 10 time steps, the numerical results agree very well with exact solutions. To further verify the heat transfer evaluation in this test, another Dirichlet boundary condition resulting in timevarying heat transfer is considered in the f ollowing. 4.3.2.2 Dirichlet Boundary condition f ( t ) = 1 exp( ) When an exponentially decaying boundary condition f ( t ) = 1 exp( ) is applied, the heat flux on the boundary becomes (Hahn and Ozisik, 2012) 1/2 1/20 /[1 ]ntlt tTqDDtD nexerfi() (4 15) where is the Gamma function and erfi is th e imaginary error function. PAGE 135 135 The results of the heat transfer rate ()wQt obtained from Eq. (43 ) are compared with the exact solutions in Fig. 4 9 for different inclination angles at = 0.1. Again, limited oscillations of the simulated wQ values are observed at small time ( N < 10) in Fig. 4 9 due to the discontinuity of the initial condition. After that very good agreement between the numerical and exact solutions is noticed in Fig. 4 9 for all cases tested. Th is test demonstrates that when the straight boundary is not perpendicular to t he discrete velocity vectors of the TLBE model and intersected by the square lattice 4 3 ), is still able to provide accurate results. 4.3.3 Two D Tr ansient H eat C onduction I nside a C ircle Following the example in Section 3.4.4, we consider here the transient sinusoidal boundary condition applied on the curved boundary of a circle with radius r0. The analytical solution for this 2D problem is given by Eq. (370), assuming no tempera ture variation in the azimuthal direction. Two error norms are defined to assess the present heat transfer evaluation E ( t*) = **0()()/wwQtQtrLBE ex and (4 16) E2 = 2/ 1** ***000011()()()()2ww wwQtQtdtQtQtdt LBE ex LBE ex ,(4 17) where t = /2 To eliminate the effect of the initial variations on the average error E2 defined in Eq. ( 4 17), the computation is performed for a sufficiently long period of time ( t > 1). The variations of E ( t *) in Eq. ( 4 16) at different grid resolution, 1/ r0, are plotted in Fig. 4 10 when the coefficient cd1 determined from Scheme 2 in Eq. ( 3 35b) is used in Eq. ( 4 1 ). After the dynamic steady state is reached the average error E2 is computed PAGE 136 136 and the results are shown in Fig. 4 11 for each of the three schemes implemented. The Stokes number St is set to be unity in all of the present simulations. On average, the decrease of the error norm E2 in Fig. 4 11 is linearly proportional to the grid resolution, 1/ r0, as shown on the log log scale for each of the three sc hemes implemented. The results demonstrate that the present heat transfer evaluation is first order accurate for curvedboundary situations. Compared with the quadratic convergence of the boundary heat flux and heat transfer for the 2D channel flow with straight walls, the degradation of the convergence order in this test is due to the first order accurate heat flux approximation on the curved boundary. It has been verified with numerical tests in Chapter 3 that the boundary heat flux obtained from Eq. ( 4 1 ) in TLBE simulations is first order accurate for Dirichlet conditions on curved boundaries because of the irregularly distributed lattice fractions intersected by the curved boundary, even though the temperature field obtained from the TLBE results is second order accurate. Thus it is concluded that the proposed heat transfer evaluation method has the same first order accuracy as the heat flux computation on curved boundaries. It is noticed that in this test the boundary heat flux in the normal direction and the surface area on the curved boundary are not required in the determination of the overall heat transfer. It is thus believed that the proposed heat transfer evaluation method is very efficient in curvedboundary situations such as thermal flows in porous media. 4.3.4 ThreeD S teady S tate C ircular P ipe F low with a Dirichlet C ondition This example simulates the thermal convection and diffusion inside a 3D circular pipe for which the temperature field has been studied in Section 3.4.5.1 for boundary co ndition implementation analysis Applying the same sinusoidal Dirichlet boundary condition the analytical solution for the temperature field is given in Eq. (375). The heat PAGE 137 137 transfer rate at the curved surface of the circular pipe based on Eq. (43 ) is com puted and the following error norm is defined to assess the present heat transfer evaluation 1/2 2 2 2()ww z w zQQ E Q LBEex Qw ex ( 4 18) Note that the boundary heat transfer in Eq. ( 4 18) is evaluated at each location along the axial z direction since the over all heat transfer over the entire surface of the circular pipe is zero. Figure 412 shows the results of 2EQw versus the grid resolution, 1/ r0, with the three boundary schemes given in Eqs. (335 a c) implemented. The characteristic Pclet number is Pe = (2 r0) U/Dt = 5 0 in all simulations Linear convergence for the heat transfer is observed in Fig. 4 12 for all three schemes; and Scheme 2 gives smaller errors than the other two. The first order accuracy of the present heat transfer ev aluation method is thus verified for 3D therm al flows with curved boundaries 4.3.4 Two D Steady State Natural Convection in a Square Enclosure with a Circular Cylinder in the Center Natural convection between an outer square enclosure with length L and an inner circular cylinder of radius r0 has been investigated in a number of previous studies (Shu and Zhu 2002, Peng et al. 2003, Kim et al. 2008, Jeong et al. 2010, Moukalled and Acharya 1996). The present study considers the basic configuration that the cylinder is located in the center of the square enclosure with r0 = 0.2 L and all the boundary walls are stationary. The straight walls of the enclosure are maintained at a constant low temperature of Tc, whereas the inner circular wall is at a constant hig h temperature Th (see Fig. 4 13). The fluid properties are assumed to be constant except that the density PAGE 138 138 in the buoyancy term is based on the Boussinesq approximation. The Prandtl number is set to be 0.71, corresponding to that of air, and the Rayleigh number, Ra varies in the range of 103 6 in present simulations. To simulate the velocity field, the 2D MRT LBE model by Lallemand and Luo (2000) = 0.5 in both directi on as shown in Fig. 4 13) so that the standard bounceback condition is used. For the inner circular boundary, the secondorder accurate noslip velocity condition treatment by Mei et al. (1999) is implemented. The LBE model for the velocity field and the TLBE model for the temperature field are coupled in the simulation. The convergence is considered to be obtained when the relative error between 50 successive lattice time steps is less than 108 for both the velocity and temperature fields. Figures 4 14 a nd 4 15 show the isotherms and streamlines, respectively, at Ra = 103, 104, 105 and 106 (the lattice grids 309 309 are used for all cases). The present results agree well qualitatively with those reported in (Shu and Zhu 2002, Peng et al. 2003, Kim et al. 2008, Jeong et al. 2010, Moukalled and Acharya 1996) To quantify the accuracy of the present simulation results, the surfaceaveraged Nusselt numbers on both the inner circular wall and the outer straight walls are examined. At s teady state, the local heat transfer rate on the boundary walls of the inner cylinder and the outer enclosure can be obtained from hcfT qhTTk nwall (4 19) PAGE 139 139 where h is the local heat transfer coefficient, kf is the thermal conductivity of the fluid, and n is the normal of the boundary walls pointing outward of the enclosure. From Eq. ( 4 19) one can determine h as ffThk kTnn1wall wall (4 20) where hcTTT and is the dimensionless temperature c ()/ The average Nusselt numbers over the inner and outer boundaries are then defined as iiiifihSNu Skn and oooofohSNu Skn (4 21) Here Si denotes half of the circumferential length of the inner circle, and So is half of the total length of the outer straight walls due to symmetry. At steady state, the average Nusselt numbers iNu and oNu should be equal to each other. Based on the present heat transfer evaluation method, the heat transfer on the inner circular surface is obtained by integrating the boundary heat flux n in the discrete velocity directions of the TLBE model with unit effective surface area. T hus the average Nusselt number iNu on the inner surface in Eq. ( 4 21 ) is obtained from i ini niiifffiNuSSnkk k211122 (4 22a) where n i is determined using Eq. (41 ), and w iQ is the total heat transfer rate on the entire circular boundary that is computed using Eq. ( 4 3 ). For the outer straight boundaries perpendicular to the discrete velocity vectors in xor ydirection, the normal heat flux is the same as the corresponding heat flux in the dis crete velocity direction, PAGE 140 140 i.e., nn and the surface area is exact as the outer boundaries are placed in the center of the lattice links. T hus the average Nusselt number oNu on the outer surface is o ono nooofffoNuSSnkk k211122 (4 22b) It should be noted no in the above is determined using Eq. ( 4 1 woQ is the heat flow rate on the entire outer surface. A grid refinement study is carried out first, an d the average Nusselt numbers iNu and oNu obtained from Eqs. ( 4 22a) and ( 4 22 b), respectively, are presented in Table 4 1 at different lattice grid sizes. The results can be considered well converged at the mesh size 359 359. Table 4 2 compares the present converged results with the averag e Nusselt number iNu evaluated on the inner surface in (Shu and Zhu 2002, Peng et al. 2003, Kim et al. 2008, Jeong et al. 2010, Moukalled and Acharya 1996) The specific methods applied in those references are also listed at the bot tom of Table 4 2. The present results in Tables 4 1 and 4 2 are obtained using Scheme 2 for both the boundary implementation and the heat flux evaluation. The present results for iNu and oNu in Tables 4 1 and 4 2 are very close to each other for all Ra values simulated, indicating that the present heat transfer evaluation method on curved boundaries is accurate compared with that obtained from straight boundaries. In addition, the present results are in good agreement with previous numerical results reported in (Shu and Zhu 2002, Peng et al. 2003, Kim et al. 2008, Jeong et al. 2010, Moukalled and Acharya PAGE 141 141 1996) It is thus verified that the present heat transfer evaluation method is accurate and efficient in thermal hydrodynamic coupling simulations involving curved geometry 4.4 Summary and Conclusions An efficient and accurate approach for evaluating local and overall heat transfer rates on curved boundaries is proposed in the thermal lattice Boltzmann equation method. The boundary heat flux in the discrete velocity directions is obtained with the temperature distribution functions at the lattice nodes close to the boundary and the given Dirichlet or mixed boundary condition. As the lattice spacing is uniform i n the Cartesian coordinate system, the effective surface area for any of the discrete heat fluxes on the boundary is constant. Based on energy conservation, the integration of the heat fluxes in the discrete velocity directions with constant effective surf ace area leads to the total heat transfer rate. There are two major benefits for the proposed heat transfer evaluation: no requirement for the normal heat flux determination and no surface area approximation. Thus the proposed heat transfer evaluation method is very efficient and convenient to implement compared with other methods that approximate the boundary heat transfer based on the simulated temperature field. Numerical tests have demonstrated that the heat transfer evaluation method is secondorder ac curate for straight walls perpendicular to one of the discrete lattice velocity vectors and first order accurate for curved boundaries PAGE 142 142 Table 41. Surfaceaveraged Nusselt numbers at different grid resolution. Average Nusselt numbers iNu oNu lattice grids Ra = 10 3 Ra = 10 4 Ra = 10 5 Ra = 10 6 159 159 3.170, 3.168 3.229, 3.226 4.932, 4.924 9.313, 9.306 209 209 3.170, 3.167 3.229, 3.225 4.925, 4.918 9.121 9.117 259 259 3.170, 3.167 3.228, 3.225 4.921, 4.915 9.041, 9.037 309 309 3.169, 3.167 3.227, 3.226 4.917, 4.913 8.997, 8.992 359 359 3.169, 3.167 3.227, 3.226 4. 916, 4.912 8.971, 8.966 Table 42. Comparison of the present s urface averaged Nusselt numbers with previous results. Average Nusselt number Nu Ra Present iNu Present oNu Shu & Z hu+ Peng et al. Kim et al. Jeong et al. Moukalled & Acharya 10 3 3.169 3.167 3.396 3.399 10 4 3.227 3.226 3.24 3.24 3.414 3.412 3.331 10 5 4.916 4.912 4.86 4.84 5.138 5.176 5.08 10 6 8.971 8.966 8.90 8.75 9.39 9.171 9.374 + Differential quad rature method with vorticity stream function formulation (Shu and Zhu 2002) ; Taylor series expansion and least squares based lattice Boltzmann method (Peng et al. 2003) ; Immersed boundary finite volume method (Kim et al. 2008) ; Immersed boundary thermal lattice Boltzmann method (Jeong et al. 2010) ; Finite volume method (Moukalled and Acharya 1996) PAGE 143 143 exterior xexfxw field boundary wall n nnn y x eexff x Figure 4 1 Illustration of the heat transfer evaluation on a curved boundary (closed circles: field nodes; closed squares: boundary nodes; open circ les: exterior nodes). (a) (b) Figure 4 2 Wall heat flux error E2_qw versus the grid resolution, 1/ H for the channel flow Dirichlet problem at (a) = 0.50, 0.25, and 0.75, and (b) = 0.50, 0.01, and 0.99. PAGE 144 144 Figure 4 3 Wall heat flux error E2_qw versus the lattice link fraction at Ny = 34 for the channel flow Dirichlet problem. Figure 4 4 Wall heat flux error E2_qw versus ( D 0.5) at N y = 66 for the channel flow Dirichlet problem. PAGE 145 145 (a) (b) Figure 4 5 Interior gradient error E2_gx (in xdirection) versus the grid resolution, 1/ H for the channel flow Dirichlet problem at (a) = 0.50, 0.25, and 0.75, and (b) = 0.50, 0.01, and 0.99. (a) (b) Figure 4 6 Interior gradient error E2_gy (in ydirection) versus the grid resolution, 1/ H for the channel flow Dirichlet problem at (a) = 0.50, 0.25, and 0.75, and (b) = 0.50, 0.01, and 0.99. PAGE 146 146 (a) (b) Figure 4 7 Interior gradient errors of (a) E2_gx, and (b) E2_gy, respectively, versus the grid resolution, 1/ H for the channel flow Neumann problem. Figure 4 8 Variations of the heat transfer rate with t ime for the Dirichlet condition T ( l =0)= t1/2 at the end of an inclined semi infinite solid (symbols: LBE results; lines: exact solutions). PAGE 147 147 Figure 4 9 Variations of the heat transfer rate with t ime for the Dirichlet condition T ( l = 0) = 1 exp( 0.1 t ) at the end of an inclined semi infinite solid (symbols: LBE results; lines: exact solutions). Figure 4 10. Variations of t he error E ( t*) defined in Eq. (41 6 ) with time at different grid resolution for the transient heat conduction inside a circle. PAGE 148 148 Figure 4 11. Heat transfer error, E2, defined in Eq. (41 7 ) versu s the grid resolution, 1/ r0, for the transient heat conduction inside a circle. Figure 4 12. Heat transfer error 2wQE_ versus the grid resolution, 1/ r0, for the circular pipe flow with a Dirichlet boundary condition. PAGE 149 149 r0 field g L L Th, Tc, u= 0,v= 0 u= 0 ,v= 0 = 0.5 = 0.5 Figure 4 13. Schematic depiction of the computational domain ( r0 = 0.2 L ) and the Dirichlet thermal and velocity boundary conditions on the inner and outer walls (a) (b) (c) (d) Figure 4 14. Isotherms in the field between the circular cylinder and the c avity walls at (a) Ra = 103, (b) Ra = 104, (c) Ra = 105, and (d) Ra = 106. (a) (b) (c) (d) Figure 4 15. Streamlines in the field between the circular cylinder and the cavity walls at (a) Ra = 103, (b) Ra = 104, (c) Ra = 105, and (d) Ra = 106. PAGE 150 * This work has been submitted to International Journal of Heat and Mass Transfer (Li et al. 2013c). 150 CHAPTER 5 MULTIPLE RELAXATION TIME LATTICE BOLTZMANN MODEL FOR THE AXISYMMETRI C CONVECTION DIFFUSION EQUATION* 5 .1 Introduction Heat and mass convectiondiffusion in axisymmetric systems such as circular cylinders and annular cavities formed by coaxial cy linders are widely encountered in engineering practices. Examples of applications include heat exchangers, desalination systems, solar energy collectors, electronic cooling devices and chemical reactors. Fluid flows and heat and mass transfer in cylindric al systems are essentially threedimensional (3D) problems, while under axisymmetric conditions they can be reduced to quasi 2 D problems, which greatly reduces the computational effort. The lattice Boltzmann equation (LBE) method has become an effective numerical method for the convectiondiffusion equation (CDE) (Flekky 1993, van der Sman and Ernst 2000, ServanCamas and Tsai 2008, Shi and Guo 2009, Yoshida and Nagaoka 2010, Ginzburg et al. 2010, Ginzburg 2013, Li et al. 2013a, Li et al. 2013b) Especial ly, the multiplerelaxationtime (MRT) lattice Boltzmann (LB) model (Yoshida and Nagaoka 2010, Li et al. 2013a, Li et al. 2013b) is robust and retains the secondorder accuracy for curved boundaries. Although 3D LB models can be directly applied to simulate the axisymmetric CDE when the curved boundary condition is properly handled, quasi 2 D LB models, which require far less computational effort and bypass the curved boundary treatment for straight pipes, are attractive alternatives for solving the axisym metric CDE as demonstrated by the LB models developed for axisymmetric athermal flows (Halliday et al. 2001, Lee et al. 2006, Reis and Phillips 2007, Reis and Phillips 2008, Zhou 2008, Chen et al. 2008, Huang and Lu 2009, Guo et al. 2009, Wang et al. 2010, Li et al. 2010, PAGE 151 151 Zhou 2011) and thermal problems (Peng et al. 2003, Huang et al. 2007, Chen et al. 2009, Li et al. 2009, Zheng et al. 2010a, Zheng et al. 2010b) The standard LB models break the computation into a collision step and a streaming step. In order for axisymmetric LB models to work using the existing LB models developed in 2 D Cartesian coordinate systems in terms of the standard collision streaming procedure, either the governing equations for momentum and energy transport in the cylindrical c oordinate system are represented in a pseudoCartesian coordinate system, or the effect of axisymmetry is modeled at the particle distribution function level (Guo et al. 2009, Wang et al. 2010, Zheng et al. 2010a, Zheng et al. 2010b) The first approach results in additional terms, including the velocity and temperature gradients, in the transformed governing equations, and they can be treated as additional source terms in the respective macroscopic transport equations. This treatment is straightforward since it is carried out at the macroscopic level. For example, the typical axisymmetric CDE is r z rr zzruurD DGtrrzrrrzz (5 1) where G is the typical source term representing possibly the effects of viscous dissipation, pressure work, volumetric heating, heat released or mass created from chemical reaction for the energy or species transport equation. In Eq. ( 5 1), ur and uz are the velocity components in the radial r and axial z directions, respectively, Dr r and Dzz are the diagonal diffusion coefficient s of the diffusion tensor Dij. It can be rewritten as PAGE 152 152 r z rr zz aauuDDSSG trzrrzz12 with r rraauDSSrrr12, (5 2) Equation ( 5 2) becomes an ordinary CDE in the pseudoCartesian coordinate system ( r z ) with three source terms, Sa 1, Sa 2, an d G of different origins. In the second approach, new LB evolution equations for the particle distribution functions must be carefully constructed so that the correct macroscopic transport, Eq. ( 5 1), can be recovered through the ChapmanEnskog expansion for the new evolution equation. When MRT LB models such as the standard ones in (Lallemand and Luo 2000, Yoshida and Nagaoka 2010) are desired, such modification at the microscopic level for the evolution equations of the distribution functions may involve redefining the transformation matrix and the equilibrium distribution functions. The first approach is thus preferred if the gradients can be efficiently obtained from the distribution functions instead of using finitedifference schemes. A review of the various LB models for axisymmetric fluid momentum equations can be found in the recent work of Zhou (2011) which also introduced an improved axisymmetric LB model based on his original version in (Zhou 2008) The main advantages of Zhous revised model co mpared with the other axisymmetric LB models are as follows: (i) the source terms represent the same additional terms in the transformed governing equations and no extra source terms are needed in the evolution equation, (ii) it preserves the original computational procedure for both the distribution functions and the macroscopic variables including density and velocities as those in the conventional LBE method, and (iii) no velocity gradients are required. It is well known PAGE 153 153 that the MRT LB models (Lallemand and Luo 2000, Yu et al. 2003) have better numerical stability and accuracy than BGK LB models, i.e., the singlerelaxationtime (SRT) collision models based on the work of Bhatnagar, Gross and Krook (1954) Thus when the fluid momentum equations are invol ved, an MRT version of Zhous axisymmetric LB model (Zhou 2011) is used in the present work. Several axisymmetric thermal LB models, which simulate the axisymmetric CDE for temperature, have also been proposed in the literature (Peng et al. 2003, Huang et al. 2007, Chen et al. 2009, Li et al. 2009, Zheng et al. 2010a, Zheng et al. 2010b) Peng et al. (2003) and Huang et al. (2007) proposed hybrid approaches in which the fluid velocities in the radial and axial directions are simulated using the athermal axi symmetric LB model by Halliday et al. (2001) while the azimuthal velocity and the temperature field are solved using finitedifference schemes. Chen et al. (2009) proposed an LB model for axisymmtric buoyancy driven thermal flows in which the vorticity st ream function formulation was used for the velocity field so that the governing equations for vorticity and temperature became CDEs with source terms. The vorticity and temperature gradients in the source terms were obtained from finitedifference calculat ions. In addition, a Poisson equation for the stream function must be solved during each time step in ( Chen et al. 2009) To avoid the finitedifference computation of the temperature gradient, i.e., Sa2 in Eq. (2), in LB models, Li et al. (2009) proposed an implicit evolution equation for the distribution functions to recover the axisymmetric energy equation. Additional body forces can also be modeled by including a forcing term in this model. The implicitness is then eliminated by introducing a new distri bution function, which includes not only the original distribution function but also the source PAGE 154 154 term. As a result, the modified evolution equation is more difficult to implement than that in standard LB models; Furthermore, the temperature calculation in ( Li et al. 2009) requires the radial velocity component and thus it deviates from the standard calculation of macroscopic variables in the LBE method. Zheng et al. (2010a) extended their athermal axisymmetric LB model (Guo et al. 2009) to simulate thermal f lows by including a D2Q4 LB model for the temperature field. The radial coordinate and the temperature distribution functions were lumped together and the source term in the evolution equation of the temperature distribution functions was properly construc ted so that the scalar source term, Sa 1, and the gradient term, Sa 2 given in Eq. ( 5 2), do not show up explicitly in the thermal LB model. It was realized that, however, due to its special construction of the source term, the extension of Zheng et al.s mo del (2010a) to include additional source terms, such as viscous dissipation and pressure work, is not straightforward. To overcome this limitation, Zheng et al. (2010b) also proposed a coupled LB model for axisymmetric thermal flows with viscous dissipatio n and pressure work starting from the continuous axisymmetric Boltzmann equation. It is noticed that the velocity distribution functions are explicitly coupled in the evolution equation of the energy distribution functions in (Zheng et al. 2010b) It is em phasized that all the existing LB models for axisymmetric thermal flows in (Chen et al. 2009, Li et al. 2009, Zheng et al. 2010a, Zheng et al. 2010b) used the BGK collision operator. A robust and straightforward MRT LB model that can simulate the general axisymmetric CDE is highly desired. Furthermore, we notice that the numerical validation regarding the order of accuracy of the axisymmetric thermal LB models has not been well established. PAGE 155 155 In this work, we propose an MRT LB model for the general axisymmetr ic CDE for scalar transport such as temperature, mass concentration and azimuthal velocity component in rotational flows. When the CDE is coupled with the hydrodynamic equations, the radial axial velocity field is solved using Zhous axisymmetric LB model (Zhou et al. 2011) with the BGK collision model replaced by an MRT based model. The axisymmetric CDE in the cylindrical coordinate system is transformed to an ordinary CDE in the Cartesian coordinate system, and the extra terms due to the coordinate transf ormation are considered as additional sources. The scalar source term, Sa 1 in Eq. ( 5 2), is obtained from the summation of the distribution functions, and the gradient term, Sa 2 in Eq. ( 5 2), is calculated from the nonequilibrium parts of the distribution functions; thus no finite difference calculations are required. The transformed CDE in the Cartesian coordinate system is then simulated with the D2Q5 MRT LB model proposed in (Yoshida and Nagaoka 2010) Compared with the existing LB models for axisymmetr ic thermal flows in (Chen et al. 2009, Li et al. 2009, Zheng et al. 2010a, Zheng et al. 2010b) the proposed model for the axisymmetric CDE has the following features: (i) the macroscopic variable and its gradient terms in the CDE, such as Sa 1 and Sa 2 in E q. ( 5 2), are naturally treated as additional source terms, and they are obtained from the distribution functions; thus no special treatment for the LB evolution equation or the equilibrium distribution functions is required; (ii) it is easier to implement when compared with Li et al.s model (2009) and the standard calculation of the macroscopic variables is preserved; (iii) the forms of the source terms can be very general and they can be directly included in the present model, (iv) viscous dissipation and pressure work can be conveniently treated as source terms in the present model; PAGE 156 156 thus the explicit coupling of the velocity distribution function in the evolution equation for the energy distribution function in (Zheng 2010b) is avoided, and (iv) rather t han using BGK collision operators as in (Chen et al. 2009, Li et al. 2009, Zheng et al. 2010a, Zheng et al. 2010b) an MRT collision operator is applied to the present model. To validate the applicability and accuracy of the proposed LB model for the axisymmetric CDE, several numerical tests are conducted in this study. Analytical solutions are available for the first four tests, and the results show that the proposed model possesses secondorder accuracy. The axisymmetric LB model for the radial axial velo city field by Zhou (2011) and the present LB model for the temperature CDE are coupled in the fifth and sixth tests, where natural convection in an annulus between two coaxial vertical cylinders and confined swirling flows in a vertical cylindrical container are studied, respectively. The numerical results from these two tests are compared with published data, and good agreement is obtained for each case. The rest of this chapter is organized as follows. In Sec. 5. 2, the D2Q5 MRT LB model for the axisymmetr ic CDE is proposed ; in which the treatment of the macroscopic scalar variable and its gradient in the source terms in the transformed CDE is provided. Numerical validation with six tests is presented in Sec. 5. 3, and Sec. 5.4 concludes the chapter 5 2 Lat tice Boltzmann Model for the Axisymmetric ConvectionDiffusion Equations The treatment of the additional terms, Sa 1, and Sa 2, in Eq. ( 5 2) due to the polar to Cartesian coordinate transformation as additional source terms is straightforward in the LBE method. While it is preferable to obtain these macroscopic quantities, especially the gradient term Sa 2, from the microscopic distribution functions rather than PAGE 157 157 using finitedifference schemes on the solved macroscopic field since the standard collision stream ing procedure is carried out on the distribution function level. As discussed in Sec. 4.2.3, Yoshida and Nagaoka (2010) gave a formula in their MRT LB model to obtain the scalar gradient based on the distribution functions. The accuracy of the gradient for mula was examined and verified in Chapter 2 with numerous computational examples. In Sec. 4.2.3, an evaluation formula for the gradient based on the moment of the local nonequilibrium components of the distribution functions is also proposed and it is ver ified that these two formulas are essentially equivalent. The nonequilibrium approach in Sec. 4.2.3 for the scalar gradient evaluation will be used in this chapter In the LBE method, t he LBE governs the evolution of the scalar distribution function, g ( x , t ), which is defined in the discrete phase space ( x ) and time t where x = r er + z ez is the spatial vector and is the particle velocity usually discretized to a small set of discrete velocities {  = 1, 2, m }. With the objective to recover the CDE ( 5 2), we propose the following evolution equation for the distribution functions g12,,Lxe xg ( 5 3 ) where gtgt ,,, xx e is the th discrete velocity vector (the discrete velocity set { e} for the twodimensional fi ve velocity (D2Q5) LB model is shown in Fig. 5 1), is the time step, L is the standard collision operator, and is the weight coefficient. The leading order solution of the CDE is obtained from m g0 ( 5 4 ) and the addition al source terms are obtained from the distribution functions as PAGE 158 158 mmrrrauuuSggrrreq100 ( 5 5 ) mm rr a rr rrD S eg egg rrr neq eq 2 1111 11 11 22 ( 5 6 ) where geq ( see definition in Eq. ( 3 7 )) and gneq are the equilibrium and non equilibrium components of the distribution function g, respectively, with gggneq eq rr is a relaxation coefficient (see Eq. (3 26)) and e is the projection of the discrete velocity e in the radial direction. The derivation and numerical validation of the gradient evaluation in Eq. ( 5 6) in terms of gneq can be found in Sec. 4.2.3. The D2Q5 LB model with a MRT collision operator (Yoshida and Nagaoka 2010) is employed to simulate the LBE (5 3). The two step procedure is summarized as collision step: gtgttteq112,,MS,,xxmxmx and ( 5 7 ) streaming step: g ,, xe x ( 5 8 ) where g is again the post collision distribution function. The details of the transformation matrix and the multiple relaxation coefficients are give n in Sec. 3.3.2. When the CDE (2) is coupled with the momentum transport equation, the axisymmetric LB model by Zhou (2011) is applied for the radial axial velocity field and the BGK collision operator is extended to an MRT version for consistency. The ext ension of Zhous BGK model (2011) to the MRT model is given in Appendix C PAGE 159 159 5 3 Numerical Validation and Discussion To validate the accuracy of the proposed LB model for the axisymmetric CDE, we present six numerical tests in this section. The velocity fiel d is specified in each of the first four tests and the analytical solutions for are available. The accuracy and convergence order of the numerical solutions are thus investigated in detail. The fifth and sixth tests consider axisymmetric thermohydrodynamic problems, where the radial axial velocity fields are simulated with an MRT version of the athermal LB model by Zhou (2011) and the presently proposed LB model is applied for the thermal field as well as the azimuthal velocity field in the sixth test. These two LB models are coupled to reach the steady states and the numerical results are compared with published results. 5. 3.1 Steady convection diffusion in an annulus In this example, we consider the axisymmetric convectiondiffusion in an infinitely lo ng annulus defined in Ri r Ro, where the variations in both the azimuthal and axial directions are neglected. With a constant volumetric flow rate Q supplied on the inner surface, the radial velocity in the annulus is ur = Q / r The convectiondiffusion inside the annulus is c haracterized by the Pclet number, iRirrPeuRDQD // For given Dirichlet conditions ioiorR,,() the exact solution for inside the annulus is PePe PePeio oiioPePeoiRrrRRrRzRRex(,) ( 5 9 ) With no axial variation only five lattice nodes are used in the z direction, and periodic boundary conditions are applied at the ends. Each of the inner and outer boundaries is placed at the center of the lattice PAGE 160 160 Dirichlet bounceback scheme, which is recovered by all three schemes in Eqs. ( 3 35 a c), is implemented. To assess the accuracy of the proposed LB model, the following relative L2 error norm is defined in the annulus rz rzE1/2222LBEex ex, / ( 5 10) Figure 5 2 shows the error norm, E2, versus the grid resolution, 1/ Ri, at Pe = 10, 20, 50 and 100, Ro = 2 Ri, rr = D = 0.75, i = 0.5, and o = 1.0. The error norm E2 increases as Pe becomes larger; and secondord er convergence is noticed for each Pe number tested. Hence the present axisymmtric LB model is secondorder accurate. 5. 3.2 Steady convection diffusion in a circular pipe Following the circular pipe test in Sections 3.4.5 and 4.3.4, the fluid flow inside t he circular pipe is assumed to be a plug flow with constant axial velocity U and zero radial velocity. The axisymmetric CDE at steady state is thus simplified to rrrr zzDUDDzrrzzrr ( 5 11) The last term on the right side of Eq. ( 5 11) is treat ed as a source term in the present model as in Eq. ( 5 6). Only isotropic diffusion is considered for simplicity so that Drr = Dzz = D and the Pclet number is Pe = (2 R0) U / D The computational domain and the square lattice in the r z plane are schematically shown in Fig. 5 3 A periodic boundary condition in the z direction is imposed so that f ( x + L ) = f ( x) is valid for both and the distribution function g. The centerline, r = 0, is always located half way between the / r r =0 = 0 is imposed due to axisymmetry. Dirichlet and Neumann boundary conditions on the surface of r = R0 will PAGE 161 161 be studied separatel y in this section. For the Dirichlet condition on the boundary with 3 35 b) is implemented as it showed smaller relative errors than the other two schemes in Eqs. ( 3 35a) and ( 3 35c ) for most cases studied in Chapter 3. T he Pclet number Pe = 20 and the relaxation coefficient D = 0.75 are used in all simulations. This problem has been studied in Chapters 3 and 4 where the D3Q7 LB model was used in the 3D domain and curved boundary condition treatments were required to pr eserve the exact geometry of the circular pipe. In the present work, in contrast, as the 2 D axisymmetric LB model is used, the effect of the curved boundary can be eliminated. 5. 3.2.1 Dirichlet boundary condition on the wall A sinusoidal Dirichlet boundar y condition is imposed on the boundary walls and the analytical solution for is given in Eq. (375). For illustrative purposes, Fig. 5 4 shows the contours of in the r z plane from the numerical simulation using Nr = 42, Nz = 83 ( Nr and Nz denote the n umbers of lattice nodes in the r and z To assess the accuracy of the present LB model, the error norm E2 defined in Eq. ( 5 10) is computed. In addition, with the analytical solution in Eq. ( 3 75 ), t he exact solutions for / r / z in both directions can be obtained. The numerical boundary flux is obtained from the technique proposed in Chapter 3 in the LBE method (see Eq. (3 41a)) with d n ddc gt gtgt 1 111 22 2 11 xxx PAGE 162 162 and the interior gradient is evaluated using the idea shown in Eq. ( 5 6). The following error norms are further defined to assess their numerical accuracy E2_qw = rR rRrrr001/222LBE ex ex/ ( 5 12) E2_gr = rrr1/2 2 2 interiornodes interiornodes LBE ex ex/ ( 5 13) E2_gz = zzz1/2 2 2 interiornodes interiornodes LBE ex ex/ ( 5 14) The results of E2, E2_qw, E2_gr and E2_gr, defined in Eqs. ( 5 10 5 12 5 14) versus the grid resolution, 1/ R0, 5 5 (a d), respectively. For each case in Fig. 5 5 second order convergence is observed. It should be noted that only first order convergence was obtained for the boundary flux and boundary heat transfer in Figs. 3 28 and 412 in Sections 3.4.5.1 and 4.3.4 respectively, due to the curved boundary effect when the D3Q7 LB model was used. In the present axisymmetric case, the curved boundary effect is eliminated and thus secondorder accuracy is preserved for the interior distribution of the boundary flux, and the interior gradients in both directions. 5. 3.2.2 Neumann boundary condition on the wall When the sinusoidal Dirichlet condition in 5.3.2.1 is replaced with a sinusoidal Neumann condition, the analytical solution for can also be solved and it is given in Eq. (3 79). PAGE 163 163 The contours of in this case are very similar to that with a Dirichlet wall boundary condition except for the magnitude; thus the contours are not shown for brevity. According to Eq. (341b) the boundary values of are computed based on the distribution functions at the adjacent lattice nodes and the boundary flux information, d Dd ddgt gtgt 1 111222 111xxx The following L2 error norm is further defined to assess the LB results E2_tw = rR rR001/222LBEex ex/ ( 5 15) Figures 5 6 (a d) show the error norms E2, E2_tw, E2_gr and E2_gr defined in Eqs. ( 5 10 5 15, 5 13 and 5 14), respectively, versus the grid resolution, 1/ R0 Secondorder accuracy is obtained for all cases in Fig. 5 6 The proposed LB model for the axisymmetric CDE with Neumann boundary conditions is also verified to be secondorder accurate. 5. 3.3 Steady heat conduction inside a sphere When there is no temperature variation in the azimuthal direct ion, the heat conduction inside a sphere is axisymmetric and it can be represented in the radial axial cylindrical coordinate system. Thus the axisymmetric LB model can be applied to simulate the heat conduction in a sphere. The computational domain and th e lattice layout are schematically shown in Fig. 5 7 The curved geometry of the sphere is x y) for those field lattice nodes, such as point P in Fig. 5 7 adjacent to the boundary. After the l PAGE 164 164 along the lattice velocity direction e is obtained for each boundary node, the Dirichlet and Neumann boundary condition treatments in Eqs. ( 3 32) and ( 3 38 ) are implemented, respectively. 5. 3.3.1 Dirichlet boundary co ndition on the wall A Dirichlet boundary condition of a Legendre polynomial Pn( = cos ) on the sphere surface, r = R0, is considered and the analytical solution for the interior temperature is (Hahn and Ozisik 2012) n nr ex 0(,)/() ( 5 16) A representative case of P424()35303/8 is studied in this work. The axis of the sphere is placed in the center of the lattice link (see Fig. 5 7 ) so that the symmetry / r r =0 = 0, which is a Neumann boundary condition, is used in all simulations. Figure 5 8 shows the error norms of the interior temperature field and the boundary heat flux E2 and E2_qw defined in Eqs. ( 5 10) and ( 5 12), respectively, versus the grid resolution 1/ R0, with the Dirichlet Schemes 13 in Eqs. ( 3 35 a c) implemented. It should be noted that the boundary heat flux in this test is evaluated along the discrete lattice velocity directions using the formula proposed in Eq. (341a) For those boundary nodes intersected by the lattice vectors along z direct ion, rRr0/ in Eq. ( 5 12) is replaced by rRz0/ For each case, quadratic and linear convergence orders are noticed in Fig. 5 9 for E2 and E2_qw, respectively, which are consistent with the results in Figs. 3 27 and 3 28 for curved boundaries. Secondorder accuracy of the temperature field is thus preserved with the present axisymmetric LB model, and the first order PAGE 165 165 accurate boundary heat flux is due to the irregularly distributed link fractions in both the radial and axial directions caused by the curved boundary (Li et al. 2013a, 2013b) 5. 3.3.2 Neumann boundary condition on the wall When a Neumann flux condition n rR nDrDnP 00 is imposed on the sphere surface the analytical solution for is the same as that in Eq. ( 5 16) and n = 4 is again studied. To implement the Neumann condition treatment in Eq. ( 3 38), the boundary flux n in each lattice velocity direction is required. The Cartesian decomposition method discussed In Chapter 3 to obtain n based on the normal flux n on a curved boundary is used in the present case. The L2 norms E2 and E2_tw defined in Eqs. ( 5 10) and ( 5 14), respectively, are computed for this test and their convergence is shown in Fig. 5 9 where on average a first order convergence is obtained for both E2 and E2_tw. As demonstrated earlier in Fig. 5 8 for Dirichlet problems with curved boundaries, the evaluation of the boundary flux n is only first order accurate even with a secondorder accurate temperature field obtained. The same formula for evaluating the boundary flux n is directly used in the Cartesian decomposition method for Neumann problems, thus some first order error terms are introduced in the boundary condition treatment. The numerical tests in this section serve as an excellent example to show that even with a secondorder LB model implemented, the overall accuracy of the LB solutions is directly affected by the relevant boundary condition treatment. PAGE 166 166 5. 3.4 Unsteady heat conduction in a circle For unsteady heat conduction inside a circle, the governing equation simplifies to rr rrD D trrrr A sinusoidal boundary condition ( t r = R0) = f ( t ) = sin( t ) is specified and the analytical solution for can be found in Eq. (370). This axisymmetric heat conduction problem has been studied in Section 3.4.4 to assess the thermal bound ary condition treatment and in Section 4.3.3 for heat transfer evaluation on curved boundaries, respectively. In both studies the computational domain is in the 2D circular plane. In this work, as the axisymmetric LB model is used, the computational domain is on the r z plane and only a few lattice grids ( Nz = 5 in the results shown below) are used in the z direction since the variations are only in the r direction. Both r = 0 and r = R0 Following Section 3.4.4, the time averaged error norm for in the interior of the domain on the r z plane is defined as rzrzE dt1/22/22 LBEex0,12 ( 5 17) The results of E2 defined in Eq. ( 5 1 7 ) versus the grid resolution, 1/ R0, at different relaxation coefficients of D are shown in Fig. 5 10. Secondorder convergence is observed for each case in Fig. 5 10 This demonstrates that the proposed LB model can be effectively applied to simulate tim e dependent axisymmetric convection diffusion problems and the secondorder accuracy is preserved. PAGE 167 167 5. 3.5 Steady natural convection in an annulus between two coaxial vertical cylinders In the previous tests, the velocity fields are specified. With the objec tive to show that the proposed LB model for the axisymmetric CDE of can be effectively coupled with the athermal LB models for axisymmetric hydrodynamic equations, we consider in this test the natural convection in an annulus formed by two coaxial vertic al cylinders. The velocity field is solved using an MRT version of Zhous athermal axisymmetric LB model (Zhou 2011) and the temperature field is simulated with the present LB model at each time step. The two models are iterated to reach their steady states. The computational domain, the thermal boundary conditions, and the schematic layout of the square lattice are sketched in Fig. 5 11. Both the radius ratio Ro/ Ri and the aspect ratio H /( Ro Ri) are equal to 2.0 in the present simulations. The vertical c ylinder surfaces are maintained at constant temperatures Ti and To, respectively with Ti > To, while the top and bottom walls are insulated. Invoking the Boussinesq approximation, the fluid properties in the annulus are assumed to be constant, except for t he density in the buoyancy term in the momentum equation. The Prandtl number is fixed at Pr = 0.7, and the characteristic Rayleigh number is defined as Ra = ( Ti To)(Ro Ri)3/ where g is the gravitational acceleration, is the thermal expansion coefficient, is the thermal diffusivity and is the kinematic viscosity of the fluid. The effect of the temperature distribution on the velocity field is tak en into account by adding a body force term z = ( T Tm) in the momentum equation, where Tm = ( Ti + To)/2 is the average temperature. It should also be noted that the present coupling between the athermal and thermal LB models is only on the macrosco pic level, and the direct coupling of the distribution functions for the velocity PAGE 168 168 and temperature is eliminated. This makes the implementation of the LBE models more straightforward and convenient. The convergence is checked by monitoring the variations of the velocity and temperature values at some representative nodes, and the relative differences between successive iterations are at least less than 108 for both the velocity and temperature fields. The streamlines and isotherms at steady states are shown in Figs. 5 12 and 5 13, respectively, for Ra = 103, 104 and 105 with Nr Nz = 202 402. The present results agree well with those reported in (Li et al. 2009, de Vahl Davis and Thomas 1969, Kumar and Kalam 1 991, Venkatachalappa et al. 2001) To quantify the numerical results, the following surfaceaveraged Nusselt numbers on the cylinder walls are defined and compared with the published results in (Li et al. 2009) and (Venkatachalappa et al. 2001) HiiiioRTNu dzHTTr0() and HoooioRTNu dzHTTr0() ( 5 18) It is worth mentioning that in this work the boundary temperature gradients in Eq. ( 5 18 ) are evaluated on the basis of the temperature distribution functions as proposed in Section 4.2.3 rather th an using any finitedifference schemes on the temperature field. And the overall heat transfer on the walls is obtained using the technique proposed in Chapter 4. Both the heat flux and heat transfer evaluations are secondorder accurate as verified in Cha pters 3 and 4 for straight boundaries A grid refinement study is first carried out and the Nusselt numbers Nui and Nuo defined in Eq. ( 5 18) obtained from different mesh sizes are listed in Table 5 1. The numerical results for Nui and Nuo evaluated on the inner and outer surfaces are almost identical to each other for all cases, demonstrating the excellent conservation property PAGE 169 169 of the LBE scheme. It is emphasized that such high consistency between Nui and Nuo can only be obtained after the numerical soluti ons are fully converged. In the present test, all the boundaries are placed in the center of the lattice links ( = 0.5), thus the second order accuracy is preserved for the LBE method in the interior, for the bounceback hydrodynamic and thermal boundar y conditions and the heat transfer evaluation on the straight boundaries. Hence more accurate results can be obtained using the Richardson extrapolation of the raw secondorder results before they are rounded in Table 5 1. In the present case, nnNu extrapolated finecoarse/1 where fine coarse2/2 is the grid resolution ratio, and n = 2 is the order of convergence of the LBE solutions. Table 5 2 compares the extrapolated Nusselt number with the published data in (Li et al. 2009) and (Venkatachalappa et al. 2001) Both the results for Nui and Nuo are presented in (Li et al. 2009) and the results from (Venkatachalappa et al. 2001) are the average Nusselt numbers on the two cylinder surfaces, i.e., Nu = ( Nui + Nuo)/2. The present results in Table 5 2 ar e in good agreement with those from (Li et al. 2009) and (Venkatachalappa et al. 2001) demonstrating that the present LB model can be effectively coupled with the athermal axisymmetric LB models and the interaction between the velocity field and the therm al field is accurately captured. 5. 3.6 Swirling thermal flows in a vertical cylinder Confined axisymmetric flows driven by a rotational velocity is considered in this test to further demonstrate the utility of the present LB model when it is coupled with t he axisymmetric LB model for the velocity field. The fluid inside the cylindrical container is assumed to be incompressible and the Boussinesq assumption is valid for the present PAGE 170 170 flow. The top disc rotates at a constant angular speed and the rest of the solid boundaries are stationary. The top and bottom discs are maintained at constant temperatures, Th and Tc, respectively ( Th > Tc), while the side walls are insulted. The swirling flows are characterized by the Reynolds number Re and the Richardson number Ri, defined as Re = R2 / and Ri = Th3/( 2), ( 5 19) where R is the radius of the cylinder, T = Th Tc, and h is the aspect ratio of the height to the radius of the cylinder, h = H / R The flow pattern and heat transfer of the swirling flows in the range of 102 Re 3 103 and 0 Ri Pr = 1.0 and cylinder aspect ratio h = 1.0 were investigated in detail by Iwatsu (2004) using the vorticity stream function formulation. Zheng et al. (2010a) also simulated the same swirling flow problem and compared their results at Re = 2000, and Ri = 0 and 1.0 with that in ( Iwatsu 2004) to validate their axisymmetric thermal LB model. In the present simulations, the radial and axial velocity components ur and uz are solved using Zhous axisymmetric LB model (Zhou 2011) with an MRT collision operator (see Appendix C ). While it is noted that the governing equation for the azimuthal velocity, u, is also a CDE (Zhou 2011, Graebel 2001) rr u uu uuuuuutrzrzrrrr222222 ( 5 20) where is the kinematic viscosity of the fluid. Eq. ( 5 20) is very similar to Eq. ( 5 1 ) when u is considered as a scalar variable. The source terms Sa 1 and Sa 2 in Eq. ( 5 2) can be modified as rauSurr122 and auSrr2 to account for the scalar PAGE 171 171 variable and its gradient terms in Eq. ( 5 20), respectively. Thus the azimuthal velocity component in rotational flows can als o be solved using the present LB model for axisymmetric CDEs. In this study, Zhous LB model (2011) for the radial and axial velocities and the present LB model for the azimuthal velocity and the temperature are coupled at each time step to reach the steady states. The following results are obtained with mesh sizes of 150 150 for ( Re Ri) = (2000, 0), and 200 200 for ( Re Ri) = (2000, 1.0) after a grid refinement study. The contours of the azimuthal veloci ty component u and the isotherms are shown in Figs. 5 14 and 5 15 for the extreme values of Ri = 0 and 1.0, respectively, at Re = 2000. The swirling flows driven by the rotating top disc are clearly seen, and thin boundary layers are observed near the top disc in the c ontours of u. At Ri = 1.0, the effect of the buoyancy force becomes significant and the isotherms in Fig. 5 15 (b) become more horizontal than that for Ri = 0, especially in the half region near the bottom. It is also noticed that the magnitude of u is v ery small and thus not shown in the lower portion of the cylinder at Ri = 1.0 in Fig. 5 15 (a). The present flow pattern and isotherm distributions are very similar to those reported in (Zheng et al. 2010a, Iwatsu 2004) To further quantify and validate the present results, the profiles of the velocity components and the temperature at r / R = 0.8 are compared with those from (Iwatsu 2004) in Figs. 5 16 (a d). It should be noted that all the velocity components are normalized by and the dimensionless temperature is T* = ( T Tm)/(Th Tc) with Tm = ( Th + Tc)/2. Good agreement is obtained for each case in Fig. 5 16 The rapid changes of the velocity and temperature profiles across the boundary layers are fully captured in the pr esent simulations. The present LB model can thus be applied to axisymmetric PAGE 172 172 rotational flows where both the temperature and azimuthal velocity fields can be simulated with the proposed LB model for the axisymmetric CDEs. 5. 4 Summary and Conclusions A robus t multi relaxationtime (MRT) lattice Boltzmann model for the axisymmetric convectiondiffusion equation (CDE) is proposed. The axsiymmetric CDE is transformed to a standard CDE in the Cartesian coordinate system and the additional terms due to the coordinate transformation are treated as source terms at the macroscopic level. The scalar gradient in the source terms is formulated using the local moment of the non equilibrium components of the distribution functions without the need for finitedifference type of calculations. When the CDE for the scalar variable is coupled with the hydrodynamic equations, the radial axial velocity field is solved with an MRT version of the athermal axisymmetric LB model (Zhou 2011) The macroscopic velocity components, inst ead of the velocity distribution functions, are coupled in the evolution equation of the distribution functions for In this manner, there is no explicit coupling between the two sets of distribution functions, making the present model much easier to implement and more robust. The secondorder accuracy of the proposed LB model is validated with a series of numerical tests. The effect of the curved boundary on the convergence order of the LB solutions is also investigated and verified. The agreement between the present LB results and published numerical computations in the literature for the natural convection problem and the confined swirling flows also demonstrates that the present model can PAGE 173 173 be effectively coupled with athermal LB models to yield accurate results for axisymmetric thermo hydrodynamic problems. PAGE 174 174 Table 5 1. Surfaceaveraged Nusselt numbers at different grid resolution. Average Nusselt numbers iNu oNu mesh Ra = 10 3 Ra = 10 4 Ra = 10 5 102 202 1.69200, 1.69200 3.21673, 3.21673 5.83096, 5.83095 127 252 1.69206, 1.69206 3.21607, 3.21607 5.81486, 5.81486 152 302 1.69210, 1.69210 3.21572, 3.21572 5.80 636, 5.80636 202 402 1.69213, 1.69213 3.21531, 3.21531 5.79802, 5.79802 Table 5 2. Comparison of the present surface averaged Nusselt numbers with published results. Average Nusselt number Nu Ra Pre sent Li et al. 2009 iNu Li et al. 2009 oNu Venkatachalappa et al. 2001 (extrapolated) (101 201) (41 81) 10 3 1.6922 10 4 3.2148 3.216 3.218 3.1629 10 5 5.7873 5.782 5.787 5.8818 PAGE 175 175 0 1 2 3 4 x1x2 Figure 51. Discrete velocity set { e} for the D2Q5 LB model. e = (0, 0), (1, 0), and (0, 1). F igure 5 2 Relative L 2 error norm E2 for the interior field of versus the grid resolution, 1/ Ri, for the 1 D annulus convectiondiffusion problem. PAGE 176 176 z r L y 0.5 y R0 eF igure 53 Schematic of the square lattice in the r z plane of a 3D circular pipe. The axis of the pipe, r = 0, is fixed at the center of the link with = 0.5, while the pipe surface, r = R0, is placed a F igure 54 Contours of in the r z plane at Nr = 42, Nz = 83, and = 0.5 for the steady convectiondiffusion in a pipe with a Dirichlet boundary condition on the wall. PAGE 177 177 (a) (b) (c) (d) F igure 55 Relative L 2 error norms versus the grid resolution, 1/ R0, for the steady convectiondiffusion in a pipe with a Dirichlet boundary condition on the wall: (a) E2 for the interior field of (b) E2_qw for the boundar y flux, (c) E2_gr for the interior r and (d) E2_gz z PAGE 178 178 (a) (b) (c) (d) F igure 56 Relative L 2 error norms versus the grid resolution, 1/ R0, for the steady convectiondiffusion in a pipe with a Neumann boundary condition on the wall: (a) E2 for the interior field of (b) E2_tw for the boundary values of (c) E2_gr r and (d) E2_gz for the interior axial gradient z PAGE 179 179 z R0 r x y P F igure 57 Schematic representation of the computational domain and the lattice layout for the heat conduction inside a sphere. F igure 58 Relative L 2 error norms E2 and E2_qw for the interior field of and the boundary flux, respectively, versus the grid resolution, 1/ R0, for the steady heat conduction inside a sphere with a Dirichlet boundary condition. PAGE 180 180 F igure 59 Relative L 2 error norms E2 and E2_tw for the interior field and boundary values of respectively, versus the grid resolution, 1/ R0, for the steady heat con duction inside a sphere with a Neumann boundary condition. F igure 510. Time averaged error norm E2 for the interior field of versus the grid resolution, 1/ R0, for the unsteady heat conduction in a circle with a temporal Dirichlet boundary condition. PAGE 181 181 g H TiTo zT = 0 zT = 0 RiRo z F igure 511. Schematic depiction of the computational domain ( Ro/ Ri = 2.0, and H /( Ro Ri) = 2.0) and the thermal boundary conditions on the walls all placed in the PAGE 182 182 (a) (b) (c) F igure 512. Streamlines in the annulus between two coaxial vertical cylinders at (a) Ra = 103, (b) Ra = 104, and (c) Ra = 105. (a) (b) (c) F igure 513. Isotherms in the annulus between two coaxial vertical cylinders at (a) Ra = 103, (b) Ra = 104, and (c) Ra = 105. PAGE 183 183 (a) (b) F igure 514. Contours of the azimuthal velocity u (a) and isotherms (b) of the confined swirling flow at Re = 2000 and Ri = 0. (a) (b) F igure 515. Contours of the azimuthal velocity u (a) and isotherms (b) of the confined swirling flow at Re = 20 00 and Ri = 1.0. PAGE 184 184 (a) (b) (c) ( d ) F igure 516. Velocity components and temperature profiles along the vertical line at r / R = 0.8 for Re = 200 0: (a) radial velocity ur, (b) azimuthal velocity u, (c) axial velocity uz, and (d) temperature. PAGE 185 * The Monte Carlo ray tracing modeling of the radiative heat transfer inside the cavity by Singh A. (Singh 2013) based on the work of Petrasch J. (2002) is greatly appreciated. 185 CHAPTER 6 ENERGY TRANSPORT IN A SOLAR THERMOCHEMICAL REACTOR FOR HYDROGEN PRODUCTION* 6 .1 Introduction Thermochemical cycles using concentrated solar energy as the heat source to produce hydrogen or synthesi s gas (syngas, mixture of hydrogen and carbon monoxide) have been demonstrated to be a promising and efficient strategy for renewable energy pursuit (Steinfeld 2002, Steinfeld 2005, Abanades et al. 2006, Kodama and Gokon 2007, Perkins and Weimer 2009, Lout zenhiser et al. 2010, Mehdizadeh et al. 2012). The two step H2O/CO2splitting cycle by a metal oxide process is one of most studied cycles in the literature (Steinfeld 2002, Steinfeld 2005, Abanades et al. 2006, Kodama and Gokon 2007, Loutzenhiser et al. 2010, Mehdizadeh et al. 2012). In the two step H2O/CO2splitting reductionoxidation (Redox) reaction, one step is known as the regeneration or thermal reduction, in which a metal oxide is reduced, releasing lattice oxygen. During the subsequent oxidation s tep, the reduced and activated material is oxidized by reacting with water/CO2 to produce hydrogen/CO. The present study focuses on hydrogen production. Denoting M as the used metal and MO the corresponding metaloxide, the reduction and oxidation steps can be represented by oxidizedreduced2MOMO+O Reduction Step ( 6 1) reduced2 oxidized2MOHOMOH Oxidation Step ( 6 2 ) The reduction step is endothermic and the oxidation step is exothermic, both steps require high temperatures and the temperature for the reduction step is normally PAGE 186 186 much higher than that for the oxidation step. For example, most of the operation temperatures required for the reduction step as listed in Abanades et al.s work (2006) is greater than 1000 oC. Due to the high temperature requirement especially during the reduction step, concentrated highflux solar energy is considered as an effective heat source for which a solar cavity receiver configuration is usually required. Klausner et al. (2011) have been developing a solar thermochemical reactor prototype for cyclically producing renewable and carbonneutral fuels Main features in their concept include: (i) a magnetically stabilized structured bed of metal/metal oxide particles to improve cyclical behavior and minimize soli d particle losses, (ii) low pressure operation to bring down reduction temperatures below 1500 oC, and (iii) a windowless dual cavity reactor concept that avoids reliability issues associated with directly irradiated, windowed cavity type reactors ( Klausner et al. 2011). Heat transfer mechanisms in the reactor include radiation, conduction, convection and chemical reaction. A fundamental understanding of the energy transport in the solar cavityreceiver and the heat flux and temperature distributions on the surfaces and in the interior of the reactors (absorbers) is essential to the optimal design of the cavityreactor Radiative heat transport inside the cavity can be modeled either by the radiosity analysis method (Wieckert et al. 2003, Mller et al. 2008) or using the Monte Carlo ray tracing (MCRT) technique (Petrasch et al. 2007, Melchior and Steinfeld 2008). The MCRT technique, which is intrinsically threedimensional, is preferred in the present work since it can be applied to simulate the radiative heat transfer in arbitrary geometric configurations and thus be used for optimal design. With the heat flux boundary PAGE 187 187 condition on the absorber surfaces obtained from the radiative heat transfer analysis, the energy transport inside the absorbers, including conduction, convection, radiation and chemical reaction in the porous bed, can be simulated. Realizing the reemission from the absorber surfaces becomes significant when the surface temperatures are high enough, this reemission should be fed back to the radi ative analysis inside the cavity. Thus these two models should be coupled to balance the energy transfer in the whole cavityreactor system. Wickert et al. (2003) applied the radiosity method to calculate the solar energy absorption efficiency for a twoca vity configuration. Mller et al. (2008) also used the radiosity method to simulate the radiative heat transfer between the cavity wall and the absorber surfaces. Melchior and Steinfeld (2008) used the Monte Carlo method to study the radiative transfer ins ide a solar cylindrical cavity with diffusely or specularly reflecting inner walls, containing one singletube or multipletube absorbers, and a selective window or a windowless aperture. The location and relative dimensions of the absorbers are optimized for maximum energy transfer efficiency or for maximum absorber surface temperature (Melchior and Steinfeld 2008). It is noticed that the heat transfer inside the reaction chamber or the absorbers is not fully simulated in both of the work by Wieckert et al (2003) and by Melchior and Steinfeld (2008). Instead, they considered the net power absorbed by the reactors as a parameter and examined the energy transfer efficiency at different values of that parameter. Maag et al. (2010a, 2010b) developed a radiation/convection/conduction coupled solar reactor heat transfer model for syngas production and methane decomposition. Their reactor model consists of a cavity receiver model and an absorber tube model, which are coupled by the net PAGE 188 188 heat sink into each tube bas ed on the outer surface temperatures of the tubes. A constant surface temperature was assumed on each outer tube in the coupled model (Maag et al. 2010a, 2010b) considering the high thermal conductivity of the graphite tubes. The thermochamical reduction s tep in the redox cycle is very sensitive to change of temperature at high temperatures. In addition, the radiative emission from the absorber surfaces, which is directly related to the surface temperature distributions, contributes to and significantly aff ects the overall radiative transfer analysis inside the cavity. The simulation and determination of the actual temperature field inside the absorbers as well as on the absorber surfaces are thus necessary in the energy transport analysis of the cavity reac tor system. A coupled radiationconvectiondiffusion model is developed in this work. The boundary heat flux distributions on the surfaces are obtained from an MCRT radiation mode, and the energy equation inside the absorbers is solved with the thermal lat tice Boltzmann equation (TLBE) method discussed in Chapters 3, 4 and 5. T hese two models are fully coupled and iteratively solved. The geometry of the curved boundary is preserved using the boundary condition treatment in Chapter 3. 6 .2 Horizontal Cavity R eceiver Solar Reactor The geometrical configuration of the present solar reactor is schematically depicted in Fig. 61. Its main features include a cylindrical cavity with its axis being horizontal and parallel to the central incident solar flux, a circular aperture through which the incident solar flux enters the cavity, and a series of tubular absorbers at the circumference of t he cavity. A stabilized porous structure or packed bed of particles is loaded inside the absorbers for chemical reaction and it i s indirectly heated by the solar PAGE 189 189 radiation. To minimize the heat loss from the outer surface of the cavity, a thick insulation layer covering the cavity is used. 6 3 Heat Transfer Analysis A comprehensive heat transfer analysis in the hightemperature sola r thermochemical reactor includes the components of radiation, conduction, convection, and chemical reaction. This study is focused on the reduction step during which concentrated solar energy is used for high temperature requirement. The University of Flo rida Solar Simulator is able to provide high flux up to 10 kW ( Erickson 2012, Klausner 2011) and it is used as a radiation source. The distribution of the solar flux arriving at the cavity aperture from the Solar Simulator is determined utilizing an MCRT b ased VEGAS code. The details about the VEGAS code can be found in (Petrasch 2002). The center of the present analysis is on the energy transfer inside the cavity reactor with a given constant input solar power and a fixed flux distribution at the aperture. An MCRT based radiation model developed by Singh (2013) is applied to simulate the radiative transfer between the cavity wall and the absorber surfaces. The details of the radiation model are referred to (Singh 2013) and some basic concepts are highlighte d here for completeness. The energy transfer inside the absorbers is modeled using the thermal lattice Boltzmann equation (TLBE) method. 6 3 .1 Radiative Heat Transfer 6 3 .1.1 Radiation between the cavity wall and the absorbers Assuming the medium filling the space between the cavity wall and the absorbers is nonpart icipating, the radiative heat flux on a surface using the MCRT technique is governed by (Modest 2003, Melchior and Steinfeld 2008) PAGE 190 190 ''' 4 '4' 'd()()()()()dddAdAAFq Arrrrr (6 3) where q( r ) is the local surface heat flux in the normal direction, r is a spatial vect or, ( r ) is the total hemispherical emissivity of the surface at r is the Stefan Boltzmann constant, T is the surface temperature, A is the surface area of the enclosure, and 'ddAdAF is the generalized radiation exchange factor between surface elements 'dA and dA In Eq. (6 3), the first term on the right hand side represents the emission from the surface at r with temperature T ( r ), an d the integrand in the second term is the fraction of the energy, originally emitted from the surface at r that eventually gets absorbed at location r If the enclosure is not closed, such as the aperture in the cavity, some artificial closing surfaces must be introduced. Since the opening aperture is irradiated by the solar simulator, it is replaced by a nonreflecting surface with zero emissivity for all angles but the solar angle. The MCRT method traces the history of a statistically meaningful random sa mple of photons from their point s of emission to their points of absorption within or exit ing from the volume (Modest 2003) The emissive power is equally distributed to the individual photon rays, which are created in random directions. The photon rays wi th concentrated solar power, solarQ enter the cavity through the aperture, and the power carried per ray is solarrayaperturerayapertureQQN (6 4) PAGE 191 191 where Nray aperture is the number of rays at the aperture. A fter encountering multiple reflections inside the cavity the photon rays either get absorbed by the absorber surface s or the cavity wall, or leave the cavity from t he aperture. Note that the distribution of incident rays at the aperture are obtained from the VEGAS code (Petrasc h 2002). As the surface temperatures of the absorbers and the inner cavity walls increase with the irradiation absorbed, the emission from the their surfaces becomes important and it should be accounted for in the radiative transfer analysis. The absorber surfaces and the cavity inner walls are thus divided into a large number of differential elements, each assumed at a uniform temperature. The emissive power per ray from the surface element Ai is thus 4 absorber rayabsorber rayabsorberrayabsorber i i iii iiQ Q NN (6 5) Both the absorber surfaces and the cavity wall s are assumed to be diffusely reflecting and emitting. 6 .3 .1. 2 Radiation is the porous structure inside the absorbers The reactive material filled in the absorbers is in the form of optically thick p orous structures thus the Rosseland diffusion approximation for the radiation inside the bed is valid (Modest 2003, Mller et al. 2008). The radiative heat transfer is then simplified to a diffusion process with an equivalent thermal conductivity 2 316 3R Rn kT (6 6) where n is the refractive index of the medium, and R is the Rosselandmean extinction coefficient. The high value of R results in attenuation of the incident radiation mostl y at PAGE 192 192 the surface layer of the bed. The radiative exchange inside the bed is thus mainly due to the local emission and absorption of thermal radiation, justifying the application of the Rosseland diffusion approximation for optically thick media (Mller et al. 2008). 6 .3 2 Convective Heat Transfer In the reduction step of the twostep H2O/CO2splitting cycling the primary gas is the oxygen released during the reaction. Considering the oxygen is removed from the system at a very low speed, the convective heat transfer is neglected in the present simulations. The same assumption was made in those previous works (Mller et al. 2008, Loutzenhiser et al. 2010). In the subsequent oxidation step, however, as relatively high H2O/CO2 flow rate is passing through the bed, the convective heat transfer is significant and should be included in the overall heat transfer analysis. 6 .3 3 Conductive Heat Transfer For the reduction step, the heat transfer inside the absorbers is simplified to a heat conduction problem. Recalling Eq. (6 6), the effective thermal conductivity for heat conduction can be approximated by (Mller et al. 2008, Loutzenhiser et al. 2010) 23eff cond cond163RRnkkkTk (6 7) where the thermal conductivity condk for the por ous medium without considering radiation can be obtained either from experimental measurement of the samples at room temperature, or from tomography based porescale numerical simulation. 6 3 4 Endothermic Chemical Reaction The heat consumption due to the endothermic dissociation of the reactive material can be represented by PAGE 193 193 2''''''chemMOM OMO(0.5)qrhhh (6 8) where '''MOr is the volumetric reaction rate of the metal oxide MO such as FeO or ZnO, h is the specific enthalpy of the reactants or products. The reaction rate is usually modeled by an Arrhenius type rate law (Mller et al. 2008, Loutzenhiser et al. 2010). Here we follow the oxidation reaction model in (Mehdizadeh et al. 2012b) and use a similar reac tion model for the reduction step '''MO0exp()aErkfxRT (6 9) where t he activation energy, Ea, for reduction is e xperiment ally determined, the preexponential factor, k0, is a constant fitted from experimental data, and f ( x) is a dimensionles s function that depends on the reaction mechanism, x is the mass fraction of the reacted material (Mehdizadeh et al. 2012b) 6 4 TLBE Simulation of Heat Transfer in Absorbers The energy conservation equation for the reactive porous structure filled in the cylindrical absorbers can be written as '''eff eff eff chem211pTTTT trrrr zz (6 10) where and pc are the density and specific heat capacity of the solid matrix, respectively, and ''' chemq is the rate of volumetric heat sink given by Eq. (68). Eq. (6 10) is a diffusion equation and it can be solved using the cubic D3Q7 TLBE model (Yoshida and Nagaoka 2010) with the curved geometry preserved using the boundary condition treatment by Li et al. (2013a, 2013b). PAGE 194 194 6 4 1 Scaling from the physical units to the LBE units It is realized that all the numerical tests in Chapters 3 5 are presented in LBE units where both the lattice spacing and the time step are unity. This serves as a clean plate on which the analysis and numerical validation are conveniently conducted. To simulate a particular real world problem, however, the parameters/variables in physical units should be properly scaled and converted into those in LBE units that are also suitable for LBE algorithm implementation. These two sets of unit systems are connected by the same dimensionless units, such as abs abs phys LBE// rRrrR (6 11a) *abs absphys LBE//zRzzR (6 11b) phys LBE// TTTTT and (6 11c) abs abs* phys LBE abs abs *rR rRRR TTT TrrTr (6 11d) T is the reference temperature difference. Introducing basic scaling factors for the characteristic length and time scales from the physical units to the LBE units as physLBE 0absabs/xL and 0physLBE0//ttT (6 12) where t0 i s the characteristic physical time and Nt is the number of time steps. The diffusion coefficient in LBE simulations then becomes 0 eff0LBE2200ttpxxTDL (6 13) PAGE 195 195 In TLBE simulations the relaxation coefficient ij is related to the diffusion coeffici ent as 212()ijij ij (see Eq. 314) in order to correctly recover the macroscopic convectiondiffusion equation (CDE). Also as demonstrated in Figs. 38, 4 3 and 44, the errors of the LBE results increases with larger values of ij especiall y for those cases with ij > 1.0. Accurate LBE simulations thus require that ij be within a reasonable region. For highdiffusioncoefficient CDEs, the diffusion coefficient should thus be properly rescaled when the TLBE method is employed. The details about the rescaling technique are given in the following section. 6 4 2 Rescaling for highdiffusioncoefficient CDE in TLBE In hydrodynamic LBE simulations, the macroscopic hydrodynamic equations are recovered from the kinetic based LBE when the mean free p ath is much smaller than the lattice grid size. As the relaxation coefficient f is defined as the ratio of the relaxation time to the lattice time step ( f = f/ ), a high f value corresponds to a large relaxation time and large particle mean free path. Hence higher f values would yield larger numerical errors in LBE solutions, as examined in (Chen et al. 1992, Noble et al. 1995, Mei et al. 2002, Chun and Ladd 2007). Similarly, in thermal LBE simulations, the deviation of the TLBE results from the solutions to the CDE increases with large relaxation coefficient D values (corresp onding to higher diffusion coefficient in the governing CDE) as reported in (Li et al. 2013a, Yoshida and Nagaoka 2010). In this work, we show that the diffusion coefficient D in a CDE can be properly rescaled to ensure moderate D values to have good numerical stability and accuracy. It should also be noted that with a gain of numerical accuracy, using smaller relaxation values of f and D in LBE simulations would require more lattice PAGE 196 196 time steps to reach steady state or time dependent solutions (Li et al. 2013a, Yoshida and Nagaoka 2010). We rescale the large diffusion coefficient D with a parameter so that the general CDE (31) becomes ()()j ijjijvDGtxxx (6 14) with t /jjvv /ijijDD and / GG (6 15) Equation (615) shows that the time, velocity and source term are rescaled accordingly. Note that the spatial vector x is not rescaled in Eq. (614) and the characteristic Fourier number and Pclect n umber remain the same after the rescaling, i.e., 22//ccFoDtLDtLFo and // PeULDULDPe with tc U and L being the respective reference time, velocity and length scales. The Dirichlet boundary condition in Eq. (32) is not affected by the rescaling, while the Neumann boundary condition in Eq. (33) should be rescaled as(Li et al. 2013d) jniij jjvnDnx (6 16) After the heat equation (610) is converted to a CDE in the LBE units, it will be rescaled according to Eq. ( 6 14) and solved. 6 4 3 Computational Procedure for the Coupled Model The computational procedure for the coupled conductionradiation model is summarized as below: PAGE 197 197 1. The incident solar power solarQ and the distributions of the rays at the aperture of the cavity is obtained from the VEGAS ray tracing code for the solar simulator (Petrasch 2002); 2. Start with initial temperature distributions on the absorber and cavity surfaces and the emission power from each surface is computed; 3. The MCRT radiation model (Singh 2013) is used to obtain the net heat flux, which is a combination of the absorbed and emitted radiative energy, on the absorber surfaces and the cavity walls; 4. The temperature distributions in the interior of the absorbers, the cavity and the insulation layer are simulated with the net heat flux boundary conditions using the TL BE model; 5. The heat sink due to chemical reaction is calculated with the temperature field in the absorbers obtained; 6. Temperature distributions on the absorber surfaces and the cavity walls are updated so that the radiative emission from those surfaces is updated; Steps (3 6) are iteratively conducted during each time step. 6 5 Simulation Results and Discussion Simulation results from three cases with different aspec t ratios of the cavity are presented in this section. The baseline parameters for those thr ee cases are summarized in Table 6 1. A constant solar power of solar,1Q = 10 kW is applied to heat up the reactor from room temperature 27 oC (300 K) to the reduction temperature at 1450 oC, then a lower solar input solar,2Q is chosen to maintain the temperature of the reactive material at about 1450 oC. Considering the absorbers are uniformly distributed along the PAGE 198 198 az imuthal direction inside the cylindrical cavity, a representative absorber is selected and the heat flux boundary condition is assumed to be the average value of the heat fluxes on all the absorbers. Therefore only the interior temperature field in one abs orber is simulated in the coupled model, which greatly reduces the computational time. Figures 6 2 (a) and (b) show the evolution of the average reactor bed t emperatures Tbed (averaged over the 3D interior of the absorbers) with time for the three cases l isted in Table 6 1. The evolution history starts from 300 K and the desirable operating temperatures, 1200 oC for the oxidation step and 1450 oC for the reduction step, are also highlighted in Fig. 6 2 (a). A 10 min reduction duration time is used and the decrease of the solar power is illustrated in Fig. 6 2 (b). Define the solar to fuel efficiency as 0solar,11solar,22pump/mcmYVHQtQtW (6 17) where m is the total mass of the bed material, Y is the total hydrogen yield per kg of bed material, Vm is a constant Vm = 22.41 dm 3 /mol, 0cH is the standard enthalpy of combustion of hydrogen 0 cH = 286 kJ/mol, ( solar,1Q t1) is the input solar energy to ramp up Tbed from 1200 oC to 1450 oC time, ( solar,2Q t2) is the input solar energy to maintain the reduction temperature of 1450 oC for 10 mins (ideally the oxidation step is exothermic and no solar input is required), and Wpump is the energy consumption of the vacuum pump. According to the benchmark data obtained from our experiments for Cobalt ferrite based reactive material, a bulk hydrogen yield of Y = 2.5 dm 3 /kg of total material is assumed, the eff iciency for those cases are thus 1 = 11.12%, 2 = 9.91%, and 3 = PAGE 199 199 7.30%. The calculations are illustrated in Table 6 2. For cerium based reactive material, faster reaction and higher production yield are noticed, thus the overall effic iency is expected to be higher To illustrate the incident flux distributions, the contours of incident heat flux on the surface of Absorbe r 1 are shown in Figs. 63 (a), (b) and (c) for Cases I, II, and III, respectively. The incident flux remains the same during the reduction step. With the absorber surface heated up by the incident flux, its temperature increases and radiative emission fro m the surface becomes important. When the average temperature of the absorbers reaches 1450 oC the net heat flux and temperature distributions on the absorber surfaces for all three cases are shown in Figs 6 4 and 6 5, respectively where the results are presented in the ( z ) coordinates (see Fig. 61) The heat flux and temperature distributions are symmetric about the axis at = 0, which is consistent with the axisymmetric arrangement of the absorbers inside the cavity. Along the direction, relatively higher surface temperature is observed for each case at the location near = 0, which is exposed to the incident radiation coming through the aperture (see Fig. 61 ). In the axial z direction, a hot spot of temperature is noticed for each case. This hot spot region corr esponds to the location on the surface that absorbs highintensity incident flux as shown in Fig. 63 The temperature profiles along the axial direction at = 0 for the three cases are compared in Fig. 6 6. To examine the interior temperature field of th e absorbers at Tave= 1450 oC we show the temperature distributions on the xz plane that contains the centerline of Absorber 1 for all three cases i n Fig. 6 7 (the x, y, z axes for Absorber 1 are parallel to PAGE 200 200 that in the xyz Cartesian coordinate system f or the cavity as shown in Fig. 6 1 with the origin shifted to the centerline of Absorber 1). The interior temperature distribution is essential to the understanding of the local chemical reaction rate. The temperature distributions in Fig. 67 implies that further effort can be devoted to minimize the temperature difference in the bed to have more uniform chemical reaction. 6 6 Summary and Conclusions A coupled radiationconduction model is developed for simulating the energy transport in a hightemperature solar thermochemical reactor. The radiative transfer inside the porous structures is approximated as a diffusion process with large effective thermal conductivity at high temperatures. A general technique for rescaling the highdiffusion convection diffusion equation (CDE) and the boundary conditions in the ther mal lattice Boltzmann equation ( T LBE) method is presented to improve numerical accuracy. Radiative transfer between the cavity wall and the absorbers is modeled with the MonteCarlo ray tracing (MCRT) technique. As the surface temperature of the absorbers is solved, actual radiative emission from the absorber surfaces is included in the MCRT model. Simulation results for a horizontal cavity receiver reactor prototype are presented to give insights into the reactor design. PAGE 201 201 Table 6 1. Baseline simu lation parameters for cavity aspect ratio Lcav, in / Dcav = 0.69, 1.01, and 1.16 for Cases I, II and III respectively Parameter Symbol Value Maximum s olar power input solar, maxQ 10 kW Absorber emissivity 0. 8 Number of absorbers n 17, 14, 13 Absorber diameter D abs 50.8 mm Absorber length L abs 246, 305, 330 mm Absorber thickness t abs 3 18 mm Aperture diameter D aperture 50.0 mm Cavity diameter D cav 355, 302, 284 mm Cavity thickness h cav 76.2 mm Insul ation thickness h ins 127 mm Inner cavity length L cav, in 246, 305, 330 mm Outer cavity length L cav, out 398, 457, 482 mm Insulation length L ins 652, 711, 736 mm Spacing between absorber and cavity axes Habs 150, 123, 115 mm Table 6 2. Solar to fuel efficiency calculations for Cases I, II and III. Case I ( L cav, in / D cav =0. 69 ) Case II ( L cav, in / D cav =1.0 1 ) Case III ( L cav, in / D cav = 1 16 ) Mass load (kg) 13.29 13.54 13.62 Benchmark H 2 yield (liter) 33.23 33.85 34.05 HHV of H 2 (kJ) 424.28 432.26 434.81 ( solar,1Q t1) 1200 1450 oC (kJ) 2030.00 2690.00 4226.67 ( solar,2Q t2) ~1450 oC (kJ) 1697.25 1582.20 1638.00 W pump (kJ) 87.53 89.17 89.70 Efficiency (Cobalt ferrite based) 11.12% 9.91% 7.30% PAGE 202 202 Cavity Aperture Insulation Absorbers 1 xy z Figure 6 1 Schematical depiction of the horizon tal solar reactor configuration: (a) cross section view, and (b) front view. (a) (b) Figure 6 2 (a) V ariations of the average bed temperature ( Tbed) with time and (b) v ariations of Tbed and solar power input ( solarQ ) with time at 1200 o Tbed 1450 oC for Cases I, II and III. PAGE 203 203 (a) (b) (c) Figure 6 3 I ncident heat flux distribution (in W/m2) on the surface of Absorber 1 (see Fig. 6 1) for (a) Cases I, (b) Case II, and (c) Case III. PAGE 204 204 (a) (b) (c) Figure 6 4 Net heat flux di stribution (in W/m2) on the surface of Absorber 1 for (a) Cases I, (b) Case II, and (c) Case III. PAGE 205 205 (a) (b) (c) Figure 6 5 Temperature distribution (in K) on the surface of Absorber 1 for (a) Cases I, (b) Case II, and (c) Case III. Figure 6 6 Com parison of t emperature profiles (in K) along ( = 0, z ) on the surface of Absorber 1 for the three cases simulated. PAGE 206 206 (a) (b) (c) Figure 6 7 Interior temperature distribution on the central xz plane of Absorber 1 for (a) Cases I, (b) Case II, and (c) Case III. PAGE 207 207 CHAPTER 7 SUMMARY AND FUTURE WORKS Heat transfer, as one of the most significant transport phenomena in fluidsolid flows, is studied in this dissertation. Various numerical tools have been developed/improved to efficiently and accurately evaluate the heat transfer in fluidsolid systems. Most effort is devoted to the boundary condition treatments and heat transfer evaluation in the thermal lattice Boltzmann equation (TLBE) method. The major contributions in the dissertation are summarized as follows. Heat transfer between colliding particles of various materials is investigated: The singular nature of the thermal field and heat flux at the edge of the contact area is elucidated through analytical solutions. A self similar solution for the thermal field during the initial period of contact is developed. And i t serves as an accurate initial condition for the entire collisional period. A two dimensional asymptotic analysis is presented for the thermal field and heat flux in the contact area in the small Fourier number limit. A closed form expression is developed for the heat transfer prediction as a function of the Fourier number, the thermal diffusivity ratio and the thermal conductivity ratio of the impacting particles. Boundary condition treatments in TLBE method are proposed: Thermal boundary condition treatments based on the bounceback scheme and spatial interpolation are pr esented for both the Dirichlet and Neumann conditions. When second order accuracy is pursued, there is an adjustable parameter for the Dirichlet problem, while the Neumann condition treatment is unique. Three particular schemes for the Dirichlet condition treatment are presented and their stability and accuracy are examined. All three Dirichlet schemes are stable and Scheme 2 gives smaller error in most cases investigated. For curved boundaries, the Dirichlet condition treatment can be applied directly and it leads to secondorder accuracy of the temperature field and first order accuracy of the boundary heat flux. PAGE 208 208 For curved boundaries, the proposed Neumann condition treatment requires the determination of heat fluxes in the discrete velocity directions o f the TLBE model. A general formula to convert the normal heat flux into the fluxes in the discrete velocity directions is derived. O nly first order accuracy is obtained for the temperature field with Neumann conditions due to the combined effect of the curved boundary and the heat flux conversion. For the special case with no tangential variation of temperature along the curved boundary, the heat flux conversion becomes exact and secondorder accuracy of the temperature field can be obtained. Heat transfer evaluation on curved boundaries in TLBE method is studied : An efficient and accurate approach for evaluating local and overall heat transfer rates on curved boundaries is proposed. The boundary heat flux in the discrete velocity directions is obtained w ith the temperature distribution functions at the lattice nodes close to the boundary and the given Dirichlet or mixed boundary condition. Based on energy conservation, the integration of the heat fluxes in the discrete velocity directions with constant e ffective surface area leads to the total heat transfer rate. There are two major benefits for the proposed heat transfer evaluation: no requirement for the normal heat flux determination and no surface area approximation. Numerical tests demonstrate that the heat transfer evaluation method is secondorder accurate for straight walls perpendicular to one of the discrete lattice velocity vectors and first order accurate for curved boundaries. A Multi relaxationtime (MRT) lattice Boltzmann (LB) model for the axisymmetric convectiondiffusion equation (CDE) is proposed: The axsiymmetric CDE is transformed to a standard CDE in the Cartesian coordinate system and the additional terms due to the coordinate transformation are treated as source terms at the macros copic level. The scalar gradient in the source terms is formulated using the local moment of the nonequilibrium components of the distribution functions without the need for finite difference type of calculations. PAGE 209 209 For thermohydrodynamic problems, the pro posed model can be directly coupled with hydrodynamic axisymmetric LB models on the macroscopic level without explicit coupling between the two sets of distribution functions, making the present model much easier to implement and more robust. A coupled radiation conduction model for energy transport in a solar reactor is developed : The radiative transfer inside the porous structures is approximated as a diffusion process with large effective thermal conductivity at high temperatures. A general technique for rescaling the highdiffusion convection diffusion equation (CDE) and the boundary conditions in the thermal lattice Boltzmann equation (TLBE) method is presented to improve numerical accuracy. Radiative transfer between the cavity wall and the absorbers i s modeled with the MonteCarl o ray tracing (MCRT) technique. T he surface temperature of the absorbers is solved and actual radiative emission from the absorber surfaces is inc luded in the MCRT model. Simulation results for a horizontal cavity receiver reac tor prototype are presented to provide insights into the reactor design. Based on the material presented in this dissertation, the following are recommended as future work s: 1. With curved boundaries, the boundary heat flux evaluation for Dirichlet problems and the interior and boundary temperature computations for Neumann problems become first order accurate in the TLBE method. This has been verified with numerous numerical tests. While a thorough analysis of the convergence order for curved boundaries remains open. 2. A detailed parametric study is needed in the radiationconduction coupled model for optimal solar thermochemical cavity reactor design. 3. A more advanced numerical model including, i) fluid flows (especially for the oxidation step), and ii) species transport due to chemical reaction, is highly desirable to serve as a design tool for the solar thermochemical reactor. PAGE 210 210 APPENDIX A VALIDITY OF THE SEMI INFINITE MEDIA ASSUMPTION During the elastic collision of two particles initially at different temperat ures, the deformations of the spheres in the direction normal to the contact area are very small compared with their respective radii according to the Hertzian contact theory (Timoshenko et al. 1970, Johnson 1987). The heat transfer is also limited to a sm all area near the contact region. Sun and Chen (1988) thus employed the semi infinite media assumption by treating the rest of the solids in the z direction (normal to the contact area) far away from the contact region as infinite media (see Eq. (215)). I t is also implied that the solid in the r direction far away from contact region is at infinity (see Eq. (2 16)). Furthermore, the surfaces near the edge of the contact area are assumed to be flat so that the surfaces outside r = R ( t ) can be assumed to be at z = 0 (see Eq. (2 14)). Zhou et al. (2008) argued that for high Fourier numbers (when Fo is close to 10), the above assumptions may not be valid. Thus they treated the two colliding spheres as deforming bodies and thus maintained the geometric fidelit y and the exactness of the boundary conditions using a finiteelement method. Thus it is instructive to examine the conditions needed for the semi infinite media assumption to be valid. There are three issues of concern: i) Is the region of influence in t he z direction much smaller than the length scale of the body in the z direction? ii) Is the region of the influence in the r direction much smaller than the length scale of the body in the r direction? PAGE 211 211 iii) Is the region of the influence in the r directio n sufficiently small so that the curvature effect of the surface near the contact region can be assumed negligible or the surface can be assumed to be flat? These issues are addressed below based on scaling analyses and the computational results obtained i n this study using the post priori method. i) For a solid with thermal diffusivity the relevant length scale of the thermal diffusion in the z direction is on the order of ()cO in which tc is the contact duration. For the rest of the solid to be treated as an infinite medium in the z direction, it is necessary t hat ici ( i = 1, 2) (A 1) Examination of the temperature contours over a wide range of Fo for particles of the same material shows that the temperature contour lines in the direction varies very little during the entire coll ision period. This confirms that 1 is the correct scale in the z direction. Since 1/ is applied in this work the domain of the influence in the direction is limited to ~ 4 5. Thus when (A 1) is satis fied, the z direction can be treated as semi infinite for each body. ii) In the r direction, the region influenced by the impact can be estimated as ()ccOR where Rc is the maximum contact radius. In the limit of zero Fo the region is limited precisely to the contact area, r Rc; for large Fourier numbers, the thermal diffusion length in the r direction is then ()cO beyond the point of contact. At Fo = 50, the isothermal lines on the sur face of particle 1 at t* = 0.5, 0.75, and 0.95 are shown in Fig. A 1 where the coordinate is transformed back to the ( r z ) coordinates. The PAGE 212 212 shrink of the contact area during the decompression process renders a decrease in the length scale of the transient contact radius, which makes the thermal diffusion from the contact zone dominant. During the whole decompression process the diffusion length in the r direction is still on the order of ()cO Thus for Eq. (216) to be valid, it requires that ciRR ( i = 1, 2), and (A 2) ici ( i = 1, 2) (A 3) In the context of elastic collision, (A 2) must be satisfied. Since (A 3) is identical to (A 1), the second condition becomes identical to the first one. iii) The spherical surface just outside the contact area can be described using the Taylor series expansion as (Timoshenko et al. 1970, Johnson 1987) 2~...2irzR ( i = 1, 2) (A 4) Thus for the spherical surfaces near the contact area and within the region of the thermal influence by the contact to be approximated by z = 0 (Eq. (2 14)), it is necessary that 22ciiRRR and 2iciiRR ( i = 1, 2) (A 5) based on (A 2) and (A 3). Th us when the assumptions of the elastic collision (A 2) is satisfied, the thermal conditions in (A 1) and (A 5) simply reduce to ici ( i = 1, 2) (A 6) Recalling the definition of Fo in Eq. (2 23a), condition (A 6) can be rewrit ten as PAGE 213 213 1/21c iR Fo R ( i = 1, 2) (A 7) or 1/5 2 1/2 2/5 1212 12 125 1 1 4imR Fo v RE (A 8) For spherical particles of the same size and same material, condition (A 8) can be simplified to 2/5 2 12E Fo (A 9) Typically E ~ O ( 1010) (N/m2) or higher and ~ O (103) (kg/m3). For particle collision velocity at 10 (m/s), the above condition requires 2/5 2 12E Fo ~ 100 for the semi infinite media assumption to be valid. For the physical parameters used by Zhou et al (2008) where = 1451.7 kg/m3, E = 193GPa, and the normal relative velocity v12 = 0.5 m/s, the condition (A 9) requires 3090 Fo which is satisfied in all the simulations described in Zhou et al.s work. In conclusion, for practical engineering problems where conditi on (A 8) is satisfied, the semi infinite media assumption is appropriate and the results of the present work can be applied. In the extreme cases where (A 8) does not hold or plastic deformation occurs, one must adopt the approach used by Zhou et al. (2008) to capture the influence of the geometry and the relevant boundary conditions. PAGE 214 214 ( a ) ( b ) ( c) Figure A 1. Isothermals on the surface of particle 1 at Fo = 50 when (a) t = 0.50, (b) t = 0.75, and (c) t = 0.95. PAGE 215 215 APPENDIX B ASYMPTOTIC ANALYSIS OF SECONDORDER ACCURATE BOUNDARY SCHEMES The asymptotic analysis and diffusive scaling in ( Yoshida and Nagaoka 2010) are applied here with the following dimensionless parameter introduced = / (B 1) where L is the reference length. For details of the asymptotic analysis of the thermal lattice Boltzmann equation (TLBE) and the diffusive scaling please refer to the work by Yoshida and Nagaoka (2010). Here we focus on the asymptotic analysis for the present boundary condition treatments; some results given by Yoshida and Nagaoka (2010) are used directly. The distribution function g and the macroscopic variable can be expanded as (the arguments of ( x t ) are dropped for simplicity): g = (0)g + (1)g + 2 (2)g + O ( 3), (B 2) = (0) + (1) + 2 (2) + O ( 3), (B 3) and the spatial vectors xf and xff in the field are related to the boundary vector xw by (see Fig. 32) xf = xw e = xw e (B 4a) xff = xw e (B 4b) Note that the spatial vector x and time t throughout this Appendix are dimensionless variables normalized by the reference length and diffusion time scales, respectively. Based on the diffusive scaling in (Yoshida and Nagaoka 2010) the present boundary condition treatments in Eqs. (332) and (338) can be rewritten in terms of xw as follows Dirich let boundary condition: PAGE 216 216 g ( xw e t + 2) = cd 1 g ( xw e t ) + cd 2 g ( xw e t ) + cd 3 g ( xw e t ) + cd 4 d and (B 5) Neumann boundary condition: g ( xw e t + 2) = cn 1 g ( xw e t ) + cn 2 g ( xw e t ) + cn 3 g ( xw e t ) + cn 4 n / U (B 6) where U is a reference speed ( U = / ) for diffusive scaling as defined in (Yoshida and Nagaoka 2010) The source term does not affect the boundary condition treatment; thus it is not considered here for simplicity. The post collision distribution functions in Eqs. ( B 5, B 6) are manipulated with the collision operator, and after inserting the expansion ( B 2) into Eqs. ( B 5, B 6) each distribution can be expressed with respect to the boundary distribution function, g(m) ( xw, t ), or g(m) ( xw, t ) (m = 0, 1, and 2) using Taylor series expansions as below g(m) ( xw e t + 2) = g(m) e jx + 2 t + 12 22e e 2 ijxx + O ( 3)) g(m) (B 7a) g(m) ( xw e t ) = g(m) ( xw e t ) + L g(m) ( xw e t ) = ( g(m) + L g(m) ) + ( e jx + 12 22e e 2ijxx + O ( 3)) ( g(m) + L g(m) ), (B 7b) g(m) ( xw e t ) = g(m) ( xw e t ) + L g(m) ( xw e t ) = ( g(m) + L g(m) ) + ( e jx + 12 22 e e 2 ijxx + O ( 3)) ( g(m) + L g(m) ), (B 7c) g(m) ( xw e t ) = g(m) ( xw e t ) + L g(m) ( xw e t ) PAGE 217 217 = ( g(m) +L g(m) jx + 12 22 e e 2ijxx + O ( 3)) ( g(m) +L g(m) ), (B 7d) where the arguments ( xw, t ) on the RHS of the equations are skipped for brevity and the collision operator is L (m)g = M1SM ( Q(0) + Q(1)) (m)g (B 8) with Q(0) = T016(,,...,)(1,1,...,1) I7X7, and (B 9a) Q(1) = uT1234561(0,,,,,,)](1,1,1,1,1,1) (B 9 b) Substitution of Eqs. (B 7a) through (B 7d) into the boundary schemes given by Eqs. (B 5) and (B 6) and rearranging the relevant terms according to the power of one obtains the required conditions to satisfy the thermal boundary conditions to secondo rder accuracy. Dirichlet boundary condition: The equation of zeroth order in is (0) (0) 3124(1)() d cgccgc (B 10) In the power of 1, the equation is (1) (1)312(1)()dcgccg + (0) (0)33 122[(dd jjggcce cccexx (B 11) According to Yoshida and Nagaoka (2010) the following relationships are noticed based on the respective Eqs. (50), (51) and (57) (0)(0)0gg (B 12) PAGE 218 218 (0)(0)(0)gg (B 13) (1)(1)(1) gg (B 14) Equations ( B 12 B 14) imply th at the following relationships must hold for the coefficients in Eqs. ( B 10) and ( B 11) ddd dd dd ddd dccc cc ccccc c312 34 33 122 41() 1 (1)(1)() 0, ( B 15) in order to satify (0) d and (B 1 6 ) (1) 0 (B 1 7 ) Noting that the Dirichlet thermal boundary condition is given by d (3 2) and recalling the expansion (0) + (1) + 2 (2) + O ( 3) in (B 3) the following is thus obtained (0) d + O ( 2). (B 1 8 ) Equation ( B 18) verifies that the Dirichlet boundary condition is satisfied up to secondorder accuracy with the leading order solution sol ved by (0) in TLBE simulations. The coefficients in Eq. ( B 15) cannot be uniquely determined, and leaving cd 1 adjustable the constraints can be rewritten as PAGE 219 219 dddddddccccccc12131412222121. ( B 19) Neumann boundary condition: The Neum ann boundary condition is given by Eq. (33). For brevity, we assume the normal direction of the physical boundary is in the same direction as the discrete velocity direction e ; thus n can be directly applied. Substitution of expansions in (B 7 a B 7d) into the boundary scheme in (B 6) gives the equation of zerothorder in as ncgccg(0) (0)312(1)()0 (B 20) In the power of 1, the equation is (1) (1)312(1)()ncgccg nnn jjggcce cccecx xU(0) (0)33 122 4[( (B 21) and proceeding further to O ( 2) yields n jjgg cgccgceccce xx(1) (1) (2) (2) 312 3 112(1)() n ij ijgg cee ccee xx xx2(0) 2(0) 22 22 3 1211 22 (B 2 2 ) With the expansion (0) + (1) + 2 (2) + O ( 3) in (B 3) one would expect that PAGE 220 220 (0)(0)ij i njDvx (B 2 3 ) so that the leading order solution of (0) satisfies the Neumann boundary condition. Recalling the relationship between ijD and ij in Eq. (3 14) the above (B 23) can be rewritten as (0)(0)1()2inijijjvxUU (B 2 4 ) Besides Eqs. (B 12 B 14), the following relationships can be obtained from Eqs. (56) and (64) in Yoshida and Nagaoka (2010), respectively (0) (1)(1)(0) i jv gg Ux (B 2 5 ) (1)(2)(2) (1)i jvggU (B 2 6 ) Comparing Eq. ( B 20) with Eq. ( B 12) and combi ni ng Eqs. ( B 13, B 21, B 24, and B 25) we obtain the requirements nnn nn nnn n nn nnccc ccccc ccc cc312 33 122 3 33 431 (1)(1)() 12(1) 1. ( B 27) The unique solution to the coefficients in Eq. ( B 27) is n n n nc c c c1 2 3 41 2 2 2 2 2 2 ( B 28) PAGE 221 221 With the above coefficients and the relationships in Eqs. ( B 14) and ( B 26) the following can be derived from Eq. ( B 22) iijijDjvxU(1)(1)1()02 (B 2 9 ) As discussed in (Yoshida and Nagaoka 2010) since the condition f or (1) is linear and homogeneous, (1) in ( B 29) has only a trivial solution of (1)0 Thus we have verified that the Neumann boundary condition is satisfied to secondorder accuracy. PAGE 222 222 A PPENDIX C EXTENSION OF ZHOUS BGK AXISYMMETRIC LATTICE BOLTZMANN MODEL (ZHOU 2011) TO AN MRT LATTICE BOLTZMANN MODEL The governing equations for the athermal incompressible axisymmetric flows in the cylindrical coordinate system can be written as jrjuuxr (C 1 ) ii iiij irjijuuuuputx2221 (C 2 ) where i j are the indices denoting the radial r and axial z directions, xi and ui are the components of the spatial vector x ( r z ) and the velocity vector u ( ur, uz), respectively, is the density, p is the pressure, and ir is the Kroneckers delta function. Substitution of the continuity equation ( C 1) into Eq. ( C 2) gives the following momentum equation (Li et al. 2010, Zhou 2011) ji i iririij irj ijji iuu u uuuupuutx221 (C 3 ) For incompressible flows in the Cartesian coordinate system, the continuity and momentum equations recovered by the standard LBE read jjux0 (C 4 ) jiiijj ijjiuuupuutx1 (C 5 ) PAGE 223 223 As implied in Zhous LB model (2011) Eq. ( C 3) can be co nsidered as a particular form of Eq. ( C 5) with explicit force terms ir iuu rrx and iriiruurr22 and another implicit source term accounting for the continuity Eq. ( C 1). The LBE for the velocity distribution function f in Zhous BGK LB model (2011) is formulated as f eq2211,,12xe x (C 6 ) in order to recover the axisymmetric flow equations. For brevity, the special case of r = 0 discussed in (Zhou 2011) is not considered here. The source or sink term is r/ (C 7 ) and Fi is a force term defined by irii irFrr22 (C 8 ) For rotational flows with azimuthal velocity component u, the force term in Eq. ( C 8) should be modified to irii ir irFrrr222 (C 9 ) Here we rearrange Eq. ( C 6) as f eq eq21121,,2xe x (C 10) It is noted in Eq. ( C 10) that the fist bracketed term on the right side, ffeq1 is the standard BGK collision operator in the LBE method. The second bracketed term, PAGE 224 224 ereq1212 and the last term, eF 2 correspond to the velocity gradient (stress) term ir iuu rrx and the force term iriiruurr22 in Eq. ( C 3), respectively. As discussed earli er, the source term, included in Eq. ( C 10) is necessary to satisfy the continuity equation ( C 1). The proof for the correspondence of each term in the LBE ( C 10) to their relevant terms recovered in the momentum Eq. ( C 3) can be found in (Zhou 2011) I n this manner, the velocity gradient can be considered as a source term and obtained from the nonequilibrium part of the velocity distribution functions as in Eq. ( C 10). This is analogous to the treatment of the scalar gradient, rrDrr in Eq. ( 5 2) as a source term obtained from the nonequilibrium distribution functions, m rr egg r eq 111 1 2 in Eq. ( 5 6) in the present LB model for the axisymmetric CDE. The BGK model by Zhou (2011) can then be conveniently extended to an MR T model with the following evolution equation for f f eq1,,MSxe xmm (C 11) where M is a transformation matrix converting the velocity distribution functions f to their moment m by m = M f m(eq) is the equilibria of the moments, and S is a matrix of relaxation coefficients. According to Eq. ( C 10), the combined source term G in Eq. ( C 11) is determined from PAGE 225 225 Geffr eq2121 12 (C 12) The standard D2Q9 MRT LB model in (Lallemand and Luo 2000) and (Yu et al. 2003) can be applied to s imulate Eq. ( C 11) and it is used in the present work. PAGE 226 226 LIST OF REFERENCES Abanades S., Charvin P., Flamant, G., and Neveu, P., 2006, Screening of W ater S plitting T hermochemical C ycles P otentially A ttractive for H ydrogen P roduction by C oncentrated S ola r E nergy, Energy, 31, pp. 28052822. Aidun, C. K., and Clausen, J. R., 2010, Lattice Boltzmann M ethod for C omplex F lows, Annu. Rev. Fluid Mech., 42 pp. 439472. Alexander F. J., Chen S., and Sterling J. D., 1993, Lattice Boltzmann T hermohydrodynam ics, Phys. Rev. E, 47 R2249. Ammar, F. B., Kaviany M., and Barber J. R., 1992, Heat T ransfer D uring I mpact, Int. J. Heat Mass Transfer, 35, pp. 14951506. Asinari P. 2008, Asymptotic A nalysis of M ultiple relaxationtime L attice Boltzmann S chemes for M ixture M odeling , Comput. Math. Appl. 55 pp. 13921407. Barber J. R. 1989, An A symptotic S olution for S hort time T ransient H eat C onduction between T wo S imilar C ontacting B odies , Int. J. Heat Mass Transfer, 32, pp. 943 949. Behrend, O., 1995, S olid F luid B oundaries in P article S uspension S imulations via the L attice Boltzmann M ethod, Phys. Rev. E 52 pp. 11641175. Bhatnagar P. L., Gross E. P., and Krook M. 1954, A M odel for C ollision P rocesses in G ases. I. S mall A mplitude P rocesses in C ha rged and N eutral O ne C omponent S ystems , Phys Rev 94, pp. 511 525. Bouzidi M., Firdaouss M., and Lallemand P. 2001, Momentum T ransfer of a BoltzmannL attice F luid with B oundaries , Phys Fluids, 13, pp. 3452 3459. Chen, H., Chen S., and Matthaeus W. H. 1992, Recovery of the Navier Stokes E quations U sing a L attice G as Boltzmann M ethod , Phys Rev A, 45 R5339. Chen, L., Kang, Q., He, Y. L., and Tao, W. Q. 2012, PoreS cale S imulation of C oupled M ultiple P hysicochemical T hermal P rocesses in M icro R eactor for H ydrogen P roduction U sing L attice Boltzmann M ethod, Int. J. Hydrogen Energ. 37, pp. 1394313957. Chen, S., and Doolen, G. D. 1998, Lattice Boltzmann M ethod for F luid F lows , Annu. Rev. Fluid Mech., 30, pp. 329364. Chen, S., Martinez D., and Mei R. 1996, On B oundary C onditions in L attice Boltzmann M ethod, Phys Fluids, 8 pp. 2527 2536. Chen, S. Tolke, J. Geller, S. and Krafczyk, M. 2008, Lattice Boltzmann M odel for I mcompressible A xisymmetric F lows, Phys. Rev. E 78, 046703. PAGE 227 227 Che n, S., Tolke, J. and Krafczyk, M. 2009, Simulation of B uoyancy D riven F lows in a V ertical C ylinder U sing a S imple L attice Boltzmann M odel, Phys. Rev. E 79 016704. Chen, Y., Ohashi H., and Akiyama, M. 1994, Thermal L attice Bhatnagar Gross Krook M odel without N onlinear D eviations in M acrodynamic E quations , Phys Rev E, 50, 2776. de Vahl Davis, G. and Thomas, R. W. 1969, Natural C onvection between C oncentric V ertical C ylinders, Phys. Fluids 12, (suppl. 2) pp. 198 207. dHumires D. 1992, Gen eralized L attice Boltzmann E quations. In R arefied G as D ynamics: T heory and S imulations , Prog. Astronaut. Aeronaut. 159, pp. 450 458. dHumires D., Ginzburg I., Krafczyk M., Lallemand, P., and Luo, L. S. 2002, Multiple relaxtiontime L attice Boltzma nn M odels in T hree D imensions, Philos. T. Roy. Soc. A, 360 pp. 437 451. DOrazio A., Succi S., and Arrighettic C. 2003, Lattice Boltzmann S imulation of O pen F lows with H eat T ransfer , Phys Fluids, 15 pp. 27782781. Faghri, A., and Zhang, Y., 2006, Transport Phenomena in Multiphase Systems Elsevier Inc. Filippova O., and Hnel D. 1998, Grid R efinement for L attice BGK M odel , J. Comp ut. Phys 147 pp. 219 228. Flekky, E. G., 1993, Lattice Bhatnagar Gross Krook M odels for M iscible F luids, Phys Rev. E 47, pp. 42474257. Garimella S. V., and Sobhan, C. B. 2003, Transport in M icrochannels a C ritical R eview , Annu. Rev Heat Transfer, 13, pp. 1 50. Ginzbourg I., and dHumires D. 1996, Local S econdorder B oundary M ethods for L attice Bolt zmann M odels , J. Stat Phys 84 pp. 927 971. Ginzburg I., and dHumires D. 2003, Multireflection B oundary C onditions for L attice Boltzmann M odels , Phys Rev E, 68, 066614. Ginzburg I., 2005, Generic B oundary C onditions for L attice Boltzmann M od els and T heir A pplication to A dvection and A nisotropic D ispersion E quations , Adv Water Resour. 28, pp. 11961216. Graebel, W. P. 2001, Engineering Fluid Mechanics Taylor & Francis, New York. PAGE 228 228 Guo, Z., Han, H. Shi, B. and Zheng, C., 2009, Theory of t he L attice Boltzmann E quation: L attice Boltzmann M odel for A xisymmetric F lows, Phys. Rev. E 79, 046708. Guo Z., and Zhao T. S. 2002, Lattice Boltzmann M odel for I ncompressible F lows through P orous M edia , Phys Rev E, 66, 036304. Guo Z., Zhao T. S ., and Shi Y. 2006, Physical S ymmetry, S patial A ccuracy, and R elaxation T ime of the L attice Boltzmann E quation for M icrogas F lows , J. Appl Phys 99, 074903. Guo Z., Zheng C., and Shi B. 2002, An E xtrapolation M ethod for B oundary C onditions in L a ttice Boltzmann M ethod , Phys Fluids, 14, pp. 2007 2010. Guo Z., Zheng C., and Shi B. 2007, Thermal L attice Boltzmann E quation for L ow Mach N umber F lows: D ecoupling M odel , Phys Rev E, 75, 036704. Hahn, D. W., and Ozisik M. N., 2012, Heat Conducti on, 3rd ed. John Wiley & Sons, New York USA Halliday, I. Hammond, L. A. Care, C. M. Good, K. and Stevens, A. 2001, Lattice Boltzmann E quation H ydrodynamics, Phys. Rev. E 64, 011208. He X., Chen S., and Doolen, G. D. 1998, A N ovel T hermal M odel for the L attice Boltzmann M ethod in I ncompressible L imit , J. Comp ut. Phys. 146 pp. 282300. He, X., and Doolen, G., 1997, Lattice Boltzmann M ethod on C urvilinear C oordinates S ystem: Flow around a C ircular C ylinder, J. Comput. Phys. 134 pp. 306 31 5. He X., Zou Q., Luo L. S., and Dembo M. 1997, Analytic S olutions of S imple F lows and A nalysis of N onslip B oundary C onditions for the L attice Boltzmann BGK M odel , J. S tat. Phys 87, pp. 115 136. Heasley, J.H. 1965, Transient H eat F low between C o ntacting S olids , Int. J. Heat Mass Transfer, 8 pp. 147154. Huang, H. and Lu, X. Y., 2009, Theoretical and N umerical S tudy of A xisymmetric L attice Boltzmann M odels, Phys. Rev. E 80 016701. Ho C. F., Chang C., Lin K. H., and Lin, C. A. 2009, Con sistent B oundary C onditions for 2D and 3D L aminar L attice Boltzmann S imulations , Comput. Model Eng. Sci. 44 137. Howard, J. R., and Sutton, A. E. 1970, An A nalogue S tudy of H eat T ransfer through P eriodically C ontacting S urfaces , Int. J. H eat Mass Tra nsfer, 3 pp. 173183. Huang H., Lee T. S., and Shu S. 2006, Thermal C urved B oundary T reatment for the T hermal L attice Boltzmann E quation , Int. J. Mod. Phys. C 17, pp. 631643. PAGE 229 229 Huang H., Lee T. S., and Shu S. 2007, Hybrid L attice Boltzmann F in ite D ifference S imulation of A xisymmetric S wirling and R otating F lows, Int. J. Numer. Meth. Fluids 53 pp. 1707 1726. Iwatsu, R., 2004, Flow P attern and He at T ransfer of S wirling Flows in Cylindrical C ontainer with R otating T op and S table T emperature G r adient, Int. J. Heat Mass Transfer 47 pp. 27552767. Jeong H. K., Yoon H. S., Ha M. Y., and Tsutahara M. 2010, An I mmersed B oundary T hermal L attice Boltzmann M ethod U sing an E quilibrium I nternal E nergy D ensty A pproach for the S imulation of F lows with H eat T ransfer , J. Comput. Phys 229 pp. 25262543. Johnson, K. L. 1987, Contact Mechanics Cambridge University Press Cambridge, U.K. Kang, Q., Lichtner, P. C., and Zhang, D., 2006, Lattice Boltzmann PoreScale Model for Multicomponent Reactive Transport in Porous Media, J. Geophys. Res., 111 B05203. Kang, Q., Lichtner, P. C., and Zhang, D., 2007, An Improved Lattice Boltzmann Model for Multicomponent Reactive Transport in Porous Media at the Pore Scale, Water Resour. Res., 43 W12S14. Kao, P H., Ren, T. F., and Yang, R. J., 2007, An I nvestigation into F ixed B ed M icroreactors U sing L attice Boltzmann M ethod S imulations, Int. J. Heat Mass Transfer 50 pp. 42434255. Kaviany M. 1988, Effect of a M oving P article on W all H eat T ransfer in a C hannel F low , Nume r. Heat Transfer, 12, pp. 111 124. Kaviany M. 1998, Heat T ransfer in P orous M edia. Handbook of Heat Transport 3rd ed., McGraw Hill, pp. 9.1 9.82. Klausner, J. F., Petrasch, J., Hahn, D. W., and Mei, R., 2011, Solar thermochemical fuel production via a novel low pressure, magnetically stabilized, nonvolatile iron oxide looping process, Research Proposal for ARPAE DOE Kim, B. S., Lee, D. S., Ha, M. Y., and Yoon, H. S., 2008, A N umerical S tudy of N atural C onvection in a S quare E nclo sure with a C ircular C ylinder at D ifferent Vertical L ocations, Int. J. Heat Mass T ransfer, 51 pp. 18881906. Kodama, T., and Gokon N. 2007, Thermochemical C ycles for H igh T emperature S olar H ydrogen P roduction , Chem Rev 107 pp. 40484077. Kostoglou M., and Konstandopoulos A. G., 2002, On the T ransient H eat T ransfer through the S olid P hase in A ggregates and D eposits , Ind. Eng. Chem. Res. 41 pp. 33173325. PAGE 230 230 Kumar, R. and Kalam, M. A. 1991, Laminar T hermal C onvection between V ertical C oaxial I sothermal C ylinders, Int. J. Heat Mass Transfer 34, pp. 513 524. Kuo L. S., and Chen, P. H. 2009, Numerical I mplementation of T hermal B oundary C onditions in the L attice Boltzmann M ethod , Int. J. Heat Mass Transfer 52 pp. 529532. Lallemand, P., and Luo L. S. 2000, T heory of the L attice Boltzmann M ethod: D ispersion, D issipation, I sotropy, Galilean I nvariance, and S tability , Phys Rev. E, 61 pp. 6546 6562. Lallemand, P., and Luo, L. S. 2003, Theory of the L attice Boltzmann M ethod: Acoustic and T hermal P roperties in T wo and T hree D imensions , Phys Rev. E, 68 036706. Ltt J., Chopard B., Succi S., and Toschi F. 2006, Numerical A nalysis of the A veraged F low F ield in a T urbulent L attice Boltzmann S imulation , Physica A, 362, pp. 6 10. Lee, T S. Huang, H. and Shu, C., 2006, An A xisymmetric Incompressible Lattice Boltzmann M odel for P ipe F low, Int. J. Mod. Phys. C 17, pp. 645 661. Li J., and Mason D. J. 2000, A C omputational I nvestigation of T ransient H eat T ransfer in P neumatic T ransport of G ranular P articles , Powder Technol 112 pp. 273282. Li L., Mei R., Klausner J. F., and Hahn D. W. 2012, Heat T ransfer between C olliding S urfaces and P articles , ASME J. Heat Trans 134 011301. Li L., Mei R., and Klausner J. F. 2013a, Boundary C onditions for T hermal L attice Boltzmann E quation M ethod, J. Comp ut. Phys 237 pp. 366 395. Li L., Mei R., and Klausner J. F. 2013b, Heat T ransfer E valuation on C urved B oundaries in T hermal L attice Boltzmann E quation M ethod, ASME J. He at Trans accepted for publication. Li L., Mei R., and Klausner J. F. 2013c, Multiple relaxationtime L attice Boltzmann M odel for the A xisymmetric C onvection D iffusion E quation , Int. J. Heat Mass Transfer, under review. Li L., Singh, A., AuYeung, N ., Mei R., Petrasch, J., and Klausner J. F. 2013d, L attice Boltzmann Simulation of High Diffusivity Problems with Application to Energy Transport in a HighTemperature Solar Thermochemical Reactor, Proceedings of the ASME 2013 Summer Heat Transfer Con ference, HT2013, July 1419, Minneapolis, MN Li, Q., He, Y. L., Tang, G. H., and Tao, W. Q., 2009, Lattice Boltzmann M odel for A xisymmetric T hermal F lows, Phys. Rev. E 80, 037702. PAGE 231 231 Li, Q., He, Y. L., Tang, G. H., and Tao, W. Q., 2010, Improved A xisymme tric L attice Boltzmann S cheme, Phys. Rev. E 81, 056707. Liu C. H., Lin K. H., Mai H. C., and Lin, C. A. 2010, Thermal B oundary C onditions for T hermal L attice Boltzmann S imulations , Comput. Math. Appl. 59 pp. 21782193. Loutzenhiser P. G., Meier A., and Steinfeld, A., 2010, Review of the T wo step H2O/CO2S plitting S olar T hermochemical C ycle B ased on Zn/ZnO R edox R eactions , Materials, 3 pp. 4922 4938. Luo, L. S., and Girimaji, S. S. 2003, Theory of the L attice Boltzmann M ethod: T wo fluid M odel for B inary M ixtures, Phys. Rev. E 67 036302. Maier, R. S., and Bernard, R. S., 2010, Lattice Boltzmann A ccuracy in P ore S cale F low S imulation, J. Comp ut Phys. 229 pp. 233 255. Manjhi, N., Verma, N., Salem, K., and Mewes, D., 2006, Lattice Boltzm ann M odelling of U nsteady state 2D C oncentration P rofiles in A dsorption B ed, Chem. Eng. Sci. 61, pp. 25102521. Mehdizadeh, A. M., Klausner J. F., Barde A., and Mei R. 2012, Enhancement of T hermochemical H ydrogen P roduction U sing an I ron silica M agnet ically S tabilized P orous S tructure , Int J. Hy dro gen Ener g. 37 pp. 89548963. Mei R., Yu D., Shyy W. and Luo, L. S., 2002, Force E valuation in the L attice Boltzmann M ethod I nvolving C urved G eometry , Phys Rev E, 65, 041203. Mei R., Luo L. S., a nd Shyy W. 1999, An A ccurate C urved B oundary T reatment in the L attice Boltzmann M ethod , J. Comput. Phys 155, pp. 307 330. Mei R., Shyy W., Yu D., and Luo, L. S. 2000, Lattice Boltzmann M ethod for 3D F lows with C urved B oundary , J. Comp ut. Phys 161 pp. 680 699. Melchior T., and Steinfeld A. 2008, Radiative T ransfer within a C ylindrical C avity with D iffusely/ S pecularly R eflecting I nner W alls C ontaining an A rray of T ubular A bsorbers , ASME J. Sol. Energ. 130 021013. Mezrhab, A., Moussaoui M. A., Jami M., Naji H., and Bouzidi M. 2010, Double MRT T hermal L attice Boltzmann M ethod for S imulating C onvective F lows , Phys Lett A, 374 pp. 34993507. Modest M. F. 2003, Radiative Heat Transfer 2nd ed Academic Press, USA Moukalled, F., and Acharya, S., 1996, Natural C onvection in the A nnulus between C oncentric H orizontal C ircular and S quare C ylinders, J. Thermophys. Heat Transfer 10 pp. 524 531. PAGE 232 232 Mller ., 2008, Transient H eat T ransfer in a D irectly I rradiated S olar C hemical R eactor for the T hermal D issociation of ZnO , Appl Therm Eng 28, pp. 524 531. Noble D. R., Chen S., Georgiadis J. G., and Buckius R. O., 19 95, A C onsistent H ydrodynamic B oundary C ondition for the L attice Boltzmann M ethod, Phys Fluids, 7 pp. 203 209. Pan C., Luo L. S., and Miller C. T. 2006, An E valuation of L attice Boltzmann S chemes for P orous M edium F low S imulation , Comput Fluids, 35 pp. 898909. Peng, Y., Chew, Y. T., and Shu, C., 2003, Simplified T hermal L attice Boltzmann M odel for I ncompressible T hermal F lows, Phys. Rev. E 67, 026701. Peng, Y., Shu, C. Chew, Y. T., and Qiu, J., 2003, Numerical I nvestigation of F lows in Czo chralski C rystal G rowth by an A xisymmetric L attice Boltzmann M ethod, J. Comput. Phys. 186 pp. 295 307. Perkins C., and Weimer A. W. 2009, Solar T hermal P roduction of R enewable H ydrogen, AIChE J 55, pp. 286 293. Petrasch J. 2002, Thermal Modeli ng of Solar Chemical Reactors: Transient Behavior, Radiative Transfer , M asters thesis, ETH Zrich. Petrasch J., Coray P., Meier A., Brack M., Haberling P., Wuillemin D., and Steinfeld A. 2007, A N ovel 50 kW 11,000 S uns H igh flux S olar S imulator B ased on an A rray of Xenon A rc L amps , ASME J. Sol. Energ. 129, pp. 405 411. Premnath K. N., and Abraham J. 2007, Three dimensional M ultirelaxationtime (MRT) L attice Boltzmann M odels for M ultiphase F low , J. Comput. Phys., 224, pp. 539559. Qian Y. H., dHumires D., and Lallemand, P. 1992, Lattice BGK M odels for Navior Stokes E quation , Europhys Lett 17, pp. 479 484. Reis, T., and Phillips, T. N., 2007, Modified Lattice Boltzmann Model for Axisymmetric Flows, Phys. Rev. E, 75, 056703. Reis, T., and Phillips, T. N., 2008, Numerical Validation of a Consistent Axisymmetric Lattice Boltzmann Model, Phys. Rev. E, 77 026703. Rong F., Guo Z., Chai Z., and Shi B. 2010, A L attice Boltzmann M odel for A xisymmetric T hermal F lows through P orous M edia, Int. J. Heat Mass Transfer, 53, pp. 55195527. ServanCamas, B. and Tsai, F. T. C. 2008, Lattice Boltzmann M ethod with T wo R elaxation T imes for A dvectionD iffusion E quation: Third O rder A nalysis and S tability A nalysis, Adv. Water Resour. 31 pp 11131126. PAGE 233 233 Shan X. 1997, Solution of Rayleigh Benard C onvection U sing a L attice Boltzmann M ethod , Phys Rev E, 55, 2780. Shi B., and Guo, Z. 2009, Lattice Boltzmann M odel for N onlinear C onvectionD iffusion E quations , Phys. Rev. E 79 016701. Sh u, C., and Zhu, Y. D., 2002, Efficient C omputation of N atural C onvection in a C oncentric A nnulus between an O uter S quare C ylinder and an I nner C ircular C ylinder, Int. J. Numer. Meth. Fluids 38 pp. 429445. Singh, A., 2013, Multi Scale, Multi Physics M odeling of ThermoChemical Iron/Iron Oxide Cycle for Fuel Production, Doctoral dissertation, University of Florida. Steinfeld A ., 2002, Solar H ydrogen P roduction via a T wo step W ater splitting T hermochemical C ycle B ased on Zn/ZnO R edox R eactions , Int. J. Hydrogen Energ 27 pp. 611619. Steinfeld A. 2005, Solar T hermochemical P roduction of H ydrogen A R eview , Sol Energy, 78, pp. 603 615. Succi S. 2001, The L attice Boltzmann E quation for F luid D ynamics and B eyond Oxford University Press, New York. Sullivan, S. P., Sani, F. M., Johns, M. L., and Gladden, L. F., 2005, Simulation of P acked B ed R eactors U sing L attice Boltzmann M ethods, Chem. Eng. Sci. 60, pp. 34053418. Sun J., and Chen M. M. 1988, A T heoretical A nalysis of H eat T ransfer D ue t o P article I mpact , Int. J. Heat Mass Transfer 31 pp. 969 975. Tang, G. H., Tao W. Q., and He Y. L., 2005, Thermal B oundary C ondition for the T hermal L attice Boltzmann E quation , Phys. Rev. E, 72, 016703. Timoshenko, S. P., and Goodier J. N. 1970, T heory of Elasticity 3rd ed McGraw Hill, New York. van der Sman, R. G. M. Ernst, M. H. 2000, ConvectionD iffusion L attice Boltzmann S cheme for I rregular L attices, J. Comput. Phys. 160, pp. 766 782. Venkatachalappa, M. Sankar, M., and Natarajan, A. A., 2001, Natural C onvection in an A nnulus between T wo R otating V ertical C ylinders, Acta Mech. 147, pp. 173 196. Verma, N., Salem, K., and Mewes, D., 2007, Simulation of Microand MacroTransport in a Packed Bed of Porous Adsorbents by Lattice Boltzmann Methods, Chem. Eng. Sci. 62 pp. 36853698. PAGE 234 234 Verma, N., and Mewes, D., 2008, Simulation of Temperature Fields in a Narrow Tubular Adsorber by Thermal Lattice Boltzmann Methods, Chem. Eng. Sci. 63 pp. 42694279. Wang, L. Guo, Z., and Zheng, C., 2010, Multi relaxationtime L attice Boltzmann M odel for A xisymmetric F lows, Comput. Fluids 39 pp. 15421548. Wang, M., Wang, J., Pan, N., and Chen, S., 2007, Mesoscopic P redictions of the E ffective T hermal C onductivity for M icroscale R andom P orous M edia, Phys. Rev. E 75, 036702. Wieckert C., Meier A., and Steinfeld, A. 2003, Indirectly I rradiated S olar R eceiver R eactors for H igh T emperature T hermochemical P rocesses , ASME J. Sol. Energ. 125 pp. 120 123. Yin, X., Le, G., and Zhang, J., 2012, Mass and M omentum T ransfer across S olid F luid B oundaries in the L attice Boltzmann M ethod, Phys. Rev. E 86, 026701. Yoshida, H., and Nagaoka, M. 2010, Multiple relaxationtime L attice Boltzmann M odel for the C onvection and A nisotropic D iffusion E quation, J. Comput. Phys., 229 pp. 77747795. Yu D., Mei R., Luo L. S., and Shyy W. 2003, Viscous F low C omputations with the M ethod of L attice Boltzmann E quation , Prog Aero sp. Sci 39, pp. 329 367. Yu H., Luo L. S., and Girimaji, S. S. 2006, LES of T urb ulent S quare Jet F low U sing an MRT L attice Boltzmann M odel , Compu t. Fluids, 35 pp. 957965. Ziegler D. P. 1993, Boundary C onditions for L attice Boltzmann S imulations , J. Stat Phys 71, pp. 11711177. Zheng H. W., Shu C., and Chew Y. T. 2006, A L attice Boltzmann M odel for M ultiphase F lows with L arge D ensity R atio , J. Comput. Phys. 218 pp. 353371. Zheng, L. Shi, B. Guo, Z., and Zheng, C., 2010a, Lattice Boltzmann E quation for A xisymmetric T hermal F lows, Comput. Fluids 39, pp. 945 952. Zh eng, L., Guo, Z., Shi, B. and Zheng, C., 2010b, Kinetic T heory B ased LBE with V iscous D issipation and P ressure W ork for Axisymmetric T hermal F lows, J. Comput. Phys. 229 pp. 58435856. Zhou, J.G., 2008, Axisymmetric L attice Boltzmann M ethod, Phys. Re v. E 78, 036701. Zhou, J.G., 2011, Axisymmetric L attice Boltzmann M ethod Revised Phys. Rev. E 84 036704. Zhou, J. H., Yu A. B., and Horio, M. 2008, Finite E lement M odeling of the T ransient H eat C onduction between C olliding P article s, Chem Eng J. 139 pp. 510516. PAGE 235 235 Zou Q., and He X. 1997, On P ressure and V elocity B oundary C onditions for the L attice Boltzmann BGK M odel , Phys Fluids, 9 pp. 15911598. PAGE 236 236 BIOGRAPHICAL SKETCH Like Li was born in 1985 in Hunan, China He received his bachelors degree in mechanical engineering from the Central South University in Changsha, China in 2007. He was then accepted by the graduat e school of Beijing University of Aeronautics and Astronautics (BUAA) with a research assistantship. In 2010, he obtained his masters degree in mechanical engineering from BUAA and flew overseas to the United States for pur suing a doctoral degree. Like Li has been studying and working as a Ph.D. student in the Department of Mechanical and Aerospace Engineering at the University of Florida since spring 2010. His current research interests include computational modeli ng of heat transfer, momentum and mass transfer in fluidsolid systems. 