How to do FEA by hand – Part 2: 2D PLATES!

Quite some time ago, I wrote a post on how to do FEA by hand. Out of necessity, I used simple beam elements, to show the process. Today, thanks to Jakubs’ help, we are continuing this journey into the 2D elements!

You are about to learn how to do FEA by hand for 2D elements. But as I always point out – be very deliberate about why you want to learn this. Somehow I get a feeling that this topic will be mostly read by students before the FEA exam.

Don’t get me wrong, there is nothing wrong with learning before an exam. Hell, I even encourage this!

But if you feel that you have to learn this to start doing FEA designs… please watch this first, and then get back afterward!

This episode is sponsored by Jakub!

I won’t lie: writing the post with doing beams FEA by hand took me a lot of learning. This is not what I do. And to be fair, I never suspected that I would write a post about 2D FEA by hand.

And this would never happen if not for the fact that Jakub wrote to me and proposed, that he can write such a post – of course, I was more than happy to publish it here!

Since this is Jakubs’ post I will be awkwardly silent today.

Let’s go – Jakub is taking the mic!

In the previous blog post Łukasz described the way to use FEM in hand calculations of simple 1D beam elements. Such examples are really good to show the basic concepts of the Finite Element Method. That’s why they usually form the first chapters of books about FEA.

We learned from them how to calculate element stiffness matrices, assemble them into global stiffness matrix, form force, and displacement vectors and then solve the linear algebraic equation to get displacement from which we can obtain additional values like stress and strain.

However when you learn the theoretical aspects of FEA and you solve such examples, at some point you will start wondering:

Where are the differential equations?

The Finite Element Method is a way to solve these equations after all! So how can we solve any problem without using partial differential equations (PDEs)? And perhaps even more importantly, how do FEA programs do this?

It’s a common misconception that these programs use PDEs directly. Actually, they usually utilize the so-called weak form (too complicated and boring to discuss in detail here) which is obtained from the PDE governing a selected physical problem.

What’s more, they don’t use this weak form directly too. To explain how to do it, let’s go back for a while to the previous post where we solved some beam examples.

As you may remember we obtained the stiffness matrix pattern from simple mechanical relationships. You can use such a direct approach for simple 1D elements. But what about the most common 2D and 3D elements?

Their stiffness matrices can’t be obtained directly and thus they are calculated in another way. That’s where PDEs become needed. When we combine the PDE, constitutive equations, and shape functions we can get the weak form (or so-called functional if we use the variational method).

Then we extract the formulas for the stiffness matrix as well as the force vector for the specific element types. The rest is similar to what we’ve seen in the previous post:

• Firstly, you calculate stiffness matrices for all elements using the formula
• Then you assemble them
• Next, you calculate the force vector
• You add boundary conditions
• And finally, you solve the equation.

This is how all FEA programs work in practice.

The task at hand…

Let’s proceed to the example to see how can we do it by hand. We will solve a simple flat plate (plane stress member) fixed on one side and subjected to a concentrated force at a single node on the other side:

We will use a “typical structural steel”, so Young’s modulus is 210 000 MPa, and Poisson’s ratio is 0.3.

The first thing that we need to do is discretize the plate. This means that we will divide it into finite elements. In this case, two triangular elements should be enough. More would make it too hard to solve without a computer because of the size of the matrices. Let’s also number the elements and nodes:

The math slowly starts

The stiffness matrix for such plane stress element looks like this:

[k] = t A [B]^T [D] [B]

Where:

t – element thickness

A – area of the element

[B] – strain-displacement matrix

[D] – material property matrix

T in the superscript indicates matrix transposition (even Łukasz knows this!).

FEA programs use integration over the volume of the element but fortunately, in this simple case with constant arguments, there’s no need to do it. You can calculate the element area using nodal coordinates. Luckily, in this example, we can use simple geometrical relations. We will use symbols as seen below:

Knowing the above numbering, we can calculate the strain-displacement matrix for a linear triangular element from:

where β and γ are shape-function coefficients. You can calculate them using nodal coordinates:

The key here is that nodes are always numbered counterclockwise.

It would be possible to solve the problem using opposite node numbering. But please note that all formulas that we use here apply to the counterclockwise approach. The same is done in FEA software.

You may wonder what to do if there are multiple elements with node numbers different than 1,2 and 3. Do we take their nodes in the order of their numbers (for example 16,17,19)?

No, we actually don’t. Instead, we choose any node as the “first one” and then go counterclockwise treating the two subsequent nodes as the second and third nodes in the aforementioned formulas, respectively.

This is called global and local node numbering.

Global node numbers are the ones given when considering the whole computational domain. They are unique for each node (16,17 and 19 for instance).

Then when we calculate matrices for each element we use local numbering. In this case, we prescribe numbers 1,2, and 3 to the nodes keeping the rule of counterclockwise numbering.

In the end, we will also need the formula for the material property matrix. For plane stress, it’s given by:

Where:

E – is the Young’s Modulus

v – is the Poisson’s ratio

The numbers are about to start!

Now that we have all the required formulas, we may proceed to the solution to our problem. Let’s start from the first element and calculate its matrices using local node numbering:

Thanks to everything above, we can finally calculate the stiffness matrix for this element:

Rows and columns in local stiffness matrices correspond to the nodal degrees of freedom in global coordinates.

In the case of the first element global node numbers, in the order given by local numbering, are 1, 3, 2 so the following Degrees of Freedom (DOFs) are included in this stiffness matrix:

1 (x for node 1), 2 (y for node 1), 5 (x for node 3), 6 (y for node 3), 3 (x for node 2), 4 (y for node 2).

Ok… I took me a bit of time to actually digest this.

If you got this, move on skipping the next chapter… below I will explain what just happened a bit more!

Łukasz adds a small explanation!

OK, let’s digest the numbering shenanigans one bite at a time. I will simply state the facts we already know, with some masterfully crafted images in the hope this will make things easier to see.

FACT 1: We have a model, that has a GLOBAL node numbering. You already saw those at the beginning, here is a short reminder (Global Node Numbers marked in green circles):

FACT 2: Each of those GLOBAL nodes has 2 Degrees of Freedom (since it’s a plane stress example). There is a potential to confuse stuff… so I will add a “G” subscript to global coordinates. The entire set of GLOBAL COORDINATES is given below:

FACT 3: Before we can solve the global problem, we will calculate each element “separately”. And to solve each element separately, we need to assign LOCAL NODE NUMBERS for each element. Those are always “1, 2, 3” and those are assigned counterclockwise. I will use Element 1 from above, and I will mark the LOCAL node numbers in red:

Before we move on… please note one thing! Global Node No. 3 (in green on the left) has a “local” number 2 (in red on the right) when only a “local” element 1 is considered.

Such numbering problems will always be an issue. Although it’s tricky only if you start numbering global nodes from “1”… since local nodes (the red ones!) will always be numbered 1, 2, and 3 in a counterclockwise direction.

FACT 4: I’m sure you know where this is going… the LOCAL NODES have their Local coordinates. We already numbered them before, but just so everything is in one place:

Of course, the above LOCAL COORDINATES have subscripts fitting local node numbers. Just as the global coordinates have subscripts fitting to global node numbers.

FACT 5: Now we are ready to once again take a look at what Jakub did with the local stiffness matrix for Element 1. Of course, the math is there if you want to review it again.

For me, the tricky part was the “rows and columns numbering”. Of course, the single-element stiffness matrix is assembled in the local system of that element, which means that we can initially number the columns this way:

But such an approach to numbering isn’t super useful (but orderly and makes sense!). It’s not useful, since every element has local nodes numbered from 1 to 3 in a counterclockwise direction. So for every element, this local stiffness matrix has the same row and column numbers!

And this is where we go “back” to global coordinates. We know what the GLOBAL Node Numbers are, so we can “map” the above matrix with rows and columns “fitting” the global system.

As you could see above, Node 1 has the same number in our example in both local and global coordinates (totally doesn’t have to be the case BTW!). Local node 2 (red) was global node 3 (green), and local node 3 (red) was global node 2 (green).

This means, that with the above we could add:

It’s not very comfortable to use subscripts all the time. This is why we will use “normal numbers” to number GLOBAL DOFs. Simply put:

x₁ will be “1”

y₁ will be “2”

x₂ will be “3”

y₂ will be “4”

x₃ will be “5”

y₃ will be “6”

x₄ will be “7”

y₄ will be “8”

The above is for GLOBAL DOFs, so for the DOFs marked in green above. If we simply substitute the symbols with the numbers, we will get the stiffness matrix Jakub showed:

Well… I hope that the numbering is clear now : )

I’m running away again, and Jakub is taking over the stage (once more!).

The journey continues

Now let’s calculate the matrices for the second element assigning local numbers again. The material property matrix will be identical to the one for the first element since they are made of the same material:

And yet again, thanks to the above we can build the stiffness matrix for the second element, as you can see below:

In the case of the second element, the following global numbered nodes belong to this element: 1,4,3. So the DOFs here are: 1 (x for node 1), 2 (y for node 1), 7 (x for node 4), 8 (y for node 4), 5 (x for node 3), 6 (y for node 3).

I hope that at this point, the green numbers aren’t surprising!

Note that Element 2 had LOCAL node number 1 in the GLOBAL node 1. This meant, that local node 2 (remember, always counterclockwise) would be in global node 4, and local node 3 will be in global node 3.

The rest works just as I described before.

Next phase – global stiffness matrix

Since we have the stiffness matrices for both elements, we can proceed to the next step – the assemblage of these matrices into a global stiffness matrix.

We can do this by just adding proper terms to each other and forming a new 8×8 matrix. For example take the 25^th element (column 2, row 5 – using GLOBAL DOF-based Numbering) from the first local stiffness matrix and add it to the corresponding element of the second local stiffness matrix to get the 25^th element of the global stiffness matrix.

This is the outcome of this operation:

The above sounds like a magic spell?

Don’t worry – you already know all that is needed. But just to be sure, let’s add a simple explanation here as well.

So far, we have two LOCAL stiffness matrixes (one for each element). We already managed to “assign” rows and columns in both those matrixes to GLOBAL coordinates (the green numbers). Both of those look like this:

And now, all we need to do is to make an 8×8 GLOBAL Stiffness Matrix. Our model in total has 8 GLOBAL DOFs. This is why our Global Matrix will be 8×8. Those will be the “green” DOFs starting from 1 to 8 both in columns and rows. Like this:

The cool thing is, that we have a Global Stiffness Matrix. The bad thing is, that it’s empty! Luckily Jakub already explained how to populate it.

All you need is to take numbers from the same row and column (GLOBAL coordinates so green numbers!) and add them together. And then put the outcome into the same row and column on a global stiffness matrix.

If you make a graphic, it’s actually simple:

Ok, I will shut up now!

Go Jakub!

Jakub takes the stage… for the 3rd time!

After doing this we should make sure that the matrix is symmetric about its diagonal. Then we can eliminate constrained DOFs (1,2,3,4) by crossing out columns and rows related to these DOFs.

Let’s calculate the displacements!

Now the equation that has to be solved is:

[K] {u} = {f}

Where:

{u} – displacement vector

{f} – external forces vector

To solve this equation we may invert the stiffness matrix:

{u}=[K]^-1 {f}

The results from the above are:

– node 3, x displacement: u₅ = 0.03 mm

– node 3, y displacement: u₆ = -0.154 mm

– node 4, x displacement: u₇ = -0.03 mm

– node 4, y displacement: u₈ = -0.153 mm

Finally, stresses and strains!

Of course, we can use these nodal displacements to calculate strains and stresses:

{ε} = [B] {u}

{σ} = [D] {ε} = [D] [B] {u}

The stresses are:

• For the first element:

• For the second element:

The results were confirmed by comparing them to the solution obtained with an open-source FEA solver CalculiX. They are in perfect agreement. Exemplary screenshots showing displacements are shown below. Those outcomes were visualized using an open-source pre/post processor for CalculiX called PrePoMax:

As you can see, the manual calculations are really laborious even in such simple cases. Even with two triangular elements, matrices become quite large. If you would like to make the calculations a little easier, you can use the Computer Algebra System (CAS) software. Free programs of such kind include SMath Studio and wxMaxima.

The great advantage of this approach is that we can learn how FEA software actually works in linear static plane stress analyses.

If you want to find more examples like this, check out Udai Borker’s book “Finite Element Analysis”. Other books about FEM that include similar problems are those written by P. Seshu, S.S. Bhavikatti, T.R. Chandrupatla, and D.L. Logan.

Łukasz is rushing to the stage yet again!

First of all, thank you Jakub for writing such a cool guide. While editing your work I really learned something! I guess this is one of the perks of having such a blog as mine!

It’s always fun when the “stars align” and the final outcome is just right! This was expected of course. But I well remember my math teacher at Uni who said in one of the lectures: “I will do the riskiest thing a lecturer can do! We will solve the same integral with the second method… hoping to get the same result!”.

A classical cruel twist at the end!

If you’ve read my first post on how to do FEA by hand, you already know what is coming!

I’m a super practical guy, so I always wonder: is this useful?

And the answer is… well, no, at least not really!

I mean, it’s always better to “know something” than “not to know it”. I totally agree. But learning anything takes time, and learning how your solver does the math seems like “doubling skills”.

To some extent it would be equivalent to me, learning how accounting works, so I can “double-check” my accountant. Technically, this is reasonable, but I hire a good accountant for a reason, right?

I know, that nowadays, you will be forced to learn things like that at Universities and courses. So let’s wonder for a second not about what you’ve learned, but what you haven’t learned here today!

These are some questions, that desperately need answers if that would be a “design case”:

Of course, both me and Jakub are perfectly aware of this! And I totally understand that the goal of this post was NOT to give you those answers… but to show you how the solver does the math.

But in a world where so many people feel that “math is all in FEA”, this may be my only chance to show you, that design is so much more than FEA-math! And this is why I’m mentioning it here!

Knowing how solver works seem cool – I would definitely like to know that in detail. But learning that would take months (if not years!) so instead, I decided to focus on the design skills I will need if I want to use my solver correctly.

I’m not an accountant in my company… I run things here, and I let my accountant do her work! I don’t feel the need to double-check the accounting at all. But I have an overall understanding of things, so if something “outrageously wrong” would pop up, I would notice it and ask my accountant. I don’t need to know the details of her work, I just need to be aware of what outcomes I should expect. And just as with the solver, this has nothing to do with “the math”, and all to do with understanding how things work. You don’t need the math for that.

If you want to learn some practical design/FEA things, sign up for my free course below!