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.