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

We analyze 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 as well as the general multirelaxation-time 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.

Figure 1

Lattices used in simulation: (a) D2Q5 lattice; (b) D2Q9 lattice.

Figure 2

Direction of propagation for plane waves on the lattice. The numerical dispersion relation strongly depends on the angle $\theta$ along which the waves are propagating.

Figure 3

The BGK- and MRT-LB numerical dispersion and attenuation curves for D2Q5 and D2Q9 lattices. The real or imaginary component of the nondimensional frequency $\left(\omega \delta x\right)/c_s\pi$, is plotted against the nondimensional wave number $\left(k\delta x\right)/\pi$. Panels (a) and (c) show numerical dispersion for propagation along $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$ respectively and panels (b) and (d) show numerical attenuation for propagation along $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$ respectively. Lattice weights used for D2Q5 are $w_0=0,w_1=1/4$. For D2Q9, conventional lattice weights $w_0=4/9,w_1=1/9,w_5=1/36$ are used (curves marked D2Q9-con). For the MRT implementation, relaxation parameters $s_p$ and $s_e$ are set to 2 for both lattices. For D2Q9, relaxation parameter $s_{\epsilon }$ is also set to 2 and $s_q$ is set to 1. Relaxation parameters for the conserved moments $s_{\rho }$ and $s_v$ can be set to any value as it does not affect the hydrodynamics. The D2Q5 dispersion curves for the BGK- and MRT-LB schemes are identical (better than D2Q9). The legend in the first plot also applies to the rest of the plots.

Figure 4

The optimum dispersion curves for the D2Q9 lattice. Here, the real component of the nondimensional frequency, $\left(\omega \delta x\right)/c_s\pi$, is plotted against the nondimensional wave number $\left(k\delta x\right)/\pi$ for propagation along $\theta =0^{\circ }$ shown in panel (a) and $\theta ={45}^{\circ }$ shown in panel (b). Lattice weights used for D2Q5 are $w_0=0,w_1=1/4$. However, lattice weights required $w_0=0,w_1=0.001,w_5=0.249$ make the D2Q9 lattice equivalent to the D2Q5 lattice rotated by ${45}^{\circ }$.

Figure 5

Comparison of the BGK-LB D2Q5 dispersion curves (normalized sound speed $\left(\omega \delta x\right)/k{c}_s\pi$ vs normalized wave number $\left(k\delta x\right)/\pi$) with second- and fourth-order FD schemes. Panels (a) and (b) show numerical dispersion at Courant number $C=0.66$ (maximum allowed value for the fourth-order FD scheme) for propagation along $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$ respectively. Panels (c) and (d) show numerical dispersion at Courant number $C=0.33$ for propagation along $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$ respectively. The BGK-LB D2Q5 curves are identical to the second-order FD curves for these propagation angles. The fourth-order FD scheme has better dispersion characteristics overall. The legend of the first plot also applies to the rest of the plots.

Figure 6

Comparison of the BGK-LB D2Q5 dispersion curves [normalized sound speed $\left(\omega \delta x\right)/k{c}_s\pi$ vs normalized wave number $\left(k\delta x\right)/\pi$] with exact and second-order FD scheme at Courant number $C=0.71$ (maximum allowed value for both the BGK-LB D2Q5 and the second-order FD scheme), for propagation along $\theta ={15}^{\circ }$ shown in panel (a) and $\theta ={30}^{\circ }$ shown in panel (b). These BGK-LB D2Q5 dispersion curves are slightly better than second-order FD. The legend of the first plot also applies to the second plot.

Figure 7

Mesh refinement scheme: (a) coarse-fine interface; (b) LBM sequence on composite mesh—interpolation sequence is based on scheme proposed in [35].

Figure 8

Time variation of model source for acoustic simulations.

Figure 9

Simulations of waves in homogeneous medium. Panels (a) and (b) show simulation using the BGK-LB scheme on D2Q5 and D2Q9 lattice respectively. Panels (c) and (d) show simulation using the MRT-LB scheme on D2Q5 and D2Q9 lattice respectively. The D2Q5 lattice with weights $w_0=0$ and $w_1=1/4$ and the D2Q9 lattice with weights $w_0=0, w_1=0.01$, and $w_5=0.24$ are used. For the MRT-LB simulation $s_p,s_e=2$ on both lattices. For the D2Q9 lattice, $s_{\epsilon }$ is also set to 2 and $s_q$ is set to 1. The grid resolution is 16 points per wavelength. The source is located at the center. Since the BGK-LB D2Q9 dispersion relation deviates significantly from exact, different wave numbers propagate at different speeds, resulting in a series of trailing waves behind main wave front. These trailing waves are eliminated in the MRT-LB simulation due to attenuation. Both the BGK-LB and the MRT-LB simulations on the D2Q5 lattice are identical and better than D2Q9.

Figure 10

The time responses of the medium at distance $r=21\lambda$ from source. Panels (a) and (b) show the time response from BGK-LB simulation of wave propagation along $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$ respectively. Panels (c) and (d) show the time response from MRT-LB simulation of wave propagation along $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$ respectively. Numerical results are compared with the analytical result. Both BGK- and MRT-LB D2Q5 simulations are identical and the schemes capture the response pulse in close approximation to the analytical result. For the BGK-LB D2Q9 simulation, the primary pulse arrives approximately at the same time as the analytical result. However, there are a series of trailing response pulses after the main pulse. The MRT-LB D2Q9 simulation removes the trailing waves due to attenuation. But overall the performance is less accurate than D2Q5.

Figure 11

Comparison of the BGK-LB D2Q5 scheme with second- and fourth-order finite-difference schemes. Detailed view of time responses at distance $r=21\lambda$ from source for wave propagation along $\theta ={15}^{\circ }$ shown in panel (a) and $\theta ={30}^{\circ }$ shown in panel (b).

Figure 12

A snapshot of the propagating wave front using grid refinement with the BGK-LB scheme on the D2Q5 lattice (a) in a homogeneous (uniform) medium and (b) in a heterogeneous (nonuniform) medium. For both the media, the left side is the coarse computational domain and the right side is the fine computational domain. A line separating the two domains is also shown. For the homogeneous medium, we see no artificial reflection at the coarse-fine interface. For the heterogeneous medium, we see a reflection at the interface because of the sound speed contrast (0.8) between the coarse and fine domains.

Figure 13

The D1Q3 lattice: Projection of the D2Q5 lattice in one dimension.