High-frequency structural-acoustic analysis using an unstructured zero-order energy FEM formulation

2.1. Energy flow equation

The EFA governing differential equation for steady state problem can be developed by considering orthogonal and incoherent waves traveling within the wave-bearing domain of interest and expressed as [11]:

where $e$ is the time- and space-averaged energy density, which characterizes the vibration behavior. $\omega$ is the circular frequency of the analysis, $c_g$ is the group speed, $\eta$ is the hysteresis damping factor, and ${\pi }_{in}$ is the injected power density. The first term on the left side of Eq. (1) represents the transmitted power, and the second term represents the dissipated power. Thus, Eq. (1) expresses an energy flow balance that involves the power flow, energy loss, and injected power density.

In typical vibrating systems, the form of Eq. (1) is slightly different. As enumerated in Table 1, ${\rho }_s$ is the structural density, $E$ and $\nu$ are the Young’s modulus and Poisson ratio, respectively, $D$ is the flexural rigidity, and $c_0$ is the sound speed. Notably, the energy densities in plates and acoustic domains are recognized as surface and bulk densities, respectively.

2.2. Power flow on coupling boundary

Notably, the governing equations enumerated in Table 1 are valid if the structure and the fluid (air) do not affect each other. However, the primary variable $e$ is discontinuous at the boundaries of connecting elements at the junctions for a structural–acoustic coupled system. These discontinuities can originate from the geometry, from multiple components connected together, from the change of material, or from the different medium interface [7], where ever a reflection/refraction of traveling wave occurs. The structural bending, longitudinal, shear motion, and acoustic motion will not affect each other unless structural-structural or structural-acoustic interfaces exist [19].

In this work, two kinds of discontinuities are considered, namely, Plate-Plate (P-P) joint, which represents a junction of plates, and Plate-Acoustic (P-A) joint, which represents an interface of plate and acoustic domain. As deduced by Dong et al. [19], the relationship between the power flow through the discontinuities and the energy density on the boundary of the discontinuities has the same form in both P-P and P-A joints, as:

where $\mathbf{I}$ is an identity matrix, ${\mathbf{c}}_g$ is a diagonal matrix formed by the group speeds of coupling waves, and $\mathbf{\tau }$ is the power transfer coefficient matrix. The dimensions of the matrices are 2×2 for P-A joints, and 3$n$×3$n$ for P-P joints, where $n$ is the number of plates which constitute the joint. In this formulation, all the three wave types in plates are considered.

2.3. P-A joint relationship

For a structural-acoustic coupling problem, the bending vibration induces the acoustic wave, and the acoustic wave also affects the bending vibration. Zhang et al. [20] pointed out that because of the effect of fluid loading on the structure, the governing equation of bending wave in Table 1 should be modified as follows:

where ${\eta }_{rad}$ is the radiation damping, which is defined as:

In Eq. (4), the parameter $\alpha$ can be calculated by:

where $f_c$ is the coincidence frequency, at which the structural bending wave number $k_b$ equals to the acoustic wave number $k$. The bending group speed should be computed with the formulation $c_{gB}=2{\left(D{\omega }^2/\alpha {\rho }_{s}h\right)}^{1/4}$, where $D$ is the flexural rigidity of the plate.

The radiation efficiency ${\sigma }_{rad}$ quantifies the interaction between structure bending and acoustic waves, the formulation proposed by Leppington [21] is used in this work:

where $r=a/b$$\text{(}r\ge \text{1)}$ is the ratio between the characteristic length $a$ and $b$ of the plate, $\mu =k_{sB}/k$ is the wave number ratio, and ${\lambda }_c=c_0/f_c$ is the wavelength of acoustic wave at the coincidence frequency $f_c$.

By the physical meaning of ${\sigma }_{rad}$, Bitsie [8] derived the power transfer coefficients between bending wave in plate and acoustic wave in air, as:

As the result of conservation law, ${\tau }_{bb}=1-{\tau }_{ba}$ and ${\tau }_{aa}=1-{\tau }_{ab}$.

Considering the P-A joint, the form of the matrices in Eq. (2) can be written as:

2.4. Power flow through the P-P joint

Langley and Heron [18] discussed how to compute the power transfer coefficients between an arbitrary number of plates connected to a common beam, using a wave dynamic stiffness matrix method. In this method, the complete set of transmission coefficients is written in the form of ${\tau }_{pr}^{ij}$, where respectively $i$ and $p$ represent the carrier plate and waveform of the incident wave, and $j$ and $r$ represent those of the generated wave. If $i=j$, $\tau$ is recognized as a reflection coefficient. Their published paper has several errors, which are enumerated in Appendix.

If considering a P-P joint, which consists of two plates denoted by $i$ and $j$, the matrices in Eq. (2) are in the form of:

3.1. Calculation of ${\mathit{Q}}_{\mathit{i},\mathit{j}}$ in 2D domain

Fig. 1(a) shows the schematic of two adjacent control volumes of the unstructured grids in a plate, whose centroids are denoted as $Z_i$ and $Z_j$, respectively, and their common boundary is $\stackrel{-}{AB}$. In Fig. 1(b), ${\mathbf{t}}_{i,j}$ and ${\mathbf{n}}_{i,j}$ are the tangential and normal unit vectors of $\stackrel{-}{AB}$, respectively, and ${\mathbf{\tau }}_{i,j}$ and ${\mathbf{\nu }}_{i,j}$ are the unit vectors of the one starting at $Z_i$ pointing to $Z_j$ and its vertical one, respectively. $\theta$ is the angle between ${\mathbf{\tau }}_{i,j}$ and ${\mathbf{t}}_{i,j}$. $d_i$ and $d_j$ are the distances from each centroid of $Z_i$ and $Z_j$ to $\stackrel{-}{AB}$ , respectively. $M$ is the intersection of $\stackrel{-}{AB}$ and $\overline{Z_i{Z}_j}$.

Fig. 1Two adjacent control volumes and vectors definition on the common boundary in 2D domain

a) Volumes $Z_i$ and $Z_j$

b) Vectors definition on $\stackrel{-}{AB}$

In 2D plates, the net power flow density $q_{i,j}$ is assumed to be uniform on the common boundary $\overline{AB}$; thus, the power flux $Q_{i,j}$ can be easily expressed as:

where $L_j$ is the length of $\overline{AB}$, and $h$ is the thickness of the plate. In this work, the triangular grid is used for the plates; thus, the number of adjacent control volume for each grid, $m$ in Eq. (12), is equal to 3, and $q_{i,j}$ is in the form of:

Note that the vectors ${\mathbf{\tau }}_{i,j}$ and ${\mathbf{t}}_{i,j}$ are non-collinear, so ${\mathbf{n}}_{i,j}$ can be decomposed into ${\mathbf{n}}_{i,j}=a{\mathbf{\tau }}_{i,j}+b{\mathbf{t}}_{i,j}$, where:

then the power flux density is resolved into these two directions, and a finite difference is used in each direction for numerical approximation, as in $V_i$ one can obtain:

and similarly in $V_j$, the opposite net flow density is:

The symbol $\left|\overline{P_1{P}_2}\right|$ indicates the distance of the two involved points, $P_1$ and $P_2$.

The power flux through the common boundary $\overline{AB}$ is continuous, which yields:

Substituting Eqs. (16), (17) to Eq. (18), $e_M$ can be easily solved as:

then substituting Eq. (19) into Eq. (16), and notice that:

One can eliminate $u_M$ and the approximation formula of $q_{i,j}$ can be written as:

From Eqs. (13) and (21), the power flux through one boundary of the element $Z_i$ is solved. This procedure is repeated $m$ times and is superposed. Then, the superposition of all the power flux out of the element $Z_i$ can be worked out, which indicates that the first term of Eq. (10) is finally solved.

3.2. Calculation of ${\mathit{Q}}_{\mathit{i},\mathit{j}}$ in 3D domain

As shown in Fig. 2, $Z_i$ and $Z_j$ denote the centroids of two adjacent tetrahedral grids in 3D domain, whose common boundary is a triangle denoted by ${∆}_{BCD}$, and $E$ is the centroid of ${∆}_{BCD}$.

First, two functions should be introduced:

$V_{abcd}$ represents the volume of a tetrahedron whose vertices are $a$, $b$, $c$ and $d$. An area vector, ${\mathbf{S}}_{a,b,c}$, is in the form of:

The use of commas emphasizes that the subscripts $a$, $b$ and $c$ in Eq. (23) are directive, which is different from those in Eq. (22). The direction of vector ${\mathbf{S}}_{a,b,c}$ is determined by right hand rule, starting from $a$ to $b$ and then to $c$.

The governing equation of acoustic domains has the same form as that of plates (listed in Table 1), but the power flow flux from $Z_i$ to $Z_j$, $Q_{i,j}$, is divided into three parts in 3D solution:

The gradient, $\nabla e\text{,}$ in each part, ${\left(\nabla e·S\right)}_{B,C,E}$ for example, is averaged in a secondary tetrahedron, $T_{Z_{i}BCE}$:

Otherwise, a similar approximation in $T_{Z_{j}BCE}$ while calculating $Q_{j,i}$ will lead to:

To guarantee the conservativeness of the power flux through the common boundary ${∆}_{BCD}$:

Substituting Eqs. (25), (26) into Eq. (27), $e_E$ can be eliminated, and results in:

where the scalar coefficients can be solved as:

and:

Similar to Eq. (28), the other two parts of power flux, ${\left(\nabla e·S\right)}_{C,D,E}$ and ${\left(\nabla e·S\right)}_{D,B,E}$, can be worked out. Once they are all determined, submit them into Eq. (24), then $Q_{i,j}$ in 3D domain can be obtained.

Fig. 2Two adjacent control volumes in 3D domain

Fig. 3Interpolation diagram for computing vertex weighted average e0 from surrounding cell-centered ei

3.3. Node interpolation

When calculating $Q_{i,j}$ in both 2D and 3D domains, the variable at the vertex, $e_0$, is unknown, as shown in Fig. 3; it is evaluated by a weighted average of the surrounding cell-centered values, $e_1$∼$e_N$, as [24]:

in which the weight functions $W_i$ are selected to ensure zero Laplacian of the linear function in the shadow region (Fig. 3). $N$ is the number of the surrounding cells.

Aiming at this, Lagrangian multipliers, ${\lambda }_x$, ${\lambda }_y$ and ${\lambda }_z$, are employed, and the weight function is expressed as:

where $\left(x_i,y_i,z_i\right)$ are the spatial locations of the centroids of the cells, and $\left(x_0,y_0,z_0\right)$ is that of the vertex. The Lagrangian multipliers are computed by the local geometric parameters as:

where:

The interpolation coefficients are computed as:

3.4. Assembling the global equation

First, $Q_{i,j}$ is initially computed for each common boundary in 2D and 3D domains, and if the boundary is on any joint, $Q_{i,j}$ is computed using Eq. (2) instead. Then, the results are summarized into Eq. (10) and all the element stiffness equations are worked out. Lastly, the corresponding coefficients of each control volume are gathered, and the global equation can be written in a general form as:

where $N_{\mathrm{d}\mathrm{o}\mathrm{f}}=3N_P+N_A$, $N_P$ and $N_A$ are the numbers of control volumes in plates and air cavities, respectively. The notations $B$, $L$, $S$, and $A$, represent the bending, longitudinal, shear, and acoustic waves, respectively. The diagonal matrices BB, LL, and so on, are assembled from $Q_{i,j}$ of each wave domain, and the other matrices are from the power transfer of the junctions. $\mathit{e}$ is the energy density vector, and $\mathbf{\Pi }$ is the vector of corresponding power injected into every element.

4.1. Seven-plate engine foundation model

Several researchers have investigated the 1:2.5 scale model of the diesel propulsion engine foundation and the hull structure of a frigate ship to study the solution of high-frequency structural vibration based on energy flow approach. Lyon [26] divided this model into 7 and 12 substructures and studied the significance of the in-plane energy in energy transmission by the SEA method. Vlahopoulos et al. [3] applied EFEM on this seven-plate structure and obtained consistent results with those of Lyon.

Fig. 4Seven-plate engine foundation model

In this paper, the material properties and the dimensions of plates of the seven-plate model are the same as those in Ref. [3]. The uEFEM⁰ model comprises 660 nodes, 1170 triangular control volumes, and 27 P-P joints, as shown in Fig. 4.

The analysis is carried out in thirteen 1/3 octaves that cover from 1.25 kHz to 20 kHz, as listed in Table 2. In each octave band, a unit bending power of 1 W is injected at the center of plate 1, whereas it is applied on the top plate in the corresponding SEA model. On all the boundaries of the system, the boundary condition is set as $\mathbf{I}=-\left(c_g^2/\eta \omega \right)\nabla e\cdot \mathbf{n}=0$.

Fig. 5 shows the comparison of the total bending energy of plates 1, 4, and 7. The bending energy ratios between plates 7 to 1 and plates 2 to 1 are presented in Figs. 7(a) and 7(b), the results with the legend of ‘EFEA’ were developed by Vlahopoulos [3], and the experimental results were provided by Lyon [26].

It can be demonstrated that:

1) The bending energy results of the incident plate (Plate 1) fit well, whereas those of other plates have slight difference, which is less than 3 dB. This finding is attributed to the plates except the incident plate; the energy is induced by the joints, and different treatments are applied with the joints in uEFEM⁰ and SEA.

2) As observed in Figs. 7(a) and 7(b), the developed uEFEM⁰ shows the best overall correlation to the experimental data.

3) When only the bending mechanism is included in the formulation, power can be transferred from plates 1 to plate 7 only from the bending motion, while in fact an important amount of power is transferred from the in-plane motion of plates 2 and 5. Thus, the SEA result with only bending wave included shows a significant error on the energy ratio between plate 7 to plate 1.

4) The uEFEM⁰ formulation is based on a FEM mesh of the structure. Therefore, the spatial distribution of the energy density can be identified easily (Fig. 6).

Fig. 5Bending energy results in plates 1, 4 and 7 (results from uEFEM0, SEA/AutoSEA)

Fig. 6Distribution of bending energy density on the seven-plate model at 1.25 kHz (result from uEFEM0, and the unit: dB, ref = 1×10-12 J/m2)

Fig. 7Bending energy ratio between plates (results from uEFEM0, SEA/AutoSEA, and literature)

a) Plates 2 to 1

b) Plates 7 to 1

4.2. Simplified passenger vehicle model

A passenger vehicle is simplified to seven structural plates with an internal air cavity. The plates are made of aluminum, with a uniform thickness $h=$5 mm, and the internal cavity is filled with air. The dimensions and mesh of the vehicle are both shown in Fig. 8. The model comprises 468 structural plate elements and 1220 acoustic air elements, with 79 P-P joints and 468 P-A joints. The injected powers, 0.25 W each in any octave, are applied at about the positions of four wheels, which are (0.6, 0.3), (0.6, 1.2), (2.5, 0.3), and (2.5, 1.2) in the coordinate system $XOY$ set on the base plate. In the SEA model, a unit power of 1 W is applied at the base plate. The side plates are examples of plates with irregular shape, which would be difficult to deal with EFEM⁰ developed by Wang [5], but easy with this presented uEFEM⁰.

The material properties are listed in Table 3, and the analysis frequency range is the same as the prior numerical example, which is listed in Table 2.

Fig. 8uEFEM0 model of simplified vehicle

Figs. 9(a) and 9(b) show the comparison of results from uEFEM⁰ and AutoSEA, including the bending energy of plates and sound pressure level (SPL) of the cavity. The spatial variation of the energy density is presented in Fig. 10. From the analysis of the simplified vehicle, the following conclusions can be summarized:

1) Similar to the behaviors of uEFEM⁰ and SEA analysis in the prior numerical example, the total bending energy of the incident plate shows a good agreement, whereas the results in other plates are slightly deviated. The result of SPL fits well for these two methods.

2) As observed in Fig. 10, the bending energy level is highest at the loading points and decreases as the distance of the loading points increases. The in-plane energy takes its peak in the side plates, which have the same direction with the loading forces. This phenomenon explains the importance of the in-plane energy in the energy transmission in such structures including ‘L’ angles. The SPL shows a significant value at the corners between the base and side plates. This result is due to a number of plate elements with relatively high bending energy that radiate energy into the air cavity at those corners.

3) The simplified passenger vehicle application demonstrates how the uEFEM⁰ formulation can be utilized as a simulation tool for structural-acoustic car applications, thus it may potentially be further developed to be an NVH tool.

Fig. 9Comparison of bending energy in several plates and sound pressure level in air cavity

a) Bending energy

b) Sound pressure level