Lattice Boltzmann model for the simulation of the wave equation in curvilinear coordinates

TL;DR

This paper develops a lattice-Boltzmann model for simulating wave equations in curvilinear coordinates, enabling accurate analysis of vibrational modes in complex geometries like cylinders and tori.

Contribution

It introduces a novel lattice-Boltzmann approach for scalar wave propagation in curvilinear coordinates, expanding the method's applicability beyond traditional rectangular systems.

Findings

Accurately reproduces theoretical vibrational modes in complex geometries

Demonstrates the model's effectiveness in cylindrical, trumpet, and toroidal shapes

Shows potential for broad applications in physics and engineering

Abstract

Since its origins, lattice-Boltzmann methods have been restricted to rectangular coordinates, a fact which jeopardises the applications to problems with cylindrical or spherical symmetries and complicates the implementations with complex geometries. However, M. Mendoza recently proposed in his doctoral thesis a general procedure (based on Christoffel symbols) to construct lattice-Boltzmann models on curvilinear coordinates, which has shown very good results for hydrodynamics on cylindrical and spherical coordinates. In this work, we construct a lattice-Boltzmann model for the propagation of scalar waves in curvilinear coordinates, and we use it to determine the vibrational modes inside cylinders, trumpets and tori. The model correctly reproduces the theoretical expectations for the vibrational modes, and exemplifies the wide range of future applications of lattice-Boltzmann models on general curvilinear coordinates.

Tables

Table 1. Table 1: First three zeros of the Eq. 81

0rder	Value
1	1,2556
2	4,0793
3	7,1390

Equations

$$\begin{array}{r}{f}_i\left(\vec{x}+{\vec{\xi }}_i{\delta }_t,{\vec{\xi }}_i,t+{\delta }_t\right)-f_i\left(\vec{x},{\vec{\xi }}_i,t\right)=\\ -\frac{{\delta }_t}{\tau }\left[f_i\left(\vec{x},{\vec{\xi }}_i,t\right)-f_i^{eq}\left(\vec{x},{\vec{\xi }}_i,t\right)\right],\end{array}$$

$$P=\sum _i{f}_i,\vec{J}=\sum _i{f}_i{\vec{\xi }}_i.$$

$$f_i^{(0)}=f_i^{(eq)},$$

$$-\frac{1}{\tau }{f}_i^{(1)}=\left[\frac{\partial }{\partial t_1}+{\vec{\xi }}_i\cdot {\vec{\nabla }}_1\right]f_i^{(0)},$$

$$\frac{\partial }{\partial t_2}{f}_i^{(0)}+\left(1-\frac{1}{2\tau }\right)\left[\frac{\partial }{\partial t_1}+{\vec{\xi }}_i\cdot {\vec{\nabla }}_1\right]f_i^{(1)}=\frac{1}{\tau }{f}_i^{(2)}$$

$$-\frac{1}{\tau }{f}_i^{(2)}=\frac{\partial }{\partial t_2}{f}_i^{(0)}.$$

$$-\frac{1}{\tau }\left(ϵf_i^{(1)}+{ϵ}^2{f}_i^{(2)}\right)=\frac{\partial }{\partial t}{f}_i^{(0)}+\vec{\nabla }\cdot \left({\vec{\xi }}_i{f}_i^{(0)}\right),$$

$$\frac{\partial P}{\partial t}+\nabla \cdot \vec{J}=0.$$

$$\frac{\partial \vec{J}}{\partial t}+\nabla \cdot {\vec{\Pi }}^{(0)}=0,$$

$${{\Pi }^{(0)}}^{\alpha \beta }=\sum _i{f}_i^{eq}{\xi }_i^{\alpha }{\xi }_i^{\beta }=c^{2}P{\delta }^{\alpha \beta },$$

$$\frac{\partial \vec{J}}{\partial t}+c^2\nabla P=0,$$

$$\frac{{\partial }^{2}P}{\partial t^2}-c^2{\nabla }^{2}P=0,$$

$$f_i^{eq}=\{\begin{array}{ll}P+\frac{c^2}{c_s^2}P\left(w_0-1\right)& \text{i}fi=0\\ \frac{w_i}{c_s^2}\left(c^{2}P+{\vec{\xi }}_i\cdot \vec{J}\right)& \text{i}fi\ne 0\end{array}.$$

$$f^{eq}\left(\vec{\xi }\right)=\frac{P}{(\sqrt{2\pi c_s^2}{)}^D}{\text{e}}^{-(\vec{J}-\vec{\xi }{)}^2/2c_s^2}.$$

$$f_{i}^{eq}=\begin{cases}P-\left(\frac{5P}{2}-\frac{3c^{2}P}{2c_{s}^{2}}\right)\left(1-w_{0}\right)+\\ \frac{3c^{2}P}{2c_{s}^{2}}\left(c^{2}-c_{s}^{2}\right)&\text{if }i=0\\ \\ w_{i}\bigg{[}P+\frac{\xi_{i}\cdot J}{c_{s}^{2}}\\ +\frac{P}{2c_{s}^{4}}\left(c^{2}-c_{s}^{2}\right)\left(\xi_{i}^{2}-3c_{s}^{2}\right)\bigg{]}&\mbox{otherwise}\end{cases}\quad,$$

$$\frac{{\partial }^{2}P}{\partial t^2}-\frac{c^2}{\sqrt{g}}{\partial }_i\left(\sqrt{g}{g}^{ik}{\partial }_{k}P\right)=0,$$

$$(\nabla P{)}^l={\partial }^{l}P=g^{lk}{\partial }_{k}P,$$

$$\nabla \cdot \vec{J}=\frac{1}{\sqrt{g}}{\partial }_i(\sqrt{g}{J}^i),$$

$${\nabla }_j{\left({\Pi }^{(0)}\right)}^{ij}=\frac{1}{\sqrt{g}}{\partial }_j\left(\sqrt{g}{\left({\Pi }^{(0)}\right)}^{ij}\right)+{\Gamma }_{jk}^i{\left({\Pi }^{(0)}\right)}^{jk}.$$

$$\frac{\partial \left(\sqrt{g}P\right)}{\partial t}+{\partial }_i(\sqrt{g}{J}^i)=0,$$

$$\frac{\partial \left(\sqrt{g}{J}^i\right)}{\partial t}+{\partial }_j\left(\sqrt{g}{\left({\Pi }^{(0)}\right)}^{ij}\right)=-{\Gamma }_{jk}^i{\left({\Pi }^{(0)}\right)}^{jk}\sqrt{g}.$$

$${\left({\Pi }^{(0)}\right)}^{ij}=c^{2}P{g}^{ij},$$

$$\frac{\partial \left(\sqrt{g}{J}^i\right)}{\partial t}=-{\partial }_j\left(\sqrt{g}{c}^{2}P{g}^{ij}\right)-{\Gamma }_{jk}^i{c}^{2}P{g}^{jk}\sqrt{g}$$

$$\frac{{\partial }^2\left(\sqrt{g}P\right)}{\partial t^2}+{\partial }_i\left(\frac{\partial \left(\sqrt{g}{J}^i\right)}{\partial t}\right)=0,$$

$$\frac{\partial^{2}\left(\sqrt{g}P\right)}{\partial t^{2}}-c^{2}\partial_{i}\bigg{(}\partial_{j}\left(\sqrt{g}Pg^{ij}\right)+\Gamma^{i}_{jk}g^{jk}\sqrt{g}P\bigg{)}=0\quad\quad.$$

$$\sqrt{g}P=\sum _l{f}_l^{eq},$$

$$\sqrt{g}{J}^{\prime i}=\sum _l{f}_l^{eq}{\xi }_l^i-\frac{1}{2}{c}^{2}P{\Gamma }_{jk}^i{g}^{jk},$$

$$\sqrt{g}{\left({\Pi }^{(0)}\right)}^{ij}=\sum _l{f}_l^{eq}{\xi }_l^i{\xi }_l^j=\sqrt{g}{c}^{2}P{g}^{ij}.$$

$$f^{eq}\left(\vec{\xi }\right)=\frac{P}{(\sqrt{2\pi c_s^2}{)}^D}{\text{e}}^{-[g_{\alpha \beta }\left(J^{\alpha }-{\xi }^{\alpha }\right)\left(J^{\beta }-{\xi }^{\beta }\right)]/2c_s^2}.$$

$$f_{i}^{eq}=\begin{cases}w_{0}\sqrt{g}P\left[\frac{1}{w_{0}}-\frac{c^{2}}{2c_{s}^{2}}g^{\alpha\alpha}\right]-\sqrt{g}P+\frac{5}{2}\sqrt{g}Pw_{0}&\text{if }i=0\\ \\ w_{i}\sqrt{g}\bigg{[}P+\frac{\vec{\xi}_{i}\cdot\vec{J^{\prime}}}{c_{s}^{2}}\\ +\frac{P}{2c_{s}^{4}}\left(g^{\alpha\beta}c^{2}-c_{s}^{2}\delta^{\alpha\beta}\right)\left(\xi_{i}^{\alpha}\xi_{i}^{\beta}-c_{s}^{2}\delta^{\alpha\beta}\right)\bigg{]}&\text{if }i\neq 0\end{cases}\,.$$

Full text

I Introduction

Lattice-Boltzmann methods (LBM) were first introduced as mesoscopical models to simulate a wide variety of processes in fluid dynamics, from wind tunnels and turbulence to porous media and general rheology Succi (2001); but later on they were extended to more general systems, like waves Chopard et al. (1997); Chopard (2009); Guangwu (2000), electrodynamics Mendoza and Muñoz (2010), or even Quantum Mechanics Palpacelli and Succi (2008) and, therefore, they can be considered nowdays as a general numerical scheme to solve differential equations that can be written as a set of conservation laws. In comparison with other numerical schemes (like finite-differences or finite-element methods) all variables needed to compute for the next time step at a single node are present in the node itself (before being moved to the neighbouring cells), making them perfect to run parallel on graphic cards. Additionally, LBM models like the one for electrodynamics, has shown to be up to five times faster than the mentioned methods even in serial calculations for the same precision order Mendoza and Muñoz (2010). Because of these advantages, LBMs have gained the interest of a wide range of research areas and industrial applications; however, since most lattice-Boltzmann models assume a homogeneous and isotropic set of velocity vectors to move the information from node to node, the computational domain has been restricted to a rectangular array of cubic cells, forcing the use of staircase approximations on curved boundaries and imposing three-dimensional simulation domains for systems that, because of axial or spherical symmetries, were essentially two-dimensional, with an exorbitant increase in computational costs. This is also the case of simulating waves with lattice Boltzmann on acoustic systems. Symphonic and traditional instruments like violins, trumpets or drums, modern auditoriums, complex geological wells in seismic prospection and hearing organs like the Cochlea have too complex geometries to be properly described in Cartesian coordinates, asking for the need of new LBM for acoustics in generalized coordinates to take advantage of their versatility and parallel nature.

The simulation of waves is a wide area of research by itself. Waves are present almost in every phenomenon, from the surface patterns on water and the propagation of sound in solids and fluids to electromagnetic potentials and gravitational waves, and not all are related with fluids. Moreover, even in the case of fluids, focusing directly on reproducing the wave equation of interest reduces the complexity and computational costs, in comparison to simulate the whole fluid mechanics problem. The first approach to simulate the wave equation by LBM was proposed in 1998 by Chopard, Luthi and Wagen Chopard et al. (1997). The model modifies the equilibrium distribution function and the macroscopic variables to reproduce the wave equation, and with simple modifications runs on almost every velocity set, including D2Q5 and D3Q7. Later, in 2000 Guangwu Guangwu (2000) proposed to redefine the first macroscopic momentum of the LBM as the temporal derivative of the wave pressure. This also leads to the pressure wave equation, but an additional integration step is needed. In addition, some approaches use the standard LBM to simulate the Navier-Stokes equations (NSE) and meassure the pressure waves present in such dynamics, but with some restrictions and at a very high computational cost. In 2011, Li and Shan Li and Shan (2011) simulated the NSE with a Multiple Relaxation Time lattice-Boltzmann and studied the decaying of pressure signals due to the viscous, thermal and acoustic damping. Their model requires very high order approximations, which make it very demanding on computational resources, and applies only for low frequencies. In 2015 Sun et. al. Sun et al. (2015) used a forced lattice-Boltzmann model to solve the NSE for the fluid inside a fully saturated porous medium, and studied the propagation of a pressure front there. The absorption coefficients derived from the simulations coincide well with the theoretical expectations, but not with the experiments, due, perhaps, to the need of extra factors like thermal dissipation and flexibility of the porous medium. In 2016, Salomons, Lohman and Zhou Salomons et al. (2016), studied various cases of outdoor acoustics in 2D with a Multiple Relaxation Time lattice-Boltzmann model for the NSE, but the restrictions of minimal viscosity in the model avoids them to obtain the almost non-dissipative behaviour found in the real air. These limitations in the use of using LMBs for the NSE to simulate acoustics favours the use of simpler schemes directly designed for the wave equation, as the one proposed by Chopard et. al.

Various approaches have been proposed to overcome the Cartesian restriction in LBMs. One of the first proposals was made by Nannelli and Succi Nannelli and Succi (1992), consisting in a finite volume solution of the kinetic Boltzmann equation where the two-dimensional cells can be irregular quadrilaterals; but this scheme introduces a new coarse-grained distribution function (averaging the standard lattice-Boltzmann distributions) that should be retrieved through additional interpolation steps. Later, in 1997, He and Doolen He and Doolen (1997) implemented a lattice-Boltzmann scheme in polar coordinates, based on an algorithm proposed by He, Luo and Dembo He et al. (1996), where that additional interpolation step is still required. In 1998 a grid refinement scheme with boundary-fittings for complicated geometries was proposed by Filippova and Hänel Filippova and Hänel (1998). In this work, smaller cells are located in regions where a higher resolution is needed. Even though this strategy reduces some computational costs, the staircase approximation does not disappear and some interpolations on the curved boundaries of complex systems are still needed. Up to now, the implementation of lattice-Boltzmann models in curved general geometries has been a hard task, not only because of the insertion of additional computing steps, but mainly because each new geometry needs a new cautious discretization scheme to ensure a correct definition of the boundary conditions. In 2010 and 2015 Li et. al.Li et al. (2010) and Reijers, Gelderblom and Toschi,Reijers et al. (2016), respectively, developed a lattice-Boltzmann scheme for fluid mechanics to be implemented on axisymmetric coordinate systems. In particular, Reijers and co-workers constructed a bi-phasic lattice-Boltzmann model capable to achieve density ratios between phases up to 1000 and used it to study the propagation of pressure waves in such coordinates, but with an excessive dissipation due to numerical viscosity. In 2012 a more general approach in two dimensions was developed by Budinsky for both the shallow water and the Navier-Stokes equationsBudinski (2012). In that scheme, the equations are written in general coordinates, and the additional geometric terms (containing the Jacobian and Jacobian spatial derivatives) are introduced as forcing terms in the collision operator. The main advantage of this model is that all information about the curvilinear coordinates is included in the equilibrium distribution function and the forcing term, and no further discretization or interpolation steps are needed. Simultaneously with the Budinsky proposal, M. Mendoza Mendoza (2012) introduced in his doctoral thesis a new strategy to built lattice-Boltzmann models for fluids on any three-dimensional curvilinear coordinate system. The strategy also reproduces in the macroscopic limit the desired equations in generalized coordinates, but using the metric tensor and Christoffel symbols instead of the Jacobian . In contrast with Budinsky’s method, the forcing terms are included both in the equilibrium functions and in the macroscopic quantities, following the procedure by Guo et. al. Guo et al. (2002) and reaching second-order accuracy. Again, the strategy does not require neither a specific discretization scheme for each problem nor additional interpolation steps. This model has been successfully used to study the Dean’s instability in ellipsoidal coordinates Debus et al. (2014), the flow through randomly curved media Mendoza et al. (2013) and, more recently, the energy dissipation due to curvature Debus et al. (2017).

In this work we extend a modified version of the model by Chopard et. al Chopard et al. (1997); Chopard (2009) to simulate acoustic waves on generalized coordinates by following the proposal of M. Mendoza Mendoza (2012). The method, that can be used on any coordinate system, was tested by simulating the normal modes inside a cylinder, a trumpet and a torus, finding second-order accuracy. Section II reviews the LBM proposed by Chopard et. al and derives an alternative form by using Hermite polynomials, comparing their performance in the simple case of a point source in two dimensions. Section III extends that alternative form to generalized coordinates, including the general recipe to build the LBM for waves on any coordinate system. The model is tested in Section IV by simulating the acoustic waves inside a cylinder, a trumpet and a torus. Finally, Section V summarizes the main results and conclusions. Videos of the simulations can be found in the supplementary material attached to this manuscript.

II LBM for waves in Cartesian coordinates

Let us start from the lattice-Boltzmann’s equation with the Bhratnagar-Gross-Krook approximation Bhatnagar et al. (1954),

[TABLE]

where $\delta_{t}$ is the time step, we choose $\delta_{t}=1$ hereafter, $\tau$ is a characteristic relaxation time and $f_{i}\left(\vec{x},\,\vec{\xi}_{i},\,t\right)$ is the distribution function assigned to the velocity vector $\vec{\xi}_{i}$ 111In lattice-Boltzmann models for fluid mechanics, this distribution function is proportional to the probability to find a molecule at position $\vec{r}$ and time $t$ with velocity $\vec{\xi}_{i}$, but it is just a system’s variable in the general case. If the new distribution function (i.e. the first term on the left) were written at the same place, instead of moving to the neighbouring cells, this equation would represent an exponential decay to equilibrium. The left-hand side represents the time evolution of the distribution function without external forces and is called the advection term, whereas the right-hand side represents the interactions among distribution functions and is called the collision term. Note that the velocity space is reduced to a discrete set of velocity vectors $\vec{\xi}_{i}$, usually with some associated weights $w_{i}$, like the set D3Q7 described in Fig. 1. The macroscopic quantities are computed in terms of the discrete velocities and the distribution functions,

[TABLE]

In the case of acoustic waves, $P$ represents the pressure and $\vec{J}=-B\frac{\partial D}{\partial t}$ is proportional to the time derivative of the mean particle displacement $D$, with $B$ the bulk modulus.

The set of conservation laws the system fulfils in the macroscopic limit can be found by performing a Chapman-Enskog expansion Chopard (2009), of the evolution law Eq. (1). The left-hand side of such equation is replaced by a Taylor expansion up to second order and, both the distribution functions $f_{i}=f_{i}^{(0)}+\epsilon f_{i}^{(1)}+\epsilon^{2}f_{i}^{(2)}$ and the differential operators $\nabla=\epsilon\nabla_{1}$ and $\frac{\partial}{\partial t}=\epsilon\frac{\partial}{\partial t_{1}}+\epsilon^{2}\frac{\partial}{\partial t_{2}}$ are expanded on the Knudsen’s number $\epsilon$, such that $\epsilon\to 0$ is the continuous limit. By replacing and equating order by order, one obtains

[TABLE]

[TABLE]

[TABLE]

where we have used Eq. 4 to rewrite Eq. 5. If one chooses $\tau=1/2$, Eq (5) reduces to.

[TABLE]

By multiplying Eq. 4 by $\epsilon$ and Eq. 6 by $\epsilon^{2}$ and adding both equations, one recovers the time and space derivatives, obtaining

[TABLE]

where we use the fact that the velocity vectors $\vec{\xi}_{i}$ are the same from cell to cell. Now, by summing over the whole set of discrete velocities $\vec{\xi}_{i}$ and taking the limit $\epsilon\to 0$, we find a first conservation law

[TABLE]

Similarly, multiplying by $\vec{\xi}_{i}$ before summing give us a second conservation law

[TABLE]

where we defined $\Pi^{(0)}=\sum_{i}f_{i}^{eq}\vec{\xi}_{i}\otimes\vec{\xi}_{i}$.

Our goal now is to define the equilibrium distribution functions in such a way that those two conservation laws combine to obtain the scalar wave equation. Note that if we handle to obtain

[TABLE]

Eq. (9) becomes

[TABLE]

Now, by taking the time derivative of Eq. (8) and replacing into Eq. (11) one obtains

[TABLE]

that is, the wave equation for the pressure $P$ in Cartesian coordinates.

Finally, it is possible find the proper functional form for the equilibrium distribution in two ways: The first one consists in proposing an Ansatz for the equilibrium distribution to obtain the desired form for ${\Pi^{(0)}}^{\alpha\beta}$ and complete it to reproduce the macroscopic quantities Eq.(2) with $f_{i}^{eq}$ instead of $f_{i}^{0}$. If the velocity vectors $\vec{\xi}_{i}$ and weights $w_{i}$ fulfil $\sum_{i}w_{i}\xi_{i}^{\alpha}\xi_{i}^{\beta}=c_{s}\delta^{\alpha\beta}$, one can propose

[TABLE]

where the additional terms are added to obtain $P=\sum_{i}f_{i}^{eq}$ and $\vec{J}=\sum_{i}f_{i}^{eq}\vec{\xi}_{i}$. This is the equilibrium function proposed by Chopard et. al. Chopard et al. (1997). The second one consists in discretizing a proposed continuous equilibrium distribution function as a truncated a Hermite polynomial series Krüger et al. (2017). Let propose the continuous form

[TABLE]

which lead us to the discrete equilibrium distribution function (see AppendixVI.4)

[TABLE]

where the $f_{0}^{eq}$ contribution was found from $f_{0}^{eq}=P-\sum_{i=1}^{Q}f_{i}^{eq}$ (hereafter, Einstein summation is assumed).

Fig. 2 shows the wave patterns obtained through both methods for a point source in two dimensions. The simulation domain was $100\times 100\times 1$ cells, the wave speed was set to $c=0.5$ cells per timestep and the source at the central cell oscillates as $P=A\sin(\omega t)$, with $\omega=2\pi/\lambda$ and $\lambda=10$ cells. One can observe that both approaches reproduce the theoretical profile $P(r)=A\mathcal{J}_{0}\left(kr\right)\sin\left(\omega t\right)$, with RMS errors of 0.0487 for the Moment Matching approach and 0.0484 for the one obtained from Hermite polynomials expansion. Although the first approach is simpler, is the second the one we will use as starting point for our LBM in generalized coordinates.

III LBM for waves in generalized coordinates

Our strategy to build a LBM in generalized coordinates is to maintain a Cartesian array of cells to store the information in the computer, each dimension in this array representing a generalized coordinate in the real space. So, the velocity vectors and weights will be the same as in a traditional cubic discretization scheme. In contrast, the macroscopic quantities, equilibrium functions and forcing terms are chosen to simulate the conservation laws in generalized coordinates, as they were just novel and complicated equations to be represented on that Cartesian scheme.

Let us start by writing the differential operators in its generalized form and considering the additional terms as forcing terms. The wave equation in general coordinates reads

[TABLE]

where $g^{kl}$ is the contravariant metric tensor of the geometry and $g$, its determinant. Here, the second term is the explicit definition of the Laplace-Beltrami operator applied on the scalar quantity $P$. The gradient of a scalar function and the divergence of vector and tensor fields are given by Lawden (2002)

[TABLE]

[TABLE]

[TABLE]

Thus, the conservation laws (Eqs. (8) and (9)) now read

[TABLE]

[TABLE]

Note that the right-hand side in Eq. (21) can be considered as a source term of the conservation law, it is, therefore, a forcing term in the lattice-Boltzmann scheme. The forcing term can be treated as usual by following the approach by Guo et. al. Guo et al. (2002) with $\tau=1/2$ as shown below (Appendix VI.1). Observe that, if we handle to find an equilibrium distribution function such that the stress tensor were

[TABLE]

Eq.(21) would be

[TABLE]

and, by taking the time derivative of Eq. (20),

[TABLE]

the left hand side of Eq. (23) could be replaced in the second term of Eq. (24) to obtain

[TABLE]

This is exactly the wave equation in general coordinates (Eq. 16) we want to reproduce.

Next, because the macroscopic fields $P$, $\vec{J}$ and $\Pi^{(0)}$ in Eq. (24) and (25) appear multiplied by $\sqrt{g}$, let us define the statistical moments of the distribution functions such that

[TABLE]

[TABLE]

[TABLE]

The additional term in the righ-hand side of the macroscopic quantity $\vec{J}^{\prime}$ comes from the forcing term of the conservation law (Eq. (21)), as described by Guo et. al. for $\tau=1/2$.

The problem reduces again to find the equilibrium distribution functions such that Eq. (26-28) hold. First, we write a continuous equation that is just the extension of Eq. (14)

[TABLE]

Up to second order, the obtained form of the discrete equilibrium distribution function for the wave equation in curvilinear coordinates is (see Appendix VI.4)

[TABLE]

This form of the equilibrium distribution function evidences the existence of a conflictive term, the factor $\left(g^{\alpha\beta}c^{2}-c_{s}^{2}\delta^{\alpha\beta}\right)$, which accounts for how different is the metric of the involved curvilinear coordinate system from the Cartesian one. If the difference is too high, this term will produce numerical instabilities, inducing negative values for the equilibrium distributions. To diminish that effect, we move that term from the equilibrium function to the forcing term in $\vec{J}$, which becomes

[TABLE]

the whole forcing term is, therefore

[TABLE]

and the equilibrium function reduces to

[TABLE]

The spatial derivative $\partial_{j}\left[\sqrt{g}c^{2}P\left(g^{ij}-\delta^{ij}\right)\right]$ in Eq (31) is easily computed as (See Appendix VI.2)

[TABLE]

where $A^{ij}=\sqrt{g}c^{2}P\left(g^{ij}-\delta^{ij}\right)$. Eq. (34) is obtained from a Taylor expansion in the discrete velocity space and, therefore, it has the same properties of the overall lattice-Boltzmann scheme Thampi et al. (2012).

Summarizing, the model consists in solving the Eq. (1) with $\tau=1/2$ using the equilibrium distribution given in Eq. (33), which in turn is computed from the updated macroscopic variables obtained from Eqs. (26), (31) and (28). Note that all the information about the curvature of the system is included in the product $\sqrt{g}$ present in the Eq. (33) and Eqs. (26) and (28), and in the forcing term of the Eq. (31), but the set of velocities $\{\vec{\xi}_{i}\}$ remains the same. So, it is possible to adapt the geometry just by changing the metric tensor and Christoffel symbols, without altering the discretization scheme in the computer.

IV Three examples: Cylinder, Trumpet and Torus

As examples, we will measure the normal modes for three cases: A cylinder, a trumpet and a torus.

IV.1 Cylinder

For the first case, cylindrical coordinates are given by the transformation

[TABLE]

as shown in Fig. 3.

The metric tensor reads:

[TABLE]

and the non-zero Christoffel symbols are

[TABLE]

It is important to keep in mind that the information inside the computer is still stored in a cubic lattice; therefore, the visualization of the actual simulated system requires an inverse mapping (Fig. 4).

To validate the method, let us start by simulating a radial cylindrical wave expanding from the $z$ axis, whose theoretical solution is given by

[TABLE]

where $\mathcal{J}_{0}$ is the 0-th order Bessel function. To perturbate the system, we impose the macroscopic quantity $P$ at the axis to be a harmonic perturbation, $P(r,\theta,z,t)_{source}=C\sin\left(\omega t\right)\delta(r-r_{\rm min}),$ where $\delta$ is the Dirac’s delta and $r_{\rm min}$ is the minimum value of $r$ to be simulated (since cylindrical coordinates have an indetermination at $r=0$). For this simulation we use a cubic storage array of $100\times 3\times 10$ cells for the coordinates $r$, $\theta$ and $z$, respectively, with periodic boundary conditions in $\theta$ and $z$. For $r=r_{max}$ we imposed open boundary conditions (in other words, the simulation time is not enough for the wave to reach such boundary). Note that because of symmetry our field does not depend on $\theta$ and $z$, and the simulation could be one-dimensional in $r$ with a mesh $L_{r}\times 1\times 1$; but we choose to have some extra cells in $\theta$ and $z$ directions just to be sure that everything will run well for future simulations. The results for this benchmark (Fig. 5) are in excellent agreement with the theoretical predictions Eq. (38). In contrast with previous schemes Nannelli and Succi (1992); He et al. (1996); Filippova and Hänel (1998), no additional interpolation steps in the evolution are required. All the values of the quantities shown here and henceforth are in lattice units. for this simulation we have chosen $c=0.6$, $\omega=0.075$, $r_{\rm min}=1$ and $r_{\rm max}=500$.

Now, let us simulate an acoustic wave inside an open-closed end pipe whose dimensions in the computational domain in directions $r,\,\theta,\,z$ are $20\times 5\times 100$ cells with $r_{\rm max}=25$ and $z_{\rm max}=100$. Fig. 6 describes the boundary conditions imposed: Free boundary conditions are imposed at the end of the pipe at $z=L$ as usual by copying the values of the neighbouring cell normal to the surface. In addition, because the axis $r=0$ is now part of the simulation space (in contrast to the previous case, where the axis was a source boundary condition), we need to take care about the coordinates’ singularity at $r=0$ by also implementing a free boundary condition at $r=r_{\rm min}=1$ (Fig.3). Bounce back boundary conditions are considered in the pipe walls at $r=r_{\rm max}$ and in the closed end at $z=0$ by interchanging the distribution functions travelling on opposite velocity vectors. Finally, some dissipation must be added to reach steady-stable oscillations after a transient. This is done just by adding a damping factor for the distribution functions at the bounce-back boundary condition after the interaction with the walls,

[TABLE]

where the indexes $in$ and $out$ means the incident and reflected directions, respectively. For all simulations we choose $D=0.65$, which was enough to reach stable normal oscillation modes after 8400 timesteps.

As a first case, let us consider that the left end oscillates as a rigid body like a moving uniform plate; therefore, the functional form of the source is $J^{\prime}_{z}(r,\theta,0,t)=C\sin\left(\omega t\right)$. This wave propagates along the $z$ direction only (as in a one-dimensional cord) and, we could expect that the resonant frequencies correspond to the open-end pipe in one dimension given by $\omega_{n}=\frac{2\pi c(2n+1)}{4L_{z}}$. However, these frequencies are shifted and some peaks are suppressed (see Fig.7) due to the appearance of radial modes that cannot be suppressed in this geometry, as will be shown below. Next, we considered that the left end is a fixed wall plus an oscillating ring in front of it. The ring’s radius is chosen to be of two cells, such that the ring fits just around the cylinder at $r=0$ that has been excluded from the simulation space. The functional form of this source is $P(r_{min},\theta,0,t)=C\sin(\omega t)$. The expected normal modes are now three-dimensional with axial symmetry. The pipe imposes a von Newmann boundary condition at the cylindrical wall,

[TABLE]

Under such conditions the characteristic frequencies are (see Appendix VI.3):

[TABLE]

In order to identify the resonant frequencies for the two different perturbations, we vary the source frequency from [math] to $0.03$ oscillations per timestep measuring the maximum pressure at a given point inside the pipe ($r=3,\,z=z_{\rm max}/10$). For these simulations, we set $c=0.1$. Fig. 7 shows the simulations results. The solid vertical lines are the expected values for the resonant frequencies obtained from Eq.(41), the dashed line corresponds to the rigid uniform plate and the solid line, to the oscillating ring. Note that in the case of the uniform ring, longitudinal modes are more defined but shifted due to the appearance of the first radial mode at $\omega=0.164$, in addition, some resonant peaks that are present for the case of the ring were suppressed in the spectrum of the plate. The resonant frequencies simulated with the oscillating ring as source are in great agreement with the expected values, and the RMS deviation between the maxima and the black vertical lines is less than 1$\%$.

Fig.8 shows a snapshot of the pressure waves inside the cylinder for the ring perturbation.

Next, let us measure the model convergence. To this aim, we choose a theoretically known quantity, in this case the third characteristic frequency of the open-closed end pipe, to be compared with its simulated value. As the resolution increases, we expect that the smaller the cell, the smaller the difference between the theoretical value and the simulated one. Such behaviour is shown in Fig. 9. The error grows with cell size as a power law with exponent 2, that is the model convergence is second order.

IV.2 Trumpet

The coordinate system Eq. 35 can be modified to include a flared end resembling a trumpet’s bell. In this work, we consider a particular shape named Bessel Horn, given by the following coordinate system

[TABLE]

where $\lambda$ is a positive real number and $\mathcal{Z}=Z_{\text{max}}-z$ with $Z_{\text{max}}$ the trumpet’s length (Fig.10).

The surfaces for $r=\text{const}$ resembles the shape of a different horn for each value of $\lambda$ (Fig. 11).

The Non-zero Christoffel symbols calculated for these coordinates in terms of $\lambda$ are given by

[TABLE]

and the metric tensor, by

[TABLE]

By using Eqs. 42, 43 and 44 it is possible to simulate the acoustical waves inside a trumpet. The simulation space and boundary conditions are the same to the ones for the pipe with an oscillating ring (Fig. 6). Fig. 12 shows a snapshot after 8400 timesteps of the pressure profile for an axisymmetric Bessel horn with $\lambda=0.2$, $r_{min}=2$, $r_{max}=150$ and $Z_{\text{max}}=750$, at a frequency $\omega=18\pi c/4z_{\rm max}$. As expected, the intensity of the pressure wave is reduced at the end of the trumpet due to the bell. Fig. 13 compares the oscillation envelope as a function of $z$ at $r=3$ cells when the ring oscillates at fixed frequency $\omega=0.005$ for two shapes: a horn and a cylindrical pipe. Note that the amplitude of the oscillations decrease near the bell when it is present (at highest values of $z$ when $\lambda\neq 0$)

Let us now study how the natural frequencies change by the inclusion of the bell. Figure 14 shows the resonant peaks for a trumpet with $\lambda=0.1$, compared with the previously obtained results for $\lambda=0$. As expected, the longitudinal normal modes for the trumpet are less defined, and the second radial modes are strongly suppressed by the bell, similarly to what is observed in real pipes Rossing et al. (2002).

IV.3 Torus

As a third benchmark, we study the pressure waves inside a torus. A torus can be seen as a circle of radius $r_{\rm max}$ whose center travels along a larger circle of radius $R$ (See Fig. 15). It can be described by the coordinate transformation

[TABLE]

The metric tensor for this coordinates is, therefore,

[TABLE]

and the non-zero Christoffel symbols are

[TABLE]

The chosen torus has $R=640$, $5.12<r<409.6$, $0<\theta<2\pi$ and $0<\phi<2\pi$. For the simulations we used 20 cells in $r$, 200 cells in $\theta$ and 30 cells in $\phi$. The boundary conditions imposed are bounce-back with damping factor $D=0.65$ in $r_{\rm max}$, free in $r_{\rm min}$ and periodic both in $\theta$ and $\phi$.

To test the LBM we run three different simulations, each one with a different source. The first one is a ring along the major circle of radius $R$, and its functional form is $P_{source}(r,\theta,\phi,t)=A\sin(\omega t)\delta(r-r_{min})$. Fig. 16 shows the evolution of the pressure wave, showing expected radial symmetry.

The second source is similar to the small ring employed in pipes and trumpets, with functional form $P_{source}(r,\theta,\phi,t)=A\sin(\omega t)\delta(r-r_{\rm min})\delta(\theta)$. This source allows the excitation of normal modes in both $r$ and $\theta$ directions. The results for different frequencies are shown in Fig. 17.

We also studied the resonant frequencies for this geometry. To this aim we located the small oscillating ring at $\theta=0$ and measured the pressure intensity at the opposite point on the major circle (i. e. at $\theta=\pi$) after 8000 time steps, when the wave has reached a steady-stable condition. Fig. 18 shows the obtained intensity for frequencies between 0 and 0.5 radians per click varying the frequency in steps of $\Delta\omega=0.001$.

Note that the larger resonant peaks are placed whitin 0.05 and 0.1 radians per click, with a cutoff frequency near to 0.25 radians per click. Although the analytical solution for the normal modes of the torus was not found in the literature, it is interesting to point out that the pattern is very similar to the one for a trumpet with an open bell Rossing et al. (2002).

To simulate more general perturbations, the third source we used is an oscillating point outside any symmetry axis, i.e. at $r=\frac{r_{max}-r_{min})}{2}$, $\theta=0$, and $\phi=\pi/2$. The functional form of the source is, therefore, $P_{source}(r,\theta,\phi)=A\sin(\omega t)\delta\left(r-\frac{r_{max}-r_{min})}{2}\right)\delta(\theta)\delta(\phi-\pi/2)$. The resulting pattern is shown in Fig. 19. We can observe that the pattern is no longer symmetric in $\phi$ and shows variations along $r$.

For this last case, the evolution inside the torus is general, since it is not restricted to any type of symmetry, showing that our LBM can be applied to non-axisymmetric systems.

V Conclusions

A lattice-Boltzmann model (LBM) for the simulation of acoustic waves in general curvilinear coordinates was developed and implemented. The model considers the additional terms related to the wave equation in curvilinear coordinates as forcing terms of the conservation laws and uses the approach of Guo et. al. to construct the proper forced lattice-Boltzmann model. The model was tested by reproducing the propagation of a cylindrical wave, by finding the normal modes inside a cylindrical pipe and a Bessel horn trumpet and by computing steady-state pressure profiles inside a Torus. Free boundary conditions were set at the open ends, and bounce-back boundary conditions at the walls, with a damping factor in the lateral ones. The simulation results show very good agreement with the theoretical expectations: both the RMS errors for the cylindrical wave and the resonant frequencies for the pipe are less than 1$\%$. In addition, the results for the trumpet show a strong suppression of the second radial modes, similar to the one observed in real trumpets. For the case of the Torus, the resonant spectrum shows an interesting similarity with the one for a trumpet and, although we did not find theoretical expressions for the normal modes, our results are in qualitative agreement with the ones expected for that geometry.

Unlike previous models, our proposed LBM is completely general and does not ask for additional interpolation steps. The geometric information of the coordinate system is completely included in the equilibrium distribution function and in the macroscopic quantities, but the system stored in the computer is still a cubic lattice. For this reason, it is not necessary a new discretization scheme, i.e. the discrete velocity set is defined as usual and the velocity vectors are the same for all cells. The model also allows to take advantage of the symmetry of the system to reduce the number of cells to reach certain precision, reducing computational requirements. For instance, our simulations of pipes and horns are actually two-dimensional in cylindrical coordinates, taking advantage of the azimutal symmetry, and therefore, they consume fewer computational resources. Those versatilities widely opens the application range of lattice-Boltzmann models to acoustic systems with very complicated geometries, like more complex musical instruments, auditoriums and concert halls or even detection organs as the Cochlea. The convergence is second order, which keeps the accuracy of the standard lattice-Boltzmann models.

Summarizing, this work introduces a lattice-Boltzmann model for the simulation of acoustic waves in general curvilinear coordinates. The manuscript has written in such a way that the interested reader can implement any desired geometry without effort, just by computing the metric tensor and Christoffel symbols and inserting them into the model. So, hopefully, it will be of great utility in future research.

VI Appendix

VI.1 Chapman-Enskog Expansion for a forced lattice-Boltzmann

In this appendix we will show how the microdynamics of the lattice-Boltzmann equation, including a forcing term, can lead to a set of conservation laws with source terms like Eqs. (20) and (21) for any lattice-Boltzmann scheme, in particular, the forces considered in this work. We will obtain Eq. (27) as the momentum correction needed to include such forcing terms in the lattice-Boltzmann scheme. The approach here described is based on the analysis made by Li, Zhou and Yan in 2016 Li et al. (2016)

Let us start by writing the lattice-Boltzmann equation with an additional forcing term discrete in the velocity space, $F_{i}$

[TABLE]

now, let us perform the Chapman-Enskog expansion with this additional term. First, we expand the left-hand side of the Eq. (48) in a Taylor series up to second order which gives

[TABLE]

where $f_{i}=f_{i}\left(\vec{x},\,\vec{\xi}_{i},\,t\right)$ and $\delta_{t}=1$. Second, we make a perturbative expansion in the small parameter $\epsilon$ of the distribution functions, the differential operators and the forcing term up to second order as follows $f_{i}=f_{i}^{(0)}+\epsilon f_{i}^{(1)}+\epsilon^{2}f_{i}^{(2)}$, $\nabla=\epsilon\nabla_{1}$, $\frac{\partial}{\partial t}=\epsilon\frac{\partial}{\partial t_{1}}+\epsilon^{2}\frac{\partial}{\partial t_{2}}$ and $F_{i}=\epsilon F^{(1)}_{i}$. By replacing the perturbative expansions in Eq. (49), we get

[TABLE]

equating order by order, we obtain for the zero-th first and second order, respectively

[TABLE]

[TABLE]

[TABLE]

by using the Eq. (52) we can write Eq.(53) in the form

[TABLE]

or, equivalently

[TABLE]

If we sum Eqs. (52) and (55) over $i$ and use the definitions of the macroscopic variables (Eq. (2)) and the Eq. (51), we obtain

[TABLE]

[TABLE]

where we have assumed $\sum_{i}F_{i}=0$ and $\sum_{i}\vec{\xi}_{i}F_{i}=\vec{\mathcal{F}}$ as suggested by Guo et. al.Guo et al. (2002). By adding the Eqs.(56) and (57) the first conservation law for the lattice-Boltzmann scheme reads

[TABLE]

with

[TABLE]

Similarly, if we multiply by $\vec{\xi_{i}}$ before summing over $i$, we obtain the conservation laws of higher order

[TABLE]

[TABLE]

by summing these two equations and taking the choice $\tau=\frac{1}{2}$ for the wave equation, we obtain

[TABLE]

Again, the functional form of the forcing term in this case have to be chosen in order to fulfil $\sum_{i}\vec{\xi}_{i}\vec{\xi}_{i}F^{(1)}_{i}=0$. Eqs. (58) and (62) are exactly the conservation laws we use for our forced lattice-Boltzmann scheme (Eqs. (20) and (21)) plus a rescaling factor of the macroscopic variables due to the geometry $\sqrt{g}$. Note, in addition, that the Eq.(59) corresponds to Eq. (27).

Finally, for the case of the LBM for the wave equation a simplification can be made on Eq. (48) when we insert the value of $\vec{J}^{\prime}$ into the equilibrium distribution function i.e. $f_{i}^{eq}\left(P,\vec{J}\right)\to f_{i}^{eq}\left(P,\vec{J}^{\prime}\right)$. Let us write the equilibrium distribution function in the form

[TABLE]

which allow us to rewrite the Eq. (48) as

[TABLE]

where

[TABLE]

and the first three moments of this force are

[TABLE]

[TABLE]

[TABLE]

Therefore, it wont have effect on the computed macroscopic fields and we can write:

[TABLE]

which is the evolution equation we use in our model.

VI.2 Fields derivatives computed using a cell configurationThampi et al. (2012)

In some cases, it is necessary to compute the derivative of any scalar or vectorial field during a lattice-Boltzmann calculation, this can be done by implementing a discretization scheme like finite differences scheme, however, it can either reduce the accuracy order of the overall scheme or increase the computational time and model complexity. However, there is an useful alternative to compute the gradients of the fields that are involved in the lattice-Boltzmann simulation.

The procedure to find these gradients is described here. Let us start from a Taylor expansion of the field in the direction of a given discrete velocity from the discrete set used in the lattice-Boltzmann calculations.

[TABLE]

now, we multiply the whole expansion by $w_{i}\vec{\xi}_{i}$

[TABLE]

if we sum over $i$, we can use the isotropy conditions (LABEL:Weights) to cancel the terms of odd order in $\vec{\xi}_{i}$

[TABLE]

and use

[TABLE]

to write

[TABLE]

finally, if we divide by $c_{s}^{2}$ we obtain the required expression to compute the gradient of a scalar field by using the same discretization scheme implemented in the lattice-Boltzmann.

[TABLE]

It can also be generalized to compute the gradient of any tensor $A^{\alpha\beta}$

[TABLE]

Note that the order of accuracy is the same of the overall lattice-Boltzmann scheme, even if the chosen discrete velocities set is only second order as in our case (D3Q7).

VI.3 Radial modes of the pipe.

In order to define correctly the boundary condition for the rigid walls of the pipe and find a theoretical expression for the vibrational modes along the radial direction, lets first write the differential equation for $r$ obtained from the variables separation of the wave equation

[TABLE]

whose solution is

[TABLE]

where $J_{m}\left(\zeta r\right)$ and $Y_{m}\left(\zeta r\right)$ are the Bessel functions of first and second kind respectively, the constants $A_{r}$ and $B_{r}$ can be obtained from the boundary conditions at $r=0$ and $r=r_{max}$. For $r=0$, since the pressure have to be finite and the Bessel Function of the second kind $Y_{m}(\zeta r)$ diverges at $r=0$, we have $B_{r}=0$. Now replacing the Eq. 78 into the Eq. 40, we find the condition

[TABLE]

which can be numerically solved for $\zeta$ to find the characteristic radial frequencies of the pipe. In order to simplify the above equation and its solution we use the property

[TABLE]

to find the function

[TABLE]

whose zeros correspond to $\zeta_{l}r_{max}$, where $\zeta_{l}$ is the $l$-th radial characteristic frequency of the pipe.

For the simulations done in this work, we assume axial symmetry of the waves and therefore we choose $m=0$, the Figure 20 shows the graphic of the Eq. 81 as a function of $\zeta r_{max}$ for $m=0$

The first three zeros of the function happen to be

Finally, the theoretical expression for the vibrational modes inside the pipe with $m=0$ considering both the radial and longitudinal vibrations is given by

[TABLE]

VI.4 Hermite expansion of the equilibrium distribution function

Let us start by considering the conservation laws from the continuous Boltzmann equation. The statistical moments of the distributions are conserved under the collision process Krüger et al. (2017); therefore, its value is the same whether they are computed from $f$ or from its equilibrium value $f^{eq}$,

[TABLE]

Our aim is to find a simpler form for this relations, but preserving the values of the macroscopic fields $P$, $\vec{J}$ and $\Pi^{\left(0\right)}$. In fluids, which are the macroscopic systems the Boltzmann equation is intended to model, the equilibrium distribution is Maxwell-Boltzmann, which is a Gaussian distribution,

[TABLE]

Where $\rho$ is the fluid density and $\vec{U}$ its velocity, $R=k_{B}/m$ is the ideal gas constant, with $k_{B}$ the Boltzmann constant and $m$, the mass of the particles; $T$, the temperature and $D$, the system dimension. The clue is given by regarding that the continuous distribution function has the same functional form of the weight function for Hermite Polynomials,

[TABLE]

which are defined by

[TABLE]

and are orthogonal with the dot product

[TABLE]

For the case of the wave equation, nevertheless, there is not a continuous distribution function that can be taken as a starting point to find its discrete version for LBM. So, we have to propose an appropriate continuous form to successfully retrieve the wave equation in the macroscopic limit. One can propose

[TABLE]

[TABLE]

for the Cartesian and generalized case respectively. Since the Hermite polynomials are orthogonal and form a basis, we can write any function $f(x)$ as a multidimensional Hermite polynomial series of the form

[TABLE]

where, the coefficients $\vec{a}^{\left(n\right)}$ can be obtained as

[TABLE]

That is valid on each direction in velocity space,

[TABLE]

Note that one of the reasons for choosing Hermite polynomials is that the series coefficients directly correspond with the statistical moments of the distribution i. e. the system macroscopic quantities.

[TABLE]

Next, the expansion can be truncated up to certain order $N$

[TABLE]

If we truncate the expansion up to second order we can write

[TABLE]

Another important feature of the Hermite polynomials is called Gauss-Hermite Cuadrature, which allow us to find exactly the value of an integral of a weighted polynomial of grade $n$, $P^{(n)}$, by considering a discrete sum over precise values $\xi_{i}$, which are actually the roots of the Hermite polynomial $\mathcal{H}^{(n)}$,

[TABLE]

where $N$ must satisfy $n\geq\left(N+1\right)/2$. We can use the Eq. (97) to easily compute the coefficients $\vec{a}^{\left(n\right)}$ (Eq. (91)),

[TABLE]

The coefficients are

[TABLE]

The corresponding equilibrium distribution function is

[TABLE]

The expression for $f_{0}^{eq}$ can be obtained from

[TABLE]

Summarizing, we have

[TABLE]