Category Archives: Math

Edge diffraction with geometrical acoustics

In this post, I'm using geometrical acoustics to calculate the sound field around a wedge, including specular reflections, shadow zones and diffraction.

The simulation

You might encounter some problems with the simulation if you're not using recent hardware/software with full WebGL support.

The simulation:

Note that when you move the source or change the frequency of the source, it will take a while for the diffraction calculations to update.

I recommend that you keep the source to the left of the wedge if you want to observe the results of the diffraction calculations, as the calculations are done starting to the right of the wedge (you'll see what I mean if you play around with the simulation).

I could make the code much more effective by optimizing things (many times more effective, probably), but I won't do any of that as it's just a proof of concept.

The theory

The specular reflections are calculated by reflecting the source relative to the wedge. The shadowing is self-explanatory.

The diffraction is calculated using the Biot-Tolstoy expressions, using the method presented by Svensson et al (JASA 1999). The simulation spans 1 meter by 1 meters. The wedge spans 1 meter to both sides of the source. The simulation represents a cut plane perpendicular to to the wedge.

The amplitude and phase of the diffracted signal is calculated from the impulse response at many points in space, using the Goertzel algorithm (for a single frequency). The algorithm is implemented server-side, using Python/Numpy.

The calculations are done with some simplifications, to keep the server happy. I haven't validated the results in the simulation thoroughly, small errors could be seen for the few points I tested (most likely due to the simplifications). But the calculation model works; I compared my Python implementation more rigorously with the examples in the paper by Svensson.

The implementation

A lot of things are calculated using the shaders which definitely should not be calculated using them. If you delve deeper into the code (esp. the shaders), you will most definitely encounter quite a few ugly things. But the simulation runs happily on my computer(s), so I'm happy.

Each time you move the source or change the frequency, the contribution of the diffracted signal is calculated for a polar grid, spanning 32 x 32 points (it extends quite a bit outside the visible view).  The origin of the grid is at the edge of the wedge. The calculations are done server-side, and fetched one data point at a time using jQuery/ajax/JSON. This makes the calculations really slow, but it's partly also because I don't want to strain the server too much. I mainly wanted to test how these types of calculations can be done using Django.

This diffraction data is passed to the shaders using a 2D texture, with two bytes ("red" and "green") representing the amplitude of the signal and one byte ("blue") representing the necessary phase data.

WebGL can handle linear interpolation for textures automatically, so the data is interpolated nicely in-between the data points.

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}(\pmb{x},t)=\frac{1}{c^2}\frac{\partial^2 \tilde{p}(\pmb{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). \pmb{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}(\pmb{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(\pmb{x}) + k^2p(\pmb{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:


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. Using local coordinates, they look as following:

Lagrange quadratic shape functions

Lagrange quadratic shape functions

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. A resulting two-dimensional shape function for the node at coordinates (0,0), for example, looks as following:

Two-dimensional Lagrange shape function

Two-dimensional Lagrange shape function

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:


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:


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.


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:




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:


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


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.

Finite element analysis of a beam

This post assumes you already have at least some basic knowledge of structural mechanics. I'm not a mathematician, so if you spot some incorrect notation, please let me know. Also, if some parts of the post are (too) unclear, let me know.

Here are a few examples of the industries where finite element analysis is used, just to give you an idea of how varied the areas where FEA can be used are: acoustics, construction, aerospace engineering, economics, electronics and health care. As with many things, there are quite a lot of details you should understand (at least on some level) before you'll be able to fully make use of all the possibilities the method offers you. In this post, I'll explore the theory behind the basics of the finite element method, or FEM, from a construction point of view. I'll examine a single element, as the same principle applies to multiple elements as well.

Some background information you might find useful when you read the text:

The beam we're going to inspect

The beam we're going to inspect

The finite element method

The finite element method (FEM) is "a numerical technique for finding approximate solutions to boundary value problems" [wiki]. It's based on discretizing the domain we are examining, in all sorts of ways. From the perspective of acoustics and dynamics, or vibration, it allows for us to model complex vibrating systems using differential equations. The following image search gives a lot of images describing the results of the method, if you're not familiar with it. A complex model is discretized into simple elements, which makes the problem easier to understand.

Here, we'll start by examining a simple piece of a beam. In this introduction, we'll assume the beam consists of only one element (which makes things easier to understand). I'll probably return to cases where we use multiple elements in another post. You'll see that things aren't really that different, once you understand how things work with a single element.

Shape functions

I wrote something about shape functions in this post. Shape functions allow for us to approximate a continuous function using discrete values and are very central throughout finite element analysis. They are also relatively easy to understand, so read the post.

In this post, I will use Hermite shape functions. It's just a fancy name for certain types of functions which fit the purpose of this post really nicely. With Hermite shape functions, the derivative (which will describe the slope of the beam) and value (which will describe the deflection of the beam) are 0 at the start and end of the function interval for all shape functions except one at a time. So (assuming you read the post i wrote about shape functions) this means that we now have shape functions for both the deflection of the beam and the slope of the beam at the boundaries of the element.

Hermite basis functions

Hermite basis functions

For example: By examining the figure to the right, you can see that the red shape function has a value of 1 at the start of the function, while the slope is 0. At the end of the function, both the value and the slope are 0. So when we multiply the red function with a scalar value, we will only affect the deflection at the beginning. We won't change the deflection or slope in any other way at the boundaries of the element, as the red function doesn't affect them in any way.

Let's introduce the local coordinate system we will use for one element. We'll see that this makes things a lot easier later on. The local coordinate system is defined using the following definitions:

  • The length of the element we're inspecting is l_e
  • The local coordinate system is given by \xi, which is defined as \xi = \frac{2x}{l_e} - 1

The meaning of these definitions is perhaps made clearer by the following image:

Local coordinates

Local coordinates

The Hermite shape functions for the local coordinate system are now as following:

Name Function
H_1 \frac{1}{4}(2-3\xi+\xi^3)
H_2 \frac{1}{4}(1-\xi-\xi^2+\xi^3)
H_3 \frac{1}{4}(2+3\xi-\xi^3)
H_4 \frac{1}{4}(-1-\xi+\xi^2+\xi^3)
Derivatives of the shape functions in global coordinates

We'll need to take the derivative of the shape functions at some point in this post. It will become clearer why later, but for now I'll just examine how it will be done.

First, we need some relationships between the local and global coordinates. This is because we want to consider the shape functions in local coordinates (the calculations will be easier if the domain of the shape function is constant). We'll see that we can consider the shape functions purely from the view of the local coordinate system, and convert them to global coordinates by multiplying them by a factor taking into account the scaling of the shape function.

This is how it's done.

When x changes the amount of \frac{l_e}{2}, \xi changes the amount of 1. Thus dx = \frac{l_e}{2}d\xi, and \frac{d}{dx} = \frac{2}{l_e}d\xi.

Using the properties defined above, let's take the derivative of a shape function.


= \frac{d}{d\xi}\frac{d\xi}{dx}H_a (the chain rule)

= \frac{d}{d\xi}\frac{2}{l_e}H_a

This means that we can derive the shape function in local coordinates, and convert the result to global coordinates by multiplication (\cdot\frac{2}{l_e}). Following the same principle, the second derivative becomes:

\frac{d^2}{dx^2}H_a = \frac{4}{l_e^2}\frac{d^2}{d\xi^2}H_a

Making sure the shape functions have unit values in global coordinates

We want the shape functions which describe the slope to have a value of 1 when describing the slope in global coordinates. At the moment, the shape functions H_2 and H_4, which will define the slope at the beginning and end of the element, both will have a value of 1 when taking the derivative in local coordinates.

Note that \frac{d}{d\xi}H_{2,4} = \frac{d}{dx}\frac{2}{l_e}H_{2,4}, as per the chain rule. To get \frac{d}{dx}\frac{2}{l_e}H_{2,4} to equal 1, we multiply the shape functions 2 and 4 by \frac{l_e}{2}. As \frac{l_e}{2}\cdot\frac{2}{l_e} = 1, the shape functions will now have unit values in global coordinates.

In the future, the shape functions will be defined as:

N_1 = H_1

N_2 =\frac{l_e}{2}H_2

N_3 = H_3

N_4 =\frac{l_e}{2}H_4

Integrating the shape functions in global coordinates

We'll also need to integrate the shape functions (actually shape functions multiplied by each other, but the principle shown here applies). Once again, we want to do the integration in local coordinates, and convert the result to global coordinates by some multiplying factor.

Let's take the integral of a single shape function and apply the principles shown in conjunction with taking the derivative of the shape function:

\int_{-l_e/2}^{l_e/2}N_a dx

= \int_{-1}^{1}N_a\frac{l_e}{2}d\xi

Converting the integrated function from local to global coordinates in this case then simply means replacing dx with \frac{l_e}{2}d\xi.

Important integrals related to multiplied shape functions

The following table shows the result of integrating multiplied shape functions from -l_e/2 to l_e/2, using the principles shown above. The importance of this table will become clearer later on. We'll do this for the second derivative of the shape functions. For example: the first row in the first column describes the result of \int_{-l_e/2}^{l_e/2} N''_1 N''_1 \mathrm{dx}, the second row the result of \int_{-l_e/2}^{l_e/2} N''_1N''_2 \mathrm{dx}, etc.

All the results can be calculated as following:

\int_{-l_e/2}^{l_e/2} N_{a,xx}N_{b,xx} \mathrm{dx}

= \int_{-1}^1 N_{a,xx} N_{b,xx}\frac{l_e}{2}\mathrm{d\xi}

= \int_{-1}^1 \frac{4}{l_e^2}N_{a,\xi\xi}\frac{4}{l_e^2}N_{b,\xi\xi}\frac{l_e}{2}\mathrm{d\xi}

= \frac{8}{l_e^3}\int_{-1}^1N_{a,\xi\xi}N_{b,\xi\xi}\mathrm{d\xi}

The resulting table is symmetric, so I'll omit the lower half.

N''_1 N''_2 N''_3 N''_4
N''_1 \frac{12}{l_e^3} \frac{6}{l_e^2} -\frac{-12}{l_e^3} \frac{6}{l_e^2}
N''_2 \frac{4}{l_e} -\frac{6}{l_e^2} \frac{2}{l_e}
N''_3 \frac{12}{l_e^3} -\frac{6}{l_e^2}
N''_4 \frac{4}{l_e}

We'll also need the following integrals of the shape functions:

\int_{-1}^1 N_1 \frac{l_e}{2}\mathrm{d\xi}=\frac{l_e}{2}
\int_{-1}^1 N_2\frac{l_e}{2}\mathrm{d\xi}=\frac{l_e^2}{12}
\int_{-1}^1 N_3 \frac{l_e}{2}\mathrm{d\xi}=\frac{l_e}{2}
\int_{-1}^1 N_4\frac{l_e}{2}\mathrm{d\xi}=\frac{l_e^2}{-12}

Important equations related to the behaviour of the beam

We'll assume the material behaves linearly in our analysis. This means that it can be described by Hooke's law, and thus the elastic modulus E is constant. As we're following the Euler-Bernoulli beam theory, we won't have to worry about any other material properties for now. We do need to worry about the second moment of area of the beam, but it's not really that central for this post. So let's just assume that both E (the elastic modulus) and I (the second moment of area) are known, and that they are constant throughout the length of the beam.

The Euler-Bernoulli equation to solve is then as following:

EI\frac{d^4u}{dx^4} - q(x) = 0

u is the deflection of the beam and q is a distributed load (perpendicular force per unit length)

The relationships between the different derivatives of the deflection might be somewhat confusing if you're not familiar with them. Anyway, the following equations which relate forces in a beam to different derivatives of the deflection also become important:

\theta(x) = \frac{du}{dx}

M(x) = -EI\frac{d^2u}{dx^2}

Q(x) = -EI\frac{d^3u}{dx^3}

The weak form of the problem

Now we'll arrive at a central concept, which I've personally found the most difficult to understand: defining the "weak" counterpart (link: Wikipedia's explanation, which does look somewhat mathematical) of the solvable equation.

Background: Boundary conditions

All problems we need to have some sort of boundary conditions (otherwise we don't have a well defined problem, and we can't really solve anything). In the case of the static beam we have the following boundary conditions (given in global coordinates here):

  • u(0) = 0 The deflection at the beginning is 0.
  • u_x(0) = 0 The slope of the beam at the beginning is 0 (the boundary doesn't allow for rotation).
  • EI u_{xxx}(L) = P The shear force at the end of the beam is P, as we'll assume that the point load is applied as a shear force (makes things simpler)
  • EI u_{xx}(L) = 0 The moment at the end of the beam is 0

The boundary conditions which define values for the deflection / shape functions directly (the first two) are called essential boundary conditions. Note that these automatically state that the shape functions defining the deflection and the slope at the beginning of the beam need to be multiplied by 0.

The rest of the boundary conditions are called natural boundary conditions. I don't know if there's some deeper meaning behind the word.

Background: The trial function and the weighting function

We'll use two types of functions to derive the weak form of the problem. We'll later see that these functions will be defined by the shape functions defined earlier.

  1. The solution for u is described by the trial function \delta. As such, we require both the value and slope for \delta to be 0 at the beginning of the beam (as per the essential boundary conditions).
  2. We apply "virtual displacements" using the weighting function, or variation w (this will become clearer later). Even virtual displacements are hindered by essential boundary conditions, so no displacement/variation in our case is possible at the beginning of the beam. The value and the slope for w should thus be 0 at the beginning of the beam.

The weak form

I've personally understood the weak form as following:

  • We assume (or require) that the trial function is the same as the solution throughout the beam, \delta(x) = u(x).
  • This means that the Euler-Bernoulli beam equation is constantly 0 when we replace the deflection u(x) with the trial function \delta(x), as they both satisfy the equation
  • If we integrate the Euler-Bernoulli beam equation throughout the length of the beam, the result will also be 0 (because integrating 0 will be 0, no matter what)
  • We can multiply the beam equation with anything, at any point on the beam, and the integrated result will still be 0 (0 times anything is still 0)
  • We will multiply the beam equation using the weighting function w(x). It's called "weighting function" because it will put more weight on the correct solution at the places where the weighting function will have larger values.
  • By proper choosing of the weighting and trial functions, we can derive a function which is much easier to solve

Alright! So what we want to do is integrate the beam function from the start of the element to the end of the element. We will also multiply the beam function with the weighting function w(x). The solvable equation is then:

\int_{-{l_e/2}}^{l_e/2} w(EI u_{,xxxx} - q) \mathrm{dx} = 0

where u_{,xxxx} stands for the fourth derivative of the deflection u. The next thing we want to do is remove the menacing fourth derivative. We do this by integrating it by parts, which results in the following equation:

|_{-{l_e/2}}^{l_e/2} w EI u_{,xxx} - \int_{-{l_e/2}}^{l_e/2} w_{,x} EI u_{,xxx} \mathrm{dx} - \int_{-{l_e/2}}^{l_e/2} w q \mathrm{dx} = 0

After integrating by parts again (the integral with the third derivative), we get the following equation:

|_{-{l_e/2}}^{l_e/2} w EI u_{,xxx} - |_{-{l_e/2}}^{l_e/2} w_{,x} EI u_{,xx} + \int_{-{l_e/2}}^{l_e/2} w_{,xx} EIu_{,xx} \mathrm{dx} - \int_{-{l_e/2}}^{l_e/2} w q \mathrm{dx} = 0

This is the weak form of the problem.

Sorting out the boundary conditions

At this point, we recognize two things:

  • The first two terms of the weak form describe the boundary values of the problem
  • EI u_{,xxx} can be replaced by -Q (the shear force), and that EI u_{xx} can be replaced by -M (the moment).

This makes the equation look as following:

-w({l_e/2}) Q({l_e/2}) + w(-{l_e/2}) Q(-{l_e/2}) -w_{,x}({l_e/2}) M({l_e/2}) + w_{,x}(-{l_e/2}) M(-{l_e/2}) + \int_{-{l_e/2}}^{l_e/2} w_{,xx} EIu_{,xx} \mathrm{dx} - \int_{-{l_e/2}}^{l_e/2} w q \mathrm{dx} = 0

Using the shape functions

We will now replace w with \sum_{A=1}^n c_A N_A, and u with the trial function \sum_{B=1}^n d_B N_B.

Let's take a look at how the equation looks now (yes, it's starting to get confusing):

-\sum_{A=1}^n c_A N_A({l_e/2}) Q({l_e/2}) +\sum_{A=1}^n c_A N_A(-{l_e/2}) Q(-{l_e/2}) -\sum_{A=1}^n c_A N'_A({l_e/2}) M({l_e/2}) +\sum_{A=1}^n c_A N'_A(-{l_e/2}) M(-{l_e/2})+ EI\int_{-{l_e/2}}^{l_e/2} \sum_{A={l_e/2}}^n c_A N''_A \sum_{B=1}^n d_B N''_B\mathrm{dx} - \int_{-{l_e/2}}^{l_e/2} \sum_{A=1}^n c_A N_A q \mathrm{dx} = 0

If you look carefully, you'll notice that all the multiplications have a common term from the variation, \sum_{A=1}^n c_A. As such, the equation can be simplified as following:

\sum_{A=1}^n c_A \cdot \left(-N_A({l_e/2}) Q({l_e/2}) + N_A(-{l_e/2})Q(-{l_e/2})-N'_A({l_e/2}) M({l_e/2}) + N'_A(-{l_e/2}) M(-{l_e/2})+EI\int_{-{l_e/2}}^{l_e/2} N''_A\sum_{B=1}^n d_B N''_B\mathrm{dx} - \int_{-{l_e/2}}^{l_e/2} N_Aq \mathrm{dx}\right) = 0

The equation has to hold for any value of c_A, because this term is associated with the arbitrary weighting function. The only way the sum is 0 with any value of c_A is if the part of the equation inside the parentheses is always 0. The value inside the parentheses needs to be 0 for all values of A, not just the sum! We now get a group of equations:

-N_A({l_e/2}) Q({l_e/2}) + N_A(-{l_e/2})Q(-{l_e/2})-N'_A({l_e/2}) M({l_e/2}) + N'_A(-{l_e/2}) M(-{l_e/2})+EI\sum_{B=1}^n d_B \int_{-{l_e/2}}^{l_e/2}N''_AN''_B\mathrm{dx} - \int_{-{l_e/2}}^{l_e/2} N_A q \mathrm{dx} = 0

for all A = 1 ... n

Examining the first equation, A = 1

We now have a group of equations where A = 1 ... 4. For clarity's sake, let's go through the case where A = 1:

-N_1({l_e/2}) Q({l_e/2}) + N_1(-{l_e/2})Q(-{l_e/2})-N'_1({l_e/2}) M({l_e/2}) + N'_1(-{l_e/2}) M(-{l_e/2})+EI\sum_{B=1}^n d_B \int_{-{l_e/2}}^{l_e/2}N''_1N''_B\mathrm{dx} - \int_{-{l_e/2}}^{l_e/2} N_1 q \mathrm{dx} = 0

Considering the terms containing the boundary conditions, only N_1(-{l_e/2}) differs from 0 for the case where A = 1. The rest will disappear. Thus we get the following equation:

N_1(-{l_e/2})Q(-{l_e/2})+EI\sum_{B=1}^n d_B \int_{-{l_e/2}}^{l_e/2}N''_1N''_B\mathrm{dx} - \int_{-1}^1 N_1 q \mathrm{dx} = 0

This is where the table with the integrals of the shape functions will come in handy. We will assume that q is constant in this post, to make things simpler. Inserting the correct integrals into the equation, we get:

N_1(-{l_e/2})Q(-{l_e/2})+EI(d_1 \frac{12}{l_e^3} + d_2\frac{6}{l_e^2} + d_3\frac{-12}{l_e^3} + d_4\frac{6}{l_e^2}) - \frac{l_e}{2}q = 0

This can also be written as:

Q(-l_e/2) + \frac{EI}{l_e^3}\begin{bmatrix} 12& 6l_e& -12 &6l_e\end{bmatrix}\begin{Bmatrix}d_1\\ d_2\\ d_3 \\ d_4\end{Bmatrix} - \frac{l_e}{2}q = 0

Writing up all the equations, A = 1 ... 4

We can now write all of the equations A = 1...4 in matrix form. By going through all of the cases in the same manner as above, we get the following system of equations:

\begin{Bmatrix}Q(-l_e/2) \\M(-l_e/2) \\-Q(l_e/2) \\-M(l_e/2) \end{Bmatrix} + \frac{EI}{l_e^3}\begin{bmatrix} 12& 6l_e& -12 &6l_e\\6l_e&4l_e^2&-6l_e&2l_e^2\\-12&-6l_e&12&-6l_e\\6l_e&2l_e^2&-6l_e&4l_e^2\end{bmatrix}\begin{Bmatrix}d_1\\ d_2\\ d_3 \\ d_4\end{Bmatrix} - \begin{Bmatrix}\frac{l_e}{2}\\\frac{l_e^2}{12}\\\frac{l_e}{2}\\-\frac{l_e^2}{12}\end{Bmatrix}q = 0

Inserting the boundary conditions

Remember that there were two kinds of boundary conditions:

  • Essential boundary conditions
  • Natural boundary conditions

The essential boundary conditions define some of the unknowns d_i directly. Stepping a few equations backwards (under "Using the shape functions"), the term c_A becomes 0 for the cases where essential boundary conditions are defined. In our case (a cantilever beam) the deflection and slope at the left node are defined as 0. This means that equations 1 and 2 become undefined, as 0 = 0 always holds true. We should remove these equations from the matrix equation, thus making the equation look like this:

\frac{EI}{l_e^3}\begin{bmatrix}-12&-6l_e&12&-6l_e\\6l_e&2l_e^2&-6l_e&4l_e^2\end{bmatrix}\begin{Bmatrix}d_1\rightarrow0\\ d_2\rightarrow0\\ d_3 \\ d_4\end{Bmatrix} - \begin{Bmatrix}\frac{l_e}{2}q+Q(l_e/2)\\-\frac{l_e^2}{12}q+M(l_e/2)\end{Bmatrix} = 0

We can now insert the rest of the boundary conditions:

\frac{EI}{l_e^3}\begin{bmatrix}12&-6l_e\\-6l_e&4l_e^2\end{bmatrix}\begin{Bmatrix} d_3 \\ d_4\end{Bmatrix} - \begin{Bmatrix}Q\\0\end{Bmatrix} = 0

Solving the equations

Now we've finally arrived at a very central form of the equation, which describes a static situation:

Ku - f = 0

What does the equation define? It states that the elastic constant times the deflection must equal the force. We're examining linear elasticity, so we can draw an analogy to Hooke's law. K represents the spring constant, u the deflection and f the force. K is called the stiffness matrix.

After all the complexity, we've finally arrived at something which is useful and relatively simple to play with. The solution can now be obtained as following:

u = K^{-1}f

which, in our case, will be:

\begin{Bmatrix} d_3 \\ d_4\end{Bmatrix} = \frac{EI}{l_e^3}\begin{bmatrix}12&-6l_e\\-6l_e&4l_e^2\end{bmatrix}^{-1}\begin{Bmatrix}Q\\0\end{Bmatrix}


\begin{Bmatrix} d_3 \\ d_4\end{Bmatrix} =\frac{l_e^3}{EI}\begin{bmatrix}\frac{1}{3}&\frac{1}{2l_e}\\\frac{1}{2l_e}&\frac{1}{l_e^2}\end{bmatrix}\begin{Bmatrix}Q\\0\end{Bmatrix}

This gives the following solutions for the deflection and slope at the end of the cantilever beam (remember that the deflection and slope are 0 at the beginning):

d_3 = \frac{Ql_e^3}{3EI} (deflection)

d_4 = \frac{Ql_e^2}{2EI} (slope)

In this specific case, the solution is identical to the analytical solution.


Finite element analysis can be thought of having a lot of abstraction layers. The math quickly becomes troublesome, but the different phases of the analysis can be separated from each other quite well. For example: in the case above, it's enough if one derives the stiffness matrix once, there's no need to do it any more. Actually, the stiffness matrix which was derived above is readily available from a multitude of sources online.

Similar methods can be utilized to calculate the mass matrix, which describes the inertia of the object in relation to different degrees of freedom. Once the mass matrix is calculated, there are a multitude of equations, based on Newton's laws, which can be used to explore the dynamic behavior of the system.

Shape functions in finite element analysis

Disclaimer: Please note that quadratic shape functions as described below aren't the way you actually usually would model a beam (I'll return to this in a later post). Still, I think they demonstrate the principle nicely.

The concept of a shape function is one of the most important concepts which needs to be understood when doing finite element analysis of continuous problems. It really is a relatively simple concept to grasp (although I've always hated it when textbooks call something simple). To clarify the meaning of shape functions, I will only show the meaning of the shape functions, and leave finite element analysis out of this post.

Quadratic 1D shape functions

Quadratic 1D shape functions

The next picture shows three shape functions, by which we can approximate the deflection of a beam. It's important to note that the shape functions have values of 1 at one node, and 0 at the other nodes. This means that when we sum up the shape functions, the value at each of the nodes is purely determined by a single shape function at a time. This is really important to understand, so take a moment to think about it.

To illustrate the principle without finite element analysis: we know how the beam will deflect (in real world cases of finite element analysis, we usually won't). This is assuming we know the material properties. This means that we also know the deflections at points i, j and k. We multiply each shape function by the deflection value at each respective node (note that one of the shape functions will be multiplied by 0). By summing the multiplied shape function back together, we will get a close approximation of the correct deflection.

u(x) \approx \sum_{A=1}^n c_A N_A(x) = c_i N_i(x) + c_j N_j(x) + c_k N_k(x)

Important properties of the shape functions

The single most important consequence of using shape functions to approximate some continuous function is that we now have a group of known functions (which are multiplied by some unknown scalar). If we ingrate a function times a constant, the constant stays the same: \int c N(x) \mathrm{dx} = c\int N(x) \mathrm{dx}. This means that we can integrate the continuous function without actually knowing what this constant is, thus converting a continuous problem to a discrete one!

Some thoughts

I think an analogy can be drawn to Fourier analysis, or something of the likes (for people involved in DSP). I won't delve into this any deeper, but there's one important point to recognize: the shape function discretization will only take into account low frequencies, leaving high frequency data out. This is important when simulating acoustics or vibration; the shape functions will only be able to carry vibrations accurately up to some frequency. This can be remedied by shrinking the elements or using higher order shape functions. For acoustic propagation in air (for example), a proper element size is often given as being the shortest wavelength examined divided by five or six.