Kinetic lattice methods are a very attractive representation of nonlinear macroscopic systems because of their inherent parallelizability on multiple processors and their avoidance of the nonlinear convective terms. By uncoupling the velocity lattice from the spatial grid, one can employ higher order (non-space-filling) isotropic lattices-lattices which greatly enhance the stable parameter regions, particularly in thermal problems. In particular, the superiority of the octagonal lattice over previous models used in 2D (hexagonal or square) and 3D (projected face-centered hypercube) is shown

Pavlo, Pavol

Vahala, George

Vahala, Linda L.

English

Old Dominion University

Old Dominion UniversityODU Digital CommonsElectrical & Computer Engineering FacultyPublications Electrical & Computer Engineering1998Higher Order Isotropic Velocity Grids in LatticeMethodsPavol PavloGeorge VahalaLinda L. VahalaOld Dominion University, lvahala@odu.eduFollow this and additional works at: https://digitalcommons.odu.edu/ece_fac_pubsPart of the Physics CommonsThis Article is brought to you for free and open access by the Electrical & Computer Engineering at ODU Digital Commons. It has been accepted forinclusion in Electrical & Computer Engineering Faculty Publications by an authorized administrator of ODU Digital Commons. For more information,please contact digitalcommons@odu.edu.Repository CitationPavlo, Pavol; Vahala, George; and Vahala, Linda L., "Higher Order Isotropic Velocity Grids in Lattice Methods" (1998). Electrical &Computer Engineering Faculty Publications. 55.https://digitalcommons.odu.edu/ece_fac_pubs/55Original Publication CitationPavlo, P., Vahala, G., & Vahala, L. (1998). Higher order isotropic velocity grids in lattice methods. Physical Review Letters, 80(18),3960-3963. doi: 10.1103/PhysRevLett.80.3960VOLUME 80, NUMBER 18 P H Y S I C A L R E V I E W L E T T E R S 4 MAY 1998Higher Order Isotropic Velocity Grids in Lattice MethodsPavol PavloInstitute of Plasma Physics, Czech Academy of Sciences, Praha 8, Czech RepublicGeorge Vahala*Department of Physics, College of William & Mary, Williamsburg, Virginia 23187Linda VahalaDepartment of Electrical & Computer Engineering, Old Dominion University, Norfolk, Virginia 23529(Received 8 December 1997)Kinetic lattice methods are a very attractive representation of nonlinear macroscopic systemsbecause of their inherent parallelizability on multiple processors and their avoidance of the nonlinearconvective terms. By uncoupling the velocity lattice from the spatial grid, one can employ higherorder (non-space-filling) isotropic lattices—lattices which greatly enhance the stable parameter regions,particularly in thermal problems. In particular, the superiority of the octagonal lattice over previousmodels used in 2D (hexagonal or square) and 3D (projected face-centered hypercube) is shown.[S0031-9007(98)05978-X]PACS numbers: 47.11.+ j, 05.50.+qLattice methods can play a major role in solving sys-tems, particularly those whose evolution involves nonlin-ear advection terms. In particular, the kinetic lattice gas [1]and lattice Boltzmann [2] methods have been implementedsuccessfully in the solution of Navier-Stokes flows, mag-netohydrodynamics, flow through porous media, turbulentflows, phase transitions, multiphase flows, and chemicallyreacting flows as well as to liquid crystals and semiclassi-cal transport in semiconductors. In many instances, latticemethods can be viewed as a maximally discretized molecu-lar dynamics—an optimized mesoscopic description be-tween the full molecular dynamics and the macroscopicconservation laws. However, nearby all applications havebeen restricted to macroscopic systems that are closed be-low the second moment (energy) level. Only a few at-tempts [3–6] have been successful in simulating thermallattice Boltzmann systems (TLBE) because inherent nu-merical instabilities severely limit the accessible parameterregimes [5]. It is the purpose of this Letter to suggest anew TLBE model that displays numerical stability over awide range of parameters—not only in 2D, but also in 3D.In their simplest form, these kinetic lattice models areideally implemented on multiparallel processors becauseall operations are purely local; the Lagrangian schemesare explicit, with complex boundaries being readilyhandled. One of the major requirements is that the under-lying discrete symmetry of the chosen lattice not mar theintended physics emulating from the continuously sym-metric macroscopic equations. Initially, the spatial and(phase space) velocity lattices were intrinsically coupledand, because of the need to be space filling, this restrictedthe lattice geometry in 2D to be either a composite squaregrid or a hexagonal grid (with sufficient number of propa-gation speeds) or the projection of a 4D face-centeredhypercube (FCHC) into 3D [7]. Recently, nonuniformspatial grids have been introduced to more accuratelyhandle wall-bounded flows [8]. This, as well as therecent use of finite difference schemes [9], allows for theecoupling of the spatial grid from the velocity lattice.However, it does not seem to have been appreciated thathis decoupling now frees us to employ higher ordersymmetric velocity lattices which result in impressivegains in stability. Indeed, here we shall examine someof the consequences of employing an octagonal velocitylattice in 2D and its generalizations to 3D.We shall utilize a Lagrangian free-streaming algorithmbut, of course, our ideas are immediately applicable to anyother scheme (finite difference, finite volume, etc.). On aparticularly chosen discrete lattice, the kinetic equation (inlattice units) takes the formNpisx 1 cpi, t 1 1d 2 Npisx, td ­ Vpisx, td,i ­ 1, . . . , bp , p ­ 1, . . . , S ,(1)wherecpi is the lattice velocity vector with free streamingin direction i at speedjcpij. bp is the number of linksat speed indexp, S is the total number of speeds inthe model, andNB is the total number of base vectors(or “bits”): NB ­Pbp . For lattice gas [1],Npi is aninteger representing the particle occupation number onsite x at time t for momentum state (pi) while, for latticeBoltzmann [2],Npi is a (continuous) particle distributionfunction. The collision operatorVpi is constrained tosatisfy (at each spatial node) particle, momentum, andenergy conservationXpiVpisx, td ­ 0 ­XpiVpisx, tdcpi­12XpiVpisx, tdc2pi . (2)3960 0031-9007y98y80(18)y3960(4)$15.00 © 1998 The American Physical SocietyVOLUME 80, NUMBER 18 P H Y S I C A L R E V I E W L E T T E R S 4 MAY 1998For variable Prandtl number flows, one must employ anextended [1,6] Bhatnagar-Gross-Krook (BGK) collisionoperatorVpi ­ 21tXjAijfNpj 2 Neqpj g,Aij ­ dij 1utbpc2pcpi ? cpj ,(3)with Aij being a circulant matrix and the relaxation ratest and u controlling the viscositym(t) and conductivitykst, ud transport coefficients [1,2].m ! 0 ast ! 0.51.The form of Neq is critical. A local H theorem [10]can be shown for exponential (Maxwellian-like)Neq,but, to enforce the nodal conservation quantities (2) andthe moment definitions of mass, momentum, and internalenergy, XpiNeqpi ­ r,XpiNeqpi cpi ­ ru,XpiNeqpi cpiacpib ­2Dredab 1 ruaub ,(4)exactly, one is forced into a polynomial representation:Neqpi ­ rfApsed 1 Bpsedcpi ? u 1 Cpsed scpi ? ud21 Dpsedu2 1 1Epsedcpi ? uu2 1 Fpsed scpi ? ud31 Gpsed scpi ? ud2u2 1 Hpsedu4g . (5)The coefficientsAp , . . . , Hp are constrained by the levelof isotropy in the chosen velocity lattice as well assatisfying (4) and those constraints needed to recover themacroscopic conservation equations in the long-time long-wavelength Chapman-Enskog limit:XpiNeqpi cpiacpibcpig ­2Dresuadbg 1 ubdga 1 ugdabd1 ruaubug , (6)XpiNeqpi cpiacpibc2pi ­4sD 1 2dD2re2dab 12Dreu2dab12sD 1 4dDreuaub 1 ruaubu2(7)(here D is the dimensionality of the spatial grid andsummation over repeated subscripts is understood).It is clear from (4)–(7), that thenth-lattice velocitymomentsTsnda,...,z ;XpTsndp,a,...,z ;XpXicpia · · · cpiz (8)play a critical role. For arbitrary regular lattices, one findsTs4dp,abgd ­ cpYabgd 1 fpsdabdgd 1 c.p.d ,Ts6dp,abgd´z ­ CpYabgd´z 1 LpsdabYgd´z 1 c.p.d1 FpsdabTs4dgd´z 1 c.p.d , (9)whereY... is the higher rank Kronecker tensor (Y... ­ 1if all indices are equal, 0 otherwise) and is anisotropic.cp , . . . , Fp are coefficients determined explicitly oncethe lattice is specified and “c.p.” denotes cyclic per-mutation of indices. For example, consider the 2DTLBE 13-bit square latticehp ­ 1, i ­ 1, . . . , 4: c1i ­s61, 0d, s0, 61dj; hp ­ 2, i ­ 1, . . . , 4: c2i ­ s61, 61dj;hp ­ 3, i ­ 1, . . . , 4: c3i ­ s62, 0d, s0, 62dj and the restparticle hp ­ 0, i ­ 1: c0i ­ s0, 0dj. This lattice doesnot yield an isotropic fourth rank lattice tensorT s4dp atany individual speedp, (9). To recover overallT s4disotropy, (8), the minimal isotropy needed to recover thecorrect conservation equations (at least to leading order),further constraints (and speeds) must be enforced on thecoefficientsAp . . . of the relaxation distribution function,(5). These extra constraints have such strong adverseeffects on the numerical stability (see Fig. 1) that this 2DTLBE 13-bit free-streaming square model is practicallyuseless for our applications. It has been noted [11] that a2D 17-bit square lattice can recover overallT s6d isotropy,but these models are even more numerically unstabledue to the extra constraints placed on the expansioncoefficients in (5).The 2D 13-bit hexagonal lattice [3] does yield isotropicT s4dp but not Ts6d tensors. As a result, the hexagonallattice exhibits better stability properties than the squarelattice [12] (see Fig. 1), but it does not exhibit sufficientisotropy to eliminate the spurious cubic deviations inthe macroscopic conservation equations [11]. Moreover[11], there does not seem to be a higher-bit hexagonalFIG. 1. Linear stability windowssrl ­ 1d for the octagonalgrid, 0.28 , e , 0.555, and the hexagonal grid,0.32 , e ,0.35 and0.49 , e , 0.505 for t ­ 0.504, and Prandtl numberPr ­ myk ­ 0.5 (i.e., u ­ 0). For this value oft, there isno stable region for the 13-bit square grid, with minrl ­ 1.18(t ­ 0.5 corresponds to infinite Reynolds number).3961VOLUME 80, NUMBER 18 P H Y S I C A L R E V I E W L E T T E R S 4 MAY 1998lattice that will enforceT s6d isotropy. The reason so muchemphasis has been placed on these 2D lattices is that theyare space filling.Since here we uncouple the velocity lattice fromthe spatial grid, the velocity lattice need no longer bespace filling. Consider the 2D 17-bit octagonal model,with lattice vectorsfcpi ­ pscospiy8, sinpiy8d, i ­0, . . . , 7, p ­ 1, 2g and rest particle. This octagonal lat-tice exhibits individual speed isotropy up to the sixth-ranklattice tensorsT s6dp [see Table I below and Eq. (9)] so thatone can recover the thermal macroscopic equations ex-actly, without spurious cubic nonlinearities [11].To proceed, we must now specify a spatial grid, and,for convenience, we consider the standard 17-bit squaregrid. If we free stream one time step using the octagonalvelocitiescpi , only those directions withi even will fall ona spatial square node. In the diagonal directions (i odd),the spatial nodal radius­ 21y2p while the octagonal free-streaming distance­ p. Thus we must resort to someinterpolation schemes to latch onto the chosen spatial gridnodes. The numerical stability [12] of this 2D octagonalTLBE model is shown in Fig. 1 for relaxation timest ­ 0.504 andu ­ 0 (i.e., for the standard BGK collisionterm. Previous numerical tests on the 2D hexagonallattice have shown that the stability is independent ofthe Prandtl numbers Pr­ myk for Pr . 0.03). Note thatthe 2D 13-bit square lattice has no stability window inthe internal energye for any t , 0.58. The hexagonallattice has limited disjoint stability windows0.32 , e ,0.35 and 0.49 , e , 0.505 for t ­ 0.504, while the2D octagonal mode has an extensive continuous stabilitywindow 0.28 , e , 0.555.Here we have employed a second order interpolationscheme for the diagonal directions (i odd). It is importantto ascertain the effect of interpolation on the transportcoefficients—especially since interpolation will alwaysbe needed in Lagrangian schemes on uncoupled lattices.Performing the standard [3] isothermal plane Poiseuilleflow determination of the viscosity and the Fourierheat flux law for the conductivity, we find excellentagreement between theoretical Chapman-Enskog transportand that determined from our interpolated octagonalTLBE simulation, Figs. 2 and 3. (A standard bounce-back scheme is employed at the boundaries.)The temperature profile in 2D Rayleigh-Benard con-vection is presented in Fig. 4 at timest ­ 5 K and t ­TABLE I. Individual speed isotropy of the octagonal lattice.Octagonal grid(17-bit) cp fp Cp Lp Fpp ­ 1Speed-1 0 1 0 0 16p ­ 2Speed-2 0 16 0 0 332FIG. 2. Comparison of the theoretical Chapman-Enskog vis-cosity nTH (solid lines) to the 17-bit octagonal TLBE simu-lation viscosity nSIM (filled squares fore ­ 0.5, and opensquares fore ­ 0.3) as a function oft at Pr­ 1 for isother-mal plane Poiseuille flow. We have verified thatnSIM is in-dependent the Prandtl number and find excellent agreementwith nTH: e.g., fore ­ 0.5 and Pr­ 1, nTH ­ 20.25 1 0.5t,while nSIM ­ 20.249 1 0.499t.200 K lattice units. The lower plates y ­ 0d is at inter-nal energye ­ 0.5, while the upper plates y ­ 129d isat e ­ 0.3. The gravitational forcerG is the2y direc-tion and is readily incorporated into the collision operator,Eq. (3), by the addition [11] of a term proportional tocpi .The Rayleigh number Ra, a ratio of the buoyancy to diffu-sion, for this run was6 3 106. A weak period-4 perturba-tion in thex direction was initially imposed on a straightline conduction profile. After 5 K iterations, Fig. 4(a),four strong convection cells are evident. At this highRayleigh number, the flow is turbulent. Byt . 30 K,turbulent effects are evident [c.f., Fig. 4(b)].Finally, we have performed some preliminary stabilitytests on our extended 3D “octagonal” 53-bit TLBE modeland compared them to that for the standard 4D FCHC-projected 34-bit model in 3D. Fort ­ 0.52, the stan-dard FCHC-projected model has no stability window inFIG. 3. The dependence of the conductivity on the Prandtlnumber att ­ 0.51. Again, excellent agreement is found be-tween the theoretical Chapman-Enskog conductivitykTH (solidcurve) and that from the 17-bit TLBE octagonal simulation(solid diamonds).3962VOLUME 80, NUMBER 18 P H Y S I C A L R E V I E W L E T T E R S 4 MAY 1998FIG. 4. Temperature profiles for 2D Rayleigh-Benard convec-tion. The lower plates y ­ 0d is at temperature ­ 0.5 whilethe upper plates y ­ 129d is ate ­ 0.3. The Rayleigh numberRa ­ 6.0 3 106. After a short time, the convection patterns[(a) t ­ 5 K] become turbulent [(b)t ­ 200 K].energye, while our 53-bit model has an extensive stableregion 0.28 , e , 0.72. Thus, not only do we have amuch greater stability window ine, but our model canenforceT s6d isotropic—while the FCHC model has onlyT s4d isotropy. More details of our 3D model will bepresented elsewhere.This work was supported by a joint U.S.-Czech grantand DOE.*Also at Institute for Computer Applications in Scienceand Engineering, NASA Langley Research Center,Hampton, Virginia 23665[1] H. Chen, C. Teixeira, and K. Molvig, Int. J. Mod. Phys. C8, 675 (1997).[2] S. Chen, S. P. Dawson, G. D. Doolen, D. R. Janecky, andA. Lawniczak, Comput. Chem. Eng.19, 617 (1995), andreferences therein.[3] F. J. Alexander, S. Chen, and J. D. Sterling, Phys. Rev. E47, R2249 (1993).[4] Y. Chen, H. Ohashi, and M. Akiyama, Phys. Rev E50,2276 (1994).[5] G. R. McNamara, A. L. Garcia, and B. J. Alder, J. Stat.Phys.81, 395 (1995).[6] G. Vahala, P. Pavlo, L. Vahala, and M. Soe, Czech.J. Phys.46, 1063 (1996).[7] U. Frisch, D. d’Humieres, B. Hasslacher, P. Lallemand,Y. Pomeau, and J. P. Rivet, Complex Systems1, 649(1987); S. Wolfram, J. Stat. Phys.45, 471 (1986).[8] X. He, L. S. Luo, and M. Dembo, J. Comput. Phys.129,357 (1996).[9] N. Cao, S. Chen, S. Jin, and D. Martinez, Phys. Rev. E55, R21 (1997).[10] M. Henon, Complex Systems1, 763 (1987); H. Chen,J. Stat. Phys.81, 347 (1995).[11] Y. Chen, Ph.D. thesis, University of Tokyo, 1994); Y. H.Qian and S. A. Orszag, Europhys. Lett.21, 255 (1993).[12] P. Pavlo, G. Vahala, L. Vahala, and M. Soe, J. Comput.Phys.139, 79 (1998).3963

Higher Order Isotropic Velocity Grids in Lattice Methods

https://digitalcommons.odu.edu/cgi/viewcontent.cgi?article=1048&context=ece_fac_pubs

Higher Order Isotropic Velocity Grids in Lattice Methods

Abstract

Similar works

Full text

Available Versions

Old Dominion University