Lattice Boltzmann model with self-tuning equation of state for coupled thermo-hydrodynamic flows

# Lattice Boltzmann model with self-tuning equation of state for coupled thermo-hydrodynamic flows

Rongzong Huang Huiying Wu Nikolaus A. Adams School of Mechanical Engineering, Shanghai Jiao Tong University, 200240 Shanghai, China Institute of Aerodynamics and Fluid Mechanics, Technical University of Munich, 85748 Garching, Germany
###### Abstract

A novel lattice Boltzmann (LB) model with self-tuning equation of state (EOS) is developed in this work for simulating coupled thermo-hydrodynamic flows. The velocity field is solved by the recently developed multiple-relaxation-time (MRT) LB equation for density distribution function (DF), by which a self-tuning EOS can be recovered. As to the temperature field, a novel MRT LB equation for total energy DF is directly developed at the discrete level. By introducing a density-DF-related term into this LB equation and devising the equilibrium moment function for total energy DF, the viscous dissipation and compression work are consistently considered, and by modifying the collision matrix, the targeted energy conservation equation is recovered without deviation term. The full coupling of thermo-hydrodynamic effects is achieved via the self-tuning EOS and the viscous dissipation and compression work. The present LB model is developed on the basis of the standard lattice, and various EOSs can be adopted in real applications. Moreover, both the Prandtl number and specific heat ratio can be arbitrarily adjusted. Furthermore, boundary condition treatment is also proposed on the basis of the judicious decomposition of DF into its equilibrium, force (source), and nonequilibrium parts. The local conservation of mass, momentum, and energy can be strictly satisfied at the boundary node. Numerical simulations of thermal Poiseuille and Couette flows are carried out with three different EOSs, and the numerical results are in good agreement with the analytical solutions. Then, natural convection in a square cavity with a large temperature difference is simulated for the Rayleigh number from up to . Good agreement between the present and previous numerical results is observed, which further validates the present LB model for coupled thermo-hydrodynamic flows.

###### keywords:
lattice Boltzmann model, coupled thermo-hydrodynamic flows, self-tuning equation of state, viscous dissipation and compression work, boundary condition treatment, standard lattice

## 1 Introduction

The lattice Boltzmann (LB) method has developed into an attractive numerical method over the past three decades for simulating complex fluid flows Aidun2010 ; Gross2010 ; Huang2016.amr and solving various partial differential equations Ginzburg2005 ; Dellar2011 ; Chai2018 . Historically, the LB method originates from the lattice gas automata (LGA) to eliminate the statistical noise McNamara1988 , and thus it inherits some distinguishing features from LGA, such as the simple algorithm (local collision and linear streaming) and the easy incorporation of microscopic interactions Shan1993 ; Ladd1994 . Afterward, it is found that the classical LB model for hydrodynamic flows can be derived from the Boltzmann-BGK equation via systematic discretization He1997a ; He1997b , and then various LB models for multiphase flows He1998.nonideal ; Luo1998 and thermo-hydrodynamic flows (i.e., thermal fluid flows) He1998.thermal ; Guo2007 have been established from the kinetic models in an a priori manner.

Since most hydrodynamic flows involve some forms of thermal effects, thermo-hydrodynamic flows are extensively encountered in nature and engineering, and the LB method for simulating thermo-hydrodynamic flows has attracted continuous attention since the early 1990s He1998.thermal ; Guo2007 ; Shan1997 ; Alexander1993 ; McNamara1997 ; Zheng2008 ; Lallemand2003 ; Mezrhab2004 ; Mezrhab2010 ; Huang2014.ibm ; Contrino2014 . However, it remains open-ended though the LB method has achieved great success in simulating isothermal fluid flows Sbragaglia2009 ; Succi2015 . Generally, the existing LB models for thermo-hydrodynamic flows can be categorized into three major groups: the multispeed approach McNamara1997 ; Zheng2008 , the double-distribution-function (DDF) approach He1998.thermal ; Guo2007 ; Shan1997 , and the hybrid approach Lallemand2003 ; Mezrhab2004 . The multispeed approach uses a single distribution function (DF) to describe the mass, momentum, and energy conservation laws, and thus it requires more discrete velocities than the standard lattice (i.e., it requires the multispeed lattice). By definition, the DDF approach consists of double DFs, with one DF for the mass and momentum conservation laws and the other DF for the energy conservation law. In the hybrid approach, the mass and momentum conservation laws are described by one DF, while the energy conservation law is described by a macroscopic governing equation that is solved via the conventional computational fluid dynamics methods. Severe numerical instability Lallemand2003 and complexity of boundary condition treatment Lee2018 are usually encountered in the multispeed approach. As to the hybrid approach, it acts as a compromised solution that deviates from the mesoscopic LB method Lallemand2003 , and the viscous dissipation is usually ignored in this approach Feng2018 ; Safari2018 . On the contrary, the DDF approach, free of the above drawbacks, is most widely studied and adopted in real applications.

Most of the existing DDF LB models for thermo-hydrodynamic flows are inherently a decoupling model, which means that the recovered equation of state (EOS) is a decoupling EOS ( is the gas constant and is the reference temperature), where the pressure is not directly related to the temperature. Consequently, these LB models are restricted to the thermo-hydrodynamic flows under the Boussinesq approximation (i.e., the decoupling thermo-hydrodynamic flows). Based on the DDF kinetic model constructed by Guo et al. Guo2007 , and by applying the discretization of velocity space presented by Shan et al. Shan2006 that can lead to the temperature-independent discrete velocities, Hung and Yang Hung2011 proposed a DDF LB model aimed at recovering the ideal-gas EOS . However, the deviation in the third-order moment of the equilibrium distribution function (EDF) for density due to the constraint of standard lattice, as previously identified by Prasianakis and Karlin Prasianakis2007 , is not considered in Hung and Yang’s model, and meanwhile, an error also exists in their derived EDF for total energy. In 2012, by introducing the correction term for the third-order moment of the EDF for density and deriving the correct EDF for total energy, Li et al. Li2012 developed a DDF LB model for simulating coupled thermo-hydrodynamic flows. The ideal-gas EOS can be recovered by Li et al.’s model, and the simulation of natural convection with a large temperature difference is reported Li2012 . Following the similar way, Feng et al. Feng2015 proposed three-dimensional DDF LB models. A correction term for the second-order moment of the EDF for total energy is further introduced by Feng et al. Feng2015 to enhance the numerical stability of the LB equation for total energy DF. Recently, the cascaded collision scheme is employed in the LB equation for density DF to enhance the numerical stability by Fei and Luo Fei2018 , while the single-relaxation-time (SRT) collision scheme is still used in the LB equation for total energy DF.

It is worth pointing out that the ideal-gas EOS is recovered by the above DDF LB models Hung2011 ; Li2012 ; Feng2015 ; Fei2018 , which indicates that these models are only applicable to the coupled thermo-hydrodynamic flows of ideal gases. Moreover, in these models, the LB equation for total energy DF is complicated due to the consideration of the viscous dissipation and compression work, and thus it is difficult to employ the multiple-relaxation-time (MRT) or cascaded collision schemes in this LB equation to enhance the numerical stability although the MRT and cascaded collision schemes have been employed in the LB equation for density DF Li2012 ; Fei2018 . Most recently, we developed an LB model with self-tuning EOS for multiphase flows Huang2018.eos . Since the recovered EOS can be self-tuned via a built-in variable, this model serves as a good and distinct starting point for developing a novel LB model for coupled thermo-hydrodynamic flows, which is the main objective of the present work. To be specific, a novel MRT LB equation for solving the energy conservation equation, with considering the viscous dissipation and compression work, is developed. Furthermore, boundary condition treatment for simulating coupled thermo-hydrodynamic flows is also proposed on the basis of the judicious decomposition of DF into three parts rather than two. The remainder of the present paper is organized as follows. In Section 2, a novel LB model for coupled thermo-hydrodynamic flows is developed. In Section 3, boundary condition treatment is proposed. Numerical validations of the present LB model are carried out in Section 4, and a brief conclusion is drawn in Section 5.

## 2 Lattice Boltzmann model

The present LB model for coupled thermo-hydrodynamic flows is developed on the basis of the recent LB model with self-tuning EOS for multiphase flows. Double DFs are involved: one is the density DF used to solve the velocity field (i.e., the mass-momentum conservation equations), and the other is the total energy DF used to solve the temperature field (i.e., the energy conservation equation). The full coupling of thermo-hydrodynamic effects is achieved via the self-tuning EOS recovered by the LB equation for density DF and the viscous dissipation and compression work considered in the LB equation for total energy DF. Both the LB equations for density and total energy DFs are based on the standard lattice. For the sake of simplicity and clarity, the two-dimensional model will be developed here, and its extension to three-dimensional model is straightforward. The standard two-dimensional nine-velocity (D2Q9) lattice is given as Qian1992

 ei=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩c(0,0)T,i=0,c(cos[(i−1)π/2],sin[(i−1)π/2])T,i=1,2,3,4,√2c(cos[(2i−1)π/4],sin[(2i−1)π/4])T,i=5,6,7,8, (1)

where the lattice speed with and being the lattice spacing and time step, respectively.

### 2.1 LB equation for density DF

The recently developed LB equation for density DF that can recover a self-tuning EOS is briefly introduced here for self-completeness. The MRT LB equation for density DF can be expressed as Huang2018.eos

 fi(x+eiδt,t+δt)=¯fi(x,t), (2a) ¯m(x,t)=m+δtFm−S(m−meq+δt2Fm)−R(I−S2)(m−meq+δt2Fm)−δxT⋅∇ρ−δxc2X⋅∇p\tiny LBE, (2b)

where Eq. (2a) is the streaming process executed in velocity space and Eq. (2b) is the collision process executed in moment space at position and time . The moment of density DF in Eq. (2b) is given as . Here, is the dimensionless transformation matrix Lallemand2000

 M=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣111111111−4−1−1−1−122224−2−2−2−21111010−101−1−110−20201−1−110010−111−1−100−20211−1−101−11−10000000001−11−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (3)

and denotes the vector . The post-collision density DF in Eq. (2a) is obtained via the inverse transformation , and the post-collision moment is computed by Eq. (2b). The last three terms on the right-hand side (RHS) of Eq. (2b) are the correction terms aimed at eliminating the additional cubic terms of velocity in the recovered momentum conservation equation Huang2018.cubic , where denotes the recovered EOS by the LB equation. The macroscopic density and velocity are defined as

 ρ=∑ifi,ρu=∑ieifi+δt2F, (4)

where is the force term. In the recent LB model for multiphase flows Huang2018.eos , is the total force due to the long-range molecular interaction, while in the present LB model for coupled thermo-hydrodynamic flows, is simply an external force, such as the gravity force.

In Eq. (2b), the equilibrium moment function for density DF is given as Huang2018.eos

 meq=[ρ,2α1ρ+2β1η+3ρ|^u|2,α2ρ+β2η−3ρ|^u|2+9ρ^u2x^u2y,ρ^ux,−ρ^ux+3ρ^ux^u2y,ρ^uy,−ρ^uy+3ρ^uy^u2x,ρ(^u2x−^u2y),ρ^ux^uy]T, (5)

where and is the built-in variable aimed at achieving a self-tuning EOS. The coefficients and are set to and , respectively, while the coefficients and are determined by Eq. (8). The discrete force term in moment space is given as

 Fm=[0,6^F⋅^u,−6^F⋅^u+9[^F^u^u^u]xxyy,^Fx,−^Fx+3[^F^u^u]xyy,^Fy,−^Fy+3[^F^u^u]xxy,2(^Fx^ux−^Fy^uy),^Fx^uy+^Fy^ux]T, (6)

where , and the square bracket and its subscript denote permutation and tensor index, respectively. For example, . To correctly recover the Newtonian viscous stress tensor, the collision matrix in moment space is modified as follows Huang2018.eos

 S=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣s0000000000seksεωe0h^uxsqωe0h^uysqωe0000sε000000000sj000000000sq000000000sj000000000sq0000002b^uxsqωp0−2b^uysqωpsp00000b^uysqωp0b^uxsqωp0sp⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (7)

where , and , , and are the coefficients. Through the Chapman-Enskog analysis, the coefficients in and should satisfy the following relations

 α2=−2α1+ϖ+11−ϖ,β2=−2β11−ϖ,k=1−ϖ,h=6ϖ(1−ϖ)1−3ϖ,b=1−ϖ1−3ϖ, (8)

where is related to the bulk viscosity.

In Eq. (2b), the last three terms, together with the high-order terms of velocity in and , are introduced to eliminate the additional cubic terms of velocity Huang2018.cubic , which are not considered in the previous DDF LB models for coupled thermo-hydrodynamic flows. The correction matrix is a matrix and it is set as Huang2018.eos

 R=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0000000000R1100000R17R180000000000000000000000000000000000000000000000R7100000R77R780R8100000R87R88⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (9a) where the nonzero elements can be determined via the Chapman-Enskog analysis as follows R11=−(9−15k−2h)se4ϖ(^u2x+^u2y),R17=−3(9−3k+2h)sp(2−se)4(2−sp)(^u2x−^u2y),R18=12(3k+h)sp(2−se)2−sp^ux^uy,R71=−(3−4b)se(2−sp)4ϖ(2−se)(^u2x−^u2y),R77=−3(3+4b)sp4(^u2x+^u2y),R78=0,R81=bse(2−sp)ϖ(2−se)^ux^uy,R87=0,R88=6bsp(^u2x+^u2y). (9b)

The correction matrix is set as Huang2018.eos

 T=(0,T1,0,0,0,0,0,T7,T8)T, (10a) whose element is a vector implying that the dimensions of T are 9×2. The nonzero elements in T can also be determined via the Chapman-Enskog analysis as follows T1=3(2−se)2[(1−k)^u3x−(2k+h)^ux^u2y(1−k)^u3y−(2k+h)^u2x^uy],T7=2−sp2[^u3x+2b^ux^u2y−^u3y−2b^u2x^uy],T8=−b(2−sp)2[^u3y+2^u2x^uy^u3x+2^ux^u2y]. (10b)

Similarly to , the correction matrix is a matrix and it is set as Huang2018.eos

 X=(0,X1,0,0,0,0,0,X7,X8)T, (11a) where the nonzero elements are given as (11b)

Here, we would like to point out that for the coupled thermo-hydrodynamic flows under the low Mach number condition, these correction terms for the additional cubic terms of velocity can be simply ignored. However, they are kept in the present work for the sake of theoretical completeness and computational accuracy.

Through the Chapman-Enskog analysis, the following mass-momentum conservation equations can be recovered Huang2018.eos

 {∂tρ+∇⋅(ρu)=0,∂t(ρu)+∇⋅(ρuu)=−∇p% \tiny LBE+F+∇⋅Π, (12)

where and are the recovered EOS and viscous stress tensor

 p\tiny LBE=c2s[(2+α1)ρ+β1η],Π=ρν[∇u+u∇−(∇⋅u)I]+ρς(∇⋅u)I, (13)

where the lattice sound speed , the kinematic viscosity , and the bulk viscosity . As seen in Eq. (13), the recovered EOS can be arbitrarily tuned via the built-in variable .

### 2.2 LB equation for total energy DF

Since the EOS recovered by the above LB equation for solving the velocity field can be self-tuned, we are now well equipped to simulate coupled thermo-hydrodynamic flows. The remaining task is to develop an LB equation for solving the temperature field, in which the viscous dissipation and compression work are consistently considered.

#### 2.2.1 Energy conservation equation

The collision term of an LB equation conserves macroscopic quantity, and the recovered macroscopic conservation equation for this quantity usually has a conservative form (see Eq. (13) as an example). On the basis of this principle, the total energy conservation equation, in which the viscous dissipation and compression work are expressed as , is a better and more natural starting point for directly developing an LB equation at the discrete level than the internal energy conservation equation, in which the viscous dissipation and compression work are expressed as . Here, is the pressure determined by the adopted EOS. To facilitate the development of an LB equation, the total energy conservation equation is reformulated as

 ∂t(ρE)+∇⋅(ρHu)=∇⋅(λ∇T+u⋅Π)+u⋅F+q, (14)

where is the total energy, is the total enthalpy, is the temperature that can be determined by the internal energy () and density , is the heat conductivity, and is the source term. In Eq. (14), the viscous dissipation combines with the conduction term to constitute the term , and the compression work combines with the convection term to constitute the term . Meanwhile, we can also combine the work done by force and the source term to constitute an equivalent source term . Thus, Eq. (14) can be viewed as a general convection-diffusion equation with source term. Here, we would like to point out that the above reformulation is consistent with the Chapman-Enskog analysis, which means that the two terms combined together are of the same order.

#### 2.2.2 Viscous stress tensor

To consider the viscous dissipation in the LB equation for total energy DF, we first recall the recovery of viscous stress tensor by the above LB equation for density DF. On the basis of the Chapman-Enskog analysis, the viscous stress tensor is of order and can be expressed as Huang2018.eos

 (15)

where is the small expansion parameter in the Chapman-Enskog analysis and is

 ~G(1)=m(1)+¯m(1)2, (16)

where and are the terms of and in their Chapman-Enskog expansions and , respectively. Here, it is worth pointing out that the post-collision moment is kept in Eq. (16) rather than being substituted by Eq. (2b). As a consequence, the post-collision moment , which is computed in the collision process of density DF, can be directly utilized to consider the viscous dissipation in the LB equation for total energy DF (see C). Moreover, from the Chapman-Enskog analysis of the LB equation for density DF, we can easily know that the terms of and satisfy

 m(0)=¯m(0)=meq. (17)

#### 2.2.3 LB equation

For the energy conservation equation given by Eq. (14), the total energy DF is introduced here, and the MRT LB equation for is devised as

 gi(x+eiδt,t+δt)=¯gi(x,t), (18a) ¯n(x,t)=n+δtQm−L(n−neq+δt2Qm)+c2Y(m+¯m2−meq), (18b)

where Eqs. (18a) and (18b) represent the streaming process in velocity space and the collision process in moment space, respectively. The moment of total energy DF in Eq. (18b) is given as , the post-collision total energy DF in Eq. (18a) is obtained via , and the post-collision moment is computed by Eq. (18b). Here, the dimensionless transformation matrix is also given by Eq. (3). On the RHS of Eq. (18b), the last density-DF-related term is introduced to consider the viscous dissipation, in which is a matrix that will be discussed and determined later. By definition, the macroscopic total energy is given as

 ρE=∑igi+δt2qe, (19)

where is the equivalent source term. Then, the total enthalpy and the temperature can be determined via the thermodynamic relations and (a function of internal energy and density ), respectively. In the present work, a simple relation , though it strictly holds only for the ideal gases, is adopted for the sake of simplicity, and more general or empirical relations can be adopted as required by specific applications. Here, is the specific heat at constant volume.

To recover the targeted energy conservation equation, as well as inspired by the ideas of our previous works on solid-liquid phase change Huang2015 ; Huang2016.amr , the equilibrium moment function for total energy DF is devised as

 neq=[ρE,−4(ρE−ρ0Cp,0T)+γ1ρ0Cp,0T,4(ρE−ρ0Cp,0T)+γ2ρ0Cp,0T,ρH^ux,−ρH^ux,ρH^uy,−ρH^uy,0,0]T, (20)

where and are the reference density and the reference specific heat at constant pressure, respectively, and and are the coefficients related to the heat conductivity. Similarly to , the discrete source term in moment space is devised as

 Qm=(qe,γ1qe,γ2qe,qe^ux,−qe^ux,qe^uy,−qe^uy,0,0)T. (21)

To avoid the deviation term caused by the convection term recovered at the order of in the diffusion term recovered at the order of , the collision matrix in moment space is modified as follows Huang2014.cde

 L=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣σ0000000000σe000000000σε000000000σjσqωj00000000σq000000000σjσqωj00000000σq000000000σp000000000σp⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (22)

where .

Since the viscous stress tensor is only related to , , and (see Eq. (15)), the matrix in the density-DF-related term, which is introduced in Eq. (18b) to consider the viscous dissipation, is set as follows

 Y=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0000000000000000000000000000Y3100000Y37Y380Y4100000Y47Y480Y5100000Y57Y580Y6100000Y67Y68000000000000000000⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (23)

where and for , , and . Through the Chapman-Enskog analysis (see A), the nonzero elements in can be determined as follows

 Y31=^ux/3,Y37=^ux,Y38=2^uy,Y51=^uy/3,Y57=−^uy,Y58=2^ux. (24)

Then, the following macroscopic conservation equation can be recovered

 (25)

Compared with Eq. (14), the heat conductivity is given as . It can be seen from Eq. (25) that the viscous dissipation and compression work are correctly considered.

Before proceeding further, some discussion on the present LB model for coupled thermo-hydrodynamic flows is in order. First, the MRT collision scheme is employed in both the LB equations for density and total energy DFs, and the collision matrix in moment space is modified to be a nondiagonal matrix rather than being set as the conventional diagonal matrix. Second, the Prandtl number can be arbitrarily adjusted. Here, is the specific heat at constant pressure, and is the dynamic viscosity. Third, the specific heat ratio can also be arbitrarily adjusted. Note that depends on the adopted EOS, and holds only for the ideal-gas EOS. Lastly, and most importantly, an arbitrary EOS (including the nonideal-gas EOS) can be prescribed, and the built-in variable is inversely calculated via with .

## 3 Boundary condition treatment

In real applications, the boundary conditions are usually given in terms of the macroscopic variables, and thus additional treatment is required to obtain the mesoscopic DFs at the boundary node. In this section, we propose the boundary condition treatment for simulating coupled thermo-hydrodynamic flows.

### 3.1 Macroscopic variables

For the velocity field, the nonslip velocity boundary condition is considered and the velocity on the boundary is directly specified. Due to the full coupling of thermo-hydrodynamic effects, the density may significantly vary near the boundary and also has a direct effect on the heat transfer process. Thus, it is important to ensure the mass conservation at the boundary node for simulating coupled thermo-hydrodynamic flows. In the present boundary condition treatment, the boundary node is exactly placed on the wall boundary, as shown in Fig. 1. The post-collision density DF hitting the wall (i.e., streaming out of the computational domain) reverses its direction as follows

 f¯i,temp(xb,t+δt)=¯fi(xb,t), (26)

where means , and the subscript “temp” implies that the density DF is temporary. After this “bounce-back” process, all the unknown density DFs at and due to the absence of adjacent nodes are now obtained. Then, the density can be computed via definition as usual (i.e., ). Note that the velocity is directly specified. Obviously, the local conservation of mass can be strictly satisfied at the boundary node. Figure 1: Schematic of boundary condition treatment with xb denoting the boundary node placed on the wall, xf and xff denoting the nearest and next-nearest fluid nodes in the normal direction, and the solid and dashed arrow lines denoting the known and unknown distribution functions after the streaming process.

As for the temperature field, the Dirichlet boundary condition with specified temperature and the Neumann boundary condition with zero heat flux (i.e., the adiabatic boundary condition) are considered. For the Dirichlet boundary condition, since the temperature is directly specified, all the involved macroscopic variables, such as the total energy , the pressure , and the total enthalpy , can be determined via the corresponding thermodynamic relations. For the Neumann boundary condition with zero heat flux, the post-collision total energy DF hitting the wall reverses its direction as follows

 g¯i,temp(xb,t+δt)=¯gi(xb,t), (27)

and thus all the unknown total enthalpy DFs at and due to the absence of adjacent nodes are temporarily obtained. Then, the total energy can be computed via definition as usual (i.e., ), and all the involved macroscopic variables, such as the temperature , the pressure , and the total enthalpy , can be determined via the corresponding thermodynamic relations.

### 3.2 Density and total energy DFs

At the boundary node, the unknown density and total energy DFs obtained via Eqs. (26) and (27) are only used to compute the macroscopic density and total energy. In the present boundary condition treatment, all the known and unknown DFs at the boundary node will be updated to make sure that the defining equations of density, velocity, and total energy (i.e., Eqs. (4) and (19)) exactly hold at the boundary node. For this purpose, we decompose the moment of DF (the same as the DF) into its equilibrium, force (source), and nonequilibrium parts, i.e.,

 m=meq−δt2Fm+mneq, (28a) n=neq−δt2Qm+nneq. (28b)

Note that the present nonequilibrium part () in Eq. (28) is different from the previous nonequilibrium part defined as () Kruger2017 when the force (source) term exists. Since the equilibrium parts ( and ) and the force (source) parts ( and ) are determined by the macroscopic variables, and can be directly computed. As to the nonequilibrium parts ( and ) at and , extrapolations are employed following the idea of the nonequilibrium-extrapolation approach Guo2002a ; Guo2002b . However, instead of simply extrapolating and , we introduce the following terms

 ~mneq=(I−S2)mneq, (29a) ~nneq=(I−L2)nneq+c2Y2(I−S2)mneq, (29b)

where is the identity matrix; and then the first- and second-order nonequilibrium extrapolations are given as

 (30a) (30b)

where and denote the nearest and next-nearest fluid nodes in the normal direction, as shown in Fig. 1. Based on our numerical tests, the first-order extrapolation has better stability but lower accuracy than the second-order extrapolation. Note that although the present collision matrices and are nondiagonal, and in Eq. (29) are still invertible, and their inverse matrices are given in B. Therefore, Eq. (30) is compatible and can be easily implemented due to the special forms of and . Moreover, different from the previous nonequilibrium-extrapolation approach Guo2002a ; Guo2002b ; Guo2007 , the present boundary condition treatment is applicable to the situation when the dynamic viscosity and thermal conductivity significantly vary with temperature and hence with space because the collision matrices are considered in the present extrapolations of nonequilibrium parts.

## 4 Validations and discussions

In this section, simulations of thermal Poiseuille and Couette flows are first carried out to validate the present LB model with self-tuning EOS for coupled thermo-hydrodynamic flows. Three different EOSs, including the decoupling EOS, the ideal-gas EOS, and the Carnahan-Starling EOS for rigid-sphere fluids Carnahan1969 , are adopted, which are given in order as follows

 p\tiny EOS=ρRgT0, (31a) p\tiny EOS=ρRgT, (31b) p\tiny EOS=ZρRgTwithZ=1+bρ/4+(bρ/4)2−(bρ/4)3(1−bρ/4)3, (31c)

where is the compressibility factor with the coefficient set to here. Then, the present LB model is applied to the simulation of natural convection in a square cavity with a large temperature difference. The ideal-gas EOS is adopted and the Rayleigh number varies from up to . In the following simulations, , , and are chosen. The relaxation parameters in satisfy , , and Huang2016.3rd , and the relaxation parameters in satisfy , , , and Huang2016.amr . Meanwhile, the ratio of bulk to kinematic viscosity is fixed at unless otherwise stated.

### 4.1 Thermal Poiseuille flow

The thermal Poiseuille flow, driven by a constant force between two parallel walls, is first simulated. Both the lower and upper walls are at rest, and the temperature of the lower and upper walls are kept at and (), respectively. The Prandtl number , the specific heat at constant pressure , and the dynamic viscosity are assumed to be constant. Thus, the analytical solutions for velocity and temperature are given as Guo2007

 uxU0=4yW(1−yW),uyU0=0, (32a) T−TcTh−Tc=PrEc3[1−(1−2yW)4]+yW, (32b)

where is the channel width, is the maximum velocity, and is the Eckert number. As seen in Eq. (32), the analytical solutions for velocity and temperature are fully determined by and . However, the analytical solution for density further depends on both the initial state and the adopted EOS. In the simulations, the density, velocity, and temperature are initialized as , , and (), respectively, and the initial pressure is determined by the adopted EOS. Thus, for the decoupling EOS (i.e., Eq. (31a)), the analytical solution for density can be easily obtained as

 ρρ0=1; (33a) for the ideal-gas EOS (i.e., Eq. (31b)), the analytical solution for density is given as ρρ0=AT0T, (33b) where the coefficient A−1=∫W0T0/Tdy; as for the Carnahan-Starling EOS (i.e., Eq. (31c)), the analytical solution for density satisfies W∫0ρdy=ρ0Wandp∞=ZρRgT, (33c)

where is the final pressure in the channel. Although an explicit expression for cannot be derived from Eq. (33c), can be easily obtained with high precision using numerical integration.

In the simulations, the lattice sound speed is set as

 cs=√∂ρp\tiny EOS∣∣∣ρ=ρ0,T=T0, (34)

and the specific heats at constant pressure and volume are fixed at

 Cp=3.5c2sT0,Cv=Cp−p\tiny EOS(ρ0,T0)ρ0T0. (35)

With this configuration, the specific heat ratio is for the decoupling and ideal-gas EOSs and for the Carnahan-Starling EOS. The simulations are carried out on a grid with lattice spacing and periodic boundary in . The lower and upper walls are treated by the present boundary condition treatment with second-order extrapolation. The basic parameters are set as , , and , and the dimensionless relaxation time for density DF, defined as , is fixed at for the kinematic viscosity .

Fig. 2 shows the velocity and temperature distributions across the channel and compares the numerical results with the analytical solutions given by Eq. (32). Two sets of and are considered here: for the first set, is fixed at and varies from to ; while for the second set, is fixed at and varies from to . As an important computational parameter, the lattice Mach number, defined as , is fixed at in the simulations. Good agreement between the numerical results and the analytical solutions can be observed in Fig. 2, which demonstrates that the effects of the viscous dissipation and compression work are successfully captured by the present LB model. From Fig. 2, we can also see that the distributions of and obtained with different EOSs are almost identical, which agrees with the aforementioned discussion. Note that the simulation with and loses numerical stability for the ideal-gas EOS. To further validate the present LB model with self-tuning EOS, comparisons of the density distributions are carried out in Fig. 3. Good agreement is observed between the numerical results obtained with different EOSs and the corresponding analytical solutions given by Eq. (33), which demonstrates that various EOSs (including the nonideal-gas EOS) can be handled by the present LB model. For the decoupling EOS, keeps constant across the channel; as for the ideal-gas and Carnahan-Starling EOSs, varies across the channel due to the full coupling of thermo-hydrodynamic effects. In the Carnahan-Starling EOS, the molecular volume is considered, which implies that the rigid-sphere fluid is less compressible than the corresponding ideal gas. Therefore, the variation in density across the channel obtained with the Carnahan-Starling EOS is smaller than that obtained with the ideal-gas EOS, as clearly shown in Fig. 3. Figure 2: Comparisons of velocity ux/U0 (left) and temperature (T−Tc)/(Th−Tc) (right) distributions between the numerical results and the analytical solutions for thermal Poiseuille flow with (a) Pr=0.71 and Ec varying from 0.1 to 100, and (b) Ec=10 and Pr varying from 0.1 to 4. Figure 3: Comparisons of density ρ/ρ0 distributions between the numerical results obtained with different EOSs and the corresponding analytical solutions for thermal Poiseuille flow with (a) Pr=0.71 and Ec=20, and (b) Ec=10 and Pr=1.

Considering that the lattice Mach number plays an important role in the LB method, we further investigate the accuracy of the present simulation with respect to . It is worth pointing out that the lattice Mach number is not only a computational parameter but also closely related to the Eckert number (or the real Mach number) due to the lattice sound speed given by Eq. (34). In the following simulations, the Prandtl number is fixed at , and the Eckert number is set as with