Time-domain modelling of wave-based room acoustics including viscothermal and relaxation effects in air

: Air absorption can be a signiﬁcant source of attenuation, which should be considered in long-duration wideband acoustics simulations. In this short contribution, a time-domain model for three-dimensional wave propagation including vis-cothermal and relaxation effects (air absorption) is developed and coupled with locally reactive impedance wall conditions through a conservative energy framework. The model is discretised with a ﬁnite-difference time-domain method, and numerical stability is established with a discrete energy balance. Numerical examples are presented to demonstrate the proposed method


Introduction
Wave-based methods for the simulation of room acoustics have seen increased interest recently (Das et al., 2021;Wilson, 2020;Oxnard, 2018;Pind et al., 2019;Saarelma et al., 2018;Savioja and Xiang, 2020;Stein et al., 2020); such methods have great potential for acoustic prediction and auralization across the entire audible range of frequencies.For fullbandwidth wave-based simulations, it is important to consider high-frequency air attenuation (Bass et al., 1995;ISO 9613-1, 1993) comprised of viscothermal and relaxation effects (Pierce, 2019).In the case of time-domain simulations [e.g., finite-difference time-domain (FDTD)] (Botteldooren, 1994;Kowalczyk and van Walstijn, 2011), ideally this should be accomplished while maintaining passivity of the discrete system (including boundary conditions) to ensure numerical stability.In related previous work, frequency-dependent damping from air attenuation effects has been efficiently simulated in the linear, free-space case through modal domain-decomposition schemes (Botts and Savioja, 2015).Also, non-linear sound propagation including viscothermal relaxation effects has been simulated with FDTD methods (Jim enez et al., 2016;Wochner et al., 2005).
In the context of room acoustics-including non-rigid boundary conditions in non-trivial geometries-a full model of sound propagation with viscothermal and relaxation effects has yet to be presented in a concretely passive continuous formulation with numerical stability ensured in the discrete case.This has been accomplished in the absence of relaxation effects, modeling Stokes' equation with FDTD methods (Webb and Bilbao, 2011) and finitevolume time-domain (FVTD) methods (Bilbao and Hamilton, 2017).Under certain conditions, air attenuation from Stokes' equation may be sufficient for audible frequencies (Hamilton et al., 2014;Pierce, 2019), but it remains of general interest to extend such approaches to include relaxation effects.While filter-based approaches for adding air absorption to room impulse responses are available (Kates and Brandewie, 2020;Saarelma and Savioja, 2016), the interest here is the more general and direct simulation of air absorption processes, which allows for moving sources and receivers and the embedding of distributed sound sources (Bilbao and Ahrens, 2020) in both near-and far-field settings.
The aim of this short contribution is to provide a suitable model for sound propagation in the context of time-domain room acoustics simulation with the inclusion of viscothermal and relaxation effects in air, along with nonrigid wall conditions while maintaining passivity as a whole.Using judicious choices of FDTD operators in the context of the simplest Cartesian scheme in three dimensions, this system may be discretised in a numerically stable manner over non-trivial geometries by ensuring that the property of passivity is transferred to discrete time.This work can be seen as a more complete presentation of a contribution treating just free-space propagation (Hamilton and Bilbao, 2020) and may also be seen as an extension of previous work by the same authors (Bilbao and Hamilton, 2017;  Hamilton et al., 2016).a) Author to whom correspondence should be addressed.b) ORCID: 0000-0001-5332-9887.

Linear acoustic equations for air
The following system describes linear sound propagation in air including viscothermal and relaxation effects under the assumption of zero mean flow and an irrotational velocity field (Pierce, 2019): Here, the dependent variables q, p, and v are the density, acoustic pressure, and particle velocity, respectively, and the thermodynamic variables s fr and T are the entropy and deviation from ambient temperature, respectively.All are functions of t 2 R and spatial coordinates x 2 X & R 3 .To accommodate the viscothermal effects of a mixture of gases, the additional variables T ¼ T ðt; xÞ represent the apparent vibration temperature associated with a -type molecule (oxygen or nitrogen).@ t represents partial differentiation with respect to time t, and r, rÁ, and r 2 are the three-dimensional gradient, divergence, and vector Laplacian operators, respectively.Under the approximation q % ð1=c 2 Þp, (1a) and (1b) reduce to a second-order viscothermal wave equation (Stokes' equation with viscous and thermal effects) and further to the lossless wave equation for l ¼ 0. The constants that appear above are q 0 , the ambient density of air; c, the frozen sound speed in the highfrequency limit; l ¼ l B þ 4 3 l, a combination of the bulk viscosity l B and shear viscosity l; T 0 , the ambient temperature; b, the coefficient of thermal expansion; C p , the specific heat at constant pressure; c, the ratio of specific heats; j, the coefficient of thermal conductivity; and C v and s , the specific heat associated with internal vibrations of a -type molecule and its relaxation time.For a gas with ratio of specific heats c, we also have b ¼ 1=T 0 and ðc À 1Þ=c 2 ¼ b=C p .
Without further approximations, (1a)-(1e) can be reduced to where ' h ¼ j=q 0 cC p and ' v ¼ ð1=q 0 cÞ l are length constants associated with thermal and viscous effects, respectively.At this point, we can make use of the approximation T % ðT 0 b=q 0 C p Þp, as employed in the linear case for the acoustic mode (Pierce, 2007(Pierce, , 2019)).Then (2a) and (2b) accompanied by (1f) simplify to where ã ¼ ðc À 1ÞðC v =C p Þ, and p ¼ q 0 C p T is a pressure associated with a -type relaxation process.The combination of (3a) and (3b) results in the equivalent a second-order system [with (3c)] (Hamilton and Bilbao, 2020), where g ¼ ' v þ ðc À 1Þ' h .Thus, we have derived a second-order wave equation including viscothermal and relaxation effects [with (3c)]. 1 It is worth noting that a similar third-order-in-time system, which appears in Pierce (2007), can be achieved through the substitution cg@ t r 2 p % ðg=cÞ@ 3 t p. Required constants for (4) include c, g, and ã for oxygen and nitrogen.These are temperature, humidity, and atmospheric-pressure dependent, for which formulas may be found in Bass et al. (1995), ISO 9613-1 (1993), and Pierce (2019).

Dispersion relation
It is worthwhile to confirm that this system gives the desired dispersion relation (Pierce, 2019).Consider a plane wave solution of the form pðt; xÞ ¼ pe jð knÁxÀxtÞ , with direction n, angular frequency x, complex wavenumber k, and complex amplitude p (and j is the imaginary unit).Similarly, let p ðt; xÞ ¼ p e jð knÁxÀxtÞ with complex amplitude p .Inserting these solutions into the second-order system and simplifying, we arrive at 2 To solve for k ¼ k þ ja with k and a real, we can assume x ( c=g and take an approximate square-root to arrive at k % x=c and the desired air attenuation coefficient (ISO 9613-1, 1993;Pierce, 2019), (6)

Energy analysis
Boundary conditions suitable for room acoustics can be derived through an energy analysis.For this purpose, we introduce q as an intermediate variable, defined by @ t q ¼ p. Then we multiply both sides of (3a) by p and use (3b) and (3c) to get Integrating over the volume X with boundary surface C of outward normal n and using Green's first identity results in The above is an energy balance of the form ðd=dtÞEðtÞ ¼ ÀQðtÞ þ BðtÞ, with internal energy EðtÞ !0, power loss QðtÞ !0, and a boundary term BðtÞ of indeterminate sign.For the energy of the system to be non-increasing (passive), it is sufficient that BðtÞ 0. One possible boundary condition is (for where Y is frequency-independent specific wall admittance with Y !0, and ðrqÞ ?¼ n Á rq, ðrpÞ ?¼ n Á rp, and v ?¼ n Á v.As such, we have a simple locally reactive immittance boundary condition with a viscothermal perturbation, which is necessary to preserve passivity.For x ( c=g, this has a normal-incidence reflection coefficient and it is guaranteed that jRj 1 by passivity of the system (with Y !0).In terms of an acoustic velocity potential W ¼ Àrv and for p % q 0 @ t W, (9) is equivalent to a boundary condition derived in Bilbao and Hamilton (2017) in the frequency-independent case for the simpler Stokes' equation.As such, an extension to (locally reactive) frequency-dependent boundary conditions [as in Bilbao and Hamilton (2017)] follows for this system as well, but such further developments are omitted here for brevity.

Definitions
In this section, we build a simple Cartesian FDTD scheme for the systems of interest.Consider the integer lattice with points at i 2 Z 3 and standard unit vectors êw for w 2 fx; y; zg.Henceforth, w is used as a direction placeholder, and it is implied that w 2 fx; y; zg (unless otherwise specified).We then consider a spatiotemporal grid function p n i u pðnD t ; iD x Þ defined at locations iD x and times nD t , with grid spacing D x , time step D t , and integer time-index n !0. Similarly, let us define p n ;i u p ðnD t ; iD x Þ.For the particle velocity vðt; xÞ ¼ ½v x ðt; xÞ; v y ðt; xÞ; v z ðt; xÞ T , it is helpful to think of the points iD x as the centres of cubic cells with side-length D x .Using the faces of the cubic cells we define a vector grid function v n i comprised of three scalar grid functions: For compactness of notation, we use superscript n in v n i , but it should be remembered that it is staggered in time from p n i (and similarly in space).The use of a staggered grid is traditional (Botteldooren, 1994) to ensure centered differences through the scheme, which will be the case unless otherwise specified.
is taken to the inner product of a vector grid function v n i , such that jv n i j is a scalar grid function.
For an arbitrary grid function u n i , we can define temporal shift, difference, and averaging operators, respectively, Similarly, along a spatial coordinate w, we can define shift and difference operators.These are, respectively, Furthermore, we can define a second-order temporal difference, a discrete gradient, and a discrete Laplacian operator, respectively, 3.2 Staggered FDTD scheme for the first-order system A staggered FDTD scheme for the first-order system (3) follows from substitution of continuous variables with corresponding discrete grid functions and partial differential operators with difference operators.There are many such choices, but one possible FDTD scheme allowing explicit operation and easily determined numerical stability conditions is the following: All differences are centered due to staggering of the grid, except for the term d rþ d tÀ p n i in (13b).The uncentered difference is required to keep the scheme explicit (avoiding the need for a linear system solution at each time step).In theory, this reduces the order of accuracy of the scheme to first-order (yet the scheme remains consistent with the target differential system).

FDTD scheme for scalar second-order system
For the second-order system, it is convenient to make use of the variable qðt; xÞ and define q n i u qððn þ 1=2ÞD t ; iD x Þ, which can be related to p n i by d tÀ q n i ¼ p n i .In this case, Eqs.(13a) and (13c) become where d tÁ ¼ l tþ d tÀ .Equation (14a) is a scheme for the second-order system that may be simulated without the use of the vector velocity grid function v n i (thereby saving on storage). 3Meanwhile, (13b) is available for an auxiliary calculation of v n i if necessary. 4If it is preferred to simulate pressure directly, applying d tÀ to (14a) gives the scheme (Hamilton and Bilbao, 2020) which also requires (13c).In the above, and in (14a), the viscothermal term is uncentered in time to keep the scheme explicit [as discussed previously in regard to the term d rþ d tÀ p n i in (13b)].

Energy analysis
We can carry out a discrete energy analysis to derive suitable boundary conditions and establish passivity of the system (numerical stability).For this, let G represent the set of points i for which iD x 2 X (interior points).We then take a spatial inner product over (14a) against the grid function d tÁ q n i and use (14b) to get a discrete analogue to (7), Then in analogy with (8), we have a discrete energy balance of the form Here, d 0 r6 q n i is d r6 q n i with directional differences d w6 q n i evaluated only when i 6 êw 2 G (when both points in difference are interior), otherwise returning zero scalar components (Hamilton, 2016).Its complement is d 0 r6 q n i ¼ d r6 q n i À d 0 r6 q n i , which, in turn, returns directional differences across the boundary, including "ghost points" (i.e., points i 6 2 G with at least one i6ê w 2 G).
For passivity, it is required that E n !0, Q n !0, and B n 0. Clearly, Q n is non-negative, but the signs of E n and B n are indeterminate.To examine this, we can use the following identity and bound (Bilbao, 2009;Hamilton, 2016): Thus, for E n !0, we can group together terms with ðd tÀ q n i Þ 2 to arrive at the condition À g, which can be seen as the global necessary stability condition for the scheme in the absence of boundary conditions.For passivity on the whole, we need to ensure that B n 0. This is ensured by the following numerical boundary conditions, mimicking (9): This boundary condition does not interfere with the global stability condition above, which itself, importantly, is not dependent on relaxation effects.This independence is ultimately a by-product of the difference operators chosen for the scheme.In structure, this FDTD scheme can be seen as a generalisation of that presented in (Bilbao and Hamilton, 2017;Webb and Bilbao, 2011) (replacing q by a scalar acoustic velocity potential W).It follows that it is possible to extend this FDTD scheme to non-Cartesian grids (Hamilton, 2016) and, further, to fully unstructured grids including frequencydependent impedance boundary conditions that are spatially varying (Bilbao and Hamilton, 2017).

Update equation
The FDTD scheme can be solved recursively in time using its corresponding update equation.To provide this, it will help to define the neighbour set N ðiÞ as the interior nearest-neighbours of i; i.e., N ðiÞ ¼ fi þ êw 2 Gg, and its cardinality is denoted jN ðiÞj.Expanding (14a) and (14b) and using the boundary condition (19) to eliminate "ghost points" yields the (two-stage) update equations, , and b ¼ a =1 þ a ; also where Setting initial conditions q 0 i ; q 1 i ; p 0 ;i ; p 1 ;i (e.g., with localised distributions) may be used to excite the scheme [details left out for brevity; see, e.g., Strikwerda (2004)].

Numerical air attenuation
To check consistency of the scheme, we consider the numerical dispersion relation through insertion of a plane wave ansatz q n i ¼ qe jð k ÁiDxÀxnDt Þ [with p n ;i ¼ p e jð k ÁiDxÀxnDt Þ ] propagating along a Cartesian-aligned axis.Inserting the ansatz into (14) and solving for the imaginary part of the wavenumber, =f kg gives the numerical air attenuation from the scheme, which can be compared to the desired attenuation (6).This comparison is shown in Fig. 1 for three air conditions [with other required constants derived from ISO 9613-1 (1993)] with the grid spacing D x chosen for 7.5 points per wavelength (PPW) at f max ¼ 80 kHz [chosen to demonstrate the full range of air attenuation effects while keeping less than 2% phase velocity errors (Hamilton, 2016)], which corresponds to D x ¼ 0:57 mm.Each specific air condition has associated constants c and g, according to ISO 9613-1 (1993), and for each air condition, D t is set to the associated stability limit of the scheme.It can be seen that the scheme reproduces the desired air attenuation in general except for the appearance of erroneous numerical attenuation near the limiting frequencies simulated by the scheme, which is expected by approximation errors in the FDTD scheme.Such deviations increase with frequency, reaching as high as 6% relative (dB) error at f max in these cases.Simulation measurements that also display consistency may be found in Hamilton and Bilbao (2020).

Complex geometry
To demonstrate the energy conservation properties of the FDTD scheme [detailed in ( 17)], we consider the non-trivial geometry shown in Fig. 2(a) with specific admittance Y ¼ 0:01 on walls.Simulations are carried out with air temperature 10 C and relative humidity 15%, using raised cosine (Hann window) initial conditions of varying spatial supports L (widths) with L 2 f0:5; 0:25; 0:125g m, centered about the source position (as pictured).The grid spacing D x is set to 1 cm (7.5 PPW at 16 kHz), and D t is set to the stability limit (D t % 18:3 ls).A two-dimensional (2D) slice of the pressure field after 4 ms is shown in Fig. 2(a).Discrete energy quantities are also calculated for the system at each time step and plotted in Fig. 2(c).Supplementary video files show the pressure field and discrete energies evolving over time for initial condition L ¼ 0.5 m (Mm.1), L ¼ 0.25 m (Mm.2), and L ¼ 0.125 m (Mm.3).It may be observed that the total energy of the simulation (internal energy plus accumulated remains constant to within 12 decimal places (variations due to machine precision), which is a result of the discrete energy balance of the scheme.

Conclusions
In this study, linear sound propagation in air including viscothermal and relaxations is coupled with locally reactive immittance wall conditions in a passive framework.A FDTD discretisation is provided with discrete energy conservation and conditions for passivity and thus numerical stability.Numerical attenuation was compared to target air attenuation curves, and simulation examples were provided to demonstrate the application of the scheme to a non-trivial room acoustic scenario, with discrete energy conserved to machine precision.The work presented here uses the simplest Cartesian scheme in three dimensions, which is known to feature staircasing errors.A generalisation to a finite-volume framework (Bilbao and Hamilton, 2017), which fixes such errors, is possible and constitutes future work, along with complex impedance and boundary conditions that may vary (spatially) over the boundary.Optimised schemes, possibly with higher than first-order accuracy, could be investigated as well as the inclusion of realistic source models (Bilbao and Ahrens, 2020) and performance on parallel computing architectures (Hamilton et al., 2016;Saarelma et al., 2018).

Fig. 2 .
Fig. 2. (a) Room setup with source location (denoted by "S") and axis units in meters.(b) Snapshot of wave propagation at 4 ms through the y-z plane (at source x-location) for the case L ¼ 0.25 m.(c) Plot of discrete energy quantities as a function of time.