Monthly Archives: July 2013

Listening to vibrations

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

[audio kaistale.com/blog/130722beam1d/kaistale_blog_beamtune.mp3]

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

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:

[audio http://kaistale.com/blog/130722beam1d/Waveform.wav]

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.

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.

I've always found that spectrograms describe situations like this nicely. Here is a spectrogram of the first 0.1 seconds of the resulting sound. The y-axis shows the frequency content, while the x-axis shows time.

Spectrogram

Spectrogram

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.

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.