Applications of Finite Element Methods (FEM)

In this article, we will explore various applications of FEM in the field of Mechanical Engineering.

Structural Analysis

In structural analysis, FEM is used to compute displacement, strain and stresses in the structure subject to external loads. A structural problem can categories based on structural geometry, material properties used and involvement of time.

• 1D, 2D & 3D

Based on the geometry, there can be 1D, 2D and 3D structural analysis problems. All structures are three dimensional in real world. However, we model them as 1D if its 1 dimension is relatively larger than other two dimension. For example, bar or truss structures are often modelled as 1D structures. Similarly, the plates and shell structures subjected to plane stress and plane strain conditions can be conveniently model as 2D structures. When, all three dimensions of a structure are comparable to each other, then, we model them as 3D structures. Also, we investigates the complex 3D stress or strain state in a structure model by modelling it as 3D structural problem. The computational cost increases from 1D to 3D structural analysis.

• Linear vs Nonlinear

Based on the material properties used, a structural analysis can be considered as linear or nonlinear one. If the material properties are within the proportional limit, it will be considered as a linear structural problem or nonlinear otherwise. Generally, we design our structure to withstand loads below its yield strength. In other words, such structures are designed to operate within the proportional limit of the materials, consequently, requiring linear structural analysis.

• Static vs Dynamic

A structural problem can be static if the state variables do not change with time, otherwise, it will be a dynamic problem. For static structure problems, only spatial discretization of the domain is required, however, for solving the dynamic problems, there is need of both spatial and temporal discretization. The temporal (time) discretization is needed to evaluate how the response of structure evolves over the passage of time. The dynamic problems are frequently solved using direct integration methods e.g. central difference method, Newmark method etc. Another method for solving linear dynamic problems is the modal method. It should be noted that the dynamic problems involves the terms of inertial force, damping force, elastic force and applied force. In case of static problems, the inertial and damping force terms become zero, thus, the static problems becomes a special case of dynamic problems that involve elastic force and applied force terms only.

For dynamic problems, the global stiffness matrix is obtained by combining the element stiffness matrices. In a similar fashion, the global mass matrix is built from the element mass matrices. However, the construction of global damping matrix is done differently. Rather than building it, we often express the global damping matrix in terms of global stiffness matrix and global mass matrix. This is one way of driving the global damping matrix, and it is named ‘proportional damping model’. It is frequently used because it is convenient as compared to other damping models (i.e. viscous damping model etc.).

Here is an example of how to drive the global mass matrix from the element mass matrix shown below. Similarly, how to express the global damping matrix in terms of global stiffness matrix and global mass matrix and how to find out the relevant constants from the vibration graph representing the relationship between the excitation frequency and the damping ratio of the structure.

Modal Analysis

This analysis is carried out to compute the natural frequencies and the mode shapes of the structure, which plays very important role to design the structures. If the frequency of an external source becomes close or equal the natural frequency of a structure, then, resonance can occur. The resonance phenomenon causes high amplitude vibration that may lead to catastrophic failure.

The modal analysis can be performed by setting assumptions (i.e. free vibration and zero damping) on the equation of motion. Free vibration means that we want computed the natural frequencies of structure by neglecting the excited force term. The damping term is set zero to conveniently solve the problem using the special mathematics technique i.e. Eigenvalue problem. It should be noted the damping force just reduces the amplitude of vibration and effects on the structure frequency can be neglected.

Thermal Analysis

It is used to compute the heat transfer problems. Heat transfer can be through conduction, convection and/or radiation. The thermal problems are solved in a similar fashion as that of the structure problems. For example, for 1D heat conduction problem can be expressed by Fourier equations, which then can be solved using FEM by building a set of simultaneous algebraic equations expressed in matrix form as shown below. The thermal problems can be of steady-state or dynamic nature. For steady-state problems, only spatial discretization is needed, however, for dynamic problems, both spatial as well as temporal discretization is carried out.

Coupled Problem Analysis

All physical problems are coupled in nature. We often decoupled or simplify those using appropriate assumptions so that we can conveniently solve them. However, when and where required, the FEM can also be applied to solve coupled problems. Coupled problems are generally complicated and require relative more computational resource. For example, coupled structural-thermal analysis is often carried to analyze the structures subjected to mechanical as well as temperature loads. For this kind of coupled problem requires both mechanical and thermal properties of the structure materials as inputs. Similarly, the FEM can also be deployed to solve fluid-structure interaction (FSI) problems. Investigation of the interaction between the structural dynamics and fluid mechanics at the interface is not only intricate to solve but also consumes a lot of computation resource. 

The coupled problem, in nutshell, require multidisciplinary knowledge and their solution can be challenging due to the presence of strong nonlinearity.

Developing Finite Element Mathematical Model (Formal Approach) - Part 2

In this Part-2 of previous article, we will be covering new concepts of Jacobian matrix and Gauss-Legendre Quadrature numerical integration method.

In Part-1, we learn how to derive the element mathematical model starting from the choice of shape functions expressed in the natural coordinate system (CS). In case of structure problems, we apply the minimization of potential energy principle to drive the element stiffness matrix presented in the integral form. For 1D bar element, it was easy to solve this integral manually, however, for complex problems we often make use of numerical integration methods.

One of the various ways to solve the integral form of the element stiffness matrix is to deploy Gauss-Legendre numerical integration method, as discuss in the following:

Gauss-Legendre (GL) numerical integration

GL numerical integration method is accurate and widely used in FEM numerical computation. The general form of this GL method is shown below. In this method, the integral form of a function is computed by summing up all terms, each obtained by multiplying a weight with function value, estimated at various integration points within the domain. Since, the method requires computation at few pre-determined integration points within the domain, rather than the integration over the entire domain, therefore, its computation is fast. The accuracy is dependent on the selected integration points selected from a Table. The more are number of integration points, higher will be the accuracy of the integration results.

It should be noted that the GL method can be applied to the functions spanning over the domain ranging [-1, 1] only. In other words, we need to express every function in the natural CS which are normalized, in order to numerical integrate them using this GL method.

This can further explained with the following example where a comparison between manual integration and GL numerical integration methods. In case of GL method, 02 integration points are selected with in the normalized domain [1, -1], but still the computed results are very close to that obtained from manual integration.

Jacobian Matrix

The Jacobin matrix can be considered as a transformation matrix which translate the given information between the physical spaces expressed in local/element coordinate system, and the natural space which is normalized space expressed in natural coordinate system. Below Figure below shows the Jacobian matrix linking the physical and natural spaces for a 2D element.

In FEM, the use of Jacobian matrix is linked the GL integration method. Whenever, we need to integrate a function expressed in local/element coordinate system, first we need to transform it in the natural CS so that we can make use of GL integration method. For transformation between the local and natural CSs, the Jacobin matrix comes in action. It should be well noted that for complex engineering problems, the element geometry in local CS is often irregular which makes it difficult to perform the integration in the local CS. That is the reason, we map the physical geometry space to an abstract/mathematical space defined using natural CS which is not only normalized but also regular, making it convenient to perform numerical integration. Here for 2D element, we can see how a Jacobian matrix looks like.

Now further expand the integration of stiffness matrix of a bar element from 1D to 3D. It will look like this:

Let’s summarize the all concepts related to developing of element mathematical model. We begin with shape functions using natural CS to express the state variables with the element domain. The natural CS is used because we need to perform Gauss-Legendre numerical integration in the later staged. After executing various mathematical steps, and applying the minimization of potential energy principle, we derive an integral form of the element stiffness matrix which we developed previously using simplified direct method. The Jacobian matrix comes in action, when we need to map state variables in physical space to the natural space. This transformation makes it convenient for us to evaluate state variable which might be challenging to compute otherwise in the physical space due to its complex geometry. For details, see below Figure.

One last thing that we must keep in mind, that the element stiffness matrix is expressed in the local/element coordinates. The evaluation of the integral form of the stiffness matrix is done in the natural coordinate system is carried out for computation convenience only. For mapping between local/element CS and natural CS is carried out using Jacobian matrix. Once, the element stiffness matrix is computed in the local CS, then, it needs to be transformed according global/system CS, before the element assembly process begins. For transformation from local CS to global CS, we require a transformation matrix which we already covered in our previous articles. Here is the complete table describing the summary of transformations/mappings.

Developing Finite Element Mathematical Model (Formal Approach) - Part 1

 In this article, we will be covering how to drive the element mathematical model using formal approach that involves the concepts of shape functions, natural coordinates, minimization of potential energy and numerical integration using Gauss-Legendre quadrature method. The article has 2 parts. Part-1 covers the how to drive element stiffness matrix. Part-2, coming next, will be focused on how to make use of numerical integration to evaluate the entries of element stiffness matrix.

Let’s dig out the formal approach by considering a simple 1D bar element with 2 displacement dofs, as show in Figure 1. There are 03 types of coordinate systems (CS) to define the element mathematical model: global CS, local/element CS and natural CS. The global CS is used for element assembly process. The local CS is used to define the element stiffness matrix, and the natural CS is used to conveniently solve for the entries of element stiffness matrix etc. The concept of using natural CS will be further explained in the following text.

In case of formal approach, we begin with choosing with a displacement function that is expressed in terms of nodal values. Such functions are named as shape functions as they define the distribution of displacement field within the element domain by interpolating nodal values. The number of shape functions depends on the number of element nodes and associated dofs. For example, 2 shape functions are needed for 1D bar element with 2 nodes, each having 1 displacement dof.

Pay attention to the fact that:

·      Shape function are expressed in terms of natural coordinate system (CS), which are normalized in nature as they cover the domain space within range [-1, 1]. Why is there need of natural CS? It is needed to compute the entries of element stiffness matrix conveniently using numerical integration method.

·    Shape functions are not arbitrary selected. They must fulfill the 02 criteria which are: (1) the value of shape function at a given node is 1 and 0 at other node(s), and (2) the sum of shape functions at any point with the element domain must be equal to 1.

Here is demonstration of how to choose shape functions and recheck them.

One the choice of shape functions is made, then, we express the displacement and strain fields within elements in terms of shape functions and natural coordinates, present them in matrix form as following.

Now, to drive the element stiffness matrix, we can deploy various methods i.e. weighted residual, minimization of potential energy etc. The minimization of potential energy principle is very commonly used in the field of engineering structures, and the same we will be using here. Using this principle, we equate the work done on the structure with the corresponding strain energy stored in the structure. After doing tedious mathematics, we can easily derive the element stiffness matrix as following.

For 1D bar element with 2 nodes, we solve this integral to get the element stiffness matrix. Also, we can nodal force and displacement vectors, as following.

Developing Finite Element Mathematical Model (Direct Method)

 In this article, we will be discussing how to develop element mathematical model from scratch using the direct method. The development of an element mathematical model is a crucial task, and it belongs to step # 2 of the procedure, requiring to perform 6 steps in total, for solving problems using FEM. For further details on FEM procedure, see previous articles on this topic.

When we are using a commercial FEM software package, we choose an appropriate element type from the element library of the software depending on the nature of our problem. Being end users, we do not need to derive and implement the element mathematical model. We just need to use these ready-made elements. However, using there elements without knowing how they are developed, is not that interesting. Is not it?

Let assume that we may not have access to costly commercial software packages, or sometimes we want to solve a problem and there is no suitable element available in the element library of an open-source or commercial packages, then, what will we do? The answer is that we simply drive the element mathematical model, implement it and solve our problem.

The development of element mathematical modelling can be broadly done either by direct formulation method or using formal approach. The direct method is based on intuition and is suitable for solving simple problems. On the contrary, the formal approach requires intensive involvement of mathematics and is suitable for complex problems. Moreover, the formal approach is attributed to various methods like the weighted residuals methods (e.g. Galerkin method etc.), and the variational methods involving the calculus of variation and the minimization of potential energy (e.g. the Rayleigh-Ritz method etc.).

Here we will be deriving the mathematical model of 1D bar/spring element having 02 nodes and 1 dof per node, using direct stiffness method. It should be kept in mind that for understanding of FEM, the development of element mathematical model using direct formulation is a good starting point.

Consider an axially loaded elastic bar member as shown in Figure 1a. The governing differential equation for the bar element is as following:

In case of FEM, we need to discretize the continuous do

main of axially load bar structure member. The discretized element having 2 nodes, where each node represent 1 nodal displacement dof, subject to axial nodal forces, is shown in Figure 1b.

Let’s express continuous displacement in terms of nodal displacements of discretized element. Applying force equilibrium at each node, we can drive force equation, as shown below. Finally, expressing the force equations in the matrix form, which is the element stiffness matrix we are looking for. How simple force equation is.

The above element stiffness matrix will be assigned to the discretized domain, followed by the element assembly process and other steps involved in performing FEA to solve the problem.

Source of Errors in a FEM based Solution

When solving an engineering problem using FEM, we begin with developing its mathematical model. The mathematical model cannot cover all details of the physical problems, therefore, we need to simplify it and include the essential aspects of physical problem only. The mathematical model involves the solution of the governing differential equation for the system to get its exact solution which is often not possible.

The next step is to solve this mathematical model using FEM. For this purpose, we discretized the model domain into finite elements. This is equivalent to transforming the governing differential equation into a set of simultaneous algebraic equations which are solved for unknown variables to get the required numerical solution with the help of computer.

In doing above procedure, various kinds of errors are involved as shown in Figure 1.

1. When we represent a physical system with a mathematical model, we make several assumption. That’s why a mathematical model cannot fully describe the all details of a physical system. This is called the modelling error. The better the assumption we make in building a mathematical model of a physical system, the lower will be the modelling error.
2. Solution of the governing differential equation of a mathematical model is exact if we are able to solve it. Often, we resume to FEM to get a numerical solution which is an approximation. Here is the involving of discretization error. It can be reduced by increasing the mesh density at the cost of computational resources. 
3. Lastly, whenever we compute a numerical solution, a third type of error is involved, called the numerical error. It is due to the limitation of computer memory to store numerical values. It depends on the specification of the computer system and can never be zero.

Figure 1. Sources of errors (i.e. modelling error, discretization error & numerical error) involved in getting FEM based solution of a physical problem.

Being engineer, we need to be familiar with these error sources. Let me further explain these errors with a following example.

Let’s we want to analysis the structure of an overhead storage tank as shown in Figure 2(a). The actual physical system is complex in nature and it would be rather difficult to develop a mathematical model that cover all details of this physical system. Remember we are engineers, not mathematicians, so we begin to simplify this problem to develop its mathematical model. Here are few assumptions that can make our lives at peace:

• Assume that the ground behaves as a rigid body.
• The inlet and outlet water piping and its connections does not contribute to the structure performance
• The effect of staircase can be neglected
• The weight water body can be distributed among the columns. 

By applying the above assumptions, we are now able to simplify the problem as shown in Figure 2b. Whenever, we build a mathematical model of the physical problem, the modelling error is involved that depends on nature of assumptions we establish to simplify the problem. The mathematical model involves the solution of governing differential equation which do not have close-form solution in the most of practical applications. But do remember that we get an exact solution if we are able to solve the differential equation.

Figure 2. An overhead water storage tank, (a) physical problem, (b) mathematical model, and (c) discretized model.

Let’s resume to FEM to solve this mathematical model to get an approximate solution. We discretize the problem domain into finite elements as shown in Figure 2(c). During discretization, there is involvement of discretization error as we are getting an approximate numerical solution. If we discretize the model very fine, then, the discretization error will be small but the computational cost will be high. Typically, we begin with a coarse mesh, then, we keep on refining the mesh density while keeping an eye on the solution. At the onset, the solution changes drastically with increase in the mesh density until it begins to converge. At this point, a further increase in the mesh size does not affect the solution, and it is said the mesh convergence is reached. So, this discretization error is the second source of error that we should keep in mind while doing FEA of a physical problem.

The last source of error is the numerical error. This error always exist whenever we make use of computer as it has limited memory to store the numerical values being computed. This error cannot be avoided in the numerical computations due finite precision involving floating point or integer values.

In nutshell, whenever we apply FEM to solve an engineering problem, we must understand that the numerical solution we are getting is infect with various types of errors i.e. modelling error, discretization error and numerical error. Consequently, we should always doubt the solution we are getting from FEM and always try to verify or validate the results whenever possible before making a critical decision.

About the Author: Dr. Khazar Hayat is a professional engineer with almost 15+ year of experience in research, design, analysis and development of products made of fiber reinforced plastics composites (FRPCs). Currently, he is working as an Associate Professor at Mechanical Engineering Department, The University of Lahore, Pakistan, can be reaching by emailing at khazarhayat@gmail.com.

Coordinates Transformation in FEM

In this article, we will be discussing how the coordinate transformation is carried out to transform the element information i.e. element stiffness, element nodal displacements and element nodal forces, in the local/element coordinate system oriented at an angle to the global coordinate system.

This is a necessary step to be perform before starting the element assembly process. We need to check the orientation of local/element coordinate system xy-CS, used for building the element model with reference to the global coordinate system XY-CS. If the xy-CS is aligned with the XY-CS, then, we can begin the assembly of elements by constructing the element connectivity table as discussed in the previous article. Otherwise, it is essential to do coordinate transformation to align the element local information to the global coordinate system before the element assembly if the xy-CS is oriented at a certain angle  to the XY-CS.

We need to transform all information of an element model i.e. nodal displacements, nodal forces and element stiffness matrix. In the below text, we discuss each of them one by one.

Transformation of Element Nodal Displacements

Let consider a linear spring element model with 2 nodes, as shown in Figure 1. Each node has 02 displacement degrees of freedom, which are represented with blue colored vectors along the element/local coordinate system xy-CS.

Please pay attention to the fact that the xy-CS is NOT aligned to the global coordinate system XY-CS, and is oriented at a certain angle. Therefore, we have to transform the blue colored nodal displacements to red colored nodal displacements which are aligned to XY-CS.

Figure 1. Transformation of element nodal displacement

Now, we discuss how to perform this transformation of the element nodal displacements. Here is the geometrical representation of the transformation between the components of nodal displacements in xy-CS and XY-CS in Figure 2.

Figure 2. Geometrical representation of the nodal displacement transformation

In terms of mathematical equations, we can write the same geometrical representation of nodal displacement transformation as following in Figure 3.

Figure 3. Mathematical representation of the nodal displacement transformation

We prefer writing all nodal displacements in the following matrix form as a matter of convenience.

Where [T] is the transformation matrix, {u'} is the nodal displacement vector in the xy-CS oriented an angle, and  {u} is the nodal displacement vector in aligned with the XY-CS.

Transformation of Element Nodal Forces

In a similar fashion for the nodal displacements, we can develop transformation relation for the nodal forces as shown below in Figure 4.

Figure 4. Transformation of element nodal forces

The above information can be also present in the matrix form as below:

Where [T] is the same transformation matrix, {f'} is the nodal displacement vector in the xy-CS oriented an angle, and  {f} is the nodal displacement vector in aligned with the XY-CS.

Transformation Matrix [T]

For conversion of nodal displacements and forces, we developed a transformation matrix [T] with entries of trigonometric sine and cosine ratios, shortly denoted by symbols l and m. This transformation matrix is a very special matrix because if we take its transpose or try to invert it, we get the same result. The matrix demonstrating such special property is called orthogonal matrix. Hence, the transformation matrix in our case is an orthogonal matrix as its transpose is equal to its inverse, as shown below.

The orthogonality of transformation matrix is very useful in numerical computation involving the matrix inversion. Whenever we need to inverse an orthogonal matrix, then, rather than taking its inverse we just take its transpose as inverting a matrix is a computationally expensive task as compare to taking its transpose which is merely the switching of its rows with columns. Please do remember that if a matrix is not an orthogonal matrix then we need takes its inverse as usual.

Transformation of Element Stiffness Matrix

So for we are able to develop the transformation relationships for the nodal displacement and load vectors using the transformation matrix [T], as shown below.

The remaining is the element stiffness matrix that also needs to be transformed from the xy-CS to the XY-CS. However, in this case use the previous transformation equations for the nodal displacements and forces, and manipulating them to derive the transformation relation for the element stiffness matrix as following.

Summary

Here is the summary of complete transformation relations for the element information in xy-CS oriented at an angle to XY-CS. 

Once all element information is fully transformed from -CS to -CS, then, we proceed to the assembly process by constructing the element connectivity table. For further details on the element assembly process, please visit following article.

About the Author: Dr. Khazar Hayat is a professional engineer with almost 15+ year of experience in research, design, analysis and development of products made of fiber reinforced plastics composites (FRPCs). Currently, he is working as an Associate Professor at Mechanical Engineering Department, The University of Lahore, Pakistan, can be reaching by emailing at khazarhayat@gmail.com.

How a FEM software package assemble elements?

One of the strengths that makes the FEM a tool of choice is its generalized way of assembling the elements in the form global stiffness [K], and the global nodal displacement vector {U} and the global nodal load vector {F}, to build a set of simultaneous algebraic equations i.e. [K]{U} = {F}, to be solved for unknown variables.

In this article, we will go through an example problem to explain how a FEM software performs the element assembly process. The problem chosen is a linear 1D bar with legs. One of the legs is subjected to an axial tensile load while others are fixed, as shown in Figure 1.

Figure 1. A linear 1D example problem.

Let’s begin with the decomposition of the problem into small segments. Say the spatial domain of the problem is discretized into 4 subdomains so that we can easily managed their computation as shown in Figure 2. Each subdomain is assigned with an element model of a linear spring type element with two nodes, where each node has one displacement dof, as shown in Figure 3. We will not go into details of how to get this element model. Just assume that we selected it from an existing element library of the FEM software and want to use it to solve our problem.

Also, do remember that the nodes in the element model shown are local in nature. In other words, they are based on local or element coordinate system which was used to build the element model. To distinguish them from the global node numbers which we will see in the coming section, we are using a pair of braces to represent these local node numbers.

Figure 2. Discretization of the example problem.

Figure 3. A linear spring element model with local nodes to be assigned.

Now, we need to specify the global node number and element number for the discretized domain. This numbering can be arbitrary. Here is Figure 4, in which the global node number and element number used current problem are shown. The global node numbers are enclosed in circles, while, element numbers are mentioned in rectangles for distinction. Let me emphasis again, you can have your own number sequence, just do not duplicate any number.

Figure 4. Specification of global node number and element number.

After that, we have to build an element connectivity table, which simply shows the connection between the local node number, enclosed in braces, and the global node number, enclosed in circles, for each element having element number enclosed in rectangle, as shown in Table 1. This table should be read as, say for an element # 1, the local node 1 of element model is linked with the global node 4. Similarly, the second local node 2 of the element #1 is linked to the global node number 2, and so on.

Table 1. Element connectivity table.

Before starting the assembly in the form of global matrix and global vectors, we need to define the size of system. The size of system depends on the number of global nodes and number of dof per node.

Figure 5. Initialization of global stiffness matrix, and global nodal displacement and force vectors.

Finally, here begins the assembly process…

For this element we have the local stiffness matrix [k], the local nodal displacement vector {u} and the local force vector {f}, having the local size of 2x2, 2x1 and 2x1, respectively, as there are two local nodes, with each having 1 dof, as shown in Figure 6.

Figure 6. Element stiffness matrix, and element nodal displacement and load vectors.

We have to place each element of [k], {u} and {f} into their respective positions in [K], {U} and {F} initialized earlier. The positions are determined using the element connectivity Table 1, which we covered earlier. Here is the mapping of local and global positions for each entry for the element # 1.

As a consequence of above connectivity, we can see the populated global stiffness matrix [K], and global nodal displacement vector {U} and global nodal load vector {F}, with the entries of element # 1 in Figure 7 as following.

Figure 7. Global stiffness matrix, and global nodal displacement and force vectors populated with entries of element # 1.

Exactly the same process is repeated for element # 2, 3 and 4, as shown in Figure 8, 9 and 10, respectively.

Figure 8. Assembly process for element # 2.

Figure 9. Assembly process for element # 3.

                                                  Figure 10. Assembly process for element # 4.

Once assembly of all element is done, the final form of global stiffness matrix [K], global nodal displacement vector {U} and global nodal load vector {F}, can be presented in the form of [K]{U}={F}, as shown in Figure 11. This is a set of simultaneous algebraic equations, a so-called presentation of the differential equation of the system, we are looking for.

Figure 11. Final assembled form of all 4 elements.

After the element assembly, we need to apply load and boundary conditions, followed by the solving for unknown global nodal displacement. Once these unknown primary/state variables are determined, we can then compute the dependent variables as per our requirement.

About the Author: Dr. Khazar Hayat is a professional engineer with almost 15+ year of experience in research, design, analysis and development of products made of fiber reinforced plastics composites (FRPCs). Currently, he is working as an Associate Professor at Mechanical Engineering Department, The University of Lahore, Pakistan, can be reaching by emailing at khazarhayat@gmail.com.