Contact Hours/Week : 3+0+2 (L+T+P)	Credits:	4.0
Total Lecture Hours : 39	CIE Marks:	50
Total Practical Hours : 26	SEE Marks:	50
Course code : N2CSE01	SEE Duration (Hours):	3.0

Course objectives:

The objective of this course is to make the students conversant with basic procedure of finite element analysis of structural problems by introducing the concept of structural discretization and different types of displacement models, different types of elements, and iso-parametric elements. The application of the technique bar, beam, truss, and frame element will be introduced and utilized for solution of 2Dimensional problems of structural analysis.

Unit 1. Introduction to FEM	8+6 hours
Approximate method of structural analysis, Concept of Rayleigh-Ritz method, Advantages and disadvantages of FEM, Basic procedure of FEM for structural problems- Process of Discretization, finite elements for 1-D, 2-D and 3-D problems, Element aspect ratio, mesh refinement Vs higher order elements, numbering of nodes to minimize band width.
Unit 2. Discretization of Structures and Displacement Models	7+5 hours
Nodal displacement parameters, Polynomial displacement functions for standard bar and beam elements triangular elements, and rectangular elements, conditions to be satisfied by displacement functions- invariance, continuity, degree of continuity of displacement functions – C0, C1 and C2 functions, convergence and compatibility, Generalized and natural coordinates.
Unit 3. Interpolation functions and stiffness matrices		8+5 hours
Lagrangian interpolation functions, Shape functions for one, two and 3-dimensional elements, Derivation of element stiffness matrices for Bar, and Beam elements. Linear static analysis of one-dimensional problem using Linear and Quadratic bar elements and beam elements.
Unit 4. Iso-parametric elements and numerical integration	8+5 hours
Concept of Iso-parametric elements, sub and super parametric elements, Convergence requirements for Iso-parametric elements, Condensation of internal nodes, Formulation of Jacobian matrix and strain displacement matrix, and element stiffness matrix, consistent load vector and numerical integration.
Unit 5. Stiffness matrices for 2-D elements and FE softwares	8+5 hours
Two dimensional problems-Derivation of element stiffness matrices and equivalent nodal force vectors for constant strain triangular elements and quadratic elements. Auto mesh generation, Computer Program for FEM – Organization – basic flowcharts, Desired features of Pre and Post Processors.

Textbooks

1.	Tirupathi R. Chandrupatla, Ashok D. Belegundu, Introduction to Finite Elements in Engineering, ‎ Cambridge University Press, 5th edition, 2021
2.	Krishnamoorthy C.S., Finite Element Analysis – Theory and Programming, 2nd Edition, Tata McGraw Hill, New Delhi, 2011.

References

1.	Cook, R.D., Malkus, D.S., and Plesta, M.E. and Robert J Witt, Concepts and Applications of Finite Element Analysis,4th Edition, Reprint by Wiley India Pvt. Ltd, New Delhi of John Wiley and Sons, Singapore Edition, 2007
2.	Desai, C.S. and T Kundu, Introductory Finite Element Method, CRC Press, London, 2001
3.	Rajasekaran S., Finite Element Analysis in Engineering Design, S. Chand and Co., New Delhi, 2006.
4.	Singirsu S. Rao, The finite element method in Engineering, Fourth edition, Elsevier Inc., New Delhi, 2005.
5.	Yang T.Y., Finite Element Structural Analysis, Prentice Hall, New Jersey, 1986.
6.	Zienkiewicz, O.C. R L Taylor and J Z Zhu, The Finite Element Method- Its Basis and Fundamentals, 7th Edition, Butterworth and Heinemann, Elsevier , London 2013
7.	K J Bathe, Finite Element Procedures, 2nd edition, Prentice Hall, Pearson Education Inc., 2014

Course outcomes: After completing this course, the students will be able to,

CO1:	Integrate his basic knowledge of structural mechanics and computational structural mechanics with advanced knowledge of FE concepts, energy methods to find solution to the problems of beams with various loading and boundary conditions.
CO2:	Apply critical thinking and original judgment to formulate displacement functions and shape functions for simple bar, beam and triangular and rectangular element using polynomial and interpolation functions in Cartesian and natural coordinates. He will be able to make creative advanced and apply lateral thinking to check for invariance, degree of continuity, convergence and compatibility requirements of shape functions, and provide solution to simple problems of trusses and frames
CO3:	Distinguish between iso-parametric, sub- and super-parametric elements with examples, check convergence criteria for displacement functions of iso-parametric elements, explain the concept of condensation of internal nodes, and formulate of Jacobian matrix, strain displacement matrix, and element stiffness matrix, consistent load vector for quadrilateral elements and apply numerical integration
CO4:	Think laterally, conceptualize and formulate and solve simple two-dimensional problems using triangular/four noded quadrilateral elements.
CO5:	Apply critical thinking and independent judgment to explain the behaviour of different types of plate and shell elements and explain their conformity, convergence and compatibility. He will be able to formulate the procedure for deriving element stiffness matrix and force vectors for four noded plate elements and specify the procedure for assemblage into global stiffness matrix, applying boundary conditions and finding solution for displacements, rotations and moments.

Lab component: Use of FE software for solution of simple problems of beams, frames, trusses, plane stress and plane strain problems, foundation modeling, Axisymmetric problems

1 INtroduction to FEM

In this chapter the following topics are covered

· Rayleigh Ritz method as an approximate method similar to FEM

· Galerkin method as an approximate method similar to FEM

· Basic procedure of FEM for structural problems

· Advantages and Disadvantages of FEM

· Types of Finite Elements

· Assembling of elements and effect of node numbering on band width

1.1 Rayleigh Ritz method as an approximate method similar to FEM

1.1.1 Two concepts:

This method uses two concepts:

1. Any unknown function can be expressed by a sum of a series of base functions and their coefficients;

2. There exists a functional (like potential energy) which needs to take a minimum value of all possible values. In structural mechanics problem, we use the functional called the Potential energy function() , which is the sum of the change in internal strain energy (U) of the body and the work done (W) by the external forces.

negative sign indicates the work is done on the body. Internal strain energy is work done by the body (due to internal forces).

1.1.2 Steps

1. Express the required displacement solution in discretized form using base functions

2. Express the internal strain energy (U) in terms of discretized solution

3. Express the Work done (W) by external forces in terms of discretized solution

4. Solve the minimization of the potential energy (), to get the unknows of discretization

1.1.3 Illustrative Example #1: Cantilever beam subjected to UDL

1. Approx. discretized solution is

The functional to be minimized is the potential energy (), which is the difference between strain energy (U) gained by the body and the work lost by the loads

Substitute the discretized solution into the functional

should be such that is minimum

5. So, the approximate solution is

So, the approximate value of free-end deflection is

Exact Value is

1.1.4 Example #4: Simply supported beam subjected to point load P using polynomial function

1. Approx. Discretized solution with trigonometric base function

The functional to be minimized is the potential energy (), which is the difference between strain energy (U) gained by the body and the work lost by the loads

Substitute the discretized solution into the functional

should be such that is minimum

Exact deflection at the centre of SS beam with conc. load

1.1.5 Example #4: Simply supported beam subjected to point load P using trigonometric functions

A diagram of a line with a red arrow pointing to the line

AI-generated content may be incorrect.

2. Approx. Discretized solution with trigonometric base function

Both these conditions are satisfied for any value of

The functional to be minimized is the potential energy (), which is the difference between strain energy (U) gained by the body and the work lost by the loads

Substitute the discretized solution into the functional

should be such that is minimum

Exact deflection at the centre of SS beam with conc. load

1.2 Galerkin method as an approximate method similar to FEM

1.2.1 Example#1:

Galerkin method- cantilever beam with concentrated load at free end-

Given beam
Approximate function
Check for admissibility
Formulation of residual
Integration of weighted residual

Exact solution	At This gives the exact deflection

1.2.2 Example#2:

Galerkin method- cantilever beam with concentrated load at free end- Triginometric

Given beam
Approximate function
Check for admissibility
Formulation of residual
Integration of weighted residual

Let
Exact solution	At At

1.2.3 Example# 3:

A bar fixed at both ends and subjected to a concentrated axial load at the mid-section.

Simply supported beam with concentrated load at mid-span

1.3 Basic procedure of FEM for structural problems

The basic/ general procedure of a Finite Element Analysis (FEA) involves

Select suitable field variables and the elements

· Discritize the continua

· Select interpolation functions

· Find the element properties

· Assemble the element properties to get the global properties

· Impose the boundary conditions

· Solve the system of equations to fet the nodal unknowns

· Additional quantities are derived from the nodal unknowns

The explanation for each these steps with an illustrative example is given below

1.3.1 Selection of the suitable variables and the element

The given problem can be modelled by using bar elements or plane stress elements. If only axial variations of displcements and stresses are sufficient, bar element modelling is enough. If the displacements and stress distribution across the cross-section are required the plane stress modelling is to be done. In this case, to reduce the computational effort, we can use the two-way symmetry. Thus, only one-quarter of the specimen can be modelled by Finite Elements.

1.3.2

Discritize the continua

As shown in the previous figure two elements are used in bar-element modelling and 18×4 = 72 elements are used in the plane stress modelling with 4-noded quadrilateral elements. The total number of nodes in the plane-stress modelling is 19×5=95 nodes. At every node of 4-noded quadrilateral elements, there will be 2 degees of freedom (), the total unknown displacements in the plane-stress case will be 95×2=190.

1.3.3 Select the interpolation functions

The interpolation function will be same as the shape functions used in the formulation of the selected elements. If the 4-noded quadrilateral elements are used, then the shape functions used to interpolate the field variables u, v will be as follows

1.3.4 Element properties

For each element the strain-displacement matrix [], and using it, the element stiffness matrices [] are determined by integrating the over each element. Numerical integration is usually used to do this integration. The [D] matrix corresponding to the plane-stress condition is used. [] matrix will be 8 x 8 matrix for 4-noded quadrilateral element.

For each element the nodal force vector shall also be calculated using consistent formulation.

1.3.5 Assembling

The element stiffness matrix [] of each element is assembled into the global stiffness matrix based on the DOF of the element. If an element has DOF numbers as , then the elements of [] will be placed in these rows and columns as shown below (only first four rows and first four columns are shown). Here are the elements of the matrix of element. When the subsequent elements are assembled the values at the overlapping locations get added up algebraically.

The size of the global stiffness matrix will be 192 × 192, where 192 is the total number of DOF. The nodal force vectors of each element are similarly assembled into the global nodal force vector, which is of size 192 ×1.

1.3.6 Applying the boundary conditions

The boundary conditions are applied by identifying the degrees of freedom whose values are known. For example, in the above problem, the degrees of freedom corresponding to the nodes on the fixed edge will be all zeros. So, those rows and columns can be eliminated from the matrix. This is the elimination method of imposing BCs. Due to symmetry the vertical displacements at the left face of the model will also be zeros. The corresponding rows of the force vector shall also be eliminated to get the reduced global force vector.

1.3.7 Solving the simultaneous equations

The reduced global stiffness matrix will be solved by multiplying its inverse with the reduced force vector, to get the global unknown displacements.

1.3.8 Post-processing

Depending on the requirements, the deformed shape, stress at interior points can be calculated and plotted.

1.4 Advantages and disadvantages of FEM

1.4.1 Advantages: FEM is applicable for numerous field problems due to its ability

· To model complex shaped bodies quite easily –geometries

· To handle several load conditions without difficulty- loading

· To handle different kinds of boundary conditions- boundary conditions

· To handle different types of materials in the domain- materials

· To handle different types of assumptions in different portions of domain- elements

· To handle different levels of refinement required at different places of the model- Mesh

· To handle nonlinear behaviour in different portions of the model- Nonlinearity

· To handle steady-unsteady, compressible & incompressible, laminar and turbulent flow problems

· To handle time dependent problems- Dynamics

1.4.2 Disadvantages:

It requires great level of expertise & discretion

1. to understand the underlying assumptions

2. to choose the appropriate elements

3. to choose the appropriate level of discretization

4. to choose the appropriate type of analysis

5. to choose other parameters of the modelling

2 Discretization of Structures and Displacement Models

2.1 Discretization of geometry

2.1.1 Aspects of Finite Element meshing

2.1.2 Types of Finite Elements

Objective: To reduce the number of unknowns in the solution to finite number of unknowns

Method: The domain is divided into finite number of elements and the variation of the field variable is assumed within each element, so that if the nodal values of the filed variables are known at the nodes, then the whole solution can be determined.

Various types of Elements:

The element is defined by the number of nodes and the type of interpolation. Based on dimensions, the elements can be classified as follows.

One-D elements:

· Bar elements : 2-noded, 3-noded, 4-noded, etc

· Beam elements: cublic elements, shear-flexible beam elements

Two-D elements:

· 3-noded triangular elements

· 6-noded triangular elements

· 4-noded quadrilateral elements

· 8-noded quadrilateral elements

Three-D elements:

· 4-noded tetrahedron elements

· 8-noded hexahedral elements

· 20-noded hexahedral elements

2.1.3 Assembling of elements

2.1.4 Numbering nodes for band-width minimization

The band width of the overall or global characterisitcs matrix depends on the node nubering scheme and the number of degrees of freedom considered per node. If we can minimize that band-wdith, the storage requirements as well as the solution time can also be minimized. Since, the number of DOF per node is generally fixed for any given type of problem, the band-width can be minimized by using a proper node numbering scheme.

Band width, [Max difference between the number DOF at the ends of any member+1]

Band width (B) = (D+1)*f

Where, D is the maximum largest difference in the node numbers, f is the number of degrees of feedom at each node.

A shortest bandwidth can be obtained simply by numbering the nodes across the shortest dimension of the body.

Case1: Numbering along the shorter direction. There are four spans (bays) and 20 storeys. So, the shorter direction is the width (span) direction, and the longer direction is the height (storey) direction.

Case2: An alternative node numbering will be considered as follows

Therefore it can be seen that when the numbering is given along the shorter direction (Case1), the band width is 15, while the band width value becomes 66, if the numbering is done along the longer direction (Case2)

Banded Matrix:

In a banded matrix, the non-zero elements are obtained or contained within the band and the outside the band elements are zero. In the matrix, all non-zero elements (coefficients) are on either side of the main diagonal joining the elements, , and these are within the band whose width is given by , called the half-band width. The half-band width is defined as the greatest number of coefficients of the matrix in any row of the matrix from and including the leading diagonal to the right hand side non-zero coefficients. The main use of banded matrix lies in reducing the space required to store the given matrix, and to minimize the solution time, the banded matrix is used to store only the non-zero elements. These can be stored completely in - matrix instead of original matrix.

2.1.5 Element aspect ratio

The shape of the element also effects the accuracy of the analysis. Defining the aspect ratio as the ratio of the largest to smallest size in an element, the conclusion of many researchers is aspect ratio should be close to unity as possible. For a two-dimesnional rectangular element the aspect ratio is conveniently defined as the length to breadth ratio.

Effect of aspect ratio on the error in Finite element solution

Fig. 2.1Effect of aspect ratio on the numerical error in Finite element solution

The plot in Fig. 2.1 shows the effect of element aspect ratio on the finite element solution. The higher the value of aspect ratio, larger will be the error in Finite element solution.

2.1.6 Higher Order elements Vs refined Mesh

Accuracy of calculation increases if higher order elements are used. Accuracy can also be increased by using more number of elements.

Aspect	Higher-Order Elements (e.g., quadratic, cubic)	Refined Mesh (More linear elements)
Accuracy per DOF	Higher accuracy with fewer elements and DOFs	May require many elements to achieve same accuracy
Computational Cost	Fewer elements but more complex integration	More elements, simpler computations
Convergence	Faster convergence for smooth problems	Slower, especially for curved geometries or gradients
Geometry Representation	Better for curved boundaries	May need excessive refinement to match curves
Stress Concentrations	May still need refinement locally	Mesh refinement more direct for stress risers
Programming Complexity	More complex shape functions and Jacobians	Easier to implement in custom codes
Solver Performance	Fewer elements but larger bandwidth matrices	Many small elements, can strain memory
Nonlinear Problems	Can be more sensitive and less robust	Linear elements can be more stable under nonlinearity

Use Higher-Order Elements when:

The solution is smooth (e.g. in structural problems without sharp gradients).
You're modeling curved geometries (e.g., beams, shells, contact surfaces).
You want accuracy with fewer elements.
You’re using commercial FEM software (e.g., Abaqus, ANSYS) that handles integration well.

Use Mesh Refinement (h-refinement) when:

The solution has singularities or high gradients (e.g., near holes, cracks, supports).
You are performing adaptive meshing.
You’re solving nonlinear, contact, or plasticity problems.
You are limited to linear element solvers or simpler code frameworks.

Real-World Example:

For a beam under bending, quadratic elements can capture curvature better than many linear elements.
But for a plate with a hole under tension, localized mesh refinement near the hole edge is often more effective.

Best Practice: Combine Both (p- and h-refinement) Often, the optimal strategy is a hybrid:

Use higher-order elements globally, and
Apply local mesh refinement in critical regions.

2.1.7 Polynomial Displacement models

The basic idea of FEM is piecewise approximation. That is the solution of a complicated problem is obtained by dividing the region of interest into small regions (finite elements) and approximating the variation of the solution over each region by a simple function. Thus, a necessary and important step is that of choosing a simple function for the solution in each element. The function used to represent the behaviour of the solution within the element are called “interpolation functions/ models or displacement models or approximating functions”

These are simple functions which are assumed to approximate the displacements for each element. Polynomial displacement function have been mostly used due to the folowing reasons:

1. It is easier to perform mathematical calculations such as integration/ differentitation

2. It is easier to formulate and computerize

3. It is possible to icrease or improve the accuracy of the results by increasing the order of the polynomial. Although trigonometric functions also possess these properties, they are seldom used in FE analysis

The general form of the polynomial displacement function for a 1D case is given by

For 2D case

For 3D case

Where are the coefficients of the polynomial, and are known as “generalized coordinates”.

In most of the practical applications the order of the polynomial in the interpolation functions is taken as one, two and three. Thus,

For first order model, or linear model,

For second order model, or quadratic model,

The greater the number of terms included in the approximation function, the greater will be the accuracy of the solution (closeness to the exact solution).

For 2D case, the number of terms (m) in a complete n^th order polynomial is given by the expression

For 3D case, the number of terms (m) in a complete n^th order polynomial is given by the expression

The terms can be easily remembered by using the Pascals’ triangle in 2D and 3D cases as shown

2.1.8 Convergence requirements of displacement functions

· Compatibility Requirments

· Rigid body modes

· Constant strain mode

· Geometric isotropy or spatial isotropy or geometric invariance

2.2 Generalized Vs Natural Coordinates

Generalized coordinates are the mathematical parameters to describe a variation. For variation of displacements, the expression can be as follows

It can be observed that the coordinates are mathematical parameters that govern the rate of variation, but do not change the type of variation.

2.3 Polynomial displacement functions for standard bar and beam elements triangular elements

2.3.1 Shape functions of bar element from first degree polynomial in local coordinates

The bar element is subjected to axial forces only. The displacement will also have only one component, i.e., the axial component. The displacement variation can be expressed as linear functions of spatial coordinates. The polynomial displacement function will be,

This displacement function should satisfy the boundary conditions at the end of the bar, i.e.,

2.3.2 Shape functions of bar element from first degree polynomial in natural coordinates

This displacement function should satisfy the boundary conditions at the end of the bar, i.e.,

Writing these boundary conditions in matrix form

On substituting in the matrix form of the displacement model

Shape function matrix relates the field variables at a point to the corresponding nodal values as,

On comparing the last two expressions, we can write

2.3.3 Shape functions of 3-noded bar element from second degree polynomial in natural coordinates

The bar element can be subjected to varying tractions along the axial directions. The displacement will have only one component, i.e., the axial component. But, the displacement can vary quadratically along the length of the member. The displacement variation can be expressed as a quadratic function of spatial coordinates. The polynomial displacement function will be,

In matrix form

This displacement function should satisfy the boundary conditions at the end of the bar, i.e.,

Writing these boundary conditions in matrix form

Substituting this, in matrix form of displacement function

2.3.4 Shape function for 2-noded beam element from cubic polynomial

The displacement function needs to be of C1 continuity. So, a cubic polynomial is to be considered for shape function. The field variable for beam element is the transverse displacement or deflection, and the nodal degrees of the freedom consists of deflection and rotation at each node. The

In matrix form

The BCs in terms of the nodal degrees of the freedom are expressed as below

BCs in matrix form

Substituting in the matrix form of displacement function

2.3.5 Shape functions for 3-noded triangular element from linear polynomial displacement

Let the displacement components at any point (x,y) are given by

IN matrix form

On applying the boundary values

All these can be combined into one matrix equation

On substituting the values of

where

The value of N₁ at any point P can be obtained as the ratio of two triangles (PBC and ABC), as shown below

2.3.6 Shape functions for 3-noded triangular element in natural coordinate system

Let the natural coordinates of nodes 1,2,3 be and shape functions be . Since the CST element is linear displacement model, the displacement function shall be

The BCs at the nodes are at node1, are

The BCs at the nodes are at node2, are

The BCs at the nodes are at node3, are

On substituting in displacement function

The displacement model for will be exactly similar. Will lead to the same shape functions

2.3.7 Shape functions for 6-noded triangular element from quadratic polynomial displacement

Consider a 6-noded triangular element as shown in

Fig. 2.2Six noded triangular element in a) global and b) natural coordinates

The displacement function for 6-noded triangular element will have 6 unknowns as given below

In matrix form

Similarly for v also it can be written.

(Alternatively, the shape functions can be derived from the displacement function in global coordinates also

)

BCs (in natural coordinates) are

Node#
1	1	0	0
2	0	1	0
3	0	0	1
4	½	½	0
5	0	½	½
6	½	0	½

On substituting the above BCs and coordinates in the matrix form

If [B] is the inverse of the above matrix

Therefore

Therefore,

Substituting in the first equation

where

similarly

2.3.8 Shape functions for 4-noded quadrilateral element from polynomial displacement function in natural coordinates

The shape functions can be derived from the displacement function in natural (and isoparametric) coordinates as follows (refer to the Pascal triangle of functions)

Substitute the boundary conditions in this displacement field

BCs are

· At Node 1 ,

· At Node 2

· At Node 3

· At Node 4

Note that the DOF are independent and so we can uncouple the BCs (i.e., we can write the expressions for and separately. On substitution of these boundary conditions we get,

Similarily

On solving these simultaneous equations

To find the inverse of 4 x 4 matrix we need to partition the matrix as follows

we can observe that there is a simplicity in [A] matrix

From the above two equations

Similarily

Adding above two equations

So, the inverse matrix of will be

where

2.4 Shape functions using Lagrangian and Hermitian polynomials

Lagragian functions give polynomial functions that have C1 continuity in the field variable. That means, only the filed variable (displacement) will be continuous.

General form of Lagrangian function

In one-dimension

For example, for 5-noded bar element

These expressions indicate the shape functions of a 5-noded bar element. The shape functions are shown below

In two-dimensions the Lagrange polynomials will be as follows

For example for a 4-noded quadrilateral element, the linear shape functions will be written in terms of the natural coordinates in stead of

2.5 Shape functions using ‘Serendipity Concept’

2.6 Concept of iso-parametric, sub-parametric and super-parametric elements

Shape functions are interpolation functions used to get the values at the intermediate points from the values at the nodes. These functions can be used both for interpolating the field variables (displacement components like ) and can also be used to interpolate the geometry () . Based on the type of these two kinds of shape functions used in a particular element type, the elements are classified into iso-parametric, sub-parametric and super-parametric elements.

2.6.1 Iso-parametric elements

If same number of nodal values are used for interpolating the displacements and geometric coordinates, then such elements are called iso-parametric elements. For example, in the 4-noded quadrilateral element shown below, the geometry is defined by four nodes and displacements are also defined by 4-nodes. Similarily, a curved geometry requires 8 nodes to define the geometry. If the same number of nodes are used to define the displacements, then it becomes an iso-parametric element.

2.6.2 Sub-parametric elements

If less number of nodal values are used for interpolating the geometric coordinates, than the number of nodal values for interpolating the displacements, then such elements are called sub-parametric elements. Such elements are useful to model high stress variations, in simple geomtry. For example, in the 8-noded quadrilateral element shown below, the geometry is defined by four nodes and displacements are also defined by 8-nodes.

2.6.3 Super-parametric elements

If higher number of nodal values are used for interpolating the geometric coordinates, than the number of nodal values for interpolating the displacements, then such elements are called super-parametric elements. Such elements are useful to model curved geomtries with not high stress variations. For example, in the 8-noded quadrilateral element shown below, the geometry is defined by eight nodes and displacements are defined only by 4-nodes.

3 Strain-displacement relationship in discretized form

3.1 Review of strain displacement relationship

In theory of elasticity (or mechanics of deformable bodies) the strain displacement relationships are given as below

1D
	Axially loaded bar problems
	Transversily loaded beams
2D
	Plane strain and plane stress problems
3D
	3D solids

3.2 Strain -displacement relationship in discretized form of FEM

The continuous displacement field is discretized in terms of nodal values, as shown below.

Accordingly, the continuous form of strain-displacement relationship can be used to express strains in terms of nodal displacements as follows

3.2.1 Two-noded bar element

The LINEAR shape functions of two-noded bar element in local coordinate system is used for discretization as follows

where

Using the strain displacement relationship we get

Where

is the strain displacement matrix for twonoded bar element

3.2.2 Two-noded beam element

The relationship between the strain (curvature) and displacement (deflection) of beams is given by

The CUBIC shape functions of two-noded beam element in local coordinate system is used for discretization as follows

where

Using the strain displacement relationship we get

Where

, we can express this in terms of natural coordinates as follows

3.2.3 Three-noded triangular element OR Constant Strain Triangle (CST ) Element

The LINEAR shape functions of three-noded triangular element in natural coordinate system is used for discretization as follows

The above two componets of the displacement field can be compactly written as follows

where

푤here

Using the strain displacement relationship we get

Where

This is the strain displacement matrix for three noded triangular (CST) element

Determine the strain-displacement matrix for the two CST elements shown below (fig. 3)

Fig. 3

The strain displacement matrix of three-noded triangular (CST) element is

Where

And

For element #1

For element #2

3.2.4 Four-noded quadrilateral element

The LINEAR shape functions of four-noded quadrilateral element in natural coordinate system is used for discretization as follows

The above two componets of the displacement field can be compactly written as follows

where

Using the strain displacement relationship we get

Where

where [J] is called as Jacobian matrix [J] which is given by the concept of isoparametric formulation

3.3 Illustrative examples

3.3.1 CST Element

Show that the strain is constant at any point in a CST element. Thus, determine the constant strain for the element shown below, when

The strain displacement matrix for a CST element is given by

It can be seen that the matrix is independent of the coordinates of the interior point. So, it represents a state of constant strain at all points inside the triangle.

Here

Therefore,

The strain components of this CST element is given by

3.3.2 Four-noded quadrilateral element

Determine the strain at the point corresponding to for the four-noded quadrilateral element shown below, when

Determination of Jacobian matrix at the given point to

Now, the strain displacement matrix is given by

Finally, the strain at the given point is

4 Gauss quadrature

4.1 Concept of Gauss Quadrature

4.2 Examples

4.2.1 Line Integrals

Evaluate the following integrals using appropriate order of Gauss quadrature. Check the results using conventional method.

4.2.2 Line integrals involving shape functions

Evaluate for a quadratic bar element, using appropriate order of Gauss quadrature. Hint: The shape functions for 3-noded quadratic element are .

Since higest degree is 2 for each shape function, will have a degree of 4.

4.2.3

5 Stiffness matrix

The stiffness of the structral system also gets discretized into a finite sized matrix. This matrix relates the nodal displacements with the nodal forces. The expression for this matrix is dervied from Variational principle of minimization of the potential energy (functional).

5.1 Variational principle of minimization of the potential energy

5.1.1 Calculus principle:

A function has extreme value (maximum or minimum) when its first derivative with respect to variables is zero. The function will have a maximum value if the second derivative is negative and has a minimum value when the second derivative is positive. The first derivative of a function is called its variance.

A function of a function is termed as functional. A functional is said to have attained a stationary value when the first variance of the functional is zero.

In structural mechanics, the functional used is the potential energy of the structure or body. It can be showed that the minimization of this potential energy functional will lead to equailibrium equation.

5.1.2 Variational method in Finite Element Analysis

Step1: Formulate the Potential energy functional in terms of nodal displacements (unknown variables) at element Level

Potential energy = Strain energy stored in the body – Work done by the forces

Strain energy expression:

The strain energy density at a point in a deformed body is given by , Strain energy in the element will be obtained by integrating the volume times strain energy density over all points in the element

In discretized form the strain at a point is related to the nodal displacement as follows

and stress at a point is related to the strain at a point through the elasticity matrix [D]

Therefore the expression for strain energy will be

Expression for Workdone by external forces:

External forces on any element can be

i. Volume-proportional forces (or body forces)

ii. Area-proportional forces (or surface forces)

iii. Concentrated forces at the non-nodal points

iv. Concentrated forces at the nodal points

The body force in x-direction at a point is equal to the product of body force density at that point and differential volume around that point. Work done by body forces at a point is given by sum of work done in each of the three directions (x-, y-, z-). Thus, the total work done by body forces in an element is

But, the displacement field is discretized (or expressed) in terms of nodal displacements using shape function matrix [N] as follows

Thus, the work done will be

Where is equivalent nodal body force vector of element ‘e’

Similarly, for surface forces the equivalent nodal surface force vector is

Similarly, for concentrated forces on non-nodal points

For concentrated forces at nodal points, the work done is calculated at global level and not at element level

Step2: Condition for minimization of potential energy functional

If the total potential enegry of the structure is to be minimized then

where

For the work term

and so on..

5.2 Stiffness matrix for 2-noded bar element

5.2.1 Shape function matrix [N]

For 2-noded bar element the shape functions are

5.2.2 The strain displacement matrix [B_e]

The shape function matrix for 2-noded bar element ‘e’ is

5.2.3 The stiffness matrix [K_e]

The stiffness matrix will be given by

For one-dimensional case of bar element, the volume integral becomes a line integral as,

5.3 Stiffness matrix for quadratic bar element

The general formula for Stiffness matrix for any element is

For prismatic bar element

For homogeneous bar element

For 3-noded bar element, the strain-displacement matrix is given by

where

Where, the relation between cartesian coordinates and natural coordinates is given by linear relationship (as the geometry of the bar is straight line)

This satisfies the following boundary conditions also

Therefore,

Thus, the strain-displacement matrix will be

The stiffnesss matrix will therefore be

5.4 Stiffness matrix for beam element

5.4.1 Shape function matrix [N]

For 2-noded beam element the shape functions are

5.4.2 The strain displacement matrix [B]

The strain-displacement matrix for 2-noded bar element ‘e’ is

5.4.3 The Stress-strain relationship [D]

The strain for beam element is the curvature and the stress will be the bending moment M. So, the stress-strain relationship will be

5.4.4 The stiffness matrix [K_e]

The stiffness matrix will be given by

For one-dimensional case of beam element, the volume integral becomes a line integral as,

5.5 Stiffness matrix for CST element

5.6 Illustrative examples of determining [K]

5.6.1 For beam element

Q1. Determine the element stiffness matrices for the two elements (AB and BC) of the beam shown below

For element AB of L=6 m,

For element AB of L=4 m,

The global stiffness matrix will be

5.6.2 For CST element

Q1.For the CST element shown below, determine the shape function at the point P, and the stiffness matrix if thickness t=1m, E=65 GPa, and μ=0.3 for plane stress conditions

Shape function at point P will be same as the natural cordinates at this point

The negative area indicates that the point is outside the triangle

Stiffness matrix calculation

The strain displacement matrix for a CST element is given by

It can be seen that the matrix is independent of the coordinates of the interior point. So, it represents a state of constant strain at all points inside the triangle.

Here

Stiffness matrix

5.6.3 Element stiffness matrix for 4-noded element rectangular element

The element stiffness matrix of a 2D element is given by

When using one-point Gauss Quadrature, the above integral can be written as

where

Therefore,

Finally

5.6.4 Example#1

Determine the stiffness matrix for the 4-noded quadrilateral element shown below using one-point Gauss quadrature. Assume plane stress conditions with E and

Where

where

Therefore,

5.7 Assembling stiffness matrix

Each node is associated with n DOF. The DOF are numbered based in the same sequence of the node numbers. So, Based on the node numbers of an element, the element DOF can be determined and the element stiffness matrix will be assembled at that location

5.7.1 Example

Assemble the stiffness matrices shown for the four CST elements shown below. Assume the element stiffness matrices to be as given in the figure.

A diagram of a square with numbers and points

Description automatically generated

Element #1 connects between nodes 1-2-4. Thus, the DOF will be 1,2,3,4,7,8

Element #1 connects between nodes 1-4-5. Thus, the DOF will be 1,2, 7,8,9,10

Element #3 connects between nodes 2-3-4. Thus, the DOF will be 3,4,5,6,7,8

On adding the above three matrices we get

5.7.2 Example#2

Q. Assemble the stiffness matrices shown for the four CST elements shown below. Assume the element stiffness matrices to be as given in the figure.

A diagram of a triangle with green lines and dots

Description automatically generated

Solution:

5.8 Consistent Load vector

5.8.1 For quadratic bar element

For a uniformly distributed selfweight load

Because

5.8.2 For CST element under body forces

The general expression for consistent load vector for any element is obtained from the principle of virtual work as follows

Under body forces,

In 3D elements (not for CST element)

In 2D, the above expression can be rewritten as

Since CST element is 2D element, and self-weight is a body force, the second expression is used

Here, is the unit weight of the element. Gravity direction is assumed to be in negative y-direction.

The shape function matrix for CST element is

Therefore,

General formula for integration in natural coordinates of CST elements is

From this we can write that

5.8.3 For CST elements under surface forces

The general expression for consistent load vector for any element is obtained from the principle of virtual work as follows

Under surface forces,

For 3D elements (not for CST element)

For 2D elements (like CST element), the above expression can be rewritten as

Under traction forces on the free boundary,

5.8.3.1 Example: CST element under uniform pressure and along edge 1

All the integrations here are to be done along the edge1 connecting nodes 2&3. Noting that along edge 1, the value of will only exist

We also know that

Using he above formula for

5.8.3.2 Example#2

Derive the equivalent nodal force vector for CST element subjected to uniform pressures

A diagram of a triangle

Description automatically generated

5.9 Solution of FEM problems of Bar and beams

5.9.1 Condensation of internal nodes

To reduce the size of the matrices involved in solving the equilibrium equations, the internal nodes can be eliminated. This elimination is mathematically done so that the effect of the nodes presence can also be taken into account. Mathematically, this is equivalent to removing the rows and columns corresponding to the degrees of freedom of the internal nodes

· Firstly, the DOF of nodes to the eliminated are to be partitioned as follows

Here indicates the rows and columns of the nodes that are retained, and are the rows and columns of the nodes that are to be condensed. The matrix is the matrix containing the rows corresponding to the retained nodes and columns corresponding to eliminated nodes. The matrix is the matrix containing the column corresponding to the retained nodes and rows corresponding to eliminated nodes.

· To eliminate the rows and columns corresponding to the eliminated rows, we need to modify the matrix to take the effect of other three matrices as follows

This modified matrix will be the condensed stiffness matrix corresponding to the retained nodes. Note that the size of the matrix will be smaller than the original stiffness matrix

5.9.2 Bar problems

For the bar shown in Fig. 2a, determine the axial displacements at nodes 2 & 3, and support reactions at nodes 1 & 4. Take E= 200 Gpa, and P= 30 kN

Element#1:

Element#2:

Element#3:

Global stiffness matrix

Applying the BCs at DOF 1 and 4

Force vector is

Unknown Displacement vector will be

5.9.3 Beam problems

Solve the continuous beam shown in Fig. 3 to determine the rotation over the continuous support B and over support C.

Element#1:

Element#2:

Global stiffness matrix

Applying the BCs at DOF 1, 2, 3 and 5

Force vector is obtained by taking equaivalent nodal forces for each element and then assembling them as follows

For element#1

For element#2

The global force vector will be

On removing the rows corresponding to the BCs 1,2,3 and 5

Unknown Displacement vector will be

6 Finite Element software

In this chapter the following topics are covered

Commerical software packages for Finite Element Analysis

Basic features of a general purpose Finite Element Software

Desired features of pre- and pos-processing

a. Auto-mesh generation

6.1 Software packages avaiable

The following different FE softwares are available

1. Structural Design Language (STRUDEL) by MIT, USA

2. NASA Structural Analysis (NASTRAN) by NASA, USA

3. Non-linear Incremental Structural Analysis (NISA) by E Ramm Germany

4. Engineering Analysis Systems (ANSYS), by Swanson, sweden

5. Structural Analysis Program (SAP), by Dr. Wilson, UC, USA

6. Advanced Dynamic Incremental Nonlinear Analysis (ADINA)

7. ABAQUS FEA software by Dessault Systems, France

8. Open Source Earthquake Engineering Software (OPENSEES), UC, USA

Improvements leads to newer versions of the software. These improvements include

Increasing the variety and capability of elements

Provision of using different types of elements at a time

Addition of various kinds of dynamic analysis

Addition various types of geometric and material non-linearities

Addition of optimization modules

Addition of user-defined modules and scripting interfaces

Addition of user-firendly interfaces

Addition of parallel processing and other modules

6.2 Basic features of general purpose FE software

Input modules

b. Geometric data

c. Sectional data

d. Material data

e. Meshing data

f. Loading data

g. Boundary conditions

h. Interactions between the various components

i. Type of analysis and the parameters governing the analysis

j. Output required

Pre-processing modules

k. Initializing the variables

l. Element level analysis to obtain element stiffness matrices and load vectors. This requires looping over each gauss point (or integration point) to obtain corresponding stiffness matrix and adding them with appropriate weights to get the element stiffness matrix

m. Assembling each element stiffness matrix in global stiffness matrix, using the element DOF

n. Assembling the element force vectors into the global force vector using element DOF

o. Imposing the boundary conditions by penalty method or elimination method

Solver module

p. The simulateous algebraic equations obtained after imposing BCs are solved by numerical methods like Gauss-siedal method, etc

Post-processing module

q. The global displacements are transformed into element level displacements to get deformed shapes

r. The gauss point level stress & strain analysis is done and plotted as contour plots, vector plots, etc

s. Error analysis is done

6.3 Desirable features of a FEA package

· Excellent pre- and post-processing modules

· Availability of all types of elements, 1D, 2D, 3D, plate and shell elements

· Ability to handle different types of load distributions like concentrated, unformly distributed, parabolic and other user-defined variations of load distribution, both in time and space.

· Ability to consider different types of loads, like gravity, internal stresses, initial strains, lack of fit, temperature stresses, centrifugal loads, traction loads

· Ability to apply all types of boundary conditions

· Ability to handle problems with huge DOF by efficient use of memory, parallel processing, frontal solution techniques, band width reduction algorithms

· Ability to performa variety of analysis like dynamic analysis, buckling analysis, frequency analysis, no-linear analysis, etc

· Ability to performa optimization algorithms

· Avaialbility of network licensing for marge number of users

· Availability of design modules for specific requirements. For example, general purpose FE packages can be customized for Civil Engineering design, Mechnical designs, as per the corresponding design standards applicable in that locality

6.3.1 Auto-mesh generation feature of FEA software

Auto mesh generation in the Finite Element Method (FEM) is a crucial step in setting up a model for simulation. It involves automatically dividing the domain into smaller, simpler regions (elements) to facilitate numerical analysis. Here’s an overview of how auto mesh generation works, common techniques, and considerations.

1. Mesh Generation Techniques

The primary goal in mesh generation is to discretize a complex geometry into elements (like triangles or quadrilaterals in 2D, and tetrahedrons or hexahedrons in 3D) to allow for accurate simulations. Techniques used include:

Structured Meshes: These consist of regular patterns, often used for simple geometries (like grids in square or rectangular domains). While easy to generate, they are less adaptable for complex shapes.
Unstructured Meshes: These consist of elements of varying size and shape, allowing for better adaptability to complex geometries. Unstructured meshes are generated using algorithms that adjust the element shapes and sizes based on geometry and desired resolution.
Hybrid Meshes: Some meshes combine structured and unstructured regions. For example, the core of a geometry may use a structured mesh, while boundaries use an unstructured one.

2. Meshing Algorithms

Various algorithms are available for automatic mesh generation, including:

Delaunay Triangulation: A popular method for generating unstructured meshes in 2D and 3D. It creates triangular or tetrahedral meshes and ensures good element quality with minimum angle optimization.
Advancing Front Method (AFM): This method starts from the boundary and "advances" inward, generating elements as it moves. It’s useful for adaptive meshing and produces well-structured elements in complex geometries.
Octree and Quadtree Partitioning: This technique divides the domain recursively into smaller quadrants (2D) or octants (3D), which can be refined based on geometric complexity. It’s particularly useful for adaptive meshing.

3. Adaptive Mesh Refinement (AMR)

Adaptive Meshing: AMR is a technique used to refine the mesh in regions where higher accuracy is required. In FEM, AMR can help by increasing the element density in areas with large gradients (like stress concentrations or boundary layers) while keeping a coarser mesh in less critical areas.
Error Estimation and Refinement: The solution is first computed with a coarse mesh, then error estimations guide the refinement process. This can be achieved either by refining the mesh globally or selectively in specific regions.

4. Meshing Considerations in FEM

Element Quality: Poor-quality elements (e.g., highly skewed or distorted shapes) can lead to inaccurate solutions and slower convergence in FEM analysis. Mesh quality is typically improved by enforcing geometric criteria during generation.
Boundary Conformity: Ensuring that the mesh conforms well to the geometry’s boundaries is essential for accurate results, particularly in cases involving complex or curved surfaces.

Resolution: Choosing the right resolution (element size) is essential; smaller elements lead to higher accuracy but increase computational cost. Adaptive meshing can help balance this by refining only the necessary regions.

Node	X	Y	Z
1	-1	-1	-1
2	1	-1	-1
3	1	1	-1
4	-1	1	-1
5	-1	-1	1
6	1	-1	1
7	1	1	1
8	-1	1	1
9	0	-1	-1
10	1	0	-1
11	0	1	-1
12	-1	0	-1
13	0	-1	1
14	1	0	1
15	0	1	1
16	-1	0	1
17	-1	-1	0
18	1	-1	0
19	1	1	0
20	-1	1	0

Node	X	Y	Z
1	-1	-1	-1
2	1	-1	-1
3	1	1	-1
4	-1	1	-1
5	-1	-1	1
6	1	-1	1
7	1	1	1
8	-1	1	1
9	0	-1	-1
10	1	0	-1
11	0	1	-1
12	-1	0	-1
13	0	-1	1
14	1	0	1
15	0	1	1
16	-1	0	1
17	-1	-1	0
18	1	-1	0
19	1	1	0
20	-1	1	0

Friday, 27 June 2025

FORMULATION OF 20-NODED HEXAHEDRAL ELEMENT

INTRODUCTION

Figure 1: 20-Noded Hexahedral Element showing nodes

SHAPE FUNCTION

ELEMENT STIFFNESS MATRIX

Strain Displacement Matrix

Constitutive Matrix

Jacobian Matrix

Numerical Integration

VERIFICATION OF CODE

Static analysis of a straight cantilever beam

Static analysis of a tapered cantilever beam

Monday, 23 June 2025

1.3.2 Discritize the continua

2.3.6 Shape functions for 3-noded triangular element in natural coordinate system

4.2.3

5.1 Variational principle of minimization of the potential energy

5.8.3.1 Example: CST element under uniform pressure and along edge 1

5.8.3.2 Example#2

1.3.2

Discritize the continua

Node	X	Y	Z
1	-1	-1	-1
2	1	-1	-1
3	1	1	-1
4	-1	1	-1
5	-1	-1	1
6	1	-1	1
7	1	1	1
8	-1	1	1
9	0	-1	-1
10	1	0	-1
11	0	1	-1
12	-1	0	-1
13	0	-1	1
14	1	0	1
15	0	1	1
16	-1	0	1
17	-1	-1	0
18	1	-1	0
19	1	1	0
20	-1	1	0