3.4 FE Formulation of Electro-Magnetic Field Problems

Chapter 3
through a surface having unitary length in z-direction with the

sides of the coil parallel to z-axis. (Bargallo, 2006). In a coil with
sides 1 and 2, the ux per unit length is given by:
= A1 A2
The ux linked with the coil is given by:
1
=
S

A dS
As per Faradays law, induced e.m.f. in the conductor can be

obtained by numerical derivation of the ux linkage. Thus, from
the value of the primary variable obtained from the eld simulation, all secondary variables are derived and hence machine
parameters like induced voltage, inductances, d and q-axes reactances, saturation factor are computed (Bianchi, 2001) based on
the fundamental electro-magnetic principles.
3.4
FE Formulation of Electro-magnetic Field

Problems
Design optimisation of electrical machinery is a challenging

job for the designer because it demands competitive price, high
power to weight ratio and high degree of reliability in operation.
This necessitates virtual prototyping and parameter estimation
of the machine under development. But, complex geometry together with presence of non-linear magnetic materials pose diculty in analysis and prediction of magnetic eld pattern in such
machines.
School of Engineering, CUSAT
59
Chapter 3
In order to predict the performance of electrical machines accurately at the design stage, magnetic eld mapping of the machine structure is necessary which requires accurate modelling
of the machine. Conventional methods are based on approximations and hence inaccurate. FEM oers an ecient and accurate
means of predicting the magnetic eld distribution in machines,
taking into account non-linearity of the iron material.
3.4.1
Formulation of Brushless Alternator Field Problem
3D analysis of the synchronous generator involves dividing

the whole structure by 3D nite elements which require heavy
processing and results in long computation time. The magnetic
eld existing in the machine structure can be approximated to be
2D since there exists a symmetry due to repetition of the physical
phenomena, plane by plane. This happens as the currents ow
in a direction parallel to the z-axis and the transverse sections to
this axis always present same geometry and the same materials,
point by point (Bargello, 2006) Considering planar symmetry of
the eld problem as shown in Figure 3.1, the eld properties of
the machine are invariant along the z-axis and the dimensionality of the problem can be reduced to two, thereby considerably reducing the computational time (Bianchi, 2005; Bargello,
2006; Pyrhonen et al., 2009). Hence, the brushless alternator eld
problem under the present study is reduced to 2D eld analysis
without compromising the accuracy of results.
Study of the 2D eld problem with planar symmetry has been
detailed as follows:
1. The current density vector has z-axis component only.
i.e.,
= [0, 0, z ]
(3.12)
60
Chapter 3
Fig. 3.1: 2D eld and boundary conditions for a salient pole synchronous machine(Source: Pyrhonen et al. (2009))
2. The magnetic vector potential A is parallel to the vector ,

thus it has z-axis component only. i.e.,
A = [0, 0, Az ]
(3.13)
As regards the divergence of the magnetic vector potential,

Coulombs gauge is adopted, i.e., divA = A = 0
3. The ux density vector B has components only in the (x,y)
plane and is derived from A using the formula:
B = curlA = A = [0, 0, Az ]

Az Az
=
,
,0
y
x
(3.14)
4. With a constant magnetic permeability, the magnetic eld

distribution along the cross section of a synchronous machine is described by the two dimensional Poissons as given
61
Chapter 3
in Equation 3.10. MVP can be determined from Poissons
equation, at each point in the problem domain if the current density is xed and the value of Az is known on the
boundary of the domain, i.e., in a part of the domain,
say 1 , the value of Az is assigned, and in the remaining
part of the boundary, say 2 , the value of the derivative of
Az normal to the boundary line is assigned.
5. Since the ux lines are tangential to the outer boundary
of the machine structure, no ux lines cross the boundary.
Hence homogeneous Dirichlets condition (Bianchi, 2005) is
assigned, xing the MVP, Az = 0 along the outer boundary.
6. The magnetic behaviour in the machine structure has periodic symmetry among pole pairs around the z-axis. The
machine under study being 4-pole, the domain under analysis is reduced to one fourth on account of the periodic
symmetry of the eld pattern, which further reduces the
computational time for eld analysis. Therefore, along the
axes of periodic symmetry, the periodic boundary conditions are applied. Since the ux density components cannot
be discontinuous across the symmetry axis, the ux lines
must be normal to the axis. Hence the Neumanns boundary condition, i.e., A
= 0, is applied along the symmetry
n
axis (Bianchi, 2005).
3.4.2
Formulation by Galerkins Method
Let A be the approximate solution of the magnetic vector

potential. Substituting the approximate solution A in Equation
(3.10), the residue R can be made available as (Silvester et al.,
1996; Bianchi, 2005; Virjoghe et al., 2012):
R=
62
1 2 A 1 2 A
+
+
x2
y 2
(3.15)
Chapter 3
The weighted integral form for 2D solution is:

R W dx dy = 0
(3.16)
where W is the weight function (Finlayson and Scriven, 1966;

Finlayson, 2013). Substituting the value of R in the weighted
integral form,

1 2 A 1 2 A
dx dy =
W
+
W dx dy (3.17)
x2
y 2
Integrating L.H.S. of Equation (3.17) by parts,

1 W A W A
1 A
+
W
dc
dx dy
x x
y y
n

c

=
W dx dy (3.18)
where n
is the outward unit vector on boundary c.
Discretisation of the problem domain in Equation (3.18) is performed using 2D nite elements. Six noded triangular elements
are used for discretisation of the problem domain which has three
nodes at the vertices of the triangle, three mid-side nodes and the
variable interpolation within the element is quadratic in x and y
(Bastos and Sadowski, 2003). For M triangular elements,

1
W e Ae W e Ae
+
dx dy
e
x x
y y
e
M

1 Ae
e
W dc =
W e dx dy
(3.19)
e
n
c
e
Geometry of the machine discretised with second order triangular
elements are shown in Figure 3.2.
63
Chapter 3
Fig. 3.2: Geometry of the machine discretised with triangular elements
The approximate solution of A at any point (x,y) in the element

is represented in terms of the nodal values of A and the shape
functions,
A1
A2
A3

e
A(x,y) = N1 N2 N3 N4 N5 N6
(3.20)
A
4
A5
A6
Three digit numbering is given for each node of the second order
triangular element, where , and are integers satisfying the
relation, + + = n, and n is the order of the interpolation
function for the element (Chari and Silvester, 1980; Hoole, 1989;
Humphries, 1997; Salon et al., 1999). Triple-subscript designation is used for the interpolation functions, so that (1 , 2 , 3 )
will denote the shape function for node P (, , ) in terms of the
area coordinates 1 , 2 , 3 . As per Silvesters method (Chari and
Salon, 2000), the interpolation function for an nth order triangular element is derived by the formula:
(1 , 2 , 3 ) = (1 ) (2 ) (3 )
64
(3.21)
Chapter 3
where
(1 ) =

n1 i + 1
i=1
, 1
(3.22)
The interpolation functions for the six noded second order triangular element, with n = 2 and three-digit notation for each of
the nodes, 200 , 020 , 002 for the three vertices and 101 , 110 , 011
for the three mid-side nodes respectively are evaluated by the
formula:
For node 1,
= 2, = 0, = 0
= 2 = 1 (21 1)
= 0 = 1
= 0 = 1
i.e.,
N1 = 200 = 2 (1 ) 0 (2 ) 0 (3 ) = 1 (21 1)
(3.23)
Shape functions for the remaining nodes are derived as:

N2 = 020 = 2 (22 1)
N3 = 002 = 3 (23 1)
N4 = 101 = 43 1
N5 = 110 = 41 2
N6 = 011 = 42 3
As per the iso-parametric interpolation formula,
N1

N2
1
1 1 1 1 1 1
x x 1 x 2 x 3 x 4 x 5 x 6 N3

y = y 1 y 2 y 3 y 4 y 5 y 6 N4
A1 A2 A3 A4 A5 A6 N5
A
N6
(3.24)
(3.25)
(3.26)
(3.27)
(3.28)
(3.29)
65
Chapter 3
in which the natural derivatives are:
41 1
N
0
=
4
1
2
0
43
0
42 1
N
0
=
2 41
43
0
43 1
N
=
0
3
42
41
(3.30)
(3.31)
(3.32)
Then, the magnetic vector potential A is expressed as:

A =
6

Ni (x, y) Ai
(3.33)
i=1
In Galerkins weighted residual method, weight function is the

shape function itself.i.e.,
N1
200
N2 020
N3 002
e

W =
(3.34)
N4 = 101
N5 110
N6
011
66
Chapter 3
The shape function logic is concerned with computation of the
partial derivatives of the shape functions with respect to x and
y at any point in the element. For this purpose, a generic scalar
function, w (1 , 2 , 3 ) is considered that is quadratically interpolated over the triangle by,
w = w 1 N 1 + w 2 N 2 + w 3 N3 + w 4 N 4 + w 5 N 5 + w 6 N 6

= w1 w2 w3 w4 w5 w6 N T
(3.35)
where w stands for 1, x, y or A . Taking partials of Equation
(3.35) with respect to x and y and applying the chain rule twice
yields,

Ni
x
Ni
y

=
Ni
1
Ni
2
Ni
3

P
(3.36)
In matrix form,

w
x
w
y
1
x
1
y
2
x
2
y
3
x
3
y
6
i=1
i
wi N
1
6
Ni
w
i
i=1
2
6
Ni
i=1 wi 3
(3.37)
Transposing both sides of Equation (3.37) while exchanging sides

yields,

6
Ni
i=1 i 1
6
Ni
i=1 i 2
6
Ni
i=1 i 3
1
x
2
x
3
x
1
y
2
y
3
y
(3.38)
Making w = 1, x, y and stacking the results row-wise,
6
Ni
i=1 1
6
Ni
i=1 2
6
Ni
i=1 3
6
6
6
Ni
Ni
Ni
x
x
x
i
i
i
i=1
i=1
i=1
1
2
3
6
6
6
Ni
Ni
Ni
y
y
y
i=1 i 1
i=1 i 2
i=1 i 3
1
x
2
x
3
x
1
y
2
y
3
y
1
x
x
x
y
x
1
y
x
y
y
y
(3.39)
67
Chapter 3
1
1
But, x
= y
= 1 and x
= y
= x
= y = 0 because x and y
x
y
y
6 x
are independent coordinates. Since i=1 Ni = 1, the entries of
the rst row of the coecient matrix are equal to a constant C,
which can be scaled to unity because the rst row of the right
hand side is null. Therefore, Equation (3.39) gets reduced to,
1
1
1
6
6

Ni
Ni
6 xi Ni
i=1 xi 2
i=1 xi 3
i=1 1
6
6
6
Ni
Ni
Ni
y
y
y
i=1 i 1
i=1 i 2
i=1 i 3
1
x
2
x
3
x
1
y
2
y
3
y
0 0
= 1 0
0 1
(3.40)
The coecient matrix of equation (3.40) is called the Jacobian
matrix denoted by J. For compactness, Equation (3.40) is written
as
1 1
x
y
1
1
1
0 0
2 2

JP = Jx1 Jx2 Jx3

(3.41)
x y = 1 0
Jy1 Jy2 Jy3
0 1
3
3
x
Solving this system gives,

1 1
x
2
x
3
x
Jy23 Jx32
= 1 Jy31 Jx13 = P
2J
Jy12 Jx21
y
2
y
3
y
(3.42)
in which Jxji = Jxj Jxi , Jyji = Jyj Jyi and

J = 12 detJ = 12 (Jx21 Jy31 Jy12 Jx13 )
Substituting into Equation (3.37):
w Ni wi
=
=
wi
x
x
2J
i=1
i=1
6
68

Ni
Ni
Ni
Jy23 +
Jy31 +
Jy12
1
2
3
(3.43)
Chapter 3

Ni
Ni
Ni
Jx32 +
Jx13 +
Jx21
1
2
3
(3.44)
The shape function derivatives are,

Ni
1 Ni
Ni
Ni
(3.45)
=
Jy23 +
Jy31 +
Jy12
x
2J 1
2
3

Ni
1 Ni
Ni
Ni
(3.46)
=
Jx32 +
Jx13 +
Jx21
y
2J 1
2
3
w Ni wi
=
=
wi
y
y
2J
i=1
i=1
6
i
can be read o Equations
in which the natural derivatives N
j
(3.30) to (3.32). Use of the P matrix dened in Equation (3.42),
yields Equations (3.45) and (3.46) in compact form as,

Ni
Ni
Ni
Ni
Ni
= 1 2 3 P
(3.47)
x
y
These dierentials are substituted in the weighted integral form

of Equation (3.18) and the second term in the equation,
1 A
W n dc = 0 as the component of A along the outward unit
c
vector n
on boundary c is zero. Also,

W dx dy = 3
0
0
0
1
1
1
0
0
0
1
1
1
(3.48)
Product of the current density, and area of the triangle, gives

the total current owing through the triangle. Hence, the governing partial dierential equation of the eld form represented
by Equation (3.18) gets reduced to the matrix form,
[S][A] = [I]
(3.49)
69

3.4 FE Formulation of Electro-Magnetic Field Problems

Uploaded by

Copyright:

Available Formats

3.4 FE Formulation of Electro-Magnetic Field Problems

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

3.4 FE Formulation of Electro-Magnetic Field Problems

Uploaded by

Copyright:

Available Formats

Chapter 3

through a surface having unitary length in z-direction with the

As per Faradays law, induced e.m.f. in the conductor can be

FE Formulation of Electro-magnetic Field

Design optimisation of electrical machinery is a challenging

Formulation of Brushless Alternator Field Problem

3D analysis of the synchronous generator involves dividing

School of Engineering, CUSAT

2. The magnetic vector potential A is parallel to the vector ,

As regards the divergence of the magnetic vector potential,

4. With a constant magnetic permeability, the magnetic eld

Formulation by Galerkins Method

Let A be the approximate solution of the magnetic vector

School of Engineering, CUSAT

where W is the weight function (Finlayson and Scriven, 1966;

Integrating L.H.S. of Equation (3.17) by parts,

Fig. 3.2: Geometry of the machine discretised with triangular elements

The approximate solution of A at any point (x,y) in the element

School of Engineering, CUSAT

Shape functions for the remaining nodes are derived as:

Then, the magnetic vector potential A is expressed as:

In Galerkins weighted residual method, weight function is the

School of Engineering, CUSAT

Transposing both sides of Equation (3.37) while exchanging sides

School of Engineering, CUSAT

JP = Jx1 Jx2 Jx3

Solving this system gives,

in which Jxji = Jxj Jxi , Jyji = Jyj Jyi and

These dierentials are substituted in the weighted integral form

Product of the current density, and area of the triangle, gives

You might also like