Stabilization of the lattice Boltzmann method by the H theorem: A numerical test Santosh Ansumali and Iliya V. Karlin* ETH-Zu ¨rich, Department of Materials, Institute of Polymers, ETH-Zentrum, Sonneggstrasse 3, ML J 27, CH-8092 Zu ¨rich, Switzerland Received 7 June 2000 For a one-dimensional benchmark shock tube problem, we implement the lattice Boltzmann method based on the H theorem I. Karlin, A. Ferrante, and H. C. O ¨ ttinger, Europhys. Lett. 47, 182 1999. Results of simulation demonstrate significant improvement of stability, as compared to realizations without explicit en- tropic estimations. PACS numbers: 47.11.+j, 05.20.Dd I. INTRODUCTION Since the invention of the lattice-gas model 1, lattice- based methods for simulations of complex hydrodynamic phenomena received much attention over the past decade. In these methods, hydrodynamic equations are not addressed by a direct discretization procedure, rather, a simple pseudopar- ticle kinetics is introduced in such a way that the hydrody- namic equations are obtained on the large space and time scale. Particularly promising is the well-known lattice Bolt- zmann method LBM2. It is based on the fully discrete velocity-space-time kinetic equation of the form, Nx +c, t +1 -Nx , t =Nx , t  . 1 Here N( x , t ) is the b-component vector of populations N i of the pseudoparticles with velocities c i , at the sites x of a lattice at discrete time t. The system of discrete velocities at any site is formed by the outgoing links of the lattice, and it also may include the zero vector. One of the most important problems related to the LBM, recognized by many authors, is the problem of numerical stability. For the LBM related to incompressible flow simu- lations, numerical instabilities preclude so far a study of high Reynolds number flow situations. Instabilities become even more annoying for compressible flows 3. It has been discussed for some time in the literature that stability of the LBM could be improved if the method could be equipped with an analog of the Boltzmann H theorem. Recently, theoretical progress in this direction has been achieved 4–8. In particular, for the isothermal LBM, the hydrodynamic fields are the density, =( 1, N), and the av- erage momentum, u =( c , N), where ( , ) denotes the standard scalar product in the b-dimensional space of popu- lation vectors. In this case, one can construct entropy func- tions in such a way that its local equilibrium implies the crucial relation for the stress tensor, c c , N eq =c s 2  +u u , up to the admissible degree of accuracy of the LBM 6. Furthermore, it has been suggested how to construct the col- lision integral based solely on the knowledge of the en- tropy function, and how to stabilize the updates on the basis of the discrete-time H theorem 6,8. In the sequel, we term the LBM based on the H theorem the entropic lattice Boltz- mann method ELBM. It is the goal of this paper to test the aforementioned the- oretical developments for a shock tube problem. This one- dimensional benchmark problem has been suggested some time ago 9for testing various ideas in the LBM. Though this model is based on a very simple three-velocity lattice, it provides a stringent test of stability. Implementation of the ELBM demonstrated a large improvement of stability in this benchmark problem. In fact, we were able to reach the values of the kinematic viscosity as low as 10 -12 without any sign of numerical instabilities. The most important part of the realization is a robust root-finding procedure which imple- ments the H theorem. II. CONSTRUCTION OF THE ELBM An advantage of the LBM in comparison to the lattice-gas method is that the Galilean invariance of the Navier-Stokes equation is easier to control in the former than in the latter. In order that this advantage should not get lost in the ELBM, entropy functions should be found for each lattice separately. We here consider the one-dimensional lattice with spacing c, and the population vector at each site x has three compo- nents, N=( N + , N 0 , N - ) ² , corresponding to velocities c + =c , c 0 =0, and c - =-c , respectively. For this model, the entropy function has been found in 6 H =N 0 lnN 0 /4+N - lnN - +N + lnN + . 2 Realizations of the ELBM based on the entropy function 2 result in the one-dimensional Navier-Stokes equation with the sound speed c s =1/3c , within the accuracy of the order ( u / c s ) 4 . Construction of the ELBM involves the following two steps these steps are independent of the choice of the entropy. First, we specify bare collision integral in such a way as to satisfy the admissibility condition, ( , 1) =( , c ) =0, ( , H ) 0, and ( N eq ) =0. Here H is the gradient of H in the space of population vectors. The choice of the bare collision integral is not unique. We here consider three cases, =g + -g - exp H , g -  -exp H , g +  , 3a *Corresponding author. FAX: +41 01 632 10 76. Email address: ikarlin@ifp.mat.ethz.ch PHYSICAL REVIEW E DECEMBER 2000 VOLUME 62, NUMBER 6 PRE 62 1063-651X/2000/626/79995/$15.00 7999 ©2000 The American Physical Society