Can You do Finite Element Analysis by Hand?

It’s the beginning of 2009. I’m sitting in one of the lectures of my Ph.D. program. Finally, we will start to learn about Finite Element Method! My company is slowly developing, and I’m really hopeful that this will allow me to spread my wings! And then professor says: “let’s start with basic matrix operations”…

You really can do Finite Element Analysis by hand. The real question is, why would you?! Even with computer aid, you would be more or less bound to linear analysis… as nonlinear stuff is actually pretty difficult to develop. But just for giggles, let’s take a look at how this could be done!

A fair warning: here be mathematics! While I try to avoid that at all cost, apparently this is a popular topic. I figured I will write on it, adding my anty-math twist to the mix! Take a look and let me know what do you think!

This is a long post (solid 30 min to read it all!), so I’ve made a content list.

This Article Content List:

• Simple stuff
• Element sitffness matrix
• Simple example
• Complex example
  • Global stiffness matrix
  • Load and displacement vectors
  • Matrix operations
  • Outcomes checking
  • Internal forces
  • Shape functions, strains and stresses
• A cruel twist!

Simple Stuff

I think that every “basic” FEA book or guide starts in the same way. Well, they are more about “how FEA works”, but it’s the same thing.

Usually, such guides focus on rod elements. Without much comment, the beginning is usually something like this:

There is a really good reason why that is the case. The above is the simplest possible element: a 1 DoF “spring”. And what I mean is not only that this is a 1D beam element (but of course this is included!). What is also important is, that our element can carry only tension or compression. The only deformation we “see” is along the length of the element.

Why this is important? Well… first of all, because such an element won’t “see” bending nor shear. It is only carrying a normal force and nothing more!

This makes the problem simple, and we can easily follow the math in such a case! Let’s take a look!

Firstly, a few basic equations you should easily recognize. The only thing I did there is, that the elongation of our rod is u2-u1 based on the drawing above. Technically, this could be defined differently because it only depends on what you would consider “positive”. However, the above setup seems to be the most “popular” (at least as far as explaining goes).

The above shows you the basic relation between stress and strain. I also explained how all of the parameters can be established in the case of our rod (Note that Young’s Modulus is a material constant).

I remember well, that in primary school I had problems in physics class. When I had no idea how to solve a task I would just write down all the equations I know. Then I would just substitute stuff between them until something worked. This approach rarely failed me… so let’s try it here!

Now, let’s get back for a second to our rod. We would love it to be in equilibrium right? This means that the forces on both ends should have equal values, but work in opposite directions. In such a case…

You can notice that I’ve decided to make the F1 force “negative”. One of them had to be (for them to work in the opposing directions), and F1 just “draw the shorter straw”. If I would do it the other way around it would still work at the end : )

The Place Where Matrixes are Born…

So far, we have played a bit with a few relatively well-known equations. This is not exactly what your FEA solver does, but we are getting closer.

It’s easy to notice, that if you would know the forces (or deformations) you could relatively easily solve the above. It’s just a set of 2 “linked” equations. But if there would be plenty of those (and there will be!) it’s much more convenient to solve those in a matrix form. Well… maybe not convenient for you, but definitely more convenient for your computer!

I don’t want to go into matrix notation here, but if you don’t know how vectors and matrixes interact check this out lesson from Khan Academy!

There are also other lessons in the subject so if you need more basics or if you want to go deeper into this check out the linear algebra portion of what Khan Academy offers.

The above 2 equations are pretty easy to write in the matrix form. You just take the forces and put them in one vector, and you do the same with corresponding deformations. The rest is just matrix multiplication really:

Now you can start to notice things. First of all, if we gather up forces and displacements into vectors there is some “stuff” left. This stuff, as you can see above will form a matrix. The terms there, clearly correspond to how “rigid” the model is, and this is why it’s called the stiffness matrix. This is how your solver “knows” what cross-sections and thicknesses you used in your model, along with correct material properties.

For now, it’s only a matrix for a single element (rod in our case), but if we would have more elements in our example we would “bash them” together into one big “global” stiffness matrix.

Keeping It Real

So far it’s simple I hope… normally I would test if the previous equations “work” but the above is not a “real case”. We just assumed that we will apply equal forces at both ends of the element to keep it stable. Obviously, according to our assumptions, this would work. Let’s try to do something more “real”:

As you can see, this would be the simplest static task you can have. There is a 2m long rod loaded on one end and supported on another end. The difference here is, that it’s actually a “solvable” task, so we can try to use FEA to solve it, and gain some confidence that it actually works!

All I did below, is simply input the “real numbers” from the example into the equation from the previous part. All the logic still applies, nothing fancy happened. You just need to remember that in a place where support is, deformations are equal to zero (hence u1=0) but there is a force applied (reaction force)… you just don’t know how big the reaction force is (this is why F1 is written as… F1 – we don’t know the value yet):

I’m pretty sure there is plenty of ways in which you could solve such a set of equations. Heck, at the Uni as a student I most likely knew at least a few myself…

… a quick google later I know that the easiest for the 2×2 matrix is to actually write the equations 😛

This sucks, as I will have to make a bigger example to use any sort of matrix operations (!). But there is also another “rule”. Somewhere it was called “elimination” but I’m not too fancy with math “names” so don’t quote me on this. Basically, if you have a “u” that is equal to “0” you can cross out the row containing the zero, as well as the corresponding column. In our case, we can cross out the first row (because u1 = 0) and the first column (because it is corresponding to the first row obviously):

This way we formed a really simple equation:

We got the u2 deformation. But at this stage, we can also check if what we are getting makes any sense at all. Since u1 = 0, u2 is the deformation of node 2, but it will also be equal to total element elongation. If you remember the equation, you know that this is what we did above. But in any case, you can google for it, and check what the total elongation of our rod would be:

Since we already solved u2, we can get back and solve F1… the only other thing that is left to solve. I got back to the original matrix equation and “multiplied out” the first equation (for F1):

Notice, that the value of the force is negative. When we started I assumed that all forces “to the right” are positive, and hence, reaction forces working in the opposite direction must be negative! It seems we can do something with FEA by hand!

A Bit More Complex, a Bit More Real

Ok, I will have to incredibly dumb down this example if we want to solve it with FEA by hand. It’s just that doing more elements, or more degree of freedom (i.e. to capture bending) is simply too much! I want this example to be relatively easy to follow!

Imagine a steel column in a hall building with a crane. There is a roof supported by the column at the top (loading the column with 120kN of compression). Lower, the crane is operating and the load is 160kN at that level. The question is: what is the maximal load at the bottom of the column, and what is the vertical deformation on the top and crane level.

Important note!

This is extremely dmubed down! Normally you would have to include eccentricities of the applied load, but those would create bending. This in turn would require additional DOF per node… and the matrixes would be too big to solve by hand in a short time!

You know… use the software when you try to solve real stuff! This is just for demonstration/fun!

All right then! Let’s take a swing at this.

The first step is pretty easy. Since elements are loaded at nodes I need to have a node at every place where the load or support is applied. However, I’ve decided to divide the 8000mm part of the column into 2 elements. It’s not even about the length of the elements (outcomes would be the same)… I simply wanted to make an “unloaded” node (2). This will allow me to point out something later.

The second part requires a small description. We need to build a global stiffness matrix!

Bashing up the Global Stiffness Matrix

In the last example, we had a very simple case. There was only one element, so not a lot of worries. I told you that if there are more elements, the stiffness matrix of all elements is “bunch up” into one big global stiffness matrix. This is something that we will have to do this time!

Last time we assumed that the “positive” direction is to the right… now let’s assume it’s upward (I guess I’m an optimist like that). As you already know the stiffness matrix of a single element looks like this:

We already know that it “works”, but note one important thing. We already assumed how the nodes are numbered, and that the numbering “increases” toward a positive direction. This is why in the column node (1) is at the bottom, while node (4) is a the top. I could change that for sure, but it would just make the thing more complex – let’s stick to the basics : )

Firstly, let’s build the individual stiffness matrixes for all elements separately. All I do is just take the “template” above, and input the variables for each element. To make the matrixes shorter I will already calculate values. This is a steel column so I will use E = 210GPa. Values in the matrixes themselves are in MN/m.

As you already know, those “small” element stiffness matrixes correspond with certain global Degrees of Freedom. (DoF). For instance, the stiffness matrix of element A corresponds with deformations and loads in nodes (1) and (2). Stiffness matrix of element C corresponds with deformations and loads in nodes (3) and (4). This can be visualized like that (using our element B as an example):

What is cool about each FEA model is, that the nodes are connecting elements together. As you can see above, element A and element B “share” the node (2). This means that both elements A and B “add” something to the “situation” in the node (2). When we will build a global stiffness matrix it’s made from “blocks”. Since we already know that matrixes of each element correspond to specific DOF, we can denote them as:

Finally, we can assemble the global stiffness matrix! Notice, that we have 4 DOF in our system so the matrix will be 4×4. Since we already know which elements contribute to which Degrees of Freedom (DOF) we can use all the above matrixes to create the global stiffness matrix:

Wow… now that was a lot of matrixes! Anyway, we know the schematic for building the global stiffness matrix (click on the image above for a larger version if needed). The trick is, that on the “overlaps” components of the matrix are added together. In such a case, we create a matrix for our task:

Yesss! We have the global stiffness matrix! How awesome is that!?

Well… on its own not too much to be honest. We still need the deformation and load vectors and to solve the whole thing. If you are following me from the start in one go… maybe a small break would do at this point? I know I need one before I will write more!

Do you want to practice or learn more about the assembly of stiffness matrix?

If so, definitely check lecture from MIT Open Courseware by the legendary professor K. J. Bathe.

Load and Displacement Vectors

Having the first part of the task, we can move on. While assembling global stiffness matrix may be a bit time-consuming and overwhelming, load and displacement vectors are super friendly! It’s been a while, so let’s take a look at how the model really looks like, and what those vectors are:

Let’s wonder for a second what just happened. As before I marked the corresponding DOF in the circle. Firstly the force “F”. For (3) and (4) loads are applied, so it’s obvious. The only “difficult” thing is, that we have to remember that we agreed that the positive direction is upward – so the applied load is negative! The curious thing happens in (2). Note, that there is no load applied in node 2, so the value is 0. However, it doesn’t mean that there is no “internal force” in the elements meeting at node 2. It only means that we are not applying an external load in that DOF!

Just as previously, the reaction force is applied in (1). We don’t know the value yet, so I marked it as R1.

Deformations are so much easier! We only know that at DOF (1) deformation is 0 (this is where the support is). Other values are unknown – we will calculate them in a second.

On a side note:

It’s good to understand that corresponding DOF number is not a node number. Here it’s like that, but in “real” FEA problems it’s not the case! We are solving a very specyfic problem where each node has only 1 DOF (along the element). Normally beams have 6 or 7 DOF per node, which means that first 6 or 7 DOF’s are still at node 1!

Numbers in the circles refere to “global DOF”. This can be translation along X in node 1 or rotation around Y in node 1 in “real” cases. Here, since we have only 1 DOF per node the node and DOF number are the same. I thought it’s a good idea to point this out!

Matrix Operations and All the Jazz!

This is the moment you’ve been waiting for (I’m sure!). We have all the things needed to make this one “final” move! Let’s gather up all the things we did so far… and wonder how to solve this!

Firstly, since u1 = 0 we will cross out the first row/column to simplify the problem a bit:

I knew that, but I admit I double-checked just in case. The matrix above is a set of “normal” equations. This means that the normal “algebra” still works. It’s like solving a set of equations in high school. You can multiply on both sides, or add/subtract equations from each other. The goal is, to have the “bottom left triangle” in matrix equal to 0 (positions {2,1}, {3,1}, and {3,2}). This is basically solved in such a form, as starting from the last equation you just calculate each parameter per matrix row (going upward). To do this I did some mathematical “harm” to the initial matrix. Let’s take a look:

STEP 1: (Row 2) = (Row 2) + 0.5 (Row 1)

This allows me to get a 0 in position {2,1}

STEP 2: (Row 3) = (Row 3) + 350/743.5 (Row 2) = (Row 3) + 0.4707 (Row 2)

I already have a 0 in position {3,1}, but I need a 0 in position {3,2} as well

STEP 3: Use what you’ve got!

With the last row having {0, 0, 185.3} it’s very easy to calculate u4. There is only one parameter in the last equation so:

Having u4, we can use the higher row to calculate u3:

Now, we have u3 and u4… so we can solve the first row for u2. At this stage I think it’s clear why I wanted to have zeros in the bottom left corner of the matrix – such a state makes solving really easy:

STEP 4: Don’t forget about what you crossed out!

We already solve all deformations, but we are missing a reaction force. I’ve crossed out the first row, as it was quicker this way. But knowing u2 it’s time to get back to the “real” first row of our system that contains R1:

Yes! We actually did it. I must say that I was writing this post step by step doing the math as I went. I was wondering if this will play out right. After all, I’ve only made a few calculation errors I had to fix, and it seems this is doable 🙂

This is how you do FEA by hand… at least part of it! Let’s see if the outcomes have any sense!

Drumroll… Outcomes Checking!

The reaction force is pretty simple. Value is positive, so reaction force is upward, working against the applied load. No complaints there!

Total load was 280kN, the reaction I’ve got is 279.4kN. It seems I’m missing 0.2% of the outcome. This is due to round up errors. I was only using 3-4 digits, and those add up. No worries though, the computer uses much more than 4 digits… I’m just lazy!

With deformations, it’s a bit more tricky, but we can check this as well.

Firstly I will calculate the shortening of the “main column”. It was represented with elements A and B in our small FEA task. Here, there is no need to treat those separately, as they are constantly loaded, with constant cross-section and material. Just remember that the total applied load to this part of the column is 120kN from the top with the additional 160kN at the start – we need to sum those up!

The total shortening of the “main column” should be equal to deformation downward (negative value in such case) at node 3. Our FEA “solver” spit out -0.71mm – a spectacular success!

Shortening of the top “smaller column” alone is:

This alone isn’t helpful for checking, but the total deformation of the column (so deformations at node 4) should be equal to the sum of the values above:

Our solver showed displacement as -1.05mm… yet another spectacular victory!

Internal Forces

Great! Above we’ve managed to solve displacements and reaction forces of the system. We already know a few useful things. But obviously not all of them. Now, let’s wonder how to calculate the internal forces in each element – those can come in handy in design right?

Firstly, let’s sum up what we already know. We managed to calculate the reaction force and displacements in our system, and in the process, we had to create the stiffness matrix for each and every element. All of the outcomes are below:

Calculation of internal forces is actually pretty simple. We will just make the equilibrium equation for each individual element, but instead of active forces, we will use internal forces in the element. Since we know the deformations already, this should work out nicely. In case you were wondering, when I will mark an internal force in subscript there will be a node number near which the value is “located”. In superscript the element letter. This way we won’t get lost here : )

As you can see, finding the internal forces isn’t so bad! We are losing some accuracy as before (in element C the force is 119kN instead of 120kN), but this is because I’m lazy and I round up numbers pretty fast!

Note, that in elements C and B we did not use the applied loads. However, the difference in internal forces, of course, is present. At node 3 element C shows 119kN while element B shows 279.4. The difference is, that at node 3 (think about it as “between” elements C and B) additional load of 160kN is applied!

Strains, Stresses and Shape Functions

Strain (and resulting stress) isn’t really a nodal value. You can, of course, provide the value in the element near the node (like we did with internal forces) but there is a twist. For strain to appear something with “initial length” must change its length. Nodes don’t have length at all… so strains are associated with elements rather than nodes.

The twist is… we don’t know yet how the elements deform! I mean we know how the beginning and the end of each element deforms, and we are using linear elements, so it’s painfully obvious. But our solver doesn’t know how the elements deform… and we need to do something about this if we want to use it!

To help our precious machine we have to “explain” what is happening in each element based on what values it has at the beginning and at the end. This is what shape functions are for!

Since we are using a linear element, the shape functions are relatively easy. This is how a typical element looks like:

I’m using “i” as the “beginning number” and “j” as the “ending number”. Obviously depending on the element different DOF will be at the beginning and end so using “real” numbers isn’t the best choice!

This is a very simple linear ramp, but we want to use coordinates in the mix, so the solver knows where it is. Also, instead of an “equation”, I will prepare a vector, you should use with the deformation vector. This way it will have an option to use the deformation vector we already have! This vector is the shape function for our element. Usually, it’s described as “N”:

Note, that for x = xi the first value is 1, while the second is 0. On the same note for x = xj the first value is 0, while the second is 1. This is because the shape functions are always “calibrated” in such a way, that at a node, only the value from that node impacts the outcome. It can also be seen that regardless of the value of x, both terms in the vector will add up to 1. Think about the shape functions, not as any physical value… it’s just a mathematical tool to “spread out and average” nodal values to the entire finite element.

Of course, on its own shape function is useless. But using it, we can describe the deformation state of the entire element. Such deformation will be a function of coordinate x because depending on the x value the deformation in the given point on the element will vary. Let’s call the “elemental” deformations U(x), while the nodal deformation will use lower case u as before. In such a case:

The above looks nice and dandy, but let’s do a quick test if it is any good. Element “B” seems to be a good candidate. It starts with node 2 (u2 = -0.355mm) and ends with node 3 (u3 = -0.71mm). It’s pretty obvious that in the middle of the element the deformation should be -0.533mm, Well, let’s see if this is also obvious to our solver:

Of course, we already know the coordinates. I will use the global ones here, assuming that x=0 is at the bottom of the column. This means that xi, being the beginning of element B will be at xi = 4m, the xj = 8m and the place we want (the middle of the element) x = 6m. With this in mind I can simply solve the above:

Nice! It seems that our solver has a new trick – it knows how to distribute displacements along finite elements. In itself, it’s not all that fancy, but there is a neat thing going on in here. And that is strain is a derivative of displacements! Now, since we have the displacement function, we can actually mathematically calculate strain… finally we are getting somewhere! Normally there would be a matrix [D] describing to the solver which derivatives to use on which displacements (in a real 3D problem). But here, we know that we will do a d/dx, since this is the only coordinate we have! Regardless, the outcome is usually denoted with [B], so let’s keep the classical notation:

Above I could “throw out” [u] from the derivative, as those are constant values of nodal deformations. Notice, that B(x) doesn’t really depend on x in our case! We are using elements with the linear distribution of deformation along the length… so the resulting strains will be constant on the entire element!

 Let’s see what strain will we get:

It’s easy to see, that finally, all the operations lead to the correct one: displacement at the end minus the displacement at the beginning divided by length. It looks like we can calculate strains as well! Following the same method, strain in element A is the same and in element B, and strain in element C is 0.011333%.

Take heart! Only stresses are missing and we will call it a day!

Luckily for us, the relation between stress and strain is primitive at best. Sure, out FEA solver should have a matrix of “potential” material parameters but in our case, it’s obvious that we will use Young’s modulus:

Finally, we did everything I wanted to cover here… and that was something! Take a minute to celebrate, but if you want to learn what I think is most important here… read on!

A Cruel Twist!

I must admit that I was surprised that I pulled this off. I never did FEA by hand, and it took me several hours to find materials in books to make this work. But I did it, and I definitely learned something new. I think that this post contains more equations than all other posts I’ve written in the last 3 years O.o’

But being me, I want to finish in a traditional way when it comes to math: is this useful?

The answer is: no, not really. I mean if I would be aiming at writing my own solver then perhaps this would be the first step. But if I would like to design that column from the example? Don’t be ridiculous – this is not how it’s done! Even if I would expand the example to 6 DOF to take bending into account this would be far from “perfect”.

Here are some questions that are important in solving this problem from a design standpoint… that has nothing to do with anything I wrote today. Answering those is far more important than learning how to operate on matrixes. After all, you can use a solver done by someone else (there are plenty, even open-source), but not a single solver will answer you those questions:

Look, I admit that this is far-fetched. I know, that the above does not cover today’s topic at all, and instead goes into “steel design” much more than into FEA. Trust me, I really know there are codes and rules to follow in such designs, and there is specific software that will aid you in such a task to some degree.

I wrote it here because I see that more and more people do FEA, but have no clue on what is going on. I mean I’m all for “FEA democratization”. For me, FEA was a wonderful tool that greatly helped my career. I wish everybody to be even more successful than I am.

But… this is still engineering for crying out loud! The fact that you know how to move a matrix around does not mean that you know how to design stuff! In fact, I would say that those are two completely different things!

If you are just starting, please (please!) remember – FEA is like a calculator, it only does math operations. Designing is far more complex than that, and requires a LOT of knowledge well outside of the math operations! You need to understand how to mesh your model, how to load is and how to support it. Finally, you need to know what analysis is the “right one” and how to interpret outcomes (and what to pay attention to). There are just so many things that are important in design, that have nothing to do with what I covered today.

I know that at Uni they mostly (if not only) talk about this stuff in a way I did in this post… but this doesn’t mean that this is the only thing that exists. Far from it… it’s not even the beginning. After all, I’ve been successfully using FEA in design for a decade now… and I learned how to solve those freaking equations a few hours ago…

Want to Learn Something Useful?

Definitely check out my FREE FEA essentials course. You can get it by subscribing below.