*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.

john hudsonHi Kai

I found your interesting blog while searching for 2D acoustic finite elements. I am starting to model domestic room acoustics.

Just a thought on 2D interpolation. The Lagrange interpolation functions on a rectangular grid you mention, eg. N1 are not particularly “flat”. For example, if I set up a 2D grid N1(x)*N1(y) at all the integer axis point and set all their amplitudes to unity then the values between the integer axis points are not unity. At worst at a point such as (x=0.5, y=0.5) has value 4*N1(0.5) ^2 = 9/4 or 7.2 dB too high and very different from unity. So if I solve the FEM for the pressure values at the integer 2D axis points, I don’t expect an accurate interpolation result in between the integer points. Have you any comments on this.

regards

John hudson

KaiPost authorHi! Thanks for your comment, let me know if I misunderstood something. I’ll go through some initial stuff first, so we can be sure we’re talking about the same thing.

When you calculate the interpolated values, you need to sum up

allthe shape functions at the given point. N1*N1 is the shape function for just a single node, there will be 9 shape functions (N1*N2, N1*N3, N2*N2, etc) which each form a seperate 2-dimensional shape function, with value 1 at one of the nodes and 0 at the others. One such shape function is shown in the image “two-dimensional lagrange shape function” in the post.You’ll get the value like this: $$\sum_{i=1}^9N_i$$, where N_i is all of the possible 2D shape functions. So, assuming the value is 1 for all the nodes (you’ll want a single value everywhere), you’ll sum up all the 2D shape functions times 1.

With the quadratic Lagrange functions mentioned in the post, you’ll get the single value 1 when summing up all the shape functions. $$\xi$$ and $$\eta$$ will cancel each other out. Which means that the interpolation will work perfectly, you’ll get the value 1 everywhere!

PS, you can try it yourself (this mess is the sum of the shape functions times 1):

$$\xi*\eta*(\xi – 1)*(\eta – 1)/4 + \xi*\eta*(\xi – 1)*(\eta + 1)/4 + \xi*\eta*(\xi + 1)*(\eta – 1)/4 + \xi*\eta*(\xi + 1)*(\eta + 1)/4 – \xi*(\xi – 1)*(\eta – 1)*(\eta + 1)/2 – \xi*(\xi + 1)*(\eta – 1)*(\eta + 1)/2 – \eta*(\xi – 1)*(\xi + 1)*(\eta – 1)/2 – \eta*(\xi – 1)*(\xi + 1)*(\eta + 1)/2 + (\xi – 1)*(\xi + 1)*(\eta – 1)*(\eta + 1)$$ will simplify to the value one.

john hudsonthanks for your rapid reply !. My mistake was to truncate the Lagrange function for x,y1. If they are allowed to have infinite extent then the problem disappears as you say. I was starting from the triangle (1st order) interpolation 1-|x|, -1<x<1, which works in 1D, but has an interpolation problem in 2D. An alternative was using 1st order hexagonal shape functions in 2D but these don't fit well into orthogonal coordinates x,y for the wave equation.