Basics of FEM: Prof. S. V. Kulkarni

Indian Institute of Technology Bombay
Basics of FEM
Prof. S. V. Kulkarni
Department of Electrical Engineering
Outline
Introduction
FEM Procedure
Discretisation of the Domain
Approximation of the Solution
Assembly of the System
Boundary Conditions and Solution of the Final System
Properties of FEM Matrices

Introduction
Introduction
Computational Electromagnetics:
Finite difference method
Finite element method (FEM)
Boundary element method (BEM)
Method of moments (MoM)
Meshless methods: Particle-in-cell, Petrov Galerkin etc.
Finite Element Analysis:

Variational Procedure
Galerkins Method

Introduction
Illustrative Example: Parallel Plate Capacitor

=10
2 = 0
|y =1, 0<x <1 = 10

P:
|
y = 0 , 0 <x <1 = 0

(0, 1)

on
(0, 0)
(1)
x
(1, 0)
=0
(1, 1)

FEM Procedure
Variational Procedure
Minimisation or Maximisation of a Functionalfunction of

functions
Functional for the example problem:
E =
1
2
||2 d
(2)
Physical Significance of functional: Energy of the system

For a single dielectric system, the energy minimisation and
the corresponding solution is not influenced by and hence it
can be dropped from the expression

FEM Procedure
Overview of Discretisation
FEM Problem domain

divided into a large number of
small elements/sub-domains
Type of element:
Geometry of problem
Shape of element:
2DTriangular, Rectangular
and 3DCubic, Tetrahedral,
Prismatic
Nodal versus Vector
Size of element: specific to

geometry requirements
(x,3 y3)
3e
e
1
(x,1 y1)
(x,2 y2)
2e

FEM Procedure
Element Size Illustration: Arbitrary Boundary

FEM Procedure
Example Problems Domain Discretisation

7
3
6
1
3
5
3
1
2
2 6
3
4
2

FEM Procedure
Approximation of the Solution
Overview of Solution Approximation
Analytical solution may not exist for many problems

Piece-wise linear polynomial approximation: Polynomial
solution over each element
Examples:
3D Elemental: e = a + bx + cy + dz + exy + gyz + hxz + ixyz
2D Elemental: e = a + bx + cy + dxy
Global approximation:
e
e
Rayleigh-Ritz method (variational procedure)
(3)

FEM Procedure
Overview of Linear System Assembly

Nodal method: 2D Linear
polynomial approximation
(x,3 y3)
3e
e1 = a + bx1 + cy1
e2 = a + bx2 + cy2
(4)
e3 = a + bx3 + cy3

1
a 1 x1 y1 e1
e 1
1

b = 1 x2 y2 e2 (5)(x,1 y1)

c
1 x3 y3 e
3
We shall, henceforth, drop the on e s for simplicity and
understand that refers to the approximate solution.
(x,2 y2)
2e

FEM Procedure
2D Linear System Assemly (1 of 2)

1 x1 y1 e1
a

e = 1 x y b = 1 x y 1 x2 y2 e2
e

3
1 x3 y3
c
(6)
Ni (x , y )ei
(7)
i =1
where,
1
[(x2 y3 x3 y2 ) + (y2 y3 )x + (x3 x2 )y ]
2
1
[(x3 y1 x1 y3 ) + (y3 y1 )x + (x1 x3 )y ]
N2 (x , y ) =
2
1
[(x1 y2 x2 y1 ) + (y1 y2 )x + (x2 x1 )y ]
N3 (x , y ) =
2
N1 (x , y ) =
(8)

FEM Procedure
2D Linear Assembly (2 of 2)
and,
= Area of the elemental triangular element =
1 x1
1
1 x2
2
1 x3
y1
y2
y3
(9)
Ni (x , y ) has the property that,

Ni (xj , yj ) = ij
(10)
1 i = j
ij =
0 i j
(11)
where,

FEM Procedure
Functional Approximation (1 of 3)
Using equations 3 and 7 in equation 2, we can express the
functional in terms of the (unknown) nodal potentials as:
E =
e
E =
1
2
1
E =
2
1
2
Ni (x , y )ei
{Ni (x , y )}ei
e
de
2
de
(13)
i =1
N1 (x , y )e1 + N2 (x , y )e2 + N3 (x , y )e3

e
(12)
i =1
3
de
(14)

FEM Procedure
For any vector a, a a = |a |2
E =
1
2
N1 (x , y )e1 + N2 (x , y )e2 + N3 (x , y )e3

e
N1 (x , y )e1 + N2 (x , y )e2 + N3 (x , y )e3 de

(15)
E =
e
1
2
ei Ni (x , y ) Nj (x , y )ej de
e
i =1 j =1
where is the elemental domain
(16)

FEM Procedure
1
E =
2
3
e
i =1 j =1
Ni (x , y )
ei
Nj (x , y )de ej
(17)
Ni Nj de
aije
(18)
e
3
Ae
i =1 j =1
e
e
e
a13
a11 a12
e
e
e
a22
a23
aije = a21
e
e
e
a31
a32
(19)
a33
A is referred to as the elemental stiffness matrix. The elemental

energy can, thus, be represented as,
E e = e T Ae e
where e =
e1 e2 e3
and E =
(20)
E e.

FEM Procedure
Ae Example Calculations (1 of 2)
For example,
N1 N1 de
e
a11
=
(21)
where, from equation 8,
N1 =
1
(x2 y3 x3 y2 ) + (y2 y3 )x + (x3 x2 )y
2
1
[(y2 y3 )x + (x3 x2 )y]
(22)
2
1
= N1 N1
de =
(y2 y3 )2 + (x3 x2 )2
42
e
a11
(23)
1
(y2 y3 )2 + (x3 x2 )2
4
(24)

FEM Procedure
Ae Example Calculations (2 of 2)
Similarly, we have,
1
[(y2 y3 )(y3 y1 ) + (x3 x2 )(x1 x3 )]
4
(25)
1
[(y2 y3 )(y1 y2 ) + (x3 x2 )(x2 x1 )]
= N1 N3 =
4
(26)
e
a12
= N1 N2 =
e
a13

FEM Procedure
Global Connectivity Matrix
e
e
e
e
e
e
e
e
1
2
3
4
5
6
7
8
1e
1
1
2
2
4
4
5
5
2e
5
2
6
3
8
5
9
6
3e
(27)

FEM Procedure
Assembly of Global Stiffness Matrix (1 of 3)
The global stiffness matrix is of the form:
A11 A12 A13

A21 A22 A23
A = A31 A32 A33

.
..
..
.
A91 A92 A93
A19
A29
A39
..
.
A99
(28)

FEM Procedure

Example elements:
7
3
5
3
2 6
3
4
2
i = 3, 6, 7, 8, 9
e2
e3
e4
A22 = a22
+ a11
+ a11
e4
A23 = a12
= A32
..
.
3
2
A1i = Ai1 = 0
e2
A12 = a12
= A21
8
e1
e2
A11 = a11
+ a11
e1
e2
e3
A55 = a22
+ a33
+ a33
2
e6
e7
e8
+ a22
+ a11
+ a11
e1
e2
A51 = a21
+ a31
= A15
and so on

FEM Procedure
Thus, the global stiffness matrix is:

A=
1
2
3
..
.
9
1
2
3
e2
a e1 + a e2
a12
0
11
11
a e2
e3
e2
e4
e4
a22 + a11 + a11 a12
21
e4
e4
0
a21
a22
..
..
.
.
0
0
0
0
(29)
..
e8
e7
a22 + a33
0

FEM Procedure
Approximate Functional Expression
Thus, equation 17 for the functional can now be expressed as,
E =
1 T
A
2
(30)
where, A is the matrix in equation 29 and is a vector of all nodal

potential values,
= 1 2 3 4 5 6 7 8 9
(31)

FEM Procedure
Gradient of Functional is Zero
E = 0
E
1
E

. 2 = 0
..

(32)
(33)
E
9
A11 1 + A12 2 + A13 3 + + A19 9 = 0

A21 1 + A22 2 + A23 3 + + A29 9 = 0
..
.
A91 1 + A92 2 + A93 3 + + A99 9 = 0
A = 0
(34)
(35)

FEM Procedure
Incoporating Boundary Conditions: Method 1
The boundary conditions are:
1 = 2 = 3 = 0
7 = 8 = 9 = 10
This specifies six of the nine unknowns in equation 34.
The remaining three unknowns can be determined by
substituting these six values into any three of the nine
equations in (34) and solving the resulting set of linear
equations.
(36)

FEM Procedure
Equation 34 now becomes the reduced order system as shown

below:
A44 4 + A45 5 + A46 6 = 10(A47 + A48 + A49 )
A54 4 + A55 5 + A56 6 = 10(A57 + A58 + A59 )
A64 4 + A65 5 + A66 6 = 10(A67 + A68 + A69 )
(37)
A33 31 = b31
(38)
Symmetric

FEM Procedure

Alternatively, the boundary conditions can be applied without
reducing the size of A as follows:
1 = 0
2 = 0
3 = 0
A41 1 + A42 2 + A43 3 + + A49 9 = 0
A51 1 + A52 2 + A53 3 + + A59 9 = 0
A61 1 + A62 2 + A63 3 + + A69 9 = 0
7 = 10
8 = 10
9 = 10
A99
Unsymmetric
91 = b91
(39)
(40)

FEM Procedure
Linear System Solution
Matrix A of the order of 102 : Direct inversion

Matrix A of the order of 103 to 104 : Direct solution techniques
such as Gaussian elimination, LU decomposition or LDU
decomposition
For linear systems of the order of 105 or higher:
Computationally efficient to use only iterative solution
techniques such as steepest descent method or the conjugate
gradient method. Further, may require the use of
preconditioners.

Properties of FEM Matrices
Important Properties of FEM Induced Matrices (A)

1
Matrix A is (generally) a large sparse matrix with a sparsity

factor of less than 1%. The sparsity factor is defined as the
total number of non-zero elements in A divided by the total
number of elements in A expressed as a percentage.
Matrix A is positive definite (that is, x T Ax > 0 x) and hence,
diagonally dominant. The positive definiteness is a direct
consequence of the physical nature of the laws governing an
electromagnetic system.
The condition number of A deteriorates as the size of A
increases. This requires the use of special techniques in the
solution of the large linear system in order to ensure that the
solutions obtained are not spurious.

Basics of FEM: Prof. S. V. Kulkarni

Uploaded by

Copyright:

Available Formats

Basics of FEM: Prof. S. V. Kulkarni

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Basics of FEM: Prof. S. V. Kulkarni

Uploaded by

Copyright:

Available Formats

Indian Institute of Technology Bombay

Indian Institute of Technology Bombay

Properties of FEM Matrices

Indian Institute of Technology Bombay

Finite Element Analysis:

Indian Institute of Technology Bombay

Illustrative Example: Parallel Plate Capacitor

|y =1, 0<x <1 = 10

Indian Institute of Technology Bombay

Minimisation or Maximisation of a Functionalfunction of

Physical Significance of functional: Energy of the system

Indian Institute of Technology Bombay

FEM Problem domain

Size of element: specific to

Indian Institute of Technology Bombay

Element Size Illustration: Arbitrary Boundary

Indian Institute of Technology Bombay

Example Problems Domain Discretisation

Indian Institute of Technology Bombay

Overview of Solution Approximation

Analytical solution may not exist for many problems

Rayleigh-Ritz method (variational procedure)

Indian Institute of Technology Bombay

Overview of Linear System Assembly

Indian Institute of Technology Bombay

2D Linear System Assemly (1 of 2)

Indian Institute of Technology Bombay

= Area of the elemental triangular element =

Ni (x , y ) has the property that,

Indian Institute of Technology Bombay

N1 (x , y )e1 + N2 (x , y )e2 + N3 (x , y )e3

Indian Institute of Technology Bombay

N1 (x , y )e1 + N2 (x , y )e2 + N3 (x , y )e3

N1 (x , y )e1 + N2 (x , y )e2 + N3 (x , y )e3 de

where is the elemental domain

Indian Institute of Technology Bombay

A is referred to as the elemental stiffness matrix. The elemental

Indian Institute of Technology Bombay

where, from equation 8,

Indian Institute of Technology Bombay

Indian Institute of Technology Bombay

Global Connectivity Matrix

Indian Institute of Technology Bombay

Assembly of Global Stiffness Matrix (1 of 3)

The global stiffness matrix is of the form:

A11 A12 A13

A = A31 A32 A33

A91 A92 A93

Indian Institute of Technology Bombay

Assembly of Global Stiffness Matrix (2 of 3)

Indian Institute of Technology Bombay

Assembly of Global Stiffness Matrix (3 of 3)

Thus, the global stiffness matrix is:

Indian Institute of Technology Bombay

Approximate Functional Expression

Thus, equation 17 for the functional can now be expressed as,

where, A is the matrix in equation 29 and is a vector of all nodal

Indian Institute of Technology Bombay