Coupling of field and circuit equations (Electrical Machine) Part 3


Solution of the system of coupled equations

Particular attention must be paid to the solution procedure of the coupled system of FEM and circuit equations. The matrix obtained for a 2D time-harmonic solution coupled with an electric circuit described with the proposed method is complex, symmetric, has no zero diagonal elements but is not hermitian.

The FEM block is positive definite and the circuit coupling block is negative definite. Therefore, the conjugate gradient (CG) method can not be used. Other suggestions are the bi-conjugate gradient (BiCG) method, the conjugate gradient method on the normal equations (CGN), other orthogonal Krylov-subspace methods and block elimination schemes (BES).

Conjugate gradient on the normalised equations

Instead of solving the system


the system


is solved. The matrix is hermitian.

Block elimination schemes

Due to the fact that the different parts of the matrix have different properties, it is advantageous to split up the matrix into several blocks. By matrix calculus, some block elimination schemes can be derived. The kind of elimination is dependent on the criterion to divide the matrix.

The system of equations can be partitioned as:


where A is a sparse matrix with the same dimension as the mesh, C is a full matrix with all circuit equations and B is a full matrix containing the coupling terms between the FEM magnetic field description and the circuit equations.

The system can be written as:


The vector Q and each column of P can be calculated by


using an iterative equation solver for sparse systems. The number of systems that has to be solved in this manner equals the amount of network unknowns plus one. This is an important reason for describing the electric circuit with as few unknowns as possible.

Another method separates the complex equations (these are the equations in the solid conductor regions and the electric circuit equations) from A. This results in dividing the system matrix in a real and a complex part.


The system of equations can be written as:


The system can be solved with the unknowns X,M,Yand N by iteration.

Two examples

In the first example, eddy currents are induced in a conducting plane, passing between two symmetric inductors (Fig. 5.88). The external circuit is shown in Fig. 5.82. In the second example an induction machine at start-up is analysed (Fig. 5.89). Stator windings and rotor bars are connected as shown in Fig. 5.90. The numbers of additional circuit equations for the different methods are shown in Table 5.8. The structure of the matrix is shown in (Fig. 5.90).

Table 5.8. Number of circuit equations.

circuit analysis










induction motor





Plot of the equipotential lines of an inductor and a conducting plane moving at 10 m/s to the right side.

Fig. 5.88. Plot of the equipotential lines of an inductor and a conducting plane moving at 10 m/s to the right side.

Equipotential plot of a 50 Hz time-harmonic solution for an induction motor.

Fig. 5.89. Equipotential plot of a 50 Hz time-harmonic solution for an induction motor.

Electric circuit of a) the stator and b) the rotor of an induction motor.

Fig. 5.90. Electric circuit of a) the stator and b) the rotor of an induction motor.

Matrix structure of a 50 Hz time-harmonic solution of the induction motor example.

Fig. 5.91. Matrix structure of a 50 Hz time-harmonic solution of the induction motor example.

Transient electromagnetic problem

The partial differential equation for a two-dimensional transient magnetic problem is:


To stay in the same notation as in section 5.7 using the relative reluctivitytmp101415_thumb[2][2]we can write:


and the functional is:


Functional within an clement

The functional within an element becomes


Evaluating the third term yields:


with the matrix entriestmp101422_thumb[2][2]



Time stepping

The system of equations in matrix-vector notation can be written by:


The time discretisation can be applied to the Galerkin approach in the time domain. The solution is only computed at discrete points in time, spaced in finite intervals A/, the time-steps. First order shape functions are chosen for A and T as functions in time.






If the Galerkin approach is applied withtmp101428_thumb[2][2]as the weighting function then


This corresponds to a standard central difference formula. By introducing a more general set of weighting functions other difference schemes can be obtained. (5.275) can be generalised by using the parametertmp101432_thumb[2][2]


Equation (5.275) can be written in a general form:


Withtmp101437_thumb[2][2]the forward difference Euler method is obtained. This approach is an explicit method because the term KA is evaluated at the beginning of time intervaltmp101438_thumb[2][2]


tmp101444_thumb[2][2]gives the backward difference fully implicit method since the term KA is evaluated at the end of the time intervaltmp101445_thumb[2][2]


The Crank-Nicolson scheme is obtained for


the Galerkin scheme of (5.275) is obtained.


the Galerkin scheme of (5.275) is obtained.

Particular attention must be paid to the choice of the value oftmp101454_thumb[2][2] because it influences the stability and the accuracy of the numerical results.

Stability and accuracy

Stability is important when time stepping is used for the solution of partial differential equations. By using the time stepping formula eq.(5.277), the-scheme is unconditionally stable fortmp101455_thumb[2][2]The stability does not prevent oscillations, but guarantees that oscillations do not grow out of control. An oscillation-free scheme is the fully implicit method withtmp101456_thumb[2][2]while for all valuestmp101457_thumb[2][2]the implicit scheme oscillates if the time-steptmp101458_thumb[2][2]is too large. The stability limit fortmp101459_thumb[2][2]has the form:




the diffusivity,tmp101468_thumb[2][2]the characteristic length of the elements and tmp101469_thumb[2][2]a numerical constant which depends on the elements used and on the choice oftmp101470_thumb[2][2]

The errortmp101471_thumb[2][2]in the approximation of the time derivative is fortmp101472_thumb[2][2] andtmp101473_thumb[2][2], i.e. for the fully implicit and the Euler explicit methods, of the order numerical constant which depends on the elements used and on the choice oftmp101479_thumb[2][2]




for the Crank-Nicolson method it is of the order:


Therefore, the choicetmp101485_thumb[2][2]is advantageous, since it corresponds to an implicit method which is unconditionally stable and gives thus second order accuracy for the time integration. On the other hand, the valuetmp101486_thumb[2][2]gives only first order accuracy for the time integration but completely avoids numerical oscillations even with large time steps. Thus a value oftmp101487_thumb[2][2]may be used, with a relatively small time step, in order to obtain an accurate solution. A value oftmp101488_thumb[2][2]can be used with a large time-step to obtain a less accurate solution to estimate the transient behaviour of the problem in principle.

Slow motion

We will discuss motion problems with a static magnetic field and with a uniform moving Cartesian geometry at relatively low speed v:


Problems considering high speed with the solution of this A-formulation will be discussed.

The first term of the differential equation is called the diffusion term. The second term is called the convection term. The right-hand-side is called the load.

Unlike the previous formulations, the variational principle corresponding to the solution of the differential equation is not known. Therefore, a weighted residual method has been applied.

Galerkin approach

The multiplication of the differential equation with the weighting functionstmp101494_thumb[2][2]yields:


The magnetic vector potential A is written in terms of the basis functionstmp101498_thumb[2][2]


In a Galerkin approach the same functions are chosen as weighting functions. The first term is then written as:


The third term can be evaluated by:


The second term is:


The Galerkin approach and the variational technique give the same result for the diffusion term and the load term. In the case of linear basis functions, the integral in the convection term becomes:


System of equations

The system of equations is given by:


This system is not symmetric because:


Appropriate iterative methods to solve the system of equations are minimal residual methods such as BiCG and GMRes.

GMRes has the problem that the memory requirement increases with the number of iteration steps. Therefore, restarted versions can be used. Petrov-Galerkin schemes such as BiCG are fast iterative equation solvers but have the disadvantage of breakdowns while solving the system.

External circuits

The current density in a moving conductor is calculated as


In the case of a stranded conductor, the current density is assumed to be constant. The voltage drop over a stranded conductor withtmp101509_thumb[2][2]turns moving at the speed v is


In the case of a solid conductor, the voltage is assumed to be constant. The current through the solid conductor is


If the currents through one or more stranded conductors or the voltages through one or more solid conductors are unknown, extra circuit equations are needed to describe the full behaviour of the model. The choice of the circuit unknowns, the construction of the extra equations and the coupling terms is done by the signal flow graph methods already described for the case of time-harmonic magnetic fields. The non-moving conductors, however, have to be treated in a slightly different way. No induced effects appear in these conductors so long as the geometry does not change while moving.


A conductive plane moves between two inductors. If the speed of the plane is zero, the solution corresponds to the static solution (Fig. 5.92). If motion is considered, the flux lines are pushed away in the direction of the motion (Fig. 5.93). If speed increases, the flux lines have less space to pass through the conducting plane. For very high speeds,

numerical problems result in unstable solutions (Fig. 5.94). In the figure, a separated closed flux line occurs. Up-winding schemes promise to surmount this difficulty. For further details please refer to the literature.

Static solution of a non-moving conductive plane.

Fig. 5.92. Static solution of a non-moving conductive plane.

Static solution of a conductive plane moving at 10 m/s to the right side.

Fig. 5.93. Static solution of a conductive plane moving at 10 m/s to the right side.

Static solution of a conductive plane moving at high speed to the right side.

Fig. 5.94. Static solution of a conductive plane moving at high speed to the right side.