Basics of Finite Element Method — Direct Stiffness Method Part 2

In Part 1 of this series of articles on direct stiffness method, we covered formation of stiffness matrix for a 1D, 2D and 3D truss element. Further, we also introduced two separate approaches — a longer approach and a shorter approach towards direct stiffness method. In the longer approach, global structural stiffness matrix was derived by combining individual element stiffness matrix and then applying displacement compatibility and nodal force equilibrium equations. This approach was rather tedious one, involving several steps. But, using the observations gathered in the longer approach, a set of general rules for the direct assembly of global structural stiffness matrix were developed. These rules were then used to establish a so called shorter approach of direct stiffness method, which forms the basis for displacement based finite element formulations. In case if you have missed the previous article, it is strongly recommended that you first go through the previous article on the basics of direct stiffness method. A link to the part 1 is below —

Basics of Finite Element Method — Direct Stiffness Method Part 1

In part 1 you have seen the derivation of an element stiffness matrix for a simple 1D finite element i.e. a linear truss element which takes only axial loads. Further, by representing individual components of a truss like structure with these simple 1D truss finite elements, it was shown how a global stiffness matrix for the entire structure could be derived. Apart from an engineering interpretation of these element stiffness matrices as representing the mechanical properties of the underlying continuum, there’s also a few important mathematical aspects which helps in the analysis of the computational results.

To understand the mathematical properties and their physical interpretation of element and global structural stiffness matrix, the example problem from article part 1 is considered here too. In fact, we would directly takeover the expression for an element stiffness matrix k and global structural stiffness matrix K, obtained in article part 1 via symbolic computation using a SymPy code.

Fig. 1 Example problem taken from article part 1 — truss like structure composed of 1D linear finite elements

To start with, a stiffness matrix, whether for an element or for entire structure, is just like any other matrix in linear algebra, which linearly transform a given displacement vector d to a corresponding force vector F in a 3D Euclidean space. Analogously, if the stiffness matrix is invertible (non-singular) with non-zero determinant, then it’s inverse can be used to map a force vector F to a corresponding displacement vector d 3D Euclidean space.

Fig. 2 Interpretation of a stiffness matrix as a linear transformation

Let us consider the element and global structural stiffness matrices from the example in article part 1 for understanding following mathematical properties in detail.

Fig. 3 Formulation of a 1D spring/truss finite element

Stiffness matrix for the above linear truss element, as derived in article part 1 is rewritten below as

Fig. 4 Stiffness matrix of an individual truss finite element — k

Similarly, the global stiffness matrix for the entire truss structure as derived in article part 1 is also rewritten here for sake of simplicity. We would use these expressions of element and global structural stiffness matrices for our analysis of their mathematical properties in next few sections.

Fig. 5 Global stiffness matrix of the truss structure problem as derived in article part 1 — K

Stiffness Matrix is Real-Symmetric

The stiffness matrix of the entire structure from example in article part 1 shown above in figure no. 5 was assembled using a SymPy code using rules of assembly as derived from the direct stiffness method.

First observation which you could make is that the global structural stiffness matrix K is also a square matrix of size n. Where n is the total number of degrees of freedom in the problem. Since the example problem from article part 1 had total 8 global degrees of freedom, global structural stiffness matrix K is a square matrix with 8 rows and 8 columns. Further, one observes that the stiffness matrix comprises of only real values, without any complex entries. In fact, in the field of structural mechanics, we normally deal only with real matrices, i.e. when all the entries belong to a set of real numbers R.

Fig. 6 Stiffness matrix — square with real entries

Additionally, it can also be seen that the global structural stiffness matrix is symmetric in nature. Hence, K.transpose = K. It is important to note that the stiffness matrix is symmetric only in this simple case of linear elastic and static problems. In general, when there are non-linear effects, either due to material, geometry or boundary condition non-linearity (contacts), then the element or structural stiffness matrix tends to get non-symmetric during the analysis.

Fig. 7 Symmetry interpretation of a stiffness matrix

Hence, for this class of problems — i.e. linear elastic and static, both structural and element stiffness matrices are real-symmetric square matrices.

Eigenvectors and Eigenvalues of Stiffness Matrix

A lot of interesting insights can be obtained when one tries to perform an eigenvalue analysis of the stiffness matrix (element or structural). A pair of eigenvector and it’s corresponding eigenvalue of a stiffness matrix could be interpreted physically as a certain deformation mode and a corresponding stiffness resisting against that deformation mode. For a square matrix of size n, there would also be n eigenvalues and n eigenvectors.

Fig. 8 Eigenvalue problem of a stiffness matrix

A zero eigenvalue for a stiffness matrix means “zero stiffness” for the corresponding deformation mode (or eigenvector). For a simple case of 1D truss element, there are 3 zero eigenvalues, each corresponding to a deformation mode which can be classified as a rigid body mode, for which no stiffness in the element is generated. In other words, internal strain energy in element will be zero corresponding to these deformation modes. The fourth non-zero eigenvalue corresponds to an elongation of 1D element with equal amount at both nodes. An overview of all four eigenvalues for a 1D truss element is shown below

Fig. 9 Eigenvalues of 1D linear truss element

It is also possible to visualize the 3 rigid body modes and one deformation mode for the 1D truss element graphically as shown in the illustration below. There are two rigid body modes corresponds to pure translation of entire element in X and Y direction respectively. The third rigid body mode corresponds to a pure rotation about an axis perpendicular to the plane containing truss element. It is interesting to see that, all eigenvectors with zero eigenvalues represent nullspace of the stiffness matrix. In other words, no forces would be generated inside the element when rigid body deformation mode is enforced upon it.

Fig. 10 Physical interpretation of eigenvalues and eigenvectors of a 1D linear truss element

Similarly, when we perform an eigenvalue analysis of the global structural stiffness matrix for example from article part 1, as shown above in figure no. 5, we get following eigenvalues -

Fig. 11 Eigenvalues of global structural stiffness matrix

Since, the stiffness matrix corresponds to a structure where no boundary conditions were enforced, we see that the stiffness matrix has got three zero eigenvalues, each corresponding to a rigid body motion in 2D space, similar to an individual 1D truss element. However, there is also a fourth zero eigenvalue in this global structural stiffness matrix, which corresponds to an internal mechanism in the structure — the element 1 is free to rotate relative to the remaining triangular frame. In other words, there is a zero stiffness associated to rotation of element 1 about node 4, hence it contributes towards an additional zero eigenvalue for the entire structure.

Fig. 12 Summary of eigenvalues for example problem from article part 1

The figure below illustrates the graphical interpretation of eigenvectors corresponding to the zero eigenvalues of global structural stiffness matrix -

Fig. 13 Zero eigenvalue deformation modes

Further, it can be seen that both element and master stiffness matrices have zero determinant. Hence, they are not invertible in the current form, which means that they are singular. This is also associated with existence of zero eigenvalue. From linear algebra we know that if an eigenvalue of a matrix becomes zero, it’s determinant also becomes zero and hence matrix becomes singular or non-invertible. So, if any eigenvalue becomes zero for stiffness matrix, it would not be possible to invert it and hence no unique solution for displacements can be obtained.

Fig. 14 Singularity of stiffness matrix

However, generally, after applying displacement boundary conditions on element or master stiffness matrix, they don’t have any zero eigenvalues, and hence they become invertible. So, it would be possible to obtain a unique solution for the displacement.

But, this is true only for linear problems! For non-linear problems, even after application of displacement boundary conditions, the stiffness matrix can become singular, or it can suffer a loss of rank. There are many instability phenomenon where such a scenario might occur — like buckling analysis or snap through analysis. In these cases, eigenvalues becomes zero, even after accounting for displacement boundary condition. The corresponding deformation mode is then referred to as an instability mode, which is not a rigid body mode, but still results in a zero eigenvalue.

Positive Definite and Positive Semi-definite Stiffness Matrix

Till now we have concluded that the stiffness matrices are real valued and symmetric (well, at least in the case of linear problems). From linear algebra we know that any real-symmetric matrix has got real eigenvalues and orthogonal eigenvectors. Similarly, stiffness matrices have also got real eigenvalues which are greater than or equal to zero. Due to this reason, stiffness matrices are at least positive semi-definite. So, stiffness matrices prior to the application of boundary conditions are positive semi-definite and non-invertible, because they have at least one zero eigenvalue.

However, after application of displacement boundary conditions, stiffness matrices becomes positive definite with real-positive eigenvalues (except in case of instability problems). Positive definite matrices are invertible, non-singular with all eigenvalues, determinant and sub-determinants greater than zero.

The advantage of working with positive definite stiffness matrix is that it’s eigenvalues, which in turn represents “stiffness” for corresponding eigenvectors (deformation modes) are real and greater than zero. This is so desired because negative stiffness does not make any sense physically and zero stiffness represents either a rigid body mode or an instability mode. In the field of structural mechanics, where we work predominantly with stiffness matrices of structures, it is helpful to work with positive definite matrices because it results in physically sensible eigenvalues (stiffness)and eigenvectors (deformation modes).

Fig. 15 Stiffness matrix as a positive definite and positive semi-definite matrix

Conclusion

In this article part 2 of basics of finite element using direct stiffness method, we covered some mathematical properties and their physical interpretation for an element and a global structural stiffness matrix. We figured out that stiffness matrix has following properties -

• Real-symmetric
• Real eigenvalues and orthogonal eigenvectors
• At least Positive semi-definite

We learnt that the stiffness matrices are singular or non-invertible prior to the application of displacement boundary condition, which means they have at least one zero eigenvalue, which represents either a rigid body movement or an internal free mechanism inside the structure. However, there are cases of instability like buckling or snap through processes, where stiffness matrix becomes singular or non-invertible even after application of displacement boundary conditions. In such cases, there exists deformation modes (eigenvectors) which are not rigid body modes, still corresponding stiffness (eigenvalue) is zero.

All these ideas and concepts actually helps a lot in verifying the results of finite element simulations and it’s computer implementation. Apart from that, they are helpful in gaining a graphical perception of otherwise cryptic numbers and matrices.

I hope you might have learnt something new in this second part. In case you have any questions regarding content of this article or would like to suggest any modifications or corrections, do drop me a message. It would be great to hear your feedback to keep motivating myself to bring more interesting articles in future!

All the best for your learning and see you in next article!