What is Finite Element Analysis (FEA)?

The finite element method (FEM) or finite element analysis (FEA) is a numerical technique used to solve problems of engineering and mathematical physics. This technique is used to find the approximate solution of engineering problems, described by partial differential equations. It is very useful for problems with complicated geometries, loadings, and material properties where an analytical solution cannot be obtained. 

These equations are essentially boundary value problems in which one or more dependent variables must satisfy a differential equation over a known domain. The domain is generally a physical area of interest. The boundary conditions are the field variables or derivatives on the boundaries of the domain. Boundary conditions may be of type displacement, temperature, support, velocity, and heat flux.

How is FEA carried out?

The domain is divided and represented into small pieces known as elements, the corner points of these elements are called nodes. Approximating functions in finite elements are determined in terms of the nodal values of the elements which are sought. In this way, a continuous physical problem is transformed into a discretized finite element problem in which nodal values are calculated. Then the global system of equations is solved and the unknown nodal values are calculated. Mathematically, FEM gives an approximate numerical solution to a described physical problem. It is not easy to determine whether the obtained solution is desirable for every problem. Analytical or experimental solutions can help a lot to verify the finite element result.

Procedure to solve an FEA problem

The solution of a general static structural finite element problem in a step-wise procedure is stated below: 

Step 1: Discretization of the domain/structure 

The very first step in FEA is to divide the domain into finite elements or subdomains. Basically, discretization is done by a preprocessor program. Mesh description has mesh attributes like node numbers, element connectivity, and element type. 

 Step 2: Selection of interpolation function or displacement model

As the displacement solution of a complex structure under load cannot be predicted accurately, we assume some suitable solution which is used to interpolate the field variable over an element. These assumed solutions are often in the form of a polynomial function.

Step 3: Finding element stiffness matrix

The stiffness matrix for an element is derived from the equilibrium equation or a suitable variational approach principle. Therefore, equations in a matrix form are established which relate the unknown nodal values to other parameters.

Step 4: Finding a global equation system

The complete structure is composed of finite elements with respective element stiffness and load vectors. All element equations are assembled in a proper manner to get a global equation system. In simple words, all local element equations for all elements are combined which are used for discretization. 

 `nbsp;                                                                                                        `{K} ϕ ⃗=P ⃗`                                                                           (1`

The global equation system is formulated in the form of equation (1). Where, `{K}` is the global stiffness matrix, `ϕ ⃗` is the nodal displacement vector and `P ⃗` is the nodal force vector for the complete structure.

Step 5: Solving the global equation system

Boundary conditions are imposed on the global equation system before solving it. In the finite element method, the global equation system is typically sparse, symmetric, and positive definite. For linear problems, the vector `ϕ ⃗` can be easily solved because the stiffness matrix remains constant. Whereas, for non-linear problems the vector `ϕ ⃗` is solved in a sequence of steps, here the stiffness matrix changes due to the nonlinearity.

Step 6: Calculating the stress and strain

Once the nodal displacements are known, the stress and strain values for elements are calculated using the structural mechanics equations (stress-strain relation, strain-displacement relation).

An example of an axially loaded step bar is explained to show the application of the above steps.

Figure 1.1 shows the axially loaded bar with force`F`acting at node 3 and fixed at node 1 end. The cross-section area, Young’s modulus, and length of each bar element are denoted by `A`, `E` and `L` respectively. A local coordinate system is defined along the axial direction by `x`. This problem is solved step by step for unknown displacements using all 6 steps shown above using the finite element method.

How do we calculate the stress and strain?

Step 1: Discretization of the structure 

The given bar is discretized into two separate elements, which are element 1(left) and element 2(right) with different cross-section areas, material properties, and length. The discretized elements are shown in figure 1.2 above. Assuming this bar to be one-dimensional, it will have only axial displacement as a degree of freedom per node. There are three nodes in the structure so it has three unknowns as ϕ1, ϕ2, and ϕ3 as shown in figure 1.2.

Step 2: Selection of interpolation function or displacement model

Figure 1.3 depicts an element with defined properties, loads, and displacement along with the boundary condition. By approximating the unknown displacement through a linear function we can write,

 

                                                                                                                  `ϕ(x)= a+bx`                                                                                                         (2)

 

Here `a` and `b` are constants. From boundary condition, the end displacements can be written as `ϕ1`(at `x=0`) and `ϕ2` (at `x=l^e` element length). Further using these boundary conditions we get,

 

                                                                                                                    `a= ϕ_1^e`                                                                                                          (3`

                                                                                                         `b= ((ϕ_2^e-ϕ_1^e))/l^e`                                                                                               (4`

 

So, equation 2 can be written as,

 

                                                                                                 `ϕ(x)= ϕ_1^e (1-(x/l^e ))+ϕ_2^e (x/l^e)`                                                                                (5)

 

Step 3: Finding element stiffness matrix

By using the principle of minimum potential energy theorem we can derive the stiffness matrix of an element. Potential energy () for a given bar can be written as,

 

                                                                                            `I="Strain energy-work done by external forces"`          

                                                                                            `I= π- W_p`                                                                                                                                 (6)

 

Strain energy for an element ‘e’ can be written as (this equation is derived from strain energy for the axially loaded elastic bar),

 

                                `π^e= A^e ∫_0^(l^e)1/2 σ^e ε^e dx=A^e ∫_0^(l^e)1/2 (E^e ε^e) ε^e dx= ( A^e E^e)/2 ∫_0^(l^e) (ε^e)^2 dx`                                     (7)

 

Further, the strain of an element can be described by equation 8 and is evaluated from equation 5,

 

                                                                                      `ε^e= (∂∅)/(∂x)= ((ϕ_2^e-ϕ_1^e))/l^e`                                                                                             (8)

 

So, substituting equation 8 in equation 7 we get strain energy for an element as,

 

`π^e= ( A^e E^e)/2 ∫_0^(l^e) {(ϕ_2^e-ϕ_1^e)/l^e }^2 dx= ( A^e E^e)/(2l^e ) {(ϕ_1^e)^2+(ϕ_2^e)^2- 2 ϕ_1^e ϕ_2^e }`                                                      (9)

 

Equation 9 can be written in the matrix form as shown below,

 

                                                                                      `π^e= 1/2 (ϕ ⃗^e)^T [K^e ] (ϕ ⃗^e)`                                                                                                (10)

 

Where,

`ϕ ⃗^e= ((ϕ_1^e),(ϕ_2^e) )`  is the nodal displacement vector for an element ‘e’

For element 1, `ϕ ⃗^e= ((ϕ_1),(ϕ_2 ))` and for element 2, `ϕ ⃗^e= ((ϕ_2),(ϕ_3))`

`K^e=( A^e E^e)/l^e [[1 -1],[-11]]= [[K_11^eK_12^e],[K_21^eK_22^e ]]``is the stiffness matrix of an element ‘e’

 

Now nodal loads are acting on each node of the bar, the work done by these nodal loads can be written as,

 

                                                                                         `W_p= ϕ_1 P_1+ ϕ_2 P_2+ ϕ_3 P_3`                                                                                              (11)

 

The `P_i` denotes the force applied at each node of an element and `ϕ_i` denotes the nodal displacement. In the above-defined problem, `P_1="reaction force" (R)`, `P_2=0` and `P_3=F`

 

Finally, the complete load vector for the given bar can be written as, `P ̅= ((P_1),(P_2),(P_3 ))`

 

Step 4: Finding a global equation system

In this step, all the element equations will be assembled together in a suitable manner to get the global equation system. The element stiffness matrix, load vector, and displacement vector all these equations will be assembled to get the overall global equilibrium equations.

Through the principle of minimum potential energy theorem,

 

                                                                                                  `(∂I)/(∂ϕ)_i =0`                                                                                                                 (12)

 

So, using equation 10 and 11, equation 6 can be written as,

 

                                                                                    `(∂I)/(∂ϕ_i) = ∂/(∂ϕ_i) (∑_(e=1)^2(π^e- W_p ))=0`                                     

                                                                                       `∑_(e=1)^2([(K^e ] ϕ ⃗^e- P ⃗^e ) =0`                                                                                         (13)

 

The summation in equation 13 denotes vector assembly in which the elements corresponding to a particular degree of freedom in different vectors are added.

The stiffness matrix for element 1 can be written as,

 

                                                               `K^1=( A^1 E^1)/l^1 [[1-1],[-11]]= [[K_11^1K_12^1],[K_21^1K_22^1 ]]`                                                             (14)

 

Similarly the stiffness matrix for element 2 can be written as,

 

                                                                 `K^2=( A^2 E^2)/l^2 [[1-1],[-11]]= [[K_11^2K_12^2],[K_21^2K_22^2]]`                                                            (15)

 

For element 1, the displacements at node 1(left) and node 2(right) are denoted by `ϕ_1` and `ϕ_2`.

For element 2, the displacements at node 1(left) and node 2(right) are denoted by `ϕ_2` and `ϕ_3`.

 

So, the global stiffness matrix can be obtained by assembling the element 1 and element 2 stiffness matrices. As, there are three unknown displacements (DOFs) `(ϕ_1,ϕ_2,ϕ_3)`, the global stiffness matrix will have a 3*3 matrix dimension. The global stiffness matrix is shown below,

                                                                                                                                    (16)

After assembling each load vector of element 1 and element 2 the global load vector can be written as,

                                                                                             `P ̅= ((P_1),(P_2),(P_3 ))=((R),(0),(F))`                                                                                     (17) 

Here, `R` is the reaction force at node 1 of element 1 and `F` is the applied force at node 2 of element 2. From equation 16 and 17, the global system of the equation can be written in the matrix form as,

                                                                                                                                    (18)

The above global equation system can be derived by potential energy theorem for a bar using equations 9, 11, and 13. A step by step procedure is shown below. The potential energy for the given bar can be written as,

 

                                                                                                          `I= π^1+ π^2- W_p`                                             

`I= ( A^1 E^1)/(2l^1 ) {ϕ_1^2+ϕ_2^2- 2 ϕ_1 ϕ_2 }+ ( A^2 E^2)/(2l^2 ) {ϕ_2^2+ϕ_3^2- 2 ϕ_2 ϕ_3 }- (ϕ_1 P_1+ ϕ_2 P_2+ ϕ_3 P_3)`                                       (19`

Through the principle of minimum potential energy theorem, we know,

                                                                                                               `(∂I)/(∂ϕ_i) =0` `"where" i=1,2,3`

So, by substituting the load vector values, equation 19 yields,

`(∂I)/(∂ϕ_1) =0=( A^1 E^1)/(2l^1 ) {ϕ_1- ϕ_2 }- R`                                                                                                                                                               (20)

`(∂I)/(∂ϕ_2) =0=( A^1 E^1)/(2l^1 ) {ϕ_2- ϕ_1 }+ ( A^2 E^2)/(2l^2 ) {ϕ_2- ϕ_3 }`                                                                                                                  (22)

`(∂I)/(∂ϕ_3) =0=( A^2 E^2)/(2l^2 ) {ϕ_3- ϕ_2 }- F`                                                                                                                                                                (23)

Equations 20, 21, and 22 can be written in matrix form which gives the same result as equation 18.

Step 5: Solving the global equation system

From the given boundary conditions, we know `ϕ_1=0`. By substituting this boundary condition in the global equation system, we can further delete the first row and first column of equation 18.

                                                                                                                                     (23)

Using equation 23, the unknown displacements `ϕ_2` and `ϕ_3` can be evaluated. We have two unknowns and two equations.

Step 6: Calculating the stress and strain

After getting unknown displacements, strain and stress for each element can be further calculated from the derived equations as shown below,

                                                                                               `ε^1= (∂∅)/(∂x)= ((ϕ_2-ϕ_1))/l^1` 

                                                                                               `ε^2= (∂∅)/(∂x)= ((ϕ_3-ϕ_2))/l^2`

 And stress for each element can be calculated as,

                                                                                                             `σ^1=ε^1 E^1`

                                                                                                              `σ^2=ε^2 E^2`