【精品论文】Threedimensional finitedifference lattice Boltzmann model and its application to invi.doc

精品论文Three-dimensional finite-difference lattice Boltzmann model and its application to inviscid compressible flows with shock waves5HE Yaling, LIU Qing, LI Qing(Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and PowerEngineering, Xian Jiaotong University, Xian 710049)Abstract: In this paper, a three-dimensional (3D) finite-difference lattice Boltzmann model for simulating compressible flows with shock waves is developed in the framework of the10double-distribution-function approach. In the model, a density distribution function is adopted to model the flow field, while a total energy distribution function is adopted to simulate the temperature field.The discrete equilibrium density and total energy distribution functions are derived from the Hermiteexpansions of the continuous equilibrium distribution functions. The discrete velocity set is obtained by choosing the abscissae of a suitable Gauss-Hermite quadrature with sufficient accuracy. In order to15capture the shock waves in compressible flows and improve the numerical accuracy and stability, an implicit-explicit finite-difference numerical technique based on the total variation diminishing flux limitation is introduced to solve the discrete kinetic equations. The model is tested by numerical simulations of some typical compressible flows with shock waves ranging from 1D to 3D. The numerical results are found to be in good agreement with the analytical solutions and/or other20numerical results reported in the literature.Keywords: compressible fluid dynamics; lattice Boltzmann model; Hermite expansion; finite- difference; shock wave; three-dimensional0Introduction25The lattice Boltzmann (LB) method, which originates from the lattice gas automata (LGA) method1, has been developed into a promising numerical method for simulating complex fluid flows and modeling complex physics in fluids2-4. Different from the conventional numerical methods, which are based on the discretization of macroscopic governing equations, the LB method is based on the kinetic theory and it simulates fluid flows by tracking the evolutions of the30particle distribution functions. Owing to its kinetic nature and distinctive computational features, such as the simple form of the governing equation, the easy implementation of complex boundary conditions, and the inherent parallelizability on multiple processors5, the LB method is a very attractive representation of the nonlinear macroscopic systems, especially for those containing shock waves or contact discontinuities.35In the past decade or so, many brilliant efforts have been made in constructing LB models for the fully compressible Euler6-11 or Navier-Stokes (NS) equations12-19. Among these models, Sun et al.13-15 developed a locally adaptive LB model. In this model, a Kronecker delta function is adopted to replace the Maxwellian equilibrium distribution function as the continuous equilibrium distribution function, and the lattice velocities are chosen according to the local flow velocity and40internal energy. So it permits the mean flow to have high Mach number, but the relaxation time is fixed at 1 and the Prandtl number is equal to the specific-heat ratio, which may limit its application. Subsequently, Qu et al.11 developed a non-free-parameter LB method to construct equilibrium distribution functions for 2D problems. In this method, the conventional Maxwellian function is replaced by a circular function and a Lagrangian interpolation polynomial is then45constructed to discretize the circular function to the discrete velocity directions. However, thePrandtl number in Qu and co-workers model is also equal to the specific-heat ratio. Following theFoundations: The Ph.D Programs Foundation of Ministry of Education of China (No. 20110201110038).Brief author introduction:HE Yaling, (1963-),female,professor,new technologies in energy storage,numerical principle and its applications in fluid flow and heat transfer process. E-mail: yalinghe- 15 -non-free-parameter LB method11, Li et al.6 developed a 3D non-free-parameter LB model for inviscid compressible flows with shock waves. In Li and co-workers model, a spherical function, which satisfies the zeroth- through third-order moments of the Maxwellian distribution function in503D space, is introduced to replace the Maxwellian function as the continuous equilibrium distribution function.In the LB community, the determination of the equilibrium distribution function is a key issue. The equilibrium distribution function can be determined by several methods. The first method, which is the conventional method, is based on the Maxwellian equilibrium distribution55function. This method is widely used in the current LB community for it is easy to be accepted and implemented. The conventional method includes two approaches, i.e., the Taylor expansion approach20-22 and the Hermite expansion approach23. In the Taylor expansion approach, the equilibrium distribution function can be derived by applying the truncated Taylor series expansionto the exponential form of the Maxwellian function in terms of the local velocity u . However, the60Taylor expansion approach is limited to low- and moderate-Mach-number flows. In the Hermite expansion approach, the equilibrium distribution function can be determined by projecting theMaxwellian function onto the tensor Hermite polynomial basis in terms of the particle velocity and up to a certain order. As pointed out by Shan et al.23, the Hermite expansion approach allowsfor simulations of high-Mach-number flows. The second method has nothing to do with the65Maxwellian function. In this method, the equilibrium distribution function can be a Kronecker function13-15, a circular function11 or a spherical function6, which satisfies the needed statisticalrelations and is used to replace the Maxwellian function as the continuous equilibrium distribution function. The main advantage of the second method is that it can be used to simulate high-Mach-number flows.70In this paper, following the Hermite expansion approach23, we will present a 3Dfinite-difference LB model for simulating compressible flows with shock waves in the framework of the DDF LB approach, which can be viewed as an extension to some previous works23, 24. In the computation, like the earlier works6, an implicit-explicit (IMEX) finite-difference scheme based on the total variation diminishing (TVD) flux limitation is adopted to solve the discrete75kinetic equations in order to improve the numerical stability and capture the shock waves with a finite number of grid points. The rest of this paper is organized as follows. Sec. 1 briefly introduces the two-relaxation-time kinetic model. Sec. 1 presents the 3D finite-difference LB model in detail. In Sec. 1, numerical simulations are carried out for some typical compressible flows with shock waves to validate the proposed model. Finally, a brief conclusion is given in Sec.804.1Two-relaxation-time kinetic modelIn this section, we will briefly introduce the two-relaxation-time kinetic model17, 24, 25. Historically, the first DDF LB model including the compression work and the viscous heat dissipation was devised by He et al.22. The key point in He and co-workers model is the use of85two sets of distribution functions: the density distribution function to solve the flow field and the internal energy distribution function to solve the temperature field. This model has attracted much attention for its good numerical stability and the adjustability of the Prandtl number. Following He and co-workers model, Guo et al.25 developed a decoupling thermal LB model for low Mach number flows. In Guo and co-workers model, a total energy distribution function is introduced to90represent the total energy which enables the DDF LB models to be simpler and more efficiency as compared with He and co-workers DDF LB model. Subsequently, Li et al.17 developed a95100coupled DDF LB method for compressible flows by combining Guo and co-workers DDF LB approach with the multispeed approach. Recently, Li et al.24 developed a coupling LB model for simulating thermal flows on standard lattices in the framework of the DDF LB approach25. In this model, the discrete equilibrium density and total energy distribution functions are derived from the Hermite expansions of the continuous equilibrium distribution functions. Both the Prandtl number and the specific-heat ratio are adjustable in the coupling LB model.According to Refs. 17, 24, 25, the kinetic equations of the two-relaxation-time kinetic model are written as follows:f + f = 1 ( f f eq ) ,(1a)t fh + h = 1 (h heq ) + u ( f f eq ) ,(1b)t h h fwheref ( x, , t )is the density distribution function for particle moving with velocity atposition x and time t ,h ( x, , t )is the total energy distribution function; fand hare themomentum and total energy relaxation times, respectively; hf is a relaxation time related to f105and has hf = h f( f h ) ;f eq is the Maxwellian equilibrium distribution function andheqis the equilibrium total energy distribution function, which are defined as follows:f eq ( ; , u,T ) = ( u )2exp ,(2a) ( 2 RT )D 2 2RT 2h ( ; , ,T ) = + (b D ) RT ( u )2 exp eq u ,(2b)2 ( 2 RT )D 2 2RT where 2 = , u 2 = u u , D is the spatial dimension and R is the gas constant.110The macroscopic density , velocity u and total energy E are defined as follows: = f ( x, , t ) d , u = f ( x, , t ) d , E = h ( x, , t ) d (3)whereE = bRT2 + u 22 , and bis a constant related to the specific-heat ratio via = cpcv = (b + 2) b , in whichcp = (b + 2) R 2andcv = bR 2are the specific-heat coefficients115at constant pressure and volume, respectively.Through the Chapman-Enskog analysis17, 24, 25, 26, the following macroscopic equations can be obtained at the NS level: + ( u ) = 0 ,(4a)t ( u )+ ( uu ) = p + ,(4b)t ( E )+ ( E + p ) u = (T ) + ( u ) ,(4c)120where p = RTtis the pressure, = h (b + 2) 2 pRis the thermal conductivity. The viscousstress tensor is given by B = u + (u )T ( 2 D ) ( u) I + ( u) I ,(5)where = f pis the dynamic viscosity,B = 2 (1 D 1 b) f pis the bulk viscosity, andI is the unit tensor. The Prandtl number of the system is defined as Pr = f h .1252Three-dimensional finite-difference lattice Boltzmann model1302.1Hermite expansions of the continuous equilibrium distribution functioneqIn the LB community, the determination of the equilibrium distribution functions is a key issue. In this subsection, we will first specify the detailed forms of the equilibrium distribution functions based on the continuous equilibrium distribution functions given by Eqs. (2). FollowingShan and co-workers approach23, the Hermite expansions of the continuous equilibriumdistribution functionsf eqandheqcan be written as:ceqf eq = ( ,T ) 1 a( n ) H ( n) ( ) ,(6a)n = 0 n!cheq = ( ,T ) 1 b ( n) H ( n) ( ) ,(6b)wheren = 0 n!135 ( , T ) =1 2 exp (7)c( 2 RT )D 2 2RT eqcc ( n ) is the weight function, = RTcand Tcis the characteristic temperature. Note thataeq ,eqb( n) andH ( n) ( )are rank-n tensors and the products on the right-hand sides of Eqs. (6) denoteeqfull contraction. The expansion coefficientsa( n ) andb( n) are given by140The first few polynomials are23:( n ) a=eqb=( n) eq f eq H ( n ) ( ) d ,(8a) heq H ( n) ( ) d .(8b)eq H (0) ( ) = 1 , H (1) ( ) = ,H ( 2) ( ) = ,H (3) ( ) = .(9)aeq Using Eqs. (8a) and (9), the leading Hermite expansion coefficients of as follows23:f eqcan be obtained145( 0) eqeq= ,a (1) = u ,a ( 2) = uu + ( 1) , a (3) = uuu + ( 1) u,(10)whereu = uRTcand = T Tc . With Eqs. (6a), (9) and (10), the third-order Hermiteexpansion of the Maxwellianf eqis given explicitly as follows23:f eq ,3 (T ) = ( , T ) 1 + u + 1( u )2 u 2 + ( 1) ( 2 D ) 2c + u ( u )2 3u 2 + 3( 1) ( 2 D 2) ,(11)6150Using Eq. (8b) and (9), the leading Hermite expansion coefficients of as24:heqcan be obtainedbeq(0) eqeq= E ,b (1) = ( E + p ) u ,155b ( 2) = ( E + 2 p ) uu + p + ( 1) E .(12)With Eqs. (6b), (9) and (12), we can obtain the following second-order Hermite expansion ofheq 24:heq , 2 (T ) = ( , T ) E1 + u + =1 ( u )2 u 2 + ( 1) ( 2 D )2c + ( , T ) p u + ( u )2 u 2 + ( 2 D ) .(13)c 2It can be proved thatf eq ,3 (T ) satisfies the zeroth- through third-order velocity moments ofthe Maxwellian distribution functionf eq 24:160 f eq ,3 (T )d f eq d = ,(14a)eq ,3 ( )eq i fT d i f d = ui ,(14b)eq ,3 ( )eq i j fT d i j f d = ui u j + pij ,(14c)eq ,3 ( )eq () i j k fT d i j k f d = ui u j uk + p ui jk + u jik + uk ij.(14d)165It can be verified that of heq 24:heq , 2 (T )satisfies the zeroth- through second-order velocity moments heq, 2 (T )d heq d = E ,(15a)eq ,2 ( )eq () i hT d i h d =p + E ui ,(15b) heq , 2 (T )d heq d = ( E + 2 p ) u u+ p ( E + RT ).(15c) i j i j i j ij1702.2Discrete velocity m