Acoustic modal analysis and harmonic analysis are of paramount importance in Automotive, Aerospace, Room acoustics applications to assess acoustic comfort. There are many methods available to conduct acoustic analysis such as Finite Element Method (FEM), Boundary Element Method (BEM), Statistical Energy Analysis (SEA) and Modal Impact / Shaker testing. In the case of low frequency applications with finite domain(interior acoustic problems), FEM found to be a more suitable choice. In FEM, the weighted residual form of the second-order wave equation in the frequency domain is discretized and acoustic pressure is approximated with shape functions in an isoparametric sense. BEM would be a good choice for computing with infinite domains whereas SEA would give reliable results for mid and high frequency range.
The following page is derived from a project report written by Ramsubramanian Pazhanisamy and Ricky Aristio for the course Modelling and Simulation in Structural Mechanics in the winter term 2019/20.
Acoustic Finite Element Method
Strong form
The three dimensional wave equation for sound waves is given by
| \begin{equation} \frac{1}{v^{2}} \frac{\partial^{2} p(x, y, z,t)}{\partial t^{2}}=\frac{\partial^{2} p(x, y, z,t)}{\partial x^{2}}+\frac{\partial^{2} p(x, y, z,t)}{\partial y^{2}}+\frac{\partial^{2} p(x, y, z,t)}{\partial z^{2}}. \end{equation} |
This wave equation can be transformed into the frequency domain
| \begin{equation} \nabla^{2} p(x, y, z)+k^{2}, p(x, y, z)=-j \rho_{0} \omega \cdot q(x, y, z) \end{equation} |
The second order acoustic wave equation in the frequency domain is also called \textit{Helmholtz equation}. In the equation k=\omega / c=2 \pi f / c is the wavenumber, c is the speed of sound and \rho_{0} is fluid density.
For a unique definition of the field pressure in V , either one of boundary condition must be specified at each position of the domain boundary \Omega=\Omega_{p} \cup \Omega_{v} \cup \Omega_{Z}:
Pressure boundary condition (Dirchlet):
p=\bar{p} \quad \text { on } \Omega_{p} Acoustic velocity boundary condition (Neumann):
v_{n}=\frac{j}{\rho_{0} \omega} \frac{\partial p}{\partial n}=\bar{v}_{n} Impedance boundary condition (Robin) :
p=\bar{Z} \cdot v_{n}=\frac{v_{n}}{\bar{A}}=\frac{j \bar{Z}}{\rho_{0} \omega} \frac{\partial p}{\partial n}=\frac{j}{\rho_{0} \omega \bar{A}} \frac{\partial p}{\partial n}
Weak form - weighted residual formulation
Multiplying the strong form with the test function and then integrating them into the whole domain results in weak form.
| \int_{V} \widetilde{p}\left(\nabla^{2} p+k^{2} p+j \rho_{0} \omega q\right) \cdot d V=0 |
Applying Gauss divergence theorem and reformulation leads to a weighted residual statement of Helmholtz equation,
| \int_{V}(\vec{\nabla} \widetilde{p} . \vec{\nabla} p) . d V-\omega^{2} \int_{V}\left(\frac{1}{c^{2}} \widetilde{p} p\right) . d V=\int_{V}\left(j \rho_{0} \omega \widetilde{p} q\right) . d V-\int_{\Omega}\left(j \rho_{0} \omega \widetilde{p} \vec{v} . \vec{n}\right) . d \Omega |
FEM Approximation
The pressure field is approximated with shape functions in an isoparametric way as follows,
| p(x, y, z) \approx \hat{p}(x, y, z)=\sum_{i=1}^{n} N_{i}^{e}(x, y, z) \cdot a_{i} \quad(x, y, z) \in V_{e} =[N] \cdot\left\{\hat{p}_{i}\right\} |
where n is the total number of nodes in the model, [N] is the shape function matrix and \left\{\hat{p}_{i}\right\} the vector of degree of freedoms.
The pressure gradient only contains normal derivatives and does not have cross derivatives because an acoustic fluid does not support shear,
| \vec{\nabla}{\tilde{p}}=[\partial][N] \cdot\left\{\tilde{p}_{i}\right\}=[B]\left\{\tilde{p}_{i}\right\} |
where [B] is called B operator matrix and does not have cross derivative terms. By discretization, the weighted residual helmhotz equation can be transformed as follows
| \left\{\widetilde{p}_{i}\right\}^{T} \cdot\left([K]+j \omega[C]-\omega^{2}[M]\right) \cdot\left\{\hat{p}_{i}\right\}=\left\{\widetilde{p}_{i}\right\}^{T} \cdot\left(\left\{Q_{i}\right\}+\left\{V_{n i}\right\}+\left\{P_{i}\right\}\right) |
each term of the above equation are explained in the following:
Acoustic stiffness matrix:
| \left\{\widetilde{p}_{i}^{e}\right\}^{T} \cdot\left[K^{e}\right]\left\{\hat{p}_{i}^{e}\right\} = \left\{\widetilde{p}_{i}^{e}\right\}^{T} \cdot\left(\int_{V_{e}}\left(\left[B^{e}\right]^{T} \cdot\left[B^{e}\right]\right) \cdot d V\right) \cdot\left\{\hat{p}_{i}^{e}\right\} |
Acoustic damping matrix: (Not considered in this implementation)
| C_{i j}=\int_{\Omega_{Z}}\left(\rho_{0} \bar{A} N_{i} N_{j}\right) d \Omega |
Acoustic mass matrix:
| -\omega^{2}\left\{\widetilde{p}_{i}\right\}^{T} \cdot[M]\left\{\hat{p}_{i}\right\} =-\omega^{2}\left\{\widetilde{p}_{i}\right\}^{T} \cdot\left(\int_{V}\left(\frac{1}{c^{2}}[N]^{T} \cdot[N]\right) \cdot d V\right) \cdot\left\{\hat{p}_{i}\right\} |
Input displacement boundary condition:
| \left\{\widetilde{p}_{i}\right\}^{T} \cdot\left\{P_{i}\right\} =\left\{\tilde{p}_{i}\right\}^{T} \cdot\left(\int_{\Omega_{u}}\left(\rho_{0} \omega^{2}[N]^{T} \vec{u} . \vec{n}\right) . d \Omega\right) |
Pressure boundary condition: (Not considered in this implementation)
| P_{i}=\int_{\Omega_{p}}\left(-j \rho_{0} \omega N_{i} \vec{v} . \vec{n}\right) . d \Omega |
Acoustic source excitation term: (Not considered in this implementation)
| \left\{\widetilde{p}_{i}\right\}^{T} \cdot\left\{Q_{i}\right\} = \left\{\tilde{p}_{i}\right\}^{T} \cdot\left(\int_{V}\left(j \rho_{\theta} \omega[N]^{T} q\right) \cdot d V\right) |
Out of all the terms mentioned above only acoustic mass matrix, acoustic stiffness matrix and input pressure boundary condition terms are considered in this scope of work, other terms are neglected.
FEM Modal Analysis formulation
| \begin{equation} \left(-\omega^{2} M + K \right) \boldsymbol{\phi}=0 \end{equation} |
where \omega denotes natural frequency(eigen value) and \phi is the mode shape(eigen vector).
FEM Harmonic Analysis formulation
| \begin{equation} \hat{\mathbf{p}}(\omega)=\left(-\omega^{2} \mathbf{M}+\mathbf{K}\right)^{-1} \hat{\mathbf{P}}(\omega) \end{equation} |
where \hat{\mathbf{P}}(\omega) is the right hand side representing the input displacement boundary condition
Results
A modal and a harmonic analysis of a cavity model using 2D and 3D acoustic finite elements that are implemented in Kratos have been performed. The result are compared to a 1D Analytical solution.
Geometric Modelling and Load Conditions
A cavity of length l = 1m is filled with air with density \rho = 1.225 kg/m^3 and speed of sound c = 340 m/s. Moreover, the 2D and 3D model of the cavity are created with height h = 0.2 m and thickness t = 0,2 m respectively.
At the end of the cavity, a time-harmonic displacement is imposed with frequency \omega and amplitude u = 0.001 m. There are three different types of the imposed displacement. For the 1D model, the displacement is only a point at the end of the cavity. For the 2D model, the displacement is assigned in two different approaches, first as a point load at the right corner of the cavity and second as a line load at the end of the cavity. Furthermore, for the 3D model, the displacement is assigned in three different approaches, as a point load at the right corner of the cavity, as a line load at the top right of the cavity and as an area load at the end of the cavity. And finally, the rest of the cavity is rigid walled.
Based on the figure above, we can distinguish three different cases. First, the blue models are comparable because they only have a longitudinal mode in one direction due to the distributed load at all directions at the end of the cavity. Second, the green models are comparable, because there are two types of modes, the longitudinal and transversal modes since we only put distributed load in one direction. Third, for the yellow model, similarly, we expect that it has three different types of modes, one longitudinal and two transverse modes.
Modal Analysis
In the following, the modal analysis result from 2D and 3D models will be presented and compared to the analytical eigenfrequency of the 1D cavity model. The figure below shows the type of modes of the 2D model.
As expected, compared to the 1D model, the 2D model has not only longitudinal mode but also transverse modes.
The figure below shows the type of modes of the 3D model.
Similarly, the 3D model has another transverse mode due to the addition of the third direction.
Now, all eigenfrequencies from the 2D and 3D cavity are compared with the eigenfrequencies of the 1D analytical cavity.
Since the 1D cavity has only longitudinal modes, we consider only longitudinal modes of the 2D and 3D cavity (neglect the transverse modes). From the table above, it shows that for each eigenfrequency in the 1D cavity, there are corresponding modes in 2D and 3D cavity and mostly those modes are also longitudinal. As the mode increases, the mode shape also increases by a half-sine wave. For example, the mode shape of the first eigenfrequency can be represented as a half-sine wave, then the second eigenfrequency can be represented as one sine wave and so on.
The following figure shows the first mode of the 2D and 3D models. It represents the longitudinal modes very well and the eigenfrequencies are close to the 1D analytical solution.
As mentioned earlier, almost all the modes here are aligned. However, there is an outlier here in the fifth mode. Although the value of the eigenfrequencies is close between 1D, 2D, and 3D model, the mode shape does not represent longitudinal modes as can be shown in the figure below.
Both 2D and 3D models have multiple mode shapes in one eigenfrequency value (two mode shapes in 2D model and three mode shapes in the 3D model) and those mode shapes are not longitudinal modes. Moreover, if we consider the pattern of the mode shapes, it still aligns with the longitudinal modes wherein fifth modes we have 2.5 wavelengths.
There is a reason why the mode shapes are coupled in 2D and 3D cavity in the fifth mode. First, we consider that in 2D and 3D models, we have modes in two and three different directions. Next, coincidentally the length of the model is five times the height and thickness. If we check the modes in each direction one by one at the fifth eigenfrequency, for instance in the 2D model, we have 2.5 sine wave in longitudinal direction and half-sine wave in the transverse direction. Thus, the superposition of those two modes will lead to the combination of the longitudinal modes and the transverse modes, since both of the waves are in phase.
Harmonic Analysis
In the following, harmonic analysis' results of two dimensional and three dimensional acoustic cavities are going to be discussed. Three load cases such as point load, line load, and area load are applied.
One dimensional axial modes case
In this subsection, results of two dimensional acoustic cavity with line load at the right edge and three dimensional acoustic cavity with area load at the right hand side are compared. In both cases only the longitudinal(axial) modes are excited.
By comparing the two figures above, it can be concluded that both models excite the same mode at a particular frequency of interest.
The comparison graph above shows all longitudinal modes which are excited by the axial perturbation displacement in the frequency range. Frequency range up to 5000 rad/s, both 2D and 3D models have good match whereas at higher frequencies 3D model behaves less stiff compared to the 2D models as can be seen from the shift of natural frequencies. The acoustic pressure amplitude cannot be compared between 2D and 3D models, because both line and area load are only applied point-wise.
Two dimensional modes case
In this subsection, results of two dimensional acoustic cavity with point load at the right edge and three dimensional acoustic cavity with line load at the right hand side are compared. Here both the longitudinal(axial) modes and transverse modes are excited.
As can be seen from the two figures above both models excite the same mode at a specific frequency of interest. In the figure below, as expected both axial and transverse modes are excited due to the transverse loading. The same argument from the previous subsection holds here regarding the shift of natural frequencies at high frequency range.
Three dimensional modes case
In this subsection, the results of three dimensional acoustic cavity with a point load at the right hand side are discussed. This harmonic perturbation excites all available modes in the frequency range of excitation. The figure below shows three dimensional mode excited at a frequency step.
In the figure above, high modal density can be observed in the frequency range between 5000 rad/s and 10000 rad/s signifying the presence of a high number of closely spaced modes.
Mesh dependency analysis
Two dimensional acoustic cavity with point load at the top-right node have been successively refined starting from 500 elements until 3920 elements, then their frequency responses are compared. At low frequency range (below 5000 rad/s), natural frequencies are not shifting as refinement goes up, signifying the sufficiency of a low number of elements. At high frequency range, the wavelength reduces hence the model requires high resolution. As expected natural frequencies are shifting in the high frequency range as the mesh is refined.
