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.

Exploring FEM procedure by manually solving an Example

While performing a simulation task using any in-house developed, an open source or a commercial FEM software package, we have to follow these standard steps: (1) discretize the domain, (2) select an appropriate element type depending on the nature of the problem or develop it if not available, and assign this element model to the discretized domain, (3) assemble the elements, (4) apply the relevant load and boundary conditions, (5) solve for unknown primary or state variables, and (6) finally compute the secondary or dependent variables if required. Another terminology often used by the FEM community to describe the above steps in a more brief manner is the pre-process (i.e. step # 1 to 4), the solution (step # 5) and the post-process (step # 6).

In this article, we will be learning how a FEM simulation is conducted behind the curtain of GUI (Graphical User Interface) of a FEM software package. For this purpose, a simple structure problem shown in Figure 1 is selected. We will be solving this problem through FEM and all relevant steps involved will be performed manually to get a deeper insight of this amazing yet a powerful method. This learning journey will be a bit intensive but it will finally pay you off if you show patience.

Problem Statement

So here is the problem, we want to solve with the FEM. It is a simple bar with cross-section A(x) varying along its length L, fixed at its left end and subject to an axial load P at its right end, as shown in Figure 1. The bar is made of an isotropic material having Young’s modulus E. We want to compute the displacement u at the free end.

Figure 1. An axially loaded bar

Here is the so-called “governing differential equation” of problem along with the relevant boundary conditions.

Even though if someone is very good at solving differential equations, still, there is a slim chance to find an exact solution at the first attempt. Majority of us are not proficient at solving differential equations, that’s why, we want to solve this problem with a simple yet generalized way i.e. FEM. You will see in the following why it is so convenient to solve problem using FEM.

Finding Approximate Solution using FEM

Step # 1 Discretization

The first step in performing FE analysis involved the decomposition of the geometry of the problem (i.e. continuous domain) into simpler parts (i.e. subdomains). This process is also called meshing. Generally, all FEM software provide various mesh control options to discretize the geometry automatically or manually depending on our needs. We can also mesh the geometry outside of a FEM software or use another a software dedicated to meshing only, and import the meshed model into the FEM environment.

For the problem on hand, we will mesh the axially loaded bar into two subdomains. We can further increase the discretization, and it will surely increase the accuracy of our approximate solution, however, we limit ourselves to two subdomains only for a very good reason. And that is we will be performing the procedure of FEM manually. So here is how the discretized model looks like (see Figure 2).

Figure 2. Discretization of domain (i.e. Step # 1)

Please have a closer look at Figure 2, each subdomain has a constant cross-section. For these prismatic subdomains we already have simple solutions with which we are very much familiar in mechanics of materials course. In FEM, we are now trying to use these simple solutions for constant cross-sections, to approximate the problem with a varying cross-section along its length.

Please also note that in the above problem, we are interested in the displacement in the x-direction only, so it is 1D problem. Besides, we are using a global coordinate system for the discretized model on hand. It is named so because complete system is built using this coordinate system.

Step # 2 Develop Element Model and Assign to Discretized Domain

In this step, we have to develop an element model based on the nature of problem under investigation. In our case, the axially loaded bar behaves just like a spring. The bar deforms when loaded, and upon removal of load, it returns to its initial undeform shape.

There are numerous ways for driving an element model e.g. direct formulation, minimum total potential energy method, weighted residual method etc. Below Figure 3 shows the derivation of the element model through direct formulation method as the element stiffness matrix is derived using simple formulae of mechanics of materials. So, we model a linear spring element with two nodes. The mathematical model of a m spring element having local nodes i,j is shown below in Figure 3. The derived element model is based on a local/element coordinate system.

Figure 3. A spring element model

The next Figure 4 shows the derivation of the spring element model. Since, our focus here in not on exploring these ways of driving element model, rather than learning FEM general procedure, so we leave this discussion to another day.

Figure 4. Derivation of a spring element model.

Note that being users of a FEM software, generally we do not need to develop these element models. We have to just select an appropriate element type from the element library of the software which contains a plenty of element covering various engineering fields. The software developers have already derived and implemented these element models.

Once we developed our element model or select an element type from the element library of the software, then, we have to assign that element model to the discretized domain. Depending on the requirement of the problem, there can be more than one element type that can be assigned to various portion of the domain, but, one portion of the domain cannot be simultaneously assigned more than one element type.

In our case, we have discretized the domain into two subdomains, and each of subdomain is assigned a linear spring element with two nodes as shown in Figure 5.

Figure 5. Assignment of element type/model to the discretized domain.

Before we move ahead, it is very important to remember that the element model we developed or element type we select from the element library of the software, is built on the local/element coordinate system. For distinction, we write the local node numbers of element in the braces.  On the contrary, the discretized model is built on the basis of global coordinate system, and therefore, we write the global node numbers of the discretized model in circle to make a clear distinction.

We do we need to make above distinction? Because the element model is built in local/element coordinate system. When we assigned the elements to the discretized model built in global coordinate system, we have to transform the relevant attributes of the element from local to global coordinate system.

Step # 3 Assembly of elements

The next step involved in a general procedure of FEM, is to assemble the elements which were previously assigned to the discretized domain. In doing so, we need to map the local node numbers of the elements to the global node numbers in the system. This is done in the form of an element connective table as shown n Figure 6. The first column of the table lists the element number of elements in the discretized domain. The next columns list the global node numbers and corresponding local node numbers for each element.

Table 1. Element connectivity table

Before, we start to assemble the elements, an empty global stiffness matrix, a global displacement vector and a global load vector, each populated with zero values, is constructed.  The size of the global matrix and the global vectors is based on the number of global nodes (i.e. 3 in our case) and number of degrees of freedom per node (i.e. 1 in our case), and is determined as shown in Figure 6.

Figure 6. Sizing of global stiffness matrix and global vectors.

With the help element connectivity table, we transform the local stiffness matrix of each element to corresponding global stiffness matrix. After the completion of this operation, we add the elements of global stiffness matrix of each element to get the global stiffness matrix of the system, as shown in Figure 7. Likewise, we construct the global nodal displacement vector and the global nodal force vector.

Figure 7. Assembly of elements

Here is the final form as shown in Figure 8, the assembled elements in the form of a set of simultaneous algebraic equations, and computers are very very good at solving such system of equations without getting bored.

Figure 8. System global stiffness matrix, nodal displacement vector and nodal force vectors.

 

Step # 4 & 5 Apply load & boundary conditions and Solve for unknown Primary variables

Before we proceed further we need to apply the load and boundary conditions. There is a tip load at global node number 3. The displacement of global nodal number 1 is zero as this node is fixed boundary condition. There is no external force on global node number 2, so it will have net zero nodal force, as shown in the upper portion of Figure 9.

After the application of load and boundary conditions, the set of simultaneous equations is solved for unknown primary variables i.e. displacement dof at each node. The computer solve the set of simultaneous equation of the form Ax = b, by inverting the matrix A and multiplying it to vector b. To inverse a matrix, it should be not a singular matrix because a singular matrix cannot be inverted.

Figure 9. Application of load and boundary conditions, and solve for unknown primary variables.

Step # 6 Compute Dependent or Secondary Variables

The last step is an optional post-processing step. It includes result visualization, data plotting, or computation of the secondary variables from the primary variables etc. The example of computed secondary variables in the current problem could be the strains which can then be converted into stresses using the constitutive relationship at the element level, as shown in Figure 10.

Figure 10. Compute dependent or secondary variables.

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.