SHIFTED FREQUENCY ANALYSIS FOR EMTP SIMULATION OF POWER SYSTEM DYNAMICS by Peng Zhang B.Sc., Shandong University, China, 1996 M.Sc., Shandong University, China, 1999 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) March 2009 © Peng Zhang, 2009 Abstract Electromagnetic Transients Program (EMTP) simulators are being widely used in power system dynamics studies. However, their capability in real time simulation of power systems is compromised due to the small time step required resulting in slow simulation speeds. This thesis proposes a Shifted Frequency Analysis (SFA) theory to accelerate EMTP solutions for simulation of power system operational dynamics. A main advantage of the SFA is that it allows the use of large time steps in the EMTP solution environment to accurately simulate dynamic frequencies within a band centered around the fundamental frequency. The thesis presents a new synchronous machine model based on the SFA theory, which uses dynamic phasor variables rather than instantaneous time domain variables. Apart from using complex numbers, discrete-time SFA synchronous machine models have the same form as the standard EMTP models. Dynamic phasors provide envelopes of the time domain waveforms and can be accurately transformed back to instantaneous time values. When the frequency spectra of the signals are close to the fundamental power frequency, the SFA model allows the use of large time steps without sacrificing accuracy. Speedups of more than fifty times over the traditional EMTP synchronous machine model were obtained for a case of mechanical torque step changes. This thesis also extends the SFA method to model induction machines in the EMTP. By analyzing the relationship between rotor and stator physical variables, a phase-coordinate model with lower number of equations is first derived. Based on this, a SFA model is proposed as a general purpose model capable of simulating both fast transients and slow dynamics in induction machines. Case study results show that the SFA model is in excess of seventy times faster than the phase-coordinate EMTP model when simulating the slow dynamics. In order to realize the advantage of SFA models in the context of the simulation of the complete electrical network, a dynamic-phasor-based EMTP simulation tool has been developed. 11 Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii Acknowledgements xi Chapter 1 Chapter 2 Introduction 1 1.1 Background 1 1.2 Motivation 3 1.3 Contributions 5 Shifted Frequency Analysis 6 2.1 Analytic Signal and Hubert Transform 6 2.1 .1 7 2.2 2.3 Chapter 3 Shifted Frequency Analysis SFA-Based Network Component Models 10 2.2.1 Equivalent Circuits for RLC in the Shifted Frequency Domain 11 2.2.2 Options between Complex Arithmetic and Real Arithmetic 15 2.2.3 Transformer Model in the Shifted Frequency Domain 17 2.2.4 Load Models 20 Numerical Accuracy Analysis 27 Synchronous Machine Modelling Based on Shifted Frequency Analysis 33 3.1 Introduction 33 3.2 Voltage-behind-Reactance Synchronous Machine Model 33 111 Table of Contents 3.3 Synchronous Machine Modelling with SFA .38 3.3.1 Synchronous Machine Model Based on SFA 38 3.3.2 Discrete Time Model 39 3.3.3 Note on the Cylindrical-Rotor Machine Model 45 3.4 Simulation Results 47 3.4 Summary 57 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis 58 4.1 Introduction 58 4.2 Equivalent-Reduction Approach to Induction Machine Modelling in EMTP. 59 4.3 4.4 Chapter 5 4.2.1 Equivalent-Reduction (ER) to Stator Quantities 60 4.2.2 Induction Machine Modelling Based on ER Technique 60 4.2.3 Simulation Results 65 Induction Machine Modelling with SFA 71 4.3.1 Induction Machine Modelling Based on SFA 71 4.3.2 Discrete Time Model 72 4.3.3 Simulation Results 73 Summary 80 EMTP Implementation 81 5.1 Introduction 81 5.2 Program Structure 82 5.3 Test Cases 86 Chapter 6 Conclusions 102 6.1 Summary of Contributions 102 6.2 Future Research 103 Bibliography Appendix A 105 Machine Parameters A. 1 111 Synchronous Machine Parameters iv 111 Table of Contents A.2 Appendix B Induction Machine Parameters . Voltage behind Reactance Induction Machine Model 112 113 B. 1 Voltage Behind Reactance Model of Induction Machine 113 B. 2 Discrete Time VBR Model 116 B. 3 Free Acceleration Simulation of a 3-hp Induction Machine 118 B. 4 Dynamic Performance of Induction Machine during Mechanical Torque Appendix C Changes 118 Test Case Data 123 V List of Tables Table 2.1 Load Increased with the Time 23 Table 3.1 CPU Times for 4s Simulation in Case C 50 Table 3.2 Accuracy Comparison between SFA Model and VBR EMTP Model in Case C 56 Table 4.1 CPU Times for Simulations 71 Table 4.2 CPU Times for Simulations 75 Table 4.3 CPU Times for Simulations 79 vi List of Figures Figure 2.1(a) M Coupled Inductances; (b) EMTP Equivalent Circuit; (c) SFA Equivalent Circuit 9 Figure 2.2 (a) Series Connection of M-phase R and L; (b) EMTP Equivalent Circuit; (c) SFA Equivalent Circuit 13 Figure 2.3 Linear Test Case 13 Figure 2.4 Current Flowing through Branch 1-2 (At = ims) 14 Figure 2.5 Linear Time Varying Test Case 14 Figure 2.6 Voltage across the Time Varying Inductance (At 0.5 ms) 15 Figure 2.7 Single-Phase Two-Winding Transformer 18 Figure 2.8 Secondary Winding Voltage (At 20 = 50 ms) Figure 2.9 Equivalent Circuit for the Exponential Load 21 Figure 2.10 Two Node Test Case 21 Figure 2.11 Simulation Results of the Node Voltage (a) SFA Solution (At = 10 ms) (b) EMTP Solution (At = 0.5 ms) 22 Figure 2.12 Simulation Results of Voltage Collapse (a) SFA Solution (At = 10 ms) (b) EMTP Solution (At = 0.5 ms) 22 Figure 2.13 SFA Simulation Results of the Voltage Collapse Process (At = 10 ms) 23 Figure 2.14 Steady State Equivalent Circuit 24 Figure 2.15 Discrete Time Equivalent Circuit 24 Figure 2.16 Motor Test Case 26 Figure 2.17 Induction Motor Terminal Voltage (At = 5 ms) 27 Figure 2.18 Real Power Absorbed by the Induction Motor (At = 5 ms) Figure 2.19 Slip of Induction Motor (At = 5 ms) 27 27 Figure 2.20 Equivalent Circuit Le for Backward Euler and Forward Euler (L=1; & 1 cycle, 3 cycles, 5 cycles) 30 Figure 2.21 Accuracy of integration rules (a) Magnitude (b) Phase vii 31 List of Figures Figure 3.1 Salient-Pole Synchronous Machine and Its Windings 34 Figure 3.2 SFA Equivalent Circuit of a Synchronous Generator 43 Figure 3.3 Simulation Results with the SFA Model for Field Voltage and Mechanical Torque Changes(At=7ms) 51 Figure 3.4 Simulation Results with the SFA Model for a Three Phase Fault (At 0.5 ms): A Stable Case (fault removed at 1 .4s) 52 Figure 3.5 Simulation Results with the SFA Model for a Three Phase Fault (At = 0.5 ms): An Unstable Case (fault removed at 1 .5s) Figure 3.6 Simulation Results with the SFA Model for a Three Phase Fault (At 53 = 0.5 ms): Zoomed-in View of a Portion of Results in Figure 3. 5 54 Figure 3.7 Simulation Results with the SFA Model for a Single-Phase-to-Ground Fault (At = 0.5 ms) 54 Figure 3.8 Simulation Results by SFA: Phase A Stator Current 55 Figure 3.9 Simulation Results by SFA: Zoomed-in View of Phase A Stator Current 55 Figure 3.10 Simulation Results by SFA: Field Current 55 Figure 3.11 Simulation Results by SFA: Electromagnetic Torque 56 Figure 3.12 Simulation Results by SFA: Rotor Speed 56 Figure 3.13 Time Domain Results Using the VBR Model (At = 10 ms) 57 Figure 4.1 Torque-Speed Characteristics during Free Acceleration of a 2250-hp Induction Machine (At = 500 us) 66 Figure 4.2 Dynamic Performance of a 3-hp Induction Machine during Free Acceleration (At =500ts) 67 Figure 4.3 Dynamic Performance of a 3-hp Induction Machine during Step Changes in Load Torque (At = 500 jis) 68 Figure 4.4 Simulation Results for a 3-Phase Fault at the Terminals of a 2250-hp Induction Machine (At = 100 us) 69 Figure 4.5 Simulation Results for a 3-Phase Fault with Different Time Steps 70 Figure 4.6 SFA Equivalent Circuit of an Induction Machine 73 Figure 4.7 Stator Current during Free Acceleration of a 2250-hp Induction Machine (At 500 us) = 74 viii List of Figures Figure 4.8 Dynamic Performance of a 3-hp Induction Machine during Free Acceleration (At =500.is) 75 Figure 4.9 Dynamic Performance of a 3-hp Induction Machine during Step Changes in Load Torque (RE At = 0.5 ms, SFA At =5 ms) 76 Figure 4.10 Simulation Results for a 3-Phase Fault at the Terminals of a 2250-hp Induction Machine (At = 500 ps) 77 Figure 4.11 Simulation Results by SFA 78 Figure 4.12 Time Domain Results from an EMTP Algorithm Implemented in MATLAB 79 Figure 5.1 Schematic Structure for the EMTP Dynamic Phasor Simulation Tool 83 Figure 5.2 The H Transmission Line Model 84 Figure 5.3 Contributions of the H-Circuit Transmission Line Model to the Nodal Admittance Matrix 85 Figure 5.4 The MATLAB Code for Inserting H-Circuit Transmission Line Model into G Matrix 86 Figure 5.5 Flow Chart for Dynamic-Phasor-Based EMTP Simulator 87 Figure 5.6 One-line Diagram of a Radial Test System 87 Figure 5.7 Phase A Voltage at the Load Node (At = 5 ms) 88 Figure 5.8 Zoomed-in View of Phase A Voltage at the Load Node (At = 5 ms) 88 Figure 5.9 Phase B Voltage at the Load Node (At 88 = 5 ms) Figure 5.10 Zoomed-in View of Phase B Voltage at the Load Node (At = 5 ms) 89 Figure 5.11 Phase C Voltage at the Load Node (At = 5 ms) 89 Figure 5.12 Zoomed-in View of Phase C Voltage at the Load Node (At = 5 ms) 89 Figure 5.13 One Line Diagram for the Test Feeder 90 Figure 5.14 Phase A Voltage at the Induction Machine Node (At = 1 ms) 91 Figure 5.15 Zoomed-in View of Phase A Voltage at the Induction Machine Node (At = 1 ms) ..91 Figure 5.16 Phase B Voltage at Load Node 6 (At = 1 ms) 91 Figure 5.17 Zoomed-in View of Phase B Voltage at Load Node 6 (At = 1 ms) 92 Figure 5.18 Phase C Voltage at Node 1 (At = 1 ms) 92 Figure 5.19 The First Benchmark Network for Subsynchronous Resonance Studies 93 Figure 5.20 Generator Terminal Voltage: Phase A (At = 1 ms) 94 Figure 5.21 Transformer High Side (Bus A) Voltage: Phase A (At = 1 ms) 95 ix List of Figures Figure 5.22 Voltage across Series Capacitor: Phase A (At = 1 ms) 95 Figure 5.23 Infinite Bus (Bus B) Voltage: Phase A (At = 1 ms) 96 Figure 5.24 The Second Benchmark Network for Subsynchronous Resonance Studies 98 Figure 5.25 Generator Terminal Voltages: Phase A (At = 1 ms) 99 Figure 5.26 Voltage across Series Capacitor: Phase A (At = 1 ms) 99 Figure 5.27 Transformer High Side (Bus 2) Voltage: Phase A (At = 1 ms) 100 Figure 5.28 Infinite Bus (Bus 1) Voltage: Phase A (At = 1 ms) 100 Figure B. 1 Free Acceleration Characteristics in Stationary Reference Frame 119 Figure B.2 Free Acceleration Characteristics in Rotor Reference Frame 120 Figure B.3 Free Acceleration Characteristics in Synchronous Reference Frame 121 Figure B.4 Dynamic Performance of a 3-hps Induction Machine during Step Changes in Load Torque 122 x Acknowledgements I would like to express my deep gratitude to my supervisors Dr. José R. MartI and Dr. Hermann W. Dommel. My admiration for their great achievements and contributions was the reason why I decided to study with them in Canada. I want to thank Dr. MartI for his consistent guidance and financial support during my study at UBC. The inspiration and warm personality of Dr. MartI and Dr. Dommel have won the author’s highest respect and love. I sincerely thank my friends and colleagues at UBC and British Columbia Transmission Corporation for their encouragement and understanding. My wife Helen, my parents Yihua and Guangfu, are sources of unconditional love and support since the beginning of this project. They are entitled to more thanks than can be expressed here. My parents-in-law Guier and Jingxiang, my brother Kun, my brother-in-law Hailu, my sister-in-law Haiyang are also important members of the team that supports me permanently. To all my family, from the depths of my heart, thank you. xi Chapter 1 Introduction Chapter 1 Introduction 1.1 Background Although a power system blackout is a low-probability event, whenever it happens, its impact on the power system and the society is catastrophic. A typical example is the Aug. 14, 2003 Blackout in the USA and Canada, which caused a total loss of 62,000 MW of loads and tripping of 531 generators [1]. Fifty million people were affected by this system collapse. During the next two years, twelve major blackouts happened in Europe, Asia and North America. A recent blackout happened in Los Angeles on September 12, 2005 affecting millions in California. In today’s deregulated electricity market environment, power system planning and operations are largely driven by economic motivations, without sufficient investments in new generation and transmission equipments. This might have contributed to the higher frequency of blackouts in recent years [2]. Because of aging infrastructures operating under stressed conditions, power system stability, including transient stability, voltage stability, and inter-area oscillations have become a major concern in the North American power grid. Traditionally in power system studies a number of simplifying assumptions are made to analyze different types of stability problems with specific tools [3]-[5]. For instance, the static techniques for long-term voltage stability analysis [6] only solve the steady-state (algebraic) response of the system. Another example is quasi-steady-state (QSS) analysis [7], which neglects the fast network transients and only considers electromechanical dynamics. Because the time varying electrical quantities in QSS are represented with phasors, these tools can only capture snapshots of the system operation ignoring the dynamics between states. Problems such as voltage stability depend strongly not only on the machine-network interaction, but also on the interaction between the electrical network and the loads. Moreover, 1 Chapter 1 Introduction severe disturbances in the system often involve frequency excursions in addition to voltage excursions [8]. Especially when the power system is approaching critical collapse points, the QSS or steady-state assumptions can result in inaccurate predictions. Tools such as the EMTP [9][10] can most closely simulate the real power system dynamics by continuously tracing the evolution of the system state in arbitrary multi-phase networks with lumped or distributed parameters. Therefore EMTP-type simulators [11] -[15] are very appealing to be used for power system stability assessment. These tools, however, require small discretization steps dictated by the need to trace the instantaneous values of all waveforms. This makes the EMTP unnecessarily slow to trace phenomena around the 60 Hz fundamental frequency. As a matter of fact, nowadays one can only simulate small-scale equivalent power systems on the real-time EMTP simulators [10] [16]. Therefore, efforts should be taken to bridge the gap between the EMTP-type simulators and a unified power system analysis tool. Reference [18] aims at combining electromagnetic and electromechanical simulations [19] by considering both the network and machine differential equations without making QSS assumptions. However, no efficient solution algorithms are provided in this reference and the presented machine models are still reduced order models in phasor form. Reference [20] provides the first systematic method for a unified steady-state and transient power system analysis tool, by combining dynamic phasor concepts with EMTP simulation. Dynamic phasor models for transmission lines and linear elements are presented in this reference. However, the machine model is still based on the dq0 model, which can present numerical stability problems for large time steps and may not be the most appropriate for a general EMTP-based dynamic phasor unified theory. Important contributions to dynamic phasor modeling of electric machines have been presented in [21], [22], [23] based on the approach of the generalized averaging method [23] which, however, is also an approximate method with higher frequency components truncated. In general, it is a challenging topic to correctly, efficiently obtain the time-domain simulation results in the neighborhood of the fundamental frequency without making QSS or other simplif4ng assumptions. Hence, the intention of this thesis is to explore an effective way to combine the EMTP solution and the dynamic phasor concept, to build dynamic phasor models for power system components, and to develop a general-purpose simulation tool to obtain dynamic phasor results from EMTP solutions. 2 Chapter 1 Introduction 1.2 Motivation It is important to have some insight into the power system signals before we can further develop an efficient and accurate method for power system dynamic simulation. A conmion observation is that, when a contingency happens, the frequencies in a power system are usually close to fundamental frequency co 3 (60 Hz). That is to say, in the dynamic situation, the electrical signals in power systems have a fundamental frequency component modulated by slower events. This fact is analogous to the situation in communication systems where the carrier frequency is very high and carries on its sidebands the much lower frequency (e.g., audio or video) information. In communications theory, the carrier frequency is modulated to send out the information and then demodulated to recover the information at the receiving end [24]. Typically the carrier in a communication system is a sinusoidally varying signal at some frequency which is much higher than the frequencies contained in the information or message. The process of modulation gives the signal a nonzero bandwidth that is usually much smaller than the carrier frequency. Thus the signal can be regarded as a narrowband process, and can be accurately modeled as the product of a bandlimited ‘message’ waveform and a sinusoidal carrier. By this analogy, one can directly draw the conclusion that usually a power system signal u(t) can be modeled by a narrowband signal which has a small bandwidth around a. In steady-state AC circuit theory, one uses complex exponentials (i.e. phasors) to represent real sinusoidal signals with the understanding that the real part of the complex exponential gives the physical time-domain quantity of the AC signal. Going beyond the steady-state analysis, the analytic signal of a real function plays the similar role for a dynamic waveform. In fact, the analytic signal z(t) of a bandpass signal u(t) can be derived as: (1.1) where UQ) is called the dynamic phasor of u(t), zQ) = uQ) + 1H[u(t)], and H[.] denotes the Hilbert transform [25][26]. The dynamic phasor UQ) is the complex envelope of u(t), and will degenerate to the phasor U if uQ) becomes a sinusoidal signal in steady state. In fact, the analytic signal of a t+ço) is Be°’ where B 3 steady-state sinusoidal function Acos(a 3 = Here the dynamic phasor Chapter 1 becomes Introduction which is the phasor ofAcos(cot+q). Similar to the phasor, the dynamic phasors can be transformed back to time domain waveforms by taking the real part of the analytic signal in (1.1) (proof is shown in Chapter 2) It can be seen in (1.1) that the spectrum of U(t) is the spectrum of the analytic signal of u(t) shifted by synchronous angular frequency o. This can be called the “frequency shifting” — property. If the power network is formulated by using dynamic phasors instead of the instantaneous time quantities, it is said to be modeled in the shifted frequency domain, as opposed to the model in the time domain. Similar idea is proposed in reference [27]. The way to build power system component models and solve the network equations in the shifted frequency domain is called the Shifted Frequency Analysis (SFA). A major advantage of the Shifted Frequency Analysis is that it allows the use of large time steps in the EMTP solution environment to accurately simulate dynamic frequencies within a band centered around a fundamental frequency. The original system is transformed into a shifted frequency system where the frequencies around the power frequency (e.g. 60Hz) become frequencies around dc (0 Hz) [28][29]. Because the time step in EMTP simulations is limited by the Nyquist frequency, low frequencies in the shifted frequency system imply the feasibility of using large integration time steps. By coupling SFA models into an EMTP-type simulator, the shifted system is then numerically integrated to obtain dynamic phasor solutions, which are more easily understood by power system operators and planners than instantaneous time domain results. At the same time, the dynamic phasor results can then be transformed back to time domain waveforms using the inverse transformation. EMTP uses the trapezoidal integration rule that is A-stable, reasonably accurate and simple. The discrete-time nodal equations formed in EMTP are elegant and can preserve the sparsity of the power system structure. All these features makes EMTP a standard simulation tool in the power industry that can model the power system at the device level. The goal of this thesis is to implement the SFA method into EMTP, which is the first practical step for building a unified power system analysis tool based on the EMTP solution. 4 Chapter 1 Introduction 1.3 Contributions A two-step strategy is employed in the research of Shifted Frequency Analysis. First, the SFA-based models for power network components are developed and tested. Second, a generalpurpose SFA-based EMTP tool is built to simulate a power network. The main contributions are presented in Chapter 2 to Chapter 5 and are briefly summarized as follows. I. The definition of dyhamic phasor and the theoretical basis of the SFA method are introduced. The basic procedures for analyzing time functions (signals) in the SFA domain are discussed, and the numerical accuracy of the discrete-time SFA solution is analyzed. Chapter 2 summarizes these theoretical works. II. A few types of power system components are modeled in the SFA domain. Chapter 2 describes the SFA-based models for the linear circuit components, transformer and some load models such as exponential load and steady state induction motor load. III. A new efficient synchronous machine model for the simulation of slow system dynamics using EMTP modelling in the Shifted Frequency Analysis framework is presented in Chapter 3. Discrete-time SFA synchronous machine models are derived which have a similar form as the EMTP component models. IV. In Chapter 4, the SFA method is extended to model the induction machines in EMTP. By analyzing the relationship between rotor and stator physical variables, a phase-coordinate model with lower number of equations is first derived. Based on this, a SFA model is proposed as a general purpose model capable of simulating both fast transients and slow dynamics. Case study results have confirmed that the SFA induction machine model is a valuable component for real time EMTP simulations. It is observed that the SFA model is in excess of 70 times faster than the EMTP phase-coordinate model when simulating dynamics with frequency spectra close to the fundamental power frequency. V. A final contribution in this thesis is a dynamic-phasor-based EMTP simulation tool. The structure of the tool and the test system results are reported in Chapter 5. 5 Chapter 2 Shifted Frequency Analysis Chapter 2 Shifted Frequency Analysis 2.1 Analytic Signal and Hubert Transform The Fourier transform F(co) of a real-valued functionf(t) is a Hermitian function, which means F(—co) = F*(w) in the frequency domain. Thus the F for negative frequencies can be expressed by F* (complex conjugate of F) for positive ones. This means that the positive frequency spectra is adequate to represent a real signal, and the negative frequency components of the Fourier transform can be discarded without loss of information. Now one can construct a function Fa(0) that contains only the non-negative frequency components of F(co) by defining F((D) + sgn(aF(ü) Fa(O)) (2.1) where sgn(o)= 1 a)>O 0 —1 o=0 o<0 The analytic signalfa(t) is defined to be the inverse Fourier transform Of Fa(CO) L(t)=f(o0)td = -- f do t F(o)e° (2.2) It is clearly seen thatLQ) is a complex function. If we define the Hubert transform off(t) as the imaginary part Offa(t), and denote it by H[fQ)] = f(t), thenL(t) can be expressed as fQ) =JQ) +jfQ) (2.3) 6 Chapter 2 Shifted Frequency Analysis From equation (2.1) and (2.3) we can get the following Fourier pair JQ) +jfQ) -> F(o) + sgn(F(w) F Because we haveJ(t) F(w) by definition, we can derive that -jsgn(co)F(co) f(t) Using the transform pair iF — -+ —jsgn(co) and taking the inverse Fourier transform of —jsgn(o).F(a), we can get H[fQ)]=j(t) = ±*fQ)= --PVfdr (2.4) Note that the Hubert transform is an improper integral, thus it is calculated using the Cauchy principal value (see the ‘PV’ in (2.4)), i.e. f(t) iimi[ j t—r -A dr + A dr]. Obviously, H[f(t)] is real. The analytic signal has no negative frequency components; moreover, the original real signal can be converted back from it by simply dropping the imaginary part. The analytic signal is thus a generalization of the phasor concept. This important implication leads to the shifted frequency analysis method as is explained in the following section. 2.1.1 Shifted Frequency Analysis If a power system signal u(t) is a bandpass one, it can then be represented [201 as. u(t) = u 1 (t) cos oit — UQ (t) sin ot where the lowpass signals u 1 (t) and (2.5) UQ (t) are the in-phase and quadrature components of the bandpass signal, respectively. The dynamic phasor U(t) for the signal u(t) is defined as U(t)=uJ(t)+juQ(t) which is the complex envelope of u(t). Assume that z(t) is the analytic signal ofu(t), then z(t)e°’ = {uQ) + jH[u(t)Je0t 7 Chapter 2 = Shifted Frequency Analysis 1 (t) cos ot {u — UQ (t) sin ot + 1 (t) sin cost + UQ (t) cos j[ (t)+fu 1 u ( 0 t) (2.6) where H[.] denotes the Hubert transform. Therefore, z(t)e°’ U(t) (2.7) Note that the dynamic phasor U(t) will degenerate to the phasor U if u(t) becomes a sinusoidal signal in steady state. Equation (2.7) shows that the spectrum of U(t) is the spectrum of the analytic signal of u(t) shifted by synchronous angular frequency co,. This is called “frequency shifting”. — Shifted frequency analysis (SFA) [24] allows the exact simulation of frequencies within a band centered around a fundamental frequency using large time steps in a discrete-time EMTP type of solution environment. The original system is transformed into a shifted frequency system where the frequencies around the power frequency (60Hz or 50Hz) become frequencies around dc (0 Hz). The shifted system is then solved using an EMTP solution. The equivalent circuit for network components in the shifted frequency domain can be derived in three steps: (i) Create the phase-coordinate differential equations of the component in the normal unshifted domain; (ii) Transform phase quantities into dynamic phasor variables according to Equation (2.7); (iii) Discretize the dynamic phasor equations using an integration method, e.g. the trapezoidal rule, and build the equivalent circuit suitable for an EMTP solution. For instance, the time domain voltage equations for M coupled inductances, shown in Figure 2.1(a), can be written as v(t)=tLt) (2.8) As illustrated in Figure 2.1(b), the standard EMTP equivalent circuit derived with the trapezoidal rule can be written as 1(t) = R v(t) 1 — (2.9) hL (t) where 8 Shifted Frequency Analysis Chapter 2 RL (2.10) =- is an MxM matrix of resistances, and hL(t)=L’v(t — At)At + hL (t — (2.11) At) is a vector of past histories. To obtain the SFA equivalent circuit (Figure 2.1(c)), we transform (2.8) using (2.7) to get V(t)= L dIQ) L 1(t) 3 +jo (2.12) where V(t) and 1(t) are the dynamic phasor vectors corresponding to the physical time vectors v(t) and i(t), respectively. We can now transform (2.12) to discrete time using, for example, the trapezoidal rule, 1(t) = V(t) 1 R — (2.13) HL (t) where the equivalent resistances RL and history terms HL (t) are expressed as (2 2L (2.14) RL =—+KL 4 2. —-—+JO) At 2 V(t 1 L — At) HL (t — — At) (2.15) -+jcO The corresponding SFA equivalent circuit is shown in Figure 2.1(c). L p p (a RL = RL 2L/At hL(t) (b) = (2/At+ jco)L HLQ) (c) Figure 2.1 (a) M Coupled Inductances; (b) EMTP Equivalent Circuit; (c) SFA Equivalent Circuit A major advantage of SFA modelling is that since the equivalent circuits are centered 9 Chapter 2 Shifted Frequency Analysis around 60 Hz (ja)L for the inductances of Figure 2.1), deviations from 60 Hz correspond numerically to very low frequencies (like deviations from 0 Hz in the unshifted domain) and a large integration step can be used in the solution. For example, respecting the accuracy limits imposed by the Nyquist frequency and the distortion introduced by the discretization rule [30], allowing a 3% error with the trapezoidal rule, the integration step can go, for example, from about three 60 Hz cycles, 5Oms, for frequency deviations of± 2Hz to as large as is for frequency deviations of± 0.1Hz. With these time steps, the SFA method provides an effective way to make the EMTP program a general purpose simulator for power system stability problems. By applying SFA in an EMTP-type simulator, time varying phasor solutions are obtained, which can be easily understood by power system operators and planners. At the same time, detailed waveform results can also be traced back from the SFA results by using the inverse transformation u(t) = Re[U(t)e.t 1 (2.16) Unless specifically noted, in this paper uppercase letters are used to represent dynamic phasors in the shifted frequency domain, while lowercase letters denote real variables in the time domain. 2.2 SFA-Based Network Component Models The component models in the SFA domain are building blocks for a SFA-based network simulator. In this section, we implement the SFA modeling technique described in Section 2.1.1 and build component models for the linear RLC components, transformer, exponential load and steady-state induction machine. SFA equivalent circuits for other network components can also be derived following the procedures in Section 2.1.1. Particularly, two important dynamic components in the power system, i.e. synchronous machine and induction machine, are modeled in SFA domain in Chapter 3 and Chapter 4. 10 Chapter 2 Shifted Frequency Analysis 2.2.1 Equivalent Circuits for RLC in the Shifted Frequency Domain The equivalent circuits for network components in the shifted frequency domain can be derived by first writing the component equations in the phase domain and then relating phase and dynamic phasor variables according to the frequency shifting transformation. A. M-Phase Resistances The time domain voltage equations for M-Phase resistances are expressed by vQ)=Ri(t) (2.17) By using (2.7), we can obtain the SFA form of equation (2.17) VQ)=RI(t) (2.18) Series resistance matrix sometimes appears as a part of the it-circuit representation of the transmission line. if an M-phase transmission line is modeled as an M-phase it-circuit, the series resistance matrix will be a full MxM matrix. The off-diagonal elements come about because the earth return is eliminated as the (M+1)th conductor [9]. B. M-Phase Coupled Inductances Equations (2.12)-(2.15) give us the shifted frequency domain equations and the discretetime equivalent circuit for L. C. M-Phase Coupled Capacitances The time domain voltage equations for the M-Phase capacitances are written as follows (2.19) The SFA form of equation (2.8) is obtained by employing the shifted frequency transformation, I (t) = C dV(t) +K V (t) (2.20) where K jwC. By discretizing (2.20) by the trapezoidal integration rule, the difference equations are obtained (2.21) I(t)GV(t)h(t) where 11 Chapter 2 Shifted Frequency Analysis Rc=(Gc)’ = hc(t)= Equation (2.21) is the EMTP equivalent circuit for C in the shifted frequency domain. D. M-Phase Coupled R-L Branches The time domain voltage equations for a series connection of M-Phase resistances and coupled inductances, shown in Figure 2.2(a), can be written as R 1(t) v(t) + L diQ) (2.22) The standard EMTP equivalent circuit (Figure 2.2(b)) for the series connection of R and L are reproduced as follows i(t) =rv(t) —h(t) where rL =R+--L=g is a matrix of resistances, and _I]v(t_At) is a vector of history terms. We now transform (2.22) using (2.7) to obtain the dynamic phasor equations V(t)= R I (t) + L dI(t) +KL 1(t) (2.23) where KL =jco L, V(t) and 1(t) are the dynamic phasor vectors corresponding to the time-domain 3 real variables v(t) and 1(t), respectively. Equation (2.23) can be transformed to discrete time by using the trapezoidal rule, 1(t) R;V(t) — (2.24) H (t) where the equivalent resistances R and past histories H (t) are expressed as R= 12 Chapter 2 Shifted Frequency Analysis _I1V(t_At)_G(R_—-L +KLJ H(t-At) The SFA equivalent circuit described in (2.24) is illustrated in Figure 2.2(c). R L (a) R=R+ [(2/At) +jcoJL r =R+(2/At)L H(t) (b) (c) Figure 2.2 (a) Series Connection of M-phase R and L; (b) EMTP Equivalent Circuit; (c) SFA Equivalent Circuit E. Case Study Two test cases are simulated to explore the feasibility of the SFA modeling method. The first test case is to simulate the switching operations in a typical linear time invariant (LTI) circuit that is shown in Figure 2.3. A voltage source is switched into the circuit at t = 0 s. After 8 ms, the switch Si opens once the current flowing through the switch crosses zero. In this circuit, 1 R = 3 2, R =50 2, L 2 1 = = 1000 mH, C 2 300 mH, L 1 = 20 jiF, C 2 = 61 iF. The voltage source has a rms value of 230 kV and a frequency of 60 Hz. 1 Ri Li 2 3 R2 L2 Figure 2.3 Linear Test Case The current following through the branch connecting node 1 and node 2 is illustrated in Figure 2.4. From the simulation result we can find 13 Chapter 2 Shifted Frequency Analysis (1) During the transient state, the SFA solutions are the envelop of the time domain solutions; (2) During the steady state, the SFA solutions degenerate to the phasor, which is the concept being widely used in power industry; (3) Time domain simulation results can be accurately traced back from the SFA solutions. We can see from Figure 2.4 that the time domain curve is identical to that obtained from the EMTP time-domain simulation. Shifted Frequency Analysis: Test Case -I- 9 Time (s) Figure 2.4 Current Flowing through Branch 1-2 (At = ims) The second test case is a linear time varying (LTV) circuit. A sinusoidal current source 1(t) = A cos(c#) is applied to a time varying inductance L cos(o t), as is shown in Figure 0 2.5. The frequency of the current source is co = 60Hz, and co =10Hz. The current source has a magnitude of 1 A, and we choose an integration step At = 500ps. The unknown variable v(t) is the voltage across the inductance. 1(t) = v(t) Acos(o0 Figure 2.5 Linear Time Varying Test Case The numerical solution to the dynamic phasor V(t) can be readily obtained by using the trapezoidal rule, which is given by 14 Chapter 2 V(t) Shifted Frequency Analysis = —V(t — At) + [f0JL cos(co L sin(a 0 t) + 0 t) + o) 0 [JoL cos(o. 0 (t — At)) + L sin(a 0 (t — At)) — .. cos(aot)]I(t) + 0 (t cos(a — At))]I(t — At) The time domain solution of voltage across the inductance can be transformed back from V(t) by v(t) = Re[V(t)e3]. time (s) Figure 2.6 Voltage across the Time Varying Inductance (At = 0.5 ms) The simulation result is shown in Figure 2.6. It is clearly seen that the dynamic phasor solution is the envelope of the time domain solution. This test case validates that the SFA method is suitable for the LTV circuit analysis. 2.2.2 Options between Complex Arithmetic and Real Arithmetic As can be seen in Section 2.2.1, the discretized power systems equations based on SFA models will be a system of linear complex equations, which can be solved by using either the complex arithmetic or the real arithmetic. Some references [311 [321 concluded that, for the complex matrices, the complex inversion of them might be up to twice as fast as the real inversion, and that the rounding error bound for complex inversion is tighter than that for real inversion in terms of Gauss elimination. Fortunately, this may no longer be a notable issue with the modem computer architectures. Therefore, the real arithmetic method can be used as an alternative to the complex-arithmetic-based system solver. This section gives some examples about building real valued equivalent circuit for the RLC elements. 15 Chapter 2 Shifted Frequency Analysis A. M-Phase Resistances The discretized equations (2.18) can be decomposed into two equation sets corresponding to the real and imaginary parts of the dynamic phasors for voltages and currents. This leads to the real valued equivalent circuit of R. lEVre(t)1 fIre(t)1 = rR’ [‘im (t)J L (2.25) R’ ]LYm (t)] B. M-Phase Coupled Inductances The decomposed difference equations for the M-Phase coupled inductances can be derived as Ire (t)1 = raL 1 [‘im (t)] lrvre (t)1 bL’ 1 ]L”i aL m (t)J — 1 LbL (t)1 rm ’L,I 1 L hLre — (2.26) (t)j where P’ (t)1 [hLIrn (t)] — — rcL-’ — LdL1 ‘-‘ (t 1 ][V. (t cL — — At)1 + At)] rLf e fThL,re (t — — At) e ][hLjm (t — At) 2 a=y = bCzD2 4E(2 2 Atk At) -—-il—-i-a = Lzxt 2 r(2 21 E(2 LAt) SJ LAt) (22 -I—-I+a) — I d II—I+(L) e— (42 21 VAt) (22 I—I+w VAt) I Il—I 4 2 S 2 +0) S 2 0)— (2 —1+0) VAt) 2 C. M-Phase Coupled Capacitances We can get the real valued equivalent circuit for C as the following 16 Chapter 2 Shifted Frequency Analysis rlre(t)1 = C O)C c --c At LIirn(t)J rvre(t)1_rhc,re(t)1 (2.27) LVirn(t)J Lhcim(t)J where 4 r hCre(t)1 17 0’)] Lh Ac — 0 — At rvre(t m (t c L’ — — At)1 At)] Ehc,re(t At) [hcim (t At) — — — D. M-Phase Coupled R-L Branch By decomposing the equivalent circuit in (2.24), we can obtain the real valued equivalent for a series RL branch rlre(t)1 rG,re _Gpj,jmjVre(t)][hp,re(t) [Gjm [h 0’) , Gre iLVim 0’)] 1 — , (t)] 1 LI — 228 where Eh,re (t)1 0’)] — — rc Dlrvre 0’ — LD c ]LV,rn 0’ — — At) At)] — — [F FlEhc reQ E ][hcim (i’ — — At) At) A=R---L At L 3 B=o = G D = —G C — RL,re jL,im E + G RL + G re AG jLre BG — Gpj BG + Gpj AG — G jL,reBGi?J,,im + G RL,re AG jm — — G jm G jm BGpjm G RL,re A GpjjB — F 2.2.3 Transformer Model in the Shifted Frequency Domain Since the SFA modeling focuses on the low frequency dynamics as opposed tO the standard EMTP modeling for simulating the fast transients, here only the low frequency transformer model is built in the SFA domain. For transients with highest frequency less than 2-3 17 Chapter 2 Shifted Frequency Analysis kllz, the transformer can be modeled as a series connection of multi-phase coupled R and L branches, in which each RL branch can represent one transformer winding [9]. Therefore, the SFA domain model of a transformer follows directly the model described in equation (2.24). Here a major issue is how to obtain the G 1 (thus the history terms) in (2.24). A traditional way for single-phase transformers is the [Z] matrix method, which first builds the coupled impedances from the transformer parameters. For example, the [Z] matrix of a two winding transformer illustrated in Figure 2.7 can be obtained as below [z]= Then, [Zj is inverted to get the admittance matrix [Y] that will be used to calculate G . This 1 method, however, has limitations such as: (i) Magnetizing impedance Zm cannot be neglected or set to be co in this model. Otherwise the model cannot work. (ii) Zm usually dominates [Z] matrix because Zm >> leakage impedances Zi and Z . However, it is the leakage impedances that largely determine the simulation results. This 2 means all data inputs are required to be very precise. Moreover, [Zf 1 is ill-conditioned because of the dominating Zm, which may negatively affect the accuracy in the simulation results. [Zj 1 Z VI :N 1 N 2 Zm j 2 Z I’ I Figure 2.7 Single-Phase Two-Winding Transformer 1 model [9] adopted in MicroTran, rather than the [Zj matrix In this section the [Lf model, is used. This model can directly formulate the [Lf’ matrix without doing inversion on the [Z] matrix. It works for any number of windings, and for single-phase as well as three-phase transformers. For instance, the [Lf 1 matrix of a two winding transformer in Figure 2.7 can be built by 18 Chapter 2 Shifted Frequency Analysis 1 N 1 = 1 [L] Ni (2.29) (N i 2 ) L 2 N L 2 N where 2 L=Li+JL An alternative to the [Lf’ model is to use an inductance matrix in series with an ideal transformer [14]. The theory behind the [Lf’ model and the [Lf’ matrices for other types of transformers can be found in [91. The equivalent conductance and history terms matrices can be calculated by using (2.24). Note that some modifications in (2.24) are needed for the usage of [Lf’, as expressed below G =R+--L + KU HRL(t)={[L’R + =[LR+(_+ J)I] L-’ 1 frD JI] [LR + + + [L-’R + +. [L-’R + + JDsJI1G — G}V(t — At)— jo i] H(t—At) The magnetizing branch is not required in this model. It can be added at the terminals when it is needed. This will allow one to model the saturation effect of the core by adding the nonlinear Lm at the terminals. As a test case, a two winding transformer with the on-load tap changer (OLTC) is simulated using the SFA model. The rated line-to-line voltage of the primary winding is 110kV. The winding ratio is N : N 1 2 = 110 : 28.4. R = 10 , X 33.2174 2. The transformer is first connected to the rated voltages and operating until the primary winding voltage suddenly drops to 0.95 p.u. at t= lOs. Five seconds later, the secondary side OLTC changes the tap ratios to 1+ 0.0125, attempting to bring the voltage back. The total simulation time is 20s. Figure 2.8 illustrates the dynamic phasor result together with the time domain result transformed back from SFA solution. In this simulation, a very large time step At = 50 ms is used. On the other hand, if 19 Chapter 2 Shifted Frequency Analysis we simulate this case using the EMTP algorithm, the time step for an accurate solution should at least satisfr At < 11(5 x 60) =3.3 ms to respect the Nyquist frequency limit. Transformer Test Case 02 — SFA Dynamic Phasor Solution zsrV_ 0 2 4 6 8 10 Time (s) 12 14. - 16 V 18 20 Figure 2.8 Secondary Winding Voltage (At = 50 ms) 2.2.4 Load Models A. Exponential Load Model An exponential load means the power consumed by the load depends exponentially on the load voltage, which can be expressed as \flp ( PQ) = (2.30) PoLJ 0 \,nq ( QQ) (2.31) = 0 where Po, Qo, and Vo are the pre-disturbance conditions of the load. For special cases such that load parameters np or nq becomes 0, 1, or 2, the load model will represent the constant power, constant current, or constant impedance load, respectively. The following equations express the discrete-time equivalent circuit (see Figure 2.9) for the exponential load. They can be directly derived from (2.30) and (2.31) by using the backward Euler rule. The backward Euler rule avoids predictions in the simulation with load models, though it may introduce some damping into the equivalent circuit. The trapezoidal rule can also 20 Chapter 2 Shifted Frequency Analysis be used but iterations have to be applied at each time step due to the nonlinearity of the exponential load model. 2 IV(t At)1 - RlOOd (t) = (2.32) P(t) V(t At)1 2 - kOd(t) = d (t) 1 RL hLld(t) (2.33) oQ(t) + = —4 At (2 I + (2.34) 10 (t) 0 L JO) = 1load (t)V(t L — Joi At P(t) = At) + 2 At 2 At — hL,d (t — At) Iv(t—At) (VQ—At) (2.35) (2.36) jiq (2.37) where I V(t)I is the magnitude of dynamic phasor V(t). vc —b Roa’ Lload —, 0 RlaRLl hLload P+jQ Figure 2.9 Equivalent Circuit for the Exponential Load In this section, a two node test case shown in Figure 2.10 is simulated to demonstrate the load modeling by the SFA method. This test case for voltage stability studies is taken from [33]. A B ZjO.25 P+jQ Figure 2.10 Two Node Test Case 21 Chapter 2 Shifted Frequency Analysis First, the system is operating with the load P= O.8p.u., and Q=O.4p.u. The load voltage at node B is illustrated in Figure 2.11. It can be seen that the SFA curve is the complex envelope of the time domain curve, to which the SFA result can be accurately transformed. Two Node Test Csse Two Node Test Case 92 90 > 88 86 84 -—0 1 2 3 4 - 0 0.2 t (a) 0.4 0.6 0.8 t (s) (a) (b) Figure 2.11 Simulation Results of the Node Voltage (a) SFA Solution (At = 10 ms) (b) EMTP Solution (At = 0.5 ms) Second, the system is operating with a heavier load of P= 1 .85p.u., Q=O.1 5p.u. From the nodal voltage shown in Figure 2.12, it can be seen that the system collapses in this case. These results in Figure 2.12 are almost the same as the results obtained in [33], where the system is found to collapse when P=1.843p.u and Q=O.15p.u. Two Node Test Case Two Node Test Case t (s) Figure 2.12 Simulation Results of Voltage Collapse (a) SFA Solution (At ms) 10 ms) (b) EMTP Solution (At = 0.5 Now the system performance is simulated when the system load keeps increasing with time. The system load is assumed to change in the way shown in Table 2.1. The SFA results of 22 ______ Chapter 2 Shifted Frequency Analysis the load node voltage are illustrated in Figure 2.13. It can be seen that the SFA model can accurately trace the voltage stability behavior of the system when the exponential dynamics of the load are taken into consideration. Table 2.1 Load Increased with the Time Time interval (s) ActiveLoad(p.u.) Reactive Load (p.u.) [0, 2] 0.7 0.3 [3, 4) 0.8 0.4 [4, 5] 0.8. 0.5 [5, 6] 0.9 0.5 [6, 7] 1.0 0.6 [7, 8] 1.85 0.15 Two Node Test Case 12C 110 100 90 1 80 70 , 60 50 40 t (s) Figure 2.13 SFA Simulation Results of the Voltage Collapse Process (At = 10 ms) B. Steady State Induction Machine Load The steady state induction motor load is modeled in this section. The machine may be first initialized at the beginning of the simulation. Usually, the induction machine parameters such as R, X, XM, XR, RR are known quantities as input data. In addition, the real power P and the terminal voltage are also known, where V can be determined by the steady state power flow calculations. Hence the unknowns for initialization are the slip s, and reactive power Q. With the core loss neglected, the equivalent circuit of the induction motor is shown in Figure 2.14. 23 Chapter 2 Shifted Frequency Analysis X -- P + jQ, XM :RIs Figure 2.14 Steady State Equivalent Circuit The formula to calculate the initial slip s and the reactive power Q can be derived [34] based on the equivalent circuit in Figure 2.14, as expressed below. 2A s= _g2 -4AC (2.38) RR Y) 2 Q=Im(V (2.39) where A=A’---R 2 V A’ + (x + = B=B’---X XM )2 B’ = 2RX C=C’---RS(XS+XM) 2 V C’ = (XSXM + XSXR )2 + XRXM + R:(xR + XM) 2 Then the initial slip can be validated by checking whether the power factor falls into a reasonable range, for instance, 0.7<cosq=cos 2 arctan <0.9 P-I If the inequality is not satisfied, an error message will be printed. The discrete time equivalent circuit in the shifted frequency domain is shown in Figure 2.15. The phasor dynamics are introduced to the steady-state equivalent circuit to capture the dynamics in the electrical part. The parameters in the SFA model can be obtained as follows. 3 Figure 2.15 Discrete Time Equivalent Circuit 23 G = R + 1 2L (2.40) + KLS 24 Chapter 2 Shifted Frequency Analysis 1 G 2L (2.41) + 1 30 = G (2.42) R sQ hjj_ (t) = — At) At —[v (t hm (t) = J/ (t — At) v (t — At)] — 2 At + KLm) —R 2 At At) — — 4L At + — At At) st At) — +---+K At hstator (t — At) (2.43) +_!+KLS + LrJ (2.44) At) — Rr s(t At) At Rr 2L At sQ At) — 2 Rr R 2L At 2L +—-—KLS KL At hm (t 2Lm —+KLfl - h, (t) = —J (t + S 1 (t h — At) (2.45) — When the system voltages change, there is a mismatch between the electromagnetic torque and mechanical torque since mechanical torque usually cannot change instantaneously. Then the acceleration equation of the rotor mass is represented as ITm dt (2.46) pdt where pf is the number of poles. Therefore the slip can be calculated using the backward Euler rule, sQ) = sQ — At) — Pj 00 2 J II (t) — (2.47) Tm (t)]At where V (t) I(t)3 (2.48) 2 Er Rr [Rm + s(t — At) + (XTH + Xr )2] If the mechanical torque is assumed to be a quadratic function of co, 25 Chapter 2 Shifted Frequency Analysis (2.49) (t—eXt) 2 =a +13s(t—/xt)+ys where 2 at=a+bas+ccos 2 y=ccos X v -v — + x — TH (X + XM )2 Rm - — R + (X + X, )2 RX + x sX (X + X) 1 +(X+X) R 2 Equation (2.49) represents constant torque if a = 0 and b = 0, and a torque linearly dependent on o. if c = 0. In this section, a test system is simulated to illustrate the capability of the SFA method for modeling induction machines in steady state. As shown in Figure 2.16, an induction motor is supplied by an infinite voltage source through a double-circuit line. One of the parallel lines is tripped off at t = 2 s. The total simulation time is 10 s. The time step for integration is 5 ms. The parameters are R 3 = 0.0092 pu, X. = 0.0717 pu, X, .0717 Pu, 0 = 6.9 kV, pj = 2, a = b = 0, c = 0.00308, V inertia constant H = 0.6s, line-to-line voltage 0 1.0064 p.u., P 4.1375 pu, Rr = .00698 pu, Xr 0.4 p.U., Rime = 0.161 2, Lime = 0.0042707 H, C = 0.209 mF. The per unit values are based on Vbe = 6.9 kV and Sbase = 10 MVA. Figure 2.16 Motor Test Case The following are the SFA solutions of the test case, including the load bus voltage, induction motor slip and real power absorbed by the motor. These results are illustrated in Figure 2.17, Figure 2.18 and Figure 2.19. As shown in Figure 2.17, the voltage has dropped from steady state to collapsing state with the motor stalled in about 6 seconds. Artificial numerical oscillations appear in the simulation results because of imperfect initialization of the test system. In order to capture detailed dynamics between the steady state and the collapse, one needs to apply the SFA method to a detailed motor model considering the rotor and stator transients. This 26 Chapter 2 Shifted Frequency Analysis is further investigated in Chapter 4, where a general-purpose induction machine model in the shifted frequency domain is proposed. 8000 7000 600( 500( 400( — 3001 2000 I000 1 3 2 4 5 Timo (s) 6 7 8 9 10 Figure 2.17 Induction Motor Terminal Voltage (At =5 ms) -x 10’ 2 4 3 5 T,mo 6) 6 7 8 9 Figure 2.18 Real Power Absorbed by the Induction Motor (At 2 3 4 = 5 ms) 5 T.oo 6) Figure 2.19 Slip of Induction Motor (At = 5 ms) 2.3 Numerical Accuracy Analysis As mentioned in Chapter 1, the theory of Shifted Frequency Analysis (SFA) needs to be implemented in a circuit analysis program such as the EMTP to achieve its advantage of using large solution step and getting correct simulation results in the neighborhood of the 60 Hz frequency. In the EMTP, a numerical discretization rule (integration rule) is used to convert the 27 Chapter 2 Shifted Frequency Analysis equivalent circuit of each network component into an equivalent discrete time model consisting of an equivalent resistance and a history term. These equivalent circuits of components are used to build the nodal equations of the whole system, which can be solved by using the numerical linear algebra. The trapezoidal and backward Euler rules are the most common choices to perform the discretization from continuous time models into discrete time models. An illustrative and effective way to analyze the numerical accuracy of an integrator (integration rule) is to examine the behavior of an inductance L with v(t) as input and 1(t) as output (or a capacitance C with 1(t) as input and v(t) as output) [30]. The magnitude and phase distortion introduced by an integrator can be quantitatively analyzed by calculating the frequency domain equivalent circuit of L which is discretized by the given integration rule. This approach is adopted in this section to analyze the numerical accuracy of an integrator in the SFA domain. To analyze frequency responses of an inductance L, first take an input that has only one frequency co (close to w) and a unity magnitude, i.e. v(t) = cos cot. The dynamic phasor of v(t) in the SFA domain is V(t) = [v(t) + jH(v(t))]e’ = ej°°’ = e°’ (2.50) Now suppose the output is 1(t) = Y(Aco)e’. (2.51) where ) (iXco) is the admittance of an L in the discrete-time SFA domain. A. Frequency Domain Equivalent Circuits for L The SFA domain equivalent circuit of L is obtained with the trapezoidal and backward Euler rule. The equivalent circuit parameters can be plotted as functions of Af (=.4), as is 2r shown in Figure 2.20. Substitute (2.50), (2.51) into the difference equation it is easy to find Y(z\co) for each discretization rule. (1) Trapezoidal rule With trapezoidal rule, V(t)= (._+foJ LJI(t)_ VQ—At) + —._!+JcoLJI(t—At) 3 Substitute (2.50), (2.51) into (2.52), 28 (2.52) Chapter 2 + e30f Shifted Frequency Analysis et + foLJYe3c0t + = jw Lej&0(t + Then 1 2L — — 1 — t e° jcoL+—• t e —1 — tan S + /Xa/Xt 2 AoizV 2 joL+jAoL• where the equivalent inductance is AoiAt 2 L AoAt 2 tan Le= (2.53) (2) Backward Euler rule The difference equation is 1(t) — I(t — At) = .4 (v(t) — LIQ)) 5 jo (2.54) Then — Ye W(t&) = — josLYee’) Therefore the admittance can be derived l/}(Aco) = ——(joi At + 1— e)= jo 5 L+ 5 A L Let + Y (Ao) = 1 , Re 2 e’°” — 2 e°’ then J(Ds1e IRe = 1 / Re(Ye (Am)) (2.55) = —ii[o Im(1’(Ao))J 29 _____________ Chapter 2 Shifted Frequency Analysis -1 z 0 f(Hz) Figure 2.20 Equivalent Circuit Le for Backward Euler and Forward Euler (L=1; t =1 cycle, 3 cycles, 5 cycles) B. Accuracy of Discretization Rules The accuracy of the discrete-time integration rules can be expressed by the following ratio H(Aw) = (Aa) H(Aco) Y(Aco) (2.56) L=1 (1) Trapezoidal rule He(CO) HQo) — — Js 2 JcO+—• At S - — ej0( —1 +1 10)5 2 ej ’1 2 At +1 (2) Backward Euler H(o) H(a) — — = At +1 5 jco — e3A0S j0) A 5 t At + 1 5 jo — 30 Chapter 2 Shifted Frequency Analysis Frequency domain accuracy analysis results are shown in Figure 2.21. The frequency deviation ranges from -3 to 3 Hz. 0 A f(H.z) (a) I I I — I I 1.5 ackward Euler A = 3 cycles. Backward Euler At5 cycles 5 ‘ 0 Cs Backward Euler 0.5 .Atlcycle t Trapezoidal Rule -2 —l 0 A f(Hz) (b) Figure 2.21 Accuracy of integration rules (a) Magnitude (b) Phase From Figure 2.20 we can find that the distortion on L is very small for fclose to 60 Hz. When Af —* 0, Le—*L and there is no distortion. The error grows as txf increases. The distortion in the equivalent inductance also increases with the time step. 31 Chapter 2 Shifted Frequency Analysis Figure 2.2 1(a) shows that the trapezoidal rule is slightly less accurate than the backward Euler rule in the shifted frequency domain. On the other hand, Figure 2.2 1(b) shows that the trapezoidal rule has no phase distortion while the backward Euler rule can cause larger phase distortions especially when Afbecomes larger. It is clearly shown in equation (2.55) that the backward Euler rule adds a fictitious resistance to the circuit, which introduces numerical damping. Therefore, it may be advantageous when the trapezoidal rule has the risk of numerical oscillations when it is used as a differentiator. Here the backward Euler rule may be used over a few integration steps to damp numerical oscillations whenever some discontinuities occur, which is the major idea of CDA [30]. Unless specially noted, in this thesis the trapezoidal rule is used for the component modelling and the system solver because SFA equivalent circuits similar to those in the EMTP can facilitate the implementation of SFA concepts in the EMTP algorithm. 32 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis 3.1 Introduction This chapter describes a new synchronous machine model based on the Shifted Frequency Analysis (SFA) method, which uses dynamic phasor variables rather than instantaneous time domain variables. Discrete-time SFA component models have a form which is similar to the EMTP component models. Dynamic phasors provide envelopes of the time domain waveforms and can be accurately transformed back to instantaneous time values. When the frequency spectra of the signals are close .to the fundamental power frequency, the SFA model allows the use of large time steps without sacrificing accuracy. This makes the SFA method particularly efficient for power system dynamics. 3.2 Voltage-behind-Reactance Synchronous Machine Model The SFA method can be applied to a number of synchronous machine models, such as the phase domain synchronous machine model of [35], the voltage-behind-reactance (VBR) model of [36], and others. As reported in [37], the VBR model is more efficient and stable than other traditional machine models used in the EMTP and it was chosen for the SFA implementation in this chapter. 33 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis bs axis as Cs, bs Figure 3.1 Salient-Pole Synchronous Machine and Its Windings The cross-section of a salient pole synchronous generator is shown in Figure 3.1. The generator has a 3-phase stator winding, a rotor field winding (fd winding), M rotor damper windings in the q-axis (kql, ..., kqM), and N rotor damper windings in the d-axis (kdl, ..., kdN). The equivalent circuit in Figure 3.1 can be used for modelling different types of synchronous machines without losing generality since straightforward changes can be made without much difficulty to adapt this model for other types of synchronous machines. In this chapter, some assumptions [38] are made in developing the mathematical model of the synchronous machine: The 3-phase stator winding is assumed to be symmetrical; The capacitances of all the windings are neglected; Each of the distributed winding can be represented by a concentrated winding; 34 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis The change in the inductance of the stator windings due to rotor position is sinusoidal and does not contain higher harmonics; Hysteresis loss is negligible while the effect of eddy current is included in the damper winding models. In Figure 3.1, the motor convention is adopted for the machine modelling, which means the stator currents are positive when flowing into the machine terminals. The damper windings are shown to represent the. paths for induced rotor currents. In a salient-pole machine, the rotor is laminated and therefore the damper winding currents are largely confined to the cage windings embedded in the rotor surface. Usually the dynamic behaviour of the salient-pole machine can be predicted accurately enough by using one equivalent damper winding kq2 in the q-axis and one damper winding kd in the d-axis. On the other hand, a cylindrical-rotor machine has a solid iron rotor with a cage-type winding embedded in the rotor surface, and the damper winding currents can flow either in the cage winding or in the solid iron. In order to accurately simulate the transient process in the cylindrical-rotor machine, it is necessary to put two equivalent damper windings kql, kq2 in the q-axis and one damper winding lcd in the d-axis. The dqO synchronous machine model is most widely used in power system electromagnetic transient simulations. However, to interface the dqO model to the network in the phases a, b, c coordinates, some electrical variables have to be predicted, which may theoretically cause numerical instability. The dqO model needs prediction because the dqO electrical variables are not solved simultaneously with abc external network variables. A phasedomain model avoids predicting the electrical quantities and therefore is more robust than the dqO model. The SFA model implemented in this Chapter is based on phase coordinates and therefore does not need predictions of the electrical variables. A detailed discussion of the prediction in EMTP machine models can be found in Chapter 8 of [9]. The original phase-coordinate synchronbus machine model consists of three stator voltage equations, M + N +1 rotor voltage equations, M + N +4 flux linkage equations, and mechanical part equations for the torque and speed. To improve the numerical efficiency and stability, a new phase-coordinate model called Voltage-behind-Reactance (VBR) model has been derived in [36]. The rotor voltage equations and flux linkages equations are manipulated so that the derivatives of the rotor flux linkages are algebraically incorporated into the stator voltage equations, resulting in an efficient voltage-behind-reactance form model. The VBR model for the 35 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis synchronous machine in Figure 3.1 includes the stator voltage equations, flux linkages equations and rotor mechanical part equations, which can be stated as follows. A. Stator Voltage Equations va (t) i abcs (t) + 5 r = P[Lbcs (0r )‘abcs (t)] + (3.1) (t) where r 5 is the stator winding resistance matrix, and (3.2) L:bcs(Or)=L+(—L:)A L,+L: L”0 = • - —s2 —--s- LIs — + L” L”a ” 0 L — 2 —--s- cos(20r) A= - 15 L +L’ C0S(28 — 2?r/3) Co5(20 + 2n/3) 2n/3) C05(20r 47r/3) cos(20r) cos(20r + 42r/3) cos(20r + 2nj3) cos(20r) C05(20 = LZq + — — Ld 3 = — L’ -1 M L:q =1/Lmq +1/LlkJ Ld 1 +1/L +1/L ‘ -1 N =1/Lmd , 0 and (As an example, for a salient-pole machine, L, Lq = (1/Lmq + 1/L, 2 )‘ Ld = (1/Lfl,d + 1/Llka + , L = (1/Lmd + 1/L 1+ L,d may be expressed as q 3 1/Lw )‘. Otherwise, L,q = (1/L,, + 1/L, 1 + 1/LIkq )’ 2 1/Lw )‘ may be used in modelling a cylindrical-rotor machine.) The subtransient voltages in equation (3.1) can be expanded V:bcs(t) = 1 [K’(O(t))J_ (3.3) v 0 36 _____ _____ Synchronous Machine Modelling Based on Shifted Frequency Analysis Chapter 3 where Vq = L j1 ) 1 L M) N” Lrnq+ Vd f= + B. “ ,2 N + fd OrLm L2L ,kqi 2 L j=1 Lmar,jj lkqi E, ‘‘fd Lmqr M + — k) + mqr,qj 2 L •r 1 L Ikdj J + N rfd (3.4) ,k 2 L lkaj A. r .2 +—L mdl 1 2 L 1=1 L,k [T,. LmdL+_2] null —+---— — 1=’ kqi j=1 (Afd N Lfd M (L mq L 2 ds (3.5) r L” Flux Linkages Equations P’kqi —-sct L . = kqj — kq) ,f = 1, 2, ..., (3.6) M Ikqj r (3.7) Ikdj Pfd —‘A. fd Lfd (3.8) A.md)+Vfd where MA. (3.9) = Lmq (-- + qs) 1 J A. L lkqj NA. (3.10) Ad =Lmd(++ldS) j=1 L, Mechanical Part Equations C. r0JJr 8 P (3.11) PWr°jTm) (3.12) f(2isqs (3.13) Aqslds) where M Aqs (3.14) L’qlqs 37 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis NA (A (3.15) ) L 1 = Lifd Discretizing the above VBR model with the trapezoidal integration results in a more flexible and robust EMTP-type model than the qd0 model [37]. The qdo model is numerically efficient and is widely used in EMTP simulators. However, the interfacing of qd0 quantities with the electrical network, which is modeled in abc coordinates with the EMTP, is not direct and requires some prediction. The SFA model proposed in this chapter is based on the VBR model. Therefore, the SFA model is still in abc coordinates and the interfacing in the EMTP or in a hybrid environment like the OVNI [39] simulator does not require prediction but interpolation (in a similar way to what is discussed in reference [40] for the coupling of different size time steps). Section 3.3 details the development of the SFA-based synchronous machine model. 3.3 Synchronous Machine Modelling with SFA 3.3.1 Synchronous Machine Model Based on SFA The following equations are based on the salient-pole generator with one field windingfd, one damper winding kd in the d-axis, and one kq2 damper winding in the q-axis. These formulas can be modified with minor efforts for synchronous machines with arbitrary numbers of windings in the d- and q-axis, following the procedures in Section 2.1.1. To construct the SFA-based stator voltage equations, we first rewrite the electrical variables using equation (2.5) ‘abcs (t) Vabcs (t) = acsi (t) cos cost = VabCSI ‘ abcsQ — (t) cos oit — (3.16) (t) sin co t 5 VabcsQ (3.17) (t) sinco t 5 Substituting equation (3.16) and equation (3.17) into equation (3.1) results in VabcsI = — — (t) cos °‘! — V abcsQ (t) 5 L’abcsJ (t) cos ot r 1 abcsQ (t) sin ot P — sin cost 5 sin o1)st (,itj + L [P1abcsl (t) cos (flat ‘abcsl (t)0) t P abcsl (t) 3 5 cos 5 tj L,’ {p[A cos ot] ‘ abcsl (t) + A cos oi acs (t)o abcsO (t) sin — . — — t] abcsO (t) 5 P[45h1 oJ — t PabcsQ (t)}+VbCSJ (t) cos t 3 Asina 38 — VbCSO (t) sin ot (3.18) Synchronous Machine Modelling Based on Shifted Frequency Analysis Chapter 3 Applying the Hilbert Transform to the left hand side (LHS) and right hand side (RHS) of (3.18), respectively, we can construct the analytic signal (3.19), as follows 5 r t = VabcS (t)e ‘abcs (t)e°’s’ +Lg [Plabcsl (t) + f üI abcsl (t)Je.b0)t —4 {j(2ü + )Kalabcs (t) + J(2Or + KaPIabcs (t) + K bP’ abcs (t)} + Here ‘abcs (t) respectively; I , atcs (t) j(28. +or) Ka = Vabcs (t) j(28. +t—2ir/3) ej(280st_42r3) Kb = e e j(29 —t) j(28r —t—2i/3) j(20. —,t+2nj3) abcs (t) (3.19) t Vabcs (t)e° . ei(20_23) e )K bI is the conjugate ofI abcs (t) K a and Kb are expressed as e j(2O +t+2r/3) UJs are the dynamic phasors of stator currents and voltages, e e — e j(28, +t) e j(28 +t+2,/3) e1(280t) j(28 +t+4,r/3) e j(28,. —at—2r/3) e e j(29 —t—4ir/3) e j(28 —at) e1(28, —ot) Performing the frequency shifting as explained in equation (2.7), the SFA based formulas for the stator voltage equations can be expressed as Vabcs (t) = 5 r ‘abcs (t) +L [PIabcs (t) + Jo + J(2COr + co 5 )L e sI abcs ‘abcs (t) + J(2COr (t)j+ — 0)s 1 )L “e j(2o,—2f) abcs (t) * + L”e 2 p1 abcs (t) + Le J(28_2o)1)pI abcs (t) + V:bcs (t) (3.20) where 1 L=—-- a 2 13 V 2 a=ef . a The SFA formulas for the subtransient voltages Vbcs (t) can be derived using the same approach described above. 3.3.2 Discrete Time Model Discretizing (3.20) with the trapezoidal rule of integration and rearranging, 39 Synchronous Machine Modelling Based on Shifted Frequency Analysis Chapter 3 1 Vabcs(t) = A Iabcs(t) 1*abcsQ)+A I ai,cs(t 2 +A 3 — 1*abcsQ 4 At) +A — At) + V:bCS(t) (3.21) —VObCS(t—At) +Vabcs(tAt) where =r +I—+ fJL 1 A 2 =I_+jc A +_+jwLneJ29 JLe1(29t)_2øo) At 4 2 =r +I——+ 3 A At J“L” 0 4 =1 ——+ja A JLej(20__2t_ +‘___+joLffe12t ° At Equation (3.21) can be further decomposed into real and imaginary parts rvabcsr [Vacsi (t)1 (t)j r +A 1 EA r 2 = , 1 [A — +A , 2 , 1 A Air ] +A , [Iabcsr (t) 2 — 21 A ‘abcs,i rVn I abcs,r (t)1 31 +A rA r 4 , +A 3 —A jTIabcsr(tAt)l 4 [v:bcS,(t)j , 4 LA , 3 , +A 31 41 A —A ][Iabcs,,Q_At)] (t—At)1 EVabcsr(tAt) L V”abcs,, (t At)] VObCSZ (t —At)] r’Vlf I + abcs,r — — L (3.22) The discrete time formulas for the subtransient voltage can be similarly derived, giving (3.23) V”abcs (t) = K(t)Ib(t)+e(t) The time-varying factors K(t) and e (t) are expanded as follows r e (t) = ii e°’ ef(o_a_,r/2) 243 e° ei(O_mst_23_2) j(8 —o,t+2/3) —,t+2,/3—/2) 2 (t [KI — At) + K [ 2 1 ij 1 t—At)1 — At)] 0 + vf (t) + K 3 K v (t 4 — At) + ] [K (° (t))] cos(0 (t At) + 2n/3)l 2 Fcos(Q (t At)) cos(O,. (t At) 22r/3) 1 (t At) K ] 6 — L srn(Or (t At)) sin(Or (t At) 2r/3) sin(Or (t At) + 22r/3)] — — — — — — — — — 0 (3.24) 40 _________ _____ ________ _______ __ Chapter 3 ________________ _________ __ _______ ________ _______ Synchronous Machine Modelling Based on Shifted Frequency Analysis rK 1 _ 9 10 K K(t)=[K’(O(t))] (3.25) ][K(Or(t))] 0 0 [ 0 where Ecos(Or) COS(Or K;(or) = __[Sin(Or) Sjfl(Or 21 — — 1 2 2r/3)l 2ir/3) COS(Or + 2r/3) Sjfl(Or + 2r/3) 1 2 1 2 I i is the Park’s transformation matrix, and L”r L” mq kq2 1 2—At---’1— ...f!.L_1] 2 Llkq 2 Llkq 2 L lkq2 L”mq I J 2 L, N A1jlLnq I 2 l _CUr(t)L 2 Llkq ] J 2 Llkq 2 Llkq r LJ I 8 2 =K K 7 •K At -- L” N 2- At 1 1L 1 Llkd [ Lrnd )J 1 L L [ LL 01 K =4 3 K ml +_I [LJ r\1 4 =K K 7 •K 8 Li L:qrkq2 L2 lkq2 K “ L” r(J2Ltq 2 L,,, l\Ikq2 +Athl_L mq 2 l )r(t) Llkq2 J N J 2 Llkq Llkq2 “ At = 6 K • 7 • 8 K K fdnd I fd I kdmd I At’_“ Lj L” K7JLr md LL [ + L Ldrffl(L” d_1 L L” (Or(t) Lr 2 L lfd lkd + L”r ,,,d kd 2lkd L I 41 rnd llf Synchronous Machine Modelling Based on Shifted Frequency Analysis Chapter 3 -‘-1 rfd 8 K I I In fdJ LJ =[ 1 L” kd - LJJ 1 L L =K + 9 K 0 10 K + 6 =K !kd Equation (3.23) is then decomposed into real and imaginary parts EV:bcs,r (t)1 [v:bCS, (t)] — — rKr — 1 LK 1 K Kr 1EI1,r (t)1 + re,r (t) 1 (t) JL1abcs,i (t)J Le, (3.26) Substituting equation (3.24) into equation (3.22) produces E”abCs. (t)1 = R LVabcs.j (t)] E’a,0 eq (t)1 +ECh.r (t) (3.27) [Iab, (t)] [eh, (t) where Req reh,r EA +A r 1 r +K 2 11 +A A 21 =[ (t)1 [ehl (t)] = 4 ,EA 3 r +A 31 + A [A 41 EC,r +A 1 , 21 1 —A 1 A A r +K] 2 — 31 A + 4 A (t — 3 A — r ][Iasi (t 4 A — At)1 r + VabcS. (t At)] LVa”bCS,I (t — — At)1 At)] — rva,r [Vabcsi (t — (t — At)1 + At)] (t) Q) 1 [e Equation (3.27) is the SFA-based formula, which has a similar form to the phase coordinates EMTP synchronous machine model of [41]. The SFA equivalent circuit is shown in Figure 3.2. 42 ____________________ _____________________ Chapter 3 _ _ _______ __ Synchronous Machine Modelling Based on Shifted Frequency Analysis Req(t) Iabc,r(t) eh,r(t) z - -. s-H vvv) .. 3-H - - ±1+ :- +0—H +0-H Figure 3.2 SFA Equivalent Circuit of a Synchronous Generator In each time step, the,flux linkages can be updated as follows L,q 2 r kq2 2 (t) = — [cos(O (t)) COS(Or (t) COS(Or (t) + 22r/3)]Re[Iab (t)et] 2r/3) — 2+At-hi_-LJ 2 L, L” 2_At1”1 mq l\ 2 LIkq + Efd (t)I — - . At)) [ U - At -- — At) — — r — I —1 [rfdLd1 LdL Lvd L”mdII kd 2+AJ )J 1 L U (3.28) — 1 L” - J [ At rkd L”md LlkdJ 2r/3) Sjfl(Or (t) + 22r/3)JRe[Iabcg (t)e] r U” LL _LN -- UJ L, L 12 L cos(Or (t At) + 2nj3)]Re[Iab (t At)e°’] 22r/3) mL”” 1 dI U kd —A-r-—— Llkd Ld .;L L” Llkq2 COS(O (t rfd ( r 2+At—-I 1——--- [sin(O (t)) sin S(Or (t) + At--”l 2 f 2 Llkq kq2( JJ 2 ‘!kq [2+At Lkd(t rkq2 U”mq ‘ L”mq I LJ2l. — 2 L!kq J,; 2 LIkq , 1 2+At!1i1 • [COS(O (t 1\ 2+At r U L” N imI 2AtLtLIl rfd’ Lu L L 4fd At L/kd) 43 2 At - L,kd L :{ r 1I (t N lmdI Jj 1 L — At)1 Lkdt—Ati Synchronous Machine Modelling Based on Shifted Frequency Analysis Chapter 3 2+ At — 1’l L --- — — Lw,) At-—-- 2 1 L L • [Sjfl(Or(t 2+ At — At .- At--1’lL _.-‘l L,,) r LZd L AtmdI L, Llkd At)) sin(Or(t At) 2ir/3) sin(OrQ At) + 22r/3)]Re[Iabcs(t At)e°’] — 1— L + + At--LL — ‘rnd - — — At —a-- !::‘.r; L L (3.29) [t1vfa(t+vfd(t_At L” -At------ 2+At---I1---- 1 L 1 L 1 L °’ 3 2Lqs (t) = L-[cos(Or (t)) cos(Or (t) 2r/3) cos(O,. (t) + 22r/3)]Re[Iab (t)e — ds [sin(O (t)) 5(0r (t) 2r/3) sin(Or (t) + 22r/3)]Re[Iabcs (t)e’ (t) = L — 1+ L 1+ L 2Lkq2 (t) (3.30) Ikq2 [fd (t) fd + lkd (3.31) The differential equations of the generator’s mechanical part also have to be discretized and can be solved together with the electrical equations. The discrete time equations of the mechanical part are expressed ‘r (t) = ‘r (t — _.f T (t) = At) + (t) + Tm (t — At) — Te (t) + Te (t — At)] (3.33) (t)Iqg (t) + 2qs (t)idS (t)] (3.34) Or(t) = Or(tAt)+[COr(t)+COr(tAt)] r (t) = or (t — At) + At[C0r (t) + U) (t — At) — ] where p is the number of poles and ‘qdOs = (3.32) K (o (t))Re[Iabcs (t)et] Refer to [9] for more detailed multi-mass mechanical part models. 44 (3.35) Synchronous Machine Modelling Based on Shifted Frequency Analysis Chapter 3 3.3.3 Note on the Cylindrical-Rotor Machine Model The shifted-frequency-domain equivalent circuit of a cylindrical-rotor synchronous machine is similar to that of the salient-pole machine. As mentioned in Section 3.2, however, the cylindrical-rotor machine has different rotor windings from the salient-pole machine. Therefore, particular equations and factors in Sections 3.3.1 and 3.3.2 have to be re-written., which can be expressed as follows. ei(_0t) e (t) — 1 [K [K eJ(G 243 e°’ e j(r K6 22n/3_2n/2) —sr÷2/3) e j(9,. —t+2a/3—r/2) rAk (t—At)1 LAt — 1 1 At)jAkb2 (t — At) + K 2 A (t—At) — (t ] Ecos(Or [sin(Or (t — — + At) vf (t 4 vf (t) + K 3 K — At)] 1 22r/3)1 At)) cos(Or (t At) 2r/3) cos(Or (t At) + At)) sm(0r (t At) 2n/3) sin(Or (t At) + 2r/3)J 0 — — — — — + abcs (t — [iç (o (t))]-’ At) — (3.36) Lqrkql 1L Lqi 1 LlkqI ,) 1 K rkq L:q 2 2 L:qrkq Llkql 2 Lq 2 Lq (L rkql 2 L:q ) 2 LqlLIkq L !kq2 Ikql 2— At-!!—I!L 1 Ljjqi L,,, — At = i”,)1 Ikql rkqiLrnq At 2+ LlkqlLlkq2 2— At-12 L Llkq2 L 12 L 1 [!_ — 1’-- Lmq 2 rkq Lrkql K — + 2 rkq L:q At-_I-_ 1 L, 1 Lijq At — ) Lmq 2 rkq 2 LlkqlLIkq 2 L,qrkq + ‘mqkq1 — — Ikq2 lkql lkql lkq2 lkql Ikq2 L” 11 L r(t) lkq2 45 lkq2 — ) At rkqlLrnq 2 LIkqlLIkq 2+ At-2 L, (i_ 2 Li*q — Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis r d’] H =K 9 K 5 i-I L”mq rkql I —At Llkql LIkq2 Llkql Llkq 2 rkql —At [cos(O,. (t)) r L — I At L rkq2 L” mq 2 Liq I ] + mq —At - -. Llkql L”rnq 2— At—-12 [UIkq L 2 Lrnq 2 rkq —At [rkq1LIql rkqlLrnq rkql —At — —At 2- 2 Llkql L,, L” LIkql At mq Lrnq 2 rkq 12 L 2nj3) cos(Or (t) + 2n/3)}Re[Iab (t)e°’ ) L,q 2 r 1 At 2_At__[___lJj 2 L 1 L, L”inIN LIql[LlkqI rkq2 Lrnq 2 rkq cos(9r (t) L”mq 2 , 1 Likql L rkql rkql 1+ L” L” 1L L, 12 12 LIkqI L At 112 L 2 Lijq 2 rkq L” 2 LIkql L, -At--I -n--i Llkql —At rkq2 , 1 L mq [cos(O,. (t At)) cos(Or (t At) — 2 , Llkq 1 L 2 — — 22r/3) cos(Or (t At) + 27r/3)]Re[Iabc$ (t At)e°’] — — Aq 0’) = Lq”- [COS(r (t)) cos(O (t) 2n/3) cos(Or (t) + 22r/3)]Re[Iabcs (t)e°’ — LZq[kl(t) + Ikql (3.37) 1+ 12 2 1 (t (3.38) Ikq2 Thus, it can be seen that separated subroutines have to be coded for salient-pole and cylindrical-rotor machines when implementing them in a computer program. If both the salient pole and the cylindrical-rotor machine could use a single model with the same number of damper windings, this would be more convenient for the code maintenance and extension. 46 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis Before we solve the above problem, let us first examine the subtransient reactance matrix in the stator equation (3.2), where 0 L = L + (— L)A. If the dynamic saliency effect, which means L L, can be neglected from the stator equations, we get L = L — = 0 and L 8 becomes a constant matrix. This will effectively improve the numerical efficiency in simulations. References [42j has proposed a method to neglect the dynamic saliency effect. First, the operational impedances are calculated based on the original parameters of the machine. Second, an artificial q-axis winding is added in the machine model such that X’ = Xd”, at the same time, the frequency response curve of the new operational impedance accurately matches that of the original Xq (s) over a frequency range , U’ where f. is a user-defined frequency. By fitting the frequency response curve of the new operational impedance, a new set of machine parameters can be obtained. Therefore, we can always build a three damper windings model to represent both salientpole machine and the cylindrical-rotor machine with acceptable frequency response. Because the dynamic saliency is eliminated, the refitted model is numerically more efficient than that with the dynamic saliency. We then need to develop only one uniform subroutine for different types of synchronous machines. Thus the problem is solved. 3.4 Simulation Results Three test cases have been simulated to illustrate the efficiency and accuracy of the SFA model. The first two cases illustrate that the SFA model is a general purpose model, which can simulate both slow and fast dynamics. The third case shows that the SFA model is numerically more efficient and stable than the phase-domain EMTP model. Throughout this chapter, both the EMTP and SFA programs for simulating the single machine dynamics are implemented and tested in the same MATLAB environment. The synchronous machine parameters for the three cases are shown in Appendix A taken from [43]. A. Simulation of Slow Dynamics The first test case consists of a salient-pole hydro turbine generator connected to an 47 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis infinite bus. It was assumed that the system is initially operating in steady state with a mechanical torque Tm = 0 and a field voltage Vj1 .0 p.u. At t = 1 .Os, the field voltage is stepped down to 0.8 p.u. Later, at t = 7.0 s, the field voltage is stepped up to 1.2 p.u. Finally, at t = 15.0 s the mechanical torque is stepped up to Tm 0.85 p.u. The total simulation time is 20 s. The simulation results with the SFA model are shown in Figure 3.3. The variables plotted in this figure are ‘as (phase A stator current), ld (field current), Te (electromagnetic torque), or (rotor angle), cor (rotor speed), P (real power) and r are Q (reactive power); all values except for 0 in per unit of the machine ratings. In this simulation, a large time step & = 7 ms is used. Note that the SFA model can produce both dynamic phasor and time domain results (refer to equation (2.7)). The uppermost subplot in Figure 3.3 shows the magnitude of dynamic phasor for phase A stator current. The dynamic phasor solution gives the envelope of the time domain results. The reconstructed instantaneous time domain values were identical to those obtained with the EMTP simulation with the synchronous machine represented with our own VBR model implementation. The EMTP and SFA programs for simulating the single machine dynamics are both developed in MATLAB script. The total CPU time for the SFA responses on a 1 .83G Hz PC was 7.0 s. B. Simulation of Fast Transients The second simulation is a three phase fault for a single machine infinite bus system. The parameters of the hydro turbine generator used for this case are also those in Appendix A. The system is initially operating with Tm = 0.85 p.u. and Vj is applied to the generator terminals. Then, at t = = 1.0 p.u. At t =1.0 s, a three-phase fault 1.4 s, the fault is removed with the simulation continuing to t = 8.0 s. The simulation is then repeated with the fault cleared at t = 1.5 s. A & of 0.Sms was used in order to capture the fast transients caused by the fault. When the fault is cleared in 0.4 s the system remains stable. However, when the fault is removed in 0.5 s the system loses synchronism. By repeating the simulation, we found that the critical fault clearing time to prevent instability is close to 0.462 s. The SFA-based model gave identical results to those of the instantaneous time domain simulation. The time domain results are from the same MATLAB program mentioned in Section A. The CPU time needed to calculate the SFA responses for the stable case and the unstable case were 6.13 s and 6.14 s, respectively. The simulation results of the three-phase fault case are illustrated in Figure 3.4-Figure 3.6. Zoomed 48 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis in views of the unstable case with fault clearing time of 0.5 s are shown in Figure 3.6, which has clearly captured the exponentially decaying of dc offsets in three-phase stator currents. Simulation results for a single-phase-to-ground fault are presented in Figure 3.7. This figure shows the stator currents zoomed-in segments for ‘a ‘as, ‘bs and with a fault occurring at t = 1.1 s on phase A. Two are embedded in the top subplot to depict the relationship between dynamic phasor and time domain results. Case studies (a) and (b) show that the SFA model is adequate as a general purpose machine model. It allows us to use large time steps to simulate slow transients (case (a)), while it can also accurately capture the fast transients if an appropriate small time step is used (cases (b)). C. Comparison between SFA and VBR Models In order to further illustrate the advantages of the SFA solution domain and of the proposed SFA phase-coordinates salient-pole synchronous generator model, a third case is simulated comparing the SFA solution with a conventional EMTP solution that uses the improved VBR phase-coordinates salient-pole synchronous generator model. In this example, a steam turbine generator is connected to an infinite, bus that supplies rated three-phase voltages. The system is running in no-load steady state with Tm = 0 and = 1.0 p.u. The torque Tm is stepped up to 1.1 p.u. att=0.1 s. The SFA model was used with different time steps At = 0.lms, At = 1 ms and At The study was repeated using the VBR model with At = 10 ms. 0.1 ms. The VBR model was carefully verified in [37] against the standard EMTP model and its results are assumed to be the reference accurate time domain results. The results in Figure 3.8-Figure 3.12 show that the proposed SFA model can accurately simulate the system dynamics with a very large time step of 10 ms. In fact, the simulations show that if the time step is made as large as 15 ms, we still get accurate and stable results. Table 3.1 indicates the CPU time needed for the VBR and SFA simulations. For the same accuracy, the VBR model needed, a At = 0.1 ms and a total CPU time of 13.266 s, while the SFA model needed a At = 10 ms and a total CPU time of 0.26 s. A zoomed in view of phase A stator current is shown in Figure 3.9. The SFA model was in this case 50 times more efficient than the VBR model. As reported in [37], for the same accuracy the EMTP VBR model can use a 50 times larger time step than the conventional dqO model [41] (e.g. 0.5 ms vs 10 its) at the expense of about 5 times the computational time per solution step. For this rotor dynamics 49 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis example, the VBR model in the EMTP solution frame was, therefore, about 10 times more efficient than the dq0 EMTP model, while the new proposed SFA model in the SFA solution frame was about 500 times faster than the dqO EMTP solution. This example illustrates the very high advantage of the SFA model for slow system dynamics. Table II also shows that SFA simulation is computationally more intensive than EMTP simulation for the same At step. This is due to the higher computational cost of operating with complex numbers (SFA) as compared to real numbers (EMTP) and also to the overhead involved in transferring between shifted frequency domain and normal time domain. These higher costs, however, are more than compensated for by the gains achieved by using a much larger integration step At. It is also interesting to compare the numerical stability of the SFA machine model in the SFA frame with that of the VBR model in the EMTP standard frame. Figure 3.13 shows that if we run an EMTP simulation with the VBR model using the time step of 10 ms that we use for the SFA model, the VBR model results are no longer numerically stable. It is known that the integration rule can only give reasonably correct results when the maximum frequency in the simulation is at least five times less than the Nyquist frequency fNyquist words, the time step At should be less than 1/(5x2x60) = = ) [30]. In other m 1/(2f 1.67 ms for an EMTP simulation of fundamental frequency dynamics. Table 3.2 shows the VBR model with a 5 ms time step has 11% numerical error, which makes the result no longer usable. When the time step used in the simulation is too large, i.e. At = 10 ms, the simulation in the unshifted time domain will become numerically unstable. On the other hand, as can be seen in Table 3.2, the SFA model with 10 ms time step provides much more accurate result than the VBR model with 1 ms time step, and at the same time is five times faster than the latter. Table 3.1 CPU Times for 4s Simulation in Case C TimeSteps 0.lms* lms lOms* *values for same accuracy ** VBR Model 13.266s* 1.328s ** numerically unstable 50 SFA Model 27.012s 2.668s 0.260s* Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis - — — EMTP SOICCOOn —SFASOkDOn - 1 2 4 6 8 10 12 14 16 18 0 2 4 I 6 I 8 I 10 12 14 16 18 Li a. 0€ a. I 20 H 0.5 0, 2 4 6 8 10 12 14 16 18 2 4 6 8 10 12 14 16 18 2 4 6 8 10 12 14 16 18 I C I C I 0 2 4 6 8 10 12 ‘‘5 2 4 6 30 €1) 20 0) €1) l0 2 1.00. 1 .DOC I .0O .1.00: 0.99€ 0.: -02 -0.6 -0.1 I_ 14 16 19 20 OZ 0.15 0,1 — 005 O .0,05 -0.1 .0.15 8 10 12 14 16 18 2C Figure 3.3 Simulation Results with the SFA Model for Field Yoltage and Mechanical Torque Changes (At =7 ms) 51 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis - EMTP SoIutloj Co CU I 150 a) a) 0) a) •0 8 D ci 2 Figure 3.4 Simulation Results with the SFA Model for a Three Phase Fault (At = 0.5 ms): A Stable Case (fault removed at 1.4s) 52 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis 15 I I I I EMTP Solution I isFASoluuoa ic 0 1 I I I I I I 2 3 4 5 6 7 2 3 4 5 6 7 5 6 7 4000 a) 20) 3000 _ 2000 a) to 1000 t I 2 3 a. $ Figure 3.5 Simulation Results with the SFA Model for a Three Phase Fault (At (fault removed at 1.5s) 53 = 0.5 ms): An Unstable Case Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis 0. Ill \ I,I Ij 1111 Il’? I II 1 — — 1 EMTPSob9o I-si C Ii 1 l 115 12 II 0 I II l II I I! Ij 1 1] II I 1 II I, I, 1.02 1.01 0. 0.99 I 1.05 1.1 1.15 1.2 .25 1.3 1.35 .4 1.45 Time (s) Figure 3.6 Simulation Results with the SFA Model for a Three Phase Fault (At = 0.5 ms): Zoomed-in View of a Portion of Results in Figure 3.5 cI. 0 Figure 3.7 Simulation Results with the SFA Model for a Single-Phase-to-Ground Fault (At = 0.5 ms) 54 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis - - StorC,rentl(t) x 10 Aoo,r.t, Tin,. Oonn,in Ren,,It —SFA5t0.1o-3 SFA8tle-3 — — 2 SFAt0.1.-3. SFAn,1o4n 15 05 0 -05 —1 -1.5 -2 1.5 05 25 35 Inn, (n) Figure 3.8 Simulation Results by SFA: Phase A Stator Current Figure 3.9 Simulation Results by SFA: Zoomed-in View of Phase A Stator Current (A) 0 P.IdCmnti 10 1.5 SFAI,r10.-3 t 05 0.5 1 1.5 2 2.5 3 3.5 fin,, (n) Figure 3.10 Simulation Results by SFA: Field Current 55 Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis Electoomagneto Torque T 0 z I—a) 0.5 1 1.5 2 1W,. (s) 2.5 3 3.5 Figure 3.11 Simulation Results by SFA: Electromagnetic Torque Rotor Speed e, Co CO e Wee(s) Figure 3.12 Simulation Results by SFA: Rotor Speed Table 3.2 Accuracy Comparison between SFA Model and VBR EMTP Model in Case C Model Time Step CPU Time VBR VBR SFA lms 5ms lOms* 1.328s 0.271s 0.260s 56 Numerical Error in Stator Current 1.02% 11.0% 0.21% Chapter 3 Synchronous Machine Modelling Based on Shifted Frequency Analysis 4: n,o () Figure 3.13 Time Domain Results Using the VBR Model (At = 10 ms) 3.4 Summary This chapter has presented a very efficient synchronous machine model for the simulation of slow system dynamics using EMTP modelling in the Shifted Frequency Analysis framework. Instead of instantaneous time domain variables, the model uses time varying complex variables (dynamic phasors). The SFA model results in an EMTP-type equivalent circuit that retains the numerical properties of the EMTP solution at much larger integration steps. The synchronous machine model implemented in this chapter for the SFA framework is based on the VBR synchronous machine model, which is a very efficient implementation of a phase-coordinates synchronous machine model. Working in phase coordinates provides a more stable machine model than the traditional dqO model by avoiding predictions at the interface with the external phase-coordinates power system network. The SFA framework reduces the size of the equation set that describes the synchronous machine behaviour for large efficiency gains. During system dynamics around the fundamental 60 Hz operating state of the system, a very large time step can be used to capture time domain results without loss of accuracy. Speedups of fifty times over traditional EMTP simulation were obtained for a case of mechanical torque step changes. This proves the feasibility of using EMTP solutions in the Shifted Frequency Domain for on-line assessment of large-scale power system dynamics. The SFA model is also capable of simulating fast transients caused by topological changes in the electrical network by using smaller time steps, but the extra overhead makes it slower than traditional EMTP for these situations. 57 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis 4.1 Introduction The intent of this chapter is to extend the SFA method to induction machine modelling. Induction machines form a large portion of power system loads and serve as some Independent Power Producer-owned generators. An accurate simulation of their dynamic behavior is therefore required for voltage stability and transient stability, etc [7] [8]. The efficient and stable SFA model for induction machines is well suited for such simulations in the SFA domain. This chapter proposes a new efficient SFA model for induction machines, by combining the SFA method with the phase-domain EMTP induction machine model. An equivalentreduction technique, which reduces the number of equations of stator, rotor, and flux linkages into only those of stator variables without sacrificing numerical accuracy, is used in deriving the SFA induction machine model. It allows the use of large time steps in simulating slow dynamics, but the SFA model can also accurately simulate the fast transients in the induction machine, which makes it a general purpose model. Section 4.2 first proposes a new phase-domain induction machine model used in EMTP. The philosophy inside the proposed equivalent-reduction (ER) approach will be described, and case study results will also be discussed in this section. Section 4.3 extends the Shifted Frequency Analysis method to the ER-based induction machine model, resulting in a very 58 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis efficient, numerically stable model that is a valuable component for real-time EMTP simulations. Section 4.4 summarizes the new contributions made in this chapter. 4.2 Equivalent-Reduction Approach to Induction Machine Modelling in EMTP Induction machines form a large portion of power system loads. Accurate EMTP simulation of their transient behavior may be required for voltage stability, transient stability, power quality analysis, and motor starting calculations, etc [8][44]. An efficient and stable EMTP model for induction machines is therefore of importance. Induction machine models being used in EMTP-type simulators are mainly dqO-domain models [41], [45]-[48]. In terms of the ways they interface with EMTP, the dqO-domain models [49] can be further classified into several groups. One group of dqO models, e.g. the universal machine model [50], utilizes the compensation technique where the external system is modeled as a three-phase Thevenin equivalent circuit and solved simultaneously with the machine equations to obtain machine variables. Then the system node voltages are updated by using linear superposition to take into consideration the machine effects. In order to use the Thevenin equivalent and linear superposition, the external system seen from the machine terminals has to be linear. Therefore, a distributed-parameter stub line has to be inserted to separate the machine model from other machines or nonlinear elements, which usually requires very small integration time steps [45][46]. The synchronous machine model in [41] can be directly converted to the induction machine model by feeding it with modified data entries [46]. However, this model requires predictions of flux linkages and mechanical variables. With this approach, the dqO model is transformed to a phase quantities model, after application of the trapezoidal rule and equations reduction to stator variables only. To keep the equivalent resistances constant and not rotor-angle dependent, the d- and q-resistances are averaged, and a correction term is added to compensate for this error. This correction term requires .prediction of the stator currents. All of these predictions require a smaller time step than actually needed, to avoid numerical instabilities. Other dqO models include the voltage-behind-subtransient-reactance (and armature resistance) model and Norton-equivalent-based model 59 [47], etc, which have similar Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis disadvantages as described above. A detailed review of dqO machine models can be found in [47]. To overcome the numerical instability and the limitation of smaller time steps in dqO models, the phase domain model has been proposed [35][51]. Based on the coupled-circuit induction machine equations [43], a phase domain model can be directly coupled into system equations; thus a simultaneous machine-network solution is achieved. This chapter proposes an efficient phase-domain EMTP model for induction machines based on the equivalent-reduction (ER) technique. By using ER, the equations of stator, rotor, and flux linkages are reduced to only those of stator variables, without sacrificing numerical accuracy. The ER model is more efficient than and as accurate as other phase domain models in the EMTP. 4.2.1 Equivalent-Reduction (ER) to Stator Quantities As seen from the terminal nodes, a power system component can be described by the terminal variables. A power system component can often be modeled as a set of linear (time-varying) differential equations (DEs). When the differential equations are discretized, the resulting circuit can be equivalenced to a simple equivalent conductance matrix and a current source vector that retain only the terminal variables. The ER method will significantly improve the simulator efficiency when the reduced nodal conductance matrix of the component is used quite often, for example, in a step-by-step EMTP simulation. The following subsection 4.2.2 derives the so-called ER induction machine model for the EMTP simulation. 4.2.2 Induction Machine Modelling Based on ER Technique A. Induction Machine Model The original voltage and flux linkage equations of the induction machine can be expressed in the arbitrary reference frame as - 60 Induction Machine Modelling Based on Shifted Frequency Analysis Chapter 4 Vqaos =rsiqaos +O•)Idqs 0 = rriqaor + P’s 1 co — “r)dqr (4.2) + Pqdor 1 FcL:’ = KLK; qaos KrLrK;’Jiqaor] [KrLK:’ [qaor] where (a (4.1) +P%qdos (4.3) is the speed of the reference frame. K and Kr are the transformation matrices that map the stator variables and rotor variables to the reference frame [43], respectively. B. Machine Modelling via ER Rewriting equations (4.1)-(4.3) into a new matrix form, excluding the 0 equations because it is already decoupled from the d and q equations, we get 11 A F[A 21 12 1qds A 11 1 FB 21 [B 12 [iqds B = 1 22 J[Iqdr•j B 22 A ][qdrJ + [ 0 (4.4) ] where 1 L + LM 11 A = L,S+LM 1 L LM 12 A = 21 A LM = 0 1 L + LM 22 A = 1 L —w(Ll$+LM) 3 —r 11 B + LM) — 0 — —‘ 0 —coLM 0 12 coLM B = 0 21 B = 0 (COCOr)LM (Ct;—cor)LM 0 0 61 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis (OrXLlr+LM) rr 22 = (0 X1tr + LM) — — rr rr The number of equations in (4.4) can be reduced by eliminating rotor equations via ER, as follows (4.5) + Bi 1qar + Vqds 2 = BIllqds where rL’ L” —aL —r 2M L r LM 1 L — + LM L r LIr+LM —r 1 +LM L rr -L r M LM 12 oL r M , 1 L +LM r 3 +LZ L”=L, LILM — L”M T + 1 L T LM Re-arranging (4.5), including the 0 equations, and applying Park’s inverse transformation to it, we get Vb rsiabcs (4.6) + Labcslabcs + Vres where T” Labcs — L —a- 1” M L+L” Is3M 3 L + — 3 — +LZ, 5 L, —- 62 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis LMr,. = 1 + LM L 5 (o)]’ [K COTLM L+L LM rr — r LM 1 + LM L, 1 L (47) + LM 0 In equation (4.6), no information for rotor electrical part is lost. In fact, the rotor information is included in Labcs and currents ‘qdr A 22’qcfr = 5 Vre which is an equivalent seen from the stator. Rotor can be represented using the following differential equations derived from (4.4) 21 + B22qdr B iIqas + q 2 A Discretizing the above equation for qdr, we can obtain the discrete-time equations (see equation (4.9) in subsection C) that can be used to update qdr at each time step. The rotor circuits have to be retained explicitly when power electronics devices are connected to them. C. Discrete-Time ER Model Discretizing (4.6) with the trapezoidal rule of integration gives us + Labcs Jlabcs (t) + Vrcs (t) + Vb (t) = ehl (4.8) (t) where = ehl(t) — 1dabcs )‘abcs (t — At)+ Vres(t — At) — Vb(t — At) After applying Kron’s reduction formula [52] to the discrete-time form of equation (4.4), we can get qar(t) = (4.9) 1 (t)Iqds (t) + k k (t) 2 Substituting (4.9) into (4.7) results in Vres = (4.10) k(t)Ib (t) + eh Q) 2 The equivalent circuit of the induction machine therefore becomes Vabos (t) = (4.11) R eq’ abcs (t) + Ch (t) where 63 Chapter 4 Req = eh (t) Induction Machine Modelling Based on Shifted Frequency Analysis 3 r = L abcs + Chl k(t) + (Equivalent resistances) (t) + Ch2 (t) (History term) The factor matrices and history terms in (4.9) and (4.10) can be expanded as follows ( AIV’ At 22 22 A _B -A (t)___J 21 21 +B (t)— At[co(t)—co(t)](L, Air —i- +LM) 2 + + LM — LM (t) = A 2 k 22 _B Q)4J.J [(A 22 21 +B l(t_At)•JIqds(t_At) 2 (t_At)4Jiqcir(t_At)] 22 22 +B +(A Air __: + Lir + At[a(t) LM — Air 2 CUr (t)Lir 2 + 1 L + + LM) LM LM LM qds 1 (t — At) 1 —-+L +LM At[CO(tIXt)COr(tAt)](Lir +LM) At[(D(tAt)COr(tAt)](Llr+LM) —-+LIr +LM 2 2 The time-varying factors k(t) and eh2(t) in (4.10) are expressed as 64 Induction Machine Modelling Based on Shifted Frequency Analysis Chapter 4 k(t) = [K(o)]_1[k3(t) (t) 2 eh = [K(e)p[jt)] 1 L (t)= 3 k Lir — L, + LM I — LMrr 1 +LM L COr(t)LM L Mr r O)Q)L r M O)r(t)LM + LM O)(t)LM (t) 1 k 1 +LM L + LM LMrr Q) = 4 k ]o)] 0 ( 5 [K L M rr 2 (t) k T lr It should be noted that, although the equations for the rotor electrical part are eliminated, the rotor currents could be updated at each simulation step via (4.9). Refer to [9] and [511 for the discrete-time equations of the mechanical part. The ER induction machine model is actually a full-order accurate phase-domain model because the information for the rotor electrical part has already been integrated into the stator part equations. Also, because Park’s transformation is implicitly used in the ER model, the elegance and simplicity in the discrete-time equivalent circuit are achieved. With a lower number of equations in the model, it is simpler and faster than traditional phase-domain models in the EMTP, without compromise in accuracy. 4.2.3 Simulation Results Three test cases are simulated to illustrate the efficiency and accuracy of the ER model. The first case illustrates free acceleration characteristics of a 3-hp induction machine and a 2250hp induction machine. The second case is a load torque change test which shows the slow dynamics of the 3-hp machine. The third case simulates the fast transients caused by a three phase fault. The induction machines parameters for the three cases are shown in Table A-TI taken from [43]. Rotor parameters have been converted to the stator side. As reported in [37], the voltage-behind-reactance (VBR) model [36] is more efficient than traditional phase-domain EMTP models and more stable than dqO-domain models; hence a discrete-time VBR induction machine model used in EMTP [53] was chosen to compare with the 65 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis ER model. To facilitate the comparison of CPU times, both machine models are deliberately implemented in the language of MATLAB® script. A. Simulation of Free Acceleration The start-up transients of the 3-hp and 2250-hp induction machines during free acceleration from stall are simulated. The total simulation times are 1 s for the 3-hp induction machine and 3 s for the 2250-hp induction machine, respectively. In this simulation, a time step &= 500 ,us is used. Figure 4.1 illustrates the torque-speed characteristics of the 2250-hp induction machine during free acceleration. The simulation results, for the 3-hp induction machine are shown in Figure 4.2. The variables plotted in this figure are ‘as (phase A stator current), ‘m (magnetizing flux linkage), Te (electromagnetic torque) and mechanical rotor speed. These machine variables are observed in the rotor reference frame; results viewed in other reference frames are omitted here. It can be seen that the time domain values obtained from the EMTP simulation with the ER induction machine model are practically identical to those with the VBR model. In Figure 4.1 and Figure 4.2, the results from the ER model and the VBR model are so close that they become indistinguishable from each other. Numerical tests also indicate that both ER and VBR model, when a large time step of 1 ms is used, can still get accurate, stable, and identical solutions. x i0 4 .j’ 2 —ER Result VBR Result - - 2 I, 0 200 400 600 800 1000 1200 1400 1600 1800 Speed (rpm) Figure 4.1 Torque-Speed Characteristics during Free Acceleration of a 2250-hp Induction Machine (At = 500 js) 66 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis 0.5 the (s) Figure 4.2 Dynamic Performance of a 3-hp Induction Machine during Free Acceleration (At = 500 Its) B. Load Torque Change Test Slow dynamics during step changes in load torque are simulated in this case. The 3-hp machine originally operates in steady state with no load. A mechanical torque Tm of 12 Nm is applied at t = 2.05 s. Then Tm is reversed to —12Nm at t = 2.5 s. The simulation time is 3 s with a time step of 500 ps. Figure 4.3 shows the phase A stator current, slip, real power, and reactive power during load torque changes. The machine variables are observed in the synchronous reference frame. Comparisons of the EMTP simulation results obtained from the ER model with those from the VBR model show again that the results are identical. 67 ___________ Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis —00 2.1 2.2 2.3 2.4 2.5 2.6 2.8 2.7 2.9 300 2000 -- ER Result VBRResuH 1005 -1005 -2005 11* I I 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 2.1 2.2 2.3 2.4 2.5 Ume (s) 2.6 2.7 2.8 2.9 3 200C 1955 1905 > O 1855 1805 1755 1702 Figure 4.3 Dynamic Performance of a 3-hp Induction Machine during Step Changes in Load Torque (At = 500 is) C. Three Phase Fault Test In this simulation of fast transients, the 2250-hp induction machine is first operating under rated conditions with a load torque equal to TB. Here TB is the base torque defined by TB = where PB is the rated power output of the machine and b 0 is the rated speed of the /p)cob 2 ( machine. A 3-phase fault is then applied at the terminals at t = 6.1 s. After 6 cycles, the fault is cleared. The simulation time is 10 s with a time step of 100 ,us. Figure 4.4 illustrates that the proposed ER model can simulate fast dynamics of the machine as accurately as the VBR model. 68 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis ER Reso1 VBR Res,1t - 5000 - -5000 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 E Figure 4.4 Simulation Results for a 3-Phase Fault at the Terminals of a 2250-hp Induction Machine (At = 100 its) Cases A to C show that the ER model is as accurate as the VBR model for the same At step. It is also interesting to explore the numerical accuracy and robustness of the ER model by using different At steps in the fast transients simulation of Case C. Results obtained from the VBR simulation [53] with a very small At = 10 ps are assumed to be the reference accurate time domain results. A zoomed-in view of phase A stator current is shown in Figure 4.5. Figure 4.5 indicates that, if we run an EMTP simulation with the ER model using a time step of 1 ms, relatively accurate results are still obtained. Even with a time step of 1.5 ms, the numerical solution remains stable. 69 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis Figure 4.5 Simulation Results for a 3-Phase Fault with Different Time Steps Table 4.1 summarizes the CPU times needed for the VBR and ER simulations in cases A to C. The timing in Table 4.1 was obtained on a 1.83G Hz PC. It can be calculated that the ER is usually 33%-36% faster than the VBR model for the same accuracy. The efficiency of the ER model is due to its concise form and the low computer burden in updating electrical variables and flux linkages. The ER model has a further advantage in addition to the faster speed. Instead of being specifically designed for the machines to achieve efficient phase domain models, the ER technique used in the induction machine modeling is a general-purpose approach. ER can also be used for building more efficient models for other types of machines, transformers and power system elements, as long as they can be mathematically represented by differential equations with mutually coupled (or decoupled) state variables. Since the ER-based model proposed in this chapter is actually a full-order model, it is different from the Model-Order-Reduction-based models. Model Order Reduction (MOR) approaches have been proposed, either for constructing a reduced-order model of the machine itself by neglecting stator flux derivatives, or for constructing a reduced-order model of the external system by using techniques such as Krylov-subspace or Fourier methods [54]-[57]. MOR methods improve the simulation speed; however, undesirable inaccurate step-by-step simulation results will be produced, especially from the latter one. 70 Induction Machine Modelling Based on Shifted Frequency Analysis Chapter 4 Table 4.1 CPU Times for Simulations Test Case Case A (3-hp induction machine, * T= is, IXt=500,us) Case A (2250-hp induction machine, us) 4 T=3s,At=500 Case B (3-hp induction machine, us) 4 T=3s,At=SOO Case C (2250-hp induction machine, T=iOs,At=500ps) Case C (2250-hp induction machine, T lOs,At50 us) 4 * VBR Model ER Model 0.875 s 0.556 s 2 601 . 1 698 2 606 . 696 8 391 . 5 562 82.796 55.031 . . T stands for the total simulation time. 4.3 Induction Machine Modelling with SFA 4.3.1 Induction Machine Modelling Based on SFA To build the dynamic-phasor-based stator equations, the stator variables are first expressed as ‘abcs (t) = V abcs (t) abcsl = (t) C05 — v abcsl (t) cos o t ‘abcsQ — (t) sin w,t (4.12) v abcso (t) sin OJJ t (4.13) Substituting (4.12) and (4.13) into (4.6), we get abcsI 1 = 5 r — (t) cos (J) t 5 [Iabcd — VbQ (t) cos cot — P’abcso (t) sin ot (t) sin co t 5 I abcsQ (t) sin th t] + ‘‘abcs [pi (t) cos cost 5 — abcsQ — abcs! 5 sin oi (t)o t 5 (4.14) 5 cos oJ (t)co t] 5 Then we apply the Hilbert transform to both sides of (4.14) and construct the analytic signal, which leads to Vabcs (t)e°’ = 5 r ‘abcs where VabL,s (t) and (t)e°’ + L abcs [Plabcs (t) + f cOSI abcs (t)] e°’ + V (t)e° t ‘abcs (t) are dynamic phasors of stator voltages respectively, and Vre, (t) is the dynamic phasor of v (t). 71 Vabcs (t) (4.15) and currents ‘abcs (t), Induction Machine Modelling Based on Shifted Frequency Analysis Chapter 4 Applying frequency shifting as described in (2.7), we can obtain the dynamic-phasor equations for the stator part as 55 r Iabcs(t)+Labcs[PIab ( t) + JüSIabcS(t)] Vb(t) (4.16) + Vres(t) 4.3.2 Discrete Time Model Discretizing (4.16) with the trapezoidal rule of integration gives us Vb(t) = [r +--+foJLb ] (4.17) Is(t) +V(t)+Ehl(t) where Ehl (t) + js Jtabcs] ‘abcs (t — At) — Vb (t — At) + ‘res (t — At) [i. + From the discrete-time equation of time equation of S 8 r (see Section 4.2), the corresponding discrete- Vres (t) (r) can be derived as Vres = k(t)Iabcs (t) + E$h (t) e1(00t) e1(0__2) 1 eo23) ej(ot232) 1 j(9—eot+22r/3) j(O—aI+2,r/3—nr/2) 1 . — k(t)Iabcs(t) + e e (t)] 4 [k (4 18) (t) are explained in Section 4.2. 4 The factors k(t) and k By substituting (4.18) into (4.17), we get the equivalent circuit of the induction machine (4.11) in the SFA domain as follows Vabcs(t) _ (4.19) Req Iabcs(t)+EhQ) where 5+ Req = r + Eh (t) = Ehl (t) + Eh 0’) + k(t) (Equivalent resistances) (History term) The SFA equivalent circuit for the induction machine is illustrated in Figure 4.6. 72 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis Figure 4.6 SFA Equivalent Circuit of an Induction Machine Discretizing the rotor mechanical part equations [9], we can get the discrete time equations of the mechanical part (4.20) O)r(t) = cr(t—At) f [. I (t) = Or (t) Or (t = — (t)lq At) + (4.21) (t) + Aq 3 (t)idS (t)] [C0r (t) + (4.22) cUr (t At)] — where p is the number of poles, K(Or(t))Re[Iabcs(t)e’] ‘qdOs Lilq ± mq = Lislqs + LM (1qs + qr) Aas and = dS 1 LIS Iqar(t) + kd = LiSidS + LM (Ids + dr) can be updated by using (4.9). The more complex mechanical part equations based on the multi-mass model can be found in [9]. 4.3.3 Simulation Results The same three test cases described in Section 4.2.3 are simulated on a 1.83 GHz PC to illustrate the efficiency and accuracy of the SFA model. The first case illustrates free acceleration characteristics of a 3-hp induction machine and a 2250-hp induction machine. The second case is a load torque change test which shows the slow dynamics of the 3-hp machine. The third case simulates the fast transients caused by a three phase fault. The following subsections A, B, and C explain that the SFA-based model is an accurate and general-purpose one. Subsection D explores the very high efficiency of the SFA model. The author has developed 73 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis a MATLAB program based on the EMTP algorithm and the ER machine model. All SFA results in the following subsections have been compared with the results from this EMTP program in the same MATLAB environment. The induction machine parameters for the three cases can be found in Appendix A. Rotor parameters have been converted to the stator side. A. Simulation of Free Acceleration The start-up transients of the 3-hp and 2250-hp induction machines during free acceleration from stall are simulated. The total simulation times are 1 s for the 3-hp induction machine and 3 s for the 2250-hp induction machine, respectively. In this simulation, a time step of & = 500 ps is used. Figure 4.7 illustrates the stator current of the 2250-hp induction machine during free acceleration. Note that the SFA model can produce both dynamic phasor and time domain results. The simulation results for the 3-hp induction machine are shown in Figure 4.8. The variables plotted in this figure are ‘as (phase A stator current), Vm (magnetizing flux linkage), Te (electromagnetic torque) and mechanical rotor speed. These machine variables are observed in the rotor reference frame; results viewed in other reference frames are omitted here. As is shown in Figure 4.8, the time domain values reconstructed from the dynamic phasor results are identical to those obtained from the EMTP algorithm using a ER model. The CPU times used for simulations are listed in Table 4.2. Figure 4.7 Stator Current during Free Acceleration of a 2250-hp Induction Machine (At 74 500 jts) Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis Time (s) Figure 4.8 Dynamic Performance of a 3-hp Induction Machine during Free Acceleration (At = 500 jts) Table 4.2 CPU Times for Simulations Free Acceleration Case SFA Model 0.728 s 3-hp induction machine, T = 1 s, tt = 500 4 us 2250-hp induction machine, T 3 s, & = 500 ,us 2.43 9 s * T stands for the total simulation time. B. ER Model 0.556 s 1.698 s Load Torque Change Test Slow dynamics during step changes in load torque are simulated in this case. The 3-hp machine originally operates in steady state with no load. A mechanical torque Tm of 12 Nm is applied at t = 2.05 s. Then Tm is reversed to —12Nm at t = 2.5 s. The simulation time is 3 s with a 75 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis time step of 0.5 ms for the ER model, and a time step of 5 ms for the SFA model. Figure 4.9 shows the phase A stator current and the slip during load torque changes. Figure 4.9 shows that simulation results from the SFA model with At = 5 ms are as accurate as those from an EMTP solution with At = 0.5 ms. The CPU time for a 3 s simulation with the SFA model is 0.266 s, while the simulation with the ER model uses 1.696 s. Therefore, to achieve the same numerical accuracy, the SFA induction machine model can be more than 5 times faster than its corresponding ER model in this case. 2.5 Time (s) Figure 4.9 Dynamic Performance of a 3-hp Induction Machine during Step Changes in Load Torque (RE At = 0.5 ms, SFA At =5 ms) C. Three Phase Fault Test In this simulation of fast transients, the 2250-hp induction machine is first operating 76 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis under rated conditions with a load torque equal to TB. A 3-phase fault is then applied at the terminals at t = 6.1 s. After 6 cycles, the fault is cleared. The simulation time is 10 s with a time step of 500 ps. Figure 4.9 illustrates that the proposed SFA model can simulate fast dynamics of the machine as accurately as the ER model. The CPU times are 7.258s and 5.562 s for SFA simulation and ER simulation, respectively. Time (s) Figure 4.10 Simulation Results for a 3-Phase Fault at the Terminals of a 2250-hp Induction Machine (At 500 J4s) D. Comparison between SFA and an EMTP Phase-Domain Model Cases A to C show that the SFA model is as accurate as EMTP phase-domain models for the same time step i.t. Now we will carefully explore the numerical accuracy and robustness of 77 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis the SFA model, by using different time steps At in the slow dynamics simulation. The 2250-hp machine originally operates in steady state with rated torque of 8900 Nm. The machine dynamics are recorded when a sudden reverse of mechanical torque to -8900 Nm is applied at t = 6 s. The total simulation time is 10 s. The SFA model was run with different time steps At = 0.lms, At = lms, At = 5ms, and At = lOms. The results obtained from an original phase-domain model (refer to p. 142-147 in [43]) with a small At = 0.1 ms are assumed to be the reference accurate time domain results. A zoomed-in view of phase A stator current is shown in Figure 4.11, which shows that the SFA model can accurately simulate the system dynamics with a large time step of 10 ms. In fact, even with a time step of 15 ms, the numerical solution with SFA remains stable. On the contrary, an EMTP simulation in the time domain with a time step of more than 4 ms became numerically unstable, as is shown in Figure 4.12. 5.8 6 6.2 6.4 6.6 Time (s) 6.8 Figure 4.11 Simulation Results by SFA 78 / 7.2 Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis $ Figure 4.12 Time Domain Results from an EMTP Algorithm Implemented in MATLAB Table 4.3 summarizes the CPU times needed for the SFA and ER simulations in this case. It can be seen that the SFA model was about 10 times faster (0.370 s vs 2.80 1 s) than the EMTP phase-domain model (ER model) to achieve a comparable accuracy. In certain situations, for instance, when simulating a system with power electronics components, an EMTP-type simulator may have to use very small time steps, for example, less than 0.1 ms. As can be seen from Table 4.3 that a SFA model with zt = lOms can be about 75 times faster than an EMTP model with & = 0.lms, for reasonable accuracy. Therefore, it is very promising to apply the SFA model in hybrid simulators, where the power electronics components as well as the network can be simulated with a small time step, while the SFA can be used for machine models with a larger time step. Further, the new SFA model has the capability to be extended to model the saturation and deep bar effects in induction machines. The SFA modeling framework is shown to be successful for machine modeling as a new method, which opens the door for more future research such as applying the SFA to wind turbine generator modeling, and developing SFA based EMTP simulators or hybrid simulators, and so on. Table 4.3 CPU Times for Simulations 2250-hp induction machine Torque Change Case T=lOs,&0.lms T=lOs,&=lms T=lOs,&5ms T= 10 s, & = 10 ms * Numerically unstable. 79 SFA Model ER Model 36.108 s 3.725s 0.738s 0.370 s 27.328s 2.801s -- Chapter 4 Induction Machine Modelling Based on Shifted Frequency Analysis 4.4 Summary (1) An Equivalent-Reduction-based induction machine model is proposed. With a lower number of model equations, it is actually a full-order phase-domain model but is simpler and faster. Park’s transformation is implicitly used in the ER model in order to maintain the elegance and simplicity in the discrete-time equivalent circuit. As a non-ad hoc approach, the ER can also be used for building more efficient models for other power system elements. (2) Based on the ER model, a SFA model is proposed as a general-purpose model capable of simulating both fast transients and slow dynamics. Case study results have confirmed the SFA induction machine model is a valuable component for real-time EMTP simulations. It is observed that the SFA model is in excess of 70 times faster than the EMTP solution with the RE model when simulating dynamics with frequency spectra close to the fundamental power frequency. 80 Chapter 5 EMTP Implementation Chapter 5 EMTP Implementation 5.1 Introduction Since a number of power system component models in the shifted frequency domain have been built in the previous chapters, it is now time to ask how they can be used in a general purpose simulation tool based on the Shifted Frequency Analysis. EMTP, which was originally developed for calculating the transient overvoltages in transmission systems, has been significantly expanded to tracing the evolution of the system states in arbitrary multi-phase power networks consisting of all types of components. With improved functionality, accuracy and numerical stability, EMTP has become a standard tool being widely used in the power industry for system planning and designing purposes. Now EMTP is seeing broader applications in power system steady state [58] and dynamics studies and will remain one of the mainstreams in power system research. The intention of this thesis is to expand EMTP to efficiently simulate the slow dynamics in power systems, and to bridge the gap between the EMTP and a unified power system analysis tool. This means that the EMTP algorithm will be adopted in the SFA-based simulation tool. What makes a difference here is that the electrical variables are described by dynamic phasors instead of instantaneous time values. 81 Chapter 5 EMTP Implementation 5.2 Program Structure The implicit trapezoidal rule of integration, which has attractive characteristics in terms of accuracy and numerical stability, is used in EMTP and is followed in the SFA-based modelling of system components. By using the trapezoidal rule, the differential equations representing all network components are converted into algebraic relationships which can be interpreted as equivalent resistances or admittances, voltages, currents, and known history terms. Then the nodal equations of the system can be written directly as follows. [G][V(t)] = [I(t)]+ [H(t)] (5.1) where [G] is the matrix of the nodal equivalent admittances [VQ)] is the vector of the nodal voltages, which are dynamic phasors [1(t)] is the vector of the nodal currents, which are dynamic phasors [11(t)] is the vector of the nodal history terms The voltages of the nodes connected to voltage sources are known quantities, therefore the corresponding equations can be eliminated. Suppose A is the index set denoting the nodes with unknown voltages, and B is the index set for nodes with known voltages, then the nodal equations can be written in a block matrix form rGM G1EvA(t)1 [GBA GBB][VB(t)J = EIA(tHAt)1 (5.2) [IB(t)+HB(t)] Thus the unknown voltages can be obtained at time t by solving the following equations GVA(t) [IA(t) + (5.3) HA(t)]—GVB(t) In this thesis, a toolbox for simulating power system dynamics in the shifted frequency domain has been developed with MATLAB. This dynamic phasor tool consists of several files, which are depicted in the schematic structure in Figure 5.1. A. Input Data File Input data file provides the data needed for the simulation of a multi-phase power network. These include transmission line data, load data, machine data, switching operations 82 Chapter 5 EMTP Implementation data, voltage sources data, etc. All these data are documented in a MATLAB script file and most of the data are provided as the MATLAB arrays. Refer to Appendix C for the input data format. System Solver Main.m PiLineParameters.m MacParameters.m Data Processing LoadParameters.m RLCParameters.m EMTP Dynamic Phasor Simulation Tool (EMTDP) Equivalent Admittance Matrix Forming Updating Admittance Matrix and History Terms Data Input File Gsys.m UpdateG.m UpdateHist.m Case_i 3Bus.m Figure 5.1 Schematic Structure for the EMTP Dynamic Phasor Simulation Tool B. Data Processing Files When the input data are read by the main program, the data processing functions are called to establish the equivalent circuit for different power system components, and initialize their history terms for the simulation at the first time step. The equivalent circuit for different components can be found in previous chapters. The initialization in the SFA simulations is based on the snapshot method. First, let the dynamic phasor program run with zero initial conditions and reaches the steady state. Then, a snapshot of.the system is taken by saving system variables and history terms at a particular time step. These system variables and history terms from the snapshot file are fed to the SFA program as the initial conditions. With system variables and history terms initialized, the SFA simulation 83 Chapter 5 EMTP Implementation will run with a ‘flat’ start, and the advantage of using large time step in simulating 60 HZ dynamics is achieved. C. File for Building G Matrix The file ‘Gsys.m’ builds the system [G] matrix in equation (5.1) and returns the [Gj and [Gp.i] to the system solver for calculating the unknown voltages. The following explains how this simulation tool inserts power system components into equivalent admittance matrix [G] by taking the it-circuit model of a transmission line as an example. Suppose there is a three-phase it-circuit for a transmission line connecting node set J ,k 1 ,k 2 ]) as depicted in Figure 5.2. As can be seen in Figure 5.2, 3 (=[ji, j2, j3}) and node set K (=[k the it-circuit model consists of a coupled RE branch and two coupled capacitances. When discretized by the trapezoidal rule, the coupled RE branch will contribute a 3 x3 equivalent admittance matrix to [G], and contributes a 3x1 history term hpj to [H]. Similarly, the coupled capacitances will also contribute a 3 x 3 equivalent admittance matrix Gc and a 3 xl history term hc. From equations (2.21) and (2.24), it can be found that GRL will be entered into [G]j,j and [G]jçK in the system admittance matrix [G], and —G j. into [G]J,K and [G]jçj. For 1 instance, G(ll) will be added to the element [G]13, Ggi(1,2) will be added to [GJJ1i2, 2 GII(l, to ) GRI(2,1) to [G]j21, and so on. This procedure is illustrated in Figure 5.3. The MATLAB implementation of this process is shown in Figure 5.4. This routine also moves the nodes connecting to the voltage sources to the bottom of the system equation set, and then extract the sub-blocks [Gj and [Gpj] from [G]. 13 C Figure 5.2 The fl Transmission Line Model 84 Chapter 5 EMTP Implementation D. File for Updating [G] and [H] The MATLAB files UpdateG.m will update the G matrix when changes happen in the network configuration, e.g. line tripping or faults. The file UpdateHist.m will calculate the history terms at each time step and return them to the solver for calculating unknown variables. E. System Solver The system solver Main.m is the core function in the simulation tool. In each time step, this routine assembles the right hand side of (5.3), performs downward operations on it, and does the back-substitution to obtain the unknown nodal voltages. After the unknowns are found, the routines for updating history terms and/or admittance matrix are called, in preparation for the calculations at the next step. The simulation run will continue until the total simulation time is reached. A schematic flow chart of the system solver is presented in Figure 5.5. 3 2 1 k 111213 +G -G --- \___ / jl 32 33 32 33 + \ hcj +Gc + 1 G k 1 2 Ic 3 k /hRL k 1 2 Ic 3 Ic / + -G +G hc K [hj [Gj Figure 5.3 Contributions of the fl-Circuit Transmission Line Model to the Nodal Admittance Matrix 85 Chapter 5 EMTP Implementation % Pi Lines Contribution to G Matrix for kk = 1: nLine ii = (LineFromBus(kk)- 1 )*3 + 1; jj = (LineToBus(kk)1)*3 + 1; Gsystem( ii: ii + 2,11: ii + 2) = Gsystem( ii: ii GRL(:,:, kk) + GC(:,:, kk); Gsystem( ii: ii + 2,jj: jj + 2) = Gsystem( ii: ii GRL(:,:, kk); Gsystem( jj: jj + 2,11: ii + 2) = Gsystem( jj: jj GRL(:,:, kk); Gsystem(jj: jj + 2,jj:jj + 2) = Gsystem(jj:jj GRL(:,:, kic) + GC(:,:, kk); end + 2,ii: ii + 2) + 2,jj: jj + 2) + 2,ii: ii + 2) + 2,jj:jj + 2) + - - + Figure 5.4 The MATLAB Code for Inserting H-Circuit Transmission Line Model into G Matrix 5.3 Test Cases Four test cases are tested in the SFA-based EMTP simulator developed in this chapter. A. Radial Transmission Line Case The first case is a radial network consisting of a voltage source, a three phase load and a double-circuit overhead line between the source and the load. The one-line diagram is shown in Figure 5.6. The transmission line is energized at t = 0. At t = 3s, one of the parallel circuits is tripped. The total simulation time is 5 seconds. A time step of 5 ms is used in this simulation. Detailed data for this case are included in Appendix C. Figure 5.12-Figure 5.12 illustrates the three phase voltages at the load node, including the zoomed-in view for these voltages in the chosen time interval from 2.9s to 3.ls. A MATLAB program based on the EMTP pi-circuit model is used to generate the time domain results for comparison. It can be seen that the EMTP solutions in the shifted frequency domain and those in the time domain are almost identical, and the dynamic phasor result is the envelop of the time domain curve. 86 __________ Chapter 5 EMTP Implementation ( Input Data Find Steady State Solution to Initialize History Terms and Variables at I = Jr Build GAA, GAB for Transient Solution 1. Factor G into LU Form Jr Start Time Step Simulation = Any Switching Modify and ‘‘ Evaluate Current and Voltage Sources, Update the Right Hand Side of(5.3) Output Results J,N I=t+Ai Perform Downward Operations on RHS of(5.3), Do Back-substitution to Solve Nodal Voltages ‘I, Update History Terms Figure 5.5 Flow Chart for Dynamic-Phasor-Based EMTP Simulator Source Busi Bus2 Figure 5.6 One-line Diagram of a Radial Test System 87 R,X Load Chapter 5 EMTP Implementation ‘C I SFA Solution ‘nfl —l 1 0.5 1.5 2 2.5 Time (s) 3 3.5 4 45 5 Figure 5.7 Phase A Voltage at the Load Node (At = 5 ms) Figure 5.8 Zoomed-in View of Phase A Voltage at the Load Node (At = S ms) ‘ r 0 0.5 1 1.5 2 2.5 Time (s) 3 3.5 Figure 5.9 Phase B Voltage at the Load Node (At 88 4 = 5 ms) 4.5 1 5 ___________ Chapter 5 EMTP Implementation i ,i I I I I 2.92 2.94 2.96 2.98 I I .1 3.02 3.04 SFA Solution 0.5 > > -0.5 —1 2.9 3.06 3.08 3 1 Time(s) (b) Figure 5.10 Zoomed-in View of Phase B Voltage at the Load Node (At = 5 ms) x 0 0.5 I 1.5 2 2.5 Time (s) 3 3.5 4 4.5 1 5 Figure 5.11 Phase C Voltage at the Load Node (At =5 ms) x I 0.5 Li 0 -0.5 —l I I SFA Solution EMTP Solution IVVVVV\IVVVVV 2.9 - 2.92 2.94 2.96 2.98 3 Time (s) 302 3.04 3.06 108 3.1 Figure 5.12 Zoomed-in View of Phase C Voltage at the Load Node (At = 5 ms) 89 ChapterS B. EMTP Implementation Distribution Network Case The second case is a 13 bus network with the configuration adapted from the IEEE test feeder [59]. The network configuration is shown in Figure 5.13. The detailed network data are documented in Appendix C. The network is energized at t connected to node 11 is initially operating at no-load. At t = = 0. The induction machine load 3s, a rated mechanical torque is applied on the induction machine. Then one circuit of the double-circuit lines connecting node 0 and node 1 is tripped at t = 6 s. The total simulation time is 8 s. A time step of 1 ms is used in this simulation. Figure 5.14-Figure 5.18 illustrate part of the three phase voltages at three different nodes. A MATLAB program using the EMTP algorithm is used to produce the time domain results. The time domain simulation results and the dynamic phasor results are both shown in each figure to illustrate the correctness of the dynamic phasor results. From the figures, it can be seen that the increase of the induction motor load causes drops in the system voltage. Loss of one of the double circuit lines further weakens the systems and causes voltage dip problems. 0 5 4 1 2 11• Figure 5.13 One Line Diagram for the Test Feeder 90 3 ChapterS EMTP Implementation — 3 3.5 4 4.5 5 5.5 Time (s) 6 SFA Solution EMTP Solution 7 6.5 7.5 Figure 5.14 Phase A Voltage at the Induction Machine Node (At = 1 ms) > Figure 5.15 Zoomed-in View of Phase A Voltage at the Induction Machine Node (At = 1 ms) I SFA Solution EMTP Solution I I I I I 3 3.5 4 45 5 5.5 I I 6 6.5 7 Time (s) Figure 5.16 Phase B Voltage at Load Node 6 (At = 1 ms) 91 7.5 I 8 Chapter 5 EMTP Implementation > Figure 5.17 Zoomed-in View of Phase B Voltage at Load Node 6 (At = 1 ms) SFA Solution EMTP Solution .4OQO I I I I 3.5 4 4.5 5 I 5.5 6 6.5 Time (s) Figure 5.18 Phase C Voltage at Node I (At = 1 C. ms) First Benchmark System for the Subsynchronous Resonance Studies This test system was prepared by an IEEE Subsynebronous Resonance Task Force [601 as a standard test case for computer programs to simulate subsynchronous resonance phenomena. The test system consists of an 892.4 MVA turbine-generator connected through a step-up transformer to a 500 kV transmission line with series capacitor compensation. The power system at the receiving end is represented by a Thevenin equivalent circuit (infinite bus behind reactance) [61]. Figure 5.19 shows the one-line diagram for the first benchmark system. In 92 Chapter 5 EMTP Implementation Figure 5.19, all data are represented in p.u. based on 892.4 MVA and 500 kV. The generator parameters of the electrical part are listed as below Xci = 1.79 p.U. Xq X’d = 0.169 p.U. X’q = 0.228 p.U. 1 X X”q = 0.2 p.u. f X”d = 0.135 p.U. T’cio = 4.3 T”do 1.71 p.U. T’qo S 0.032 S T”qo = 0.85 S 0.05 5 Ra= 0 0.13 p.U. 60 Hz it(O) = 1 f.U. The original purpose of this case was to simulate the interaction between the mechanical torque placed on the generator turbines and the electrical torque related to the power network, and the resulting shaft torsional oscillations. Accordingly, a detailed multi-mass modcl of the mechanical shaft is adopted in [601, [611, where the generator shaft system is modeled by 6 masses including 4 turbine sections HP, IP, LPA and LPB, 1 generator and 1 exciter. This thesis focuses on the feasibility of applying the SFA method, and thus uses a simpler single-mass representation of the mechanical part. The total inertia constant of the turbine and generator is 2.894 seconds. The SFA synchronous machine model in Chapter 3 is used in this test case. It should be noted that the SFA synchronous machine model of Chapter 3 can be extended to model the multi-mass shaft system by adding differential equations for individual spring-masses to the differential equations of the shaft system. The self and mutual damping effect can also be easily included in the shaft system equations. The detailed network parameters can be found in [601 and [611. 26kV/500kV A Generator 0.14 =0.14 0 X = X = 0.50 R 0.02 1 1 =l.56 0 =0.50 X 0 R B 1 X = 0 06 — I — Xc-O.371 1 Infinite Bus X=0.04 0 = 0.04 Figure 5.19 The First Benchmark Network for Subsynchronous Resonance Studies 93 Chapter 5 EMTP Implementation In this test, a three-phase fault occurs at bus B at t = 3 s. After 4 cycles, the fault is cleared. The total simulation time is 3.5 s. A time step of 1 ms is used in this simulation. The whole system is represented using the SFA models proposed in the previous chapters. The generator terminal voltages, the voltages at bus A and bus B, and the voltages across the series capacitors are monitored during the SFA simulation. Figure 5.20 to Figure 5.23 illustrate the dynamic phasor results for the system voltages at different locations as well as the time domain results transformed back from the corresponding dynamic phasors. The dynamic phasor results show that, after the fault is applied and cleared, there are low frequency oscillations happening at the generator terminals, on the transformer sides, and across the series capacitor banks. ‘U ‘U C C ‘U F C C ‘U C ‘U time (s) Figure 5.20 Generator Terminal Voltage: Phase A (At 94 = 1 ms) Chapter 5 EMTP Implementation 0 0 .0 I Figure 5.21 Transformer High Side (Bus A) Voltage: Phase A (At = 1 ms) 3 time (s) Figure 5.22 Voltage across Series Capacitor: Phase A (At = 1 ms) These oscillations are an electrical phenomenon because the shaft system is modeled as a single mass. The slow oscillations are caused by resonances with the series capacitor, which are excited by the fault and the switching operations. The resonance mode(s) is determined by the inherent characteristics of the power network. In fact, the modes or natural frequencies can be quantitatively found by either performing a frequency scan (steady-state solutions over a frequency range) or by calculating the eigenvalues of the admittance matrix of the electrical 95 Chapter 5 EMTP Implementation network. Similar to the network resonance analysis, if the shaft system is represented by the detailed multi-mass model, we can also determine the torsional natural frequencies (eigenvalues) and mode shapes (eigenvectors) by applying the modal analysis to the shaft differential equations. If the complement of the natural frequency of the network is close to one of the torsional frequencies of the shaft system, torsional oscillations will be excited [331. This is the mechanism of the subsynchronous resonance. The dynamic phasor results obtained from the SFA simulations show that the oscillations with frequencies around 30 Hz have been excited in this benchmark system. On the other hand, the natural frequencies of the 6-mass shaft system were found to be 15.71 Hz, 20.21 Hz, 25.55 Hz, and 32.28 Hz [601. That means the complement of the natural frequency of the series-compensated network is close to the torsional natural frequencies, which would result in a subsynchronous resonance in this system and would build up torsional oscillations on the shaft. This has been verified by EMTP simulations [611. In summary, the SFA simulations are able to capture the slow oscillations in the benchmark system and can be used for subsynchronous resonance studies once the detailed shaft model is incorporated into the SFA synchronous machine models. U, 0 > U, 0 0 to time (s) Figure 5.23 Infinite Bus (Bus B) Voltage: Phase A (At D. = 1 ms) Second Benchmark System for the Subsynchronous Resonance Studies Eight years after the first subsynchronous resonance benchmark was published, the IEEE Subsynchronous Resonance Working Group proposed a second benchmark [621. The second 96 EMTP Implementation Chapter 5 benchmark is a case consisting of a 600 MVA generator connected through a step-up transformer to two parallel lines, one of which is series compensated. The compensation rate of the line is 55%. The power system at the receiving end is represented by a Thevenin equivalent circuit (infinite bus behind impedance). Figure 5.24 shows the one-line diagram for the second benchmark system. All data in Figure 5.24 are represented in p.u. on a 100 MVA, 500kV base. The mechanical part of the round rotor turbine generator is again represented as a single-mass model with a total inertia of 2.683 seconds in the SFA simulation. In [621, a more detailed multimass representation is used with 4 masses. The parameters of the electrical part are based on the machine ratings, which are given below. Xd = 1.65 p.U. Xq = 1.59 p.u. X’d = 0.25 p.u. X’q X”d = 0.20 p.U. X”q = 0.20 p.u. f T’do = 4.5 T’qo = 0.55 if(O) Tdo = S 0.04 S = T”qo 0.46 p.U. = S 0.09 s Ra = 0.0045 Xi = 3.U. 0.14 p.U 60 Hz = 1 3.U. number of poles = 2; The detailed benchmark network parameters can be found in [62] and [61]. In this test, a three-phase fault is applied at the high voltage side of the transformer connecting to bus 2 at t = 3 s. The fault clearing time is 1 cycle, and the total simulation time is 3.5 s. The time step used in this simulation is 1 ms. Figure 5.25 to Figure 5.28 show the dynamic phasor results for the system voltages at different locations together with the time domain results transformed back from the corresponding dynamic phasors. Obviously, oscillations slower than the fundamental frequency are excited in the system after the fault occurs, which can be seen from the dynamic phasors of the generator terminal voltage, the transformer high side voltage and the voltage across the series capacitor. Note that bus 1 does not see large deviations in the voltage because it is electrically close to the infinite external system, and also because the series capacitor is acting as a highpass filter that blocks the slow oscillations in the system. The test results on the second benchmark system further verify that the SFA method can capture the so called ‘parallel resonance’ phenomenon [62] in the meshed power system. Note that relatively small time-step is used for the simulation of the subsynchronous resonance because large frequency deviation occurs after the fault is applied. The reason why a small time step has to be used is that large frequency deviation defines a wide bandwidth in the shifted frequency domain, which in turn requires smaller time steps to respect the Nyquist 97 Chapter 5 EMTP Implementation frequency limit in the shifted frequency domain. This indicates that a variable time step scheme is a future research direction to realize full potential of the SFA method. The time step can be reduced when system states are changing rapidly in order to achieve better accuracy in the SFA simulation. On the other hand, when system dynamics slowdown, large time steps can be used to avoid unnecessarily long computational time while still achieving reasonable accuracy. Generator R=0.0002 X=0.02 Bus 2 1 = 0.0067 R R = 0 0.0186 R = 0.0074 1 0.022 0 R Xi = 0.0739 =0.21 0 X Xi = 0.08 =0.24 0 X Bus C Xiijne Bus 1 Ri =R =0.0014 0 1 X = = 0.03 0 X Infinite Bus Figure 5.24 The Second Benchmark Network for Subsynchronous Resonance Studies 98 EMTP Implementation Chapter 5 V C: -C V C: C > I time (s) Figure 5.25 Generator Termiual Voltages: Phase A (At = 1 ms) V C: C: -C V C C U Ct U C: V V C’) 3 time (s) Figure 5.26 Voltage across Series Capacitor: Phase A (At = 1 ms) 99 Chapter 5 EMTP Implementation 0 •0 CD I I 2.5 time (s) Figure 5.27 Transformer High Side (Bus 2) Voltage: Phase A (At = 1 ms) > 0 ci, 0 Figure 5.28 Infinite Bus (Bus 1) Voltage: Phase A (At = 1 ms) The test results from all the above cases indicate that the dynamic phasor is a generalization of the phasor concept, which can represent the dynamic waveform in power systems, without loss of important information. The SFA method with its implementation in the EMTP environment can integrate the differential equations of the power system in the SFA domain and can produce dynamic phasors for electrical variables, which are visually clear and 100 Chapter 5 EMTP Implementation easy to follow for power engineers. With the SFA method, power engineers may gain better insight into the EMTP simulation results, and we would expect broader applications of the EMTP in power system steady state and dynamics studies, beyond the fast transient simulations. 101 Chapter 6 Conclusions Chapter 6 Conclusions 6.1 Summary of Contributions The goal of this thesis is to extend EMTP functionality for power system dynamic simulation, especially for simulating dynamics with frequency spectra close to the fundamental power frequency. This has been accomplished by developing the Shifted Frequency Analysis (SFA) method, modeling system components with SFA, and accelerating the EMTP simulations for dynamics around 60Hz. A series of contributions made in this thesis are the following. I. The theory of Shifted Frequency Analysis is proposed with the help of Hubert transform and analytic signal concept. Numerical accuracy analysis is performed for the discretetime SFA simulation. II. Linear circuit components, transformer, exponential load and steady-state induction motor are modeled in the shifted frequency domain. III. An efficient SFA-based synchronous machine model for the simulation of slow system dynamics is developed. This model is a general-purpose model that can be used for evaluating the dynamic performance of both the salient-pole and the cylindrical-rotor machines. IV. The SFA method is extended to model the induction machines. This thesis proposes a new phase-domain induction machine model based on the equivalent-reduction (ER) approach. The ER model has a concise discrete-time equivalent circuit that can be directly incorporated into EMTP-type simulators. Based on the ER model, a SFA induction machine model is proposed as a general purpose model capable of simulating both fast transients and slow dynamics. V. An EMTP simulation tool based on the SFA is developed. Simulation results 102 Chapter 6 Conclusions validate that the SFA method is capable to efficiently simulate power system fundamental frequency dynamics. This is the first practical accomplishment to build a unified power system analysis tool based on the EMTP solution. 6.2 Future Research The accomplishment achieved in this thesis will inspire researchers to apply SFA and associated techniques to power system analysis and other potential areas. The success of the applications of SFA may be achieved in different aspects such as: Apply the SFA to the modeling and analysis of renewable energy resources, dispersed generation plants, and industrial power systems. In the foreseeable future, more and more independent-power-producer-owned (IPP-owned) generation units will be connected into the power transmission and distribution system due to a deregulated electricity market. Many of these plants will be cleaner, and many of them will be renewable energy sources such as biomass, solar, wind, geothermal, small hydro, ocean energy, and so on. The IPP interconnection impact studies, which identify the system constraints brought about by the integration of IPPs, and determine network upgrades and remedial action schemes, are therefore critical for the reliable and safe operation of the TPPs and associated power systems. By applying the SFA technique in the IPP modeling, it may lead to accurate and efficient solutions in the IPP impact study, which is hardly achievable with the current phasor tools. II. Develop new component models and controller models, and improve the SFA-based EMTP solver. In addition, the future SFA simulator will adopt more object-oriented design such that each module for new components can simply ‘plug into’ the simulation engine without modifying the core codes. The robustness and efflciency of the core codes will also be improved by using up-to-date sparsity techniques. III. The equivalent-reduction idea can be used for building more efficient models for various types of machines, transformers and power system equipments, as long as they can be mathematically represented by differential equations with mutually coupled state variables. These new models can then be further explored in the shifted frequency domain. 103 Chapter 6 IV. Conclusions Transform the SFA-based simulation tool into parallel programs. The UBC Object Virtual Network Integrator (OVNI) [631 is a real time parallel simulator, which uses PC clusters [64] as hardware and Muti-Area Thevenin Equivalent (MATE) algorithm as solution engine [65]-[69]. MATE partitions the power system into subsystems and solves them in parallel. By changing the real-valued component models to the dynamic phasor models, a new OVNI simulator based on SFA can be implemented. In the future, one may expect this new simulator to serve as a distributed simulation tool for the supervision and control of self-healing power infrastructures. V. Further investigate the basic theory of the SFA method. This will still be an interesting field in future research. There are two limitations in the current SFA method. First, SFA simulation is computationally more expensive than EMTP simulation for the same integration step At. This is due to the higher computational cost of operating with complex numbers in SFA as compared to real numbers in EMTP and also to the computational cost involved in transferring between shifted frequency domain and time domain. Second, the aliasing effect may occur when simulating very fast transients in a system [24]. This may introduce error or distortion in SFA simulation result. Some new theories such as discrete-time analytic signal, Hilbert-Huang transform, and new antialiasing techniques are likely to bring theoretical breakthroughs and may lead to the next generation SFA method, which may be more flexible and accurate than the current SFA. VI. It may be worthwhile to look into transformations to d,q,0-quantities as an alternative to dynamic phasors [70], particularly for cases where the three-phase impedances are balanced, and where the faults are symmetrical three-phase faults. VII. Last but not least, the applications of the SFA method is not limited to the EMTP solution. The Shifted Frequency Analysis can also be introduced into any other circuit simulators, for instances, some state-space-based circuit and/or electric machinery solvers, SPICE, and so on. It is equally suitable to be implemented in the variable time step simulators such as Simulink®. It can also be used in hybrid simulations, and hardware-in-loop simulators. All in all, the SFA is believed to have opened new doors of opportunity for research relating to dynamic simulations, and further theoretical work and applications are justified. 104 Bibliography Bibliography [1] U.S.-Canada Power System Outage Task Force, Final Report on the August 14, 2003 Blackout in the United States and Canada: Causes and Recommendations, Apr. 2004. [2] G. T. Heydt, C. C. Liu, A. G. Phadke, V. Vittal, “Solution for the Crisis in Electric Power Supply,” IEEE Computer Applications in Power, vol. 14, no. 3, pp. 22-30, Jul. 2001. [3] IEEE/CIGRE Joint Task Force on Stability Terms and Definitions, “Definition and Classification of Power System Stability,” IEEE Transactions on Power Systems, vol. 19, no. 3,pp. 1387- 1401,Aug.2004. [4] W. Li, Risk Assessment of Power Systems: Models, Methods, and Applications, New York: IEEE Press and Wiley & Sons, 2005. [5] S. R. Aumuri, L. R. Malone, V. Burtnyk, “Representation of Single-Pole Open Conditions in Stability Studies”, IEEE Transactions on Power Systems, vol. 6, no. 1, pp 9-15, Feb. 1991. [6] V. Ajjarapu, C. Christy, “The Continuation Power Flow: a Tool for Steady State Voltage Stability Analysis”, IEEE Transactions on Power Systems, vol. 7, no. 1, pp 416-423, February 1992. [7] T. Van Cutsem, C.Voumas, Voltage Stability of Electric Power Systems, Kluwer Academic Publishers, 1998. [8] IEEE/PES Power System Stability Subcommittee, Voltage Stability Assessment: Concepts, Practices and Tools, Final Document, Aug. 2002. [9] H. W. Dommel, EMTP Theory Book, 2nd edition, Vancouver, BC: MicroTran Power System Analysis Corporation, 1996. [10] H. W. Dommel, “Digital Computer Solution of Electromagnetic Transients in Single- and Muti-Phase Networks,” IEEE Transactions on Power Apparatus and Systems, vol. 88, no. 2, pp. 734-771, Apr. 1969. 105 Bibliography [11] MicroTran Power System Analysis Corporation, Micro Tran Reference Manual: Transients Analysis Program for Power and Power Electronic Circuits, Vancouver, BC, Canada, 2002. [121 Manitoba HVDC Research Center, User’s Guide of PSCAD, Winnipeg, Manitoba, Canada, 2003. [13] Canada/America EMTP Users Group, ATP Rule Book, 1998. [14] J. Mahseredjian, S. Dennetiere, L. Dube, B. Khodabakhchian, and L. Gerin-Lajoie, “On a New Approach for the Simulation of Transients in Power Systems,” Electric Power Systems Research, vol. 77, no. 11, pp. 15 14-1520, 2007. [15] DIgSILENT GmbH, DlgSlLENTPowerFactory V13 User Manual, 2005. [16] RTDS Technologies, Real Time Digital Simulator User Manual, 2003 [17] C. Henville, R. Folkers, A. Hiebert and R. Wierckx, “Dynamic Simulations Challenge Protection Perfonnance,” Proceedings of the 30th Western Protective Relay Conference, Spokane, WA, Oct. 2003. [18] V. Venkatasubramanian, “Tools for Dynamic Analysis of the General Large Power System Using Time-Varying Phasors,” International Journal ofElectrical Power and Energy Systems, vol. 16, no. 6, pp. 365-376, Dec. 1994. [19] M. Stubbe, A. Bihain, J. Deuse, J. C. Baader, “STAG-A New Unified Software Program for the Study of the Dynamic Behaviour of Electrical Power Systems”, IEEE Transactions on Power Systems, vol. 4, no. 1, pp 129-138, Feb. 1989. [20] S. Henschel, “Analysis of Electromagnetic and Electromechanical Power System Transients with Dynamic Phasors,” Ph. D. dissertation, The University of British Columbia, Vancouver, Canada, 1999. [21] A. M. Stankoviá, T. Aydin, “Analysis of Asymmetrical Faults in Power Systems Using Dynamic Phasors,” IEEE Transactions on Power Systems, vol. 15, no. 3, pp. 1062-1068, Aug. 2000. [22] A. M. Stankoviá, S. R. Sanders, T. Aydin, “Dynamic Phasors in Modeling and Analysis of Unbalanced Polyphase AC Machines,” IEEE Transactions on Energy Conversion, vol. 17, no. l,pp. 107-113, Mar. 2002. 106 Bibliography A. M. Stankovié, B. C. Lesieutre, T. Aydin, “Modeling and Analysis of Single-Phase [23] Induction Machines with Dynamic Phasors,” IEEE Transactions on Power Systems, vol. 14, no. 1, pp. 9-14, Feb. 1999. [24] J. R. Marti, “Shifted Frequency Analysis (SFA) for EMTP Simulation of Fundamental Frequency Power System Dynamics”, Internal Report, Power System Laboratory, The University of British Columbia, Vancouver, Canada, Apr. 2005. 5. L. Hahn, Hubert Transforms in Signal Processing, Norwood, MA: Artech House, [25] 1996. A.V. Oppenheim and A. S. Wilisky, Signal and Systems, Upper Saddle River, NJ: [26] Prentice-Hall, 1996. R. Shintaku and K. Strunz, “Branch Companion Modeling for Diverse Simulation of [27] Electromagnetic and Electromechanical Transients,” IPST’OS, Montreal, Canada, Jun. 2005. P. Zhang, J. R. Marti, H. W. Dommel, “Synchronous Machine Modelling Based on [28] Shifted Frequency Analysis,” IEEE Transactions on Power Systems, vol. 22, no. 3, pp. 11391147, Aug. 2007. P. Zhang, J. R. Marti, H. W. Dommel, “Induction Machine Modelling Based on Shifted [29] Frequency Analysis,” IEEE Transactions on Power Systems, Accepted for Publication. J. R. MartI, J. Lin, “Suppression of Numerical Oscillations in the EMTP,” IEEE [30] Transactions on Power Systems, vol. 4, no. 2, pp. 739-747, May 1989. L. W. Ehrlich, “Complex Matrix Inversion Versus Real,” Communications of the ACM [31] vol. 13, no. 9, pp. 561-562, Sept. 1970. L. Tornheim, “Inversion of a Complex Matrix,” Communications of the ACM, vol. 4, [32] issue 9, pp. 398, Sept. 1961. [33] P. Kundur, Power System Stability and Control, New York: McGraw-Hill, 1994. [34] BPA Transient Stability Program User’s Guide, Bonneville Power Administration, 1987. [35] J. R. MartI and K. W. Louie, “Phase-Domain Synchronous Generator Model Including Saturation Effects,” IEEE Transactions on Power Systems, vol. 12, no. 1, pp. 222—229, Feb. 1997. [36] S. D. Pekarek, 0. Wasynczuk, and H. J. Hegner, “An Efficient and Accurate Model for the Simulation and Analysis of Synchronous Machine/Converter Systems,” IEEE Transactions on Energy Conversion, vol.13, no. 1, pp. 42-48, Mar. 1998. 107 Bibliography [37] L. Wang and J. Jatskevich, “A Voltage-behind-Reactance Synchronous Machine Model for the EMTP-Type Solution”, IEEE Transactions on Power Systems, vol. 21, no. 4, pp 1539—1549, Nov. 2006. [38] J. Machowski, J. Bialek, and J. Bumby, Power System Dynamics and Stability, J. Wiley & Sons, 1997. [39] J. R. Marti, L. R. Linares, J. Calviflo, H. W. Dommel, J. Lin, “OVNI: An Object Approach to Real-Time Power System Simulators, “Proceedings of the 1998 International Conference on Power System Technology, Powercon’98, Beijing, China, August 18-21, 1998, pp. 977-981. [40] F. Moreira, J. R. MartI, “Latency Techniques for Time-Domain Power System Transients Simulation,” IEEE Transactions on Power Systems, vol. 20, no. 1, pp. 246-253, Feb. 2005. [41] V. Brandwajn, “Synchronous Generator Models for the Analysis of Electromagnetic Transients”, Ph.D. Thesis, The University of British Columbia, Vancouver, BC, Canada, 1977. [42] 5. D. Pekarek, E. A. Walters, “An Accurate Method of Neglecting the Dynamic Saliency of Synchronous Machines in Power Electronic Based Systems,” IEEE Transactions on Energy Conversion, vol. 14, no.4, pp. l1771184, Dec. 1999. [43] P. Krause, 0. Wasynczuk, and S. Sudhoff, Analysis of Electric Machinery and Drive Systems, 2nd ed., New York: IEEE Press, 2002. [44] IEEE Recommended Practice for Power Systems Analysis, IEEE Standard 399-1997, Sept. 1997. [45] J. Lin, J. Mahseredjian and S. Lefebvre, “Improvement of Synchronous Machine Saturation Simulation in the EMTP,” Proceedings of IEEE Region 10 Conference, vol 5, pp. 127-132, Oct. 1993. [46] R. Hung, H. W. Dommel, “Synchronous Machine Models for Simulation of Induction Motor Transients,” IEEE Transactions on Power Systems, vol. 11, no. 2, pp 833—838, May 1996. [47] IEEE Tutorial on Digital Simulation of Electrical Transient Phenomena, IEEE Catalog No. 81 EHO173-5-PWR, 1981. 108 Bibliography [48] G. J. Rogers, D. Shirmohammadi, “Induction Machine Modelling for Electromagnetic Transient Program,” IEEE Transactions on Energy Conversion, vol. 2, no. 4, pp 622—628, Dec. 1987. [49] R. H. Park, “Two-Reaction Theory of Synchronous Machines—Generalized Methods of Analysis—Part I,” AIEE Trans., vol. 48, pp. 716—727, July 1929. [50] H. K. Lauw and W. S. Meyer, “Universal Machine Modeling for the Presentation of Rotating Electric Machinery in an Electromagnetic Transients Program,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-lOl, no. 6, pp 1342—1351, Jun. 1982. [51] J. R. MartI and T. 0. Meyers, “Phase-Domain Induction Motor Model for Power System Simulators,” in Proc. IEEE Western Canada Conf Exhibition, pp. 276-282, 1995. [52] G. Kron, Tensor Analysis of Networks, New York: J. Wiley & Sons; London: Chapman &Hall, 1939. [53] P. Zhang, “Symmetrical Induction Machine Simulation Based on the Voltage-behind- Reactance Model”, Project Report, Power System Laboratory, The University of British Columbia, Vancouver, Canada, Dec. 2005. [54] D. Chaniotis and M. A. Pai, “Model Reduction in Power Systems Using Krylov Subspace Methods,” IEEE Transactions on Power Systems, vol. 20, no. 2, pp 888—894, May 2005. [55] J. Roychowdhury, “Reduced-Order Modeling of Time-Varying Systems”, IEEE Transactions on Circuits and Systems—Il: Analog and Digital Signal Processing, vol. 46, no. 10, pp. 1273—1288, Oct. 1999. [56] R. W. Freund, “Krylov-subspace Methods for Reduced-order Modeling in Circuit Simulation” Journal of Computational and Applied Mathematics, vol. 123, pp. 395—421, 2000. [57] A. C. Antoulas, D. C. Sorensen, and S. Gugercin, “A Survey of Model Reduction Methods for Large-Scale Systems,” Contemporary Mathematics, vol. 280, pp 193—219, 2001. [58] J. C. Fraga, “Tunable Difference Equations for the Time Domain Simulation of Power System Operating States”, Ph.D. Thesis, The University of British Columbia, Vancouver, BC, Canada, 2003. [59] W. H. Kersting, “Radial Distribution Test Feeders,” IEEE Transactions On Power Systems, vol. 6, no. 3, pp. 975-985, Aug. 1991. 109 Bibliography [60] IEEE Subsyncbronous Resonance Task Force, “First Benchmark Model for Computer Simulation of Subsynchronous Resonance,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-96, no. 5, pp. 1565-1572, September/October 1977. [611 H. W. Dommel, Case Studies for Electromagnetic Transients, 3rd edition, Vancouver, BC, Mar. 2004. [62] IEEE Subsynchronous Resonance Working Group, “Second Benchmark Model for Computer Simulation of Subsynchronous Resonance,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-104, no. 5, pp. 1057-1066, May 1985. [63] J. R. MartI, L. R. Linares, J. A. Hollman, F. A. Moreira, “OVNI: Integrated software/Hardware Solution for Real-time Simulation of Large Power Systems,” in Proceedings of the PSCCO2, Sevilla, Spain, June, 2002 [64] J. A. Hollman, J. R. MartI, “Real Time Network Simulation with PC-Clusters,” IEEE Transactions on Power Systems, vol. 18, no. 2, pp. 563-569, May 2003. [65] M. Armstrong, J. R. MartI, L. R. Linares, and P. Kundur, “Multilevel MATE for Efficient Simultaneous Solution of Control Systems and Nonlinearities in the OVNI Simulator,” IEEE Transactions on Power Systems, vol. 21,110. 3, pp. 1250-1259, Aug. 2006. [66] P. Zhang, J. R. MartI and H. W. Dommel, “Network Partitioning for Real-time Power System Simulation,” IPST’05, Montreal, Canada, Jun. 2005. [67] J. R. MartI, L. R. Linares, “Real-time EMTP-Based Transients Simulation,” IEEE Transactions on Power Systems, vol. 9, no. 3, pp. 1309-13 17, Aug. 1994. [68] T. De Rybel, J. Hoilman, J. R. MartI, “OVNI-NET: a Flexible Cluster Interconnect for the New OVNI Real-Time Simulator,” 15th Power Systems Computation Conference, Liege, Belgium, Aug. 22-26, 2005. [69] F. A. Moreira, J. R. MartI, L. C. Zanetta, Jr., L. R. Linares, “Multirate Simulations With Simultaneous-Solution Using Direct Integration Methods in a Partitioned Network Environment,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 53, no. 12, pp. 2765-2778, Dec. 2006. [70] Note from H. W. Dommel, ‘Dynamic Phasors Versus dqo-Quantities,’ Oct. 25, 2008. 110 Appendix A Machine Parameters Appendix A Machine Parameters A.1 Synchronous Machine Parameters TABLE A-I SmcHRoNous MACHINE PARAM1TERs Parameters Hydro Turbine Generator Steam Turbine Generator 3 r 0.00234 2 0.1478(2 0.5911(2 1.0467(2 0.0005 (2 0.2523(2 0.01675 (2 0.1267(2 0.01736(2 0.1970(2 Js 7 3.51x10 2 325 MVA 0.00243 2 0.1538(2 1.457(2 1.457(2 0.00075 (2 0.1145(2 0.00144(2 0.6578 (2 0.0068 1 (2 0.07602(2 0.0108(2 0.06577(2 Ps 4 6.58x10 2 835 MVA 20 KV 26 KV 64 2 X, X Xd rfd X 1 r Xlkql 12 r 1 X rij X 1 J Rating Line-to-line voltage Poles - - 111 Appendix A Machine Parameters A.2 Induction Machine Parameters TABLE A-Il INDUCTION MACHINE PARAMi’rERs Parameters r x’s XM Xir rr J TB Line-to-line voltage Poles 3-hp Induction Machine 0.435 2 0.754 2 26.132 0.754 2 0.816Q 2 0.089kgm 11.9Nm 2250-hp Induction Machine 0.0292 0.2262 13.04 0.226 0.022Q 63.87 kgm 2 Nm 3 8.9x10 220 V 2.3 kV 4 4 112 Appendix B Voltage behind Reactance Induction Machine Model Appendix B Voltage behind Reactance Induction Machine Model This appendix faithfully reproduces reference [53], a project report for UBC course EECE 549. The only difference here is that all equation numbers and figure numbers are added with a ‘B.’. This report implemented the Voltage behind Reactance (VBR) model for the symmetrical induction machine simulations. The discrete-time VBR model has the similar form as the EMTP-type equivalent circuits. The simulation results validate the accuracy and efficiency of this new induction machine model. B. 1 Voltage Behind Reactance Model of Induction Machine The equations of the induction machine can be expressed in the arbitrary reference frame as rslqs + Vqs VdS 1 ds CO)Lqs = + P)Lq (B.1) +pAd$ (B.2) + o = rriqr + (o — OJr o ( 0)r ,. 0 o = rrior + pA, = rridr (B.3) , 0 pA — )adr + Pqr (B.4) )Aqr 1dr 2 +P (B.5) (B.6) where each variable and parameter have been converted.to stator side, and co = 0 for stationary reference frame, co = co,. for rotor reference frame and co = coe for 113 Appendix B Voltage behind Reactance Induction Machine Model synchronous reference frame. Flux linkages equations are expressed as 05 A = Lislqs + Amq = Lis ‘qs + LM qr) (B.7) = SIdS 1 L + Ad = L,SidS + LM (‘di + ld) (B.8) = 05 i 5 L, (‘qs + (B.9) Aqr LirIqr + Amq (B. 10) Adr = Lirl& + A (B. 11) 0 i 1 =L ‘0r (B.12) Rearrange and manipulate (B.7) (B.12), we can obtain - A,q = L [qs Amd1M[1ds + (B. 13) +.z’,J (B.14) where =(+J’ LL Substitute (B.13), (B.14) into (B.7), (B.8), the q-axis and d-axis flux linkages can be rewritten as 5 = L”1q Aq 5 (B. 15) + A =L”idS +A’ (B.16) where L”—L is -‘-L”M — (B.17) (B.18) Therefore (B.1) and (B.2) can be written as Vqs = + o(L”idS + A) + P( ”lqs 1 + (B.19) A:) 114 Appendix B = rSidS Voltage behind Reactance Induction Machine Model — (D(L”qs (B.20) + p(L”i + + By manipulating (B.1)-(B.18), we can obtain LZIqs (B.21) _Cor)2] )—(Co = +LZidS _ar)+(0_C0r)2qr] 1 (B.22) =r[ Substitute (B.21), (B.22) into (B.19), (B.20), the voltage behind reactance form of the voltage equations can be expressed = rsiqs Vq VdS = rSidS + — OJL”idS + P(L”lqs) + (B.23) (DL”lqs + P(L”lqs) + v (B.24) where + ‘‘ = lqg — — —(co — v=—coA:+(;—Adr)+i+(co—cor), (B.25) (B.26) Applying the inverse transformation to (B.23), (B.24), we can get the phase domain equations V (t) = (t) + p[L”i, (t)] + v(t) (B.27) where v(t) = [K(o)]’ v 0 +±L” L Is3M LbCS = —— L” _L 3 3 +—L 1 L L,+LZ 115 Appendix B Voltage behind Reactance Induction Machine Model B. 2 Discrete Time VBR Model Based on the trapezoidal rule, equation (B.27) can be discretized and rearranged as abcs 7 ‘ (t) = + --L 5 )iabcs(t) + vbC$(t) (B.28) + ek (t) where Ch (t) — = 5 L } abcs — At) + (t VbCS — At) — (t Vabcs (t — At) (B.29) Rewrite (B.25) and (B.26), we can get the matrix form ofthesubtransient voltages = fa 11 12 a [a21 a22 (t)J [v (t)1 + Ea13 0 ]Ldr (t)] [ ji (t)1 0 (B.30) a3 ][i (t)] where a =‘‘I—i L LL r 1 12 a 21 = a 12 a 22 = a a 11 Lrr L 23 = a a 13 a 13 =Wr(t)jj By rearranging and manipulating (B.4),( B.5), (B.10) and (B.11), the rotor flux linkage equations can be obtained (t) pH [A (t)J E11 21 [b (t)1 22 ][A.. b 12 b (t)] + rb! 3 [0 0 jIqs (t) (B.31) 23 ][i (t) b where b‘ =--1-—1 , 1 L 21 b = 12 —b 12 b = —[a(t) 22 b = 11 b — 0)r 13 =--LZ b (t)] 23 b = 13 b Discretizing (B.3 1) by using the Trapezoidal rule, then r[— — E °‘)1 (t)J 2 11 b At 2b1 (t)At L’dr + 2 — At 11 b 11 b —b 12 (t)At1 2 + At 2 — At 22 J [b b , (t — At)At 2 —b 12 (t)At1’ Ebi 0 At 3 2—b 22 At] [ 0 23 b At] Jj••1qS 116 — At) J[2ar — At) (t)1 + (t — At) (t)J [i (t — At) — 1b2 (t At)At1E (t 2 + At 22 b (t 1q5 (B.32) Appendix B Voltage behind Reactance Induction Machine Model Substitute (B.32) into (B.30), we can get the discrete time equations [ j r1 =K(t)[ vd(t) 1 (i’) (t) 2 j+K ldS(t) where K (t) = Hii K 2 (t) — — 2— At 11 b 12 a 21 [a 22 a 2 (t)At J[— b raii 12 a 11 b ir 2 At ]L b21 (t)At [ + Eaii , 2 [a 22 a — At 1 2 b, j 22 a 21 (t)At b 12 a — ][— — 2 12 b — Ebi 1 (t)Atl A 3 t 22 b At] L 0 1 23 At] b 0 + E’u13 [ (t)Atlr 23 a 0 At)&1r (t JL’dr (t 2 + At 11 b 1b2 (t 1b2 2 At 22 J[b b 2 + At 22 b , (t At)At 2 — — — — 2 b, — (t)At1’ P13 At 22 b 2 — At] [ 0 0 jqs 0’ — — — At) At) 23 (t — At) b At][idS Applying the inverse transformation to (B.33), we can transfer the qd0 variables to the phase domain VbCS (B.34) (t) = K (t)Iabcs (t) + ChS (t) where K(t) = [K 5 (0(t))]-’ [K (t) [ ,EK ehS(t)=[KS(0(t))] L 01 j[K. (0(t))] 0’) Finally, by substituting (B.34) into (B.28), the equivalent circuit for the stator voltage equations can be written as ‘abcs (t) = G eqVabcs (t) — h(t) (B.35) where Geq K(t)] 5 + [r + h(t) =Geq[Ch(t) + ehS(t)] 117 Appendix B Voltage behind Reactance Induction Machine Model B. 3 Free Acceleration Simulation of a 3-hp Induction Machine The first test case is to simulate the start-up transients of a 3-hp induction machine. The parameters of the machine are P=4 V 1 220V rr = 0.816c2 J r 0.4352 X, 0.754 2 XM= 26.13 2 X r 1 O.754 0.089 kg.m . 2 The simulation results viewed from different reference frames are shown in Fig. B. 1- B.3. (1) Figure B.1: Stationary Reference Frame Because the initially position of the stationary qd0 reference frame is zero, the f = fqs. Therefore the Iqs and Vqs are the same as i and v, respectively. The rotor variables are varying at 60 HZ after referred to the stationary reference frame. (2) Figure B.2: Rotor Reference Frame At the very beginning, when the speed of the machine is slow, the variables in the rotor reference frame are similar to those in the stationary reference frame. Then the stator variables will vary at slip frequency when the rotor speed is increasing. After the speed reach the synchronous speed, all variables in rotor reference frames becomes constant. (3) Figure B.3: Synchronous Reference Frame At stall the variables in the synchronous reference are varying at 60 HZ. The stator and rotor variables will become constant when the rotor speed reaches synchronous speed. Note that the zero position of synchronous reference frame that we chose is zero, therefore the magnitude of Vqs is identical to this of v and Vds is zero. B. 4 Dynamic Performance of Induction Machine during Mechanical Torque Changes In this case, the machine originally operates in the steady state. At t = 0.05s, a constant mechanical torque of 12 N.m 2 is applied. Then the mechanical torque is reversed to -12N.m 2 at t = 0.5s. The simulation time is ls. Figure B.4 illustrates the simulation results with the VBR model. 118 Voltage behind Reactance Induction Machine Model Appendix B 0’ V .0 E 200 100 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 >0• S 0 S 0 -100 Figure B.1 Free Acceleration Characteristics in Stationary Reference Frame 119 1 ______________________________________________________ Appendix B Voltage behind Reactance Induction Machine Model I i.J.. I .: 0 06 07 08 09 0.5 0.6 0.7 0.8 0.9 0.6 0.7 0.8 0.9 I I I 07 08 09 1 OC - I 30C 0 0.1 0.2 0.3 0.4 I \\45 ,nn 2O0OjOi2Oi3 04 05 06 Figure B.2 Free Acceleration Characteristics in Rotor Reference Frame 120 Appendix B Voltage behind Reactance Induction Machine Model 100 I 0 0.1 2C 0.2 I I I 0.3 0.4 0.5 0.6 I I I 0.7 0.8 0.9 I C’ > -20 -4C I I I I I I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 I J 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 05 06 07 08 09 -I (VI 0 I 50 0 -50 100O 300 1 I Figure B.3 0I2 0I3 0I4 0i5 08 0I7 08 I I I I I I I 09 I Free Acceleration Characteristics in Synchronous Reference Frame 121 Appendix B Voltage behind Reactance Induction Machine Model 21 22 22.12.2 2.3 2.4 2.5 2.6 2.7 2.8 I I 2.6 2.7 2.8 2.9 2.9 3 j 500C -500C 2 2.1 I I I I 2.2 2.3 2.4 2.5 220C I 3 I iGoc 500C —--- — 2 2.1 I I 2.2 2.3 2.4 I I I I I 2.5 2.6 2.7 2.8 2.9 Figure BA Dynamic Performance of a 3-hps Induction Machine during Step Changes in Load Torque 122 Appendix C Test Case Data Appendix C Test Case Data %*********** ************ *** *** ********** ************ ******** *************** % Test Case 1 % AnElectromagnetic Transient Program with Dynamic Phasor Solution % % % = % % = Version 1.0 = Last revised: 24 Feb.2007 = Peng Zhang = Copyright 2008 © Peng Zhang = % nBus=2; % P1 Line data nLine=1; % Three phase line LineFromBus(1)=1; LineToBus(1)2; I Ri ohm I Xl ohm I 1R21 X2lC2l.... Cl muF LineParameters =[0.0868455 0.0298305 0.0288883 0.2025449 0.0847210 0.0719161 0.00274 -0.0007 0.0298305 0.0887966 0.0298305 0.0847210 0.1961452 0.0847210 -0.0007 -0.00034; 0.00296 -0.00071 0.0288883 0.0298305 0.0868455 0.0719161 0.0847210 0.2025449 -0.00034 -0.00071 0.00274 1*10.56 % Load data nLoadl;% Three phase load LoadBus(1)=1; rLoad(1,1)=(12.47*1000/sqrt(3))P2/(2000*1000); % phase A resistance xLoad(1,1)=(12.47* 1000/sqrt(3)) ’2/(0.328684105 17886306346562595373367*2000*1000); % phase A reactance t rLoad(1,2)=(12.47*1000/sqrt(3)y2/(2000* 1000); % phase B resistance xLoad(1,2)=(12.47*1000/sqrt(3)y2/(0.32868410517886306346562595373367*2000*1000); % phase B reactance 2/(2000* 1000);% phase C resistance t rLoad(1,3)(12.47*1000/sqrt(3)) xLoad(1,3)(12.47* 1000/sqrt(3)y2/(0.32868410517886306346562595373367*2000*1000);% phase C reactance % Voltage source data nVSourcel; VSourceBus(1)2; Vabcs(1,1) Vabcs(1,2) 10182;%Vamp = 10182 * exp(-jay * 2 * p / 3); 123 Appendix C Vabcs(1,3) 10182 Test Case Data * exp(jay* 2 * pi/3); T10= 3; LineTripped = 1; return; %************************************************************************** % Test Case 2 % An Electromagnetic Transient Program with Dynamic Phasor Solution % % % = % % = Version 1.0 = Last revised: 24 Feb.2007 = Peng Zhang = Copyright 2008 © Peng Zhang = % simT = 8; % total simulation time ja’sqrt(-l); ws2*pi*60; nBusl2; % PT Line data nLine=1 1; % Three phase line LineFromBus(l)12; LineToBus(1)=1; LineFromBus(2)1; LineToBus(2)4; LineFromBus(3)4; LineToBus(3)5; LineFromBus(4)=1; LineToBus(4)=2; LineFromBus(5)2; LineToBus(5)=3; LineFromBus(6)=l; LineToBus(6)6; LineFromBus(7)6; LineToBus(7)8; LineFromBus(8)8; LineToBus(8)9; LineFromBus(9)8; LineToBus(9)l0; LineFromBus(10)=6; LineToBus(10)1 1; LineFromBus(1 1)6; 124 Appendix C Test Case Data LineToBus(1 1)7; I I Ri LineParameters(1:3,1:9)=[0.3465 0.0 1.0179 0.0 0.3375 0.0 0.0 0.0 0.0 0.3414 1R21X21C21.... ClmuF 6.2998/(2*pi*60) _1.2958/(2*pi*60) 1.2595/(2*pi*60); 0.5017 0.4236 0.5017 1.0478 0.3849 1.2958/(2*pi*60) 0.4236 0.3849 1.0348 1.2595/(2*pi*60) 1.2595/(2*pi*60) 5.9597/(2*pi*60) 1.2595/(2*pi*60); 5.6386/(2*pi*60)] 0.0 1.3569 0.4591 0.4591 4.6658/(2*pi*60) 1.3294 0.0 0.4591 1.3471 0.4591 0.8999/(2*pi*60) 4.7097/(2*pi*60) 0.8999/(2*pi*60); 0.4591 1.3569 .0.8999/(2*pi*60) .0.8999/(2*pi*60) 4.6658/(2*pi*60)] LineParameters(1:3,10:18) =[1.3238 0.0 0.0 I Xl 0.0 0.0 LineParameters(1:3,19:27) =[1.3238 0.0 1.3238 0.4591 1.3294 0.0 0.0 0.0 LineParameters(1:3,28:36) [0.7526 0.0 0.8999I(2*pi*60); 1.3569 0.4591 0.4591 4.66581(2*pi*60) 0.8999/(2*pi*60) 0.8999/(2*pi*60); 0.0 0.0 0.8999/(2*pi*60) 4.7097/(2*pi*60) 0.8999/(2*pi*60); 0.4591 1.3471 0.4591 0.8999/(2*pi*60) 1.3238 0.4591 0.4591 1.3569 0.8999/(2*pi*60) 0.8999I(2*pi*60) 4.6658/(2*pi*60)] 0.0 1.1814 5.6990/(2*pi*60) 1.0817/(2*pi*60) 1.6905/(2*pi*60); 0.4236 0.5017 0.0 0.7475 0.0 0.4236 1.1983 0.3849 1.0817/(2*pi*60) 0.0 0.0 0.5017 0.3849 1.2112 1.6905I(2*pi*60) 0.6588/(2*pi*60) 0.7436 5.1795/(2*pi*60) 0.6588/(2*pi*60); 5.4246/(2*pi*60)] LineParameters(1:3,37:45) =[O.3807 0.0000 0.0000 0.6922 0.0000 0.0000 0.0001 0.0000 0.0000; 0.0000 0.3807 0.0000 0.0000 0.6922 0.0000 0.0000 0.000 1 0.0000 0.0000 0.0000 0.3807 0.0000 0.0000 0.6922 0.0000 0.0000 0.0001 LineParameters(1:3,46:54) [ 0.3465 0.0 0.0 0.0 0.3375 0.0 0.0 0.0 LineParameters(1:3,55:63) =[1.3238 0.0 1.0179 0.5017 0.3414 0.4236 6.2998/(2*pi*60) 1.2595/(2*pi*60) 1.2595I(2*pi*60); 0.5017 0.4236 1.0478 0.3849 1.2595I(2*pi*60) 5.9597/(2*pi*60) 1.2595/(2*pi*60); 1.0348 1.2595/(2*pi*60) 1.2595/(2*pi*60) 0.3849 0.4591 5.6386/(2*pi*60)] 0.0 0.0 1.3238 0.0 0.4591 0.4591 0.4591 0.8999/(2*pi*60) 4.70971(2*pi*60) 0.8999/(2*pi*60); 0.0 0.0 0.4591 0.4591 1.3471 0.8999I(2*pi*60) 0.8999I(2*pi*60) 4.7097/(2*pi*60)] 1.3294 0.4591 4.6658/(2*pi*60) 0.8999/(2*pi*60) 0.8999/(2*pi*60); 1.3569 LineParameters(1 :3,64:72) =[1.3292 0.0000 0.0000 1.3475 0.0000 0.0000 4.5 193/(2*pi*60) 0.0000 1.3292 0.0000 0.0000 1.3475 0.0000 0.0000 0.0000 0.0000 1.3292 0.0000 0.0000 1.3475 0.0000 0.0000 4.5193/(2*pi*60) 0.0000 0.0000; 0.0000; 4.5193/(2*pi*60)] \ LineParameters(1:3,73:81)=[1.3425 0.0000 0.0000 0.5124 0.0000 0.0000 88.9912/(2*pi*60) 0.0000 0.0000; 0.0000 1.3425 0.0000 0.0000 0.5124 0.0000 0.0000 88.9912/(2*pi*60) 0.0000; 0.0000 0.0000 1.3425 0.0000 0.0000 0.5124 0.0000 0.0000 88.9912/(2*pi*60)] Lineparameters(1:3,82:90) =[ 0.3465 0.0 0.0 0.0 6.2998/(2*pi*60) .1.9958/(2*pi*60) ..1.2595/(2*pi*60); 0.0 1.0179 0.5017 0.4236 0.3375 0.0 0.5017 1.0478 0.3849 1.9958/(2*pi*60) 5.9597/(2*pi*60) 0.7417/(2*pi*60); 0.4236 0.3849 1.0348 -1 .2595/(2*pi*60) .0.7417/(2*pi*60) -0.0143 96.8897/(2*pi*60) 0.0 LineParameters(1:3,91:99) =[0.7982 0.0 0.3414 0.0 0.4463 0.0328 0.0000 0.0000; 0.0 0.7891 0.0 0.0328 0.4041 0.0328 0.0000 96.88971(2*pi*60) 0.0000; 0.0 0.0 0.7982 -0.0143 0.0328 0.4463 0.0000 0.0000 % Load data nLoad=8;% Three phase load LoadBus(1)’3; LoadBus(2)=4; LoadBus(3)5; LoadBus(4)=10; LoadBus(5)6; LoadBus(6)=7; LoadBus(7)1 1; LoadBus(8)=9; rLoad—zeros(8,3); 125 96.8897/(2*pi*60)] 5.6386/(2*pi*60)] Test Case Data Appendix C rLoad(:,:)1e6; xLoad=zeros(8,3); xLoad(:,:)=1e6; rLoad(5,1)=(4. 16* 1000/sqrt(3))s2/(385* 1000); % phase A resistance xLoad(5,1)=(4. 16*1000/sqrt(3))2/(220* 1000); % phase A reactance rLoad(5,2)=(4. 16* l000/sqrt(3))/s2/(385* 1000); % phase B resistance xLoad(5,2)=(4.16*1000/sqrt(3))f2/(220*1000); % phase B reactance rLoad(5,3)=(4.16* l000/sqrt(3))f2/(385* 1000); % phase C resistance xLoad(5,3)=(4. 16*100O/sqrt(3))l2/(220*1000); % phase C reactance % Induction machine data IndMacBus=1 1; Ttorq=3; % torque applied I IndMacParameters = S I VI IP rs J I I Lls % Voltage source data nVSourcel; VSourceBus(1)=12; Vabcs(1,1) = Vabcs(1,2) = LM I iT I Llr [500 /(1.341e-3) 4600 4 11.06 0.262*4 1 .206*4/ws 54.02*4/ws 0.187*4 1 .206*4Iws]; (4.16*1000/sqrt(3))*sqrt(2) ; % Vamp (4.16*l000Isqrt(3))*sqrt(2) * exp(-jay Vabcs(1,3) = (4.16*l000Isqrt(3))*sqrt(2) * exp(jay * 2 * pj / 3); * 2 * pi /3); return; 126
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Shifted frequency analysis for EMTP simulation of power...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Shifted frequency analysis for EMTP simulation of power system dynamics Zhang, Peng 2009
pdf
Page Metadata
Item Metadata
Title | Shifted frequency analysis for EMTP simulation of power system dynamics |
Creator |
Zhang, Peng |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | Electromagnetic Transients Program (EMTP) simulators are being widely used in power system dynamics studies. However, their capability in real time simulation of power systems is compromised due to the small time step required resulting in slow simulation speeds. This thesis proposes a Shifted Frequency Analysis (SFA) theory to accelerate EMTP solutions for simulation of power system operational dynamics. A main advantage of the SFA is that it allows the use of large time steps in the EMTP solution environment to accurately simulate dynamic frequencies within a band centered around the fundamental frequency. The thesis presents a new synchronous machine model based on the SFA theory, which uses dynamic phasor variables rather than instantaneous time domain variables. Apart from using complex numbers, discrete-time SFA synchronous machine models have the same form as the standard EMTP models. Dynamic phasors provide envelopes of the time domain waveforms and can be accurately transformed back to instantaneous time values. When the frequency spectra of the signals are close to the fundamental power frequency, the SFA model allows the use of large time steps without sacrificing accuracy. Speedups of more than fifty times over the traditional EMTP synchronous machine model were obtained for a case of mechanical torque step changes. This thesis also extends the SFA method to model induction machines in the EMTP. By analyzing the relationship between rotor and stator physical variables, a phase-coordinate model with lower number of equations is first derived. Based on this, a SFA model is proposed as a general purpose model capable of simulating both fast transients and slow dynamics in induction machines. Case study results show that the SFA model is in excess of seventy times faster than the phase-coordinate EMTP model when simulating the slow dynamics. In order to realize the advantage of SFA models in the context of the simulation of the complete electrical network, a dynamic-phasor-based EMTP simulation tool has been developed. |
Extent | 3788957 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067234 |
URI | http://hdl.handle.net/2429/7727 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2009_spring_zhang_peng.pdf [ 3.61MB ]
- Metadata
- JSON: 24-1.0067234.json
- JSON-LD: 24-1.0067234-ld.json
- RDF/XML (Pretty): 24-1.0067234-rdf.xml
- RDF/JSON: 24-1.0067234-rdf.json
- Turtle: 24-1.0067234-turtle.txt
- N-Triples: 24-1.0067234-rdf-ntriples.txt
- Original Record: 24-1.0067234-source.json
- Full Text
- 24-1.0067234-fulltext.txt
- Citation
- 24-1.0067234.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0067234/manifest