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

# Mass-spring systems should be understood before they are used

Mass-spring systems are often used to model the behavior of an isolated vibrating system.

The mass-spring viewpoint should not be utilized in cases where the system can not be reliably modeled as a mass-spring system. The importance of the static compression of the isolating material, due to gravity, is also overemphasized and often misunderstood.

Some important points:

• The mass isolates vibrations by resisting change of momentum. Usually: the more mass, the better. The weight (due to gravity) has very little to do with the isolation, it’s the mass that matters!
• The spring’s main purpose is to carry static forces (for example the weight). It will also try to resist the displacements of the mass, which means that it will transfer forces to the foundations. From the perspective of isolating vibrations, the spring should be as loose as possible, to allow for the mass to vibrate freely.
• The mass-spring combination only works on frequencies higher than the resonance frequency.
• Ideally, there would be no “spring”, but it’s impossible in practice.

# Simplifying rotating parts

Let’s start from the very basic reasoning behind why the mass-spring perspective is so common when designing vibration isolation. Consider case a. We have a machine, with a part rotating around some axis.

Case b shows the varying forces such a movement creates. Perpendicular forces exist relative to the x-axis and y-axis. The varying force in the direction of the x-axis, together with the resisting force at the support, also create a moment of force (kind of like attempting to tip the machine over). The vertical component of the force often creates a moment of force, as well. As does accelerating or decelerating angular motion.

In practical cases, it has been shown that it’s often enough to consider the sinusoidal force shown in case c (but it’s not always enough!).

Note that this model assumes both the machine and the support to be perfectly rigid. Note also that the force required to keep the rotating mass in bay will rise along with the frequency ($$\omega^2$$)! I think I will return to this in another post.

# The equilibrium of the system

So far we have simplified the system to a sinusoidally varying vertical force, acting on a mass on a spring.

We will describe the system using these parameters:

• The mass of the machine
• A varying force directed on the machine
• The varying displacement of the machine
• The spring constant, which we’ll assume to be constant (linearly elastic)

We now have five different sources for forces:

1. $$F_{m} = ma = m\ddot y$$, caused by the mass/momentum of the machine which resists acceleration
2. $$F_{k} = k y$$, caused by the compression of the spring
3. $$F_{c} = cv = c\dot y$$, caused by various energy losses which are directly proportional to the speed of the displacement
4. $$F_{e}$$, external forces which we feed into the system.
5. $$F_{g}$$, the force of gravity pulling down the mass

Using the principle of equilibrium (Newton), the sum of these forces must be 0.

$$F_k + F_c + F_m + F_e + F_g = 0$$

$$F_g$$ can be discarded in our calculations. I’ll soon explain why.

# Calculating the natural frequency of the system

We know that such a system oscillates naturally at a specific frequency. To calculate this frequency, we assume three things:

1. The motion is sinusoidal
2. There are no external forces ($$F_e$$)
3. The damping factor ($$F_c$$) approaches zero, allowing for the system to oscillate freely

By assuming the motion follows the function $$y(t) = A\sin(\omega t)$$, we get the following equation:

$$F_k + F_m + F_g = k(y(t) + y_0) + m\ddot y(t) + mg = k A\sin(\omega t) – m A \omega^2\sin(\omega t) = 0$$,

where $$y_0$$ is the compression caused by the pull of gravity.

$$k y_0$$ has to be equal to $$m g$$ for the system to balance out (Hooke’s law), which means that we can throw the effect of gravity out of the equation. Intuitively this can be explained as following: as the effects of gravity and the static compression of the string ($$k y_0$$) constantly cancel each other out, only the vibrating part of the displacement ($$y(t)$$) makes any difference. The force caused by the vibrating part of the displacement is still directly proportional to the spring constant. Gravity (and as such the static compression of the spring) makes no difference. Gravity only affects the situation if the spring constant changes because of it, not because it compresses the spring or adds a constant downward force to the mass. This holds true for all constant forces, so constant forces should never be included in a vibrating mass-spring system to avoid confusion! They should only be used to calculate the correct spring stiffness, which is almost always very nearly linear when the vibrations are small.

The forces will balance out when $$\omega = \sqrt{\frac{k}{m}}$$, or $$f_0 = \frac{1}{2 \pi}\sqrt{\frac{k}{m}}$$.

Note that this is the combined natural frequency of the spring and the mass! This might seem obvious now, but when we replace the mass and spring constant with something else, it might not be that obvious. The compression of the spring, using Hooke’s law, will be $$mg = k\delta x$$. Inserting this into the formula will give $$f_0 = \frac{1}{2 \pi}\sqrt{\frac{g}{\delta x}}$$. This formula is often used for practical purposes, but effectively hides the real parameters the natural frequency is based on!

# Calculating the forced response of an undamped system

We will consider the steady state response, meaning the response of the system after being excited at a specific frequency for a long time (when everything balances out). We will assume the mass will move at the same frequency as the force exciting it.

The equilibrium formula shown earlier will now look like this:

$$m \ddot y(t) + k y(t) = F(t)$$

Assuming the machine will be excited according to the function $$F(t) = F\sin(\omega t)$$ (ignoring everything else for now) and the machine moving according to the function $$y(t) = A\sin(\omega t)$$ we get the following formula:

$$-m A \omega^2 sin(\omega t) + k A sin(\omega t) = F_0\sin(\omega t)$$

Solving for A (the amplitude of the motion caused by the force), we get the following formula:

$$A = \frac{F_0}{-m \omega^2 + k}$$

Note that this describes the displacement of the mass. I think the following concept is really important to understand; the force transmitted into the foundations is directly proportional to the displacement. As the spring constant approaches zero (let’s imagine the mass floats in space, for example), the system will still vibrate with an amplitude of $$\frac{F_0}{m \omega^2}$$. Increasing the mass will decrease the amplitude of the vibration, also when there is no spring.

To get the force transmitted into the foundations of the machine, we need to multiply the displacement with the spring constant $$k$$ (according to Hooke’s law). From this, we get the relation of the transmitted force to the original force $$\frac{F_{tr}}{F_0}$$:

$$\frac{F_{tr}}{F_0} = \bigg\vert\frac{1}{\frac{-m}{k}\omega^2 + 1}\bigg\vert$$

The formula nicely shows us a few things:

• The resonant frequency is at $$\omega = \sqrt{\frac{k}{m}}$$, which is the same as earlier
• The effect of decreasing the stiffness of the spring
• The effect of increasing the mass of the machine

## Thoughts

Please note that this has been a significant simplification of the situation. Ideally, the simplification only works for cases where the vibrating force is purely vertical, and acts in the exact direction of the centre of mass of the object.