A new approach to generate di ff use sound pressure fields

Cédric Van hoorickx 1 and Edwin P.B. Reynders 2

KU Leuven, Department of Civil Engineering, Structural Mechanics Section Kasteelpark Arenberg 40, box 2448, B-3001 Leuven, Belgium

ABSTRACT If the acoustic wavelength is small compared to the characteristic size of an enclosed space, the sound field in that space is often modelled as di ff use. A di ff use sound field is conventionally defined as a zero-mean circularly-symmetric complex Gaussian random field. A more recent, generalized definition is that of a sound field having mode shapes that are di ff use in the conventional sense, and eigenfrequencies that conform to the Gaussian Orthogonal Ensemble. Such a generalized di ff use sound field can represent a random ensemble of sound fields that share gross features such as modal density and reverberation time, but otherwise have any possible arrangement of local wave scattering features. In this contribution, realizations of a conventional di ff use sound field or, equivalently, of the mode shapes of a generalized di ff use sound field, are generated in a Monte Carlo framework. As a discrete decomposition is numerically expensive when the sound pressures at many locations are of interest, a fast analytical decomposition based on prolate spheroidal wave functions is developed. The approach is numerically validated by comparison with a detailed room model, where random wave scatterers are explicitly modeled, and good correspondence is observed. Applications involving correlated sound sources and sound-structure interaction are presented.

1. INTRODUCTION

At high frequencies, the sound field in an enclosed space is often modelled as di ff use. Conceptually, a di ff use field model represents a random ensemble of sound fields that share gross features such as modal density and total absorption, but otherwise have any possible arrangement of local wave scattering features such as boundaries and internal objects. The conventional mathematical definition of a di ff use field is that of a random field, composed of plane waves coming from all directions, with statistically independent and uniformly distributed phases and zero-mean, uncorrelated amplitudes. It follows immediately from the central limit theorem that such a conventional di ff use field is a circularly-symmetric complex Gaussian field [1]. Numerical and experimental studies have demonstrated that this conventional mathematical definition of a di ff use field is in line with the mean energetic quantities of the random ensemble that it is supposed to represent, but that this is not the case for higher order statistics such as the variance [2]. This is due to the fact that the modal behaviour of the enclosed space is neglected in the conventional definition. Accounting for this behaviour leads to the following definition: a

1 cedric.vanhoorickx@kuleuven.be

2 edwin.reynders@kuleuven.be

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

Figure 1: Three di ff use field realizations of a rectangular room of size 4 . 15 m × 5 . 09 m × 4 . 12 m at 250 Hz in a fluid with ρ a = 1 . 2 kg / m 3 and c a = 343 m / s.

generalized di ff use sound field is generated by an enclosure that has mode shapes that are di ff use in the conventional sense [3], and a squared eigenfrequency spacings distribution that conforms to the Gaussian orthogonal ensemble (GOE) eigenvalue spacings distribution [4, 5]. This generalized di ff use field model is known to be valid for all vibro-acoustic systems for which the uncertainty of the local field quantity due to the small random wave scattering elements, is large [6]. When a component in a vibro-acoustic system is modelled as di ff use, usually only its total energy is considered and additional assumptions are made, e.g., regarding the coupling to other components. The limitations of this approach can be overcome in a Monte Carlo analysis, where samples of a di ff use wave field are generated numerically by direct application of its definition [7, 8]. Since a conventional di ff use field (hence also the mode shapes of a generalized di ff use field) is Gaussian, it is completely determined by the spatial correlation function. If the sound pressures in a finite number of positions are of interest, they can be collected into a random vector. Monte Carlo samples of the random vector can then be generated via an eigenvalue decomposition of the correlation matrix. As an example, Figure 1 displays three realizations of the di ff use sound pressure field in a rectangular room that have been generated in this way. However, this approach becomes computationally expensive when sound field realizations at a large number of points are required, which is the case for many problems in room acoustics and vibro-acoustics, as the size of the correlation matrix grows with the square of the number of points. An alternative method is therefore developed. The basic idea is to perform the eigenvalue decomposition of the spatial correlation function without any spatial discretization, resulting in eigenvalues and continuous eigenfunctions, rather than discrete eigenvectors. Such a decomposition is known as a Karhunen-Loève decomposition. Its computational cost is independent of the number of locations at which the random field realizations need to be generated. The spatial correlation function of a three-dimensional conventional di ff use sound field is a sinc function [3, 9]. When the locations of interest are situated on a line, the Karhunen-Loève decomposition problem relates to the solution of the Helmholtz equation in spheroidal coordinates and can be described in terms of prolate spheroidal wave functions (PSWFs). In this work, the PSWF theory is invoked for obtaining the Karhunen-Loève decomposition of the spatial correlation function, such that Monte Carlo samples of the di ff use sound pressures at a large number of locations can be e ffi ciently generated. First, the special case is studied where the locations of interest are situated on a line. Subsequently, the general case, where the points of interest are at arbitrary locations in the three-dimensional sound field, is tackled by introducing an approximation of the spatial correlation function that renders the eigenproblem separable along perpendicular axes. Finally, the approach is validated in applications where the sound field needs to be evaluated at many locations, such as the evaluation of the sound field in a room with a loudspeaker array or the interaction between a sound field and a vibrating surface.

2. MONTE CARLO SIMULATION OF A CONVENTIONAL DIFFUSE SOUND FIELD

A conventional di ff use field is a superposition of many elementary plane waves coming from all directions. The amplitudes and phases of the elementary contributions are statistically independent

and the phases are normally distributed in [ − π, π ]. Let us consider a harmonic sound field at frequency ω with complex pressure amplitude p ( x ) such that the pressure at position x and time t is given by Re  p ( x ) exp(i ω t ) . Assuming small fluid damping, it can be shown that such a di ff use pressure field is circularly-symmetric complex Gaussian with the following covariance function [1]:

C p x , x ′
= Cov p ( x ) , p ( x ′ )
= K p j 0 k || x − x ′ ||
, (1)

where j 0 ( x ) = sin( x ) / x denotes the spherical Bessel function of the first kind and order zero, k = 2 π

λ the wavenumber corresponding to the wavelength λ , and K p a scaling factor which depends on the excitation. For a room with an acoustic source with a known average input power E [ W in ], this equals

K p = 2 ρ a c 2 a ωη V E [ W in ] , (2)

with V the volume of the acoustic enclosure, η the damping loss factor, and ρ a and c a the air density and sound speed, respectively. Monte Carlo samples of the complex pressure amplitude p ( x ) in a domain of interest Ω ′ can be obtained from a truncated Karhunen-Loève decomposition [10]. In this approach, the covariance function is decomposed as C p x , x ′
= X

n λ n f n ( x ) f n ( x ′ ) , (3)

where λ n and f n ( x ) are the eigenvalues and eigenfunctions, respectively, of the following Fredholm equation of the second kind with kernel C p ( x , x ′ ): Z

Ω ′ C p ( x , x ′ ) f n ( x ′ ) d x ′ = λ n f n ( x ) . (4)

The pressure field can then be expressed as:

p ( x ) = X

p

λ n f n ( x ) ξ C n (5)

n

with ξ C n =  ξ R n + i ξ I n  / √

2 where ξ R n and ξ I n denote independent, standard normal random variables that can be realized with a Gaussian random number generator. The eigenproblem in Equation 4 can be solved numerically for any covariance function by discretizing the domain of interest Ω ′ . A drawback of this discrete Karhunen-Loève decomposition is that the computational cost increases if the modeshape component at many correlated points is required [11]. Some problems in which this is the case are: the computation of the detailed pressure field in the vicinity of one or more point sources, of the di ff use sound pressure loading onto a mechanical structure, and of the sound radiation from a vibrating structure. For the latter two, we need the pressure distribution at the entire interface between structure and room. In this contribution, we therefore present an alternative approach using an analytical solution of the eigenproblem.

3. DIFFUSE PRESSURE FIELD ON A LINE

If the domain of interest Ω ′ is a line, i.e. Ω ′ = [ − a , a ], the eigenproblem (Equation 4) transforms into

Z a

sin ( k ( x − x ′ ))

K p

x − x ′ f n ( x ′ )d x ′ = λ n f n ( x ) . (6)

k

− a

From the solution of this eigenproblem, yielding eigenvalues λ n and eigenfunctions f n ( x ), the pressure is obtained using Equation 5. Substitution of x = au and x ′ = au ′ in Equation 6 yields Z 1

sin ( ak ( u − u ′ ))

π ( u − u ′ ) ψ n ( u ′ )d u ′ = µ n ψ n ( u ) , (7)

− 1

where

ψ n ( u ) = ψ n  x

 = f n ( au ) = f n ( x ) and µ n = k π K p λ n . (8)

a

It has been demonstrated that the eigenfunctions ψ n of Equation 7 correspond to the solutions of the following Sturm-Liouville di ff erential equation with eigenvalues χ n [12]:

 1 − u 2  d ψ n

!  χ n − a 2 k 2 u 2  ψ n = 0 . (9)

d d u

d u

This di ff erential equation has exactly the same form as the one that defines the radial and angular prolate spheroidal wave functions R mn ( ak , ξ ) and S mn ( ak , η ) for the special case m = 0. These di ff erential equations are encountered when the Helmholtz equation ∇ 2 Θ + k 2 Θ = 0 is separated in spheroidal coordinates ξ ∈ [1 , ∞ ), η ∈ ( − 1 , 1), and ϕ ∈ [0 , 2 π ), i.e. by considering the solution to be Θ ( ξ, η, ϕ ) = R mn ( ak , ξ ) S mn ( ak , η ) Φ m ( ϕ ) where 2 a is the interfocal distance of the elliptical cross section [13]. The eigenfunctions of Equation 9 (and Equation 7) therefore equal the angular prolate spheroidal wave functions for m = 0 [12], and the eigenfunctions f n ( x ) in Equation 6 are

f n ( x ) = ψ n  x

 = S 0 n  ak , x

 . (10)

a

a

As an example, Figure 2 displays the first three eigenfunctions f n ( x ) along a line with length 2 a = 4 . 15 m at a frequency of 400 Hz. The eigenfunctions obtained from the prolate spheroidal wave function (PSWF) solution, Equation 10, agree perfectly with the ones obtained with a discrete Karhunen-Loève decomposition for 416 points equally distributed along the line.

0.2

0.2

0.2

0.5

Eigenfunction 0

Eigenfunction 1

Eigenfunction 2

0.4

0.1

0.1

0.1

Eigenvalue

0.3

0

0

0

0.2

-0.1

-0.1

-0.1

0.1

-0.2

-0.2

-0.2

0

0 1 2 3 4 x coordinate [m]

0 1 2 3 4 x coordinate [m]

0 1 2 3 4 x coordinate [m]

0 10 20 30 40 Order n

Figure 2: Discrete decomposition (blue solid line) and PSWF (red dashed line) solution of the first three eigenfunctions f n ( x ) and the eigenvalues λ n at 400 Hz.

It has furthermore been demonstrated that the eigenvalues µ n in Equation 7 and therefore the corresponding eigenvalues λ n of Equation 6 can be related to the radial prolate spheroidal wave functions [12]:

λ n = π K p

k µ n = 2 aK p [ R 0 n ( ak , 1)] 2 , (11)

where R mn ( ak , ξ ) is the radial prolate spheroidal wave function normalized such that for large ξ

R mn ( ak , ξ ) → 1 ak ξ cos ak ξ − n + 1

2 π ! . (12)

Figure 2 also displays the eigenvalues λ n along a line with length 2 a = 4 . 15 m at a frequency of 400 Hz for K p = 1 Pa 2 . The eigenvalues obtained from the PSWF solution, Equation 11, are compared with the ones obtained with a discrete Karhunen-Loève decomposition. Only the first eigenvalues are nonzero, showing that if the response in multiple degrees of freedom is searched, the eigenproblem to be solved for the PSWF solution is smaller than the one for the discrete decomposition. After computing the eigenfunctions and eigenvalues by evaluating Equation 10 and Equation 11, the pressure field can be finally obtained from Equation 5.

4. DIFFUSE PRESSURE FIELD ON A SURFACE OR VOLUME

4.1. Theory In what follows, a generalization towards multiple dimensions is achieved. We propose an approximation of the correlation function for a domain of interest Ω ′ in multiple dimensions as a product of correlation functions along multiple perpendicular lines. For a flat surface, the correlation function can be approximated by:

C p x , x ′
= K p j 0 ( k || x − x ′ || ) ≈ K p j 0 ( k | y − y ′ | )j 0 ( k | z − z ′ | ) , (13)

where x = { y , z } T and x ′ = { y ′ , z ′ } T . Only for large values of both k | y − y ′ | and k | z − z ′ | or for | y − y ′ | , | z − z ′ | ≫ λ with λ the wavelength, the approximation error becomes significant. The correlation function itself, however, then becomes small. Indeed, if x and x ′ are far apart, the sound pressures are approximately uncorrelated, justifying the proposed approximation. For a rectangular domain Ω ′ = [ − a y , a y ] × [ − a z , a z ], Equation 4 then reduces to:

Z a y

Z a z

sin( k ∆ y )

k ∆ y sin( k ∆ z )

k ∆ z f n ( y ′ , z ′ )d y ′ d z ′ = λ n f n ( y , z ) . (14)

K p

− a y

− a z

This eigenproblem is separable in y and z :

a y k , y

! S 0 n z

a z k , z

! , (15)

f n ( y , z ) = f n y , n z ( y , z ) = f n y ( y ) f n z ( z ) = S 0 n y

a y

a z

λ n = λ n y , n z = λ n y λ n z = 4 a y a z K p h R 0 n y ( a y k , 1) i 2  R 0 n z ( a z k , 1)  2 . (16)

Although these expressions are obtained for a rectangular solution domain, the shape of the domain is of no importance as we are interested in di ff use pressure components, which are the same regardless the chosen (sub)domain. This means that for a non-rectangular domain, the solutions for a circumscribed rectangle or bounding box can be used. Furthermore, with the procedures used, i.e. decomposing the solution into the ones of the eigenproblems along several lines, also the modeshapes in three dimensions can be obtained, i.e. for rectangular cuboids.

4.2. Application: di ff use sound loading Pressure field realizations are required when analyzing the loading onto a structure by a di ff use sound pressure field. This pressure field is computed here assuming light fluid loading, for which the fluid- structure coupling is weak [14]. The out-of-plane displacements of the structure at the positions x ′ are obtained through a truncated modal decomposition. Considering the structure as a perfectly reflecting boundary, the mean squared sound pressure is twice the mean squared sound pressure in the center of the room [15]. Because of this, the covariance function (Equation 1) should be multiplied by two for points located on the reflecting structural element [16]. As an example structure, consider a flat steel plate with dimensions 0 . 5 m × 0 . 5 m and thickness 5 mm. The material properties of the plate are: a mass density of 7850 kg / m 3 , a Young’s modulus of 210 GPa, a Poisson’s ratio of 0 . 25 and damping loss factor of 0 . 1. The steel plate is coupled to a room of dimensions 4 . 15 m × 5 . 09 m × 4 . 12 m at x = 4 . 15 m. The room has a damping loss factor of 0 . 01 and contains a point monopole at position { 1 . 9 m , 2 . 3 m , 1 . 1 m } T with a unit harmonic volume acceleration a p . The room has a random acoustic mass distribution in the sense that a total of 50 point air pockets (or acoustic point masses) are distributed at random and statistically independent locations within the room. The air pockets are highly idealized models for small wave scatterers in the room, in the same way as point masses would be highly idealized models for small wave scatterers on a plate. Each air pocket has 0 . 3 % of the total acoustic mass V / c 2 a of the room.

The random air pockets are explicitly modeled in a detailed MC approach, while they are implicitly accounted for in the PSWF approach by modelling the room as being di ff use. The detailed MC results are achieved with an assumed-modes model as described in Appendix A of [17], in which the hard- walled room modes are taken as basis functions. Realizations of the random sound pressure are computed by fixing the positions of the point air pockets with a random number generator. This process is repeated to generate a Monte Carlo ensemble of 1000 realizations of the sound field in the room and the displacement field in the plate.

10 -4

10 -6

Mean total room energy [J]

Mean total plate energy [J]

10 -5

10 -8

10 -6

10 -10

10 -7

10 -12

125 250 500 Frequency [Hz]

125 250 500 Frequency [Hz]

(a)

(b)

10 -8

10 -10

Variance total room energy [J 2 ]

Variance total plate energy [J 2 ]

10 -10

10 -15

10 -12

10 -20

10 -14

10 -16

10 -25

125 250 500 Frequency [Hz]

125 250 500 Frequency [Hz]

(c)

(d)

Figure 3: (a,b) Mean and (c,d) variance of the total energy in (a,c) the room and (b,d) the plate. Results for the detailed MC approach are plotted with blue solid lines, for the conventional di ff use approach in black dash-dotted lines and for the GOE di ff use approach in red dashed lines.

Figure 3 displays the mean and variance of the total energy in the room and the plate computed with the detailed MC and the presented PSWF approach. For the mean energies, a good correspondence between both is obtained. Although the correspondence for the variance of the total plate energy is good, this is not the case for the variance of the total room energy: its variance is substantially underestimated. This is because the sound field is produced by a point source, and therefore many eigenmodes contribute to the response. In a single frequency treatment of the di ff use field, a circularly-symmetric complex Gaussian pressure is assumed, implying that the standard deviation of the squared pressure modulus equals its mean value. Previous studies have shown that this is generally not the case and have illustrated that an accurate variance is obtained when assuming Gaussian Orthogonal Ensemble (GOE) statistics [18].

5. MONTE CARLO SIMULATION OF A GENERALIZED DIFFUSE SOUND FIELD

5.1. Generation of modeshapes in a di ff use field To obtain the correct variance of energetic quantities, a generalized definition of a di ff use field needs to be used, relying on the fact that the statistics of the local eigenvalue spacings are those of the GOE, while the modeshapes are Gaussian random fields. This has led to the development of the GOE-MC method [7, 8]. We will now apply the previously discussed Karhunen-Loève decomposition with prolate spheroidal wave functions to obtain these di ff use modeshapes. For the generation of natural frequencies with a correct eigenvalue density, we refer the reader to [8].

Adopting a di ff use wave field model, the modeshape of a system component in high-frequency regime at a given location consists of a summation of independent plane waves with the same mean amplitude and uncorrelated phases, coming from all directions with equal probability. It then follows from the central limit theorem that the modeshapes are zero-mean, Gaussian random fields. For di ff usely reflecting boundaries, the modeshapes ϕ p i are statistically homogeneous, i.e., the statistics of the modeshape components are independent of their position. The covariance function for an acoustic enclosure then has the form [3,9]:

C ϕ i x , x ′
= Cov  ϕ p i ( x ) , ϕ p i ( x ′ )  = K ϕ j 0 ( k p i || x − x ′ || ) , (17)

where j 0 ( x ) = sin( x ) / x is the spherical Bessel function of the first kind and order zero, k p i = ω p i

c a denotes the wavenumber corresponding to the eigenfrequency ω p i of mode i , and K ϕ is a factor which can be determined from the modeshape normalization condition. For a homogeneous three-dimensional acoustic enclosure Ω with volume V , the normalization condition reads Z

1 c 2 a ϕ 2 p i ( x )d x = 1 ⇔ K ϕ = c 2 a V . (18)

As the correlation function of the modeshapes (Equation 17) and the pressure field itself (Equation 1) have the same form, they can both be solved with a continuous Karhunen-Loève decomposition with the eigenvalues and eigenfunctions obtained from prolate spheroidal wave functions. The main di ff erence is that the modeshapes are real Gaussian (instead of circularly- symmetric complex Gaussian) and therefore ξ C n in Equation 5 should be replaced by real standard normal random variables.

5.2. Application: di ff use sound field loading Let us return to the problem of section 2, i.e. the loading of a plate by a di ff use sound field. The pressure is obtained through modal decomposition into its di ff use eigenmodes. Considering the plate as a perfectly reflecting boundary, the covariance function in Equation 17 is multiplied by two for points located on the reflecting structural element [16]. Figure 3 displays the mean and variance of the total room and plate energy. The variance of the total room energy now agrees well with the ones predicted by the detailed MC approach. The uncertainty is a bit higher for the di ff use approach than for the detailed MC approach at low frequencies, where the sensitivity of the local vibration field to the presence of the small acoustic point masses is relatively low and the di ff use field model is not yet entirely valid. This example illustrates that the GOE approach is capable to correctly predict the energetic variance in a di ff use field. The normalization factor of the modeshape covariance function is furthermore independent from the mean energy in the room, in contrast to the pressure covariance function, for which a power balance is needed.

6. CONCLUSIONS

Di ff use sound field realizations from a numerical solution of the eigenvalue decomposition of the covariance matrix become very expensive for a large number of correlated degrees of freedom. An analytical approach based on prolate spheroidal wave functions has therefore been presented to generate di ff use sound fields. The analytical approach yields the exact solution for correlated degrees of freedom situated on a line, while an approximate solution can be obtained for degrees of freedom situated on a surface or in a volume. The approach for obtaining di ff use sound fields can be generalized in a GOE-MC framework, where now the di ff use mode shapes are generated using prolate spheroidal wave functions. The approach is numerically validated by comparison with a detailed room model, where random wave scatterers are explicitly modelled as acoustic point masses with random positions, and good correspondence is observed.

ACKNOWLEDGEMENTS

The research presented in this paper has been performed within the frame of the VirBAcous project (project ID 714591) “Virtual building acoustics: a robust and e ffi cient analysis and optimization framework for noise transmission reduction” funded by the European Research Council in the form of an ERC Starting Grant. The financial support is gratefully acknowledged.

REFERENCES

[1] K.J. Ebeling. Statistical properties of random wave fields. In W.P. Mason and R.N. Thurston, editors, Physical acoustics Vol. XVII , pages 233–310. Academic Press, Orlando, FL, 1984. [2] R.S. Langley and A.W.M. Brown. The ensemble statistics of the energy of a random system subjected to harmonic excitation. J Sound Vib , 275(3–5):823–846, 2004. [3] M. V. Berry. Regular and irregular semiclassical wavefunctions. J Phys A-Math Gen , 10(12):2083–2091, 1977. [4] R.L. Weaver. On the ensemble variance of reverberation room transmission functions, the e ff ect of spectral rigidity. J Sound Vib , 130(3):487–491, 1989. [5] M.L. Mehta. Random Matrices . Elsevier, San Diego, CA, 3rd edition, 2004. [6] C. Ellegaard, K. Schaadt, and P. Bertelsen. Acoustic chaos. Physica Scripta , T90:223–230, 2001. [7] E. Reynders, J. Legault, and R.S. Langley. An e ffi cient probabilistic approach to vibro-acoustic analysis based on the Gaussian orthogonal ensemble. J Acoust Soc Am , 136(1):201–212, 2014. [8] C. Van hoorickx and E.P.B. Reynders. Gaussian orthogonal ensemble modeling of built- up systems containing general di ff use components and parametric uncertainty. J Sound Vib , 501(116045):1–18, 2021. Open access. [9] R.K. Cook, R.V. Waterhouse, R.D. Berendt, S. Edelman, and M.C. Thompson Jr. Measurement of correlation coe ffi cients in reverberant sound fields. J Acoust Soc Am , 27(6):1072–1077, 1955. [10] R.G. Ghanem and P.D. Spanos. Stochastic finite elements: a spectral approach . Dover Publications, New York, 2nd edition, 2003. [11] M. Eiermann, O.G. Ernst, and E. Ullmann. Computational aspects of the stochastic finite element method. Comput Vis Sci , 10(1):3–15, 2007. [12] D. Slepian and H.O. Pollak. Prolate spheroidal wave functions, Fourier analysis and uncertainty – I. AT & T Tech J , 40(1):43–63, 1960. [13] E. Skudrzyk. The foundations of acoustics: basic mathematics and basic acoustics . Springer, Vienna, Austria, 1971. [14] N. Atalla and F. Sgard. Finite element and boundary element methods in structural acoustics and vibration . CRC Press, Boca Raton, FL, 2015. [15] R.V. Waterhouse. Interference patterns in reverberant sound fields. J Acoust Soc Am , 27(2):247– 258, 1955. [16] J.H. Rindel. Sound insulation in buildings . CRC Press, Boca Raton, FL, 2018. [17] E.P.B. Reynders, P. Wang, C. Van hoorickx, and G. Lombaert. Prediction and uncertainty quantification of structure-borne sound radiation into a di ff use field. J Sound Vib , 463:114984, 2019. [18] R.S. Langley and V. Cotoni. The ensemble statistics of the vibrational energy density of a random system subjected to single point harmonic excitation. J Acoust Soc Am , 118(5):3064– 3076, 2005.