The behavior of any physical system can be defined mathematically. Solving these mathematical equations allows for the analysis of the system’s behavior. Depending on the complexity of the equation, finding solutions directly may not be straightforward. Finite element analysis is an approximation technique to transform the mathematical representation of a physical system into a set of computer solvable algebraic equations and calculate the numerical solution at a set of discrete locations.
There are several types of mechanical problems, such as static, dynamic, vibration, shock, thermal and corrosion. Static problems deal with structures attaining equilibrium under the applied loads. These problems involve analyzing forces and moments acting on a structure to ensure that it remains in state of balance without any motion. Principles such as equilibrium equations, Newton’s laws of motion and summation of forces and moments are used to solve these static problems. In this article let’s understand the steps needed to solve static problems using FEA by hand.
FEA Procedure Overview
The FEA procedure can be categorized into 3 parts as shown in the flow chart below.
In this article, let’s discuss the general steps for defining a finite element problem and calculating a numerical solution using simple 1D elastic bar problem.
Problem Formulation and discretization
Defining The Physical Problem
It is very important to determine the physicality of the problem at hand first. For this, we need to answer certain questions, such as whether the problem is 1D, 2D or 3D. Is it a solid mechanics, continuum mechanics, heat transfer, computational fluid dynamics problem? Is it a static or dynamic problem? Is it a single or multi- body problem? Does it require plastic and damage properties for material?
In this article, we are considering a simple solid bar with fixed cross-sectional area (A). The bar is elastic with Youngs modulus (E). The left-hand side of the bar is fixed in all degrees of freedom, and it is subjected to compressive pressure (t) on the right-hand side. The stresses, displacement and strains in the bar are unknown variables. This can be reduced to a 1D problem as the deformation and stresses are almost uniform across the cross section of the rod and they only vary along its axis. This is a quasi-static problem as we are calculating its deformation when the bar attains equilibrium under the influence of applied loads.
Mathematical Model
Mathematical equations relating input and output variables need to be established in this step. This can be obtained using fundamental laws such as equilibrium equations, elasticity law, heat transfer law and diffusion law.
For the current problem at hand, we need to obtain a relationship between the applied force and elastic deformation of the rod. For this, we use the equation of equilibrium and Hooke’s material law to arrive at the mathematical model. In this step, we also need to identify the boundary conditions and include them in the mathematical model. For the 1D problem, the boundary conditions are the body force b and the pressure t applied at x = 6 in. The mathematical model for the 1D bar problem can be obtained as,
Deriving Weak Form
In this step, optimization requirements can be applied to the strong form equation established in the previous step to obtain weak form. There are different variational principles as shown in the figure below to derive the weak form.
For the 1D bar problem, the weak form is derived to be,
“Weak” refers to the fact that the continuity requirement of the displacement variable u is weakened in this newly derived weak form equation. For the strong form differential equation to hold everywhere on the bar, the first derivation of u (du/dx) has to be continuous, and the second derivative of u(d^2u/dx^2) needs to exist. Whereas for the weak form integral equation to hold everywhere on the bar, u and its variation (δu) should be continuous, and the first derivatives can have discontinuity at finite number of locations and thus the continuity requirement is weakened in this step.
Discretization
In this step, we break the computational domain of the problem to a set of smaller segments. Each segment is termed as an element. This process of discretizing the computational domain is called meshing. The meshing process is demonstrated here by diving the 1D bar into 6 elements (1 in length each). The adjacent elements are connected through a common node. The elements and nodes are numbered for tracking purposes. The meshing process for 1D problems is pretty straight forward. This becomes more involved for 2D and 3D domains with complex geometries. This can be achieved numerically by using algorithms.
Calculating The Solution
Defining Unknown Variables For The Elements
In this step, the unknown variables such as displacement are approximated over the elements in the computational domain. Let’s assume the unknown displacement in the 1D bar problem at nodes 1, 2…6 is u1, u2 …u6 respectively. These displacements are termed as nodal displacements. The displacements over each of the elements can be assumed either to be linear, quadratic, or other higher polynomials as shown in the figure below.
A general polynomial function in x for u can be written as,
Let’s assume linear approximation for the displacement over each element in our model. Calculating the linear polynomial coefficients gives linear shape functions for each element as shown below.
The shape functions and their first derivative values must be calculated for all the elements in this step. For more information on shape functions and their properties please refer to our previous blog here.
Weak Form Discretization
In this step, the global integral in the weak form can be broken down into element integrals:
The unknown displacement function and its derivates must be replaced by finite element polynomial approximations. For the first element, these polynomial approximations are derived to be,
Substituting the finite element approximations into the weak form for all elements gives the elemental stiffness matrix and force vectors.
Global Equation Systems
The elemental stiffness matrices [k] and force vectors {f} should be assembled in the global matrix and vector as shown in the figure below. Each element and its nodal DOF are given a unique global index, which forms the basis for building global matrices. Note that, based on mesh connectivity, some of the elemental stiffness terms will get added in the global matrix.
Loads and Boundary Conditions
There are two types of boundary conditions (BCs): natural and essential. Essential boundary conditions are specified directly to the solution matrices (u, f) while the natural boundary conditions describe phenomenon that affect solution matrices. Examples of essential BCs are fixed displacement or nodal forces and examples of natural BCs are body force, fluxes, and rate changes. The natural BCs are already included in the global vectors during formulation of element matrices, but the essential BC should be applied in this step. There are different approaches like direct substitution method, penalty method and Lagrange multipliers to apply essential BCs. In the 1D bar problem, the essential boundary condition is u|x=0 =0. This has to be applied in this step.
Solving the global equation system
The overall global system of equations (KU =F) with all the boundary conditions applied, has to be solved in this step for solution vector {u}. With increase in size and complexity of the problem domain, solving this equation system becomes tricky and computationally expensive. There are a number of direct and iterative methods available to solve this equation system as shown in the table below.
| Iterative methods | Direct methods |
| SOR (successive over relaxation)Conjugate gradient methodMinimal residual methodGauss-Siedel | Gaussian eliminationGauss-Jordan eliminationLU decomposition |
Post-Processing
The nodal DOF are only obtained by solving the global system of equations. But we often need to calculate other physical quantities such as strain and stress in the element and visualize them in the contour plots or graphs. These calculations are performed in this step.
In the 1D bar problem, we can calculate strain and stress in each element by,
These quantities can be plotted on the graph, or they can be visualized using contour plots or animations.
Final Thoughts
Hopefully this article has given an overview on steps involved in formulating and solving a simple 1D FEA problem. The complexity of meshing process, global matrices and vectors increases with increase in complexity and dimensionality of the computational domains. Thankfully we have commercial FEA software like Abaqus to ease this problem formulation process. Using this software, the repetitive tasks like element formulations, assembling global matrices and solving equations are automated, and solutions can be obtained much faster. Hence the commercial FEA softwares play a vital role in modern engineering design and analysis, facilitating the development of innovative and reliable products.
If you are interested in solving mechanical problems using FEA simulations, feel free to reach us. Our Engineers have many years of experience in commercial FEA, and we would love to talk to you!
