Erfan Mostajeran, Navid Amiri, Seyyedmilad Ebrahimi, Juri Jatskevich
©SHUTTERSTOCK.COM/PHOTO SMILE
Electrical machines are extensively used in our everyday life. On the one hand, this can be seen in the rapid growth of generation from renewable energy sources such as wind, hydropower, tidal, etc. On the other hand, at the energy utilization and consumption end visible to most people, we are also witnessing revolutionary changes in many sectors, such as the electrification of all types of transportation, i.e., electric vehicles, industry-wide initiatives for more-electric aircraft and more-electric ships, military vehicles and defense systems, industrial automation, industrial and personal robots, medical devices and instruments, flying drones, electronic toys, and the multitude of household appliances and devices, all of which are designed and built with electric motors of various types and sizes.
In many applications, electrical motors and generators are controlled using switching power electronic converters that enable efficient and flexible operation and optimize energy conversion. Engineers and researchers develop computer models of electrical machines and conduct countless offline and real-time simulations to design and operate those systems and devices. Such simulations may be run on a small system with just a few machines and their converters and controllers or on a large system together with many other components and subsystems. Just as electrical motors and generators enable many products and technologies that improve human life, models of such electrical machines for studies in electromagnetic transient (EMT) simulators also make it possible to envision, design, and analyze/operate such products and devices effectively.
The models of electrical machines may be broadly categorized depending on their underlying purpose and level of fidelity. The first group includes the finite-element analysis (FEA)- or magnetic equivalent circuit (MEC)-based models, where the dimensions of equations may reach tens to hundreds of thousands, which are determined by the discretization of the machine’s geometry and the flux paths into the finite-element mesh. For the MEC models, this would be the reluctance segments whose number may still be significant on the order of several hundreds or thousands. Although these models can provide excellent accuracy necessary for design-oriented studies and optimization on an individual machine while considering the geometry and material properties to a significant level of detail, they are computationally costly.
The second group includes the so-called lumped-parameter coupled-circuit models, where the number of electrical windings in the machine determines the size of the equations. Since the total number of windings is typically small (e.g., one per phase on the stator/rotor, and if needed, a field winding plus damper windings, etc.), the resulting models have a much smaller number of equations and may be efficiently implemented in most EMT simulators for simultaneous studies of machines with power electronic converters as well as power systems with many generators and motors. It may be estimated that for one FEA-type model, thousands of lumped-parameter coupled-circuit models are implemented and simulated every day worldwide using various offline and real-time EMT simulators. Therefore, this article focuses on the lumped-parameter coupled-circuit models of electrical machines and their interfacing methods in EMT simulators.
The energy conversion process in electrical machines is one of the purest, cleanest, and most efficient known to humankind and is based on Maxwell’s equations for electromagnetic fields and Lorentz force. Electrical machines perform energy conversion based on the interaction between the magnetic fields and current-carrying conductors to generate torque. Generally, an electrical machine consists of a stator and a rotor circuit, each generating its own magnetic field through external sources, such as excitation, permanent magnets (PMs), magnetization through an external magnetic field, or induction. Different types of electrical machines exist that differ depending on the method by which the stator and rotor generate their magnetic fields. These electrical machines are commonly rotating devices (i.e., their rotor is in a circular motion), although linear motors/generators also exist where the motion is in a linear direction. Nevertheless, the modeling and analysis of both rotating and linear electrical machines are greatly similar, where one is carried out using circular coordinates and angles, and the other is carried out in linear coordinates and displacements. Electrical machines that are responsible for generating most of the electricity today and converting energy to mechanical work in various applications mainly include rotating synchronous and induction machines. Therefore, this article will mainly discuss the modeling of these machines in EMT simulation programs.
The lumped-parameter coupled-circuit models of electrical machines may be obtained directly by considering the magnetically coupled windings as rotor-position-dependent inductors and Kirchhoff’s circuit laws that consider the actual connection of all coils and conductors. This approach gave rise to the so-called coupled-circuit-phase-domain (CCPD) models that consider each stator (and rotor) phase winding in their respective physical variables and phase coordinates. Figure 1 illustrates examples of three-phase induction, field-wound synchronous, and PM synchronous machines (PMSMs) and their simplified coupled-circuit representations. According to this approach, a typical three-phase induction machine may be represented by three inductors for the stator and three inductors for the rotor, all of which are magnetically coupled, and the coupling between the stator and rotor windings also changes with the rotor position. The squirrel-cage induction motor may be easily represented by shorting the rotor windings in Figure 1(a), and the doubly-fed induction machine may be implemented if the rotor circuit terminals are not shorted and made available instead. The number of equivalent rotor windings in Figure 1(a) may be increased accordingly to represent more complicated rotors, such as squirrel-cage deep-bar conductors or two layers of conductors.
Figure 1. Typical three-phase rotating ac machines and their lumped-parameter coupled-circuit representation. (a) An induction machine. (b) A field-wound synchronous machine (SM) with equivalent damper windings on the rotor. (c) A PM synchronous machine (PMSM).
The conventional wound rotor SM may be represented similarly by using three coupled inductors for the stator windings and a field winding for the rotor excitation, as depicted in Figure 1(b). Additionally, the rotor may have a specially designed damper cage, which may be similar in construction to the squirrel-cage rotor conductors of an induction motor and act to dampen the electromechanical oscillations of the rotor with respect to the stator field. This damper structure may also be represented by equivalent windings placed on the direct and quadrature rotor axes (i.e., the d- and q-axes) of the SM. Depending on the modeling accuracy and fidelity, one or several equivalent damper windings may be assumed in the CCPD model of the SM, as depicted in Figure 1(b). Generally, one or two damper windings are considered sufficient for typical transient studies, and more equivalent damper windings may be included to better represent the stator electrical transients in the higher frequency range.
PMSMs are also gaining wide acceptance in power applications such as motors and generators ranging from several hundred watts to tens of megawatts. The CCPD model for such machines may be obtained by considering the stator phases as coupled inductors with a PM rotor, as illustrated in Figure 1(c). Similar to the conventional SMs of Figure 1(b), the rotor structural/magnetic saliency (as well as the damping effect of the rotor iron structure) may be considered by assigning appropriate values to the inductances and possibly adding some damper windings to the equivalent circuit representation.
The total number of equations describing the electrical subsystem of the machine is determined by the total number of equivalent windings considered in the coupled-circuit model. The mechanical subsystem is typically characterized by Newton’s laws of motion relating to the rotor’s mass or inertia, the angular position and speed, and all torques acting on the machine. Single-mass or multi-mass systems may be assumed depending on the modeling needs. The common properties of the CCPD models of the rotating machines depicted in Figure 1 are that their inductances will be variable and dependent on the rotor position and magnetic saturation phenomena if it is modeled. Thus, if implemented directly as it is defined, these models will generally require the recalculation of all the relevant inductances at every time step of the simulation, which adds to the numerical cost when considering these models in various EMT simulators. For a long time, this computational challenge prevented the wide use of the CCPD models directly and led to the development of the two-axis models as an alternative. However, with the development of more powerful computational hardware and efficient numerical techniques, the CCPD models are now being considered again in some EMT simulators.
The structure of the rotors of typical field-wound SMs determines the direct magnetic axis (the d-axis) and the orthogonal quadrature axis (the q-axis). This distinct geometrical rotor property (as well as the computational challenge with the variable inductances in the CCPD models) led the analysis of electric machinery from the early 20th century to express all machines’ electrical equations using transformed variables projected into the rotating qd system of coordinates. This is known as R. H. Park’s transformation and two-reaction theory, which transforms the original physical variables and equations of the CCPD models into the rotating qd reference frame, which is naturally fixed on the rotor for the SMs. Depending on the convention used, i.e., which axis leads the other, and the assumed direction of rotation, one may find qd- and dq-axes used equivalently in the literature. In the 1960/1970s, P.C. Krause generalized the reference frame theory for electrical machinery by relating all commonly used reference frames, showing their unique properties. This powerful tool allowed the expression of most ac machine models in the qd reference frame for analysis, control, and simulation purposes. The zero sequence is also included to uniquely transform balanced and unbalanced three-phase quantities, and thus, qd0 coordinates. The transformed qd0 equivalent circuits of the typical induction machines, field-wound SMs, and PMSMs are shown in Figure 2.
Figure 2. The qd0 equivalent circuit models of typical three-phase rotating ac machines. (a) An induction machine. (b) A field-wound SM with equivalent damper windings on the rotor. (c) A PMSM.
Figure 2(a) depicts the qd0 model of induction machines that assumes a shorted squirrel-cage type rotor. The zero-sequence circuit is also shown on the stator side to account for possible zero-sequence current if the machine stator windings are assumed to be Y-connected with an accessible neutral point. If the neutral point is not used, which is often the case, then the zero-sequence circuit can be removed from consideration. The doubly-fed induction machine model is obtained by opening up the rotor circuits and making their terminals available for external connections. The SM model is shown in Figure 2(b), where one field winding (on the d-axis) and three equivalent damper windings (one on the d-axis and two on the q-axis) are assumed on the rotor side. The stator sides of the induction and SM qd0 models look similar, except for the point that induction machines are typically round and result in equal magnetizing inductances in the q- and d-axes, but those inductances may be different for the SMs depending on their rotor construction (i.e., round or salient-pole rotor type).
The PMSM qd model may be obtained by simplifying the equivalent circuit to that in Figure 2(c), where the rotor is represented by a source to account for the PM rotor flux. Similar to the conventional SMs of Figure 2(b), the rotor structural/magnetic saliency (as well as the damping effect of the rotor iron structure) may be considered by appropriate values for the inductances (and possibly adding some damper windings to the equivalent circuit).
An essential property of the transformed equivalent circuits in Figure 2 is that since the q- and d-axes are perpendicular, the magnetic coupling is removed and replaced by the so-called coupled speed-voltage sources. This, in turn, results in constant (i.e., rotor-position-independent) inductances that describe the magnetically linear qd0 machine models. This property also greatly simplifies the computations needed for executing such models in computer simulations, which made the qd0 models by far the most widely used even today. As the three-phase ac side is typically symmetric and balanced, the overall number of equations may also be reduced from a typical three per stator side (as in the CCPD model) to two per stator side (one for q-axis and one for d-axis), assuming that the ac side does not have a zero sequence. These features make the qd0 models very compact and computationally efficient as standalone components.
However, when interfacing with external systems, the qd0 models require special consideration since their internal calculations are carried out in transformed qd variables, but all the external systems are typically represented in the actual physical variables and abc phase coordinates. The challenge introduced by this interface may also differ depending on the EMT simulation program considered.
EMT simulation programs are typically the preferred choice when studying electrical power systems in the time domain with high accuracy, e.g., detailed short circuit transient analysis, overvoltage transients, controller tuning and optimization, motor startup transients, generator synchronization, the coordination and interaction of converter-based resources, etc. The EMT simulation programs commonly used in the power industry may be broadly categorized based on their underlying solution methods into state-variable-based (SVB) and nodal-analysis- (or modified-nodal-analysis)-based EMT programs (EMTP type). Each category also includes offline and real-time industry-grade simulators. Examples of widely used SVB programs include MATLAB/Simulink, Simscape Electrical toolbox, and PLECS for offline simulations as well as Opal-RT and Typhoon HIL for real-time simulations. Examples of EMTP packages include EMTP-RV, ATP-EMTP, MicroTran, and the widely used PSCAD/EMTDC for offline simulations and RTDS for real-time simulations.
Some programs include modules based on both solution approaches and allow their interfacing, such as DIgSILENT. Both families of SVB and EMTP simulators have their advantages and limitations (due to solving the electrical circuits differently). These programs are widely used in industry and academia. Generally, the models of electrical machines are crucial components for most power systems studies.
The SVB EMT simulators represent electrical system dynamics based on ordinary differential equations (ODEs). The program formulates the ODEs in state-space form or as a collection of state equations based on the state variables of the system, e.g., the currents of inductors and voltages of capacitors as the states of the energy storage circuit components. Electrical systems often contain linear and nonlinear subsystems (e.g., circuit elements, controllers, etc.). Models of such subsystems may be separated into linear and nonlinear blocks, as illustrated in Figure 3. Each subsystem may accept inputs from external sources and exchange the coupling variables with other subsystems.
Figure 3. System partitioning into linear and nonlinear models in SVB EMT simulators.
Once the network topology is determined, a state-space model may be formed for the entire linear part of the circuit (for which the standard A, B, C, and D matrices may be created). However, the nonlinear or varying parts are typically solved in their general form without forming fixed matrices to minimize additional calculations. The time-domain solution is then calculated numerically by integrating the state-space equations using either fixed- or variable-time-step ODE solvers embedded in the program. The SVB EMT simulators commonly offer a variety of ODE solvers and user-defined settings to optimize the simulation performance and accuracy. However, due to computational cost and scalability challenges, such programs are typically used for simulations of smaller-scale electrical power systems, e.g., for control designs of electric drives, interconnection studies of small-scale renewable generation with few electrical machines, etc.
Implementing models of electrical machines in SVB EMT simulators is based on numerically integrating the corresponding voltage equations (either in the qd0 form or in the CCPD form), assuming either winding currents or fluxes as the state variables. Due to their inductive nature, when formulating a proper state model that is solved by numerical integration (as opposed to differentiation), these models require the machine’s terminal voltages as the inputs and compute the windings’ currents as the outputs. Thus, the natural interface of the machine model is formed by a voltage-controlled current source. This voltage-controlled current source requirement also brings challenges when interfacing the classical qd0 machine models with the external circuits. In most EMT simulators, the components and systems outside of the electrical machines are modeled and implemented in the actual physical variables (i.e., in the abc phase coordinates for the conventional three-phase systems). Depending on the composition of the external circuit to which the machine is connected, the interface may be compatible or incompatible.
Figure 4(a) shows an example of the interface, where the qd0 machine model is interfaced with an external network that receives currents as input and computes voltages as output of its submodel. The terminal voltages ${v}_{abc}$ are transformed from the abc phase coordinates into the qd0 coordinates and fed to the qd0 model as ${v}_{qd0}$, together with mechanical torque ${T}_{m}$, speed ${\omega}_{r}$, and position ${\theta}_{r}$ from the mechanical subsystem. The machine model computes the phase currents ${i}_{qd0}$ that are then transferred to the abc physical coordinates and fed back to the external subsystem as ${i}_{abc}.$ This interface is compatible since the machine model requires phase terminal voltages as input and computes winding currents as output. Examples of networks that receive current and output voltage are circuits composed of capacitors, resistors, voltage sources, etc.
Figure 4. Interfacing of the qd0 model of ac electrical machines in SVB EMT simulators. (a) A compatible interface with the external capacitive circuit. (b) An incompatible interface with the external inductive circuit. (c) A compatible interface with the external inductive circuit using artificial snubber resistors.
However, in most practical applications, electrical machines are connected to other network elements that are also inductive, such as lines, cables, transformers, or power electronic converters with a topology that changes depending on the state of its switching devices. The typical numerical solution of such networks also requires voltage as the input, and it computes current as the output, just like the qd0 machine model. Examples of networks that require input voltages and output currents are inductive circuits, nonlinear models/elements, current sources, or qd0 models of other electrical machines. An example of such an interface is shown in Figure 4(b), where the qd0 model is assumed to be connected to an inductive external network. This interface with the qd0 machine model is incompatible because of the input–output variable mismatch. Specifically, the voltages ${v}_{abc}$, which are required as the input to the machine model and the external circuit, remain unknown.
Sometimes, this problem may be solved by invoking iterations on the subsystems involved, which increases the computational cost and may have a convergence problem. The iterative approach may be used only for offline EMT simulators because such iterations are not possible in real time. However, in most SVB EMT simulators, this problem is solved using artificial snubber circuits placed at the interfacing terminals, as illustrated in Figure 4(c). Therein, with the addition of the snubber resistors ${r}_{x}$, the output currents of the ac machine (denoted as ${i}_{abc1})$ and the external network (denoted as ${i}_{abc2})$ can flow into the snubber. This enables establishing the interfacing voltages ${v}_{abc}$ required by both subsystems to make the interface compatible. Such snubbers may be realized using simple resistors or capacitors as their equations permit calculating the terminal voltages based on the input currents.
In nodal-analysis-based EMT programs (i.e., EMTP-type programs), the equations of individual circuit elements are discretized at the component level using the trapezoidal integration rule with a fixed simulation time step $\Delta{t}$. For example, the voltage-current relationships for basic RLC components in continuous and discrete-time representation are demonstrated in Figure 5.
Figure 5. A representation of basic RLC components for the EMTP nodal solution. (a) Voltage-current relationships in the continuous-time domain. (b) Voltage-current relationships in the discrete-time domain based on trapezoidal integration rule. (c) Equivalent circuits in the discrete-time domain.
After discretization, each component is represented as a conductance in parallel with a so-called history current source (which is equal to zero for resistors). It is noted that the history terms depend on the solution at the previous time step (i.e., ${t}{-}{\Delta}{t}$). Nonlinear and varying components can also be modeled using different techniques, e.g., piecewise-linear representations, current injections with a time-step delay, compensation-based methods, etc. All the resultant voltage-current relationships are organized in the nodal form as ${\bf{Gv}}_{n} = {\bf{i}}_{h}$. Here, ${\bf{G}}$ is the conductance matrix of the entire network, ${\bf{v}}_{n}$ is the vector of nodal voltages, and ${\bf{i}}_{h}$ is the vector containing the sum of the history and independent current sources injected at the nodes. In this formulation, each component is represented by the entries in the network conductance matrix ${\bf{G}}$ and in the voltage and current vectors ${\bf{v}}_{n}$ and ${\bf{i}}_{h}$.
Figure 6(a) illustrates the interface of a three-phase ac machine at the nodes a, b, and c to the network in EMTP. The discretized ac machine model is interfaced using the Norton equivalent circuit shown in Figure 6(b), with ${\bf{G}}_{abc}$ (the equivalent conductance matrix of the subject electrical machine) as the interfacing submatrix. Connecting the machine model to the external network at nodes a, b, and c is visualized in Figure 6(c) by writing the values of the machine submatrix ${\bf{G}}_{abc}$ to the appropriate entries in the overall ${\bf{G}}$ matrix as well as the ${\bf{v}}_{n}$ and ${\bf{i}}_{h}$ vectors.
Figure 6. The interfacing of ac machine models in EMTP programs. (a) A three-phase machine is connected to the external network at interfacing nodes a, b, and c. (b) The Norton equivalent interfacing circuit of the machine discretized model. (c) The nodal equation of the discretized network and addition of equivalent submatrix/subvectors of the ac machine model. (d) The modified augmented nodal-analysis formulation.
The nodal equation ${\bf{Gv}}_{n} = {\bf{i}}_{h}$ is solved at every time step ${\Delta}{t}$ for the unknown voltages in ${\bf{v}}_{n}$ using efficient numerical techniques (which avoid the direct matrix inversion), after which the currents of branches are calculated. This is known as H. Dommel’s algorithm, which enabled the efficient and scalable EMT simulation of power networks. It is important to note that for a network with time-invariant components, the ${\bf{G}}$ matrix will be constant, and the program needs to calculate it only once. Otherwise, the ${\bf{G}}$ matrix must be updated whenever the values of the circuit components change; thus, solving ${\bf{Gv}}_{n} = {\bf{i}}_{h}$ will require more computations.
Some offline EMT simulation programs (e.g., EMTP-RV) use the so-called modified augmented nodal analysis (MANA). Therein, the network equations are organized in the general form shown in Figure 6(d), where the vectors ${\bf{v}}_{n}$ and ${\bf{i}}_{h}$ as well as the ${\bf{G}}$ matrix correspond to their classical nodal-analysis counterparts, while there are several other matrices and vectors that are augmented to include circuit components that cannot be modeled as a conductance, e.g., ideal voltage and current sources, current-dependent circuit elements, and/or ideal switches and transformers. In MANA, the voltage/current relations of such elements are considered by the row, column, and diagonal block submatrices ${\bf{A}}_{r}$, ${\bf{A}}_{c}$, and ${\bf{A}}_{d}$ and by the vectors of unknowns ${\bf{x}}$ and knowns ${\bf{b}}$, respectively, as shown in Figure 6(d). For example, for ideal voltage (or current) sources, ${\bf{A}}_{r}$ and ${\bf{A}}_{c}$ would be incidence matrices having either 0 or ${\pm}{1}$ values possessing information about the connection of sources between nodes, ${\bf{b}}$ would include the values of the sources, and ${\bf{x}}$ is the currents (or voltages) of the sources.
For the ideal transformers, these submatrices would also include the turn ratios as the scaling factors. Ideal and nonideal switches can also be included using the new submatrices ${\bf{A}}_{r}$, ${\bf{A}}_{c}$, and ${\bf{A}}_{d}$ without changing the rank of the overall matrix. Although the MANA approach presents advantages compared to the classical nodal analysis, it may require iterations for the solution of nonlinear devices and electrical machines. Generally, it is preferred to solve the nodal equations noniteratively in real-time EMT simulators.
For implementation in EMT simulators, the models of ac machines are also discretized using the trapezoidal rule with a time step ${\Delta}{t}$. Moreover, since the classical qd0 models are derived in transformed variables and coordinates, obtaining their equivalent discretized Norton equivalent interfacing circuit with ${\bf{G}}_{abc}$ is not straightforward, and several approaches are possible, generally requiring some form of approximation. Different solutions have been proposed in the literature to overcome the incompatibility issues of the qd0 models in the EMTP solution. One approach to achieve the interfacing is using the stator voltages from the previous time step (i.e., at ${t}{-}{\Delta}{t}$) and modeling the machine as a Norton equivalent current source in parallel with artificial snubbers, as illustrated in Figure 7(a). This approach results in a fixed interfacing conductance ${\bf{G}}_{snub}$, and it is used in programs such as PSCAD. Another approach may be to average the unequal parameters of the discretized qd0 model to obtain a constant ${\bf{G}}_{abc}$ in abc coordinates and to use some form of prediction to compensate for the error. This approach is illustrated in Figure 7(b) and is used in MicroTran.
Figure 7. The interfacing of qd0 models of electrical machines in the EMTP solution. (a) Using snubbers and stator voltages from the previous time step. (b) Using the equivalent conductance matrix and some approximations.
Since the transformations between the qd0 and abc interfacing variables are nonlinear, it is also possible to iteratively solve the machine-network interfacing within the MANA formulation in offline EMT simulators. The so-called compensation-based approach may also be applied, where a Thevenin equivalent of the external network in qd0 coordinates is obtained, and the machine is interfaced and solved with the external grid in qd0 coordinates. This technique can be found in the universal-machine model in ATP. All these methods achieve the desired interfacing of the qd0 models in EMT simulators but possess different numerical properties, and their accuracy generally quickly degrades when the simulation time step is increased. Because of that, some researchers turned to discretizing the CCPD machine models and using them in EMT simulators directly. As it turned out, the discretized CCPD model naturally has a constant ${\bf{G}}_{abc}$ matrix for a magnetically linear induction machine. For SMs, however, this matrix is rotor-position-dependent.
To address the challenges of interfacing the traditional qd0 models of ac machines in EMT simulation programs, the voltage-behind-reactance (VBR) models have been developed more recently. The basic idea of these types of models is to separate the stator equivalent interfacing circuit (in physical variables and abc phase coordinates) from the rotor subsystem (which uses the transformed qd variables). S. D. Pekarek first pioneered this idea in the late 1990s. This approach allows interfacing the machine model using the VBR equivalent circuit, as shown in Figure 8. Depending on the considered ac machine and the VBR model type, the interfacing inductances ${\bf{L}}^{''}_{abc}$ may be variable, i.e., dependent on the rotor position and magnetic saturation, or made constant (with some approximations), resulting in constant-parameter VBR (CPVBR) models that are particularly advantageous for EMT simulation programs.
Figure 8. The interfacing of the VBR models of ac electrical machines in SVB EMT simulators.
Due to its VBR structure, in SVB EMT simulators, this model is interfaced using current-controlled voltage sources and resistor and inductor elements, as depicted in Figure 8. The interfacing VBR circuit is solved together with the external circuit. The machine’s stator phase currents ${\bf{i}}_{abc}$ (which are also the outputs of the external circuit in this system of Figure 8) are inputs to the rotor subsystem model of the ac machine. Other inputs include mechanical torque ${T}_{m}$, rotor speed ${\omega}_{r}$, and position ${\theta}_{r}$, all of which come from the mechanical subsystem. The VBR model computes the so-called subtransient back electromotive force (emf) voltages ${\bf{e}}^{''}_{qd0}$, which are then transformed to the abc coordinates as ${\bf{e}}^{''}_{abc}$ and used as inputs for the controlled voltage sources behind the reactance in the interfacing circuit, as depicted in Figure 8. This allows connecting the VBR models to arbitrary configurations of external circuits directly without artificial snubbers, as was needed in Figure 4(c).
For numerically efficient simulations, it is desirable to have constant (i.e., time-invariant) series impedances in the VBR models. Otherwise, the program must calculate the values of the interfacing branch parameters at every time step of the simulation, which can significantly increase the computational cost. For this purpose, a family of CPVBR models has been developed for induction and synchronous machines. Figure 9 shows a unified constant-parameter decoupled VBR RL-branch equivalent circuit for interfacing ac induction and synchronous machine models with external electrical networks. This circuit includes an optional zero-sequence branch, which is needed if the neutral point of the Y-connected stator winding is actually used in the studies; otherwise, this branch is simply removed. For round and symmetrical induction machines, this circuit is derived without approximations and is composed of constant and decoupled inductors ${L}_{D}$ and resistors ${r}_{D}.$ For SMs, some approximations may be required to arrive at constant-parameter and decoupled interfacing elements with explicit noniterative solutions, but the same interfacing circuit shown in Figure 9 is achieved.
Figure 9. A unified CPVBR RL-branch equivalent circuit for interfacing induction and synchronous machine models with external electrical network circuits.
A simple system composed of an induction motor fed from a source with a small inductive impedance is considered in Figure 10 to give the reader an idea about the impact of the interfacing of different models. Two models are considered to simulate a startup acceleration transient of the energized motor. The traditional qd0 model is implemented based on the equivalent circuits presented in Figure 2(a) and is interfaced with the inductive source according to Figure 4(c) using large resistive snubbers. The VBR model is implemented using the interfacing circuit of Figure 8 and is directly connected to the inductive source with no snubbers. Choosing the values for the snubbers is always a compromise decision between the modeling accuracy and the numerical performance of the simulation. Generally, one wants to select the snubbers as large as possible to make the snubber current very small, and thus, to make the model more accurate. However, this increased accuracy comes at the cost of making the model numerically stiff and possibly leading to smaller/shorter time step sizes (more time steps) and longer computation times for offline simulations. Snubbers may also limit the largest permissible time-step size for the real-time simulators.
Figure 10. A single induction motor fed from a voltage source with an inductive impedance emulating line cable. Snubber resistors are used to interface the conventional qd0 model with an inductive source.
Figure 11(a) shows the rotor phase a current obtained with the qd0 and VBR models compared to the reference solution. It is more visible in the magnified view in Figure 11(b) in which, even with the appropriately selected snubber resistors, the qd0 model has a visibly higher error and also requires more simulation time steps, making this model computationally more expensive. At the same time, the directly-interfaced VBR model uses much larger time-step sizes (i.e., fewer steps) and produces the solution points that land precisely on the reference trajectory. So, the demonstrated benefits of the VBR model include achieving better accuracy and permitting large time steps.
Figure 11. (a) The rotor current during startup acceleration transient as predicted by the qd0 and VBR models. (b) Its magnified view showing the difference between the models. Ref.: reference solution.
Including magnetic saturation is often considered to improve the fidelity of ac machine models in commercial EMT simulation packages. The saturation characteristic can be easily obtained from no-load tests, and it captures the nonlinear relationship between the magnetizing current ${i}_{m}$ and the main magnetizing flux magnitude ${\lambda}_{m}$ (typically the fundamental component). This nonlinear characteristic is depicted in Figure 12(a), and in the EMT simulators, it is typically represented using several piecewise linear segments. For induction machines, the saturation is typically applied to the main flux path, and since they are round and symmetrical, the saturation is incorporated into the qd0 models by considering the projections of the saturated flux vector ${\lambda}_{m}$ onto the q- and d-axes (therefore taking into account the so-called qd-axes cross-saturation), as depicted in Figure 12(b).
Figure 12. Modeling saturation in three-phase ac machines. (a) A typical nonlinear magnetizing flux–current relationship. (b) Projections of saturated magnetizing flux onto q- and d-axes in a round rotor machine. (c) Using the saliency factor to convert a salient-rotor machine to an equivalent round machine for representing magnetizing flux saturation.
For SMs, depending on the rotor construction, the saturation is sometimes applied to the d-axis only, to the q- and d-axes independently, or to the main magnetizing flux ${\lambda}_{m}$. Experiments and FEA suggest that saturation of the main flux magnitude results in better accuracy. For the qd0 model of the salient-pole SMs, the single-saturation factor or the two-saturation factor approaches can be utilized depending on the availability of the q-axis and intermediate directions of the saturation characteristic. However, in practice, only the d-axis saturation characteristic can be easily obtained experimentally from the open circuit tests. Therefore, a single saliency factor approach is often considered sufficient for the general-purpose models in EMT simulators. The saliency factor is defined as the square root of the ratio of the q- and d-axes magnetizing inductances ${L}_{mq}$ and ${L}_{md}$, as ${SF} = {\sqrt{{L}_{mq} / {L}_{md}}}$, and it is assumed to be constant for unsaturated and saturated conditions. Using this approach, the saturation is then implemented in the qd0 models by rescaling the magnetizing flux and current projections of the q-axis by the factor ${SF}$, as illustrated in Figure 12(c), which effectively converts the original salient rotor into an equivalent round rotor.
The described saturation approach assumes that the magnetizing flux wave is sinusoidally distributed in the machine’s air-gap and that this entire sinusoidal wave is being reduced due to saturation. The magnitude ${\lambda}_{m}$ of this flux wave follows the saturation characteristic presented in Figure 12(a). However, it has been shown in the literature that the flux distribution can become considerably distorted as a result of different saturation levels in different parts/zones of the machine, which creates spatial harmonics in the distribution of the magnetizing flux. A common cause is a deeper saturation in stator/rotor teeth/slots along the magnetization flux path compared to the sides. As a result, the magnetic flux density under saturation will have a flat-top waveform containing odd harmonics, as depicted in Figure 13(a). Therefore, the air-gap flux spatial harmonics can also be considered for higher modeling fidelity. Figure 13(b) demonstrates the fundamental component ${\lambda}_{m1}$ and the third harmonic ${\lambda}_{m3}$ for a typical induction machine as a function of the magnetizing current ${\boldsymbol{i}}_{m}$. Thus, it is possible to consider the air-gap flux harmonics to increase the accuracy of representing saturation in ac machine models in EMT simulators.
Figure 13. Considering air-gap flux harmonics when modeling saturation in ac machines. (a) A flat-top magnetizing flux waveform and its fundamental and third harmonic components. (b) The saturation characteristic of the magnitude of the fundamental and third harmonics flux versus the magnetizing current.
An induction machine example is considered here to illustrate the representation of saturation. Recent developments in advanced machine drives have given rise to robust sensorless control strategies based on monitoring certain harmonics in the machine’s voltages and currents resulting from magnetic saturation to pinpoint the magnetic flux vector position. For simulations of such drive systems in EMT simulators, it may be important to consider saturation with air-gap flux harmonics. While a wide range of harmonics can be present, the third harmonic is typically considered due to its dominant effect. Magnetic saturation may be implemented in most models, including CCPD, qd0, and variable- and constant-parameter VBR models, as illustrated in Figure 14. In all the considered models, the third harmonic ${\lambda}_{m3}$ may also be considered as an option, on top of the saturation of the fundamental magnetizing flux ${\lambda}_{m1}$. Each model considered in Figure 14 has unique interfacing properties, impacting the computational performance and numerical accuracy.
Figure 14. The implementation of the saturable (a) CCPD, (b) qd0, (c) VBR, and (d) CPVBR models of induction machines and their interface with an external network.
In the CCPD model, all the stator and rotor windings are represented as rotor-position-dependent and coupled inductive branches, as illustrated in Figure 14(a), that are solved together with the external network. The saturation is considered appropriately by defining all inductances in terms of the magnetizing current ${i}_{m}$ and the rotor position ${\theta}_{r}$. The computational cost of this model will be high due to the variable nature of all inductances.
In the qd0 model, the saturation may be represented using a flux correction method that follows the saturation characteristic of Figure 12(a) and (b) for the fundamental magnetizing flux component ${\lambda}_{m1}$. The implementation and interfacing of this model are illustrated in Figure 14(b). If the air-gap flux harmonic is considered, additional subcircuits, including the zero-sequence circuit, have to be added on top of the traditional equivalent circuit depicted in Figure 2(a) to include the effect of ${\lambda}_{m3}$. Cross-coupling between the fundamental component ${\lambda}_{m1}$ and its third harmonic ${\lambda}_{m3}$ is also considered. Snubber resistors are typically required for interfacing with this model, as depicted in Figure 14(b).
The VBR model may also be implemented with the saturation of the main magnetizing flux ${\lambda}_{m1}$ and an option of including the third harmonic ${\lambda}_{m3}$ of the air-gap flux. A direct extension from the saturable qd0 model results in a VBR model with saturation-level-dependent interfacing inductances. This property increases the computational cost, similar to the saturable CCPD model, although with only three variable coupled inductors instead of six. With some approximations, a CPVBR model may also be derived. These models and their interface are depicted in Figure 14(c) and (d).
An experimental study with a 0.25-hp squirrel-cage induction motor is presented here to demonstrate the performance of the saturable models. The experimental setup is illustrated in Figure 15(a). The measured and simulated steady-state voltages and currents for different levels of applied voltage are shown in Figure 15(b). It can be observed in Figure 15(b) that, as the voltage magnitude increases from 1 per unit (p.u.) to 1.2 p.u. and then to 1.4 p.u., the waveforms become more distorted, which is captured by the considered models. Here, all the considered models use small time steps and produce results consistent with the experiments.
Figure 15. Validation of induction machine models with saturation. (a) The experimental setup. (b) Steady-state voltage and current waveforms under different levels of saturation due to the increased voltage level.
A startup transient is simulated next, with the results presented in Figure 16. As can be seen in Figures 15 and 16, the currents are visibly distorted due to harmonics, which all models capture well along with the mechanical variables. However, Figure 16(b) shows the difference among the solutions produced by the considered models, particularly that the qd0 model requires very small time steps (i.e., a large number of steps) while being slightly off in terms of accuracy due to its snubber interface. It was also found that the CPVBR model may be about three times more efficient per step compared to the variable VBR model while achieving similar accuracy.
Figure 16. The induction machine startup transient predicted by the qd0, VBR, and CPVBR models under saturation. (a) The stator phase current and (d) its magnified view. (b) The electromagnetic torque and (e) its magnified view. (c) The mechanical speed.
The role of the interfacing of ac machine models in EMTP may be illustrated by considering a single-machine infinite-bus system with a large (835-MVA) steam turbine generator subjected to a three-phase fault at the terminals. The saturation is not considered here to focus on the interface only. The considered models include variable-parameter CCPD and VBR models and two versions of the qd0 model (qd0_1 is iterative using EMTP-RV, and the second model is qd0_2 with prediction using MicroTran). It is important to mention that when using time steps on the order of 50 µs to 100 µs, all models produce reasonable results visually close to each other, with up to 2% numerical error with respect to the reference solution. Therefore, the simulation step size is increased to 1 ms to exemplify the difference among the considered models.
The result of this transient in the stator phase a current computed by the considered models is shown in Figure 17. As seen in this figure, the results of different models are visibly different. The CCPD and VBR models achieve a direct simultaneous solution, and their results are closest to the reference even at such a large time step. However, the two qd0 models produce solution points that significantly deviate. It was also confirmed that the qd0 model interfaced using snubbers and delayed stator voltages is not numerically stable with this large time step and is therefore not included here. It is concluded that the variable-parameter CCPD and VBR models possess better accuracy, although their computational cost per step may be high due to the variable equivalent conductance matrix.
Figure 17. The three-phase fault transient study of an 835-MVA steam turbine generator computed by several models. (a) The stator phase a current. (b) A magnified view showing the differences among the models.
The advantageous interfacing properties of the CCPD and VBR models in EMTP brought more attention to these modeling approaches. However, when these models are derived for synchronous machines, their interfacing conductance matrix may become rotor-position-dependent. Moreover, considering the nonlinear saturation characteristic as a piecewise-linear function additionally makes the interfacing conductance matrix saturation-segment/level-dependent.
This property inevitably adds computational cost to implementing the CCPD and VBR models and hinders their application, particularly in real-time simulators where it is highly desirable to have a constant overall network conductance matrix for efficient and noniterative solutions. Fortunately, it is also possible to derive a constant interfacing conductance matrix using careful approximations. The procedure for developing constant-parameter CCPD and VBR machine models for EMTP is summarized in steps illustrated in Figure 18. Starting from the discretized Thevenin equivalent circuit, the first step is to separate the interfacing resistance matrix of the model into two parts—one constant part and one varying (rotor-position- and saturation-dependent) part. Next, the varying part is moved into the history voltage sources, which also requires predicting/approximating the stator currents. Finally, the new Thevenin equivalent circuit is converted to the Norton equivalent with a constant conductance matrix and appropriately modified history terms.
Figure 18. Steps for obtaining the CPVBR and CCPD models for EMTP. (I) Decomposing the original resistance matrix into constant and varying parts. (II) Moving the varying part into the history terms using predicted stator currents. (III) Obtaining the constant-parameter Norton equivalent circuit for interfacing in EMTP simulators.
The merits of possessing a constant-parameter interfacing circuit are demonstrated on an IEEE 39-bus system consisting of 18 synchronous generators, as shown in Figure 19. The system has been implemented in PSCAD, where all generators have been implemented using CCPD, VBR, and CPVBR models with the saturation of the main magnetizing flux, as presented in Figure 12(c). A transient study of a three-phase-to-ground fault at bus 25 is considered. The 2-norm error of the fault current contributed by the generators G17 and G18 connected to bus 37 is considered and reported in Table 1. The simulation times for conducting a 0.5-s study in PSCAD are provided in Table 2. Table 1 shows that all the considered models produce very similar accuracy for time steps of 100 µs and 500 µs, which confirms the approximation steps in Figure 18. However, Table 2 shows that the CPVBR model significantly outperforms the variable-parameter CCPD and VBR models by a factor of 4.7 in terms of computational efficiency. This computational gain is demonstrated due to the constant overall conductance matrix of the multimachine system, testifying to the advantage of having constant-parameter machine models.
Figure 19. The single-line diagram of the modified IEEE 39-bus system comprising 18 synchronous generators.
Table 1. Average 2-norm error of fault current contributed by synchronous generators connected to bus 37 for different time-step sizes expressed in percentage (%).
Table 2. CPU time required in PSCAD for simulation of 0.5-s short circuit transient study.
For EMT simulators, it is crucial to have computationally efficient and accurate models of ac electrical machines. Such models should be capable of representing essential electromagnetic and mechanical transient phenomena for studies of individual machines with their controllers and power electronic converters and for system-level studies with many machines participating as generators and motors. The typical machine models include a representation of magnetic saturation, which may be important for a wide range of studies. In addition to the traditional two-axis qd0 models, the CCPD and VBR models may offer interfacing advantages, translating into numerical robustness and good accuracy even with large time-step sizes. When such models are implemented in EMT simulators, it is also important that they are realized in the most efficient way, requiring the smallest possible number of calculations. Then, such ac machine models may enable extending the application of modern offline and real-time EMT simulators to larger power grids while optimizing the use of computational resources.
Electrical machines with more than three (e.g., five, six, nine, twelve, etc.) electrical phases are becoming attractive for some applications, offering increased reliability, fault-tolerant operation, and reduced losses. Models of such machines in EMT simulations are also getting attention. Further improvement of the models’ fidelity, numerical accuracy, and computational efficiency represents ongoing research with significant potential impact and benefit to many simulation software vendors, utilities, and consulting companies.
H. W. Dommel, EMTP Theory Book. Vancouver, BC, Canada: MicroTran Power System Analysis, May 1992.
P. C. Krause, O. Wasynczuk, S. D. Sudhoff, and S. Pekarek, Analysis of Electric Machinery and Drive Systems, 3rd ed. Piscataway, NJ, USA: IEEE Press, 2013.
L. Wang et al., “Methods of interfacing rotating machine models in transient simulation programs,” IEEE Trans. Power Del., vol. 25, no. 2, pp. 891–903, Apr. 2010, doi: 10.1109/TPWRD.2009.2039809.
E. Levi, “Saturation modelling in d-q axis models of salient pole synchronous machines,” IEEE Trans. Energy Convers., vol. 14, no. 1, pp. 44–50, Mar. 1999, doi: 10.1109/60.749146.
S. D. Pekarek, O. Wasynczuk, and H. J. Hegner, “An efficient and accurate model for the simulation and analysis of synchronous machine/converter systems,” IEEE Trans. Energy Convers., vol. 13, no. 1, pp. 42–48, Mar. 1998, doi: 10.1109/60.658202.
Y. Xia et al., “An efficient phase domain synchronous machine model with constant equivalent admittance matrix,” IEEE Trans. Power Del., vol. 34, no. 3, pp. 929–940, Jun. 2019, doi: 10.1109/TPWRD.2019.2897612.
N. Amiri, S. Ebrahimi, and J. Jatskevich, “Saturable voltage-behind-reactance models of induction machines including air-gap flux harmonics,” IEEE Trans. Energy Convers., vol. 38, no. 3, pp. 2022–2033, Sep. 2023, doi: 10.1109/TEC.2023.3259385.
Erfan Mostajeran (emostajeran@epeconsulting.com) is with Electric Power Engineers LLC, Austin, TX 78738 USA.
Navid Amiri (navida@ece.ubc.ca) is with the Electrical and Computer Engineering Department, University of British Columbia, Vancouver V6T 1Z4, Canada.
Seyyedmilad Ebrahimi (ebrahimi@ece.ubc.ca) is with the Electrical and Computer Engineering Department, University of British Columbia, Vancouver V6T 1Z4, Canada.
Juri Jatskevich (jurij@ece.ubc.ca) is with the Electrical and Computer Engineering Department, University of British Columbia, Vancouver V6T 1Z4, Canada.
Digital Object Identifier 10.1109/MELE.2023.3320508
2325-5897/23©2023IEEE