Applied Science, Faculty of
Electrical and Computer Engineering, Department of
DSpace
UBCV
Ramesh, Rahul Raman
Master of Applied Science - MASc
University of British Columbia
Time-domain simulations have become the preferred approach in modern terrestrial and shipboard power systems due to the increasing integration of high-frequency power electronic devices. These developments necessitate smaller transient analysis windows, further driving innovation in power systems simulation tools and computational technology. First introduced and pioneered by Prof. Herman Dommel, the nodal-analysis-based electromagnetic transient programs (EMT or EMTP-Type) lead the field. Despite the existing simulation programs based on EMTP-type solution algorithms, several critical areas remain active in the research and development of modern and future offline and real-time EMT simulators. For instance, synchronous machines (SMs) models present significant computational challenges due to their complex dynamics and extensive calculations, especially when coupled with switching power converters. Recent advancements have focused on improving computational efficiency in SM simulations. These efforts involve reformulating the SM model using constant parameter matrices, significantly reducing the computational burden. This thesis builds upon this foundation by advancing the constant parameter voltage-behind-reactance (CP-VBR) saturable machine model to support an arbitrary number of damper windings while implemented using a multirate (MR) algorithm. The proposed MR CP-VBR SM model separates fast and slow dynamics, enabling further reduction of computational costs by lowering the sampling rate for slower subsystems. Through extensive computer studies, this research demonstrates that the proposed MR CP-VBR SM model offers superior performance compared to existing models across diverse scenarios. Including magnetic saturation and an arbitrary number of damper windings make the proposed model particularly suitable for modern shipboard power systems where SM are significant components and simulation accuracy and speed are critical.
MULTIRATE CONSTANT PARAMETER GENERALIZED SYNCHRONOUS MACHINE MODEL FOR SIMULATION OF ELECTROMAGNETIC TRANSIENTS IN SHIPBOARD POWER SYSTEMS by Rahul Raman Ramesh B.Tech., SRM Institute of Science and Technology, 2022 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April, 2025 ©Rahul Raman Ramesh, 2025 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled: MULTIRATE CONSTANT PARAMETER GENERALIZED SYNCHRONOUS MACHINE MODEL FOR SIMULATION OF ELECTROMAGNETIC TRANSIENTS IN SHIPBOARD POWER SYSTEMS submitted by Rahul Raman Ramesh in partial fulfilment of the requirements for the degree of Master of Applied Science in Electrical and Computer Engineering Examining Committee: Dr. Juri Jatskevich, Professor, Electrical and Computer Engineering, UBC Supervisor Dr. José Martí, Professor, Electrical and Computer Engineering, UBC Supervisory Committee Member Dr. Christine Chen, Associate Professor, Electrical and Computer Engineering, UBC Supervisory Committee Member iii Abstract Time-domain simulations have become the preferred approach in modern terrestrial and shipboard power systems due to the increasing integration of high-frequency power electronic devices. These developments necessitate smaller transient analysis windows, further driving innovation in power systems simulation tools and computational technology. First introduced and pioneered by Prof. Herman Dommel, the nodal-analysis-based electromagnetic transient programs (EMT or EMTP-Type) lead the field. Despite the existing simulation programs based on EMTP-type solution algorithms, several critical areas remain active in the research and development of modern and future offline and real-time EMT simulators. For instance, synchronous machines (SMs) models present significant computational challenges due to their complex dynamics and extensive calculations, especially when coupled with switching power converters. Recent advancements have focused on improving computational efficiency in SM simulations. These efforts involve reformulating the SM model using constant parameter matrices, significantly reducing the computational burden. This thesis builds upon this foundation by advancing the constant parameter voltage-behind-reactance (CP-VBR) saturable machine model to support an arbitrary number of damper windings while implemented using a multirate (MR) algorithm. The proposed MR CP-VBR SM model separates fast and slow dynamics, enabling further reduction of computational costs by lowering the sampling rate for slower subsystems. Through extensive computer studies, this research demonstrates that the proposed MR CP-VBR SM model offers superior performance compared to existing models across diverse scenarios. Including magnetic saturation and an arbitrary number of damper windings make the proposed model particularly suitable for modern shipboard power systems where SM are significant components and simulation accuracy and speed are critical. iv Lay Summary Synchronous machine models are essential when conducting studies of large power networks, integrating renewable technologies into the grid, and designing shipboard power systems, where electrical machinery plays a significant role. Traditional models of synchronous machines are computationally expensive, making it hard to conduct simulations of practical systems. This thesis extends the previous work and presents new models of synchronous machines that consider magnetic saturation and MR to achieve better computational performance. The proposed models are extensively validated in the nodal-analysis-based electromagnetic transient (EMT) program PSCAD in multiple scenarios, demonstrating superior computational benefits and accuracy. The new models are envisioned to benefit many practitioner engineers working on studies with renewable energy, power grids, and shipboard power systems. v Preface The research outcomes presented in this thesis have been partly published and/or are under preparation for future publications. Unless otherwise specified, I was directly responsible for deriving mathematical equations and models, implementing these models into computer programs, analyzing the results, and writing and editing the first drafts of all papers published from this thesis. My supervisor, Dr. Juri Jatskevich, provided crucial guidance and leadership throughout all stages of this work. He also provided constructive feedback and edits on the writing and editing of all papers. Additionally, my co-authors contributed significantly to this research by providing critical input and revising manuscripts for clarity and accuracy. The list of articles from this research and further contributions from co-authors is provided below. Chapter 2 is based on the following conference paper: R. R. Ramesh, T. Vahabzadeh, and J. Jatskevich, âConstant-Parameter VBR Synchronous Machine Model for Studies of Unbalanced Faults in PSCAD,â IEEE 23rd International Symposium, (INFOTEH-JAHORINA, 2024), Jahorina, 20 - 22 Mar. 2024, 6 pages. Chapters 3 and 4 combine the models and the studies presented in the following papers and ongoing work: R. R. Ramesh, E. Mostajeran, S. Ebrahimi, and J. Jatskevich, âMultirate VBR Model of Synchronous Machines for EMTP-Type Simulation Programs,â IEEE International Conference on Modelling, Simulation and Intelligent Computing (MoSICom 2023), Dubai, UAE, 7 - 9 Dec. 2023, 6 pages. doi: 10.1109/MoSICom59118.2023.10458759. vi R. R. Ramesh, A. Safavizadeh, S. Ebrahimi, and J. Jatskevich, âEfficient Modeling of Synchronous Machines for Transient Simulation of Marine Electrical Power Systems,â IEEE International Conference on Green Industrial Electronics and Sustainable Technologies (GIEST-2024), NIT Manipur, India, 25 - 26 Oct. 2024, 6 pages. Best Paper Award. Generative Aritificical Intelligence (GEN AI) was not used in the process of drafting and preparation of this thesis or any of the research work presented in this thesis. vii Table of Contents Abstract ............................................................................................................................. iii Lay Summary ................................................................................................................... iv Preface .................................................................................................................................v Table of Contents ............................................................................................................ vii List of Figures .....................................................................................................................x List of Symbols and Notations ...................................................................................... xiii List of Abbreviations .................................................................................................... xvii Acknowledgements ...................................................................................................... xviii Dedication .........................................................................................................................xx Chapter 1: Introduction ....................................................................................................1 1.1 Motivation ........................................................................................................... 1 1.2 Literature Review................................................................................................ 1 1.3 Research Objectives ............................................................................................ 2 1.3.1 Objective I:...................................................................................................... 3 1.3.2 Objective II: .................................................................................................... 3 Chapter 2: Interfacing of SR CP-VBR Synchronous Machine in EMTP-Type Solutions ..............................................................................................................................4 2.1 State Space VBR Model ..................................................................................... 5 2.2 Main Flux Saturation in EMTP-Type Programs ................................................. 7 2.3 SM EMTP-Type VBR Model ............................................................................. 8 2.4 EMTP-Type CP-VBR Model ........................................................................... 10 viii 2.5 Interfacing of the Machine Model .................................................................... 14 2.6 Computer Study ................................................................................................ 17 Chapter 3: MR CP-VBR Machine models ....................................................................21 3.1 EMTP-Type Solver Implementation ................................................................. 21 3.2 MR Algorithm ................................................................................................... 22 3.2.1 MR counter ................................................................................................... 24 3.2.2 Interfacing of multistep algorithms ............................................................... 24 3.3 MR-based CP-VBR SM Models ....................................................................... 25 3.3.1 MR-TT CPVBR SM model .......................................................................... 26 3.3.1.1 Slow subsystem ...................................................................................... 26 3.3.1.2 Fast subsystem........................................................................................ 27 3.3.2 The MR-TBE CP-VBR SM model ............................................................... 29 3.3.2.1 Slow subsystem ...................................................................................... 29 3.3.2.2 Fast subsystem........................................................................................ 30 3.4 MR CP-VBR Implementation........................................................................... 31 Chapter 4: MR Model Studies ........................................................................................35 4.1 Preliminary SMIB Study................................................................................... 35 4.2 PSCAD/EMTDC Implementation of a Microgrid ............................................ 39 4.2.1 Time-domain simulations ............................................................................. 41 4.2.2 Numerical accuracy ...................................................................................... 43 4.2.3 Computational efficiency .............................................................................. 44 4.3 Extended PSCAD/EMTDC Marine Propulsion System ................................... 44 4.3.1 MR numerical accuracy ................................................................................ 49 ix 4.3.2 Computational efficiency .............................................................................. 51 Chapter 5: Conclusion .....................................................................................................53 5.1 Contributions..................................................................................................... 53 5.2 Anticipated Impact ............................................................................................ 55 5.3 Future Work ...................................................................................................... 55 Bibliography .....................................................................................................................57 Appendices ........................................................................................................................59 Appendix A : Coefficients introduced in Chapter 1 ..................................................... 59 Appendix B : Machine and System Parameters for the study in section 1.6 [2] .......... 62 Appendix C : Additional coefficients introduced in Chapter 2 .................................... 63 Appendix D : Machine and System Parameters for the studies in Chapter 4 ............... 64 D.1 Hydro turbine generator used in section 4.1[2-Table 5.10-1] ....................... 64 D.2 Marine microgrid system parameters for section 4.2.................................... 64 D.3 Extended marine microgrid system parameters for section 4.3 .................... 64 x List of Figures Figure 2.1 The generator connection circuit for fault detection using neutral to ground current [19].......................................................................................................................................................... 5 Figure 2.2 Piecewise linear saturation curve of d-axis saturation. ................................................. 7 Figure 2.3 Phasor representation of main-flux saturation, with projections along the qd-axes...... 8 Figure 2.4 Physical component to EMTP-Type component representation: (a) physical component, (b) branch-defined component, and (c) node-defined component. ............................................... 14 Figure 2.5 Expansion of the conductance matrix: (a) solidly grounded component, and (b) component with separate neutral node. ......................................................................................... 15 Figure 2.6 A single-machine system with generator grounding. .................................................. 17 Figure 2.7 A single-machine system implemented in PSCAD: (a) using the built-in qd0 model, and (b) using the proposed generalized CP-VBR model. ............................................................. 17 Figure 2.8 Transient response of the generator to a single-line-to-ground fault: (a) stator current, (b) neutral-to-ground current, and (c) electromagnetic torque. .................................................... 19 Figure 2.9 Magnified views from Figure 2.8: (a) stator current, (b) neutral-ground current, and (c) electromagnetic torque. ................................................................................................................. 19 Figure 2.10 Comparison of norm-2 error for the solution trajectory of stator phase current ias and neutral current ing for the PSCADâs built-in qd0 and the proposed CP-VBR models at different timestep sizes. ............................................................................................................................... 20 Figure 3.1 Simplified program flow in an EMTP-type simulator, considering a constant network conductance matrix. ...................................................................................................................... 22 Figure 3.2 MR solution algorithm implementation in an EMTP-type solver. .............................. 23 xi Figure 3.3 A typical MR simulation timeline depicting commonly used terms for an MR ratio of three............................................................................................................................................... 24 Figure 3.4 Slow transient prediction and accurate fast transient area computation depicted graphically..................................................................................................................................... 25 Figure 4.1 Phase-A current response during a three-phase fault: (a) overall transient, and (b) magnified view.............................................................................................................................. 37 Figure 4.2 Rotor speed response during a three-phase fault: (a) overall transient, and (b) magnified view. .............................................................................................................................................. 38 Figure 4.3 Cumulative 2-norm error for the stator current and rotor speed trajectories as obtained by the proposed MR CP-VBR model with different raios M. ...................................................... 38 Figure 4.4 CPU time for the considered transient study simulated using different MR ratios M. 38 Figure 4.5 Considered shipboard electrical system comprised multiple synchronous generators supplying AC and DC loads. ........................................................................................................ 40 Figure 4.6 Part of the shipâs electrical system depicted in Figure 4.5 implemented in PSCAD/EMTDC using the custom-defined MR synchronous machine model. .......................... 40 Figure 4.7 Transient response of generators G1 and G3 variables during the operation of circuit-breakers at the tie lines, opened at t=0.05s and re-closed at t=0.10s, for different values of hf and M: (a) phase-a stator current of G1, (b) phase-a stator current of G3, (c) electrical rotor speed of G1, and (d) electrical rotor speed of G3. ...................................................................................... 41 Figure 4.8 Zoomed in sections from Figure 4.7: (a) phase-a stator current of G1, (b) phase-a stator current of G3, (c) electrical rotor speed of G1, and (d) electrical rotor speed of G3. ..................... 42 Figure 4.9 Calculated 2-Norm error for: (a) phase-a stator current of G3, and (b) electrical rotor speed of G3, at various fast timesteps and MR ratios. .................................................................. 43 xii Figure 4.10 Normalized execution time for the custom model in Figure 4.6 for different values of MR ratio. ....................................................................................................................................... 44 Figure 4.11 Overview of the marine microgrid transient study. ................................................... 45 Figure 4.12 PSCAD implementation of the marine microgrid shown in Figure 4.11. ................. 46 Figure 4.13 Switched stator phase-a voltage of propulsion motor M1 from the marine micro-grid study at an SR: (a) overall transient, (b) pre-command zoomed-in transient, (c) post-command zoomed-in transient. ...................................................................................................................... 47 Figure 4.14 Switched stator phase-a currents of propulsion motor M1 from the marine micro-grid study at an SR: (a) overall transient, (b) pre-command zoomed-in transient, (c) post-command zoomed-in transient. ...................................................................................................................... 47 Figure 4.15 Electrical rotor speed of propulsion motor M1 from the marine microgrid study at a SR. ................................................................................................................................................. 48 Figure 4.16 Zoomed in sections from Figure 14 and Figure 15 with M=15 from propulsion motor 1: (a) phase-a current and (b)electrical rotor speed. ..................................................................... 50 Figure 4.17 Relative 2-norm percentage error of variables of propulsion motor 1: (a) phase-a current, (b) q-axis magnetizing flux-linkage, and (c) electrical rotor speed. ................................ 50 Figure 4.18 Normalized average per-step time across multiple MR ratios. ................................. 51 xiii List of Symbols and Notations The conventions and notations in this thesis are taken from the subject literature [5]â[11]. Specifically, matrices are denoted by bold upper-case letters and vectors by bold lower-case letters. The motor sign convention is used for the stator currents. The rotor reference frame is considered for abc ï® qd0 transformation. The matrices sK and 1sï­K represent the Parkâs and inverse Parkâs transformation, respectively [5]. The stator variables for compactness are grouped into vectors as ï ïTabcs as bs csf f fï½f , 0 0Tqd s qs ds sf f fï© ï¹ï½ ï« ï»f , where f may represent flux linkages, voltage or current. The transformation from abc ï® qd0 has the following form 0qd s s abcsï½K ff . The synchronous machine is considered to have m q-axis and n d-axis damper windings along with one field winding. Similar to the stator variables, the rotor variables are represented by vectors as 1Tqr kq kqmf fï© ï¹ï½ ï« ï»f ï , 1Tdr fd kd kdnf f fï© ï¹ï½ ï« ï»f ï , qrqdrdrï© ï¹ï½ ïª ïºï« ï»fff. Only the commonly used variables are provided in this section. Other variables are defined in this thesis on a need-to-know basis, as they appear in the text. xiv r cï´0 A zero matrix of size r cï´ r cï´1 A matrix of 1 of size r cï´ , ,as bs cse e eï¢ï¢ ï¢ï¢ ï¢ï¢ Sub-transient phase a, b, and c voltages ,qs dse eï¢ï¢ ï¢ï¢ q-axis and d-axis sub-transient voltages h SR discrete timestep fh MR fast transient discrete timestep sh MR slow transient discrete timestep , ,as bs csi i i Stator phase a, b, and c currents fdi Field winding current kdxi d-axis damper winding current ï» ï½1, ,x nï ï kqxi q-axis damper winding current ï» ï½1, ,x mï ï 0, ,qs ds si i i q-axis, d-axis, and zero-sequence stator currents rI Identity matrix of size r rï´ hJ History current injection vector lsL Stator leakage inductance lfdL Field winding leakage inductance lkdxL d-axis damper winding leakage inductance ï» ï½1, ,x nï ï xv lkqxL q-axis damper winding leakage inductance ï» ï½1, ,x mï ï , mdumquL L Unsaturated q-axis and d-axis magnetizing inductance , mdsmqsL L Saturated q-axis and d-axis magnetizing inductance M MR ratio p Heavisideâs operator P Number of poles sr Stator resistance kdxr d-axis damper winding resistance ï» ï½1, ,x nï ï kqxr q-axis damper winding resistance ï» ï½1, ,x mï ï SF Saliency factor eT Generated electromagnetic torque mT Applied mechanical torque t Time , ,as bs csv v v Stator phase a, b, and c voltages fdv Field winding voltage 0, ,qs ds sv v v q-axis, d-axis, and zero-sequence stator voltages kdxv d-axis damper winding voltages ï» ï½1, ,x nï ï kqxv q-axis damper winding voltages ï» ï½1, ,x mï ï xvi ï¢ Flux vector angle ï¨ ï©xï² 2-norm relative percentage error of x rï± Electrical rotor position fdï¬ Field winding flux linkage kdxï¬ d-axis damper winding flux linkages ï» ï½1, ,x nï ï kqxï¬ q-axis damper winding flux linkages ï» ï½1, ,x mï ï ,mq mdï¬ ï¬ q-axis and d-axis magnetizing flux linkages ï¢ï· Time derivative of flux vector angle rï· Electrical speed of the rotor z MR counter xvii List of Abbreviations AC Alternating current BE Backward Euler CP-PD Constant parameter phase domain (model) CPU Central processing unit CP Constant parameter (interfacing method) CP-VBR Constant parameter voltage-behind-reactance (model) d-axis Direct axis DC Direct current EMT Electromagnetic transients EMTP Electromagnetic transients program EMTP-Type Nodal-analysis-based programs based on the original EMTP MR Multirate (algorithm) PD Phase domain (model) q-axis Quadrature axis QD Synchronous machine defined only using the qd-axis variables (model) SM Synchronous machine SR Singlerate (algorithm) VBR Voltage-behind-reactance (model) VP Variable parameter (interfacing) xviii Acknowledgements My deepest gratitude goes to my beloved supervisor, Dr. Juri Jastkevich. He has been more than just a supervisor, a confidant, a guide, and a genuinely caring person. He is a brilliant toolsmith and an even better artisan when it comes to his work. I am genuinely grateful to have been under his supervision during my formative years and to have assisted in bringing his ideas into existence. This research was funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada under the Collaborative Research and Development (CRD) grant entitled âResearch and Development of Advanced Techniques for EMT Simulation of Modern Electrical Energy Systems,â led by Dr. Jatskevich, and NSERC Alliance grant entitled âIntegrated Design and Operation of Hybrid-electric Shipboard Systems,â led by Dr. Chen. I feel incredibly privileged and humbled to have had the opportunity to work in this vital area and contribute to the development of state-of-the-art models and simulation tools used by industries worldwide. I would also like to thank the examining committee members, Dr. José Martí and Dr. Christine Chen, for their time and effort in reviewing my thesis and the invaluable comments and encouragement they provided to me over these years. You made me feel that no question was ever too small, answering each one with patience and generosity. It has been a true pleasure learning from you! The UBCâs Electric Power and Energy Systems group gave me the time and space to grow at my own pace. I am especially grateful to my fellow colleagues and mentors, Mr. Arash Safavizadeh and Dr. Milad Ebrihami, who taught me the importance of a strong work ethic and often helped guide me through my moments of disorganization. To my friends from the lab xix Ekamjot, Ziliang (Bruce), Abhay, Imran, and Shadman, thank you for making life enjoyable and memorable. Hereâs to many more years of friendship ahead! I also want to thank Erfan Mostejeran and Taleb Vahabzadeh for helping me build my first model and proofreading my work. A heartfelt thanks go to the Department of Electrical and Computer Engineering, staff from the administrative team and IT, and the custodial staff whose constant support and kindness created an incredibly welcoming environment that made completing this thesis possible. The presence of a few people in my life was crucial to my journey these past three years, Amritpal Singh Dhir and Muskan Puri. I cannot thank you enough for standing by me during some of my most challenging times. You held me up when I was slipping and reminded me of my goals when I had lost sight of them. To my brother, Rajesh Sankar Ramesh, thank you for being my anchor. I have been able to pursue my work peacefully because of the sense of security you gave me. Competing with your immense knowledge sparked my own desire to keep learning and growing, a competition that still drives me today. Finally, I want to thank my parents, Rajani Ramesh and Ramesh Sankaran. None of this would have been possible without your sacrifices. Your financial and emotional support when I embarked on this journey made everything possible. You gave up so much for my future; I will forever be grateful for that. They say it takes a village to raise a child, and it has taken many peopleâs effort, belief, and love to bring me to where I am today. xx Dedication To Rudvi and her Future 1 Chapter 1: Introduction 1.1 Motivation With the evolution and modernization of power systems, driven significantly by advancements in power electronics and renewable integration, there has been an increased demand for detailed and computationally efficient simulation models of synchronous machines (SM). Accurate transient analysis through electromagnetic transient programs (EMTP) is critical for designing, planning, and operating complex electrical systems such as utility grids and marine electric propulsion systems. Current state-of-the-art nodal-analysis-based EMTP simulators (EMTP-type) are extensively used due to their accuracy in capturing transient phenomena [1]. However, synchronous machine models traditionally implemented in these programs, primarily based on conventional qd0 models, are known to have numerical and computational inefficiencies when used at larger simulation steps. These inefficiencies become particularly significant in systems with multiple synchronous machines and power electronic devices, requiring extensive computational resources for transient studies, especially when sampled at smaller time steps. 1.2 Literature Review Historically, electromagnetic transient simulations have relied heavily on the qd0 models [2] for the SM, despite their interface and computational overhead limitations. As computational demands grew, researchers introduced alternative models such as the Phase-Domain (PD) [2] and Voltage-Behind-Reactance (VBR) formulations [3], known to provide better numerical stability and efficiency for larger time steps. Several commercial EMT-type simulators have become standards in both academia and industry due to their robustness, computational efficiency, and accuracy. These commercial simulators can be broadly categorized based on their underlying numerical methods: nodal-analysis-based (EMTP-type) and state-variable-based (SVB). State-variable-based (SVB) simulators such as MATLAB/Simulink (Simscape Electrical) [4], PLECS [5], Opal-RT [6], and Typhoon HIL [7] utilize direct numerical integration of system equations. These simulators provide flexibility in handling nonlinearities, power electronics, and detailed 2 control system modelling but face significant computational inefficiencies for large-scale power network simulations due to the dense matrices and smaller time step requirements. Unlike state-variable-based (SVB) simulators, electromagnetic transient program (EMTP)-type nodal-analysis simulators such as EMTP-RV [8], PSCAD/EMTDC [12], ATP/EMTP [9], and RTDS/RSCAD [10] solve network equations using nodal admittance matrices with trapezoidal numerical integration [11]. These simulators provide significant computational advantages in large-scale power network simulations due to their efficient handling of sparse matrices and larger permissible time steps, especially with the recent improvement with hardware accelerated solutions[12]. However, they typically offer less flexibility in directly modelling nonlinearities, advanced control systems, and detailed switching-level power electronic devices. They often require specialized interfacing techniques or customized models to integrate these elements accurately. The VBR model, initially introduced for state-variable-based simulation platforms, has been effectively adapted to nodal-analysis-based EMTP programs [13] due to advantageous interfacing capabilities. Subsequently, constant-parameter VBR (CP-VBR) [14] models were developed, significantly enhancing computational efficiency by providing a constant interface conductance matrix. Moreover, the concept of multirate (MR) simulations has emerged [15], [16], [17], exploiting the disparity in dynamic speeds between various subsystems. Initial MR implementations in machines targeted induction machines [18] showcasing promising computational gains, and this research is expected to expand upon this literature by combining constant interfacing and multirate algorithms. 1.3 Research Objectives This thesis aims to address the gaps in existing SM modelling methodologies for EMTP-type simulators, focusing on improving transient studiesâ computational efficiency and accuracy. Specifically, this research targets the following objectives: 3 1.3.1 Objective I: Development and implementation of a generalized saturable Constant-Parameter VBR (CP-VBR) synchronous machine model for EMTP-type simulators: ï· Include support for arbitrary numbers of damper windings and incorporate the floating neutral point to accurately represent realistic system configurations and fault conditions. ï· To validate the new model against existing built-in EMTP-type simulator models, demonstrating improvements in computational efficiency and numerical accuracy, especially under unbalanced fault conditions. 1.3.2 Objective II: Enhancement of the CP-VBR synchronous machine model through multirate (MR) algorithms: ï· Implementation of multirate solutions that distinctly separate fast stator and slow rotor/mechanical subsystems. This involves heterogeneous integration methods and various prediction algorithms to accurately model slow subsystem dynamics while significantly reducing computational overhead. ï· To demonstrate considerable computational advantages and accuracy improvements, test the proposed MR CP-VBR models across different scenarios, including single-machine-infinite-bus systems and complex marine power systems. In achieving these objectives, this research contributes toward more efficient computational methods in electromagnetic transient simulations. The advanced modelling techniques developed here enable faster simulations of larger, more complex electrical systems, enhancing both the feasibility and accuracy of EMTP analyses for modern power systems, thereby bridging mathematical modelling and practical engineering implementations. . 4 Chapter 2: Interfacing of SR CP-VBR Synchronous Machine in EMTP-Type Solutions Synchronous machine modelling starts with the representation of their state space equations, then adapting the model to EMTP-Type solutions by discretizing the states using trapezoidal; hence, depending on the type of state-space equations used, we can name it as VBR, PD, or QD model of synchronous machines. In this thesis, we work with the VBR implementation of the synchronous machines due to their easier abc interfacing than the qd0 models. This chapter primarily reviews the literature on the CP-VBR formulation of the synchronous machine by first exploring the state space representation of the VBR model, going through the concepts of main flux saturation, and finally, deriving the SR CP-VBR SM for EMTP-Type solutions. This chapter also extends the previous work in [14] and makes a contribution by generalizing the saturable constant-parameter voltage-behind-reactance (CP-VBR) model with an arbitrary number of qd-axes damper windings and accommodate the generator stator winding neutral circuit connection according to Figure 2.1 and the IEEE Guide [19]. Generally, the interfacing conductance is derived as a 3-by-3 matrix, which assumes a grounded neutral point of the machineâs stator winding. However, in many practical configurations of synchronous generators, the generator neutral point is connected to the ground through a particular circuit with impedance (or even a transformer) to reduce the transient over-voltages, reduce fault currents, and detect the ground and/or unbalanced faults in the system [20]. According to the IEEE Guide for the Application of Neutral Grounding, the monitoring of the neutral-to-ground current may be done using a special distribution transformer, as depicted in Figure 2.1. 5 Figure 2.1 The generator connection circuit for fault detection using neutral to ground current [19]. The transformer used in the neutral connection, according to Figure 2.1, typically has a lower voltage rating than the generatorâs rated voltage. The impedance connected to the secondary of the transformer depends on the application and type of protection used, and it could range from low to high impedance, inductive or capacitive, etc. Therefore, it is essential to consider generatorsâ exact neutral point connection circuits in system studies that may encounter unbalanced faults and protection operation. 2.1 State Space VBR Model To derive the subject model, the VBR SM model [16] in state space form is presented briefly. The network interfacing stator voltages are represented as ï¨ ï© ï¨ ï©abcs abcs r abcs abcs r abcs abcsj jp ï± ï± ï¢ï¢ï© ï¹ï½ ï« ï«ï« ï» iR eL iv , (2.1) where ï¨ ï©abcs rj ï±L , and ï¨ ï©abcs rj ï±R are the time and saturation-dependent stator resistance and inductance matrices, respectively, given by ï¨ ï© ï¨ ï© ï¨ ï©10 00 00 0lsabcs r s r ls s rjmqj jlsmdL LL LLï± ï± ï±ï­ï¢ï¢ï© ï¹ï«ïª ïºï¢ï¢ï½ ï«ïª ïºïª ïºï« ï»L K K , (2.2) 6 ï¨ ï© ï¨ ï© ï¨ ï©10 00 00 0sabcs r s r s rjmqsj jmd sr rr rrï± ï± ï±ï­ï¢ï¢ï© ï¹ï«ïª ïºï¢ï¢ï½ ï«ïª ïºïª ïºï« ï»K KR , (2.3) where the sub-transient qd-axis magnetizing inductance jmqLï¢ï¢ , jmdLï¢ï¢ , and resistances jmqrï¢ï¢ , jmdrï¢ï¢ are provided in Appendix-A as (A.2), and (A.7). The sub-transient emf in (2.1) is represented as ï¨ ï©13 2abcs s r qdsï±ï­ï´ï¢ï¢ ï¢ï¢ï½Ke e , (2.4) where qdsï¢ï¢e is the sub transient emf in the qd-frame, and expanded as ï¨ ï© ï¨ ï©1 2 3, .j j jqds r qdr r resqd qdrï¢ï· ï· ï·ï¢ï¢ ï½ ï« ï«e λA vA λ A (2.5) The coefficients 1jA , 2jA , and 3jA are provided in Appendix-A as (A.9), (A.10), and (A.11). The state equation for the rotor flux linkages qdrλ is given by 1 2 1 ,j j jqdr qdr qds resqd qdrp ï«ï«ï½ ï­Î» v B B λ C λi (2.6) where the coefficient matrices 1jB , 2jB , and 1jC are expanded in Appendix-A as (A.12) and (A.15). Finally, the mechanical subsystemâs rotor speed and angle are represented as ï¨ ï© ,2r e mPp T TJï· ï½ ï­ (2.7) and ,r rpï± ï·ï½ (2.8) respectively. The electromagnetic torque eT is given by ï¨ ï©3 ,4e md qs mq dsPT i iï¬ ï¬ï½ ï­ (2.9) where 7 3TT jmqd mq md qs ds qdrj jmq mdL i L iï¬ ï¬ ï¢ï¢ ï¢ï¢ï© ï¹ï© ï¹ï½ ï½ ï«ï« ï» ï« ï»Î» A λ . (2.10) 2.2 Main Flux Saturation in EMTP-Type Programs In order to model the saturation, the usual practice in EMTP-Type machine models [21], [13], [22], [14], [23] is to represent the main-flux saturation as a piecewise linear function of the flux linkage mï¬ at various d-axis magnetizing currents. This curve can be experimentally achieved by performing an open circuit test on the SM, and results in distinct magnetizing inductance DL for each saturation segment Figure 2.2. A superscript j in ï¨ ï© jï represents the saturation-dependence of a variable at the jth saturation segment, mi represents the magnetizing current, and resï¬ representing the residual flux linkage. To accommodate the q-axis saturation, a main-flux saturation approach is used, assuming the rotor to be isotropic with a constant saliency factor SF defined by the assumption that the ratio between the q and d-axis inductances is constant, 1 1/ / ,j jmq md mq mdSF L L L Lï½ ï½ (2.11) for unsaturated 1j ï½ or any given saturated operating region j of the SM. With the above considerations, we can now define the resulting magnetic flux linkage as ï¨ ï©22 / ,m md mq SFï¬ ï¬ ï¬ï½ ï« (2.12) Figure 2.2 Piecewise linear saturation curve of d-axis saturation [24]. 8 Figure 2.3 Phasor representation of main-flux saturation, with projections along the qd-axes. for every magnetizing current mi given by ï¨ ï©22 .m md mqi i SFiï½ ï« (2.13) The illustrated phasor of the main-flux saturation, along with its respective projections along the qd-axes, is depicted in Figure 2.3 and can be derived from the equations (2.11)â(2.13) as 2 TTmqd mq m j jD Dd S L LFï¬ ï¬ ï© ï¹ï© ï¹ï½ ï½ï« ï» ï« ï»Î» , (2.14) and the residual flux linkages are given by cos .sinj jresq res resresqd mqdjresd mresSFï¬ ï¬ ï¢ ï¬ï¬ ï¬ï¬ ï¢ï© ï¹ï© ï¹ï½ ï½ ï½ïª ïºïª ïºï« ï» ï« ï»Î» λ (2.15) where ï¢ is the saturation angle between the scaled q-axis flux linkage and the main flux linkage mï¬ where j jm D m resL iï¬ ï¬ï«ï½ . (2.16) 2.3 SM EMTP-Type VBR Model In this section, the derivation of a VBR SM model for the purpose of EMTP-type simulation is discussed in detail [13]. Since the network in an EMTP-Type program is almost always discretized using the trapezoidal integration rule, this chapter also assumes a trapezoidal discretized stator, rotor flux linkages, and mechanical equations for the SR SMâs. 9 The mechanical subsystem can directly be discretized from (2.7) and (2.8) as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï©4r e m e mP ht T t T t T t h T t hJï· ïï½ ï­ ï« ï­ ï­ ï­ , (2.17) and ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï© ,2r r rht t t hï± ï· ï·ï½ ï« ï­ (2.18) respectively. To derive the SR VBR SM, it is obvious from (2.1) that the discrete solution for the stator interfacing equation requires a discrete sub-transient equation ï¨ ï©abcs tï¢ï¢e , hence requiring ï¨ ï©qds tï¢ï¢e (2.5). To achieve this, the state space equation of the rotor flux linkages (2.6) is firstly discretized using trapezoidal integration over h, with the resulting terms rearranged as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï» ï½12 21 ,j j jqdr qdr qds resqd ht h t t tï­ï½ ï«ï© ï¹ ï«ï« ï«ï»Î» B λiC v B λ (2.19) where the coefficient 2jC is expanded in Appendix-A, as (A.16), and the history flux is combined as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï» ï¨ ï© ï¨ ï© ï¨ ï©ï½31 1 1 12 1 2j j j jh qd e q rr qd r d ds sqt h t h t h t h h t hï­ ï­ ï­ ï­ï© ï¹ï½ ï­ ï« ï­ï« ï«ï­ï» ï« ï­Î» v B λiC B λ C , (2.20) where jâ1 depicts the operating saturation segment during the previous time step, the saturation and time step dependent coefficient 3jC is given in Appendix-A, as (A.16). Now substituting (2.19) instead of the state qdrλin the sub-transient equation, the new discrete sub-transient ï¨ ï©qds tï¢ï¢e can be expressed as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï©1 2 3, , , ,rj j jr qds r qdr r resqdqds htt h t h t h tï¢ï· ï· ï· ï·ï¢ï¢ï¢ï¢ ï¢ï¢ï½ ï« ï« ï«eÎ»e D i D v D eï±ï´ï´ï´ï´ï´ï´ï´ï´ï²ï´ï´ï´ï´ï´ï´ï´ï´ï´ï³, (2.21) where the stator and rotor-dependent ( rï¢ï¢e ) terms are separated to later derive the CP-VBR model. The coefficients 1jD , 2jD , and 3jD are expanded in Appendix-A, as (A.18), (A.19), and (A.20), respectively. The sub-transient history term hï¢ï¢e in (2.21) is given by 10 ï¨ ï© ï¨ ï© ï¨ ï©1 r hh jt tï·ï¢ï¢ ï½e A λ . (2.22) Finally, the stator interfacing equation (2.1) is discretized with a trapezoidal rule across the time step h and combined with (2.22) ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 2VP VPabcs abcs abcs s r r sht tt t tï±ï­ï´ï¢ï¢ ï¢ï¢ï½ ï« ï«i Kv R e e , (2.23) where the variable interface resistance PabcsVR is defined as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©1 13 2 2 32 ,jabcs abcs r abcs r s r rVP j jr st hhï± ï± ï± ï· ï±ï­ï´ ï´ï½ ï« ï«R R L K D K , (2.24) the history voltage is given as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 2VP t hsh habcs abcs s r rVabcsPt t h t ht t h ï±ï­ ï­ï´ï¢ ï½ ï­ï¢ ï¢ï¢ï« ï­ ï­ ï­i Ke R e v , (2.25) the history resistance matrix can be inferred as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©1 1 12113 2 32 ,t h t h t h j t h t hhabcs abcs r abcs r s r rjs rVP jt hhï± ï± ï± ï· ï±ï­ ï­ï­ ï­ ï­ ï­ ï­ ï­ï´ ï´ï­ï½ ï­ ï«R LR K D K . (2.26) Now, the final EMTP-Type Norton interface equation for the VBR model can be given as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©VPabcs abcs abcsV hPt t tt ï«ï½ v Ji G , (2.27) where the variable conductance for a fixed time step h is given by ï¨ ï© ï¨ ï© 1abcs abP sVP V ct tï­ï© ï¹ï½ ï« ï»G R , (2.28) and the history current injection VPhJ is given by ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 2VP VPh abcs s r r sPhV t tt t ï±ï­ï´ï© ï¹ï½ ïª ïºï¢ ï¢ï¢ï« ï»ï¢ï«J G K e e . (2.29) 2.4 EMTP-Type CP-VBR Model In well-optimized EMTP-type programs, the inversion of the network conductance matrix PabcsVGheavily determines the overall simulation speed. Due to this, some commercial software (RTDS/RSCAD) 11 [10] does not allow for a variable interfacing conductance. Hence, to implement the SM model on these platforms and to enable faster inversion in software that accepts both constant and variable interface conductance (PSCAD/EMTDC ), a CP-VBR SM model was derived in [14]. In this section, the derivation of the CP-VBR SM model is introduced. To achieve constant interfacing, the variable resistance matrix in (2.24) is divided into a constant interfacing resistance abcsR and a variable difference resistance jqdsïR , which accounts for the difference between the variable and the constant parameter implementations ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 2 2 3,abcs abcs qsPr s rV jdst tï± ï±ï­ï´ ï´ï« ïï½R R K KR (2.30) where it can be observed that a smaller jqdsïR would result in a more accurate network conductance at any given time. To determine the optimal value for abcsR and jqdsïR , the individual resistances in (2.26) are analyzed as ï¨ ï© ï¨ ï©ï¨ ï©ï¨ ï© ï¨ ï©ï¨ ï©ï¨ ï©ï¨ ï© ï¨ ï©ï¨ ï©ï¨ ï©1 113 33 2 2 31 12 1,1 , 1,22 ,2, 2,1 2, 2j jrabcj jmq mqj js s ls s r s rj jr md mdr L h hht r Lh h r L hhï·ï± ï±ï·ï­ï´ï´ ï´ï© ï¹ï¢ï¢ ï¢ï¢ï« ï«ïª ïºï¦ ï¶ï½ ï« ï« ïª ïºï§ ï·ï¨ ï¸ ïª ïºï¢ï¢ ï¢ï¢ï« ï«ïª ïºï« ï»RDKIDKD D (2.31) For a fixed time step h, this is compactly represented as ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï© ï¨ ï©311 1221323 22 3122 .rabcs s ls s r s rrj jVPj jr rt r Lr rhï·ï± ï±ï·ï­ï´ï´ ï´ï© ï¹ï¦ ï¶ï½ ï« ï« ïª ïºï§ ï·ï¨ ï¸ ï« ï»KI KR (2.32) The variability of the conductance in the transformation process occurs due to three distinct factors. Firstly, the off-diagonal elements in the variable part ï¨ ï©21 rjr ï· and ï¨ ï©12 rjr ï· lead to rotor-dependent resistance after transformation. Since these elements cannot be approximated, they are separated from the constant part to ensure accuracy. Secondly, the variability introduced by the saturation segment can be approximated by 12 using constants from the unsaturated segment, i.e. j=u=1. The differences are then combined into the difference resistance mentioned above, improving accuracy. Finally, the variability resulting from mismatches in the diagonal entries is approximated by averaging the diagonal values. Applying this approach, the variability in (2.32) can be rearranged as ï¨ ï© ï¨ ï©ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©11 12 11 2221 21 13 2 2 3 3 22 2 3,2r csj j u uonsts r s r r st rjf qdrsj jr r r r tr rï·ï± ï± ï± ï±ï·ï­ ï­ï´ ï´ ï´ ï´ï© ï¹ ï¦ ï¶ï«ï½ ï«ïª ïº ï§ ï·ï¨ ï¸ï« ï»ïK K KR RK (2.33) where the constant transformation consttfR is expanded as, ï¨ ï© ï¨ ï©311 221121322112 2 202 ,202u uu utfconsts r su u rr rr rr rï± ï±ï­ï´ ï´ï© ï¹ï«ïª ïºï¦ ï¶ï« ïª ïºï½ï§ ï· ï«ïª ïºï¨ ï¸ïª ïºï« ï»K KR (2.34) with the constant result given in Appendix-A, as (A.21), and the difference matrix as ï¨ ï©ï¨ ï©11 12, 11 222 2221.2rrj j u uj CPqds j jr r r rr rï·ï·ïï© ï¹ ï«ï½ ï­ïª ïºï« ï»IR (2.35) From(2.32)â(2.35) we can now write the variable conductance as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©311 22 13 33 2 22 ,2aCb sPcconstabcs s ls su uVP jtf r rqds sr rt r L thï± ï±ï­ï´ï´ ï´ï¦ ï¶ï«ï¦ ï¶ï½ ï« ï« ï«ï§ ï·ï§ ï·ï¨ ï¸ ï¨ ï¸ïRR R K KI Rï±ï´ï´ï´ï´ï´ï²ï´ï´ï´ï´ï´ï´ï³ (2.36) Substituting (2.36) into the stator interfacing equation (2.23) yields ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 2,CPabcs abcs abcs sh s drCP jqds sr qt t tt t tï±ï­ï´ï¢ï¢ ï¢ï¢ï© ï¹ï½ ï« ï« ï«ïï« ï»v iR Ri Ke e (2.37) the stator history subtransient ï¨ ï©VPsh tï¢ï¢e in (2.23) must be updated to ï¨ ï©CPsh tï¢ï¢e in (2.37) with changes similar to its non-history counterpart, or else this imbalance would result in a diverging solution, following this, ï¨ ï©CPsh tï¢ï¢ecan now be expanded to ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 2,CP jCP t hsh habcs abcs s r r ahqds q sds bct t h th ht t t hï±ï­ ï­ï´ï© ï¹ï¢ï¢ ï¢ï¢ï½ ï« ï­ ï«ï ï­ï­ ï­ ï­ï« ï»e R e Ri K i vï¥ (2.38) 13 where the constant resistance and the variable difference resistance are derived similarly to their non-history counterpart and are given by 3 11 2232 ,2u uconsCP h hthabcs s l fs tr rr Lh ï´ï¦ ï¶ï«ï¦ ï¶ï½ ï­ ï« ï§ ï·ï§ ï·ï¨ ï¸ ï¨ ï¸R I R (2.39) and ï¨ ï©ï¨ ï©ï¨ ï©11 12, 11 2221 222,2t hrtj j u uhj CP hjhhhqds jr hr r r rtr rï·ï·ï­ï­ïï© ï¹ ï«ïª ïºï½ ï­ïª ïºï« ï»IR (2.40) where 11jhr , and 22jhr are 11 11 22 224 4, ,j jmqj j j j mdh hL Lr r r rh hï¢ï¢ ï¢ï¢ï½ ï­ ï½ ï­ , (2.41) Finally, when the stator equation (2.37) is rearranged for EMTP-type Norton-based interfacing, it results in ï¨ ï© ï¨ ï© ï¨ ï©CPabcs abcs abcs hCPt tt ï«ï½ v Ji G , (2.42) where, 1abcs ab sCP CPcï­ï© ï¹ï½ ï« ï»G R , (2.43) and the history current injection hJ is given by ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï» ï½13 2.CPh abcs s r sr shCP jqd qdst tt tï±ï­ï´ï¢ï¢ ï¢ï¢ï© ï¹ï«ï ï«ï« ï»ï½J e eG K R i (2.44) As seen in (2.44), the history current injection depends on the stator current at the current time step. To eliminate this dependency, the current in the qd reference frame is predicted using a linear prediction algorithm, resulting in a history current injection that is independent of current time step parameters: ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï» ï½13 2.CP jqds qdsCP CPh abcs s r r sht tt tï±ï­ï´ï© ï¹ï¢ï¢ ï¢ï¢ï«ï«ï½ ï ï«ï»J e R eG K iï¥ (2.45) 14 With equations (2.42)â(2.44), we have now derived the CP-VBR model of the SM. This derivation serves as the foundation for the implementation of the CP-VBR model in the following section and sets the stage for the MR theory and model discussed later. 2.5 Interfacing of the Machine Model Any component designed for an EMTP-type simulator can be represented as a combination of a discrete conductance G and history current source anhJ shown in Figure 2.4. The definition of the conductance can either be as a branch-based component Figure 2.4(b), where the values are described across the component, or a node-based component Figure 2.4(c) with the parameters defined with respect to a common ground. Figure 2.4 Physical component to EMTP-type component representation: (a) physical component, (b) branch-defined component, and (c) node-defined component. The interfacing of the existing synchronous machine models does not conform to an accurate physical representation. This is due to the assumption that the neutral is assumed to be solidly grounded, with no access to it. To interface a machine model in an EMTP-type program, the final discretized equation is put into the nodal form [1], hï« ï½Gv i i . (2.46) 15 This three-phase conductance is usually defined across the terminals of the generator and the ground in a node-based component representation. Figure 2.5 Expansion of the conductance matrix: (a) solidly grounded component, and (b) component with separate neutral node. From the final interfacing equation of CP-VBR (2.42) 1aa ab acabcs abcs baca cbCP CPbb cccbg g gg g gg g gï­ï© ï¹ïª ïºï© ï¹ï½ ï½ ï½ï« ï» ïª ïºïª ïºï« ï»G G R , (2.47) and history injection currents are expanded as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© .Tag bg cghP hh hC t j t j t j tï© ï¹ï½ ï« ï»J (2.48) The 3 3ï´ conductance matrix in (38) implies an interfacing of the machine model using its Norton equivalent, as depicted in Figure 2.5(a), which also assumes a grounded neutral point (because the nodal voltages are computed with respect to the ground). This 3 3ï´ conductance matrix and history injection currents are placed into the overall network nodal equation. This interface enables a simultaneous solution of the machine equations with the rest of the network in an EMTP solution. The traditional way of achieving this in a node-defined component is by expanding the conductance matrix to include the neutral node Figure 2.5. 16 To generalize the model interfacing and account for a floating neutral, the stator voltages of the machine are now defined between nodes abc and n, as shown in Figure 2.5(b). Redefining the voltage vector vabcs to include the neutral node, we get ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© Tabcn a b c nt t t tv v tv vï½ ï© ï¹ï« ï»v . (2.49) Also, the history current vector is expanded to include the current flowing through the neutral node as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© .Tneutral ag bg cg ngh h hh ht j t j t j t j tï© ï¹ï½ ï« ï»J , (2.50) where ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï©ng ag bg cgh h h hj tj j jt t tï½ ï­ ï« ï« . (2.51) The expanded voltage and current equations now require a 4 4ï´ matrix nconstG for the proposed nodal interface with neutral, where 3 3 11 3 1 13aa abCP Tneutralca anba bb bc bnabcs xncca cb cc cnna nonstn nb n nx nc ng g g gg g g gg g g ggg g g gï´ ï´ï´ ï´ï© ï¹ïª ïºï© ï¹ ïª ïºï½ ï½ïª ïº ïª ïºï« ï» ïª ïºï« ï»G gGg. (2.52) Here, applying algebraic manipulation to (2.47) using (2.49) and (2.51), provides us with the values for additional entries in the extended conductance matrix [28] as follows ï¨ ï©ï¨ ï©ï¨ ï©ï¨ ï©ï¨ ï©ï¨ ï©,an aa ab ac na aa ba cabn bb bc nb bb cbcn cb cc nc bc ccba aca a cbg g g g g g g gg g g g g g g gg g g g g g g gï½ ï­ ï« ï« ï½ ï­ ï« ï«ï¬ ï¬ï¯ ï¯ï½ ï­ ï« ï« ï½ ï­ ï« ï«ï­ ï­ï¯ ï¯ï½ ï­ ï« ï« ï½ ï­ ï« ï«ï® ï®, (2.53) and ï¨ ï© ï¨ ï© ï¨ ï©nn aa ab bc ba bb bc ca cb ccg g g g g g g g g gï½ ï« ï« ï« ï« ï« ï« ï« ï« . (2.54) Finally, the Norton interface equation for the CP-VBR is constructed in the required form as ï¨ ï© ï¨ ï© ï¨ ï©neutral neutralconst abcn h abcnt t tï« ï½G v j i . (2.55) 17 2.6 Computer Study Figure 2.6 A single-machine system with generator grounding. Figure 2.7 A single-machine system implemented in PSCAD: (a) using the built-in qd0 model, and (b) using the proposed generalized CP-VBR model. A single-machine system composed of a synchronous generator with a grounding inductance Lng, short line impedances L1, L2, and an ideal ï/Y transformer feeding an infinite bus, as depicted in Figure 2.6, is considered to demonstrate the proposed model. This system has been implemented in PSCAD/EMTDC, in one case using its own classical built-in qd0 model, as shown in Figure 2.7(a); and in another case using the proposed generalized CP-VBR model, as depicted in Figure 2.7(b). The synchronous machine and system parameters are provided in Appendix B. The magnetic saturation has not been considered here due to the focus on the model interfacing. The PSCAD implementations of both cases are shown in Figure 2.7. 18 In the following study, the system is assumed to initially operate in a steady state defined by a mechanical torque Tm = 1.289Ã106 Nm and field excitation voltage Vfd = 14.876 V. Then, at t = 0.06s, a fault is introduced at phase-a to ground, and the simulation is continued to run till t = 0.2s. Both models are run using time steps of 1µs, 50µs, 100µs, 200µs, 250µs and 500µs, respectively, to evaluate their respective accuracy. The transient responses of the stator current ias, the current through the grounding reactance ing, and the electromagnetic torque Te by the subject models are shown in Figures 2.8 and 2.9. The reference solution is obtained using either of the models simulated using a very small time step of 1µs, at which both models converge to the same solution. As observed in Figure 2.8, the current flowing from the neutral to the ground remains zero until the fault. When the fault is applied, we observe a surge in the current through the neutral, implying an unbalanced fault. Figure 2.8 also shows that when the simulation time step is sufficiently small, all models give similar solutions that converge to the reference solution. To better see the difference, a magnified view of the fragments of Figure 2.8 are also plotted in Figure 2.9, which reveals that at larger time step sizes, such as 500µs, the conventional qd0 model is visibly less accurate compared to the proposed CP-VBR model. To further investigate the accuracy of the considered models, the 2-norm cumulative error [27] has been calculated for the solution trajectories of the stator phase current ias and the grounding current ing, for different time step sizes, and the results are plotted in Figure 2.10. As observed in figure 2.10, the traditional built-in qd0 model very quickly loses its accuracy, and the error in ias reaches 10% for the time steps 500µs. Its error in the grounding current ing remains reasonable under 1%. 19 Figure 2.8 Transient response of the generator to a single-line-to-ground fault: (a) stator current, (b) neutral-to-ground current, and (c) electromagnetic torque. Figure 2.9 Magnified views from Figure 2.8: (a) stator current, (b) neutral-ground current, and (c) electromagnetic torque. At the same time, the proposed CP-VBR model retains a fairly good accuracy with an error below 1% for both phase and grounding currents. -20020030600 0.05 0.1 0.15 0.2Time (s)1.21.31.4(a)(c)(b)see Fig 2.9(a)see Fig 2.9(b)see Fig 2.9(c)0.0675 0.068 0.06856262.50.064 0.0645 0.065 0.0655 0.066Time (s)1.321.331.34ReferenceBuilt-in qd0 (500 s)Proposed CP-VBR (500 s)0.073 0.0735 0.074 0.0745 0.075 0.07550123ReferenceBuilt-in qd0 (500 s)Proposed CP-VBR (500 s)ReferenceBuilt-in qd0 (500 s)Proposed CP-VBR (500 s)(a)(b)(c)20 Figure 2.10 Comparison of norm-2 error for the solution trajectory of stator phase current ias and neutral current ing for the PSCADâs built-in qd0 and the proposed CP-VBR models at different time step sizes. 21 Chapter 3: MR CP-VBR Machine models In this chapter, we will first explore the implementation of a standard component model in EMTP-Type programs, then expand it to provide on how MR models are interfaced with SR-based solvers, moving on to the theory and conception of MR algorithm, and finally derive the MR CP-VBR SM Machine model, with different time steps and integration methods. 3.1 EMTP-Type Solver Implementation Implementing the proposed machine models requires careful attention. While the literature addresses this topic, it seldom provides enough detail for direct implementation in commercial software. In this section, the procedure for interfacing the model in PSCAD is covered comprehensively. To the authorâs knowledge, following this section will facilitate more straightforward implementation of the model in similar EMTP-type simulators. This discussion also aims to provide foundational knowledge of EMTP-type modelling, which is typically only available in software technical documentation and/or standard example libraries provided by the developers. A thorough understanding of the simulation environment is essential for developing effective and integrable models. A typical EMTP-type solver follows a consistent algorithmic framework [25], although the specific methods may vary based on the solverâs complexity. At its core, the implementation process remains consistent across various simulators. Before the simulation begins, the solver program calls an initialization function. This function defines the initial operating conditions for all components in the network and precomputes the constant interfacing conductance to minimize computational overhead during each time step. In the case of the proposed algorithm, this function also calculates the saturation-dependent coefficients, constant coefficients, and the constant interfacing conductance. The three primary functions are ï· Pre-network solution function, 22 ï· Network solution function, ï· and post-network solution function. The pre-network solution function computes each componentâs history, current injection, and variable interfacing conductance. These values are aggregated to form the overall network conductance and network injection sources, and the system proceeds to solve for the unknown nodal voltages. Figure 3.1 Simplified program flow in an EMTP-Type simulator, considering a constant network conductance matrix. Once the nodal voltages for the current time step are determined, the solver proceeds to the post-network solution function, where time-dependent non-network variables are computed. 3.2 MR Algorithm Keen readers might have noticed that the solution algorithm presented in Figure 3.1 does not account for an MR-based simulation. Unfortunately, current state-of-the-art simulators are not inherently designed to support MR models and multistep simulations. For example, RTDS/RSCAD features supersteps and substeps but lacks the necessary interfacing to implement MR algorithms fully. On the other hand, PSCAD entirely neglects MR simulation and only supports parallel simulation, which is of minimal benefit when exploring MR models. The basic idea behind MR implementation is to algorithmically separate variables that exhibit different time constants. This separation can be leveraged to our advantage by dividing the model into two 23 sections, each executing at different intervals. While the focus of this thesis is not on an extensive implementation of MR algorithms but rather on their application, we assume only two different time steps: 1. Fast time step hf Used for network variables (terminal voltages and currents) 2. Slow time step hs Used for mechanical and rotor transients, as well as for saturation segment changes. Following the notion of a simplified MR implementation, the MR ratio (M) between the slow and fast time steps is considered to be a natural number ,sfhMhï½ ïï¥ (3.1) a rudimentary depiction of the MR algorithm in any EMTP-type solvers for an MR ratio of three is given in Figure 3.2. Figure 3.2 MR solution algorithm implementation in an EMTP-type solver. 24 Figure 3.3 A typical MR simulation timeline depicting commonly used terms for an MR ratio of three. 3.2.1 MR counter The chosen commercial simulator does not inherently support MR simulations. A counter z is defined within the algorithm. This counter calculates the distance from the previous major time step Figure 2.3 and computed every time step of the simulation. The distance z is determined using the current simulation time t, the network simulation time step hf, and the MR ratio defined by the user for the component. Most EMTP-type simulators provide easy access to the current simulation time and hf, making this implementation straightforward. Mathematically, the counter z is computed using the mod operator as follows: ï¨ ï©mod / , ,fz t h Mï½ (3.2) i.e. the reminder of ï¨ ï©/ ft h over M. 3.2.2 Interfacing of multistep algorithms The dependent variables between fast and slow variables must be decoupled to ensure proper computation. Due to differing time constants, the slow transients can be safely assumed to vary linearly between major steps. Hence, linear prediction is done using the distance z on the slow transient variables as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï©2 ,s s sf f fzt t h ft h t hMï½ ï­ ï« ï­ ï­ ï­ï¥ (3.3) Figure 3.4 depicts (3.3) graphically where the prediction is done for an MR ratio of three with the prediction for the two minor steps z = {1, 2} or t =.{hf, 2hf } 25 Figure 3.4 Slow transient prediction and accurate fast transient area computation depicted graphically. In contrast, the slow transients calculated during the major step sometimes require the area under the fast transient curve. Accurate computation of this area is essential for precise slow transient solutions. Since the minor step computation provides all stator currentâs values for the minor step, the area under a fast-transient current curve can be more accurately represented using an expanded trapezoidal rule ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©112sst Mt hf fiff ihfd t f ht h h f tï´ ï´ï­ï½ï­ï» ï« ï­ ï«ï© ï¹ï« ï­ï» ï¥ï² . (3.4) This can also be graphically verified from Figure 3.4 (striped area) to be more accurate than the trapezoidal integration over the time step hs represented by the solid area. 3.3 MR-based CP-VBR SM Models Before we proceed with the model derivation, it is essential to distinguish the slow and fast transient variables. In the slow transient variables, we include the rotor and mechanical subsystems. Meanwhile, the stator electrical network variables are lumped into the fast transient. We can, in theory, implement a three-step system with mechanical transients and rotor transients with different ratios of execution to preserve the simplicity of the models and their respective studies; a single mass mechanical system is considered, in which case combining the rotor and mechanical system would suffice. 26 3.3.1 MR-TT CPVBR SM model In this section, we will derive an MR homogeneous CP-VBR SM model, where the fast and slow transient variables are discretized using the trapezoidal integration rule, but with varying time step hf and hs, respectively. To differentiate variables in this model from the others, all the variables corresponding to this model have a superscript ï¨ ï©TTï . 3.3.1.1 Slow subsystem To derive the MR-TT CP-VBR SM, we start by discretizing the rotor-flux state space equation (2.6) using the trapezoidal rule, but over the slow time step hs, ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï» ï¨ ï©ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï½2 2 1 111 12stTT TT jsqdr qdr qdr qdr s qdst hj j j TT j TTresqd resqd s qdr qdr st t t t h dt t h hht tï´ ï´ï­ï­ ï­ï© ï¹ï­ ï« ï­ï« ï»ï© ï¹ ï© ï¹ï« ï­ ï« ï« ï­ ï»ï½ï« ï» ï«ï« ï«ï²Î» λ v v BB λ Biλ C λ C λ. (3.5) As discussed previously, the area under the fast transient variable qdsi is expanded using (3.4) and rearranged to get ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©2 211 ,1TT j j jqdr s qdr qds resqd hTTt h t t tMï­ ï¬ ï¼ï½ ï« ï«ï­ ï½ï©ï« ï» ï¾ï«ï®ï¹Î» B λ λiC v B (3.6) where ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï© ï¨ ï©ï½1121 1112 131 2j j jh qdr s qds s qds ressMTTfjkTTqq srddt h t h t h t kh t hMh t hï­ ï­ï­ï½ï­ï­ï© ï¹ï¬ï© ï¹ï½ ï­ ï« ï­ ï« ï­ ï­ï­ ïªï«ï«ï»ï«ïºï« ï» ï®ï­ï¥i iλ C v B B λC λ. (3.7) A general assumption of no saturation segment change during the minor step is considered here, similar to the literature [18]. Hence, j-1 represents the operating saturation segment during the previous major time step t-hs rather than any intermediatory fast transient time steps{t-hf, t-2hf, â¦ , t-(M-1)hf }. Substituting the flux linkages (3.6) into the sub-transient equation to find the sub-transient voltage at current time (2.5) we get 27 ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï©2 311 , , , ,TTrj j jr s qds r s qdr r s resqTdT TTqds htt h t h t h tM ï¢ï· ï· ï· ï·ï¢ï¢ï¢ï¢ ï¢ï¢ï½ ï« ï« ï«eÎ»De D i v D eï±ï´ï´ï´ï´ï´ï´ï´ï´ï´ï²ï´ï´ï´ï´ï´ï´ï´ï´ï´ï³ (3.8) where ï¨ ï© ï¨ ï© ï¨ ï©1jTT TTh r ht tï·ï¢ï¢ ï½ A λe . (3.9) The mechanical subsystem in (2.7) and (2.8), are subsequently discretized as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ,4TT TT sr r s e m e s m sP ht t h T t T t T t h T t hJï· ï·ï½ ï­ ï« ï­ ï« ï­ ï­ ï­ï© ï¹ï« ï» (3.10) and ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï© ,2TT TT TT TTsr r s r r sht t h t t hï± ï± ï· ï·ï½ ï­ ï« ï­ ï­ (3.11) respectively, using trapezoidal discretization. 3.3.1.2 Fast subsystem The derivation from here closely follows the CP-VBR modelâs stator voltage derivation. To avoid repetition from Chapter 1, and to invite the reader into component modelling, the derivation of the final results below is left as an exercise with the following guide. ï· Derive the VP VBR stator interfacing equations using the slow transient variables derived in the previous section, and following the methodology provided in (2.23)â(2.26). ï· Use the interfacing equation from the previous step, and divide the interfacing resistance into a constant TabcsTR and a difference , TjqdsTïR part following the methodology in (2.30)â(2.36). ï· Rearrange the terms to form the nodal equations, and form the stator history, along with the Norton equivalent current source TThJ , following the derivation from (2.37)â(2.45). The final stator interfacing equations for the MR-TT CP-VBR SM model are provided below, ï¨ ï© ï¨ ï© ï¨ ï©TTabcs abcs abcs hTTt tt ï«ï½ v Ji G , (3.12) where 28 1abcs ab sTT TTcï­ï© ï¹ï½ ï« ï»G R , (3.13) and ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï» ï½,13 2.TT TT TT TTh abcs s r r shTT jqds qdst tt t tï±ï­ï´ï© ï¹ï¢ï¢ ï¢ï¢ï« ï ï«ï« ï»ï½J e R eG K iï¥ (3.14) The constant resistance TabcsTR , variable difference resistance , TjqdsTïR , and the stator history TTshï¢ï¢e are expanded as 3, ,11 2322 ,2u TT unTTTT co stabc l ts sffsr rr Lh ï´ï¦ ï¶ï§ ï·ï¨ ï¸ï¦ ï¶ ï«ï½ ï« ï«ï§ ï·ï§ ï·ï¨ ï¸R I R (3.15) where ï¨ ï©ï¨ ï© ï¨ ï©, ,11 12, ,2 211 2201 , ,20jmqjmqj TT j TTsj TT j TT jj mdfr jrrfmdLrhr rhr r MLrhï· ï·ï·ï¢ï¢ï© ï¹ï¢ï¢ ï«ïª ïºï© ï¹ ïª ïºï½ ï«ïª ïº ïª ïºï¢ï¢ï« ï» ïª ïºï¢ï¢ ï«ïª ïºï« ï»D (3.16) ï¨ ï©ï¨ ï©, , , ,11 12, 11 22, ,21 222,2rrj TT j TT u TT u TTj TTqds j TT j TTr r r rr rï·ï·ïï© ï¹ ï«ï½ ï­ïª ïºï« ï»IR (3.17) and ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©2,13,ft hTT TTsh habcs abcs f s r rTT j TThqds qf f bcsd a fst h ht h tt tt hï±ï­ï­ï´ï© ï¹ï¢ï¢ ï¢ï¢ï½ ï« ï­ ï«ï ï­ ï­ï« ï»ï­ ï­i Ke R e R viï¥ (3.18) where , ,11 23232 ,2u TT u TTTT h htffconsthabcs s lsr rr Lh ï´ï¦ ï¶ ï«ï½ ï­ ï«ï§ ï·ï§ ï·ï¨ï§ ï·ï¸ï¦ ï¶ï¨ ï¸R I R (3.19) , , , ,11 11 22 224 4, ,j jmqj TT j TT j TT j TT mdh hL Lr r r rh hï¢ï¢ ï¢ï¢ï½ ï­ ï½ ï­ (3.20) 29 ï¨ ï©ï¨ ï©, ,, ,11 12, 11 22, ,21 222.2ffTt hrtj TT j TTu T u TThj TT h hhq j hrds TT j TThr r r rr rï·ï·ï­ï­ïï© ï¹ ï«ïª ïºï½ ï­ïª ïºïª ïºï« ï»IR (3.21) 3.3.2 The MR-TBE CP-VBR SM model A heterogeneous CP-VBR SM model with the slow transients discretized using BE over hs and the fast transients using trapezoidal at hf is derived in this section, further separating the rotor and stator subsystem in terms of the integration system. The variables following this derivation method are labelled with the superscript ï¨ ï©TBEï§ . It is highly recommended that the reader goes through section 2.3.1 before proceeding with this section, as common assumptions and a guide to deriving the below-mentioned equations are provided in detail. 3.3.2.1 Slow subsystem For the MR-TBE CP-VBR SM, we first discretize the rotor-flux state space equation (2.6) using the BE rule over the slow time step hs, resulting in ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 21,12TBE j j jqdr s qdr q hTds resBEqdtMt h t tï­ï© ï¹ ï«ï«ï½ ï»ï¬ ï¼ï« ï«ï­ ï½ï® ï¾Î» C v B λiB λ (3.22) where ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©4 1 11111 2 ,12MTBE TBqdrj jh s ssfkqds s qdst h t h t t khMhhï­ï­ ï­ï½ï¬ ï¼ï© ï¹ï© ï¹ï½ ï­ ï« ï­ ï« ï­ï­ ï½ïª ïºï« ï» ï« ï»ï® ï¾ï¥Î» C λ iB i (3.23) Substituting the flux linkages (3.22) into the sub-transient equation (2.5) to find the sub-transient voltage at the current time we get ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï¨ ï©4 5 61 , , , ,2TBErj jTBE TBEqd sjr s qd r s qdr r s res ds hqtt h t h t h tM ï¢ï· ï· ï· ï·ï¢ï¢ï¢ï¢ ï¢ï¢ï½ ï« ï« ï«ee D i D v D λ eï±ï´ï´ï´ï´ï´ï´ï´ï´ï´ï²ï´ï´ï´ï´ï´ï´ï´ï´ï´ï³ (3.24) where ï¨ ï© ï¨ ï© ï¨ ï©1j rTBE TBEhh t tï·ï¢ï¢ ï½ λe A , (3.25) 30 the coefficients 4jD , 5jD , and 6jD are provided in Appendix-C as (C.1), (C.2), and (C.3). The mechanical subsystem in (2.7) and (2.8), are subsequently discretized as ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ,2TBE TBE sr r s e mP ht t h T t T tJï· ï·ï½ ï­ ï« ï­ï© ï¹ï« ï» (3.26) and ï¨ ï© ï¨ ï© ï¨ ï©TBE TBE TBEr r ss rt t h h tï± ï± ï·ï½ ï­ ï« , (3.27) respectively, using trapezoidal discretization. 3.3.2.2 Fast subsystem The final stator interfacing equations after the CP-VBR formulation and rearrangement for the MR-TT CP-VBR SM model are provided below, ï¨ ï© ï¨ ï© ï¨ ï©TBEabcs abcs abcs hTBEt tt ï«ï½ v Ji G , (3.28) where 1abcs abTBE TBcsEï­ï© ï¹ï½ ï« ï»G R , (3.29) and ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©ï» ï½3,12TBE jTBE TBE TBE TBEh qa db r q s scs s r sd ht t tt tï±ï­ï´ï© ï¹ï¢ï¢ ï¢ï¢ï«ï«ï½ ï ï«ï»J e R eG K iï¥ . (3.30) The constant resistance EabcsTBR , variable difference resistance , EjqdsTBïR , and the stator history TBEshï¢ï¢e are expanded as 3, ,11 2232 ,2u TBE u TBETBEtconstabcs s l fsfr rr Lh ï´ï¦ ï¶ ï«ï½ ï« ï«ï§ ï·ï§ ï·ï¨ ï¸ï¦ ï¶ï§ ï·ï¨ ï¸R I R (3.31) where 31 ï¨ ï©ï¨ ï© ï¨ ï©, ,11 124, ,21 22201 , ,220jmqjmqj TBE j TBEsj TBE j TBE jfr jrfj mdmdrLrhr rhr r MLrhï· ï·ï·ï¢ï¢ï© ï¹ï¢ï¢ ï«ïª ïºï© ï¹ ïª ïºï½ ï«ïª ïº ïª ïºï¢ï¢ï« ï» ïª ïºï¢ï¢ ï«ïª ïºï« ï»D (3.32) ï¨ ï©ï¨ ï©, , , ,11 12, 11 22, ,21 222.2j TBE j TBE u TBE u TBEj TBEqds j TBE jrTBErr r r rr rï·ï·ïï© ï¹ ï«ï½ ï­ïª ïºï« ï»IR (3.33) and ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©13 2, ,ft hTBE TBEsh habcs abcs f s r rTBE j TBEhqds qf ds f abcs ft t ht h t t t hhï±ï­ï­ï´ï© ï¹ï¢ï¢ ï¢ï¢ï½ ï« ï­ ï«ï ï­ ï­ï« ï»ï­ ï­i Ke R e R viï¥ (3.34) where , ,11 223 32 ,2u TBE u TBETBE h hts fconsthabcs s lfr rr Lh ï´ï¦ ï¶ ï«ï½ ï­ ï«ï§ ï·ï¦ï§ ï·ï¨ï§ï¸ï¶ï·ï¨ ï¸R I R (3.35) , , , ,11 11 22 224 4, ,j jmqj TBE j TBE j TBE j TBE mdh hL Lr r r rh hï¢ï¢ ï¢ï¢ï½ ï­ ï½ ï­ (3.36) ï¨ ï©ï¨ ï©, ,, ,11 12, 11 22, ,21 222.2ffj TBE j TBEu TBE u TBEhj TBE h hhqds j TBE j ThBEthtrhrr r r rr rï·ï·ï­ï­ïï© ï¹ ï«ïª ïºï½ ï­ïª ïºïª ïºï« ï»IR (3.37) 3.4 MR CP-VBR Implementation In this section, we aim to integrate the model derived from the previous sections with the EMTP-type solver. From the final interface equations in the previous sections (3.12)â(3.21) and (3.28)â(3.37), we can see that the history current is dependent on two slow transient variables, namely rï¢ï¢e and rï· , and as discussed in the MR algorithm, we would have to predict these variables during minor steps using linear predictor defined in (3.3) to get the current step slow subsystemâs solution. To easily represent the final algorithm and visually separate the different transient calculations, algorithms that are executed only during the minor step are represented with blue, the algorithms executed only in the major step are represented by purple, and the steps executed during both in black. 32 âââStart of Pre-network solution subroutine âââ 1) Compute the step z using (3.2). 2) Predict ï¨ ï©r tï· rotor speed in the current step using the values from previous major time steps. 3) Compute the rotor angle ï¨ ï©r tï± using the predicted value of ï¨ ï©r tï· and the value of ï¨ ï©r st hï· ï­ with appropriate integration methods. 4) Predict the magnetizing main flux linkage ï¨ ï©mqd tλ and ï¨ ï©tï¢ï· using a linear predictor. 5) Calculate ï¨ ï©m tï¬ and ï¨ ï©resqd tλ using the predicted values from step two as (2.12) and (2.15) in the current step. 6) Predict the rotor transients ï¨ ï©r tï¢ï¢e , in the current step using the values from previous major time steps. 7) ï¨ ï©r tï¢ï¢e is calculated using the rotor history terms as given in a. (3.8) for MR-TT model, b. (3.24) for the MR-TBE model. 8) Calculate the difference resistance ï¨ ï©jqds tïR using a. (3.17) for MR-TT model, b. (3.33) for MR-TBE model, and multiply with the linearly predicted stator qd current ï¨ ï©qds tiï¥ . 9) Compute the stator history voltage term ï¨ ï©sh tï¢ï¢e using a. (3.18) for MR-TT model, b. (3.34) for MR-TBE model, and finally compute the interfacing history current injection ï¨ ï©h tJ for the nodal solution using a. (3.14) for MR-TT model, 33 b. (3.30) for MR-TBE model. âââEnd of pre-network solution subroutine âââ âââStart of network solution subroutine âââ 10) Solve for the nodal voltages ï¨ ï©abcs tv , usually computed by the commercial solver during each time step. âââEnd of network solution subroutine âââ âââStart of post-network solution subroutine âââ 11) Transform the stator abc current from the network solution into its qd counterpart ï¨ ï©qds ti . 12) Calculate the rotor flux linkages ï¨ ï©qdr tλ using a. (3.6) for the MR-TT model, b. (3.22) for the MR-TBE model, for the current time step. 13) Compute and correct the magnetizing flux linkage ï¨ ï©mqd tλ using (2.10) and recalculate ï¨ ï©m tï¬using(2.12). 14) Check for saturation segment change according to ï¨ ï©m tï¬ if the operating segment does not change, skip to step 16. 15) With the changed segmentâs coefficients, recompute variables in steps 12 and step 13, and recalculateï¨ ï©r tï¢ï¢e with the updated coefficients using a. (3.8) for the MR-TT model, b. (3.24) for the MR-TBE model. 34 16) Compute ï¨ ï©tï¢ as ï¨ ï© ï¨ ï©ï¨ ï©1cos /mq mt tï¬ ï¬ï­ and correct ï¨ ï©tï¢ï· using ï¨ ï© ï¨ ï© /s st t hhï¢ ï¢ï­ ï­ï© ï¹ï« ï» . 17) Finally, calculate the generated electrical torque eT and calculate the rotor speed ï¨ ï©r tï· and angle ï¨ ï©r tï±using (2.9), a. (3.10), and (3.11) for the MR-TT model, b. (3.26), and (3.27) for the MR-TBE model, respectively. âââEnd of post-network solution subroutine âââ From the above visual separation, we can immediately notice the benefit of implementing the MR algorithm; the implementation skips most of the rotor calculations in the major step, especially during the post-network solution algorithm, the frequency of which is determined by the MR ratio M. The reduced computation, in turn, speeds up the overall time spent by the synchronous machine model per time step. A note on magnetically linear MR CP-VBR SM implementation As the magnetically linear model is a simplification of the above-derived saturable model, a few simplifications of the above implementation would result in the required model. It is assumed that the SM continually operates in the segment defined by j = u or j = 1, significantly reducing the number of coefficients needed to be calculated, as well as forcing ï¨ ï© 2 1resqd t ï´ï½ 0λ . Given no change in the saturation segment takes place, steps 4, 5, and 14-16 in the implantation can be ignored. 35 Chapter 4: MR Model Studies The studies mentioned here are intended to take the reader through an incremental journey of an overall MR model implemented in a state-of-the-art simulator. Firstly, the feasibility and the ideation stages are carried out in MATLAB with a simple machine study using the magnetically linear MR-TT CP-VBR SM model. Given MATLABâs robust computation package and the need to not compile the algorithm, an interpreted software is a great step in the ideation and verification phase of the models. Then, a simple MR-TT magnetically linear CP-VBR SM is designed for PSCAD/EMTDC. The MR model was implemented in C due to its robust preprocessor directives, which provide a better platform to program models for various configurations of the SMs. The model in C is further interfaced with PSCAD using Fortran-C interfacing. The study implemented in this platform consists of a microgrid modelled after a marine microgrid, where the operation of the on-board generators is analyzed using the subject model. Finally, the MR (MR-TT and MR-TBE) CP-VBR SM models with main flux saturation were developed and interfaced in PSCAD in a similar manner to the previous case, with an extended study of the marine grid, now supplying a propulsion motor also represented using the proposed SM. 4.1 Preliminary SMIB Study To demonstrate the performance of the proposed MR CP-VBR model, a single-machine-infinite-bus system with a hydro turbine generator with parameters given in Appendix D.1 is being considered here. To obtain a reference solution, the qd0 model of the same synchronous machine is implemented in Simulink with a fixed time step of 1 sï­ and ODE4 solver. The discretized MR CP-VBR model is programmed with MATLAB with a fast time step at 50 sï­ . The ratio between the slow and fast time steps M is increased to investigate the numerical performance and accuracy. 36 In the following study, the generator is assumed to operate in a steady state with an input torque of 0.5 baseT . At 0.06t ï½ s, a three-phase fault is applied to the terminals, and the simulation is run until 0.2t ï½s. The transient response of the phase a stator current is shown in Figure 4.1. To better see the performance of the MR CP-VBR model with different values of M, the magnified view is also shown in Figure 4.1(b). The transient observed in the rotor speed is plotted in Figure 4.2, and the magnified view is shown in Figure 4.2(b) for better visual comparison. Figures 4.1(b) and 4.2(b) show that the solutions obtained with 5M ï½ and lower are visibly undistinguishable from the reference solution. The solutions obtained with 14M ï½ and higher are starting to visibly deviate, indicating a reduction in simulation accuracy. To further evaluate the accuracy of the proposed MR CP-VBR model under various ratios M, the cumulative norm-2 error [26] (expressed in %) relative to the reference solution trajectory (for the duration of the entire transient) has been calculated and plotted in Figure 4.3. As shown in Figure 4.3, the error generally increases for the fast variables (stator current) and the slow variables (rotor speed) with the increase of ratio M. However, the error increases insignificantly for the solutions obtained with 5M ï½ and lower. Even for the very high ratio 20M ï½ , the error in the stator current remains below 4%, and for the speed below 0.003% (which is insignificant). Based on this evaluation, it can be concluded that with the MR ratio M up to 10, one may obtain reasonable results without sacrificing the numerical accuracy. 37 Figure 4.1 Phase-A current response during a three-phase fault: (a) overall transient, and (b) magnified view. Increasing the MR ratio M is also expected to reduce the computational cost of the model. To investigate the computational gains, the CPU time taken by the proposed MR CP-VBR model for the subject transient fault study has been plotted in Figure 4.4. As shown in Figure 4.4, the CPU time taken by the model is significantly reduced by increasing the ratio M up to 5, and after 10, the reduction becomes less significant. Also, the CPU time can be easily reduced by a factor of 2 or more when the time step ratio M is increased up to 10, with only a slight increase in numerical error. 0 0.05 0.1 0.15 0.2-500501000.1619 0.16195 0.162 0.16205 0.1621-20-19-18-1738 Figure 4.2 Rotor speed response during a three-phase fault: (a) overall transient, and (b) magnified view. Figure 4.3 Cumulative 2-norm error for the stator current and rotor speed trajectories as obtained by the proposed MR CP-VBR model with different raios M. Figure 4.4 CPU time for the considered transient study simulated using different MR ratios M. 39 4.2 PSCAD/EMTDC Implementation of a Microgrid Here, the performance of the proposed MR CP-VBR synchronous machine model is evaluated for the simulation of a shipboard electrical system that includes multiple synchronous machines implemented in PSCAD. The overall system under consideration is shown in Figure 4.5, which consists of five synchronous generators, transformers, full-bridge rectifiers, AC and DC loads, and two propulsion induction motors fed by inverters through the DC buses. In this system, the synchronous generators G1 to G5 (where G3 is an auxiliary generator that may be activated or deactivated depending on the loading condition) supply the AC loads bus through transformers. The generator and system parameters are provided in Appendix D.2. A simplified version of the system depicted in Figure 4.5 has been implemented in PSCAD/EMTDC using the custom-defined library, as illustrated in Figure 4.6. As seen in Figure 4.6, the circuit is simplified for the studies presented in this section by modelling the propulsion motors as a load connected to each DC bus. The proposed MR synchronous machine model has been programmed using language C and interfaced with the rest of the PSCADâs network through the conductance matrix and history injection currents (see [34] for detailed implementation). For the computer studies considered in this section, the reference solution is obtained by running the same model using an SR time step of 1 µs (i.e., hs=hf=1µs, and M=1). The numerical errors converge to almost zero at such a small time step. All studies are done on a computer with an 11th Gen Intel(R) Core(TM) i5-11500 2.70 GHz processor.40 Figure 4.5 Considered shipboard electrical system comprised multiple synchronous generators supplying AC and DC loads. Figure 4.6 Part of the shipâs electrical system depicted in Figure 4.5 implemented in PSCAD/EMTDC using the custom-defined MR synchronous machine model. 41 4.2.1 Time-domain simulations In the following study, the system is assumed to operate in a steady state defined by the real output power of all the generators set to Pg1=Pg2=Pg3=Pg4= Pg5= 5.0 MW, supplying the equivalent AC load with Pac= 13.0MW and the equivalent DC load with Pdc=12.0 MW. The field voltage is set to 2.245 p.u for generators G1, G2, and G3, and the nominal field voltage of 1.123 p.u for generators G4 and G5. It is assumed that at t=0.05s, the auxiliary generator G3 is disconnected by the circuit breakers at the tie lines. At t=0.1s, the auxiliary generator G3 is reconnected. The simulation is run until t=0.5s and is carried out for various timestep sizes of hf , including 1, 10, 50, 100, and 200µs, which corresponds to the MR ratio M of 1, 2, 5, and 10, respectively. The transient responses of the system variables for the synchronous generators G1 and G3 are shown in Figure 4.7, as obtained with different values of hf and M. Figure 4.7 Transient response of generators G1 and G3 variables during the operation of circuit-breakers at the tie lines, opened at t=0.05s and re-closed at t=0.10s, for different values of hf and M: (a) phase-a stator current of G1, (b) phase-a stator current of G3, (c) electrical rotor speed of G1, and (d) electrical rotor speed of G3. 42 Figure 4.8 Zoomed in sections from Figure 4.7: (a) phase-a stator current of G1, (b) phase-a stator current of G3, (c) electrical rotor speed of G1, and (d) electrical rotor speed of G3. For stator transients, the phase-a current of the generator G1 and the auxiliary generator G3 are shown in Figures 4.7(a) and 4.7(b), respectively. For stator transients, the phase-a current of the generator G1 and the auxiliary generator G3 are shown in Figures 4.7(a) and 4.7(b), As seen in Figure 4.7(b), at t=0.05s, the stator current supplied by generator G3 (iasg3) drops to zero when the tie lines are disconnected, i.e., the output power of generator G3 drops to zero. The reduced power following the disconnection of G3 forces the other four generators to supply more power, causing their stator currents to increase. This increase is observed in Figure 4.7(a) for the stator current of generator G1 (i.e., iasg1). When the tie lines circuit breakers are re-closed at t=0.1s, generators G1 and G3 return to their normal operating condition after a short transient period. The magnified views are also provided in Figure 4.8 to show the detailed deviations in magnitude and phase resulting from the implementation of the MR algorithm at various time step sizes and MR ratios. As seen in 43 Figures 4.7(a), 4.7(b), 4.8(a), and 4.8(b), using a small time step with large MR ratios (i.e., hf=1µs and M=10), produces results that are nearly identical to the reference trajectory with hf=1µs and M=1. However, the accuracy remains acceptable even with larger time steps and MR ratios (e.g., hf=200µs and M=10). The same trend can be observed for the transient response of electrical rotor speed 1grï· and 3grï· for generators G1 and G3, respectively, as illustrated in Figures 4.7(c), 4.7(d), 4.8(c), and 4.8(d). As observed, due to the disconnection of G3, and thus the increased load for the remaining generators at t=0.05s, the speed 1grï· reduces while 3grï· it increases. 4.2.2 Numerical accuracy Figure 4.9 Calculated 2-Norm error for: (a) phase-a stator current of G3, and (b) electrical rotor speed of G3, at various fast time steps and MR ratios. To benchmark the numerical accuracy, the cumulative 2-norm error in the trajectory of phase-a stator current iasg3 and rotor speed 3grï· of generator G3 is calculated with respect to the reference solution (SR 1 02.505.001 50 100 150 20000.030.0644 µs time step) and shown in Figs. 5(a) and 5(b) respectively. The errors are calculated for various time step sizes using different MR ratios, i.e., M=1, 2, 5, and 10. As seen in Figs. 5(a) and 5(b), even at large time steps and large ratios, the error is less than 5% for the phase-a stator current and less than 0.05% for the rotor speed, respectively. This proves that the proposed MR model can provide accurate results at large time steps (i.e., hf=200µs) and large MR ratios (i.e., M=10). 4.2.3 Computational efficiency Figure 4.10 Normalized execution time for the custom model in Figure 4.6 for different values of MR ratio. The computational performance of the proposed MR model is analyzed at a typical simulation time step of hf=100µs with MR ratios M to be 1, 2, 5, and 10, where the SR is considered the reference solution. Figure 4.10 illustrates the normalized single machine simulation time for various MR ratio M values with a base of 1.16ms. As shown in Figure 4.10, increasing the MR ratio reduces the machine modelâs simulation time because the slow subsystemâs computations decrease. In particular, even at a ratio of two (M=2), the simulation time decreases by 12%, and a further increase in the ratio to M=5 shows a performance improvement close to 20%, demonstrating a significant improvement in computational efficiency. 4.3 Extended PSCAD/EMTDC Marine Propulsion System The ship micro-grid is considered again, extending the simplified grid in the previous section. A simplified view of the system is given in Figure 4.11, and its implementation in PSCAD is shown in Figure 4.12. Power electronic converters are typically used to drive propulsion motors operating at high frequencies 45 in marine vessels. The propulsion motor considered here is a synchronous machine [27]. The inverters operate at a switching frequency of 1 kHz, and due to the high-frequency components, detailed models require a small timestep for simulation (âð = 10 μs). This scenario presents an ideal case where MR algorithms can enhance even general-purpose simulators targeting high-frequency transients. The specifications for the identical synchronous generators, motors, and the series reactance of the transformers are provided in Appendix B. The generators supply two identical AC loads of 7.5 MW modelled as resistors at nominal voltage, along with the propulsion system. The propulsion system consists of a DC bus established by a full-bridge rectifier, with DC filters provided in Appendix D.3. An inverter, controlled by elementary volts-per-hertz control [28] (open loop control) Figure 4.12, powers the propulsion motors. The steady-state operation of the system is established with an excitation of 1.791 p.u. and a torque input of 0.7813 p.u. or 4 MW for the generators, while the propulsion synchronous machines operate with 2.8781p.u. of excitation and 0.8928 p.u. or 5 MW of load torque. The control set point for the inverter is maintained at 1000 rpm during the steady state. Figure 4.11 Overview of the marine microgrid transient study. 46 Figure 4.12 PSCAD implementation of the marine system shown in Figure 4.11. G1G2G4 G5M1M2G3AC LoadsN : N1 2Rectifier Inverter47 Figure 4.13 Switched stator phase-a voltage of propulsion motor M1 from the marine micro-grid study at an SR: (a) overall transient, (b) pre-command zoomed-in transient, and (c) post-command zoomed-in transient. Figure 4.14 Switched stator phase-a currents of propulsion motor M1 from the marine micro-grid study at an SR: (a) overall transient, (b) pre-command zoomed-in transient, and (c) post-command zoomed-in transient. 48 Figure 4.15 Electrical rotor speed of propulsion motor M1 from the marine microgrid study at an SR. To induce a transient in the system, the command speed to both the propulsion motors is increased to 1200 rpm at t = 0.3 seconds. This forces the entire network into a transient state; the simulation is carried out at hf =10 μs and till t = 2.5 seconds. The switched stator voltage waveform of phase-a of the motor M1 is plotted in Figure 4.13(a), with the steady state on either side of the control represented by Figure 4.13(b) and Figure 4.13(c). From the zoomed-in section, we are able to observe the two-level voltage waveform that is the result of the inverter with sine triangle modulation (STM). Further differences in the fundamental frequency of the switched waveforms can be observed between the zoomed-in sections in Figure 4.13(b) and Figure 4.13(c), corresponding to the two set speeds of 1000 rpm and 1200 rpm. This trend continues for the stator current plot in Figure 4.14 of the motor, except for the reduced magnitude of switching due to the inherent inductance of the machine, also evident from the phase difference between the quantities. Using STM, we expect sinusoidal switched currents to follow the same fundamental frequency as 49 the voltage waveforms; this is observed in both the zoomed-in sections of the current in Figure 4.14(b) and Figure 4.14(c). Finally, Figure 4.15 shows the electrical rotor speed transient starting close to 50 Hz, corresponding to 1000 rpm, and ending close to 60 Hz, corresponding to 1200 rpm. The slope of this increase is set by the custom slew rate block implemented in PSCAD, with a slope of 175. 4.3.1 MR numerical accuracy A zoomed-in section of the current and speed from Figure 4.14 and Figure 4.15, respectively, are plotted in Figure 4.16 as Figure 4.16(a) and Figure 4.16(b). Also, the solutions from the MR-TT and MR-TBE models are plotted for an MR ratio of 15, subject to a slow transient of 150 μs. Visually, both the stator and the rotor transients follow closely to the SR solution taken as reference. For a more accurate numerical analysis, the relative 2-norm error is considered, ï¨ ï© ï¨ ï©ï¨ ï©norm100normreferenarctual nceefere cef fffï­ï½ ï´ï² , (3.38) where f may represent any variable across a specific time period. The error is computed for Phase-A stator current 1masi from the stator subsystem, magnetizing q-axis rotor flux-linkages 1mqmï¬ from the rotor, and the electrical rotor speed 1rmï· from the mechanical subsystem are plotted in Figure 4.17 as Figure 4.17(a), Figure 4.17(b), and Figure 4.17(c) respectively, across various MR ratios. Observing the error trends in Figure 4.17, we can clearly see that the error in the stator transient variables supersedes the error in the rotor and mechanical variables. In fact, the percentage error in the mechanical subsystem is magnitudes smaller than both the rotor and stator subsystem. 50 Figure 4.16 Zoomed in sections from Figure 14 and Figure 15 with M=15 from propulsion motor 1: (a) phase-a current, and (b) electrical rotor speed. Figure 4.17 Relative 2-norm percentage error of variables of propulsion motor 1: (a) phase-a current, (b) q-axis magnetizing flux-linkage, and (c) electrical rotor speed. 1 2 4 6 8 10 12 14 15051015(iasm1 ) (%)1 2 4 6 8 10 12 14 1500.511.5(mqm1 ) (%)1 2 4 6 8 10 12 14 1500.01(rm1 ) (%)MR-TT MR-TB(a)(b)(c)(b) (a) 51 4.3.2 Computational efficiency A basic understanding of C compilers and their functions is assumed for the efficiency analysis. The code is compiled using GCC (x86_64-win32-seh-rev0, Built by MinGW-W64 project) 8.1.0 [29] on an Intel 11th Gen Intel(R) Core(TM) i5-11500@2.7GHz processor, with the -O3 optimization flag for compiler-optimized variant of the code. To measure the efficiency of the compiled models, a per-machine per-step efficiency metric is selected. The internal hardware clock of the processor is accessed via the Windows API [30] to ensure high-precision timing measurements. Additionally, the saturation segment check in the solution algorithm has been disabled, forcing the program to compute a segment change during every step. This approach highlights both the benefit of the algorithm and the upper bound for the computational time per step of the proposed models (worst-case scenarios). Figure 4.18 Normalized average per-step time across multiple MR ratios. The results of the timing studies are summarized in Figure 4.18 as normalized values with respect to a base of the MR-TT model with M=1 (SR operation). We can infer that the trapezoidal-52 only model MR-TT has greater time per step than the mixed integration method MR-TBE across all M. The trend is primarily due to the reduced number of computations carried out in the rotor and mechanical subsystems due to the inherent lower calculations of BE discretization. On the other hand, we can observe the trend of increase in efficiency in both the proposed models as the MR ratio increases. This inference is in line with all the previous literature on MR. In fact, the maximum efficiency plateauing about M = 10 can also be observed here and in the literature. The maximum efficiency in the MR-TT model is about 20%, with about a 5% further increment for the MR-TBE model, proving the superior computational efficiency of MR machine models. 53 Chapter 5: Conclusion 5.1 Contributions This research focuses on developing numerically efficient multirate (MR) models for synchronous machines, specifically for interfacing with EMTP-type simulation programs. Given the recent surge in the use of power electronic converters in modern power systems and shipboard systems, users of these programs are increasingly required to adopt small time steps to achieve accurate simulation results when interfacing with such converters. Although significant contributions have been made to improve the efficiency of synchronous machine models through constant interfacing techniques, the need for small time steps still forces these models to perform redundant calculations for the slow transients of the rotor. The proposed MR CP-VBR model addresses this issue by implementing multirate algorithms that separate the fast and slow dynamics, thereby reducing the computational burden by increasing the sampling period of the slow subsystem. Chapter 2 addresses Objective I of this thesis, which was achieved by generalizing the synchronous machine model to allow access to the neutral node and structuring the equations to accommodate an arbitrary number of damper windings. The resulting model is both application-oriented and flexible, as it formulates the CP-VBR synchronous machine (SM) model in this manner. A simple unbalanced fault case study was conducted to validate the model against the built-in PSCAD synchronous machine model. It was found that the CP-VBR SM performed comparably to the built-in model at small time steps while providing improved results at larger time steps. 54 Objective II is addressed in Chapters 3 and 4, which was accomplished by integrating the MR algorithm into the CP-VBR synchronous machine model initially developed in Chapter 2. The modelâs performance was systematically evaluated through incremental studies across different platforms and complexity levels, starting with a preliminary single-machine-infinite-bus (SMIB) transient fault study implemented in MATLAB. This preliminary stage demonstrated significant computational benefits, including CPU time savings of around 50% in an interpreted environment (MATLAB) while maintaining numerical accuracy with negligible deviations from the reference trajectories for MR ratios up to 10. Following this initial assessment, the MR algorithm was implemented in a more sophisticated compiled environment using PSCAD/EMTDC, where it was applied to simulate a realistic shipboard microgrid power system composed of multiple synchronous generators, transformers, AC/DC loads, and rectifiers. In the PSCAD environment, the MR modelâdeveloped using C and interfaced through Fortranâunderwent further validation by analyzing generator dynamics during load shedding and reconnection events. These simulations confirmed that the numerical accuracy remained within acceptable limits even at larger time step ratios (e.g., M = 10), with errors below 5% for currents and 0.05% for rotor speeds. Appreciable computational time reductions (around 25%â30%) were also achieved. The final study expanded the PSCAD model to simulate an advanced marine propulsion microgrid system that included high-frequency switching inverters and synchronous propulsion motors. This detailed study further affirmed that the MR model delivers substantial computational efficiency improvements, achieving a simulation time reduction of approximately 20% at moderate MR ratios. Importantly, these computational advantages were attained without 55 compromising transient accuracy, demonstrating the MR modelâs suitability for large-scale, complex transient simulations. 5.2 Anticipated Impact With the increasing prevalence of power electronics-based components in terrestrial and marine power systems, average value models [30], [31] may be considered for analysis of dynamics and simulation of transients. However, for detailed simulations, the proposed models of SMs demonstrate strong applicability and readiness for industry deployment. The presented implementation of SM models incorporates advanced features such as configurable damper windings, main flux saturation, and a numerically stable and robust interface with the rest of the network solution in a common EMT simulation program such as PSCAD, enhancing the wide adaptation of the proposed models. The primary contribution of this thesis, the SM MR CP-VBR model, is anticipated to yield a significant impact by reducing both computational load and engineering effort in simulation workflows in simulation studies of modern power grids and electrified ships [33]. Furthermore, the model architecture is conducive to real-time implementation, requiring minimal restructuring.. 5.3 Future Work The research presented in this thesis establishes a robust foundation for efficient modelling and simulation of SMs using MR algorithms interfaced with EMTP-type simulation tools. However, several promising directions remain open for future exploration, two of which are highlighted below: Firstly, the multirate implementation introduced in this work could be extended to other constant-parameter directly-interfaced synchronous machine models, notably the constant-parameter phase-domain (CP-PD) synchronous machine model. Investigating MR approaches 56 within CP-PD [23], [34] synchronous machines may reveal additional computational efficiencies and further enhance numerical performance. Initial studies could explore whether similar or improved computational savings can be achieved, as demonstrated by the MR CP-VBR model. Secondly, considering the increasing trend towards multiphase synchronous machines driven by advancements in power electronics in electric ship applications and renewable energy, the advantages of MR algorithms in simulations present a significant research opportunity. Previous chapters of this thesis have established transparent methodologies for single-rate and multirate simulations of conventional three-phase SMs. Future studies could adopt these methodologies to construct and validate MR-compatible models of multiphase machines, specifically addressing challenges related to modelling accuracy, computational efficiency, and numerical stability across different electrical phase configurations (e.g., 5-phase and 6-phase systems) [29]. Incorporating harmonic flux components could increase the accuracy of SM models, which may be required for high-fidelity simulations of ship propulsion drives to investigate the machine-converter interaction and power quality impacts on the rest of the shipâs electrical loads. 57 Bibliography [1] H. W. Dommel, EMTP Theory Book. Microtran Power System Analysis Corporation, 1992. [2] P. Krause, O. Wasynczuk, S. Sudhoff, and S. Pekarek, Analysis of Electric Machinery and Drive Systems, Third Edition., vol. 3. Wiley-IEEE Press, 2013. [3] 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. [4] âSimulink Dynamic System Simulation Software Users Manual.â MathWorks Inc., 2023. [Online]. Available: www.mathworks.com [5] PLECS User Manual. Plexim GmbH, 2023. [6] OPAL-RT User Documentation. OPAL-RT Technologies, 2023. [7] Typhoon Hardware-in-the-Loop (HIL) User Manual. Typhoon HIL, Inc., 2023. [8] EMTP-RV Software Package. Powersys Solutions, 2023. [9] ATPDraw (Alternative Transient Program) User Manual. EEUG, 2023. [10] RSCAD/RTDS Hardware-in-the-Loop Simulator. RTDS Technologies, Inc., 2023. [11] A. U. M and G. Chen, First Course in Numerical Methods. Society for Industrial and Applied Mathematics, 2011. [12] Y. Song, Y. Chen, S. Huang, Y. Xu, Z. Yu, and J. R. Martí, âFully GPU-based electromagnetic transient simulation considering large-scale control systems for system-level studies,â IET Gener. Transm. Distrib., vol. 11, no. 11, pp. 2840â2851, 2017. [13] L. Wang and J. Jatskevich, âMagnetically-Saturable Voltage-Behind-Reactance Synchronous Machine Model for EMTP-Type Solution,â IEEE Trans. Power Syst., vol. 26, no. 4, pp. 2355â2363, Nov. 2011. [14] E. Mostajeran, N. Amiri, and J. Jatskevich, âConstant-Parameter Voltage-Behind-Reactance Synchronous Machine Models Considering Main Flux Saturation for EMTP-Type Programs,â IEEE Trans. Energy Convers., pp. 1â12, 2023. [15] M. L. Crow and J. G. Chen, âThe multirate method for simulation of power system dynamics,â IEEE Trans. Power Syst., vol. 9, no. 3, pp. 1684â1690, Aug. 1994. [16] S. D. Pekarek et al., âAn efficient multirate Simulation technique for power-electronic-based systems,â IEEE Trans. Power Syst., vol. 19, no. 1, pp. 399â409, Feb. 2004. [17] V. A. Galván, J. R. Martí, H. W. Dommel, J. A. Gutiérrez-Robles, and J. L. Naredo, âMATE multirate modelling of power electronic converters with mixed integration rules,â in 2016 Power Systems Computation Conference (PSCC), Jun. 2016, pp. 1â7. [18] F. Therrien and J. Jatskevich, âMultirate EMTP-Type Induction Machine Models,â IEEE Trans. Energy Convers., vol. 31, no. 3, pp. 1142â1152, Sep. 2016. [19] âIEEE Guide for the Application of Neutral Grounding in Electrical Utility Systems. Part II - Grounding of Synchronous Generator Systems,â IEEE. [20] D. Braun and G. S. Koeppl, âIntermittent Line-to-Ground Faults in Generator Stator Windings and Consequences on Neutral Grounding,â IEEE Trans. Power Deliv., vol. 25, no. 2, pp. 876â881, Apr. 2010. [21] E. Levi, âA unified approach to main flux saturation modelling in D-Q axis models of induction machines,â IEEE Trans. Energy Convers., vol. 10, no. 3, pp. 455â461, Sep. 1995. 58 [22] F. Therrien, L. Wang, M. Chapariha, and J. Jatskevich, âConstant-Parameter Interfacing of Induction Machine Models Considering Main Flux Saturation in EMTP-Type Programs,â IEEE Trans. Energy Convers., vol. 31, no. 1, pp. 12â26, Mar. 2016. [23] E. Mostajeran and J. Jatskevich, âPhase-Domain Synchronous Machine Model with Main Flux Saturation for EMTP-type Solution,â in 2023 22nd International Symposium INFOTEH-JAHORINA (INFOTEH), Mar. 2023, pp. 1â6. [24] J. R. Martí, âLecture from EECE 560,â Department of Electrical and Computer Engineering, University of British Columbia. Accessed: Apr. 16, 2025. [Online]. Available: https://ece.ubc.ca/courses/eece-560/ [25] M. Xiong et al., âParaEMT: An Open Source, Parallelizable, and HPC-Compatible EMT Simulator for Large-Scale IBR-Rich Power Grids,â IEEE Trans. Power Deliv., vol. 39, no. 2, pp. 911â921, Apr. 2024. [26] G. H. Golub and C. F. V. Loan, Matrix computations. The Johns Hopkins University Press, 2013. [27] J. L. Kirtley, A. Banerjee, and S. Englebretson, âMotors for Ship Propulsion,â Proc. IEEE, vol. 103, no. 12, pp. 2320â2332, Dec. 2015. [28] P. Shrawane, âInduction motor open- and closed-loop control using Volt-per-Hertz controller,â in 12th IEEE International Power Electronics Congress, Aug. 2010, pp. 97â101. [29] âGCC, the GNU Compiler Collection - GNU Project.â Accessed: Feb. 03, 2025. [Online]. Available: https://gcc.gnu.org/ [30] âMinGW/MinGW/include/windows.h at master · jonasstrandstedt/MinGW,â GitHub. Accessed: Feb. 03, 2025. [Online]. Available: https://github.com/jonasstrandstedt/MinGW/blob/master/MinGW/include/windows.h [31] S. Ebrahimi, T. Vahabzadeh, and J. Jatskevich, âNumerically Efficient Average-Value Model for Voltage-Source Converters in Nodal-Based Programs,â IEEE Open J. Power Electron., vol. 5, pp. 93â105, 2024. [32] S. Ebrahimi and J. Jatskevich, âAverage-Value Model for Voltage-Source Converters with Direct Interfacing in EMTP-Type Solution,â IEEE Trans. Energy Convers., pp. 1â4, 2022. [33] I. and C. C. Housing, âUlkatcho First Nation to house largest off-grid solar project in Canada.â Accessed: Feb. 01, 2025. [Online]. Available: https://www.canada.ca/en/housing-infrastructure-communities/news/2024/04/ulkatcho-first-nation-to-house-largest-off-grid-solar-project-in-canada.html [34] L. Wang and J. Jatskevich, âA Phase-Domain Synchronous Machine Model With Constant Equivalent Conductance Matrix for EMTP-Type Solution,â IEEE Trans. Energy Convers., vol. 28, no. 1, pp. 191â202, Mar. 2013. [35] âAvK DIG156 | STAMFORD | AvK.â Accessed: Apr. 11, 2025. [Online]. Available: https://www.stamford-avk.com/alternators/avk/dig156?brand=13246&&model=765 59 Appendices Appendix A : Coefficients introduced in Chapter 1 Sub-transient qd-axis magnetizing inductance in (2.2) are defined as 1 1, ,T Tlqr lkq lkqm ldr lfd lkd lkdnL L L L Lï© ï¹ ï© ï¹ï½ ï½ï« ï» ï« ï»l lï ï (A.1) ï¨ ï© ï¨ ï©1 11 11 1 1 1, ,dm ojk kmqj jmq m jlqr md ldrL LL k L kï­ ï­ï½ ï½ï© ï¹ ï© ï¹ï¢ï¢ ï¢ï¢ï½ ï« ï½ ï«ïª ïº ïª ïºïª ïº ï« ï»ï« ï»ï¥ ï¥l l (A.2) Sub-transient qd-axis magnetizing resistance in (2.3) are defined as 1 1, ,T Tqr kq kqm d fd kd kdnrr r r r rï© ï¹ ï© ï¹ï½ ï½ï« ï» ï« ï»r rï ï (A.3) ï¨ ï© ï¨ ï©ï¨ ï© ï¨ ï©ï¨ ï©ï¨ ï©10 1101, ,n on oj jj jmq qr md drq dlqrklk d krkL k L kk kk kï½ ï½ï½ï½ï¹ï¢ï¢ ï¹ï¢ï¢ï¹ ï¹ï½ ï½ïº ïºï» ï»ïº ï»ï»r ra al l (A.4) ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©11 1111, ,j jmq mdq dlqr ldkrn on oj jkkkL Lk kk kï½ ï½ï½ï½ï¹ï¢ï¢ ï¹ï¢ï¢ï¹ ï¹ï½ ï½ïº ïºï» ï»ïº ï»ï»a al l (A.5) ï¨ ï© ï¨ ï©ï¨ ï© ï¨ ï©ï¨ ï©ï¨ ï©22 2 21 11 1, ,n on oj jj jmq qr md drq dlqr l rk kk kqL k L kk kk kï½ ï½ï½ ï½ï¹ ï¹ï¢ï¢ ï¢ï¢ï¹ ï¹ï½ ï½ïº ïºï» ï»ïº ïºï» ï»r ra al l (A.6) 0 1 0 1, ,j j j j j jmT Tq q d dmq dr rï¢ï¢ ï¢ï¢ï© ï¹ ï© ï¹ï½ ï½ï« ï» ï« ï»a a a a (A.7) Coefficients in (2.5) are defined as 0011 10 1 ,Tjm q m m om oTjo m ojj qjd om o dï´ ï´ï´ï´ ï´ ï´ï© ï¹ï­ï© ï¹ ïª ïºï½ ï½ïªï»ï¹ïº ïª ïºï« ï» ï­ïªï©ï´ ï«ï© ï¹ï´ ï« ï»ï» ïºï«1 I 0AaA0 1 a00 A I (A.8) 60 ï¨ ï© 2 0 111 2 0,Tj j jqTq djj j jq drr dTr Tï·ï·ï·ï© ï¹ï© ï¹ ï© ï¹ï« ï» ï« ï»ïª ïºï½ ïª ïºï© ï¹ ï© ï¹ïª ïºï« ï» ï« ï»ï« ï»ï­a A aAa a A (A.9) ï¨ ï©2 , ,mqr mdmq md mqrr mq md mmq mjmqjjdd mdr L SFLL L LL L rL L SF Lï¢ï¢ï¢ï·ï·ï· ï·ï· ï·ï¢ï¢ ï¢ï¢ï© ï¹ï¢ï¢ï­ïª ïºïª ïºï½ ïª ïºï¢ï¢ ï¢ï¢ ï¢ï¢ïª ïºï­ ï«ïª ïºï« ï»A (A.10) 131,m oojj qjdmï´ï´ï© ï¹ï½ ïª ïºï« ï»a 0A0 a (A.11) Coefficients in (2.6) are defined as 11 2 101 01 / 0, ,0 1 /mqq mmdn djj j jjLLï´ï´ï© ï¹ ï© ï¹ï½ ï½ïª ïº ïª ïºï« ï»ï« ï»a 0B B B0 a (A.12) 111 11 11,11j jkq mq kq mqlkq lkq lkq lkqmj jkqm mq kqm mqlkqm lkq lkqm lkqmjqr L r LL L L Lr L r LL L L Lï¢ï¢ ï¢ï¢ï­ ï­ï¢ï¢ ï¢ï¢ï­ï§ï© ï¹ï¦ ï¶ ï¦ ï¶ïª ïºï§ ï· ï§ ï·ï· ï§ ï·ïª ïºï¨ ï¸ ï¨ ï¸ïª ïºï½ ïª ïºïª ïºï¦ ï¶ ï¦ ï¶ïª ïºï§ ï· ï§ ï·ï§ ï· ï§ ï·ïª ïºï¨ ï¸ ï¨ ï¸ï« ï»ï­Cïï ï ïï (A.13) 11 1 11 1 11 111j j jfd md md mdlfd lfd lfd lkd lfd lkdnj j jkd md kd md kd mdlkd l dfd ff lkd lkd lkd lkdnjkdn mdjlkdn lddfdr r rL L LL L L L L Lr L r L r LL L L L L Lr LL Lï¢ï¢ ï¢ï¢ ï¢ï¢ï­ ï­ ï­ï¢ï¢ ï¢ï¢ ï¢ï¢ï­ ï­ ï­ï¢ï¢ï­ ï­ï¦ ï¶ ï¦ ï¶ ï¦ ï¶ï§ ï· ï§ ï· ï§ ï·ï§ ï· ï¨ ï¸ ï¨ ï¸ï¨ ï¸ï¦ ï¶ ï¦ ï¶ ï¦ ï¶ï§ ï· ï§ ï· ï§ ï·ï§ ï·ï½ ï¨ ï¸ ï¨ ï¸ï¨ ï¸ï¦ ï¶ï§ ï·ï§ ï·ï¨ ï¸Cïïï ï ï ï11j jkqn kqnmd mdlkqn lkd lkqn lkdnr rL LL L L Lï© ï¹ïª ïºïª ïºïª ïºïª ïºïª ïºïª ïºïª ïºïª ïºïª ïºï¦ ï¶ ï¦ ï¶ï§ ï· ï§ ï·ïª ïºï¨ ï¸ ï¨ ïºï¢ï¢ ï¢ï­ï»ï¢ï¸ïªï«ï (A.14) 111dmojj qjomï´ï´ï© ï¹ï½ ïª ïºï« ï»C 0C0 C (A.15) 61 Coefficients in the discrete rotor flux-linkage equation (2.19), and (2.20) ï¨ ï© ï¨ ï© 12 1 3 12 2, ,jm o m oj j jh hh hï­ï« ï«ï½ ï« ï½ ï­C I C C I C (A.16) ï¨ ï©4 11 ,m oj jh h ï«ï½ ï«C I C (A.17) Coefficients in the qd discrete sub-transient emf (2.21) are expanded as ï¨ ï© ï¨ ï© ï¨ ï© 11 1 2 1, ,j j j jr rh hï· ï·ï­ï© ï¹ï½ ï« ï»C BD A (A.18) ï¨ ï© ï¨ ï© ï¨ ï©2211 3, ,j j j jr rh hï· ï·ï­ï© ï¹ï½ ï«ï« ï»C AD A (A.19) ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©11 22 23 , , , ,j j j j jr r rh hï¢ ï¢ï· ï· ï· ï· ï·ï­ï© ï¹ï½ ï«ï« ï» AA C BD (A.20) The constant transformation of a 2 2ï´ diagonal qd resistance matrix with equal diagonal entries ï¨ ï©/ 2 / 22 / 2 / 2 .3/ 2 / 2consttfa a aa a a aa a aï­ ï­ï© ï¹ïª ïºï½ ï­ ï­ïª ïºïª ïºï­ ï­ï« ï»R (A.21) 62 Appendix B : Machine and System Parameters for the study in section 1.6 [2] 6 231 1 2 2270 MVA, 18kV, 4, 0.0062762 10 kgm 60 Hz, 0.0019 ,0.1860 , 1.7580 , 1.8540 , 0.0011 , 0.1297 , 0.0457 ,0.0011 , 0.0099 , 0.1151 , 0.6139 , 0.054,ll b sls mq md fd lfd kdlkd kq lkq kq lkqS V P J f rX X X r X rX r X r Xï¦ ï½ ï½ ï½ ï½ ï´ ï½ ï½ ïï½ ï ï½ ï ï½ ï ï½ ï ï½ ï ï½ ïï½ ï ï½ ï ï½ ï ï½ï½ ï1 2 .0.031 3,0. H9549 8 m , 0.3820 mm HH,ngL L L ï½ïï½ ï½ 63 Appendix C : Additional coefficients introduced in Chapter 2 Coefficients corresponding to the sub-transient qd voltages for the heterogeneous MR-TBE model (3.24) ï¨ ï© ï¨ ï© ï¨ ï© 14141, ,j j j jr rh hï· ï·ï­ï© ï¹ï½ ï« ï»C BD A (C.1) ï¨ ï© ï¨ ï© ï¨ ï©5 411 3, ,j j j jr rh hï· ï·ï­ï© ï¹ï½ ï«ï« ï»C AD A (C.2) ï¨ ï© ï¨ ï© ï¨ ï© ï¨ ï©11 2 26 4, , , ,j j j j jr r rh hï¢ ï¢ï· ï· ï· ï· ï·ï­ï© ï¹ï½ ï«ï« ï» AA C BD (C.3) 64 Appendix D : Machine and System Parameters for the studies in Chapter 4 D.1 Hydro turbine generator used in section 4.1[2-Table 5.10-1] 6 22 2325 MVA, 20 kV, 64, 60 Hz, 35.1 10 0.00234 , 0.1478 , 0.0005 ,0.2523 , 0.01675 , 0.1267 , 0.01736 , 0.1970 , 0.4433 ,0.8989 ., s ls fdlfd kq lkq kd lkd mqmdP J Js r X rX r X r X XXï½ ï½ ï ï½ ï ï½ ï ï½ ïï½ ï ï½ ï ï½ ï ï½ ï ï½ ï ï½ ïï½ ï D.2 Marine microgrid system parameters for section 4.2 Generatorâs data1 [35]: 231 16.40 MVA, 6.6 kV, 6, 710kgm 60 Hz, 0.03646 , 0.714656 ,5.3429 , 11.4685 , 0.0102 , 1.2449 , 0.1379 , 1.2834 ,0.6745 , 0.0577 .,ll b s lsmq md fd lfd kd lkdlkq kqS V P J f r XX X r X r XX rï¦ ï½ ï½ ï½ ï½ ï½ ï½ ï ï½ ïï½ ï ï½ ï ï½ ï ï½ ï ï½ï½ï ï½ ïï½ ï ï Network data: Transformers: N1:N2 = 2:1, and ST = 100MVA Equivalent DC Load: Lf = 2 mH, Cf = 5mF, Rdc = 3.375ï. D.3 Extended marine microgrid system parameters for section 4.3 Generatorâs saturation data (peak values) [21]: mi (A) 269.38 566.78 790.97 910.68 996.97 1085.9 1146.3 1273.4 1411.4 1547.3 mï¬ (Wb) 8.195 17.073 23.523 26.786 28.986 31.035 32.249 34.374 36.346 38.167 1 Also used for generators in section 4.3 65 Propulsion Motorâs data [35]: 21 11 15.6 , 0.0125 0.20423.1278 0.0035 0.0473 0.3557 0.36670.0198, 0.1927.MVA, 3.3kV, 6 poles, 1000 rpm, 50 Hz, 710 Js , ,, 0.682, , , , ,s lsmd fd kd lfd lkdkq lkqJ XX SF r r XrXr Xï½ ï½ïï½ ï½ ï½ï½ï½ ï ïï ï½ ï ï ï½ ï½ï½ï Propulsion Motorâs saturation data (peak values) [21]: mi (A) 471.42 991.88 1384.2 1593.7 1744.7 1900.4 2006.1 2228.5 2469.9 2707.8 mï¬ (Wb) 04.917 10.244 14.113 16.071 17.391 18.621 19.349 20.624 21.808 22.900 Network data: Transformer: 0.4356TX ï½ ï. DC Filters: mH, and mF5 =15f fCL ï½ .
Thesis/Dissertation
eng
Electrical and Computer Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Jatskevich, Juri
Multirate constant parameter generalized synchronous machine model for simulation of electromagnetic transients in shipboard power systems
Text