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

1 Introduction

In helicopters, cabin and cockpit heating is performed by taking hot air from engine compression stage. This heat is mixed to ambient or external air into a Heat and Ventilation Air Conditioning (HVAC) system through an injector in order to reach the temperature set point. The device produces a high frequency hot jet noise propagating through the air distribution system before it radiates into the cabin and cockpit. The whole architecture of the HVAC can be rather complex, see for instance the illustration of a typical cockpit air conditioning system in [1] (note this does not include the heating and ventilating part). HVAC systems are also encountered in many applications such as buildings and in the transport industry. In an early phase of designing HVAC duct networks, it is desirable to have a simulation tool for the prediction of the noise radiated from the terminal diffusers. For this purpose, HVAC systems are usually described as a network of several sub-systems joined by duct-like elements. Typical sub-systems are, for instance, bends and sudden expansion/contraction, T-junctions, duct or expansion chambers with acoustically treated walls and terminal diffusers. Each sub-element can act as an obstacle for the acoustic wave (i.e. the passive part) and also as a source of noise (i.e. the active part) which generally stems from the presence of a turbulent flow in the system. The development of simulation softwares has been the subject of active research in the last decades and we can notably cite [2, 3]. These can be regarded as “1-D” simulations tools using either the transfer matrix [3] or the scattering matrix formalism [2] and are therefore limited to the low-frequency range whereby only the plane wave mode is allowed to propagate in the ducts. In order to alleviate these limitations, a new active two-port scattering matrix formalism has been developed recently using the acoustic power as states variables instead of the acoustic pressure [4]. The theory relies on the algorithm for network equations developed earlier by Glav and Åbom [5] whereas scattering matrices as well as the source terms for each sub-element are constructed e.g. from VDI 2081 [6] and ASHRAE standards [7]. The method, although restricted to a specific choice of acoustic sub-elements, offers an efficient prediction tool for HVAC noise but ignores the effect of reflections. In the wake of Cumming’s work [8], the concept is extended by the same authors by constructing the scattering matrix via ray tracing algorithms [9]. This allows to include both reflection and transmission and to treat more accurately a variety of elements and good results are obtained as long as the “high-frequency assumption” remains valid, i.e. the number of propagating modes must be sufficiently large (say at least 10 as shown in [8]).

Many HVAC systems, such as those installed in helicopters, comprise non-standard multi-ports duct components involving complex geometries with or without absorbing liners. Moreover, the frequency spectrum of the source can largely exceed the plane wave cut-off frequency. Recall that this condition is reached as soon as the wavelength becomes smaller than the diameter of the duct. Thus, for a large part of the frequency range of interest, waves are likely to propagate in a multi-modal context. The propagation of high-order modes is accounted for in many research papers dealing with silencers of cylindrical shapes using the Mode Matching Methods [10, 11]. The technique relies on the decomposition of the acoustic waves in the silencer as a series of duct modes and the application of appropriate continuity conditions at both ends of the silencer. The whole system of equations which includes incident and reflected waves in the inlet and outlet ducts as well as in the silencer can be condensed in order to build the scattering matrix. Depending on the configuration, modes can be obtained analytically by exploiting the circular symmetry of the silencer [11] or can be calculated numerically when the silencer cross-section is elliptical [10]. However, in order to treat components with arbitrary geometries, the modal description is not valid anymore and the use of classical 3D discretization techniques such as the Finite Element Method (FEM) or the Boundary Element Method (BEM) is unavoidable [12].

The idea of connecting finite element matrices with modal representation of waves gave rise to the so-called “hybrid methods” which have been the research topic of many papers. In this regard, we can cite the work of Astley [13] who presented a general formulation for the computation of acoustic waves in exterior domains. Later, Kirby developed the technique for the propagation in acoustic waveguides [14]. An interesting feature that these hybrid methods just mentioned have in common is that finite element stiffness and mass matrices are not modified and could be built via any commercial software. Besides, the global matrix is symmetric (and possibly complex-valued if absorbing materials are present) regardless of the shape of the component. Again, the global matrix can be condensed after elimination of the internal degrees of freedom to obtain the scattering matrix.

It is the aim of the present paper to develop a scattering matrix formalism for the acoustic power in duct networks. The approach bears resemblance with [4] and [9] except that the matrix is derived from a general multi-port multi-modal scattering matrix (one application of the concept can be found in [15]) constructed from the 3D finite element discretization of each component of the network, thus allowing to take into account arbitrary geometries, with possibly the presence of absorbing materials, this is explained in Sections 2 and 3. In Section 4, the problem is then reformulated in terms of acoustic power by considering a given distribution of uncorrelated incident modes. One notes in passing that a somewhat similar idea is presented in [16] except a diffuse acoustic field loading is applied and the modal behavior is not considered. In the last section, the method is applied to the case of a simplified HVAC comprising Y-junctions, baffle-type silencers, and terminal diffusers. We should point out that the effect of air flow, either on the production or the propagation of acoustic waves (for the latter this not a major issue since the flow velocity is usually small in most HVAC systems) is neglected in this work and this could be the subject for further work.

2 Theory of the scattering matrices for N-ports acoustic systems

The theory of scattering matrices for acoustic systems, as illustrated in Figure 1 is briefly reminded here. Each port labeled l = 1, …, N consists of a straight duct with rigid walls and we introduce the associated local coordinate system (x_l, y_l, z_l) where x_l refers to the duct axis at port l. We call M_l the number of propagating acoustic duct modes and the acoustic pressure can be written as the modal series (convention exp(jωt)):

$$p_l\left(x_l,y_l,z_l\right)=\sum _{i=1}^{M_l}\mathrm{}{\varphi }_{i,l}\left(y_l,z_l\right)\left(A_i^l\exp \left(-j{k}_{i,l}{x}_l\right)+B_i^l\exp \left(j{k}_{i,l}{x}_l\right)\right).$$(1)

Figure 1 Illustration of an acoustic system Տ with N ports.

Here, it is understood that ϕ_i,l stands for the mode shape, k_i,l the (real-valued and positive) wavenumber of the duct mode propagating along the x_l direction, which is deliberately defined as pointing toward the acoustic system, and $A_i^l$ (resp. $B_i^l$) are the amplitudes of the incoming (resp. outgoing) waves. Here, mode shapes are function of the transverse coordinates (y_l, z_l) and satisfy classical orthogonality properties and are normalized such that

$${\int }_{{\mathrm{\Gamma }}_l}\mathrm{}{\varphi }_{i,l}\left(y_l,z_l\right){\varphi }_{i^{\text{'}},l}\left(y_l,z_l\right)\mathrm{d}\mathrm{\Gamma }={\delta }_{i,i^{\text{'}}},$$(2)

with dΓ = dy_ldz_l.

It is convenient for the subsequent analysis, to recast the previous equation in the compact form (here, x_l = 0 corresponds to the entrance of the system):

$$p_l={\underline{\mathsf{\Phi }}}_l^{\mathrm{T}}\left({\mathsf{A}}_l+{\mathsf{B}}_l\right), \hspace{1em}\mathrm{at} x_l=0,$$(3)

where the column vector ${\underline{\mathsf{\Phi }}}_l$ contain the mode shape functions at port l. Similarly, the acoustic velocity v_l, calculated along the x_l direction can be recovered via the Euler equation:

$$j{\rho }_l\omega v_l=-\frac{\mathrm{\partial }p}{\mathrm{\partial }{x}_l},\hspace{1em}\mathrm{at} x_l=0,$$(4)

and, using the orthogonality properties (2), it has the alternative form

$$v_l=\frac{1}{\omega }{\underline{\mathsf{\Phi }}}_l^{\mathrm{T}}{\mathsf{\Psi }}_l\left({\mathsf{A}}_l-{\mathsf{B}}_l\right)\hspace{1em}\mathrm{where} {\mathsf{\Psi }}_l=\frac{{\mathsf{k}}_l}{{\rho }_l}$$(5)

and ${\mathsf{k}}_l=\mathrm{diag}(\dots , k_{i,l}, \dots )$ is a diagonal matrix containing the duct acoustic wavenumbers. By exploiting the fact that the physical mechanisms taking place in the acoustic system are linear and time-invariant, we call $\mathsf{S}$ the scattering matrix which links incoming and outgoing waves:

$$\left(\begin{array}{l}{\mathsf{B}}_1\\ ⋮\\ {\mathsf{B}}_N\end{array}\right)=\mathsf{S}\left(\begin{array}{l}{\mathsf{A}}_1\\ ⋮\\ {\mathsf{A}}_N\end{array}\right).$$(6)

The matrix is composed of block diagonal matrices ${\mathsf{R}}_l={\mathsf{S}}_{l,l}$ which correspond to reflection of waves at port l and block matrices ${\mathsf{S}}_{l,l^{\text{'}}}$ account for the transmission from port l′ to port l:

$$\mathsf{S}=\left(\begin{array}{llllll}{\mathsf{S}}_{\mathrm{1,1}}& {\mathsf{S}}_{\mathrm{1,2}}& & {\mathsf{S}}_{1,l}& & {\mathsf{S}}_{1,N}\\ {\mathsf{S}}_{\mathrm{2,1}}& {\mathsf{S}}_{\mathrm{2,2}}& & {\mathsf{S}}_{2,}& & {\mathsf{S}}_{2,N}\\ ⋮& ⋮& \ddots & & & ⋮\\ {\mathsf{S}}_{l,1}& {\mathsf{S}}_{l,2}& & {\mathsf{S}}_{l,l}& & ⋮\\ ⋮& ⋮& & & \ddots & ⋮\\ {\mathsf{S}}_{N,1}& {\mathsf{S}}_{N,2}& & & & {\mathsf{S}}_{N,N}\end{array}\right).$$(7)

At port l, the total acoustic power is

$$W_l=\frac{1}{2}{\int }_{{\mathrm{\Gamma }}_l}\mathrm{}\mathrm{Re}\left(p_l{v}_l^{\mathrm{*}}\right)\mathrm{d}\mathrm{\Gamma }=W_l^{+}-W_l^{-},$$(8)

where acoustic power due to incoming and outgoing waves are obtained from

$$W_l^{+}=\frac{1}{2\omega }{\mathsf{A}}_l^{†}{\mathsf{\Psi }}_l{\mathsf{A}}_l\hspace{1em}\mathrm{and}\hspace{1em}{W}_l^{-}=\frac{1}{2\omega }{\mathsf{B}}_l^{†}{\mathsf{\Psi }}_l{\mathsf{B}}_l.$$(9)

For purposes of subsequent analysis, it is convenient to rewrite the last term using incoming waves and the scattering matrix of the system:

$$W_l^{-}=\frac{1}{2\omega }\sum _{l^{\text{'}}=1}^N\mathrm{}\sum _{l^{″}=1}^N\mathrm{}{\mathsf{A}}_{l^{\text{'}}}^{†}{\mathrm{S}}_{l,l^{\text{'}}}^{†}{\mathsf{\Psi }}_l{\mathsf{S}}_{l,l^{″}}{\mathsf{A}}_{l^{″}}.$$(10)

3 Construction of the scattering matrix using the finite element method

3.1 Formulation

The construction of the scattering matrix is achieved via the finite element discretization of the acoustic system. The latter is modeled as a three-dimensional cavity Ω of arbitrary shape filled with air which includes the presence of bulk reacting sound absorbing materials which can be treated using classical equivalent fluid models (see Fig. 2). A locally reacting wall with complex-valued impedance Z(ω), noted Γ_Z is also considered in the formulation. In the cavity, the acoustic pressure p must obey the following wave equation (see [17] for instance):

$$\nabla \cdot \left(\frac{1}{\rho }\nabla p\right)+\frac{{\omega }^2}{K}p=0.$$(11)

Figure 2 Acoustic system with N ports.

Here, the fluid density ρ and the bulk modulus K = ρc² must be viewed as functions of angular frequency ω and spatial coordinates (x, y, z).

After integration by parts, the weak formulation becomes (q stands for the test function):

$${\int }_{\mathrm{\Omega }}\mathrm{}\frac{\nabla p\cdot \nabla q}{\rho }\mathrm{d}\mathrm{\Omega }-{\omega }^2{\int }_{\mathrm{\Omega }}\mathrm{}\frac{\mathrm{pq}}{K}\mathrm{d}\mathrm{\Omega }+\mathrm{j\omega }{\int }_{{\mathrm{\Gamma }}_Z}\mathrm{}\frac{\mathrm{pq}}{Z}\mathrm{d}\mathrm{\Gamma }=\mathrm{j\omega }\sum _{l=1}^N\mathrm{}{\int }_{{\mathrm{\Gamma }}_l}\mathrm{}q{v}_l\mathrm{d}\mathrm{\Gamma }.$$(12)

To ease the demonstration, the acoustic pressure in the cavity is written with interpolating functions N_i, which results from the assembly process, and

$$p\left(x,y,z\right)=\sum _{i=1}^{N_E}\mathrm{}{N}_i\left(x,y,z\right)p_i={\mathsf{N}}^{\mathrm{T}}\mathsf{p},$$(13)

where N_E correspond to the number of nodes of the finite element mesh (here classical T10 quadratic elements are used) and the column vector $\mathsf{p}$ contains the value of the acoustic pressure at nodes. At port l, the same notation is used for the quantities of interest so we put

$$p_l={\mathsf{N}}_l^{\mathrm{T}}{\mathsf{p}}_l,\hspace{1em}{v}_l={\mathsf{N}}_l^{\mathrm{T}}{\mathsf{v}}_l,\hspace{1em} \mathrm{and}\hspace{1em} {\underline{\mathsf{\Phi }}}_l={\mathsf{\Phi }}_l^{\mathrm{T}}{\mathsf{N}}_l.$$(14)

Applying the previous notation to the variational form (12) leads to the algebraic system

$$\mathsf{Gp}=\sum _{l=1}^N\mathrm{}{\mathrm{\Gamma }}_l\left({\mathsf{B}}_l-{\mathsf{A}}_l\right)$$(15)

where

$$\mathsf{G}={\int }_{\mathrm{\Omega }}\mathrm{}\frac{\nabla \mathsf{N}\cdot \nabla {\mathsf{N}}^{\mathrm{T}}}{\rho }\mathrm{d}\mathrm{\Omega }-{\omega }^2{\int }_{\mathrm{\Omega }}\mathrm{}\frac{{\mathsf{NN}}^{\mathrm{T}}}{K}\mathrm{d}\mathrm{\Omega }+\mathrm{j\omega }{\int }_{{\mathrm{\Gamma }}_Z}\mathrm{}\frac{\mathsf{N}{\mathsf{N}}^{\mathrm{T}}}{Z}\mathrm{d}{\mathrm{\Gamma }}_Z,$$(16)

and

$${\mathsf{\Gamma }}_l=j{\int }_{{\mathrm{\Gamma }}_l}\mathrm{}{\mathsf{N}}_l{\mathsf{N}}_l^{\mathrm{T}}{\mathsf{\Phi }}_l{\mathsf{\Psi }}_l \mathrm{d}\mathrm{\Gamma }.$$(17)

Applying standard orthogonality properties and by multiplying by the quantity $j{k}_{i,l}/{\rho }_l$ yields

$${\mathsf{\Gamma }}_l^{\mathrm{T}}{\mathsf{p}}_l=j{\mathsf{\Psi }}_l\left({\mathsf{A}}_l+{\mathsf{B}}_l\right).$$(18)

Once assembled, the algebraic system takes the form of symmetric matrix (the vector ${\mathsf{p}}_{\mathrm{\Omega }}$ contains the value of the pressure in the domain Ω except at the ports):

$$\left(\begin{array}{lllllll}{\mathsf{G}}_{11}& & {\mathsf{G}}_{1N}& {\mathsf{G}}_{1\mathrm{\Omega }}& {\mathsf{\Gamma }}_1& 0& 0\\ ⋮& \ddots & ⋮& ⋮& 0& \ddots & 0\\ {\mathsf{G}}_{N1}& & {\mathsf{G}}_{\mathrm{NN}}& {\mathsf{G}}_{N\mathrm{\Omega }}& 0& 0& {\mathsf{\Gamma }}_N\\ {\mathsf{G}}_{\mathrm{\Omega }1}& & {\mathsf{G}}_{\mathrm{\Omega }N}& {\mathsf{G}}_{\mathrm{\Omega }\mathrm{\Omega }}& 0& 0& 0\\ {\mathsf{\Gamma }}_1^{\mathrm{T}}& 0& 0& 0& -j{\mathsf{\Psi }}_1& 0& 0\\ 0& \ddots & 0& 0& 0& \ddots & 0\\ 0& 0& {\mathsf{\Gamma }}_N^{\mathrm{T}}& 0& 0& 0& -j{\mathsf{\Psi }}_N\end{array}\right)\left(\begin{array}{l}{\mathsf{p}}_1\\ ⋮\\ {\mathsf{p}}_N\\ {\mathsf{p}}_{\mathrm{\Omega }}\\ {\mathsf{B}}_1\\ ⋮\\ {\mathsf{B}}_N\end{array}\right)=\left(\begin{array}{lll}{\mathsf{\Gamma }}_1& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\mathsf{\Gamma }}_N\\ 0& 0& 0\\ j{\mathsf{\Psi }}_1& 0& 0\\ 0& \ddots & 0\\ 0& 0& j{\mathrm{\Psi }}_N\end{array}\right)\left(\begin{array}{l}{\mathsf{A}}_1\\ ⋮\\ {\mathsf{A}}_N\end{array}\right).$$(19)

The scattering matrix is constructed by solving the previous system for each incoming mode with unit amplitude.

3.2 Numerical examples

In this section, two acoustic systems are considered, one concerns a parallel baffle-type silencer, the second is a 1-port system which simulates the radiation of acoustic waves from a flanged duct. These will serve to construct a simplified HVAC system as shown in the next section. The dissipative silencer is placed in a rigid rectangular duct of size 10 cm × 20 cm. It consists of a parallelepiped made of glass wool supported on both sides by two rigid plates. Parameters such as the density and bulk modulus are obtained thanks to the Johnson-Champoux-Allard model and properties of the glass wool are taken from [18]: tortuosity α_∞ = 1, resistivity σ = 14,066 N s m⁻⁴, porosity Φ = 0.954, viscous and thermal characteristics Λ = 91.2 μm and Λ′ = 182.4 μm. Air parameters are: air mass density ρ = 1.213 kg m⁻³, speed of sound c = 342.2 m s⁻¹, dynamic viscosity μ = 1.84 × 10⁻⁵ kg m⁻¹ s⁻¹, heat capacity ratio γ = 1.4, and Prandtl number Pr = 0.71 (see Fig. 3).

Figure 3 FEM model of baffle silencer housed in a rectangular duct.

Coefficients of the scattering matrix are compared to measurements carried out by R. Binois [18, 19]. Results of Figure 4 show the transmitted plane wave $B_1^2$ at port 2 for a given incident plane wave at port 1, with $A_1^1=1$. The reflected coefficients $B_1^1,B_1^2$, and $B_1^3$ are also shown. Note the first transverse mode is never excited due to symmetry reasons ($B_2^1\approx 0$) and a peak in the reflection coefficient for the plane wave is reached at the cut-off frequency, around 1714 Hz, which corresponds to the emergence of the symmetric second transverse duct mode. Computed results show good agreements in the frequency range of interest (up to 4000 Hz) and this confirms that the finite element discretization of the silencer is of sufficient quality.

Figure 4 Coefficients of the scattering matrix, measured (gray) and computed (black).

The construction of the scattering matrix as described earlier can be adapted to model infinite or semi-infinite media. For this purpose, the radiation of acoustic waves can be modeled using the Perfectly Matched Layer (PML). The latter, placed on the periphery of the computational domain, can be regarded as an artificial porous material which absorbs waves as they radiate away. The precise description of the method can be found in [12] and this will not be repeated here. In order to illustrate the method, we compute the radiation of acoustic waves from a flanged rigid duct of rectangular cross-section with dimensions 16 cm × 26.7 cm. The computational domain is shown in Figure 5. The entrance of the system, of rectangular shape, is identified in blue color.

Figure 5 The terminal diffuser modeled as a 1-port system.

We are interested in the Transmission Loss (TL) of this 1-port system for an incident plane wave of unit amplitude, $A_1^1=1$. Two scenarios are investigated: one with the presence of an outlet grille generally encountered in terminal diffusers and one without grille. The TL of the terminal diffuser can be estimated according to VDI 2081 [6] and is given by the following formula

$$\mathrm{TL}=10\log (1+(\frac{c}{4\pi f}{)}^2\frac{\mathrm{\Omega }}{S})+m(0.04283\log \left(f\sqrt{S}\right)-0.0303),$$(20)

where f = ω/2π, S is the outlet area, Ω = 2π and m is the length to height ratio. The computed TL is calculated as

$$\mathrm{TL}=-10\log \left(1-\frac{W_1^{-}}{W_1^{+}}\right).$$(21)

Results, from Figure 6, clearly show very good agreements (note that, up to 3 duct acoustic modes with coefficients $B_1^1,B_1^2,$ and $B_1^3$ are involved in the definition of the reflected power $W_1^{-}$). The “no grille” case shows noticeable deviations at low frequency (remind that the first cut-off frequency appears at 640 Hz).

Figure 6 Transmission loss of a terminal diffuser.

4 Scattering matrix for the acoustic power flow

4.1 Formulation and assumptions

The previous formalism allows an exact description of the scattering of waves taking place in each sub-system and in principle, this can be used to simulate the propagation of waves in arbitrary networks. In view of simplifying this, we can exploit the fact that, in many situations, many random excitation mechanisms are responsible for broadband noise. For incoherent excitation we can treat the incident mode amplitudes as uncorrelated random variables and it is assumed that

$$\langle A_{i,l}{A}_{i^{\text{'}},l^{\text{'}}}^{\mathrm{*}}\rangle =0\hspace{1em}i\ne i^{\text{'}},l\ne l^{\text{'}}.$$(22)

We call

$$\langle W_{l^{\text{'}}}^{+}\rangle =\frac{1}{2\omega }\langle {\mathsf{A}}_{l^{\text{'}}}^{†}{\mathsf{\Psi }}_{l^{\text{'}}}{\mathsf{A}}_{l^{\text{'}}}\rangle$$(23)

the incident acoustic power at port l′. We can calculate the radiated power at port l,

$$\langle W_l^{-}\rangle =\frac{1}{2\omega }\sum _{l^{\mathrm{\prime }}=1}^N\mathrm{}\langle {\mathsf{A}}_{l^{\mathrm{\prime }}}^{†}\mathrm{diag}\left({\mathsf{S}}_{l,l^{\mathrm{\prime }}}^{†}{\mathsf{\Psi }}_l{\mathsf{S}}_{l,l^{\mathrm{\prime }}}\right){\mathsf{A}}_{l^{\mathrm{\prime }}}\rangle ,$$(24)

where property (22) has been used. The concept of scattering matrix for the acoustic power flow starts by expressing the previous relation as

$$\langle W_l^{-}\rangle =\sum _{l^{\text{'}}=1}^N\mathrm{}{\widehat{S}}_{l{l}^{\text{'}}}\langle W_{l^{\text{'}}}^{+}\rangle$$(25)

where, by definition

$${\widehat{S}}_{l{l}^{\text{'}}}=\frac{\langle W_l^{-}\rangle }{\langle W_{l^{\text{'}}}^{+}\rangle }=\frac{\langle {\mathsf{A}}_{l^{\text{'}}}^{†}\mathrm{diag}\left({\mathsf{S}}_{l,l^{\text{'}}}^{†}{\mathsf{\Psi }}_l{\mathsf{S}}_{l,l^{\text{'}}}\right){\mathsf{A}}_{l^{\text{'}}}\rangle }{\langle {\mathsf{A}}_{l^{\text{'}}}^{†}{\mathsf{\Psi }}_{l^{\text{'}}}{\mathsf{A}}_{l^{\text{'}}}\rangle }.$$(26)

At this point, we should note that similar derivations can be found in [9] using the concept of acoustic exergy and this allows taking into account the effect of the flow in the analysis. The values of these scattering coefficients depend on the modal distribution and some additional assumptions must be adopted to simplify the analysis. For this purpose, we should adopt the power-law model

$$\langle |A_{i,l}{|}^2\rangle =Q(k_{i,l}{)}^{-\alpha }$$(27)

which relevance is discussed in [20]. Here, the strength Q is arbitrary and the exponent α depends on the nature of the source distribution model. Here α = 1 and α = 2 correspond, respectively, to equal energy per mode and a uniform distribution of incoherent monopoles and α = 0 to equal energy density per mode [20]. By construction, scattering coefficients for the acoustic power flow are positive and real-valued and satisfy the inequality

$$\sum _{l=1}^N\mathrm{}{\widehat{S}}_{l{l}^{\text{'}}}\le 1,$$(28)

and it is equal to unity in the no-dissipative case. In order formalize the theory for a acoustic network composed of n_s sub-systems, we call ${\widehat{\mathsf{S}}}^s$ the power scattering matrix associated with the sub-system number s ∈ [0, n_s], i.e. ${\widehat{S}}_{l{l}^{\text{'}}}^s=({\widehat{\mathsf{S}}}^s{)}_{l{l}^{\text{'}}}$.

4.2 Comparison between a single system and two sub-systems

The accuracy of the method is assessed in the case of two chained 2-ports sub-systems where each sub-system consists of an expansion chamber with or without porous material (see Fig. 7).

Figure 7 Illustration of an acoustic system with two ports (the FE mesh is represented). Scenario 1 (left) is referred to as “Monobloc” (in this example the length of the central duct is 10 cm) and scenario 2 (right) is referred to as “Chained”. Note the entrance of the system is identified by the blue color.

The main body length of the expansion chamber of cylindrical shape is set to L = 40 cm while inner and outer radius r_a and r_b are set to 4 cm and 6 cm, respectively. The chamber is either filled with air (purely reactive) or with a rock mineral wool placed between r_a and r_b (dissipative case). We are interested in computing the TL due to the association of two similar expansion chambers connected to each other.

$$\mathrm{TL}=-10\log \left(\frac{\langle W_t\rangle }{\langle W_i\rangle }\right),$$(29)

where 〈W_t〉 and 〈W_i〉 are the transmitted and incident acoustic power. Here, the simple case where 〈W_i〉 is constant is considered, this means that the losses can be evaluated directly from the knowledge of the scattering coefficients and $\mathrm{TL}=-10\log {\widehat{S}}_{21}$. Results called “Monobloc” correspond to the scenario 1 (see Fig. 7) where the system composed of the succession of two expansion chambers is modeled as single acoustic system. Results are compared with the “Chained” scenario whereby each expansion chamber is treated independently. Here, the scattering matrix for the whole system is calculated using the Redheffer star product [21]:

$$\widehat{\mathsf{S}}={\widehat{\mathsf{S}}}^1\otimes {\widehat{\mathsf{S}}}^2,$$(30)

which is equivalent to solving equation (25) for each sub-system. Figure 8 shows the TL due to the association of two rigid expansion chambers. Because this acoustic system which is purely reactive, is prone to resonance and phase effects between the two sub-systems, the Monobloc has been constructed by considering three different lengths for the central duct, here 10 cm, 20 cm, and 40 cm. Results are compared with the “Chained” scenario which does not depend on this geometric parameter. It emerges that all curves collapse showing a remarkable consistency as long as the frequency is sufficiently high. In order to interpret this more clearly, the number of propagating modes, call it M, in the connecting duct is shown in Table 1.

Figure 8 Transmission loss due to the association of two rigid expansion chambers (note the power flow model is constructed using equal energy per mode).

Results at low frequency are typical of the plane wave regime showing noticeable deviations whereas they become more consistent above 4000 Hz, which means that at least 5 modes propagate. Clearly, resonance effects which results in high losses due to the reflection of waves can only be partially captured by the “Chained” scenario. Once integrated over 1/3-octave bands¹, discrepancies are somehow smoothed out except in the plane wave regime and this also due to the fact that frequency bands are smaller and therefore very sensitive to frequency shifts (see Fig. 9).

Figure 9 Transmission loss due to the association of 2 rigid expansion chambers. Results are averaged over the 1/3-octave frequency band.

The effect of the modal distribution is now shown in Figure 10. Here, the scattering matrix for the power flow is constructed using the modal distribution law (27) with α = 0 and α = 1 as these values are representative of many configurations (see for instance [20] and references therein and in [22, 23]). Below the first cut-off frequency of the circular duct (around 2500 Hz), results are independent of the modal distribution as expected. It can be observed that differences between the two scenarios, though discernible, are modest and this remains the case even above cut-off. Above cut-off, results deviate slightly depending on the modal distribution but differences do not exceed 2 dB.

Figure 10 Transmission loss due to the association of two rigid expansion chambers.

Figure 11 shows the TL due to the association of one rigid expansion chamber and one filled with absorbing material (which properties are taken from the numerical example given in the Sect. 3.2). Deviations are only noticeable between 3000 Hz and 4000 Hz and they do not exceed 3dB which is reasonable when compared to the value of the TL, here above 45 dB of attenuation.

Figure 11 Transmission loss due to the association of one rigid expansion chamber and one filled with absorbing material.

5 Formalism for duct networks, case of a simplified HVAC

Two-port networks which are connected in series can be modeled using formula (30) which allows to compute rapidly the acoustic power scattering matrix of the whole system. In order to deal with more complex networks involving, for instance, subsystems having more than two ports, we need a more general formalism. This can be done in the spirit of Glav and Åbom’s work [5]. The approach proposes a strategy for the assembly of the scattering matrices taking into account the interactions between acoustic sub-systems. To illustrate this, the HVAC duct network of Figure 12 is chosen as an example.

Figure 12 Illustration of a simplified HVAC duct network with five acoustic sub-systems.

One can recognize the baffle-type silencer already encountered in Section 3 (sub-system 2) and the two terminal diffusers modeled as 1-port systems (labeled 4 and 5). Element 3 represents a typical Y-junction and the element 1 is a 1-port in which is placed the acoustic source. Systems are connected to each other via acoustic waveguides (rigid ducts) of arbitrary length. The starting point is the definition of the global scattering matrix $\widehat{\mathsf{S}}$ which results in the concatenation all scattering matrices involved. It is a block diagonal matrix having the following form:

$$\widehat{\mathsf{S}}=\left(\begin{array}{llllllll}{\widehat{S}}_{11}^1& 0& 0& 0& 0& 0& 0& 0\\ 0& {\widehat{S}}_{11}^2& {\widehat{S}}_{12}^2& 0& 0& 0& 0& 0\\ 0& {\widehat{S}}_{21}^2& {\widehat{S}}_{22}^2& 0& 0& 0& 0& 0\\ 0& 0& 0& {\widehat{S}}_{11}^3& {\widehat{S}}_{12}^3& {\widehat{S}}_{13}^3& 0& 0\\ 0& 0& 0& {\widehat{S}}_{21}^3& {\widehat{S}}_{22}^3& {\widehat{S}}_{23}^3& 0& 0\\ 0& 0& 0& {\widehat{S}}_{31}^3& {\widehat{S}}_{32}^3& {\widehat{S}}_{33}^3& 0& 0\\ 0& 0& 0& 0& 0& 0& {\widehat{S}}_{11}^4& 0\\ 0& 0& 0& 0& 0& 0& 0& {\widehat{S}}_{11}^5\end{array}\right).$$(31)

We now define the column vectors ${\mathsf{W}}^{-}$ and ${\mathsf{W}}^{+}$ containing all incoming and outgoing waves in the network. The presence of acoustic sources, generated for instance by the turbulent flow or a vibrating structure, can be easily accounted for by including a source vector ${\mathsf{W}}^S$ in the analysis, so the scattering mechanisms (for the acoustic power) in the network is described by the set of equations

$${\mathsf{W}}^{-}=\widehat{\mathsf{S}}{\mathsf{W}}^{+}+{\mathsf{W}}^S.$$(32)

The system of equations can be solved by realizing that incoming and outgoing acoustic power are equal at each interface. This means that

$${\mathsf{W}}^{-}=\mathsf{P}{\mathsf{W}}^{+},$$(33)where $\mathsf{P}$ is the permutation matrix

$$\mathsf{P}=\left(\begin{array}{llllllll}0& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right).$$(34)

The solution is easily found:

$${\mathsf{W}}^{+}=(\mathrm{P}-\widehat{\mathsf{S}}{)}^{-1}{\mathsf{W}}^S\hspace{1em}\mathrm{and}\hspace{1em}{\mathsf{W}}^{-}=\mathsf{P}(\mathsf{P}-\widehat{\mathsf{S}}{)}^{-1}{\mathsf{W}}^S.$$(35)

In this example, we consider that element 1 generates some acoustic power with a given sound power spectrum so we put ${\mathsf{W}}^S=(\langle W_1^S\rangle ,\mathrm{0,0},\mathrm{0,0},\mathrm{0,0},0{)}^T$. We are interested in the radiated acoustic power at ports 4 and 5. The TL is calculated via

$$\mathrm{TL}=-10\log \left(\frac{{\int }_{f_{\min }}^{f_{\max }}\mathrm{}\langle W^{\mathrm{rad}.}\rangle \mathrm{d}f}{{\int }_{f_{\min }}^{f_{\max }}\mathrm{}\langle W_1^S\rangle \mathrm{d}f}\right),$$(36)where 〈W^rad.〉 is the radiated power at both terminal diffusers, evaluated as follows

$$\langle W^{\mathrm{rad}.}\rangle =\langle W_4^{+}\rangle -\langle W_4^{-}\rangle +\langle W_5^{+}\rangle -\langle W_5^{-}\rangle .$$(37)

Results of Figure 13 show no more than 2 dB variations depending on the modal distribution and results below cut-off frequency of the silencer 2 (around 857 Hz) are identical as expected. These predictions also depend on the acoustic power spectrum of the source which is assumed constant in our calculations. Here, the impedance of the source is chosen so that ${\widehat{S}}_{11}^1=0$ which means that the source is placed in a semi-infinite duct and part of the acoustic energy is allowed to radiate upstream. The opposite scenario where there is total reflection, i.e. ${\widehat{S}}_{11}^1=1$ is also worth investigating. Results reported in Table 2 show that this condition can increase the radiated power by 0.5 dB to more than 2 dB depending on the frequency range. More importantly, the variations are comparable irrespective of the modal distribution law and this shows, at least in this example, that the method is robust and reliable for engineering purposes.

Figure 13 Transmission loss of the HVAC duct network.

We end this section by pointing out that the hypothesis that incoming waves obey a similar power-law model (27) might seem restrictive. As discussed in [9], dissipative elements are sensitive to the incident field and are also prone to strongly attenuate high-order modes. This can have an impact on the accuracy of the method when a large number of elements are connected in cascade [9]. To remedy this, scattered coefficients for the power (26) could be calculated using iterative algorithms (as shown in [11] in the case of a single expansion chamber with absorbing materials) by treating each element successively from the source to the exit of the duct system. This could be a topic for further studies.

6 Conclusion

The numerical method presented in this paper, enables to simulate the propagation of medium to high frequency sound waves, i.e. above plane wave cut-off frequency, in HVAC duct networks. The main idea is to construct a multi-port scattering matrix describing each component of the network using the acoustic power associated with incoming and outgoing waves, as a state variable instead of the acoustic pressure. The simplified model, relies on simplifying assumptions which are representative of broadband random noise: it is assumed that (i) duct acoustic modes are uncorrelated and (ii) the modal distribution is supposed to satisfy simple power-law models. Numerical tests show that the TL, once integrated over 1/3-octave bands, are not very sensitive to the power-law model. This method allows to tackle complex duct networks at a modest computational cost while taking into account the exact geometry of each multi-port element of the network. In practice, the method is shown to be an efficient and reliable tool for the study of dissipative systems when the frequency is sufficiently large, i.e. the propagation of waves in ducts must be multi-modal. It would be interesting to compare this with ray tracing algorithms based methods as advocated in [9] and we leave this for further studies. There are good reasons to believe that the new approach presented in this paper could be tailored to more specific configurations whereby propagating waves comprises both broadband noise and tonal components. In the latter case, scattering matrices could be constructed with a more appropriate model taking into account precisely the modal distribution and partial correlation between modes. The association of a large number of elements connected in cascade could be treated in an iterative manner as discussed earlier.

Conflict of interest

The authors declare no conflict of interest.

Funding information

The authors would like to acknowledge the financial support from the ANRT-Airbus Helicopters grant.

References

All Tables

All Figures

	Figure 1 Illustration of an acoustic system Տ with N ports.
In the text

	Figure 2 Acoustic system with N ports.
In the text

	Figure 3 FEM model of baffle silencer housed in a rectangular duct.
In the text

	Figure 4 Coefficients of the scattering matrix, measured (gray) and computed (black).
In the text

	Figure 5 The terminal diffuser modeled as a 1-port system.
In the text

	Figure 6 Transmission loss of a terminal diffuser.
In the text

	Figure 7 Illustration of an acoustic system with two ports (the FE mesh is represented). Scenario 1 (left) is referred to as “Monobloc” (in this example the length of the central duct is 10 cm) and scenario 2 (right) is referred to as “Chained”. Note the entrance of the system is identified by the blue color.
In the text

	Figure 8 Transmission loss due to the association of two rigid expansion chambers (note the power flow model is constructed using equal energy per mode).
In the text

	Figure 9 Transmission loss due to the association of 2 rigid expansion chambers. Results are averaged over the 1/3-octave frequency band.
In the text

	Figure 10 Transmission loss due to the association of two rigid expansion chambers.
In the text

	Figure 11 Transmission loss due to the association of one rigid expansion chamber and one filled with absorbing material.
In the text

	Figure 12 Illustration of a simplified HVAC duct network with five acoustic sub-systems.
In the text

	Figure 13 Transmission loss of the HVAC duct network.
In the text