Investigation of radiation damping in sandwich structures using finite and boundary element methods and a nonlinear eigensolver

The exposure of human beings to vibration and noise can have implications ranging from annoyance to health damage. Hence, researchers of various fields, such as material scientists and control engineers, are concerned with the development of passive and active damping devices as well as the exploitation of material-inherent damping. This is particularly important for lightweight structures, which—generally speaking—are either stiff and weakly damped or exhibit high damping but rather poor elastic properties.^1–3 

However, an often neglected contribution to the overall damping of structures is the dissipation of vibrational energy due to sound radiation. While acoustic radiation damping is a rather insignificant aspect in many bulky engineering applications, it is the primary energy dissipating mechanism for stiff lightweight structures with large radiating surfaces. It follows that attempts to reduce the vibrational response of these lightweight structures by additional mechanical damping can only be successful if the extent of mechanical damping is comparable or larger than the extent of radiation damping.⁴ Therefore, engineers are in need of reliable and flexible methods for the quantification of radiation damping in an early stage of the design process.

However, due to the coupled nature of the problem, involving the behaviors of both structure and surrounding fluid, radiation damping is not generally amenable to analytical quantifications. Early theoretical methods predict the modal radiation damping of rectangular plates^5,6 and cylindrical shells.⁷ Expressions for frequency averaged radiation damping are also derived, assuming that a sufficiently large number of modes contributes to the vibration of the plate.⁸ These methods are all based on theoretical expressions of the radiation resistance^9–11 or theoretical expressions of the acoustic impedance of the plates.¹² They are only valid for homogeneous plates that are confined in an acoustically rigid baffle prohibiting flow between the two sides of the plate. Later, correction factors are proposed to account for unbaffled plates¹³ with arbitrary boundary conditions.¹⁴ However, their applicability to more complex geometric and material configurations can hardly be judged.

Sandwich structures, consisting of two thin and stiff face sheets enclosing a thick, lightweight and often anisotropic core, account for such complex configurations. While sandwich structures excel at the ratio of bending stiffness to mass, they exhibit relatively high flexural wave speeds compared to those of solid plates with equivalent mechanical properties. In consequence, coincidence between bending and acoustic waves occurs at relatively low frequencies. Moreover, due to the anisotropy, sandwich panels do not only exhibit a single critical frequency, but rather a range of frequencies in which coincidence occurs, thus giving rise to efficient sound radiation and hence high acoustic radiation damping in a wide frequency range.

Analytical expressions for the flexural vibration of sandwich panels can be derived from Hamilton's principle.¹⁵ Experimental^16,17 and numerical approaches^18–20 have been followed to investigate the vibroacoustic behavior of sandwich panels with respect to different core and face materials, lay-ups, and geometric configurations—see also the review by D'Alessandro et al.²¹ Besides vibroacoustic studies, many researchers have also made efforts to quantify material-inherent damping of sandwich panels to enhance damping by means of viscoelastic treatments.² Most of the experimental studies are conducted in air and hence the thereby obtained loss factors include the effects of acoustic radiation damping.

Clarkson and Brown deduced the radiation loss factors of a honeycomb sandwich platform by means of reference measurements inside a vacuum chamber.²² Zhou and Crocker determined radiation damping of sandwich plates clamped between two reverberation chambers based on principles of energy flow.¹⁹ Apart from these two articles, however, little published data on actual values for radiation damping of sandwich structures exist, although radiation damping can account for the major share in the overall damping and therefore undermine the effectiveness of additional mechanical damping.

In this paper, we employ a numerical framework based on the finite element method (FEM) and the boundary element method (BEM) in order to better understand the phenomenon of acoustic radiation damping. The structural and acoustic responses are fully coupled to enable the modeling of a mutual structural acoustic interaction as it occurs in many sandwich structures. The cores are represented by three-dimensional solid finite elements in order to capture local bending deformations of the individual face sheets that cause sound radiation in addition to the global bending deformations. Using this framework, we contribute in the following aspects to gain a deeper insight into the phenomenon of acoustic radiation damping:

• First, we study the extent of radiation damping for three sandwich panels subject to different boundary conditions in both acoustic full- and halfspaces. The harmonic radiation loss factors are obtained by relating the radiated sound power to the vibrational energy of the structure. The panels are excited by point forces as well as diffuse acoustic fields. The results indicate a strong influence of boundary conditions and excitations in the low frequency range, where the responses are mainly determined by modal behavior.

• Second, using a nonlinear eigensolver based on contour integration, we perform modal analyses of the air-loaded sandwich panels to deduce their modal radiation loss factors and eigenfrequencies. The latter are lowered compared to the in vacuo eigenfrequencies due to the effect of added mass and damping. The modal radiation loss factors, which are inherent properties of the structural acoustic system, agree well with the harmonic loss factors at the respective eigenfrequencies. Furthermore, we propose a more effective strategy for checking and filtering the eigenvalues when using contour integration and also provide guidance in choosing the solver-specific parameters.

• Last, we compare our numerically obtained estimates for radiation damping to theoretical expressions and experimental results found in the literature, generally yielding a good agreement. However, the comparison also indicates the limited applicability of commonly used theoretical expressions when coincidence occurs in a frequency range where the modes are still well separated. Finally, we discuss experimental quantification of radiation damping and associated sources of error such as the reinjection of acoustic energy and the reliability of reverberation room measurements in the low frequency range.

We consider the fully coupled structural acoustic interaction in order to determine the vibratory response of sandwich structures. Under the assumption of a harmonic time dependency $e−iωt$⁠, the equations of linear elasticity and acoustics are discretized using FEM²³ and direct collocation BEM.²⁴ The resulting systems of equations read

and

Therein, u and p are the vectors of unknown displacement and sound pressure values at the nodes, respectively. The stiffness and mass matrices of the structure are denoted with K and M, respectively. The structure is excited by external forces $fs$ as well as fluid forces $ff$⁠. The latter act by virtue of the acoustic field. Structure-inherent damping is not considered in this work and, hence, acoustic radiation damping is the only dissipative mechanism occurring. Further, H and G are the frequency dependent boundary element (BE) matrices, relating the structural particle velocity $vf$ to the sound pressure. Acoustic excitation is taken into account by the incident sound pressure field $pi$ and the corresponding incident particle velocity $vfi$⁠. The angular frequency is defined as $ω=2πf$⁠, and $i$ denotes the imaginary unit.

Since we are particularly interested in applications that exhibit considerable levels of radiation damping, the influence of the acoustic field on the structural response is not generally negligible. Consequently, it is not sufficient to determine the in vacuo response of the structure by solving Eq. (1), and subsequently evaluate the acoustic field using Eq. (2) in a post-processing step. Instead, Eqs. (1) and (2) are mutually coupled on the sound radiating surface, i.e.,

The mesh coupling is established by the coupling matrices $Csf$ and $Cfs$⁠, relating the displacement and pressure degrees of freedom (DOFs).²⁵ Since structural acoustic interaction is mainly relevant for thin-walled lightweight structures, most researchers rely on shell finite elements for modeling the structural subdomain.^18,20,25 While the sandwich panels considered in this work can certainly be modeled using layered shell formulations as well, nevertheless, we will follow a different approach involving three-dimensional solid finite elements for the representation of the thick core. Additionally, shell elements are employed for the thin face sheets. While this approach leads to more DOFs compared to the use of layered shell elements, it enables us to capture local bending deformations of the individual face sheets, which otherwise would not be possible. These local bending deformations of the face sheets—also known as symmetric motion—involve thickness deformations of the core. They cause sound radiation in addition to the global bending deformations (anti-symmetric motion). These two types of lamb waves,²⁶ which are shown in Fig. 1, coexist in sandwich panels.

Regarding the boundary conditions of the panel, both the freely suspended and the simply supported cases are considered in this work. The simply supported conditions are modeled by constraining the displacement DOFs of the face sheet edges, as is schematically shown in Fig. 2.

The approach for modeling the vibroacoustic behavior of a sandwich panel using finite and boundary elements is schematically depicted in Fig. 3. The shell finite elements that represent the face sheets are defined with an offset of half the shell thickness. In this way, their nodes coincide with the outer nodes of the solid elements representing the core. These nodes on the top and bottom surfaces are also the ones that are coupled to the nodes of the boundary element mesh. In the case of an unbaffled panel—i.e., a panel with free edges where acoustic short-circuiting occurs—a single boundary element mesh with a closed surface is used. Otherwise, when the panel is confined in an acoustically rigid baffle, the two independent acoustic subdomains on each side of the panel are modeled using a halfspace formulation with a modified Green's function.²⁷ In this way, the dissipation of vibrational energy due to sound radiation is considered simultaneously on both sides of the baffled panel. The boundary element meshes corresponding to a baffled panel are shown in Fig. 4. Note that it is also possible to only model half of the baffled panel along with a single acoustic halfspace. However, this approach would require separate computations and subsequent superpositions of the symmetric and antisymmetric responses (cf. Fig. 1). While this would result in fewer DOFs, nevertheless, we use a full model involving two acoustic halfspaces for the sake of convenience.

In the case of a single acoustic subdomain, the global system of equations containing the coupling conditions emerges as

If the panel is confined in a baffle, the global system comprises three subdomains. Assuming that just one side of the panel is excited by an incident sound field (which actually resembles the situation in a window test rig), the resulting monolithic equation is given as

where $(·)(I)$ and $(·)(II)$ denote the acoustic halfspaces on the respective sides of the panel. In the case of geometrical symmetry with respect to the plane of the baffle, $H(I)=H(II)$ and $G(I)=G(II)$ hold, and consequently, the numerical integration for assembling the BE matrices needs to be performed only once.

The complex sound power P in linear time-harmonic acoustics can be obtained from

where $vf$ denotes the fluid particle velocity, and $(·)*$ is the conjugate complex. In the discrete setting, the sound power is evaluated as a post-processing step. The nodal values for the sound pressure are related to the particle velocity via Eq. (2), and the integration of their interpolation functions results in the boundary mass matrix $Θ$⁠. Finally, substitution by the acoustic impedance matrix $Z=(H−1G)TΘ$ yields the complex sound power in the discrete setting

Only the real part $Re(·)$ of the above expression contributes to the radiation to the far-field and hence to the structural damping due to sound radiation. The latter is quantified by relating the radiated sound power to the power corresponding to the total energy of the vibrating structure.²⁸ The time-averaged total vibrational energy equals twice the time-averaged kinetic energy, or equivalently, twice the time-averaged potential energy. For harmonic problems, these energy quantities can be determined from the structural response via²⁹ 

and

where the first term in Eq. (9) corresponds to the energy due to the elastic strain, and the second term is the work done by external forces. The evaluation of the kinetic energy $Ek$ requires knowledge of the additional mass $Mf$ due to acoustic loading. This frequency dependent mass contribution could be approximated by the second order term of a Taylor expansion of the acoustic impedance matrix Z.³⁰ However, for our purposes it is more convenient to simply use the potential energy $Ep$ to quantify the radiation loss factor. Hence, the radiation loss factor is expressed by²⁸ 

Note that the kinetic energy $Ek$ could be equally used to evaluate the radiation loss factor. Recent results^31,32 show that spurious numerical damping could lead to an overestimation of damping phenomena when studying them with BEM. However, the occurrence of numerical damping does not seem to be an issue in exterior acoustics.

Modal analyses provide useful information on the properties of the system, such as the eigenfrequencies of the fluid-loaded structure. In this work, in particular, it serves as an alternative way to quantify the extent of radiation damping. The modal radiation loss factors can be deduced from the complex eigenvalues of the structural acoustic system. At resonance, these modal loss factors are expected to agree with the continuous radiation loss factor defined in Eq. (10).

The purely structural equation subject to acoustic loading is obtained by forming the Schur complement of Eq. (4) and thereby omitting the pressure DOFs,³³ i.e.,

in which $iωCsfH−1GCfs$ can be interpreted as the effect of fluid loading. Note that the Schur complement of Eq. (5) can be obtained in a similar manner. By setting the right-hand side to zero, we arrive at the definition of the structural acoustic eigenvalue problem (EVP)

with the fluid-loaded structural mode v and the complex eigenfrequency $ω̃$⁠. The EVP in Eq. (12) is nonlinear since the BE matrices H and G implicitly depend on the frequency.

Several methods have been proposed for the solution of Eq. (12) during the last years. Peters et al.³⁰ employed a truncated Taylor series to approximate the frequency dependent matrices, and the resultant polynomial EVP is solved using symmetric linearization. In a subsequent work, the computational effort associated with the linearized EVP is addressed by means of Krylov subspace model order reduction of the structural subproblem.³⁴ However, the success of this method strongly depends on the convergence radius of the Taylor approximation and the decay of the coefficients of the polynomial approximation. As a remedy, the frequency range of interest needs to be subdivided, whereas proper choices of these sub-frequency ranges can hardly be made a priori.

Therefore, in recent years alternative approaches for the solutions of nonlinear EVPs have been proposed, which can be classified as contour integral methods.^35–38 Using contour integration, a nonlinear EVP is converted to a generalized EVP of reduced dimension that exhibits identical eigenvalues inside a predefined region in the complex plane. Contour integral methods are particularly appealing because of their general applicability and suitability for the execution on distributed parallel computers.

While we assume that the other contour integral methods would also fulfill our purpose of investigating air-loaded modes and radiation damping of sandwich panels, we choose to use the block Sakurai Sugiura method (block SS)^36,39 in this work. A comparison of different eigenvalue solvers is beyond the scope of this work. Moreover, we note that the focus of our contribution is not the further development of existing methods but rather its application in the context of air-loaded elastic structures. Most of following content on block SS can be also found in the papers by Asakura et al.³⁶ and Zheng et al.^36,39 Since we nevertheless propose a more effective strategy for checking and filtering the eigenvalues—which is crucial when using contour integral methods—the procedure is briefly outlined in what follows.

Block SS is a direct method, and it essentially works by replacing the nonlinear EVP in Eq. (12) by the generalized EVP

with the eigenpair $(ψ,λ)$⁠. The block Hankel matrices $H1,H2∈ℂKL×KL$ are defined as

where K and L are user-specified positive integers. The proper choice of K and L will be discussed in detail later on. The moments $Ml∈ℂL×L$ are computed from

where U and V contain randomly chosen source vectors as columns, and $(⋅)H$ denotes the Hermitian transpose. The original system matrix B is evaluated at the complex frequency parameter σ. The latter is defined along C—a closed non-self-intersecting continuous loop in the complex plane. Once the reduced EVP in Eq. (13) is solved, the fluid-loaded structural mode can be recovered from

The corresponding eigenvalue λ equals the complex eigenfrequency $ω̃$ of the original system in Eq. (12). The block matrix $S=[S0,…,SK−1]$ is also obtained by contour integration via

The range of obtained eigenvalues is enclosed by the contour C along which the integrals in Eqs. (15) and (17) are evaluated. This contour needs to be predefined by the user. In the context of fluid-loaded structures, a suitable choice of the contour is an ellipse that has its major axis aligned with the real axis. The two vertices on the real axis correspond to the upper and lower limits (⁠$fmax, fmin$⁠) of the frequency range of interest. A suitable ellipse is shown in Fig. 5 and can be expressed by

where $γ=(fmax+fmin)/2$ and $ρ=(fmax−fmin)/2$⁠. The factor ζ defines the shape of the ellipse and should be chosen according to the expected ratio of imaginary and real parts of the eigenvalues. Generally, the ellipse should be wide and short (⁠$ζ<1$⁠)—especially in the case of weak to moderately strong structural acoustic coupling. With the definition of the contour at hand, the integrals in Eqs. (15) and (17) are approximated using the N-point trapezoidal rule, i.e.,

where N denotes the number of integration points on the contour, and $θj=2π(j−1)/N, j=1,…,N.$ Using the approximated moments $M̂l$⁠, the (approximated) Hankel matrices $Ĥ1$ and $Ĥ2$ can be assembled according to Eq. (14). Finally, the corresponding generalized EVP $Ĥ1ψ̂j=λ̂jĤ2ψ̂j$ is solved and the complex eigenfrequencies $ω̃j$⁠, as well as the fluid-loaded modes $vj, j=1,…,KL$⁠, are recovered from

The main challenge with the use of block SS lies in the choice of the following parameters: the degree of moments K, the number of source vectors L, and the number of integration points N. These parameters are related to the computational effort of the method, as well as to the completeness of the determined eigenvalues, i.e., whether all eigenvalues lying inside the contour are found.

Given a fixed number of integration points N, Sakurai et al.⁴⁰ suggest to set the degree of moments to $K=N/4$ as a good compromise between accuracy and numerical efficiency of the algorithm. With L = 1, we have the original SS method,³⁵ while the choice L > 1 results in the block SS method,³⁶ which achieves higher accuracy at a similar numerical cost.⁴¹ In this work, we found that the eigensolutions improved with an increasing number of L up to roughly L = 10. Higher values than that did not change the results anymore as long as the product KL was large enough. This product defines the dimension of the subspace and hence corresponds to the number of eigenvalues that are obtained from the reduced system. Therefore, KL should be at least as large as the expected number of eigenvalues inside the contour. This number can be estimated a priori by an empirical formula.⁴⁰ However, in our case of sandwich structures interacting with air, we will rather rely on the knowledge of the number of in vacuo eigenfrequencies inside the given frequency range $[fmin, fmax]$⁠. More specifically, we solve the purely structural EVP

with the in vacuo modes $vdry$ and eigenfrequencies $ωdry$ prior to the solution of the structural acoustic EVP in Eq. (12). Then, assuming that in a given frequency range, the number of structural acoustic eigenvalues roughly equals the number $ndry$ of eigenfrequencies, we set the dimension of the reduced EVP such that $KL>ndry$⁠.

The number of integration points N determines the number of linear systems of equations $B(σ)X=V$ that have to be solved for evaluating Eq. (19), accounting for the main computational effort. In cases where the algorithm is executed in a parallel computing environment, N is chosen according to the available computing nodes. Sakurai et al.⁴⁰ note that a large N is not necessary for an accurate quadrature and suggest, e.g., N = 16 or 32. Whereas the results of Zheng et al.³⁹ confirm this suggestion, they also show that iteratively increasing N is a suitable way for checking whether all eigenvalues inside the contour are found and also for distinguishing them from spurious eigenvalues. The latter mainly occur due to the projection.

However, in the context of weak to moderately strong structural acoustic interaction, a more effective strategy for checking and filtering the eigenvalues is available based on the modal assurance criterion (MAC).⁴² Assuming that the modes of a structure subject to light fluid-loading are similar to the in vacuo modes, the validity of a complex eigenvalue $ω̃j$ can be simply tested by checking the occurrence of its associated mode $vj$ in the range of in vacuo modes, i.e.,

The values of $MAC$ range from 0, indicating no correspondence between the two modes, to 1, representing a consistent correspondence. In the case of air-loaded sandwich panels, we can expect to find an in vacuo mode $vdry$ for each actual fluid-loaded mode $vj$ that satisfies $MAC=0.95 … 1.00$⁠. Otherwise, $vj$ and the corresponding eigenvalue $ω̃j$ are identified as spurious. On the other hand, if we are interested in the complex eigenfrequencies corresponding to particular modes, a criterion based on Eq. (23) is useful to check if that eigenfrequency has been found by block SS. If not, the accuracy of the projection needs to be improved by increasing the number of integration points N—of course given that the underlying elliptic contour actually encloses that eigenvalue. Once all desired modes are found, the accuracy of eigenpairs $(ω̃,v)$ can be assessed by inserting them into Eq. (12) and computing the backward error

Furthermore, the accuracy of eigenpairs can be subsequently improved by repeating the contour integration (19) with additional integration points placed in between the previous ones. In this way, the additional numerical effort is limited to the computations required for the new integration points, while the intermediate solutions corresponding to the previous integration points can be reused. Thus, a strategy with a gradually increasing number of N is only marginally more expensive than a single execution of the procedure with the final (i.e., largest) number of N.

In the following, we will study the vibroacoustic behavior of three honeycomb sandwich panels in air with particular focus on acoustic radiation damping. The results obtained using the presented numerical framework are then compared to theoretical expressions as well as to experimental results available in the literature.

Panel A consists of two plywood face sheets enclosing a paper honeycomb core. The vibroacoustic behavior of a similar panel is experimentally investigated in the pioneering work by Moore,¹⁶ which has also served as a benchmark for other researchers in the past.^21,43 Panels B and C are made of plane weave fabric-reinforced graphite composite face sheets and a polyurethane foam-filled honeycomb core.¹⁹ The dimensions and the material properties of the panels A, B, and C are presented in Table I. The freely suspended boundary condition is denoted with “-free” and the simply supported case with “-SS.” For example, panel A subject to simply supported boundary conditions will be referred to as panel A-SS in what follows. Unless otherwise stated, the freely suspended panels are excited by a point force of $Fz=1 N$ at the corner node (⁠$x=y=0, z=h/2$⁠), and the simply supported panels likewise at the center node (⁠$x=lx/2, y=ly/2, z=h/2$⁠). Moreover, excitation by diffuse acoustic fields will also be considered in Sec. III D. The definition of the coordinate system is depicted in Fig. 6.

Eight-noded quadrilateral shell finite elements based on the Reissner-Mindlin theory are employed for the modeling of the face sheets. The cores of the panels are modeled using 20-noded hexahedral solid finite elements. The respective stiffness and mass matrices are extracted from ANSYS.⁴⁴ For all three panels, a uniform finite element mesh with 48 and 96 elements along the in-plane directions and 2 elements in the thickness direction is used (recall that all 3 panels have an aspect ratio of roughly 2:1). The finite element (FE) meshes result in 240 000 displacement DOFs. This corresponds to three quadratic elements per bending wave length of panel A in the frequency range up to $2000 Hz$⁠. For panels B and C, this mesh results in at least 6 and 13 elements per bending wave length, respectively, in the considered frequency ranges.

Quadrilateral boundary elements with bilinear interpolation functions are used for the discretization of the surrounding acoustic field. These elements have their DOF-carrying nodes inside the element, rather than on the element edges, resulting in a sound pressure interpolation that is discontinuous across element boundaries.⁴⁵ The corresponding BE matrices H and G are extracted from the non-commercial software AKUSTA.⁴⁵ A treatment for the nonuniqueness problem that occurs in exterior BE formulations is not required since the respective first irregular frequencies of the panels are beyond the frequency range of interest, e.g., the first spurious mode of panel A occurs at approximately $27 kHz$⁠.

Regarding the size of the boundary element mesh, there are no guidelines available in the literature for coupled structural acoustic radiation problems. Therefore, in this work, the adequate mesh size for the acoustic field is chosen based on a convergence study. Figure 7 shows the relative difference in radiated sound power of panel A-SS for different BE mesh sizes. These meshes do not conform with the above defined structural finite element mesh. Consequently, the relative difference shown in Fig. 7 is related to the discretization error of the acoustic field, as well as to the error introduced by the mesh coupling. It is calculated from $ϵrel,P=|P−Pref|/Pref$⁠, where $Pref$ denotes the reference sound power for a mesh with 48 and 96 elements along the in-plane directions corresponding to 19 584 pressure DOFs for each acoustic halfspace. This reference BE mesh conforms with the above defined FE mesh. As we expected, the relative difference displayed in Fig. 7 decreases monotonically as we refine the mesh. Moreover, the relative difference is of the same order of magnitude at four different frequency points, although some of those frequencies lie well above the coincidence frequency range of the panel (see Sec. III E for discussion on coincidence).

Based on this convergence study, we choose the BE mesh with 24 and 48 elements along the in-plane directions for all upcoming simulations. This mesh has 5148 pressure DOFs for each halfspace and, similarly, 10 296 DOFs for the acoustic fullspace corresponding to the unbaffled panels. Compared to the reference mesh that has four times more DOFs, it results in a relative difference of less than 3.3% (0.14 dB) at all considered frequency points. Finally, in order to assess the influence of the discretization error on the eigenfrequencies, the modal analysis that is presented in Sec. III B was repeated with the reference mesh, resulting in maximum relative differences of 0.17% in the imaginary and 0.07% in the real part of the eigenvalues.

The modal analysis scheme presented in Sec. II is now applied to panel A to obtain its air-loaded modes and associated eigenfrequecies. From a preceding in vacuo analysis, we expect panel A-Free to have 15 eigenfrequencies in the frequency range $[fmin=10 Hz,fmax=500 Hz]$⁠. Trivially, panel A-Free also exhibits six rigid body modes which are not affected by the acoustic loading. The in vacuo analysis of panel A-SS yields 11 eigenfrequencies in the frequency range $[fmin=100 Hz,fmax=600 Hz]$⁠. These bounds are also chosen for the definition of the respective ellipses in Eq. (18). The number of integration points N was gradually increased until the accuracy of the eigenpairs stagnated, resulting in N = 32. Moreover, K = 8 moments, L = 15 source vectors, and an aspect ratio of $ζ=0.1$ were chosen for both panels. In the considered examples, we found that the eigensolutions are relatively insensitive to the choice of K and L as long as the resulting subspace was large enough and $L≥10$⁠. After solving the generalized EVP, the modes corresponding to the eigenvalues lying inside the contour are checked using MAC as given in Eq. (23).

The eigenfrequencies of the air-loaded panel correspond to the real part of the eigenvalue, i.e., $fj=Re(ω̃j)/(2π)$⁠. They are given in Tables II and III, along with the in vacuo eigenfrequencies $fdry,j$ and their relative difference $Δj=|fdry,j−fj|/fdry,j$⁠. As we expect, the eigenfrequencies of the air-loaded panels are generally lowered due to the effect of added mass and radiation damping. The actual extent of the frequency shift depends on the shape of the associated mode. For example, the eigenfrequency f₁ of panel A-SS associated with the fundamental bending mode is significantly lowered (⁠$Δ1=−2.39 %$⁠), while the eigenfrequency f₁₁ of panel A-Free that belongs to an in-plane mode is almost unaffected by the air-loading (⁠$Δ11=−0.03 %$⁠).

Regarding the numerical accuracy of the frequency shifts Δ_j, we distinguish between the discretization error and the accuracy of the eigensolver. As mentioned in Sec. III A, the discretization error of the wet eigenfrequencies of panel A-SS is at most 0.07%. Although this might seem sufficiently accurate at first sight, it nevertheless needs to be set in relation to the frequency shifts. In the presented examples, the discretization error is mostly 1 order of magnitude smaller than the computed frequency shifts Δ_j in both panels. The accuracy of the eigensolver is assessed by computing the relative residuals of the air-loaded eigenpairs by using Eq. (24). They are given in Tables II and III and show that errors of order $O(10−5)$ are achieved. This verifies the accuracy of the presented modal analysis scheme.

Expressing the sound radiation by means of modal contributions is a popular procedure to accelerate active control applications. For instance, a few (orthogonal) surface velocity patterns are usually sufficient to approximate the total radiated sound power at a certain frequency. These patterns, also known as acoustic radiation modes,⁴⁶ are the eigenvectors of the acoustic impedance matrix and computed using BEM or analytical methods.⁴⁷ 

Alternatively, we can express the radiated sound power in terms of the fluid-loaded modes of the structure. This requires orthonormalization⁴⁸ of the modal basis obtained from block SS such that

holds with the orthonormal modal matrix V containing the fluid-loaded modes $v̂j$ as columns. Due to orthonormalization with respect to the frequency dependent matrix $B(ω)$⁠, these modes are now also frequency dependent. Given a structural force excitation $fs$⁠, the structural displacement in Eq. (11) can be expressed by exploiting the condition (25), i.e.,

For an individual mode $v̂j$⁠, the modal particle displacement $dj$ at the acoustic nodes can be written as³⁰ 

Inserting Eq. (27) into Eq. (7) yields the complex modal sound power contributions

Figure 8 displays the real parts of the sound power resulting from point force excitations of panels A-SS and A-Free. Here, panel A-SS is not excited at the center but at the node (⁠$x=0.61 m,y=1.754 m,z=0.038 m$⁠) to ensure that a larger number of modes participate in the response.

As expected, several resonances occur in the considered frequency ranges. In the case of panel A-SS, the peaks are noticeably rounded. Given that structure-inherent damping is not modeled here, this clearly indicates the effect of energy dissipation by sound radiation. While the first couple of resonances of A-Free exhibit sharp maxima, the effect of acoustic radiation damping also comes into play in the higher frequency range.

Besides the sound power that is obtained from a harmonic analysis, Fig. 8 also displays the diagonal modal contributions $Re(Pjj)$⁠, as well as the superposition of all (diagonal and off-diagonal) contributions, i.e., $∑j∑kRe(Pjk)$⁠. The total modal superposition agrees well with the harmonic sound power except in the higher frequency range, where the modal basis is obviously not sufficient. The relative differences between the harmonic sound power and the total modal superposition are given in Fig. 9.

At the resonances, the diagonal values almost exclusively contribute to the radiated sound power. The rigid body modes of panel A-Free are also included in the modal basis, and determine the sound radiation up to around 50 Hz. At around 100 Hz, they even exceed the total radiated sound power, which is an indication that off-diagonal sound power contributions with negative signs occur, i.e., $Re(Pjk)<0, j≠k$⁠.

The modal displacements $dj$ are generally not orthogonal with respect to the acoustic impedance matrix Z, and cross-coupling between two modes can occur. Therefore, the off-diagonal sound power contributions are not necessarily zero. In fact, the occurrence of these off-diagonal contributions can be interpreted in that the spatial distribution of inertial forces of the structure is different from the spatial distribution of the inertial forces due to the acoustic loading. We notice that despite the relatively weak structural acoustic interaction between the sandwich panels and the air, these off-diagonal values significantly contribute to the overall sound power radiation.

By relating the power loss due to far-field sound radiation to the vibrational energy of the structure, the radiation loss factor quantifies the extent of acoustic radiation damping. The harmonic radiation loss factor as given in Eq. (10) is the result of a frequency-wise response analysis and generally depends on the excitation.

Figure 10 displays the radiation loss factor for the panels A-SS and A-Free subject to point-force excitation. While both panels show qualitatively similar behaviors with an increase of radiation damping toward the coincidence region and a subsequent plateauing, significant differences in the magnitudes are observed in the low frequency range. Panel A-SS already exhibits considerable radiation damping in the low frequency range by virtue of the fundamental bending mode that exhibits a monopole radiation characteristic. In contrast, the effect of acoustic short-circuiting in conjunction with freely moving edges of panel A-Free lead to much lower radiation loss factors in the low frequency range. At the higher frequency range, when the panels contain a few bending waves, the effect of boundary conditions becomes insignificant, and both panels exhibit similar radiation loss factors. Values of $ηr>0.01$ across a wide frequency range indicate the relevance of acoustic radiation damping in honeycomb sandwich structures—particularly when considering that material-inherent loss factors are typically of the same order of magnitude.

Besides the boundary conditions, the excitation can have a significant influence on radiation damping at low frequencies as well. This is reflected in Fig. 11, which is a close-up of Fig. 10 in the low frequency range. Additionally, it shows averaged loss factors of the panels subject to 100 randomly located point forces as well as loss factors for diffuse field excitation. The responses to a diffuse incident field were computed as the mean values of 50 simulations, where the excitation in each simulation was given by the summation of 1145 random incident plane waves arriving from uniformly distributed directions in space. For a detailed description of this procedure, we refer to the appendix of a paper by Rafaely.⁴⁹ Regarding the extent of radiation damping, we notice that the diffuse field excitation leads to higher loss factors in the low frequency range in both panels. This can be explained by the spatially uniform distribution of the incident pressure fields in the low frequency range that almost act like a plane wave excitation. For panel A-Free around 125 Hz in particular, this leads to monopole radiation characteristic which is not achieved by point excitation.

In addition to the harmonic radiation loss factor, modal loss factors can be obtained characterizing the radiation damping of each individual fluid-loaded structural mode. These modal radiation loss factors are properties of the structural acoustic system and hence independent of the excitation. At the complex eigenfrequency $ω̃j$⁠, the modal radiation loss factor is defined as^4,50

in which $Im(ω̃j)$ is negative due to the choice of the time dependency $e−iωt$⁠. The modal loss factors are given in Fig. 11 for panels A-SS and A-Free at their respective eigenfrequencies. In the case of point excitation, the harmonic loss factors deviate from some of the modal loss factors, indicating that the respective modes are not (or not exclusively) excited in the harmonic analysis. This is particularly obvious for the in-plane mode of panel A-Free that occurs at $f11=389 Hz$ and exhibits only marginal radiation damping (⁠$η11=2×10−4$⁠). On the other hand, the harmonic loss factor corresponding to the diffuse field excitation coincides with all modal loss factors for both panels without exception. Summing up, both the modal and the harmonic radiation loss factors provide useful measures to characterize the extent of radiation damping.

So far, we have placed our attention only on numerical methods and how they can be employed to study the radiation damping of sandwich structures. In the following, we will compare our numerical results to commonly used theoretical expressions and experimental results available in the literature.

Theoretical expressions for radiation damping rely on approximations of the radiation resistance R_r,^8,11 i.e.,

where $mpanel$ denotes the mass of the panel. Several authors have derived the radiation resistance of simply supported, baffled plates based on the concept of power flow,^8,9 and correction factors have been proposed for taking the effect of acoustic short-circuiting into account.¹³ All of these expressions assume a multi-modal radiation of the panel and thus are not applicable in the low frequency range. A summary and discussion of these expressions for radiation resistance can be found in a publication by Renji and Nair¹¹ in which the authors also point out that some of the expressions in the above-mentioned literature have inconsistent factors. Here, we use the expression as given in Eq. (11) of the paper by Renji and Nair.¹¹ 

The accuracy of theoretical radiation resistance estimates, in turn, depends on the prediction of the critical frequency. For composite panels with symmetric cross-ply laminates, the critical frequency $fc$ under consideration of transverse shear effects can be estimated from⁵¹ 

where $ρs$ denotes the surface density of the panel, c is the speed of sound, and $α=(D12+2D66)/D$ with $D=D11D22$⁠. The flexural rigidity values D₁₁, D₂₂, D₁₂, D₆₆, and the shear rigidity N can be obtained from the properties of the face sheets and the core based on laminate theory.⁵² Using Eq. (31) yields critical frequency estimates of 146 Hz for panel A, 122 Hz for panel B, and 780 Hz for panel C.

The third-octave band-averaged theoretical estimates for the radiation loss factor resulting from Eq. (30) are given in Figs. 12 and 13, along with the numerically obtained results for panels A-SS, B-SS, and C-SS. The critical frequencies cannot be identified directly from the numerical results becayse this would require determining the bending wave content by means of a spatial Fourier transform. However, from Figs. 12 and 13 we can observe that the radiation loss factor exhibits a plateau in the higher frequency range. This plateauing that generally occurs above the critical frequency is in accordance with the theoretical results, while the actual radiation loss factors at the critical frequencies are significantly overestimated by the theoretical expressions. In the higher frequency range, where the radiation loss factors level out, the theoretical and numerical results of all three panels are in good agreement. This indicates that enough modes contribute to the radiation of the panel in this frequency range so that the theoretical expressions are valid. At lower frequencies, however, the radiation damping of panels B-SS and C-SS strongly depends on which modes are excited. Moreover, the modes are still widely separated while significant radiation damping already occurs. However, as already mentioned above, the theoretical expressions assume a multi-modal radiation of the panel. This leads to the conclusion that they are not sufficient to comprehensively assess the radiation damping of such sandwich structures, at least when coincidence and, thus, efficient sound radiation already occur at low frequencies.

Furthermore, experimentally determined loss factors of similar panels are taken from the literature, and they are also compared to our numerical results in what follows. Panel A-SS was tested my Moore¹⁶ in a window between two reverberation rooms. It was excited by a loudspeaker in the sending room, and the sound pressure, as well as the space averaged mean square accelerations of the panel, were measured in the receiving room. By relating the radiated sound power to the vibrational level, Moore¹⁶ obtained third-octave band-averaged radiation efficiencies of panel A-SS (see Fig. 5.8 in Ref. 16). The associated radiation loss factor can be obtained from the radiation efficiency $σ̂r$ via

where ρ denotes the density of air. Zhou and Crocker¹⁹ have conducted similar measurements to obtain the third-octave band-averaged radiation loss factors of the panels B-SS and C-SS. The panels were clamped in a window between two reverberation rooms and excited by a shaker. Since the loss factors resulting from the above-mentioned experiments are associated with the sound radiation of only one side of the panel, they are multiplied by a factor of 2 in order to compare them to our numerical results.

The experimentally determined radiation loss factors of panels A-SS, B-SS, and C-SS are displayed in Figs. 12 and 13, along with the numerical results. In general, the experimental and numerical results agree well for all three panels. Above the critical frequency, the experimental loss factors exhibit a similar leveling-off as the numerical results, although in this frequency range the theoretical estimates provide better agreements with the numerical results. Conversely, around the critical frequency, where the theoretical expressions significantly overestimate radiation damping, the experimental values provide better agreements with the numerical results than the theoretical estimates do. In the subcritical range, however, the experimentally obtained loss factors fall significantly below the numerical ones. Two explanations are possible for this discrepancy.

The first explanation is related to the boundary conditions of the panels. When testing panels in a window test rig, they are typically clamped between the two walls of adjacent rooms. This clamping, however, is far from ideally rigid and will, to a certain extent, always exhibit compliance. It is clear that the boundary conditions of a particular window test rig can hardly be reproduced in simulations, and as a compromise, simply supported boundary conditions were imposed on all face sheet edges, as shown in Fig. 2. The difference to the actual boundary conditions in the test rig will certainly have an influence on the low frequency results.

The second possible explanation for the discrepancy between the experimentally and numerically obtained loss factors is related to the reliability of reverberation room measurements in the low frequency range, in general. The low modal density of the room results in a nonuniform sound pressure field, and therefore microphone-based measurements are subject to high uncertainties. This issue could be addressed by a recently proposed experimental procedure in which mobility measurements are combined with a numerically obtained acoustic impedance matrix to compute the acoustic response.⁵³ In this way, experimental estimates of the radiation loss factor that only depend on the properties of the panel could be obtained.

In addition to the simply supported panels, the freely suspended, non-baffled panels B-Free and C-Free are also studied, and their respective numerical and experimental radiation loss factors are displayed in Fig. 14. Zhou and Crocker¹⁹ obtained the experimental values by exciting the freely hanging panels with a shaker and measuring both the sound pressure in the reverberation room and the mean vibrational velocity of the panel. While the numerical and experimental loss factors qualitatively show similar behaviors with an increase toward the coincidence region and a subsequent plateau, the actual magnitudes differ significantly in the higher frequency range. There, the numerically determined loss factors are higher throughout than the experimental ones.

This deviation could be related to the reinjection of energy due to reflections in the reverberation rooms. While the numerical models assume that all the radiated sound energy disappears in the far-field, in fact, part of the acoustic energy in a reverberant room is transferred back to the panel and, hence, serves as an excitation in addition to the mechanical excitation by the shaker. This line of reasoning becomes clear when considering how Zhou and Crocker¹⁹ obtained the radiation loss factor estimates from their experimental data. For this purpose, recall the power balance inside a reverberation room, i.e.,

where $Eroom$ and $Epanel$ denote the total mean energies of the acoustic field and the panel, respectively. The dissipation loss factor of the reverberation room $ηroom$ includes the sound power absorption of the walls due to air. The coupling loss factor $ηc$ reflects the transfer of acoustic energy from the room to the panel—similar to the radiation loss factor $ηr$ that quantifies the energy transfer from the panel to the acoustic field in the room. The energy quantities are given as

where $p¯$ and $v¯$ are the experimentally obtained averaged values for the sound pressure and the vibrational velocity. Further, $Vroom$ denotes the volume of the reverberation room. The coupling loss factor $ηc$ in Eq. (33) is defined based on considerations of statistical energy analysis (SEA) from the reciprocal relationship¹¹ 

with the modal density of the panel $npanel$ and the modal density of the room $nroom$⁠. Combining Eqs. (30), (33), (34), and (35), we arrive at an expression for the radiation resistance that reads

Indeed, Zhou and Crocker used another, simplified expression in order to deduce the radiation resistance from their experimental data. Their expression [cf. Eq. (21) in Ref. 19] reads

A comparison to Eq. (36) reveals that Eq. (37) misses the second term in the denominator and as a consequence, neglects the reinjection of acoustic energy into the panel. This could lead to an underestimation of the radiation resistance, which would explain the deviation between the numerical and experimental estimates for radiation damping that we observe in Fig. 14 in the higher frequency range. In future experiments, the second term in the denominator of Eq. (36) could be evaluated to assess its impact on the actual radiation damping values.

Finally, we note that the above-mentioned paper¹⁹ is not the only publication of experimental results where the reinjection of acoustic energy is left unconsidered. In fact, in the pioneering work²² on radiation damping of sandwich panels, Clarkson and Brown deduce radiation damping by means of reference measurements in a vacuum chamber. However, the measurements in air are conducted inside the (disabled) vacuum chamber as well, which clearly leads to reflection at the inner walls of the chamber, and therefore to unwanted reinjection of acoustic energy.

Using a fully coupled FEM/BEM formulation, we have systematically studied the acoustic radiation damping of sandwich structures. The extent of radiation damping is quantified by the harmonic radiation loss factor relating the radiated sound power to the structure-inherent power. Besides harmonic response analyses, modal analyses of sandwich panels interacting with the surrounding air have also been performed. The underlying nonlinear EVP has been solved using a projection method based on contour integration, resulting in the complex eigenfrequencies and modes of the air-loaded structure. Spurious eigenvalues that arise due to the contour integration are identified by checking the occurrence of the associated air-loaded modes in the range of in vacuo modes. This criterion is also used to check the completeness of the eigenvalue solution. The final eigenvalues provide the shifted (wet) eigenfrequencies, as well as the modal radiation loss factors, of the sandwich structure.

The numerical framework has been applied to three honeycomb sandwich panels subject to various boundary conditions and excitations. The reduction in eigenfrequencies of more than 2% compared to the in vacuo eigenfrequencies clearly indicates the relevance of air-loading when studying the vibration of sandwich structures. Moreover, radiation loss factors of 8% in the coincidence region and larger than 1% across wide frequency ranges demonstrate that the phenomenon of acoustic radiation damping significantly contributes to the overall damping. Furthermore, it is observed that the simply supported, baffled panels exhibit significantly larger radiation damping than the freely suspended, unbaffled panels in the lower frequency range. At higher frequencies, the effects of boundary conditions and excitations are insignificant. The modal radiation loss factors quantifying the radiation damping of each individual structural mode show excellent agreement with the harmonic radiation loss factors corresponding to the diffuse field excitation.

The comparison of the numerical results to the theoretical expressions for radiation damping yields good agreement above the critical frequency. However, commonly used theoretical expressions overestimate the radiation damping at the critical frequency, and they are also inaccurate in the lower frequency range, where the modes of the panel are widely separated and the response of the panel depends on the excitation. Given that sandwich structures exhibit high radiation damping already in the low frequency range, this deficiency of the theoretical expressions underlies the importance of numerical quantification of radiation damping.

Furthermore, we have compared our numerical results to experimentally obtained radiation loss factors found in the literature. While they qualitatively show similar behaviors with an increase toward the coincidence region and subsequent plateaus, we have also observed some significant deviations. In the case of the baffled panels, the deviation in the low frequency range could be explained by the effect of boundary conditions and also by the low modal density of reverberation rooms. A recently proposed procedure based on mobility measurements could resolve the latter issue.⁵³ The deviation in the loss factors of the unbaffled panels in the high frequency range could have its origin in the reinjection of acoustic energy into the panel when testing them in reverberation rooms. While the numerical models assume that all the radiated acoustic energy disappears in the far-field, in fact, part of the acoustic energy in a reverberant room serves as an excitation in addition to the mechanical excitation. The latter should be taken into consideration when experimentally determining radiation damping.

Future research will address the choice of nonlinear eigensolvers for computing air-loaded eigenfrequencies and modes. While the block SS method used in this paper is computationally efficient and achieved errors of order $O(10−5)$ in the considered examples, it also has some disadvantages: Ill-conditioning of the Hankel matrices could result in inaccurate eigenpairs,⁵⁴ and the choice of input parameters could represent a daunting task for the engineer. Other contour integral methods^37,54 or iterative eigensolvers,⁵⁵ in conjunction with in vacuo modes as initial guesses, could prove to be more suitable. Proper benchmarking with regard to nonlinear FEM-BEM EVPs is certainly an issue for future research.

This work was conducted with the support of the German Research Foundation (DFG Project No. MA 2395/15-1) in the context of the priority program 1897 “Calm, Smooth and Smart—Novel Approaches for Influencing Vibrations by Means of Deliberately Introduced Dissipation.”