UFDC Home  myUFDC Home  Help 
Title Page  
Dedication  
Acknowledgement  
Table of Contents  
List of Tables  
List of Figures  
Abstract  
Introduction to multiscale computational...  
Continuum model: Navierstokes...  
Error assessment of the lattice...  
Lbe method for immiscible twophase...  
Summary and future work  
References  
Biographical sketch 



Table of Contents  
Title Page
Page i Page ii Dedication Page iii Acknowledgement Page iv Table of Contents Page v Page vi List of Tables Page vii List of Figures Page viii Page ix Page x Abstract Page xi Page xii Introduction to multiscale computational dynamics with interfaces Page 1 Page 2 Page 3 Page 4 Continuum model: Navierstokes equations with multiscales Page 5 Page 6 Page 7 Page 8 Page 9 Page 10 Page 11 Page 12 Page 13 Page 14 Page 15 Page 16 Page 17 Page 18 Page 19 Page 20 Page 21 Page 22 Page 23 Page 24 Page 25 Page 26 Page 27 Page 28 Page 29 Page 30 Page 31 Page 32 Page 33 Page 34 Page 35 Page 36 Page 37 Page 38 Page 39 Page 40 Page 41 Page 42 Page 43 Page 44 Page 45 Error assessment of the lattice Boltzmann method for variable viscosity flows Page 46 Page 47 Page 48 Page 49 Page 50 Page 51 Page 52 Page 53 Page 54 Page 55 Page 56 Page 57 Page 58 Page 59 Page 60 Page 61 Page 62 Page 63 Page 64 Page 65 Page 66 Page 67 Page 68 Page 69 Page 70 Page 71 Page 72 Lbe method for immiscible twophase flow computation Page 73 Page 74 Page 75 Page 76 Page 77 Page 78 Page 79 Page 80 Page 81 Page 82 Page 83 Page 84 Page 85 Page 86 Page 87 Page 88 Page 89 Page 90 Page 91 Page 92 Page 93 Page 94 Page 95 Page 96 Page 97 Page 98 Page 99 Page 100 Page 101 Page 102 Page 103 Page 104 Page 105 Page 106 Page 107 Page 108 Page 109 Page 110 Page 111 Page 112 Page 113 Page 114 Page 115 Page 116 Page 117 Page 118 Page 119 Page 120 Page 121 Page 122 Page 123 Page 124 Page 125 Page 126 Summary and future work Page 127 Page 128 References Page 129 Page 130 Page 131 Page 132 Page 133 Page 134 Page 135 Page 136 Page 137 Biographical sketch Page 138 

Full Text  
MULTISCALE COMPUTATIONAL FLUID DYNAMICS WITH INTERFACES By JIANGHUI CHAO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006 Copyright 2006 by Jianghui Chao To my parents, my wife and my son. ACKNOWLEDGMENTS I would like to express my sincere gratitude to Drs. Wei Shyy and Renwei Mei for providing me the opportunity and flexibility to perform this work. I can not thank them enough for being so patient and understanding for years and pushing me to learn more. I would like to thank Drs. Siddharth Thakur and Alireza Haghighat for agreeing to serve in my thesis committee. Since the surroundings dictate the quality of life in several ways, I thank my research group members for helping with academic aspects while providing memorable company for the past five years. I thank my wife who has been the best friend for more than a decade and has been extremely patient and understanding during all these years, and my son who has been giving me lots of happiness. Finally I thank my parents who have been behind me every step of the way providing their unconditional support. TABLE OF CONTENTS A C K N O W L E D G M E N T S ..............................................................................................................4 L IST O F T A B L E S ......................................................................................................... ........ .. 7 LIST OF FIGURES ............................................. ............ ...........................8 A B S T R A C T .......................................................................................................... ..................... 1 1 CHAPTER 1 INTRODUCTION TO MULTISCALE COMPUTATIONAL DYNAMICS WITH IN T E R F A C E S .......................................................................................................................... 1 1.1 Motivation ............................................................................. 1 1 .2 O bje ctiv e s ................................................................................................ ..................... 3 1.3 Structure of the D issertation .................... ...............................................................4...... 2 CONTINUUM MODEL: NAVIERSTOKES EQUATIONS WITH MULTISCALES ........ 5 2.1 Introduction to Computational Techniques for Solving NavierStokes Equations ............5 2.2 Incom pressible V iscous Flow Solvers........................................................... ...............8... 2 .2 .1 In tro d u ctio n ....................................................................................................... .. 8 2.2.2 A artificial Com pressible M ethod (A CM ) ........................................ ..................... 9 2.2.3 SemiImplicit Method for PressureLinked Equations (SIMPLE) ...........................9 2.3 Computational Issues on Multiscale Computations ............................. ..................... 12 2 .3 .1 Stiffn ess ........................................ ..... ............................................................ 12 2.3.2 Grid Requirements for Multiscale Problems................................................... 14 2.3.3 M ethods to R educe Stiffness ................................................................ ............... 15 2.4 A Special Computation Example: ThermoMEMS Computation Results....................19 2.4.1 Introduction to Thermal Anemometry for Fluid Velocity and Skin Friction M e a su rem en t............................................................................................................... 19 2.4.2 G governing E quations ..................................................................... ................ 23 2 .4 .3 N um erical Schem es .. ...................................................................... ................ 25 2.4.4 C om putational Stiffness ......................................... ........................ ................ 25 2.4.5 G eom etry and G rid L ayout....................................... ...................... ................ 25 2 .4 .6 R results and D iscu ssion .......................................... ......................... ................ 30 2.4.7 Sum m ary and C conclusion: ......................................... ...................................... 36 3 ERROR ASSESSMENT OF THE LATTICE BOLTZMANN METHOD FOR VARIABLE VISCO SITY FLOW S ............................................................. ................ 46 3.1 Introduction to Lattice Boltzmann Method and Wall Boundary Condition..................48 3.1.1 LBE BGK M ethod ... .. ................................. ......................................... 48 3.1.2 B oundary C conditions ...................................................................... ................ 5 1 3.2 Error Assessment of the LBE Method due to Variable Viscosity via a Fully D developed C channel F low ............................................................................... .............. 55 3.2.1 Fullydeveloped Laminar Channel Flow with Variable Viscosity......................55 3.2.2. The Lattice Boltzmann Equation Treatment ....................................................58 3 .2 .3 A sse ssm en t ............................................................................................................. 5 9 4 LBE METHOD FOR IMMISCIBLE TWOPHASE FLOW COMPUTATION ...................73 4.1 Overview of Immiscible TwoPhase Flow Computation ...................... ..................... 73 4.2 Literature Review on LBE Method for TwoPhase Flow Computation........................76 4.3 He et al.'s Isothermal LBE Model for TwoPhase Flow.............................................80 4.3.1 The Boltzmann Equation for NonIdea Fluids ..................................................80 4.3.2 Lattice Boltzmann Scheme for Multiphase Flow in the Near Incompressible L im it.................. ... .......................................................................... ......... 8 5 4.4 Code V alidation: RayleighTaylor Instability ............................................. ................ 86 4.4.1 Linear Analysis of RayleighTaylor Instability ................................................87 4.4.2 RayleighTaylor Instability: LBE Results......................................... ................ 88 4.4.2.1 A SingleM ode G row th R ate................................................... ................ 89 4.4.2.2 Flow Field at R e=2048 and At=0.5......................................... ................ 90 4.4.2.3 Flow Field at R e=256 and A t=0.5........................................... ................ 91 4.5 LeeLin's Implicit LBE Twophase Model and He et al.'s Model with Surface T en sion ............... ....... .. .... ......... ..................................................................... . 9 1 4.5.1 LeeLin's Scheme for Large Density Ratio.......................................................91 4.5.2 Diffusion of Lee & Lin's Implicit LBE model..................................................94 4.5.3 M odeling the Surface Tension........................................................... ................ 96 4.6.1 Surface T ension C alculation .................. ............................................ ................ 97 4.6.2 A Filterbased Technique for Surface Tension .................................................98 4.6.3 Volume Conservation and Mass Conservation ....... ................. ...................100 4 .7 N um erical Sim ulations ................................................. ............................................ 104 4.7.1 Stationary B ubble .............................. ............................................ 104 4.7.2 Capillary Wave .................... .. ........... ............................... 107 4 .7.3 R ising B ubble .............. .. .................. .................. ................... ...... .. ............ 108 4.8 Sum m ary and C conclusion. ................................................................... ............... 110 5 SUMMARY AND FUTURE WORK .......................................................... 127 5 .1 Su m m ary ....................................................................................................... ....... .. 12 7 5 .2 F u tu re W o rk ................................................................................................................. ... 12 8 L IST O F R E F E R E N C E S ....................................................... ................................................ 129 B IO G R A PH IC A L SK E T C H .................................................... ............................................. 138 LIST OF TABLES Table page 41 Effect of grid resolution on computed pressure drop for Laplace number 100. The data were taken after 100 time steps and the density and viscosity ratios were set to 100 and 10 respectively .... .. .............................. ......... ........ ............... 112 42 Effect of Laplace number on the ratio of numerical pressure drop to theoretical pressure drop: Bubble diameter is 40 lattices. The data were taken at non dimensional time = 100. Density ratio is 100 and dynamics viscosity ratio is 10.........112 43 Effect of density ratio on pressure drop: Bubble diameter is 40 lattice units. Viscosity ratio was set to 10 for Laplace number = 100 and the data were taken at non dimensional time = 100 ..... ................ .......... .......................... ... 112 44 Effect of viscosity ratio on pressure drop: Bubble diameter is 40 lattice units. Laplace number = 100 and density ratio is 100. The data were taken at non dim ensional tim e = 100 ............................................................. ...... .. ...... ..... 112 LIST OF FIGURES Figure page 21 Geometry of channel flow with solid substrate ............... ....................................37 22 R educeddom ain geom etry ...................................................................... ................ 37 23 Relative reductions of residuals of singlephase and conjugate cases with specified sensor tem perature .............. .................................... ........ ....8.... .38 24 Relative reductions of residuals of conjugate cases with specified sensor temperature solved by singlegrid and m ultigrid solvers ................................................ ................ 39 25 Shear stress comparisons for the coarse grid and refine grid computations...................40 26 Relative reductions of Pressure and temperature residuals on the refined grids ...............41 27 Tem perature contours w ith Gr = 0.5 ............................................................. ............... 42 28 Tem perature contours w ith Gr = 10 ............................................................. ................ 42 29 Tem perature contours w ith Gr = 100 ........................................................... ................ 42 210 Temperature distribution on a crosssection which originates from the middle point of the sensor to the top boundary of the reduceddomain.............................................43 211 u velocity profiles of conjugate cases on the crosssections.........................................44 212 W all shear stress distribution ........................................... .......................... ................ 45 213 Shear stress variation ........................................................................................... 45 31 Boundary nodes and their neighbors using the square lattice.......................................67 32 A boundary cell using the hexagonal (FHP) lattice (Noble et al. 1995).........................67 33 Absolute L2 norm errors of LBE with Noble's scheme ...............................................67 34 A 2D 9velocity lattice (D 2Q 9) m odel ......................................................... ................ 68 35 Absolute error of a fullydeveloped channel flow using Inamuro et. al.'s scheme ............68 36 Two set of viscosity distributions used in this study .................................... ................ 68 37 The exact velocity profiles of the channel flows with different boundary layer thicknesses due to different viscosity distributions. ..................................... ................ 69 38 Square lattice distribution in channel flow simulation .................................................69 39 Com prison of the LBE velocity profiles .................................................... ................ 70 310 Dependence of the relative L2 norm error on the lattice size h in the fullydeveloped channel flow w ith variable viscosity............................................................. ................ 71 311 Comparison of ELBE = uLBE with E = u, uFD for H=200, and 3 = 0 .000 5, 0 = 0 .0 102 ................................................. ............................................. 72 41 The growth rate a (measured in units of (g2 Iv)1/3) of a disturbance vs. its wave num bers k ............................................................................. ....... ... ..................... 113 42 The grow th rate plot ............. .. .................. .................. ......................... .. ............ .. 113 43 Evolution of the fluid interface from a single mode perturbation for At=0.5 and R e = 2 0 4 8 ................................................................................................ .................... 1 14 44 He et al's results of the evolution of the fluid interface from a single mode perturb action .......................... ............................................... 115 45 Density profiles across the bubble and spike fronts at three different time steps..........116 46 Evolution of the fluid interface from a single mode perturbation. ..............................117 47 Evolution of index function of a stationary droplet with zero velocity .........................117 48 Evolution of the capillary number with for the stationary bubble simulation with density ratio 2 and dynamic viscosity ratio 2............... .........................118 49 Density profiles and pressure profiles at t=0 and t=2000 of the stationary bubble sim u nation ...................................................................................................... ....... .. 1 18 410 Surface tension profile for a stationary bubble computed from HeChenZhang m ethod .................................................................................................. .. 119 411 Surface tension calculated by Kim's formulation...... ........................................ 119 412 Rising bubble with Eo=0.971, Mo=1.26e3, density ratio=100, viscosity ratio=10........119 413 The velocity profiles and stream line of the spherical bubble at t=20 in 412...............120 414 Time evolution of the volume (A) and velocity (B) of the rising bubble with volume correction in Figure 413B .......................................................................... ............... 120 415 The contours of space derivative of hydropressure and density of a rising bubble with Eo=0.971, M=1.26e3, Re=5.19, density ratio=100, viscosity ratio=10 .......................121 416 Stationary bubble computation. (a) Spurious currents of magnitude Ca=5.03 x 103, (b) comparison between computed and theoretical pressure jump............................... 121 417 Density profiles of the stationary bubble with diameter 40, density ratio 100, viscosity ratio 10, and L a= 100 ..................................... ........................ ............... 122 418 Maximum spurious velocities for these two grid resolutions..................................122 419 Initial interface profile for a capillary wave simulation........................ ...................123 420 Time evolution of the amplitude of a capillary wave with density ratio 100 ................123 421 Shape diagram of Clift et al. (1978)...... ........... ........... ...................... 123 422 Computed bubble shapes. (A) Cylindrical, (B) Ellipsoidal, (C) Dimpledellipsoidal.....124 423 Time evolutions of rising bubbles (A) Cylindrical, (B) Ellipsoidal, (C) Dimpled ellip soid al...................................................................................................... ......... 12 4 424 Time evolutions of bubble velocities ....... ........ ........ ..................... 124 425 Density profiles of rising bubbles (A) Cylindricall, (B) Ellipsoidal, (C) Dimpled ellip soid al...................................................................................................... ....... .. 12 5 426 Time evolutions of a rising bubbles with Eo=97.1, M=0.971 and Re=31.2................. 126 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MULTISCALE COMPUTATIONAL FLUID DYNAMICS WITH INTERFACES December 2006 Chair: Wei Shyy Cochair: Renwei Mei Major Department: Mechanical and Aerospace Engineering Fluid flow and heat transfer problems involving interfaces typically contain property jumps and length and time scale disparities, resulting in computational stiffness, demanding resolution requirements due to nonlinearity and multiple physical mechanisms. In this dissertation, both continuum and kinetic approaches are investigated to address the multiscale thermofluid problems, especially those involving interfaces. A continuum model based on the NavierStokes fluid model and the Fourier law for heat transfer was employed to investigate the conjugate heat transfer problem. The problem is motivated by recent advancement in the microelectro mechanical systems (MEMS) devices for shear stress measurement. Due to the length scale disparity and large solidfluid thermal conductivity ratio, a twolevel computation is used to examine the relevant physical mechanisms and their influences on wall shear stress. The substantial variations in transport properties between the fluid and solid phases and their interplay in regard to heat transfer and nearwall fluid flow structures are investigated. It is demonstrated that for the stateoftheart sensor design, the buoyancy effect can noticeably affect the accuracy of the shear stress measurement. A lattice Boltzmann equation (LBE) model derived from the kinetic consideration is investigated to address (i) error behavior due to variable viscosity, and (ii) the interfacial fluid dynamics of some twophase flow problems. It is shown that the boundary treatment error does not have a significant interaction with the truncation error associated with variable viscosity, and the LBE model closely matches the NavierStokes model for fluid flows with large viscosity variation. Next, an improved interface lattice Boltzmann model is developed for twophase interfacial fluid dynamics with large property variations. In order to suppress numerical instabilities associated with the presence of the interface, the following approaches were investigated: (i) a new surface tension formulation originated from the diffusion interface method was used to remove unphysical pressure wiggles across interfaces; (ii) a filter scheme was used in numerical gradient calculations to maintain monotonic property variations; (iii) a volume correction procedure was devised. The performance of the improved LBE model was evaluated using the RayleighTaylor instability problem, stationary bubble under force equilibrium, capillary waves, and rising bubbles. The computational results demonstrated that this LBE two phase model is more robust than those reported in the literature, capable of treating larger density ratio up to order 0(102 ), while confining the interface thickness within 56 grids. CHAPTER 1 INTRODUCTION TO MULTISCALE COMPUTATIONAL DYNAMICS WITH INTERFACES 1.1 Motivation Many fluid mechanics and heat transfer problems involve multiple scale phenomena, such as turbulence, multiphase flows, and conjugate heat transfer. The presence of the disparity in length, time, and velocity scales is caused by the presence of different competing mechanisms, such as convection, diffusion, chemical reaction, body forces, and surface tension. Furthermore, these mechanisms are often coupled and nonlinear (Shyy et al. 1997b). Mathematically, issues related to multiscale problems can be illustrated by a system of ordinary differential equations whose characteristic matrix has a large range of eigenvalues. These large range eigenvalues physically represent different growth and decaying rates of different competing mechanisms. Numerically they make computations very sensitive to the stability of the numerical schemes (Lambert 1980; Shyy 1994; Lomax et al. 2001; Chapra and Canale 2002; Deuflhard and Bornemann 2002). If the linear algebraic equations discretized from the original stiff ODEs are written as (Miranker 1981) (I hfA)= b (1.1) where I denotes a mdimensional identity matrix, /7 is a constant, h is integration step size, A is the matrix that has large range of eigen values, x is solution vector, and b is source term vectors that can include initial/boundary conditions and other non homogeneous sources, the condition number of matrix A is defined as (Chapra et al. 2002) Cond(A)=AL A 1 A A axn (1.2) where Amax and Amin are the smallest and largest eigen values of A, respectively. From Eq. (1.2), it can be seen that the numerical solution process can be illconditioned (Miranker 1981), which brings in computational difficulties in both stability and accuracy. These kinds of problems are thus called stiff problems. In order to understand the physical processes involving different scales, separate scaling procedures for each physical regime and the coupling need to be devised. However, different multiscale problems involve different mechanisms and thus require different treatments. Generally, two classes of multiscale problems exist. The first one is that the small scales are restricted in distinct regions in space. Thus, the problems involved in this class require the matching and patching of the solutions in large and small scale domains. In the second class, the disparate scales coexit over the entire flow field and interact with each other. Resolving or modeling such disparate scales is the task of this class of multiscale problems (Shyy et al. 1997a). At present, there are two major computational fluid dynamics solvers for modeling incompressible, viscous fluid flows. One is the wellestablished NavierStokes equations solver. Another is the Lattice Boltzmann equation (LBE) solver (He et al. 1997; Chen et al. 1998; Succi 2001; Yu et al. 2003a). For multiscale problems in fluid dynamics and heat transfer, if the continuity assumption is valid, all scales can be represented by NavierStokes equations. The issue in this type of problem is how to resolve all the relevant scales based on NavierStokes equations. As an alternative computational fluid dynamics solver, the LBE method offers a mesoscale framework for fluid dynamics computations. It recovers the macroscopic fluid flow solution based on averaging of the particle distribution functions, obtained by solving the simplified form of the Boltzmann equation. Despite some difficulties in compressible flows and heat transfer, the LBE method has been successfully used for isothermal complex flows. For interfacial fluid dynamics, however, the numerical stability issues still exist from the large property jumps and scale disparities in normal and tangential directions of interfaces. A study on how to stabilize the LBE twophase computations for flows with large density ratio and prevent interface thickness from diffusion is still needed. In this dissertation, we will focus on the interfacial problems solved by Navier Stokes equations and lattice Boltzmann equation since interfaces always result in multi scale processes (Shyy 1994; Shyy et al. 1997c) due to the sharp changes in the material properties across the interface. Because of the physical nature and numerical issues of multiscale problems, they are very challenging and nontrivial. 1.2 Objectives One of the main objectives of the present research is to develop a twolevel computational method based on the NavierStokes equations solver to simulate a fluid flow and heat transfer problem involving a fluidsolid interface. Another main objective is to simulate a liquidgas twophase flow problem by using the lattice Boltzmann equation (LBE) method. Specifically, the objectives are listed below: To investigate the multiscale issues involved in the fluid flow and heat transfer surrounding a MEMSbased thermal shear stress sensor; To develop an efficient, accurate numerical method to simulate the fluid flow and heat transfer surrounding the MEMSbased thermal shear stress sensor; To investigate the performance of the lattice Boltzmann equation method when it is applied to fluid flows with rapidly varying viscosity over a thin boundary layer that is a typical multiscale phenomenon; To develop a lattice Boltzmann method for immiscible liquidgas twophase flows with large density ratio. 1.3 Structure of the Dissertation The numerical techniques for solving the NavierStokes equations are briefly presented in Chapter 2. The multiscale and stiffness issues involved in the models represented by NavierStokes equations are presented. Finally, in Chapter 2, a twolevel numerical technique is developed to simulate a multiscale physical problemthe fluid flow and conjugate heat transfer surrounding a MEMSbased thermal shear stress sensor. In Chapter 3 the method of Lattice Boltzmann equation (LBE) is briefly described, followed by the error assessment of the LBE method for variable viscosity. In Chapter 4, an improved interface lattice Boltzmann model is developed for two phase interfacial fluid dynamics with large property variations. Some numerical approaches are investigated to suppress numerical instability in the computations for large density ratio flows and to prevent interface thickness from diffusion. This dissertation ends with Chapter 5, which contains a summary and conclusion of the contributions of the present work, and the recommendations for future research efforts. CHAPTER 2 CONTINUUM MODEL: NAVIERSTOKES EQUATIONS WITH MULTISCALES 2.1 Introduction to Computational Techniques for Solving NavierStokes Equations Before addressing the multiscale flows described by the NavierStokes equations, it is necessary to introduce the different numerical methods for fluid flow and heat transfer since different governing equation features require different numerical methods, and different numerical methods have different advantages on different flow and heat transfer problems. Generally, NavierStokes equations solvers can be classified into two groups: one is the pressurebased formulation, and another is the densitybased formulation. The pressurebased formulation is often used for flows without discontinuities, such as incompressible flows without shock waves. In contrast, the densitybased formulation is often used for flows with discontinuities, such as supersonic/hypersonic flows (Chung 2002). The algorithms for compressible and incompressible flows are different, which can be understood from the role of pressure in compressible and incompressible flows (Moukalled et al. 2001). If fluid compressibility is zero (in other words, Mach number is zero), the thermodynamic pressure loses its sense for incompressible flows. Pressure gradient becomes a balance force to viscous, inertial and other body/surface forces. Therefore, no state equations exist for pressure. In addition, because density does not change locally for incompressible flows, continuity equation can no longer be considered as the governing equation for density. Thus no explicit governing equation for pressure exists for incompressible flow. Also, for incompressible flow, the speed of pressure wave is infinity. If the pressure oscillation propagating at a finite speed could not be removed, then velocity obtained from the momentum equation could not satisfy continuity equation (Patankar 1980; Blosch et al. 1993; Shyy 1994). Therefore, for incompressible flow, mass conservation is satisfied through the driving action of pressure gradients on velocity. For compressible flow the continuity equation is never a function of velocity alone but also includes the density variation. For compressible flow with different Mach number, the role of pressure is also different. In the hypersonic limit, velocity variation is relatively smaller than the velocity itself. Thus, pressure does not influence the mass conservation very much through the velocity variation by the momentum conservation as compared to through the density variation by state equation. Therefore, pressure can be considered to act on density through the state equation, and then the mass conservation is kept through the continuity equation. If Mach number is not as high as that of a hypersonic flow, the pressure has dual roles on mass conservation through the state equation and the momentum equation. For subsonic flow, mass conservation is more readily satisfied by the pressure acting on the velocity through the momentum equation than by the pressure acting on the density through the state equation. However, for supersonic flow, mass conservation is more readily satisfied by the pressure acting on the density through the state equation than the pressure acting on the velocity through the momentum equation (Moukalled et al. 2001). The different roles that pressure plays in mass conservation can be thought of as the velocity scale disparity in the physical processes. Due to the different roles of pressure on the mass conservation for incompressible flow and compressible flow, the pressurebased methods are originally designed for solving the incompressible flow and the densitybased methods are originally designed for solving the compressible flow. Also due to the dual roles of pressure, pressurebased methods can be extended to solve a compressible flow (Shyy 1994) and the concepts of densitybased Navier Stokes solvers were borrowed to solve incompressible flows, such as the artificial compressible method (Chorin 1967) and preconditioning method (Turkel 1999) if the stiffness in the convection coefficient matrixes could be removed or reduced appropriately. If a single method is needed to simulate fluid flow at all speeds, the pressurebased method and preconditioning method are two popular frontiers. The artificial compressible method is seldom used for all speed flow solvers because of the stiff solution matrix which degrades convergence rate (Moukalled et al. 2001). Among the incompressible solvers, there is another popular kinetic method, called the lattice Boltzmann equation (LBE) method. In this method, only one variable, the distribution function, is solved. Distribution function represents the mesoscopic information of the fluid flow. NavierStokes equations can be recovered from lattice Boltzmann equation on the condition of nearly incompressible assumption. Macroscopic variables can be obtained from the moment integration of distribution function (Chen et al. 1998). The governing equation, namely the lattice Boltzmann equation, is linear (the nonlinear feature is hidden in the equilibrium distribution function) and whole computational process in one time step is localized. These give the LBE method an excellent capability to parallelize the computations. Also due to the simple boundary condition schemes that the LBE has, the LBE method is very easy to use for fluid flows with random complex geometries (Chen et al. 1998; Succi 2001). In this dissertation, incompressible viscous flow solvers are used to examine the multi scale computational fluid dynamics problems with interfaces. Compressible NavierStokes solvers could be found in typical computational fluid dynamics textbooks (Hirsch 1990; Fletcher 1991; Shyy 1994; Anderson 1995; Anderson et al. 1997; Chung 2002; Ferziger et al. 2002). 2.2 Incompressible Viscous Flow Solvers 2.2.1 Introduction Incompressible flow is a physical model of general fluid flow on the condition of a very small Mach number. The time scales of the convection velocities are much larger than the time scale of sound speed as the Mach number becomes very small. If Mach number equals to zero, new issues are brought in for incompressible flow computations. As introduced in the previous section, the pressure waves have infinity propagating speed in incompressible flows, which leads more elliptic features to the governing equations, especially for the pressure fields (Kiris et al. 2003). Because of this characteristic of the governing equations, the computational schemes have to be designed to couple the continuity and momentum equations through pressure in order to keep the pressure field from oscillating, by which the conservation of mass can be preserved as the sound speed becomes much higher than the convection velocity components. In dealing with incompressible flows described by the NavierStokes equations, two major approaches are generally used: the primitive variable methods and the vortex methods. The purpose of vortex methods is to remove the difficulty of solving accurately the pressure field in incompressible flows. However, if some complexities, such as chemical reaction, are needed to be added to the fluid flows, the primitive variable methods can easily accommodate these requirements by adding additional equations to the stack (Blosch et al. 1993). The primitive variable approach includes the artificial compressibility method (ACM) (Chorin 1967) and the pressure correction methods (PCM). Pressure correction methods include the marker and cell (MAC) method, the semiimplicit method for pressure linked equations (SIMPLE) (Patankar 1980), and the pressure implicit with splitting of operators (PISO) (Issa 1985). 2.2.2 Artificial Compressible Method (ACM) The incompressible NavierStokes system of equations are written in nondimensionalized form as Continuity Equation: = 0 (2.1) Momentum Equation: au, au, 8p 1 2U, +u +2 (2.2) at 9x x 9x, Re ax2 In the ACM method, an artificial compressibility term is added into continuity equation as following S+ L = 0 (2.3) where p is the artificial density, which is equal to the product of artificial compressibility factor / and pressure, ~P = p8l (2.4) When steady state is reached, the artificial density derivative with respect to time vanishes. How to choose the artificial factor/ is the key for ACM (Kwak et al. 1986). ,/ has to be maintained low enough (close to the convective velocity) to avoid stiffness associated with a large range of eigen value magnitudes. But it has to be kept high enough such that the speed of sound can be as large as possible to achieve the incompressible requirement (Chung 2002). 2.2.3 SemiImplicit Method for PressureLinked Equations (SIMPLE) The SIMPLE method was designed originally to solve the steadystate fluid flows. In this method, the velocity and pressure are solved sequentially in one iteration (which is called the outer iteration, the corresponding inner iteration is the one for solving algebraic equations). How to solve the pressure field is the key in this method. In one outer iteration, if pressure is known (which could be given as a guess for initial condition or obtained from previous iteration), velocity could be solved from the momentum equations. The velocity obtained from the momentum equation might not satisfy continuum equation. This requires correcting the pressure, and further correcting the velocity to make the continuum equation satisfied in this outer iteration. If we denote u v and ) as the velocity and pressure which satisfy the momentum equations, the predictorcorrector procedure with successive pressure correction steps is given as p =P+ (2.5) where p is the actual pressure, and p is the pressure correction. The actual velocity components in twodimensions are u =+u (2.6) v = +v (2.7) Since u, v, and p satisfy momentum equation, u, v and p also should satisfy the momentum equation, we have ae(e +u)= ab(unb u)+b+[(P,+p,)(PE +EA, (2.8) ael = Zanbnb + b +[p E ]Ae (2.9) The subscriptions in these equations are same as those in Partankar's book (Patankar 1980). Subtracting (2.9) from (2.8), we have aeue = anbu +(p PE )A (2.10) From (2.10), it is clear that the velocity corrections on any computational point have two components, one is induced from the difference of pressure correction on the neighboring points and the other comes from the velocity corrections of its neighboring points. The former source of velocity correction, the difference of pressure correction on neighbor points, is the major factor. If the neighbor velocity correction impact is neglected (this is why this method is called semi implicit method), the velocity correction equation (2.10) could be written as aeue = ( p)p eA (2.11) u = A pp p = d( pp p) (2.12) ae Likewise, the velocity correction component in v direction is given as v A (pP , p = d (pP N p1) (2.13) Substituting (2.12) and (2.13) into (2.6) and (2.7), and then substituting (2.6) and (2.7) into the discretized continuum equation, the pressure correction equation could be obtained as app' aPEp +awp' +aNpN + asp' + b (2.14) aE PedAy aW = p dAy aN =pdnAx as = p~dAx ap =aE +aW +aN +aS b (p= pp ) AxAy (pv) A At w e A + Since the source term b represents the mass imbalance based on the velocity in the previous iteration, if its value is small enough, mass conservation could be taken to be satisfied. Therefore, the value of b could be the criterion for stopping iteration. Neglecting the velocity correction from neighboring points in the velocity correction equation does not influence the final steadystate solution. In this sense, the iteration procedure for the pressure can be simplified such that it requires only a few iterations at each time step. SIMPLE method has been modified to several versions. People call them as SIMPLE family methods. These modified methods accelerate solving the pressure correction equation which is elliptic type equation. The detailed descriptions of these methods could be found in many CFD textbooks (Patankar 1980; Fletcher 1991; Shyy 1994; Ferziger et al. 2002). 2.3 Computational Issues on Multiscale Computations 2.3.1 Stiffness Multiscale computations in fluid dynamics and heat transfer means that widely different time or length scales phenomena need to be captured with expected accuracy. As described in Chapter 1, different scales arouse from a wide acting range of different forces in dynamics. These forces can interplay with each other such that the small scale perturbations could amplify and propagate into large scale regions. Even though, for some cases, the small scale impacts can be restricted in small local regions without amplification and propagation, the local phenomena restricted by the large scale dynamic mechanisms are sometimes interested, such as the heat and momentum transfer mechanisms surrounding a heating element of hotwire/film whose length scales are generally much smaller than the measured large scale fluid flows. These kinds of measurement devices are preferred because of their small sizes such that large scale fluid flows are disturbed as little as possible, and the measurement could be more accurate. Numerically, multiscale problems are relative to stiff issue if the physical problems are modeled with sufficient accuracy by a coupled set of ODEs and these ODEs have unique and bounded solutions (Lambert 1980; Lomax et al. 2001; Deuflhard et al. 2002). Generally, the ODEs can be expressed as di= Aii j(t) (2.15) dt where ii is the dependent variable vector, t is the independent variable. If A is independent of t, Eq. (2.15) is linear, otherwise nonlinear. The second term on the right hand side is a forcing function which is determined by the inherent features of dynamic system (Kinsler et al. 2000). The difference between the dynamic scales in physical space is represented by the difference in the magnitude of the eigenvalues in the eigenspace of A. The forcing function can have its own scales. Usually its scale could be adequately resolved by the chosen computational step size. Thus, it is reasonable to assume that the forcing function has no effect on the numerical stability of the homogeneous part in one computational step. However, for nonlinear problem, the matrix A is dependent on t. If the time scale of (t) is needed to resolve, it has a limit on the discretized size of t which has an impact on the eigenvalues of A because the entries of A depend on the discretized size of t. This can be observed when simulating complex flow phenomena such as turbulence, possibly in combination with chemical reactions. In these processes, additional source terms can be added into equations. The nature of these source terms results in an increasing stiffness of these problems, which can reduce the convergence rate (Steelant et al. 1994). The major feature of stiffness is the interaction of computational accuracy and stability. Let us assume the independent variable t in Eq. (2.15) is time. The ODE's could be solved by using timemarching methods. In one time step, the integration with respect to t is an approximation in eigenspace that is different for every eigenvector. Numerically, the eigenvectors corresponding to the small eigenvalue could be well resolved. In contrast, the eigen vectors corresponding to the large eigen values are not (Shyy 1994; Lomax et al. 2001). The approximation to the eigenvectors associated with large eigenvalues requires stable numerical schemes for the global computation procedures and the accuracy represented by these eigen vectors has to be kept from ruining the complete solution. The wide range eigenvalues of multiscale problems bring string requirements on numerical schemes in terms of accuracy and stability. Accuracy and stability are always two contrary requirements for numerical schemes. As a result, the physical nature of multiscale problems and the numerical nature of accuracy and stability make multiscale problems nontrivial. Besides, grid requirement is another factor making multiscale problems harder to resolve. 2.3.2 Grid Requirements for Multiscale Problems If the physical models accurately represent the physical processes, and the numerical schemes can be appropriately selected to capture the physics, the grid arrangement will determine the accuracy of solutions. Multiscale phenomena in fluid flow and heat transfer are generally characterized by a couple of regions in which the spatial gradients of the dependent variables are much higher than in other regions. Consequently, higher spatial grid resolutions are required in those high gradient regions, and a lower grid resolution in others. This grid clustering can strongly affect the eigenvalues of the matrix A since A is calculated directly in terms of the grid sizes. This grid effect on the eigenvalues of the matrix A can be illustrated by the diffusion equation shown below (Lomax et al. 2001). The diffusion equation is au 2u v (2.16) at ax2 Using the threepoint central differencing scheme to represent the second order derivative in space leads to the following ODE diffusion model S B(1, 2,1)i+bc) (2.17) dt AX2 with Dirichlet boundary conditions folded into the (bc) vector, where B(1,2,1)is a banded matrix 2 1 1 2 1 B(1,2,1)= '. (2.18) 1 2 1 1 2 The eigenvalues of the matrix B are 4v mF 1 A = sin m =1,2,...,M (2.19) The stiffness ratio for this diffusion equation is ( 4v "(AX'D 2 Cr 2 4 (2.20) C, M, n I 4v AX2 Ax2 From this ratio, it is clearly seen that the smaller the grid is, the stiffer the numerical equations would be (Atkinson 1989; Lomax et al. 2001). 2.3.3 Methods to Reduce Stiffness In order to achieve the expected accuracy for multiscale problems, the truncation errors should be kept as small as possible. Reducing the grid size is a way to achieve a smaller truncation error if the adopted numerical schemes are consistent and stable. However, multiscale problems often involve large gradients in space or/and in time. For large time derivatives problems, one might suppose that the adaptive time stepsize routines might offer a solution; that is, using small time steps during the rapid transients and large time steps otherwise. But, that is not the case in numerical practices. The reason is that the stability requirement still necessitates very small steps over the entire solution (Chapra et al. 2002). Therefore, Astable methods (generally they are implicit methods) are preferred (Lambert 1980; Shyy 1994). For large variable gradients in space due to multilengthscales, the errors arising from insufficient grid resolution is a major contribution among the many types of errors included in the overall solution (Celik et al. 2004). In such problems, grids should be clustered at locations where large gradients exist so that the truncation errors could be generally kept as small as possible and computational expenses could be reduced (Fletcher 1991; Anderson 1995; Ferziger et al. 2002). Besides clustering grids at the known locations where large gradients exist, adaptive meshing algorithms could be used to solve more complicated problems in which the locations of large gradients might not be known in advance (Bank 1990; Anderson 1995; Park et al. 2004). However, the difficulty with adaptive meshing method lies in that finding the solution of the linear system, as the grid points, that are added for grid refinement are globally combined with the regular mesh (Park et al. 2004). The resulting linear system may also lose the banded structure or the positive definiteness of the uniform grids (Golub et al. 1996). If the physical geometries are very complicated and the length scales have wide distributions in the computational domain, the intensity of clustering grid might be very severe. The aspect ratios of some grid cells might be very large in some subregions. Such computational stencils have poor grid quality since the order of truncation errors for such clustered grid can be lower because some terms in truncation errors could not be cancelled off due to nonuniform grid sizes. As a result, if the variation of aspect ratio of computational stencils could not be kept as small as expected (generally 10% variation is accepted for practical computations suggested by Ferziger et al. (2002)), the truncation error could be larger than acceptable. Thus, for such multi scale problems with large space gradients, further improvement for reducing truncation error is necessary. One remedy is to use multilevel technique (Shyy et al. 1997a; Appukuttan et al. 2003b). This technique includes a global computation made for providing global information in the whole computational domain, which may not resolve the small scale phenomena because of the low resolution in the vicinity of the small scale process regions, and then a reduced domain computation, which includes the small scale process regions, is carried out on a higher resolution gird with the boundary conditions extracted from the global computations. In this method, the length scale of the reduced domain is smaller than the global domain and larger than that of the small scale processes. The global computation can provide accurate boundary conditions for the reduced domain computation if the small scale processes are mainly restricted to the interested regions, and their variations do not affect the flow and heat processes in the regions far away from them. In such cases, the reduced domain boundaries could be selected at such locations where no great impact comes from small scale regions. By this way, the length scale disparity could be reduced. Good grid qualities for both global and reducedregion computations can be acquired. If small scale processes could be amplified and propagated over the global domain or if the physical geometries are very complex, a single global grid of multilevelgrid technique may not be capable of accurately capturing the global phenomena. As a result it can not provide accurate boundary conditions for the reducedregion computations. For such problems Domain Decomposition Methods (DDM) or multiblock techniques can be used as alternative method to reduce the error due to grid resolutions (Shyy et al. 1997b). Unlike in multilevel method that has two separate computing processes, in multiblock methods, the computational domain is partitioned into several blocks according to different physical scales in different blocks and the computational variables are solved simultaneously (without block loop (Thakur et al. 2002)) or iteratively (with block loop (Shyy et al. 1997b)). By reasonable domain partition, multiblock techniques can reduce the topological complexity. Each block can be generated independently. Thus, grids can have expected high resolutions in the required regions. Also, grid lines do not need to be continuous across block interfaces. This makes the local grid refinement and adaptive redistribution easy to accommodate different physical length scales that exist in different regions. The poor quality of single grids on such complex problems is thus improved by block arrangement. Compared to the multigridlayout method, multiblock method can deal with more complicated problems, but, involves more work than multilevel method, in terms of specifying the block interface boundary conditions, data structure, interface region size effect on data exchange and convergence rate (Shyy 1994; Shyy et al. 1997b; Thakur et al. 2002). For very complex physical and geometrical problems, parallel computational techniques using on multiblock grids can improve convergence rate (Golub et al. 1993; Shyy et al. 1997b). Parallel computing techniques depend on both hardware architecture (parallel architecture) and software (parallel algorithms). The above techniques are very useful for multiscale problem computations. Block partition is determined in terms of different scale phenomenon regions, which is generally reflected by local dimensionless parameters in NavierStokes equations, such as Reynolds number. Grid characteristics (Grid stretching and clustering) also yield to different physical scale distributions. In pressurecorrection algorithms, the inner loop/iteration (solving algebraic equations) has a big impact on the outer loop/iteration that couples the unknown variables through conservation equations. The convergence rate of the inner loop can influence the outer loop substantially as the number of grid points increases (Shyy et al. 1997b). It is wellknown that the eigenvalues of the algebraic equations of multiscale problems might have a wide range of values (Lomax et al. 2001). Typical iteration methods, such as pointJacobi iteration on a single grid, do not have fast convergence rate for such equations. While, multigrid methods have the capability to accelerate solving these kinds of equations (Shyy 1994; Shyy et al. 1997b; Tao 2001; Chung 2002). The large range of eigenvalues can be represented by the large range of frequency in terms of grid sizes in Fourier's analysis. High frequency error components damp off faster on fine grids. Multigrid methods resolve different frequency error components on multi level grids so that the low frequency error components on fine grids become high frequency error components on coarse grids and they can be damped off faster on coarse grids. Because of this uniform (in terms of frequency) solution error treatment, multigrid methods are highly efficient to solve physical problems with disparate multiscales (Shyy et al. 1997b). 2.4 A Special Computation Example: ThermoMEMS Computation Results In this section, a multiscale problem with a sharp fluidsolid interface described by the NavierStokes is investigated. This problem is a conjugate heat transfer type resulting from a Microelectromechanical Systems (MEMS)based on thermal shear stress. In this problem, multiscales come from (1) length scale disparity of characteristic length scale of a channel flow and the size of the MEMSbased thermal shear stress sensor, and (2) large conductivity ratio of fluid to solid substrate on which the MEMSbased sensor is deposited. Before addressing the multiscale issues in this simulation, an introduction about thermal anemometry and its issues are presented below. 2.4.1 Introduction to Thermal Anemometry for Fluid Velocity and Skin Friction Measurement Thermal anemometry has been widely used to measure the fluid velocity and the skin friction for decades (Winter 1977; Perry 1982; Bruum 1995). All thermal sensors are temperatureresistive transducers that essentially measure heattransfer rate. They are indirect sensors and thus require an empirical or theoretical correlation, valid for very specific conditions, to relate the measured Joulean heating rate to the flow parameter of interest (i.e., velocity or wall shear stress). All indirect thermal shearstress sensors possess several limitations when used for quantitative measurements. Specifically, they are limited by the difficulty in obtaining a unique calibration or relationship between heat transfer and flow parameters of interest. In addition, the frequencydependent conductive heat transfer into the support substrate (i.e., prongs for a hot wire and supporting substrate for a hot film) reduces the sensitivity and complicates the dynamic response. Thermal sensors are also sensitive to the mean temperature variations which can be difficult to correct. Finally, the thermal sensor is heated well above the ambient temperature during the operation, which can potentially influence the nearwall flow structure and introduce measurement error. A numerical simulation can help interpret the sensor output and correct the calibration formula. For example, Durst et al. (Durst et al. 2002) reported a numerical study of conjugate heat transfer effect on the nearwall hotwire correction. They found that the influence of the local length scale (or Reynolds number), Y+, on the heat transfer process adjacent to hotwire is very important, especially when the thermal conductivity of the wall is very small, such as that of glass or perspex. In the above, Y+ is defined as Y+ = uY (2.21) V where u,= is the friction velocity, and Y is the wiretowall distance. Shi et al. (Shi et al. 2003) conducted twodimensional numerical simulations of forced convection from a micro cylinder in a laminar crossflow. Their numerical results showed that the heat diffusion from the wire is pretty large in the case of small wiretowall distance (Y<3). The reason is the modification of thermal boundary condition at the fluidwall interface caused by property variations between phases. Thermalbased skin friction sensors, which relate the convection from a thin heated film deposited on a substrate to the local wall shear stress, are limited by the heat conduction into the supporting substrate as well as the flow perturbations due to local fluid heating (Naughton et al. 2002; Sheplak et al. 2002; Appukuttan et al. 2003a). Recently, researchers have developed silicon micromachined sensors because of the possibility of improved thermal isolation of the sensing element from the substrate by depositing the sensor on a thin membrane stretched over a vacuum cavity (Ho et al. 1998; Sheplak et al. 2002). While the vacuum cavity structure greatly improved the performance in terms of sensitivity with respect to conventional sensors, MEMS based thermal sensors still suffer from the inherent limitations of localized flow heating and heat conduction into the supporting membrane (Naughton et al. 2002; Sheplak et al. 2002). The inherent small size of MEMSbased thermal sensors invalidates the classical 1/3 power law of the hotfilm theory, i.e., Nu oc rz3 (Lin et al. 2000), which assumes that the thermal boundary layer resides entirely within the linear region of the velocity profile and that the boundary layer approximation holds for the energy equation (Lin et al. 2000). Because of the small size of the MEMSbased sensors, diffusion can substantially affect the fluid flow and heat transfer in the vicinity of the MEMSbased thermal sensors (Ho et al. 1998). Lin (Lin et al. 2000) thus included 2D heat conduction of fluid flow with unidirectional convection, resulting in the following energy equation in the fluid phase: OT (02 T 02 T u = a + (2.22) Ox 2 2y ) In addition, onedimensional heat conduction in a solid membrane was employed to account for the conjugate heat transfer effect. Yoshino et al. (Yoshino et al. 2003) conducted a numerical analysis of frequency response of micro thermal flow sensor. In their simulations, no buoyancy force was taken into account and a linear fluid velocity was given. The heat conduction in the solid substrate was considered in their numerical simulations. They found an optimum diaphragm size to frequency response for their thermal sensor, which is 200300m long. Recently, mixed convection induced by MEMSbased thermal shear stress sensor was studied by Appukuttan et al. (Appukuttan et al. 2003a). They discussed the flow and heat transfer surrounding a thermal shear stress sensor embedded on a wall of a channel. Buoyancy force effects induced by the thermal sensor on the shear stress were examined for different Grashof numbers. As they pointed out, buoyancy force has little impact on the wholedomain flow structure in the channel, while its impact on shear stress measurement can be noticeably observed because the sensor performance depends on the nearwall velocity gradients. In the work of Appukuttan et al's, heat conduction in the solid substrate on which the thermal sensor is deposited was not considered. However, as Naughton and Sheplak stated in their review paper (Naughton et al. 2002), the heat conduction in the substrate can noticeably influence the heat transfer process. In the above numerical studies on MEMSbased thermal sensors (Lin et al. 2000; Yoshino et al. 2003), velocity profiles were given. Only energy equations were solved based on the assumed velocity profiles. Buoyancy force was not accounted for in the numerical simulations. Although onedimensional (Lin et al. 2000) and twodimensional (Yoshino et al. 2003) conjugate heat transfer was considered in the energy equations, the impact of the solid heat transfer on the fluid velocity distribution through the buoyancy forces was absent. In the hotwire simulations of Durst et al (Durst et al. 2002) and Shi et al (Shi et al. 2003), the temperatures of hotwire heating elements were specified and conjugate heat transfer was involved, in order to identify the correction on the velocity profile. However it was found that the buoyancy force has very little effect on the velocity profile (Durst et al. 2002). Hence the buoyancy terms were neglected in the corresponding momentum equations. Sheplak et al (Sheplak et al. 2002) proposed that sensor measurement is quite sensitive to the ambient temperature variation. In this dissertation, we will investigate the conjugate heat transfer around the solid substrate and the surrounding fluid. The effect of the buoyancy force which couples the momentum and energy equations is investigated because it plays an important role in establishing the velocity gradient at the fluidsolid surface. Two kinds of sensor boundaries are used, namely, specified temperature and specified heat flux. In both cases a singlephase (involving the gas phase only), and a conjugate heat transfer with coordinated thermal boundary conditions are considered to highlight the individual heat transfer modes. The effect of heat transfer on the IVIMEMSbased wall shear stress sensor will be addressed. 2.4.2 Governing Equations In the previous work (Appukuttan et al. 2003a), 2D steady state NavierStokes equations with the Boussinesq approximation for buoyancy force were solved. In order to capture the buoyancy effect, velocities are scaled by the buoyancy force term (Shyy 1994). The following dimensionless variables are employed to normalize the governing equations: __ u v T* TTo p u U v T T T p gf.ATH gfATH Tv ... To P x (2.23) P x =_ y* =_ (2.23) p .gfATH H H The nondimensionalized NavierStokes equations subject to the above scaling references are: u momentum u  +v =  + +p* I+2U* (2.24) S* y* Ox* 7 ax*2 y *2 v* momentum S* v v *P 1 2 2 * u* +v* = + + +T (2.25) &* y* y* .7r\Ox 2 y* ) Energy equation ._ _T .* OT 1 2 T* 2 T* + +v _+ (2.26) Ox" y" Pr xGr Ox*2 ay*2) where Pr and r gfA TH3 The last term in Eq. (2.25) represents the dimensionless a V buoyancy force. The current formulation corresponds to a horizontal channel in which the gravity is perpendicular to the streamwise velocity. In the solid domain, heat conduction is the only transport phenomena. Without a heat source in the solid (the sensor is a heat source but treated via the boundary condition due to its very small thickness), the governing equation is, k a2+ = 0 (2.27) where ks is the dimensionless solid thermal conductivity; here, in dimensionless form, it is the solidfluid thermal conductivity ratio. It is noted that the dimensionless fluid viscosity and 1 1 thermal conductivity are and  respectively. For convenience, all asterisks are Gr PrVGr dropped hereafter unless specifically mentioned. In the present study, the Grashof number is defined based on the length (L) of the thermal sensor, even though the governing equations are normalized by the channel height (H). For the MEMSbased sensor cases, it is better to describe the buoyancy effect by the sensorlength based Grashof number Gr, which can be expressed as Gri = Gr j (2.28) In the present study, Prandtl number is fixed as 0.71. 2.4.3 Numerical Schemes The wellestablished pressurebased approach with finite volume formulation is adopted to solve the governing equations. A secondorder upwind scheme for the convective terms and a secondorder central difference scheme for pressure and diffusion terms are used. Detailed information can be found in (Shyy 1994; Shyy et al. 1997b; Thakur et al. 2002). Code validation and numerical accuracy assessment were performed by Appukuttan et al. (Appukuttan et al. 2003a). 2.4.4 Computational Stiffness Computational stiffness arises from the large difference in length and time scales, caused by variations in transport properties and sensortochannel sizes (Chapra et al. 2002). In this study, the dimensionless solid thermal conductivity is 1200, corresponding approximately to the ratio between the MEMSbased thermal sensor materials and air. 2.4.5 Geometry and Grid Layout As shown in Figure 21, the 2D channel is 25cm long and 1cm high. A sensor is located on the bottom wall of the channel, and is 20cm downstream from the channel inlet. The 2D sensor is 200[pm in length, which is much smaller than the characteristic length scale of the channel flow. A solid substrate is below the bottom wall of the channel. The solid substrate and flow channel are symmetric with respect to the channel bottom wall. The grid layout is related to the computational method which is used to overcome the scale disparity. If a single uniform grid is used to capture the entire physical flow and heat transfer process for the present problem, the grid would be very fine in the vicinity of the sensor and coarse away from the sensor. This kind of grid can introduce a noticeable computational error if the grid aspect ratios in some parts of the computational domain are too large (Ferziger et al. 2002). In order to circumvent the length scale disparity, twolevel layout grids similar to those used by Appukuttan et al. (Appukuttan et al. 2003a) are adopted. First, a wholedomain computation (Figure 21) is carried out, which generally does not provide the resolution necessary for the velocity field near the sensor. Based on the wholedomain solution, a reduced domain computation focusing on the sensor region (Figure 22) with a much finer grid is performed. In the reduceddomain computation, the boundary conditions are obtained from the wholedomain computation by using a bilinear interpolation. Unless otherwise stated, the wholedomain grid is meshed nonuniformly with finer grid near the vicinity of the sensor. The fluid region has 342x36 grid points. The solid region has 342x 12 grid points. Such arrangement is necessary because the temperature gradient in the fluid region is larger than that in the solid region due to convective effects. Since the solidfluid thermal conductivity ratio is large, this grid arrangement is also helpful for obtaining the numerical convergence. The relationship of convergence with grid size and thermal conductivity ratio was presented by Shyy and Burke (Shyy et al. 1994). As will be presented later, the grid refinement study supports the adequacy of the present grid system. The reduceddomain computation is employed on a uniform grid consisting of 162x202 points which has much better resolution than the wholedomain grid and is expected to capture the detailed heat and fluid flow structure surrounding the sensor. The issue of grid resolution will be discussed while presenting the results. The boundary conditions for the wholedomain computational domain (Figure 21) are as follows. (a) Boundary conditions for the momentum equations: Inlet: parabolic velocity profile is determined by Reynolds number and Grashof number. (The average inlet velocity is determined by Reynolds number, u = Rev which H is in the dimensional form. If this velocity is normalized by the velocity scale in Eq.(2.23), we have =Re. This velocity value is used to specify the parabolic inlet aveu 7Gr7 velocity.) Channel walls: noslip Outlet: velocity extrapolation (b) Boundary conditions for the energy equation: Inlet of channel: T = (c) Boundaries other than inlet and outlet in Figure 21: 2T S= 0 (2.29) an2 where n is the direction normal to the boundaries. This temperature boundary condition is adopted from the energy equations because the temperature and heat flux on these boundaries are unknown a priori. The boundary conditions of the reduceddomain (Figure 22) are as follows. (a) Inlet: velocity and temperature are interpolated from the solutions of the whole domain computation by bilinear interpolation. (b) Outlet: velocity and temperature are extrapolated from inner grid points. (c) Boundaries of the solid region (except fluidsolid interface): temperature is interpolated from the solution of the wholedomain computation by bilinear interpolation. At the solidfluid interface (except the sensor surface on which the temperature or the heat flux is imposed), the interface conductivity k, (Patankar 1980) is employed. This artificial conductivity arises from the nonhomogeneity of the materials on both sides of the interface. Let us consider two neighboring computational nodes next to the interface, one in the fluid region represented byf and the other in the solid region represented by s. The temperatures on these two points are T7 and T,, respectively. Since no computational node exists on the interface, the heat flux through the interface could be represent as TfT q=k 'f q S (2.30) whereCo y is the normal distance from the node f to the solidfluid interface, and y, is the normal distance from the point s to the solidfluid interface. In Eq. (2.30), k, is the effective, interfacial conductivities. If onedimensional heat conduction resistances between the nodes f and s are considered, the heat flux through the interface can also be expressed as T Ts i y (3 (2.31) kf ky Combination of Eqs. (2.30)(2.3 1) leads ke = y + y (2.32) 3Yf 9y + + kf k The interfacial conductivity ke in Eq. (2.32) is the harmonic mean of the fluid and solid thermal conductivities. Between the sensor and the solid substrate, adiabatic boundary condition is used to account for the vacuum cavity beneath the thermal sensor. The thickness of the vacuum cavity is neglected because it is very small (Sheplak et al. 2002). On the fluid side of the sensor surface, either one of the two thermal boundary conditions is adopted, namely, (a) specified sensor temperature, or (b) specified sensor heat flux. The specified sensor temperature condition is the same as the sensor temperature boundary condition of the singlephase computation previously reported in (Appukuttan et al. 2003a). The specified sensor heat fluxes are calculated from the corresponding single phase cases with the same Grashof number. The values of Grashof numbers are selected based on the operational temperature range of the MEMSbased thermal shear stress sensor (Sheplak et al. 2002) used in the Interdisciplinary Microsystems Group in the University of Florida, which is about 20400C. For air, the Grashof numbers are in the range 0.050.8 based on the sensor length. Therefore, Gi, =0.08 and 0.5 are chosen. Numerically G, =10 and 100 are also used to examine larger Grashof numbers limit, which correspond small kinematic viscosity fluids. For the case with Gi, =0.08, only shear stress distribution and variation are shown in the next section since the flow structure variations in these cases are too small to be observed, while the shear stress variations are noticeable. 2.4.6 Results and Discussion Computations are conducted for a horizontal channel flow in which gravity is transverse to the streamwise direction. The Reynolds number is 500 with the length scale based on the channel height. Before discussing computational results, it is necessary to first discuss the numerical convergence because of the numerical stiffness encountered in the present investigation. The convergence might be problematic if iterative methods are used. To quantify the convergence behavior, residuals of the governing equations are used. Poor convergence is indicated by a lower rate of decrease of residuals. Furthermore, since small values of residuals do not necessarily represent small errors in the solutions for stiff problems, at least a threeorder of magnitude drop in the residuals is needed for judging the convergence. Since the momentum and the energy equations are coupled by the buoyancy force, the stiffness in one equation also affects the accuracy and convergence of other equations. When the thermal conductivity ratio is fixed as a constant, the stiffness arising from the large thermal conductivity ratio is similar for all the cases with different Grashof numbers. This can be seen from the energy equations for the fluid and solid regions. Figure 23 shows the residual reduction histories of u, v, pressure and temperature of the singlephase and conjugate cases using a logarithm scale with specified sensor temperature andGr, =0.5. Both are solved using the algebraic multigrid technique, specifying the same number of inner loop iterations and relaxation factors. The residuals of singlephase and conjugate cases are normalized based on the same normalization values. The four residual history plots show that the singlephase case converges faster than the corresponding conjugate case. These different convergence rates mainly result from the stiffness of thermal conductivity ratio of the conjugate case. To demonstrate the impact of the solution technique, Figure 24 shows u, v, pressure and temperature residual reduction history plots of a conjugate case with Gr = 0.5, using both singlegrid and multigrid solvers. The multigrid solver shows noticeably faster convergence than the single grid solver, indicating that the multigrid solver performs better for solving stiff problems than the singlegrid solver. This phenomenon is also observed in twophase flow problems that generally have multiscale issues and significant property jumps (Francois et al. 2004). In order to assess the resolution requirement for the current problem, grid refinement has been conducted for both the wholedomain computation and the reduceddomain computation. For the wholedomain computation, the fluid region has 514x53 nodes, compared to the baseline grid of 342x36 nodes. The solid region has 514x18 nodes, refined from the baseline grid of 342x 12 nodes. For the reduceddomain computation, the refined grid is still uniform, which has 242x302 nodes in comparison to the baseline grid of 162x202 nodes. The velocity and temperature fields computed on the coarse baseline grid are only marginally different from those on the refined grid. Figure 25 shows that the difference between the shear stresses computed on these coarse and refined grids is small. As expected, the convergence rates of the computation on the refined grid are noticeably slower than those on the coarse grid (see Figure 26), where the same multigrid strategies have been adopted. Hence the baseline (coarse) grid is adopted in the rest of this study. Figure 27 Figure 29 show temperature contours for the singlephase and conjugate cases for different Grashof numbers. In these temperature contour plots of the conjugate cases, similar to singlephase cases, the warmedup regions enhance with Grashof numbers, showing stronger buoyancy effect for higher Grashof numbers. Comparing with single phase temperature contours, near the sensor in the fluid region, temperature distributions of the conjugate cases are very similar to those of the single phase cases. Similarly distorted elliptical contours can be observed for both the singlephase and conjugate heat transfer cases with different sensor boundary conditions. This can be seen from the temperature contour line with the value 0.0625. In the outer fluid region, the singlephase temperature contours have different patterns than the conjugate cases. The temperature contours of single phase cases originate from the sensor and expand asymptotically to the main stream, while the temperature contours of the conjugate cases in that region no longer connect with the sensor since the heat conduction effect of the solid substrate raises the upstream and downstream fluid temperature, and thus, the effective heating length is enlarged due to this substrate heating. At the solidfluid interface, temperature contours have abrupt kinks, reflecting the effect of the large solidfluid thermal conductivity ratio. In the leadingedge regions of the conjugate cases, the temperature contours turn in different directions at the solidfluid interface. Near the sensor, temperature contours from the solid region turn left into the fluid region. Beyond a small upstream distance from the sensor, the temperature contours turn right from the solid into the fluid region. In the vicinity of the sensor at the leading edge, heat goes upstream in the fluid by diffusion near the sensor which heats up the solid substrate. Since the heat diffusion effect in the fluid region at the leadingedge is restricted very close to the sensor by the fluid flow, the temperature gradient of the fluid in that region is larger than that of the solid. Therefore, beyond a small upstream distance from the sensor, the temperature of the solid is higher than that of the fluid at the solidfluid interface. This leads to the fact that the solid substrate heats up the fluid due to the heat conduction in the solid region and the large solidfluid thermal conductivity ratio. Then the heat transferred from the solid to the fluid is carried downstream by fluid convection. At the trailing edge, fluid temperature is always higher than that of the solid since the heat diffusion in the fluid enforces heat convection. As a result, the temperature contours at the interface downstream of the sensor have similar shapes. In the solid region, heat flux is determined by solid thermal conductivity, solidfluid interface thermal conductivity and convection effect of fluid at the interface. The interface heat resistance results from both conduction and convection at the interface. The interface resistance of conduction can be indicated by the interface conductivity (Eq.(2.32)). Because the solidfluid thermal conductivity ratio is fixed, the interface thermal conductivity can be rewritten as 3)Af + g2 kf k,= k s + k (2.33) Therefore, the dimensionless interface conductivity varies with dimensionless fluid 1 conductivity which is proportional to In other words, the resistance from conduction increases with increasing Grashof number. The interface heat resistance from the convection can also be indicated by Grashof number since the Reynolds number is fixed as 500 for the present study. The resistance from convection decreases with increasing Grashof number because larger Grashof numbers denote larger buoyancy forces which enhance the convection effect. Thus, the two parts of the interface heat resistances have different tendencies of variation with increasing Grashof number. Dimensionless solid thermal conductivity also varies with dimensionless fluid conductivity in a similar way owing to fixed solidfluid thermal conductivity ratio. Hence, convective heat transfer enhancement in the fluid region reduces the conduction effect at the solidfluid interface, and vice versa. Due to the interface resistances, the major part of the heat generated by the thermal sensor is carried away by convection. Therefore, the wall has a function to isolate fluid and solid, which allows relatively small part of heat to transport into the solid substrate. This is the reason why, near the sensor, the temperature contours of conjugate cases do not vary so much from the singlephase cases. Comparing the conjugate cases with specified sensor temperature, temperature contours of G, = 0.5 case are very similar to those of G, = 10. For the conjugate cases with specified sensor heat flux, at the same locations, temperatures of the case with Gi = 10 are lower than those ofGr = 0.5. While for the conjugate cases withGr = 100, temperatures in the solid are higher than the other cases. This might result from the different variation trends of conduction effect and convection effects at the interface with regard to Grashof numbers. This opposite effect of Grashof number on the heat conduction and convection at the fluid solid interface is shown more clearly in Figure 210. For the cases of the specified sensor temperature, Figure 210 shows the temperature distributions on a crosssection, which is perpendicular to the sensor surface and originates from the middle point of the sensor to the top boundary of the reduceddomain. Due to the opposite effect of Grashof number on the heat transfer on the fluidsolid interface, temperature distributions of the cases with G, = 0.5 and Gi = 10 are very similar. For the case with Gr = 0.08, the temperature of the fluid far away from the sensor is the lowest of these four cases, which mainly results from the lower buoyancy force effect on heat convection in the fluid region. For the case withGr = 100, temperature is higher in the fluid region far away from the sensor and smaller in the region closer to the sensor. It can be explained that, for the case withGr, = 100, the convection effect is larger in the far field region and the diffusion effect is smaller in the region near the solid surface. Figure 211 presents the uvelocity profiles of conjugate cases with different sensor boundary conditions. Due to the wall isolation effect, the uvelocity profiles of conjugate cases do not have noticeable differences from those of the singlephase results (Appukuttan et al. 2003a). This is also because of the isolation of the wall described previously. However, the wall effect on the shear stress can be clearly observed in Figure 212 and Figure 213. In Figure 213 shear stress variations are calculated based on the analytic shear stress of a developed channel flow without buoyancy force acting on the fluid. The analytical shear stress du U isp A l. If it is normalized by p, and U, is the average velocity of the channel flow, the dy H analytic shear stress becomes d(u* /IU ) = dy =6 (2.34) Comparing singlephase shear stresses (Appukuttan et al. 2003a), for the small Grashof number cases, the shear stresses of conjugate cases with the specified sensor heat flux have smaller variation than those of the specified sensor temperature cases; for the larger Grashof numbers, their shear stress variations are higher than those of specified temperature cases. Therefore, different energy boundary conditions for the sensor lead to different shear stress distributions and variations for conjugate cases. This difference should be caused by the buoyancy force. Different sensor temperature conditions cause different temperature distributions surrounding the sensor, which generate different buoyancy forces. Therefore, if a thermal sensor works in the constant temperature mode, its shear stress behavior is different from those of thermal sensors working in constant current mode or constant voltage mode. 2.4.7 Summary and Conclusion: The effects of conjugate heat transfer on shear stress variations of a MEMSbased thermal shear stress sensor are investigated. In the very close vicinity of the sensor, temperature contours of conjugate cases with different sensor temperature conditions are similar to those of the single phase cases, while the temperature contours in the region away from the sensor are different from the singlephase cases. For conjugate cases, the effective heating length is enlarged due to the substrate heating, which was emphasized in (Sheplak et al. 2002). Even though the velocity distributions of both singlephase and conjugate cases are very similar, shear stress distributions for the conjugate cases have observable deviations from singlephase cases. Therefore, though buoyancy force does not change velocity profiles significantly, it noticeably influences velocity gradients near the thermal sensor, which introduces errors in the shear stress measurements. Inlet 1 Fluid Outlet 1 Substrate / Focused region boundaries 20 Figure 21. Geometry of channel flow with solid substrate 0.32 inlet inle( 13 2 0.2 Fluid outlet / 0.2 Substrate sensOr  0.1 Figure 22. Reduceddomain geometry. 1, 2 and 3 represent the sections on which velocity profiles are plotted in Figure 211  singlephase  conjugate  singlephase  conjugate U ' 2 2 I " 3   s e N 33 3 N 3 4 4  5 5 . 0 1000 2000 3000 0 1000 2000 3000 a iteration number iteration number b 0 1 1 s1 inl1 epsinglephase singlephase  conjugate 2 conjugate 2  2 3 7 7 5 7 N 6 8 7 S 9 8 10 0 1000 2000 3000 0 500 1000 1500 2000 2500 3000 C iteration number iteration number d Figure 23. Relative reductions of residuals of singlephase and conjugate cases with specified sensor temperature, GCr = 0.5. (a) u residual, (b) v residual, (c) pressure residual and (d) temperature residual 0 0  smgle grid solver  multi grid solver 25 0 1'II 201 3000 4 5000 6000 00 200 00 4000 00 100 0 1 smgle grid solver smgle grid solver 2 2  32 3  3 4  5 5 3 4  4 \ 5 0 1000 2000 3000 4000 5000 6000 0 1000 2000 3000 4000 5000 6000 a iteration number iteration number b 0 0 1 solved bysingle grid solver solvers,1 Gr = 0.5. (a) u residual, (b)d solver (c) pressure residual angrid solved) temperature multi grid solver 2 3 3 ~4 1 5\ \ ~6\ 7 \8 0 1000 2000 3000 4000 5000 6000 0 1000 2000 3000 4000 5000 6000 C iteration number iteration number d Figure 24. Relative reductions of residuals of conjugate cases with specified sensor temperature solved by singlegrid and multigrid solvers, Gr, = 0.5. (a) u residual, (b) v residual, (c) pressure residual and (d) temperature residual. 6.31  coarse grid a 6.25  fine grid 6.2 6.15 6.1 6.05 6 5.95 5.9 0.1 0.2 0.3 x Figure 25. Shear stress comparisons for the coarse grid and refine grid computations (Gr, =0.5 and Re= 500)  coarse grid fine grid  coarse gnd fine grid 3 3 N 4  1 5 \ I "4 \ 6   \ 2 5 \ o N 9 N 77 0 1000 2000 3000 4000 0 1000 2000 3000 4000 iteration number iteration number The wholedomain grid 0 1 coarse grid 1  coarse grid 2  fine grid  fine grid 2 3 0 \ 3 \ 4  4 5 o/ = 5 6  S6 S1 7 9  0 2000 4000 6000 0 1000 2000 3000 4000 5000 6000 7000 iteration number iteration number The reduceddomain grid Figure 26. Relative reductions of Pressure and temperature residuals on the refined grids ( Gr = 0.5 ) 025 02  0 15 01 9 / 01 00177 0 0 5 0 0 62 5 I Z " S01 02 03 0 1 02 03 b 0 01 02 03 0 375 a 0 b 1 c Figure 27. Temperature contours with Gr, = 0.5 (a) Singlephase cases with specified sensor temperature; (b) Conjugate cases with specified sensor temperature; (c) Conjugate cases with specified sensor heat flux. 01 1 127 / 0 [5 00625 S125 015 015 0012 0 0 020 01 02 03 02 1 02 03 a 0 b 0 l 02 03 Figure 28. Temperature contours withGr5 =10 (a) Singlephase cases with specified sensor temperature; (b) Conjugate cases with specified sensor temperature; (c) Conjugate cases with specified sensor heat flux. '5 3E )5 00625 oiS^,'26 02 a 00625 S 0125 000205 00202 00196 01 02 03 x C Figure 29. Temperature contours withG, = 100 (a) Singlephase cases with specified sensor temperature; (b) Conjugate cases with specified sensor temperature; (c) Conjugate cases with specified sensor heat flux. S0 1 Gr s=0 08 0.15 Gr s=0 5 Gr s=10  Gr s=100 a 0.1 0.05 10 10 10 T Figure 210. Temperature distribution on a crosssection which originates from the middle point of the sensor to the top boundary of the reduceddomain (Specified sensor temperature cases) 02 02 0175 0 S Un1 0175 S1: 015 Secton3 015 on 0125 0 125 A 01 01 0 075 0 075 005 005 0025 0025 025 05 075 025 05 075 U4(GH L)/Re U (Grf/L)/R0 (a) Specified sensor temperature (b) Specified sensor heat flux Gr, =0.5 02 02 ] 0175 s 0175 S 015 015 0125 0125 A 01 01 0075 0075 * 005 005 0025 0 025 025 75 025 05 075 UL (G^H/I/)R (a) Specified sensor temperature (b) Specified sensor heat flux Gr, = 10 02 02  0175 0 175  015 015 0 125 0 125 * A 01 A 01 0 075 0075  005 005 * 0025 0025  0 025 075 0 025 075 (a) Specified sensor temperature (b) Specified sensor heat flux Gr, = 100 Figure 211. u velocity profiles of conjugate cases on the crosssections just before (section 1), after (section 2) and in the region of the sensor (section 3) at Re=500 631 6 2 61 59 0 05 0 1 <I n n n; Figure 212. Wall shear stress distribution: (a) Specified sensor temperature (b) Specified sensor heat flux , I II I I I I I 0 02 04 Grashof Number (Gr ) 0 01 02 03 04 05 06 Grashof Number (Gro) Figure 213. Shear stress variation: (a) Specified sensor temperature (b) Specified sensor heat flux ' ' 01 015 02 025 40 CHAPTER 3 ERROR ASSESSMENT OF THE LATTICE BOLTZMANN METHOD FOR VARIABLE VISCOSITY FLOWS As an alternative CFD tool, lattice Boltzmann Equation (LBE) method has been developing for about two decades. From scale classification, the LBE is a mesoscale method for fluid dynamics computations. NavierStokes equations for compressible flows can be recovered from the LBE with the incompressible constraint. In other words, the LBE describes fluid flows with very low Mach number (Chen and Doolen, 1998; Succi, 2001; Yu et al. 2003a). Although there are some difficulties for the LBE method in compressible flows and thermofluid flows, it has been successfully used for isothermal complex fluid flows, such as interfacial dynamics (He et al. 1999), turbulent flows (Yu et al. 2005) and porous media fluid flows (Chen et al. 1991). Since large velocity gradients often accompany substantial viscosity variations due to either shearthinning effects or turbulence models, errors associated with the variable viscosity model need to be addressed before exploring the accuracy for these problems. With variable viscosity such as for turbulent flows treated with the eddy viscosity models (Filippova and Hanel, 1998), or for fluids with flowdependent properties, additional truncation error appears in the macroscopic equations derived from the LBE (Hou et al. 1995). With the help of the truncation error, computational accuracy for those flows can be further understood. The truncation error behavior of the LBE with constant viscosity has been studied by Holdych et al. (Holdych et al. 2004). However, the truncation error for variable viscosity problems has not been studied systematically in the literature. Like other computational models, the LBE method needs to address the wall boundary conditions (Mei, et al., 1999). In particular, the noslip wall boundary condition in the LBE method is based on the bounceback concept (Chen and Doolen, 1998; Succi, 2001). Much work has been done to extend this simple scheme to arbitrary curvature boundaries with second order accuracy (Mei et al., 1999; Ladd, 1994; Yu et al., 2003; Bouzidi et al., 2001). These methods can handle geometrical complexities easily; however, they cannot exactly recover the noslip boundary condition (He et al., 1997) at the mesoscopic level (Mei, et al., 1999). This inconsistency gives rise to the boundary condition error. Although the noslip boundary condition of the LBE has been derived from hydrodynamic conditions on walls by Noble et al. (1995) and Inamuro et al. (1995), these exact noslip boundary schemes are incapable of handling geometrically complex boundaries because they require the computational nodes to coincide with the physical boundaries. The original bounceback scheme and other improved treatments have been successful in constant viscosity laminar flows with expected accuracy (Chen and Doolen, 1998; Yu et al., 2003; Mei et al., 1999; Mei et al., 2000). However, for variable viscosity problems, the boundary condition error for noslip walls may act on the truncation error. Thus, besides the truncation error behavior itself, whether the secondorder accuracy of the bounceback scheme can be maintained in the presence of the variable viscosity is another open question. In this Chapter, these issues are investigated via a fullydeveloped laminar channel flow with a specified variable viscosity. With the help of the finite difference analysis, the truncation error of the LBE with variable viscosity is investigated. Two different specified viscosity distributions, which lead to different boundary layer characteristics, are employed to examine the errors associated with the variable viscosity, in the LBE model. The exact solution of this channel flow exists and can be used for error analysis. In the following, the LBE method and the boundary conditions will first be reviewed. We then present the variable viscosity laminar channel flow equations along with the exact solution. Based upon the exact solution and the finite difference analysis, the error behaviors due to the variable viscosity and the boundary condition schemes will be assessed. To separate the error associated with the variable viscosity from that with the boundary treatment, both Noble et al.'s scheme (Noble et al. 1995), which is exact but restricted to straight boundaries, and bounce backonlink (BBL) scheme, which is not exact but can handle irregular geometries, will be employed. 3.1 Introduction to Lattice Boltzmann Method and Wall Boundary Condition 3.1.1 LBE BGK Method The Lattice Boltzmann equation with the singerrelaxationtime (SRT) BhatnagarGross Krook (BGK) model can be written as (He and Luo, 1997a) fa (x, +e.9t, t +9t) f (x,, t) = x (x,,t) f eq (x,,t)] (2.35) where f, denotes f(x, e, t), which is the distribution function in the direction of the ath discrete velocity ea, f(eq) is the corresponding equilibrium distribution function in the discrete velocity space, r =A/ t, Ais the relaxation time, and x) represents computational nodes in physical space. The most popular lattice model for simulating twodimensional flows is the ninevelocity square lattice model, which is often referred to as the 2D 9velocity (D2Q9) model (Qian et al., 1992). In this lattice model, the discrete velocities (ea) are as follows: eo = 0, e =c(cos((a1);r/4),sin((a1);r/4)) for a=1,3,5,7 (2.36) e = ,c (cos((a1);r/4),sin((a1);r/4)) for a=2,4,6,8 where c = 8x /St. The corresponding equilibrium distribution functions are defined as 3f(eq) W 3 u+9 U)2 3 U f[Ieq) 2 er 2+ (e 2C.21 u I a= 0,1,...,8 (2.37) c/ 2c 2c where w, is the weight coefficient given by 1 W1 = W3 = 5 = W7 = W2 = 4 9 1 W6 8 =W 36 The density and momentum fluxes can be obtained from the moments of the distribution function as 8 a=0 8 a=0 8O f(eq) a=0 f (eq) a=0 (2.39) The kinematic viscosity associated with the D2Q9 lattice model can be expressed as V= )T c23t I 2 S (2.40) where c is the speed of sound, which is equal to c/l3 for D2Q9 lattice model. The corresponding equation of state is p = pc2. For a hexagonal lattice model (Noble et al. 1995), the discretized velocities are e = 0, e= c(cos((a1)2;r/6),sin((a1)2;r/6)) for a= 1,2,...,6 (2.41) The corresponding equilibrium distribution functions are f (eq) do _ (u.u) c (2.42) f (eq) =p d pD +pD (D +2) 2 pD S + 2(e . u) + (e u) (u . u) for a= 1,2,...,6 b c2b 2c4b 2c2b Where for the 2D hexagonal lattice model, the dimension rank is D = 2, the number of lattice direction is b = 6, the average rest particle density is do = p/2. The kinematic viscosity for 2D hexagonal lattice model is =2r 1 3x2 8 8t (2.43) 4 w = 9 (2.38) The LBE is usually solved in two steps: Collision step: (xt+3t)= f(x t) [f(x,t) f (xt)] (2.44) Streaming step: f. (x + e. t, + 8t)= (xx,, + 9t) (2.45) where f, represents the postcollision state. In order to get the NavierStokes equations, the ChapmanEnskog expansion can be used with time and space being rescaled as t = t, t2 = '2t, X, = EX, a a 2 a a a (2.46) t 8t8 8t2 at at, at2 a 1 &: and the distribution function is expanded as f = fl0 f(1) 2 (2) O( 3). (2.47) From Eq. (2.46) and (2.47), it is clear that the ChampanEnskog expansion is essentially a multi scale expansion (Frisch et al. 1987). In the incompressible flow limit, that is <<1, the Navier C Stokes equations can be recovered from the LBE (Eq. (2.35)) with the leading truncation error of O (3x)2 (u)3 I(Hou et al. 1995): =0, (2.48) chu chu 1 p a +u 1 = +V2u (2.49) at 8x, p 8x 3.1.2 Boundary Conditions In contrast to NavierStokes equations, the lattice Boltzmann equation has no physical wall boundary conditions for distribution function f,. The distribution functions are unknown and have to be constructed from inside information. Symmetric and periodic boundary conditions can be constructed without any ambiguities. Outlet boundary conditions are generally extrapolated from internal nodes (Yu et al. 2003a). Inlet boundary conditions can be extended directly from the boundary treatments for solid wall. Since a flow with a large Reynolds number usually has boundary layer phenomenon near walls, which is a typical multiscale process (Batchelor 1967) and sensitive to the wall boundary condition, the noslip wall boundary condition is thus the most important one among other boundary conditions in fluid dynamics. Historically, not only the models of the LBE evolved directly from its predecessor, the latticegas automata (LGA) (Frisch et al. 1986), but also the wall boundary conditions, which is namely the bounceback scheme (Lavallee et al. 1991; Chen et al. 1998; Succi 2001). This boundary condition method has a very attractive feature: it is very easy to implement computationally, especially for complex geometries. In comparison with the NavierStokes solvers, this amazing feature offers the LBE method a distinguished capability to handle flows with geometry complex (Chen et al. 1998; Mei et al. 1999; Mei et al. 2000). For the noslip wall boundary condition, if computational boundary nodes are located on walls, the LBE method with the bounceback scheme has firstorder accuracy (Ziegler 1993). A slip wall velocity exits (He et al. 1997; He et al. 1997c) and increases with relaxation time (Noble et al. 1995). This bounceback scheme with the computational nodes on walls can be referred as the bouncebackonnode (BBN) scheme. While, if the boundary is shifted into the fluid by one half mesh unit, i.e. placing the wall between the computational nodes, secondorder accuracy can be achieved by the bounceback scheme (Ziegler 1993; Ladd 1994a; Ladd 1994b; Yu et al. 2003a). This bounceback scheme with wall positions locating between computational nodes can be referred as the bouncebackonlink (BBL) scheme (Mei et al. 1999). In the BBL scheme, after a particle distribution function f streams from a fluid node at xf to a boundary node at xb along the direction of e, the particle distribution function f, reflects back to the node xf along the direction ofe, (= e), f = fa, (2.50) which is illustrated in Figure 31. The BBL scheme has secondorder accuracy (Ziegler 1993; Ladd 1994a; Ladd 1994b), even through the slip wall velocity can not be exactly removed. The BBL scheme can be improved to arbitrary curved wall boundaries with secondorder accuracy, such as Filippova and Hanel's scheme (Filippova et al. 1998; Fillippova et al. 1998). Mei, Luo and Shyy's scheme (designed as MLS's scheme (Lai et al. 2001)) improved the numerical stability of Filippova and Hanel's scheme (Mei et al. 1999). Bouzidi et al. (Bouzidi et al. 2001) presented a simpler bounceback boundary condition for a wall located at an arbitrary position and no requirements for constructing fictitious fluid points inside walls (Bouzidi et al. 2001; Yu et al. 2003a). Yu et al (Yu et al. 2003b) presented a unified boundary condition based on Bouzidi et al. and Mei et al.'s works. Although the exact noslip velocity boundary condition could not be achieved by using the BBL and the improved BBLs, they are very attractive for flows with complex geometries (He et al. 1997c). Besides the bounceback schemes, there are some other schemes for noslip wall boundary condition, such as Noble et al.'s hydrodynamic boundary condition scheme for hexagonal lattice (Noble et al. 1995), and Inamuro's scheme for square lattice (Inamuro et al. 1995). Noble et. al's scheme was derived for the hexagonal lattice from the macroscopic mass and momentum at a noslip boundary without any assumptions. This scheme is illustrated in Figure 32. On the boundary nodes, the following expressions of the unknown distribution functions can be obtained by solving density and mass flux from the boundary nodes: f2 + A3 = P(fO + fl + f4 + f + f6 ) f2 f3 = 2pu (2f, 2f4 f5+ f6) (2.51) f2 + f3 =(2/ )pv+(f5 +f6) Through Eq.(2.51), the unknown distribution functions on boundary nodes can be obtained from the distribution functions on the nodes of fluid side. This scheme can provide the exact noslip wall boundary condition, which can be exhibited through simulations of a fullydeveloped channel flow with constant viscosity. The absolute L2 norm error of this simulation defined as H [ULB (y) u.,t (y)] dy\ E2 = [IILBM H f (2.52) is shown in Figure 33, which has the order of the roundoff error. This demonstrates the exact noslip boundary condition can be attained via Noble et al.'s scheme. Inamuro's scheme was developed for the 2D square lattice model (D2Q9). A lattice cell is shown in Figure 34. The definitions of the density and velocity on the noslip wall are same as those on the internal nodes: 8 8 8 8 p = i/f = e.e pi = 4cf = fq a=0 a=0 a=0 a=0 In Figure 34, the boundary node is on the wall. Thus, the distribution functions f f, f, fs, f6, f8 are known; while f2, f, f4 are unknown and they need to be constructed from the known distribution functions. Because the slip wall velocity exists, Inamuro et. al. assume the unknown distribution functions f2, f3, f4 have similar forms as other known distribution functions except that there is a counter velocity which is added to the slip wall velocity to force the wall velocity as zero. The unknown distribution /2, f, and /4 are expressed as: f2 = P 1+3(UW+u +vW)+ (u,+u +v) _ [(u +u) +v2, (2.53) f3= p{1+3vw + 3[(u+ +)2v+ (2.54) f4 P 1+3 uwu, +vw )+ UW u +vW u +L +V)2 (2.55) where uw, vw are the slip wall velocity, u is the counter velocity, and p is the density change due to the counter velocity. The velocity on the wall in the horizontal directionuw and vertical direction v, are calculated as u =fe2 (A 4 +f8 6 +f 5f) =(f2 4 +A) (256) Pwa Pw Pw vw= Z f/e,= (=f2+f4+f3 6 8 7)= (f2+f4+f+ B) (2.57) pw a Pw Pw where A = f +f / and B =1ff f1 In the above equations, the unknowns are u', p and pw. With some algebraic manipulations, the unknowns can be obtained u = [6 u 3uwvw (2.58) 1+ 3v, p p'6 p v B 1+3v, +3v2 v1 [( (2.59) S vW (2.60) It can be seen that u', p and p, are obtained from the known distribution functions and wall velocity. Substitutingu', p and p back into(2.53), (2.54)and(2.55), the unknown distribution functions on the wall can be obtained. Inamuro et. al.'s scheme could enforce the noslip wall boundary condition. However, it could not remove all errors from the boundaries. This can be observed from simulations of a fullydeveloped channel flow with a constant viscosity. The absolute L2norm error is shown in Figure 35. It is clearly shown that the absoluteL2 norm error is larger than the roundoff error (O(1015)). The error in Inamuro et. al.'s scheme might come from the counter velocity assumption which is the only hypothesis in this scheme. Whether the unknown distribution functions really have the form (Eq. (2.53), (2.54) and (2.55)) or not is unknown. Through the comparison between Noble et al.'s scheme and Inamuro et al.'s scheme, it is clear that Noble et. al.'s scheme can exactly realize noslip boundary condition, while Inamuro et. al.'s scheme does not. Noble et al.'s scheme is thus adopted in this chapter. 3.2 Error Assessment of the LBE Method due to Variable Viscosity via a Fully Developed Channel Flow 3.2.1 Fullydeveloped Laminar Channel Flow with Variable Viscosity In this section a 2D fullydeveloped steady laminar channel flow with variable viscosity is adopted to assess the error behavior of the LBE with variable viscosity. The governing equation for this flow in a channel of height H is + I d o 0 and u = 0, du =0 (2.61) p dx dyv dy) o dy yH2 The viscosity vo,,, is modeled as Total V0 VSt, S2 ( )2 0 < = < 0.5, (2.62) vo 02(1+13/3) [1+(1/70)031/3' H where vo is a constant and the symmetry gives vo,, for 0.5 parameters a and 0 control the profile shapes for viscosity and velocity. Employing the boundary conditions and integrating Eq. (2.61) twice yield u H2dp j 0.5)( 4 ,J) for 0<4<0.5 (2.63) pvdxo + +,5 Let a = /02 and r73 +ar/2 + = (r7+a)(1r2 +br +c). It is easily seen that a3 aa2 =6, b =aa and c =, where the value of a can be obtained numerically in terms of a and a Integrating Eq. a (2.63) results in the exact solution of this laminar channel flow, that is, (q2 q)Dq+Aln(ql/a+1)+B [In(q2 +bq+c)ln(c)] H2 dp 2 2 .6 uM(*) = (2.64) a[(1/2+a)a2 +3] B a(1+2aa) where D= a, A= a[( ,)a2 +] B = 3(+2a) 8 and Q= (Ab+Ba). a2 + 2c a(a2 + 2c) In order to explore the truncation error due to variable viscosity under the condition of large velocity gradient over a short distance, the parameters 3 and 0 in Eq. (2.62) need to be chosen carefully so that boundarylayer like velocity profile can be obtained in Eq. (2.64). The first term in Eq. (2.64) in the curly bracket corresponds to the parabolic velocity profile with the constant viscosity v0. Other terms are due to variable part of the viscosity vt For small values of a and 0, the fourth term [in(r72 +b7+c)ln(c)] does not change dramatically across the whole channel. The second term Dr is a negative linear term and the third term Aln(Q7/a+1) is positive across the channel, the sum of these two terms varies slowly in the near wall region. The last term in the near wall region varies rapidly because, for small 8 and 0, under further assumption of 0 >> 52 (2.65) it is asymptotically equal to I Bl b2 [tan b _tan1 b 1l0 tan ' (2.66) I 2 2 W {2 0 Thus, with small values of 8 and 0, the tan (...) term in Eq. (2.64) results in the boundary layer type of behavior. For comparison purpose, two sets of viscosity distributions, which satisfy Eq. (2.65) for small 3 and small 0, are used in this study: Casel: (o,3,0)= (0.004167,0.0005,0.0102) CaseHl: (vo, ,0) =(0.008333,0.002,0.0289) The corresponding viscosity profiles are shown in Figure 36. The resulting velocity distributions corresponding to these two viscosity distributions are shown in Figure 37. The boundary layer effect from the tan' (...) terms in Eq. (2.64) or Eq. (2.66) for Case I is shown in the inset of Figure 37. Eq. (2.66) also gives a guideline for estimating the grid resolution required for resolving the boundary layers. Since tan1 (1)/tan 1(oc) = 1/2, it is seen that over a distance of 7 = 0, the velocity reaches 50% of the maximum given by Eq. (2.66). Thus, for dimensionless grid size h=1/H that is close to or larger than 0, numerical solutions will not have sufficient resolution for this thin layer. 3.2.2. The Lattice Boltzmann Equation Treatment In LBE method, the viscosity is associated with the relaxation time r by Eq. (2.40). The variable viscosity in LBE can still be realized via a spatially varying relaxation time, and the total viscosity at a given location can be expressed as 2r(x, y)1 total (x, y) = (xy) for square grid 6 (2.67) total (x, y) = for hexagonal grid 8 Thus, the relaxation time can be represented by the local fluid viscosity as r (x, 6v (x, y)l for square grid 2 (2.68) (x, y) = 8vtrtl (x, y) +1 for hexagonal grid 2 The computational setup for BBL scheme is shown in Figure 38. The first computational node is half a lattice away from the channel wall. For Noble et al.'s scheme, the setup is the same except that the grid is hexagonal and the computational boundary nodes are on the channel walls. For comparison purpose, the computations on hexagonal grid and square grid should have the same grid resolution across the channel height. For the same grid resolution the channel height of the square grid setup is 2//3 times that of the hexagonal grid setup. In order to make absolute comparison among the velocity profiles, the pressure gradients are accordingly adjusted in each computation so that H2 dp remains the same. Periodic boundary condition is used at the left and dx right side boundaries. The constant pressure gradient in Eq. (2.61) is treated as a body force, and is added to the distribution functions after the collision step. The error analysis for the computations with both boundary conditions is carried out by examining the difference between the exact solution and the computational results at each y location. 3.2.3 Assessment Since Case I exhibits thinner boundary layers and sharper nearwall velocity gradients than Case II, we first discuss Case I. Figure 39 compares three velocity profiles for H=50 (h=1/H=0.02): a) exact solution; b) LBE solution using the standard BBL boundary condition and square lattice formulation; c) LBE solution using Nobel et al's exact boundary condition and hexagonal lattice formulation. For the parameters considered, H=50 (h=0.02 >0 = 0.0102) does not resolve the thin boundary layer near the wall, as clearly shown in the inset of Figure 39. It is noted that due to the exactness of the boundary condition used, the hexagonal formulation gives only a slight overshoot in the velocity one full lattice away from the wall (at h/O 2 ); however, it results in an overshoot for the rest of the lattices in the channel. The LBE solution using square lattice formulation suffers from both the inaccuracy of the bounceback boundary condition and the insufficient near wall resolution and gives lower velocity throughout the entire channel. Both LBE velocity profiles have errors of similar magnitudes. This implies that when the velocity in the boundary layer is not sufficiently resolved, the exact boundary condition with the hexagonal lattice formulation would not offer any advantage compared to the approximate BBL condition using the square lattice formulation. For h ranging from 0.004 to 0.1 (H: 10 250 in lattice unit) the square lattice formulation with the BBL boundary scheme shows different over/under predictions for the velocity profiles from the hexagonal formulation comparing with Noble et al.'s scheme. While the results are not shown here for brevity, it suffices to note that the magnitudes of the overpredictions in the hexagonal lattice case are comparable to those of the underpredictions in the square lattice case. These velocity profiles indicate that the errors of the LBE solution using both boundary condition schemes are comparable. Thus, the relative L2  norm error of the LBE solution, defined as { [I0 ,LM (y)u_,, (y)]2 dY E U 1y 2 (2.69) ru2o (Y)dy _0 is examined for both boundary condition schemes. The LBE computations with these two sets of viscosity distributions are carried out by using both Noble et al's scheme and the BBL scheme for grid size h ranging from 0.004 to 0.1. Their relative L2 norm errors in the velocity profiles with respect to grid resolution are shown in Figure 310a. As expected, the relative L2 norm error curve of Case I shifts up with respect to that of Case II because Case II has a thicker boundary layer which implies better computational resolution than in Case I for the same h (= 1/H). For sufficiently high grid resolution (h<0.01, points AD shown in Figure 310a), both Noble et al's scheme and the BBL scheme yield the asymptotic second order accuracy, which is consistent with the truncation error analysis for both square and hexagonal lattice schemes. For the present fully developed 2D channel flow with low speed, the velocity field satisfies Vu = 0 and O(u3), resulting in very small modeling error. Thus the errors in LBE computations are mainly from the truncation error (due to variable viscosity) and the boundary condition treatment. With sufficient grid resolution (points AD in Figure 310a), the relative L, norm errors of BBL scheme are about 15% larger than those of Noble et al.'s scheme. Since Noble et al.'s scheme does not contain boundary conditioninduced error, this 15% difference in error in the BBL scheme results from the boundary condition treatment. Thus comparing the results using Noble et al.'s scheme and the BBL scheme, it can be inferred that in the presence of the strong velocity profile variation the truncation error contributes a significant part of the overall error. Because a substantial part of the overall error is from the truncation error as opposed to the boundary condition error with highly variable viscosity, the truncation error associated with the variable viscosity thus deserves close attention. However, the truncation error could have very complex form, even for the problems with constant viscosity (Holdych et al. 2004). It is recognized that there is a close relation between the LBE and the finite difference form of the momentum transport represented by the NavierStokes equation (He et al. 1997; He et al. 2002; Mei et al. 1998). For example, He et al.(1997) showed that the square lattice formulations for the particle distribution functions in a 2D pressure driven channel flow with constant viscosity, after averaging, leads to a second order, centraldifference formulation for the axial velocity. For the present channel flow problem with the large variation in viscosity, such derivation could not be easily obtained. However, if the relative L2norm error of the finite difference method still behaves similarly to that of the LBE, it is reasonable to expect that the truncation error behavior of the LBE is similar to that of the finite differencebased macroscopic model (in the present case, the NavierStokes equation). This hypothesis will be assessed first. In the NavierStokes model, the governing equation (Eq.(2.61)) is discretized by standard central difference scheme (with Ay = 1 lattice unit): 0 = d 1(A)2 }+1/2 (+1 +) v1/2 (u 1). (2.70) To compare the errors of the finite difference and LBE schemes, two finite difference solutions on different grid arrangements are obtained for the velocity profile. The first finite difference solution has the same grid arrangement as in the hexagonal lattice so that the first fluid node is one full mesh away from the wall. The second finite difference solution has the same grid arrangement as in the square lattice with the BBL scheme so that the first fluid node locates half a mesh away from the wall. This second finite difference solution requires an approximation for the velocity condition at the walls; a linear extrapolation is used in conjunction with the noslip condition at the wall. For completeness, a second order extrapolation is also used to approximate the derivative at the wall in solving Eq. (2.61); the error is consistently larger than the linear extrapolation and hence the results are not presented. The viscosity distribution of Case I is chosen for the finite difference computations for Case I gives sharper nearwall velocity gradients than Case II. The results of the relative L2 norm errors over a range of h=0.004 to 0.1 from these two finite difference solutions are shown in Figure 310b. The results labeled as "FDH" refer to the finite difference solution obtained on the hexagonal grid and the results labeled as "FDS" refer to the finite difference solution obtained on the square grid. For small values of h (h<0.01), the O(h2) asymptotic behavior in error is clearly visible, which is similar to the LBE cases. This asymptotic error behavior is expected because the LBE scheme has second order accuracy and the boundary conditions are, depending on the specific scheme chosen, either exact or second order accurate, and the finite difference scheme with the central difference discretization gives global second order accuracy. As h increases close to or greater than 0.01, both the finite difference and the LBE errors increase faster than O(h2). Since the velocity profile given by the exact solution has a thin layer of thickness 0= 0.0102, the error starts to increase more rapidly for h>0.01 when the resolution of the thin layer becomes inadequate. For 0.01 the numerical resolution is insufficient, all the error curves are quite close, indicating that the need for a better resolution of the boundary layer exceeds the need for a better boundary condition. One also notes that when h exceeds the boundary layer thickness, truncation errors associated with high order velocity gradients may not be small. Comparing the relative L2 norm errors of the finite difference solution with those of the LBE solution in Figure 31Ob, it is seen that the relative L2 norm error of"FDH" follows that of the LBE with Noble et al.'s scheme closely, and the relative L2 norm error "FDS" follows that of the LBE with the BBL scheme closely. This consistent behavior of the relative L2 norm errors between both the LBE method and the finite difference method suggests that one could obtain an insight on the truncation error of the LBE schemes by studying the truncation error of its finite difference counterpart. For both LBE and finite difference schemes, their modified equations associated with the corresponding truncation errors can be represented as L(ELBE) = T.E. of the LBE method (2.71) L(EFD) = T.E. of the finite difference method (2.72) where L is the differential operator d vota ELBE = t uLBE and ED = u,t u,. Using dx dx Taylor series expansion the truncation error in Eq. (2.72) is 1 04U _I VtotI_ 2 tota 2 I V total C2 T.E.= v +1 a _+ a (A)2 +H.O.T.(2.73) 12 o 4 6 ao] o]3 8 a]2 Or2 24 r73 ar] Although the truncation error of the LBE is unknown, comparing the solution errors, EFD and ELBE, on the lefthandsides of Eq. (2.71) and (2.72) can offer insight into the truncation errors on the righthandsides of these equations while this indirect comparison avoids the tedious derivation of the truncation error of the LBE. The value of EFD can be computed by subtracting uFD from the exact solution after the velocity uFD is solved from Eq. (2.70). To highlight the behavior of the leading term of the truncation error, the values of EFD are approximated by ED obtained from solving Eq. (2.72) with small h. For ELBE it can only be determined directly by subtracting ULBE from the exact solution. Eq. (2.72) can be solved, with the boundary conditions ED (0) = ED (1) = 0, by using the same central difference scheme given by Eq. (2.70), and replacing the entire T.E. by the leading term of the T.E. provided that the resolution is sufficient and the highorder terms in the T.E. can be neglected. Caution must be exercised in interpreting the results from the numerical solution of Eq. (2.72). In this equation, even if the leading term of the T.E. on the right hand side can be analytically evaluated, the variation of the viscosity on the left hand side is the source of the T.E. in the first place. When solving for E(77) the effect of the T.E. associated with the ODE (Eq. (2.72)) now is further compounded by the variation of the viscosity. The finite difference computation for Eq. (2.72) is carried out on a grid with the same grid arrangement and size as in the hexagonal lattice case so that the exact velocity boundary condition can be used. For Case I, the thin wall layer has a length scale 0=0.0102 and can be adequately resolved using Aq = 1/H = 1/200, which is confirmed by the comparison between ED and E in Figure 311. The curve representing ED is almost on top of the curve representing E Figure 311 also shows the variation of LBE. across the channel. Overall, E from Eq. (2.72) is smaller than LBEL across the whole channel which implies that these two methods indeed have quantitatively different truncation errors. However if we adjust the scale for E* by a factor of 2, these two curves lie almost on top of each other, as shown in the inset of Figure 311. The matching of the error curves between two solutions demonstrate that the LBE truncation error behaves very similarly to that observed in the finite difference scheme in the presence of a strong variation in the viscosity. This is in agreement with the observation on the relative L2 norm error behavior shown in Figure 310b. Based on the foregoing discussions, it is seen that BBL scheme performs similarly as Nobel's scheme in the presence of strong velocity profile variation. In view of the simplicity and flexibility of the BBL scheme over the Noble et al.'s scheme and comparable performance in accuracy, it is attractive to use the BBL scheme (or its extended version such as Bouzidi et al.(2001) and Yu et al. (2003a)) for handling problems of substantial velocity profile variations caused either by complex geometry or variable viscosity. For complex 3D flows with curved boundaries, NavierStokes solvers currently have more flexibility on grid arrangement. For example, bodyfitted coordinates and grid stretching can be more easily implemented to improve grid resolution near boundaries than in the LBE solvers. Although the recent developments in LBE method such as multiblock (Yu et al. 2003) and composite grid (Fillippova and Hanel 1998) techniques can alleviate the difficulty in grid arrangement in LBE simulations to certain extent, further research efforts in LBE are needed. As the present study has indicated the boundary condition induced error may not be dominant, there is a great potential in improving the overall capability of the LBE solvers for complex flows with very strong velocity variation by focusing future research efforts on extending the grid flexibility of LBE schemes. 3.3 Summary and Conclusion In this chapter, the error behavior of the lattice Boltzmann equation (LBE) method for a flow with strong variation in viscosity is studied. The variable viscosity in the LBE method is modeled through a variable relaxation time. Solutions for a laminar channel flow with a specified variable viscosity are obtained using both the LBE method and the finite difference method to examine the truncation error behavior of the LBE method for flows with strong varying viscosity. The effect of the boundary condition error of the bouncebackonlink (BBL) scheme on the overall error is investigated via the comparison of the error behavior of the BBL scheme and that of Noble et al.'s scheme. The results show that with rapid viscosity variation the boundary condition error of the BBL scheme does not induce noticeable, additional errors, and the overall error of such flows is dominated by the truncation error itself. The results also show that in the presence of strong variable viscosity the truncation error behavior of the LBE solution is consistent with that of the finite difference solution to NavierStokes solution. Xf Xf Wall Figure 31. Boundary nodes and their neighbors using the square lattice 3 \ / 2 2 Fluid (f) 1 Boundary \Wall 5 6 (w) Figure 32. A boundary cell using the hexagonal (FHP) lattice (Noble et al. 1995) 20 40 60 80 100 . Absolute L2 norm errors of LBE wi channel flow with constant viscosity th Noble's scheme for fullydeveloped laminar eh e_ ee e. e,+ Rx a x 10 10 10 10 10 10 10 10 10 Figure 33 40 60 80 100 H Absolute error of a fullydeveloped channel flow using Inamuro et. al.'s scheme 4 3 7 8 A 2D 9velocity lattice (D2Q9) model 6 Figure 34. 6E11  5E11  4E11  .2 3E11  0 2j 2E 11 Figure 35. 20 I I I I I Figure 36. Two set of viscosity distributions used in this study 0.8 se ii. 025  0.6 0 2 Case I  Case II 015 tan terms S in mEq (3 35) (Case ) 0. , 0.4  00 5  0 0 0002 0004 / . 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 u Figure 37. The exact velocity profiles of the channel flows with different boundary layer thicknesses due to different viscosity distributions. The parameters are: (v;8; O) = (0.004167; 0.0005; 0.0102) for Case I and (vo0;;0) =(0.008333;0.002;0.0289) for Case II 1 Outlet: Periodic boundary condition Figure 38. Square lattice distribution in channel flow simulation j=Ny j=Ny1 Inlet: Periodic boundary condition j=2 j=i 0 0006 Exact Exact u Hex 00005 u_Square 0 0004  o 00003  0 0002000 0 01 02 03 04 05 11 Figure 39. Comparison of the LBE velocity profiles using square lattice (with bounceback on the link boundary condition) and hexagonal lattice (with Nobel et al's exact boundary condition) with the exact solution at H=50 lattice units. p 0 C BBL scheme (Case I) Noble et al.' s scheme (Case I) BBL scheme (Case II) Noble et al.' s scheme (Case II) I I I I I I I I I I I I II li li 0.02 h 0.04 0.06 0.08 0.1 0.12 a 10' ** 21O 102 10  2  FDS (Case I) 1 FDH (Case I)  BBL scheme (Case I) 10  Noble et al.' s scheme (Case I) 10 I I i I I I I 0.02 0.04 0.06 0.08 0.1 0.12 h b Figure 310. Dependence of the relative L2 norm error on the lattice size h in the fully developed channel flow with variable viscosity. The viscosity parameters are: (3,0) = (0.0005,0.0102) for Case I, (0.002,0.0289) for Case II. (a) The LBE with the boundary conditions of Noble et al's scheme and BBL scheme for both Case I and Case II. (b) The finite difference and the LBE with boundary conditions of Noble et al's scheme and BBL scheme for Case I. 8E06 E E G E 6E06 4E06 6E06  E_ 3E06 4E06 :4E06 2E06 2E06 E06 2E06 0 0 0 002 004 006 008 01 012 014 '1 0 0.1 0.2 0.3 0.4 0.5 T] Figure 311. Comparison of ELBE = c, uLBE with E = ue,, uF for H=200, and 8 = 0.0005,0 = 0.0102, E*, is the numerical approximation of ED obtained from solving Eq. (2.72). CHAPTER 4 LBE METHOD FOR IMMISCIBLE TWOPHASE FLOW COMPUTATION The characteristic feature of immiscible two/multiphase flows is the interphase boundary (interface) between different phases. The structure of any interphase boundaries is on the mesoscopic length scale level (Pismen 2001). The steep gradients of material properties in the normal direction of an interface result from the fact that the normal length scale approaches molecular level, while the length scale along the interface is still on the macroscopic level (Rowlinson et al. 1982a; Carey 1992). As described in Chapter 1 these scale disparities make the computations more difficult. In this chapter, a LBE immiscible twophase model is developed based on the framework of He et al.'s model (1999). Several numerical techniques are implemented to the numerical stability so that twophase flows with high liquidgas density ratio can be simulated. This chapter begins with an overview of the immiscible twophase flow computations and a literature review on the LBE twophase models, followed by an introduction to He et al.'s model (He et al. 1999). The RayleighTaylor instability is simulated and computational results are compared with analytical ones to validate the code. After this, the numerical issues in LeeLin's implicit LBE model and He et al.'s model are addressed. A filter based LBE model along with a new surface force formulation and a volumecorrection step is proposed. The stability and accuracy of this new model are tested by computing flows around a stationary bubbles, capillary waves and rising bubbles. The results show that the present model is capable of handling the flows with large density ratio up to 0(102). 4.1 Overview of Immiscible TwoPhase Flow Computation For macroscopic computations of immiscible twophase flow the numerical methods can be classified into two categories. One is the interfacetracking method, in which interfaces have to be tracked numerically to separate materials of two different phases or properties (Shyy et al. 1996). The other is the continuous interface method in which interfaces are not tracked explicitly (Anderson et al. 1998). Two approaches exist in terms of the grids used for computations. One is the moving grid (Lagrangian) approach. The other is the fixed grid (Eulerian) approach (Shyy et al. 1996; Shyy 2004). These two approaches can also be combined to form a mixed approach that is called the EulerianLagrangian approach. In the Lagrangian approach the interface is tracked explicitly by a bodyfitted grid that deforms with the interface. No modeling is necessary to define the interface or its effect on the flow field. However, the bodyfitted grid has to be regenerated as the interface deforms. This requires a set of equations that need to be solved for the grid regeneration at each time step. Moreover, Lagrangian approach is not easy to use when fluid flows involve complex topological changes of interfaces. The Eulerian approach is more robust in general but it needs elaborate procedures to deduce the interface location based on the volume fraction information. Typical Eulerian methods are the levelset (LS) method and the volume of fluid (VOF) method. Interface is often constructed after the field solution obtained. Thus the interface construction decouples from the field equation solver. Since the grid is fixed in the Eulerian approach, the grid generation required by the Lagrangian approach is obviated. Another advantage of the Eulerian approach over the Lagrangian approach is that the topological changes of interfaces are automatically taken into account in the Eulerian approach, which makes the Eulerian approach more suited for flows with complex interface deformations. However, the Eulerian approach suffers several disadvantages. The interface shape is smeared off. If high resolution information of interface is desired, an adaptive local grid refinement may be needed, which is another complicated process. Thus, some physics may be lost if the complicated grid refinement is not available, such as the uniqueness of the shape interpretation, and continuity and smoothness of the interface (Shyy 2004). The EulerianLagrangian approach combines the advantages of both Lagrangian and Eulerian approaches. In this approach, the grid is fixed as the Eulerian approach. Grid refinement is not needed. The interface is explicitly tracked by some marker points which are the Lagrangian components. This explicit interfacetracking can provide detailed interface information as the Lagrangian approach (Shyy et al. 1996). The Eulerian, Lagrangian and EulerianLagrangian approaches are numerical methods for twophase flow computations with interfaces. Physically, the interface can be modeled as a finite thickness line (2D) and surface (3D), or a zero thickness line (2D) and surface (3D). They are called the continuous interface method (CIM) and the sharp interface method (SIM), respectively. At the macroscopic modeling level, the SIM is consistent with the concept of the continuum mechanics (Shyy 2004). Interface is tracked explicitly with some computational cells which are cut by the interface. Thus interface smearing is not involved. In contrast to the SIM, the CIM does not track interfaces explicitly. Of these two methods, the SIM and CIM, if detailed interface information is desired, the SIM is a better choice since it can provide a secondorder accuracy while the CIM can only provide a firstorder accuracy in general. However, the CIM has some advantages which the SIM lacks. For instance, the CIM can be easily adopted for three dimensional flows with topological changes of interfaces; the CIM is particularly useful for flows with phase changes; and the CIM is especially appropriate for some problems that are currently tough for the SIM, such as contact line dynamics (Pismen 2001; Lee et al. 2005). More detailed information for moving interface computations based on the macroscopic computation point of view can be found elsewhere (Shyy et al. 1996, Shyy 2004) The method of lattice Boltzmann equation is essentially a kinetic equation solver over mesoscopic length scale (Succi 2001). Physically the essential length scale of this method is closer to microscopic level than macroscopic continuumbased methods. Thus it is relatively easier for the LBE method to incorporate microscopic modeling for twophase flows with interfaces. As a result, the LBE method has numerous advantages that the molecular dynamics method has, which are especially useful for the simulations of phase interfaces of nonideal gases (Shan, et al, 1993; He, et al, 1999) or binary fluids and near wall treatment at microfluid level (Nie, et al, 2002). Phase segregation and interfacial dynamics can be achieved naturally by incorporating intermolecular interactions. Because the interface spreads over several lattices in the LBE method (He et al. 1999) most of the twophase LBE methods can be considered as CIM. As one of the continuous interface methods, the LBE method has all the advantages of the CIM for twophase flows. However, the nature of the modeling on the mesoscopic level may allow the LBE method to incorporate the intermolecular interaction more easily than the conventional CIM. Therefore, the LBE method may be a good alternative for twophase flows simulations (Fang et al. 2001; Lee et al. 2005). 4.2 Literature Review on LBE Method for TwoPhase Flow Computation The first latticebased twophase flow model was the lattice gas color particle model proposed by Rothman and Keller (Rothman et al. 1988). This model has not been widely used in practice since it suffers from the defects of the lattice gas method, such as the computational noise. The LBE method overcomes the natural defects of the lattice gas method. Gunstensen et al. (Gunstensen et al. 1991; Gunstensen et al. 1993) proposed the first LBE multiphase model developed from Rothman and Keller's lattice gas model. Although Gunstensen et al's model possesses the essential features of the LBE method, i.e., the Galilean invariance and statistical noise free, it cannot simulate multiphase flows with different densities and viscosities. Grunau et al. (Grunau et al. 1993) improved Gunstensen et al.'s model by using a singletime relaxation approximation and a proper equilibrium distribution function. This model can be used for flows with different densities and viscosities. Since these models were developed from Rothman and Keller's pioneering lattice gas color particle model, all require a "recolor" step to maintain the interfaces. Shan and Chen proposed a new twophase LBE model hereinafter referred as (SC model) in which the nonlocal microscopic particle interaction is incorporated (Shan et al. 1993; Hou et al. 1997). The molecular interaction introduced into the ideal gas LBE model is represented by an interaction potential which models the multiphase separation and dynamics. Although the SC model has performed much better for the multiphase flows than the RK model (Hou et al. 1997), it possesses some anisotropy (He et al. 1998). The freeenergy LBE model developed by Swift et al. (1995) uses an equilibrium distribution which makes the pressure tensor consistent with that of the freeenergy function of nonuniform fluids. However the Galilean invariance is not satisfied. Although the aforementioned models are based on different physical considerations, He et al. pointed out that they all originate from the kinetic theory which can be represented by the continuous Boltzmann equation with certain approximations (He et al. 1998). Subsequently He et al. proposed an improved LBE scheme using the singletimerelaxation approximation to simulate of multiphase flows in the incompressible limit (He et al. 1999a). In this new model, two distribution functions corresponding to two evolution equations are employed. One acts as an index function to track interfaces implicitly between different phases. The other is an evolution equation for pressure distribution which satisfies the mass and momentum conservation on the macroscopic level. Molecular interactions, such as the molecular exclusion volume effect and the intermolecular attraction, are incorporated into this model. These interactions maintain the phase separation by the mechanical instability on the supernodal curve of the phase diagram (He et al. 1999). Compared to the SC model, this thermodynamicsbased concept has much more flexibility for maintaining the phase separation. Furthermore, the numerical stability of He et al.'s model is improved by reducing the numerical errors from 0(1) to 0(u) during the calculation of the molecular interactions. One common deficiency in all the aforementioned models is that they cannot be used for the multiphase flows with large density ratio (Lee et al. 2005). RK model, SC model, and Swift model can only be used for the twophase flows with density ratio less than 2. He et al.'s improved model brings the density ratio to around 10. Chen et al. (1998) discretized the Boltzmann equation with a total variation diminishing scheme. In Chen et al.'s method, the secondorder RungeKutta timematching was used for the discretized Boltzmann equation (DBE), while the collision term and the intermolecular interaction term were treated explicitly as source terms. Although the use of TVD scheme improved the stability of the LBE method for the large density ratio case, the accuracy is reduced in the presence of large gradients because of the additional numerical dissipation of the TVD scheme. Inamuro et al. (2004) proposed another LBE model for the twophase flows with large density ratio. They also used two particle distribution function equations: one as an index function to indicate the phase separation, the other for the momentum evolution but without the pressure gradient. Velocity and pressure are coupled by an additional pressure Poisson equation to satisfy the incompressibility requirement. The procedure is analogous to the fraction step method (Ferziger 2003) which is used to solve the macroscopic NavierStokes equation. It is well known that it is time consuming to solve the Poisson equation (Shyy, 1997c). Although Inamuro et al.'s were able to simulate multiphase flows with large density ratio using their scheme, such as rising bubbles in liquid, the natural simplicity of LBE method is lost due to the necessity of solving the pressure Poisson equation. Lee and Lin developed an implicit twophase LBE model which can simulate large density/viscosity ratio problems without using the iteration procedure (Lee and Lin 2005). Lee and Lin's model is similar to He et al.'s model but with a different pressure updating process. In He et al.'s model, the pressure in calculating the intermolecular interaction term is updated from equation of state (EOS). Lee and Lin argued that this pressure updating process can lead to the spurious pressure fluctuations at the interfaces due to the EOS. Lee and Lin proposed a pressure evolution equation which can overcome this difficulty and at the same time allow for large density difference across the phase interfaces when phase change occurs due to pressurization and depressurization (Lee et al. 2003). Using this model, Lee and Lin simulated 1D advection equation with a source term, a stationary droplet, an oscillating droplet and a droplet splashing on a thin film at a density ratio of 1000 with varying Reynolds numbers. However, very little information on the details of the flow field near the large gradient interfacial region was reported. In the following sections, He et al.'s LBE model for twophase flows will be reviewed. It will be applied to the RayleighTaylor instability problem for a code validation. The interface behavior of Lee and Lin's implicit LBE is then studied. To maintain the interface thickness as sharp as the initial interface thickness and improve numerical stability, a filterbased twophase LBE model is developed based on He et al.'s model. Stationary liquid flow around a static bubble, the capillary wave caused by two superposed fluids with same viscosity, and flow induced by a rising deforming bubble are computed at density ratio of 0(102) to asses the efficacy of the proposed twophase LBE model. 4.3 He et al.'s Isothermal LBE Model for TwoPhase Flow 4.3.1 The Boltzmann Equation for NonIdea Fluids The Boltzmann equation for nonidealfluid flows intermolecular interaction force is (He et al. 1998) f Vf f feq (+ u) (F +G) fq (2.1) +( Vf_= + Ieq (2.1) 8t A pRT The singlerelaxationtime BGK model is used in this equation. The derivation of Eq. (2.1) is based on the assumption that the distribution function gradient with respect to the particle velocity is equal to the equilibrium distribution function gradient with respect to the particle velocity (He et al. 1998). Compared to the standard LBE model for ideal gas flows, the last term on the RHS of Eq. (2.1) represents the effects of the intermolecular interaction F and gravity G. In Eq.(2.1), R is the gas constant and 1 is the relaxation time. The Maxwellian distribution for an equilibrium state is feq = ) ex (e R u)2 (2.2) (27RT)D/2 2RT where D is the dimension of the space. The intermolecular attractive force originates from the van der Waals theory (Rowlinson et al. 1982b). This intermolecular force can be expressed as F,, = pV 2ap+ AV2p) (2.3) The parameters a and K result from the intermolecular pairwise potential U 1I a = U (r)dr, K = r2U (r)dr, (2.4) 2 r>d 6 r>d where d is the effective molecular diameter, dr is the effective differential volume. Besides the intermolecular force, He et al. also included the exclusionvolume effect into the intermolecular interaction FE, = bp2RTzVln(p2z) (2.5) which accounts for the extra volume effect of nonideal fluids (Chapman et al. 1970). The 2 d3 parameter b is given by b = 3 where m is the single molecular mass. The parameter, in (2.5) 3m is the increase in the collision probability due to the increasing fluid density. This parameter can be expressed as a polynomial in terms of the product of the density p and the parameter b (p)= l+bp +0.2869(bp)2 +0.1103(bp)3 +... (2.6) The overall intermolecular interaction is thus the sum of the intermolecular force and the exclusionvolume effect F = FM, + FV (2.7) The intermolecular interaction can be recast as F = V V + F (2.8) where F, is associated with the surface tension, F = KpVV2p The potential term VVi is associated with the pressure gradient as = p pRT (2.9) Thus the term VVi depends on the equation of state (EOS) and it in turn determines the phase segregation. From the van der Waals theory, if the temperature is lower than the critical point of a fluid, the phase segregation can be generated by the molecular attraction. The critical temperature is determined by the equation of state. For the van der Waals equation of state pRT ap2 (2.10) 1bp the critical temperature is T =a If the hard sphere model is used for the repulsion 27bR interaction between molecules, the CarnahanStarling equation of state is obtained (Carnahan and Starling 1969, 1972) 1+bp +(bp 2bp3 p = pRT 4ap2 (2.11) (1^ and the critical temperature is Tr = 0.3773a, where b is the covolume of spheres. During the bR computation, the parameters a and b should be selected to satisfy the condition that the temperature is lower than the critical temperature. He et al. pointed out that vIy(p) in Eq. (2.8) is usually very large at the interfaces. As a result, the calculation of VI,(p) can involve to a substantial numerical error, which lends the numerical scheme unstable. In order to reduce this sensitivity, He et al. introduced another variable g to replace the distribution function g= fRT + (p)F (0) (2.12) where F(0)is the value of F(u) evaluated at u=0, and F(u) (exp ( ) The (2T/1/)' 2RT derivative of g is then g g =Vg =RT +F(O)) (2.13) Dt 9t Dt Dt Since y(p) is a function of density, we can rewrite DV Dt DyI V/ + = +u.V+V. u.V = d +(yu).Vi (2.14) Dt at at dt For incompressible flow dy =y p was then treated as zero in the wok of He et al (1999). dt dp dt Substituting Eqs. (2.1), (2.12) and (2.14) into (2.13) the evolution equation of g can be written as g+ +.Vg gg +( u){[F(u)(F1 +G)(F(u)F(O))VV(p)] (2.15) The equilibrium distribution function of g is geq = pRTF(u)+V(p)F(0). (2.16) From g, the macroscopic variables pressure p, instead of the density p, and velocity u can be calculated p= fgdf, pRTu = fgd. (2.17) Note the function V/ is calculated from the hydrodynamic pressure as (p)= p, pRT (2.18) where the hydrodynamic pressure p, is obtained from Eq. (2.17), instead of equation of state (Eq. (2.11)) as was typically the case for single phase flow. The sensitivity of V/(p) in Eq. (2.15) is alleviated by multiplying (F(u)F(0)) which is proportional to the Mach number and is small in the incompressible limit. For multiphase flows, besides pressure and velocity which are given by (2.17), density is another important variable which describes the different phase distributions. Density is only a single function of pressure for isothermal singlephase flows, but it cannot be uniquely determined by the pressure alone for isothermal multiphase flows. Therefore, like other diffusive interface methods for multiphase incompressible flows, an index function evolution equation is needed for tracking different phases and then the density or composition of fluids could be obtained from this index function. In the LBE models for twophase flows, Eq. (2.1) can serve as the index function evolution equation. Since the gravity does not affect the mass conservation, He et al. dropped the gravity in the index function evolution equation. It is given by, = +Vf= +( )V )feq (2.19) Dt 8t 2 pRT It is noted that / is the function of the index function 0, instead of density p Thus it is defined as (0) = p, ORT (2.20) where the thermodynamic pressure p, is the thermodynamic pressure calculated from the equation of state Eq. (2.11), while the density is replaced by the index function 0 1+ + p, = ORT 3 a02 (2.21) 4 To recapitulate, the evolution equations forf and g are D.f Vf f + (U) VI(0) feq (2.22) Dt 8t 2 pRT g+ Vg =  +( u).[F(u)(F1 +G)(F(u)F(O))VV((p)] (2.23) where V/ () = p, ORT, and /(p) = p, pRT. The equilibrium distributions forfand g are f q = F (u) (2.24) geq = pRTF(u)+V(p)F(O), (2.25) where F(u)= 1 exp L u )2. The macroscopic variables can be obtained from the (2RT)D 2 2RT moments off and g S= ffde, p = gde, pRTu = f gde. (2.26) Density and viscosity can be calculated by the linear interpolation from 0 P () = P, + P,), (2.27) where p, and p, are the material densities of light and heavy fluids, respectively; v, and V, are the kinematical viscosities of light and heavy fluids. The above formulation is equivalent to the following macroscopic equations (Zhang et al. 2000): + .(u)= V p(p) Vp()V (2.29) 1 p +V.u=0 (2.30) pRT Ot S +]t(V) =Vp+ V[pv(Vu+uV)]+V2p+G (2.31) 4.3.2 Lattice Boltzmann Scheme for Multiphase Flow in the Near Incompressible Limit The lattice Boltzmann equation can be obtained by discretizing the Boltzmann equations described in the previous section. The discretization procedure is same as given in Chapter 3 for single phase flows. After discretization in the phase space, the discretized Boltzmann equations are a a+e .V (eu) ) (2.32) Qt A pRT +e Vg a q +(e u).[F (u)(F +G)(Fa (u)F (0))VV(p)] (2.33) To improve the accuracy and maintain the explicit computational scheme, He et al. further introduced two new variables, which are 7 = (e u)V (0) (u),, (2.34) 2RT a ga (e. u).[F (u)(F] +G)(F,(u)Fa(O))Vw(p)]a1 (2.35) 2g 7c The new evolution equations in terms of these two new variables are (x+e?, ,t + ,) (x,t) ( (x,t) {(xxt)(2r_ 1) 1){( u).y() F (2.36) T 2z RT X +(x +e6,,t + 8 )' (x,t) = t) _'q(xt) +(2z (e, u). [F, (u)(F + G) (F (u) F (0))VV((p)] (2.37) where is still the nondimensional relaxation time, r=/,. The macroscopic variables can be calculated in terms of the new distribution functions as S=Z (2.38) p = uVV/(p),, (2.39) pRTu= e, + RT(F +G)6,. (2.40) In the D2Q9 model, the kinetic viscosity is independent of surface force and is related to the non dimensional relaxation time as v = ( 0.5)RT8,. (2.41) Zhang et al. (2000), and McCracken et al. (2005) have used the following integral relationship to analytically relate surface tension a with the coefficient ic : cr=kl(a) = k dz (2.42) where z is a direction normal to a flat interface. 4.4 Code Validation: RayleighTaylor Instability A twodimensional D2Q9 LBE (Eq. (2.36)) code has been developed to solve Eqns (2.36) (2.41). To validate the code the RayleighTaylor instability is simulated and the results are compared with analytical solutions. The RayleighTaylor instability problem is a very challenging computational task because it possesses great complexities including strong non linearity associated with growth of the secondary KelvinHelmholtz instability and turbulent mixing at later stages. On the other hand, the RayleighTaylor instability has been systematically studied in the past (Chandrasekhar 1981), and the results from previous studies can be used to validate the present LBE simulation. 4.4.1 Linear Analysis of RayleighTaylor Instability The RayleighTaylor (RT) instability occurs when a heavier fluid of density p2 superposes over a lighter fluid of density p and the lighter one accelerates the heavier one in the direction normal to the plane interface between these two fluids. The heavier fluid moves downward producing spikes and lighter fluid grows upward producing bubbles. The linear analysis is usually used to predict the linear instability. In this approach, the perturbation analysis is used to linearize the nonlinear convection terms. The resulting pressure and density are functions of the vertical coordinate, say z, only (Chandrasekhar 1981). Chandrasekhar gave the relation between the growth rate of the disturbance as a function of wave number: 2iAt+ k2 +1 }(f +fq,k)4kf +4k(fl v2)[(f2q flq2)+k(f, f)] (2.43) 4k3 + (fv f)2) (qk)(q k)=0, The Atwood number is defined as At= 2 p (2.44) P1 + p2 where a is growth rate, g is gravity, k iswave number, v is viscosity, and P 2 1 Aq2 2 +a f= A, f2k2+, q2 = k2 . PA 2 Al + A V, V 2 4.4.2 RayleighTaylor Instability: LBE Results The twodimensional RayleighTaylor instability with singlemode is simulated using the same physical parameters as in He et al.'s work (He et al. 1999a). Although He et al. presented the multimode RayleighTaylor instability results, the singlemode RayleighTaylor instability is preferred to serve as the benchmark for the present numerical simulation. The computational domain is bounded by the top and bottom walls on which the noslip boundary condition is applied. The periodic boundary condition is applied on the two lateral boundaries. Surface tension is neglected in this simulation and the kinetic viscosities of the two fluids are assumed to be same. The function V/(0)= p,, RT is calculated from the pressure which satisfies the CarnahanStarling equation of state, i.e., Eq. (2.11). In the CarnahanStarling equation, the parameter a is chosen to be bc2, which induces sufficient molecular interaction to separate different phases. He et al. (1999), Zhang et al. (2000) and McCracken et al. (2005) reported that a =bc2= 4 is sufficient to separate phases. With these values for a and b, the function I(a) in Eq. (2.42), that is I(a)= f (c \dz, can be numerically determined for Canahan Starling fluid (Zhang et al. 2000) as 0."1518(aa5 (2.45) I(a)= 05 (2.45) l+3.385(aa )5 where a, = 3.53374 is the critical value of CarnahanStarling equation of state, below which a fluid cannot separate into different phases. With these parameters, the range of the index function can be theoretically or numerically determined by the equation of state. For CanahanStarling equation of state, Zhang et al. (2000) obtained the range of the index function as 0.023810.2508. In order to make meaningful comparisons two dimensionless numbers, Reynolds number and Atwood number, are used for this problem. Reynolds number is defined as 