Finite Element Method Programming - Introduction

The problem with finite element analysis learning is that everybody is starting it from theory and after you are completely brainwashed, they start doing practical calculations. The theory for today FEM is so wide, that when you go for this super interesting practical case and you finally will apply your FE knowledge, suddenly ... it occurs you just do not remember theory... 

My approach here will be a little bit different. I will start explanation from a simple practical case, and I will be coming back to theory during the explanation. In this post, I will show 1D bar element and the method of basic programming it in C#. I assume you want to learn C# and FEM so I will try to say about both of them simultaneously.

First thing first! The task is to calculate 1D bar, like on this figure:

For simplification, we assumed: this is simple bar element, only compression or tension is allowed, there is no gravity, the deformation can be only according to the one direction (X). This kind of finite element can be used for the bars in a truss.

To check if algorithms are working well, we always need some benchmark. It is necessary to check even the simplest codes with the different method! In this case, we will go into the theory of elasticity to solve tension element. Before we start the reminder of definitions. If we have an element, it has some basic length l0. When we apply the tension force P, the element will stretch to the actual length l.  

Most of the building materials have got the elastic phase of this deformation, it means that deformations will grow up linear, as the load is growing up linear. In this elastic phase we can find the relation between the stress and strain, which is, of course, Young modulus E. Young modulus is defining the stiffness of the material, but graphically can be explained as the gradient of the stress-strain curve (the bigger E number is the stiffer the material is).   

What is, in this case, the stress and strain? Stress 𝞼 is just force P divided by area A. In general, stress is a physical quantity to measure the state of the material and normally we need tensor to described it, but right now P/A is enough. Strain 𝞮, is used to measure the deformation of the material. The formulas with units you will find here:

Ok, so we have all definitions to solve this primitive element. Let's do some calculation. First, we have to clarify what we are looking for? We are looking for displacement of the node on the right. The displacement is the difference between the basic length and the actual length of the element. If you look into the strain formula, you will easily find this variable there. Strain in tension element will be expressed by the ratio between displacement and basic length. This is freaky good information, because if I will put all the data into Young modulus equation, suddenly! Out of nowhere! we will have only one (1) unknown. And... It is our displacement! WOW! This was super fast quick and simple. Good Job Me.

With the values from the first figure you can count that 10 m long steel bar with 10x10cm cross section loaded with 1 kN will deform... just 0.0047619m, it is u=5mm. It is time to see what FEM will propose to us. In almost every book about basic FEM you will find the bench of basic equations like this:

How you are building element stiffness matrix will be a separate post, but what is important now is that this stiffness matrix is just relevant for bar element. To explain what is what in this equations we have to introduce some nomenclature. 

Remember that the solution (in our case displacement) in the FEM is only precise in nodes. For the rest of the structure, we have to introduce approximation. The displacement and force vector will be than in size of 2. We have 2 nodes we have 2 displacements and 2 forces. Easy! If we will have 3 nodes we will have how many displacements? ... of course 3. If we will have 230213 nodes we will have 230213 displacements in the displacement vector. These calculations are of course true when we have 1D task with only displacements allowed according to the one direction. 

Ok, so we know now enough to solve our task. We have to just modified equilibrium equation to have unknowns on one side and rest on the other (right one). To do like this:

We have to inverse K matrix (element stiffness matrix), but wait a minute! If  I will just try to inverse K matrix something will go wrong: 

The determinant for my stiffness matrix is equal to zero?!?! Do you think, what I am thinking right now? Somebody is trying to sell me an old stinky fish for really big money... But we are engineers, we need to solve a problem, the mathematicians will come afterwards and figure it out the methods for this... We need to solve the equation system, so maybe the boundary conditions will help us. So before we inverse matrix, let's try to specify and apply the boundary conditions to K matrix:

   

Why we wiped out our first column and first raw? We know, that the displacement in the first node (u1) is equal to 0 (it is fixed), and we know that K matrix should be symmetrical. If you do not understand why K should be symmetrical, you can also ask yourself: why u1 should influence u2, if we know that u1 is fixed so displacement u1 for sure will be zero, null, nothing, 0? This is why we are cleaning raws and columns responsible for fixed nodes. Let's try to inverse matrix now:

Done! WOW! That was extremely hard! Now we can start to find displacements. It is more than simply, that we have to multiply inverse K matrix and force vector. The result is:
  u1 = 0
  u2 = ( P2∙l )/( E∙A )
Since l is the basic length of the element (length before loading) we can say that it is the same result as an analytical solution... AMAZING! But try to calm down for a second, because we  still have to learn how to program this thing:

The C# code (VS 2015 community) from the movie:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;

namespace Truss_FEM
{
    class Program
    {
        static void Main(string[] args)
        {
            // The heading of the code
            Console.WriteLine("FEM 1D program for bar elements");

            // The variables
            double emodulus = 1; // 210 GPa = 210 000 000 kN/m2
            double npoisson = 0; // no units
            double barlength = 10; // [m] meters
            double elementarea = 10; // [m2] meter x meter 

            // The bar stiffness matrix
            Matrix<double> ke = DenseMatrix.OfArray(new double[,] {
            {1,-1},
            {-1,1} });

            ke = emodulus * elementarea / barlength * ke;

            Console.WriteLine("Element stiffness matrix:");
            Console.Write(ke);

            // The boundary conditions (BC)
            ke[0, 0] = 1;
            ke[1, 0] = 0;
            ke[0, 1] = 0;

            Console.WriteLine("Element stiffness matrix(after applying BC):");
            Console.Write(ke);

            // Inverse stiffness matrix
            Matrix<double> ike = ke.Inverse() ;

            Console.WriteLine("Inverse element stiffness matrix:");
            Console.Write(ike);

            // The force vector
            Vector<double> Re = DenseVector.OfArray(new double[] {0 , 1});

            Console.WriteLine("Force vector:");
            Console.Write(Re);

            // The general equation Kr=R, deformation vector re
            Vector<double> re = ike*Re;

            Console.WriteLine("Deformation vector:");
            Console.WriteLine(re);

            // The elongated beam
            double newbarlength = barlength + re[1];

            Console.WriteLine("new bar length="+newbarlength);

            // The finish line of the code
            Console.ReadKey();
        }
    }
}