# Finite element analysis of 2D acoustics

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:

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.

# Listening to vibrations

Update: here’s a small tune I made to demonstrate what the resulting sample would sound like as an instrument:

I’ve always been fascinated by how things, especially sound waves, look like in slow motion. By examining something in slow motion, somehow even complex phenomena start to seem intelligible.

In this post, I’m examining the behavior of a very simple instrument in slow motion. The physics are very much the same as when examining vibrations in structures, but instruments are (very much) nicer to listen to. The principles shown here are naturally applicable to larger scale phenomena that occur in structures of different kind (buildings, bridges, etc).

## The instrument

The instrument will consist of a cantilever beam made out of steel, with a length of 20 cm. The cross-sectional dimensions of the beam are 2 cm x 2 cm. The instrument will thus resemble something like a simple vibraphone / glockenspiel / large music box.

We’ll use the same principle to get the sound of the beam as an electric guitar uses to pick up the sound of a string, by placing a (virtual) pickup at a position 3/4:th to the left of the tip of the beam (or 1/4:th from the left of the beam, whichever suits your fancy).

What will happen if we hit the very tip of the beam, with a very sharp, impulse-like force (even sharper and quicker than what’s possible with the hammer in the picture)?

## The result

First, it should be noted that the same principle is used as when an electric guitar picks up the sound of a string. The pickup transforms the deflection at a certain point on the string (or in this case, beam) directly to sound. The resulting sound for the beam is as following:

To see what the beam looks like directly after being hit by the sharp impulse, examine the resulting deformation of the beam as a function of time here (click on the blue area to load the content). The sound is generated from the movement of the circle in the direction of the y-axis.

http://kaistale.com/blog/130722beam1d/animation.html

Some things to notice when watching  the time lapse of the deformation:

• There’s a transient part at the very beginning (slightly visible in the 0.0015 s time span) which attenuates very quickly. This is the part of the response where standing waves (resonances) have not yet formed.
• The two lowest natural frequencies (resonances) can be seen clearly; one at 464 Hz and another at 2910 Hz. The third natural frequency, 8150 Hz, can be seen at the very beginning of the response.
• For this setup, the higher frequencies attenuate quickly. In the end, only the lowest natural frequency, 464 Hz, can be heard. This gives the distinct pitch you can hear in the sound sample above.

## The theory

I’m using finite element analysis to examine the behavior of the beam. The theory is the same as in the previous post, but I’ve also calculated the mass matrix for the beam. I’ve used 25 elements for the beam, thus solving a 50-degree of freedom system. Note that I’m using the Euler-Bernoulli beam theory, so some simplifications are made.

The damping was done using Rayleigh damping. For those familiar with this type of damping, the value for $$\alpha$$ was 1e-05 and the value for $$\beta$$ was 1.5e-06. I chose these values simply on the basis of listening to which damping values sounded better than others, not much thought were put into them. They seem surprisingly small as compared to other values I encountered online, maybe someone more acquainted with Rayleigh damping could offer me their opinion on this?

I used the Newmark algorithm for time-stepping, with a value of 0.5 for both $$\beta_1$$ and $$\beta_2$$.

## The methods

I used Python for the calculations. Python is ideal for such an endeavor, and free! SymPy provided me with the tools I needed to solve the necessary equations, while NumPy did the calculations for me.

I saved the resulting calculations as a binary file, which only contained the necessary information. For example: only one byte / element describes the deformation at each time step in the animation seen above, as 8 bits is more than enough in this case.

I wrote a script that loads up these binary files (of which 3 can be seen above). The deformations can be examined as a function of time. I’m quite happy about how it turned out, even though some parts of the code are quite hacky-ish. You can see the code for the javascript here. Note that this code only displays the data from the binary files, the binary files have been calculated using Python.

## Thoughts for the future

Perhaps the parameters could be tweaked a little, so the sound would be closer to the sound from a music box (which is the instrument I think the model most closely resembles). Also, the Timoshenko beam model would be nice to try, for comparison. All in all, this post was a nice exercise in programming Python, JavaScript and doing finite element analysis in the time domain.