Still having doubt’s?
Awab Sir-89 76 104646
MODULE 04 COMPUTATIONAL ELECTROMAGNETICS Finite Difference Method (FDM): Neumann type and mixed boundary conditions, Iterative solution of finite difference equations, solutions using band matrix method Finite Element Method (FEM): Triangular mesh configuration, Finite element discretization, Element governing equations, assembling all equations and solving resulting equations Method of Moment (MOM): Field calculations of conducting wire, parallel conducting wires and complicated geometries Text Book Bhag Guru and Huseyin Hiziroglu, “Electromagnetic field theory fundamentals”, Cambridge University Press, 2nd Edition, 2010. Reference Book Matthew N.D. SADIKU, Elements of Electromagnetics, Oxford International Student 4th Edition, 2007 EXPECTED QUESTIONS There will be questions on the basic methods used and method of applying boundary condition. Procedure (steps) for solving the method should be known. No numerical where iterations are required would be asked in the exam. INTRODUCTION When the complexities of theoretical formulas make analytic solution intractable, we resort to non analytic methods, which include 1. Graphical methods 2. Experimental methods 3. Analog methods 4. Numerical methods Graphical, experimental, and analog methods are applicable to solving relatively few problems. Numerical methods have come into prominence and become more attractive with the advent of fast digital computers. The three most commonly used simple numerical techniques in EM are 1. Moment method 2. Finite difference method 3. Finite element method. Most EM problems involve either partial differential equations or integral equations. Partial differential equations are usually solved using the finite difference method or the finite element method; integral equations are solved conveniently using the moment method. Although numerical methods give approximate solutions, the solutions are sufficiently accurate for engineering purposes. www.awabsir.com
Page 1
Still having doubt’s?
Awab Sir-89 76 104646
Consider the Laplace equation which is given as follows: And a source free equation given as
Where, u is the electrostatic potential. Now if we take an example of a parallel plate capacitor as shown in the figure below. When we neglect the fringing field we get the following equation. Where, V is the closed form analytical solution (can see the effects of varying any quantity on the RHS by significant change in the LHS). For solving the equation 3 we do not need the use of numerical techniques.
However if we now take into account the fringing field the solution for every point x & y such that equation 4 is satisfied is not manually possible
To find the field intensity at any point we need to use numerical techniques. www.awabsir.com
Page 2
Still having doubt’s?
Awab Sir-89 76 104646
Step 1: Divide the given problem domain into sub-domains It is a tough job to approximate the potential for the entire domain at a glance. Therefore any domain in which the field is to be calculated is divided into small elements. We use sub domain approximation instead of whole domain approximation. Considering a one dimensional function
Step 2: Approximate the potential for each element The approximate potential for an element can be given as OR Where, a, b and c are constants. Considering the first equation for each element, the potential distribution will be approximated as a straight line for every element. Step 3: Find the potential u for every element in terms of end point potentials Assuming The above equation can also be written as We can write
We can write equation 1 and 2 in matrix form as follows Rearranging the terms we get
Substituting equation 5 in equation 1 we get
www.awabsir.com
Page 3
Still having doubt’s?
Awab Sir-89 76 104646
Similarly we can find the electrostatic potential for each element.
Step 4: Find the energy for every element The energy for a capacitor is given as follows:
The field is distributed such that the energy is minimized. The electric field intensity is related to the electrostatic potential as follows
Step 5: Find the total energy The total energy is the summation of the individual electrostatic energy of every element in the domain. Step 6: Obtain the general solution The field within the domain is distributed such that the energy is minimized. For minimum energy the differentiation of electrostatic energy with respect to the electrostatic potential is equated to zero for every element.
Solving this we get a matrix The matrix [K] is a function of geometry and the material properties. The curly brackets denote column matrix. Equation 9 does not lead to a unique solution. For a unique solution we have to apply boundary conditions. Step 7: Obtain unique solution We need to apply boundary conditions to equation 9 to obtain a unique solution. For example we assume u1=1 V and u4=5 V. On applying boundary conditions the RHS becomes a non zero matrix and a unique solution can be obtained.
THE FINITE ELEMENT METHOD
www.awabsir.com
Page 4
Still having doubt’s?
Awab Sir-89 76 104646
The finite element analysis of any problem involves basically four steps: 1. Discretizing the solution region into a finite number of sub regions or elements 2. Deriving governing equations for a typical element 3. Assembling of all elements in the solution region 4. Solving the system of equations obtained. A. Finite Element Discretization
We divide the solution region into a number of finite elements as illustrated in the figure below, where the region is subdivided into four non overlapping elements (two triangular and two quadrilateral) and seven nodes. We seek an approximation for the potential Ve within an element e and then inter-relate the potential distributions in various elements such that the potential is continuous across inter-element boundaries. The approximate solution for the whole region is
Where, N is the number of triangular elements into which the solution region is divided. The most common form of approximation for Ve within an element is polynomial approximation, namely
for a triangular element and
for a quadrilateral element. The potential Ve in general is nonzero within element e but zero outside e. It is difficult to approximate the boundary of the solution region with quadrilateral elements; such elements are useful for problems whose boundaries are sufficiently regular. In view of this, we prefer to use triangular elements throughout our analysis in this section. Notice that our assumption of linear variation of potential within the triangular element as in eq. (2) is the same as assuming that the electric field is uniform within the element; that is, www.awabsir.com
Page 5
Still having doubt’s?
Awab Sir-89 76 104646
B. Element Governing Equations
The potential Ve1, Ve2, and Ve3 at nodes 1, 2, and 3, respectively, are obtained using eq. (2); that is,
We can obtain the values of a, b and c. Substituting these values in equation 2 we get
OR
Where
www.awabsir.com
Page 6
Still having doubt’s?
Awab Sir-89 76 104646
And A is the area of the element e
The value of A is positive if the nodes are numbered counterclockwise. Note that eq. (5) gives the potential at any point (x, y) within the element provided that the potentials at the vertices are known. This is unlike the situation in finite difference analysis where the potential is known at the grid points only. Also note that α, are linear interpolation functions. They are called the element shape functions. The shape functions α1 and α2 for example, are illustrated in the figure below
The energy per unit length can be given as
Where
And www.awabsir.com
Page 7
Still having doubt’s?
Awab Sir-89 76 104646
The matrix C(e) is usually called the element coefficient matrix. The matrix element Cij(e) of the coefficient matrix may be regarded as the coupling between nodes i and j.
C. Assembling of all Elements Having considered a typical element, the next step is to assemble all such elements in the solution region. The energy associated with the assemblage of all elements in the mesh is
Where n is the number of nodes, N is the number of elements, and [C] is called the overall or global coefficient matrix, which is the assemblage of individual element coefficient matrices. The properties of matrix [C] are 1. It is symmetric (Cij = Cji) just as the element coefficient matrix. 2. Since Cij = 0 if no coupling exists between nodes i and j, it is evident that for a large number of elements [C] becomes sparse and banded. 3. It is singular. D. Solving the Resulting Equations From variational calculus, it is known that Laplace's (or Poisson's) equation is satisfied when the total energy in the solution region is minimum. Thus we require that the partial derivatives of W with respect to each nodal value of the potential be zero; that is,
To find the solution we can use either the iteration method or the band matrix method. The finite element method (FEM) has a number of advantages over the finite difference method (FDM) and the method of moments (MOM). First, the FEM can easily handle complex solution www.awabsir.com
Page 8
Still having doubt’s?
Awab Sir-89 76 104646
region. Second, the generality of FEM makes it possible to construct a general-purpose program for solving a wide range of problems. A single program can be used to solve different problems (described by the same partial differential equations) with different solution regions and different boundary conditions; only the input data to the problem need be changed. However, FEM has its own drawbacks. It is harder to understand and program than FDM and MOM. It also requires preparing input data, a process that could be tedious.
NOTE: Refer the example solved in class.
THE FINITE DIFFERENCE METHOD Boundary Conditions A unique solution can be obtained only with a specified set of boundary conditions. There are basically three kinds of boundary conditions: 1. Dirichlet type of boundary 2. Neumann type of boundary 3. Mixed boundary conditions
1. Dirichlet Boundary Condition Consider a region s bounded by a curve l. If we want to determine the potential distribution V in region s such that the potential along l is V=g. Where, g is prespecified continuous potential function. Then the condition along the boundary l is known as Dirichlet Boundary condition. www.awabsir.com
Page 9
Still having doubt’s?
Awab Sir-89 76 104646
2. Neumann Boundary Condition Neumann boundary condition is mathematically represented as
Where, the conditions along the boundary are such that the normal derivative of the potential function at the boundary is specified as a continuous function. 3. Mixed Boundary Condition There are problems having the Dirichlet condition and Neumann condition along l1 and l2 portions of l respectively. This is defined as mixed boundary condition.
A problem is uniquely defined by three things: 1. A partial differential equation such as Laplace's or Poisson's equations. 2. A solution region. 3. Boundary and/or initial conditions. A finite difference solution to Poisson's or Laplace's equation, for example, proceeds in three steps: 1. Dividing the solution region into a grid of nodes. 2. Approximating the differential equation and boundary conditions by a set of linear algebraic equations (called difference equations) on grid points within the solution region. 3. Solving this set of algebraic equations. Step 1: Suppose we intend to apply the finite difference method to determine the electric potential in a region, shown in the figure below. The solution region is divided into rectangular meshes with grid points or nodes as shown. A node on the boundary of the region where the potential is specified is called a fixed node (fixed by the problem) and interior points in the region are called free points (free in that the potential is unknown).
www.awabsir.com
Page 10
Still having doubt’s?
Awab Sir-89 76 104646
Step 2: Our objective is to obtain the finite difference approximation to Poisson's equation and use this to determine the potentials at all the free points. Step 3: To apply the following equation, to a given problem, one of the following two methods is commonly used. ) ……. (1)
A. Iteration Method We start by setting initial values of the potentials at the free nodes equal to zero or to any reasonable guessed value. Keeping the potentials at the fixed nodes unchanged at all times, we apply eq. (1) to every free node in turn until the potentials at all free nodes are calculated. The potentials obtained at the end of this first iteration are not accurate but just approximate. To increase the accuracy of the potentials, we repeat the calculation at every free node using old values to determine new ones. The iterative or repeated modification of the potential at each free node is continued until a prescribed degree of accuracy is achieved or until the old and the new values at each node are satisfactorily close. B. Band Matrix Method Equation (1) applied to all free nodes results in a set of simultaneous equations of the form
Where: [A] is a sparse matrix (i.e., one having many zero terms), [V] consists of the unknown potentials at the free nodes, and [B] is another column matrix formed by the known potentials at the fixed nodes. www.awabsir.com
Page 11
Still having doubt’s?
Awab Sir-89 76 104646
Matrix [A] is also banded in that its nonzero terms appear clustered near the main diagonal because only nearest neighboring nodes affect the potential at each node. The sparse, band matrix is easily inverted to determine [V]. Thus we obtain the potentials at the free nodes from matrix [V] as
The concept of FDM can be extended to Poisson's, Laplace's, or wave equations in other coordinate systems. The accuracy of the method depends on the fineness of the grid and the amount of time spent in refining the potentials. We can reduce computer time and increase the accuracy and convergence rate by the method of successive over relaxation, by making reasonable guesses at initial values, by taking advantage of symmetry if possible, by making the mesh size as small as possible, and by using more complex finite difference molecules. One limitation of the finite difference method is that interpolation of some kind must be used to determine solutions at points not on the grid. One obvious way to overcome this is to use a finer grid, but this would require a greater number of computations and a larger amount of computer storage.
Example 1: Determine the electrostatic potential distribution inside the region shown in the figure below, with boundary conditions as specified. y
v = 100 V 3
v=0
v=0
0
v=0
3
x
The boundaries range between
www.awabsir.com
Page 12
Still having doubt’s? At
x=0 y=0 x=0
Awab Sir-89 76 104646
0
zero potentials
i.e. the potentials along these boundaries is constant satisfies Dirichlet‟s condition. At, y = 3
0
the potential is 100v which is also a constant and thus, this boundary also satisfies Dirichlet‟s condition. Note: The small separation „δ‟ are desirable as the upper horizontal boundary is maintained at a different potential than the others. To use FDM, we divide the region into square meshes with h = 1 y
v = 100 (1, 3)
(2, 3)
(0, 3)
(3, 3) (1, 2) 1
(2, 2) 2
(0, 2)
(3, 2)
v=0
v=0 (0, 1)
(1, 1)
(2, 1)
3
(3, 1) 4 x
(0, 0)
(1, 0)
(2, 0)
(3, 0)
v=0
we have to determine the potential at nodes (1, 1), (2, 1), (1, 2) and (2, 2) the potential at nodes (1, 3) and (2, 3) is given as 100 and the potential at all other nodes is zero From the general equation V1 + V2 + V3 + V4 – 4V0 = 0 V0 = ¼ (V1 + V2 + V3 + V4)
V1
we get, V (1, 1) = V3 = ¼ (0 + 0 + V1 + V4) V (2, 1) = V4 = ¼ (0 + 0 + V3 + V2) V (2, 2) = V2 = ¼ (0 + 100 + V1 + V4) V (1, 2) = V1 = ¼ (0 + 100 + V2 + V3) www.awabsir.com
V0 V2
V4
Page 13
Still having doubt’s?
Awab Sir-89 76 104646
by rearranging the equation we get, V3
4V1 – V2 – V3 = 100 - V1 + 4V2 – V4 = 100 - V1 + 4V3 – V4 = 0 - V2 – V3 + 4V4 = 0
V4 =
In matrix form, =
This is standard form for system of linear equation AV = b where, A is square matrix V = A-1b V4 = (V2 + V3) 0.25 substitute in other three equations 4V1 – V2 – V3 = 100 ------------------------- (1) - V1 + 4V2 – 0.25(V2 + V3) = 100 - V1 + 3.75V2 – 0.25V3 = 100 -------------- (2) - V1 + 4V3 – V4 = 0 - V1 + 4V3 – 0.25(V2 + V3) = 0 - V1 – 0.25V2 + 3.75V3 = 0 ----------------- (3) Solving the simultaneous equation (1), (2) & (3) we get V1 = 37.5, V2 = 37.5, V3 = 12.5, and V4 = 12.5. Example 2: Solve the node potentials in the geometry using SOR method. y v = 100 3
v=0
v=0
www.awabsir.com
Page 14 x 0
v=0
3
Still having doubt’s?
Awab Sir-89 76 104646
Solution, Dividing into mesh, 100
100
100
1
100
2
0
0
0
0 3
4
0
0 0
0
Let us consider an initial guess of 50V at nodes and an error criteria of 0.1 and acceleration factor of 1 Iteration 1 =
+
= 50 + |
– – –
–
+
–4
)
(100 + 0 + 50 + 50 - 200) = 50 (100 + 0 + 50 + 50 - 200) = 50 (0 + 0 + 50 + 50 - 200) = 25
| = |50 – 25| = 25
= 50 + |
+
| = |50 – 50| = 0
= 50 + |
+
| = |50 – 50| = 0
= 50 + |
(
(50 + 25 + 0 + 0 - 200) = 18.75
| = |50 – 18.75| = 31.25
Iteration 2 = 50 + 0.25 (100 + 0 + 25 + 50 - 200) = 43.75 | – | = |50 – 43.75| = 6.25 = 50 + 0.25 (100 + 43.75 + 18.75 + 0 - 200) = 40.625 | – | = 9.375 = 25 + 0.25 (43.75 + 0 + 0 + 18.75 - 100) = 15.625 | – | = 9.375 = 18.75 + 0.25 (40.625 + 15.625 + 0 + 0 - 75) = 14.06 | – | = 4.68
100
www.awabsir.com
100
100
100
50 40.625 38.28 37.69
Page 15 0
Still having doubt’s?
1 0
50 43.75 39.06 37.89
Awab Sir-89 76 104646
2
37.59 37.53 25 15.625 13.28 0
3
0
12.69 12.55 12.52
0
37.55 37.52 18.75 14.06 12.89 4
0 12.59 12.53 12.51
0
0
Iteration 3 | |
= 43.75 + 0.25 (100 + 0 + 15.625 + 40.625 – 4 | = 4.6875 = 38.28 = 13.28 | = 2.345 | | = 2.345
Iteration 4 = 37.89 | | = 1.17 Iteration 5 = 37.59 | | = 0.295 Iteration 6 = 37.53 | | = 0.065 We get,
|
= 37.69 | = 0.585
|
= 37.55 | = 0.145
|
= 37.52 | = 0.035
= 37.53 V , = 12.52 V and
43.75) = 39.0625
|
|
= 12.69 | = 0.585
|
= 12.55 | = 0.145
|
= 12.52 | = 0.035
= 12.89 | = 1.17
|
= 12.59 | = 0.295
|
= 12.53 | = 0.065
|
= 12.51 | = 0.02
= 37.52 V , = 12.51 V
The finite difference method represents the solution region by an array of grid points; its application becomes difficult with problems having irregularly shaped boundaries. Such problems can be handled more easily using the finite element method. METHOD OF MOMENTS
www.awabsir.com
Page 16
Still having doubt’s?
Awab Sir-89 76 104646
Like the finite difference method, the moment method or the method of moments (MOM) has the advantage of being conceptually simple. While the finite difference method is used in solving differential equations, the moment method is commonly used in solving integral equations. MoM uses integral method. The advantage of this method is that the order of the problem is reduced by one. For example a parallel plate capacitor is a 3-D domain. However, we will be working only on the surface of the capacitor plates, it becomes reduced to 2-D domain.
The potential at any point on the plate is a function of the charge distribution. Step 1: The charge can be given as
Step 2: To determine the charge distribution, we divide the plate section into smaller rectangular elements. The charge in any section is concentrated at the centre of the sectionThe potential V at the centre of any section is given by
Step 3: Assuming there is uniform charge distribution on each subsection we get This equation can be rearranged to obtain Where, [B]: column matrix defining the potentials www.awabsir.com
Page 17
Still having doubt’s?
Awab Sir-89 76 104646
[A]: square matrix In MoM, the potential at any point is the function of potential distribution at all points, this was not done in FDM and FEM. E.g. Using the moment method find the capacitance of the parallel plate capacitor shown in the figure below a = 1m, b = 1m, d = 1m, r = 1 z b a Vj
ρ1
+1V
Rij
d Vi
y -1V
x
Solution: Let the potential difference between the plates be V0 = 2v the top plate is maintained at +1v and bottom plate at -1v. we have to find C =
=
_______ (1)
to find charge Q = _______ (2) Use MOM to determine Divide P1 into n subsections , …….. , ……. The potential Vi at the center of any section is
and P2 into n subsections
Vi =
Vi =
www.awabsir.com
________ (3)
Page 18
Still having doubt’s?
Awab Sir-89 76 104646
Assuming there is uniform charge distribution on each subsection equation (3) can be written as Vi = where,
Aij =
Thus,
V1 =
________ (4)
=1
: : Vn =
=1 Vn + 1 = :
= -1
: V2n = = -1 From the set of 2n simultaneous equation with 2n unknown charge densities we get, [A] [ ] = [B] ______________ (5) hence, by matrix inversion [ ] = [A] -1 [B]
where, [B] is column matrix defining potentials [A] is square matrix containing elements Aij (xj , yj , zj )
z
y1
y2 y
x1
www.awabsir.com
Page 19
Still having doubt’s?
Awab Sir-89 76 104646
x2 (xi , yi , zi )
x
To determine Aij, we consider the two subsections I and j shown in the figure above (the subsections could be on the same plate or different plate). Aij =
____________(6)
where, Rij = [(xj - xi)2 + (yj - yi)2 + (zj - zi)2]1/2 x2 – x1 =
= y2 – y1
we get, Aij =
and
Aij =
=
ln (1 +
i
)=
j
(0.8814)
Using these equation values C can be determined.
Compare FDM, FEM and MoM
Basic principle
Advantage
Finite Difference Finite Element Method Method of Moments Method (FDM) (FEM) (MoM) Based on 1. Energy based Based on integral method differentiation i.e. (energy the differential minimization) equation is 2. Weighted residual converted to a (reducing the error) difference equation. 1. Simplest 1. Computationally 1. More accurate
www.awabsir.com
Page 20
Still having doubt’s?
Disadvantage
www.awabsir.com
Awab Sir-89 76 104646
method 2. Taylor series based
easier than MoM 2. Can be applied to unisotropic media. 3. [K] matrix is sparse
1. Need to have uniform rectangular sections (not possible for real life structures) 2. Outdated
1. Difficult as compared to FDM
(errors effectively tend to cancel each other) 2. Ideally suited for open boundary conditions 3. Popular for antennas 1. Mathematically complex
Page 21
