Coupling of A Kirchhoff Plate with An Acoustic Cavity with The Wave Based Method

Gündoğan, Elif Ali Sorkun Unbekannter Benutzer (ge32nis)

The wave based method (WBM) is a numerical method that uses weighted wave functions to approximate steady-state field response in boundary value problems. WBM is useful, particularly for mid to high frequency vibro-acoustic problems, where element based methods such as FEM require large amount of elements to capture the short wavelength [1]. To obtain the frequency response in the wave-based method, a harmonic analysis is performed, where the excitation is modeled as a sinusoidal function with a specific frequency. The differential equations of the wave-based method are then solved for this frequency, and the solution is used to compute the response of the structure. An example problem is to select, scale, and add wave functions to the following pressure contour on the left such that the resultant satisfies not only the governing equation but also the prescribed boundary conditions as shown on the right.

The purpose of this project is to extend an existing MATLAB script that solves acoustic problems using WBM to account for all typical boundary conditions in a vibro-acoustic problem, including a Kirchhoff plate boundary, which requires a coupling between the acoustic cavity and plate. Subsequently, the script is used to solve several vibro-acoustic problems, including a plate-cavity system with a flexible plate and an acoustic cavity. The results are then compared with corresponding finite element method (FEM) solutions, which highlight the efficiency of the WBM.

Theoretical Basis

This section discusses mathematical models and idealisations used to formulate governing equations for both the acoustic part as well as the structural part. Only general derivations are included while detailed derivations are omitted and deferred to the Appendix.

Helmholtz Equation

Acoustic responses in a fluid are usually regarded as small perturbations to an ambient reference state, which is assumed to be constant [2]. Consequently, pressure, density, and velocity can be decomposed into their ambient and perturbation values as shown in the following illustration and expressions.

$\begin{aligned} p \left( \mathbf{x}, t \right) &= p_{0} + p' \left( \mathbf{x}, t \right) \\ \rho \left( \mathbf{x}, t \right) &= \rho_{0} + \rho' \left( \mathbf{x}, t \right) \\ \mathbf{v} \left( \mathbf{x}, t \right) &= \mathbf{v}_{0} + \mathbf{v}' \left( \mathbf{x}, t \right) \end{aligned}$

Decomposition of pressure, velocity, and density to their ambient and perturbation simplifies mass and momentum conservation equations since the ambient terms are constants and therefore their derivatives are zero. Furthermore, small perturbations justifies linearisation of both equations. The two simplified equations, without external ambient flow, are given below.

$\begin{aligned} \dot{\rho}' &= - \nabla \cdot \left( \rho_{0} \mathbf{v}' \right) + \rho_{0}q' \\ \rho_{0} \dot{\mathbf{v}}' &= - \nabla p' \end{aligned}$

Because there are three unknowns i.e. pressure, velocity, and density, aforementioned two conservation equations are not sufficient and an additional equation is necessary. Typically, the adiabatic relationship between pressure and density is used. Again, linearisation of the relationship is justified by the small perturbations. This third equation is given below.

$\begin{array}{l}\displaystyle p' = \rho ' \cdot \frac{\gamma p_{0}}{\rho_{0}}\end{array}$

Differentiation of the mass conservation equation with respect to time, followed by substitutions of the other two equations yields the governing equation for the acoustic part. The equation as well as its Fourier transformation are given below. The latter is used to model the acoustic domain in this project.

$\begin{array}{l}\displaystyle -\rho_{0} \dot{q}' = \nabla^{2} p' -\frac{1}{c^{2}} \ddot{p}' \phantom{x} \xrightarrow{\phantom{x}\mathscr{F}(\phantom{x})\phantom{x}} \phantom{x} -j\rho_{0}\omega Q = \nabla^{2} P + \frac{\omega^{2}}{c^{2}} P\end{array}$

Plate Equation

According to Kirchhoff’s plate theory, the in-plane and out-of-plane motions in a thin, flat plate are decoupled and are governed by separate dynamic equations [2]. For the Kirchoff plate, the following assumptions are used:

Plate material is elastic and its elastic quantities are governed by Hooke's law.
Plate material is homogenous and isotropic, with elasticity modulus, density, and Poisson ratio given by $\begin{array}{l}E,\rho_{s},\nu\end{array}$ respectively.
Plate thickness $\begin{array}{l}t\end{array}$ is small in comparison to other dimension of the plate.
Plate displacement is small compared to plate thickness.
Sections normal to mid-plane remain normal to it after deformation.
Normal stress in direction normal to plate surface is negligible.

The plate in coupled vibro-acoustic systems is loaded by the acoustic pressure and external mechanical load as shown in the following figure.

The steady state dynamic equation for this plate is given by the following equation. To take damping into account, a visco-elastic damping model is assumed. Thus, the damping effect is modelled by defining a complex elasticity modulus; $\begin{array}{l}E^{*} = E(1+j\eta)\end{array}$ , where $\begin{array}{l}\eta\end{array}$ is the material loss factor. In addition, the term $\begin{array}{l}\frac{1}{D}F(x,y)\end{array}$ is added to the right hand side of the equation, to take external mechanical force into account. The derivation of this formula is quite long and is deferred to the Appendix.

$\begin{array}{l}\displaystyle \nabla^{4}w(x,y) - k_{b}^{4} w = \frac{1}{D}P(x,y) + \frac{1}{D}F(x,y), \phantom{xx} D = \frac{E(1+j\eta)t^{3}}{12\left(1-\nu^{2}\right)}; \phantom{x} k_{b} = \sqrt[4]{\frac{\rho_{s}t\omega^{2}}{D}}\end{array}$

Boundary Conditions

For interior coupled vibro-acoustic systems, the following types of boundary conditions may occur. They are also shown in the following diagram.

Imposed pressure (Dirichlet boundary condition)
Imposed normal velocity (Neumann boundary condition)
Normal impedance boundary condition (mixed boundary condition)
Normal displacement continuity at the fluid-structure interface

The relations between values at the boundaries for these are given by the following equations, respectively.

$\begin{array}{l}\begin{array}{rmll} P(\mathbf{r}) &=& \bar{p}(\mathbf{r}) & \forall \mathbf{r} \in \Omega_{p} \\ v_{n}(\mathbf{r}) &=& \bar{v}_{n}(\mathbf{r}) & \forall \mathbf{r} \in \Omega_{v} \\ v_{n}(\mathbf{r}) &=& P(\mathbf{r}) / \bar{Z}(\mathbf{r}) & \forall \mathbf{r} \in \Omega_{z} \\ v_{n}(\mathbf{r}) &=& j\omega w(\mathbf{r}) & \forall \mathbf{r} \in \Omega_{s} \end{array}\end{array}$

Design of Numerical Solution

The first half of this section discuss the construction of approximate solution for both acoustic pressure and plate displacement. The second half is dedicated to define the "best" approximate solution in fulfilling the prescribed boundary conditions, as well as reformulation of the problem into a system of linear equations to be solved.

Approximate Solution to Acoustic Pressure

The Helmholtz equation for point acoustic source is given below, where $\begin{array}{l}(x_{q},y_{q})\end{array}$ denote the position of the acoustic source.

$\begin{array}{l}\displaystyle -j\omega |Q|\cdot \delta \left(x-x_{q},y-y_{q}\right) = \nabla^{2} P(x,y) - k^{2} P(x,y), \phantom{xx} k = \omega / c\end{array}$

Using indirect Trefftz method, an approximate solution is constructed from its particular solution and a finite set of its homogenous solutions. An analytical particular solution is available and given by a scaled Hankel function of the second kind with order zero. Hankel function of the second kind itself is a combination of Bessel function of first and second kind with the same order. The particular solution is given by the following.

$\begin{array}{l}\displaystyle P_{q}(x,y) = \frac{\rho_{0}\omega |Q|}{4} \cdot H_{0}^{(2)} \left( k \cdot \sqrt{\left(x-x_{q}\right)^{2}+\left(y-y_{q}\right)^{2}} \right), \phantom{xx} H_{0}^{(2)}(z) = J_{0}(z) - jY_{0}(z)\end{array}$

While it is beyond the scope of this project to find all homogenous solutions of the equation, it is known that functions of the following forms are in the homogenous solution space.

$\begin{aligned} \Phi_{a1,i} (x,y) &= \cos{ \left( k_{xa1,i}x \right) } \cdot \exp{ \left[ -jk_{ya1,i} \cdot \left( y - f_{ya1,i}L_{y} \right) \right] }, \phantom{xx} k_{xa1,i}^{2}+k_{ya1,i}^{2} = k^{2} \\ \Phi_{a2,i} (x,y) &= \cos{ \left( k_{ya2,i}y \right) } \cdot \exp{ \left[ -jk_{xa2,i} \cdot \left( x - f_{xa2,i}L_{x} \right) \right] }, \phantom{xx} k_{xa2,i}^{2}+k_{ya2,i}^{2} = k^{2} \end{aligned}$

The pressure can then be approximated as a sum of the particular solution and a linear combination of some homogenous solutions as shown below. Through this construction, the approximate solution automatically fulfil the Helmholtz equation. Subsequently, it is only necessary to scale each homogenous solution such that the resulting approximate solution fulfils the prescribed boundary conditions the "best".

$\begin{array}{l}\displaystyle P(x,y) \approx P_{q}(x,y) + \sum_{i=1}^{n_{a}}\hat{p}_{1,i}\Phi_{a1,i}^{*}(x,y) + \sum_{i=1}^{n_{a}}\hat{p}_{2,i}\Phi_{a2,i}^{*}(x,y)\end{array}$

Approximate Solution to Plate Displacement

The plate equation for mechanical point load and pressure load is given below, where x' is a local coordinate system with its origin at one end of the plate.

$\begin{array}{l}\displaystyle \frac{d^{4}w \left(x'\right)}{dx'^{4}} - k_{b}^{4} w \left(x'\right) = \frac{|F|}{D} \cdot \delta\left(x'-x'_{F}\right) + \frac{1}{D}P\left(x'\right)\end{array}$

In the same manner as acoustic pressure, an approximate solution to plate displacement is constructed from its particular solutions and a finite set of its homogenous solutions. Analytical particular solution to the mechanical point load is available and is given below.

$\begin{array}{l}\displaystyle w_{F}\left(x'\right) = -\frac{j|F|}{4Dk_{b}^{3}} \left[ \exp{\left( -jk_{b}\left|x'-x'_{F}\right| \right)} - j\exp{\left( -k_{b}\left|x'-x'_{F}\right| \right)} \right]\end{array}$

The pressure load contains several homogenous solution to the Helmholtz equation, each has an associated particular solution for the plate equation. Analytical particular solutions are available and each of them assume the following form.

$\begin{array}{l}\displaystyle w_{ai,j}\left(x'\right) = -\frac{j}{4Dk_{b}^{3}} \int_{0}^{L} \Phi_{ai,j} \left(\xi\right) \cdot \left[ \exp{\left(-jk_{b}\left|x'-\xi\right|\right)} - j\exp{\left(-k_{b}\left|x'-\xi\right|\right)} \right] \phantom{.}d\xi\end{array}$

Although slightly more complicated, analytical particular solution for pressure's particular solution is also available and given below. The Hankel function of second kind with order zero appear again here.

$\begin{aligned} w_{q}\left(x'\right) =& -\frac{j\rho_{0}\omega|Q|}{16Dk_{b}^{3}} \int_{0}^{L}H_{0}^{(2)} \left[ k \cdot \sqrt{\left(x_{0}+\xi\cdot\cos{\alpha}-x_{q}\right)^{2}+\left(y_{0}+\xi\cdot\sin{\alpha}-y_{q}\right)^{2}} \right] \cdot \exp{\left(-jk_{b}\left|x'-\xi\right|\right)}\phantom{.}d\xi \\ & -\frac{\phantom{j}\rho_{0}\omega|Q|}{16Dk_{b}^{3}} \int_{0}^{L}H_{0}^{(2)} \left[ k \cdot \sqrt{\left(x_{0}+\xi\cdot\cos{\alpha}-x_{q}\right)^{2}+\left(y_{0}+\xi\cdot\sin{\alpha}-y_{q}\right)^{2}} \right] \cdot \exp{\left(-\phantom{j}k_{b}\left|x'-\xi\right|\right)}\phantom{.}d\xi \end{aligned}$

Unlike the Helmholtz equation that has infinitely many homogenous solutions, the plate equation is a linear ODE with order four. Consequently, it only has four homogenous solutions. The four homogenous solutions are given by the following.

$\begin{array}{l}\displaystyle \Psi_{s}\left(x'\right) = \exp{\left(-j^{s}k_{b}x'\right)}, \phantom{x} \forall \phantom{x} s\in\{1,2,3,4\}\end{array}$

The plate displacement can then be approximated as a sum of the particular solutions and linear combination of the four homogenous solutions as shown below. Through this construction, the approximate solution automatically fulfil the plate equation. Subsequently, it is only necessary to scale the four homogenous solution such that the resulting approximate solution fulfil the prescribed boundary conditions the "best".

$\begin{array}{l}\displaystyle w\left(x'\right) \approx w_{F} \left(x'\right) + w_{q} \left(x'\right) + \sum_{i=1}^{n_{a}}\hat{p}_{a1,i}w_{a1,i}^{*} \left(x'\right) + \sum_{i=1}^{n_{a}}\hat{p}_{a2,i}w_{a2,i}^{*} \left(x'\right) + \sum_{i=1}^{4}\hat{w}_{s,i}\Psi_{i}^{*} \left(x'\right)\end{array}$

Weak Formulations of The Boundary Conditions

Approximation of acoustic pressure and plate displacement do not fulfil the strong form of boundary conditions. Instead, there are going to be deviations from desired pressure or velocity at boundary. To quantify these deviations, four residual equations are defined in the followings, one for each type of boundary conditions.

$\begin{array}{l}\begin{array}{lll} R_{p}(x,y) & = P(x,y) - \bar{p}(x,y) \phantom{\frac{P(x,y)}{\bar{Z}(x,y)}} & (x,y)\in\Omega_{p} \\ R_{v}(x,y) & = v_{n}(x,y) - \bar{v}_{n}(x,y) \phantom{\frac{P(x,y)}{\bar{Z}(x,y)}} & (x,y)\in\Omega_{v} \\ R_{z}(x,y) & = v_{n}(x.y) - \frac{P(x,y)}{\bar{Z}(x,y)} & (x,y)\in\Omega_{z} \\ R_{s}(x,y) & = v_{n}(x,y) - j\omega w(x,y) \phantom{\frac{P(x,y)}{\bar{Z}(x,y)}} & (x,y)\in\Omega_{s} \end{array}\end{array}$

To obtain the best approximation, homogenous solutions are scaled such that:

Residuals in the form of pressure are orthogonal to the span of acoustic velocity homogenous solutions
Residuals in the form of velocity are orthogonal to the span of acoustic pressure homogenous solutions

This orthogonalisation / projection yields the following integral form of weak formulation of boundary conditions.

$\begin{array}{l}\displaystyle \int_{\Omega_{p}} \tilde{v}_{n}^{*}R_{p}\phantom{.}d\Omega + \int_{\Omega_{v}} \tilde{p}^{*}R_{v}\phantom{.}d\Omega + \int_{\Omega_{z}} \tilde{p}^{*}R_{z}\phantom{.}d\Omega + \int_{\Omega_{s}} \tilde{p}^{*}R_{s}\phantom{.}d\Omega = 0\end{array}$

Problem Formulation as System of Linear Equations

To take advantages of MATLAB's capability, all of the equations need to be converted to matrix equations. Subsequently, the problem is reformulated into a system of linear equations. Collecting the homogenous solutions into a vector, the acoustic pressure approximation and pressure test function can be written in the following matrix notation.

$\begin{array}{l}\displaystyle P = P_{q} + \mathbf{\Phi}_{a}^{H}\hat{\mathbf{p}}_{a} \phantom{xx} \text{and} \phantom{xx} \tilde{p}^{*} = \tilde{\mathbf{p}}_{a}^{H}\mathbf{\Phi}_{a}\end{array}$

Subsequently, the acoustic velocity approximation and velocity test function can also be rewritten in the following matrix notation.

$\begin{array}{l}\displaystyle v_{n} = \frac{j}{\rho_{0}\omega} \left(\nabla P_{q}\right) \hat{\mathbf{n}} + \frac{j}{\rho_{0}\omega} \left[ \left( \nabla \mathbf{\Phi}_{a} \right) \hat{\mathbf{n}}\right]^{H}\hat{\mathbf{p}}_{a} \phantom{xx} \text{and} \phantom{xx} \tilde{v}_{n}^{*} = -\frac{j}{\rho_{0}\omega} \tilde{\mathbf{p}}_{a}^{H} \left[ \left( \nabla \mathbf{\Phi}_{a} \right) \hat{\mathbf{n}}\right]\end{array}$

Last, the plate displacement approximation is also rewritten in the following matrix notation.

$\begin{array}{l}\displaystyle w = w_{F} + w_{q} + \mathbf{w}_{a}^{H}\hat{\mathbf{p}}_{a} + \hat{\mathbf{\Psi}}^{H}\hat{\mathbf{w}}_{s}\end{array}$

Substitution into the weak formulation of boundary conditions and plate's end condition, followed by rearrangement of terms and assembly of terms into matrix notation yields the following system of linear equations. The "stiffness" matrix consists of terms that contain product of homogenous solutions. The difference between particular solutions and desired boundary values appear as the "force" term in the system of equations, while the coefficients associated with each homogenous solution appear as the "DOF".

$\begin{array}{l}\begin{bmatrix} \mathbf{K}_{as} & \mathbf{K}_{aa} \\ \mathbf{K}_{ss} & \mathbf{K}_{sa} \\ \end{bmatrix} \begin{pmatrix} \hat{\mathbf{w}}_{s} \\ \hat{\mathbf{p}}_{a} \end{pmatrix} = \begin{pmatrix} \mathbf{F}_{a} \\ \mathbf{F}_{s} \end{pmatrix}\end{array}$

As shown below, particular solutions in general don't fulfil the desired boundary conditions. One can then interpret the system of equations as follow; difference between particular solutions and boundary condition forces us to add homogenous solutions, and the problem becomes finding the right coefficients, $\begin{array}{l}\hat{\mathbf{w}}_{s}}\end{array}$ and $\begin{array}{l}\hat{\mathbf{p}}_{a}}\end{array}$ , to compensate this difference.

Implementation

Structure of Scripts

The project is implemented completely in a MATLAB live script. This approach was chosen because the code is intended for educational purpose. Consequently, readability takes higher priority over generality or performance. The script is written in a structure similar to the report, with each equation followed immediately by its MATLAB implementation. A snippet of the script is shown below.

Matrices Computations and Numerical Integrations

The bulk of the program deals with boundary integrations to compute matrices for the system of linear equations. Through use of three dimensional array, with each of the third indices associated with different quadrature point, it is possible to perform page - wise operations e.g. page matrix multiplication and sum all the elements over the third indices. This is equivalent to looping through all quadrature points and sum the integrands from each point. The computations of the matrices is summarised by the following flowchart.

After solving for the coefficients, the pressure and velocity field, as well as the plate displacement is calculated in a post-processing step using the approximate solution formulas provided in the Design of Numerical Solution section.

Results

This section contain several results from wave-based methods. In the first half of the section, results for a system with all four boundary conditions is shown. Each boundary condition is checked and its convergence rate is examined.Then the performance is compared against finite element method, particularly convergence rate and size of resulting system of equations.

Numerical Solution Example

As an example, consider an acoustic cavity with all four types of boundary conditions present. There is a harmonic point source inside the cavity and the plate is loaded with mechanical point load with the same frequency but with opposite phase. The boundary conditions and loads are shown in the following figure.

Pressure Field Approximation

The resulting pressure field is shown in the following figure.

As can be seen in the figure, the pressure along $\begin{array}{l}y=0\end{array}$ is zero, as prescribed by the boundary condition. On the contrary, the pressure along $\begin{array}{l}y=1\end{array}$ is largely negative. This is because a large outgoing velocity is prescribed by one of the other boundaries, so a large negative pressure is needed to allow air to enter the acoustic cavity from the impedance boundary. The dimple indicates the location of negative acoustic point source.

Velocity Field Approximation

The resulting velocity vector field is shown in the following figure. There is a uniform outward flux on the left side, as prescribed by the velocity boundary condition.

Plate Displacement Approximation

The resulting plate displacement is shown in the following figure. As can be seen from the figure, the displacements and derivatives are zero at the ends, as prescribed by the plate's end conditions.

Boundary Conditions Check

To check the structural boundary conditions, the acoustic normal velocity at the boundary is plotted against the plate velocity below. As can be seen from the figure, the two velocities match each other and satisfy the structural boundary condition.

To check the impedance boundary conditions, the acoustic normal velocity at the boundary is plotted against the acoustic pressure below. As can be seen from the figure, the points lie on the line $\begin{array}{l}y = x / 428\end{array}$ and therefore satisfy the impedance boundary condition.

Convergence for Each Type of Boundary

Residuals from numerical solutions with various number of DOFs are checked. The square of residuals are integrated along the boundaries and plotted in the following figures.

As can be seen from the graphs, all residuals diminish as more and more DOFs are used. One particular observation from the plate boundary condition is that the residual only start to diminish when the number of DOFs reach a certain value. This is because, when using low number of DOFs, all the spatial wavenumber of acoustic homogenous solutions lie below the spatial wavenumber of plate's homogenous solution. Consequently, good approximation of the plate BC requires a certain minimum number of DOFs.

Plate Displacement vs FEM

The steady state displacement of the plate boundary is compared with result from a FE model consisting of linear beam elements. The pressure along the plate calculated with WBM is used as load in the FE model. The problem used in this part is identical to the previous one, but without the mechanical load and with the harmonic load frequency reduced from 50Hz to 3Hz because of available hardware limitation to run the FE simulation. The displacement from both WBM and FEM are shown in the following figure.

As the figure shown, the displacements from both method agree quite well. To show that displacement from FE model converges to the displacement from WBM, the FE simulation is run using various number of element. The minimum displacements are plotted against number of DOFs in the following figure.

As shown in the figure, the minimum displacement from FE model tend to converges to displacement computed by the WBM as the number of DOFs is increased. The mean square error is calculated by taking the FE solution with the highest number of DOF as the reference solution. The mean square error of the WBM model is 0.323%.

Acoustic Pressure vs FEM

In this section, a simple problem is considered. Inside the acoustic cavity, there is a harmonic point source. Pressure boundary conditions are imposed to all four sides of the acoustic cavity. The problem is depicted in the following figure.

For finite element method, a simple bilinear regular mesh is used. The number of DOFs for both wave-based method and finite element methods are varied and the convergence of pressure at two points are checked. The convergence plots are given below.

As we can see from the graphs, finite element method requires higher number of DOFs. The difference is about 4-5 order of magnitude. This is because FEM approximate gradient and derivatives with lower order polynomial, which leads to locking and lower convergence rate. On the other hand, in WBM the span of homogenous solutions and their derivatives are the same.

While higher order elements and anti-locking techniques are available, it is noteworthy that really low frequency is used in the example due to hardware limitations. As frequency increase, significantly more elements are needed, since the element size needs to be in the order of the acoustic wavelength. This emphasises the advantage of wave-based methods over element-based methods in mid to high range of frequencies.

Size of the System of Equations

High number of DOFs results in large systems matrices. However, the stiffness matrix from finite element method is sparse and real, while the matrix from wave-based method is dense and complex. The sparsity patterns for stiffness matrices from both methods are shown below.

The sparsity of matrix from finite element method is because the elements only interact with their respective neighbours. Meanwhile, the homogenous solutions used in wave-based method are not necessarily orthogonal with each other. The magnitude of the matrix elements from wave-based method is depicted by the following figure.

Conclusion

The given MATLAB script's functionalities have been extended by the implementation of impedance and structural boundary conditions. The implementations have been shown to produce solutions that fulfil prescribed boundary conditions well. Pressure, velocity, and impedance boundary conditions converges rapidly, at convergence rate inversely proportional to the number of DOFs raised to the power of three or more. Structural boundary condition also converges at similar rate, provided that there are wave functions with wavenumber higher than plate's homogenous solutions'.

The WBM has been shown to produce smaller system of equations compared to FEM, particularly for harmonic acoustic load frequencies in the audible range. However, WBM produces dense and complex matrix, while FEM produces sparse and real matrix. Consequently, FEM and/or other element-based methods are still preferred over WBM in low frequency range.

At the time this report is written, the script is limited to one acoustic cavity and one plate. In addition, only plate clamped at both ends have been implemented. Therefore, future works may include extension of current script and/or implementation of additional functionalities to address these.

References

Genechten, B. (2006). On the Coupling of Wave Based Models with Modally Reduced Finite Element Models for Structural-Acoustic Analysis. Proceedings of ISMA2006: International Conference on Noise and Vibration Engineering.
Desmet, W. (1998). A Wave Based Prediction Technique for Coupled Vibro-Acoustic Analysis [Dissertation]. Katholieke Universiteit Leuven.
Howard, C. Q., & Cazzolato, B. S. (2015). Acoustic Analyses Using MATLAB® and ANSYS®.
Liu, L. (2021). Simulation and testing of a coupled plate-cavity system targeted for vehicle interior noise analysis and control [Dissertation]. Politecnico di Milano.
Du, J. T., Li, W. L., Xu, H. A., & Liu, Z. G. (2012). Vibro-acoustic analysis of a rectangular cavity bounded by a flexible panel with elastically restrained edges. The Journal of the Acoustical Society of America, 131(4), 2799–2810.

Appendix

Contains more detailed derivations of equations used in the project.

Helmholtz Equation

Using decomposition of density and velocity into their ambient and perturbation parts, the conservation of mass equation is given by the following. Most of the terms are dropped because ambient terms are assumed to be constant with respect to both time and space. Particularly, ambient velocity is assumed to be zero. Furthermore, the perturbation terms are assumed much smaller than ambient terms, so non-linear perturbation terms are also dropped.

$\begin{aligned} 0 &= - \dot{\rho}_{0} - \dot{\rho} ' + \left( \rho_{0} + \rho ' \right) q' - \nabla \cdot \left[ \left( \rho_{0} + \rho ' \right) \left( \mathbf{v}_{0} + \mathbf{v} ' \right) \right] \\ &= - \dot{\rho} ' + \rho_{0} q ' - \rho_{0} \nabla \cdot \mathbf{v} ' - \dot{\rho}_{0} + \rho ' q' - \mathbf{v}_{0} \nabla \cdot \rho ' - \nabla \cdot \left( \rho ' \mathbf{v} ' \right) \\ &= - \dot{\rho} ' + \rho_{0} q ' - \rho_{0} \nabla \cdot \mathbf{v} ' \end{aligned}$

Using decomposition of density, velocity, and pressure into their ambient and perturbation parts, the conservation of linear momentum is given by the following. The equation is simplified in a similar manner with the conservation of mass equation.

$\begin{aligned} 0 &= - \left( \rho_{0} + \rho ' \right) \left[ \dot{\mathbf{v}}_{0} + \dot{\mathbf{v}} ' + \nabla \left( \mathbf{v}_{0} + \mathbf{v} ' \right) \cdot \left( \mathbf{v}_{0} + \mathbf{v} ' \right) \right] - \nabla \left( p_{0} + p' \right) \\ &= - \rho_{0} \dot{\mathbf{v}} ' - \nabla p ' - \rho_{0} \left[ \dot{\mathbf{v}}_{0} + \nabla \left( \mathbf{v}_{0} + \mathbf{v} ' \right) \cdot \left( \mathbf{v}_{0} + \mathbf{v} ' \right) \right] - \rho ' \left[ \dot{\mathbf{v}}_{0} + \dot{\mathbf{v}} ' + \nabla \left( \mathbf{v}_{0} + \mathbf{v} ' \right) \cdot \left( \mathbf{v}_{0} + \mathbf{v} ' \right) \right] - \nabla p_{0} \\ &= - \rho_{0} \dot{\mathbf{v}} ' - \nabla p ' \end{aligned}$

Using adiabatic relation, pressure is only a function of density. Taylor series expansion of pressure around ambient state is given by the following. Using small perturbation assumption, higher order terms are dropped out.

$\begin{aligned} p' &= \rho ' \left. \frac{dp}{d\rho} \right|_{p_{0}.\rho_{0}} + \left( \rho ' \right)^{2} \left. \frac{d^{2}p}{d\rho^{2}} \right|_{p_{0}.\rho_{0}} + ... \\ &= \rho ' \left. \frac{dp}{d\rho} \right|_{p_{0}.\rho_{0}} \\ &= \rho ' \frac{\gamma p_{0}}{\rho_{0}} \end{aligned}$

Mass conservation equation is differentiated with respect to time, followed by substitution of linear momentum conservation equation and the simplified adiabatic relation between pressure and density. Subsequently, Fourier transformation is applied to the result to obtain the Helmholtz equation.

$\begin{aligned} -\rho_{0} \dot{q}' &= - \nabla \cdot \left( \rho \dot{\mathbf{v}}' \right) -\ddot{\rho}' \\ &= \phantom{-}\nabla \cdot \nabla p' -\frac{d \rho '}{d p'} \cdot \ddot{p}' \\ &= \phantom{-}\nabla^{2} p' -\left( \frac{\gamma p_{0}}{\rho_{0}} \right)^{-1} \cdot \ddot{p}' \\ &= \phantom{-}\nabla^{2} p' -\frac{1}{c^{2}} \ddot{p}' \end{aligned}$

$\begin{array}{l}\displaystyle -\rho_{0} \dot{q}' = \nabla^{2} p' -\frac{1}{c^{2}} \ddot{p}' \phantom{x} \xrightarrow{\phantom{x}\mathscr{F}(\phantom{x})\phantom{x}} \phantom{x} -j\rho_{0}\omega Q = \nabla^{2} P + \frac{\omega^{2}}{c^{2}} P\end{array}$

Weak Formulation of The Boundary Conditions

Approximate of acoustic pressure and plate displacement are substituted into each term of the weak formulation of all boundary conditions to generate a system of equations. Subsequently, matrix notation is used to form a concise equation. The result for prescribed pressure boundary condition is shown below.

$\begin{aligned} \int_{\Omega_{p}} \tilde{v}_{n}^{*} R_{p} \phantom{.} d\Omega &= \int_{\Omega_{p}} \tilde{v}_{n}^{*} \left( P - \bar{p} \right) \phantom{.} d\Omega \\ &= - \tilde{\mathbf{p}}_{a}^{H} \frac{j}{\rho_{0}\omega} \int_{\Omega_{p}} \left[ \left( \nabla\mathbf{\Phi}_{a} \right) \hat{\mathbf{n}} \right] \left( P_{q} + \mathbf{\Phi}^{H}\hat{\mathbf{p}}_{a} - \bar{p} \right) \phantom{.} d\Omega \\ &= \tilde{\mathbf{p}}_{a}^{H} \mathbf{K}_{aap} \hat{\mathbf{p}}_{a} - \tilde{\mathbf{p}}_{a}^{H} \mathbf{F}_{ap} \phantom{\int_{\Omega_{p}}} \\ \\ \mathbf{K}_{aap} &= \frac{-j}{\rho_{0}\omega} \int_{\Omega_{p}} \left[ \left( \nabla\mathbf{\Phi}_{a} \right) \hat{\mathbf{n}} \right] \mathbf{\Phi}^{H} \phantom{.} d\Omega \\ \mathbf{F}_{ap} &= \frac{j}{\rho_{0}\omega} \int_{\Omega_{p}} \left[ \left( \nabla\mathbf{\Phi}_{a} \right) \hat{\mathbf{n}} \right] \left( P_{q} - \bar{p} \right) \phantom{.} d\Omega \end{aligned}$

The result for prescribed velocity boundary condition is shown below.

$\begin{aligned} \int_{\Omega_{v}} \tilde{p}^{*} R_{v} \phantom{.} d\Omega &= \int_{\Omega_{v}} \tilde{p}^{*} \left( v_{n} - \bar{v}_{n} \right) \phantom{.} d\Omega \\ &= \tilde{\mathbf{p}}_{a}^{H} \frac{j}{\rho_{0}\omega} \int_{\Omega_{v}} \mathbf{\Phi}_{a} \left( \left( \nabla P_{q} \right) \hat{\mathbf{n}} + \left[ \left( \nabla \mathbf{\Phi}_{a} \right) \hat{\mathbf{n}} \right]^{H}\hat{\mathbf{p}}_{a} - \bar{v}_{n} \right) \phantom{.} d\Omega \\ &= \tilde{\mathbf{p}}_{a}^{H} \mathbf{K}_{aav} \hat{\mathbf{p}}_{a} - \tilde{\mathbf{p}}_{a}^{H} \mathbf{F}_{av} \phantom{\int_{\Omega_{v}}} \\ \\ \mathbf{K}_{aav} &= \frac{j}{\rho_{0}\omega} \int_{\Omega_{v}} \mathbf{\Phi}_{a}\left[ \left( \nabla \mathbf{\Phi}_{a} \right) \hat{\mathbf{n}} \right]^{H} \phantom{.} d\Omega \\ \mathbf{F}_{av} &= \frac{j}{\rho_{0}\omega} \int_{\Omega_{v}} \mathbf{\Phi}_{a} \left( \bar{v}_{n} - \left( \nabla P_{q} \right) \hat{\mathbf{n}} \right) \phantom{.} d\Omega \end{aligned}$

The result for impedance boundary condition is shown below.

$\begin{aligned} \int_{\Omega_{z}} \tilde{p}^{*} R_{z} \phantom{.} d\Omega &= \int_{\Omega_{z}} \tilde{p}^{*} \left( v_{n} - \frac{1}{Z}P \right) \phantom{.} d\Omega \\ &= \tilde{\mathbf{p}}_{a}^{H} \int_{\Omega_{z}} \mathbf{\Phi}_{a} \left( \frac{j}{\rho_{0}\omega} \left( \nabla P_{q} \right)\hat{\mathbf{n}} + \frac{j}{\rho_{0}\omega} \left[ \left( \nabla\mathbf{\Phi}_{a} \right)\hat{\mathbf{n}} \right]^{H}\hat{\mathbf{p}}_{a} - \frac{1}{Z} \left( P_{q} - \mathbf{\Phi}_{a}^{H} \hat{\mathbf{p}}_{a} \right) \right) \phantom{.} d\Omega \\ &= \tilde{\mathbf{p}}_{a}^{H} \mathbf{K}_{aaz} \hat{\mathbf{p}}_{a} - \tilde{\mathbf{p}}_{a}^{H} \mathbf{F}_{az} \phantom{\int_{\Omega_{z}}} \\ \\ \mathbf{K}_{aaz} &= \int_{\Omega_{z}} \frac{j}{\rho_{0}\omega} \mathbf{\Phi}_{a} \left[ \left( \nabla\mathbf{\Phi}_{a} \right)\hat{\mathbf{n}} \right]^{H} - \mathbf{\Phi}_{a} \mathbf{\Phi}_{a}^{H} \phantom{.} d\Omega \\ \mathbf{F}_{az} &= \int_{\Omega_{z}} \frac{1}{Z} P_{q} - \frac{j}{\rho_{0}\omega} \left( \nabla P_{q} \right)\hat{\mathbf{n}} \phantom{.} d\Omega \end{aligned}$

The result for structural boundary condition is shown below.

$\begin{aligned} \int_{\Omega_{s}} \tilde{p}^{*} R_{s} \phantom{.} d\Omega &= \int_{\Omega_{s}} \tilde{p}^{*} \left( v_{n} - j\omega w \right) \phantom{.} d\Omega \\ &= \tilde{\mathbf{p}}_{a}^{H} \int_{\Omega_{s}} \mathbf{\Phi}_{a} \left( \frac{j}{\rho_{0}\omega} \left(\nabla \mathbf{P}_{q}\right)\hat{\mathbf{n}} + \frac{j}{\rho_{0}\omega}\left[ \left( \nabla \mathbf{\Phi}_{a} \right)\hat{\mathbf{n}} \right]^{H}\hat{\mathbf{p}}_{a} - j\omega \left[ w_{F} + w_{q} + \mathbf{w}_{a}^{H}\hat{\mathbf{p}}_{a} + \hat{\mathbf{\Psi}}^{H} \hat{\mathbf{w}}_{s} \right] \right) \phantom{.} d\Omega \\ &= \tilde{\mathbf{p}}_{a}^{H} \mathbf{K}_{aas} \hat{\mathbf{p}}_{a} + \tilde{\mathbf{p}}_{a}^{H} \mathbf{K}_{as} \hat{\mathbf{w}}_{s} - \tilde{\mathbf{p}}_{a}^{H} \mathbf{F}_{as} \phantom{\int_{\Omega_{s}}} \\ \\ \mathbf{K}_{aas} &= \int_{\Omega_{s}} \frac{j}{\rho_{0}\omega} \mathbf{\Phi}_{a} \left[ \left( \nabla \mathbf{\Phi}_{a} \right)\hat{\mathbf{n}} \right]^{H} - j\omega \mathbf{\Phi}_{a} \mathbf{w}_{a}^{H} \phantom{.} d\Omega \\ \mathbf{K}_{as} &= -j\omega \int_{\Omega_{s}} \mathbf{\Phi}_{a} \hat{\mathbf{\Psi}}^{H} \phantom{.} d\Omega \\ \mathbf{F}_{as} &= \int_{\Omega_{s}} \mathbf{\Phi}_{a} \left( j\omega w_{F} + j\omega w_{q} - \frac{j}{\rho_{0}\omega} \left( \nabla P_{q} \right) \hat{\mathbf{n}} \right) \phantom{.} d\Omega \end{aligned}$

The results are all summed to form the weak formulation of the boundary condition below. The term $\begin{array}{l}\tilde{\mathbf{p}}_{a}^{H}\end{array}$ is arbitrary and thus can be dropped.

$\begin{aligned} \mathbf{F}_{a} &= \begin{bmatrix} \mathbf{K}_{aa} & \mathbf{K}_{as} \end{bmatrix} \begin{pmatrix} \hat{\mathbf{p}}_{a} \\ \hat{\mathbf{w}}_{s} \end{pmatrix} \\ \\ \mathbf{K}_{aa} &= \mathbf{K}_{aap} + \mathbf{K}_{aav} + \mathbf{K}_{aaz} + \mathbf{K}_{aas} \\ \mathbf{F}_{a} &= \mathbf{F}_{ap} + \mathbf{F}_{av} + \mathbf{F}_{az} + \mathbf{F}_{as} \end{aligned}$

The system of equations is underdetermined because it is missing the plate's boundary conditions. The plate's boundary condition is given below.

$\begin{array}{l}\begin{array}{mlmlmlml} 0 =& w_{F}(0) &+& w_{q}(0) &+& \mathbf{w}_{a}^{H}(0)\hat{\mathbf{p}}_{a} &+& \hat{\mathbf{\Psi}}^{H}(0)\hat{\mathbf{w}} \\ 0 =& w_{F}(L) &+& w_{q}(L) &+& \mathbf{w}_{a}^{H}(L)\hat{\mathbf{p}}_{a} &+& \hat{\mathbf{\Psi}}^{H}(L)\hat{\mathbf{w}} \\ 0 =& \frac{d}{dx'}w_{F}(0) &+& \frac{d}{dx'}w_{q}(0) &+& \frac{d}{dx'}\mathbf{w}_{a}^{H}(0)\hat{\mathbf{p}}_{a} &+& \frac{d}{dx'}\hat{\mathbf{\Psi}}^{H}(0)\hat{\mathbf{w}} \\ 0 =& \frac{d}w_{F}(L){dx'} &+& \frac{d}{dx'}w_{q}(L) &+& \frac{d}{dx'}\mathbf{w}_{a}^{H}(L)\hat{\mathbf{p}}_{a} &+& \frac{d}{dx'}\hat{\mathbf{\Psi}}^{H}(L)\hat{\mathbf{w}} \end{array}\end{array}$

Written in matrix notation

$\begin{array}{l}\begin{bmatrix} \mathbf{K}_{ss} & \mathbf{K}_{sa} \end{bmatrix} \begin{pmatrix} \hat{\mathbf{w}} \\ \hat{\mathbf{p}}_{a} \end{pmatrix} = \mathbf{F}_{s} \\ \\ \begin{aligned} \mathbf{K}_{ss} =& \begin{bmatrix} \mathbf{\Psi}(0) & \mathbf{\Psi}(L) & \frac{d}{dx'}\mathbf{\Psi}(0) & \frac{d}{dx'}\mathbf{\Psi}(L) \end{bmatrix}^{H} \\ \mathbf{K}_{sa} =& \begin{bmatrix} \mathbf{w}_{a}(0) & \mathbf{w}_{a}(L) & \frac{d}{dx'}\mathbf{w}_{a}(0) & \frac{d}{dx'}\mathbf{w}_{a}(L) \end{bmatrix}^{H} \\ \\ \mathbf{F}_{s} =& \begin{pmatrix} -w_{F}(0) - w_{q}(0) \\ -w_{F}(L) - w_{q}(L) \\ -\frac{d}{dx'}w_{F}(0) - \frac{d}{dx'}w_{q}(0) \\ -\frac{d}{dx'}w_{F}(L) - \frac{d}{dx'}w_{q}(L) \end{pmatrix} \end{aligned}\end{array}$

The two equations are combined to form the system of equations to be solved.

$\begin{array}{l}\begin{bmatrix} \mathbf{K}_{as} & \mathbf{K}_{aa} \\ \mathbf{K}_{ss} & \mathbf{K}_{sa} \end{bmatrix} \begin{pmatrix} \hat{\mathbf{w}} \\ \hat{\mathbf{p}}_{a} \end{pmatrix} = \begin{pmatrix} \mathbf{F}_{a} \\ \mathbf{F}_{s} \end{pmatrix}\end{array}$

Seitenhierarchie

Table of Contents

Theoretical Basis

Helmholtz Equation

Plate Equation

Boundary Conditions

Design of Numerical Solution

Approximate Solution to Acoustic Pressure

Approximate Solution to Plate Displacement

Weak Formulations of The Boundary Conditions

Problem Formulation as System of Linear Equations

Implementation

Structure of Scripts

Matrices Computations and Numerical Integrations

Results

Numerical Solution Example

Pressure Field Approximation

Velocity Field Approximation

Plate Displacement Approximation

Boundary Conditions Check

Convergence for Each Type of Boundary

Plate Displacement vs FEM

Acoustic Pressure vs FEM

Size of the System of Equations

Conclusion

References

Appendix

Helmholtz Equation

Weak Formulation of The Boundary Conditions