*This is a purely mathematical post covering some of the theoretical basis behind the matrices used in finite element analysis for acoustics in fluids. In coming posts, I’ll cover some of the more practical uses for the method.*

A lot of the content presented here is based on the book “Computational Acoustics of Noise Propagation in Fluids – Finite and Boundary Element Methods”. I’ll do the examination in 2D, to make things as simple as possible.

## The Helmholtz equation

The linearized Euler equation is the most basic theoretical basis I will consider, even though even this equation can be derived using continuum mechanics. The linearized Euler equation is as following:

$$\Delta\tilde{p}(x,t)=\frac{1}{c^2}\frac{\partial^2 \tilde{p}(x,t)}{\partial t^2}$$

where $$\Delta$$ is the Laplacian and $$\tilde{p}$$ refers to the local perturbation of the pressure (i.e. small disturbance around the equilibrium that causes sound). $$x$$ is a position vector. The pressure is dependent on both position and time. If the problem is assumed to be time harmonic, i.e. the time variation is sinusoidal, $$\tilde{p}(x,t)$$ can be replaced by $$p(\pmb{x})e^{-i\omega t}$$. Placing this into the Euler equation gives the Helmholtz-equation:

$$\Delta p(x) + k^2p(x)=0$$

where $$k$$, the wave number, equals $$\omega/c_0$$. Note that the speed of sound can in some cases vary inside the medium. In these cases, where diffraction is to be considered, an additional factor $$n(x) = c_0/c(x)$$ (*the index of refraction*) should be used to multiply the term on the right.

## The boundary conditions

The boundary conditions can be divided into the following types of conditions:

- Dirichlet boundary conditions, often called a
*soft boundary*, which in acoustics means that the pressure at the boundary is specified - Neumann boundary conditions, often called a
*hard boundary*, which in acoustics means that the mean particle velocity normal to the boundary is specified - Robin boundary conditions, which allows for specification of the boundary
*admittance*.

The Robin boundary condition allows for the specification of both hard and soft boundaries, and it can be defined on the boundary of the problem as following:

$$v_f(x)-v_s(\pmb{x})=Y(\pmb{x})p(\pmb{x})$$

where $$Y$$ represents the boundary admittance. $$v_s$$ denotes the normal velocity of the boundary and $$v_f$$ denotes the normal fluid velocity of the acoustic domain. What exactly does this state? It states that the difference between the speed of the fluid and the speed of the structure equals the impedance times the pressure at the boundary. $$v_f$$ can be calculated from the normal derivative of the sound pressure as following:

$$v_f(\pmb{x}) = \frac{1}{i\omega\rho_0}\frac{\partial p(\pmb{x})}{\partial n(\pmb{x})}$$

## The weak formulation

The weak formulation is calculated in precisely the same way as was the case with the Euler-Bernoulli beam equation, in a previous blog post. In this case, the weak formulation is calculated from the following equation (the $$(\pmb{x})$$-parts are omitted for clarity), where $$w$$ represents the *weighting function*:

$$\int_\Omega w[\Delta p+k^2 p]d\Omega=0$$

After integrating by parts (the first term, the one with the Laplace operator):

$$\int_\Gamma w\nabla p \pmb{n}d\Gamma-\int_\Omega[\nabla p \nabla w-k^2wp]d\Omega=0$$

The first term can now be rewritten:

$$\int_\Gamma w \frac{\partial p}{\partial n}d\Gamma-\int_\Omega[\nabla p \nabla w-k^2wp]d\Omega=0$$

After inserting the earlier equation for $$v_f$$:

$$\int_\Gamma wi\omega\rho_0v_f d\Gamma-\int_\Omega[\nabla p \nabla w-k^2wp]d\Omega=0$$

Finally, let’s insert the Robin boundary condition defined earlier (replacing $$v_f$$ with $$Yp + v_s$$):

$$\int_\Gamma wi\omega\rho_0[Yp + v_s] d\Gamma-\int_\Omega[\nabla p \nabla w-k^2wp]d\Omega=0$$

Substituting $$\omega = kc$$, we get:

$$ikc\rho_0\int_\Gamma wYp d\Gamma+ikc\rho_0\int_\Gamma wv_s d\Gamma-\int_\Omega[\nabla w \nabla p-k^2wp]d\Omega=0$$

## The shape functions

I’ll choose Lagrange quadratic shape functions for this endeavor.

$$N_1$$ | $$N_2$$ | $$N_3$$ |
---|---|---|

$$\frac{(\xi – 0)(\xi-1)}{(-1-0)(-1-1)}$$ | $$\frac{(\xi+1)(\xi-1)}{(0+1)(0 – 1)}$$ | $$\frac{(\xi+1)(\xi-0)}{(1+1)(1 – 0)}$$ |

As these functions are one-dimensional, we need to multiply them by each other in order to get two-dimensional shape functions. To get a 2D shape function with the value 1 at some specific node $$(i,j)$$, we need to do the following multiplication:

$$\phi_i(\xi,\eta)=N_i(\xi)\cdot N_j(\eta)$$

with $$\xi$$ and $$\eta$$ denoting local coordinates inside the element.

We do the same multiplication for all the nodes 1 … 9. Thus we have 9 shape functions describing the evenly spaced values in a 2D square. The node numbering doesn’t really matter, as long as it’s something that is agreed upon.

Note that we can now get the continuous values of the pressure inside the 2D square using the shape functions as following:

$$\sum_{i=1}^N\phi_i(\pmb{x})p_i=\pmb{\phi^T}(\pmb{x})\pmb{p}$$

## Integrating the shape functions

Later on, it will become apparent that we need to integrate the 2D quadratic shape functions. There will be two cases of integration:

- Multiples of shape functions across the fluid domain
- Multiples of shape functions on the boundary of the problem

The derivatives and integrals are, as always, performed in global coordinates. As the shape functions are defined in local coordinates, we need a way of transforming a derivation/integration performed in one coordinate system to the other. It turns out that this can be done using the *jacobian matrix*. It’s basically the same thing as the chain rule, but in matrix form. The jacobian matrix $$\pmb{J}$$ in our case is defined as following:

$$\begin{Bmatrix}\frac{\partial N_a}{\partial \xi}\\\frac{\partial N_a}{\partial \eta}\end{Bmatrix}=\begin{bmatrix}\frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi}\\\frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta}\end{bmatrix}\begin{Bmatrix}\frac{\partial N_a}{\partial x}\\\frac{\partial N_a}{\partial y}\end{Bmatrix}=\pmb{J}\begin{Bmatrix}\frac{\partial N_a}{\partial x}\\\frac{\partial N_a}{\partial y}\end{Bmatrix}$$

To convert the derivative the other way around, we invert the jacobian matrix, like this:

$$\begin{Bmatrix}\frac{\partial N_a}{\partial x}\\\frac{\partial N_a}{\partial y}\end{Bmatrix}=\pmb{J}^{-1}\begin{Bmatrix}\frac{\partial N_a}{\partial \xi}\\\frac{\partial N_a}{\partial \eta}\end{Bmatrix}$$

For integration, the following holds true:

$$dx\,dy=\begin{vmatrix}\pmb{J}\end{vmatrix}d\xi\,d\eta$$

In this case, I won’t write up all the possible integrals. Suffice to say, all the necessary integrals and derivatives can be calculated using these equations.

## Discretization

The continuous problem can now be discretized using the shape functions. I’ll use the following notation for the different discretizations, adopted from the book I mentioned in the beginning of the post:

$$p(\pmb{x})=\pmb{\phi^T}(\pmb{x})\pmb{p}$$

$$v_s(\pmb{x})=\pmb{\overline{\phi}^T}(\pmb{x})\pmb{v}$$

$$Y(\pmb{x})=\pmb{\tilde{\phi}^T}(\pmb{x})\pmb{Y}$$

Inserting the shape functions into the weak formulation (again, omitting the $$(\pmb{x})$$-parts), we now get:

$$ikc\rho_0\int_\Gamma w\pmb{\tilde{\phi}^T}\pmb{Y}\pmb{\phi^T}\pmb{p} d\Gamma+ikc\rho_0\int_\Gamma w\pmb{\overline{\phi}^T}\pmb{v} d\Gamma-\int_\Omega[\nabla w \nabla \pmb{\phi^T}\pmb{p}-k^2w\pmb{\phi^T}\pmb{p}]d\Omega=0$$

The weighting functions are once again arbitrary. We’ll use the Galerkin discretization, which here means that the weighting functions $$w$$ are also replaced with 2D Lagrange quadratic shape functions, $$\phi_i$$. Thus we get N equations (I tried to explain this more thoroughly in this post, which covers a more simple case of finite element analysis). The equations A = 1 … N can be written as following:

$$ikc\rho_0\int_\Gamma \phi_A[\sum_{k=1}^N\tilde{\phi}_kY_k][\sum_{m=1}^N\phi_mp_m]d\Gamma+ikc\rho_0\int_\Gamma \phi_A[\sum_{j=1}^N\overline{\phi}_jv_j]d\Gamma-\int_\Omega[\nabla\phi_A[\sum_{j=1}^N\nabla\phi_jp_j]-k^2\phi_A[\sum_{j=1}^N\phi_jp_j]]d\Omega=0$$

Okay, it’s starting to get confusing. But now we can finally start to simplify things. Let’s write the equation in the following matrix form:

$$(\pmb{K}-ik\pmb{C}-k^2\pmb{M})\pmb{p}=ikc\rho_0\pmb{\Theta}\pmb{v_s}$$

with the following matrices. For the purposes of this post, the impedance will be assumed to be constant on the boundary, which makes the damping matrix simpler to calculate. We’ll also assume that Lagrange quadratic shape functions are used consistently.

$$\pmb{K} = \int_\Omega\nabla\pmb{\phi}\nabla\pmb{\phi}^Td\Omega$$

$$\pmb{C} = \rho_0cY\int_\Gamma\pmb{\phi}\pmb{\phi}^Td\Gamma$$

$$\pmb{M} =\int_\Omega\pmb{\phi}\pmb{\phi}^Td\Omega$$

$$\pmb{\Theta} =\int_\Gamma\pmb{\phi}\pmb{\phi}^Td\Gamma$$

## Conclusion

This post by no means went through everything one should know to be able to do finite element analysis of acoustics in 2D. Still, it covered the theoretical basis behind what I will cover in the following posts.