domingo, 17 de julio de 2011


Finite Element Analysis is a mathematical tool very extended among engineers. However, after more than a year researching on the topic of computer simulation, where FEA plays such an important role, I haven't yet found a satisfactory explanation on how they really really work...
Hopefully by means of this post I get to clarify some of the topics, trying to remove excessive algebraic and mathematical verbosity that is so annoyingly present everywhere.

The main background of FEM is that of structural engineering in the 60s. Engineers are very practical people, so initially they devised a system which allowed them to set algebraic equations where the relations between different points or "nodes" of their structures could be set.
These algebraic relations were further proven to be of many different types, not only structural, so also thermal relations could be formulated and many others: all one needed to analyze was a proper discretizacion of the space in the form of a mesh with nodes to relate to each other.

For our particular case of  structural engineering (the one that matters for my PhD thesis),  I have tried to illustrate the procedure in its three main steps, so the main ideas come up in a graphic taste:

STEP 1: Discretization
In this basic yet classic example, I have chosen to divide my sample structure in only 5 elements where nodes are clearly identifiable in the meetings between beams. There are four nodes.
Degrees of Freedom
Each node will have 6DOF (Degrees Of Freedom): three for linear displacement on each axis (X,Y,Z) and three for rotational around each axis (X,Y,Z), because we are working on 3D. Many examples available provide the more "simple" 2D situation, but in my opinion this only complicates things further.

A sample structure (mesh) and its topology table

Once the nodes are located and there is a network of how they relate to each other, we can consider we have a mesh. In our example, the correlation is depicted in the table: The elements serve to "link" nodes to each other. The table establishes a "topology" for the nodes.

STEP 2: Element characterization and shape functions
In this step is where FEM formulation and literature get really really awkward and nasty. In fact, this is the core of everything and where FEM differenciates from other ways of solving PDEs.

Different types of Finite Elements
In reality, and despite its mathematical complexity (also unneccesary to be explained so much in detail in my opinion), what we are looking for is a way of characterizing the material properties of the element. For such purpose, the method requires that the behavior of those links among nodes obeys some formula. This formula is the actual Shape Function. In fact, the shape function can be any mathematical formula that helps us to interpolate what happens wherever there are no points to define the mesh. This "ghost" entity that appears between nodes is in fact the Finite Element. In practical terms, as engineers we are more interested in the implementation of the FEM, not so much on its formulation, so what is important to understand is that for different shape functions we obtain different element matrices.
2DOF Beam element matrix
3DOF Beam Element matrix
3DOF Timoshenko Beam element matrix
 6DOF Timoshenko Beam element matrix
Depending on the chosen formulation we have different degrees of interpolation and hence presumably higher or lower degrees of precision. Also depending on the chosen formulation we might have different ways of locating and relating our nodes to each other.

For our example we can choose any of the formulations provided in literature (above are the most common used in structural engineering). It is important to note the internal structure of these element matrices, which are symmetrical and clearly divided into parts each corresponding to the nodes that reside on the element's boundaries (2 nodes in the case of a beam - 4 quadrants in the matrix).

STEP 3: Matrix assembly and solution

Because the relations between nodes need to be accomplished all at the same time, we have to set all the equations in such a manner that they compose an algebraic system of equations. The matrix equation we want to solve (at least in statics) is as follows: 

[F]= [Kg]·[u]

Where [F] is the vector of applied external forces, [Kg] is the system's global stiffness matrix and [u] is the vector of displacements resulting from the application of the forces. The size of [Kg], [F] and [u] is that of the number of DOFs times the number of nodes, being [Kg] squared and [F] and [u] unidimensional.
In order to generate the Force vector, all we need to do is collect the applied forces (linear and moment forces) and sort them according to the index of the node they are applied to.
For the global stiffness matrix, it is necessary a bit more laborious procedure by means of which we iterate throughout each element's particular stiffness matrix. Out of each one of those, we get only the part that corresponds to the position of the node we are storing in the matrix, and add it to the possible concurrent data that comes from other elements. Warning: before entering in the global stiffness matrix, we must convert local coordinates to global coordinates!

To get the displacement vector, it is needed to first enforce the constraints and then solve the resulting algebraic system of equations. To do this there are two classical approaches: Penalty Method and Lagrange Multipliers method, but we wont enter into this here...

Global Matrix Assembly
Afterwards, all we need to do is use any available algebraic equation solver (LU decomposition is one of the most extended), and obtain the solution of the system.

18 comentarios:

  1. Hello sir. Im a student of mechanical engineering. Im having FEM this semester. I dint understand shape function so i googled for it and i found your blog.

    Why did you consider elements 4 and 5 when they actually dont exist ?

    "This "ghost" entity that appears between nodes is in fact the Finite Element" I did not understand what you meant by ghost entity.

    1. Hi Pruthvi,
      As you might have guessed from the tone of the articles published here I am trying to avoid the excess of terminology that -in my opinion- often blurs the physical notions and concepts. By "ghost entity" I mean an artificially introduced idea (in this case an interpolation function) that in reality -again in my opinion- is not actually supporting the physics underneath but some kind of mathematical "rigor".
      That being said, I emplace you to the next article I am currently preparing where I will try to make the Rayleigh-Ritz methodology understandable for people. This is often recognized as the embryo of the currently method known as FEM, and surely will help you out in your struggle with shape functions.

      Well, good luck and thanks for reading and the feedback :)!

      Regarding elements 4 and 5, they indeed do exist in the system pictured, though their Young's and inertia modulus (E and I) have different values from the others. Their corresponding nodal matrices will be assembled in step 3 and will affect to the global solution as any other element.

  2. Thank you for replying sir. Since the time i last commented my knowledge in FEM has increased a bit(infinitesimally) :) .

    How can elements 4 and 5 have youngs modulus and other quantities ? Its either that i did not understand or the shadow in the picture misses elements 4 and 5. Are nodes 1 and 3 connected by a solid ??

    1. Hi again Pruthvi,

      The correct answer is that of the renderer ignoring such thin bars (it was made with Sketchup...).
      Indeed for the effects of computation elements 4 and 5 have Youngs modulus like any other element. This is material dependent. What could be zeroed (or better, made very little=0,00000001) in the matrix equations are the values of its inertia modulus, who oppose to the rotational efforts (bending, torsion), and its area, who opposes to longitudinal deformations.

  3. Hi- just wanted to let you know that I'm a graduate student studying FEM for my research, and have been unable to find a description of shape functions this clear- most of the literature is far too math-intensive and verbose! Thank you very much for posting this. -Neha

    1. Thanks Neha!
      I had exactly the same problem..:)
      I am very happy it helps other people


  4. Hi,
    What is the reason behind the rule that "sum of shape functions should be one at a particular node and zero at all other nodes"?....
    Need explanation...please

  5. Hi Shashidar K,

    I assume you are following chapter 18 of the introductory course to FEM from the Colorado university program ( believe it is a formidable source for learning, though this chapter can be a bit confusing if one tackles it without having everything absolutely clear. I would suggest to go through chapters 11 and 12, where the notion of shape function is explained in an excellent manner for a very easy example...then visit chapter 15.2.4...

    Anyway, when they expose the four necessary conditions the "shape functions" have to meet they are simply talking about how the polynomial must be, preparing the reader for what is to come:
    "(A) Interpolation condition. Takes a unit value at node i, and is zero at all other nodes."

    In chapter 15.2.4 it is explained how this "shape function" i.e. the polynomial that is used to interpolate the behavior of the phenomenon (elasticity, electric flux, thermal flux,...) between the three points of a triangular surface is constructed.
    In English, basically the requirement says that in the nodes 100% of the behavior must be met, slowly reducing towards the other corners until becoming zero.
    Using the example of our beam elements in this article, imagine how between the two nodes communicated by our beam element the tension is biggest at the node where a force would be applied and slowly decreases as we are getting further from it, becoming zero in the other end.

    I hope this will be of help!

    Thanks for reading :)

  6. Great explanation :)
    I think you eliminating the maths rigour really made it more meaningful.
    Just one question, What fills the element matrices? E.g. for element 1 for the N1, N2 * N1, N2 matix what are the actual numbers in the matrix.
    Also, is putting in the different matrices values in the global matrix not changing its coordiantes , or is this an extra step?

    1. Hi,

      The element matrices are filled in with the values AE/L, 12EI/L³, etc...which are obtained from the derivation of the shape functions for each degree of freedom (rotation in x of first node, deflection in y, etc).
      These values are the relation (generally an interpolation polynomial) between the force or moment applied in one node of the beam and the displacement or rotation they create through the whole beam.
      For different polynomial types (hermitian, Lagrangian, etc) there is different approaches on how to fill in the element's matrix. The ones that I provide in the article are the most common for beam elements. After derivation, we can obtain the effect IN the nodes that define the beam element, omitting the rest of the beam.

      Regarding your second question, before assembly it is necessary to change the coordinate system. It is a common error not to multiply the local element's stiffness matrix by the transformation matrix. There are different algorithms to create this transformation matrices.

      Thanks for reading! :)

  7. Respected Sir, I am a researcher and trying to model with Finite Element Method. Please tell me the formulation of the rotor disk element.

  8. Oho...dear PMBG Adasque..I am afraid the question you are making is far beyond my capabilities...My research is oriented towards buildings and structural dynamics. Rotor blades and turbines belong to the discipline of aeronautics which unfortunately I am not so familiar with (though I never discard it...).
    I really wish I would be able to help you, but all I can do is wish you luck in the research.
    Maybe the project QBlade of the Institute of Fluid Dynamics and Technical Acoustics of the TU Berlin, which is Open Source, can be of help for you (

    So..that being said,

    Good luck!

  9. Hello Sir,

    I want more information about shape functions like does selection of shape functions influence the key characteristics of an element
    2.the steps involved in evaluating the shape function required for carrying out structural finite element analysis
    3.why natural co-ordinates are used in this methodology.

    where can i get all these information in detail?


  10. Hello Meet Bhatt,

    I would refer you to the Introduction to Finite Element Method of the
    Department of Aerospace Engineering Sciences of the University of Colorado at Boulder (, chapters 16, 17 and 18.

    I believe the precise answer to your questions should be in those documents.
    Bear in mind that the shape function is nothing but a polynomial (or set of them) used to approximate the expected behavior of the phenomenon (elastic deflection of a beam, heat transfer through a surface or volume, etc), "connecting" different nodes whose coordinates are degrees of freedom in a matrix system of equations. These coordinates can be referred to a global (0,0) point or to a local (natural) initial point of the element. This simplifies the process of coding and programming very much.
    Also, for structural analysis, most of the time we just pass with beams and plane elements (maybe sometimes also cables), so in general the more sophisticated elements are useful in other disciplines (mechanical engineering, robotics,..).

    I hope you enjoy the reading


  11. I really appreciate this blog.
    What exactly is the purpose of adding an "internal node" at the center of, say, a quadrilateral element? Texts say that the 3-node element is "too stiff in bending," which is just plain vague. I've seen it depicted as four triangular elements arranged together to form a quadrilateral, but how does this fix the stiffness problem?

  12. hi sir,
    could you explain little more briefly about shape function and stiffness matrix!

  13. I still dont see where the shape functions come into play with the calculations?