## Abstract

The lattice Boltzmann method (LBM) has been proved to be a useful tool in many areas of computational fluid dynamics, including computational aero-acoustics (CAA). However, for historical reasons, its applications in CAA have been largely restricted to simulations of isothermal (Newtonian) sound waves. As the recent kinetic theory-based reformulation establishes a theoretical framework in which LBM can be extended to recover the full Navier–Stokes–Fourier (NS) equations and beyond, in this paper, we show that, at least at the low-frequency limit (sound frequency much less than molecular collision frequency), adiabatic sound waves can be accurately simulated by the LBM provided that the lattice and the distribution function ensure adequate recovery of the full NS equations.

## 1. Introduction

Historically, the lattice Boltzmann method (LBM) was derived as the kinetic description of the lattice cellular automaton of Frisch *et al.* [1] with simplified collision models, e.g. the single-relaxation-time model of Bhatnagar *et al.* [2] introduced [3–6] to recover the Navier–Stokes–Fourier (NS) equations in the small-Mach-number (near-incompressible) limit. Owing to its kinetic nature, the LBM intrinsically includes sound waves in its macroscopic dynamics and has recently been found to be a competitive method in aero-acoustic applications by Crouse *et al.* [7] and Marié *et al.* [8].

Besides the conceptual near-incompressible assumption, two obstacles, both related to the *a posteriori* manner in which the lattice BGK (Bhatnagar–Gross–Krook) models were derived, were met when the LBM was applied to adiabatic acoustics and thermodynamic problems in general. First, as shown by Qian & Orszag [9], the original LBM model has an error term cubic in the Mach number that was only fully removed recently [10,11]. Second, inclusion of energy conservation in the microscopic dynamics is found to have caused numerical instability. Although it is known from experience that stability improves as more velocities are included, it was not known at the time how this can be achieved systematically [12]. Lallemand & Luo [13] alleged that the difficulties with the thermal equation are intrinsic to the LBM and that the energy-conserving LBM is not suitable for numerical simulation. Li *et al.* [14] resorted to a finite-difference implementation to circumvent the difficulties. Buick & Cosgrove [15] used additional body force terms to achieve a variable speed of sound. More recently, a more general compressible LBM model was developed using additional velocities by Chen *et al.* [16].

The LBM was shown early in its development to be a velocity-space discretization of the Boltzmann BGK equation [17–20], although the implications were not immediately comprehended by all authors. As suggested by Grad [21], by noticing that the distribution function only directly contributes to the hydrodynamics through its leading moments, one can describe a fluid system by the leading moments instead of the entire distribution. Furthermore, by using the one-to-one correspondence established between the moments and discrete function values of a polynomial by the Gauss–Hermite quadratures, one can formulate a lattice BGK (LBGK) equation as a projection of the continuum BGK equation in a truncated Hilbert space spanned by the leading Hermite polynomials. The truncated kinetic equation will have the same hydrodynamics if the relevant moments are preserved in the equilibrium distribution and the underlying lattice forms a sufficiently accurate Gauss–Hermite quadrature in the velocity space [20]. In particular, the full NS equation is recovered beyond the small-Mach-number limit if the equilibrium distribution is taken to be the third-order Hermite expansion of the Maxwellian and the lattice forms a quadrature that is at least sixth-order accurate. Nie *et al.* [11] demonstrated that the cubic error term can be removed and the Galilean invariance fully restored once these conditions are met.

In the present work, we present further numerical validations that adiabatic sound can be adequately simulated using a high-order LBM. The paper is organized as follows: In §2, the kinetic theory behind the LBM is briefly reviewed and results relevant to computational aero-acoustics (CAA) outlined. A sequence of LBM models with increasing orders of accuracy are also introduced. Numerical simulations of typical benchmark problems are presented in §3, and conclusions and further discussion offered in §4.

## 2. Lattice Boltzmann hydrodynamics and acoustics

Following Shan *et al*. [20], the *D*-dimensional LBGK equation is the discretized version of the following continuum BGK equation, which describes the evolution of the single-particle distribution, *f*, in phase space (** x**,

**): 2.1 Here**

*ξ***,**

*x***and**

*ξ**t*are, respectively, the spatial coordinates, microscopic velocity and time, and

*τ*is the relaxation time. The external body force terms are ignored for simplicity. Finally,

*f*

^{(eq)}is the following dimensionless Maxwell–Boltzmann distribution: 2.2 where

*ρ*and

**are the macroscopic density and fluid velocity, and**

*u**θ*is the temperature, which is related to the internal energy density per mass,

*ϵ*, by

*ϵ*=

*Dθ*/2. Defining the intrinsic velocity

**≡**

*c***−**

*ξ***, then**

*u**ρ*,

**and**

*u**ϵ*are all moments of the distribution: 2.3 By direct integration, they satisfy the conservation equations of mass, momentum and energy: 2.4a 2.4b 2.4c Here d/d

*t*≡∂/∂

*t*+

**⋅∇ is the Lagrangian derivative, and**

*u***and**

*P***are, respectively, the molecular pressure tensor and heat flux, which are also moments of the distribution: 2.5 Clearly,**

*q**P*

_{ii}=2

*ρϵ*=

*Dρθ*. The hydrostatic pressure is then given by

*p*≡

*P*

_{ii}/

*D*=

*ρθ*, which is the equation of state of an ideal gas. Denoting the dimensionless entropy per mass as

*s*, from the fundamental thermodynamic equation, we have 2.6 The heat capacities at constant volume and pressure are thus

*c*

_{v}=

*D*/2 and

*c*

_{p}=(

*D*+2)/2, and the ratio of specific heats

*γ*≡

*c*

_{p}/

*c*

_{v}=(

*D*+2)/

*D*. Writing

**=**

*P**p*

**+**

*δ***, where**

*σ***is the traceless stress tensor, equation (2.4c) can also be written as 2.7**

*σ*The hydrodynamic equations (2.1) are obtained by the Chapman–Enskog asymptotic approximations [22]. At the zeroth order, we have *f*=*f*^{(eq)}, which yields ** σ**=

**=0. Equations (2.4) become Euler’s equations. At the next order, we have 2.8 where**

*q**μ*and

*λ*are, respectively, the dynamic viscosity and thermal conductivity: 2.9 The bulk viscosity is zero.

In the case that *θ* is a constant, only equations (2.4a) and (2.4b) are valid. As Dellar [23] pointed out, the stress tensor is
2.10
corresponding to a shear viscosity of *μ* and a bulk viscosity of *ζ*=2*μ*/*D*.

Consider a small planar-wave perturbation of the form , where *ω* is the frequency and ** k** is the wavevector. Linear analysis on equations (2.4) yields the four well-known hydrodynamic modes, i.e. the viscous, thermal and two acoustic modes. The viscous mode is independent of the rest and has a simple dispersion relation of

*ω*=−

*νk*

^{2}, where

*ν*=

*μ*/

*ρ*is the kinematic shear viscosity. The other three modes are coupled together. We define the acoustic Reynolds number

*Re*=

*c*

_{s}/

*νk*, the Prandtl number

*Pr*=

*ν*/

*κ*and the Péclet number

*Pe*=

*PrRe*, where

*κ*=

*λ*/

*ρc*

_{p}is the thermometric heat conductivity and

*c*

_{s}is the sound speed at the long-wavelength limit. With full thermodynamics, the dispersion relations for the viscous, thermal and acoustic modes have, respectively, the following dimensionless forms: 2.11a 2.11b 2.11c where . For isothermal sound, the dispersion relations of the viscous and acoustic modes are 2.12 where . Note that, with the BGK collision model,

*κ*=

*ν*and

*Pr*=1. It is easy to see from equation (2.11) that adiabatic sound has formally the same attenuation rate as isothermal sound despite the different sound speed. To truly simulate adiabatic sound waves at non-unity Prandtl numbers, we employ the general multiple-relaxation-time model of Shan & Chen [24] in which the third moments of the distribution are allowed to relax independent of the second moments.

As shown previously by Shan *et al.* [20] and Nie *et al.* [11], the LBGK equations approximate equation (2.1) by expanding it in Hermite polynomials and evaluating on a set of discrete velocities that forms a quadrature in ** ξ**-space. The order of expansion,

*N*, and the algebraic precision of the quadrature,

*Q*, jointly determine the order of the highest hydrodynamic moments, and thus the corresponding physics, which are accurately retained in the discrete system. To recover the correct oscillatory behaviour of the sound waves, isothermal or adiabatic, it is sufficient to retain the dynamics of the second moments. The popular D

*x*Q

*x*models of Qian

*et al.*[4], all employing a second- or third-order expansion and a fifth-order lattice, are adequate for that purpose. The attenuation of isothermal sound is solely determined by viscous dissipation, which the D

*x*Q

*x*models are also capable of simulating correctly at the small-Mach-number limit. As the simulation of the molecular transport of energy and of momentum at finite Mach number calls for accurate recovery of the dynamics of the third moment, a third-order expansion of the Maxwellian and a sixth-order accurate lattice are required for adiabatic sound attenuation to agree with the continuum results above at small Mach number and constant flow temperature. To fully recover the attenuation of adiabatic sound waves without the Mach number and temperature limitation, a fourth-order expansion and a ninth-order lattice are required.

Although both high-order on- and off-lattice LBM models are known, to avoid the additional interpolation error of the off-lattice models, we focus on the benchmark of on-lattice models in this study. We follow the D*x*Q*y* convention of naming the lattices by their dimensionality, *x*, and number of velocities, *y*. A total of eight models are evaluated in simulations of linear mode dynamics. They are the conventional second-order models D2Q9, D3Q15 and D3Q19, third-order models D2Q17, D2Q21 and D3Q39, and fourth-order models D2Q37 and D3Q121, among which D2Q21 and D2Q37 are, respectively, the two-dimensional projections of D3Q39 and D3Q121. Here, the order of a model is the order of the highest moments retained. The velocity sets of D2Q17, D2Q21 and D2Q17 are illustrated in figure 1, and more details about these models can be found in Shan [25].

## 3. Numerical verifications

The linear modes given in equations (2.11)–(2.12) are first verified by measuring the decay rate and/or oscillation frequency in simulations performed on a periodic lattice represented by {(*x*_{1},…,*x*_{D}):*x*_{i}=1,…,*L*_{i}}, where *L*_{i} is the number of grids in the *i*th dimension. Decomposing the velocity field into components parallel and perpendicular to the wavevector, denoted by *u*_{∥} and *u*_{⊥}, respectively, and defining ** R**≡(

*ρ*,

*u*

_{∥},

*u*

_{⊥},

*θ*), the linear waveform combining all modes can be written as 3.1 where the subscript

_{0}represents quantities of the base flow;

*ϵ*is the amplitude of the perturbation, and

**′ is a superposition of all eigenmodes with the dynamics of 3.2 where , and are, respectively, the eigenvectors of the viscous, thermal and standing acoustic modes;**

*R**η*

_{v},

*η*

_{t}and

*η*

_{s}are the decay rates;

*a*

_{v},

*a*

_{t}and

*a*

_{s}are the initial amplitudes; and

*ω*is the angular frequency of the acoustic wave. The eigenvectors of the thermal and acoustic modes correspond, respectively, to homogeneous pressure and entropy in the initial condition at leading order. The wavevector is given by

*k*

_{i}=2

*π*/

*cL*

_{i}, where

*c*is the lattice constant. Waves arbitrarily oriented with respect to the grid can be simulated by choosing the grid dimensions. Whereas the viscous mode is independent of the rest and an exact solution to the Navier–Stokes equation for any

*ϵ*, the thermal and acoustic modes are coupled approximate solutions only accurate at small amplitudes (

*ϵ*≪1). To set up the initial condition for the thermal and acoustic modes, we choose

*a*

_{t}=1,

*a*

_{s}=0 and

*a*

_{t}=0,

*a*

_{s}=1, respectively. Nevertheless, owing to mode coupling and the approximate nature of the solution, the thermal and acoustic modes could be excited both at the base wavelength and at higher spatial harmonics even if they are absent initially. To accurately extract the decay rate and the oscillation frequency, a spatial discrete Fourier transform is performed on the density field at each time step to obtain the time history of the amplitude, which is then fitted in the time domain using the Levenberg–Marquardt algorithm according to equation (3.2) to obtain the decay rates and frequency.

We first measured the frequency and attenuation of the standing acoustic mode using all the models in the previous section for both isothermal and thermal cases and relaxation time in the range 0.001≤*τ*≤5. The simulations were performed on a grid of 100×1×1, with the sound wave propagating along the *x*-direction. Plotted in figure 2 is the ratio between the measured sound speed and the theoretical predictions of equations (2.11) and (2.12) against the acoustic Reynolds number, which is proportional to the inverse sound frequency against the molecular collision frequency. To be noted first is that, at large Reynolds number, corresponding to the Navier–Stokes limit where the sound frequency is much lower than the collision frequency, all models predict sound speed in excellent agreement with linear theory. This is not surprising, as only the hydrodynamics at the level of the Euler equation is needed to predict the correct sound speed. At the sound frequency approaches the collision frequency, all models show an increase in sound speed in agreement with the experiments of Greenspan [26] and the direct simulation Monte Carlo (DSMC) method of Hadjiconstantinou & Garcia [27]. To be noted is that, whereas in the isothermal case results using all models collapse quite well onto a single curve as a function of Reynolds number, in the thermal case the results of the second-order models show a deviation among themselves and from those of the higher-order models. The results of the high-order models agree with each other in the thermal case.

Similarly shown in figure 3 are the measured sound attenuation against theoretical predictions. At large Reynolds numbers, all models give the correct sound attenuation in the isothermal case whereas only the high-order models give the correct sound attenuation in the adiabatic case. This is also expected, as sound attenuation is a consequence of viscous dissipation for isothermal sound but a combined effect of both viscous and thermal dissipation for adiabatic sound. Correct thermal transport is therefore required to have the right sound attenuations. At lower Reynolds numbers, the high-order models show a decrease in sound attenuation in both the isothermal and adiabatic cases, which is again in agreement with experiments and DSMC simulations. On the contrary, as also shown by Dellar [23], the second-order models show an increase in sound attenuation for isothermal sound.

We then investigate the isotropy of the thermal modes, as Lallemand & Luo [13] alleged that energy-conserving LBM models are *not* suitable for numerical simulations owing to a universal mode coupling between the viscous and thermal modes. The simulations are conducted with the grid dimensions so chosen that the wavevector forms the desired angle with the *x*-axis and has the specified norm. The same viscosity (*ν*=0.01) and Prandtl number (*Pr*=0.71) as in Lallemand & Luo [13] are used in all the simulations, and the wavenumber is *π*/15*c*, where the model-dependent lattice constant *c* ranges between 1.19 and 1.73. Shown in figure 4*a* is the ratio between the measured thermal diffusivity and its theoretical value against the angle that the wavevector forms with respect to the *x*-axis, and in figure 4*b* the time histories of the perturbation amplitudes of the viscous and thermal modes as simulated using D2Q21 are shown. Except for the small initial noise, the time histories agree with the linear theory prediction given in equation (3.2). The oscillation in the thermal mode is due to the coupling with the acoustic modes. Both isotropy and agreement with theory are good and no abnormal ‘mode coupling’ is observed.

As a typical benchmark problem in computational acoustics, the propagation of a two-dimensional Gaussian acoustic pulse, both isothermal and adiabatic, is simulated. The initial perturbation to the otherwise static flow field has the form of and ** u**′=0, where , with

*b*being the half-width of the Gaussian pulse and

*r*the distance to the pulse centre. In addition, for adiabatic sound, the temperature is set to

*θ*=

*θ*

_{0}(

*ρ*/

*ρ*

_{0})

^{2/D}, corresponding to a homogeneous initial entropy. The analytical inviscid solution is found to be 3.3 where

*j*

_{0}(

*ξr*) is the zeroth-order Bessel function of the first kind. Simulations are performed on a two-dimensional 192×192 grid using D2Q9 for isothermal sound and D2Q21 for adiabatic sound. Viscosity is set to the very small value of

*ν*=10

^{−6}and

*ϵ*=10

^{−4}. Shown in figure 5 is the pressure on all grid points plotted against

*r*at a few typical time steps. Simulation results in both the isothermal and adiabatic cases are in excellent agreement with the analytical solutions.

## 4. Conclusions and discussion

We present initial numerical evidence that adiabatic sound can be adequately simulated by the LBM provided sufficient hydrodynamic moments are retained in velocity-space discretization and the sound frequency is much smaller than the molecular collision frequency. Previous difficulties in simulating adiabatic sound are due to insufficient accuracy rather than to intrinsic limitations of the LBM. At higher frequencies where kinetic effects start to emerge, LBM results agree at least qualitatively with experiments and DSMC simulations. More detailed studies at higher frequencies, in variable temperature environments, and convergence are necessary and will be presented in future publications.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society