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}{dx}H_a

= \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}

or:

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

Epilogue

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.

Leave a Reply

Your email address will not be published. Required fields are marked *

Time limit is exhausted. Please reload CAPTCHA.