Jean Mahseredjian, Mohammed Naidjate, Mehdi Ouafi, Juan Antonio Ocampo Wilches
©SHUTTERSTOCK.COM/ROMAN ZAIETS
This article introduces a comprehensive and unified environment that enables the exploration of power systems across various modeling levels, encompassing load-flow analysis and extending to intricate time-domain simulations featuring inverter-based resources (IBRs).
Initializations of IBR circuits and parallel computations are presented for practical system problems. In the second part of this article, it is demonstrated that the time-domain accuracy level can be further increased through detailed and integrated modeling of magnetic circuits in the proposed unified environment.
The Electromagnetic Transients Program (EMTP) is designed for the simulation of electromagnetic transients. Since its computational engine is circuit based, it allows simulations for phenomenon frequencies ranging from dc to 100 MHz and more. It is possible to expand into a large spectrum of power system analysis applications.
The main concept is the usage of an expandable system of equations (ESE) that accommodates various circuit models through generic relations. The ESE is based on modified augmented nodal analysis (MANA) created for the first time for the development of EMTP code. It allows formulating a generic system of equations stated as \[{\bf{Ax}} = {\bf{b}} \tag{1} \] where the matrix ${\bf{A}}$ is the matrix of coefficients, ${\bf{x}}$ is the vector of unknowns, and ${\bf{b}}$ is the vector of known values. The matrix ${\bf{A}}$ consists of the classic nodal admittance matrix but augments it with various other model equations, such as ideal transformers, ideal or nonideal switches, independent source functions, dependent source functions, and any other type of equation derived for a model.
The time-domain solution of (1) is based on converting differential equations of models into the well-known companion circuit representation through numerical integration. For power electronics converters, the ESE allows the modeling of switching devices, such as insulated-gate bipolar transistors (IGBTs) or diodes, through ideal switch models or detailed nonlinear functions. Conversely, controlled sources can be employed to furnish the converters with average value models, depending on the modulation index. In the time-domain, the augmented part of ${\bf{A}}$ also contains the relations for nonlinear functions, representing, for example, a nonlinear inductance function for the transformer magnetization branch or IGBT nonlinear characteristic. The matrix ${\bf{A}}$ is time variant and must be iteratively updated to account for switch position changes or changes in nonlinear function slopes. In fact, ${\bf{A}}$ is the Jacobian matrix used in Newton’s method.
The ESE formulation method is used, furthermore, for performing multiphase unbalanced load-flow simulations with phasors. Such simulations allow the study of operating conditions in transmission, distribution, and microgrid networks. The ESE for load-flow can easily accommodate traditional load-flow constraints and integrate more recent constraint (such as droop functions, for example) requirements for power electronics converters. In the load-flow solution, (1) is a complex equation, with the augmented part of ${\bf{A}}$ now containing load-flow constraint equations. Once again, ${\bf{A}}$ becomes the Jacobian matrix and must be updated iteratively until convergence for a given tolerance. It has been shown that this formulation is more efficient and converges much faster than traditional mismatch-based load-flow solvers and various other methods available in the literature. Its advantages are remarkable for ac/dc systems, including microgrids in connected or islanded modes of operation. It can accommodate multifrequency conditions in islanded mode. This formulation accommodates arbitrary balanced or unbalanced multiwire networks. The time-domain models are seamlessly compatible with load-flow equations, facilitating a direct transition into the steady-state solution that serves as the initialization for the time-domain models. Contrary to the load-flow solution, the steady-state solution of (1) uses lumped circuits for all components based on initialization from load-flow solution phasors.
The multiphase load-flow solver is used to initialize IBRs for achieving a quick steady state for time-domain waveforms from which various analyses can be performed. This is an essential aspect not only for dramatically reducing computing times but also for finding the correct and stable steady-state solution. A unified environment is realized for modeling IBRs in the steady state and time-domain, with all related control systems. The computational burden of time-domain simulations with IBRs is alleviated through parallel processing methods.
The simulation of very large grids with IBRs is a challenging problem. At first, it is necessary to use advanced graphical interfaces to present system schematics and enter data for components. A schematic editor must be used to simplify the tasks of the study engineer. Hierarchical arrangements are necessary to achieve high-level visualization. Detailed component model information is encapsulated within multilevel subcircuits, providing a structured approach for data entry through masks.
A sample studied system is presented in Figure 1. The grid-forming converters of this system are initialized from the shown load-flow solution. IBRs can effectively benefit from load-flow PV or PQ constraints. It is possible to also include droop control. Such a system includes all components, with detailed control systems for conducting performance studies with traditional generation sources and modern power electronics components, such as grid-forming converters. Scenario variation tools allow comparing the performances of synchronous generators replaced by grid-forming converters. Sample simulation results are shown in Figure 2 for a three-phase-to-ground fault (see the fault location in Figure 1) appearing at 1 s and disappearing at 1.0833 s. It is apparent that the grid-forming inverter performs better than the synchronous generator for the same fault condition.
Figure 1. A sample test system integrating grid-forming converters with conventional generation.
Figure 2. A performance analysis with synchronous generators replaced by grid-forming inverters. (a) Grid-forming inverter 3, with per-unit active power, replacing synchronous generator 3. (b) Synchronous generator 3, with per-unit active power.
Enhanced unification is achieved by including all the protection relays in Figure 1. Such complete system representation is in line with the digital twin concept since all necessary components are included to create an electronic copy of the studied grid.
A much larger grid model is presented in Figure 3. It contains 10 wind parks (WPs) with aggregated models. In such a system, the WPs can be simulated using different numerical integration time-steps by decoupling them through propagation delay-based transmission line models. This approach is sufficiently accurate due to the interpolation process in the history buffer of transmission lines or cables. This concept is illustrated in Figure 4, where the network is simulated using a 50-µs time-step, with average value models for all IBRs, and only one of the WPs is simulated using a time-step of 10 µs, with detailed modeling of its converter (IGBT nonlinear models). A phase-A-to-phase-B-to-ground fault (see the fault location in Figure 3) is occurring at 2 s, with unsuccessful fault clearing followed by line tripping at 2.25 s. The results are compared to 10-µs usage for the entire system to validate accuracy. Such a setup allows reducing by five times the computing time on two CPUs.
Figure 3. A transmission grid with 10 wind parks.
Figure 4. The active power of a WP, using two-time-step simulation.
Initialization can be further improved through more detailed initialization of IBR control blocks. The basic top-level procedure is given in Figure 5 for a WP. It is applicable also to other types of IBRs. In this version, the ideal voltage source phasor is found from the load-flow solution. It is used in the steady-state solution and then transposed into the time-domain. It stays connected for a duration of ${T}_{i}$ (typically less than 200 ms). This duration allows the IBR to initialize rapidly and reach quasi-steady-state conditions when the source switch is disconnected to interface the IBR directly with its network. The IBR contents must be also initialized automatically, including the initialization of control system blocks and the collector grid.
Figure 5. The WP initialization procedure.
As demonstrated in Figure 6, all the IBRs of Figure 3 achieve steady-state conditions before 200 ms. A phase-A-to-ground fault (see the location in Figure 3) occurs at 1 s, followed by line opening for clearing the fault and then reclosing to recover the system and return to a stable steady state. The entire system is simulated with average value models for WP converters and a time-step of 50 µs.
Figure 6. The active powers of WPs in Figure 3.
The needs for modeling and simulating large-scale power systems in EMT mode are constantly increasing. Recently, the Chilean system operator (Coordinador Eléctrico Nacional) performed EMT simulations for the Chilean power grid. This work was part of a road map to enable the decommissioning of fossil fuel-based power plants and reaching a 100% share of renewable energies as of 2030. This accelerated energy transition scenario must meet the enabling conditions to prepare the grid to integrate new technologies and execute the necessary investments in renewable generation and storage. Network reliability being a key issue, the integration of IBRs requires detailed modeling and analysis for the large-scale Chilean grid. The current EMT (EMTP) model contents are detailed as follows:
The above IBR models are generic in preparation for studies of the eventual integration of manufacturer models. The generic models are adequate for testing various performance scenarios and operating conditions.
The computing times for such a system can reach 120 s for a 1-s simulation interval, for a numerical integration time-step of 50 µs when using a single CPU on a laptop computer. This speed can be improved significantly by using several CPUs in parallel. With 60 CPUs, the computing time for 1 s drops to 13 s. A sample simulation waveform for a fault applied at 1 s is shown in Figure 7 to validate the parallel solution.
Figure 7. The active power at one WP at the connection point.
More recently, the ESE approach has been used to integrate the capability to model in detail magnetic circuits. It is presented in the following section.
Computational electromagnetics aims to develop high-performance numerical models and provides engineers with tools to conceive devices that meet technical requirements and to analyze them in various scenarios. Simulation allows testing different situations that may be challenging to achieve experimentally due to safety concerns, high costs, lack of equipment, or time constraints. Simulations can replicate extreme cases, transient states, and faulty conditions, which deepens the understanding of complex situations and helps prevent or mitigate them.
Historically, the evolution of simulation tools has been closely intertwined with the progress of computer computational capabilities. It has transitioned from basic analytical models and expert-derived formulas to very complex numerical models and, further, to multiphysics simulations. Presently, breakthroughs in cloud computing, edge computing, and artificial intelligence are igniting ambitions toward real-world simulation. This involves modeling the subject device while considering its ultimate environment. Evidence of this current trend is the rising popularity of the concept of digital twins in both academic and industrial spheres. This concept serves to shape the next generation of power grids, characterized by their minimized environmental footprint but also well suited to withstand the evolving impacts of climate changes.
In numerous instances, the device of interest operates within a large system with interactions among multiple components. Thus, it is crucial to adopt a system approach that considers the device’s environment. For instance, when studying faults affecting transformers, it is insufficient to analyze only their impacts on the transformers themselves. Effects on other power grid equipment must also be evaluated. Similarly, transients occurring in the power grid, such as lightning, ferroresonance, and geomagnetically induced currents (GICs), can significantly affect transformers and should be considered during design and operation phases.
While detailed simulators that use classical numerical methods, such as finite elements, are powerful modeling tools, they often consume considerable computing time and focus on devices individually, rendering them incompatible with complete simulations. Several issues are encountered, including very long delays to reach the steady state.
A possible solution is to use a circuit-based simulation tool, such as EMTP. EMTP is a comprehensive simulation software package. Its extensive library empowers users to generate digital replicas of the studied systems. The ESE-based kernel of EMTP incorporates the ability to model the magnetic circuits of devices, such as transformers and rotating electrical machines, using lumped circuit-based methods.
Circuit-based modeling methods offer several distinct advantages, the most important one being seamless integration within complex power systems networks. This integration allows harnessing the benefits of a system approach, enabling a more comprehensive analysis of devices within large systems. One noteworthy advantage is the combination of a magnetic circuit meshing approach with circuit methods, based on a device’s geometry. This fusion is highly valuable for dimensioning and optimization studies, facilitating precise analysis and fine-tuning of designs. Moreover, the mesh approach ensures a generic aspect of the modeling process, leading to the automatic generation of models. This automation streamlines the modeling procedure and enhances efficiency.
Furthermore, it is essential to highlight that the precision achieved by circuit-based methods is comparable to that of classical numerical methods, particularly when the meshing resolution is equivalent. This means that by appropriately adjusting the meshing parameters, circuit-based models can yield results on par with traditional methods, ensuring accuracy and reliability in the simulation process. This attribute further solidifies the value of circuit-based modeling as a viable alternative capable of providing fast, reliable, and precise results while benefiting from seamless integration with circuit-based software.
Diverse circuit-based methods for magnetic circuit modeling have been proposed in the literature, among which the partial element equivalent circuit method and the transmission line matrix method are prominent examples. This presentation adopts the equivalent circuit-based modeling approach for magnetic circuits, including all windings and interfacing with surrounding electric circuits. This choice is reinforced by compelling reasons, including robustness, simplicity of implementation, and excellent suitability for low- and midfrequency applications, which aligns perfectly with the characteristics of most electrical network components.
The following presents an application focused on transformers. Transformers exhibit nonlinear characteristics and operate at diverse power grid levels, subject to various transients. Three circuit-based methods can be selected:
The ESE of EMTP accommodates all dependent sources needed for the integration of reluctance network models. Elements such as current-controlled voltage sources (CCVSs), current derivative-controlled voltage sources (CDCVS), current controlled current sources (CCCS), and voltage-controlled voltage sources (VCVSs) can be seamlessly integrated within this framework. These elements, when combined, effectively operate as mutators, gyrators, or ideal transformers. This inherent versatility ensures a simultaneous and strong coupling between magnetic and electrical models. The ESE framework also offers flexibility by including nonlinear components solved iteratively for the highest accuracy and stability. This adaptability aids in optimizing convergence, particularly when dealing with huge numbers of nonlinear elements.
For all three analogies above, the magnetic circuit model is generated through an automated procedure. Initially, the transformer topology undergoes discretization, effectively partitioning the entire domain into a mesh structure, as exemplified in Figure 8. A three-phase and three-column transformer is modeled in this figure. Each area (domain) within the transformer can have a distinct mesh size, contributing to the establishment of an adaptive and tailored mesh, as in traditional finite-element modeling (FEM). Areas are determined based on two physical considerations. First, there is the inherent nature of the materials involved (ferromagnetic, dielectric, conductive, and so on). Second, the presence or absence of magnetomotive forces further influences the demarcation of these regions.
Figure 8. The transformer geometry mesh and examples of cells used to represent an element of the mesh for each of the three circuit-based methods.
Subsequently, to formulate the circuit-based equivalent model corresponding to the meshed system, elementary circuit “cells” are systematically associated with each individual mesh element; the code dynamically selects the appropriate type of cell, depending on the diverse characteristics exhibited by each element’s area. A novel library has been developed, consisting of an exhaustive array of cells to encompass all conceivable configurations that may arise in electromagnetic device modeling. The parameters of these cells are derived from the geometric and physical attributes of the individual mesh element. By effectively interconnecting these cells, a cohesive and comprehensive equivalent circuit model is established in the EMTP graphical user interface (see Figure 9).
Figure 9. A detailed magnetic circuit model of a three-phase transformer generated automatically in EMTP.
In essence, this automated approach allows for the systematic generation of models encapsulated as subcircuit devices. Upon double-clicking the device, users are granted access to a detailed representation of the model, faithfully aligned with the employed mesh configuration, as demonstrated in Figure 9.
This hierarchical structure closely mirrors the topological design of the studied device. The users can access internal magnetic circuit details in addition to surrounding complex power grid models.
The results stemming from the three circuit-based models, generated within EMTP, are juxtaposed against FEM outcomes acquired via ANSYS Electromagnetics. These evaluations pertain to a three-phase transformer (15 MVA, 25.63 to 4.53 kV). In Figure 10, the voltage current characteristics extracted from both open-circuit and short-circuit tests are presented. In the open-circuit test, the secondary is unloaded, and the primary winding is excited across varying voltage levels. The corresponding primary winding current is determined for each voltage level. Subsequently, in the short-circuit test, the secondary winding is short-circuited, and the primary winding is excited using different voltage levels. The secondary winding current is then deduced for each voltage condition.
Figure 10. The transformer V-I characteristic for (a) open-circuit and (b) short-circuit tests.
The observed maximum discrepancies with the FEM approach amount to 1.5% for the open-circuit test and 3.6% for the short-circuit. In this simulation, the circuit-based model utilizes a mesh including 2,848 elements, with the inclusion of over 2,200 nonlinear components. In contrast, the FEM mesh involves 3,292 elements. The simulations were conducted over a 32-ms interval, employing a time-step of 1 ${\mu}$s. The circuit-based model demonstrates a speed advantage, operating approximately three times faster than FEM, particularly under conditions characterized by significant nonlinearity and saturation effects. The good alignment between the outcomes of the circuit-based model and FEM results underscores the accuracy and efficiency of the proposed approach in reproducing electromagnetic behavior. One more important advantage is the capability to initialize systems dramatically faster than in the FEM approach.
Another validation scenario encompassing an electromagnetic transient is the study of inrush current arising from transformer energization. This investigation is presented in Figure 11. Once more, the maximum relative error observed between the circuit-based method and FEM remains below 3%.
Figure 11. The transformer energization transient.
This modeling approach not only facilitates the determination of global quantities but also possesses the capacity to capture localized behaviors. For instance, it can effectively provide the induction level at specific points within the transformer over time. Furthermore, the mesh data can be exploited with circuit solutions to generate field maps and animations.
To illustrate this capability, an internal fault scenario occurring in phase A of the primary winding is simulated. It represents a short circuit between internal winding turns appearing at 50 ms. The global quantities in Figure 12 enable studying the transient stemming from this fault. Complementing this, Figure 13 provides field maps that depict the distribution of induction levels at the instances of 52.5 and 60 ms. These maps offer a visual representation that readily reveals potential hot spots arising due to system malfunction.
Figure 12. The transient due to a short circuit between primary winding turns that occurs at 50 ms: the (a) primary current for the three phases and (b) active and reactive powers.
Figure 13. The distribution of induction levels, in teslas, at (a) 52.5 and (b) 60 ms for a fault case.
An additional advantage of these models lies in their ability to harness the EMTP environment for the purpose of conducting study scenarios necessitating a systemwide approach. Notably, GICs provoked by solar storms exemplify such scenarios. This phenomenon spreads over vast distances. Consequently, simulating such scenarios using conventional numerical methods becomes either prohibitively costly in computing times or practically unattainable.
The proposed circuit-based modeling approach offers a tangible solution. It not only facilitates the simulation of such complex scenarios within reasonable computing times but also provides access to internal constraints. This holds considerable value, not only for designers seeking to fortify transformer resilience but also for system operators tasked with implementing safeguards and protections to preemptively mitigate the impact of such disturbances.
Figure 14 illustrates a simplified configuration adopted to emulate a GIC scenario. The unloaded transformer is connected to the network. The GIC disturbance is injected in the system through the transformer neutral. The transformer is tested with two connectivity modes: Y-ground primary and Y-grounded secondary and Y-grounded primary and delta secondary. Figure 15 clearly demonstrates the significant differences in GIC impacts on the Y and delta configurations. This is explained by the fact that the delta connection does not allow the third harmonic. The shown reactive power oscillation exhibits slower stabilization compared to actual behavior, mainly due to the absence of iron losses in the current circuit-based model. These losses could be added as a parallel resistor to each winding for better accuracy.
Figure 14. The configuration for studying the GIC effect on a transformer.
Figure 15. The disturbance due to a GIC in the case of Yy0 and Yd11 connections: the (a) primary current and (b) reactive power.
The circuit-based model enables estimating iron losses through postprocessing. These estimations encompass both a comprehensive global value and a meticulous evaluation at the granularity of individual elements within the mesh. This process discerns the contributions of hysteresis losses and eddy current losses. This capability is demonstrated in Figure 16, where the aggregate losses over a single period, at a voltage level of 50% of the nominal supply, are estimated at 11.1 kW. Within this estimation, 23% of the losses are attributed to eddy current losses, while the remaining 77% arise from hysteresis losses. Moreover, the comprehensive distribution of these losses is pictured through the colormap.
Figure 16. The estimation of the distribution of iron losses, in watts, for a healthy case.
The presented unified environment is circuit based and allows simulating power systems, from load-flow solutions to fast electromagnetic transients. The load-flow solution plays a pivotal role in initializing simulated networks, enabling rapid attainment of steady-state conditions, even in the presence of IBRs, such as WPs.
Furthermore, the presented environment integrates the capability to model in detail magnetic circuits, with the noteworthy advantage of combining magnetic circuit meshing based on a device’s geometry.
A. Haddadi and J. Mahseredjian, “Power system test cases for EMT-type simulation studies,” CIGRE Brochure, Paris, France, Rep. WG C4.503, Aug. 2018.
N. Rashidirad, J. Mahseredjian, I. Kocar, U. Karaagac, and O. Saad, “MANA-based load-flow solution for islanded AC microgrids,” IEEE Trans. Smart Grid, vol. 14, no. 2, pp. 889–898, Mar. 2023, doi: 10.1109/TSG.2022.3199762.
J. Mahseredjian, U. Karaagac, S. Dennetiere, and H. Saad, “Simulation of electromagnetic transients with EMTP-RV,” Numer. Anal. Power Syst. Transients Dyn. Inst. Eng. Technol., vol. 78, pp. 103–134, Jan. 2015.
U. Karaagac, J. Mahseredjian, S. Jensen, R. Gagnon, M. Fecteau, and I. Kocar, “Safe operation of DFIG-based wind parks in series-compensated systems,” IEEE Trans. Power Del., vol. 33, no. 2, pp. 709–718, Apr. 2018, doi: 10.1109/TPWRD.2017.2689792.
M. Ouafi, J. Mahseredjian, J. Peralta, H. Gras, S. Dennetière, and B. Bruned, “Parallelization of EMT simulations for integration of inverter-based resources,” Electric Power Syst. Res., vol. 223, Oct. 2023, Art. no. 109641, doi: 10.1016/j.epsr.2023.109641.
W. Nzale, J. Mahseredjian, X. Fu, I. Kocar, and C. Dufour, “Improving numerical accuracy in time-domain simulation for power electronics circuits,” IEEE Open Access J. Power Energy, vol. 8, pp. 157–165, Apr. 2021, doi: 10.1109/OAJPE.2021.3072369.
H. Xue, A. Ametani, and J. Mahseredjian, “Very fast transients in a 500 kV Gas-insulated substation,” IEEE Trans. Power Del., vol. 34, no. 2, pp. 627–637, Apr. 2019, doi: 10.1109/TPWRD.2018.2874331.
M. Naïdjate, N. Bracikowski, M. Hecquet, M. Fratila, M. M. Duro, and J. P. Ducreux, “An intelligent reluctance network model for the study of large power and distribution transformers,” in Proc. 6th Int. Adv. Res. Workshop Transformers (ARWtr), Cordoba, Spain, 2019, pp. 89–92, doi: 10.23919/ARWtr.2019.8930174.
S. R. Pordanjani, M. Naïdjate, N. Bracikowski, M. Fratila, J. Mahseredjian, and A. Rezaei-Zare, “Electromagnetic modeling of transformers in EMT-type software by a circuit-based method,” IEEE Trans. Power Del., vol. 37, no. 6, pp. 5402–5413, Dec. 2022, doi: 10.1109/TPWRD.2022.3177137.
S. R. Pordanjani, J. Mahseredjian, M. Naïdjate, N. Bracikowski, M. Fratila, and A. Rezaei-Zare, “Electromagnetic modeling of inductors in EMT-type software by three circuit-based methods,” Electric Power Syst. Res., vol. 211, Oct. 2022, Art. no. 108304, doi: 10.1016/j.epsr.2022.108304.
Jean Mahseredjian (jeanm@polymtl.ca) is with Polytechnique Montréal, Montréal, QC H3T 1J4, Canada.
Mohammed Naidjate (mohammed.naidjate@polymtl.ca) is with Polytechnique Montréal, Montréal, QC H3T 1J4, Canada.
Mehdi Ouafi (mehdi.ouafi@rte-france.com) is with RTE, 69330 Lyon, France.
Juan Antonio Ocampo Wilches (juan-antonio.ocampo-wilches@polymtl.ca) is with Polytechnique Montréal, Montréal, QC H3T 1J4, Canada.
Digital Object Identifier 10.1109/MELE.2023.3320511
2325-5897/23©2023IEEE