ಓಂ ಶ್ರ ೀ ಗಣೇಶಯನಮಃ

PRN 5002: FINITE ELEMENT ANALYSIS

 Introduction

The finite element method (FEM) also referred as finite element analysis (FEA), has

become a powerful tool for the numerical solution of a wide range of engineering problems

Example:

• Deformation and stress analysis of automotive, aircraft, building and bridge structures.

• Field analysis of heat flux, fluid flow, magnetic flux, seepage and other flow problems.

• With advances in CAD systems, several alternative configurations can be tested on a

computer before the first prototype is built.

• FEM discretizes a large problem into smaller and simpler parts, called finite elements.

• The material properties and governing relationships are considered over these elements

and expressed in terms of unknown values at element corners, that results in a set of

elemental equations.

• The elemental equations are then assembled into a global system of equations that

models the entire problem.

• FEM then uses calculus variational methods to solve the global system of equations.

• Solution of these equations gives us the approximate behavior of continuum.

• Stress analysis
𝝐 = 𝑫/ 𝝈

• Temperature analysis
𝝐=α∆𝑻
 Introduction

V – Volume
S – Surface
T – Traction
f – Distributed force per unit
volume

Fig . Three-dimensional body

• V – Volume
• S – Surface
• T – Traction
• f – Distributed force per unit volume
𝑇
𝑓 = 𝑓𝑥, 𝑓𝑦, 𝑓𝑧
• Deformation of a point ‘x’ is given by
three components of its displacement
𝑇
𝑢 = 𝑢, 𝑣, 𝑤

• Surface traction T given by

𝑇
𝑇 = 𝑇𝑥, 𝑇𝑦, 𝑇𝑧

• P – Load acting at a point ‘i ’ represented by three components

𝑇
𝑃𝑖 = 𝑃𝑥, 𝑃𝑦, 𝑃𝑧
The stresses acting on the
elemental volume 𝑑𝑉 are shown
in Fig.

Fig. Equilibrium of elemental volume

For equilibrium of elemental volume shown in above Fig. the forces on the faces derived by
multiplying the stresses by corresponding areas.

With σ 𝐹𝑥 = 0, σ 𝐹𝑦 = 0, σ 𝐹𝑧 = 0 and 𝑑𝑉 = 𝑑𝑥, 𝑑𝑦, 𝑑𝑧 the equilibrium equations are

derived as
𝜕𝜎𝑥 𝜕𝜏𝑥𝑦 𝜕𝜏𝑥𝑧
+ + + 𝑓𝑥 = 0
𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕𝜏𝑥𝑦 𝜕𝜎𝑦 𝜕𝜏𝑦𝑧

+ + + 𝑓𝑦 = 0
𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕𝜏𝑥𝑧 𝜕𝜏𝑦𝑧 𝜕𝜎𝑧

+ + + 𝑓𝑧 = 0
𝜕𝑥 𝜕𝑦 𝜕𝑧

Where 𝜎𝑥 , 𝜎𝑦 , 𝜎𝑧 are normal stresses and 𝜏𝑥𝑦 , 𝜏𝑥𝑧 , 𝜏𝑦𝑧 are shear stresses and are represented
in vector form as;
𝑇
σ = 𝜎𝑥 , 𝜎𝑦 , 𝜎𝑧 , 𝜏𝑥𝑦 , 𝜏𝑥𝑧 , 𝜏𝑦𝑧
 Strain Displacement Relations

Strains are represented in vector form as

𝑇
𝜖 = 𝜖𝑥 , 𝜖𝑦 , 𝜖𝑧 , 𝛾𝑥𝑦 , 𝛾𝑥𝑧 , 𝛾𝑦𝑧

where 𝜖𝑥 , 𝜖𝑦 , 𝜖𝑧 are normal strains and 𝛾𝑥𝑦 , 𝛾𝑥𝑧 , 𝛾𝑦𝑧 are shear strains.

The strain relations for small deformations given as;

The below Fig. gives the deformation of dx-dy face for small deformations.
Strain Displacement Relations

Fig. Deformed element surface

Stress-Strain Relations
For an elemental cube inside the body, Hook’s law gives (for linear elastic material)

where shear modulus

(or modulus of rigidity)
G is given by;
Stress-Strain Relations are expressed as;

𝝈𝑥 𝝐𝑥
𝝈𝑦 𝝐𝑦
𝝈𝑧 𝝐𝑧
τ𝑥𝑦 γ𝑥𝑦
τ𝑥𝑧 γ𝑥𝑧
τ𝑦𝑧 γ𝑦𝑧

𝝈 = 𝑫𝝐
where D is the material matrix given by;
 𝝈 = 𝑫𝝐
where D is the material matrix given by;

One dimension: In one dimension, we have normal stress 𝜎 along x and the corresponding
normal strain 𝜖. Stress-strain relations are simply;

𝝈 = 𝑬𝝐
Two dimensions:
In two dimensions, the problems are modelled as plane stress and plane strain.

Plane Stress:
A thin planar body subjected to in-plane loading on its edge surface is said to be in plane
stress.
A ring press fitted on a shaft (Fig.) is an example.

Here stresses 𝜎𝑧 , 𝜏𝑥𝑧 and 𝜏𝑦𝑧

are set as zero.
Then Hook’s law relations given as;

And stress-strain relations (for Plane stress) are given as;

which is used as

𝝈 = 𝑫𝝐
Plane Strain:
If a long body of uniform cross section is subjected to transverse loading along its length, a
small thickness in the loaded area (as shown in Fig.) cab be treated as subjected to plane
strain.

Here 𝜖𝑧 , 𝛾𝑥𝑧 , 𝛾𝑦𝑧 are taken as zero. Stress 𝜎𝑧 may not be zero in this case.
The stress-strain relations (for Plain strain) can be written as;

which is used as
𝝈 = 𝑫𝝐

where 𝑫 is a (3 x 3) matrix, which relates three stress with three strains.

• Our problem is to determine the
displacement u of the body shown in Fig.,
by satisfying the equilibrium equations.

• Note that stresses are related to strains,

which, in turn, are related to displacements.

• This leads to requiring solution of second-

order partial differential equations.

• Solution of this set of equations is generally

referred as an exact solution.
• Such exact solutions are available for simple geometries and loading conditions.

• For problems of complex geometries and general boundary and loading conditions,
obtaining such solutions is an almost impossible task.

• Approximate solution methods usually employ potential energy or variational methods.

Potential Energy and Equilibrium; The Rayleigh-Ritz Method

The total potential energy (Π) of an elastic body, is defined as the sum of total strain energy
(U) and the work potential:

Π = Strain energy (U) + Work potential (WP)

For linear elastic materials, the strain energy per unit volume in the body is

𝟏 𝑻
𝝈 𝝐
𝟐

For the elastic body shown in Fig. 1 the total strain energy (U) is given by
The work potential WP is given by

The total potential for the general elastic body

shown in Fig. is

This potential energy Π can be used for finding an approximate solution.

The Rayleigh-Ritz method involves the construction of an assumed displacement field as

polynomial functions.
Example:
Fig. shows the a system of springs.

The total potential energy is given by

Since
𝛿1 = 𝑞1 − 𝑞2
𝛿2 = 𝑞2
𝛿3 = 𝑞3 − 𝑞2

𝛿4 = −𝑞3

we have

For equilibrium of this three degrees of freedom system, we need to minimize with respect
to 𝑞1 , 𝑞2 , 𝑞3 .
These three equations are given by

These equilibrium equations cab be put in the form of

Kq = F as follows:
From this we can write

Which is precisely the set of equations.

Example 2:

The potential energy equation for the linear elastic one dimensional rod (shown in above
Fig.), with body force neglected is given by;

Strain energy

work potential = Point load X Displacement

Let us consider a polynomial function

This must satisfy u = 0 at x = 0 and u = 0 at x = 2.

Thus,

Hence,

u=
 Then

where
Note here that an exact solution is
obtained if piecewise polynomial
interpolation is used in the
construction of u.

The finite element method

provides a systematic way of
constructing the polynomials.
 Galerkin’s Method
• Galerkin’s method uses the set of governing equations in the development of an integral
from.

• It is usually represented as one of the weighted residual methods.

• For the one-dimensional rod considered in the previous example, the governing equation
is the differential equation.
Galerkin’s Method

For three-dimensional stress analysis the Galerkin’s equation is given as (which is

“variational form” or “weak form” for three-dimensional stress analysis)

Where 𝝓 is an arbitrary displacement consistent with the specified boundary conditions of

u.
Galerkin’s Method
Example 3: Consider the problem of Example 2 and solve it by Galerkin’s approach.

The equilibrium equation is given by

Multiplying this differential equation by 𝝓 , and integrating by parts, we get

Where 𝝓 is zero at x = 0 and x = 2.
𝒅𝒖
𝑬𝑨 is the tension in the rod, which takes a jump of magnitude 2 at x = 1.
𝒅𝒙
Thus

Now we use the same polynomial for u and 𝝓.

If u1 and 𝝓1 are the values at x = 1, we have
Substituting these and E = 1, A = 1 in the previous integrals yields

This is to be satisfied for every 𝝓1.

Then we get
u1 = 0.75
Assignment
1.

2.
3.
 Matrix Algebra and Gaussian Elimination
MATRIX ALGEBRA
The study of matrices here is largely motivated form the need to solve systems of
simultaneous equations of the form

where 𝒙𝟏 , 𝒙𝟐 ,......, 𝒙𝒏 are the unknowns. The above equations can be conveniently
expressed in matrix form as
𝐀𝐱 = 𝐛
 𝐀𝐱 = 𝐛

where A is a square matrix of dimensions (n x n). And x and b are vectors of dimensions (n
x 1), given as

The matrix A is also denoted as [A]. An element located at the ith row and jth column of A
is denoted by aij.
The product of a (m x n) matrix A and an (n x p) matrix B results in an (m x p) matrix. That
is,

It should be noted that AB ≠ BA; in fact, BA may not even be defined, since the number of
columns of B may not equal the number of rows of A.
Transposition
If A = [aij], then the transpose of A, denoted as AT,
is given by AT = [aji]. The rows of A are the columns of AT.
For example, if

then

In general, if A is of dimension (m x n), then AT is of dimension (n x m).

The transpose of a product is given as the product of the transpose in reverse order:
Differentiation and Integration

The components of a matrix do not have to be scalars; they may also be functions. For
example,

In this regard, matrices may be differentiated and integrated. The derivative (or integral)
of a matrix is simply the derivative (or integral) of each component of the matrix. Thus,
Differentiation and Integration
If A be as (n x n) matrix of constants, and 𝑥 = [𝑥1 , 𝑥2 , … … . , 𝑥𝑝 ]𝑇 be a column vector of n
variables. The then derivative of Ax with respect to variable 𝑥𝑝 is given by

The derivative of Ax with respect to 𝑥𝑝 yields the pth column of A.

Square Matrix

A matrix whose number of rows equals the number of columns is called a square matrix.

Diagonal Matrix

A diagonal matrix is a square matrix with nonzero elements only along the principal
diagonal.
For example,
Identity Matrix

The identity (or unit) matrix is a diagonal matrix with 1’s along the principal diagonal.
For example,

If I is of dimension (n x n) and x is an (n x 1) vector, then

Ix = x
Symmetric Matrix
A symmetric matrix is a square matrix whose elements satisfy

Or equivalently,

That is, elements located symmetrically with respect to the principle diagonal are equal.

For example,
Upper Triangular Matrix

An upper triangular matrix is one whose elements below the principal diagonal are all
zero.

For example,
Determinant of a Matrix

The determinant of a square matrix A is a scalar quantity denoted as det A. The

determinants of a (2 x 2) and a (3 x 3) matrix are given by the method of cofactors as
follows.
Matrix Inversion

Consider a square matrix A. If det A ≠ 0, then A has an inverse, denoted by A-1. The inverse
satisfies the relations

If det A ≠ 0, then we say that A is nonsingular.

If det A = 0, then we say that A is singular.

The inverse of a square matrix is given as

For example, the inverse of a (2 x 2) matrix A is given by

Matrix Inversion
Find the adjoint of the matrix A=

We will first evaluate the cofactor of every element

Therefore
Quadratic Forms

Let A be an (n x n) matrix and x be an (n x 1) vector. Then, the scalar quantity

Is called a quadratic form, since upon expansion we obtain the quadratic expression
Example: The quantity

can be expressed in the matrix form as

Eigenvalues and Eigenvectors

Let A be a (2 x 2) matrix

Let’s multiply A by the vector y=

i.e Ay=

Ay=2y → Ay=λy
where
2 is called as an eigenvalue (λ) and y is called as an eigenvector
• This vector (y) is special! When we multiply this vector by the matrix A, we get the vector
back, although it is scaled by a number.

• This vector (y) is called an eigenvector of A and the scaling number is

an eigenvalue associated with the eigenvector.

• Remember what happened when we multiplied the matrix A with an eigenvector of y?

• The eigenvector was unchanged except for a scaling factor.

• If we had multiplied any other vector by A, the vector would have changed.

• It's only these special eigenvectors that remain the same. An equation summarizing this is
Ay = λy where λ is the eigenvalue associated with the eigenvector y.
To find a nonzero eigenvector and the corresponding eigenvalues λ,
we rewrite

I represents the identity matrix, which has 1 along the main diagonal and 0 everywhere else.

λy is the same as λIy.

Then,
Ay = λy is the same as Ay = λIy.
So,

We want to find eigenvectors y and eigenvalues λ. This last equation is true if y = 0 but having
a vector equal to zero is not interesting.
We want to find eigenvectors y and eigenvalues λ. This last equation is true if y = 0 but
having a vector equal to zero is not interesting.

We want solutions for y not equal to 0.

Nonzero solution for y will occur when (A - λI) is a singular matrix or

det (A - λI) = 0 (i.e., )

This equation is called as the characteristic equation

Example:
λ 0
Consider the matrix λI =
0 λ

The characteristic equation is

det(A - λI) = 0 →

Which yields

Solving the above equation we get

The above are the two eigenvalues

To find the eigenvectors
𝐲 𝟏 = 𝑦11 , 𝑦21 𝑇

Corresponding to the eigenvalue λ1, we substitute λ1= 3

λ 0 3 0
(A - λI)y = 0 where and λI = =
0 λ 0 3

This yields the equation

We may now normalize the eigenvector by making 𝐲 𝟏 a unit vector.

This is done by setting 𝑦21 = 1, resulting in

𝐲 𝟏 = 𝑦11 , 𝑦21 𝑇 = [2.236, 1]

Dividing 𝐲 𝟏 by its length i.e. 2.2362 + 12 yields

𝐲 𝟏 = 0.913, 0.408 𝑇

Now 𝐲 𝟐 is obtained in a similar manner by substituting

λ2 = 9 in to (A - λI)y = 0
After normalization
𝐲 𝟐 = −0.408, 0.913 𝑇
 4 −2.236 0.913
𝐀𝐲 =
−2.236 8 0.408

4 x 0.913 + (−2.236 x 0.408)

−2.236 x 0.913 + (8 x 0.408)

3.652 − 0.9122
−2.0414 + 3.264

2.739
1.224

0.913 2.739
𝐀𝐲 = 𝛌𝐲 = 3 =
0.408 1.224
Positive Definite Matrix

A symmetric matrix is said to be positive definite if all its eigenvalues are strictly positive

(grater then zero).

In the previous example, the symmetric matrix

had eigenvalues λ1 = 3 > 0 and λ2 = 9 > 0 and, hence, is positive definite matrix.
 Gaussian Elimination
Consider a linear system of simultaneous equations in matrix form

Ax = b
Where A is (n x n) and b and x are (n x 1).

The unique solution for x as

x = 𝐀−𝟏 b

• Here the construction of 𝐀−𝟏 , by the cofactor approach, is computationally expensive and
subjected to round-off errors.

• Instead, an elimination scheme is better.

• The powerful Gaussian elimination approach for Ax = b is discussed here.
Gaussian Elimination is the method used for solving simultaneous equation by successively
eliminating unknowns.

Consider the simultaneous equations

The equations are labeled as I, II, III.

It is required to eliminate 𝒙𝟏 from II and III.

We have (from Eq. I) 𝒙𝟏 = +2 𝒙𝟐 - 6 𝒙𝟑

Substituting for 𝒙𝟏 into Eqs. II and III yields

Notice the zeroes below the main diagonal in column 1, representing the fact that 𝒙𝟏 has
been eliminated from Eqs. II and III.

The superscript (1) on the labels in the above equations denotes the fact that the equations
have been modified once.

Now proceed to eliminate 𝒙𝟐 from Eq. III.

1
For this, we subtract times II from III.
6
The resulting system is
The coefficient matrix on the left side of above equations is upper triangular.
1
This solution now is virtually complete, since the last equation yields 𝒙𝟑 =
5
4
Upon substituting 𝒙𝟑 into the second equation yields 𝒙𝟐 =
5
1
Then 𝒙𝟏 = from the first equation.
5
1 4 1
i.e. 𝒙𝟏 = , 𝒙𝟐 = , and 𝒙𝟑 =
5 5 5

This process of obtaining the unknowns in reverse order is called back-substitution.

These operations can be expressed more concisely in matrix form as follow:

Working with the augmented matrix [A, b], the Gaussian elimination process is

→ →

Which, upon back-substitution, yields

1 4 1
𝒙𝟑 = , 𝒙𝟐 = , 𝒙𝟏 = ,
5 5 5
General Algorithm for Gaussian Elimination
Let the simultaneous equations

Which can be restated as

Let us consider the start of step 1, with A and b written as follows:

The idea at step 1 is to use equation 1 (the first row) in eliminating 𝒙𝟏 from remaining
equations. Denote the step number as a superscript set in parentheses.

The reduction is carried out for all the elements in the shaded area.
The elements in rows 2 to n of the first columns are zeros since 𝒙𝟏 is eliminated.
At the start of the step 2, thus we have

The start of the step k

After (n – 1) steps, we get

The superscripts are for the convenience of presentation. Drop the superscripts for
convenience, and the back-substitution process is give by
And then

This completes the Guass elimination algorithm.

Symmetric Banded Matrices

In a banded matrix, all of the nonzero elements are contained within a band; outside of the
band all elements are zero.
The stiffness matrix that will come across in this FEA is a symmetric and banded matrix.

Consider an (n x n) symmetric banded matrix

where nbw is called the

half-band width
Since only the nonzero elements need to be stored, the elements of this matrix are
compactly stored as
(n x nbw) matrix

• The principal diagonal or 1st

(main) diagonal of previous
matrix is the first column of this
matrix.
• In general, the pth diagonal of
previous matrix is stored as the
pth column of this matrix.
Skyline Solution

If there are zeroes at the top of a column, only the elements starting from the first nonzero
value need to be stored.

The line separating the top zeroes from the first nonzero element is called the skyline.

For example:
For efficiency, only the active columns need be stored.
These can be stored in a column vector A and a diagonal pointer vector ID as
Assignment
1)

2)
3)