{'fraction_non_alphanumeric': 0.04465164222461697, 'fraction_numerical': 0.021740148994986646, 'mean_word_length': 4.2727272727272725, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quartic Parameters for Acoustic Applications of Lattice Boltzmann Scheme

Using the Taylor expansion method, we show that it is possible to improve the lattice Boltzmann method for acoustic applications. We derive a formal expansion of the eigenvalues of the discrete approximation and fit the parameters of the scheme to enforce fourth order accuracy. The corresponding discrete equations are solved with the help of symbolic manipulation. The solutions are explicited in the case of D3Q27 lattice Boltzmann scheme. Various numerical tests support the coherence of this approach.


1) Introduction
• The lattice Boltzmann equation methodology is a general framework for approximating problems that are modelled by partial differential equations under conservative form arising in physics. It has been first proposed for fluid dynamics in the context of cellular automata (see e.g. Hardy et al [17], Wolfram [36], d'Humières et al [6]). The classical lattice Boltzmann scheme (McNamara-Zanetti [27], Higuera et al [20], Qian et al [29], d'Humières [5]) has been first developed for nonlinear fluid problems. It has also the capability to approximate thermal flows. (Chen et al [4]), magnetohydrodynamics (Dellar [7]) or coupling between fluid and structures (see e.g. Kwon [22]).
• We have proposed in [8,9] the Taylor expansion method to analyze formally the d'Humières lattice Boltzmann scheme [5] when the mesh size tends to zero. We have then replaced the (much more) formal Chapman-Enskog expansion methodology used in [5] by a simple Taylor expansion relative to the grid spacing. In this way we can obtain the modified equations of the scheme (in the sense of Lerat-Peyret [25] and Warming-Hyett [34]; see also [16,2,33]) at an arbitrary order in the general nonlinear case. We made more explicit in [10] the result at third order accuracy in all generality. We have observed at this occasion that a serious difficulty with the lattice Boltzmann scheme lies in the fact that the equivalent mass conservation equation contains an a priori non-null third order term. We have also proposed an algorithm to derive the modified equation in the case of linearized problems [11]. We have applied this methodology to derive the so-called "quartic parameters" to enhance the accuracy of the lattice Boltzmann scheme to simulate shear waves [11]. An application of this result is used by Leriche et al [26] for computing with spectral method and lattice Boltzmann scheme the eigenmodes of the Stokes problem in a cube. This Taylor expansion method can also be used for a detailed analysis of boundary conditions. In [12] we have shown that the "magic boundary parameters" of Ginzburg and Adler [14] depend in fact upon the detailed choice of the parameters of the lattice Boltzmann scheme.
• In this contribution, we consider linearized athermal acoustics in two and three space dimensions. In Section 2, we recall the essential properties of the d'Humières scheme. Then in Section 3, we use the method of formal expansion to expand the discrete eigenmodes of the acoustic system of partial differential equations as the mesh size tends to zero. In Section 4, we increase the accuracy of the lattice Boltzmann scheme with the development of "quartic" parameters and develop also a weaker approach by enforcing isotropy at fourth order accuracy. Applications in two and three space dimensions (D2Q13 and D3Q27 schemes) are presented in Sections 5 and 6 respectively. We detail our version of the D3Q27 scheme and the explicit formulae for the determination of the quartic parameters.

2) d'Humières lattice Boltzmann scheme
• In the framework proposed by d'Humières [5], the lattice Boltzmann scheme uses a regular lattice L with typical mesh size ∆x and for each node x in L, a discrete set of (J + 1) velocities V is given. In this contribution, the set V does not depend on the vertex x. We introduce a velocity scale λ and the time step ∆t is obtained according to the so-called "acoustic scaling" : For x ∈ L and v j ∈ V, the point x + v j ∆t is supposed to be a new vertex of the lattice. The unknown of this lattice Boltzmann scheme is the particle distribution f j (x, t) for x ∈ L, 0 ≤ j ≤ J and discrete values of time t. The numerical scheme proceeds into two major steps.
• The first step describes the relaxation f −→ f * of particle distribution f towards a locally defined equilibrium. It is local in space and nonlinear in general. In this paper devoted to acoustics we consider only linear contributions. It is defined with the help of a fixed invertible matrix M. The moments m k are defined through a simple linear relation Note that this very simple linear hypothesis is not satisfied in the scheme proposed by Geier et al [13]. The first d + 1 moments (where d is the space dimension, equal to d = 2 or d = 3 in the present applications) are the total density ρ and the momentum q α : We denote by W the vector (in R N ) composed of the density and the components of the momentum. These moments are supposed to be at equilibrium: The relaxation evolution m −→ m * due to linearized collisions is local and trivial for the vector W of conserved quantities, as remarked in (4). For the other moments, a distribution of equilibrium moments G(•) is given as a (linear for the present study of acoustic waves) function of the vector W : For k ≥ N, the evolution m −→ m * is supposed to be uncoupled: It is parametrized by the so-called "relaxation rate parameters" s k . For stability reasons (see e.g. Lallemand and Luo [23]), we have the inequalities When the components m * k (x, t) are computed for each x ∈ L at discrete time t, the distribution f * j (x, t) after relaxation is reconstructed by inversion of relation (2): • The second step is just a free advective evolution without collision through characteristics: • The asymptotic analysis of cellular automata provides evidence supporting asymptotic partial differential equations with transport coefficients related to the induced parameter defined by the so-called Hénon's relation [18] (10) The lattice Boltzmann scheme is classically considered as second order accurate (see e.g. [23]). We describe in Section 6 the D3Q27 (d = 3, J = 26) numerical scheme that we used for acoustic application.

3) Formal development of eigenmodes
• With the method described in [11], it is possible to derive the equivalent system of (d + 1) equations associated to a lattice Boltzmann scheme. The system of linearized conservation equations issued from this analysis can be written under the form (11) A(∆x, ∂) • W = 0 which uses a compact notation stating that A(∆x, ∂) is a (d + 1) × (d + 1) (4 × 4 for three-dimensional applications, 3 × 3 for two-dimensional models) matrix of differential operators of high order acting on the conserved variables W .
We seek the eigenmodes of the operator A(∆x, ∂), id est the eigenvalues λ j (∆x, ∂) and the eigenvectors r j (∆x, ∂) such that (12) A We introduce the diagonal matrix Λ(∆x, ∂) composed of the eigenvalues λ j (∆x, ∂) and the square matrix R(∆x, ∂) composed of the eigenvectors. Then relation (12) can be written under the following synthetic form: • Moreover, the operator A(∆x, ∂) in equation (11) is a formal series in the variable ∆x that is truncated at order 4 in this contribution: We can apply in this case the so-called stationary perturbation theory for linear operators (see e.g. [3] for an elementary introduction). First for ∆x = 0, the operator A 0 (∂) is exactly the perfect linear acoustic model. For d = 3 (the case d = 2 is a simple consequence of this particular analysis and we omit the details), we have Introducing the Laplacian operator ∆ ≡ ∂ 2 x + ∂ 2 y + ∂ 2 z , a system R 0 (∂) of reference eigenvectors can be given under the form It is elementary to observe the relation and the associated matrix Λ 0 (∂) of order zero is simply given by We observe that even if the eigenvalue ∂ t is degenerated, i.e. when the eigenvalues are multiple, the family of eigenvectors proposed in (16) is complete and defines a basis. This allows us to show the existence of two acoustic waves and two shear waves (see e.g. [24]).
• The parameter ∆x is supposed to be infinitesimal and we introduce a formal expansion of the eigenvalues with diagonal matrices Λ j (∂): and relative perturbations Q j (∂) of the eigenvectors: We adopt in this work the scaling condition for eigenvectors proposed e.g. by Cohen-Tannoudji et al [3]. The component of r j (∆x, ∂) relative to the corresponding unperturbed eigenvector is exactly equal to unity. In other words, all the diagonal terms of the Q j (∂) matrices are null: We insert the expansions (20) and (21) into the eigenmode condition (13). We introduce the expression A j (∂) for the perturbation relative to the unperturbed eigenbasis, id est Then we identify each term relative to ∆x and we obtain in this manner: .
Explicit forms of the eigenvalues at each order may be found without difficulty using symbolic manipulation.
• We introduce a reference density ρ 0 , a reference sound velocity c 0 , the shear viscosity µ and the bulk viscosity ζ. It is well known (see e.g. Landau-Lifshitz [24]) that the four compressible isothermal acoustic linearized equations take the form: (15) when µ = ζ = 0. This system has two families of eigenvalues.
• On the one hand, the acoustic waves λ ± are given by and are parametrized by the attenuation of the sound waves γ given according to We can expand the previous acoustic wave eigenvalues in powers of γ when the attenuation is sufficiently small compared to the frequency of the sound waves. Then we have the expansion The previous (relatively unusual) notations with matrices of operators are exactly the one used when implementing without difficulty (due to linearity) our approach with a symbolic manipulation software. As is well known, the pseudo-differential operators like √ ∆ and ∆ 3/2 are defined by their action on Fourier transforms (we refer e.g. to the book of Hörmander [21]). It is also possible to introduce a Fourier decomposition on harmonic waves of the type exp i (ω t − k•x) . Then we have the usual change of notation: The corresponding dispersion relation associated with the annulation of the eigenvalues λ ± given in (29) allows the introduction of the parameter Γ ℓ according to (30) − • On the other hand, the shear wave λ 0 define a multiple eigenmode that is simply given by With the spectral notation, relation (31) can be written with the help of the attenuation Γ t of the shear waves under the form As previously, a family of two independent eigenvectors with two different "polarizations" corresponds to this shear eigenvalue. Moreover, the perturbation ∆x (14) is diagonal on the basis of unperturbed eigenvectors proposed at relation (16). Due to well-established arguments (see e.g. the textbook of Quantum Mechanics of Cohen-Tannoudji et al [3] or e.g. the sections "multiple roots" and "degenerate roots" in the book of Hinch [19]), the expansion (20) (21) is completely valid even in this (relatively) complicated degenerate case.

4) Increasing the accuracy of shear and acoustic waves
• When we use a lattice Boltzmann scheme, the viscosity µ admits a representation of the type (c.f. the D3Q27 scheme): and similarly the bulk viscosity ζ is given for the D3Q27 scheme (to fix the ideas) according to with the notations of Section 6. We can identify the developments (19) on the one hand and (29)(31) on the other hand because, due to (33), (34) and (28), the viscosities γ and µ are proportional to ∆x. The usual implementation of lattice Boltzmann scheme consists in adjusting the eigenvalues Λ 0 (∂) and Λ 1 (∂) in order to fit the viscosities.
• When we identify the developments (19) and (31) at the order 4 for the shear waves, we recover the "quartic parameters" obtained previously [11], in particular for the D2Q9 [23] and D3Q19 [32] schemes. In these cases, the shear viscosity µ follows the "quartic" condition If we add to the previous study the fitting of the acoustic waves (29), there is no solution of the equations for the two previous schemes D2Q9 and D3Q19. We have then considered the previous conditions for extended stencils such as the D2Q13 [30,35] and D3Q27 (see e.g. [28]) lattice Boltzmann schemes. The results for quartic coefficients are displayed in Section 6 when fitting all the waves (29) and (31) of the linear problem.
• In order to take advantage of the flexibility of the d'Humières scheme, we have developed an "isotropic" condition, not as strong as the quartic condition. With this isotropic condition, we do not impose anymore the conditions that the numerical waves (19) fit exactly the physical ones (29) and (31) but suggest that the numerical waves are isotropic.
In other words, the operator matrices Λ 2 (∂) (for the order 3) and Λ 3 (∂) (for the order 4) are differential operators that are constrained in order to depend only on the radial variable r (r 2 ≡ x 2 +y 2 +z 2 ) and on the operator ∂/∂r . When this condition is imposed, it is possible to fit isotropic waves for the four previous schemes. In Table 1 and 2, we display the number of (a priori non-independent) discrete equations that must be satisfied by the parameters of the lattice Boltzmann schemes when we impose isotropy for the shear waves, isotropy for the acoustic waves, exact fitting of the shear waves and exact fitting of the acoustic waves. All these discrete equations have been solved exactly with the help of symbolic manipulation. To fix the ideas, we give in Section 6 the parametrization of the quartic shear and acoustic waves for the D3Q27 lattice Boltzmann scheme (solution of 18 discrete equations, according to Table 2).
isotropic shear isotropic acoustic exact shear exact acoustic order 3 0 1 0 2 order 4 1 1 2 2 Table 1. Fitting "isotropic" and "exact" shear and acoustic waves at various orders of accuracy. Number of equations for the D2Q9 and D2Q13 lattice Boltzmann schemes. Note that the same numbers apply to both lattices.
isotropic shear isotropic acoustic exact shear exact acoustic order 3 0 1 0 2 order 4 10 2 13 3 Table 2. Fitting "isotropic" and "exact" shear and acoustic waves at various orders of accuracy. Number of equations for the D3Q19 and D3Q27 lattice Boltzmann schemes. Note that the same numbers apply to both lattices.

5) Linearized acoustics with the D2Q13 scheme
• Two numerical "experiments" have been performed: emission from the center ( Figures  1, 2 and 3) and emission from the boundary of the circle (see Figures 1, 4, 5, 6). In both cases, a first order anti-bounce-back boundary condition (see e.g. [1]) is implemented to impose a given density (pressure) on the boundary: either a constant when sound is emitted at the center, or a sinusoidal time dependent value for the emission by the boundary. In each case, we measure the value of the density at a fixed position in space vs time (see Figure 6) or we analyze the pressure field at a given time (see
Due to the large number of moments, we detail in this subsection the way the matrix M parametrizing the transformation (2) is obtained. First, velocities v α j for 0 ≤ j ≤ J ≡ 26 and 1 ≤ α ≤ 3 are given according to relation (36). The first four conserved moments ρ and q α are determined according to (3). The generation of other moments is done according to the tensorial nature of the variety of moments that can be constructed, as analyzed by Rubinstein and Luo [31]: scalar fields are naturally coupled with one another, idem for vector fields, and so on. So components of kinetic energy are introduced: The entire set of second order tensors is completed according to Flux of the energy and of the square of energy: Square and cube of the energy: Product of XX and WW by the energy: Extra-diagonal second order moments of the energy: two different components times energy and permutations with a natural permutation convention and the third order totally antisymmetric tensor XY Z: A first (non-orthogonal) matrix M is obtained by using (3) Then matrix M is orthogonalized from relations (3) The coefficients g iℓ are computed recursively in order to enforce orthogonality: • The equilibrium values of the moments define the equilibrium distribution. The first moments ρ and q α are in equilibrium and have respectively a scalar and a vectorial structure. The equilibrium values of the other moments respect this structure. We have e eq = θ λ 2 ρ XX eq = W W eq = XY eq = Y Z eq = ZX eq = 0 ϕ eq α = c 1 λ 2 q α ψ eq α = c 2 λ 4 q α ε eq = β λ 4 ρ e eq 3 = ξ λ 6 ρ XX eq e = W W eq e = XY eq e = Y Z eq e = ZX eq e = 0 τ eq α = c 3 λ 2 q α XY Z eq = 0 .
We have the following relation between the parameter θ and the sound velocity c 0 : The relaxation parameters for these moments are denoted respectively by the following values e eq : s e XX eq , W W eq , XY eq , Y Z eq , ZX eq : s x ϕ eq α : s ϕ ψ eq α : s ψ ε eq : s ε e eq 3 : s ξ XX eq e , W W eq e : s γ XY eq e , Y Z eq e , ZX eq e : s χ τ eq α : s τ XY Z eq : s ω .
With these relaxation parameters, the shear viscosity is isotropic when • The quartic parameters for athermal acoustics are solution of the 13 + 2 + 3 = 18 discrete equations obtained at Table 2. They are displayed in this subsection. (51) We observe that the presence of the factor (σ x − σ e ) in the denominator of the algebraic expression (54) of σ ε makes this expression incompatible with a BGK or TRT [15] type hypothesis. The same remark holds for relations (55) and (56) relative to the parameters σ γ and σ χ respectively. (55) x σ e + 8064 σ 4    Of course, we used the previous numbers with the usual truncation allowed by the 32 bit or 64 bit floating point arithmetic. We observe that for these particular parameters, due to the link (33) between the shear viscosity µ and the parameter σ x , to relation (34) between the bulk viscosity ζ and the parameter σ e and to formula (28) making explicit the attenuation γ of sound waves, we have    µ = 0.013 λ ∆x ζ = 0.0920492 λ ∆x γ = 0.0546913 λ ∆x .
• The results of simple periodic waves are displayed in Figures 7, 8 and 9. The numerical results confirm that parameters proposed in (59) and (60) allow one to get fourth order (relative) accuracy. Periodic wave with the lattice Boltzmann scheme D3Q27. Error for shear eigenvalue Γ t defined in (32). For quartic parameters, we have Γ t − ν |k| 2 = O(|k| 6 ).  Periodic wave with the lattice Boltzmann scheme D3Q27. Error for the imaginary part of the acoustic eigenvalue Γ ℓ defined in (30). For quartic parameters, we have fifth order accuracy: Im Γ ℓ − c 0 |k| 1 − γ 2 |k| 2 /(2 c 2 0 ) = O(|k| 5 ) with γ detailed at relation (28). The numerical experiments have been done for a three-dimensional converging acoustic wave. A pulsating sphere of radius R = 46.08 is embedded in a 95 3 cube. The simulations used only 64 bit or even 32 bit data with CUDA: (using a Nvidia 9800GT card, one full update of D3Q27 takes approximately 11 nanoseconds per site). The "linear anti-bounceback" numerical boundary conditions (in the sense given in Bouzidi et al. [1]) are used to impose a boundary density on the sphere oscillating with a period T = 10. Since the sound velocity is around 0.62 (see (59)), there are around six mesh points by wave length. In Figure 10, the scatter plot of density vs radius after 82 iterations is displayed. The results are of excellent quality with the use of quartic parameters. Nevertheless, since the attenuation of the sound wave is relatively important for the above parameters (γ ≈ 0.0547 as recalled in (61)), the range of propagation is relatively limited (five wave lengths in our case, see Figure 10). For the less stringent "isotropic" case, i.e. to fix the ideas we can find situations of "TRT type" satisfying (63) σ φ = σ ψ = σ τ = σ ω σ e = σ x = σ ǫ = σ ξ = σ γ = σ χ with much lower attenuation of sound. An example is given in figure 10 with parameters of the scheme following relations (62) and (63) with σ x = 0.006, that is shear viscosity 0.002 and sound attenuation 0.002. These results show that choosing the goal in terms of accuracy or small attenuation will influence the choice of the parameters in the lattice Boltzmann simulation.

7) Conclusion
• We have proposed to use the Taylor expansion method in conjunction with the use of symbolic manipulation to increase the accuracy of the lattice Boltzmann scheme in the multiple time relaxation approach proposed by D. d'Humières for linear acoustic waves. The result is the use of the previous scheme with very particular "quartic" parameters. We have obtained a family of such parameters for the D3Q27 numerical scheme. Numerical experiments confirm the predictions of the theoretical analysis. The catch is to make sure that situations are stable. This problem cannot be solved in practice with developments around small wave numbers. In the absence of efficient techniques to find stable situations with small viscosities and/or sound attenuation, we tried at random a number of sets of the free parameters and have used the best of them for the explicit results shown above. We note that the less stringent "isotropic" condition is compatible with small viscosities which may be quite useful for practical simulations. Comparison of various stencils is under way and will be presented elsewhere.