Finite Element Method Application in Electrical Engineering

Abstract 

This paper deals with Finite Element Method (FEM) of analyzing different types of physical systems. Mathematical systems of equations that describe physical systems are often too complex for analytical solving. FEM is a numerical method used for solving these systems. It breaks up the system into the group of elementary elements, called “finite elements”. Then, it analyses each finite element separately. After obtaining information from calculations for each finite element, the results are once brought back together in a single system of equations in order to get data of the system.

In the study of FEM, two systems, a mechanical and its analogue electrical system were analyzed. The FEM is grounded on the principle  of the superposition, and can be used on various types of systems. It has great application in the engineering, while solving complex real – life problems which are usually approximated using numerical method of solving. The FEM method is one of the most precise and straightforward numerical methods.

Authors: Denis Radoncic, and Keisha Baxter

Linear Spring System

MechanicalSystem1

We examined the linear spring system from the figure above. There are six springs in the system. Each spring, along with two nodes between which it lies, is considered to be a finite element, one increment of the system. We inspected forces applied to and distributed on each element and displacements of nodes due to these forces.

Linear spring system is taken to be in a state of equilibrium. It follows that forces applied to nodes at ends of each element must be of equal intensities but opposite directions.

Equilibrium conditions for each element:

    \[ (I)   :   f1 + f2 = 0   => f2 = -f1 \]

    \[ (II)  :   f2 + f3 = 0   => f3 = -f2 \]

    \[ (III) :   f3 + f4 = 0   => f4 = -f3 \]

    \[ (IV)  :   f4 + f5 = 0   => f5 = -f4 \]

    \[ (V)   :   f5 + f6 = 0   => f6 = -f5 \]

    \[ (VI)  :   f2 + f4 = 0   => f4 = -f2 \]

    \[ (VII) :   f2 + f5 = 0   => f5 = -f2 \]

Each force applied to a node causes a displacement of a node. The deformation of the element is determined by displacements of two nodes of the element. Greater the force, greater the deformation (F ~ x). With the proportionality coefficient k (spring constant), the equation is constructed:

    \[ F = kx \]

Deformation is determined by the displacement of nodes, which is taken to be in the same direction for all nodes. Thus:

    \[ x1 = U1 - U2 \]

    \[ x2 = U2 - U3 \]

    \[ x3 = U3 - U4 \]

    \[ x4 = U4 - U5 \]

    \[ x5 = U5 - U6 \]

    \[ x6 = U2 - U4 \]

    \[ x7 = U2 - U5 \]

 

 

Single Element Analysis

Note: fi(j) – component of the force applied to node i that deforms spring j.

element (I) :

    \[ f1(1) = k(U1 - U2)    =>  f2(1) = k(-U1 + U2) \]

Matrix form:

    \[ \begin{bmatrix} k & -k \\ -k & k \end{bmatrix} * \begin{bmatrix} U1 \\ U2 \end{bmatrix} = \begin{bmatrix} f1(1) \\ f2(1) \end{bmatrix} \]

element (II) :

    \[ f2(2) = 2k(U2 - U3)    =>  f3(2) = 2k(-U2 + U3) \]

Matrix form:

    \[ \begin{bmatrix} 2k & -2k \\ -2k &2 k \end{bmatrix} * \begin{bmatrix} U2 \\ U3 \end{bmatrix} = \begin{bmatrix} f2(2) \\ f3(2) \end{bmatrix} \]

element (III) :

    \[ f3(3) = k(U3 - U4)    =>  f4(3) = k(-U3 + U4) \]

Matrix form:

    \[ \begin{bmatrix} k & -k \\ -k & k \end{bmatrix} * \begin{bmatrix} U3 \\ U4 \end{bmatrix} = \begin{bmatrix} f3(3) \\ f4(3) \end{bmatrix} \]

element (IV) :

    \[ f4(4) = 2k(U4 - U5)    =>  f5(4) = 2k(-U4 + U5) \]

Matrix form:

    \[ \begin{bmatrix} 2k & -2k \\ -2k &2 k \end{bmatrix} * \begin{bmatrix} U4 \\ U5 \end{bmatrix} = \begin{bmatrix} f4(4) \\ f5(4) \end{bmatrix} \]

element (V) :

    \[ f5(5) = 2k(U5 - U6)    =>  f6(5) = 2k(-U5 + U6) \]

Matrix form:

    \[ \begin{bmatrix} 2k & -2k \\ -2k &2 k \end{bmatrix} * \begin{bmatrix} U5 \\ U6 \end{bmatrix} = \begin{bmatrix} f5(5) \\ f6(5) \end{bmatrix} \]

element (VI) :

    \[ f2(6) = k(U2 - U4)    =>  f4(6) = k(-U2 + U4) \]

Matrix form:

    \[ \begin{bmatrix} k & -k \\ -k &k \end{bmatrix} * \begin{bmatrix} U2 \\ U4 \end{bmatrix} = \begin{bmatrix} f5(5) \\ f6(5) \end{bmatrix} \]

element (VII) :

    \[ f2(7) = 3k(U2 - U5)    =>  f5(7) = 3k(-U2 + U5) \]

Matrix form:

    \[ \begin{bmatrix} 3k & -3k \\ -3k &3k \end{bmatrix} * \begin{bmatrix} U2 \\ U5 \end{bmatrix} = \begin{bmatrix} f2(7) \\ f5(7) \end{bmatrix} \]

Note:

A 2×2 matrix from the equation is called “stiffness matrix” [Ke]. Stiffnes matrix of a single element is symmetric (Kij = Kji) and singular ([Ke]^-1 does not exist since det[Ke]=0). Size of the stiffness matrix is determined by the number of nodal displacement (degrees of freedom). A linear spring system with N nodal displacements has a stiffness matrix of size NxN. In our system we have six nodal displacements, which delivers a 6×6 corresponding matrix.

From singularity of stiffness matrix:

    \[ \begin{bmatrix} U1 \\U2 \end{bmatrix} = \begin{bmatrix} Ke \end{bmatrix}^-1 = \begin{bmatrix} f1(1) \\ f2(1) \end{bmatrix} \]

As [Ke]^(-1) does not exist, it is impossible to find individual displacements of nodes from their stiffness matrices, but only difference between those displacements (distortion). To find displacements of singular nodes, we need to establish a matrix equation for the whole system.

Incorporating single element results into whole system analysis

Now we need to add up all the systems of equations of single elements to create a corresponding system of equations for entire linear system of springs. Since we have 6 degrees of freedom, this system will fit into the 6×6 matrix. From the illustration and equations of each element of the system we can find out how is each applied force distributed.

As an example, force applied to node to 2 is distributed between elements 1, 2, 6 and 7, and so f2 appears in equations of these elements. Hence, the total force applied to node 2 is is:

    \[ f2 = f2(1) + f2(2) + f2(6) + f2(7) = kU1 + 2kU2 + kU2 + 3kU2 = 7kU2 \]

Examining equations with stiffness matrices we can obtain distribution of each force to each node.

The corresponding system matrix:

    \[ M= \begin{bmatrix} k & -k & 0 & 0 & 0 & 0 \\ -k & 7k & -2k & -k & -3k & 0 \\ 0 & -2k & 3k & -k & 0 & 0 \\ 0 & -k & -k & 4k & -2k & 0 \\ 0 & -3k & 0 & -2k & 7k & 0 \\ 0 & 0 & 0 & 0 & -2k & 2k  \end{bmatrix} \]

Main system equation:

    \[ \begin{bmatrix} k & -k & 0 & 0 & 0 & 0 \\ -k & 7k & -2k & -k & -3k & 0 \\ 0 & -2k & 3k & -k & 0 & 0 \\ 0 & -k & -k & 4k & -2k & 0 \\ 0 & -3k & 0 & -2k & 7k & 0 \\ 0 & 0 & 0 & 0 & -2k & 2k  \end{bmatrix} * \begin{bmatrix} U1 \\ U2 \\ U3 \\ U4 \\ U5 \\ U6 \end{bmatrix} = \begin{bmatrix}  f1 \\ f2 \\ f3 \\ f4 \\ f5 \\ f6   \end{bmatrix} \]

From the figure:

Node 1 is fixed. This means that displacement U1=0.

The force applied to node 6 is 0, since no weight is attached to it.

Incorporating statements from above, the main system equation is reduced to this form:

    \[ \begin{bmatrix} 7k & -2k & -k & -3k & 0 \\ -2k & 3k & -k & 0 & 0 \\ -k & -k & 4k & -2k & 0 \\ -3k & 0 & -2k & 7k & 0 \\ 0 & 0 & 0 & -2k & 2k  \end{bmatrix} * \begin{bmatrix} U2 \\ U3 \\ U4 \\ U5 \\ U6 \end{bmatrix} = \begin{bmatrix}   f2 \\ f3 \\ f4 \\ f5 \\ 0   \end{bmatrix} \]

While we extract:

    \[ f1 = -kU2 \]

For calculation purposes, applied forces’ and spring’s constant values are assumed. Forces applied to nodes 2, 3, 4, 5 are weights of attached masses. (W=Mg) Their values are taken to be 10N. Spring constant is taken to be 100 N/m.

MATLAB code:


f2=10;

f3=10;

f4=10;

f5=10;

f6=0;

k=100;

M=[7*k,-2*k,-k,-3*k,0;-2*k,3*k,-k,0,0;
-k,-k,4*k,-2*k,0;-3*k,0,-2*k,7*k,0;0,0,0,-2*k,2*k];

B=[f2;f3;f4;f5;f6];

U=linsolve(M,B);

U1=0

U2=U(1)

U3=U(2)

U4=U(3)

U5=U(4)

U6=U(5)

f1=-k*U(1)

Note:

Just like with the mechanical system above, the matrix M1 can be obtained in another way too. First we can expand system of equations for each element to the 6×6 size for stiffness matrix, and 6×1 size for column vectors. (k is taken to be 10 N/m)

For example, element (I):

    \[ EL1= \begin{bmatrix} k & -k & 0 & 0 & 0 & 0 \\ -k & k & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0  \end{bmatrix} * \begin{bmatrix} U1 \\ U2 \\ U3 \\ U4 \\ U5 \\ U6 \end{bmatrix} = \begin{bmatrix}  f1 \\ f2 \\ f3 \\ f4 \\ f5 \\ f6   \end{bmatrix} \]

Matlab code:


EL1=[k,-k;-k,k];

D1=zeros(2,4);

G2=zeros(4,6)

EL_fixed=cat(2,EL1,D1);

EL1=[EL_fixed;G2]

Element (II):

    \[ EL2= \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2k & -2k & 0 & 0 & 0 \\ 0 & -2k & 2k & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0  \end{bmatrix} * \begin{bmatrix} U1 \\ U2 \\ U3 \\ U4 \\ U5 \\ U6 \end{bmatrix} = \begin{bmatrix}  f1 \\ f2 \\ f3 \\ f4 \\ f5 \\ f6   \end{bmatrix} \]

MATLAB code:


EL2=zeros(6);

EL2(2,2)=2*k;

EL2(2,3)=-2*k;

EL2(3,2)=2*k;

EL2(3,3)=-2*k;

EL2

Note: Element (II) lies between nodes 2 and 3, which are exposed to forces f2 and f3, and which make displacements of U2 and U3. This is why non-zero elements of its stiffness matrix are placed on (2,2), (2,3), (3,2) and (3,3) positions.

In this way we can find all 6 systems of equations, add them together, and get a final system MU=B like we did in the first method.

Electrical System

ElectricalSystem

The finite element system can be applied to other systems as well.  In the figure below, there is an electrical circuits, an analogue of the mechanical system already studied.

Analogies:

Hook’s law (F=kx) finds it’s analogy in Ohm’s law i = (1/R)V.

Force                        ~  Current                     f  ~  i

Spring stiffness      ~  1/Resistance            k  ~  1/R

Displacement         ~  Voltage                      U  ~  V

Fixed node              ~  Grounded node

As in Spring system analysis, we can make a system of 2 equations for each element.

Note: In ij(k), j represents a node, and k represents an element to which a component of current i is going.

element(I) :

    \[ \frac{1}{R}(V1-V2) = i1(1)     =>   i2(1) = 1/R(-V1+V2) \]

Matrix form:

    \[ \begin{bmatrix} 1/R & -1/R \\ -1/R & 1/R \end{bmatrix} * \begin{bmatrix} V1 \\ V2 \end{bmatrix} = \begin{bmatrix} i1(1) \\ i2(2) \end{bmatrix} \]

element (II) :

    \[ i2(2) = 2/R(V2 - V3)    =>  i3(2) = 2/R(-V2 + V3) \]

Matrix form:

    \[ \begin{bmatrix} 2/R & -2/R \\ -2/R & 2/R \end{bmatrix} * \begin{bmatrix} V2 \\ V3 \end{bmatrix} = \begin{bmatrix} i2(2) \\ i3(2) \end{bmatrix} \]

element (III) :

    \[ i3(3) = 1/R(V3 - V4)    =>  i4(3) = 1/R(-V3 + V4) \]

Matrix form:

    \[ \begin{bmatrix} 1/R & -1/R \\ -1/R & 1/R \end{bmatrix} * \begin{bmatrix} V3 \\ V4 \end{bmatrix} = \begin{bmatrix} i3(3) \\ i4(3) \end{bmatrix} \]

element (IV) :

    \[ i4(4) = R/2(V4 - V5)    =>  i5(4) = R/2(-V4 + V5) \]

Matrix form:

    \[ \begin{bmatrix} 2/R & -2/R \\ -2/R & 2/R \end{bmatrix} * \begin{bmatrix} V4\\ V5 \end{bmatrix} = \begin{bmatrix} i4(4) \\ i5(4) \end{bmatrix} \]

element (V) :

    \[ i5(5) = 2/R(V5 - V6)    =>  i6(5) = 2/R(-V5 + V6) \]

Matrix form:

    \[ \begin{bmatrix} 2/R & -2/R \\ -2/R & 2/R \end{bmatrix} * \begin{bmatrix} V5 \\ V6 \end{bmatrix} = \begin{bmatrix} i5(5) \\ i6(5) \end{bmatrix} \]

element (VI) :

    \[ i2(6) = 1/R(V2 - V4)    =>  i4(6) = 1/R(-V2 + V4) \]

Matrix form:

    \[ \begin{bmatrix} 1/R & -1/R \\ -1/R & 1/R \end{bmatrix} * \begin{bmatrix} V2 \\ V4 \end{bmatrix} = \begin{bmatrix} i2(6) \\ i4(6) \end{bmatrix} \]

element (VII) :

    \[ i2(7) = 3/R(V2 - V5)    =>  i5(7) = 3/R(-V2 + V5) \]

Matrix form:

    \[ \begin{bmatrix} 3/R & -3/R \\ -3/R & 3/R \end{bmatrix} * \begin{bmatrix} V2 \\ V5 \end{bmatrix} = \begin{bmatrix} i2(7) \\ i5(7) \end{bmatrix} \]

Using the same technique as with Mechanical system, we get a corresponding system matrix.  Since there is no current generator attached to node 6, i6=0. The corresponding system matrix:

    \[ M1= \begin{bmatrix} 1/R & -1/R & 0 & 0 & 0 & 0 \\ -1/R & 7/R & -2/R & -1/R & -3/R & 0 \\ 0 & -2/R & 3/R & -1/R & 0 & 0 \\ 0 & -1/R & -1/R & 4/R & -2/R & 0 \\ 0 & -3/R & 0 & -2/R & 7/R & 0 \\ 0 & 0 & 0 & 0 & -2/R & 2/R  \end{bmatrix} \]

Main system equation:

    \[ \begin{bmatrix} 1/R & -1/R & 0 & 0 & 0 & 0 \\ -1/R & 7/R & -2/R & -1/R & -3/R & 0 \\ 0 & -2/R & 3/R & -1/R & 0 & 0 \\ 0 & -1/R & -1/R & 4/R & -2/R & 0 \\ 0 & -3/R & 0 & -2/R & 7/R & 0 \\ 0 & 0 & 0 & 0 & -2/R & 2/R  \end{bmatrix} * \begin{bmatrix} V1 \\ V2 \\ V3 \\ V4 \\ V5 \\ V6 \end{bmatrix} = \begin{bmatrix}  i1 \\ i2 \\ i3 \\ i4 \\ i5 \\ 0   \end{bmatrix} \]

Just like force f2 in mechanical analogous system, a current applied to the second node i2 is distributed between elements 1,2,6, and 7. So, in the main system equation, in the (2,2) position (which belongs to a row that corresponds to i2), we get 7/R in this way:

    \[ i2 = i2(1)+i2(2)+i2(6)+i2(7) \]

    \[ i2 = (1/R)V2 + (2/R)V2 + (1/R) V2 + (3/R)V2 \]

    \[ i2 = (7/R)V2 \]

In this way we fill out the whole matrix. As the node one is grounded, its voltage will equal to zero. This gives:

    \[ i1=(-1/R)V2 \]

When we abstract this equation out of the system, we can get smaller system of equations:

    \[ \begin{bmatrix} 7/R & -2/R & -1/R & -3/R & 0 \\  -2/R & 3/R & -1/R & 0 & 0 \\  -1/R & -1/R & 4/R & -2/R & 0 \\  -3/R & 0 & -2/R & 7/R & 0 \\ 0 & 0 & 0 & -2/R & 2/R  \end{bmatrix} * \begin{bmatrix} V1 \\ V2 \\ V3 \\ V4 \\ V5 \\ V6 \end{bmatrix} = \begin{bmatrix}  i1 \\ i2 \\ i3 \\ i4 \\ i5 \\ 0   \end{bmatrix} \]

For the calculation purposes, values of R and all currents are assumed. Value of R is assumed to be 10 ohms, and those of all applied currents 1A.

MATLAB code:

i2=1;

i3=1;

i4=1;

i5=1;

i6=0;

R=10;

M1=[7/R,-2/R,-1/R,-3/R,0;-2/R,3/R,-1/R,0,0;
-1/R,-1/R,4/R,-2/R,0;-3/R,0,-2/R,7/R,0;0,0,0,-2/R,2/R];

B1=[i2;i3;i4;i5;i6];

V=linsolve(M1,B1);

V1=0

V2=V(1)

V3=V(2)

V4=V(3)

V5=V(4)

V6=V(5)

i1=(-1/R)*V(1)

References:

Kwon, Y., & Bang, H. (2000). The finite element method using MATLAB (2nd ed.).
Boca Raton, FL: CRC Press.