Numerical analysis of the lattice Boltzmann method for simulation of linear acoustic waves

TL;DR

This paper analyzes a lattice Boltzmann method for simulating linear acoustic waves, comparing different collision operators and lattice models, and demonstrating its efficiency and accuracy in heterogeneous media.

Contribution

It introduces a detailed dispersion relation analysis for LB models and implements a grid-refinement algorithm, showing LB's competitiveness with finite-difference schemes.

Findings

D2Q5 lattice is most suitable for linear acoustics

LB scheme performance is comparable to finite-difference methods

LB method is efficient for complex geometries and parallel computing

Abstract

We analyse a linear lattice Boltzmann (LB) formulation for simulation of linear acoustic wave propagation in heterogeneous media. We employ the single-relaxation-time Bhatnagar-Gross-Krook (BGK) as well as the general multi-relaxation-time (MRT) collision operators. By calculating the dispersion relation for various 2D lattices, we show that the D2Q5 lattice is the most suitable model for the linear acoustic problem. We also implement a grid-refinement algorithm for the LB scheme to simulate waves propagating in a heterogeneous medium with velocity contrasts. Our results show that the LB scheme performance is comparable to the classical second-order finite-difference schemes. Given its efficiency for parallel computation, the LB method can be a cost effective tool for the simulation of linear acoustic waves in complex geometries and multiphase media.

Figures

a](#S3.F7.sf1)) of few grid points where the coarse domain extends into the fine domain and vice versa. Streaming distributions coming from either domains, instead of stopping at the interface, proceed in the overlap region on respective grids. The overlap regions consist nodes which form sponge layers Colonius ([2004](#bib.bib37)). At a sponge node, during each collision, the distribution function is reduced by a fraction of its value. Thus, the sponge layers gradually absorb the wave as it propagates in the overlap region. This results in negligible reflections back into the physical domain.

a](#S5.F12.sf1) shows simulations in a homogeneous medium with grid refinement. The source is located on the left i.e. coarse-grid. The BGK-LB scheme on the D2Q5 lattice is used. Waves propagate smoothly across the interface without any artificial reflection. The only difference between the left and the right is the resolution, which is $8$ points per wavelength on the coarse side and $12.8$ points per wavelength on the fine side. Note that the simulated wavefront is broader in the coarse domain compared to [Fig. 9a](#S5.F9.sf1). This is because resolution in the coarse domain is lowered by a

Tables

Table 1. Table 1: Lattice velocities

Lattice	Lattice velocities
D2Q5	$c_0=(0,0),$ $c_{1,3},c_{2,4}=(\pm 1,0),(0,\pm 1).$
D2Q9	$c_0=(0,0),$ $c_{1,3},c_{2,4}=(\pm 1,0),(0,\pm 1),$ $c_{5,7,6,8}=(\pm 1,\pm 1).$

Equations

$$\frac{1}{c_s^2(\text{x})}{\partial }_t^{2}p(\text{x},t)-{\nabla }^{2}p(\text{x},t)={\partial }_{t}S(\text{x},t).$$

$$\partial_{t}p(\textbf{x},t)+c_{s}^{2}(\textbf{x})\nabla\cdot\big{[}\rho_{0}(\textbf{x})\textbf{v}(\textbf{x},t)\big{]}=c_{s}^{2}(\textbf{x})S(\textbf{x},t),$$

$$\partial_{t}\big{[}\rho_{0}(\textbf{x})\textbf{v}(\textbf{x},t)\big{]}+\nabla p(\textbf{x},t)=0.$$

$$p\left(\text{x},t\right)=c_s^2\left(\text{x}\right)\rho \left(\text{x},t\right).$$

$${\partial }_{t}g(\text{x},\text{c},t)$$

$$=$$

$$g_i(\text{x}+{\text{c}}_i\delta t,t$$

$$=$$

$$g_i^{eq}(\text{x},t)=$$

$$+\frac{{\rho }_0(\text{x})}{2c_s^2(\text{x})}[{\text{c}}_i\cdot \text{v}(\text{x},t){]}^2$$

$$\displaystyle-\frac{\rho_{0}(\textbf{x})}{2}||\textbf{v}(\textbf{x},t)||^{2}\Big{\}},$$

$$\rho (\text{x},t)=\sum _i{g}_i(\text{x},t)=\sum _i{g}_i^{eq}(\text{x},t),$$

$${\rho }_0(\text{x})\text{v}(\text{x},t)=\sum _i{\text{c}}_i{g}_i(\text{x},t)=\sum _i{\text{c}}_i{g}_i^{eq}(\text{x},t).$$

$$\text{g}(\text{x}+\text{c}\delta t,$$

$$=$$

$$\text{m}=\{\rho ,{\rho }_0{v}_x,{\rho }_0{v}_y,e,p_{xx}\},$$

$$\text{m}=\{\rho ,{\rho }_0{v}_x,{\rho }_0{v}_y,e,p_{xx},ϵ,q_x,q_y,p_{xy}\}.$$

$$\text{g}(\text{x}+\text{c}$$

$$=$$

$$g_i(\text{x}+{\text{c}}_i\delta t,$$

$$=$$

$$+w_{i}S(\text{x},t),$$

$$g_{i}^{eq}(\textbf{x},t)=\frac{w_{i}}{c_{s}^{2}(\textbf{x})}\Big{[}p(\textbf{x},t)+\rho_{0}(\textbf{x})v_{\alpha}(\textbf{x},t)c_{i\alpha}\Big{]}.$$

$$\nu (\text{x})=c_s^2(\text{x})\left(\tau -\frac{1}{2}\right)\delta t.$$

$$\begin{array}{rl}\sum _i{w}_i(\text{x})& =1,\\ \sum _i{w}_i(\text{x})c_{i\alpha }{c}_{i\beta }& =c_s^2(\text{x}){\delta }_{\alpha \beta }.\end{array}$$

$$w_0(\text{x})=1-\eta (\text{x}{)}^2,$$

$$e^{(0)}=$$

$$=$$

$$p_{xx}^{(0)}=g_1^{(0)}-g_2^{(0)}+g_3^{(0)}-g_4^{(0)}=0.$$

$$g_{i}^{eq}=w_{i}\Big{\{}\sum_{j}g_{j}+\big{(}1/c_{s}^{2}\big{)}\big{[}c_{ix}(g_{1}-g_{3})+c_{iy}(g_{2}-g_{4})\big{]}\Big{\}},$$

$$g_i^{eq}=$$

$$\displaystyle+c_{iy}(g_{2}-g_{4}+g_{5}+g_{6}-g_{7}-g_{8})\big{]}\Big{\}}.$$

$$\text{g}(\text{x},t+\delta t)$$

$$=$$

$$\mathrm{A}=\left[\begin{array}{ccccc}1-w_0& -w_0& -w_0& -w_0& -w_0\\ -w_1& 1-w_1(1+1/c_s^2)& -w_1& -w_1(1-1/c_s^2)& -w_1\\ -w_2& -w_2& 1-w_2(1+1/c_s^2)& -w_2& -w_2(1-1/c_s^2)\\ -w_3& -w_3(1-1/c_s^2)& -w_3& 1-w_3(1+1/c_s^2)& -w_3\\ -w_4& -w_4& -w_4(1-1/c_s^2)& -w_4& 1-w_4(1+1/c_s^2)\end{array}\right].$$

Full text

Lattice Boltzmann method \sepseismic wave \sepdispersion analysis \sepCourant number \sepgrid-refinement

I Introduction

Over the past three decades, the lattice Boltzmann (LB) method has been established alongside conventional CFD methods as an efficient scheme for the numerical solution of partial differential equations. A variety of complex flow problems Benzi et al. (1992); Aidun and Clausen (2010) have been tackled successfully, even within the simplest Bhatnagar-Gross-Krook (BGK) framework Bhatnagar et al. (1954). It has proven particularly useful for applications involving flows through porous and multi-phase media Chen and Doolen (1998); Illner et al. (2007); Succi (2001); Perlekar et al. (2014). The method has been extended well beyond hydrodynamics to solve the governing equations of magnetohydrodynamics Chen et al. (1991), acoustic and electromagnetic wave propagation Chopard and Luthi (1999); Hanasoge et al. (2011), fracturing in solids Chopard and Droz (1998), the Schrödinger equation Palpacelli and Succi (2007) etc. With a generalized multi-relaxation-time (MRT) collision operator d’Humières (1992, 2002); Lallemand and Luo (2000), the LB scheme can be fine tuned to more accurately model the physical problem at hand.

The LB method has several distinct advantages over conventional CFD approaches Chen and Doolen (1998). The LB equation comprises a single first-order differential equation, which, in the asymptotic limit (Chapman-Enskog expansion procedure Chapman and Cowling (1970); Wolf-Gladrow (2000)), produces the macroscopic equations for mass and momentum conservation. Only immediate neighbour lattice sites interact in LB, so that the computation is highly suitable for a parallel-computing implementation. Imposing boundary conditions in LB is relatively simple, which makes it suitable to handle complex geometries.

In the current work, we discuss an application of LB to model linear wave propagation in seismology. For simplicity and clarity in understanding numerical performance, we consider a seismic wave equation of only acoustic waves (P-waves) with applications in e.g. helioseismology Gizon et al. (2010). The propagation medium for these seismic waves can be highly heterogeneous, comprising a mixture of multi-phase components in a porous environment. The geometry of the physical domain of interest may be complicated. The distinct success of the LB method for multi-phase and porous flows combined with its efficiency in handling complex boundaries makes it useful to investigate how LB simulations of acoustic waves compare with standard numerical schemes — such as finite-differences — employed in exploration seismology.

LB schemes have been analysed in the context of simulation of acoustic waves in prior literature. In Marié et al. (2009), Marié et al. investigate sources of errors in LB dispersion as compared to Navier-Stokes. In Xu and Sagaut (2011), Xu et al. propose an optimization strategy to minimize dispersion errors in the MRT-LB scheme. In both these studies, a conventional LB scheme with a non-linear equilibrium distribution, which yields the Navier-Stokes equation, i.e., on performing the Chapman-Enskog expansion, is used Wolf-Gladrow (2000). Hence, the resulting wave equation retains some non-linearity as well as viscous dissipation and does not faithfully represent the linear inviscid wave equation 1. Instead, a straightforward way to simulate the inviscid linear wave equation is to use the linear LB scheme proposed by Chopard Chopard and Luthi (1999). The linear LB scheme utilises the linear equilibrium distribution function 15 and appropriately recovers only the linear inviscid part of the wave equation (see Appendix B). In this work, we analyse the linear LB scheme for the simulation of acoustic waves. Quantifying the numerical ability of a scheme is best done on a fully linear equation, e.g. Lele (1992). Indeed, our work is relevant in this regard since it highlights unexpected strengths and weaknesses of the LB method.

Few authors have considered the linear LB scheme for the simulation of waves. Viggen Viggen (2014, 2013a, 2013b) provide a detailed account of the application of LB scheme for acoustic waves. In particular, Viggen (2014) presents a derivation of the dispersion relation for the discrete-velocity Boltzmann equation using the linear equilibrium distribution function and compares it with the Navier-Stokes equation. As we will show later, the linear equilibrium distribution function can be written as a linear combination of single-particle distribution functions. Thus we can directly obtain the numerical dispersion relations without additional linearisation procedures. We derive dispersion relations for the MRT-LB scheme on 2D lattices — D2Q5 and D2Q9. We also study the numerical anisotropy of the dispersion relation on these 2D lattices. The dispersion relations for the BGK-LB scheme are obtained as a special case when all the relaxation parameters are identical. (similar to finite-difference schemes) The numerical dispersion relation, as well as the stability of the LB scheme, are sensitive to the Courant number Robertsson and Blanch (2011). A heterogeneous medium may display large variations in sound speed, and hence the local Courant number. Thus we study LB dispersion relations at different Courant numbers and compare it with second- and fourth-order finite-difference schemes. Significant changes in sound speed result in corresponding variations in wavelength. A grid-refinement algorithm (based on Dupuis and Chopard (2003)) is presented to simulate waves with uniform grid resolution, i.e., the number of grid points that resolve the wave. Although our analysis is performed on a strictly linear problem, the numerical limits of LB that we have identified are also relevant for standard LB schemes.

In Section II, we briefly discuss the LB methodology and its application to seismic wave propagation. We describe the mathematical model for seismic waves and linear LB formulation. We also discuss the parameters which are important for performance of the LB scheme for the simulation of waves. In Section III, we derive the dispersion relation for the LB scheme and present LB dispersion curves on 2D lattices. We compare LB dispersion curves with exact as well as finite-difference schemes at different Courant numbers. In Section IV, we give a detailed account of the grid-refinement algorithm used to model waves in heterogeneous media. Finally, in Section V, we present results from simulations in homogeneous and heterogeneous media.

II Seismic waves

II.1 Macroscopic equations

Seismic P-waves in acoustic media can be described by a linear pressure wave equation with a source term

[TABLE]

This second-order equation derives from two coupled first-order equations for continuity

[TABLE]

and conservation of momentum

[TABLE]

Here $p(\textbf{x},t)$ is pressure fluctuation, $\textbf{v}(\textbf{x},t)$ is velocity fluctuation, $S(\textbf{x},t)$ is the scalar source of pressure fluctuations, and $c_{s}(\textbf{x})$ and $\rho_{0}(\textbf{x})$ are the prescribed temporally stationary sound speed and density of the background medium. The functional form of $c_{s}(\textbf{x})$ and $\rho_{0}(\textbf{x})$ are dictated by the heterogeneity of the medium and may thus be complicated. The pressure and the density fluctuations are related by a linearised ideal gas equation of state

[TABLE]

The sound speed $c_{s}(\textbf{x})$ can be modelled using an isentropic bulk modulus for the medium Krüger et al. (2017).

II.2 The LB model

Kinetic theory describes the dynamics of a gas in terms of the distribution function, which is the probability density of finding a gas particle in a differential phase-space volume. The distribution function evolves as particles move (or in LB terminology, stream) and collide with each other. In the BGK-Boltzmann transport equation, collision is modelled through a linear operator with a single relaxation time

[TABLE]

where $g(\textbf{x},\textbf{c},t)$ is the single-particle distribution function Huang (1987), $\tau$ is the relaxation time for particle collisions and $g^{eq}(\textbf{x},\textbf{c},t)$ is the equilibrium Maxwell-Boltzmann distribution.

A second-order discretisation of the BGK-Boltzmann transport equation 5 gives the BGK-LB equation Reider and Sterling (1995); He and Luo (1997)

[TABLE]

Both position space and velocity space are discretised; i.e., only a finite number of microscopic velocities are allowed. This is achieved by mapping the computational domain onto a lattice. In 6, $\textbf{c}_{i}$ is the $i^{th}$ microscopic velocity on the lattice and $g_{i}(\textbf{x},t)$ is the corresponding single-particle distribution function. For the LB scheme, the equilibrium Maxwell-Boltzmann distribution is truncated at the second order in velocity He and Luo

[TABLE]

where $w_{i}$ is the lattice weight for the $i^{th}$ microscopic velocity and $c_{s}(\textbf{x})$ is the lattice sound speed (see Section II.3). Here, the density and momentum fluctuations are obtained by taking zeroth and first microscopic velocity moments of the distribution function over the lattice velocity space,

[TABLE]

[TABLE]

The pressure fluctuations $p(\textbf{x},t)$ is then obtained by using the equation of state 4.

Depending on the problem, various 2D and 3D lattices are used in LB. The general nomenclature for lattices is DnQm where n is the dimensions of position space and m is the number of microscopic velocities available at each lattice site. For simulating 2D acoustic waves, we have used two 2D lattices - D2Q5 and D2Q9, with 5 and 9 microscopic velocities respectively. These lattices are shown in Fig. 1 and the corresponding lattice velocities are given in Table 1. In both lattices, the microscopic velocity $c_{0}$ corresponds to the particle at rest.

In the MRT-LB scheme, the different components of the distribution function g may relax to the equilibrium distribution function $\textbf{g}^{eq}$ with different relaxation parameters d’Humières (1992). The general LB equation is written for the distribution function column vector $\textbf{g}=\{g_{i}\}$ as,

[TABLE]

Here, the column vector m consists of conserved and non-conserved velocity moments of the distribution function. The orthogonal matrix $\rm{M}$ transforms the distribution function vector into the moment vector as $\textbf{m}=\rm{M}\textbf{g}$. The conserved moments density and momentum fluctuations are given by equations 8 and 9 respectively. The non-conserved moments are at higher order in velocity and are calculated using the microscopic velocity set on a given lattice d’Humières (1992, 2002); Lallemand and Luo (2000). For the D2Q5 lattice, the moment column vector is

[TABLE]

and for the D2Q9 lattice,

[TABLE]

Second-order velocity moments $e$, $p_{xx}$ and $p_{xy}$ correspond to the energy, diagonal and off-diagonal components of stress tensor respectively and $\epsilon$, $q_{x}$ and $q_{y}$ are higher-order velocity moments on lattice Lallemand and Luo (2000).

The equilibrium distribution function vector $\textbf{g}^{eq}$ transforms into the equilibrium moment vector $\textbf{m}^{eq}$ in moment space. Transformation matrices M for D2Q5 and D2Q9 lattices are specified in Appendix A. $\rm{S}$ is the diagonal relaxation matrix consisting of the inverse relaxation-times for the conserved as well non-conserved moments. For the D2Q5 lattice, $\textup{S}=\textup{diag}(s_{\rho},s_{v},s_{v},s_{e},s_{p})$ and for the D2Q9 lattice, $\textup{S}=\textup{diag}(s_{\rho},s_{v},s_{v},s_{e},s_{p},s_{\epsilon},s_{q},s_{q},s_{p})$ (where diag refers to diagonal matrix). The various inverse relaxation-times — $s$ — are related to the macroscopic properties of the fluid which govern the hydrodynamics and kinetics (Section II.3). The values of these parameters can be fine tuned to suit the macroscopic dynamics of the physical system and also to improve the numerical stability of the model Lallemand and Luo (2000). The generalised LB equation 10 can be expanded to obtain the MRT-LB equation in terms of the distribution column vector alone,

[TABLE]

Here, $\rm{C}=\rm{M}^{-1}{S}{M}$ is the collision matrix. When all the inverse relaxation-times are identical, we recover the BGK-LB equation 6.

The linear wave equation 1 can be modelled by introducing a source term in the MRT-LB equation 13 Viggen (2013b). Writing component wise, the corresponding equation is

[TABLE]

where $w_{i}$ are the lattice weights. Since the first-order macroscopic equations 2 and 3 are linear in fluctuating quantities $p(\textbf{x},t)$ and $\textbf{v}(\textbf{x},t)$, we truncate the equilibrium distribution at the linear term in $\textbf{v}(\textbf{x},t)$. Thus the numerical model of linear acoustic waves involves the following linear equilibrium-distribution function

[TABLE]

The MRT-LB equation 14 along with the linear equilibrium distribution 15 are collectively referred to as the linear LB scheme. The solution of 14 using the equilibrium distribution function 15, with proper choice of relaxation times, yields the numerical solution of the linear acoustic wave equation for pressure 1 as well as velocity. We discuss relaxation times and other LB parameters in the next section.

II.3 The LB parameters

The LB equation 14 represents the governing equation for linear acoustic waves 1 provided we correctly adjust the relaxation parameters — $s$ for the MRT scheme and $\tau$ for the BGK scheme. The systematic procedure to obtain macroscopic conservation equations 2 and 3 from the LB equation is the Chapman-Enskog analysis Chapman and Cowling (1970); Wolf-Gladrow (2000). For the BGK-LB scheme, the Chapman-Enskog expansion gives the kinematic viscosity on the lattice as Wolf-Gladrow (2000)

[TABLE]

The wave equation 1 has no dissipation term i.e. kinematic viscosity is zero. This is achieved by setting the collision relaxation time $\tau=1/2$. This particular choice of $\tau$ thus recovers macroscopic equations 2 and 3 describing the linear acoustic wave Chopard and Droz (1998).

For the MRT-LB scheme, the Chapman-Enskog analysis is carried out in moments space for each of the conserved and non-conserved moments. For the D2Q9 lattice, the analysis shows that the kinematic and bulk viscosities on the lattice may be set to zero with the choice $s_{p}=2$ and $s_{e}=2$ respectively Lallemand and Luo (2000). For the D2Q5 lattice, setting $s_{p}=2$ and either $s_{e}=2$ or the lattice sound speed $c_{s}^{2}=1/2$, recovers linear macroscopic equations 2 and 3 (see Appendix B). Other relaxation parameters (i.e. apart from $s_{p}$ , $s_{e}$) do not affect the hydrodynamics of the problem. However, in the case of the D2Q9 lattice, any choice except $s_{e}=2$ and $s_{\epsilon}=2$ causes numerical instability. The origin of these instabilities may be revealed by an analysis for the linear MRT-LB scheme similar to one performed in Dellar (2003).

In order to obtain the full Navier-Stokes equation from the LB equation, it is essential that microscopic velocity moments of lattice weights up to the fourth order are identical to that of the Maxwell-Boltzmann distribution with zero-mean velocity Wolf-Gladrow (2000); Frisch et al. (1986). However, the Chapman-Enskog expansion shows that, for the linear LB scheme that we use, velocity moments of lattice weights up to the second order are relevant for obtaining the macroscopic conservation Eqs. 2 and 3 Chopard and Luthi (1999). The constraints on lattice weights for the linear LB scheme are (all odd-order moments vanish)

[TABLE]

In the Chapman-Enskog expansion, the second constraint above determines proportionality between density and the pressure term in the momentum equation (see equations 43 and 47) and hence the local sound speed $c_{s}(\textbf{x})$. The Chapman-Enskog expansion also shows that the local lattice sound speed may be controlled by spatially adjusting the rest-particle lattice weight $w_{0}$ Chopard and Luthi (1999). This is achieved by setting

[TABLE]

where $\eta(\textbf{x})=c_{s}(\textbf{x})/c_{s\max}\leq 1$ , where $c_{s\max}$ is the maximum sound speed in the medium. In addition, because of lattice symmetry, not all lattice weights are different. For instance, in the D2Q5 lattice, we must have $w_{1}=w_{2}=w_{3}=w_{4}$. Also, lattice weights cannot be negative. Thus, given the value of $\eta(\textbf{x})$, Eqs. 17 and 18 determine the lattice parameters — weights $w_{i}(\textbf{x})$ and sound speed $c_{s}(\textbf{x})$.

Using 15, lattice weight constraints 17 and transformation matrix M (see Appendix A) for LB lattices, second- and higher-order velocity moments of the linear equilibrium distribution can be calculated. Thus, for the D2Q5 lattice, we get

[TABLE]

and

[TABLE]

Similarly, higher-order velocity moments of the linear equilibrium distribution on the D2Q9 lattice are calculated.

III Dispersion analysis

III.1 Grid Resolution (R)

In a numerical simulation, $\delta t\sim\delta x$ (acoustic scaling) and $\delta x\propto N^{-1/3}$ in 3D where N is the total number of grid-points. Hence, total number of time steps scale with $N^{1/3}$. The cost of a time step scales linearly with the total number of grid-points $N$. Hence, the cost of the numerical simulation scales as $\mathcal{O}\left(N^{4/3}\right)$ in 3D. Therefore, an important performance criterion of a numerical technique used to model wave propagation is the accuracy of the solution at low grid resolution. A natural way to define the grid resolution is to consider the number of grid points required to resolve the characteristic wavelength in the problem. Thus, the grid resolution is given by $R=(\lambda_{c}/\delta x)$, where $\lambda_{c}$ is the characteristic wavelength. The characteristic wavelength can, for instance, be the shortest wavelength in the problem.

III.2 LB dispersion relation

The dispersion relation captures the connection between the vector wave-number k and the frequency $\omega$ of the wave as captured by the numerical scheme. For waves travelling in homogeneous, dissipation-free media, $\omega$ and $k$ are proportional to each other in all directions i.e. there is no dispersion. However, a numerical scheme used to simulate these waves may introduce artificial dispersion and attenuation. We analyse the dispersion characteristics of the MRT-LB scheme for D2Q5 and D2Q9 lattices by studying the response of the scheme to plane waves in homogeneous media i.e $c_{s}(\textbf{x}),\rho_{0}(\textbf{x})=\textrm{constant}$.

We can express the equilibrium distributions $\textbf{g}^{eq}$ as a linear combination of the distribution functions g using 8, 9 and 15. This relation for the D2Q5 lattice is given by

[TABLE]

and for the D2Q9 lattice,

[TABLE]

Using these relations, equation 13 can be rewritten as

[TABLE]

where for the D2Q5 lattice

[TABLE]

Similarly, matrix A for the D2Q9 lattice can be constructed.

The distribution function is initialized such that macroscopic variables $p(\textbf{x},t)$ and $\textbf{v}(\textbf{x},t)$ evolve in time as a plane wave. This is achieved by initializing the $j^{th}$ component of the distribution function (in Fourier space) as

[TABLE]

Note that the coefficient of these distribution function plane waves is chosen such that we recover initial pressure and velocity plane waves consistent with equations 2 and 3 from equations 8 and 9. Substituting 24 in the Fourier transform of 23 gives us the dispersion relation for the MRT-LB scheme with a linear equilibrium distribution;

[TABLE]

where $\textbf{g}(\textbf{k},\omega)$ is the distribution function column vector in Fourier space.

The equation 25 is solved numerically to obtain the linear MRT-LB dispersion relation. The eigenvector corresponding to the propagating mode is, in general, a linear combination of hydrodynamic variables pressure and velocity in moments space. In Fig. 3, dispersion curves for the BGK- and MRT-LB schemes on D2Q5 and D2Q9 lattices are plotted for propagation along $\theta=0\degree\textup{and}~{}\theta=45\degree$ (see Fig. 2). The corresponding numerical attenuation (imaginary part of frequency $\omega$) is also plotted. The non-dimensional frequency (real and imaginary) $\omega^{*}=(\omega\delta x)/c_{s}\pi$ is plotted as a function of the non-dimensional wave-number $k^{*}=(k\delta x)/\pi$. Note that the non-dimensional wave-number $k^{*}$ also corresponds to grid resolution $R$ (Section III.1), with $k^{*}=1$ corresponding to the Nyquist limit of $2$ grid points per wavelength. Thus the working grid resolution for the LB scheme can be deduced from the dispersion curve by marking a point on the curve where it starts to deviate significantly from the exact curve. Lattice weights are determined from constraint equations 17 and 18 with the condition that $\eta(\textbf{x})=1$ corresponding to a homogeneous medium. For the D2Q5 lattice, we obtain $w_{0}=0,w_{1}=w_{2}=w_{3}=w_{4}=1/4$. For the D2Q9 lattice, 17 and 18 result in an under-determined set of equations for lattice weights. Thus, there is freedom to choose any combination of lattice weights for D2Q9 (but satisfying 17 and 18). For the curve in Fig. 3 (D2Q9-con), we choose conventional lattice weights for D2Q9: $w_{0}=4/9,w_{1}=w_{2}=w_{3}=w_{4}=1/9,w_{5}=w_{6}=w_{7}=w_{8}=1/36$.

Dispersion curves for the BGK- and MRT-LB scheme with $s_{p}=2$ and $s_{e}=2$ on the D2Q5 lattice are identical. Further, because the values of relaxation parameters for conserved moments, $s_{\rho}$ and $s_{v}$, are not relevant, the BGK- and MRT-LB D2Q5 schemes we use are identical. Additionally, these schemes are attenuation free, as desired for our problem. The D2Q5 dispersion is exact for propagation along $\theta=45\degree$. The BGK-LB scheme on the D2Q9 lattice is also attenuation free. However, the dispersion curve deviates significantly from the exact, especially for propagation along $\theta=45\degree$. The MRT-LB scheme on the D2Q9 lattice gives better dispersion characteristics at the cost of (small) attenuation. Also, these dispersion curves at best match D2Q5 counterparts. Thus the BGK-LB scheme on the D2Q5 lattice has the best dispersion characteristics.

Here we highlight the subtle distinction between free and forced waves which may be relevant for the dispersion analysis. Free waves have real wavenumber and complex frequency, whereas forced waves have complex wavenumber and real frequency. Hence free and forced waves respond differently to dissipation in the system and their dispersion characteristics are also affected Krüger et al. (2017). In our work, we are interested in comparing LB numerical schemes with standard finite-difference schemes. Hence, we are considering acoustic wave propagation in dissipation-free media for which the distinction between free and forced waves is not relevant Lele (1992).

III.3 Optimized lattice weights for the D2Q9 lattice

For the D2Q5 lattice, equations 17 and 18 yield a complete system of equations for lattice parameters. However, for the D2Q9 lattice, we obtain an under-determined system of equations. Since the D2Q9 lattice weights cannot be uniquely determined, we vary the lattice weights within these constraints. For each set of weights for the D2Q9 lattice, we calculate the numerical dispersion in equation 25. Performing this exercise, the lattice weights for the D2Q9 lattice which yield the best dispersion curves are obtained as $w_{0}=0,w_{1}=\varepsilon,w_{5}=1/4-\varepsilon$ where $\varepsilon$ is a very small positive number. Note that $\varepsilon=0$ causes instability, as the lattice sound speed $c_{s}$ becomes unity (see Appendix D). The optimum curves are plotted in Fig. 4 (D2Q9-opt) and compared with the D2Q5 as well D2Q9 curves with conventional weights. We see that dispersion curves of D2Q9-opt match best with the exact curves for propagation along $\theta=0\degree$. However, this particular choice of lattice weights makes D2Q9 roughly equivalent to D2Q5 rotated by $45\degree$. Note that the D2Q9 lattice weights $w_{0}=0,w_{1}=1/4,w_{5}=0$ gives the D2Q5 lattice and hence identical dispersion relations to the D2Q5 scheme.

III.4 Comparison with finite-difference schemes

The discussion hitherto suggests that BGK-LB with D2Q5 is the best LB scheme for the simulation of linear acoustic waves. We compare this scheme with standard schemes used in seismology — the second- and fourth-order finite-difference (FD) schemes. The stability of FD schemes as well as the LB scheme depends on the local Courant number which is given by Robertsson and Blanch (2011); Kupershtokh (2010)

[TABLE]

where $c_{s}=c_{s}(\textbf{x})$ is the local sound speed. From equation 26, we see that the Courant number is ratio between the physical information propagation speed and the numerical information propagation speed. For the BGK-LB D2Q5 scheme, the numerical information can propagate only along the axes with speed $\delta x/\delta t=1$ and hence it propagates with speed $1/\sqrt{2}\approx 0.71$ along the diagonal. Hence, the maximum Courant number $C_{max}\approx 0.71$ which is also the maximum lattice sound speed (see Appendix D). For the second-order and the fourth-order FD scheme used here (see Appendix C), the maximum allowed Courant number $C_{max}$ is approximately equal to 0.71 and 0.66 respectively. Considering these limits, we study the dispersion relation for different values of the Courant number.

In Fig. 5, the BGK-LB D2Q5 dispersion relation is compared with the second- and the fourth-order FD scheme at $C=0.66$ and $C=0.33$, for propagation along $\theta=0\degree$ and $\theta=45\degree$. These plots have normalised sound speed on the y-axis rather than normalised frequency. For these directions, the dispersion curves of the BGK-LB D2Q5 and the second-order FD scheme are identical (see Appendix E). As shown in Fig. 6, the BGK-LB D2Q5 and the second-order FD dispersion curves for the intermediate propagation angles are also identical, except at very low grid resolution where the BGK-LB D2Q5 is slightly better. The fourth-order FD scheme has better dispersion characteristics which further improve with decreasing Courant number. However, for propagation along $45\degree$ at high Courant numbers (close to $C_{max}=0.66$ for the fourth-order FD scheme), the BGK-LB D2Q5 (and the second-order FD) is more accurate.

IV Grid Refinement

Seismic waves travel through heterogeneous media where, background density and sound speed may vary significantly across the domain. The characteristic wavelength of seismic waves is directly proportional to the sound speed and the accuracy at low grid resolution (Section III.1) is an important performance criterion for numerical schemes modelling seismic waves. It is therefore desirable to keep the grid resolution, $R=(\lambda_{c}/\delta x)$, constant across the computational domain. If we model heterogeneous media as a patchwork of domains of homogeneous media, we will need a grid-refinement scheme for the LB simulation to maintain grid resolution across a varying sound speed.

The underlying requirement for the LB grid-refinement algorithm is to facilitate streaming of distribution functions across different patches of locally homogeneous media. With the grid-refinement scheme, lattice spacing $\delta x$ is set to vary depending on the sound speed. To maintain the microscopic velocities, corresponding changes in the time step $\delta t$ are required Dupuis and Chopard (2003); Lin and Lai (2000). The LB algorithm naturally splits into two steps — streaming and collision. So it is convenient to maintain a time-step ratio in multiples of $2$ between neighbouring patches of locally homogeneous media. Accordingly, the ratio of lattice spacing between two neighbouring patches of homogeneous media should also be in multiples of $2$.

Simulations in heterogeneous media were carried out using the LB grid-refinement scheme based on an algorithm suggested by Dupuis et al. Dupuis and Chopard (2003). In order to understand the method, let us consider a heterogeneous medium consisting of two regions of homogeneous media connected at the interface (Fig. 7a). The region on the left has greater sound speed and thus larger characteristic wavelength than the region on the right. To maintain the grid resolution, we need a coarse grid on the left and a fine grid on the right. Let the ratio of lattice spacing between the two domains be $n=(\delta x_{c}/\delta x_{f})$ where $c$ stands for coarse and $f$ stands for fine. Note that $n$ is a multiple of $2$. During the interchange between the coarse and fine domains, distribution functions should be appropriately scaled. The distribution function can be composed as a sum of equilibrium and non-equilibrium parts. Since the equilibrium distribution is a function of local macroscopic variables $p$ and v, it should remain unaltered during the interchange. The scaling of the non-equilibrium part is determined through Chapman-Enskog analysis Dupuis and Chopard (2003). Thus we have

[TABLE]

where $g_{i}^{eq,c}$ and $g_{i}^{eq,f}$ are the equilibrium distribution functions for the coarse and the fine domain respectively and

[TABLE]

where $\hat{g}_{i}$ are spatially and temporally interpolated distribution functions on the coarse-grid. Since $\tau=1/2$ (BGK-LB scheme), the kinematic viscosity 16 remains zero even if the time step $\delta t_{l}$ is varying across computational domain. As illustrated in Fig. 7b, we follow the following procedure for time evolution on a composite coarse fine grid:

• •

Initialise coarse (C) and fine (F) grid distribution functions to equilibrium values at $t=0$.

• •

Advance the coarse-grid distributions to $t=2$ with collision and streaming.

• •

If the wave is travelling from the coarse to the fine domain, calculate the fine-grid distribution functions at the interface at $t=1$ (thick line) using 28 with spatial and temporal interpolation of the coarse-grid distribution functions at $t=0$ and $t=2$. Else, fine-grid distribution functions are unchanged.

• •

Evolve fine-grid distribution functions from $t=1$ to $t=2$ through collision and streaming.

• •

If the wave is travelling from the coarse to the fine domain, correct fine-grid distribution functions at the interface according to 28 using coarse-grid distribution functions only at $t=2$. Else, if the wave is travelling from the fine to the coarse domain, correct coarse-grid distribution functions at $t=2$ according to 28 using fine-grid distribution functions at $t=2$.

• •

Evaluate and register macroscopic variables at $t=2$ on all grid nodes.

• •

Evolve the fine-grid distribution functions from $t=2$ to $t=3$.

The procedure is repeated for the desired total evolution time.

Attendant to simulating on non-uniform grids are spurious reflections that occur at the interface Bažant and Celep (1982). Indeed, the low order of accuracy of LB combined with the abrupt halving in grid spacing at the coarse-fine grid interface creates spurious reflections. Here we mitigate these spurious reflections by specifying the direction of wave propagation at the coarse-fine grid interface.

Another type of spurious reflection arises when distribution functions from the coarse region or the fine region hit the domain boundary at the interface. To counter it, we introduce an overlap region (see Fig. 7a) of few grid points where the coarse domain extends into the fine domain and vice versa. Streaming distributions coming from either domains, instead of stopping at the interface, proceed in the overlap region on respective grids. The overlap regions consist nodes which form sponge layers Colonius (2004). At a sponge node, during each collision, the distribution function is reduced by a fraction of its value. Thus, the sponge layers gradually absorb the wave as it propagates in the overlap region. This results in negligible reflections back into the physical domain. However, physical reflections are not affected as accurate boundary conditions are still maintained at the interface. Macroscopic variables are evaluated using distribution functions at legitimate coarse and fine nodes.

V Numerical experiments

V.1 Waves in homogeneous media

The findings of the dispersion-relation analysis are reinforced through simulation tests in homogeneous media. For test simulations, the source used to mimic the seismic disturbance is

[TABLE]

where

[TABLE]

and where central frequency $f_{c}=10$ Hz and wave velocity $v_{s}(\textbf{x})=4$ km/s. The source is located at a point $\textbf{x}_{s}$ in the domain. The temporal variation of the source is shown in Fig. 8. For homogeneous media, we set the mean density $\rho_{0}(\textbf{x})=1$ kg/m3. The domain size for the simulation is $60\lambda_{c}$, $\lambda_{c}$ being the central wavelength and simulations are carried out with grid resolution varying from $16$ grid points per $\lambda_{c}$ to $4$ grid points per $\lambda_{c}$. For the source term given by 29, wave equation 1 admits an analytical solution (in Fourier space)

[TABLE]

where $H_{n}^{(2)}$ is the Hankel function of the second kind. The solution in time domain is obtained by taking the inverse Fourier transform.

The wave is excited by the source 29 at a point $\textbf{x}_{s}$, typically at the center. The resulting waves propagate in all directions. Fig. 9 shows a snapshot of the propagating wavefront in the homogeneous medium. The LB simulations are carried out using the 2D lattices D2Q5 and D2Q9. In the BGK-LB D2Q5 and D2Q9 simulations, we use relaxation time $\tau=1/2$ (see section II.3). In the MRT-LB D2Q5 simulation we set only $s_{p}$ and $s_{e}$ equal to 2 and values of the other relaxation parameters are not relevant. In the MRT-LB D2Q9 simulation, $s_{\epsilon}$ is also set to 2 and $s_{q}$ is set to 1. For the D2Q5 lattice, simulations with the MRT-LB scheme and the BGK-LB scheme are identical. For both simulations, we see a distinct wavefront with a compact shape at $\theta=45\degree$. For the BGK-LB D2Q9 simulation, we see a trail of waves following the main wavefront as the lattice responds differently to waves of different spatial wavenumbers. With the MRT-LB D2Q9 simulation, this trail of waves is eliminated as dispersion performance improves (see Fig. 3). However, the main wavefront in D2Q5 simulation is more isotropic and less dispersed than in D2Q9 simulation.

These results are compared to the exact analytical solution 31 by studying the time response at different points in the medium. In Fig. 10, the time response curve at two points in the medium obtained using the BGK- and MRT-LB scheme on D2Q5 and D2Q9 lattices is compared with the exact result. The overall results for D2Q5 are better than for D2Q9, as expected from the dispersion analysis.

These results suggest that the BGK-LB D2Q5 is the best choice for simulating linear acoustic waves. A detailed look at the time response pulse and a comparison with second- and fourth-order finite-difference schemes for propagation along $\theta=15\degree$ and $\theta=30\degree$ is presented in Fig. 11. We see that BGK-LB D2Q5 and the second-order finite-difference scheme results are comparable, consistent with the dispersion curves in Fig. 6. Overall however, the fourth-order finite-difference scheme is more accurate than the LB scheme. The simulations at low grid resolution show similar trends.

V.2 Waves in heterogeneous media

The simulations in a heterogeneous medium are carried out using the same source 29. The heterogeneous medium we consider comprises two homogeneous media with different sound speeds, joined together to form an interface. The ratio of sound speeds in the two domains, i.e velocity contrast, is 0.8. Accordingly, the grid on the left is coarse and the grid on the right is fine. We implement the grid-refinement algorithm discussed in Section IV. The results are shown in Fig. 12.

Fig. 12a shows simulations in a homogeneous medium with grid refinement. The source is located on the left i.e. coarse-grid. The BGK-LB scheme on the D2Q5 lattice is used. Waves propagate smoothly across the interface without any artificial reflection. The only difference between the left and the right is the resolution, which is $8$ points per wavelength on the coarse side and $12.8$ points per wavelength on the fine side. Note that the simulated wavefront is broader in the coarse domain compared to Fig. 9a. This is because resolution in the coarse domain is lowered by a factor of two compared to the resolution in the computational domain in Fig. 9a. Fig. 12b shows the simulation in the heterogeneous medium. Waves from the coarse domain suffer a reflection when they cross the medium boundary. On the fine side, waves propagate with a reduced speed.

VI Conclusions

We have demonstrated a successful formulation of linear acoustic wave propagation using the BGK- and MRT-LB frameworks with a linear equilibrium distribution function, first proposed in Chopard and Luthi (1999). Similar to analysis performed in Lallemand and Luo (2000) for the conventional LB scheme, we develop a formalism to calculate the dispersion relations for the linear BGK- and MRT-LB scheme. With our formalism, it is possible to compare the performance of various LB lattices for simulation of linear acoustic waves. Our formalism is also useful in comparing the dispersion relation of the LB schemes with standard finite-difference schemes for a given Courant number. The LB dispersion relations as well as numerical simulations are in reasonable agreement with theoretical results.

Our dispersion analysis establishes that the fourth-order finite-difference scheme is better than any LB numerical scheme for the simulation of the linear acoustic waves. As with the finite-difference schemes Lele (1992), the numerical dispersion of the LB scheme is also anisotropic. Our dispersion analysis shows that for the linear LB scheme, the numerical dispersion is most inaccurate if the wave is propagating along one of the directions of LB particle streaming. Hence, the dispersion relations for the D2Q9 lattice with the BGK and the MRT schemes are worse than for the D2Q5 lattice. In general, the dispersion performance of the LB scheme is degraded as we increase the number of streaming directions on the lattice. For the D2Q5 lattice however, the BGK and the MRT dispersion performance is identical even with only two relaxation parameters $s_{p}$ and $s_{e}$ set equal to 2 in the MRT scheme, a requirement necessary to recover macroscopic equations 2 and 3. Thus, the BGK-LB scheme on the D2Q5 lattice is the best suited LB numerical model for seismic waves. Also, the BGK-LB D2Q5 scheme — which is a second-order scheme — is comparable with the second-order finite-difference scheme.

To keep the computational cost of simulation low, numerical scheme used should be sufficiently accurate even with low grid resolution i.e. the number of grid points per wavelength. For heterogeneous media, the wavelength in the domain changes as a function of the local sound speed. Depending on the change in wavelength, the grid size needs to be altered to maintain the number of grid points that resolve the wave i.e. the grid resolution. This requires a grid-refinement algorithm for the LB scheme to smoothly accommodate the change in the grid size across the computational domain. With the grid-refinement algorithm presented here, it is possible to maintain a uniform grid resolution and accuracy of the simulation over the entire computational domain. We have successfully implemented the algorithm to simulate waves in heterogeneous media with velocity contrasts using the BGK-LB scheme on D2Q5 lattice. The techniques used here can be easily extended to 3D wave propagation problems in much more complicated environments. Overall, the LB scheme is therefore a promising tool for faster, cost-effective simulation of waves in seismology.

Acknowledgement

We would like to thank TIFR for their support and the use of the SEISMO computer cluster which was used to perform the calculations presented here. SMH acknowledges support from Ramanujan fellowship SB/S2/RJN-73 2013, the Max-Planck Partner Group Program and the NYUAD Center for Space Science. We would also like to thank L-S Luo and J. O. Blanch for useful conversations. We thank the anonymous reviewers for their valuable comments and suggestions.

Appendix A Transformation matrix M

For the D2Q5 lattice,

[TABLE]

and for the D2Q9 lattice,

[TABLE]

Appendix B Chapman-Enskog analysis for the linear MRT-LB scheme on D2Q5 lattice

In the Chapman-Enskog multiscale expansion procedure Chapman and Cowling (1970); Wolf-Gladrow (2000), the distribution function and the time and spatial derivatives are expanded in terms of the small expansion parameter $\epsilon$ (Knudsen number) He and Luo .

[TABLE]

[TABLE]

[TABLE]

Using these expansions in 14 (without the scalar source term) and retaining terms only upto $\mathcal{O}(\epsilon^{2})$ we obtain following equations (order by order in $\epsilon$):

[TABLE]

[TABLE]

[TABLE]

Transforming these equations into the moment space gives

[TABLE]

[TABLE]

[TABLE]

where $\textbf{{E}}=\textup{M}~{}\textup{diag}(\textbf{c}_{0},\textbf{c}_{1},\textbf{c}_{2},...)~{}\textup{M}^{-1}$(diag refers to diagonal matrix) and I is the Identity matrix. Equation 40 can be simplified using 39 as

[TABLE]

In the case of the D2Q5 lattice we have

[TABLE]

also

[TABLE]

[TABLE]

Using 15, components of the equilibrium moment vector are obtained as

[TABLE]

For the conserved moments — fluctuating density and momentum — $\rho^{(k)},\rho_{0}v_{x}^{(k)},\rho_{0}v_{x}^{(k)}=0$ for $k>0$. Writing $\rho=\rho(\textbf{x},t)$, $\textbf{v}=\textbf{v}(\textbf{x},t)$, $\rho_{0}=\rho_{0}(\textbf{x})$ and $c_{s}=c_{s}(\textbf{x})$ and substituting these values in 39, we obtain at the $\mathcal{O}(\epsilon)$ following set of equations for each component of m,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Similarly using 41, we obtain at $\mathcal{O}(\epsilon^{2})$, for the conserved components of m (using $e^{(1)}$ and $p_{xx}^{(1)}$ from 44 and 45 respectively)

[TABLE]

[TABLE]

42 and 46 add to give the continuity equation 2 (without the scalar source term). Since the D2Q5 lattice cannot satisfy isotropy of the fourth order tensor Frisch et al. (1986); Wolf-Gladrow (2000), the exact recovery of the Navier-Stokes equation is not possible from 43 and 47. However, by setting $s_{p}=2$ and either $s_{e}=2$ or $c_{s}^{2}=1/2$, the linear conservation of momentum 3, is recovered.

Appendix C Finite-difference schemes used for comparison

The second-order accurate (O[2,2]) finite-difference solution of the linear wave equation is obtained by using centered differences for the spatial and temporal dervatives,

[TABLE]

The fourth-order accurate (O[2,4]) solution is obtained by solving coupled Eqs. 2 and 3 simultaneously. These equations are of the general form

[TABLE]

We use a 4-point stencil for spatial derivatives as specified below.

[TABLE]

[TABLE]

Note that we can also derive equation 48 from equation 49 by taking second-order centered differences and eliminating $v$. We stop the simulation well before the waves reach boundary of the computational domain.

Appendix D The Courant number for the LB scheme

The local Courant number given by 26 depends on the local sound speed $c_{s}$. For the LB scheme, the dependence of $C$ on lattice parameters can be deduced from 17. In particular, for the D2Q5 lattice, the relation between $C$ and the rest particle weight $w_{0}$ is given by

[TABLE]

The rest particle weight $w_{0}$ is modulated in order to obtain the locally varying lattice sound speed for heterogeneous media. However, it is essential to choose $w_{0}$ such that $C$ is always real and less than unity which is always true for the D2Q5 lattice. For homogeneous media, it is convenient to set $w_{0}=0$ which gives maximum allowed Courant number $C_{max}=1/\sqrt{2}$ for the LB D2Q5 scheme. In case of D2Q9, there is no definite relation between $C$ and $w_{0}$. Hence, $w_{0}$ needs to be chosen carefully.

Appendix E The analytical 1-D dispersion relation for the MRT-LB D2Q5 scheme

The solution of the eigenvalue equation 25 gives the dispersion relation for the MRT-LB scheme. For the D2Q5 lattice, the characteristic equation is a fifth-order polynomial and hence writing down the general analytical solution is not possible. However, one can obtain the analytical solution for the propagation along the x-axis ($\theta=0\degree$) by considering projection of the D2Q5 lattice on the x-axis i.e. the D1Q3 lattice (see Fig. 13).

In the D1Q3 lattice, the lattice velocities are $c_{0}=(0,0)$ and $c_{1,2}=(\pm 1,0)$. Using the constraints on the lattice weights (equation 17) and the lattice symmetry we have, $w_{1}=w_{2}=w$ and $w_{0}=1-2w$. Also the lattice sound speed $c_{s}=2w$ (homogeneous media). The matrix

[TABLE]

transforms distribution functions to the moments space variables $\textbf{m}=\{\rho,\rho_{0}v_{x},e\}$. For the D1Q3 lattice, the matrix A in equation 23 is

[TABLE]

Substituting this in the eigenvalue equation 25 yields the characteristic equation

[TABLE]

where $\lambda=\exp(-i\omega)$. Clearly, one root is $\lambda=-1$ and hence $\omega=\pi$ which corresponds to a non-propagating mode. The other roots are

[TABLE]

With the choice $\left(1-2c_{s}^{2}\sin^{2}(k/2)\right)=\cos{\theta}$, $\lambda=\exp{(\pm i\theta)}$ and eigenvalues for the propagating modes

[TABLE]

This is identical to the 1-D dispersion relation for the second-order Fd scheme in 48.