EDP Sciences logo
Submit your paper
Search
Menu
All issues Volume 8 (2024) Acta Acust., 8 (2024) 75 Full HTML
Open Access
Issue
Acta Acust.
Volume 8, 2024
Article Number 75
Number of page(s) 13
Section Computational and Numerical Acoustics
DOI https://doi.org/10.1051/aacus/2024071
Published online 13 December 2024

© The Author(s), Published by EDP Sciences, 2024

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Room acoustic simulation techniques have been extensively studied in recent years, with two main categories emerging: geometric and wave-based methods [1]. Geometric approaches, such as the ray tracing method, the image source model or the beam tracing method [2], assume sound propagation as a ray, i.e., wave properties are ignored. In contrast, wave-based methods numerically approximate the solution of the governing acoustic partial differential equations [3], which also provide solutions for complex wave effects such as diffraction and scattering. Various numerical discretization techniques have been established, such as the finite difference time domain method (FDTD) [4, 5], the finite element method (FEM) [6, 7], the boundary element method (BEM) [8, 9], the finite volume method (FVM) [10, 11], the discontinuous Galerkin finite element method (DGFEM) [12] and the spectral element method (SEM) [13].

Recently, the interest in wave-based methods has increased [14–17], since the high computational requirements can be met by increasing computational resources, while geometrical acoustic methods are still unable to adequately model wave-related phenomena such as the behavior of acoustic obstacles or reflectors of finite size [18]. In the case of room acoustic simulations the governing equations are usually written as the second-order wave equation [5, 19, 20]. Alternatively, an equivalent system of first-order linear acoustic equations (LAE) is employed [4, 21]. Recently, also the full nonlinear Euler equations (NLEE) have been used to numerically simulate indoor and outdoor sound propagation [22, 23]. The NLEE imply nonlinear effects and flow, which are of interest in several fields where high sound pressure levels or thermodynamic properties such as air-conditioning flow or heat distribution have to be considered, e.g., in open air sound reinforcement. Also in musical acoustics nonlinear effects are involved in the generation and propagation of sound in musical instruments [24], e.g., in brass instruments [25, 26] or in string instruments [27, 28].

The numerical accuracy and efficiency required for these applications can be achieved through FDTD discretization. This approach is particularly advantageous because the numerical mathematics used is common to the engineering community, and high-order and high-quality simulations can be achieved while maintaining numerical efficiency and scalability for high-performance computing. The governing equations are discretized in space and time. For the spatial derivatives of the first-order partial differentials (LAE, NLEE), explicit and implicit schemes are available that require only a few points per wavelength (PPW) while ensuring good dissipation and dispersion properties [29, 30]. Explicit integration methods are typically used for time discretization. Multi-step methods such as the Adams-Bashfourth method [31] or one-step methods such as the Runge-Kutta method [32, 33] are commonly used. When using nonlinear governing equations, the suppression of unphysical waves, which are generated by the FDTD method, can be achieved by modifying of the governing equations in terms of Selective Artificial Damping (SAD) [34] or by filtering [35].

The boundary conditions are a crucial part of the quality [1] and a challenging part of FDTD simulations. Characteristic boundary formulations allow for reflective or non-reflective boundary treatment [36]. Perfectly matched layers (PML) [37] and sponge layers [38] are used for improved non-reflective boundary conditions. In addition, there are several approaches to modeling complex surface impedances. Most models used in room acoustics distinguish between the complexity of the absorbing material, i.e., the absorptivity of the material [14]. These include local reaction (LR) models [39] and extended reaction (ER) models [40]. In general, absorption behavior varies with the angle of incidence. Recent focus has been on volume penalization methods [41], which can be interpreted as a subset of the more general equivalent fluid method [42–44]. The methods show good performance for LR and ER absorption behavior and are easy to implement in NLEE and LAE [45–48]. A review of recent trends in immersed boundary methods is given in [49].

This paper discusses the applicability of the aforementioned FDTD methods for room acoustic simulations using NLEE and LAE, with focus on the governing equations, spatial and temporal discretization as well as boundary conditions. To demonstrate the computational feasibility and performance of this approach, three benchmark cases are simulated, including linear cases with analytical solutions and room acoustic measurements as references as well as a nonlinear case with high sound pressure levels.

The main message of this paper is that, in addition to the commonly used second-order wave equation, the (linearized or non-linearized) first-order Euler equations can be used to simulate acoustic phenomena. In particular, the latter ensures that complex factors such as non-linearities and inhomogeneous environmental conditions, including wind flow or temperature gradients, can be taken into account, which can be particularly relevant in scenarios such as large air-conditioned halls and outdoor sound reinforcement.

The manuscript is organized as follows: Section 2 reviews the governing equations, spatial and temporal discretization schemes, filtering schemes, and boundary conditions in detail. Depending on the application, recommendations for the optimal scheme are given. In Section 3 spatial discretization and filtering schemes are analyzed in a 3D cubic room with rigid boundaries to prove their impact on the transmission behavior. In Section 4, the reflection scene 1 of the BRAS database [18] is compared with the experimental results. The difference between non-linear and linear equations in the context of high sound pressure levels is addressed in Section 5. Finally, the results are summarized and discussed in Section 6.

2 Methods

2.1 Governing equations

Depending on the intended application and the effects to be studied, the governing equations must be chosen. The full compressible Navier-Stokes equations [50] include fluid dynamics and acoustics. For purely acoustic applications, they can be reduced to the full nonlinear Euler equations (NLEE) [51]

∂t(ϱϱujpγ-1)+∂xi(ϱuiϱuiuj+pδijuipγγ-1)-ui∂xi(00p)=(000)$$ {\partial }_t\left(\begin{array}{l}\varrho \\ \varrho {u}_j\\ \frac{p}{\gamma -1}\\ \end{array}\right)+{\partial }_{{x}_i}\left(\begin{array}{l}\varrho {u}_i\\ \varrho {u}_i{u}_j+p{\delta }_{{ij}}\\ \frac{{u}_i{p\gamma }}{\gamma -1}\end{array}\right)-{u}_i{\partial }_{{x}_i}\left(\begin{array}{l}0\\ 0\\ p\end{array}\right)=\left(\begin{array}{l}0\\ 0\\ 0\end{array}\right) $$(1)

which describe all effects of sound propagation including nonlinear effects. Here, the viscous stresses are neglected. In equation (1) γ is the heat capacity ratio, uj is the velocity in the xj-direction, ϱ is the density, p is the pressure, and Δij is the Kronecker delta. Therein ∂t = ∂/∂t and ∂xi=∂/∂xi$ {\partial }_{{x}_i}=\partial /\partial {x}_i$ holds. For the sake of brevity an index notation is used. Einstein’s summation convention applies for i, j = [1–3]. Note that in addition to acoustics, laminar flows, and thermodynamic quantities such as temperature gradients can be studied with this set of equations.

If nonlinear effects can be neglected, the full linear Euler equations (LEE) can be derived from the NLEE using the decomposition approach ϱ = ϱ0 + ϱ′, uj = u0,j + uj′ and p = p0 + p’. The variables ρ0, u0, p0 denote the base flow quantities, i.e., the state in the absence of acoustic perturbations. The LEE are given in [52]. Considering no flow, i.e., u0 = 0 and a constant atmospheric pressure, the LEE reduce to linear acoustic equations (LAE)

∂t(p'u'j)+∂xj(ϱ0c02u'j1ϱ0p')=(00),$$ {\partial }_t\left(\begin{array}{l}p\prime\\ {u\prime}_j\end{array}\right)+{\partial }_{{x}_j}\left(\begin{array}{l}{\varrho }_0{c}_0^2{u\prime}_j\\ \frac{1}{{\varrho }_0}p\prime\end{array}\right)=\left(\begin{array}{l}0\\ 0\end{array}\right), $$(2)

with the acoustic relation

c0=p'ϱ'.$$ {c}_0=\sqrt{\frac{p\prime}{\varrho \prime}. $$(3)

2.2 Spatial discretization

In FDTD, the spatial derivatives are approximated by finite differences, which can be derived by combining different Taylor expansions or deriving a suitable interpolation stencil. The stencil schemes used in this manuscript are given in Appendix. Besides the stability, two other quality criteria are important for a good approximation: the order of consistency [53], which should be chosen as high as possible, and the dispersion behavior. Please note, that the following analysis considers first-order partial differential equations as in equations (1) and (2).

In order to analyze the stencil properties with respect to different wavenumbers, the dispersion error can be analyzed in one dimension using the approach

qm∼exp(jkxm)$$ {q}_m\sim \mathrm{exp}(\mathrm{j}k{x}_m) $$(4)

where qm denotes the field variable at point xm = mΔx, with m = [1, …, M], M the number of equidistantly distributed grid points and Δx as spatial step size. Furthermore, k is the physical wavenumber, and the imaginary unit j2 = −1. The effective wavenumber k′ describes the numerical computed wave number. Hence, ideally the relation of k and k′ should be unity. Using a simple second-order symmetric stencil

∂xqm=qm+1-qm-12Δx+O(Δx2).$$ {\partial }_x{q}_m=\frac{{q}_{m+1}-{q}_{m-1}}{2\Delta x}+\mathcal{O}\left(\Delta {x}^2\right). $$(5)

Equation (4) results in

k'Δx=sin(kΔx).$$ k\prime\Delta x=\mathrm{sin}(k\Delta x). $$(6)

The effective wavenumber k′ is given by a sine. Only small wavenumbers are treated in good agreement with the analytical solution. The error leads to wrong phase velocities and thus to numerical dispersion. The dispersion error is given by the relation k′/k modifying the actual phase velocity vph′=k'/k⋅vph$ {v}_{\mathrm{ph}}^\mathrm{\prime}=k\prime/k\cdot {v}_{\mathrm{ph}}$ [29]. An extension to more dimensions is presented in [29]. An extensive discussion is given in [54].

Higher-order schemes requires increasing the number of points in the stencil. In general, it improves the dispersion behavior. Instead of using additional points to increase the order of convergence, they can also be used to improve the dispersion behavior directly. The best-known family of stencils is provided by the dispersion-relation-preserving (DRP) finite-difference scheme [30] (Eq. (A4)). Another possibility to increase the order of consistency and to improve the dispersion behavior while keeping the stencil width small is to use compact schemes. The general idea is to include spatial derivatives of neighboring points on the left hand side of the deviation stencil, such as in equation (5), e.g.,

β+1∂xqm+1+∂xqm+β-1∂xqm-1=…$$ {\beta }_{+1}{\partial }_x{q}_{m+1}+{\partial }_x{q}_m+{\beta }_{-1}{\partial }_x{q}_{m-1}=\dots $$(7)

where β±1 are additional real valued weights, see [29] for further information. However, this leads to an implicit definition of the derivative. In matrix-vector notation, the derivative is expressed as Bnm∂xqm = Anmqm, with n, m = [1, …, M] and M the number of grid points. The matrix A acts on the discrete field values as before while B, typically being tri- or pentadiagonal, contains suitable weights.

A straightforward approach to applying the compact or implicit derivative is to invert the matrix B and compute the derivative by

∂xq=B-1Aq=Dq.$$ {\partial }_xq={B}^{-1}{Aq}={Dq}. $$(8)

In general, the inverse of B is dense, nearly full, so evaluating equation (8) is computationally expensive. Instead, the linear system Bq′ = Aq = b is solved within each application of the derivative using an appropriate method, such as a lower-upper (LU) decomposition [55] or the Thomas algorithm [56]. Optimized compact schemes are available [29], see Figure 1a (O6 imp LELE) (Eq. (A7)).

thumbnail Figure 1

(a) Harmonic analysis of the dispersion relation c′/c and (b) the error |k′Δx − kΔx|real of different finite-difference-stencils for first order partial differential equations. (c) CPU-time with respect to the error. For the analysis, the vector length was chosen to match PPW resolution and dispersion error from (b). The short forms OX denote the order, exp or imp, explicit or implicit / compact derivatives (Eq. (8)). The error of 10−3 is highlighted by a black dotted line as it is used for discussion the text.

In room acoustic applications, the optimal spatial derivative is defined by the following requirements: spectral-like dispersion behavior using only a few points, and good performance in terms of computation time.

For an analysis of the computational cost, the difference between analytical and numerical dispersion behavior is translated into the required number of points per wavelength (PPW), see Figure 1b for different derivative stencils. Figure 1c measures the CPU-time required for each derivative and compares it to the dispersion behavior error. Note, that the results depend on the computer architecture used. Here, the analysis was performed on a ×86 computer.

It is found that explicit derivatives have the best performance, although they require more PPW to achieve comparable quality. Thus, explicit derivatives are advantageous when the derivatives dominate the total CPU time, as it is the case with direct computation. If the numerical simulation is dominated by other means, this may no longer be the case.

For simulations where the computation time is heavily influenced by I/O-operations, such as adjoint-based gradient simulations which require the storage of and access to the entire data of the computational field [57], a reduced number of grid points is desired. Then, implicit schemes are favorable, especially the optimized scheme [29].

The validity of the previous analyses is formally restricted to 1D problems, and a direct extension to 2D or 3D is challenging. Since a detailed discussion exceeds the scope of this manuscript, the reader is referred to [29, 58, 59] for further details. In addition, the dispersion relation is not solely influenced by spatial discretization. Temporal discretization also plays a role in determining both dispersion and dissipation.

2.3 Temporal discretization

In general, two different types of temporal integration schemes can be distinguished: explicit and implicit schemes. Explicit schemes are characterized by the fact that only already known (state) information [tn, tn−1, …] is used to compute the next time step tn+1, whereas implicit schemes include the future state. For purely linear problems, this requires solving a linear system. For nonlinear problems, such as fluid dynamics, iterative methods, such as Newton methods, are required, which usually require more computational resources. However, implicit time integration schemes have superior stability properties [60]. Regardless, implicit methods typically allow a coarser time increment compared to explicit methods before the resulting errors become unacceptable. Therefore, implicit methods are advantageous when the short-time dynamics of an analyzed configuration are not of interest. Thus, for acoustic simulations, explicit methods will provide the best performance.

For explicit schemes in particular, a distinction is made between one-step and multi-step methods. The former needs only the current state information at time tn to compute the next time step tn+1, while multi-step methods use several time steps in the past. Both approaches differ in their stability behavior, which leads to different maximum time increments [61]. Since one-step methods compute substeps between tn and tn+1 they need to call multiple times the right-hand-side (rhs) of the governing equations, see Figure 2b.

thumbnail Figure 2

Analysis of selected explicit time integration schemes. The short forms OX denote the order, RK the Runge-Kutta schemes and AB the Adams-Bashfourth schemes. (a) Stability regions for selected explicit time integration schemes. (b) Maximum stable time increment for different methods i normalized with respect to the largest value in blue. Maximum time increment additionally normalized with respect to the number of rhs-calls (evaluations of right-hand-sides of the governing equations) needed for time integration indicating the numerical efficiency in orange.

Various stability analyses are known, e.g., the classical von Neumann analysis [53] and the matrix method [62]. Here, absolute stability is discussed in terms of a linear model problem without any filtering

 q̇(t)=λ̃q(t),q∈R.$$ \enspace \begin{array}{cc}\dot{q}(t)=\stackrel{\tilde }{\lambda }q(t),& q\in \mathbb{R}.\end{array} $$(9)

(Absolute) stability is found in regions where

|qn+1||qn|≤1$$ \frac{|{q}^{n+1}|}{|{q}^n|}\le 1 $$(10)

holds. By analyzing different schemes in terms of this relation and plotting the contour |qn+1|/|qn| = 1, which defines the stability limit, the schemes can be evaluated and an optimal scheme can be chosen.

Since the fastest dynamics, e.g., sound waves are of interest in room acoustic simulations, the focus is on explicit time integration methods. The requirements are a stable scheme for the intended time step, or at least a sufficient stability region, with as few as possible calls of the governing equations/rhs.

To analyze the numerical efficiency a representative rhs, mimicking transport and viscous effects using a second-order symmetric spatial derivative matrix D, is analyzed. The eigenvalues of the effective 1D-operator

rhs=-λD+μ(D⋅D),λ,μ∈R;D∈Rn×n$$ \begin{array}{cc}\mathrm{rhs}=-{\lambda D}+\mu \left(D\cdot D\right),& \lambda,\mu \in \mathbb{R};D\in {\mathbb{R}}^{n\times n}\end{array} $$(11)

are scaled by increasing time step Δt as long as the resulting spectrum is within the stability region to determine the maximum stable time step, see Figure 2a. The parameter values are set to λ = 343 m/s and μ = 0, corresponding to sound propagation in air without absorption effects. If viscous effects are to be taken into account, μ must be selected accordingly in the model equation for the analysis.

Figure 2b shows that the Adams-Bashfourth 3rd order scheme (O3 exp AB) and the Runge-Kutta 4th order scheme (O4 exp RK) have the best performance. However, the latter is the method of choice because it allows a significantly larger time step while providing a higher order of consistency. The 4th-order Runge-Kutta is available in several realizations, e.g., a low-storage implementation [32]. The time discretization impacts the disperion relation. For example the standard Runge-Kutta method of fourth order is known to be dissipative. A detailed analysis exceeds the focus of this contribution. However, it should be noted that dispersion- and dissipation-optimized methods are available for the previously examined procedures [30, 33, 63, 64], and their application is recommended. Furthermore, it must be noted that the chosen spatial discretization influences time integration, and both should not be considered separately. However, a detailed analysis of the system under investigation would exceed the scope of this manuscript, and is therefore omitted for the sake of brevity.

2.4 Filtering and damping

Acceptable dispersion behavior is limited to high resolutions, i.e., points per wavelength, or low wave numbers, while high wave numbers are inadequately treated. In particular, for wave numbers kΔx ≳ 2, parasite waves occur that contaminate the simulation. This effect cannot be avoided by increasing the spatial discretization, since the parasitic waves are eigensolutions of the finite difference method.

In order to obtain time-stable acoustic simulations artificial damping [30] or filtering is necessary, since the Euler equations do not provide any inherent physical mechanism to avoid the appearance or damping of parasitic waves. Note, that this is only required for the nonlinear equations. In the following, filtering is discussed.

Filters directly modify the considered field, in the sense that the unwanted wavelengths are removed rather than attenuated. The explicit filter approach is given by

q̅=∑ngnqm+n,$$ \bar{q}=\sum_n {g}_n{q}_{m+n}, $$(12)

where qm denotes the unfiltered field at grid point m with m = [1, …, M] and M the number of grid points and q̅$ \bar{q}$ the filtered one. The weights gn define the filter stencil, which is usually chosen to be symmetric g−n = gn with n = [−N, …, N] and N the filter stencil width.

The stencil coefficients are set so that the highest grid mode kΔx = π is removed. Again, there are several stencils available that differ in number of points and order of convergence. In general, the filter approach is less dispersive compared to a selective damping approach and easier to control because the field is filtered outside of the rhs.

Similar to the implicit or compact spatial derivatives, implicit filters also exist [35]. The approach in equation (12) is extended by

∑ohoq̅m+o=∑ngnqm+n$$ \sum_o {h}_o{\bar{q}}_{m+o}=\sum_n {g}_n{q}_{m+n} $$(13)

with ho as additional coefficients, o = [−O, … O] and O being the stencil width. Again, the implicit approach can be written in matrix-vector notation Hq̅=Gq$ H\bar{q}={Gq}$. Different sets of coefficients are given in [35]. This formulation requires solving a linear system for each vector to be filtered. However, implicit filtering has been found to be less dispersive than explicit filtering with the same stencil width [35].

The choice of the filter depends on the dispersion characteristic given by the spatial discretization. The frequency responses of explicit and implicit filters are shown in Figure 3. The coefficients ho and gn of equation (13) are taken from [35]. A low order filter has a lower cut-off frequency than higher order filters. Implicit filters have a higher cut-off frequency than explicit filters with the same stencil width. In general, if a low order spatial discretization is used, filters with a lower cut-off frequency must be incorporated. On the other hand, high order spatial discretization requires high order filtering. In particular, when implicit spatial derivatives are used and the framework for solving the resulting linear system is already in place, implicit filters are the method of choice.

thumbnail Figure 3

Frequency response of different explicit (exp) and implicit (imp) filters. Cut-off frequencies are at 0.5$ \sqrt{0.5}$ marked by the black dashed line.

2.5 Boundary conditions/layer

Boundary conditions are decisive for the quality of acoustic simulations [1]. Two types of boundary conditions are of particular interest: non-reflective and impedance boundary conditions, with the latter especially involving porous, acoustically damping wall elements.

Non-reflecting boundary conditions are not highly relevant for interior room acoustical problems but become relevant in architectural acoustics. They can be realized by different approaches such as the characteristic boundary condition (CB) [36, 65] or the perfectly matched layer (PML) technique [37, 66].

Impedance boundary methods are subject of current research in the context of finite difference methods [67, 68]. A recent trend is to use a volume penalization approach to model complex impedances [46, 48].

In the following, CB and a seleceted immersed boundary approach are discussed in more detail.

2.5.1 Characteristic boundary conditions

The method of characteristics includes a generalized separation of variables. The basic idea is to decompose the flow field at the boundary into characteristic waves or Riemann-invariants respectively and to set incoming waves to a certain value. The requirements are quasi-linearity and hyperbolicity of the governing partial differential equations, which are satisfied by the Euler equations [69].

For the sake of clarity the approach is first considered in one dimension. Writing the Euler equations in 1-D in a quasi-linear form results in ∂tq + X∂xq = 0, with X as governing operator. Since X can be diagonalized, it can be rewritten with a local basis of eigenvectors L as

X=LλL-1.$$ X={L\lambda }{L}^{-1}. $$(14)

The eigenvalues of X are λs = u and λ±= u ± c, where λs is related to the so-called entropy wave and λ± describes the acoustic waves. The characteristic variable is defined as w = L−1q. By choosing a suitable normalization of the left ei and the right ei eigenvectors

X=(   ese+e-   )(λs   λ+   λ-)( es  e+  e- )$$ X=\left(\begin{array}{lll}\enspace & \enspace & \enspace \\ {e}_s& {e}_{+}& {e}_{-}\\ \enspace & \enspace & \enspace \end{array}\right)\left(\begin{array}{lll}{\lambda }_s& \enspace & \enspace \\ \enspace & {\lambda }_{+}& \enspace \\ \enspace & \enspace & {\lambda }_{-}\end{array}\right)\left(\begin{array}{lll}\enspace & {e}_s& \enspace \\ \enspace & {e}^{+}& \enspace \\ \enspace & {e}_{-}& \enspace \end{array}\right) $$(15)

eiej = Δij (bi-orthonormality) holds.

For a non-reflecting boundary condition the following decomposition is applied: The actual state at the boundary is written in terms of a disturbance Δq = q − q0 with a suitable reference state q0 (e.g. initial condition) and decomposed by Δq = Δrses + Δr+e+ + Δr−e−. Using the bi-orthonormal relation one finds Δrs = esΔq, Δr+ = e+Δq, Δr− = e−Δq. For non-reflecting boundary conditions, assuming a subsonic flow in positive direction, the ingoing waves must be set to zero, e.g. at the left boundary, λs=!0$ {\lambda }_s\stackrel{!}{=}0$ and λ+=!0$ {\lambda }_{+}\stackrel{!}{=}0$ and at the right boundary condition λ-=!0$ {\lambda }_{-}\stackrel{!}{=}0$. Thus, the disturbances become at the left boundary Δq̃=0+0+Δr-e-$ \Delta \mathop{q}\limits^\tilde=0+0+\Delta {r}_{-}{e}_{-}$ and Δq̃=Δrses+Δr+e++0$ \Delta \mathop{q}\limits^\tilde=\Delta {r}_s{e}_s+\Delta {r}_{+}{e}_{+}+0$ at the right one. The boundary value becomes finally q̃=q0+Δq̃$ \mathop{q}\limits^\tilde={q}_0+\Delta \mathop{q}\limits^\tilde$.

In two- and three-dimensional systems, the waves are separated by introducing a direction kxi$ {k}_{{x}_i}$, which is typically chosen orthogonal to the boundary surface or in the direction of the base flow [70]. The eigenvalues in the more-dimensional case are λ = u (multiple eigenvalue), λ+ = u + c and λ− = u − c with u=kxiui$ u={k}_{{x}_i}{u}_i$.

To improve the quality of characteristic boundary conditions they can be applied not only to the boundary layer but also to neighbour layers via the zonal approach [71].

To reduce remaining spurious reflections of the characteristic boundary conditions, a sponge layer can be employed. It essentially acts as an absorbing condition, designed to reduce outgoing waves and prevent them from reflecting back into the computational domain. By doing so, the sponge layer helps to simulate the behavior of an open or infinite domain. A detailed discussion is provided in [38]. Guidance on its use with additional grid stretching can be found in [72].

2.5.2 Volume penalization

To facilitate the representation of walls with complex impedance, here modeled with porous (acoustic damping) materials, an immersed volume approach can be utilized. A pioneering work on the immersed boundary (IB) methods was done by Peskin [41]. An IB method for compressible flows was presented by de Palma et al.[73]. Liu and Vasilyev [74] included a linear friction relative to the object (Darcy’s law) and an effective volume fraction remaining for the fluid. Based on the work of Kevlahan et al. [75] for the shallow water equation, Reiss [45] showed the fluid dynamical and first acoustic properties of the method. The volume penalized NLEE are

φ∂t(ϱϱujpγ-1)+∂xi(φϱuiφϱuiujφuipγγ-1)+φ∂xi(0pδij0)-ui∂xi(00φp)=(0φχ(ujt-uj)0).$$ \begin{array}{c}\phi {\partial }_t\left(\begin{array}{l}\varrho \\ \varrho {u}_j\\ \frac{p}{\gamma -1}\\ \end{array}\right)+{\partial }_{{x}_i}\left(\begin{array}{l}{\phi \varrho }{u}_i\\ {\phi \varrho }{u}_i{u}_j\\ \phi \frac{{u}_i{p\gamma }}{\gamma -1}\end{array}\right)+\phi {\partial }_{{x}_i}\left(\begin{array}{l}0\\ p{\delta }_{{ij}}\\ 0\end{array}\right)\\ -{u}_i{\partial }_{{x}_i}\left(\begin{array}{l}0\\ 0\\ {\phi p}\end{array}\right)=\left(\begin{array}{l}0\\ {\phi \chi }({u}_j^{\mathrm{t}}-{u}_j)\\ 0\end{array}\right).\end{array} $$(16)

The dimensionless effective volume φ represents a volume fraction φ(xi) = Vfluid/Vtotal due to the presence of a porous medium, ranging from 0 to 1. The effective volume φ corresponds to the porosity. When φ = 0, the Euler equations degenerate, necessitating a small but finite value for φ in simulations to avoid this. Conversely, with φ = 1, the unmodified Euler equations are restored. Negative or non-physical values (>1) are excluded. A time-dependent φ is addressed in [45, 76] and not considered here. Further discussion on φ can be found in [45, 76].

In addition, the momentum equations expand via a penalization φχ(ujt-uj)$ {\phi \chi }({u}_j^{\mathrm{t}}-{u}_j)$, with χ(xi) controlling spatial location and amplitude which is different from zero when (acoustic) material is presented. This is referred to as the Darcy term, damping the (particle) velocity due to the presence of a porous material if ujt=0$ {u}_j^{\mathrm{t}}=0$ as done here. χ’s dimension is Ns/m4, with χ → ∞ causing the Euler equations to degenerate.

Both parameters of the volume penalization (φ and χ) are typically tabulated for acoustically damping materials or can be inquired about from the respective manufacturers of the acoustic damping materials. φ corresponds to the porosity, and χ corresponds to the length-specific flow resistance.

Solid/semi-permeable boundaries can be modeled using large χ values yielding larger negative rhs operator eigenvalues, affecting simulation stability or necessitating smaller time steps for explicit time marching. However, also a vanishing φ (e.g., 10−6 to 10−4) models rigid walls. A detailed description and examples mimicking rigid and impedance boundary conditions are given in [45, 46].

Please note, the volume penalization approach can be interpreted as a simplified variant of the equivalent fluid method in porous materials.

Equation (16) can be simplified in the same way as in Section 2.1 to linear acoustic equations:

∂t(p′u′j)+ϱ0c02φ∂xj(φu′j0)+1ϱ0∂xj(0p′)=(0-1ϱ0χu'j).$$ {\partial }_t\left(\begin{array}{l}p\mathrm{\prime}\\ {u\mathrm{\prime}_j\end{array}\right)+\frac{{\varrho }_0{c}_0^2}{\phi }{\partial }_{{x}_j}\left(\begin{array}{l}\phi {u\mathrm{\prime}_j\\ 0\end{array}\right)+\frac{1}{{\varrho }_0}{\partial }_{{x}_j}\left(\begin{array}{l}0\\ p\mathrm{\prime}\end{array}\right)=\left(\begin{array}{l}0\\ -\frac{1}{{\varrho }_0}\chi {{u\prime}_j\end{array}\right). $$(17)

A formulation of the linear acoustic equations using only the Darcy term with applications to room acoustic problems is given in [47, 48].

2.6 Classification and comparison

The Finite-Difference Time-Domain (FDTD) method discretizes the entire problem domain, making it particularly suitable for well-defined, bounded regions and location-dependent material properties (e.g., volume penalization). While it can be computationally intensive for large domains, as it scales spatially with the number of grid points (N) in three dimensions (N3) and temporally with the number of time steps, it remains very efficient for time-dependent problems. This is because the discretization is sparse, avoiding the need to solve large linear systems. In addition, its memory requirement is low. The FDTD approach is highly parallelizable, making it particularly efficient on parallel computing platforms. These properties make it a very suitable approach for large-scale acoustic simulations in comparison to other commonly used approaches like BEM, which reduces the problem by one dimension, resulting in fewer equations but necessitating the solution of a dense system of equations, which requires substantial memory and computation time. The BEM method is efficient for problems with a small surface-to-volume ratio but is less scalable and only marginally parallelizable. Furthermore, it is not suitable for non-linear problems [77].

3 3D cubic room with rigid boundaries

To evaluate the spatial derivatives and filter schemes a 3D cubic room of 1 × 1 × 1 m with rigid boundaries is considered. This common test case is of interest because an analytical solution exists ([21], Eqs. 32–37]. The spatial resolution is set to 50 points distributed equidistantly in each direction, i.e., 1.25 · 105 DOFs and Δx = 0.02 m. The time-stepping method is an explicit 4th-order Runge-Kutta scheme. The sampling frequency is set to fs = 50 kHz and the governing equations are the NLEE. Note that for this test case, the LAE are exactly equivalent if only small acoustic amplitudes are considered. The initial conditions at rest result in a speed of sound of c0 = 343 m/s resulting from ϱ0 = 1.2057 kg/m3 and p0 = 101325 Pa. The domain is excited by a Gaussian pulse

p′(xj,t=0)=exp(-∑j=1,2,3(xj-xj,src)2σ2)$$ p\mathrm{\prime}({x}_j,t=0)=\mathrm{exp}\left(-\frac{\sum_{j=\mathrm{1,2},3} ({x}_j-{x}_{j,\mathrm{src}}{)}^2}{{\sigma }^2}\right) $$(18)

with a spatial variance of σ = 8Δx in the first time step t = 0. In the case of NLEE the adiabatic relation is used to initialize the density

ϱ′(xj,t=0)=ϱ0(p′(xj,t=0)p0)1κ,$$ \varrho \mathrm{\prime}({x}_j,t=0)={\varrho }_0{\left(\frac{p\mathrm{\prime}({x}_j,t=0)}{{p}_0}\right)}^{\frac{1}{\kappa }}, $$(19)

where the adiabatic exponent in air is κ = 1.4. For all spatial schemes, the points near the boundary are treated with second order accuracy. The boundaries are treated as fully reflecting using the characteristic approach. All 3D FDTD simulations in this paper were performed using the CRUNA code [78].

3.1 Spatial derivatives

Figure 4 presents the linear frequency response for the 2nd- and 4th-order explicit and 4th-order implicit spatial derivatives. For the case shown at the bottom the simulation is initialized with a high amplitude leading to non-linear effects.

thumbnail Figure 4

Frequency response of 2nd and 4th order explicit and 4th-order implicit spatial derivatives. The receiver is located at xr = [0.7, 0.8, 0.7]T m and the source is located at xs = [0.2, 0.4, 0.3]T m. The analytical solution is shown in orange dashed and the numerical solution in blue. The gray area marks regions with a theoretical error (k’Δx − kΔx)real > 10−3. The lower non-linear case shows results for high amplitudes, at 0.75 m distance from the source the peak magnitude of the blue line is at ≈146 dBSPL and in light blue dashed at ≈152 dBSPL.

Assuming an error of 10−3 in the dispersion relation as the maximum acceptable, we expect good results for the 2nd-order explicit for frequencies discretized with up to 32 PPW, the 4th-order explicit 14 PPW, and the 4th-order implicit 10 PPW, see Section 2.3. This holds for the 2nd-order explicit derivative. Both 4th-order derivatives show good agreement with the acoustic modal frequencies but differ in the amplitude already below 20 PPW and −20 dB. The frequency shift of the individual modes up to 1.4 kHz is analyzed in Figure 5. The deviation shows good agreement with the theoretical error (k′Δx − kΔx)real > 10−3. For all schemes, the amplitude error increases with the frequency. Analyzing the non-linear case, deviations from the analytical solution are observed at lower frequencies with increasing amplitude, corresponding to [79]. At high frequencies, additional room modes appear. Therefore, it is necessary to include the nonlinearities when simulating high sound pressure levels, as is usually the case in sound reinforcement scenarios in rooms.

thumbnail Figure 5

Errors in simulated modal frequencies for the cube-shaped room with rigid boundaries.

A simulation of the room using the volume penalization approach can be found in [46].

3.2 Filtering schemes

Filtering is needed if the nonlinear governing equations are applied. To analyze the effect of filtering the same setup is used with the implicit 4th-order spatial derivative scheme. Results for an explicit 4th-order filter scheme and implicit 4th- and 8th-order filter schemes are shown in Figure 6.

thumbnail Figure 6

Frequency response of 4th-order explicit, 4th- and 8th-order implicit filter schemes. The receiver is at xr = [0.7, 0.8, 0.7]T m and the source is at xs = [0.2, 0.4, 0.3]T m. The analytical solution is shown in orange dashed and the numerical solution in blue.

According to Figure 3 no significant effects are expected for frequencies discretized with more than 16 points. However, there are amplitude deviations, even before the expected minimum resolution. This is true for all filter stencils examined. The frequencies of the spatial impulse response are not affected.

Therefore, it is advisable to choose high order schemes with an equal or higher order compared to the spatial derivative. Otherwise, correctly simulated frequencies may be filtered.

4 Sound reflections with different boundary conditions

To demonstrate the suitability of the NLEE in combination with the volume penalization for a real room acoustic environment, FDTD simulations of a sound reflection on surfaces with different boundary conditions are compared with experimental results from the BRAS database [80] as well as with results from a ray tracing simulation [18]. The FDTD simulation is performed with a 6th-order implicit spatial derivative scheme, a standard 4th-order Runge-Kutta scheme, and a 6th-order spatial filtering scheme. The spatial discretization is realized by 780 × 472 × 472 equidistant points with a spacing of Δxi = 0.013 m. Hence, the physical size of the computational domain is 10.14 × 6.14 × 6.14 m. The time step is set to Δt = 2.0833 ⋅ 10−5 s, which results in a CFL condition of ≈0.56. The initial conditions give a speed of sound of c0 = 343 m/s and the same Gaussian pulse as in Section 3 is used. The boundaries of the computational domain are treated as non-reflecting using the characteristic approach to be aligned with the experimental setup of the BRAS scenario conducted in an semi-anechoic room. Additionally, a quadratic sponge layer with a length of 0.3 m was applied. Immersed boundaries are defined by tanh-functions similar to [46] and used to approximate the rigid floor and the diffuser by setting φ = 2.5 ⋅ 10−4 and χ = 0 Ns/m4 in equation (16). Please note, that for a fully reflecting surface no Darcy-term is needed, if the acoustic boundary layer is not of interest, see [46]. The absorber was approximated by χ = 6000 Ns/m4 and φ = 1 since the fluid resistivity of the Rockfon absorber (Sonar G) in the measurement is ≈6000 Pa s/m2. The setup is shown in Figure 7a–7c. The experimental data and the simulation results were normalized to 0 dB in case (Fig. 7a). The same normalization value was applied to cases (Figs. 7b and 7c).

thumbnail Figure 7

Results of the FDTD simulation (blue) compared to experimental results of the BRAS database, scene RS1 (red dashed line) and simulation results of the BRAS E solver (yellow). In case (a), the blue dashed vertical lines show the analytical comb filter frequencies and the gray dashed vertical lines show the number of PPWs. (a) Rigid, (b) Absorber, (c) Diffuser.

The FDTD simulation in case (Fig. 7a) shows a perfect agreement with the analytical comb filter frequencies and in all cases a very good agreement with the experimental data in the range between 200 Hz and 2.2 kHz. Furthermore, the simulation results of the FDTD solver show a better agreement with the experimental data than the ray tracing solver BRAS E, especially in case (Fig. 7c), proving the necessity of wave-based methods when analyzing wave effects such as diffraction. The results indicate that the approach is suitable for accurately reproducing the investigated boundary conditions or room acoustic configurations. A comprehensive validation of the approach is presented in [46].

5 Non-linear behavior in a sound reinforcement scenario

This case demonstrates the difference between linear and nonlinear equations in a sound reinforcement scenario. It simulates a shoebox room with a physical size of 15 × 8 × 6 m. A sketch of the room is shown in Figure 8. The walls are assumed to be rigid with φ = 10−4 and χ = 0 Ns/m4. Two centered absorber plates with a size of 11 × 4 × 0.05 m are attached to the side walls and one absorber plate with a size of 6 × 4 × 0.05 m is attached to the back wall. They have a fluid resistivity of χ = 6000 Ns/m4 and an effective volume of φ = 1. Five sound sources were used at typical positions of a line source array computed by PALC [81], see Table 1. A peak magnitude of Lp̂,1m=[147,130,96]$ {L}_{\widehat{p},1\mathrm{m}}=[\mathrm{147,130,96}]$ dBSPL was measured at a distance of 1 m from the center of the array to mimic a powerful but not untypical line source array [82]. The same numerical setup was used as in Section 4. Figure 9 shows third octave smoothed differences of the frequency responses of the linear and the nonlinear equations at the receiver position x = [8.0; 5.0; 1.5] m, corresponding to a typical audience distance. The results of the 147 dBSPL case show that the magnitude difference reaches 8 dBSPL before 2500 Hz. Also the 130 dBSPL case exhibits a perceptually relevant difference of 3 dBSPL before 3000 Hz. At lower amplitudes like 96 dBSPL the differences disappear. Similar to Section 3.1, the results correspond to [79]. It can be concluded that nonlinearities have a remarkable influence at high SPLs, e.g., at typical audience positions in sound reinforcement scenarios. Similar high SPLs can also be reached by acoustic instruments [83].

thumbnail Figure 8

Sketch of the shoebox room. Edges of the room are the black lines, gray areas are porous absorber. Red dots are the source positions, see Table 1 and the blue dot marks the receiver position at x = [8.0; 5.0; 1.5] m.

thumbnail Figure 9

Differences of the frequency response of linear and nonlinear equations at the measurement position x = [8.0; 5.0; 1.5] m.

Table 1

Positions of the sound sources.

6 Discussion

This paper presents in detail FDTD simulation techniques for first-order Euler partial differential equations used in a room acoustic context. This includes spatial and temporal discretization, filtering schemes, and boundary conditions. Depending on the requirements of the simulation, e.g., the shortest wavelength, recommendations for the optimal choice of simulation techniques are given.

It has been shown that explicit spatial discretization schemes give the best performance for direct simulations. Otherwise, implicit schemes are the method of choice, e.g., for adjoint-based simulations, as they allow a lower spatial resolution and thus a larger time increment. The optimal temporal discretization is typically achieved through one-step methods, such as the explicit 4th-order Runge-Kutta scheme. Only when a very high temporal discretization is chosen, linear multi-step methods become advantageous, e.g. the 3rd-order Adams-Bashfourth scheme. The choice of an appropriate spatial filter order depends on the spatial discretization scheme. In general, it is recommended that the order of the filter be higher than the order of the spatial discretization scheme.

The dispersion errors of three spatial discretization and filtering schemes were investigated in a 3D cubic room with rigid boundaries. The theoretical results were confirmed in the linear case. The nonlinear case shows increasing deviations with amplitude and frequency. A set of the aforementioned techniques was used to simulate the BRAS scene 1 (SR1). Rigid and absorbing boundaries, as well as a rigid diffuser, were approximated by a volume penalization method. The results show a very good agreement with the experimental data. Finally, an artificial room was employed to demonstrate the difference between linear and nonlinear equations. When typical high sound pressure levels of line source arrays in sound reinforcement applications are reached, linear equations are unable to accurately approximate the solution correctly with increasing frequency.

In summary, it has been shown that the proposed methods, mostly known from computational aeroacoustics, are applicable to room acoustic problems. Non-linear effects must be taken into account when high sound pressure levels are reached, such as in sound reinforcement scenarios. In particular, the volume penalization method is able to mimic impedance boundaries without the need for grid adjustment or other numerical effort. In the future, additional scenes from the BRAS dataset will be used to further validate the method. In addition, it is necessary to conduct measurements of the setup as described in Section 5 in order to verify the results.

Acknowledgments

The authors gratefully acknowledge financial support funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 532453154. We acknowledge support by the Open Access Publication Fund of TU Berlin.

Conflicts of interest

The authors declare no conflict of interest.

Data availability statement

Data are available on request from the authors. All FDTD simulations were performed using the CRUNA code [78] (https://github.com/cruna-toolkit/).

Appendix

Spatial derivative schemes

Second order explicit

∂ x q m = q m + 1 - q m - 1 2 Δ x $$ {\partial }_x{q}_m=\frac{{q}_{m+1}-{q}_{m-1}}{2\Delta x} $$(A1)

Fourth order explicit

∂xqm=∑i=12αi(qm+i-qm-i)12Δx$$ {\partial }_x{q}_m=\sum_{i=1}^2 \frac{{\alpha }_i({q}_{m+i}-{q}_{m-i})}{12\Delta x} $$(A2)

with α1 = 8 and α2 = −1.

Sixth order explicit

∂xqm=∑i=13αi(qm+i-qm-i)60Δx$$ {\partial }_x{q}_m=\sum_{i=1}^3 \frac{{\alpha }_i({q}_{m+i}-{q}_{m-i})}{60\Delta x} $$(A3)

with α1 = 45, α2 = −9 and α3 = 1.

Fourth order DRP explicit

∂xqm=∑i=13αi(qm+i-qm-i)Δx$$ {\partial }_x{q}_m=\sum_{i=1}^3 \frac{{\alpha }_i({q}_{m+i}-{q}_{m-i})}{\Delta x} $$(A4)

with α1 = 0.79926643, α2 = −0.18941314 and α3 = 0.02651995.

Fourth order implicit

β(∂xqm+1+∂xqm-1)+∂xqm=α(qm+1-qm-1)2Δx$$ \beta ({\partial }_x{q}_{m+1}+{\partial }_x{q}_{m-1})+{\partial }_x{q}_m=\frac{\alpha ({q}_{m+1}-{q}_{m-1})}{2\Delta x} $$(A5)

with α = 1.5 and β = 0.25.

Sixth order implicit

β(∂xqm+1+∂xqm-1)+∂xqm=∑i=12αi(qm+i-qm-i)Δx$$ \beta ({\partial }_x{q}_{m+1}+{\partial }_x{q}_{m-1})+{\partial }_x{q}_m=\sum_{i=1}^2 \frac{{\alpha }_i({q}_{m+i}-{q}_{m-i})}{\Delta x} $$(A6)

with α1 = 7/9, α2 = 1/36 and β = 1/3.

Optimized sixth order implicit (Lele)

∑l=12βl(∂xqm+l+∂xqm-l)+∂xqm=∑i=13αi(qm+i-qm-i)2iΔx$$ \sum_{l=1}^2 {\beta }_l({\partial }_x{q}_{m+l}+{\partial }_x{q}_{m-l})+{\partial }_x{q}_m=\sum_{i=1}^3 \frac{{\alpha }_i({q}_{m+i}-{q}_{m-i})}{2i\Delta x} $$(A7)

with α1 = 1.3025166, α2 = 0.99355, α3 = 0.03750245, β1 = 0.5771439 and β2 = 0.0896406.

References

  1. M. Vorländer: Computer simulations in room acoustics: Concepts and uncertainties, Journal of the Acoustical Society of America 133, 3 (2013) 1203–1213. [CrossRef] [PubMed] [Google Scholar]
  2. L. Savioja, U.P. Svensson: Overview of geometrical room acoustic modeling techniques, Journal of the Acoustical Society of America 138, 2 (2015) 708–730. [Google Scholar]
  3. D. Takeuchi, K. Yatabe, Y. Oikawa: Source directivity approximation for finite difference time-domain simulation by estimating initial value, Journal of the Acoustical Society of America 145, 4 (2019) 2638–2649. [CrossRef] [PubMed] [Google Scholar]
  4. D. Botteldooren: Finite-difference time-domain simulation of low-frequency room acoustic problems, Journal of the Acoustical Society of America 98, 6 (1995) 3302–3308. [CrossRef] [Google Scholar]
  5. K. Kowalczyk, M. van Walstijn: Room acoustics simulation using 3D compact explicit FDTD schemes, IEEE Transactions on Audio, Speech, and Language Processing 19, 1 (2011) 34–46. [CrossRef] [Google Scholar]
  6. A. Craggs: A finite element method for free vibration of air in ducts and rooms with absorbing walls, Journal of Sound and Vibration 173, 4 (1994) 568–576. [CrossRef] [Google Scholar]
  7. A.G. Prinn: A review of finite element methods for room acoustics, Acoustics 5, 2 (2023) 367–395. [CrossRef] [Google Scholar]
  8. A. Hargreaves, T.J. Cox: A transient boundary element method model of Schroeder diffuser scattering using well mouth impedance, Journal of the Acoustical Society of America 124, 5 (2008) 2942–2951. [CrossRef] [PubMed] [Google Scholar]
  9. N.A. Gumerov, R. Duraiswami: Fast multipole accelerated boundary element methods for room acoustics, Journal of the Acoustical Society of America 150, 3 (2021) 1707–1720. [CrossRef] [PubMed] [Google Scholar]
  10. S. Bilbao: Modeling of complex geometries and boundary conditions in finite difference / finite volume time domain room acoustics simulation, IEEE Transactions on Audio, Speech, and Language Processing 21, 7 (2013) 1524–1533. [CrossRef] [Google Scholar]
  11. S. Bilbao, B. Hamilton: Wave-based room acoustics simulation: explicit / implicit finite volume modeling of viscothermal losses and frequencydependent boundaries, Journal of the Audio Engineering Society 65, 1/2 (2017) 78–89. [CrossRef] [Google Scholar]
  12. H. Wang, I. Sihar, R. Pagan Munoz, M. Hornikx: Room acoustics modelling in the time-domain with the nodal discontinuous Galerkin method, Journal of the Acoustical Society of America 145, 4 (2019) 2605–2663. [Google Scholar]
  13. F. Pind, A.P. Engsig-Karup, C.H. Jeong, J.S. Hesthaven, J.S. Mejling, J. Strømann-Andersen: Time domain room acoustic simulations using the spectral element method, Journal of the Acoustical Society of America 145, 6 (2019) 3299–3310. [CrossRef] [PubMed] [Google Scholar]
  14. F. Pind, C.H. Jeong, A.P. Engsig-Karup, J.S. Hesthaven, J. Strømann-Andersen: Timedomain room acoustic simulations with extendedreacting porous absorbers using the discontinuous Galerkin method, Journal of the Acoustical Society of America 148, 5 (2020) 2851–2863. [CrossRef] [PubMed] [Google Scholar]
  15. B. Hamilton, S. Bilbao: Time-domain modeling of wave-based room acoustics including viscothermal and relaxation effects in air, JASA Express Letters 1, 9 (2021) 092401. [CrossRef] [PubMed] [Google Scholar]
  16. G. Fratoni, B. Hamilton, D. D’Orazio: Feasibility of a finite-difference time-domain model in large-scale acoustic simulations, Journal of the Acoustical Society of America 152, 1 (2022) 330–341. [CrossRef] [PubMed] [Google Scholar]
  17. Y. Li, J. Meyer, T. Lokki, J. Cuenca, O. Atak, W. Desmet: Benchmarking of finite difference time-domain method and fast multipole boundary element method for room acoustics, Applied Acoustics 191 (2022) 1–13. [Google Scholar]
  18. F. Brinkmann, L. Aspöck, D. Ackermann, S. Lepa, M. Vorländer, S. Weinzierl: A round robin on room acoustical simulation and auralization, Journal of the Acoustical Society of America 145, 4 (2019) 2746–2760. [CrossRef] [PubMed] [Google Scholar]
  19. J. van Mourik, D. Murphy: Explicit higherorder FDTD schemes for 3D room acoustic simulation, IEEE Transactions on Audio, Speech, and Language Processing 22, 12 (2014) 2003–2011. [CrossRef] [Google Scholar]
  20. B. Hamilton, S. Bilbao: FDTD methods for 3D room acoustics simulation with high order accuracy in space and time, IEEE Transactions on Audio, Speech, and Language Processing 25, 11 (2017) 2112–2124. [CrossRef] [Google Scholar]
  21. S. Sakamoto: Phase-error analysis of high-order finite difference time domain scheme and its influence on calculation results of impulse response in closed sound field, Acoustical Science and Technology 28, 5 (2007) 295–309. [CrossRef] [Google Scholar]
  22. A. Emmanuelli, D. Dragna, S. Ollivier, P. Blanc-Benon: Characterization of topographic effects on sonic boom reflection by resolution of the Euler equations, Journal of the Acoustical Society of America 149, 4 (2021) 2437–2450. [CrossRef] [PubMed] [Google Scholar]
  23. L. Stein, F. Straube, J. Sesterhenn, S. Weinzierl, M. Lemke: Adjoint-based optimization of sound reinforcement including non-uniform flow, Journal of the Acoustical Society of America 146, 3 (2019) 1774–1785. [CrossRef] [PubMed] [Google Scholar]
  24. N.H. Fletcher: The nonlinear physics of musical instruments, Reports on Progress in Physics 62, 5 (1999) 723–764. [CrossRef] [Google Scholar]
  25. S. Bilbao, J. Chick: Finite difference time domain simulation for the brass instrument bore, Journal of the Acoustical Society of America 134, 5 (2013) 3860–3871. [CrossRef] [PubMed] [Google Scholar]
  26. H. Berjamin, B. Lombard, C. Vergez, E. Cottanceau: Time-domain numerical modeling of brass instruments including nonlinear wave propagation, viscothermal losses, and lips vibration, Acta Acustica united with Acustica 103, 1 (2017) 117–131. [CrossRef] [Google Scholar]
  27. A. Chaigne, A. Askenfelt: Numerical simulations of piano strings. I. A physical model for a struck string using finite difference methods, Journal of the Acoustical Society of America 95, 2 (1994) 1112–1118. [CrossRef] [Google Scholar]
  28. J. Bensa, S. Bilbao, R. Kronland-Martinet, J.O. Smith III: The simulation of piano string vibration: From physical models to finite difference schemes and digital waveguides, Journal of the Acoustical Society of America 114, 2 (2003) 1095–1107. [CrossRef] [PubMed] [Google Scholar]
  29. S.K. Lele: Compact finite difference schemes with spectral-like resolution, Journal of Computational Physics 103, 1 (1992) 16–42. [CrossRef] [Google Scholar]
  30. C.K.W. Tam, J.C. Webb: Dispersion relation-preserving finite difference schemes for computational acoustics, Journal of Computational Physics 107 (1993) 262–281. [CrossRef] [Google Scholar]
  31. E. Hairer, S.P. Nørsett, G. Wanner: Solving ordinary differential equations I: nonstiff problems. Springer series in computational mathematics, 2nd edn., Springer, Berlin, Heidelberg, 2000. [Google Scholar]
  32. M. Calvo, J.M. Franco, L. Randez: A new minimum storage Runge-Kutta scheme for computational acoustics, Journal of Sound and Vibration 201, 1 (2004) 1–12. [Google Scholar]
  33. F.Q. Hu, M.Y. Hussaini, J.L. Manthey: Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics, Journal of Computational Physics 124, 1 (1996) 177–191. [CrossRef] [Google Scholar]
  34. C.K.W. Tam, J.C. Webb, Z. Dong: A study of the short wave components in computational acoustics, Journal of Computational Acoustics 1 (1993) 1–30. [CrossRef] [Google Scholar]
  35. D.V. Gaitonde, M.R. Visbal: Pade-type higher-order boundary filters for the Navier-Stokes equations, AIAA Journal 38, 11 (2000) 2103–2112. [CrossRef] [Google Scholar]
  36. T.J. Poinsot, S.K. Lele: Boundary conditions for direct simulations of compressible viscous flows, Journal of Computational Physics 101, 1 (1992) 104–129. [CrossRef] [Google Scholar]
  37. F.Q. Hu: Development of PML absorbing boundary conditions for computational aeroacoustics: A progress review, Computers & Fluids 37 (2008) 336–348. [CrossRef] [Google Scholar]
  38. A. Mani: Analysis and optimization of numerical sponge layers as a nonreflective boundary treatment, Journal of Computational Physics 231 (2012) 704–716. [CrossRef] [Google Scholar]
  39. K. Attenborough: Acoustical characteristics of porous materials, Physics Reports 82, 3 (1982) 179–227. [CrossRef] [Google Scholar]
  40. M. Aretz, M. Vorländer: Efficient modelling of absorbing boundaries in room acoustic fe simulations, Acta Acustica united with Acustica 96, 6 (2010) 1042–1050. [CrossRef] [Google Scholar]
  41. C.S. Peskin: Flow patterns around heart valves: a numerical method, Journal of Computational Physics 10, 2 (1972) 252–271. [CrossRef] [Google Scholar]
  42. Y. Cai, S. van Ophem, W. Desmet, E. Deckers: Model order reduction of time-domain vibro-acoustic finite element simulations with nonlocally reacting absorbers, Computer Methods in Applied Mechanics and Engineering 416 (2023) 116345. [CrossRef] [Google Scholar]
  43. I. Moufid, D. Matignon, R. Roncen, E. Piot: Energy analysis and discretization of the time-domain equivalent fluid model for wave propagation in rigid porous media, Journal of Computational Physics 451 (2022) 110888. [CrossRef] [Google Scholar]
  44. A. Alomar, D. Dragna, M.A. Galland: Time-domain simulations of sound propagation in a flow duct with extended-reacting liners, Journal of Sound and Vibration 507 (2021) 116137. [CrossRef] [Google Scholar]
  45. J. Reiss: Pressure-tight and non-stiff volume penalization for compressible flows, Journal of Scientific Computing 90, 86 (2022) 1–29. [CrossRef] [Google Scholar]
  46. M. Lemke, J. Reiss: Approximate acoustic boundary conditions in the time-domain using volume penalization, Journal of the Acoustical Society of America 153, 2 (2023) 1219–1228. [CrossRef] [PubMed] [Google Scholar]
  47. S. Bilbao: Immersed boundary methods in wave-based virtual acoustics, Journal of the Acoustical Society of America 151, 3 (2022) 1627–1638. [CrossRef] [PubMed] [Google Scholar]
  48. S. Bilbao: Modeling impedance boundary conditions and acoustic barriers using the immersed boundary method: The one-dimensional case, Journal of the Acoustical Society of America 153, 4 (2023) 2023–2036. [CrossRef] [PubMed] [Google Scholar]
  49. W.X. Huang, F.B. Tian: Recent trends and progress in the immersed boundary method, Proceedings of the Institution of Mechanical Engineers. Part C: Journal of Mechanical Engineering Science 233, 23–24 (2019) 7617–7636. [CrossRef] [Google Scholar]
  50. L.D. Landau, E.M. Lifshitz: Fluid mechanics. Course of theoretical physics, vol. 6, 2nd edn., Butterworth-Heinemann, London, , 1987. [Google Scholar]
  51. E.F. Toro: Riemann solvers and numerical methods for fluid dynamics: a practical introduction, 3rd edn., Springer Science & Business Media, Berlin, Heidelberg, 2009. [CrossRef] [Google Scholar]
  52. K. Van den Abeele, J. Ramboer, G. Ghorbaniasl, C. Lacor: Numerical solution of the linearized euler equations using compact schemes, in: H. Deconinck, E. Dick (eds), Computational fluid dynamics 2006, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, pp. 837–842. [CrossRef] [Google Scholar]
  53. C. Hirsch: Numerical computation of internal and external flows, volume 1: fundamentals of numerical discretization, John Wiley and Sons, New York, 1988 [Google Scholar]
  54. L.N. Trefethen: Group velocity in finite difference schemes, SIAM Review 24, 2 (1982) 113–136. [CrossRef] [Google Scholar]
  55. J.R. Bunch, J.E. Hopcroft: Triangular factorization and inversion by fast matrix multiplication, Mathematics of Computation 28, 125 (1974) 231–236. [CrossRef] [Google Scholar]
  56. S.D. Conte, C. de Boor: Elementary numerical analysis, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017. [CrossRef] [Google Scholar]
  57. M.B. Giles, N.A. Pierce: An introduction to the adjoint approach to design, Flow, Turbulence and Combustion 65, 3 (2000) 393–415. [CrossRef] [Google Scholar]
  58. K.L. Shlager, J.G. Maloney, S.L. Ray, A.F. Peterson: Relative accuracy of several finite difference time-domain methods in two and three dimensions, IEEE Transactions on Antennas and Propagation 41, 12 (1993) 1732–1737. [CrossRef] [Google Scholar]
  59. B. Finkelstein, R. Kastner: Finite difference time domain dispersion reduction schemes, Journal of Computational Physics 221, 1 (2007) 422–438. [CrossRef] [Google Scholar]
  60. A. Iserles: A first course in the numerical analysis of differential equations (2nd edn.). Cambridge Texts in Applied Mathematics, Cambridge University Press, 2008. [Google Scholar]
  61. J.C. Butcher: Numerical methods for ordinary differential equations in the 20th century, Journal of Computational and Applied Mathematics 125, 1 (2000) 1–29. [CrossRef] [Google Scholar]
  62. A. Hindmarsh, P. Gresho, D. Griffiths: The stability of explicit Euler time-integration for certain finite difference approximations of the multidimensional advection-diffusion equation, International Journal for Numerical Methods in Fluids 9, 4 (1984) 853–897. [CrossRef] [Google Scholar]
  63. F. Nazari, A. Mohammadian, M. Charron: High-order low-dissipation low-dispersion diagonally implicit Runge-Kutta schemes, Journal of Computational Physics 286 (2015) 38–48. [CrossRef] [Google Scholar]
  64. J. Berland, C. Bogey, C. Bailly: Low-dissipation and low-dispersion fourth-order Runge–Kutta algorithm, Computers & Fluids 35 (2006) 1459–1463. [CrossRef] [Google Scholar]
  65. K.W. Thompson: Time dependent boundary conditions for hyperbolic systems, Journal of Computational Physics 68, 1 (1987) 1–24. [CrossRef] [Google Scholar]
  66. F. Pled, C. Desceliers: Review and recent developments on the perfectly matched layer (PML) method for the numerical modeling and simulation of elastic wave propagation in unbounded domains, Archives of Computational Methods in Engineering 29, 1 (2022) 471–518. [CrossRef] [Google Scholar]
  67. R. Troian, D. Dragna, C. Bailly, M.A. Galland: Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow, Journal of Sound and Vibration 392 (2016) 200–216. [Google Scholar]
  68. X. Li, X. Li, C. Tam: Construction and validation of a broadband time domain impedance boundary condition, in: 17th AIAA/CEAS Aeroacoustics Conference (32nd AIAA Aeroacoustics Conference), Portland, Oregon, 5–8 June, 2011. [Google Scholar]
  69. H. Deconinek, C. Hirsch, J. Peuteman: Characteristic decomposition methods for the multidimensional euler equations, in: F.G. Zhuang, Y.L. Zhu (eds), Tenth International Conference on Numerical Methods in Fluid Dynamics, Springer, Berlin, Heidelberg, 1986, pp. 216–221. [CrossRef] [Google Scholar]
  70. E. Albin, Y. D’Angelo, L. Vervisch: Flow streamline based Navier–Stokes characteristic boundary conditions: Modeling for transverse and corner outflows, Computers & Fluids 51, 1 (2011) 115–126. [CrossRef] [Google Scholar]
  71. R.D. Sandberg, N.D. Sandham: Nonreflecting zonal characteristic boundary condition for direct numerical simulation of aerodynamic sound, AIAA Journal 44, 2 (2006) 402–405. [CrossRef] [Google Scholar]
  72. C. Bogey, C. Bailly, D. Juvé: Numerical simulation of sound generated by vortex pairing in a mixing layer, AIAA Journal 38, 12 (2000) 2210–2218. [CrossRef] [Google Scholar]
  73. P. de Palma, M. de Tullio, G. Pascazio, M. Napolitano: An immersed-boundary method for compressible viscous flows, Computers & Fluids 35, 7 (2006) 693–702. [CrossRef] [Google Scholar]
  74. Q. Liu, O.V. Vasilyev: A brinkman penalization method for compressible flows in complex geometries, Journal of Computational Physics 227, 2 (2007) 946–966. [CrossRef] [Google Scholar]
  75. N.K.R. Kevlahan, T. Dubos, M. Aechtner: Adaptive wavelet simulation of global ocean dynamics using a new brinkman volume penalization, Geoscientific Model Development 8, 12 (2015) 3891–3909. [CrossRef] [Google Scholar]
  76. F. Kemm, E. Gaburro, F. Thein, M. Dumbser: A simple diffuse interface approach for compressible flows around moving solids of arbitrary shape based on a reduced Baer–Nunziato model, Computers & Fluids 204 (2020) 104536. [CrossRef] [Google Scholar]
  77. J.T. Katsikadelis: Boundary elements. Theory and applications, Elsevier, Oxford, 2002. [Google Scholar]
  78. M. Lemke, L. Stein, A. Hölter, Y. Schubert: cruna-toolkit, 2024. Available at https://github.com/cruna-toolkit. [Google Scholar]
  79. D.A. Webster, D.T. Blackstock: Finiteamplitude saturation of plane sound waves in air, Journal of the Acoustical Society of America 62, 3 (1977) 518–523. [CrossRef] [Google Scholar]
  80. F. Brinkmann, L. Aspöck, D. Ackermann, R. Opdam, M. Vorländer, S. Weinzierl: A benchmark for room acoustical simulation. Concept and database, Applied Acoustics 176 (2021) 107867. [Google Scholar]
  81. F. Straube, F. Schultz, D. Albanes Bonillo, S. Weinzierl: An analytical approach for optimizing the curving of line source arrays, Journal of the Audio Engineering Society 66, 1/2 (2018) 4–20. [CrossRef] [Google Scholar]
  82. d&b audio: 2024. https://www.dbaudio.com/global/en/products/all/series/sl-series/gsl8/#tab-technicaldata, accessed: 2024–09-20 [Google Scholar]
  83. S. Weinzierl, S. Lepa, F. Schultz, E. Detzner, H. von Coler, G. Behler: Sound power and timbre as cues for the dynamic strength of orchestral instruments, Journal of the Acoustical Society of America 144, 3 (2018) 1347–1355. [CrossRef] [PubMed] [Google Scholar]

Cite this article as: Hölter A. Weinzierl S. & Lemke M. 2024. Finite difference time domain discretization for room acoustic simulation based on the non-linear Euler equations. Acta Acustica, 8, 75. https://doi.org/10.1051/aacus/2024071.

All Tables

Table 1

Positions of the sound sources.

In the text

All Figures

thumbnail Figure 1

(a) Harmonic analysis of the dispersion relation c′/c and (b) the error |k′Δx − kΔx|real of different finite-difference-stencils for first order partial differential equations. (c) CPU-time with respect to the error. For the analysis, the vector length was chosen to match PPW resolution and dispersion error from (b). The short forms OX denote the order, exp or imp, explicit or implicit / compact derivatives (Eq. (8)). The error of 10−3 is highlighted by a black dotted line as it is used for discussion the text.
In the text
thumbnail Figure 2

Analysis of selected explicit time integration schemes. The short forms OX denote the order, RK the Runge-Kutta schemes and AB the Adams-Bashfourth schemes. (a) Stability regions for selected explicit time integration schemes. (b) Maximum stable time increment for different methods i normalized with respect to the largest value in blue. Maximum time increment additionally normalized with respect to the number of rhs-calls (evaluations of right-hand-sides of the governing equations) needed for time integration indicating the numerical efficiency in orange.
In the text
thumbnail Figure 3

Frequency response of different explicit (exp) and implicit (imp) filters. Cut-off frequencies are at 0.5$ \sqrt{0.5}$ marked by the black dashed line.
In the text
thumbnail Figure 4

Frequency response of 2nd and 4th order explicit and 4th-order implicit spatial derivatives. The receiver is located at xr = [0.7, 0.8, 0.7]T m and the source is located at xs = [0.2, 0.4, 0.3]T m. The analytical solution is shown in orange dashed and the numerical solution in blue. The gray area marks regions with a theoretical error (k’Δx − kΔx)real > 10−3. The lower non-linear case shows results for high amplitudes, at 0.75 m distance from the source the peak magnitude of the blue line is at ≈146 dBSPL and in light blue dashed at ≈152 dBSPL.
In the text
thumbnail Figure 5

Errors in simulated modal frequencies for the cube-shaped room with rigid boundaries.
In the text
thumbnail Figure 6

Frequency response of 4th-order explicit, 4th- and 8th-order implicit filter schemes. The receiver is at xr = [0.7, 0.8, 0.7]T m and the source is at xs = [0.2, 0.4, 0.3]T m. The analytical solution is shown in orange dashed and the numerical solution in blue.
In the text
thumbnail Figure 7

Results of the FDTD simulation (blue) compared to experimental results of the BRAS database, scene RS1 (red dashed line) and simulation results of the BRAS E solver (yellow). In case (a), the blue dashed vertical lines show the analytical comb filter frequencies and the gray dashed vertical lines show the number of PPWs. (a) Rigid, (b) Absorber, (c) Diffuser.
In the text
thumbnail Figure 8

Sketch of the shoebox room. Edges of the room are the black lines, gray areas are porous absorber. Red dots are the source positions, see Table 1 and the blue dot marks the receiver position at x = [8.0; 5.0; 1.5] m.
In the text
thumbnail Figure 9

Differences of the frequency response of linear and nonlinear equations at the measurement position x = [8.0; 5.0; 1.5] m.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.