Główna
Table of contents
Table of contents
Florin Iov
0 /
0
How much do you like this book?
What’s the quality of the file?
Download the book for quality assessment
What’s the quality of the downloaded files?
Rok:
2004
Język:
english
Plik:
PDF, 2,32 MB
Ściągnij (pdf, 2,32 MB)
 Otwórz w przeglądarce
 Checking other formats...
 Konwertować w EPUB
 Konwertować w FB2
 Konwertować w MOBI
 Konwertować w TXT
 Konwertować w RTF
 Przekonwertowany plik może różnić się od oryginału. Jeśli to możliwe, lepiej pobierz plik w oryginalnym formacie.
 Najpierw zaloguj się na swoje konto

Czy potrzebna jest pomoc? Zobacz instrukcje dotyczące wysyłania książki do Kindle
Plik zostanie dostarczony na Twój email w ciągu 15 minut.
Plik zostanie dostarczony do Twojego Kindle w ciągu 15 minut.
Please note: you need to verify every book you want to send to your Kindle. Check your mailbox for the verification email from Amazon Kindle.
Please note: you need to verify every book you want to send to your Kindle. Check your mailbox for the verification email from Amazon Kindle.
0 comments
Możesz zostawić recenzję książki i podzielić się swoimi doświadczeniami. Inni czytelnicy będą zainteresowani Twoją opinią na temat przeczytanych książek. Niezależnie od tego, czy książka ci się podoba, czy nie, jeśli powiesz im szczerze i szczegółowo, ludzie będą mogli znaleźć dla siebie nowe książki, które ich zainteresują.
1

2

Wind Turbine Blockset in Matlab/Simulink General Overview and Description of the Models Florin Iov, Anca Daniela Hansen, Poul Sørensen, Frede Blaabjerg Aalborg University March 2004 22 Wind Turbine Blockset in Matlab Simulink Abstract. This report presents a new developed Matlab/Simulink Toolbox for wind turbine applications. This toolbox has been developed during the research project “Simulation Platform to model, optimize and design wind turbines” and it has been used as a general developer tool for other three simulation tools: Saber, DIgSILENT, HAWC. The report provides first a quick overview over Matlab issues and then explains the structure of the developed toolbox. The attention in the report is mainly drawn to the description of the most important mathematical models, which have been developed in the Toolbox. Then, some simulation results using the developed models are shown. Finally, some general conclusions regarding this new developed Toolbox as well as some directions for future work are made. ISBN 8789179463 Institute of Energy Technology Aalborg University Copyright © by Aalborg University Print: UNI.PRINT Aalborg University March 2004 2 33 Wind Turbine Blockset in Matlab Simulink Table of contents Table of contents ........................................................................................................................3 Preface ........................................................................................................................................5 Chapter 1 Introduction ............................................................................................................6 Chapter 2 Matlab/Simulink® – An overview .........................................................................7 Chapter 3 Wind Turbine Blockset...........................................................................................9 3.1 General structure ........................................................................................................9 3.1.1 Mech; anical Components Library .....................................................................10 3.1.2 Electrical Machinery Library............................................................................11 3.1.3 Power Converters Library ................................................................................12 3.1.4 Common Bocks Library ...................................................................................13 3.1.5 Transformations Library...................................................................................14 3.1.6 Measurements Library......................................................................................14 3.1.7 Control Library .................................................................................................15 3.2 Simulation speed in Matlab/Simulink ......................................................................16 3.3 Summary.....................................................................................................................7 Chapter 4 Description of the models.....................................................................................21 4.1 Mechanical component models ................................................................................21 4.1.1 Wind model ......................................................................................................21 4.1.2 Wind turbine rotor model .................................................................................23 4.1.3 Drive trains models...........................................................................................25 4.1.3.1 Threemass model ........................................................................................25 4.1.3.2 Twomass model ..........................................................................................26 4.1.3.3 Onemass model ...........................................................................................27 4.2 Electrical machinery models ....................................................................................28 4.2.1 Induction machine models................................................................................28 4.2.1.1 ABC/abc model of induction machine .........................................................28 4.2.1.2 Model of induction machine in dqo arbitrary reference frame.....................33 Fluxes as state variables ...........................................................................................34 Currents as state variables ........................................................................................35 Electromagnetic torque.............................................................................................36 4.2.1.3 Reduced order model....................................................................................37 4.2.1.4 Steadystate model........................................................................................38 4.2.1.5 Modelling deepbar effect ............................................................................41 4.2.1.6 Modelling saturation.....................................................................................45 4.2.1.7 Modelling iron losses ...................................................................................46 4.2.2 Synchronous Machine Model...........................................................................48 4.2.2.1 Dynamic Equations of Synchronous Machine .............................................48 Untransformed model ...............................................................................................48 Transformed model with quadraturephase stator windings ....................................49 Commutator model ...................................................................................................50 Electromagnetic torque.............................................................................................52 4.2.2.2 Simulink model ............................................................................................52 4.2.3 Permanentmagnet synchronous machine ........................................................53 3 44 Wind Turbine Blockset in Matlab Simulink 4.3 Power converter models........................................................................................... 56 4.3.1 AC Controllers (Softstarters).......................................................................... 56 4.3.1.1 Star connected 3phase load......................................................................... 58 4.3.1.2 Delta connected 3phase load ...................................................................... 60 4.3.1.3 Branchdelta connected load........................................................................ 61 4.3.1.4 Simulink implementation of the softstarter ................................................ 62 4.3.1.5 RMS model for softstarter .......................................................................... 62 4.3.2 Voltage Source Converters .............................................................................. 64 4.3.2.1 Star connected generator.............................................................................. 65 4.3.2.2 Delta connected generator............................................................................ 66 4.3.2.3 Simulink implementation of the voltage source converter .......................... 67 4.3.3 DClink Circuit ................................................................................................ 68 4.4 Modulation strategies............................................................................................... 68 4.4.1 PWM sinusoidal modulator using sinetriangle with third harmonic .............. 69 4.4.2 Space Vector Modulation................................................................................. 70 4.5 Grid component models ........................................................................................... 72 4.5.1 Grid model ....................................................................................................... 72 4.5.2 Circuit breaker model....................................................................................... 73 4.5.3 Capacitor bank ................................................................................................. 75 4.5.4 Cable model ..................................................................................................... 76 4.5.5 Transformer model........................................................................................... 78 4.5.5.1 3phase 2winding transformer .................................................................... 78 4.5.5.2 Determination of ideal transformer matrix .................................................. 79 4.5.5.3 Parameter estimation for linear transformer model ..................................... 80 Determination of parameters for lossless threephase transformer.......................... 80 Inclusion of loss terms ............................................................................................. 81 4.5.5.4 Dynamic equations of 3 phase 2 windings transformer............................... 82 4.5.6 Phase Locked Loop.......................................................................................... 83 4.6 Control blocks – Control of PQ for DFIG ............................................................... 84 Statorflux oriented control.............................................................................................. 84 4.7 Summary .................................................................................................................. 87 Chapter 5 Simulation results................................................................................................. 88 5.1 Connection of a noload transformer to the grid...................................................... 88 5.2 Disconnection of a noload transformer................................................................... 89 5.3 Startup sequence of a fixedspeed wind turbine ..................................................... 90 5.4 Variable speed/pitch wind turbine with DFIG......................................................... 92 5.5 Summary .................................................................................................................. 21 Conclusions and Future work .................................................................................................. 94 References................................................................................................................................ 97 Appendix A An example of a “C” SFunction ..................................................................... 101 Appendix B Other developed models ................................................................................... 105 Harmonic voltage source and grid angle detection............................................................ 105 Threephase diode bridge rectifier ..................................................................................... 107 Synchronous PWM modulator........................................................................................... 107 ABC to Magnitude and phase transformation.................................................................... 108 Average wind calculator .................................................................................................... 108 4 55 Wind Turbine Blockset in Matlab Simulink Preface This report describes the Wind Turbine Blockset developed in Matlab/Simulink during the project “A Simulation Platform to Model, Optimize and Design Wind Turbines”. The project has been funded by the Danish Energy Agency contract #ENS1363/010013, and it was carried out in cooperation between Aalborg University and RISØ National Laboratory. 5 66 Wind Turbine Blockset in Matlab Simulink Chapter 1 Introduction In the last years Matlab/Simulink has become the most used software for modelling and simulation of dynamic systems. Wind turbine systems are an example of such dynamic systems, containing subsystems with different ranges of the time constants: wind, turbine, generator, power electronics, transformer and grid. Matlab/Simulink provides a powerful graphical interface for building and verifying new mathematical models as well as new control strategies for the wind turbines systems. Then, using a dSPACE prototyper these new control strategies can be easily implemented and tested in a HardwareIntheLoop structure. This report presents a new developed Matlab/Simulink Toolbox for wind turbine applications. This toolbox has been developed during the research project “Simulation Platform to model, optimize and design wind turbines”. Matlab/Simulink has been used in this Simulation Platform as a general developer tool for other three simulation tools, namely: Saber, DIgSILENT and HAWC. The report provides first a quick overview over Matlab issues and then explains the structure of the new developed toolbox. A special attention has been on increasing the simulation speed for all the models developed in this new Toolbox. Some aspects regarding the implementation methods in Simulink for a dynamic model are presented. The attention in the report is mainly drawn to the description of the most important mathematical models, which have been developed in the Toolbox. Some simulation results using the developed models are shown. During the work carried out in developing of this Toolbox a large amount of references including books, papers and other sources have been used. In the report only the most relevant references are mentioned for each particular topic. However, the reference list presented in report includes almost all references in alphabetical order used during the work. 6 77 Wind Turbine Blockset in Matlab Simulink Chapter 2 Matlab/Simulink® – An overview MATLAB® is a highperformance language for technical computing. It integrates computation, visualization, and programming in an easytouse environment where problems and solutions are expressed in a familiar mathematical notation [78]. Typical uses include: • Math and computation; • Algorithm development; • Data acquisition; • Modelling, simulation, and prototyping; • Data analysis, exploration, and visualization; • Scientific and engineering graphics; • Application development, including graphical user interface building. The MATLAB system consists of five main parts [78]: • Development Environment. This is the set of tools and facilities that help you use MATLAB functions and files. Many of these tools are graphical user interfaces. It includes the MATLAB desktop and Command Window, a command history, an editor and debugger, and browsers for viewing help, the workspace, files, and the search path. • MATLAB Mathematical Function Library. This is a vast collection of computational algorithms ranging from elementary functions like sum, sine, cosine, and complex arithmetic, to more sophisticated functions like matrix inverse, matrix eigenvalues, Bessel functions, and fast Fourier transforms. • MATLAB Language. This is a highlevel matrix/array language with control flow statements, functions, data structures, input/output, and objectoriented programming features. It allows both "programming in the small" to rapidly create quick and dirty throwaway programs, and "programming in the large" to create complete large and complex application programs. • Graphics. MATLAB has extensive facilities for displaying vectors and matrices as graphs, as well as annotating and printing these graphs. It includes highlevel functions for twodimensional and threedimensional data visualization, image processing, animation, and presentation graphics. It also includes lowlevel functions that allow you to fully customize the appearance of graphics as well as to build complete graphical user interfaces on your MATLAB applications. • MATLAB Application Program Interface (API). This is a library that allows to write C and Fortran programs that interact with MATLAB. It includes facilities for calling routines from MATLAB (dynamic linking), calling MATLAB as a computational engine, and for reading and writing MATfiles. 7 88 Wind Turbine Blockset in Matlab Simulink In the last few years, Simulink® has become the most widely used software package in academia and industry for modelling and simulating dynamic systems [78]. Simulink® is a graphical software package for modelling, simulating, and analyzing dynamic systems and it is based on Matlab®. It supports linear and nonlinear systems, modelled in continuous time, sampled time, or a hybrid of the two. Systems can also be multirate, i.e., have different parts that are sampled or updated at different rates. For modelling, Simulink provides a graphical user interface (GUI) for building models as block diagrams, using clickanddrag mouse operations. With this interface, the desired dynamic systems can be easily built. Simulink includes a comprehensive block library of sinks, sources, linear and nonlinear components, and connectors. Using SFunctions it is also possible to customize and create userdefined blocks. Models are hierarchical, so the models can be built using both topdown and bottomup approaches. The system can be viewed at a high level, then doubleclick blocks to go down through the levels to see increasing levels of model detail. This approach provides insight into how a model is organized and how its parts interact. After defining a model, it can be simulated, using a choice of integration methods, either from the Simulink menus or by entering commands in the MATLAB Command Window. The menus are particularly convenient for interactive work, while the commandline approach is very useful for running a batch of simulations (e.g. Monte Carlo simulations or sweeping a parameter across a range of values). Using scopes and other display blocks, the simulation results can be analyzed while the simulation is running. In addition, the parameters can be changed during the simulation for "what if" exploration. The simulation results can be put in the MATLAB workspace for postprocessing and visualization. Model analysis tools include linearization and trimming tools, which can be accessed from the MATLAB command line, plus the many tools in MATLAB and its application toolboxes. Since MATLAB and Simulink are integrated, the models can be simulated, analyzed, and revisited in either environment at any point. An Sfunction is a computer language description of a Simulink block. Sfunctions can be written in MATLAB, C, C++, Ada, or Fortran. C, C++, Ada, and Fortran Sfunctions are compiled as MEXfiles using the mex utility [78]. As with other MEXfiles, they are dynamically linked into MATLAB when needed. Sfunctions use a special calling syntax that enables the user to interact with Simulink's equation solvers. This interaction is very similar to the interaction that takes place between the solvers and builtin Simulink blocks. The form of an Sfunction is general and can accommodate continuous, discrete, and hybrid systems. Sfunctions allow the user to add new models to the Simulink builtin models. These new blocks can be written in MATLAB®, C, C++, Fortran, or Ada. An algorithm or model can be implemented in an Sfunction and then customize a user interface by using masking. 8 99 Wind Turbine Blockset in Matlab Simulink Chapter 3 Wind Turbine Blockset In this chapter the structure of a new Matlab/Simulink Toolbox for wind turbine applications is presented. The content of the main libraries from this toolbox is briefly shown. Then, some considerations about the simulation speed for the developed models are done. 3.1 General structure A new Matlab/Simulink Toolbox for wind turbine applications has been developed during the project. In order to analyze the dynamic and/or steady state behaviour of a wind turbine, the basic components of a wind turbine have been modelled and structured in seven libraries: Mechanical Components, Electrical Machinery, Power Converters, Common Models, Transformations, Measurements and Control as shown in Fig. 3.1. Fig. 3.1. Structure of the new Matlab/Simulink Toolbox for wind turbine simulations These models are built based on Simulink blocks as well as using C SFunctions. Some of the features of this developed Toolbox can be summarized as follows: • All the developed models use basically only Simulink blocks. No particular Toolbox is necessary to simulate a wind turbine concept. So, the user should have installed on the PC only Matlab and Simulink. However, some models uses particular blocks from other Toolboxes from Matlab e.g. Power Spectra Density calculation block; • It uses the matrix support in order to minimize the number of blocks and connection lines; • All models which involves a great number of differential equations (e.g. electrical machines, drivetrains and transformer) are also available as ‘C’ SFunctions for highspeed simulations; • In order to be able to use different drivetrain models the equation of motion is not included in the electrical machine models; 9 1010 Wind Turbine Blockset in Matlab Simulink • Special models have been developed for induction generators (abc/abc models), which can be used in some particular fault analysis e.g. unsymmetrical, unbalanced faults with unbalanced loads. In this case cannot be used the standard dqo model. Moreover, this model is used to simulate the softstarter fed induction machines, which are used in fixedspeed wind turbines. A short description of the libraries illustrated in Fig. 3.1 is presented in the following paragraphs. 3.1.1 Mechanical Components Library The first library Mechanical Components contains: wind model, aerodynamic models of the wind turbine rotor, different types of drive train models as shown in Fig. 3.2. Fig. 3.2. Mechanical Components Library. The wind model has been developed at RISØ National Laboratory based on the Kaimal spectra. The wind speed is calculated as an averaging of the fixedpoint wind speed over the whole rotor, and it takes the tower shadow and the rotational turbulences into account [71]. Since one of the main components in this model is the normally distributed white noise generator some investigations have been done in order to obtain the same wind time series in all considered simulation tools. It has been found that the builtin white noise generator from different simulation tools uses different algorithms and therefore different wind time series are obtained. A new normally distributed white noise generator has been implemented using a ‘C’ SFunction based on the Ziggurat Algorithm developed by G. Marsaglia [45]. 10 1111 Wind Turbine Blockset in Matlab Simulink The aerodynamic models of the wind turbine rotor are based on the torque coefficient CQ or power coefficient CP lookup table. As mentioned before, the equation of motion is not included in the machine model. Different types of drive train models are developed e.g. onemass model, twomass model of the drive train with torsional torques. Since the last model involves 4 differential equations an alternative implementation using C SFunction is also available. 3.1.2 Electrical Machinery Library The Electrical Machinery Library contains models with different level of detailing for electrical machines used currently in wind turbine systems as shown in Fig. 3.3. The standard dq models as well as the new models for squirrelcage induction machine, wound rotor induction machine, salientpoles synchronous machine and permanent magnet synchronous machines have been included in this library. Fig. 3.3. Electric Machinery library based on Simulink blocks. Currently, all commercial simulation software uses the dq or dqoreference frames in modelling the electrical machines [83][85]. However, some particular operating modes e.g. softstarterfed induction machines requires an abc/abc model [67], [69], [73], [74]. Therefore, the models of the electrical machines are written in statevariable form in the dqoreference frame as well as in the abc/abc natural reference frame. Some special features, as deepbar effect in the rotor of the induction machines, are also included in models. Further extensions of the models e.g. temperature, saturation, etc. can easily be added. 11 1212 Wind Turbine Blockset in Matlab Simulink Another important feature of the abc/abc machine models is that the dynamic equations are written in perphase quantities. Therefore, the desired windingconnection can be obtained e.g. star, delta or branchdelta as in the case of softstartersfed induction machines. Besides the complete dynamic models of induction machine some reduced order models have been developed and implemented in the advanced aero elastic tool HAWC [35]. In the first stage the focus has been in modelling of the induction machine and some special models including deepbar effect, reduced order models (neglecting the stator transients) and steady state models have been developed. Next step is to extend the library with the corresponding models for the synchronous machine both field winding and permanent magnet. Since all these electrical machines models involves a relatively big number of derivates for the state variables two methods have been chosen in order to implement the models in Simulink: using Simulink blocks and C SFunction. In §3.2 it will be shown how the implementation method of the dynamic equations affects the simulation speed. In Fig. 3.4 is shown the electrical machine library based on C SFunctions. An example of a C SFunction for the dqmodel of the induction machine is shown in Appendix A. Fig. 3.4. Electric machinery C SFunction based models. However, contrary to C SFunction the implementation method using Simulink blocks is much more open and accessible for the user to “see” the model and eventually to make modifications on it. 3.1.3 Power Converters Library The third library contains models for power converters based on switching functions. At the moment the following models are available: 3phase diode bridge rectifier, voltage source converter (VCS), softstarters and different modulation strategies for power converters, as shown in Fig. 3.5. 12 1313 Wind Turbine Blockset in Matlab Simulink Fig. 3.5. Content of the Power Converters Library. Using two VSC blocks connected via a DClink circuit a backtoback converter topology can be obtained. The connection type of the softstarter can be selected from the mask interface, e.g. star connection, delta or branch delta connection. Two types of modulation strategy for the Voltage Source Converter are available, namely sinusoidal PWM and Space Vector PWM both average and switching models. The switching models include minimum pulse width, deadtime compensation and pulse dropping. 3.1.4 Common Bocks Library This library contains models related with the grid connection of the wind turbines as shown in Fig. 3.6. Some relevant models from this library are: • Voltage source with a harmonic content according with the EN6100022 standard is included in this library; • Transformer model, which takes the core geometry as well as the iron losses into account, has been developed. The model is implemented both using Simulink blocks and C SFunction; • Distribution line model is based on the dynamic equations of the short line without distortion, which is a relatively simple model with good results in the case of a distribution line; • Switch model (onoff) of the circuit breaker, which takes different opening time moments for poles into account, has been developed; • Capacitor bank model. 13 1414 Wind Turbine Blockset in Matlab Simulink Fig. 3.6. Common Blocks Library. 3.1.5 Transformations Library As in a dynamic model of a wind turbine there are several transformations for the signals, a special library, which contains the main transformations, has been implemented as shown in Fig. 3.7. The basic transformations are from abc to dqo arbitrary reference and the opposite one, as well as from a dq signal to the complex representation of it. In order to take the connection type into account, some transformations from a threephase system in phase quantities to a threephase one but in linetoline quantities have been implemented. In order to find the magnitude and the phase for a threephase system the corresponding block can be used. Fig. 3.7. Content of the Transformations Library. Some of these models e.g. ABC >Magnitude&Phase uses special blocks from DSP Blockset. 3.1.6 Measurements Library This library contains some special blocks as: calculation of the period for a sinusoidal variable, calculation of the grid angle using a Phase Locked Loop system, calculation of the 14 1515 Wind Turbine Blockset in Matlab Simulink average wind for a given time interval, different modes of active and reactive power calculation (monophase, dqreference frame, complex), as shown in Fig. 3.8. Fig. 3.8. Measurements Library. Some of these models e.g. Average Wind uses special blocks from DSP Blockset. 3.1.7 Control Library This library contains models related with the control of the wind turbine system. Since the main component of a control scheme is the PI Controller, a model with an anti windup structure has been implemented in Simulink. The Maximum Power Point Tracker (MPPT) block is based on a lookup table obtained from the wind turbine characteristics. A simplified control for active and reactive for a doubly fed induction generator is developed. This control algorithm can be used with a reduced order model of the generator Fig. 3.9. Content of the Control Blocks Library. This library will be extended with other types of controllers as for example the active stall controller and capacitor bank controller used in fixedspeed wind turbines. 15 1616 Wind Turbine Blockset in Matlab Simulink 3.2 Simulation speed in Matlab/Simulink Using the available blocks from Matlab/Simulink three methods have been considered to investigate the simulation speed. Each method implements the same set of dynamic equations for the induction machine in different ways. The general structure for these models is shown in Fig. 3.10. Fig. 3.10. Simulink model of induction machine in Simulink. The first two models, Model A and Model B, have the same structure as shown in Fig. 3.11. Fig. 3.11. Implementation of induction machine model in dqrotor reference frame used in Model A and Model B. The only difference between Model A and Model B consists in the implementation of the statespace equations with currents as statevariables. Model A implements the dynamic equations using Simulink blocks as Sum, Product, Gain, GoTo, From and Integrator as shown in Fig. 3.12. The same equations are implemented in Model B by using few blocks, e.g. Function and Integrator blocks. The third model  Model C – shows in Fig. 3.12 is based on an Sfunction written in “C” language according with the Simulink format for this type of function. The coordinate transformations abc/dq and dq/abc are also included in this function as well as the calculation of the reference frame angle. So, the Sfunction called “scim” implements five differential equations. 16 1717 Wind Turbine Blockset in Matlab Simulink Fig. 3.12. Simulink implementation of dynamic equations for Model A. Fig. 3.13. Simulink implementation of dynamic equations for Model B. Fig. 3.14. Simulink implementation of dynamic equations for Model C. 17 1818 Wind Turbine Blockset in Matlab Simulink Using Simulink Performance Tools Toolbox® (Profiler) the considered models have been compared regarding the simulation speed. The Simulink Profiler collects performance data while simulating a model and generates a simulation report. This report is useful to find out how much time Simulink spends on executing each block from a model and hence where to focus the model optimization efforts. The simulation mechanism in Simulink involves several callbacks to some specific Simulink Functions at each time step. Fig. 3.15 shows the pseudocode, which summarizes the execution of a dynamic system in Simulink. Fig. 3.15. Pseudocode for execution of a dynamic system model in Simulink According to this conceptual model, Simulink executes a dynamic model by invoking the functions presented in Table 3.1 zero, one, or more times, depending on the function and the model. The Profiler measures the time required to execute each invocation of these functions and generates a report at the end of the model that details how much time was spent in each function. Table 3.1 Function invoked by Simulink during simulation. Purpose Level Simulate the model. This toplevel function invokes the other functions required to simulate the model. The time spent in this System sim function is the total time required to simulate the model. Set up the model for simulation System ModelInitialize Execute the model by invoking the output, update, integrate, etc., functions for each block at each time step from the start to the end of System ModelExecute simulation. Compute the outputs of a block at the current time step. Block Output Update a block’s state at the current time step. Block Update Compute a block’s continuous states by integrating the state Block Integrate derivatives at the current time step. Compute a block’s output at a minor time step. Block MinorOutput Compute a block’s state derivatives at a minor time step. Block MinorDeriv Block MinorZeroCrossings Compute a block’s zero crossing values at a minor time step. Free memory and perform any other endofsimulation cleanup. System ModelTerminate Function 18 1919 Wind Turbine Blockset in Matlab Simulink The considered models have been tested using Simulink Profiler for a simulation time of 0.3 s with a simulation step of 1 ms. Ode5 fixed step solver has been used. Table 3.2. Simulink Profiler Summary Report for the considered models. Parameter Model A Model B Model C Total recorded time 13.35 sec 6.98 sec 3.02 sec Number of block methods 71 40 15 Number of internal methods 9 9 9 The parameters illustrated in Table 3.2 are: • Total recorded time – total time required to simulate the model using Simulink Profiler; • Number of block methods – total number of invocations of blocklevel functions (e.g. Output()); • Number of internal methods – total number of invocations of systemlevel functions (e.g. ModelExecute). The number of blocks used in a model has a big influence on the total simulation time as shown in Table 2.2. Model A, which uses much more blocks than Model B, has the biggest recorded time even if the number of Integrator blocks is the same. Fig. 3.16 shows the simulation time in percent for functions invoked to simulate the entire model. The functions for the constitutive blocks from the models have been omitted. Function Integrate is the most time consuming procedure during the simulation around 7080 % from the total simulation time. So, the total simulation time for a Simulink model is affected first by the number of blocks used in the model and then by the number of Integrator blocks. Fig. 3.16. Total time spent executing all invocations of specific functions for the considered models as a percentage of the total simulation time for Model A. 19 2020 Wind Turbine Blockset in Matlab Simulink A comparative study regarding execution time for a C Sfunction block and one Integrator block, which compute the rotor speed in this case, for the Model C is presented in Fig. 3.17. Fig. 3.17. Average time for execution of a C SFunction and for an Integrator block. It has been considered the average time (selftime) required to execute Update, Output and MinorDeriv functions for these blocks as a percent from the sum of the selftime for these functions for the Integrator block. It is obvious that the total time required for all invocations of functions related to the Sfunction block is smaller than the corresponding time for one Integrator block even if this Sfunction implements five differential equations. Therefore, using a Sfunction which implements all the differential equations from the considered dynamic system, the simulation speed can be increased with at least a factor of two. Moreover, the numerical stability increases using an SFunction. So, in modelling of large dynamic systems, which involves a big number of Integrator blocks it is desirable to use C SFunction concept in order to increase the simulation speed and improve the numerical stability. 3.3 Summary In this chapter the structure of the new Matlab/Simulink Toolbox for wind turbine applications developed during the Simulation Platform Project is presented. The content of the main libraries is briefly presented and then some aspects regarding the simulation speed for the developed models are shown. 20 2121 Wind Turbine Blockset in Matlab Simulink Chapter 4 Description of the models In this chapter a detailed description of the models from Wind Turbine Blockset is presented. The mathematical descriptions for the models as well as their Simulink implementation are shown. 4.1 Mechanical component models 4.1.1 Wind model The wind model has been developed at RISØ National Laboratory based on the Kaimal spectra. The wind speed is calculated as an average value of the fixedpoint wind speed over the whole rotor, and it takes the tower shadow and the rotational turbulences into account [71]. A main component in this model is the normally distributed white noise generator. Therefore, in order to obtain the same wind time series in all considered simulation tools used in the Simulation Platform some investigations have been done. It has been found that the builtin white noise generator from different simulation tools uses a different algorithm and thus a different wind time series is obtained. A new normally distributed white noise generator has been implemented using a ‘C’ SFunction based on the Ziggurat Algorithm developed by G. Marsaglia [45]. Currently, the wind model is available in two basic versions. The first one uses the normally distributed white noise generator from Matlab/Simulink, while the second one is based on the new developed normally distributed white noise generator. The general structure of these Simulink models is shown in Fig. 4.1. Fig. 4.1. Simulink implementation of the wind model. 21 2222 Wind Turbine Blockset in Matlab Simulink The parameters defined in the block’s mask are: rotor diameter of the wind turbine, average wind speed, length scale, turbulence intensity and sample time as shown in Fig. 4.2. Fig. 4.2. Mask interface for wind model. In order to validate the new white noise generator model some comparisons have been done. A wind time series for 3600 sec with 0.05 sec sample time, an average wind speed of 10 m/sec and 12% turbulence intensity has been generated for both models as shown in Fig. 4.3. Fig. 4.3. Wind time series for the considered models (ZA – Ziggurat Algorithm, SB – builtin Simulink Block). The 20 binwidth histograms for these wind time series are presented in Fig. 4.4. 22 2323 Wind Turbine Blockset in Matlab Simulink Fig. 4.4. Wind histogram for the considered models (ZA – Ziggurat Algorithm, SB – builtin Simulink Block). Finally, the power spectra density for both wind time series has been calculated as shown in Fig. 4.5. Fig. 4.5. Power spectra density for the considered models (ZA – Ziggurat Algorithm, SB – builtin Simulink Block). Even if the wind time series for the considered models are not identically as shown in Fig. 4.3 the power spectra density is approximately the same, which is an important issue. 4.1.2 Wind turbine rotor model The aerodynamic model of the wind turbine rotor is based on the torque coefficient CQ or the power coefficient CP lookup table. The torque coefficient CQ is used to determine the aerodynamic torque directly by using: Twt = 0.5πρR 3 v∞2 CQ (4.1.1) 23 2424 Wind Turbine Blockset in Matlab Simulink where: ρ is the air density, R is the blade radius, v∝ is the wind speed and CQ is the torque coefficient. Alternatively the aerodynamic torque can be determined using the power coefficient CP based on: Twt = 0.5πρR 2 v3∞ C P (4.1.2) Important to underline that both coefficients CP and CQ can be function of the tip speed ratio λ for passivestall wind turbines or function of tip speed ratio λ and pitch angle θ for active stall and variable pitch/speed wind turbines. The Simulink model for the variable pitch wind turbine rotor is presented in Fig. 4.6. Fig. 4.6. Simulink model of the wind turbine rotor. The parameters from the mask interface are: blade radius, air density, cutin and cutout wind speeds, as shown in Fig. 4.7. Fig. 4.7. Mask interface for the wind model in Simulink. Using a reduced lookup table in respect with the variation of the torque coefficient/power coefficient with the pitch angle ( −10o ≤ θ ≤ 10o ) the model for the activestall wind turbine is obtained. Imposing a zero value for the pitch angle θ the model for passive stall wind turbine can be obtained. 24 2525 Wind Turbine Blockset in Matlab Simulink 4.1.3 Drive trains models Starting from the threemass model for a wind turbine drive train the twomass model and the onemass model are derived. 4.1.3.1 Threemass model The equivalent model of a wind turbine drive train is presented in Fig. 4.8. The masses correspond to a large mass of the wind turbine rotor, masses for the gearbox wheels and a mass for generator respectively. Fig. 4.8. Twomass model of a wind turbine drive train. Taking into account the stiffness and the damping factors for both shafts the dynamic equations can be written as [66]: dΩ wtr + D wtr Ω wtr + k swtr (θ wtr − θ1 ) dt dΩ1 + D wtr Ω1 + k swtr (θ1 − θwtr ) T1 = J1 dt dΩ 2 dΩ 2 + Dgen + k sgen (θ2 − θgen ) T2 = J 2 dt dt dΩgen −Tgen = J gen + Dgen Ωgen + k sgen (θgen − θ2 ) dt dθwtr dθ1 = Ω wtr , = Ω1 dt dt dθgen dθ2 = Ω2 , = Ω gen dt dt Twtr = J wtr (4.1.3) where: Twtr – wind turbine torque; Jwtr – wind turbine moment of inertia; Ωwtr – wind turbine mechanical speed; kswtr – spring constant indicating the torsional stiffness of the shaft on wind turbine part;Tgen – generator torque; Jgen – generator moment of inertia; Ωgen – generator mechanical speed; ksgen – spring constant indicating the torsional stiffness of the shaft on generator part; T1 – torque that goes in the gearbox; T2 – torque out from the gearbox; 1 T2 = ⋅ T1 , Kgear – the gearbox ratio; Ω 2 = k gear ⋅ Ω1 k gear 25 2626 Wind Turbine Blockset in Matlab Simulink 4.1.3.2 Twomass model The threemass model can be reduced to a twomass model by considering an equivalent system with an equivalent stiffness and damping factor. The moment of inertia for the shafts and the gearbox wheels can be neglected because they are small compared with the moment of inertia of the wind turbine or generator. Therefore the resultant model is essentially a two mass model connected by a flexible shaft. Only the gearbox ratio has influence on the new equivalent system. The dynamic equations can be written in two points: on the wind turbine side with the influence of generator component through the gearbox and on the generator side respectively. The equivalent system on the generator side is shown in Fig. 4.9. Fig. 4.9. Equivalent diagram of the wind turbine drive train on the generator side. The dynamic equations of the drivetrain written on the generator side are: T ' wtr = J ' wtr dΩ' wtr + D'e (Ω' wtr − Ω gen ) + k 'se (θ' wrt − θgen ) dt dθ' wtr = Ω ' wtr dt dΩ gen −Tgen = J gen + D'e (Ω gen − Ω ' wtr ) + k 'se (θgen − θ' wtr ) dt dθgen = Ωgen dt (4.1.4) where: the equivalent stiffness is given by: 1 1 1 = + ' k wtr k se k gen 2 k gear (4.1.5) and the equivalent moment of inertia for the rotor is: J ' wtr = 1 k 2 ⋅ J wtr gear The Simulink implementation of (4.1.4) is shown in Fig. 4.10. 26 (4.1.6) 2727 Wind Turbine Blockset in Matlab Simulink Fig. 4.10. Simulink implementation of a twomass model for the wind turbine drive train. The parameters from the mask interface are: moment of inertia for machine and wind turbine rotor, equivalent stiffness and damping coefficients for the shafts, gearbox ratio and initial conditions for the state variables (speeds and angles for the low speed and high speed shafts), as shown in Fig. 4.11. Fig. 4.11. Mask interface for the twomass model of the drive train. 4.1.3.3 Onemass model When the stiffness and the damping factor are neglected further reduction of the drive train model can be obtained. The onemass model is obtained in this case, which is described by: ' Tgen − Twtr = J ech dΩ gen dt (4.1.7) 27 2828 Wind Turbine Blockset in Matlab Simulink J wtr and the equivalent wind turbine 2 k gear where: the equivalent moment of inertia is J ech = J gen + ' torque is Twtr = Twtr 2 k gear The Simulnik implementation of (4.1.7) is presented in Fig. 4.12. Fig. 4.12. Onemass Simulink model for the wind turbine drive train. 4.2 Electrical machinery models In this paragraph the dynamic equations and the Simulink implementation for the generators used in wind turbine applications are shown. 4.2.1 Induction machine models In this paragraph the dynamic equations of induction machine in the natural reference frame abc as well as in the dqo arbitrary reference frame are presented [34]. These equations are written using both fluxes and currents as statevariables. 4.2.1.1 ABC/abc model of induction machine The basic equations of the machine (assuming sinusoidal mmf and neglecting saturation and losses in the core) by using the diagram shown in Fig. 4.13 are [9], [23], [34], [54]: sA stator phase A axis iA rotor phase a axis ωr ra vA θr ia va vb ib vC rb iC sC ic rc vc vB iB sB Fig. 4.13. Simplified diagram of induction machine with rotor and stator windings. 28 2929 Wind Turbine Blockset in Matlab Simulink [ V ] = [ R ] ⋅ [ I] + [ Ψ ] = [ L] ⋅ [ I] d [Ψ ] dt (4.2.1) (4.2.2) where: • [ V ] = [ vA v B vC va v b vc ]  are the voltages applied to each stator and rotor t phases; • [ I ] = [i A i B i C i a i b i c ]  are the currents in each stator and rotor phases; • [ Ψ ] = [ ϕA t ϕB ϕC ϕa ϕb ϕc ]  are the fluxes linked with each stator and rotor t phases; • [R] is the resistance matrix and [L] is the inductance matrix. Substituting (4.2.2) in (4.2.1) and arrange the equation in statevariable form results in: d [ I] d [ L] −1 −1 = [ L ] ⋅ − [ R ] − ⋅ [ I] + [ L] ⋅ [ V ] dt dt (4.2.3) Most terms of the inductance matrix [L] are functions of the rotor position θr. The derivative of the inductance matrix from (4.2.3) can be written as: d [ L ] d [ L] dθr d [ L ] = ⋅ = ⋅ ωr dt dθr dt dθr (4.2.4) where ωr is the rotor speed. Grouping (4.2.3) and (4.2.4) results in: d [ I] d [ L] −1 −1 = [ L ] ⋅ − [ R ] − ωr ⋅ ⋅ [ I] + [ L] ⋅ [ V ] dt dθr (4.2.5) which is the state space form of the dynamic equations for the induction machine. The inductance matrix is defined as: −0.5M s −0.5M s M sr f1 M sr f 2 M sr f 3 Ls −0.5M −0.5M s M sr f 3 Ls M sr f1 M sr f 2 s −0.5M −0.5M L M sr f 2 M sr f 3 M sr f1 [ L] = M f s M f s M sf −0.5M r −0.5M r Lr sr 3 sr 2 sr 1 M sr f 2 −0.5M r M sr f1 M sr f 3 −0.5M r Lr M sr f 2 M sr f1 −0.5M r −0.5M r L r M sr f 3 where: Ls = L σs + M s  stator selfinductance; L r = Lσr + 1.5M r L σs , Lσr Ms , M r M sr (4.2.6)  rotor selfinductance;  stator and rotor leakage inductance;  stator and rotor mutual inductance;  statorrotor mutual inductance. 29 3030 Wind Turbine Blockset in Matlab Simulink The coefficients f1, f2, f3 are: 2π 2π f1 = cos θr , f 2 = cos θr + , f 3 = cos θr − 3 3 (4.2.7) All these inductances can easily be measured in a wound rotor induction machine. For squirrelcage induction machines, these inductances can normally be estimated from noload and locked rotor tests. Referring the rotor parameters to stator the mutual inductances become equals: M s = M r = M sr (4.2.8) If it is taken into account that the machines windings forms an balanced/symmetrical system and: i A + i B + i C = 0 and i a + i b + i c = 0 (4.2.9) the inductance matrix [ L ] can be rewritten as: L's 0 0 M sr f1 M sr f 2 M sr f 3 L's 0 M sr f 3 M sr f1 M sr f 2 0 0s 0 L's M sr f 2 M sr f 3 M sr f1 L = [ ] ' 0 0 M sr f1 M sr f 3 M sr f 2 L r M f M f M f 0 L' r 0 sr 1 sr 3 sr 2 0 0 L'r M sr f 3 M sr f 2 M sr f1 (4.2.10) where: L's = Lσs + 1.5M s and L' r = Lσr + 1.5M s . From (4.2.10) the results are: 0 0 0 g1 g 2 g 3 0 0 0 g g g 3 1 2 0s 0 0 g 2 g 3 g1 d [ L] = − M sr dθr g1 g 3 g 2 0 0 0 g 2 g1 g 3 0 0 0 g 3 g 2 g1 0 0 0 and 2π 2π g1 = sin θr , g 2 = sin θr + , g 3 = sin θr − 3 3 (4.2.11) (4.2.12) The main problem in implementing (4.2.5) is to calculate the inverse of the inductance matrix at each simulation step [23]. A solution is to use modern software tools, which have these capabilities (e.g. Matlab). However, an analytical expression for this matrix is useful since the coefficients depend on the electrical angle θr. One approach is to apply a Cholesky decomposition since the inductance matrix is “symmetric positive definite”[23]. Starting from the general form (4.2.10) the inverse of the inductance matrix can be written as: 30 3131 Wind Turbine Blockset in Matlab Simulink b11 b 21 b −1 [ L] = [ B] = b31 41 b51 b61 b12 b13 b14 b15 b16 b 23 b 24 b 25 b 26 b33 b34 b35 b36 b 43 b 44 b 45 b 46 b53 b54 b55 b56 b62 b 63 b 64 b65 b 66 b 22 b32 b 42 b52 (4.2.13) Using a symbolic mathematical tool e.g. Matlab, Maple or Mathcad the coefficients from (4.2.13) can be determined easily. First, some coefficients are calculated as: K1 = Ls L r M sr2 3 4 K2 = 9 K1 − 4 K1 − K3 = 3 4 − 9 K1 − 4 1 M sr2 K4 = 9 K1 − 4 − (4.2.14) and then all other elements of the matrix B: b11 = b 22 = b33 = K2 Ls b 44 = b55 = b66 = K2 Lr b12 = b13 = b 21 = b 23 = b31 = b32 = K3 Ls b 45 = b 46 = b54 = b56 = b64 = b65 = K3 Lr (4.2.15) b14 = b 25 = b36 = b 41 = b52 = b63 = K 4 f1 b15 = b 26 = b34 = b 43 = b51 = b 62 = K 4 f 2 b16 = b 24 = b35 = b 42 = b53 = b61 = K 4 f 3 The electromagnetic torque can be obtained starting from the voltage equations using currents as state variables d [ V ] = [ R ][ I] + ([ L][ I]) = [ R ][ I] + dt d [ L] d [ I] [ I] + [ L] dt dt (4.2.16) and multiplying it with the transpose of the currents matrix: [ I] [ V ] = [ I] [ R ][ I] + [ I] [ Ψ ] = [ I] [ R ][ I] + [ I] t Alternatively using t t d dt t t d [ L] dt [ I] + [ I] [ L] t d [ I] dt (4.2.17) d [ L ] d [ L ] dθr d [ L ] = ⋅ = ⋅ ωr , where ωr is the electrical speed of the rotor dt dθr dt dθr results in: [ I] [ V ] = [ I] [ R ][ I] + ωr [ I] t t t d [ L] d [ I] t [ I] + [ I] [ L] dθr dt (4.2.18) The terms from (4.2.18) have the following meaning: 31 3232 Wind Turbine Blockset in Matlab Simulink • Pi = [ I] [ V ]  instantaneous power; • Pcopper = [ I ] [ R ][ I ]  copper losses in the machine windings; • d [ I]  magnetic power stored in the machine (due to the variation in dt time of the magnetic energy); • Pm = ωr [ I ] t t Pmag = [ I ] [ L ] t t d [ L] [ I]  mechanical power. dθr The electromagnetic torque is then: Te = Pm P t d [ L] = p m = p [ I] [ I] Ωr ωr dθr (4.2.19) where: p is the number of pole pairs and Ωr is the mechanical speed of the rotor. Substituting (4.2.11) in (4.2.19) the electromagnetic torque as a function of currents is: Te = − pM sr {( i A i a + i Bi b + i Ci c ) sin θ + ( i A i b + i Bi c + i Ci a ) sin ( θ + 2π / 3) + ( i A i c + i Bi a + i Ci b ) sin ( θ − 2π / 3)} (4.2.20) Based on (4.2.5) and (4.2.20) the dynamic model of the induction machine can easily be implemented in Matlab/Simulink, as shown in Fig. 4.14. Fig. 4.14. Simulink implementation of the ABC/abc model for a squirrelcage induction machine. In Fig. 4.15 the mask interface for this model is shown. The input parameters in this interface are based on the standard data sheet: resistance and leakage inductance for the stator and rotor 32 3333 Wind Turbine Blockset in Matlab Simulink windings, magnetizing inductance, number of pole pair and initial conditions for the statevariables (currents and electrical angle of the rotor). Moreover, the operating mode, motor or generator, as well as the connection type for the stator windings, star or delta, can be selected. Fig. 4.15. Mask interface for the abc/abc model of the squirrelcage induction machine. A similar model and mask interface are available for wound rotor induction machine 4.2.1.2 Model of induction machine in dqo arbitrary reference frame Some of the machine inductances are functions of the rotor speed, whereupon the coefficients of the statespace equations (voltage equations), which describe the behaviour of the induction machine, are timevarying (except when the rotor is at standstill). A change of variables is often used to reduce the complexity of these statespace equations. There are several changes of variables, which are used but there is just one general transformation [34]. This general transformation refers the machine variable to a frame of reference, which rotates at an arbitrary angular velocity ωg. In this reference frame the machine windings are replaced with some equivalent windings as shown in Fig. 4.16. Fig. 4.16. Induction machine windings in the dqoarbitrary reference frame. Based on (4.2.1) and using the general transformation, the voltage equations in the arbitrary reference frame can be written as: 33 3434 Wind Turbine Blockset in Matlab Simulink 0 0 vsd R s Ψ sd 0 −ωg 0 isd v i Ψ ωg 0 0 0 0 Rs 0 sq sq sq 0 0 vso iso d Ψ so 0 0 0 Rs = + + − ( ωg − ωr ) 0 0 0 0 Rr v rd i rd dt Ψ rd v rq i rq Ψ rq 0 Rr 0 0 0 ( ωg − ωr ) 0 Ψ v i R ro ro ro 0 0 0 r 0 0 0 Ψ sd 0 Ψ sq 0 Ψ so (4.2.21) 0 Ψ rd Ψ 0 rq Ψ ro 0 where: ωr is the electrical speed of the machine and ωg is the speed of the general reference frame. The compact form of (4.2.21) is: d dt [ V ] = [ R ][ I] + [ Ψ ] + [Ω][ Ψ ] (4.2.22) The relation between the linkage fluxes and the currents is given by: Ψ sd Ls Ψ 0 sq Ψ so 0 = Ψ rd L m Ψ rq 0 Ψ ro 0 0 0 Lm 0 Ls 0 0 Lm 0 Lσs 0 0 0 Lr Lm 0 0 0 0 0 0 Lr 0 0 0 isd 0 isq 0 iso 0 i rd 0 i rq Lσr i ro (4.2.23) or in compact form: [ψ ] = [ L][ I] (4.2.24) where: 3 M sr is the mutual inductance in the general reference frame; 2 • Lm = • Ls = Lσs + L m is the stator self inductance; • L r = Lσr + L m is the rotor self inductance. The voltage equations are written again in terms of currents and flux linkages. Clearly, these variables are related based on the matrix inductance [ L ] and both cannot be independent or state variables. Fluxes as state variables If the fluxes are selected as state variables the currents can be expressed in compact form by: [ I] = [ L] [ Ψ ] −1 The inverse of the inductance matrix is: 34 (4.2.25) 3535 Wind Turbine Blockset in Matlab Simulink 0 0 −Lm 0 0 Lr 0 Lr 0 0 − L m 0 D 0 0 0 0 0 Lσr 1 −1 [ L] = 0 Ls 0 0 D −L m 0 0 Ls 0 0 −Lm 0 D 0 0 0 0 0 Lσs (4.2.26) where D = Ls L r − L2m . After some mathematical manipulation, the statespace form of the dynamic equations (4.2.21) is: d −1 [ Ψ ] = − [ R ][ L] [ Ψ ] + [ V ] dt (4.2.27) or in an expanded form: RL R s Lr −ωg 0 − s m 0 D D R s Lr R s Lm ω − 0 0 g Ψ sd D D Ψ R s sq 0 0 0 0 Lσs d Ψ so = − R r Ls dt Ψ rd − R r Lm − ( ωg − ωr ) 0 0 Ψ rq D D R r Ls RL Ψ ro 0 − r m 0 ( ωg − ωr ) D D 0 0 0 0 0 0 0 Ψ sd vsd Ψ v 0 sq sq Ψ so vso + Ψ rd v rd 0 Ψ rq v rq 0 Ψ ro v ro Rr Lσr (4.2.28) Currents as state variables If the currents are selected as state variables the fluxes will be replaced based on (4.2.23) Assuming that all parameters are constant (saturation and other effects are not taken into account) (4.2.21) become: d [ I] −1 −1 −1 = − [ L ] [ R ] + [ L ] [ Ω ][ L ] [ I ] + [ L ] [ V ] dt ( ) (4.2.29) So the statespace form of the dynamic equations using currents as state variables is: d [ I] = [ A ][ I ] + [ B][ V ] dt (4.2.30) where: 35 3636 Wind Turbine Blockset in Matlab Simulink RL LL L2 R r Lm Lr Lm − s r + s r ωg − m ( ωg − ωr ) 0 ωr 0 D D D D D R s Lr L2m Lr Lm R r Lm Ls L r 0 0 − ω + ω − ω − − ω ( ) r r D g D g D D D RD 0 0 0 0 0 − s Lσs [ A] = R s Lm Ls L m R r Ls Ls L r L2m − ω − − ω + ω − ω 0 0 ( ) g g g r D D D D D 2 Ls L m R s Lm Ls L r R r Ls Lm ω ω − ω − ω − 0 0 ( g r) g g D D D D D RD − r 0 0 0 0 0 Lσr 0 0 −L m 0 0 Lr 0 Lr 0 0 − L m 0 D 0 0 0 0 0 L σr 1 −1 [ B] = [ L] = 0 Ls 0 0 D −Lm 0 0 Ls 0 0 −L m 0 D 0 0 0 0 0 L σs Electromagnetic torque The electromagnetic torque can be obtained starting from (4.2.22) and multiplying it from the left with the transpose of the currents vector: [ I] [ V ] = [ I] [ R ][ I] + [ I] t t t d t [ Ψ ] + [ I] [Ω][ Ψ ] dt (4.2.31) where • Pi = [ I] [ V ] is the instantaneous power; • Pcopper = [ I ] [ R ][ I ] are the copper losses in the machine windings; • d [ψ] is the magnetic power stored in machine (due to the variation in time dt of the magnetic energy); • Pm = [ I ] [ Ω ][ Ψ ] is the mechanical power. t t Pmag = [ I ] t t The electromagnetic torque is then: Te = Pm P 3 = p m = p ( Ψ sd isq − Ψ sq isd ) Ωr ωr 2 (4.2.32) where: p is the number of pole pairs and Ωr is the mechanical speed of the rotor. Other equivalent expressions for the electromagnetic torque of an induction machine are: 36 3737 Wind Turbine Blockset in Matlab Simulink • Te = 3 p ( Ψ rq i rd − Ψ rd i rq ) 2 • Te = 3 pL m ( isq i rd − isd i rq ) 2 • Te = 3 Lm p ( Ψ sq Ψ rd − Ψ sd Ψ rq ) 2 D The above equations may be somewhat misleading since they seem to imply that the leakage inductances are involved in the energy conversion process. This, however, is not the case. Even though the flux linkages contain the leakage inductances, they are eliminated by the algebra. The Simulink implementation and the mask interface for this model have the same structure as for the abc/abc model (see Fig. 4.14 and Fig. 4.15). 4.2.1.3 Reduced order model The theory of neglecting electric transients is presented in [34]. For balanced steadystate operation, the variables in synchronous reference frame are constants. Hence, the electric transients of the stator can be neglected by setting the derivative of the stator fluxes in the d and q axis to zero in [34] as: R s Lm R s L rr −ω − 0 e D • D ψ v sd 0 0 0 0 sd R s L rr R s L m ψ sd ωe − 0 v ψ sq D D sq = 0 0 0 0 ⋅ ψ sq + (4.2.33) ⋅ ψ v rd 0 0 1 0 ψ rd R r L m R L rd r ss 0 − ( ωe − ωr ) − ψ rq D D v ψ 0 0 0 1 rq rq R r Lss RL 0 − r m ( ωe − ωr ) D D In order to avoid the computational problems due to algebraic loops (4.2.33) should be written in state space form only in terms of the rotor fluxes as: τrm ατrm σr − σs ψ rd ψ = rq ατrm L m + 1 ωe − ωr − τsr • where: τsr = ατrm ατrm ατrm vsd L m + 1 ωe − ωr 1 0 2 v τ sr τsr ψ rd τsr sq (4.2.34) ψ + ατ ατ v rm τ rq − 2 rm 0 1 rd ατrm σr − rm τ sr τsr v rq σs R s L rr RL L L , τrm = r m , σ r = m , σs = m , α = D L rr Lss D 1 ω 1+ e τsr 2 . Based on rotor fluxes and input voltages the stator fluxes can be obtained using: 37 3838 Wind Turbine Blockset in Matlab Simulink α α α L ασ ω ωe r m e 2 τsr τsr τ sr Ψ sd ψ v rd + sd Ψ = α vsq ψ rq α sq − α L ω τ m e ασ r − τ2 ωe τ sr sr sr (4.2.35) Then the electromagnetic torque is calculated based on: Te = 3 Lm p⋅ (ψ sq ⋅ ψ rd − ψ rq ⋅ ψ sd ) 2 D (4.2.36) However, the reduced order model, which neglects the stator transients, can be extracted from the full order model given by (4.2.30) written in synchronous reference frame using Control System Toolbox from Matlab. This toolbox is a very powerful tool for reduced order modeling. Based on some specific functions as modred, ssbal, minreal, and balreal the desired state variables can be eliminated from the state space equations and the reduced order models are obtained. The Simulink model for this machine is based on (4.2.34)(4.2.36) as shown in Fig. 4.17. Fig. 4.17. Simulink implementation of the reduced order model for the induction machine. The mask interface for this model has the same structure as for the abc/abc model (see Fig. 4.15). The only difference consists in the initial conditions for the state variables, which are the fluxes in this case. 4.2.1.4 Steadystate model The voltage equations, which describe the balanced steadystate operation of an induction machine, can be obtained in several ways [34]. For balanced steady state conditions the d and q variables are sinusoidal in all reference frames except the synchronously rotating reference frame wherein they are constants. In the following paragraphs the prime index will be used for the rotor equations in order to highlight that these are related to the stator circuit of the machine. 38 3939 Wind Turbine Blockset in Matlab Simulink Starting from (4.2.21) written in the synchronously rotating reference frame ωg = ωe and neglecting the time rate changing of all flux linkages and then employ the relationships: ~ 2 Fas = Fsd − jFsq ~ ' ar (4.2.37) 2 F = F − jF ' rd ' rq the steadystate voltage equations for induction machine can be written as follows: ( V s = ( R s + jX σs ) Is + jX m Is + I r V = ( R r + jsX ' r ' σr )I ' r ' ( ) + jsX m Is + I ' r ) (4.2.38) ωe − ωr 2πf s − pΩ r pΩ r = = 1− ; f s is the frequency of stator voltage; Ω r is the ωe 2πf s 2πf s mechanical rotor speed; p is the number of pole pair. where: slip s = Equation (4.2.38) suggests the perphase equivalent circuit shown in Fig. 4.18. Fig. 4.18. Equivalent circuit for steady state operation of a symmetrical induction machine. The steadystate voltage equations (4.2.38) can be arranged as follows: V s = ( R s + jX s ) Is + jX m I r ' V r = ( R 'r + jsX 'r ) I r + jsX m I r ' ' ' (4.2.39) where: X s = X σs + X m is the self reactance of the stator windings; X 'r = X 'σr + X m is the self reactance of the rotor windings. Or in matrix form: V s R s + jX ss '= V r jsX m jX m Is ⋅ R + jsX 'rr I'r ' r (4.2.40) The compact form of (4.2.40) is: V = Z ⋅ I (4.2.41) 39 4040 Wind Turbine Blockset in Matlab Simulink Based on (4.2.41) the stator and rotor currents are obtained from: −1 I = Z ⋅ V (4.2.42) So, the induction machine can be treated as a complex quadripole as shown in Fig. 4.19 where the complex impedance is a function of the slip. Fig. 4.19. Equivalent perphase quadripole for the induction machine used in the steadystate analysis. Recalling (4.2.31) the electromagnetic torque is given by: Te = * ' 3 pX m Re jIs I r 2 (4.2.43) * where Is is the conjugate of Is . The active power is given by: * 3 Px = ⋅ Re al(V x ⋅ I x ) 2 (4.2.44) The reactive power is given by: * 3 Q x = ⋅ Im ag(V x ⋅ I x ) 2 (4.2.45) These equations can be used both for stator and rotor windings by replacing the indices x with s and r respectively. Using this approach the induction machine can be completely characterized in terms of steadystate values of the stator and the rotor currents, phase angle for both currents, electromagnetic torque, active and reactive power for both stator and rotor circuit. Moreover, an analysis, both in motor and generator mode of operation, can be performed. The Simulink implementation of the steadystate model, both for a squirrelcage machine and wound rotor one, is shown in Fig. 4.20. Fig. 4.20. Simulink implementation for steadystate model of induction machine. This model uses some particular blocks from the DSP Blockset (Simulink) such as: Matrix Multiplication, Horizontal and Vertical Concatenation, and LU Inverse (Matrix inverse using 40 4141 Wind Turbine Blockset in Matlab Simulink LU factorization). It has been chosen to use these blocks because the model also implements the rotor deepbar effect. Using these blocks the model can be implemented in matrix form with the minimal number of blocks. However, the model can also be implemented using standard blocks from Simulink. 4.2.1.5 Modelling deepbar effect In larger machines like MW generators rotor deepbar effect causes the equivalent rotor circuit resistance and leakage reactance to vary significantly with the slip [3], [46]. Fig. 4.21 illustrates the leakage flux paths for a squirrelcage induction machine with deepbars. The parts of each bar extending deeply into the rotor iron have higher leakage inductances than those parts of the bar crosssection near the air gap, because more of the leakage flux links the deeper parts of the bar. One may think of each bar as being composed of several layers of equal cross section, and thus of equal resistance, but with inductances increasing with depth. Fig. 4.21. Deepbar effect for a squirrel cage induction machine. Under running conditions, the slip is quite small and the frequency of the rotor currents is only 12 Hz. As a result, the leakage reactance is neglected, and the current distribution is essentially uniform throughout each rotor bar. The effective resistance of the rotor is that of all layers in parallel. This low resistance makes for low slip and high efficiency at full load. As the motor is overloaded, however, slip and frequency increase. Thus, the value of rotor leakage reactance in the circuit model becomes smaller and the value of rotor resistance becomes greater. Another important need is to include the leakage path saturation effect into a reduced model. The main saturation leakage path is at the tooth tips over the closed or nearly closed slots [46]. However, large machines usually have large slots openings, and hence the leakage path saturation effect is not significant for large machines [3]. Numerous approaches regarding the modelling of deepbar effect have been presented in literature; therefore, all these models are complicated and require complex mathematics [3], [41], [42]. In [3] is presented a simplified approach for the modelling of deep bar effect, which gives accurate results. The equivalent rotor winding resistance and leakage reactance are modelled as follows: 41 4242 Wind Turbine Blockset in Matlab Simulink R 'r = R 'r _ DC + (R 'r _ 50Hz − R 'r _ DC ) ⋅ s X 'σr = X 'σr _ DC + (X 'σr _ 50Hz − X 'σr _ DC ) ⋅ s (4.2.46) where: • R 'r _ DC is the equivalent rotor winding resistance at zero frequency (DC value or at • synchronous speed) ; R 'r _ 50Hz is the equivalent rotor winding resistance at 50Hz (at standstill the slip s=1); • X 'σr _ DC is the equivalent stator leakage reactance at zero frequency (DC value or at • synchronous speed); X 'σr _ 50Hz is the equivalent stator leakage reactance at 50Hz (at standstill the slip s=1). Therefore, this approach assumes a linear dependency of rotor resistance and rotor leakage reactance against the entire range of slip. Moreover, it requires extra parameters determined from tests at synchronous speed. Usually, the parameters from data sheets are given for rated operating point, near the synchronous speed and the method cannot be applied. Based only on the parameters from datasheet, the equivalent rotor winding resistance and leakage inductance at any operating point can be modelled as: R 'r _ dbe R 'r , s ≤ s m = ' 1 − K rsm R r ( K r − 1) s + 1 − s , s > s m m X σ' r _ dbe X σ' r , s ≤ s m = ' 1 − K xsm X σr ( K x − 1) s + 1 − s , s > s m m (4.2.47) where: • R 'r and X 'σr are the values the from datasheet given for rated operating point • K r is the resistance deepbar factor; • K x is the leakage reactance deepbar factor; • sm = R ' r R s2 + X s2 is the breakdown slip in motoring mode. ( X 2m − Xs X'r ) + R s2 X'2r This approach assumes that the parameters of induction machine are constant for the values of the slip less or equal with the breakdown slip. Resistance and leakage reactance deepbar factors can be determined based on rated voltage and current, starting current and phase angle (from shortcircuit test at rated current). If the phase angle is not available, a value around 80o can be used as a good approximation in the first step of iteration. Then it can be corrected so that the calculated starting torque will fit the value from datasheet or test reports. 42 4343 Wind Turbine Blockset in Matlab Simulink The deepbar factors for rotor resistance and rotor leakage reactance are calculated as: Kr = 1 R 'r U ph _ rated 1 R − s 2 β I ph _ rated 1 + tan ϕsc (4.2.48) 1 K x = ' ( R s + K r R 'r ) ⋅ tan ϕsc − X σs X σr where: • U ph _ rated is the rated voltage per phase; • I ph _ rated is the rated current per phase; • β is the relative starting current (from data sheet); • R s is the stator resistance perphase (from data sheet); • X σs is the stator leakage reactance perphase (from data sheet); • R 'r is the rotor resistance perphase for the rated operating point (from data sheet); • X 'σr is the rotor leakage reactance perphase for the rated operating point (from data sheet); • ϕsc is the phase angle at standstill (s=1): from shortcircuit test at rated current or around 800 for large induction machines. Fig. 4.22a shows the rotor resistance and the rotor leakage reactance against slip based on (4.2.47) and (4.2.48), while Fig. 4.22b shows the electromagnetic torque with and without deepbar effect. a) b) Fig. 4.22. Deepbar effect for a 2 MW squirrelcage induction machine as a function of slip: a) rotor resistance rotor leakage reactance and b) electromagnetic torque. The same approach can be used for doubly fed induction machine since the manufacturers performs the standard tests as for the squirrelcage induction machines. The Simulink model of induction machine, which implements the variation of the rotor parameters with the slip, is shown in Fig. 4.23. 43 4444 Wind Turbine Blockset in Matlab Simulink Fig. 4.23. Simulink model of the induction machine in dq rotor reference frame including deepbar effect. The mask interface for this model is shown in Fig. 4.24. Fig. 4.24. Mask interface for the induction machine in dq rotor reference frame including deepbar effect. In order to include the deepbar effect some additional parameters have been added in the mask interface: rated voltage and current, base frequency, number of pole pairs, starting current, phase angle at standstill. All these parameters are usually available in standard datasheets. 44 4545 Wind Turbine Blockset in Matlab Simulink 4.2.1.6 Modelling saturation The flux equations highlight the terms, which are dependent on magnetic saturation and therefore introduce nonlinearity. It will be assumed that the saturation of the leakage fluxes into a large induction machine can be neglected [33], [46] thus, the stator and rotor leakage inductances are constant. The magnetizing inductance is a function of magnetization current as shown in Fig. 4.25. Fig. 4.25. Magnetizing flux versus magnetizing current. Assuming that the machine has magnetic, electric and geometric cylindrical symmetry, then the magnetic saturation caused by the main magnetic field is independent by the direction of this field and depends only on the absolute value of it. The magnetizing flux can be written as: d ( L m ⋅ i m ) di m d d ⋅ ( Ψ m ) = ( Lm ⋅ im ) = dt dt di m dt (4.2.49) The magnetization current is independent on the reference frame. Based on magnetization curve shown in Fig. 4.25 two inductances can be defined as follows: Lm = Ψm = tgα1 = f1 (i m ) im (4.2.50) Ld = d ( L m ⋅ i m ) dΨ m = = tgα 2 = f 2 (i m ) di m di m (4.2.51) where: Lm is the magnetizing inductances and Ld is the dynamic inductance. Fig. 4.26 shows some typical dependencies of these two inductances with the magnetizing current. Fig. 4.26. Magnetizing and dynamic inductance versus magnetizing current. 45 4646 Wind Turbine Blockset in Matlab Simulink Based on the noload curve for an induction machine these two inductances can be determined as a function of the noload current. Then, using a mathematical approximation e.g. a function or a lookup table, the desired values for these inductances for a given current can be found. Taking into account the above mentioned, the phasor voltage equations of the induction machine in the general reference frame including the saturation of the main path can be written as follow: ( ) ( ) d di m is + L d + jωg Lσs is + L m i m dt dt d di m v r = R r i r + Lσr i r + Ld + j ( ωg − ωr ) Lσr i r + L m i m dt dt i m = is + i r vs = R s is + Lσs ( ) ( ) =f (i ) ( ) (4.2.52) L m = f1 i m Ld 4.2.1.7 2 m Modelling iron losses Considerable attention has been paid during the last years in modelling of iron losses into an induction machine. A simple search in IEEE database will produce around 120 results. All these papers deals with the iron loss modelling for an inverterfed induction motor and take into account the hysteresis losses as well as the eddycurrent losses in a wide range of speed and frequency. In wind turbine applications usually a large squirrelcage induction generator is directly connected to grid and the rotor speed is almost fixed. When a doubly fed machine is used, again the stator windings are directly connected to the grid and the variable speed operation is achieved using a power converter in the rotor side. Usually, the speed varies in a range of ± 30% from the synchronous speed. In data sheets the noload power, current and power factor are provided for the machine rated voltage as well as the electrical parameters for the rated operating point. So, it is convenient to use only this information in the case of a large induction machine, both squirrelcage and wound rotor. A classical approach in iron losses modelling for an induction machine involves a resistance inserted in parallel with the magnetizing reactance in the equivalent perphase diagram, as shown in Fig. 4.27. Fig. 4.27. Equivalent diagram perphase for induction machine including iron losses resistance. 46 4747 Wind Turbine Blockset in Matlab Simulink One major drawback of this approach is that the magnetizing current is divided into two components: an active component and a reactive component. So, it is difficult to use this equivalent diagram for a dynamic model as well as for a steady state analysis. A modified equivalent scheme, called Γ has been developed. Unfortunately, the currents in this equivalent scheme are multiplied with some factors and the dynamic modelling become difficult. A simple approach assumes the iron losses resistance in series with the magnetizing reactance as shown in Fig. 4.28. Fig. 4.28. Equivalent T diagram of the induction machine with iron losses. Where the parameters of the magnetizing branch can be determined based on parameters from Fig. 4.27 as follows: Ro = R Fe ⋅ X 2m  equivalent series resistance; R 2Fe + X 2m Xo = R 2Fe ⋅ X m  equivalent series reactance. R 2Fe + X 2m Since the machine manufacturers give the electrical parameters for both schemes presented in Fig. 4.27 and Fig. 4.28 in the data sheet, or at least for the first one it is easier to use this approach. The electrical parameters for the parallel mode can be transformed through some simple mathematical calculations to the series mode equivalent circuit. Recalling (4.2.1) the iron losses can be introduced into the dynamic equations written in the abc/abc reference frame as follows: d dt [ V ] = [ R ][ I] + [ R o ][ I] + [ Ψ ] (4.2.53) where the iron losses resistance matrix is defined in a similar manner as (4.2.6) by: 1 0 0 [R 0 ] = R o ⋅ 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 = R o ⋅ [ To ] 0 0 1 (4.2.54) In this condition the inductance matrix (4.2.10) will be rewritten as: 47 4848 Wind Turbine Blockset in Matlab Simulink 0 0 M o f1 M o f 2 M o f 3 Ls 0 Ls 0 M o f3 M o f1 M o f 2 0 0 L Mf Mf Mf [ L] = M f M f M sf Lo 2 0o 3 0o 1 o 3 o 2 r o1 M o f 2 M o f1 M o f 3 0 Lr 0 0 L r M o f 3 M o f 2 M o f1 0 (4.2.55) where: • • • 3 Ls = Lσs + M o 2 3 L r = Lσr + M o 2 Mo  stator selfinductance;  rotor selfinductance;  magnetizing inductance for the series circuit. Finally the state space equations with fluxes as state variables are: d −1 [ Ψ ] = − ([ R ] + [ R 0 ] ) [ L ] [ Ψ ] + [ V ] dt (4.2.56) or using currents as state variables: d [ I] d [ L] −1 −1 = [ L ] ⋅ − ([ R ] + [ R 0 ]) − ωr ⋅ ⋅ [ I] + [ L] ⋅ [ V ] dt dθ r (4.2.57) A similar approach can be used to introduce the iron loss resistance into the dynamic equations written in arbitrary reference frame. 4.2.2 Synchronous Machine Model In the following the dynamic equations for salientpole synchronous machines are presented. 4.2.2.1 Dynamic Equations of Synchronous Machine The dynamic equations of synchronous machine are derived for a salient pole machine with damping bars, which is the most complicated model for this type of machine [7], [9], [21], [33], [34], [46], [75]. This model can be also used for a nonsalient pole machine assuming identical values for inductances in d and qaxis. Untransformed model In synchronous machines with salient poles the damping bars can be represented by a directand a quadratureaxis damping winding, which are shortcircuited windings. Assume that on the stator of the machine there is a symmetrical threephase sinusoidal distributed stator winding and on the rotor there are the field winding and direct and quadratureaxis damper winding. Schematic of the machine is shown in Fig. 4.29. 48 4949 Wind Turbine Blockset in Matlab Simulink sQ sB rq sB ωr Q sA D sA sD θr f rd sC sC Fig. 4.29. Schematic of the threephase synchronous machine. The phasor form of the stator voltage equation in the stationary reference frame is: u s A, B, C = R s ⋅ i s A, B, C + dψ s A, B, C (4.2.58) dt where R s is the resistance of the stator winding. The voltage equations of the field winding and the damper windings in the reference frame fixed to the rotor are as follows: dψ f dt dψ D u D = R D ⋅ iD + dt dψ Q u Q = R Q ⋅ iQ + dt u f = R f ⋅ if + (4.2.59) Transformed model with quadraturephase stator windings Replacing the threephase stator winding by an equivalent quadraturephase stator winding, the number of differential equations can be reduced. Schematic of the corresponding model is shown in Fig. 4.30. oaxis so sq qaxis Q sd D f ωr daxis Fig. 4.30. Schematic of the transformed model for salientpole synchronous machine. The stator voltage equations are: 49 5050 Wind Turbine Blockset in Matlab Simulink u sD = R s isD + u sQ = R s isQ + dΨ sD dt dΨ sQ (4.2.60) dt The voltage equations of the field winding and the damper windings in the reference frame fixed to the rotor are the same. Commutator model It can be shown that when the stator quantities of the salient pole machine are expressed in the dqo reference frame fixed to the rotor, the components of the stator flux linkages expressed in the rotor reference frame will not contain the rotor angle. This corresponds to the commutator model of the salientpole machine. This approach is useful for studying the unsymmetric behaviour of the synchronous machine. It follows that in the rotor reference frame, under linear conditions, the salient pole machine can be described by a system of voltage equations with constant coefficients and, in general, the only changing quantity is the rotor speed in these equations. Thus, the stator equations in phasor form in the rotor reference frame are given by: u s = R s is + dΨ s + jωr Ψ s dt (4.2.61) where ω r is the rotor speed. The stator flux linkages in the dqo rotor reference frame can be expressed in terms of the currents of the machine as: 3 3 3 M Di D Mf if + ψ sd = L σs + M 0 + L 2 ⋅ i sd + 2 2 2 3 3 M Qi Q ψ sq = L σs + M 0 − L 2 ⋅ i sq + 2 2 ψ so = (L σs − 2M 0 ) ⋅ i so (4.2.62) We note: 3 L d = L σs + M 0 + L 2 2  directaxis synchronous inductance 3 L q = L σs + M 0 − L 2 2  quadratureaxis synchronous inductance; L o = L σs − 2M 0  homopolar inductance; M 'f = winding; 50 3 Mf 2  the mutual inductance between the field winding and the stator 5151 Wind Turbine Blockset in Matlab Simulink 3 MD 2 axis damper winding; M 'D =  the mutual inductance between the stator winding and the direct 3 MQ  the mutual inductance between the stator winding and the 2 quadratureaxis damper winding . M 'Q = The synchronous inductances can be expressed as follows: L d = L σs + L md L q = L σs + L mq where L σs  the stator leakage inductance; 3 L md = M 0 + L 2  the magnetizing inductance in directaxis 2 3 L mq = M 0 − L 2  the magnetizing inductances in quadratureaxis, 2 Splitting (4.2.61) into real and imaginary part yields the following dqo equations for the stator are obtained: dψ sd − ωr ⋅ ψ sq dt dψ sq u sq = R s ⋅ i sq + − ωr ⋅ ψ sd dt dψ u so = R s ⋅ i so + so dt u sd = R s ⋅ i sd + (4.2.63) where ψ sd = L d i sd + M 'f i f + M 'Di D ψ sq = L q i sq + M 'Qi Q ψ so = L o i so The voltage equations of the field winding and the direct and quadratureaxis damper windings can be obtained in the reference frame fixed to the rotor as: dψ f dt dψ D u D = R D ⋅ iD + dt dψ Q u Q = R Q ⋅ iQ + dt u f = R f ⋅ if + (4.2.64) 51 5252 Wind Turbine Blockset in Matlab Simulink where Rf, RD and RQ are the resistances of the field winding and direct and quadratureaxis damper windings respectively, and the flux linkages of the field winding and the damper windings are obtained as: ψ f = L f i f + M fDi D + M 'f i sd ψ D = L Di D + M fD i f + M 'D i sd (4.2.65) ψ Q = L Qi Q + M 'Q i sq where Lf  self inductance of the field winding; L D , L Q  selfinductances of the damper winding in the direct and quadratureaxis; M fD  mutual inductance between field winding and direct damper winding Thus the matrix form of the voltage equations of the commutator model are obtained as follows: ω r M 'Q i sd 0 pM 'f pM 'D u sd R s + pL d − ω r L q u ω r M 'f ω r M 'D 0 pM 'Q i sq sq ω r L d R s + pL q i u so 0 0 R s + pL o 0 0 0 ⋅ so (4.2.66) = ' if 0 0 R f + pL f pM f 0 u f pM fD i uD pM 'D 0 0 pM fD R D + pL D 0 D ' u Q 0 pM Q 0 0 0 R Q + pL Q i Q In case of symmetrical and balanced operation of the synchronous machine the homopolar component of the stator current is zero, so the homopolar component of the stator flux is zero. Electromagnetic torque The general expression of the electromagnetic torque is: Te = ( 3 P Ψ s × is 2 ) (4.2.67) Taking into account the expressions of the stator flux and stator current in the rotor reference frame the electromagnetic torque can be written as follows: ( ) [( ) ( ) ] 3 3 Te = P ψ sd i sq − ψ sq i sd = P L d − L q i sd i sq − M 'Q i Qi sd + M 'f i f + M 'Di D i sq (4.2.68) 2 2 4.2.2.2 Simulink model The dynamic equations (4.2.66) and (4.2.68) are implemented in Matlab/Simulink, as shown in Fig. 4.31. 52 5353 Wind Turbine Blockset in Matlab Simulink Fig. 4.31. Simulink implementation of synchronous machine model. In Fig. 4.32 is shown the mask interface for this model in Simulink. Fig. 4.32. Mask interface for synchronous machine model in Simulink. Since the model is written for a salient pole synchronous machine, the mask parameters include different values for the rotor resistances and leakage inductances in d and qaxis as well as for the magnetizing inductances. An extra input are the parameters for the field winding. 4.2.3 Permanentmagnet synchronous machine The analysis of permanent magnet (PM) synchronous machines can be made using an quadrature equivalent circuit where the damping windings are replaced with two equivalent windings D, Q in direct and quadrature axis respectively and the permanent magnet is replaced with an equivalent superconductor winding placed in the directaxis (Fig. 4.33). The current through the equivalent winding of the permanent magnet ( I f ) is constant in all mode of operation. 53 5454 Wind Turbine Blockset in Matlab Simulink sQ sB rq sB ωr Q sA D f sA sD θr rd sC sC Fig. 4.33. Schematic of PM synchronous machine with damper winding. The voltage equations of the PM synchronous machines in dqo rotor reference frame can be derivate from the voltage equation of synchronous machine (2.26) taking into account the above assertions. oaxis so sq qaxis Q sd D f ωr daxis Fig. 4.34. Schematic of the transformed model of PM synchronous machine in the dqo rotor reference frame. The voltage equations for PM synchronous machine can be written as: dψ sd − ωr ψ sq dt dψ u sq = R s isq + sq + ωr ψ sd dt dψ so u so = R s iso + dt dψ D u D = R Di D + dt dψ Q u Q = R QiQ + dt u sd = R s isd + (4.2.69) The flux linkage equations can be written as: ψ sd = L d i sd + M 'f I f + M 'D i D ψ sq = L q i sq + M 'Q i Q ψ so = L 0 i so ψ f = L f I f = ct. ψ D = L D i D + M fD I f + M 'D i sd ψ Q = L Q i Q + M 'Q i sq 54 (4.2.70) 5555 Wind Turbine Blockset in Matlab Simulink The matrix form of the voltage equations for PM synchronous machines in dqo rotor reference frame is: u sd R s + pL d − ωr L q − ωr M 'Q i sd 0 0 pM 'D u ' ' ' i ω + ω ω L R pL 0 M M pM sq r d s q r f r D Q sq i u so 0 0 R s + pL 0 0 0 0 ⋅ so (4.2.71) = 0 If 0 0 0 0 0 0 i 0 = u D pM ' 0 0 0 R D + pL D 0 D D ' = 0 u 0 pM Q 0 0 0 R Q + pL Q i Q Q The term pM 'f corresponds to the e.m.f. induced by the rotating magnet in the equivalent quadratureaxis stator winding. The general expression of the electromagnetic torque is given by: Te = 3 3 P(ψ s × i s ) = P(ψ sd i sq − ψ sq i sd ) 2 2 (4.2.72) The Simulink implementation of (4.2.71) (4.2.72) is shown in Fig. 4.35. Fig. 4.35. Simulink model for permanent magnet synchronous machine. The model take into account the iron losses as well as the parameters variation with the operating temperature. The mask interface for this model is presented in Fig. 4.36. Notice that the operating temperature is an input parameter in the model and there is also the possibility to take the iron losses into account. 55 5656 Wind Turbine Blockset in Matlab Simulink Fig. 4.36. Mask interface for the permanent magnet synchronous machine in Simulink. 4.3 Power converter models In this paragraph some aspects regarding the modelling of the power converter topologies for wind turbine applications will be presented. The focus is on the modelling of different softstarterfed induction machine topologies as well as on the voltage source converters modelling using the switching function concept. The influence of the load connection is taken into account in both cases. 4.3.1 AC Controllers (Softstarters) Many authors deal with the theory of operation for threephase ACcontrollers (softstarters), [67], [69], and [73]. However, only few of them present in details the theory of operation with a 3phase resistiveinductive load [67] and [73]. Especially the influence of the load connection type and the steadystate analysis are not presented in detail in the literature. In [70] a comparison for full and halfwave controlled load both for star and delta connection is carried out, unfortunately, the branch delta connection is not treated. An analysis of variablevoltage thyristor controlled induction motors in terms of specific operation modes, harmonic content and performances for a star connected machine is presented in [26]. In [74] is presented a hybrid abcdqo model for the induction machine, unfortunately this model has a complicated mathematical model and it can be used only for a star connection of the stator windings. A generalized approach in modelling of the power convertersfed induction machine, which involves a timedomain static network, is presented in [22]. In order to eliminate the computation errors this model uses some “dummy” shunt resistances. However, the proposed method can only be used for a star or deltaconnected machine. Moreover these approaches do not take into account the deepbar effect, which is present into a large induction machine. 56 5757 Wind Turbine Blockset in Matlab Simulink Therefore in analysis of a softstarterfed induction machine is very convenient to use the abc/abc model, since it is based on perphase quantities. At present many wind turbines, up to 2.3 MW, are based on the “Danish concept” in which a squirrelcage induction generator is directly connected to the grid as shown in Fig. 4.37. Fig. 4.37. Block diagram of directly gridconnected wind turbine. The scheme comprises the wind turbine rotor, linked via a gearbox to a generator, which through an electrical interface is connected to the grid. A control system is necessary to assure a proper operation of the wind turbine under all conditions. The electrical interface consists of: a softstarter, a capacitor bank and a transformer, which makes the connection with the medium voltage grid. The capacitor bank is used to control the power factor of the generator output. Softstarters are used only during the startup sequence of the generator in order to limit the inrush currents and the starting torque transients in the drive train. There are many configurations of softstarters, which feed an induction machine. However there are three topologies interesting for wind turbine applications, as presented in Fig. 4.38. Fig. 4.38. Possible configurations of softstarterfed induction machine in wind turbine applications: a) star connection, b) delta connection and c) branchdelta connection. The star and delta configurations have basically the same layout for the semiconductors (SCRs), the difference consists in the machine winding connection. There are two antiparallel thyristors on each phase, each one conducting on the positive cycle of the applied voltage. For a star connection the applied per phase voltages depend on the onstate of the SCRs on each phase. So, the softstarter can operate only when two or three SCRs are conducting in the same moment. The similar considerations can be done for a delta connection. However a softstarter with a branch delta connection, will operate only with one SCR conducting at a given moment. 57 5858 Wind Turbine Blockset in Matlab Simulink In wind turbine applications mainly the delta connection for the induction machine are used because the current rating of the stator windings can be reduced, and the third harmonic in the line currents is eliminated in this case. Depending on the firing angle α, three different modes of operation of the softstarter can be distinguished when a star or delta connected resistive load is used [67]: • Mode 1  0 ≤ α < 600 • Mode 2  600 ≤ α < 900  two SCRs are conducting; • Mode 3  900 ≤ α < 1500  none or two SCRs are conducting.  two or three SCRs are conducting (in either direction); where α is the firing angle for the softstarter. When a resistiveinductive load is used the analysis of the controller is difficult, since the operation modes depend on the extinction angle ξ and the limit angle αlim, both dependent on the phase angle ϕ. Mode 2 of operation, characterized by rapid changes of the output current, is not possible due to the load inductance. The ranges of the two remaining operation modes are: • Mode 1: ϕ ≤ α < α lim  two or three SCRs are conducting; • Mode 3: α lim ≤ α < 150o  none or two SCRs are conducting. The limit angle can be determined numerically from: 4π π − sin α lim − ϕ − 3tan ϕ −1 3 2e = π − sin ( α lim − ϕ ) 2 − e 3tan ϕ (4.3.1) Switching functions SwA, SwB, and SwC in two levels can be introduced in modelling of the SCRs and defined as equal to one when a given thyristor is conducting and equal to zero otherwise as shown in Fig. 4.39. Fig. 4.39. Modelling two anti parallel SCRs using a switching function. These switching functions can be determined based on phase voltage and phase (or line) current for each particular topology. 4.3.1.1 Star connected 3phase load A softstarterfed 3phase resistiveinductive load in star connection is shown in Fig. 4.40. 58 5959 Wind Turbine Blockset in Matlab Simulink Fig. 4.40. Softstarterfed 3phase RL load in star connection. Analyzing this topology the applied perphase voltages at the machine terminals can be written as: v UX SwA −SwB −SwC v AN v = 1 −S VY 2 wA SwB −SwC ⋅ v BN v WZ −SwA −SwB SwC v CN (4.3.2) where: [ v UX v VY v WZ ] are the phase voltages of the machine; [ v AN v BN vCN ] are the phaseT T toground supply voltages and SwA , SwB , SwC are the switching functions for each phase. Based on the phasetoground voltages and the phase (line) currents the switching functions are computed. The considered load has a power factor of 0.85. Using (4.3.1) a limit angle around 95o will be the result. In order to highlight the two possible operation modes for the softstarters two values for the firing angle have been used. A value of 70o for the firing angle will correspond to Mode 1 and 110o to Mode 3 respectively. The perphase voltage and the phase current for the considered load are shown in Fig. 4.41 for different firing angles. a) b) Fig. 4.41. Typical waveforms for voltage and currents for a 3phase star connected load: a) Mode 1and b) Mode 3 of operation. From Fig. 4.41 it can be observed that when the firing angle α is greater than the limit angle αlim the phase currents occurs in discontinuous, nonsinusoidal, alternating pulses. Therefore, a high harmonic content for the phase currents will be the result. 59 6060 Wind Turbine Blockset in Matlab Simulink 4.3.1.2 Delta connected 3phase load Fig. 4.42 shows a threephase resistiveinductive load in delta connection. Fig. 4.42. Softstarterfed 3phase RL load in delta connection. It can easily be demonstrated that the voltages applied to a deltaconnected load are given by: v UX SwASwB −0.5 ⋅ SwA −0.5 ⋅ SwB v AB v = −0.5 ⋅ S − ⋅ S S 0.5 S wC wB wC wB ⋅ v BC VY v WZ −0.5 ⋅ SwC −0.5 ⋅ SwA SwCSwA vCA (4.3.3) where: [ v AB v BC vCA ] are the linetoline supply voltages (or phase voltages for the load). T For a delta connection the switching functions are determined based on phasetoground voltages and line currents. Considering the same load as in the precedent case the applied perphase voltage, the phase and the line currents for different firing angles are shown in Fig. 4.43. a) b) Fig. 4.43. Typical waveforms for voltage and currents for a 3phase delta connected load: a) Mode 1 and b) Mode 3 of operation. Again, it can be observed that when the firing angle α is greater than the limit angle αlim the line currents occurs in discontinuous, nonsinusoidal, alternating pulses. Therefore, for this topology a high harmonic content for both line and phase currents will result especially in Mode 3. As a result the electromagnetic torque will have a high harmonic content. 60 6161 Wind Turbine Blockset in Matlab Simulink 4.3.1.3 Branchdelta connected load A softstarterfed 3phase resistiveinductive load in branchdelta connection is presented in Fig. 4.44. Fig. 4.44. Softstarterfed 3phase RL load in branchdelta connection. The output voltages of the softstarter in this case are: v UX SwA 0 0 v AB v = 0 S 0 ⋅ v BC wB VY v WZ 0 0 SwC v CA (4.3.4) For a branchdelta connection the switching functions are determined based on the line voltages and line currents. Again, a 3phase resistiveinductive load as in precedent cases has been considered. The perphase voltage, the phase and the line current are shown in Fig. 4.45 for the considered firing angles. a) b) Fig. 4.45. Typical waveforms for voltage and currents for a 3phase branch delta connected load: a) Mode 1 and b) Mode 3 of operation. Using a branchdelta connection the line current waveform is approximately the same for the entire range of possible firing angle. Therefore, the harmonic content in the line / phase currents as well as in the electromagnetic torque is not so high comparing with the other two topologies. 61 6262 Wind Turbine Blockset in Matlab Simulink 4.3.1.4 Simulink implementation of the softstarter The softstarter model including the connection type as well as the algorithm for generation of the switching function has been implemented in Simulink as shown in Fig. 4.46. Fig. 4.46. Simulink model of the softstarter. The mask interface for this model is shown in Fig. 4.47. The input parameters for this model are the fundamental frequency of the grid voltages and the initial states for the switches. The selection of the connection type is also an option in this mask. Fig. 4.47. Mask interface for the softstarter model in Simulink. 4.3.1.5 RMS model for softstarter In steadystate analysis of the softstarterfed induction machine or when reduced order models of the machine are used, an RMS model for the AC controller should be used. The expression of the softstarter output voltage is a function of firing angle, phase angle and socalled extinction angle. Since, a SCR will not permit the flow of reverse current, the conduction of it will end at a point ξ, which is called extinction angle or cutoff angle. The extinction angle can be obtained solving the following transcendental equation [73]: sin ( ξ − ϕ ) − sin ( α − ϕ ) e − cot ϕ( ξ−α ) = 0 where: α is the firing angle and ϕ is the phase angle (power factor). 62 (4.3.5) 6363 Wind Turbine Blockset in Matlab Simulink Only an iterative solution of (4.3.5) is possible. This yields the set of characteristics shown in Fig. 4.48. Fig. 4.48. Extinction angle versus firing angle for a series RL load with different power factors. When α = ϕ , which represents sinusoidal operation, (4.3.5) reduces to sin ( ξ − ϕ ) = 0 , which gives ξ = π + ϕ . The values ξ for sinusoidal operation are seen to lie on the dashed linear characteristic of Fig. 4.48. For a purely inductive load, therefore, the variation of ξ with α is linear, as shown in the ϕ = 90o characteristic from Fig. 4.48. Considering a single RL load as shown in Fig. 4.49 Fig. 4.49. Singlephase controller with a resistive inductive load. the applied voltage at the load terminals is given by [73]: VUX = VAB 1 ( ( ξ − α ) + 0.5sin 2α − 0.5sin 2ξ ) 2π (4.3.6) Based on (4.3.6) the variation of the fundamental voltage with firing angle is shown in Fig. 4.50 for several fixed values of phaseangle. Equation (4.3.6) is valid only for a branchdelta connection. Based on the input voltage and the phase current in complex form and a lookup table for the extinction angle the output voltage perphase can be obtained. Unfortunately, for the star or the delta connection due to the complexity of the conducting mechanism is very difficult to establish an analytical expression for the RMS ACcontroller output voltage. The Simulink implementation of (4.3.6) is based on a lookup table for the extinction angle as shown in Fig. 4.51. 63 6464 Wind Turbine Blockset in Matlab Simulink Fig. 4.50. Fundamental voltage versus firing angle for singlephase controller with a series RL load. Fig. 4.51. Simulink implementation of the RMS model for the softstarter. 4.3.2 Voltage Source Converters Currently, in wind turbine applications the backtoback voltage source converter (VSC) is mainly used. A block diagram of this power converter is shown in Fig. 4.52. Fig. 4.52. Structure of the backtoback voltage source converter. This topology comprises a double conversion from AC to DC and then from DC to AC. Both converters can operate in rectifier or inverter mode and therefore a bidirectional power flow can be achieved. A voltage source converter can be implemented in several ways: sixstep, pulse amplitude modulated (PAM) or pulse width modulated (PWM). Moreover, the implementation of a 64 6565 Wind Turbine Blockset in Matlab Simulink PWM VSC may be realized by three methods: harmonic elimination, “sinusoidal” PWM or space vector strategy (SVPWM). Independent on the implementation method the converter can be seen as a black box with some inputoutput characteristics as a function of the control strategy. Again, the switching function concept can be used in modelling of these topologies. So, in the following paragraphs the inputoutput characteristics as a function of the switch