COMSOL 6.4 - Wave-Based Time-Domain Room Acoustics with Frequency-Dependent Impedance

Wave-Based Time-Domain Room Acoustics with Frequency-Dependent Impedance

The use of wave-based techniques for room acoustic simulations has spread in the last years due to the increase in computational performance as well as the development of new numerical methods. The challenge of including realistic impedance conditions at walls is traditionally solved in the frequency domain. Recent research has focused on the implementation of frequency-dependent impedance in the time-domain models using partial fraction representation of the frequency-dependent data (see Ref. 1, Ref. 2, Ref. 3).

This tutorial shows how to get partial fraction representation of frequency-dependent wall impedance data via the Partial Fraction Fit function and how to use the results to set up the built-in impedance boundary conditions for room acoustic response simulation in the time domain. The model uses the Pressure Acoustics, Time Explicit interface to simulate the propagation of sound. The physics interface is based on the discontinuous Galerkin (dG-FEM) method which uses a matrix free approach and a time explicit solver. The method is very memory efficient and well suited for distributed computing on a cluster architecture.

Model Definition

The model studies the response of a small 50.5 m3 room shown in Figure 1 on the left. The loudspeaker to the left of the TV emits an acoustic signal given as a sine wave modulated by a Gaussian envelope at the center frequency f0 = 700 Hz. As the signal propagates from the speaker, the acoustic pressure is measured at four listening points located equidistantly on a line drawn from the speaker to the sofa, as shown in Figure 1 on the right.

Figure 1: Room geometry (left) and location of listening points (right).

The sound absorption properties of the carpet, ceiling, sofa, and walls are modeled using the simplified local reacting approximation which states that the response at a certain point of the surface depends on the sound pressure at that point. This behavior is described by an impedance boundary condition imposed on the surface. When not purely resistive, as, for example, for porous materials, the boundary impedance will depend on the frequency, resulting in the following boundary condition in the frequency domain

(1)

where Z is the frequency-dependent specific impedance and Y is the admittance. The time-domain equivalent of Equation 1 involves a convolution instead of the multiplication

(2)

and thus requires the knowledge of Y(t) in the time domain and thus the computation of the inverse Fourier (or Laplace) transform, which, obviously, cannot be done analytically for a general Y(ω) obtained from various poroacoustic models or measurements. Moreover, it is not enough to know the values of the impedance (admittance) for a number of frequencies or for a frequency range to translate the boundary condition from the frequency to the time domain. Therefore, Y(ω) has to be extended to the whole complex plane.

An approximate analytical representation of Y(t) can be retrieved from a rational approximation of Y(ω) that is defined on the whole complex plane:

(3)

Indeed, each fraction term on the right-hand side of Equation 3 corresponds to an exponential decay in the time domain

(4)

where L−1 is the inverse Laplace transform and H(t) is the Heaviside step function.

However, the approximation given by Equation 3 has to fulfill three conditions in order for the time-domain boundary condition to be physical:

•

Causality, Y(ω) is analytic and nonzero in ℑ(ω) > 0;

•

Reality,

, that is, ℜ(Y(ω)) is even and ℑ(Y(ω)) is odd;

•

Passivity, ℜ(Y(ω)) > 0 for all real ω.

The causality and reality conditions are fulfilled if Y∞ is real; the residues, Ak, and poles, αk, are either real or come in complex-conjugate pairs; and ℜ(αk) < 0 (for the exponentials in Equation 4 to decay as the time increases). That is

(5)

where NR and NC are the numbers of pure real poles and complex-conjugate pole pairs, respectively.

The form given in Equation 5 is used to reduce the evaluation of the integral in Equation 2 to system of auxiliary ordinary differential equations (ODEs) for memory variables. This approach is referred to as ADE method (see Ref. 1, Ref. 2, Ref. 3). The system of ODEs is automatically created and solved when you set up an Impedance boundary condition with a option.

Results and Discussion

The real and imaginary parts of the fitted frequency dependent admittance data are depicted in Figure 2. As seen, the sofa surface is more absorptive that the others, especially at the higher frequencies.

Figure 3 shows the acoustic pressure recorded at the listening points in blue, green, red, and cyan colors as the point moves away from the source. The values are normalized to the maximum pressure at the source. As seen from the plot, the pressure magnitude decreases with the time.

The history of the signal propagation through the room is shown in Figure 4. The signal emitted from the loudspeaker reaches listening point 1 at t = 4T0 (T0 = 1/f0) and listening point 4 at t = 11T0. Then multiple reflections occur while the signal magnitude goes down at t = 18T0 and t = 25T0.

Figure 2: Original admittance data and partial fraction expansions for the carpet, ceiling, sofa, and wall (from the top, left to right).

Figure 3: Acoustic pressure at listening points.

Figure 4: Normalized acoustic pressure at t = 4T0, 11T0, 18T0, and 25T0 (from the top, left to right).

Notes About the COMSOL Implementation

Partial Fraction Fit of Admittance Data

The partial fraction approximation (expansion) is obtained for the carpet, ceiling, sofa, and walls admittance through the built-in Partial Fraction Fit function within the frequency range from 50 Hz to 1.5 kHz. The form given by Equation 5 ensures that the reality condition is fulfilled. This is not always the case for the causality and passivity conditions. However, the Partial Fraction Fit function has the necessary tools that can be used to make the result fulfill the causality and reality conditions, thus making the time-domain impedance boundary condition physical.

The Partial Fraction Fit function fits the input data within a given tolerance that is found in the section (default 10−3). A higher tolerance error results in a higher degree, N, of the polynomials in Equation 3 and therefore a larger number of terms in the expansion Equation 5. The number of poles/residues in the expansion in turns corresponds to the number of ODEs solved for the memory variables. The default tolerance provides a balance between the approximation accuracy and the number of terms in the expansion, thus the computation costs when the expansion is used in an impedance boundary condition. On the other hand, increasing the tolerance may yield better results.

For example, the first Partial Fraction Fit function present in the model fits the admittance of the carpet. The default tolerance results in one real-valued pole and one complex-valued pole pair. From a distance the result looks good, but a closer look reveals that the passivity condition is violated at low frequencies: the real part of Y becomes negative between 50 and 90 Hz. A tighter tolerance of 10−5 yields two extra terms with real-valued poles and a partial fraction expansion that is passive within the given frequency range (see Figure 5).

Figure 5: Real part of the carpet admittance fit obtained with tolerance of 10-3 and 10-5.

The violation of the causality condition takes place when one or more poles in the partial fraction expansion have positive real part (that is, they are unstable), which can usually happen to the pure real poles. If an unstable pole has a relatively large residue, click the button to flip it onto the left half plane and circumvent the issue. The residues will be recomputed and updated automatically. Otherwise, enable the option and run the fit. The poles with relatively small residues will be removed during the cleanup process. This procedure is carried out when fitting the ceiling admittance data with the second Partial Fraction Fit function.

A different situation appears with the third Partial Fraction Fit function that is used to build a partial fraction expansion for the sofa admittance. Fitting the data with the default tolerance yields an approximation that fulfills all three conditions. However, the absolute value of the real-valued pole is much larger than that of the complex-valued pole. This results in a stiff ODE, which affects the stability of the time-integration scheme. Indeed, the Pressure Acoustics, Time Explicit physics interface used in this tutorial relies on explicit time-integration schemes. The time step is deduced solely from the minimum mesh element size (see the section below) and the speed of sound, and it can be too large to achieve a stable solution of the ODE.

The influence of such poles is often localized and therefore they can be discarded with the following update of the residues by clicking the button. For the default tolerance, the approximation becomes worse after discarding the real-valued pole (see Figure 6 on the left). Increasing the tolerance results in an extra complex-valued pole pair, and the result visually remains the same after the real-valued pole has been removed.

Figure 6: Sofa admittance approximation computed with the tolerance of 10-3 (left) and
10-5 (right) after discarding the real-valued pole.

The same procedure should be applied when the expansion contains spurious poles (also known as Froissart doublets). Those are the poles with residues close to zero, which is equivalent to very near pole and zero pairs in the rational approximation given by Equation 3.

Mesh and Time Explicit Solver

Solving wave propagation problems in the time domain has some requirements on both spatial and temporal resolution of the wave pattern. The mesh has to be fine enough to resolve the frequency content of the signal. For the quartic discretization used by default in the Nonlinear Pressure Acoustic, Time Explicit physics interface, the proper accuracy is achieved when the maximum mesh element size does not exceed 2/3 of the minimum wavelength. This tutorial studies the propagation of a broadband signal shown together with its frequency content in Figure 7. As seen, the frequency content is not limited by f0 = 700 Hz. Therefore, the mesh should resolve smaller wavelengths to achieve accurate results; at least up to λ0/2 = c0/(2f0). This yields hmax ≤ λ0/3.

Figure 7: Source signal (left) and its frequency content (right).

The Pressure Acoustic, Time Explicit interface is based on dG-FEM and uses an explicit time integration schemes. The time step is supposed to obey the CFL condition to ensure the stability of the time integration method. That is, Δt ≤ hmin/c0, where hmin is the minimum mesh element size. The solver automatically deduces the time step from the mesh and the speed of sound.

References

1. H. Wang, M. Cosnefroy, and M. Hornikx, “An arbitrary high-order discontinuous Galerkin method with local time-stepping for linear acoustic wave propagation,” J. Acoust. Soc. Am., vol. 149, p. 569, 2021; doi.org/10.1121/10.0003340.

2. F. Pind, A.P. Eising-Karup, C-H Jeong, J.S. Hesthaven, and J. Strømann-Andersen, “Time-domain room acoustic simulations with extended-reacting porous absorbers using the discontinuous Galerkin method,” J. Acoust. Soc. Am., vol. 148, p. 2851, 2020; doi.org/10.1121/10.0002448.

3. H. Wang and M. Hornikx, “Time-domain impedance boundary condition modeling with the discontinuous Galerkin method for room acoustics simulations,” J. Acoust. Soc. Am., vol. 147, p. 2534, 2020; doi.org/10.1121/10.0001128.

Acoustics_Module/Building_and_Room_Acoustics/wave_based_room

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Name	Expression	Value	Description
f0	700[Hz]	700 Hz	Signal center frequency
T0	1/f0	0.0014286 s	Signal period at center frequency
c0	343[m/s]	343 m/s	Speed of sound
lam0	c0/f0	0.49 m	Signal wavelength at center frequency

Interface	Output	Selection
Pressure Acoustics, Time Explicit (pate)	Selection