Prediction of acoustic transmission in heavily damped system using hybrid Ray-Tracing-SEA method

Classic statistical energy analysis mainly deals with the energy transmission of system with relatively low damping. With the application of passive damping treatments, one of the fundamental assumptions in SEA, i.e., a diffused field, tends to fail. The energy attenuation along transmission path becomes so significant that there may exist large energy level difference within one structural component. In light of this, this study proposes to use a hybrid ray-tracing-SEA method to predict the acoustic transmission in heavily damped system. Heavily damped structural components are treated as “coupling elements” instead of normal “subsystems.” The energy transmission from one structural element to a connected structural element through the edge can be represented by using certain number of point sources and assuming each point source radiates certain number of acoustic rays. By tracing the travelling history of each ray, the energy attenuation along the travelling path can be achieved. With the information of energy input and energy attenuation, the equivalent coupling loss factors can be computed. By rebuilding a hybrid raytracing-SEA model, the energy level differences between different subsystems can be determined. Numerical validation of the ray tracing algorithm is conducted by comparing the calculated coupling loss factor with normal SEA method. Numerical study of a oneroom system is given. The room is assumed to consist of six homogeneous concrete plates and the analysis assumes no coupling between in-plan and outplane waves. Comparisons between classic and hybrid method show that when a small number of the structural components are heavily damped, classic SEA gives similar results with the hybrid method because the prediction errors tend to cancel with each other and the transmission is dominated by paths that are not damped. With the increasing number of damped elements, SEA tends to underestimate the energy level difference.


INTRODUCTION
For acoustic and vibrational energy transmission prediction, there generally exist two approaches, i.e., the deterministic approach and the statistical approach. The deterministic approaches such as finite element method (FEM) and boundary element method (BEM) have been demonstrated to work well for simple system at low frequency range. The statistical approach, such as statistical energy analysis (SEA), is known to be reliable at mid-and high-frequency range. There are highand low-frequency limitations for those two approaches respectively. In order to achieve a reliable broad band simulation result, researchers proposed some hybrid deterministic-statistical approaches [1][2][3][4] . The general idea is to use deterministic method to simulate the system behavior at low-frequency range where the system exhibits deterministic manner with low modal density and to use statistical method to predict the system behavior at mid-and high-frequency range where the system responses in a statistical manner with high modal density and statistical overlap. These methods only concern normal system with low damping and can be seen as extensions of SEA modeling.
Statistical energy analysis has been widely used in predicting mid-and high-frequency vibrational energy transmission in complex systems such as buildings, automobiles, ships, and satellites 5 . One of the fundamental assumptions in SEA is a diffused field, i.e., independent sound waves with the same strength arriving at a receiver from all directions 6 , or in other words, energy is evenly distributed within a subsystem. This assumption tends to fail when the damping level (referred as the internal loss factor in SEA) increases. With higher damping level, the energy attenuation within a physical enclosure (i.e., ceiling, wall or floor) becomes increasingly significant that there might exist large energy level difference. Similar findings have been reported in large walls and long spaces such as corridors and train cabins 1 . Passive damping treatments, i.e., free layer, constrained layer, granular damping etc., are generally used in acoustic and vibrational control. The loss factor of well-designed passive damping treatment can reach the level as high as 0.21 to 0.53 [7][8] . The reliability of conventional SEA modeling as well as the hybrid deterministic-SEA method is questioned for systems consisting of elements with such high damping level.
Ray tracing method was first developed for computational image synthesis by tracking the travelling history of different rays [9][10] . Its application was later extended to sound field simulation [11][12] , acoustic visualization 13 , and auralization [13][14] . Recent studies have shown its application in improving the accuracy of simulated results in periodic box-like structures 15 and long spaces 16 . The former study was limited to lightly damped system, i.e., the internal loss factors (ILFs) of the numerical model examined are 0.01 and 0.005. The latter study only concerned the energy transmission in space with a point excitation. Yan and Wilson 17 examined the energy transmission along heavily damped plate using ray tracing method and proposed to use the concept of effective length ratio as an indicator of heavily damped system. However, it only concerned very simple system which is rare in real applications.
In this article, a hybrid ray-tracing-SEA method is proposed to examine the acoustic energy transmission within heavily damped system. It uses ray tracing to account for the behavior of heavily damped components within the framework of SEA. The energy transmission along these damped elements can be determined using ray trace algorithm. This information can then be substituted back into SEA modeling to work out the energy level difference between the source and the receiver. This method is suitable for vibroacoustic prediction of heavily damped system at mid-and high-frequency range.
Introduction of the theory and rationale of ray tracing will first be given. Numerical validation of the ray tracing algorithm will be followed. Comparisons and discussions between classic SEA and the hybrid approaches are given for a concrete one-room system to understand the implication on performance.

THE RAY TRACING ALGORITHM
One may use certain number of point sources to simulate a line source on the edge of a plate assuming each point source radiates certain number of acoustic rays as indicated in Fig. 1. By tracking the travelling history and taking account of the energy attenuation along the path and on the edge of each ray, the combined results can give the energy transmission along each edge. The point sources are assumed to have the same amount of energy and to be evenly distributed along the edge. Each point source is assumed to radiate same number of acoustic rays that are also evenly distributed. Homogeneous thin structure and specular reflections at edges are assumed.
Assuming a heavily damped plate with length a and width b is placed in a 2D Cartesian coordinate system as shown in Fig. 1, for a given point source (x0, y0), the equation of the acoustic ray generated with a start angle θ is given by: where 0 < x0 < a, y0 = 0. The start angle θ is always de-fined in this equation as the angle between the x-axis and the ray measured in an anti-clockwise direction. The first intersection can be calculated by solving four sets of equations: If the solution (x1,y1) satisfies that and then solution (x1, y1) is accepted as the first intersection.
The second intersection can simply be found by conducting the same procedure where (x1,y1) is used as the new start point and the new start angle is (180 -θ). The intersections at the four corners of the plate are excluded in this study for simplicity and only specular reflection is considered. Once the coordinates of the intersections are determined, the distance the ray travelled from the origin to the Nth intersection can be obtained using: The incident energy is derived using wave propagation theory. The complex wavenumber for pure bending waves in a homogeneous isotropic plate is: where kB ⊥ represents the real bending wavenumber and the imaginary bending wavenumber kB ⫫ = 0:25kB ⊥ .  is the loss factor. For plane waves propagating in the positive x-direction, the velocity can be expressed as 7 : Making use of the complex wave number results in This implies exponential decay of propagating plane waves. The amplitude of particle velocity decays by a factor of e -k B ⫫ . For bending waves, kB ⫫ = =2B ; thus, the amplitude decays by a factor of . The bending wave energy can be expressed as EB = mu 2 (x, t). Hence, the amplitude of energy decays by a factor of . Within a distance Δx, there occurs an energy reduction in level by: When the bending wave starts with u(x , t ) and travels with a distance L, the travelling time is . Thus, the ratio of the initial bending wave energy, i.e., E0, to the energy when the wave arrives x = x0 + L, i.e., EL, is: This gives: where the real part of bending wave number is the bending wave speed in plates the bending stiffness for plates, as given by Craik 5 .
Each time a ray is incident with the boundary, certain amount of energy is transmitted if the boundary is connected to another plate. The proportion of the transmitted energy to the incident energy is called the transmission co-efficient, i.e., t. It depends on the incident angle of the ray, the nature of the joint and the material properties of both subsystems. Considering the relationship between the incident angle θiaN and the reflection angle θraN of the Nth intersection (as shown in Fig. 1), the transmission coefficient (transmission efficiency) for two perpendicular plates (commonly referred to as a corner joint) is 5 : where  and  are introduced to simplify the expression and m1′ and m2′ are the surface densities for each plate, respectively. When the two plates share the same properties,  =  = 1, thus: The incident energy of the Nth intersection is then calculated by: where E0 is the original total energy contained in each ray. Here, it is assumed that the damped subsystem is connected to four similar subsystems and they share similar properties except the internal loss factor. The averaged transmission coefficient at the boundary can then be determined using: where IEi is the incident energy at that specific boundary and i is the angle dependent transmission coefficient of that specific ray at the boundary. The coupling loss factors (CLFs) between subsystem i and j can then be calculated through 7 :

NUMERICAL VALIDATION OF THE RAY TRACING ALGORITHM
The model shown in Figure 1 is assumed to have five identical plates, sharing same material and geometric properties ( Table  1). The plates are connected perpendicular to each other. One hundred evenly distributed point sources are used to simulate a line source on the edge between plate 1 and plate 2. Each point source is assumed to radiate 180 acoustic rays. The travelling history of each ray is traced up to 50 reflections. If one assumes that there is no internal loss in plate 2, by using the ray tracing algorithm, the CLFs at 500 Hz can be calculated using Eqn. (16) and Eqn. (17) to be: The corresponding CLFs calculated by SEA are: Those two approaches yield very similar answers. Similar comparisons have been undertaken for 1/3 octave bands with center frequencies from 20 Hz to 20,000 Hz and have shown same findings, i.e., the difference between those two approaches is not frequency dependent. This comparison serves as a theoretical validation of the proposed ray tracing code. When the ILF of plate 2 gets higher, classic SEA tends to fail while this hybrid SEA still works.

SIMPLE ONE-ROOM SYSTEM
Assuming that there is a room with a dimension of 3m × 4m × 2.5 m. All structural elements are assumed to be made from concrete and share the same thickness and material properties, i.e., h = 0.15 m,  = 2300 kg/m 3 , cL = 3686 m/s and ILF = 0.015. The reverberation time for the room is assumed to be 0.5 s. The energy level difference between subsystem 1 and subsystem 6 can be predicted using SEA if one assumes that only subsystem 1 has a power input. A bending only model for SEA modeling is used for a direct comparison with the hybrid approach because the ray trace code only considers bending transmission. In addition, comparison between bending-only and three-wave models is conducted to prove that for the one-room system in this study, in-plane waves do not have significant contributions in the transmission in the frequency range considered.

Inherent (Case 1)
The system shown in Fig. 2 may be described using the SEA model and the classic SEA matrix can be written for the system as: where = 42 = 53 = 24 = 35 = 16 = 0. All the loss factors can be readily calculated using formulas suggested by Craik 5 .
For given power input W1, one may predict the energy level difference between subsystem 1 and subsystem 6 as:

Damping on Subsystem 2 (Case 2)
If one applies damping treatment on subsystem 2, i.e., ILF = 0.3, subsystem 2 is no longer modal and needs to be taken out of the SEA model (here it will be referred to as "coupling element 2" to distinguish it from a subsystem capable of supporting modes). However, energy transmissions through paths 1-2-3, 1-2-5, and 1-2-6 still exist (Fig. 3). One may use ray tracing code to determine the equivalent CLFs (between subsystem 1 and 3 via 2, 1 and 5 via 2, and 1 and 6 via 2) and substitute them into SEA matrix to work out the energy distribution. In the diagram of the systems under investigation, normal SEA subsystems are represented using solid rectangles with numbers. When a subsystem is heavily damped, it is treated as a coupling element, represented by circle with a number. The coupling loss factor between normal SEA subsystems for path 1-2-3 is 12 and 23, represented with solid arrow (Fig. 4a). When component 2 is heavily damped, it can no longer be represented using a subsystem, and the coupling between subsystem 1 and 3 is denoted as 123, which means the equivalent CLF from subsystem 1 to 3 via coupling element 2. This equivalent CLF is represented with dashed arrow (Fig. 4b). If there are two coupling elements in the path, the equivalent CLF is denoted as 1234, which means that the transmission occurs from subsystems 1 to 4 via coupling element 2 and 3 (Fig. 4c). These notations are intended for use in this article only and seek to aid the explanation of how heavily damped elements are accumulated within the SEA framework.
In the model presented in Fig. 2, for path 1-2-i, the energy transmitted from coupling element 2 to subsystem i has the relationship: where 12i is the equivalent CLF from subsystem 1 to sub system i via coupling element 2. The energy transmitted from subsystem 1 to coupling element 2 can be written as: Eqn. (20) and Eqn. (21) give: The ratio in the denominator can be readily determined from ray tracing code. One may note that W12/W2i is always greater than unity, which means that the equivalent CLF is always smaller than the CLF between the source subsystem and the heavily damped subsystem. Therefore, the new CLF from subsystem 1 to subsystem i can be expressed as: This means that the CLF from the source subsystem to the other subsystem contains two parts: one is the direct coupling from the source to the receiver, and the other considers the transmission via the coupling element. The general form of the new structural CLF from subsystem i to subsystem j, where it is assumed that the coupling element is 2, is thus: The SEA matrix when subsystem 2 is heavily damped is now:

Damping on Subsystems 2 and 4 (Case 3)
When both component 2 and component 4 are damped, i.e., ILF = 0.3, they are both considered as coupling elements. Similar analysis can be completed to give the new CLFs. The new structural coupling between subsystem i to subsystem j, where it is assumed that the coupling elements are 2 and 4, is: Therefore, the SEA matrix when both components 2 and 4 are heavily damped is:

Damping on Subsystems 2, 3, 4 and 5 (Case 4)
When components 2, 3, 4 and 5 are all damped, i.e., ILF = 0.3, they are all considered as coupling elements. In this case, the SEA model only has three subsystems, i.e., subsystems 1, 6 and 7. The couplings between the structural subsystems and the room remain unaltered. The determination of new CLF from subsystem 1 to subsystem 6 is complicated. The equivalent CLFs for paths such as path 1-2-6 can be calculated easily. However, for paths such as path 1-2-5-6, the equivalent CLFs are hard to determine. The ray trace code can predict the transmission and equivalent CLF for path 1-2-5, but the estimation of the equivalent CLF for path 2-5-6 is not that reliable. This is because, when component 2 is no longer modal, the energy trans-mission from component 2 to 5 is less likely to be a line source and hence the approximation using discrete acoustic rays with same initial energy tends to fail. This happens every time two coupling elements are physically connected.

RESULTS AND DISCUSSION
Predictions for the energy level difference between subsystems 1 and 6 for different damping arrangements are shown in Fig. 5. Higher energy level difference can generally be observed when adding damping treatment to the system. The results from SEA and hybrid ray-tracing-SEA predictions have similar results for case 2 where only component 2 is damped and case 3 where both components 2 and 4 are damped. When components 2, 3, 4 and 5 are all damped (case 4), the hybrid method tends to give higher level difference, especially at high frequencies. One may also find that firstorder flanking transmission in case 4 makes a big contribution in the transmission. Considering direct paths only results in overprediction of the energy level difference.
In order to further understand the reason why SEA and the hybrid approach give similar results in case 2 and case 3, the contribution of specific paths has been analyzed. SEA underpredicts the transmission across the perpendicular edges and overpredicts the transmission across the parallel edge. These errors may cancel with each other in a oneroom system. If one assumes that subsystem 1 is the source, the energy in subsystem n due to transmission along path 1-2-3-4….n can be given as 5 In case 2, the same analysis can be performed for paths 1-2-3, 1-2-4, 1-2-6 and path 1-2-1. In pure SEA modeling, one has: The total energy contribution from subsystem 1 via subsystem 2 to other structural subsystems is thus: If one excludes subsystem 2 from the SEA model and uses equivalent CLF, one has: Similarly, the total energy contribution is: The difference between EHybrid and ESEA is thus: The same analysis can be performed for paths 3-2-1, 3-2-6, 3-2-5 and 3-2-3. The result is shown in Fig. 6. If the energy enters subsystem 2 through the coupling between subsystems 1 and 2, SEA underpredicts the overall energy that leaves subsystem 2. If the energy enters subsystem 2 through the coupling between subsystems 3 and 2, SEA overpredicts the total energy that leaves subsystem 2. Those two errors have similar values and tend to cancel with each other. As a result, SEA gives similar prediction result as the hybrid method. Same curve can be obtained for case 3 because of geometry similarity.
The result for case 4 is shown in Fig. 7. SEA tends to predict higher transmission from subsystem 1 to subsystem 6. In this case, overprediction of the transmission via the parallel edge dominates the influence of underestimating the transmission via the perpendicular edges and results in an overprediction of the total transmission. This is the reason why SEA shows a smaller energy level difference to the hybrid method in Fig. 5 and the divergence between two modeling methods tends to increase at higher frequencies. To further understand the prediction difference when the coupling elements are physically connected, the same analysis is done assuming the case where components 2 and 3 are damped, and another case where components 2, 3 and 4 are damped. Again, only direct paths and first-order flanking paths are considered. The result is shown in Fig. 8. One may find that the difference between SEA and the hybrid prediction is small when only some of the flanking structures are heavily damped. The difference becomes more significant if one increases the frequency or the number of damped components. Significant difference occurs when all of the flanking structures (i.e., subsystems 2, 3, 4 and 5) are heavily damped, i.e., all the paths except the airborne path 1-7-6 have been damped. This finding suggests that when all of the structural paths are damped, SEA tends to overpredict the transmission. Hence, a greater benefit of applying damping treatment can be expected.
A plot of the energy level difference between subsystems 1 and 3 when subsystem 5 is damped and undamped due to specific path 1-5-6 is shown in Fig. 9 to help explain the result in Fig. 8. There exists significant difference between SEA and hybrid prediction for specific path when subsystem 5 is damped. When not all of the structural paths are damped, the transmission between subsystems 1 and 6 is dominant by the undamped paths so the difference between SEA and ray trace prediction is not obvious. This is another reason why SEA and hybrid methods show no significant difference in cases 1 and 2 as shown in Fig. 5. When all of the structural paths are damped, the estimation difference becomes much more obvious.

CONCLUSIONS
A definition of heavily damped subsystem is proposed as the subsystem with effective length ratio smaller than 0.5. Implication of ray trace in SEA framework has been done for a one-room system. The result shows that the prediction from pure SEA for damped system is relatively reliable when a small number of the structural components are heavily damped because the prediction errors tend to cancel with each other and the transmission is dominated by paths that are not damped. With the increasing number of paths being damped, SEA tends to underestimate the energy level difference, which implies that a greater benefit may be observed when applying damping treatments globally. The above findings are valid for simple one-room system and are believed to reflect the system behavior for more complicated systems since the underlying principle does not change. Further validation of these findings for more complicated systems is far more complicated and would require further numerical and experimental investigation. This study only considers specular reflection within homogenous plates, further inclusion of refraction, curvaton and sandwich-type plates, which may increase both accuracy and applicable range of the prediction results.