## Abstract

Acoustic wave propagation in lattice Boltzmann Bhatnagar–Gross–Krook simulations may be analysed using a linearization method. This method has been used in the past to study the propagation of waves that are viscously damped in time, and is here extended to also study waves that are viscously damped in space. Its validity is verified against simulations, and the results are compared with theoretical expressions. It is found in the infinite resolution limit *k*→0 that the absorption coefficients and phase differences between density and velocity waves match theoretical expressions for small values of *ωτ*_{ν}, the characteristic number for viscous acoustic damping. However, the phase velocities and amplitude ratios between the waves increase incorrectly with (*ωτ*_{ν})^{2}, and agree with theory only in the inviscid limit *k*→0, *ωτ*_{ν}→0. The actual behaviour of simulated plane waves in the infinite resolution limit is quantified.

## 1. Introduction

In a seminal paper, Lighthill showed how flow fields may act as acoustic sources [1], thereby planting the seed for the scientific field now known as aeroacoustics. The related field of computational aeroacoustics (CAA) has been studied since the 1980s, and two main approaches of CAA have been identified.

**The hybrid approach,**where the flow and acoustic fields are found separately. The flow field is analysed to find the acoustic source strength of the flow, which is then used to compute the acoustic field in the surrounding domain.**The direct approach,**where the acoustic field comes out as a natural part of a numerical solution to the compressible Navier–Stokes equation.

One important weakness of the hybrid approach is its inability to simulate the feedback of the acoustic field on the flow field. This means that it is not usable in cases where this feedback is critical to the process of sound generation, such as in the singing risers problem currently studied in the natural gas industry [2]. It is also difficult to capture complex geometries with this method.

The direct approach has no such problems, but traditional numerical simulations of the compressible Navier–Stokes equation are far more complex and demanding than hybrid methods. It is therefore our hope that the lattice Boltzmann method (LBM) can be used for direct CAA. Compared with traditional methods, the LBM is straightforward to implement. It is also an explicit method, which can be parallelized with a nearly linear speedup in number of processors.

There exists a large body of work which validates the LBM as a useful tool for simulations of incompressible flow [3]. However, it has been shown mathematically that the LBM may also be used for solving the compressible Navier–Stokes equation [4], albeit with spurious terms that limit its applicability to small Mach numbers and cause incomplete Galilean invariance.

To ascertain whether the LBM may become a useful tool for CAA, we need a better picture of the acoustic behaviour of the method. This paper will focus on viscous damping of acoustic waves. In §2, the theory of this behaviour is described. Section 3 describes a method of linearization analysis that can be applied to study this behaviour in the LBM. The method is validated in §4, and is then used to compare the behaviour of LB-simulated waves with theoretically expected behaviour.

## 2. Theoretical background

When performing the Chapman–Enskog analysis on the LBM, certain assumptions must be made to retrieve the compressible Navier–Stokes equation. First, one assumes that the Mach number is sufficiently low that the aforementioned spurious terms can be neglected. Second, one assumes isothermal compression, so that the absolute pressure and density *p* and *ρ* are linked through the thermodynamic speed of sound as *p*=*c*^{2}_{s}*ρ*. These assumptions fit well with those made in linear acoustics: low Mach number and adiabatic compression, so that the off-equilibrium pressure and density *p*′ and *ρ*′ are linked as *p*′=*c*^{2}_{s}*ρ*′.

One result of the Chapman–Enskog analysis is that the kinematic shear and bulk viscosities for the Bhatnagar–Gross–Krook (BGK) collision operator are given through the relaxation time *τ* as, respectively,
2.1a
and
2.1b
This constant ratio between bulk and shear viscosities makes it impossible to perform simulations of fluids with correct transport coefficients without resorting either to extended LBMs (e.g. [4]), or to another type of collision operator such as the multiple relaxation time (MRT) operator (e.g. [5]). Both these choices would allow us to set bulk and shear viscosity independently.

### (a) Solutions to the lossy wave equation

Conservation of mass and momentum is expressed through the continuity and compressible Navier–Stokes equations, which relate density *ρ*, pressure *p*, and particle velocity **u**,
2.2a
and
2.2b
Assuming low Mach number, adiabatic compression, and sufficient distance to boundaries, it may be shown similarly to [6] that this implies a lossy wave equation,
2.3
where is the viscous relaxation time. Inserting from equation (2.1), we find for the LBM that *τ*_{ν}=2*τ*−1. From the extended analysis of Hamilton & Morfey [7], we see that equation (2.3) is correct up to second order in Mach number and *ωτ*_{ν} if the effects of thermal conduction and molecular relaxation are neglected.

In one dimension, this equation implies general damped wave solutions in pressure, density and particle velocity,
2.4a
2.4b
and
2.4c
where a hat implies a complex quantity, and and are complex angular frequency and complex wavenumber, respectively. *α*_{t} and *α*_{x} are the temporal and spatial absorption coefficients, while the real parts *ω*′ and *k*′ are linked through a phase velocity, *c*_{p}=*ω*′/*k*′. By comparison, we have in lossless media that , and *c*_{p}=*ω*/*k*=*c*_{s}.

Two important special cases are temporal damping (when ) and spatial damping (when ). In both cases, properties of the waves can be found by inserting equation (2.4a) into equation (2.3). For spatially damped waves, we thus find
2.5
whereas for temporally damped waves we find
2.6
These expressions only depend on the dimensionless quantity *ωτ*_{ν}, which may be seen as a characteristic dimensionless number for the effect of viscous acoustic damping. Spatially damped waves occur whenever a wave is radiated from an acoustic source in an absorbing medium, and is the type commonly encountered in nature.

Early numerical studies of wave propagation with the LBM focused on temporally damped waves, propagating a single wavelength in a periodic system [8]. Later numerical studies have taken spatial damping into account as well [9].

### (b) Linearization of conservation equations

By inserting the damped wave solutions in equation (2.4) into the conservation equations in equation (2.2) under the linearizing assumptions that the wave amplitudes , and are very small, we may reduce these equations to relations between the complex wave amplitudes, 2.7a and 2.7b The absolute values of these ratios indicate the amplitude ratios between the waves, while the complex phases describe the phase difference between them. For temporal and spatial damping, we may use expressions from equations (2.5) and (2.6).

## 3. Linearization analysis of lattice Boltzmann

The acoustic behaviour of LB may be analysed by performing a similar linearization [4,5]. We assume that we have a rest state of density *ρ*_{0} with a fluctuation on top. We assume this to be on a general complex one-dimensional wave form,
3.1
If this wave is very small, we may linearize the equilibrium distribution,
3.2
giving
3.3
In this small-amplitude limit, the aforementioned errors disappear.

While this behaviour could be examined for a number of different lattices, we will follow [10] and limit ourselves to the D1Q3 lattice, which is a one-dimensional projection of a number of common lattices (D2Q9, D3Q15, D3Q19 and D3Q27). Thus, plane waves in D1Q3 are equivalent with plane waves along any main axis of these other lattices. The D1Q3 lattice is given by velocities [*c*_{−},*c*_{0},*c*_{+}]=[−1,0,1] and weights [*t*_{−},*t*_{0},*t*_{+}]=[1/6,2/3,1/6].

With density and momentum given as
3.4a
and
3.4b
we find a post-collision distribution
3.5
If we for simplicity assume *x* and *t* such that , we find
3.6
which results in the eigenvalue problem as found in [4],
3.7

Through assuming known and finding the eigenvalue numerically, we may study the behaviour of temporally damped waves in the LBM. Alternatively, we may assume known and find the value of that gives the correct eigenvalue e^{iω} through a nested binary search in *k*′ and *α*_{x}. This allows us to study the behaviour of spatially damped waves also.

The eigenvector itself may be used to determine the relative complex amplitudes and of the density and velocity waves. It may be found from equations (2.4) and (3.4) that and , so that 3.8 In this way, we may compare the ratios of the LB amplitudes with theoretical ones from equation (2.7a).

## 4. Results

We may validate the linearization method by comparing with actual LB simulations. A spatially damped plane wave may be created in a simple fashion in a D1Q3 lattice by pinning density and velocity to small sinusoidal functions around an equilibrium, setting , in the first node and always relaxing it to equilibrium. We may then measure the phase velocities and absorption coefficients of the resulting waves.

The comparison is performed in figure 1, where we also compare with theoretical values. We see that the linearization analysis for spatially damped waves accurately predicts the values measured from simulations, validating this method as a tool for predicting the behaviour of spatially damped waves. Validation for temporally damped waves has already been performed [4]. From this point, we will therefore use the linearization analysis to study the behaviour of LB waves.

The theoretical values in equations (2.5) and (2.6) change only with the dimensionless parameter *ωτ*_{ν}. Our analysis shows that the LB phase velocity has a small numerical dispersion which for low *k* is linear in *k*^{2}. This *k*^{2} error also affects the spatial absorption coefficient, while the temporal absorption coefficient is hardly influenced in this way, being instead mainly affected by a very small *k*^{4} error.

In addition, the values for the phase velocity and absorption coefficients do not agree even when the numerical resolution goes to infinity (the *k*→0 limit), which means that the LBM with the BGK operator does not appear to be a consistent method for simulating waves. We will now examine how the waves behave in this limit.

Figure 2 compares the behaviour of LB waves in this limit with theoretical behaviour. We find an agreement with theory and Chapman–Enskog analysis in the *k*→0, *ωτ*_{ν}→0 limit for both phase velocity and absorption. The phase velocity is a linear function of (*ωτ*_{ν})^{2} as expected, but the slope is too steep in both cases. The absorption is correct apart from some small errors in (*ωτ*_{ν})^{3}.

Finally, we compare the ratio , as given from theory (equation (2.7a)) and as predicted through the linearization analysis (equation (3.8)). From figure 3, we see that the amplitudes of the waves are linear functions of (*ωτ*_{ν})^{2}, but that the slopes of spatially and temporally damped waves are again both steeper than they should be. The phase shift is largely correct, with a small (*ωτ*_{ν})^{3} error for spatially damped waves and an insignificant (*ωτ*_{ν})^{5} error for temporally damped waves. Again, the results agree with theory in the *k*→0, *ωτ*_{ν}→0 limit.

We have seen that the waves’ properties are quite well-behaved, although not necessarily correct. Table 1 compares the LB and theoretical properties of the waves we have studied in the *k*→0 limit, with terms of order or higher neglected. We find that while the normalized absorption coefficients and phase shifts are consistent with theory in this limit, the phase velocities and amplitude ratios do not match except in the inviscid limit. In fact, both of these increase at rates which are (*ωτ*_{ν})^{2}/4 greater than they should be.

## 5. Conclusion

We have examined the behaviour of viscously damped LB waves using the BGK collision operator. The linearization analysis was shown to be capable of analysing not only temporally damped waves, but also spatially damped waves.

We found that while the absorption coefficients and phase shifts match theory well, the phase velocities and amplitude ratios do not, except in the limit of an inviscid lossless fluid. The results are valid for plane waves propagating along a main axis in any lattice for which D1Q3 is a projection. The BGK operator allows no freedom to affect this, so addressing this problem would require an alternative method. On the other hand, the actual behaviour in the infinite resolution limit *k*→0 has been quantified in table 1, and the regularity of the deviations makes it quite possible to predict LB phase velocities and absorption coefficients *a priori* [11].

Since there is no independent pressure wave in the LBM, it is not possible to perform comparisons with equation (2.7b). In fact, the Chapman–Enskog analysis for the LBM assumes that *p*=*c*^{2}_{s}*ρ*, which may be shown from equation (2.7) to be incorrect for damped acoustic waves. That this incorrect assumption underlies the LBM may be a reason for the systematic deviations from theory which we have seen here.

However, an extended energy-conserving LBM has been proposed [12], and has been validated for a number of compressible flow problems [13]. Perhaps this method would give results more consistent with theory.

The deviations are still quite small, and have been found to be comparable to the ones found for high-order finite-difference schemes [14]. The LBM has also been found to be faster than high-order methods of similar accuracy [14]. This indicates that the LBM is a promising tool for CAA.

In this paper, we have only looked at the effects of viscosity on damped acoustic waves. In real fluids, there are very significant contributions to this damping from thermal conduction and molecular relaxation [6,15]. It remains to be seen how the LBM may be adapted to handle all three damping mechanisms.

## Footnotes

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

- This journal is © 2011 The Royal Society