Applied Science, Faculty of
Electrical and Computer Engineering, Department of
DSpace
UBCV
Mostajeran, Erfan
Master of Applied Science - MASc
University of British Columbia
In today’s wider adaptation of renewable energy resources, computer simulations and power system transient studies play a vital and enabling role. For conducting the power system studies, the nodal-analysis-based electromagnetic transient program (EMT or EMTP)-type simulators have been widely accepted and used by researchers and engineers as the industry standard tools, wherein the synchronous machines are commonly represented by the general-purpose lumped-parameter models. The traditional machine models are expressed in qd0 coordinates, which represents a challenge when interfacing with external networks that are represented in abc physical variables and coordinates, often making the machine models the computational bottleneck of computer simulations. This thesis advances the state-of-the-art and develops computationally-efficient and accurate models of synchronous machines for EMTP-type solutions. To enhance the modelling fidelity, the main flux magnetic saturation is considered. Two groups of models are proposed based on voltage-behind-reactance (VBR) and phase-domain (PD) methodologies, respectively. The previously established models of this type possess rotor-position- and saturation-dependent interfacing circuit, which requires costly refactorization of the entire network conductance matrix at each time step. The new VBR and PD models are proposed that have a constant interfacing conductance matrix. Extensive computer studies demonstrate the superiority of the proposed models in terms of computational efficiency compared to other state-of-the-art models. It is envisioned that the new models will become adopted by many standard offline and real-time EMTP simulators.
CONSTANT-PARAMETER DIRECT INTERFACING OF SATURABLE SYNCHRONOUS MACHINE MODELS FOR ELECTROMAGNETIC TRANSIENT PROGRAMS by Erfan Mostajeran B.Sc., Isfahan University of Technology, 2020 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 2023 © Erfan Mostajeran, 2023 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled: Constant-Parameter Direct Interfacing of Saturable Synchronous Machine Models for Electromagnetic Transient Programs submitted by Erfan Mostajeran in partial fulfilment of the requirements for the degree of Master of Applied Science in The Faculty of Graduate and Postdoctoral Studies Examining Committee: Dr. Juri Jatskevich, Professor, Electrical and Computer Engineering, UBC Supervisor Dr. Hermann Dommel, Professors Emeriti, Electrical and Computer Engineering, UBC Supervisory Committee Member Dr. Christine Chen, Associate Professor, Electrical and Computer Engineering, UBC Supervisory Committee Member iii Abstract In today’s wider adaptation of renewable energy resources, computer simulations and power system transient studies play a vital and enabling role. For conducting the power system studies, the nodal-analysis-based electromagnetic transient program (EMT or EMTP)-type simulators have been widely accepted and used by researchers and engineers as the industry standard tools, wherein the synchronous machines are commonly represented by the general-purpose lumped-parameter models. The traditional machine models are expressed in qd0 coordinates, which represents a challenge when interfacing with external networks that are represented in abc physical variables and coordinates, often making the machine models the computational bottleneck of computer simulations. This thesis advances the state-of-the-art and develops computationally-efficient and accurate models of synchronous machines for EMTP-type solutions. To enhance the modelling fidelity, the main flux magnetic saturation is considered. Two groups of models are proposed based on voltage-behind-reactance (VBR) and phase-domain (PD) methodologies, respectively. The previously established models of this type possess rotor-position- and saturation-dependent interfacing circuit, which requires costly refactorization of the entire network conductance matrix at each time step. The new VBR and PD models are proposed that have a constant interfacing conductance matrix. Extensive computer studies demonstrate the superiority of the proposed models in terms of computational efficiency compared to other state-of-the-art models. It is envisioned that the new models will become adopted by many standard offline and real-time EMTP simulators. iv Lay Summary Conventional synchronous machine models and power systems transient studies play a significant role in transitioning the current power grids to future green power grids with renewable energy resources. However, the implementation and interfacing of models of electrical machines in widely-used electromagnetic transient simulation programs represent a computational bottleneck. This thesis proposes new models of synchronous machines that are computationally more efficient compared to the existing state-of-the-art models. The proposed models have been validated through extensive computer studies using single and multi-machine power systems and demonstrated to offer significant computational benefits. It is envisioned that the new models will be adopted by many common offline and real-time electromagnetic transient program-type simulators and thus become useful tools for researchers and engineers-practitioners. v Preface The majority of research results presented in this thesis have been presented and published in conference proceedings and have been or will be submitted to scientific journals for peer review. Unless stated otherwise, I was responsible for deriving the mathematical equations and models, implementing models in computer programs, analyzing the results, as well as writing and editing the first draft of all papers that came out of this thesis. My supervisor, Dr. Juri Jatskevich, has provided vital leadership and guidance through all the aforementioned steps. He also provided constructive comments in the writing and editing of all papers. The other co-authors have also contributed to this research by providing critical feedback and editorial comments, and by revising the manuscripts. The articles resulting from this work and additional contributions by the co-authors are listed below: Chapter 2 is based on the following conference paper: E. Mostajeran, S. Ebrahimi, A. Safavizadeh, and J. Jatskevich, “Efficient simulation of variable-speed diesel-engine generators using constant-parameter voltage-behind-reactance formulation,” in Proc. IEEE Int. IOT. Electron. Mech. (IEMTRONICS), Toronto, ON, Canada, Jun. 1‒4, 2022, pp. 1‒6. Chapter 3 is based on the following journal paper: 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. (under review), pp. 1‒10, submitted on Feb. 22, 2023: TEC-00189-2023. vi Chapter 4 is based on the following conference paper: E. Mostajeran and J. Jatskevich, “Phase-domain synchronous machine model with main flux saturation for EMTP-type solution,” in Proc. IEEE Int. Symp. INFOTEH-JAHORINA (INFOTEH), East Sarajevo, Bosnia and Herzegovina, Mar. 15‒17, 2023, pp. 1‒6. Chapter 5 is based on the following journal article: E. Mostajeran and J. Jatskevich, “Constant-parameter phase-domain synchronous machine models considering main flux saturation for EMTP-type solution,” IEEE Trans. Energy Convers. (to be submitted), pp. 1‒10. vii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary ............................................................................................................................... iv Preface .............................................................................................................................................v Table of Contents ........................................................................................................................ vii List of Tables ................................................................................................................................ xi List of Figures .............................................................................................................................. xii List of Symbols .............................................................................................................................xv List of Abbreviations ................................................................................................................. xix Acknowledgements ......................................................................................................................xx Dedication .................................................................................................................................. xxii Introduction ................................................................................................................1 1.1 Motivation ....................................................................................................................... 1 1.2 Literature Review............................................................................................................ 4 1.3 Research Objectives ........................................................................................................ 7 1.3.1 Objective I: Develop a constant-parameter VBR synchronous machine model with main flux saturation for EMTP-solution ................................................................................. 8 1.3.2 Objective II: Develop a constant-parameter PD synchronous machine model with main flux saturation for EMTP-solution ................................................................................. 9 Modelling Synchronous Machines for EMT Simulations.....................................10 2.1 Model Derivation .......................................................................................................... 12 2.1.1 PD Model .................................................................................................................. 13 2.1.2 qd0 Model ................................................................................................................. 15 viii 2.1.3 VBR Model ............................................................................................................... 18 2.1.4 CPVBR Model .......................................................................................................... 21 2.2 Computer Studies .......................................................................................................... 24 2.2.1 Case Study System .................................................................................................... 25 2.2.2 Transient Simulation ................................................................................................. 27 2.2.3 Assessment of Numerical Performance .................................................................... 30 Constant-Parameter VBR Synchronous Machine Models Considering Main Flux Saturation for EMTP-Type Solution .................................................................................33 3.1 Representation of Magnetic Saturation ......................................................................... 34 3.2 Variable-Parameter SVB VBR Model .......................................................................... 36 3.3 Variable-Parameter EMTP-Type VBR Model ............................................................. 38 3.3.1 Model Derivation ...................................................................................................... 38 3.3.2 Simulation Steps ....................................................................................................... 40 3.4 Constant-Parameter EMTP-Type VBR Models ........................................................... 42 3.4.1 Constant-Parameter Formulation I (CPVBR-I Model) ............................................. 42 3.4.2 Constant-Parameter Formulation II with Current Correction (CPVBR-II Model) ... 46 3.4.3 Simulation Steps ....................................................................................................... 48 3.5 Computational Efficiency Assessment ......................................................................... 50 3.6 Computer Studies .......................................................................................................... 52 3.6.1 Single-Machine-Infinite-Bus (SMIB) System .......................................................... 53 3.6.2 Multimachine System ............................................................................................... 57 Phase-Domain Synchronous Machine Model with Main Flux Saturation for EMTP-Type Solution ...................................................................................................................61 ix 4.1 Model Derivation .......................................................................................................... 62 4.1.1 State-Space Model .................................................................................................... 62 4.1.2 Discrete Model .......................................................................................................... 63 4.1.3 Simulation Steps ....................................................................................................... 64 4.2 Computer Studies .......................................................................................................... 66 4.2.1 Change of Stator Voltage .......................................................................................... 66 4.2.2 Change of Mechanical Torque .................................................................................. 68 4.2.3 Change of Stator Voltage Under Full Load .............................................................. 71 Constant-Parameter PD Synchronous Machine Models Considering Main Flux Saturation for EMTP-Type Solution .........................................................................................73 5.1 Constant-Parameter Formulation .................................................................................. 74 5.1.1 Constant-Parameter Resistance Matrix Formulation ................................................ 74 5.1.2 Formulation of Interfacing Circuit with Constant Conductance Matrix ................... 79 5.2 Augmented Current Correction Scheme ....................................................................... 81 5.3 Efficient Computation of History Terms ...................................................................... 83 5.3.1 Models without Correction Scheme (CPPD-I Models) ............................................ 84 5.3.2 Models with Correction Scheme (CPPD-II Models) ................................................ 85 5.4 Execution Steps ............................................................................................................. 85 5.5 Analysis of Computational Efficiency .......................................................................... 87 5.6 Computer Studies .......................................................................................................... 90 5.6.1 Single-Machine Infinite-Bus System ........................................................................ 91 5.6.2 Multimachine System ............................................................................................... 95 Contributions, Anticipated Impact, and Future Work.........................................99 x 6.1 Contributions................................................................................................................. 99 6.2 Anticipated Impact ...................................................................................................... 103 6.3 Future Work ................................................................................................................ 104 Bibliography ...............................................................................................................................106 Appendix A: Reference Frame Transformation Matrices Used in Chapter 2 .....................115 Appendix B: Parameters for the Case Study of Subsection 2.2.1 ..........................................116 Appendix C: Additional Scalers and Matrices Introduced in Chapter 3 .............................117 Appendix D: Parameters for the Case Study of Subsection 3.6.1, and Sections 4.2 and 5.6......................................................................................................................................................120 Appendix E: Parameters for the Case Study of Subsection 3.6.2 ..........................................121 Appendix F: Additional Scalers and Matrices Introduced in Chapter 4 .............................122 Appendix G: Additional Scalers and Matrices Introduced in Chapter 5 .............................123 xi List of Tables Table 1.1 Classification of synchronous machine models with/without flux saturation based on interfacing with an external network (RPD: rotor-position-dependent; SSD: saturation-segment-dependent) ....................................................................................................................................... 8 Table 2.1 Comparison between the computational performance of subject models .................... 31 Table 3.1 Flops and nonlinear functions count per time step ....................................................... 51 Table 3.2 CPU time per step of the PD, VBR, CPVBR-I and CPVBR-II models ....................... 52 Table 3.3 Average 2-norm error of fault current contributed by the two synchronous generators for different time-step sizes ................................................................................................................ 59 Table 3.4 CPU per time step required for simulation of the IEEE 9-bus system with different models normalized by that of the VBR model.............................................................................. 59 Table 5.1 Genesis of the Interfacing Resistance Matrix of the PD Model in qd0 Coordinates and the Sequential Steps for Constructing Constant-Parameter Resistance Matrix ............................ 77 Table 5.2 Flops and Nonlinear Functions Count Per Time Step .................................................. 88 Table 5.3 CPU Time Per Step of the Models in Table 5.3 ........................................................... 90 Table 5.4 Average 2-Norm Error of Fault Current Contributed by the Set of the Synchronous Generators Connected to Bus 37 for Different Time-Step Sizes Expressed in Percentage (%) ... 97 Table 5.5 CPU Time Required for Simulation of the Modified IEEE 39-Bus System with Different Models for 0.5s Short-Circuit Transient Study ............................................................................. 98 Table 6.1 Classification of synchronous machine models with/without flux saturation based on interfacing with an external network (RPD: rotor-position-dependent; SSD: saturation-segment-dependent) ................................................................................................................................... 103 xii List of Figures Figure 2.1 Equivalent circuit of a conventional synchronous machine expressed in physical variables. ....................................................................................................................................... 12 Figure 2.2 Equivalent circuit of a conventional synchronous machine represented in qd0 coordinates. ................................................................................................................................... 17 Figure 2.3 Implementation of variable-parameter VBR formulation of the synchronous machine model............................................................................................................................................. 21 Figure 2.4 Implementation of the CPVBR model using a low-pass filter Hv (s) to relax the algebraic loop. .............................................................................................................................................. 24 Figure 2.5 Variable-speed DEG-rectifier system consisting of voltage and speed control loops supplying a dc load. ...................................................................................................................... 25 Figure 2.6 Transient response of several variables when a load rejection of 0.25 pu and step down in reference speed occur at t = 0.5s: (a) rectifier dc voltage, (b) dc load power, (c) electromagnetic torque, (d) rotor speed, and (e) stator phase-a current. ................................................................. 28 Figure 2.7 Magnified view during the moment of transient: (a) electromagnetic torque in Figure 2.6(c), and (b) stator phase-a current in Figure 2.6(e). ................................................................. 29 Figure 2.8 Magnified view of variables in Figure 2.7: (a) electromagnetic torque, and (b) stator phase-a current. ............................................................................................................................. 30 Figure 3.1 Piecewise-linear representation of the main flux saturation characteristic. ................ 35 Figure 3.2 Projections of the saturated magnetizing and residual fluxes onto the d- and re-scaled q-axes. ........................................................................................................................................... 36 Figure 3.3 Transient responses of the subject models to a single-phase-to-ground fault: (a) phase a current ias, (b) phase b current ibs, (c) phase c current ics, (d) electromagnetic torque Te, (e) field xiii current ifd, and (f) main flux λm. .................................................................................................... 54 Figure 3.4 Magnified view of variables in Figure 3.3: (a) phase c stator current ics, (b) electromagnetic torque Te, and (c) field current ifd. ...................................................................... 55 Figure 3.5 2-norm errors of stator currents ε (iabcs) and electromagnetic torque ε (Te) as a function of time-step size ∆t for the single-phase-to-ground fault case study. ........................................... 56 Figure 3.6 One-line diagram of the IEEE 3-generator 9-bus system. ........................................... 57 Figure 3.7 Three-phase-to-ground fault study in IEEE 9-bus system with the time-step of 500 µs: (a) phase a current of the primary side of T1 bus 1 (b) magnified view of phase a current of the primary side of T1. ........................................................................................................................ 58 Figure 4.1 Simulation results obtained by the subject models for the stator voltage step-up transient for the selected variables: (a) q-axis stator current qsi , (b) d-axis stator current dsi , and (c) electromagnetic torque eT . ............................................................................................................ 67 Figure 4.2 Reference transient trajectories of the mechanical torque step-up case study obtained for the selected variables: (a) q-axis magnetizing flux λmq , (b) d-axis magnetizing flux λmd , and (c) main flux angle  . ................................................................................................................... 69 Figure 4.3 Simulation results obtained by the subject models for the mechanical torque step-up transient for the selected variables: (a) q-axis stator current qsi , (b) d-axis stator current dsi , (c) electromagnetic torque eT , and (d) rotor angle δ. ........................................................................ 70 Figure 4.4 Simulation results obtained by the subject qd0 and PD models for the step-up stator voltage transient under full load for the stator phase a current asi : (a) overall transient, and (b) magnified view. ............................................................................................................................. 72 xiv Figure 5.1 Norton equivalent representation of the interfacing circuit of synchronous machine model for EMTP-type solution where G is the equivalent conductance matrix, and j is the equivalent injection current vector. ............................................................................................... 81 Figure 5.2 Transient response of several variables to SLG fault and its removal: (a) phase a current ias, (b) phase b current ibs, (c) phase c current ics, (d) electromagnetic torque Te, and (e) main flux λm. .................................................................................................................................................. 92 Figure 5.3 Magnified views from Figure 5.2 obtained for models executed with 1 ms time step: (a) phase c stator current ics, and (b) electromagnetic torque Te. ....................................................... 93 Figure 5.4 The 2-norm error of stator currents ε(iabcs) and electromagnetic torque ε(Te) as a function of time step solver for the subject SLG fault. ............................................................................... 94 Figure 5.5 One-line diagram of the modified IEEE-39 bus system comprising 18 synchronous generators (G1 is modeled as an ideal source). ............................................................................. 95 Figure 5.6 The subject three-phase-to-ground fault study obtained with ∆t = 500μs for phase c fault current contributed by the generators connected to bus 37 iBus37,c: (a) the overall transient and (b) the magnified view. ....................................................................................................................... 96 xv List of Symbols In this thesis, scalers are written using italic fonts (e.g., ias), vectors are denoted by bold lowercase characters (e.g., iabcs), and matrices by bold uppercase characters (e.g., R). For notational compactness, stator variables are often written grouped using the following vector notation:  , ,abcs as bs csf f f=f , 0 0, ,qd s qs ds sf f f =  f , ,qds qs dsf f =  f , where f may represent v, i, and λ (these variables are defined below). The rotor variables are defined as 1 2, , ,qdr kq kq fd kdf f f f =  f , where f may denote v, i, and λ. For a simpler presentation of machine variables, all rotor variables are referred to the stator side by considering the appropriate turns ratio. All other vectors and matrices are defined explicitly throughout the thesis. Only the basic variables are aggregated in this section. All other variables are defined within the thesis. ase , bse , cse Substransient phase a, b, and c voltages qse , dse q-axis and d-axis substransient voltages xfdE Re-scaled field winding voltage asi , bsi , csi Stator phase a, b, and c currents xvi qsi , dsi , 0 si q-axis, d-axis, and zero-sequence stator currents fdi Field winding current lkdi d-axis damper winding current kqzi , 1,2z = q-axis damper winding current I Identity matrix J Combined machine-load moment of inertia hJ Current injection vector Kp, Ki Proportional-integral controller parameters sK , 1s−K Park’s and inverse of Park’s transformation matrix for stator variables lsL Stator leakage inductance Llfd Field winding leakage inductance Llkd The d-axis rotor damper winding leakage inductance lkqzL , 1,2z = The q-axis damper winding leakage inductances mquL , mduL Unsaturated q-axis and d-axis magnetizing inductances mqsL , mdsL Saturated q-axis and d-axis magnetizing inductances p Heaviside’s operator P Number of poles sr Stator resistance kdr The d-axis damper winding resistance kqzr , 1,2z = The q-axis damper winding resistances xvii SF Saliency factor eT Electromagnetic torque mT Mechanical torque t Time asv , bsv , csv Stator phase a, b, and c voltages qsv , dsv , 0 sv q-axis, d-axis, and zero-sequence stator voltages fdv Field winding voltage kdv d-axis damper winding voltage kqzv , 1,2z = q-axis damper winding voltage β Flux vector angle ∆t Simulation step size ( )x 2-norm relative error of x θr Electrical rotor position δ Rotor angle λas , λbs , λ cs Stator phase a, b, and c flux linkages λqs , λds , 0λ s q-axis, d-axis, and zero-sequence stator flux linkages λmq , λmd q-axis and d-axis magnetizing flux linkages λlfd Field winding flux linkage λ lkd d-axis damper winding flux linkage λlkqz , 1,2z = q-axis damper winding flux linkages xviii  Time derivate of flux vector angle ωr Electrical speed of the rotor xix List of Abbreviations ac Alternating current CPU Central processing unit d-axis Direct axis DAE Differential-algebraic equation dc Direct current DER Distributed energy resource EMT Electromagnetic transients EMTP Electromagnetic transients program EMTP-type Nodal-analysis-based programs based on the original EMTP LPF Low-pass filter OCT Open circuit test ODE Ordinary differential equation PD Phase-domain (model) PI Proportional-integral q-axis Quadrature axis SG Synchronous generator SM Synchronous machine SVB State-variable-based (program) VBR Voltage-behind-reactance (model) xx Acknowledgements First and foremost, I would like to express my sincere gratitude to my supervisor Dr. Juri Jatskevich. His visionary and inspirational leadership, profound wisdom, dedication and passion for research, and constructive feedback were essential for realizing this thesis and its contributions. I would like to profoundly thank Dr. Jatskevich for his patience in answering my endless questions and his unlimited support during the most challenging days of my studies at UBC, especially during the COVID-19 pandemic at the beginning of my program. He has taught me to become a better student, a more curious researcher, a more confident engineer, and above all, a better human being. The funding for this research was made through 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” in collaboration with Manitoba Hydro International and RTDS Technologies. I feel very privileged and humbeled having the opportunity to work in this vital area and contribute to the development of state-of-the-art simulation tools used by industries worldwide. I would like to express my appreciation to the examining committee members, Dr. Hermann Dommel and Dr. Christine Chen, for their outstanding dedication and valuable time spent on reading this lengthy thesis and providing invaluable feedback and comment on how to improve it. I am also thankful to the faculty, staff and fellow students at UBC, Department of Electrical and Computer Engineering, who have inspired me to continue my work and graduate study. Specifically, I am deeply thankful to my colleagues and friends in the UBC Electric Power and Energy Systems research group. I want to express my appreciation to Dr. Seyyedmilad xxi Ebrahimi and Dr. Navid Amiri for engaging discussions and for their valuable comments and feedback both in research and editing my papers. I would also like to thank Arash Safavizadeh, Parastoo Hosseinian, Jinhe Zhou, Taleb Vahabzadeh, Sheraz Baig, Levi Bieber, Gloria Wang, and of course, many others for their support and for making my time in our lab fun and enjoyable. Particularly, I want to thank Arash Safavizadeh for his listening to my endless complaints, our many technical discussions regarding the papers we co-authored, and overall for his genuine and deep friendship during all these years at UBC. I also want to acknowledge the support of my dearest friends Arash, Mehdi, Adel, Pouria, and Amir. I am forever thankful to you for always being there for me and helping me go through the highs and downs of my life. I also owe a debt of gratitude to my lovely parents Ramezan and Mozhgan, for their unconditional support throughout all years of my life and studies, both financially and morally. I cannot be thankful enough to them for inspiring and motivating me to always strive for good and be a better version of myself. I am also immensely grateful to my two brothers, Mehdi - for being a role model and helping me navigate the next chapters of my life, and Meisam - for always rooting and caring for me even when we are thousands of miles away. xxii Dedication To my lovely Mom and Dad 1 Introduction 1.1 Motivation Electrical power systems worldwide are undergoing rapid and fundamental changes, where renewable energy resources are replacing conventional fossil-fuel-based generators and power plants. Contrary to traditional power plants, distributed energy resources (DERs) comprise a high share of generation in the future grids and are mostly interfaced with the grid through power electronic converters. This advanced and intricate interfacing tends to alter the long-known characteristics of the traditional grids by introducing fast and complex transients and interactions with the rest of the network. In this regard, running many computer simulations is essential to ensure a reliable, safe, stable, and optimum operation of today’s grids and their transition to future green and sustainable power grids. Consequently, more than ever, the need for developing numerically efficient and accurate computer models of power systems components and respective simulation tools is pronounced. In many power system computer studies, synchronous machines have a dominant presence because of their significant impact on the realization of the smooth transition from fossil-fuel-based to modern sustainable power systems. Specifically, on the one hand, synchronous machines are still the most significant source of electrical energy conversion in many of today’s conventional, bulk, and generation-centralized power systems [1]. For instance, in today’s the United States power grid, they contributed to 78.4% of the overall electrical energy generated in 2022 [2]. On the other hand, they are also mainly responsible for providing inertia and frequency support required in future grids. Consequently, due to their significant presence, efficient modelling of synchronous machines can significantly impact the fast-growing number of 2 conducted power system studies in terms of numerical accuracy and efficiency. The input mechanical energy in synchronous machines is converted to the output electrical energy through a magnetic field [3]. From far below to near nominal operating points, the machine experiences linear to minor-to-moderate saturation [1], [3] yielding a linear to the moderately-nonlinear relationship between flux and current. However, during network disturbances such as faults or overvoltages, the machine may deeply enter the saturation region and have a more pronounced nonlinear relationship between flux and current. This can significantly impact the machine dynamics during transients and steady-state. Therefore, it is necessary to include the flux saturation in the synchronous machine models for ultimate fidelity and consistency with experimental results. Among the available tools for conducting power system studies, the electromagnetic transient (EMT) programs [4], [5] offer great capabilities and fidelity [6] for studying modern and future power grids with fast and complex transients. The EMT simulation programs may be categorized into two groups based on the underlying modelling approach: 1) state-variable-based (SVB) programs [5], [7]‒[9] and 2) nodal-analysis-based EMT programs which are the derivatives of the original EMTP program [4] (thereafter referred to EMTP-type programs) [10]‒[14]. In the SVB simulators, a set of differential-algebraic equations (DAEs) from all component models are formed in state-space form, which are then discretized at the system-level, and solved numerically using variable- or fixed-time-step ordinary differential equation (ODE) solvers [15], [16]. Contrary to the EMTP-type programs, in the SVB packages, the system of ODEs or the equivalent system “matrix” are independent of the time step and discretization method. As a result, numerically efficient variable-time step ODE solvers can be used without the need to run the costly state matrix generation algorithm at each time step [17]. Consequently, the SVB simulation programs have 3 been successfully used for studies of smaller systems with many control elements and/or power electronic converters, where having the ability to use a variable time step may be beneficial. Conversely, in EMTP-type programs, the circuit equations are discretized at the branch/component level using numerical integration rules (e.g., implicit Trapezoidal rule), and the resulting difference equations are then combined to form a set of nodal equations to be solved at each time step [4]. This modelling approach yields a time-step-dependent conductance matrix that is subject to computationally expensive refactorization/re-triangularization whenever the time step or system parameters change. Hence, most EMTP-type programs [10]‒[13] use a fixed time step for the entire simulation interval, resulting in significant computational efficiency and fast computer studies. The high efficiency of EMTP-type simulators makes them an ideal candidate for large-scale practical networks. This thesis mainly focuses on deriving mathematical models of synchronous machines with high numerical performance for their implementation in EMTP-type simulation packages. In the EMTP-type programs, in addition to adopting a fixed time step, further measures are required to realize a constant conductance matrix in more complex networks. Specifically, although having a fixed time step makes the fundamental elements (i.e., inductors and capacitors) to become constant, the existence of saturable and/or position-dependent components (as is the case in rotating electrical machines [3]) yields a time-variant conductance matrix, which necessitates its re-formation and refactorization at each time step. Avoiding frequent refactorization of the entire conductance matrix has been reported to improve the simulation speed by a factor of 2 to 4 in practical case studies [18], [19]. Hence, the computational saving can effectively reduce the CPU utilization and the number of required CPU cores for conducting a given computer study which provides high incentives for deriving models that yield to constant 4 conductance matrix. In EMTP-type simulators, which are the focus of this thesis, the so-called general-purpose models are widely used to represent synchronous machines due to their simplicity, computational efficiency and acceptable accuracy for system-level studies [3]. Although different variations of the general-purpose models have been developed and implemented in the literature and are algebraically-equivalent [20], they show distinctive numerical performance and characteristics [21]. A specific machine model can significantly impact the accuracy, efficiency, stability, and reliability of the entire simulation [21], especially when it comes to interfacing with the external network. Specifically, some models yield a simultaneous solution of the machine-network’s variables compared to other alternatives, which do not and hence, have less accuracy. However, these highly accurate models suffer from time-variant conductance matrix contrary to their less accurate counterparts, which benefit from having a constant conductance matrix. Therefore, there exists a gap for models that can achieve both a constant interfacing conductance matrix and high accuracy, and it is the goal of this thesis. 1.2 Literature Review Three main groups of general-purpose synchronous machines models, which are derived explicitly for EMTP environment, are classical qd0, coupled-circuit phase-domain (PD), and voltage-behind-reactance (VBR) models [21], which will be discussed in great detail in the following. The classical qd0 models have been widely used as the default option for representing synchronous machines in EMTP packages. Although offering simplicity and computational efficiency, the qd0 models suffer from interfacing incompatibility with the external network, 5 which is typically expressed in abc phase coordinates and physical variables [21]. To address the interfacing issue of qd0 models, several solutions have been presented in the literature [4], [21]‒[26]. Specifically, by prediction of several mechanical and electrical variables (such as speed-voltage terms), a direct interface through a Thevenin equivalent resistance matrix and a history voltage source can be formed to enable a simultaneous solution of the machine-network variables [4], [18], [22]. Alternatively, using the values of stator voltages from the previous time step, the machine model can also be indirectly interfaced with the external network using the so-called Norton equivalent current source in parallel with artificial snubbers [23], [24], as available in PSCAD/EMTDC [10]. Additionally, the qd0 models can be interfaced with a reduced-order Thevenin equivalence of the rest of the system, which results in direct interfacing with the external grid [21], [25], [26]. While such solutions address the interfacing incompatibility, they introduce other side-effects such as numerical errors and possibly even instability at large time steps [21]. Therefore, significant attention has been focused on developing interface-compatible machine models that are computationally efficient, numerically stable and accurate even for larger time-step sizes. For improving the numerical accuracy/stability upon the classical qd0 models, the linear phase-domain (PD) models were advocated in [27], [28], wherein the machine equations are expressed and discretized in physical abc coordinates. The linear voltage-behind-reactance (VBR) models were proposed in [18] ,[27], by considering the rotor flux linkages as the state variables (contrary to the selection of rotor currents in PD models) and discretizing the machine differential equations to obtain the discrete models. Both linear PD and VBR models are directly interfaced with the external network through a varying-parameter (i.e., rotor-position-dependent) interfacing circuit. This necessitates the refactorization of the entire network conductance matrix at each time-6 step, and hence, hinders the computational performance of the PD and VBR models [21]. In this regard, constant-parameter PD models have been proposed in [29] and [30], assuming linear magnetic characteristic, and have reported 23% and 27% improved computational efficiency over the linear VBR model [15] in their corresponding multimachine case study, respectively. Incorporating saturation in synchronous machine models has been shown in the literature to enhance the accuracy of system-level studies [31], [32]. In the case of the salient-pole machines, considering magnetic saturation is more challenging since the typical saturation characteristic obtained from the no-load open-circuit test (OCT) is only valid for the magnetic path along the d-axis of the rotor. Hence, the machine saturation has been initially considered only along the d-axis of the rotor [1], [33], as it is the option in PSCAD/EMTDC [10] and RTDS/RSCAD [14]. For greater accuracy, the saturation along the q- and d- axes were advocated to be considered separately [34], [35], as adopted in EMTP-RV [11] and ATP [13]. For better accuracy and consistency with the experimental results [36]‒[39], the main flux saturation (which also accounts for the coupling between the q- and d- axes) can be augmented into the machine models. To achieve this, several single-saturation factor [31], [32], [40] and two-saturation factor [37], [38] approaches have been proposed in the literature. In EMTP-type programs, the saturation characteristic is typically approximated using piecewise-linear segments. Adopting this methodology, the changes in saturation level result in varying incremental inductances in the machine model. Hence, except for the saturable Norton equivalent qd0 model [23], [24], which is interfaced with a time step delay [21], all other saturable qd0 [41], PD [42], and VBR [43] models have variable-parameter interfacing circuit (i.e., rotor-position-dependent and/or saturation-dependent conductance matrices), which requires refactorization/triangulation of the network’s conductance matrix at each time step and increases 7 the computational burden of these models. To enhance the computational performance of the saturable PD and VBR models, partially constant (i.e., rotor-position-independent and saturation-level-dependent) PD and VBR models have been proposed in [44] and [45], respectively. Along each segment of the saturation characteristic, these models have a constant-parameter interfacing conductance matrix, which is achieved by the prediction of the stator currents in q- and d- axes. However, the saturation level changes still necessitate updating the machine interfacing conductance matrix, which leads to deteriorating the numerical performance of these saturable models [44], [45]. 1.3 Research Objectives The ultimate goal of this research is to enhance the numerical efficiency of the EMTP-type simulation packages without compromising their numerical accuracy, stability, robustness, and reliability, which is achieved by introducing numerically-efficient and improved synchronous machine models. Notably, this research holds significant value in light of the fast-growing number of computer studies in the power systems industry, where it can contribute by saving significant computational resources, i.e., reducing the CPU time of many offline studies and improving the utilization of necessary CPU cores in real-time simulators. The properties of the existing general-purpose lumped-parameter EMTP-type synchronous machine models discussed in Section 1.2 are summarized in Table 1.1. As presented in Table 1.1, there is no saturable synchronous machine model with a fully constant-parameter direct interfacing circuit that considers main flux saturation. The other existing state-of-the-art machine models, to some extent, deteriorate the numerical accuracy and/or efficiency of the subject computer studies. 8 Table 1.1 Classification of synchronous machine models with/without flux saturation based on interfacing with an external network (RPD: rotor-position-dependent; SSD: saturation-segment-dependent) Models Interfacing type Interfacing conductance matrix/parameter Magnetically-linear Magnetically-saturable Classical qd0 models [4], [8], [22], [41] Direct Constant Partially Constant (only SSD) [23], [24] Indirect Constant Constant PD models [18] ,[27] Direct Variable (RPD) - [28], [42] Direct Variable (RPD) Variable (RPD/SSD) [29] ,[30] Direct Constant - [44] Direct Constant Partially Constant (only SSD) VBR models [18] ,[27] Direct Variable (RPD) - [43] Direct Variable (RPD) Variable (RPD/SSD) [45] Direct Constant Partially Constant (only SSD) The main focus of this research is to formulate new machine models that include main flux saturation and demonstrate improved numerical performance. This aims to enhance the numerical efficiency of EMTP-type simulators while preserving their numerical accuracy and stability. Therefore, this thesis considers the two most promising state-of-the-art approaches, namely the VBR and PD formulations, as the framework to achieve the main research objectives. 1.3.1 Objective I: Develop a constant-parameter VBR synchronous machine model with main flux saturation for EMTP-solution It is shown in Table 1.1 that all saturable models with direct interfacing capabilities have rotor-position- and/or saturation-segment-dependent interfacing circuit. This thesis’s first objective is to propose new EMTP-type saturable constant-parameter VBR (CPVBR) models for synchronous machines. Specifically, two CPVBR models (CPVBR-I and CVPVBR-II) are derived based on the saturable variable-parameter VBR model proposed in [43], which considers the main flux saturation for predicting both steady state and transients with high accuracy and consistency 9 with experimental results. The new models are formulated to have a true constant-parameter interfacing circuit (i.e., both rotor-position- and saturation-level-independent) which boosts their computational efficiency and makes them desirable for real-time applications. Lastly, a current correction scheme in the CPVBR-I model is further augmented to improve the numerical accuracy, resulting in another version of the proposed model, the CPVBR-II. 1.3.2 Objective II: Develop a constant-parameter PD synchronous machine model with main flux saturation for EMTP-solution As the continuation of the research path set forth in the first objective, several efficient constant-parameter PD synchronous machine models with main flux saturation, referred to as the CPPD-I/II-q, CPPD-I/II-d, and CPPD-I/II-qd, with improved computational efficiency are derived based on the PD formulation. In this regard, the genesis of the interfacing resistance matrix in the variable-parameter PD model is analyzed, based on which, approaches for constructing a constant-parameter interfacing circuit are proposed. Although the CPPD models share the desirable feature of having a constant-parameter interfacing circuit with the CPVBR models, they offer additional computational gain due to their diagonal matrices and more efficient calculation of history terms. This further leads to reducing the number of required mathematical operations to run these models at each simulation time step. 10 Modelling Synchronous Machines for EMT Simulations This section reviews the fundamental modelling of synchronous machines in SVB simulation packages to lay the foundation for deriving the proposed saturable EMTP models in Chapters 3 to 5. To provide the reader with a perspective on the computational gain that can be harvested from having a constant-parameter interfacing circuit, the SVB CPVBR model in [46] is also recalled here. Firstly, the state-space PD model is directly obtained from the magnetically-coupled stator and rotor equivalent windings [3]. Hence, due to being modelled in physical coordinates and variables, the PD model can be directly interfaced with an external network both in rotor and stator subcircuits. However, the inductances of the machine are rotor-position-dependent (time-variable) [3], which limits the PD model’s applicability since most power system simulators do not allow the inclusion of variable-parameter inductors [7]. Even if allowed, the variable-parameter interfacing circuit hinders computational efficiency [15]. Additionally, the size of the inductance matrix in the PD model is another drawback, as it depends on the number of rotor windings and can further increase the computational burden [3]. The classical qd0 model is derived from applying the reference frame theory to the PD model, resulting in constant and decoupled inductances in qd planes [3], [47]. The formulation in qd0 reference frame is highly efficient since variables in fundamental frequency become constant in steady-state. The machine is modelled and interfaced with a noninductive external network by equivalent current sources. However, in the case of an inductive external grid, the machine cannot be directly interfaced with the network due to its incompatible voltage-input-current-output relationship [21]. Augmenting resistive/capacitive snubbers [7] and time-step relaxation [21] are 11 the most common solutions for overcoming the interfacing problem in qd0 models, respectively. However, these approaches create an undesirable trade-off between accuracy and numerical efficiency. Specifically, in the first approach, the higher value of snubbers is chosen, the higher accuracy will be achieved but at the cost of a stiffer system and poorer computational efficiency [48]. Similar logic also applies to the second approach as smaller time steps may be necessary to maintain numerical stability and accuracy [21]. Next, to tighten the gap between qd0 and PD models, the VBR model is derived based on considering the rotor flux linkages, and it directly interfaces the machine in the stator subcircuit [49]. Unlike the PD model, the computational burden of the VBR model is marginally affected by the number of rotor windings, and the size of the machine’s matrix remains unchanged [49]. However, the interfacing circuit of the VBR model is still rotor-position-dependent [49] and prohibits its application in most of the commercially-available SVB simulators [21]. Lastly, to close the gap between the constant-parameter qd0, and variable-parameter VBR and PD models, the constant-parameter VBR (CPVBR) models are proposed, which can result in both direct interfacing in the stator and/or rotor subcircuits [46], [50]. These groups of models are derived by moving the variable part of the interfacing circuit to the so-called subtransient voltage facilitated by using a low-pass filter or a high-pass filter and solving algebraic loops [46], [50]. The CPVBR models are shown to achieve fast computational efficiency without compromising numerical accuracy compared to qd0 models with interfacing snubbers [51], [52]. 12 Figure 2.1 Equivalent circuit of a conventional synchronous machine expressed in physical variables. 2.1 Model Derivation Herein, a conventional three-phase synchronous machine with two damper windings on the q-axis and one damper winding and field winding on the d-axis is considered, as illustrated in Figure 2.1. A linear magnetic circuit is assumed here. The saturable SVB models of VBR and PD are introduced in Chapters 3 and 4, respectively. For consistency with [3], the motoring convention is considered here, meaning the direction of the stator current is into the machine windings, as shown in Figure 2.1. Lastly, throughout this thesis, the synchronous machine, its rotor, and all attached mechanical loads are considered as one rigid body to simplify the presented formulations without losing generality. Consequently, the dynamics of the mechanical subsystem are defined as [3] ( )2r e mPp T TJ = −  , (2.1) rspλ bsrsrspλ asiasicsibsvasvbsvcsvfd                 Stator Rotorpλ cs13 r rp = , (2.2) where the Heaviside’s operator is used for denoting the time derivative, i.e., /p d dt= ; Tm is the mechanical torque; Te is the electromagnetic torque; ωr and θr are the rotor electrical angular speed and position, respectively; J is the combined inertia of the rotor and its mechanical loads; and P is the number of poles. 2.1.1 PD Model Based on Figure 2.1, the voltage and flux equations of the PD model are compactly expressed as [3] 00abcs abcs abcssqdr qdr qdrrp      = +            v i λRv i λR, (2.3) ( ) ( )( )abcs abcss r sr rqdr qdrrs r r     =         λ iL Lλ iL L, (2.4) where Rs and Rr are the stator and rotor resistance matrices given as  , ,s s s sr r r=R , (2.5) 1 2, , ,r kq kq fd kdr r r r =  R , (2.6) where rs is the stator resistance; rkq1 and rkq2 are the resistances of the q-axis damper windings; and rfd and rkd are the resistances of the field winding and d-axis damper winding, respectively. In (2.4), the inductance matrix comprising self and mutual inductances are defined as [3] 14 ( )2 2cos(2 ) cos(2 ) cos(2 )2 3 2 32 2cos(2 ) cos(2 ) cos(2 )2 3 3 22 2cos(2 ) cos(2 ) cos(2 )2 3 2 3A Als A B r B r B rA As r B r ls A B r B rA AB r B r ls A B rL LL L L L LL LL L L L LL LL L L L L           + − − − − − − +   = − − − + − + − −   − − + − − + − −  L , (2.7) ( )cos( ) cos( ) sin( ) sin( )2 2 2 2cos( ) cos( ) sin( ) sin( )3 3 3 32 2 2 2cos( ) cos( ) sin( ) sin( )3 3 3 3mq r mq r md r md rsr r mq r mq r md r md rmq r mq r md r md rL L L LL L L LL L L L                    = − − − −   + + + + L , (2.8) ( ) ( )23Trs r sr r =   L L , (2.9) 120 00 00 00 0lkq mq mqmq lkq mqrlfd md mdmd lkd mdL L LL L LL L LL L L+  + = + + L , (2.10) where ( )s rL is the stator self-inductance matrix; ( )sr rL and ( )rs rL are the mutual-inductance matrix between stator and rotor; rL is the rotor self-inductance matrix; Lls is the stator leakage inductance; Llkq1 and Llkq2 are the leakage inductances of q-axis damper windings; Llfd is the field winding leakage inductance; Llkd is the leakage inductance of d-axis damper winding; Lmq and Lmd are the q- and d-axes magnetizing inductacnes, respectively; and the equivalent inductances in (2.7) are defined as 15 3mq mdAL LL+= and 3md mqBL LL−= . (2.11) Finally, the electromagnetic torque is calculated as [3] ( ) ( )12 2T Te abcs s r abcs abcs sr r qdrr rPT      = +     i L i i L i . (2.12) The PD model derived in (2.3)‒(2.12) has 7×7 inductance matrix and can be directly interfaced with the external network in both stator and rotor subcircuits. The inductance matrices in (2.7)‒(2.9) are time-variable as they are rotor-position-dependent. 2.1.2 qd0 Model Instead of expressing the stator variables in the synchronous reference frame, they can be transferred into the rotor reference frame to make inductance matrices in (2.4) rotor-position-independent. This can be realized by applying Park’s transformation matrix ( )s rK [see Appendix A] into (2.3) and yields [3] λ λqs s qs r ds qsv r i p= + + , (2.13) λ λds s ds r qs dsv r i p= − + , (2.14) 0 0 0λs s s sv r i p= + , (2.15) 1 1 1 1λkq kq kq kqv r i p= + , (2.16) 2 2 2 2λkq kq kq kqv r i p= + , (2.17) λfd fd fd fdv r i p= + , (2.18) λkd kd kd kdv r i p= + , (2.19) 16 where vqs, vds, and v0s are the q-axis, d-axis, and zero-sequence stator voltages, respectively; iqs, ids, and i0s are the q-axis, d-axis, and zero-sequence stator currents, respectively; λqs, λds, and λ0s are the q-axis, d-axis, and zero-sequence stator flux linkages, respectively; vkq1, vkq2, vfd, and vkd are the q-axis damper windings, field winding, and d-axis damper winding voltages, respectively; ikq1, ikq2, ifd, and ikd are the q-axis damper windings, field winding, and d-axis damper winding currents, respectively; λkq1, λkq2, λfd, and λkd are the q-axis damper windings, field winding, and d-axis damper winding flux linkages, respectively. It should be mentioned that the damper windings in the synchronous machine are all short-circuited, and hence, the voltages in (2.16)‒(2.17) and (2.19) are zero. Next, the flux linkages in (2.13)‒(2.19) are defined as [3] λ λqs ls qs mqL i= + , (2.20) λ λds ls ds mdL i= + , (2.21) 0 0λ s ls sL i= , (2.22) 1 1 1λ λkq lkq kq mqL i= + , (2.23) 2 2 2λ λkq lkq kq mqL i= + , (2.24) λ λfd lfd fd mdL i= + , (2.25) λ λkd lkd kd mdL i= + , (2.26) where the magnetizing fluxes are given as ( )1 2λmq mq qs kq kqL i i i= + + , (2.27) ( )λmd md ds fd kdL i i i= + + . (2.28) 17 Figure 2.2 Equivalent circuit of a conventional synchronous machine represented in qd0 coordinates. In (2.20)‒(2.28), imq , imd, λmq, and λmd are the q-axis and d-axis magnetizing currents and flux linkages. The electromagnetic torque in (2.12) can be significantly simplified and calculated as [3] ( )3λ λ4e md qs mq dsPT i i= − . (2.29) Finally, the equivalent circuit of the synchronous machine expressed in qd0 coordinates is represented in Figure 2.2. It should be noted that the formulation in the qd0 model considers stator and rotor flux linkages as state variables, currents as outputs, and stator voltages. rs Lls rλ qsLmdlkdrkdidsvdsω Llfdrfd vfd ifdL ikdrs Lls rλ dsLmqlkq1rkq1iqsvqsω Llkq2rkq2 ikq2L ikq1rs Llsi0sv0s18 2.1.3 VBR Model To set the stage for the CPVBR formulation, the derivation of the variable-parameter VBR model is presented first. The VBR formulation was first introduced in [49] by considering the rotor flux linkages as state variables leading to a current-input-voltage-output model. The full-order VBR model has the following stator interface [49], [53] ( )abcs s abcs abcs r abcs abcp  = + +  v R i L i e , (2.30) where ( )abcs rL is the rotor-position-dependent subtransient inductance matrix and is given as in [53], ( )( )2 22 (2 ) (2 )3 32 2cos(2 ) (2 ) cos(2 )3 32 2cos(2 ) cos(2 ) (2 )3 3s r m r m rabcs r m r s r m rm r m r s rL L LL L LL L L             − +      = − +     + −  L , (2.31) with ( ) ( )coss ls a bL L L L  = + − , (2.32) ( ) ( )cos2am bLL L = − − , (2.33) ( )3mq mdaL LL + = , and ( )3md mqbL LL − = . (2.34) In (2.34), the so-called subtransient magnetizing inductances are defined as 11 21 1 1mqmq lkq lkqLL L L−  = + +   , (2.35) 19 11 1 1mdmd lfd lkdLL L L−  = + +   . (2.36) The subtransient voltages abce are obtained as [53], [46] ( )1 , , 0Tabc s r q de e−   =  e K , (2.37) ( )221λ λ λmq kqzq r d mq kqzz lkqzL reL=  = + −    , (2.38) ( ) ( )2 2λ λ λ λ λmd fdmd kd mdd r q md kd fd md fdlkd lfd lfdL rL r Le vL L L  = − + − + + − , (2.39) where 1s−K is the inverse of Park’s transform defined in Appendix A; the subtransient flux linkages are defined as 1 21 2λ λλkq kqq mqlkq lkqLL L  = +   , (2.40) λ λλfd kdd mdlfd lkdLL L  = +   , (2.41) and the magnetizing fluxes are obtained as λ = λmq mq qs qL i + , (2.42) λ = λmd md ds dL i + . (2.43) The state derivates of the VBR model are represented as [49], [53] ( )11 11λ λ λkqkq kq mqlkqrpL= − − , (2.44) 20 ( )22 22λ λ λkqkq kq mqlkqrpL= − − , (2.45) ( )λ λ λfdfd fd md fdlfdrp vL= − − + , (2.46) ( )λ λ λkdkd kd mdlkdrpL= − − . (2.47) According to (2.30), the VBR model is interfaced with the external network via three time-variable (i.e., rotor-position-dependent) RL branches in series with subtransient voltage sources as illustrated in Figure 2.3. Contrary to the qd0 model, the VBR model does not need artificial snubbers or time-relaxation to be connected to the inductive network. As can be seen from (2.31), the size of the machine inductance matrix is reduced to three from seven in the PD model, which boosts the computational efficiency of the VBR model. However, similar to the PD model, the VBR model still suffers from the varying-parameter interfacing circuit. Lastly, it should be mentioned that all PD, qd0, and VBR models are algebraically equivalent since they can be derived from each other with algebraic manipulation and without any approximation. Hence, the mechanical subsystem of the VBR model is similar to that of PD and qd0 models, and the electromagnetic torque can also be calculated using (2.29) herein. 21 Figure 2.3 Implementation of variable-parameter VBR formulation of the synchronous machine model. 2.1.4 CPVBR Model Multiple versions of the CPVBR models were proposed in [46] to achieve a compatible and constant-parameter interface. The main idea of all such models is to transfer the varying variables in the interfacing circuit to the subtransient voltages (i.e., abce ) and represent the interfacing branches with constant RL elements. As discussed in [46], after moving the varying parameters to the subtransient voltages, the interfacing voltage equation for the decoupled constant-interface circuit [assuming a wye (Y) connection in the stator windings and having access to the neutral] can be written as 0abcs D abcs D abcs ng abcr L p L p = + + +v i i i e , (2.48) rsrsrsExternal networkiasKseabceaebecrL (θ )abcs rRotor State-Space Subsystem(2.31), (2.37)‒(2.3 ), (2.44)‒(2.47) ibsicsKsiqdsvqdsvabcsiabcsω rθrθvasvbsvcsStator Interfacing CircuitvfdL (θ )abcs rn22 where D sr r= , (2.49) D ls mdL L L= + , (2.50) 013mdL L= − , (2.51) and ing is the current of the neutral path of the machine. The subtransient voltages in (2.48) are defined as ( )1 0Tabcs s r q de e−   =  e K , (2.52) ( ) ( )221λ λ λmq kqzq r d q d qs mq kqzz lkqzL re p L L iL=   = + − + −    , (2.53) ( ) ( ) ( )2 2λ λ λ λ λmd fdmd kd mdd r q r q d qs md kd fd md fdlkd lfd lfdL rL r Le L L i vL L L     = − − − + − + + − , (2.54) where q ls mqL L L = + , (2.55) d ls mdL L L = + . (2.56) As noted in (2.53), this formulation needs the value of the derivative of stator current at the present time, which has to be approximated by either a high-pass filter or other numerical methods to prevent the formulation of a nonproper state model [46]. As demonstrated in [46], this approximation leads to a more significant error compared to the other implementations, where alternatively, the current derivative term is firstly obtained based on system inputs and states as [46] 23 ( ) ( )2211λ λ λmq kqzqs qs s qs r d ds r d mq kqzz lkqzqL rpi v r i L iLL =   = − − − − −        . (2.57) Then, the new q-axis subtransient voltage can be written by substituting (2.57) in (2.53) as [46] ( )( ) ( ) ( )221λ λ λmq kqz q ddq r d q d ds mq kqz qs s qsz lkqzq qL r L LLe L L i v r iLL L=    −    = − − + − + −        . (2.58) Using (2.58), the need for approximating the derivate of stator current is eliminated. However, the presence of the stator voltage term vqs in (2.58) forms an algebraic feedthrough in the model. Solving this algebraic loop requires iterative solutions of nonlinear equations and hence increases the computational cost of the model. Relaxation of this algebraic loop can be achieved by using a low-pass filter as [46] ( ) 00vpH ss p=+, (2.59) where ‒p0 is the pole of the filter. The selection of the filter pole is a trade-off between the accuracy and stiffness of the system. As the magnitude of the pole of the low-pass filter increases, the accuracy improves at the cost of making the system stiffer. The systematic approach presented in [54] can be adopted to tune the pole of the filter for the desired accuracy. It should be noted that the rest of the formulations of the CPVBR model are similar to that of variable-parameter VBR, including the mechanical subsystem (2.1)‒(2.2), (2.29) and rotor flux linkages dynamics (2.44)‒(2.47). A block-diagram of the implementation of the CPVBR model is depicted in Figure 2.4. As seen, controlled voltage sources in series with decoupled and constant series resistances and inductances are used for interfacing with the external network. 24 Figure 2.4 Implementation of the CPVBR model using a low-pass filter Hv (s) to relax the algebraic loop. 2.2 Computer Studies Here, the numerical performance of the CPVBR model is compared with the models available in the PLECS library [8]. Presently, PLECS has three built-in models for synchronous machines, i.e., qd0, stator VBR (SVBR) [55], and full VBR (FVBR) models [56]. As stated in [56], the two SVBR and FVBR models are algebraically equivalent, and both result in a variable-parameter interface, which hinders their computational performance. The only difference between these two models is that the field sub-circuit in SVBR needs artificial snubbers to interface with an inductive network, whereas the FVBR model has a direct interface capability both in stator and field windings. rDrDrDExternal networkiasKseabceaebecrRotor State-Space Subsystem(2.44)‒(2.47),(2.52), (2.54), (2.58) ibsicsKsiqdsvds vabcsiabcsω rθ rθvasvbsvcsStator Interfacing CircuitvfdgLDLDLDL0Hvv qs rθingn25 Figure 2.5 Variable-speed DEG-rectifier system consisting of voltage and speed control loops supplying a dc load. 2.2.1 Case Study System For the purpose of evaluating the numerical performance of the CPVBR models, the diesel-engine generator (DEG)-rectifier system shown in Figure 2.5 is implemented in MATLAB/Simulink using the PLECS toolbox. Here, a three-phase salient wound-field synchronous generator (WFSG) is considered, which supplies a resistive load through a six-pulse diode rectifier. All the parameters of the system are summarized in Appendix B. The voltage control loop consists of a PI controller and an actuator (modelled as a transfer function), which regulates the rectifier dc voltage at 1.0 pu by adjusting the voltage applied to the field winding. Also, to limit the voltage ripple introduced by the diode-bridge, a large LC filter is used at the output of the rectifier, as shown in Figure 2.5. To realize the variable speed operation of the generator and obtain the efficient power-speed trajectory, a considerable amount of data on the engine and the generator, including saturation, is needed [57]. Here, without loss of generality, it is assumed that the reference speed is changed from two arbitrary values at the moment of load Efficient power-speedregimeRsnub(only for qd0) rloadLdcCdcvabcsvdcvdcidcabcsiPI(T s+1)11e-T s2 SGω m*ω m PdcLPFPI(T s+1)13vdc*vfdTmTeActuatorEngineActuatorVoltage Control loopSpeed Control Loopisnub26 change. In this regard, the regulation of the rotor speed around the set-point is achieved by the speed control loop consisting of another PI controller and an actuator which is modelled as a transfer function. A simplified diesel engine is also used, similar to [58], which is modelled as a time-delay unit, as shown in Figure 2.5. Three fictitious snubbers are also required across the machine terminal in the qd0 model. As the machine terminal voltages are regulated at the nominal value, the snubbers draw the same amount of current at all operating points. Therefore, as the loading of the machine decreases from the nominal operating point, a larger portion of stator currents iabcs will be formed by snubbers currents isnub. Hence, in different operating points, different values of snubbers are required to maintain the stator current error in an acceptable range (i.e., ≈1%). In this study, the snubbers are chosen to be 350 pu to ensure that the stator current error in the qd0 model stays in the proximity of 1% for the considered loading conditions. Also, the pole of the low-pass filter required for the CPVBR model is chosen to be ‒4300, corresponding to the switching frequency of 360 Hz because of the six-pulse rectification [54]. This way, the solutions error of the CPVBR model remains in the acceptable range of 1%. For consistency, all the subject models are run using MATLAB’s stiff solver ode23tb with maximum and minimum time steps of 1 ms and 0.01 µs, respectively. The absolute and relative errors are chosen to be 10-4. The simulation studies are performed on a personal computer with Intel® CoreTM i7-10750H CPU @ 2.60 GHz processor. Since the FVBR formulation includes no approximations/snubbers in its implementation, it is chosen as the reference model. The reference solution is obtained with the maximum time-step of 1µs using the same ode23tb solver. In similar studies for comparing different models of ac machines [59], [60] in the presence of power-electronic converters, two different approaches are commonly adopted. In one approach, 27 the machine's speed is considered an input to the system, e.g., in wind turbine studies. This way, the same speed is considered for all models, which eliminates the need for solving the mechanical subsystem and excludes the impact of its error on simulation results [59]. The other alternative is the consideration of a grid-connected application where the grid provides a point of reference for the system [60]. These considerations cause all the results from different models to be synchronized, and a fair comparison can be drawn between all implementations by calculating the cumulative 2-norm error [46]. A similar approach is adopted here by inserting the speed profile obtained by the reference model into all subject models and, consequently, eliminating the impact of the error introduced by the mechanical subsystem in the simulation results. This approach is justifiable since the major differences between all the subject models originate from the properties of the electrical subsystem; and neglecting the impact of the error of the mechanical subsystem does not affect the final conclusions about the numerical performance of the subject models. 2.2.2 Transient Simulation For the considered transient study, the system in Figure 2.5 is simulated for 2.5 seconds. The simulation results for the selected variables are provided in Figure 2.6. As shown in Figure 2.6(a)–(d), the system initially operates in the steady state where the generator supplies a load of 0.5 pu at the speed of 0.78 pu while the dc voltage is maintained constant at 1.0 pu. Then at t = 0.5s, a 0.25 pu load rejection occurs [as shown in Figure 2.6(b)], which leads to an immediate 0.08 pu over-voltage across the load, and it is regulated back to 1.0 pu by the control loop, as shown in Figure 2.6(a). At the moment of load change, the reference speed is also stepped down to 0.68 pu, which is assumed to be the fuel-efficient speed according to the new loading condition. 28 Figure 2.6 Transient response of several variables when a load rejection of 0.25 pu and step down in reference speed occur at t = 0.5s: (a) rectifier dc voltage, (b) dc load power, (c) electromagnetic torque, (d) rotor speed, and (e) stator phase-a current. 29 Figure 2.7 Magnified view during the moment of transient: (a) electromagnetic torque in Figure 2.6(c), and (b) stator phase-a current in Figure 2.6(e). The governor adjusts the rotor speed to the new set point according to Figure 2.6(d). As can be seen in Figure 2.6(c), the electromagnetic torque profile follows the same transient as the dc power profile and decreases suddenly at the moment of load rejection. The stator phase-a current is also affected by the decrease in the load, as depicted in Figure 2.6(e). To provide a more detailed snapshot of the considered transient about the moment of load rejection, the magnified views of electromagnetic torque and stator phase-a current are demonstrated in Figure 2.7. To qualitatively compare the results of the subject models with respect to the reference model, further magnetified views of the electromagnetic torque and the stator current are shown in Figure 2.8. It is observed in Figure 2.8 that the VBR models and the CPVBR model predicted the current and torque with great accuracy, contrary to the qd0 model, which has a visibly large error due to the interfacing snubbers. 30 Figure 2.8 Magnified view of variables in Figure 2.7: (a) electromagnetic torque, and (b) stator phase-a current. 2.2.3 Assessment of Numerical Performance For a quantitative comparison between the accuracy of the subject models, the cumulative 2-norm errors of the stator phase-a current and the electromagnetic torque for the considered transient study are calculated and provided in Table 2.1. As observed in Table 2.1, both FVBR and SVBR obtained the torque and current with a very small error (~0.03%), resulting from having a direct interface with the network. The qd0 model has around 1% error both in torque and current due to having fictitious snubbers. Meanwhile, the CPVBR model obtained an acceptable (i.e., less than 1%) error in torque (0.13%) and in stator current (0.15%). 31 Table 2.1 Comparison between the computational performance of subject models Model 2-norm error No. of steps (ode23tbode45) CPU time (ode23tbode45) Torque (Te) Stator current (ias) qd0 1.05% 1.04% 53,918 5.71 s 1,484,733 103.78 s FVBR 0.03% 0.03% 41,928 10.92 s 31,947 8.09 s SVBR 0.02% 0.03% 42,808 6.48 s 28,717 6.22 s CPVBR 0.13% 0.15% 43,102 3.56 s 35,162 3.25 s The comparison of the computational performance between the subject models using ode23tb (stiff) and ode45 (non-stiff) solvers is also included in Table 2.1. As seen in Table 2.1, the CPVBR model achieves the fastest simulations compared to other models. Specifically, when using the stiff solver ode23tb, the CPVBR model allows 38%, 67%, and 45% faster simulation compared to the qd0 (3.56s vs. 5.71s), FVBR (3.56s vs. 10.92s), and SVBR (3.56s vs. 6.48s) models, respectively. It is worth mentioning that, as a result of having constant inductances, the qd0 model was faster than the FVBR and SVBR models, which have variable inductances in their interface, although the qd0 model needs more time steps (i.e., 53,918 vs. 41,928 and 42,808, respectively). It is noted that the number of steps used by the CPVBR model is in the same order as the VBR models (i.e., 43,102). Using the non-stiff solver ode45, the qd0 model takes significantly more time (103.78s) and computational resources (1,484,733 time steps) due to having large snubbers across the machine terminals to maintain acceptable accuracy. Meanwhile, the CPVBR outperforms all other models by achieving 97%, 60%, and 48% efficiency in CPU time compared to the qd0 (3.25s vs. 103.78s), FVBR (3.25s vs. 8.09s), and SVBR models (3.25s 32 vs. 6.22s), respectively. It should be noted that the number of time steps of the two VBR models (31,947 for FVBR and 28,717 for SVBR) and the CPVBR model (35,162) models are significantly fewer than the time steps used by the qd0 model. 33 Constant-Parameter VBR Synchronous Machine Models Considering Main Flux Saturation for EMTP-Type Solution Classical qd0 models are commonly used to represent synchronous machine models in EMTP-type simulators. Due to having incompatible interfacing with the external network, which is expressed in abc coordinates, several solutions have been proposed in the literature [4], [21]‒[24]. Although offering simplicity and low computational cost, all qd0 models are still susceptible to numerical inaccuracy and instability at large time-step sizes [21]. Alternatively, PD and VBR models offer significantly improved numerical accuracy but at the cost of having a variable interfacing circuit (i.e., rotor-position-dependent). This necessitates refactorization of the entire network’s conductance matrix at each time step hindering the computational efficiency of these models. The incorporation of magnetic saturation in EMTP models can be straightforwardly realized by representing the machine’s saturation characteristic by a piecewise-linear function. While improving the modelling fidelity, it makes the already rotor-position-dependent interfacing circuit of the PD and VBR models to be also saturation-segment-dependent, which in turn, further deteriorates the computational performance of these varying-parameter models. Two new EMTP-type saturable constant-parameter VBR (CPVBR) synchronous machines models are proposed in this Chapter. The models, CPVBR-I and CVPVBR-II, are derived based on the saturable variable-parameter VBR model proposed in [43], which considers the main flux saturation. The new models are formulated by moving the rotor-position- and the saturation-dependent portions of the original variable interfacing circuit to the history terms. This formulation 34 results in an actual constant-parameter (i.e., rotor-position- and saturation-independent) interfacing circuit, boosting the computational efficiency of the proposed models compared to other state-of-the-art models. A current correction scheme in the CPVBR-I model is further augmented to improve the numerical accuracy, resulting in another version of the proposed model, CPVBR-II. Lastly, the provided extensive computer studies on single-machine and multimachine systems demonstrate the outstanding computational superiority of the proposed CPVBR models over PD and VBR models. It has also been shown that the CPVBR models maintain a similar numerical accuracy to PD and VBR models even at large time steps. 3.1 Representation of Magnetic Saturation Herein, it is assumed that the d-axis saturation characteristic can be obtained from a typical open-circuit test (OCT). Since the q-axis saturation characteristic is often not available, the single-saturation factor approach [31], [32] is considered herein. This approach incorporates both the q-axis saturation and cross saturation based only on the available d-axis saturation curve, which is achieved by re-scaling the q-axis and converting the anisotropic salient-pole machine into an equivalent isotropic machine as 22λλ λmqm mdSF = +  , (3.1) ( )22m mq mdi SFi i= + , (3.2) where λmq, λmd, imq, and imd are the q- and d-axis magnetizing fluxes and currents, respectively. In (3.1)‒(3.2), SF is the saliency scaling factor applied to q-axis. Considering qd axes saturate to the 35 Figure 3.1 Piecewise-linear representation of the main flux saturation characteristic. same level [31], [32], SF is defined based on the unsaturated and saturated magnetizing inductances (denoted here by Lmqu, Lmdu, Lmqs, Lmds, respectively) and is computed as mqu mqsmdu mdsL LSFL L= = . (3.3) For computational efficiency in EMTP-type programs, it is preferred to use a piecewise-linear function with several segments to represent the magnetic saturation characteristic as illustrated in Figure 3.1 [43]‒[45]. Accordingly, the magnetizing flux is expressed as , ,λ λm D j m res jL i= + , (3.4) where LD,j and λres,j are the dynamic/incremental inductance and the residual flux of the jth segment, respectively, as shown in Figure 3.1. Combing (3.1)–(3.4), the projections of the main flux onto the qd axes are illustrated in Figure 3.2 and are calculated as [43] ( ),λ λmq resqD j mqL SFiSF SF= + , (3.5) ,λ λmd D j md resdL i= + , (3.6) λ λ λ res,(j+1)res,jres,1res,2λ λ mimLD,jLD,(j+1)L D,2L D,136 Figure 3.2 Projections of the saturated magnetizing and residual fluxes onto the d- and re-scaled q-axes. ( ), ,λλ λ cos λλmqresq res j res jmSFSF= = , (3.7) ( ), ,λλ λ sin λλmdresd res j res jm= = , (3.8) where β is the angle of the main flux vector with respect to the q-axis; and λresq and λresd are the projections of residual flux λres,j onto the q- and d-axis, respectively. 3.2 Variable-Parameter SVB VBR Model To set the stage for the derivation of the proposed CPVBR models, the saturable state-space VBR model with four rotor windings (i.e., two damper windings on the q-axis; and one field and one damper winding on the d-axis) and main flux saturation is briefly reviewed [43], [61]. The machine-network interfacing can be formulated as [43] ( ),abcs s abcs abcs j r abcs abcsr p   = + + v i L i e , (3.9) where the definition of ( ),abcs j rL can be found in [43], and rs is the stator resistance. The q- and d-axis substransient magnetizing inductances of the jth segment (required for calculating resqλ SFLD,j (SFi )mqresdλ mqλ mdλ mλ LD,ji mdβ SF37 ( ),abcs j rL [43]), ,mq jL and ,md jL , are defined as [43] 1, 2, 1 21 1 1mq jD j lkq lkqLL SF L L−  = + +   , (3.10) 1,,1 1 1md jD j lfd lkdLL L L−  = + +   , (3.11) where Llkq1, Llkq2, Llfd, and Llkd are the leakage inductances of q-axis damper windings, field winding, and d-axis damper winding, respectively. The subtransient voltages in (3.9) are defined as ( ) ( )1 0TTabcs s r qds−   =  e K e , (3.12) with 1, 2, 3,qds j qds j qdr j resqd = + +e A i A λ A λ , (3.13) λ , λTresqd resq resd =  λ , (3.14) and 1s−K is the inverse of Park’s transformation matrix [3]. The rotor flux state equation is [43] 1, 2, 3, 4qdr j qds j qdr j resqd fdp v= + + +λ B i B λ B λ B , (3.15) where vfd is the rotor field voltage input. Next, the magnetizing flux linkages determining the saturation level are formulated as 1 2, 21 2 ,λ λ λλkq kq resqmq mq j qsLlkq Llkq D jL iL L L SF = + + +   , (3.16) 38 ,,λ λ λλfd kd resdmd md j dsLlfd Llkd D jL iL L L = + + +   . (3.17) The coefficient terms A1,j, A2,j, A3,j, B1,j, B2,j, B3,j, and B4 are defined in (C.1)–(C.13) in Appendix C. It should be mentioned that the derivative of the main flux angle (denoted by ωβ (t)) is required in A3,j and has to be numerically approximated at each time step [46], [61]. 3.3 Variable-Parameter EMTP-Type VBR Model 3.3.1 Model Derivation For EMTP-type programs, the state-space equations are discretized by applying the implicit Trapezoidal integral rule with a fixed step size of ∆t. The time dependency of matrices and scalers (which comes from the dependence on the rotor position and/or speed) is denoted by “(t)” and “(t ‒ ∆t)” for the current and previous time steps, respectively. The saturation-segment dependency is denoted by the subscript “j”. Applying the Trapezoidal rule to (3.9), the discretized machine interface equation at the current time step t is found as [43] ( ) ( )( ) ( ) ( ) ( )3 3 ,abcs s abcs j abcs abcs sht r k t t t t  = + + +v I L i e e , (3.18) where k = 2/∆t, and the stator equivalent history voltage source is given as ( ) ( )( ) ( ) ( ) ( )3 3 ,sh s abcs j abcs abcs abcst r k t t t t t t t t  = − − − + − − −e I L i e v . (3.19) Discretizing the rotor state-space equation (3.15) yields ( ) ( ) ( )( ) ( ) ( ) ( )1, 2, 3, 4,qdr j qds qds j qdr j resqd j fdt t t t t t t v t= + − + − + +λ C i i C λ C λ C . (3.20) The definitions C1,j, C2,j, C3,j, and C4,j are provided in (C.14)–(C.16) in Appendix C. The 39 subtransient voltages (3.13) are then re-formulated by substituting (3.20) into (3.13) as [43] ( ) ( ) ( ) ( )qds j qds qdrt t t t = +e D i h , (3.21) where the history term vector hqdr (t) is given as ( ) ( ) ( ) ( ) ( )1, 2,qdr j qds j qdrt t t t t t t= − + −h E i E λ ( ) ( ) ( ) ( )3, 4,j resqd j fdt t t v t+ +E λ E , (3.22) with definitions of the terms Dj, E1,j, E2,j, E3,j, and E4,j are provided in (C.17)–(C.19) in Appendix C. Converting (3.21)–(3.22) into abc coordinates, we have [43] ( ) ( ) ( ) ( )abcs j abcs rt t t t = +e M i e , (3.23) where ( ) ( )( )( )2 111 2 0jj s stt t t− =   D 0M K K0, (3.24) and ( ) ( ) ( )1 [ 0]T Tr s qdrt t t−=e K h . (3.25) Substituting (3.23) into (3.18), the interfacing equation becomes ( ) ( ) ( ) ( )eq,vbr vbrabcs j abcs ht t t t= +v R i e , (3.26) where the equivalent resistance matrix ( )eq,vbrj tR is defined as ( ) ( ) ( )eq, 3 3 ,vbrj s abcs j jt r k t t = + +R I L M . (3.27) In (3.26), the three-phase history voltage source ( )vbrh te is given by ( ) ( ) ( )vbrh r sht t t= +e e e . (3.28) 40 The existence of both saturation-segment- and rotor-position-dependent matrices ( )j tM and ( ),abcs j tL in (3.27) leads to a model with the variable-parameter interfacing circuit. Finally, since all discrete models are interfaced with the system solver as a Norton equivalence, the Norton equivalent of (3.26) for the VBR model can be expressed as ( ) ( ) ( ) ( )eq, eqvbr vbrabcs j abcst t t t= −i G v j , (3.29) with the Norton admittance as ( ) ( )1eq, eq,vbr vbrj jt t− =  G R , and the injection currents as ( ) ( ) ( )eq eq,vbr vbr vbrj ht t t=j G e . (3.30) 3.3.2 Simulation Steps The necessary steps to simulate the VBR model given the known solution at the previous time step t‒∆t are as follows: 1) Predict the rotor speed ωr (t), the main flux angle derivative ωβ (t), and the magnetizing flux linkages λmq (t) and λmd (t) using linear two-point extrapolation as [4] ( ) ( ) ( )2 2f t f t t f t t= − − −  . (3.31) 2) Calculate the electrical rotor angle θr (t) based on the trapezoidal integration of the predicted rotor speed ( )r t using ( ) ( ) ( ) ( )( )1r r r rt t t k t t t   −= − + + − . (3.32) 3) Compute the main flux λm (t) using (3.1) based on the predicted values of ( )λmq t and ( )λmd t . 41 4) Calculate the residual fluxes λresqd (t) using (3.7), (3.8), ( )λm t , ( )λmq t and ( )λmd t . 5) Compute hqdr (t) using (3.22), calculate er (t) using (3.25), evaluate esh (t) using (3.19), and combine all to form ( )vbrh te using (3.28). 6) Evaluate the Thevenin resistance matrix ( )eq,vbrj tR using (3.27) and calculate its inverse ( )eq,vbrj tG to be put into the network’s conductance matrix. 7) Compute the equivalent Norton current source ( )eqvbr tj using (3.30). Then, the calculated injection current source is inserted into the system solver. 8) Obtain the solution of the network (i.e., solving the network for unknown voltages). 9) Receiving the stator voltages vabcs (t) from the system solver, compute stator currents iabcs(t) using (3.29) and convert them to qd coordinates. 10) Compute the rotor fluxes λqdr (t), the magnetizing flux linkages λmq (t) and λmd (t), and the main flux λm (t) using (3.20), (3.16), (3.17), and (3.1), respectively. 11) Calculate the main flux angle β (t) using either (3.7) or (3.8) and its approximated time derivative ωβ (t) as follows ( )( ) ( ) ( )3 4 22t t t t ttt  − − + − =. (3.33) 12) Compare λm (t) with the limits of the current saturation segment. If required, update the saturation segment/level and the relevant variables as specified in step 13. Otherwise, skip to step 14. 13) Update j and re-calculate λresqd (t), hqdr (t), er (t), λqdr (t), λmq (t), λmd (t), λm (t), β (t), ωβ (t), ( )j tM , and ( ),abcs j tL , respectively. 42 14) Calculate ( )abcs te using (3.23). 15) Compute the electromagnetic torque Te (t) using (2.33), the electrical rotor speed ωr (t) and position θr (t) by applying the Trapezoidal integral rule to mechanical subsystem dynamics (2.5)‒(2.6) as ( ) ( ) ( ) ( ) ( ) ( ) ( )( )1 0.5r r e e m mt t t k P J T t T t t T t T t t −= − +  + − − − − , (3.34) ( ) ( ) ( ) ( )( )1r r r rt t t k t t t   −= − + + − . (3.35) 3.4 Constant-Parameter EMTP-Type VBR Models 3.4.1 Constant-Parameter Formulation I (CPVBR-I Model) To avoid the refactorization/retriangulation of the network’s conductance matrix, the equivalent resistance matrix in (3.27) should first become rotor-position- and then saturation-segment-independent. In this regard, the subtransient inductance matrix ( ),abcs j tL in (3.27) can be reformulated as [27], [43] ( ) ( ) ( ),1, ,0 00 00 0mq j lsabcs j s md j ls slsL Lt t L L tL− +   = +   L K K . (3.36) Substituting (3.36) and (3.24) into (3.27), the resistance matrix ( )eq,vbrj tR can be rewritten as ( ) ( ) ( ) ( )1 0eq, eq,vbr qdj s j st t t t−=R K R K , (3.37) where ( )0eq,qdj tR is composed of a constant matrix and a rotor-speed- and saturation-segment-43 dependent matrix ( )0var,qdj tR as ( ) ( ) ( )0 0eq, 3 3 var,qd qdj s ls jt r kL t= + +R I R , (3.38) ( )( )( ), 11, 12,0var, 21, , 22,000 0 0mq j j j rqdj j r md j jkL D D tt D t kL D +  = +   R . (3.39) As explained in [4], [22], to eliminate the rotor-position dependency, the matrix (3.39) should be reformulated to be diagonal where the first and second diagonal entries are equal. This way, the resulting matrix in abc coordinates becomes rotor-position-independent when substituted in (3.37). As discussed in [4], [22] and later adopted in [44] and [45], averaging the diagonal entries in (3.39) for making a constant interfacing matrix leads to marginal numerical error compared to other possible alternatives [30]. Following the same logic here, the variable resistance matrix ( )0var,qdj tR in (3.39) can be reformulated as ( ) ( )( )( )12,0var, 21,0diag , , 0 00 0 0qj j rqd dj j j j r jR D tt R R D t R   = +     R , (3.40) with ( ) ( ), , 11, 22,0.5 0.5j mq j md j j jR k L L D D = + + + , (3.41) ( ) ( ), , 11, 22,0.5 0.5j jq d mq j md j j jR R k L L D D = − = − + − . (3.42) Consequently, the term ( )diag , , 0j jR R   in (3.40) leads to a rotor-position-independent resistance matrix when substituted into (3.37). However, the resulting matrix is still saturation-segment-dependent, similar to the partially-constant models in [44] and [45]. 44 To overcome the saturation-level dependency, as argued for the CPVBR model of a squirrel-cage induction machine in [62], the variation of the dynamic inductance LD,j has negligible impact on the values of ,mq jL and ,md jL defined in (3.10) and (3.11), respectively. This is justifiable because the values of magnetizing inductances are significantly larger than the values of leakage inductances [62]. Since ,mq jL and ,md jL are the only saturation-dependent terms in ( )0var,qdj tR , then it can be concluded that saturation level changes have a marginal effect on the magnitude of ( )0var, .qdj tR Based on this, ( )0var,qdj tR is reformulated as ( )  ( )( )( )12,0var, 1 1 21,0diag , , 0 00 0 0qj j rqd dj j r jR D tt R R D t R  = +    R , (3.43) where R1 is the evaluation of (3.41) for the unsaturated segment at j = 1; and the definitions of qjR and djR are given as 1q qj j jR R R R = + − , 1d dj j jR R R R = + − . (3.44) Substituting (3.43) in (3.38), the resulting matrix in (3.37) is formulated as ( ) ( ) ( ) ( )1 0eq, eq-const var,vbr cpvbr qdj s j st t t t−= + R R K R K , (3.45) with 2 1 2 1 211 2 2 1 2eq-const 11 2 1 2 212 2 2 2 223 2 2 2 2 22 2 2 2 2a a a a aacpvbr a a a a aaa a a a aaR R R R RRR R R R RRR R R R RR + − + − +   = − + + − +   − + − + +  R , (3.46) 45 ( )( )( )12,0var, 21,000 0 0qj j rqd dj j r jR D tt D t R   =    R , (3.47) 1 1a s lsR r kL R= + + , 2a s lsR r kL= + . (3.48) After insertion of (3.45) into (3.26), eq-constcpvbrR can be isolated by moving the variable part (i.e., saturation-segment and rotor-position/speed-dependent term ( ) ( ) ( )1 0varqds j st t t− K R K ) into the history term voltage source by prediction of stator current at each time step [62]. As a result, the interfacing equation of the first constant-parameter VBR model (hereafter referred to as CPVBR-I) is obtained as ( ) ( ) ( )eq-constcpvbr cpvbrabcs abcs ht t t= +v R i e , (3.49) where ( ) ( ) ( ) ( ) ( ) ( )1 0var,cpvbr vbr qdh h s j s abcst t t t t t−= + e e K R K i . (3.50) In (3.50), the stator current prediction ( )abcs ti is used since the value at the current time step is yet to be computed by the solver but is required in (3.50) [1], [7], [33]. Thereafter, the symbol “⁓” denotes the predicted values. In this regard, a two-point linear prediction is adopted as ( ) ( ) ( )0 0 02 2qd s qd s qd st t t t t= − − − i i i , (3.51) which, according to [4], gives comparable results to the three-point prediction considered in [30], [44], [45], [62]. The CPVBR-I model can then be implemented using the Norton equivalent of (3.49), which is defined as 46 ( ) ( ) ( )eq-const eqcpvbr cpvbrabcs abcst t t= −i G v j , (3.52) where 1eq-const eq-constcpvbr cpvbr− =  G R , and the injection current source is ( ) ( )eq eq-constcpvbr cpvbr cpvbrht t=j G e . (3.53) To enhance the computational efficiency, the voltage history term in (3.50) can be further simplified as ( ) ( ) ( ) ( ) ( )1 0var, 0cpvbr vbr qdh h s j qd st t t t t−= + e e K R i , (3.54) eliminating the need to convert the predicted qd0 variables to abc coordinates, contrary to the formulation proposed in [62], which predicts the stator current in qd0 coordinates and then converts it back to abc variables. Furthermore, as mentioned in [62], the stator history term ( )sh te in (3.19) should be reformulated to reflect the actual operating point from previous time step as ( ) ( ) ( )h-constsh abcs abcst t t t t= − − −e R i v ( ) ( ) ( ) ( )1 0-var, 0qds h j qd s rt t t t t t t t−+ −  − − + −K R i e , (3.55) where the definitions of -consthR and 0-var,qdh jR are provided in (C.20)–(C.23) in Appendix C. 3.4.2 Constant-Parameter Formulation II with Current Correction (CPVBR-II Model) To enhance the accuracy of the CPVBR-I model, the impact of stator currents prediction (3.51) can be alleviated by considering the variable conductance matrix ( )eq,vbrj tG for calculating the stator currents instead of using the constant conductance matrix eq-constcpvbrG that has been inserted 47 into the network solver [30]. Specifically, the stator currents in (3.29) can be calculated using (3.30) and (3.54) as ( ) ( ) ( ) ( ) ( ) ( ) ( )( )1 0eq, var, 0vbr cpvbr qdabcs j abcs h s j qd st t t t t t t−= − + i G v e K R i . (3.56) This way, the stator currents are obtained without prediction error [30]. Since the direct calculation of ( ),abcs j tL and ( )j tM and the inversion of ( )eq,vbrj tR required to calculate ( )eq,vbrj tG at each time step is costly, ( )eq,vbrj tG should be analytically calculated and simplified as much as possible in advance. Specifically, the variable matrix ( )eq,vbrj tG can be calculated using (3.37) and the property of inversion of the product of matrices as ( ) ( ) ( ) ( )1 0eq, eq,vbr qdj s j st t t t−=G K G K , (3.57) where ( )0eq,qdj tG is the inversion of ( )0eq,qdj tR and is found as ( )( )( )1, 12,0eq, 21, 2,3det det 0det det 00 0j j rqdj j r jg D tt D t gg−  = −   G , (3.58) ( )21, 2, 12, 21,det j j j j rg g D D t= − , (3.59) with definitions of 1, jg , 2, jg , and 3g provided in (C.24)–(C.25) in Appendix C. The saturation-segment-dependent terms in (3.58) and (3.59) can be pre-computed for all segments in Figure 3.1 to improve the overall computational efficiency. Therefore, the costly inversion operation is entirely avoided, contrary to the original formulation in [30] for the linear constant-parameter PD model. After that, this CPVBR model with stator current correction is referred to as CPVBR-II. 48 Lastly, since the variable conductance matrix ( )eq,vbrj tG is used for calculating the stator currents in the CPVBR-II model, the VBR stator history term (3.19) should be used here. For boosting computational efficiency, the stator history term (3.19) can be reformulated similarly to (3.55) as ( ) ( ) ( )h-constsh abcs abcst t t t t= − − −e R i v ( ) ( ) ( ) ( )1 0-var, 0qds h j qd s rt t t t t t t t−+ −  − − + −K R i e . (3.60) 3.4.3 Simulation Steps During initialization, all constant and saturation-segment-dependent matrices and scalers are pre-calculated. Considering this stage as step zero, the steps for the simulation of the CPVBR models are as follows [62]: 1) Predict the rotor speed ωr (t), the main flux angle derivative ωβ (t), and the magnetizing flux linkages λmq (t) and λmd (t) using linear two-point extrapolation, i.e., using (3.31) [4]. 2) Calculate the rotor angle θr (t) based on the trapezoidal integration of the predicted rotor speed ( )r t using (3.32). 3) Predict the stator currents ( )0qd s ti using (3.51). 4) Compute the main flux λm (t) using (3.1) based on the predicted values of ( )λmq t and ( )λmd t . 5) Calculate the residual fluxes λresqd (t) using (3.7), (3.8), ( )λm t , ( )λmq t and ( )λmd t . 6) Compute hqdr (t) using (3.22), calculate er (t) using (3.25), and evaluate esh (t), 49 a. using (3.55) for the CPVBR-I model, b. using (3.60) for the CPVBR-II model. Next, combine all history terms to form ( )vbrh te using (3.28). 7) Compute ( )0var,qdj tR using (3.47) and ( )r t and form the history term ( )cpvbrh te using (3.54). 8) Compute the equivalent Norton current source ( )eqcpvbr tj using (3.53). Then, the calculated injection current source is inserted into the system solver. 9) Obtain the solution of the network (i.e., solving the network for unknown voltages). 10) Receiving the stator voltages vabcs (t) from the system solver, compute stator currents iabcs(t) a. using (3.52) in the CPVBR-I model, b. using (3.56) in the CPVBR-II model. 11) Compute the rotor fluxes λqdr (t), the magnetizing flux linkages λmq (t) and λmd (t), and the main flux λm (t) using (3.20), (3.16), (3.17), and (3.1), respectively. 12) Calculate the main flux angle β (t) using either (3.7) or (3.8) and its approximated time derivative ωβ (t) using (3.33). 13) Compare λm (t) with the limits of the current saturation segment. If required, update the saturation segment/level and the relevant variables as specified in step 14). Otherwise, skip to step 15). 14) Update j and re-calculate λresqd (t), hqdr (t), er (t), λqdr (t), λmq (t), λmd (t), λm (t), β (t), and ωβ(t), respectively. 15) Calculate ( )0-var,qdh j tR using (C.22). 16) Compute the electromagnetic torque Te (t), the rotor speed ωr (t), and angle θr (t) using 50 (2.33), (3.34), and (3.35), respectively. 3.5 Computational Efficiency Assessment It is a common practice to evaluate the computational efficiency of the discretized models based on the number of floating point operations (flop) and nonlinear functions (e.g., trigonometric functions (trig)) per time-step [27], [29], [30]. One flop is considered one addition, subtraction, multiplication, and division of two numbers [63]. Since the direct evaluation of the trigonometric functions is expensive, the method proposed in [29] is adopted here. The VBR formulations include other nonlinear functions such as square root (sqrt) and inverse trigonometric (inv trig) functions, as in Sections 3.4.1 and 3.4.2. The computational cost of the new CPVBR models is compared to that of the saturable PD [42] (which will be discussed in Section 4.2 of this thesis) and VBR [43] models. Herein, a typical salient-pole synchronous machine is considered with three rotor windings (i.e., one q-axis damper winding, one d-axis damper winding, and one field winding). The machine parameters are summarized in Appendix D and in [64]. The number of flops and nonlinear functions for all subject models is summarized in Table 3.1, wherein all possible savings have been exploited, including pre-computing all constant coefficients and matrices and utilization of the symmetry. The flops are roughly assigned to different terms/equations. Due to modelling the saturation characteristic as a piecewise-linear function, additional flops are required to update the saturation-segment-dependent variables whenever the saturation level changes. This extra number of flops is also provided in Table 3.1 under the name of the saturation subroutine (sat. sub.) term. The total number of flops is reported in Table 3.1 in case of saturation level changes and considering the cost of nonlinear functions. 51 Table 3.1 Flops and nonlinear functions count per time step PD model VBR model CPVBR-I model/ CPVBR-II model Term flops nl. func. Term flops nl. func. Term flops nl. func. ( )( )( ),,,s jsr jrs jtttLLL 38 2(1) ( )( )( ),1abcs jssttt−LKK 32 2(1) ( )( )1sstt−KK 14 2(1) eqpdR 9 0 eqvbrR 59 0 eqcpvbrR 0 0 pdhe 124 1(2) vbrhe 115 2(3) cpvbrhe 177 2(3) eqpdj 68 0 eqvbrj 68 0 eqcpvbrj 12 0 iabcs 18 0 iabcs 18 0 iabcs 15/ 81 0 iqdr 25 0 λqdr 26 0 λqdr 26 0 λm 11 1(2) λm 16 1(2) λm 16 1(2) Te 4 0 Te 4 0 Te 4 0 Sat. sub. 66 0 Sat. sub. 157 2(3) Sat. sub. 101 2(3) Sub-Total 363 4 Sub-Total 495 7 Sub-Total 365/ 431 7 Total 397 - Total 551 - Total 421/487 - (1) trigonometric functions (approx. 15 flops), (2) square root function (approx. 2 flops), and (3) inverse trigonometric (approx. 10 flops) and square root functions. Table 3.1 demonstrates that the VBR model requires 551 flops; the PD model requires 397 flops; the CPVBR-I model requires 421 flops; and the CPVBR-II model requires 487 flops, respectively. The improvement in the CPVBR-I model with respect to the VBR model originates from a more efficient evaluation of abcsi , eq ,cpvbrj and the saturation subroutine realized by having an actual constant-parameter interfacing conductance matrix and effective utilization of resulting symmetrical matrices. Although the augmented current correction feature in the CPVBR-II model increases the computational burden, the CPVBR-II still requires fewer flops than the VBR model. As shown in Table 3.1, the PD model requires the least number of flops per time step, mainly due to the saturation modelling and diagonal structure of its matrices. This motivates the formulation 52 Table 3.2 CPU time per step of the PD, VBR, CPVBR-I and CPVBR-II models Model PD VBR CPVBR-I CPVBR-II CPU Time (µs) 0.374 0.433 0.387 0.412 of efficient constant-parameter PD models that will be introduced in the next Chapter. With that, the CPVBR-I model and, to a lesser extent, the CPVBR-II model still demonstrate comparable computational efficiency. All models have also been implemented using standard C language in Visual Studio Code to analyze the computational efficiency in a standalone execution without interference from the simulation environment. The CPU times per step of the subject models are summarized in Table 3.2. As presented in Table 3.2, the PD model requires the smallest CPU per step, consistent with the number of flops reported in Table 3.1. The CPU time of the CPVBR-I and CPVBR-II models are approximately 3.5% and 10.2% greater than that of the PD model. Lastly, the VBR model is the most computationally expensive by being 15.8% more expensive than the PD model. 3.6 Computer Studies Two case studies are considered. In the first study, a single-machine-infinite-bus (SMIB) system is implemented. In the second study, all subject models are implemented in C and interfaced in PSCAD/EMTDC [10]. These models in PSCAD are then used to implement the IEEE 9-bus system. For consistency, all studies are executed on a PC with Intel® CoreTM i7-10750H CPU @ 2.60 GHz and Windows 10 OS. The 2-norm (cumulative) error [65] of the solution trajectory over the entire simulation interval is considered to examine the accuracy of the subject models. To evaluate the accuracy for 53 three-phase quantities, the average of single-phase errors is computed and reported [62]. 3.6.1 Single-Machine-Infinite-Bus (SMIB) System This study allows a detailed investigation of the numerical errors of each subject model without interference from the external network while choosing larger time-step sizes. The same synchronous machine from Section 3.5 is considered here with the machine parameters provided in Appendix D. Initially, the machine is assumed to operate with a fixed input mechanical torque Tm = 23.45×106 Nm and a fixed excitation vfd = 14.5 V [3]. This selection of inputs makes the machine operate in a deeply saturated region, increasing the magnitude of 0var,qdjR in the CPVBR models. At t = 0.04 s, a single-phase-to-ground fault is applied at the machine terminal and is removed after 200 ms. This study places the subject models under severe stress in terms of asymmetrical operation and changes in the saturation regions. Since the analytical solution is not available, the reference solution is obtained by a saturable state-space qd0 model [64] implemented in MATLAB/Simulink. The reference model is run with MATLAB’s solver ode4 with a very small fixed time step of 1µs to obtain numerically highly accurate results. The transient solutions are shown in Figure 3.3. All subject models are run at time-step sizes ranging from 50 µs to 1 ms. For a qualitative comparison, the magnified views of ics and ifd following the fault and Te just after removal of the fault, obtained by the subject models with the time step of 1 ms, are illustrated in Figure 3.4. As shown in Figure 3.4(a)–(c), the VBR model and, to a lesser extent, the PD model predict all variables with satisfactory accuracy. This comes at the cost of inversion of their conductance matrix at each time step. The CPVBR-I model is observed to be slightly less accurate than the VBR and PD models. With the current correction feature, the CPVBR-II model produces more accurate results than the CPVBR-I model, close to the VBR 54 Figure 3.3 Transient responses of the subject models to a single-phase-to-ground fault: (a) phase a current ias, (b) phase b current ibs, (c) phase c current ics, (d) electromagnetic torque Te, (e) field current ifd, and (f) main flux λm. model [e.g., see Figure 3.4(b)]. Also, as can be seen in Figure 3.4(c), all models obtain highly accurate results for the field current. 55 Figure 3.4 Magnified view of variables in Figure 3.3: (a) phase c stator current ics, (b) electromagnetic torque Te, and (c) field current ifd. To provide a more detailed perspective into the accuracy of the subject models, the average of single-phase stator current 2-norm errors ε (iabcs) and the 2-norm error of electromagnetic torque ε (Te) are illustrated in Figure 3.5 as a function of time-step size ∆t. As shown in Figure 3.5, since 56 Figure 3.5 2-norm errors of stator currents ε (iabcs) and electromagnetic torque ε (Te) as a function of time-step size ∆t for the single-phase-to-ground fault case study. all subject models are algebraically equivalent, and they asymptotically converge to the reference solution as the time-step size decreases. It is noted in Figure 3.5 that the CPVBR-I model predicts results with an accuracy close to the VBR and PD models up to the time step of 500 µs, illustrating that the current prediction error is less than the discretization error. At time steps above 500 µs, the prediction error starts to dominate the discretization error (i.e., the higher slope of error increase vs the common trend followed by the PD and VBR models). Also, as presented in Figure 3.5, the CPVBR-II model produces highly accurate results similar to the VBR model showing that the current correction scheme greatly minimizes the current prediction error (i.e., the consistent slope of error increase vs the common trend followed by the PD and VBR models). 57 Figure 3.6 One-line diagram of the IEEE 3-generator 9-bus system. 3.6.2 Multimachine System The standard IEEE 9-bus system shown in Figure 3.6 is considered next. This system is implemented in PSCAD/EMTDC [10]. To prevent the error that the drift of rotor angles may produce, causing the solution of electrical variables to drift at different time-step sizes, the swing generator connected to bus 1 [see Figure 3.6] is replaced by an ideal three-phase source with a fixed frequency. The parameters of generators G2 and G3 are obtained based on [66] and provided in Appendix E. All machines include main flux saturation with characteristics obtained using [67]. The rest of the parameters, including network, transformers, lines (modelled as lumped-parameter PI sections), and loads (modelled as constant impedance), can be found in [33]. The nominal voltages of the primary of transformers T1 and T3 are changed to 24 kV and 15.5 kV to match the data for the synchronous machines. With the slack generator operating at 24.96 kV, the system is initially in a steady state operating point specified in [33]. A three-phase-to-ground fault occurs at bus 1 [see Figure 3.6] at t = 0.025 s, and is cleared after 200 ms. For comparison, the reference solution is obtained using G2 G3ST2 T3T1321469875abc-giT1a58 Figure 3.7 Three-phase-to-ground fault study in IEEE 9-bus system with the time-step of 500 µs: (a) phase a current of the primary side of T1 bus 1 (b) magnified view of phase a current of the primary side of T1. the PD model with a small time step of 10 µs, similar to [30] and [62]. The simulation results computed with a time step of 500 µs for phase a of the primary side of T1, iT1a, are illustrated in Figure 3.6(a). A zoomed-in view of the current immediately after fault removal is presented in Figure 3.6(b). As demonstrated in Figure 3.6(b), all models produce similar results with a small numerical error, consistent with the results reported in Section 3.6.1. For a quantitative accuracy analysis of each model, the average 2-norm error of single-phase currents contributed by the two generators iT1 is reported in Table 3.3. As presented in Table 3.3, the proposed CPVBR models obtain comparable accuracy to their variable-parameter counter- 59 Table 3.3 Average 2-norm error of fault current contributed by the two synchronous generators for different time-step sizes Time-step size PD VBR CPVBR-I CPVBR-II 100 µs 0.9741 % 0.9747 % 0.9742 % 0.9747 % 500 µs 3.7550 % 3.7055 % 3.6941 % 3.7064 % Table 3.4 CPU per time step required for simulation of the IEEE 9-bus system with different models normalized by that of the VBR model Models qd0 PD VBR CPVBR-I CPVBR-II CPU Time (pu) 0.6552 0.9944 1 0.5803 0.5998 parts testifying to the efficacy of the method in Section 3.4 to make constant-parameter models. It can also be observed from Table 3.3 that since the numerical error associated with solving the network equations dominates the simulation results, the current correction scheme in CPVBR-II is less effective here compared to the SMIB study, where the numerical errors are only related to the machine model and not the network. Next, the CPU time per step of each model has been normalized by that of the VBR model and is presented in Table 3.4. As shown in Table 3.4, the CPVBR-I model outperforms the PD and VBR models by a factor of 1.71 and 1.72, respectively. It is noted that the CPVBR-II model also achieves faster simulation time by being only 3.4% slower than the CPVBR-I model. The excellent performance of the CPVBR models illustrates the computational efficiency gained from having an actual constant-parameter interfacing circuit. Finally, the efficiency of the proposed CPVBR models is also compared to the PSCAD’s built-in qd0 model (although it applies saturation to the d-axis only). As reported in Table 3.4, the 60 proposed CPVBR-I and CPVBR-II models are slightly more efficient (about 1.13 and 1.09 times faster, respectively) despite having a more complicated representation of saturation, but also a constant-parameter interfacing conductance matrix (similar to PSCAD’s qd0 model). 61 Phase-Domain Synchronous Machine Model with Main Flux Saturation for EMTP-Type Solution To set the stage for the derivation of constant-parameter PD models (CPPD) in Chapter 5, the saturable varying-parameter PD model considering main flux saturation is proposed in this chapter. The new EMTP PD model is achieved by applying the implicit Trapezoidal rule to the state-space PD model, which directly interfaces with the external network through a variable interfacing circuit (i.e., both rotor-position- and saturation-segment-dependent). The magnetic saturation is included as laid out in Section 3.1. Extensive computer studies in a single-machine-infinite-bus (SMIB) system are provided for demonstrating the differences among the common approaches for representing magnetic saturation in synchronous machines (i.e., neglecting saturation, considering only d-axis saturation, and considering main flux saturation). It is shown both neglecting saturation and d-axis saturation can lead to noticeable errors in the prediction of internal machine variables, including rotor angle. Furthermore, similar to Chapter 2, the impact of having indirect interfacing vs direct interfacing is also investigated but in the EMTP environment. It is demonstrated that the numerical accuracy of the PSCAD’s built-in classical qd0 model deteriorates significantly as the time-step increases due to having a one-step delay in its formulation. In contrast, the accuracy of the proposed PD model is highly preserved. 62 4.1 Model Derivation 4.1.1 State-Space Model To lay out the foundation for deriving the discrete PD model, the state-space model is briefly reviewed here. The voltage equation of the machine considering the flux linkages as state variables is formulated as [28], [44] 00abcs abcs abcssqdr qdr qdrrp      = +            v i λRv i λR, (4.1) where Rs and Rr are the stator and rotor matrices defined (2.1) and (2.2), respectively; and the flux linkages in (4.1) are calculated as ( ) ( )( ), ,, ,abcs abcss j r sr j r ressqdr qdrrs j r r j resr       = +           λ iL L λλ iL L λ, (4.2) with stator and rotor residual flux vector defined as 1ress s resqd−=λ K λ , 1 λ , λ , 0Tresqd s resq resd−  =  λ K , (4.3) λ , λ , λTresr resq resd resd =  λ , (4.4) where 1s−K is the inverse of Park’s transformation [3] (see Appendix A) and the definitions of matrices ( ),s j rL , ( ),sr j rL , ( ),rs j rL , and ,r jL can be found in Appendix F. Lastly, the qd-axes magnetizing fluxes required for determining saturation level are obtained as ( )2 , 1λ λmq D j qs kq resqSF L i i= + + , (4.5) 63 ( ),λ λmd D j ds kd fd resdL i i i= + + + . (4.6) 4.1.2 Discrete Model Applying the implicit trapezoidal integral rule with the fixed time step of t to (4.1), the discrete PD model can be obtained as ( ) ( )( ) ( ) ( ) ( ) ( ) ( ), , pdabcs s s j abcs ress sr j qdr sht k t t k t k t t t= + + + +v R L i λ L i e , (4.7) ( ) ( ) ( ) ( ) ( ) ( ) ( ), ,qdr r r j qdr resr rs j abcs rht k t k t k t t t= + + + +v R L i λ L i e , (4.8) where 2kt=, and the history voltage terms are given as ( ) ( )( ) ( ) ( ),pdsh s s j abcs resst k t t t t k t t= − − − − −e R L i λ ( ) ( ) ( ),sr j qdr abcsk t t t t t t− − − − −L i v , (4.9) ( ) ( ) ( ) ( ),rh r r j qdr resrt k t t k t t= − − − −e R L i λ ( ) ( ) ( ),rs j abcs qdrk t t t t t t− − − − −L i v . (4.10) Given known rotor voltages at the current time step t, the rotor currents of the machine are calculated from (4.8) as ( ) ( ) ( ) ( )( ( ) ( ) ( ))1, ,qdr r r j qdr resr rs j abcs rht k t k t k t t t−= + − − −i R L v λ L i e . (4.11) Furthermore, substituting (4.9) into (4.7), the compacted form of the PD model can be formulated as ( ) ( ) ( ) ( )eq,pd pdabcs j abcs ht t t t= +v R i e , (4.12) where the equivalent resistance matrix ( )eq,pdj tR is given as 64 ( ) ( ) ( )( ) ( )12eq, , , , ,pdj s s j sr j r r j rs jt k t k t k t−= + − +R R L L R L L , (4.13) and the full history term ( )h te as ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( )1, ,pd pdh sr j r r j qdr resr rh ress sht k t k t k t t k t t−= + − − + +e L R L v λ e λ e . (4.14) Since in the nodal-analysis-based EMTP programs, the machine is interfaced with the external network as a current source with a parallel admittance matrix, the Norton equivalent of (4.12) is computed as ( ) ( ) ( ) ( )eq, eqpd pdabcs j abcst t t t= −i G v j , (4.15) ( ) ( )1eq, eq,pd pdj jt t− =  G R , (4.16) ( ) ( ) ( )eq eq,pd pdj ht t t=j G e , (4.17) where ( )eq,pdj tG the is the equivalent Norton admittance/conductance matrix and ( )eqpd tj is the injection current source. 4.1.3 Simulation Steps 1) Predict the rotor speed ωr (t) and the magnetizing flux linkages λmq (t) and λmd (t) using linear two-point extrapolation, i.e., using (3.31) [4]. 2) Calculate the rotor angle θr (t) based on the trapezoidal integration of the predicted rotor speed ( )r t using (3.32). 3) Compute the main flux λm (t) using (3.1) based on the predicted values of ( )λmq t and ( )λmd t . 65 4) Calculate the residual flux vector ( )resqd tλ using (3.7), (3.8), ( )λm t , ( )λmq t and ( )λ .md tThen, combine all and evaluate the stator and rotor residual flux vectors ( )ress tλ and ( )resr tλ using (4.3) and (4.4), respectively. 5) Compute stator history term ( )pdsh te using (4.9), calculate rotor history term ( )pdrh te using (4.10), and combine all history terms to form ( )pdh te using (4.14). 6) Compute ( )eq,pdj tR using (4.13). 7) Compute ( )eq,pdj tG using (4.16) and calculate the equivalent Norton current source ( )eqpd tj using (4.17). Then, the calculated injection current source is inserted into the system solver. 8) Obtain the solution of the network (i.e., solving the network for unknown voltages). 9) Receiving the stator voltages vabcs (t) from the system solver, compute stator currents iabcs(t) using (4.15). 10) Compute the rotor currents ( )qdr ti using (4.11). Next, evaluate the magnetizing flux linkages ( )λmq t and ( )λmd t , and the main flux ( )λm t using (4.5), (4.6), and (3.1), respectively. 11) Compare ( )λm t with the limits of the current saturation segment. If required, update the saturation segment/level and the relevant variables as specified in step 12). Otherwise, skip to step 13). 12) Update j and re-calculate ( )resqd tλ , ( )ress tλ , ( )resr tλ , ( )qdr ti , ( )λmq t , ( )λmd t , respectively. 13) Compute the electromagnetic torque Te (t), the rotor speed ωr (t), and angle θr (t) using (2.33), (3.34), and (3.35), respectively. 66 4.2 Computer Studies For the purpose of validating and benchmarking the new saturable PD model, the proposed model from Section 4.1 has been implemented in C and interfaced inside PSCAD with its computational engine. The linear PD model [27] and the saturable VBR model with main flux saturation [43] have also been implemented inside PSCAD in the same way. For comparison purposes, the built-in PSCAD’s qd0 model with only d-axis saturation has also been considered, and the saturable VBR model in [43] is chosen as the reference model. To highlight the differences among the considered methods for representing saturation (i.e., neglecting saturation, d-axis saturation only, and main flux saturation) and the accuracy of the subject models, a simple single-machine infinite-bus (SMIB) system is assumed in this paper. Without loss of generality, a salient-pole synchronous machine with three rotor windings (i.e., one q-axis damper winding, one d-axis damper winding, and one field winding) is chosen where the machine parameters and the saturation characteristic are summarized in Appendix D. Lastly, for consistency, the reference model is run with 10 µs ensuring highly accurate solution. 4.2.1 Change of Stator Voltage In this study, the SMIB system is operating at no-load under nominal excitation voltage (i.e., 1 1xfd fd md fdE r X v−= = pu [3]). To have all models operate from the same steady-state operating point, a new magnetizing inductance obtained from the modified airgap line with 2.3318mdL = mH is used in the magnetically-linear PD model [31], [43]. To mitigate the impact of discretization error on the accuracy of simulation results, all subject models are run with a small time-step size of 50 µs guaranteeing a negligible numerical error. For the subject transient study, at t = 0.04 s, 67 Figure 4.1 Simulation results obtained by the subject models for the stator voltage step-up transient for the selected variables: (a) q-axis stator current qsi , (b) d-axis stator current dsi , and (c) electromagnetic torque eT . 68 the stator voltage steps up from 0.7 pu to 1.05 pu pushing the machine deep into the saturation region. The simulation results obtained by the subject models for the q-axis stator current iqs, the d-axis stator current ids, and electromagnetic torque Te are provided in Figure 4.1. As can be seen from Figure 4.1, the proposed PD model (PD w/ Main Flux Sat.) matches the reference model in all selected variables. This is expected since both the PD and VBR models are algebraically equivalent, and they converge to the same solution at small time-step sizes. Furthermore, since the machine is operating in idle mode, the deviations of the rotor angle from zero is negligible during the transient, which causes the main flux angle β and, consequently, the projection of the main flux vector on the q-axis to be close to zero [see Figure 3.2]. Hence, only the d-axis path of the flux in the machine is subjected to saturation under the no-load transients. Accordingly, as also demonstrated in Figure 4.1, both the saturable PD model (PD w/ Main Flux Sat.) and the qd0 model with the d-axis saturation (qd0 w/ d-axis Sat.) predict identical results. However, the linear model predicts the d-axis current, which is subjected to saturation, with noticeable error, as seen in Figure 4.1(b). Lastly, the transient responses of the electromagnetic torque predicted by all models are identical as expected and shown in Figure 4.1(c). 4.2.2 Change of Mechanical Torque The following step-up torque-change transient study evaluates the efficacy of the proposed PD model for an accurate representation of the q-axis saturation and cross-saturation. In this regard, the machine starts from an idle steady state under nominal excitation and terminal voltages. At t = 0.04 s, the input mechanical torque of the machine is stepped up to 1.2 pu, making the machine to operate in the generating mode. Similar to the previous study, to draw a fair comparison between the linear and the saturable models, the modified airgap line with 2.2014mdL = mH has 69 Figure 4.2 Reference transient trajectories of the mechanical torque step-up case study obtained for the selected variables: (a) q-axis magnetizing flux λmq , (b) d-axis magnetizing flux λmd , and (c) main flux angle  . been considered for the linear PD model, and all the models are run with the time step of 50 µs. The predicted transient trajectories of the subject case study are illustrated in Figures 4.2 and 4.3. As shown in Figure 4.2(a)‒(b), following the change in the mechanical torque, the q-axis magnetizing flux λmq undergoes a drastic decrease. In contrast, the d-axis magnetizing flux λmd drops less significantly. This observation suggests that the main flux angle β for the loaded machine is no longer 90 degrees [see Figure 3.2], which is in accordance with Figure 4.2(c). Similar to the previous study, the proposed PD model (PD w/ Main Flux Sat.) obtains highly accurate results identical to the reference model, as depicted in Figure 4.3. However, having a nonzero projection of the main flux vector onto the q-axis, the d-axis saturable and linear models predict noticeably different trajectories, as shown in Figure 4.3. Specifically, these two models obtain very similar yet inaccurate results for the qd-axes stator currents, as shown in Figure 4.3(a)‒ 70 Figure 4.3 Simulation results obtained by the subject models for the mechanical torque step-up transient for the selected variables: (a) q-axis stator current qsi , (b) d-axis stator current dsi , (c) electromagnetic torque eT , and (d) rotor angle δ. 71 (b). Although the predicted trajectory of the electromagnetic torque is approximately similar for all models (except for a part of the transient), the models without the q-axis saturation (i.e., qd0 w/ d-axis Sat. and Linear PD) predicted rotor angle response with considerable error both in transients and steady state. This observation confirms the necessity of including the main flux saturation into the machine models for accurate and reliable computer studies. 4.2.3 Change of Stator Voltage Under Full Load A step change in the stator voltage has been considered here to further investigate the impact of direct/indirect interfacing on the accuracy of machine models, especially at large time-step sizes. The proposed fully saturable PD model is benchmarked against the PSCAD’s built-in model, running at two different time steps of 50 µs and 1000 µs. In the subject study, the machine starts under nominal conditions with Te =1 pu and Exfd = 1 pu. At t = 0.04 s, the terminal voltage increases from from 1 pu to 1.05 pu. The simulated trajectory of the phase a current ias is shown in Figure 4.4. As illustrated in Figure 4.4(a)‒(b), the results obtained by the PD model are accurate and close to the reference model even at the very large time step of 1000 µs. This performance is attributed to having a direct interfacing circuit enabling a simultaneous solution of the entire network with the machine’s equations. However, the performance of the qd0 model drastically deteriorates as the time-step size increases, i.e., see Figure 4.4. This is also expected due to the existing intrinsic one-step delay in the interfacing of qd0 model in PSCAD [10]. Observations in Figure 4.4 verify the capability of the proposed PD model to have robust numerical performance even at large step sizes. 72 Figure 4.4 Simulation results obtained by the subject qd0 and PD models for the step-up stator voltage transient under full load for the stator phase a current asi : (a) overall transient, and (b) magnified view. 73 Constant-Parameter PD Synchronous Machine Models Considering Main Flux Saturation for EMTP-Type Solution Based on the proposed PD model from previous Chapter, a group of saturable constant-parameter PD models (CPPD-q-I/II, CPPD-d-I/II, and CPPD-qd-I/II) with main flux saturation are proposed in this Chapter. The new CPPD models are considered as an extension of the CPVBR models that are introduced in Chapter 3. In this regard, the genesis of the interfacing resistance matrix in the variable-parameter PD model (introduced in Chapter 4) is analyzed in qd0 coordinates, based on which approaches for constructing a constant-parameter interfacing circuit including appropriate reformulation of history terms are proposed. The introduced systematic procedure can be straightforwardly extended for formulation of synchronous machines with any number of phases. Furthermore, similar to the CPVBR-II model in Chapter 3, an effective current-correction scheme is also augmented in the CPPD-q model making this model the enhanced extension of the linear CPPD model proposed in [30]. Considerable effort has also been invested to boost the computational efficiency of the CPPD models by reducing the number of flops required by these models at each time step. Lastly, extensive computer studies are carried out in MATLAB and PSCAD/EMTDC to validate the numerical performance of the proposed models in single-machine and multimachine systems. It has been shown that the new CPPD models can achieve significant computational efficiency compared to the original VBR and PD models and, more importantly, compared to the new CPVBR models (both of which are interfaced with external network with constant-parameter 74 interfacing circuit). 5.1 Constant-Parameter Formulation Recalling the original PD formulation from Chapter 4, to avoid refactorization of the entire network conductance matrix at each time step, the interfacing resistance matrix ( )eq,pdj tR in the PD model should be constant. However, as shown in (4.13), ( )eq,pdj tR is both a function of rotor position θr and saturation segment j. In this regard, the genesis of ( )eq,pdj tR in qd0 coordinates is first analyzed, and based on that, different approaches for separating ( )eq,pdj tR into constant and variable terms are introduced in Subsection 5.1.1. Then, adopting those approaches, different constant-parameter models may be obtained based on the original PD model, as presented in Subsection 5.1.2. 5.1.1 Constant-Parameter Resistance Matrix Formulation Step I) According to [3], the self and mutual inductance matrices constituting (4.13) have the following relations in the qd0 coordinates as ( ) ( ) ( )1, ,rs j s s j st t t−=L K L K , (5.1) ( ) ( )1, ,rsr j s sr jt t−=L K L , (5.2) ( ) ( ), ,32rrs j rs j st t=L L K , (5.3) where the saturation-segment-dependent sparse matrices ,rs jL , ,rsr jL , and ,rrs jL are defined in (G.1)‒(G.2) in Appendix G. Employing (5.1)‒(5.3) in (4.13), and factoring ( )1s t−K and ( ) ,s tK75 (4.13) can be reformulated as a product of three matrices [30], ( ) ( ) ( )1 0eq,pd qdj s j st t t−=R K R K , (5.4) ( )0 , , 0diag , ,qdj q j d jR R R =  R , (5.5) where , ,q j s ls q jR r kL Z= + + , , ,d j s ls d jR r kL Z = + + , 0 s lsR r kL= + . (5.6) In (5.6), the so-called equivalent q- and d-axis subtransient impedances ,q jZ  and ,d jZ  are defined as [29] ( )11 1 1, , 1 2q j mq j lkq lkqZ Z Z Z−− − − = + + , (5.7) ( )11 1 1, ,d j md j lfd lkdZ Z Z Z−− − − = + + , (5.8) with 2, ,mq j D jZ kSF L= , , ,md j D jZ kL= , l lZ r kL  = + , (5.9) where = 1, 2, , kq kq fd kd . It is noted that the rotor-position-dependency in ( )eq,pdj tR only comes from the transformation matrices ( ) ,s tK and ( )1s t−K , whereas the saturation-dependency only comes from the middle matrix 0qdjR in (5.4), which happens to have a diagonal structure. As demonstrated in [4], [22], the resulting matrix from the product in (5.4) can become rotor-position-independent if the middle matrix has a unique diagonal structure as ( ) ( )1:0 020 030 0abcs st t v        −=      = =         MM K K , (5.10) 76 with 2  = + , and 2 2  = − + . (5.11) Thereafter, the middle matrix in the product of (5.10) is named Mα. As observed in (5.10), the first and the second diagonal entries in Mα are equal by construction, resulting in Mabc becoming time-invariant. Therefore, 0qdjR in (5.4) has to be reformulated such that it has two equal diagonal entries similar to Mα, which will be realized in the following Steps II to IV as demonstrated in Table 5.1. Step II) After obtaining 0qdjR in Step I [see (5.5) and (T1) in Table 5.1], three approaches are possible for forming the special matrix Mα out of 0qdjR by setting α in (5.10) equal to: 1) the equivalent q-axis resistance Rq,j; 2) the equivalent d-axis resistance Rd,j; and 3) the average of the equivalent qd-axes resistances ( ), ,0.5 q j d jR R+ . Adopting these approaches, 0qdjR can be rewritten in each case as a summation of two diagonal matrices, with the first term possessing the unique structure of Mα in (5.10). These first terms in Table 5.1 [see (T2)‒(T4)] corresponding to the three approaches are referred to as qjR , djR , and qdjR , respectively. Accordingly, in each approach, the other diagonal matrix (i.e., var,qjR , var,djR , and var,qdjR , respectively) is obtained as the difference between 0qdjR and the specially-arranged matrices qjR , djR , and qdjR [see (T2)‒(T4) in Table 5.1]. Although substituting the new matrices qjR , djR , and qdjR into (5.4) yields rotor-position independency, the resulting matrices in abc coordinates are still saturation-segment-dependent. Step III) Considering (5.5)‒(5.9), the only saturation-segment-dependent terms in 0qdjR are ,q jZ  and ,d jZ  , which are defined in (5.7) and (5.8), respectively, as the parallel combination of the so-called equivalent magnetizing impedances, Zmq,j and Zmd,j, and the rotor windings impedances 77 Table 5.1 Genesis of the Interfacing Resistance Matrix of the PD Model in qd0 Coordinates and the Sequential Steps for Constructing Constant-Parameter Resistance Matrix Step I Derivation of the initial diagonal matrix with unequal entries Step II, Separating diagonal part with equal entries and the remaining part Step III, Computing part for first sat. segment j = 1 Step IV, Combing the saturation-dependent and variable terms Approach 1: Considering the first diagonal entry, q-axis 0qdj =R ( ), , 0diag , ,q j d jR R R   (T1) ( )0 , , 0diagqdj q j q jR R R =  R ( ), , var,diag 0 0 q qd j q j j jR R + − = +  R R (T2) 1q q qj j= +R R R (T5) var, var,q q qj j j =  +R R R (T8) 01 var,qd q qj j= + R R R (T9) Step II, Separating diagonal part with equal entries and the remaining part Step III, Computing part for first sat. segment j = 1 Step IV, Combing the saturation-dependent and variable terms Approach 2: Considering the second diagonal entry, d-axis ( )0 , , 0diagqdj d j d jR R R =  R ( ), , var,diag 0 0 d dq j d j j jR R + − = +  R R (T3) 1d d dj j= + R R R (T6) var, var,d d dj j j =  +R R R (T10) 01 var,qd d dj j= + R R R (T11) Step II, Separating diagonal part with equal entries and the remaining part Step III, Computing part for first sat. segment j = 1 Step IV, Combing the saturation-dependent and variable terms Approach 3: Considering the average of the first and second diagonal entries, averaging q- and d-axes ( ) ( )( )0 , , , , 0diag 0.5 0.5qdj q j d j q j d jR R R R R = + + R( ) ( )( ), , , , var,diag 0.5 0.5 0 qd qdq j d j d j q j j jR R R R + − − = +  R R (T4) 1qd qd qdj j= + R R R (T7) var, var,qd qd qdj j j =  +R R R (T12) 01 var,qd qd qdj j= + R R R (T13) 78 Zlγ. Since the incremental inductance LD,j is much larger than the leakage inductances in (5.9) [68], even when the machine saturates, the values of ,q jZ  and ,d jZ  , and hence the magnitude of 0qdjR are insignificantly affected by saturation. As a result, for each approach in Table 5.1, the matrices qjR , djR , and qdjR obtained in Step II, can be further reformulated as a summation of constant and considerably smaller saturation-segment-dependent matrices, as presented in (T5)‒(T7) in Table 5.1 [see Step III]. Therein, the constant matrices denoted by 1qR , 1dR , and 1qdR correspond to the evaluation of qjR, djR , and qdjR using the 1st unsaturated segment at j = 1 in Figure 3.1, respectively. Step IV) Combining the equations in Steps II and III for each approach and grouping the saturation-segment-dependent terms as a separate matrix [see qjR in (T8), djR in (T10), and qdjR in (T12) in Table 5.1], the original resistance matrix 0qdjR can be ultimately decomposed into constant and variable matrices as presented in (T9), (T11), and (T13), respectively. Finally, after undergoing the four steps of Table 5.1 and substituting (T9), (T11), and (T13) from Table I into (5.4), the variable resistance matrix ( )eq,pdj tR in (5.4) is rewritten as ( ) ( ) ( )1eq, eq-const var,pdj s j st t t −= + R R K R K , (5.12) where the superscript γ = {q, d, qd} denotes the approach/choice in Table 5.1. In (5.12), the true constant (i.e., both rotor-position- and saturation-segment-independent) matrix eq-constR is the conversion of 1R into abc coordinates and has the same structure as Mabc in (5.10), where β equals R0 and α takes different values for different approaches in Table 5.1. Specifically, for Approaches 1, 2, and 3, α equals 1qR , 1dR , and ( )1 10.5 q dR R+ , respectively. 79 5.1.2 Formulation of Interfacing Circuit with Constant Conductance Matrix Substituting (5.12) in (4.12), a constant-parameter interfacing circuit for each approach is derived by moving the resulting term ( ) ( ) ( )1 var,s j s abcst t t− K R K i into the history term as ( ) ( ) ( )eq-constabcs abcs ht t t = +v R i e , (5.13) ( ) ( ) ( ) ( ) ( )1 var,pdh s j s abcs ht t t t t −=  +e K R K i e . (5.14) Hereafter, the superscripts γ = {q, d, qd} is used to denote each Approach in Table 5.1, i.e. referring to CPPD-q, CPPD-d, and CPPD-qd models, for Approach 1, 2, and 3, respectively. In the evaluation of the history term in (5.14), the stator current at the current time step t is required, but it is yet to be computed by the network solver. Therefore, unless an iterative network solver is used, the stator currents should be approximated at the beginning of each time step to evaluate (5.14). Herein, similar to the formulation of the CPVBR models in Chapter 3, the two-point linear extrapolation in qd coordinates in (3.51) is adopted here [68]. To avoid converting (3.51) into abc coordinates, the history term (5.14) is simplified as ( ) ( ) ( ) ( )1 var, 0pdh s j qd s ht t t t −=  +e K R i e . (5.15) In applying the Trapezoidal integral rule to the state-space PD model in (4.1), the discretization has been achieved only along a given jth saturation segment/level. This requires that the coefficients in (5.13) and (5.14) be from the same jth saturation segment. Therefore, in (5.13)‒(5.14), the coefficients of ( )abcs ti and ( )abcs t t−i have to be from the first saturation segment (since eq,constR has been evaluated for the unsaturated segment j = 1), and the coefficients of ( )abcs ti and ( )abcs t t−i have to be saturation-dependent. In this regard, substituting (4.11) in (4.9) from 80 Chapter 4, and factoring out the terms containing ( )abcs t t−i , we have ( ) ( ) ( ) ( ) ( ), ,pd pdsh h j abcs ress sr jt t t t t k t t k t t= − − − − − − e R i λ L ( ) ( )( ( ) ( )) ( )1,r r j qdr resr rh abcsk t t k t t t t t t−+ − − − − − − −R L v λ e v , (5.16) with ( ) ( ) ( )( ) ( )12, , , , ,pdh j s s j sr j r r j rs jt t k t t k t t k t t−− = − − + − + −R R L L R L L . (5.17) Since (5.17) and (4.13) share a similar structure, ( ),pdh j tR can also be separated into a constant and a variable term based on the procedure in Steps I to IV in Subsection 5.1.1. Consequently, the history term (5.16) can be reformulated for each Approach in Table 5.1 as ( ) ( ) ( ) ( ) ( ) ( )1const var, 0 ,pdsh h abcs s h j qd s ress sr jt t t t t t t k t t k t t −− −= − + −  − − − − − e R i K R i λ L ( ) ( )( ( ) ( )) ( )1,r r j qdr resr rh abcsk t t k t t t t t t−+ − − − − − − −R L v λ e v , (5.18) where the definitions of consth−R and var,h j−R are provided in (G.3)‒(G.9) in Appendix G. This way, the coefficients in front of ( )abcs ti in (5.13) and ( )abcs t t−i in (5.18), and the coefficients in front of ( )abcs ti in (5.14) and ( )abcs t t−i in (5.18) become consistent in terms of saturation level. Consequently, in addition to a unique interfacing resistance/conductance matrix, each CPPD model also has a unique stator history term defined in (5.18). Furthermore, for implementation in EMTP environment, the Norton equivalent of (5.13) is calculated as ( ) ( ) ( )eq-constabcs abcs ht t t = −i G v j , (5.19) 81 Figure 5.1 Norton equivalent representation of the interfacing circuit of synchronous machine model for EMTP-type solution where G is the equivalent conductance matrix, and j is the equivalent injection current vector. where 1eq-const eq-const  − =  G R and the history injection current source is obtained as ( ) ( )eq-consth ht t  =j G e , (5.20) with γ = {q, d, qd}. The equivalent circuit of (5.19) is shown in Figure 5.1, wherein the generic matrix G and vector j are equal to eq-constG and ( )h tj for each model, respectively. It is noted that the original PD and CPPD models share the same equivalent circuit as in Figure 5.1; however, the former has a variable, and the latter has a constant interfacing conductance matrix. In the formulation of CPPD-q/d/qd models, no terms have been neglected, and the only approximation is the prediction of stator currents. Hence, if the prediction of stator current in (3.51) is exact (e.g. in steady-state), the CPPD and PD models will yield the same solutions. 5.2 Augmented Current Correction Scheme A current correction feature has been recently introduced in [30] to improve the accuracy of the linear CPPD model. Later, it was adopted for deriving the state-of-the-art saturable CPVBR-vasvbsvcsG j82 II model [68]. The accuracy of CPPD-q/d/qd models can also be enhanced here by considering (4.15), which uses the variable conductance matrix ( )eq,pdj tG to calculate the machine’s stator currents. Combing (4.15)‒(4.17) and (5.15), the stator currents for each model are calculated as ( ) ( ) ( ) ( )( ( ) ( ))1eq, var, 0pdabcs j abcs h s j qd st t t t t t −= − + i G v e K R i . (5.21) Hereafter, the models that do not adopt the correction scheme in (5.21) are referred to as CPPD-I, while the models that adopt the correction scheme (5.21) are referred to as CPPD-II. This naming convention is used in this paper in addition to letters that signify the Approach for making the interfacing matrix in Table 5.1 (i.e., q, d, qd for Approaches 1, 2, and 3 in Table I, respectively). The direct calculation of ( )eq,pdj tG in (5.21) (i.e., forming ( )eq,pdj tR using (4.13) and then calculating its inverse) is costly. Therefore, the calculation of ( )eq,pdj tG should be optimized, wherein all constants and saturation-dependent terms should be calculated in advance. Specifically, ( )eq,pdj tG can be analytically calculated considering (5.5) and the property of inverse of the product of matrices as ( ) ( ) ( ) ( )11 0eq, eq,pd pd qdj j s j st t t t−− = = G R K G K , (5.22) ( )0 , , 0diagqdj q j d jG G G =  G , (5.23) where 1, ,q j q jG R− =   , 1, ,d j d jG R− =   , and 10 0G R−= . Furthermore, since ( )eq,pdj tG has a similar form as ( )eq,pdj tR [compare (5.22) and (5.5)], it can be computed more efficiently as ( ) ( )eq, var,pdj j jt t= +G G G , (5.24) with 83 1, 2, 2,2, 1, 2,2, 2, 1,23j j jj j j jj j jG G GG G GG G G  =    G , (5.25) and ( ) ( )( )( )( )var, , ,2 2 23 31 2cos 2 2 23 3 322 2 23 3r r rj q j d j r r rr r rt G G              − +                = −  − −               + +           G , (5.26) where 1, jG and 2, jG are defined in (G.10) in Appendix G. All saturation-segment-dependent terms in (5.25) and (5.26) can be computed during the initialization step reducing the computational burden of calculating ( )eq,pdj tG at each time step. Hence, the costly inversion of ( )eq,pdj tR is avoided, contrary to the original formulation [30]. Finally, based on Subsection 5.1.2, for having a consistent saturation level in (5.21), the stator history term ( )pdsh te in CPPD-II models should be calculated using (4.9), instead of (5.18), which entirely contains saturation-dependent coefficients similar to ( )eq,pdj tG used in (5.21). 5.3 Efficient Computation of History Terms To further enhance the computational efficiency of the proposed models, the calculations of history terms should also be optimized. In this respect, substituting (4.10) in (4.14) and applying (5.2)‒(5.3), and inserting the resulting equation into (5.15), the full history term he in (5.15) can be reformulated by factoring out ( )1s t−K as [30] 84 ( ) ( ) ( ) ( ) ( )( )( )1 a, rh s j rh resr resrt t t k t t t −= − − −e K M e λ λ ( ) ( ) ( ) ( )f, 0 var, 0j qd s j qd s resqd sht t t k t t + − + + +R i R i λ e , (5.27) with a, ,rj sr j jk=M L H , ( )1,j r r jk−= +H R L , (5.28) ( )12f, , , ,r rj sr j r r j rs jk k−= +R L R L L , (5.29) ( ) ( ) ( ) ( ) ( ),rrh qdr qdr r r j qdrt t t t k t t= + − − − −e v v R L i . (5.30) Inserting (4.10) in (4.11) and using (5.30), the rotor currents in (4.11) in Chapter 4 can be more efficiently calculated as [30] ( ) ( )( ( ) ( )( ) ( ) ( )( )), 0 0r rqdr j rh rs j qd s qd s resr resrt t k t t t k t t t= − − − − − −i H e L i i λ λ . (5.31) 5.3.1 Models without Correction Scheme (CPPD-I Models) Furthermore, in the CPPD-I models, by adding ( ), 0rsr j qd sL i t t − in (5.18), applying (5.2)‒(5.3) and (5.29) to the resulting equation, and grouping the terms containing ( )1s t t− −K , the stator history term ( )pdsh te in (5.18) can be effectively rearranged as ( ) ( ) ( ) ( ) ( )1const var, 0 f, 0sh h abcs s h j qd s j qd st t t t t t t t t  −− −= − + −  − − −e R i K R i R i ( ) ( ) ( ),rsr j qdr resqd abcsk t t k t t t t− − − − − −L i λ v . (5.32) 85 5.3.2 Models with Correction Scheme (CPPD-II Models) Additionally, for the CPPD-II models, the stator history voltage term can also be efficiently reformulated by combing (4.2) and (4. ) and substituting from (5.1)‒(5.2) as [30] ( ) ( ) ( ) ( )sh s abcs abcs abcst t t k t t t t = − − − − −e R i λ v , (5.33) where the stator flux vector λabcs (t) is computed as ( ) ( ) ( )1 0abcs s qd st t t−=λ K λ , (5.34) with ( ) ( ) ( ) ( )0 , 0 ,r rqd s s j qd s sr j qdr resqdt i t i t t= + +λ L L λ . (5.35) 5.4 Execution Steps Assuming the solution from t ‒ ∆t is available, the required steps to calculate the solution at the current time step t are as follows: 1) Predict rotor speed ωr (t) and the magnetizing flux linkages λmq (t) and λmd (t) using the linear two-point extrapolation in (3.31). 2) Compute rotor angle θr (t) based on the trapezoidal integration of the predicted rotor speed ( )r t using (3.32). 3) Predict stator currents 0 ( )qd s ti using (3.51). 4) Compute the main flux ( )λm t using (3.1) based on the predicted values of ( )λmq t and ( )λmd t . Then, calculate the residual flux vector ( )resqd tλ using (3.7), (3.8), ( )λm t , ( )λ ,mq tand ( )λmd t . 86 5) Calculate ( )rrh te using (5.30). 6) Calculate the stator history voltage term ( )sh te a. using (5.32) for the CPPD-q/d/qd-I models, b. using (5.33) for the CPPD-q/d/qd-II models. 7) Calculate the full history term ( )h te using (5.27). 8) Compute the injection currents ( )h tj to be inserted to the network solver using (5.20). 9) Solve the network solution for unknown voltages. 10) Using the solved voltages, obtain the stator currents ( )abcs ti a. using (46) for the CPPD-q/d/qd-I models, b. using (48) for the CPPD-q/d/qd-II models, and then convert them to qd0 coordinates. 11) Calculate the rotor currents ( )qdr ti using (5.31). 12) Compute the magnetizing flux linkages ( )λmq t and ( )λmd t using (4.5)‒(4.6) and then calculate ( )λm t using (3.1). 13) Evaluate if changing the operating point on the saturation characteristic is required. If the saturation level changes, go to step 14). Otherwise, jump to step 15). 14) Update j and recalculate ( )λresq t , ( )λresd t , and ( )qdr ti as ( ) ( ) ( )(( , 0r rqdr j rh rs j qd st t k t= −i H e L i ( )) ( ) ( )( )),( 1) 0rrs j qd s resr resrt t k t t t−− − − − −L i λ λ , (5.36) 87 Then, update the values of ( )λmq t and ( )λmd t , and var, .jR Only for the CPPD-q/d/qd-I models, update the value of var, .h j−R Also, only for the CPPD-q/d/qd-II models, re-evaluate jG , ,q jG , and ,d jG . 15) Only for the CPPD-q/d/qd-II models, calculate the stator flux linkages ( )0qd s tλ using (3.35) and convert them to ( )abcs tλ using (3.34). 16) Compute the electromagnetic torque Te (t), the rotor speed ωr (t), and angle θr (t) using (2.33), (3.34), and (3.35), respectively. 5.5 Analysis of Computational Efficiency The number of floating operating operations (flop) and nonlinear functions [i.e., trigonometric (trig) and square root (sqrt) functions] required by each model per time step is considered here to assess the computational efficiency of the proposed CPPD models. Similar to [27], [29], [30], [68], one addition, subtraction, multiplication, and division is considered herein as one flop. The computational cost of nonlinear functions is a complier-specific quantity and can differ significantly in different compilers and operating systems. Hence, similar to Chapter 3, it is assumed that, on average, evaluating one trigonometric and one square root function requires 15 and 2 flops, respectively [68]. This section considers a salient-pole synchronous machine with three rotor windings (i.e., one q-axis damper winding, one field winding, and one d-axis damper winding). The machine parameters and saturation characteristic are summarized in Appendix D. The number of flops and nonlinear functions for executing one step of CPPD models is summarized in Table 5.2, where the PD model [42] and the CPVBR-I model [68] are also compared. Table 5.2 shows that the CPPD-88 I model (without) and CPPD-II model (with the current correction scheme) are presented in general without emphasizing the Approach in Table 5.1 taken for making their constant-parameter interfacing matrix. This is because all Approaches in Table 5.1 have in overall the same (similar) computational cost. As stated in [68], since the saturation characteristic in the EMTP models is approximated with a piecewise-linear function, the changes in saturation level necessitate updating the saturation-dependent variables. Hence, an additional number of flops are required, as reported in Table 5.2 under the name of saturation subroutine (sat. sub.) term. It is also assumed that all saturation-dependent terms and constants are pre-computed, and the matrices’ symmetry has been exploited to boost the run-time computational efficiency. Table 5.2 Flops and Nonlinear Functions Count Per Time Step PD model CPVBR-I model CPPD-I models/CPPD-II models Term flops nl. func. Term flops nl. func. Term flops nl. func. Ls,j (t) Lsr,j (t) Lrs,j (t) 38 2(1) ( )( )1sstt−KK 14 2(1) ( )( )1sstt−KK 14 2(1) eqpdR 9 0 eqcpvbrR 0 0 eq,cosntR 0 0 pdhe 124 1(2) vbrhe 177 2(3) he 115/112 1(2) eqpdj 68 0 eqvbrj 12 0 eqj 12 0 iabcs 18 0 iabcs 15 0 iabcs 15/39 0 iqdr 25 0 λqdr 26 0 iqdr 20 0 λm 11 1(2) λm 16 1(2) λm 11 1(2) Te 4 0 Te 4 0 Te 4 0 Sat. sub. 66 0 Sat. sub. 101 2(3) Sat. sub. 33 0 Sub Total 363 4 Sub Total 365 7 Sub Total 224/245 4 Total 397 - Total 421 - Total 258/279 - (1) trigonometric functions, (2) square root function, and (3) inverse trigonometric (approx. 10 flops from Chapter 3 [68]) and square root functions. 89 Table 5.2 shows that, including the cost of nonlinear functions, the PD model requires 397 flops; the CPVBR-I model requires 421 flops; the CPPD-I models require 258 flops; and the CPPD-II models require 279 flops. The computational gain of the CPPD models with respect to the CPVBR-I model comes from the diagonal structure of matrices in the PD-based models. This property significantly decreases the number of flops inside the history terms he and allows a more efficient saturation subroutine. The CPPD models also outperform the PD model in terms of computational cost. The CPPD models do not require the inversion of the interfacing resistance matrix at each time step and costly evaluation of time-variable matrices (e.g., Ls,j(t), Lsr,j(t), Lrs,j(t)). The CPPP-I models require fewer flops compared to the CPPD-II models, which have a variable conductance matrix for calculating the stator currents. It is also noted that although in the formulation of the CPPD-I models, a time-variable term is added to/subtracted from the stator history terms (5.32), it can be considerably simplified using (4.11) and (5.29) such that it adds only a few flops to the overall cost of the history terms. All these computational savings combined make the CPPD-I models more efficient among all subject models. The computational cost of the subject models is also evaluated when these models are implemented in standard C. This standalone implementation objectively assesses each model’s cost without interference from the simulation environment (e.g. PSCAD and MATLAB). The computational cost per time step is summarized in Table 5.3, where the reported CPU times are obtained on a PC with Intel® CoreTM i7-10750H CPU @ 2.60 GHz and Windows 10 OS. Consistent with Table 5.2, only the CPPD-I and CPPD-II models are compared in general in Table 5.3 since the current correction scheme is the only differentiating factor impacting the overall cost. As shown in Table 5.3, and in accordance with Table 5.2, the CPPD-I models are the cheapest, whereas the CPVBR-I is the most expensive model being 16.2% slower than the former. 90 Table 5.3 CPU Time Per Step of the Models in Table 5.3 Models PD CPVBR-I CPPD-I CPPD-II CPU Time (µs) 0.374 0.387 0.333 0.349 The CPPD-II models are only 4.8% more expensive than the CPPD-I models, where a variable conductance matrix is used for the internal calculation of stator currents. 5.6 Computer Studies To investigate the numerical accuracy and computational efficiency of CPPD models, out of the possible six models introduced in Sections 5.1 and 5.2 (i.e., CPPD-q-I/II, CPPD-d-I/II and CPPD-qd-I/II), only the CPPD-qd-I and CPPD-q-II are considered here because of two reasons. First, as shown in the derivation of classical qd0 models in [4], [22], constructing a constant-parameter interfacing resistance matrix by averaging the diagonal entries (i.e., similar to Approach 3 in Table I) leads to better numerical accuracy compared to other alternatives (i.e., Approach 1 and 2 in Table I). Hence, only the CPPD-qd-I model (which does not include the current correction scheme) is considered here. Second, the CPPD-q-II and CPPD-d-II models with the current correction scheme share similar numerical performance. Therefore, only the CPPD-q-II is considered herein, which can also be regarded as the enhanced, saturable version of the linear CPPD model [30]. Two case studies are presented to investigate the accuracy and computational efficiency of the subject models. In the first study, a simple single-machine-infinite-bus (SMIB) system is considered and implemented in MATLAB to study a short-circuit transient of a salient-pole synchronous machine with main flux saturation. For comparison purposes, saturable VBR [43], 91 PD [42] (introduced in Chapter 4), and CPVBR-I (introduced in Chapter 3) [68] models are also implemented in the same environment. In the second study, all subject models are implemented in PSCAD/EMTDC and interfaced with its computational engine. Therein, a modified IEEE 39-bus system is considered where the number of synchronous machines has been increased from 10 to 18 to better assess the computational difference between the subject models. For the evaluation of numerical accuracy, the 2-norm (cumulative) error over the entire simulation interval is assumed as the accuracy index [65]. To evaluate accuracy for three-phase quantities, the average of single-phase errors is calculated and reported [68], [62]. 5.6.1 Single-Machine Infinite-Bus System The same salient-pole synchronous machine with parameters in Appendix D is considered here. A single-line-to-ground (SLG) fault is considered, which places the machine under severe asymmetrical conditions. The machine starts from a steady state under nominal excitation, terminal voltages, and nominal mechanical input torque. At t = 0.03s, a SLG fault occurs at the machine terminal and is removed after 200 ms. This study pushes the magnetizing flux in and out of the saturation region, which tests the models’ accuracy. The reference solution is obtained using a saturable state-space qd0 model [64] implemented in Simulink and executed using fixed step solver od4 with a very small time step of 1µs for high accuracy. The transient response to the SLG fault and its subsequent removal observed in stator currents iabcs, electromagnetic torque Te, and machine flux λm are demonstrated in Figure 5.2. As shown in Figure 5.2, the machine starts from saturation region (based on the saturation characteristic in Appendix D), and during the transient, it undergoes abrupt changes in the saturation levels [see Figure 5.2(e)]. For comparison purposes, all subject models are run with diff- 92 Figure 5.2 Transient response of several variables to SLG fault and its removal: (a) phase a current ias, (b) phase b current ibs, (c) phase c current ics, (d) electromagnetic torque Te, and (e) main flux λm. 93 Figure 5.3 Magnified views from Figure 5.2 obtained for models executed with 1 ms time step: (a) phase c stator current ics, and (b) electromagnetic torque Te. -erent time-step sizes ranging from 50 µs to 1ms. The magnified views of the stator current phase c after the fault removal and the electromagnetic torque after the fault occurrence are illustrated in Figure 5.3. As seen in Figure 5.3, the variable-parameter VBR and PD models predict the transient with high accuracy. Due to the current correction scheme in CPPD-q-II model, it obtains similar accuracy to the PD model testifying to the efficacy of the scheme in Section 5.2 for mitigating the current prediction error in (3.51). The CPVBR-I and CPPD-qd-I models obtain less accurate yet comparable results to other models, with the former having slightly better performance. 94 Figure 5.4 The 2-norm error of stator currents ε(iabcs) and electromagnetic torque ε(Te) as a function of time step solver for the subject SLG fault. To provide a quantitative analysis of accuracy, the 2-norm errors of stator currents ε(iabcs) and electromagnetic torque ε(Te) are calculated for different time steps and shown in Figure 5.4. As seen in Figure 5.4, all models converge to the same solution as the time step decreases, which is expected. Consistent with Figure 5.3, the VBR and, to a lesser extent, the PD and CPPD-q-II models achieve the highest accuracy among the subject models. It is also noted that, up to 500 µs, the prediction error is insignificant as the discretization error still dominates the overall numerical error. Only for time-step sizes larger than 500 µs, the slopes of the error trajectory in CPVBR-I and CPPD-qd-I models become steeper than the slopes of variable-parameter VBR and PD models. This may be because the current prediction error in these models starts to dominate the discretization error. 95 Figure 5.5 One-line diagram of the modified IEEE-39 bus system comprising 18 synchronous generators (G1 is modeled as an ideal source). 5.6.2 Multimachine System To further evaluate the numerical properties of the proposed CPPD models, a modified IEEE 39-bus system [69] is considered herein, as illustrated in Figure 5.5. To prevent the numerical error stemming from phase shift in the solutions, the generator connected to bus 39 (G1) is modelled as an ideal three-phase source with a fixed frequency [68]. The number of remaining machines (i.e., 9 in the original network [69]) has been increased to 18, as illustrated in Figure 5.5. 23039253718172726 28 29383451624151416789311211 131032343320363522212319G1G19G14, G15, G16G10, G11G12, G13G6, G7 G8, G9G4, G5G2, G3Bus37iG17, G1896 Figure 5.6 The subject three-phase-to-ground fault study obtained with ∆t = 500μs for phase c fault current contributed by the generators connected to bus 37 iBus37,c: (a) the overall transient and (b) the magnified view. Specifically, each of the generators in the original network has been replaced with one or more salient-pole machines from Section 5.6.1 to match the required MVA at each bus best. The system is initially at a steady-state operating point specified in [69]. At t = 0.04s, a three-phase-to-ground fault occurs at bus 25 [see Figure 5.5], and it is removed after 200 ms. For comparison purposes, the reference solution is obtained by the saturable PD model [42] with ∆t = 10 μs, which has been considered as an adequately small time step [30], [62], [68]. The solution trajectory for the fault current phase a contributed by the set of generators connected to bus 37 (i.e., G17 to G18) obtained by the subject models with ∆t = 500 μs is demonstrated in Figure 5.6(a). 97 Table 5.4 Average 2-Norm Error of Fault Current Contributed by the Set of the Synchronous Generators Connected to Bus 37 for Different Time-Step Sizes Expressed in Percentage (%) Time-step size PD VBR CPVBR-I CPPD-qd-I CPPD-q-II 100 µs 1.1228 1.1289 1.1307 1.1256 1.1237 500 µs 6.2430 5.9157 5.9291 6.1933 6.2324 The magnified view of the same current immediately after the fault removal is shown in Figure 5.6(b). As seen in Figure 5.6(a)‒(b), all subject models achieve similar accuracy, testifying to the efficacy of the proposed models. For a more detailed comparison, the average of the 2-norm single-phase errors of the fault current contributed by the generators G17 and G18 are provided in Table 5.4. Based on Table 5.4, with a small time step of ∆t = 100μs, all models offer a similar accuracy. As the time step increases, the PD-based models’ error starts to deviate slightly from that of the VBR-based model [27]. Additionally, since the error associated with solving the network equations dominates the simulation results (not the error associated with the machine models), the current correction feature augmented in CPPD-q-II is less visible here. Lastly, the computational performance of the subject models is evaluated here for the 0.5 s transient study. The CPU times for this study obtained with ∆t = 100 μs are presented in Table 5.5. As seen in Table 5.5, having a constant-parameter interfacing circuit and not refactorizing the network’s conductance matrix at each time step leads to around 80% faster simulation speed (e.g., 2.127 s for the variable-parameter PD vs 0.426 s for the CPPD-qd-I model). Consistent with the findings in Table 5.3, the VBR-based models are more expensive than the PD-based models, both in variable-parameter (i.e., 2.282 s vs 2.127 s) and the constant-parameter implementations (i.e., 98 Table 5.5 CPU Time Required for Simulation of the Modified IEEE 39-Bus System with Different Models for 0.5s Short-Circuit Transient Study Models PD VBR CPVBR-I CPPD-qd-I CPPD-q-II CPU Time (s) 2.127 2.282 0.466 0.426 0.433 0.466 s vs 0.433 s and 0.426 s). As a result, it is concluded that the CPPD-qd-I model achieves the fastest simulation speed (by being 3% faster than the CPPD-q-II model), making it a perfect candidate for simulating practical systems in offline and real-time simulators. 99 Contributions, Anticipated Impact, and Future Work 6.1 Contributions The research presented in this thesis has been focused on developing numerically efficient models of synchronous machines for their implementation in EMTP-type programs. Synchronous machine models have a significant presence in all power system studies conducted by the industry to ensure optimum, safe, and reliable operation of existing and future grids. Since synchronous machines always operate at some level of saturation, modelling saturation is necessary to obtain high-fidelity results. From a numerical perspective in EMTP-type programs, the interfacing circuit of synchronous machine model should have the following properties: 1) have a direct interface with the program, i.e., have a simultaneous solution of machine-network’s equations; and 2) have a constant-parameter interface so that the entire network’s conductance matrix becomes time-invariant and does not need to be reformulated after the initialization stage. Although several algebraically-equivalent state-of-the-art models have been proposed in the literature (e.g., qd0, PD, and VBR), their interfacing circuit is either rotor-position- and/or saturation-dependent. This thesis fulfills this gap by proposing two groups of synchronous machine models based on VBR and PD formulations that achieve better computational efficiency and comparable accuracy to other state-of-the-art models. The main features of the new models and the contributions of each research objective are summarized below. In Chapter 3, Objective I of this thesis is achieved by proposing two new synchronous machine models (CPVBR-I and CPVBR-II) with main flux saturation by re-formulating the original varying-parameter saturable VBR model to have a fully constant-parameter interfacing matrix, which requires prediction of the stator current at each time step. To compensate for this 100 prediction, a novel current-correction scheme is augmented in the first model (CPVBR-I) and a new model is formulated, referred to as the CPVBR-II in Subsection 3.4.2. Counting the number of flops and nonlinear functions in Section 3.5, the CPVBR models require fewer flops compared to the VBR model (i.e., 421 for the CPVBR-I and 487 for the CPVBR-II models vs 551 for the VBR model). Therein, it was also shown that the CPVBR models offer comparable costs to the PD model, which requires 397 flops. In terms of CPU time per step, the CPVBR-I and CPVBR-II models require 0.387 µs and 0.412 µs, thus outperforming the VBR model with 0.433 µs, and being marginally slower than the PD model requiring 0.374 µs (as demonstrated in Section 3.5). In the single-machine case study in Subsection 3.6.1, it was shown the CPVBR-I model achieves a similar level of accuracy compared to the variable-parameter PD and VBR models up to a time-step size of 500 µs. Additionally, the CPVBR-II achieved comparable accuracy to the VBR model over a wide range of time steps, thanks to its current-correction feature. A three-phase short-circuit study on the IEEE 9-bus system consisting of three synchronous generators was implemented in PSCAD to validate the proposed models and their computational efficiency. As illustrated in Subsection 3.6.2, the CPVBR-I achieves the fastest response by a factor of 1.72, 1.71, 1.03, and 1.13 compared to the VBR, PD, CPVBR-II, and the PSCAD’s built-in qd0 model, where the last two models also have a constant-parameter interfacing circuit. This outstanding performance makes the new CPVBR models particularly appealing to the real-time simulation environment where the network matrix refactorization is not desirable/applicable. Objective II is addressed in Chapters 4 and 5 by introducing the variable-parameter PD model in Chapter 4, followed by the proposed formulations in Chapter 5 for deriving the CPPD models. It is shown that the original PD model with main flux saturation has similar accuracy to the state-of-the-art saturable VBR model and outperforms the PSCAD’s built-in qd0 model in 101 terms of numerical and saturation modelling accuracy, especially at large time-step sizes (e.g., up to 1 ms) [see Section 4.2.3]. Furthermore, based on the genesis of the interfacing circuit of the PD model in qd0 coordinates, three different approaches for making a constant-parameter interfacing matrix are introduced in Section 5.1, each resulting in a constant-parameter model, but with a different set of interfacing matrix and history terms, referred to as the CPPD-q, CPPD-d, and the CPPD-qd, respectively. Similar to the CPVBR models, the CPPD models also require the prediction of stator current at each time step. Hence, following the research path of Chapter 3, a current-correction scheme is also introduced in Section 5.2 that can be augmented in any of the proposed CPPD models to compensate for the prediction of the current error. As reported in Section 5.5, the CPPD models without the correction scheme (CPPD-q/d/qd-I models) require the least number of flops per time step (258 flops), compared to the original PD model with 397 flops, and the CPVBR-I model from Chapter 3 with 421 flops. It is also reported that the inclusion of the current correction scheme adds only 21 flops, and results in new models referred to as the CPPD-/q/d/qd-II. Very importantly, even with 21 added flops, the new CPPD-/q/d/qd-II models are demonstrated to be more efficient than the PD and CPVBR-I models. In terms of CPU per time step, the CPPD-q/d/qd-I models achieve the fastest CPU time compared to the CPPD-q/d/qd-II, the CPVBR-I, and the PD models (0.333 µs vs 0.349 µs, 0.387 µs, and 0.374 µs, respectively). The short-circuit study in the single-machine study in Section 5.6.1 demonstrates that the proposed CPPD-qd-I model predicts the subject transient with high accuracy similar to the variable-parameter VBR and PD models up to 500 µs, which is an excellent performance. It is also verified that the inclusion of the current correction scheme in CPPD-q-II can effectively minimize the numerical error, resulting in a numerical performance very similar to the original variable-parameter PD model. Lastly, in the multimachine case study, a modified version of IEEE 39-bus 102 system comprising 18 synchronous generators has been implemented in PSCAD/EMTDC as described in Subsection 5.6.2. It is illustrated that the CPPD-qd-I and CPPD-q-II models significantly outperform the variable-parameter VBR and PD models in terms of simulation speed by a factor of 4.9. The CPPD-qd-I and CPPD-q-II models also achieve faster simulation speeds than the CPVBR-I model by a factor of 1.08 and 1.09, respectively. It can be concluded that the CPPD-qd-I model achieves the fastest simulation speed while preserving its numerical accuracy in system-level case studies. To provide the reader with a perspective of the contributions of this thesis, Table 6.1 summarizes the numerical properties of the state-of-the-art synchronous machine models while comparing them to the models developed in Chapters 3 to 5. As demonstrated in Table 6.1, this thesis fulfills the gap in research and literature by introducing new models with main flux saturation that can achieve both direct and constant-parameter interfacing circuits in EMTP-type programs. 103 Table 6.1 Classification of synchronous machine models with/without flux saturation based on interfacing with an external network (RPD: rotor-position-dependent; SSD: saturation-segment-dependent) Models Interfacing type Interfacing conductance matrix/parameter Magnetically-linear Main flux saturation Classical qd0 models [4], [8], [22], [41] Direct Constant Partially Constant (only SSD) [23], [24] Indirect Constant Constant PD models [18] ,[27] Direct Variable (RPD) - [28], [42](1) Direct Variable (RPD) Variable (RPD/SSD) [29] ,[30] Direct Constant - [44] Direct Constant Partially Constant (only SSD) Proposed in Chapter 3 Direct Constant Constant VBR models [18] ,[27] Direct Variable (RPD) - [43] Direct Variable (RPD) Variable (RPD/SSD) [45] Direct Constant Partially Constant (only SSD) Proposed in Chapter 5 Direct Constant Constant (1) it is introduced in Chapter 4. 6.2 Anticipated Impact The constant interfacing circuit and very small number of floating point operations per time step make the proposed models desirable alternatives. Therefore, it is envisioned that the newly proposed CPPD machine models will replace the existing machine models in commercial offline EMTP-type programs such as PSCAD/EMTDC and real-time simulators such as RTDS. The proposed models are also specially formulated to have a similar set of interfacing variables and parameters as other alternative or existing models while providing significant computational benefits. These benefits may result in extending the range of applications of the EMT simulators to larger systems that can be handled using the existing or smaller computational hardware resources. The new models and enhanced EMT tools for the power industry will also translate into tangible benefits. Companies that are performing a large number of offline transient studies will 104 be able to perform the same number of studies but in a fraction of the time, be able to extend the size of systems and the level of details that are considered, increase the fidelity and accuracy of the studies, etc. Companies that are performing real-time transient studies will be able to utilize the simulation hardware better. 6.3 Future Work The research presented in this thesis may be continued in several directions. Due to limited space, only four possible directions that may be considered and/or currently undertaken by other students in our research group are mentioned here: Firstly, several constant-parameter modelling approaches have been developed for single-rate simulations in Chapters 3 and 5 of this thesis. Extending these formulations to multi-rate applications is also desirable, where instead of one time-step size (i.e., single-rate) for the entire model, different time steps are used for the electrical and mechanical subsystems. This can further help to boost the computational savings of the proposed CPVBR and CPPD models [70]. To derive these new models, one would start by considering the multi-rate CPVBR induction machine model [71]. Then, the methodologies laid out in Chapters 3 and 5 can be adopted to obtain constant-parameter variations. Secondly, the conventional three-phase synchronous machines are considered in this thesis due to their dominant presence in today’s power grids. However, with the latest advancements in power electronics and renewable energy technologies, multiphase synchronous machines (i.e., machines with more than three phases) are becoming considered in applications with power electronic converters. So far, efficient state-variable-based VBR and CPVBR models of 5 and 6 phases synchronous machines have been proposed in the literature [60], [72]. However, the EMTP-105 type version of those models is yet to be derived, and the methodologies introduced in Chapters 3 and 5 can be further utilized and extended to machines with more than three electrical phases. Thirdly, the synchronous machine models presented in the thesis only considered the fundamental component of the airgap flux (i.e., main flux saturation), which is acceptable for many system-level studies. However, to increase the fidelity of general-purpose machine models, airgap flux harmonics may also be considered. As a continuation of the CPVBR SVB model of the five-phase synchronous machine [73], the future research task will develop the CPVBR and CPPD models of three-phase and multiphase machines considering the airgap flux harmonics. Lastly, machine models considered in the thesis are mainly valid for operation around the power systems’ nominal frequency and up to several kHz only. However, with the advent of renewable energy resources connected through power electronic converters, higher frequency components are also becoming important during transients. Also, in machine-drive applications (such as electric cars), the machine can operate at several times the nominal grid frequency. To have the required accuracy in computer studies, the new models should also be able to predict high-frequency phenomena. In this regard, a universal wide-band induction machine model that can predict both high- and low-frequency range efficiently is currently being developed in our research group for the SVB packages using the general-purpose VBR-type models. The research approach laid out in Chapters 3 and 5 can be extended for the formulation of the wind-band constant-parameter VBR and PD models with a validity of up to several MHz. 106 Bibliography [1] P. S. Kundur and O. P. Malik, Power system stability and control. McGraw-Hill Education, 2022. [2] US Energy Information Adminstration. http://www.eia.gov. [3] P. C. Krause, O. Wasynczuk, S. D. Sudhoff, and S. D. Pekarek, Analysis of electric machinery and drive systems, 3rd ed. John Wiley & Sons, 2013. [4] H. W. Dommel, EMTP Theory Book. Vancouver, BC, Canada: Microtran Power System Analysis Corporation, May 1992. [5] “RT-LAB”, OPAL-RT TECHNOLOGIES, Montreal, QC, Canada. [6] J. Mahseredjian, V. Dinavahi, and J. A. Martinez, “An overview of simulation tools for electromagnetic transients in power systems,” in Proc. IEEE Power Energy Soc. Gen. Meeting, Tampa, FL, USA, Jun. 24-28 2007, pp. 1‒6. [7] MATLAB/Simulink: Dynamic System Simulation Software, User's Manual, MathWorks Inc., Natick, MA, USA, 2023. [Online]. Available: http://www.mathworks.com. [8] Piecewise Linear Electrical Circuit Simulation (PLECS), User's Manual, Version 4.7.1, Plexim GmbH, Zurich, Switzerland, 2023. [Online]. Available: http://www.plexim.com. [9] Simscape Electrical (SimPowerSystems): Model and Simulate Electrical Power Systems, User’s Guide, MathWorks Inc., Natick, MA, USA, 2023. [Online]. Available: http://www.mathworks.com. [10] PSCAD/EMTDC, Ver. 5.0.1, On-Line Help Manitoba HVDC Research Centre and RTDS Technologies Inc., Winnipeg, MB, Canada, 2023. [Online]. Available: http://www.pscad.com. 107 [11] Electromagnetic Transient Program, EMTP-RV, Ver. 4.2.1, POWERSYS Solutions Inc., Le Puy-Sainte-R´eparade, France, 2021. [Online]. Available: http://www.emtp.com. [12] MicroTran Reference Manual, MicroTran Power System Analysis Corp., Vancouver, BC, Canada. 1997. [Online]. Available: http://www.microtran.com. [13] Alternative Transients Programs, ATP-EMTP, ATP User Group, 2007. [Online]. Available: http://emtp.org. [14] RTDS Simulator Software, RTDS Technologies Inc., Winnipeg, MB, Canada, 2023. [Online]. Available: http://www.rtds.com. [15] O. Wasynczuk and S. Sudhoff, “Automated state model generation algorithm for power circuits and systems,” IEEE Trans. Power Syst., vol. 11, no. 4, p. 1 51‒1 56, Nov. 1 6. [16] J. V. Jatskevich, “A state selection algorithm for the automated state model generator,” Ph.D. dissertation, Purdue University West Lafayette, IN, USA, 1999. [17] C. Dufour, J. Mahseredjian, and J. Bélanger, “A combined state-space nodal method for the simulation of power system transients,” IEEE Trans. Power Del., vol. 26, no. 2, p. 28‒ 35, Apr. 2011. [18] U. Karaagac, J. Mahseredjian, O. Saad, and S. Dennetière, “Synchronous machine modeling precision and efficiency in electromagnetic transients,” IEEE Trans. Power Del., vol. 26, no. 2, pp. 1072–1082, Apr. 2011. [19] X. Cao, A. Kurita, H. Mitsuma, Y. Tada, and H. Okamoto, “Improvements of numerical stability of electromagnetic transient simulation by use of phase‐domain synchronous machine models,” Elect. Eng. Jpn., vol. 128, no. 3, p. 53‒62, Apr. 1999. 108 [20] L. Wang, J. Jatskevich, and H. W. Dommel, “Re-examination of synchronous machine modeling techniques for electromagnetic transient simulations,” IEEE Trans. Power Syst., vol. 22, no. 3, p. 1221‒1230, Aug. 2007. [21] L. Wang et al., “Methods of interfacing rotating machine models in transient simulation programs,” IEEE Trans. Power Del., vol. 25, no. 2, pp. 891–903, Apr. 2010. [22] V. Brandwajn, “Synchronous generator models for the simulation of electromagnetic transients,” Ph.D. disssertation, Dept. Elect. Eng., Univ. British Columbia, Vancouver, Canada, 1977. [23] D. Woodford, A. Gole, and R. Menzies, “Digital simulation of DC links and AC machines,” IEEE Trans. Power App. Syst., no. 6, pp. 1616–1623, Jun. 1983. [24] A. Gole, R. Menzies, H. Turanli, and D. Woodford, “Improved interfacing of electrical machine models to electromagnetic transients programs,” IEEE Trans. Power App. Syst., no. 9, pp. 2446–2451, Sep. 1984. [25] H. W. Dommel, “Nonlinear and time-varying elements in digital simulation of electromagnetic transients,” IEEE Trans. Power App. Syst., no. 6, p. 2561‒2567, Nov./Dec. 1971. [26] H. K. Lauw and W. S. Meyer, “Universal machine modeling for the representation of rotating electric machinery in an electromagnetic transients program,” IEEE Trans. Power App. Syst., no. 6, pp. 1342–1351, Jun. 1982. [27] L. Wang and J. Jatskevich, “A voltage-behind-reactance synchronous machine model for the EMTP-type solution,” IEEE Trans. Power Syst., vol. 21, no. 4, pp. 1539–1549, Nov. 2006. [28] J. R. Marti and K. W. Louie, “A phase-domain synchronous generator model including saturation effects,” IEEE Trans. Power Syst., vol. 12, no. 1, pp. 222–229, Feb. 1997. 109 [29] 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. [30] Y. Xia, Y. Chen, Y. Song, S. Huang, Z. Tan, and K. Strunz, “An efficient phase domain synchronous machine model with constant equivalent admittance matrix,” IEEE Trans. Power Del., vol. 34, no. 3, pp. 929–940, Jun. 2019. [31] E. Levi, “State-space d-q axis models of saturated salient pole synchronous machines,” IEE Proc. Electr. Power Appl., vol. 145, no. 3, pp. 206–216, Mar. 1998. [32] E. Levi, “Saturation modelling in d-q axis models of salient pole synchronous machines,” IEEE Trans. Energy Convers., vol. 14, no. 1, pp. 44–50, Mar. 1999. [33] P. M. Anderson and A. A. Fouad, Power system control and stability. John Wiley & Sons, 2008. [34] R. Harley, D. Limebeer, and E. Chirricozzi, “Comparative study of saturation methods in synchronous machine models,” in IEE Proc. B—Elect. Power Appl., Jan. 1980, vol. 1, no. 127, pp. 1–7. [35] F. De Mello and L. Hannett, “Representation of saturation in synchronous machines,” IEEE Trans. Power Syst., vol. 1, no. 4, pp. 8–14, Nov. 1986. [36] P. Vas, K. Hallenius, and J. Brown, “Cross-saturation in smooth-air-gap electrical machines,” IEEE Trans. Energy Convers., no. 1, pp. 103–112, Mar. 1986. [37] A. El-Serafi, A. Abdallah, M. El-Sherbiny, and E. Badawy, “Experimental study of the saturation and the cross-magnetizing phenomenon in saturated synchronous machines,” IEEE Trans. Energy Convers., vol. 3, no. 4, pp. 815–823, Dec. 1988. 110 [38] S.-A. Tahan and I. Kamwa, “A two-factor saturation model for synchronous machines with multiple rotor circuits,” IEEE Trans. Energy Convers., vol. 10, no. 4, pp. 609–616, Dec. 1995. [39] E. Levi and V. A. Levi, “Impact of dynamic cross-saturation on accuracy of saturated synchronous machine models,” IEEE Trans. Energy Convers., vol. 15, no. 2, pp. 224–230, Jun. 2000. [40] R. Wamkeue, I. Kamwa, F. Baetscher, and J. El Hayek, “A new and efficient approach for analysis of a saturated synchronous generator under the load rejection test,” Elect. Power Compon. Syst., vol. 34, no. 5, pp. 539–563, 2006. [41] V. Brandwajn, “Representation of magnetic saturation in the synchronous machine model in an electro-magnetic transients program,” IEEE Trans. Power App. Syst., no. 5, pp. 1996–2002, Sep./Oct. 1980. [42] E. Mostajeran and J. Jatskevich, “Phase-Domain Synchronous Machine Model with Main Flux Saturation for EMTP-type Solution,” in Proc. Int. Symp. INFOTEH-JAHORINA (INFOTEH), East Sarajevo, Bosnia and Herzogovina, Mar. 15‒17, 2023, pp. 1‒6. [43] 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. [44] U. Karaagac, J. Mahseredjian, and O. Saad, “An efficient synchronous machine model for electromagnetic transients,” IEEE Trans. Power Del., vol. 26, no. 4, pp. 2456–2465, Oct. 2011. 111 [45] U. Karaagac, J. Mahseredjian, I. Kocar, and O. Saad, “An efficient voltage-behind-reactance formulation-based synchronous machine model for electromagnetic transients,” IEEE Trans. Power Del., vol. 28, no. 3, pp. 1788–1795, Jul. 2013. [46] M. Chapariha, F. Therrien, J. Jatskevich, and H. W. Dommel, “Explicit formulations for constant-parameter voltage-behind-reactance interfacing of synchronous machine models,” IEEE Trans. Energy Convers., vol. 28, no. 4, pp. 1053–1063, Dec. 2013. [47] R. H. Park, “Two-reaction theory of synchronous machines generalized method of analysis-part I,” Trans. Amer. Instit. Electr. Eng., vol. 48, no. 3, p. 716‒727, Jul. 1 2 . [48] U. M. A. a. R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Philadelphia, PA, USA: SIAM, 1998. [49] S. Pekarek, O. Wasynczuk, and H. 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. [50] M. Chapariha, F. Therrien, J. Jatskevich, and H. W. Dommel, “Constant-parameter circuit-based models of synchronous machines,” IEEE Trans. Energy Convers., vol. 30, no. 2, p. 441‒452, Jun. 2015. [51] M. Chapariha, F. Therrien, J. Jatskevich, and H. W. Dommel, “Detailed modeling of telecom AC-DC backup generation power systems using voltage-behind-reactance synchronous machine model,” in Proc. IEEE Int. Telecomm. Energy Conf. (INTELEC), Vancouver, BC, Canada, Sep. 28‒Oct. 2 2014: IEEE, p. 1‒6. [52] M. Chapariha, L. Wang, J. Jatskevich, H. W. Dommel, and S. D. Pekarek, “Constant-Parameter RL-Branch Equivalent Circuit for Interfacing AC Machine Models in State-112 Variable-Based Simulation Packages,” IEEE Trans. Energy Convers., vol. 27, no. 3, p. 634‒645, Sep. 2012. [53] W. Zhu, S. Pekarek, J. Jatskevich, O. Wasynczuk, and D. Delisle, “A model-in-the-loop interface to emulate source dynamics in a zonal DC distribution system,” IEEE Trans. Power Electron., vol. 20, no. 2, p. 438‒445, Mar. 2005. [54] F. Therrien, M. Chapariha, and J. Jatskevich, “Pole selection procedure for explicit constant-parameter synchronous machine models,” IEEE Trans. Energy Convers., vol. 29, no. 3, p. 7 0‒7 2, Jun. 2014. [55] D. C. Aliprantis, O. Wasynczuk, and C. D. R. Valdez, “A voltage-behind-reactance synchronous machine model with saturation and arbitrary rotor network representation,” IEEE Trans. Energy Convers., vol. 23, no. 2, p. 4 ‒508, May. 2008. [56] A. M. Cramer, B. P. Loop, and D. C. Aliprantis, “Synchronous machine model with voltage-behind-reactance formulation of stator and field windings,” IEEE Trans. Energy Convers., vol. 27, no. 2, p. 3 1‒402, Mar. 2012. [57] M. Mobarra, M. Rezkallah, and A. Ilinca, “Variable speed diesel generators: Performance and characteristic comparison,” Energies, vol. 15, no. 2, p. 592, Jan. 2022. [58] T. Theubou, R. Wamkeue, and I. Kamwa, “Dynamic model of diesel generator set for hybrid wind-diesel small grids applications,” in Proc. IEEE Canadian Conf. Elec. Comp. Eng. (CCECE), Oct. 2012, p. 1‒4. [59] S. Ebrahimi, N. Amiri, L. Wang, and J. Jatskevich, “Efficient modeling of six-phase PM synchronous machine-rectifier systems in state-variable-based simulation programs,” IEEE Trans. Energy Convers., vol. 33, no. 3, p. 1557‒1570, Sep. 2018. 113 [60] N. Amiri, S. Ebrahimi, J. Jatskevich, and H. W. Dommel, “Saturable and decoupled constant-parameter vbr model for six-phase synchronous machines in state-variable simulation programs,” IEEE Trans. Energy Convers., vol. 34, no. 4, pp. 1868–1880, Dec. 2019. [61] F. Therrien, M. Chapariha, and J. Jatskevich, “Constant‐parameter synchronous machine model including main flux saturation,” IET Elect. Power Appl., vol. 10, no. 6, pp. 477–487, Jul. 2016. [62] 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. [63] C. Van Loan, Computational frameworks for the fast Fourier transform. Philadelphia, PA, USA: SIAM, 1992. [64] F. Therrien, L. Wang, J. Jatskevich, and O. Wasynczuk, “Efficient explicit representation of AC machines main flux saturation in state-variable-based transient simulation packages,” IEEE Trans. Energy Convers., vol. 28, no. 2, pp. 380–393, Jun. 2013. [65] G. H. Golub and C. F. Van Loan, Matrix computations, 4th ed. Baltimore, MD, USA: Johns Hopkins Univ. Press, 2013. [66] P. Demetriou, M. Asprou, J. Quiros-Tortos, and E. Kyriakides, “Dynamic IEEE test systems for transient analysis,” IEEE Syst. J., vol. 11, no. 4, pp. 2108–2117, Jul. 2015. [67] K. Corzine, B. T. Kuhn, S. Sudhoff, and H. Hegner, “An improved method for incorporating magnetic saturation in the qd synchronous machine model,” IEEE Trans. Energy Convers., vol. 13, no. 3, pp. 270–275, Sep. 1998. [68] E. Mostajeran, N. Amiri, and J. Jatskevich, “Constant-parameter voltage-behind-reactance synchronous machine models considering main flux saturation for EMTP-type programs,” 114 IEEE Trans. Energy Convers. (under review), pp. 1‒10, submitted on Feb. 22, 2023: TEC-00189-2023. [69] M. A. Pai, Energy function analysis for power system stability. Norwell, MA, USA: Kluwer, 1989. [70] I. Kamwa and P. Viarouge, “On equivalent circuit structures for empirical modeling of turbine-generators," IEEE Trans. Energy Convers., vol. , no. 3, p. 57 ‒5 2, Sep. 1 4. [71] F. Therrien and J. Jatskevich, “Multirate emtp-type induction machine models,” IEEE Trans. Energy Convers., vol. 31, no. 3, p. 1142‒1152, Sep. 2016. [72] N. Amiri, S. Ebrahimi, Y. Huang, J. Jatskevich, and S. D. Pekarek, “Constant-parameter voltage-behind-reactance modeling of five-phase synchronous machines with air-gap flux harmonics,” IEEE Trans. Energy Convers., vol. 35, no. 1, p. 11 ‒12 , Sep. 201 . [73] N. Amiri, S. Ebrahimi, Y. Huang, J. Jatskevich, and S. D. Pekarek, “Constant-parameter voltage-behind-reactance modeling of five-phase synchronous machines with air-gap flux harmonics,” IEEE Trans. Energy Convers., vol. 35, no. 1, p. 11 ‒12 , Mar. 2020. 115 Appendix A: Reference Frame Transformation Matrices Used in Chapter 2 The Park’s and inverse of Park’s transformation matrices used for derivation of qd0 models are defined respectively as ( )( )( )2 2cos cos cos3 32 2 2sin sin sin3 3 31 1 12 2 2r r rs r r r r           − +            = − +           K , (A.1) ( )( ) ( )1cos sin 12 2cos sin 13 32 2cos sin 13 3r rs r r rr r      −       = − −            + +        K . (A.2) 116 Appendix B: Parameters for the Case Study of Subsection 2.2.1 Three-phase WFSG parameters [9]: Rated ac voltage: 2.4 kV, rated dc voltage: 3.24 kV, rated power 3.125 MVA, poles: 4, rated speed: 1800 rpm, inertia time-constant: 1.07s. 𝑟𝑠 = 0.0036 pu, 𝑋𝑙𝑠 = 0.052 pu, 𝑋𝑚𝑑 = 1.508 pu, 𝑋𝑚𝑞 = 1.008 pu, 𝑟𝑘𝑞 = 0.0616 pu, 𝑋𝑙𝑘𝑞 = 0.1532 pu, 𝑟𝑓𝑑 = 0.00245 pu, 𝑋𝑙𝑓𝑑 = 0.2911 pu, 𝑟𝑘𝑑 = 0.0159 pu, 𝑋𝑙𝑘𝑑 = 0.2563 pu. DC filter parameters: 𝐿𝑓 = 1 mH, 𝐶𝑓 = 10 mF. Voltage controller parameters [20]: 𝐾𝑝 = 500, 𝐾𝐼 =10, 𝑇3 = 0.02s. Speed controller parameters [4]: 𝐾𝑝 = 6.09, 𝐾𝐼 =5.9, 𝑇1 =0.1s, 𝑇2 =0.1s. 117 Appendix C: Additional Scalers and Matrices Introduced in Chapter 3 The matrices in (3.13) are defined as 1 22 21, , ,2 2 2 21 2diag ,kq kq fd kdj mq j md jlkq lkq lfd lkdr r r rL LL L L L       = + +              A , (C.1) , ,1, 2,2,, ,1, 2,1 2r md j r md jq j q jlfd lkdjr mq j r mq jd j d jlkq lkqL La aL LL La aL L      =   − −  A , (C.2) 2, 1 , , 21, 2 21 1 1 21mq j kq mq j mq j kqq jlkq lkq lkq lkqL r L L raL L L L   = − +   , (C.3) 2, 2 , , 12, 2 22 2 2 11mq j kq mq j mq j kqq jlkq lkq lkq lkqL r L L raL L L L   = − +   , (C.4) 2, , ,1, 2 21md j fd md j md j kdd jlfd lfd lfd lkdL r L L raL L L L   = − +   , (C.5) 2, , ,2, 2 21md j kd md j md j fdd jlkd lkd lkd lfdL r L L raL L L L   = − +  , (C.6) 1, ,,23,, ,,1,2(1,1)1(2,2)j mq jr md jjmd j r mq jD jjLLSF SFL LLSF SF   −  =  −  AAA. (C.7) The matrices in (3.15) are defined as 1 , 2 11,2 1 1 ,q jjd j =   B 0B0 B, 2 , 2 22,2 2 2 ,q jjd j =   B 0B0 B, (C.8) 118 1 21 , ,1 2Tkq kqq j mq jlkq lkqr rLL L =    B , 1 , ,Tfd kdd j md jlfd lkdr rLL L =    B , (C.9) 1 , 1 ,1 1 1 22 ,2 , 2 ,1 2 2 211kq mq j kq mq jlkq lkq lkq lkqq jkq mq j kq mq jlkq lkq lkq lkqr L r LL L L Lr L r LL L L L   −     =     −     B , (C.10) , ,2 ,, ,11fd md j fd md jlfd lfd lfd lkdd jkd md j md jkdlfd lkd lkd lkdr L r LL L L Lr L LrL L L L   −     =     −    B , (C.11) 3 , 2 13,2 1 3 ,q jjd j =   B 0B0 B,  4 0 0 1 0T=B , (C.12) ( )123 , , 1 ,q j D j q jL SF−=B B , 13 , , 1 ,d j D j d jL−=B B . (C.13) The matrices in (3.20) are defined as ( )110, 4 4 2,j jk−−= −C I B , (C.14) ( )1 11, 0, 1, 2, 0, 4 4 2,,j j j j j jk k− −= = +C C B C C I B , (C.15) 1 13, 0, 3, 4, 0, 4,2 , 2j j j j j jk k− −= =C C B C C B . (C.16) The matrices in (3.21) and (3.22) are defined as 11, 12,1, 2, 1,21, 22,j j rj j j jj r jD DD D = = +  D A A C , (C.17) 1, 2, 1, 2, 2, 2,,j j j j j j= =E A C E A C , (C.18) 119 3, 3, 2, 3, 4, 2, 4,,j j j j j j j= + =E A A C E A C . (C.19) The matrices in (3.55) are defined 2 1 2 1 211 2 2 1 2h-const 11 2 1 2 212 2 2 2 223 2 2 2 2 22 2 2 2 2h h h h hhh h h h hhh h h h hhR R R R RRR R R R RRR R R R RR + − + − +   = − + + − +   − + − + +  R , (C.20) ( ) ( )1 ,1 ,1 11,1 22,10.5 0.5h s ls mq mdR r kL k L L D D = − − + + + , 2h s lsR r kL= − , (C.21) ( )( )( ), 12,0-var, 21, ,000 0 0qh j j rqd dh j j r h jR D tt D t R   =    R , (C.22) , , 11, 1qh j mq j j hR kL D R = − + − , , , 22, 1dh j md j j hR kL D R = − + − . (C.23) The scalers in (3.58) are defined as ( )1, , 22,j s ls md j jg r k L L D= + + + , ( )2, , 11,j s ls mq j jg r k L L D= + + + , (C.24) ( )13 s lsg r kL−= + . (C.25) 120 Appendix D: Parameters for the Case Study of Subsection 3.6.1, and Sections 4.2 and 5.6 Salient-pole synchronous machine parameters [3]: 325 MVA, 20 kV, 64 poles, 60 Hz, J = 35.1·106 J·s2, rs = 0.00234 Ω, Xls = 0.1478 Ω, rfd = 0.0005 Ω, Xlfd = 0.2523 Ω, rrkq2 = 0.01675 Ω, Xlkq2 = 0.1267 Ω, rkd = 0.01736 Ω, Xlkd = 0.1 70 Ω, Xmd = 0.8 8 Ω, Xmq = 0.4433 Ω. Machine saturation data* (peak values) [43], [64] λm (Wb) 10.8 22.5 31.0 35.3 38.2 mi (kA) 4.5143 9.4981 13.255 15.261 16.707 λm (Wb) 40.9 42.5 45.3 47.9 50.3 mi (kA) 18.198 19.210 21.340 23.652 25.930 * In Section 4.2, the second saturation operating point is not considered due to the limitation on the maximum number of data entry points (i.e., 10 including the origin) that can be inserted into the PSCAD’s built-in model. 121 Appendix E: Parameters for the Case Study of Subsection 3.6.2 Synchronous generator G2 parameters [66]: 270 MVA, 18 kV, 4 poles, 60 Hz, J = 0.063·106 J·s2, rs = 0.001 Ω, Xls = 0.1860 Ω, rfd = 0.0011 Ω, Xlfd = 0.1297 Ω, rrkq1 = 0.0099 Ω, Xlkq1 = 0.1151 Ω, rrkq2 = 0.6139 Ω, Xlkq2 = 0.0540 Ω, rkd = 0.0457 Ω, Xlkd = 0.0512 Ω, Xmd = 1.8540 Ω, Xmq = 1.7580 Ω. Machine saturation curve (peak values) [67]: λm (Wb) 7.80 23.39 31.97 34.31 36.65 mi (kA) 1.605 4.838 6.652 7.162 7.694 λm (Wb) 38.21 39.37 40.15 41.32 mi (kA) 8.087 8.599 9.942 12.230 Synchronous generator G3 parameters [66]: 125 MVA, 15.5 kV, 4 poles, 60 Hz, J = 0.034·106 J·s2, rs = 0.0077 Ω, Xls = 0.14 Ω, rfd = 0.00071 Ω, Xlfd = 0.2014 Ω, rrkq1 = 0.0131 Ω, Xlkq1 = 0.3 31 Ω, rrkq2 = 0.0186 Ω, Xlkq2 = 0.15 6 Ω, rkd = 0.0356 Ω, Xlkd = 0.2583 Ω, Xmd = 2.1 4 Ω, Xmq = 2.07 6 Ω. Machine saturation curve (peak values) [67]: λm (Wb) 6.71 20.14 27.53 29.54 31.56 mi (kA) 1.167 3.516 4.832 5.201 5.586 λm (Wb) 32.90 33.91 34.24 34.91 mi (kA) 5.869 6.259 6.954 8.575 122 Appendix F: Additional Scalers and Matrices Introduced in Chapter 4 The matrices in (4.2) are defined as ( ), ,, , , ,,, ,, ,, , , ,2 2cos(2 ) cos(2 ) cos(2 )2 3 2 32 2cos(2 ) cos(2 ) cos(2 )2 3 3 22 2cos(2 ) cos(2 ) cos(2 )2 3 2 3A j A jls A j B j r B j r B j rA jAs j r B r ls A B r B j rA j A jB j r B j r ls A j B j rL LL L L L LLLL L L L LL LL L L L L          + − − − − − − += − − − + − + − −− − + − − + − −L          , (F.1) ( )2,,13D jA jL SFL+= , ( )2,,13D jA jL SFL−= , (F.2) ( )2 22 2, ,2 2cos( ) cos( ) sin( ) sin( )2 2 2 2cos( ) cos( ) sin( ) sin( )3 3 3 32 2 2 2cos(2 ) cos( ) sin( ) sin( )3 3 3 3r r r rsr j r D j r r r rr r r rSF SFL SF SFSF SF                    = − − − −   + + + + L , (F.3) ( ) ( ), ,23Trs j r sr j r =L L , (F.4) 2 21 , ,2 2, 2 ,,, ,, ,0 00 00 00 0lkq D j D jD j lkq D jr jlfd D j D jD j lkd D jL SF L SF LSF L L SF LL L LL L L + + = + +  L . (F.5) 123 Appendix G: Additional Scalers and Matrices Introduced in Chapter 5 The matrices in (5.1)‒(5.3) are defined as ( )2, , ,diag , ,rs j ls D j ls D j lsL L SF L L L = + + L , (G.1) 2 2, ,, , ,0 00 00 0 0 0D j D jrsr j D j D jL SF L SFL L  =    L , , ,Tr rrs j sr j =  L L . (G.2) The matrices in (5.18) are defined as 1 2 2const 2 1 22 2 1a a ah a a aa a aR R RR R RR R R       −  =    R , where = , , q d qd , (G.3) 1 ,123qa s ls qR r kL Z = − − , 2 ,113qa qR Z = , (G.4) 1 ,123da s ls dR r kL Z = − − , 2 ,113da dR Z = , (G.5) ( )1 ,1 ,113qda s ls q dR r kL Z Z = − − + , ( )2 ,1 ,116qda q dR Z Z = + , (G.6) ( )var, 1, 2,diag 0h j b j b jR R  −   =  R , where = , , q d qd , (G.7) 1, ,1 ,qb j q q jR Z Z = − , 2, ,1 ,qb j q d jR Z Z = − , 1, ,1 ,db j d q jR Z Z = − , 2, ,1 ,db j d d jR Z Z = − , (G.8) ( )1, ,1 ,1 ,0.5qda j q d q jR Z Z Z  = + − , ( )2, ,1 ,1 ,0.5qda j q d d jR Z Z Z  = + − . (G.9) The Scalers in (5.25) are defined as ( )1, , , 00.5j q j d jG G G G= + − , ( )2, , , 00.25 2j q j d jG G G G= − + + . (G.10)
Thesis/Dissertation
eng
Electrical and Computer Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Jatskevich, Juri
Constant-parameter direct interfacing of saturable synchronous machine models for electromagnetic transient programs
Text