Basics of Finite Element Method — Direct Stiffness Method Part 1
The term FEM (Finite Element Method) has gained a lot of traction in past few decades, specially in the field of virtual product development which involves creating mathematical models of a real system and using numerical methods to analyse its response for a variety of real load-case scenarios. Whether you are from a mechanical, aeronautical, civil, environmental or a nuclear engineering background and if you are involved in simulation and analysis of relevant processes and structures, you might have heard of FEM as a mathematical tool to obtain an approximate solution for the underlying governing differential equations, for which an exact analytical solution is otherwise not possible.
The aim of this series of articles is to impart the reader a basic understanding of what FEM is all about, and how it could be used to solve problems in the field of structural mechanics in particular. If you are a beginner in this field and want to start learning about FEM, this series will help you bring different basic aspects of FEM into perspective and build a solid foundation for the same, which will be helpful when you take more advanced courses related to FEM in particular and structural mechanics in general. But, if you are already an advanced user of FEM for your Job or Academic Research, this series could still be of some help for you to connect the dots. These articles will cover some aspects which you might have missed earlier, or simply provide a review of the basic concepts and theories you might have already learned in your studies. For these series of articles, I am following material from my master course at university of Stuttgart — COMMAS.
To bring things into perspective, I have organised the first article into following subsections:
- Theory of Direct Stiffness Method (DSM)
- Derivation of 1D Truss Element Stiffness Equation
- Analysis of 1D Truss Stiffness Matrix
- Extension of 1D Truss Stiffness Matrix to 2D & 3D Truss Stiffness Matrix
- Transformation of Element Equations from Local to Global CS
- Direct Stiffness Method using Example — Longer Approach
- Direct Stiffness Method using Example — Shorter/Direct Approach
- Solution of Structural Stiffness Equation
- Further Motivation
So, lets dig into each of these topics right away!
Theory of Direct Stiffness Method (DSM)
If you have studied or used FEM before, then you might have heard about the term Direct Stiffness Method (DSM). The concept of Direct Stiffness Method is used to introduce basics of Finite Element Method. So, first, let us try to understand what direct stiffness method is and what its advantages are in context of the finite element method (FEM) —
Direct stiffness method is the most primitive way to understand the basic concepts of FEM in general, and specially the displacement based Finite Elements which are most commonly used in solving problems related to the field of structural mechanics. In other words, direct stiffness method is the most simple (or, direct) form of the finite element method which is used to solve discrete problems and is considered as a good starting point to understand underlying idea of finite element method. Direct stiffness method can be used to analyse structures which are composed of discrete components and where each structural component can itself be considered as an ‘element’ (in analogy to an ‘element’ in the finite element method).
The concept of ‘discretization’ is however trivial in direct stiffness method, because the structure is assumed to be already discretized and the underlying properties (e.g. stiffness matrix) of each individual element is assumed to be known in advance, that too in an exact form.
Hence, in direct stiffness method, there is no discretization of the structure in a real sense and also there is no approximation of the underlying material properties. Using the already known structural properties (stiffness matrix) we are just assembling each component to get the structural property for the entire structure, using certain rules of displacement compatibility and force equilibrium. You will understand this idea more clearly when we go through all the steps of direct stiffness method using an example in the later sections of this article.
In direct stiffness method, entire structure is broken down into sub-components (think of “elements”) and for each component the displacement values at its two ends (nodes) is considered as the primary unknown variable of the problem. This is also a motivation for the displacement based finite element method, where the primary unknown of the problem is displacement values at finite element nodes. We represent properties of underlying continuum of each sub-component or element via a so called ‘stiffness matrix’. A stiffness matrix basically represents the mechanical properties of the underlying structure. It relates the nodal displacements to nodal forces via a linear transformation. Using the direct stiffness method, the rules for assembly of global/master structural stiffness matrix and load vectors are established, which could also be used in a general displacement based finite element method applied to more complicated domains, as we will see in the later articles of this series.
But, before going into detailed steps involved in direct stiffness method, we first need to establish the fundamental building blocks, i.e. stiffness matrix and load vector for a simple ‘line’ element. In the simplest form, a line element could be a linear elastic spring with stiffness k, which can take loads only in the axial direction. Alternatively, a simple line element could also be just a simple truss/bar element in 1D, 2D and 3D with a local coordinate system attached to the line element itself. Such a simple truss element can also take only axial loads, they are assumed to have a certain length (L), cross section property (Cross-sectional area, A) and material property (E, Young’s Modulus).
In other words, the mechanical properties of such a truss element could also be represented by its equivalent stiffness value k. This stiffness term combines truss’ geometrical and material properties in a single expression. From any undergrad solid mechanics text book you can find that this stiffness term k for a continuous 1D truss can be expressed as k = AE/L. If, however, the structure consist of discrete ‘springs’ , then the mechanical properties of underlying continuum could be directly represented via a discrete spring stiffness k.
Derivation of 1D Truss Element Stiffness Equation
The 1D finite element under consideration has only one displacement degree of freedom per node, with one node at each end. Hence, a total of 2 degrees of freedom per element. The element represents a 1D truss or discrete spring member which is assumed to have a stiffness k. If the underlying continuum of the truss member is assumed to be linear elastic, and if the element undergoes only small deformations, then it could be easily derived that stiffness k = AE/L. Each end of this element is called as a ‘node’ which experiences a load in the axial direction, which in this case is also the direction of the displacement at each node.
To derive the above equations, two important conditions were made to be fulfilled: Force equilibrium at each node and material response (constitutive equation) of 1D Truss. The force equilibrium states that the resultant force at each node should be zero, i.e. net external force = net internal force at each node. Secondly, we consider that the constitutive relation for the linear elastic spring or truss, which states that the internal force in the element is proportional to the net elongation of the element, where the proportionality constant is the stiffness of the underlying continuum.
After minor rearrangement of the nodal equilibrium equations, it is possible to represent them in terms of a matrix notation Ku = F. Here, the elemental stiffness matrix could be interpreted as a linear transformation matrix, which linearly transforms nodal displacement vector, u of an element onto corresponding nodal force vector, F. Hence, stiffness matrix of an element represents it’s mechanical properties and relates nodal displacements to nodal forces. Stiffness matrix plays a central role in identifying mechanical response on a structure.
Analysis of a 1D Truss Stiffness Matrix
It can be observed from the stiffness matrix of an individual element, that it is symmetric in nature, i.e. k_12 = k_21. The attribute that stiffness matrix is symmetric comes from the Maxwell’s Reciprocal Theorem which states that for any linear elastic body, displacement produced at any point A due to certain load applied at point B should be equal to displacement produced at point B when same load is applied at Point A.
Since we have assumed the truss member to be linear elastic, the Maxwell’s Reciprocal Theorem applies here and hence the stiffness matrix is symmetric. This, however, is not always true! The stiffness matric tends to get un-symmetric when material behaves in-elastically or has local instability, for example in problems involving damage and failure.
In element stiffness matrix, any general stiffness component k_ij can be thought off as a stiffness term (or a material parameter) which relates the displacement in direction of degree of freedom j with the force in direction of degree of freedom i. In other words, the amount of force required in the direction of d.o.f i to produce a unit displacement in the direction of d.o.f j is represented as k_ij. Thus, according to the Maxwell’s Reciprocal Theorem, k_ij = k_ji must be true for 1D linear elastic truss problem under consideration.
In a more general sense, the stiffness term k_ij from a global structural stiffness matrix K, can also be interpreted as a ‘coupling’ term which ‘couples’ global DOF i with global DOF j. Here, both i and j varies from 1 till the total number of degrees of freedom in an element or in the entire structure. In other words, if there’s a direct load transfer from DOF i to j or vice-versa, then the corresponding stiffness term in stiffness matrix is given as k_ij or k_ji. However, if DOFs i and j are independent from each other and there is no direct load transfer between them, then the corresponding stiffness or coupling term would also be zero. This concept would be more clear when a global stiffness matrix is derived using the direct stiffness method.
Further, it can also be observed that the summation of all rows and columns of this stiffness matrix results to 0. This has to be true because if the element suffers a rigid body movement, i.e. u1 = u2, then there should not be any force generated inside the spring element and hence nodal force vector would also be zero, which leads to the condition that sum of all rows of stiffness matrix is zero. And because the stiffness matrix is symmetric in this case, it follows that sum of all column entries are also zero.
Extension of 1D Truss Stiffness Matrix to 2D & 3D Truss Stiffness Matrix
Till now we have derived the element stiffness matrix for 1D truss element with single degree of freedom at each node. Now, the stiffness matrix for 1D Truss bar with one degree of freedom per node can be extended one step further to also represent a similar 1D Truss bar but with two degrees of freedom per node— one longitudinal (in axial direction) and other transverse displacement at each node. Because the truss element can take loads only in the longitudinal directions and have no resistance for transverse loads, the corresponding stiffness terms are also zero.
From the above stiffness matrix it could also be interpreted that longitudinal force at node 1 i.e. f1 is not coupled to transverse displacement at node 1 as well as at node 2, and transverse load at node 1 i.e. f2 is not coupled to longitudinal displacement at node 1 and 2. Due to this reason the corresponding stiffness terms/coupling terms in the stiffness matrix are also zero. It means that no amount of longitudinal force f1 at node 1 could produce non-zero transverse displacement u2 and u4. The same analogy also holds true for node 2.
Similar expression for an element stiffness equation can also be derived for a linear truss element in 3D with three translational degrees of freedom per node. Hence, resulting stiffness matrix would be a square matrix of size 6, because now there would be 6 DOF in total.
Transformation of Element Equations from Local to Global CS
We have seen till now the element stiffness equation in local element coordinate system, which was essentially parallel to the global coordinate system. In this configuration, line of action or longitudinal axis was parallel to global X-Axis. But, following a basic coordinate transformation, element stiffness equation for an element in arbitrary orientation in space can also be derived. To derive the coordinate transformation rules, lets consider a linear truss element in 2D space which is oriented at an angle theta with respect to the global coordinate system.
As it was done before, the conditions of displacement compatibility and equilibrium of nodal forces can be established here also to relate global and local values of nodal displacements and nodal forces. First, expressing the global nodal displacements in term of its components in local elemental coordinate system leads us to following relations:
Above relations could be rearranged in a Vector and Matrix notation as follows:
Similarly, a transformation relation could also be written for nodal force vector:
Combining transformation relations for both nodal displacements and nodal forces together with the element stiffness equation in local coordinate system, the element stiffness equation in global coordinate system could be derived as follows:
After doing a few basic matrix operations, the global element stiffness equation can be expressed in terms of cosines and sines of element orientation angle theta
Direct Stiffness Method using Example — Longer Approach
Till now we have established the basic concepts needed to derive stiffness matrix for a simple displacement based finite element i.e. a 1D, 2D or 3D truss element, we can now work towards the establishment of direct stiffness method using an example problem to explain all the steps involved.
Consider the sample structure as shown in the figure below, where each of its member is assumed to be of same material (hence, same strength of underlying continuum) and have similar cross-section properties. The left side of the structure is constrained (no displacement or rotations allowed) and a vertical point load is applied to the top right corner of the structure. The goal is to find the deformed configuration of the structure and forces in each individual structural member. I have taken this problem from the course on Computational Structural Mechanics taken by Prof. Bischoff in Fall 2013 during COMMAS course.
Now, the direct stiffness method involves following steps to establish a global stiffness matrix, which will describe the behaviour of entire system and not just an individual element.
Step 1: First step is breakdown of the entire structure into individual sub-components or ‘finite elements’. Each element is assigned a unique ID. Next, identify the ends of each finite element as a ‘node’, also with a unique Node ID. Since, the primary unknown variable here is displacement and as we are using displacement based finite element, the unknown degrees of freedom here would be displacements at nodes in longitudinal and transverse direction. Assign unknown global nodal displacements for the entire structure. For the sake of ease in computations, start degrees of freedom from left hand side where the structure is fixed and where the global displacement values are already known (=0). We also consider external forces acting on the structure on nodes in global coordinate system. So, in total for current problem we have 16 variables involved, i.e. 2 displacement and 2 forces per node and total 4 nodes in total. Out of these 16 variables, 8 are already known to us. Hence, the total number of unknown variables remaining in this system are 8.
Step 2: Next, create a topology matrix which would establish a relationship between elements and their respective global degrees of freedom. Each column of this topology matrix represents an element and each row represent the global degrees of freedom associated with that element. In practical applications, number of elements are much larger as compared to number of degree of freedom per element. Hence, the topology matrix has got large number of columns. But, what we have here is a special case in which number of elements and number of degrees of freedom per element are both 4.
Similar to the topology matrix, it’s also helpful to create a vector for element orientation in global CS with angles of each element w.r.t global CS and a separate vector for individual element stiffness property. These vectors are helpful in computer implementation of the direct stiffness matrix.
Step 3: In this step, we use the knowledge which we acquired in previous sections to write stiffness equations for each element in the local elemental CS, then transform them into global CS using a unique transformation matrix for each element as derived earlier. In the figure below, element stiffness equations for element #1 is derived in local CS. Since, the local element CS aligns with the global CS (i.e. theta = 0), no transformation is needed for this element.
It is also possible to express above matrix equation in an “expanded” form, accounting for the global degrees of freedom of the entire structure and their corresponding stiffness terms. In the figure below, stiffness equations for element #1 are expressed in terms of global degrees of freedom.
As you can notice in the above expanded form, only the contributions from element #1 towards the global structural stiffness matrix is considered. All other terms are currently zero as we are considering only element #1 here.
Using the similar approach, stiffness equations for all remaining elements with respect to global CS can be established. Further, they could also be expanded to the global degrees of freedom system, to express their contributions towards global structural stiffness matrix.
So till now, we have established individual element stiffness equations in the global CS. We even established an expanded version of the stiffness equation, which represents contribution of individual finite element towards the global structural matrix.
Step 4 (Longer Approach): After establishing individual element equations, now in this important step, we combine these equations to create a global structural stiffness matrix. Since, the entire structure is assumed to be built up of individual finite elements with their own stiffness equation, adding stiffness equations of all finite elements will represent the response of the entire structure.
Let us first start with the addition of nodal force vectors. Since, the nodal force vector for each element is already represented in the global coordinate system in previous steps, we can add the force vectors without having to perform any further coordinate transformation.
The above addition of vectors looks intimidating at first, but we will simplify it further based upon certain conditions. Firstly, let us consider the most fundamental equilibrium condition in finite element method i.e. the equilibrium of forces at each finite element node. The condition states that
The net force on any structural node would be zero. The net force incorporates both externally applied forces and internal forces in each element. In simpler terms, net external force should be in equilibrium with net internal force on any finite element node.
To understand the implications of the above statement, let us try to establish equilibrium of forces at each node of the problem at hand. As described in the picture below, at each node an equilibrium is established between external nodal forces F and internal forces P in each finite element. For the current problem, we already know that the external forces on nodes 1 and 2 are reaction forces from the constraints and they are treated as unknowns which needs to be solved for. Further, it is also known that there is no external force acting at node 3 and only a vertical load acting at node 4. All of these information is considered to simplify the nodal force balance equations shown below.
After applying the relations established in the picture above, the resultant global nodal force vector can be represented as follows
Next, we are interested to obtain a system of equation which describes the mechanical behaviour of the entire structure, i.e. we want to identify a global stiffness matrix K, which linearly transforms the global nodal displacement vector D to the global nodal force vector F as follows
{F} = [K]{D}
Since we have broken down entire structure into multiple sub-components or elements, global equation of motion could be written as a summation of individual element equations. The left part of the equation, i.e. resultant global nodal force vector has already been obtained in previous section. It is also possible to add stiffness terms on the right hand side for all individual elements as below
Now, the condition of displacement compatibility is used to combine all the element global displacement vectors into a single master global displacement vector for the entire structure. Compatibility of displacement states that
Two or more elements connected at a common node will have common degrees of freedom at that node. In other words, all elements connected to a common node will have same displacement value at that node or the common node at element will not move apart.
The global nodal degrees of freedom for the entire structure can be expressed as
Then, using the condition of compatibility at above nodes we obtain the following relations for nodal displacements
Using the above relations from compatibility condition, the addition of element stiffness equations can be rearranged as follows
Now, the individual element stiffness matrix in global coordinate system, as derived previously could be inserted into the above summation equation to obtain a master structural stiffness matrix which relates global nodal forces to global nodal displacements of the entire structure
Finally, the global structure equation of motion can be written as below, which is essentially a linear system of equations which be solved either be linear solver or iterative solver, depending upon the size of the problem.
It is interesting to note that in the master stiffness matrix, there are non-zero stiffness terms at only those positions where the corresponding degrees of freedom are directly coupled with each other, or, if there is a direct load transfer between those degrees of freedom. So, terms k_13 = k_14 = 0, because degree of freedom 1 is not directly coupled in any way to degree of freedoms 3 and 4, or in other words there’s no direct load transfer between nodes 1 and 2. Similar observation could also be made for nodes 2 and 3.
Further, it could also be observed, that the nodes where multiple elements are connected together, gets stiffness contributions from all connected elements in the concerned degrees of freedom. For example, in the direction of global degree of freedom 7, stiffness contributions are coming from elements 1 (k1) and 2 (k2/2) only. Hence, k_77 = k1 + k2/2. Similarly, the global degree of freedom 8 is getting stiffness contributions from elements 2 (k2/2) and 4 (k4) only, hence, k_88 = k2/2 + k4.
Direct Stiffness Method using Example — Shorter Approach
As you could see, it was a rather tedious process to convert each elemental stiffness equations first from local to global CS, expanding them to incorporate other global degrees of freedom and finally using compatibility of displacement and nodal force equilibrium to obtain master stiffness matrix. The final result was very intuitive form of master stiffness matrix which has non-zero stiffness terms for only those degrees of freedom which are coupled and zero for those which are not coupled.
Based upon the steps followed above and the results obtained, following rules for direct assembly of a master stiffness matrix can be obtained:
Rule 1: All the diagonal terms, k_ii consists of sum of direct stiffness of all the elements contributing directly to the degree of freedom i.
Rule 2: All the off-diagonal terms, k_ij consists of sum of indirect stiffness of all the elements connecting degree of freedom i to j. Add a negative sign to the resultant indirect stiffness.
Rule 3: Put a zero for all the other degrees of freedom which do not interact with each other, or which are not coupled to each other in any way.
Using these simple rules, it is possible to assemble a master stiffness matrix for a structure with any level of difficulty.
The global structural matrix assembled above for the entire structure is same as the one obtained in previous section using more tedious process of fulfilling displacement and force compatibility conditions. Hence, it is proven here that it is possible to use rules of assembly of stiffness matrix to obtain global structural stiffness matrix under direct stiffness method framework.
Solution of Structural Stiffness Equation
Once the global structural stiffness matrix is assembled, the resulting system of equations could be solved for the system unknowns — nodal displacements and reaction forces. However, if one tries to find a solution to the system of equations directly, one will face the issue that there are no unique solutions to the problem. The reason being that the assembled global structural stiffness matrix is not invertible, or it’s a singular matrix as it’s determinant is zero.
In order to find a unique solution to the problem, we apply boundary conditions, or Dirichlet Boundary Condition to be specific. Using this boundary condition, displacements at nodes are specified and system of equations is reduced to a smaller system of equations with invertible stiffness matrix.
The reduced system of equations is further simplified by inserting the expressions for individual stiffness terms
Further, the system of equations could be solved by using direct or iterative solvers for linear system of equations to obtain unknown nodal displacement values
The obtained nodal displacement values could be inserted back into the original system of equation to obtain unknown reaction forces
Finally, the unknown reaction forces at nodes are
Thus, the nodal displacements and reaction forces on the structure are obtained as follows
In the subsequent steps, which are left to the reader to try themselves, this global nodal displacement vector can be used to obtain nodal displacement of each element in local CS and thus forces in each element could be obtained using constitutive relation for each element.
Summary
In this article we covered following topics —
- Motivation behind direct stiffness method
- Derivation of stiffness matrix and equations for 1D, 2D and 3D linear truss element
- Analysis of truss element stiffness matrix
- Transformation of equations from local to global coordinate system using orthogonal transformation matrix Q
- Explanation of direct stiffness method using an example — longer approach using compatibility and equilibrium equations
- Establishment of rules from direct stiffness method and application to an example problem — shorter method
- Application of boundary conditions and solution to the linear system of equation
Further Topics
The purpose of this article was to introduce beginners in the field of finite element method to the domain of direct stiffness method, which forms the basis for displacement based finite element methods. We have reviewed the process of obtaining truss element stiffness matrix, combining individual element stiffness matrix using certain rules to obtain stiffness matrix for the entire structure. However, there are still many things which needs to be discussed in further detail, as far as direct stiffness method is concerned. In the next article, I would like to emphasize on mathematical interpretation of an elemental and global structural stiffness matrix, which provides many useful insights into the structural behaviour. Some of these aspects are
- Eigenvectors and eigenvalues of stiffness matrix
- Positive definite and positive semi-definiteness of stiffness matrix
- Singularity of stiffness matrix — structural stability and instability
- Rigid body modes and internal mechanisms in a structure
I hope the readers might have got a good introduction of the topic, or a decent review of the topic for those who are already familiar with the field of finite element method. In case you have any comments, suggestions, corrections regarding this article or the filed of finite element in general, do leave a comment or send me a message. I would love to hear back from you guys!
For a deeper understanding of direct stiffness method and mathematical interpretation of element/structural stiffness matrices, don’t forget to read the Part II of this series of article!
