Page d'accueil Table of contents

Table of contents

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?
Année:
2004
Langue:
english
Fichier:
PDF, 2,32 MB
Télécharger (pdf, 2,32 MB)
0 comments
 

To post a review, please sign in or sign up
Vous pouvez laisser un commentaire et partager votre expérience. D'autres lecteurs seront intéressés de connaitre votre opinion sur les livres lus. Qu'un livre vous plaise ou non, si vous partagez honnêtement votre opinion à son sujet, les autres pourront découvrir de nouveaux livres qui pourraient les intéresser.
1

Sex Is Still A Laughing Matter

Année:
1972
Langue:
english
Fichier:
PDF, 19,57 MB
0 / 0
2

Electrical Installation Guide 2016

Année:
2015
Langue:
english
Fichier:
PDF, 19,14 MB
0 / 0
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 87-89179-46-3
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 Three-mass model ........................................................................................25
4.1.3.2 Two-mass model ..........................................................................................26
4.1.3.3 One-mass 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 Steady-state model........................................................................................38
4.2.1.5 Modelling deep-bar 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 quadrature-phase stator windings ....................................49
Commutator model ...................................................................................................50
Electromagnetic torque.............................................................................................52
4.2.2.2 Simulink model ............................................................................................52
4.2.3
Permanent-magnet synchronous machine ........................................................53
3

44
Wind Turbine Blockset in Matlab Simulink

4.3
Power converter models........................................................................................... 56
4.3.1
AC Controllers (Soft-starters).......................................................................... 56
4.3.1.1 Star connected 3-phase load......................................................................... 58
4.3.1.2 Delta connected 3-phase load ...................................................................... 60
4.3.1.3 Branch-delta connected load........................................................................ 61
4.3.1.4 Simulink implementation of the soft-starter ................................................ 62
4.3.1.5 RMS model for soft-starter .......................................................................... 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
DC-link Circuit ................................................................................................ 68
4.4
Modulation strategies............................................................................................... 68
4.4.1
PWM sinusoidal modulator using sine-triangle 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 3-phase 2-winding 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 three-phase 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
Stator-flux oriented control.............................................................................................. 84
4.7
Summary .................................................................................................................. 87
Chapter 5 Simulation results................................................................................................. 88
5.1
Connection of a no-load transformer to the grid...................................................... 88
5.2
Disconnection of a no-load transformer................................................................... 89
5.3
Start-up sequence of a fixed-speed 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” S-Function ..................................................................... 101
Appendix B Other developed models ................................................................................... 105
Harmonic voltage source and grid angle detection............................................................ 105
Three-phase 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 #ENS-1363/01-0013, 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 Hardware-In-the-Loop 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 high-performance language for technical computing. It integrates
computation, visualization, and programming in an easy-to-use 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 high-level matrix/array language with control flow
statements, functions, data structures, input/output, and object-oriented programming
features. It allows both "programming in the small" to rapidly create quick and dirty
throw-away 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 high-level
functions for two-dimensional and three-dimensional data visualization, image
processing, animation, and presentation graphics. It also includes low-level 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 MAT-files.
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 clickand-drag 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 S-Functions it is also possible to customize and create
user-defined blocks. Models are hierarchical, so the models can be built using both top-down
and bottom-up approaches. The system can be viewed at a high level, then double-click
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 command-line 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 S-function is a computer language description of a Simulink block. S-functions can be
written in MATLAB, C, C++, Ada, or Fortran. C, C++, Ada, and Fortran S-functions are
compiled as MEX-files using the mex utility [78]. As with other MEX-files, they are
dynamically linked into MATLAB when needed. S-functions 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 built-in Simulink blocks. The form of
an S-function is general and can accommodate continuous, discrete, and hybrid systems. Sfunctions allow the user to add new models to the Simulink built-in models. These new blocks
can be written in MATLAB®, C, C++, Fortran, or Ada. An algorithm or model can be
implemented in an S-function 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 S-Functions.
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, drive-trains and transformer) are also available as ‘C’ S-Functions for highspeed simulations;

•

In order to be able to use different drive-train 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 soft-starter fed induction machines, which are used
in fixed-speed 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 fixed-point 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 built-in 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’ S-Function 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 look-up 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. one-mass model, two-mass model of the drive
train with torsional torques. Since the last model involves 4 differential equations an
alternative implementation using C S-Function 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 squirrel-cage induction machine, wound rotor
induction machine, salient-poles 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 dqo-reference frames in
modelling the electrical machines [83]-[85]. However, some particular operating modes e.g.
soft-starter-fed induction machines requires an abc/abc model [67], [69], [73], [74].
Therefore, the models of the electrical machines are written in state-variable form in the dqoreference frame as well as in the abc/abc natural reference frame. Some special features, as
deep-bar 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 per-phase quantities. Therefore, the desired winding-connection can be obtained
e.g. star, delta or branch-delta as in the case of soft-starters-fed 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 deep-bar 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 S-Function. 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 S-Functions. An example of a C
S-Function for the dq-model of the induction machine is shown in Appendix A.

Fig. 3.4. Electric machinery C S-Function based models.

However, contrary to C S-Function 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: 3-phase diode bridge rectifier, voltage source
converter (VCS), soft-starters 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 DC-link circuit a back-to-back converter topology can
be obtained.
The connection type of the soft-starter 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, dead-time 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 EN61000-2-2 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 S-Function;

•

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 (on-off) 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 three-phase system in phase quantities to a
three-phase one but in line-to-line quantities have been implemented. In order to find the
magnitude and the phase for a three-phase 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 (mono-phase, dq-reference 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
wind-up structure has been implemented in Simulink. The Maximum Power Point Tracker
(MPPT) block is based on a look-up 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 fixed-speed 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 dq-rotor reference frame used
in Model A and Model B.

The only difference between Model A and Model B consists in the implementation of the
state-space equations with currents as state-variables.
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 S-function 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 S-function 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 pseudo-code, which summarizes the
execution of a dynamic system in Simulink.

Fig. 3.15. Pseudo-code 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 top-level 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 end-of-simulation 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 block-level functions
(e.g. Output());

•

Number of internal methods – total number of invocations of system-level 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 70-80 % 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 S-function 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 S-Function and for an Integrator block.

It has been considered the average time (self-time) required to execute Update, Output and
MinorDeriv functions for these blocks as a percent from the sum of the self-time for these
functions for the Integrator block. It is obvious that the total time required for all invocations
of functions related to the S-function block is smaller than the corresponding time for one
Integrator block even if this S-function implements five differential equations.
Therefore, using a S-function 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 S-Function.
So, in modelling of large dynamic systems, which involves a big number of Integrator blocks
it is desirable to use C S-Function 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 fixed-point 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 built-in
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’ S-Function 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 – built-in Simulink Block).

The 20 bin-width 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 – built-in 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 – built-in 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 look-up 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 passive-stall 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, cut-in and cut-out wind
speeds, as shown in Fig. 4.7.

Fig. 4.7. Mask interface for the wind model in Simulink.

Using a reduced look-up table in respect with the variation of the torque coefficient/power
coefficient with the pitch angle ( −10o ≤ θ ≤ 10o ) the model for the active-stall 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 three-mass model for a wind turbine drive train the two-mass model and the
one-mass model are derived.
4.1.3.1

Three-mass 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. Two-mass 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

Two-mass model

The three-mass model can be reduced to a two-mass 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 drive-train 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 two-mass 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 two-mass model of the drive train.

4.1.3.3
One-mass model
When the stiffness and the damping factor are neglected further reduction of the drive train
model can be obtained. The one-mass 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. One-mass 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 state-variables.
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 state-variable 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 self-inductance;

L r = Lσr + 1.5M r
L σs , Lσr
Ms , M r
M sr

(4.2.6)

- rotor self-inductance;
- stator and rotor leakage inductance;
- stator and rotor mutual inductance;
- stator-rotor 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
squirrel-cage induction machines, these inductances can normally be estimated from no-load
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 squirrel-cage 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 squirrel-cage 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 state-space equations (voltage equations), which describe the behaviour of the induction
machine, are time-varying (except when the rotor is at stand-still). A change of variables is
often used to reduce the complexity of these state-space 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 dqo-arbitrary 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 state-space 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 state-space 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 steady-state
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

Steady-state model

The voltage equations, which describe the balanced steady-state 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 steady-state 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 per-phase equivalent circuit shown in Fig. 4.18.

Fig. 4.18. Equivalent circuit for steady state operation of a symmetrical induction machine.

The steady-state 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 quadri-pole as shown in Fig. 4.19
where the complex impedance is a function of the slip.

Fig. 4.19. Equivalent per-phase quadri-pole for the induction machine used in the steady-state 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
steady-state 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 steady-state model, both for a squirrel-cage machine and
wound rotor one, is shown in Fig. 4.20.

Fig. 4.20. Simulink implementation for steady-state 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 deep-bar 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 deep-bar effect

In larger machines like MW generators rotor deep-bar 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 squirrel-cage induction machine with deep-bars. The
parts of each bar extending deeply into the rotor iron have higher leakage inductances than
those parts of the bar cross-section 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. Deep-bar effect for a squirrel cage induction machine.

Under running conditions, the slip is quite small and the frequency of the rotor currents is
only 1-2 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 deep-bar 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 data-sheet given for rated operating point

•

K r is the resistance deep-bar factor;

•

K x is the leakage reactance deep-bar 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 deep-bar factors can be determined based on rated voltage
and current, starting current and phase angle (from short-circuit 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 deep-bar 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 per-phase (from data sheet);

•

X σs is the stator leakage reactance per-phase (from data sheet);

•

R 'r is the rotor resistance per-phase for the rated operating point (from data sheet);

•

X 'σr is the rotor leakage reactance per-phase for the rated operating point (from data
sheet);

•

ϕsc is the phase angle at standstill (s=1): from short-circuit 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
deep-bar effect.

a)

b)

Fig. 4.22. Deep-bar effect for a 2 MW squirrel-cage 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 squirrel-cage 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 deep-bar 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 deep-bar effect.

In order to include the deep-bar 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 non-linearity. 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 no-load curve for an induction machine these two inductances can be
determined as a function of the no-load current. Then, using a mathematical approximation
e.g. a function or a look-up 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 inverter-fed induction motor and take
into account the hysteresis losses as well as the eddy-current losses in a wide range of speed
and frequency. In wind turbine applications usually a large squirrel-cage 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 no-load 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 squirrel-cage 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 per-phase diagram, as
shown in Fig. 4.27.

Fig. 4.27. Equivalent diagram per-phase 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 self-inductance;
- rotor self-inductance;
- 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 salient-pole 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 non-salient pole machine assuming
identical values for inductances in d- and q-axis.
Untransformed model

In synchronous machines with salient poles the damping bars can be represented by a directand a quadrature-axis damping winding, which are short-circuited windings. Assume that on
the stator of the machine there is a symmetrical three-phase sinusoidal distributed stator
winding and on the rotor there are the field winding and direct- and quadrature-axis 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 three-phase 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 quadrature-phase stator windings

Replacing the three-phase stator winding by an equivalent quadrature-phase stator winding,
the number of differential equations can be reduced. Schematic of the corresponding model is
shown in Fig. 4.30.
o-axis
so

sq
q-axis
Q
sd
D

f

ωr
d-axis

Fig. 4.30. Schematic of the transformed model for salient-pole 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
d-q-o 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 salient-pole 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 d-q-o 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 


- direct-axis synchronous inductance

3 

L q =  L σs + M 0 − L 2 
2 


- quadrature-axis 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
quadrature-axis 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 direct-axis
2
3
L mq = M 0 − L 2 - the magnetizing inductances in quadrature-axis,
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 quadrature-axis 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 quadrature-axis
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 - self-inductances of the damper winding in the direct- and quadrature-axis;
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 homo-polar
component of the stator current is zero, so the homo-polar 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 q-axis as
well as for the magnetizing inductances. An extra input are the parameters for the field
winding.
4.2.3 Permanent-magnet 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 direct-axis (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.
o-axis
so

sq
q-axis
Q
sd
D

f

ωr
d-axis

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
quadrature-axis 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 soft-starter-fed 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 (Soft-starters)

Many authors deal with the theory of operation for three-phase AC-controllers (soft-starters),
[67], [69], and [73]. However, only few of them present in details the theory of operation with
a 3-phase resistive-inductive load [67] and [73]. Especially the influence of the load
connection type and the steady-state analysis are not presented in detail in the literature.
In [70] a comparison for full and half-wave controlled load both for star and delta connection
is carried out, unfortunately, the branch delta connection is not treated. An analysis of
variable-voltage 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 abc-dqo 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 converters-fed induction
machine, which involves a time-domain 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 delta-connected machine. Moreover these
approaches do not take into account the deep-bar effect, which is present into a large
induction machine.
56

5757
Wind Turbine Blockset in Matlab Simulink

Therefore in analysis of a soft-starter-fed induction machine is very convenient to use the
abc/abc model, since it is based on per-phase quantities.
At present many wind turbines, up to 2.3 MW, are based on the “Danish concept” in which a
squirrel-cage induction generator is directly connected to the grid as shown in Fig. 4.37.

Fig. 4.37. Block diagram of directly grid-connected 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 soft-starter, 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.
Soft-starters are used only during the start-up 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 soft-starters, 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 soft-starter-fed induction machine in wind turbine applications:
a) star connection, b) delta connection and c) branch-delta 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 on-state of the
SCRs on each phase. So, the soft-starter 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 soft-starter 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 soft-starter 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 soft-starter.
When a resistive-inductive 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 3-phase load

A soft-starter-fed 3-phase resistive-inductive load in star connection is shown in Fig. 4.40.

58

5959
Wind Turbine Blockset in Matlab Simulink

Fig. 4.40. Soft-starter-fed 3-phase R-L load in star connection.

Analyzing this topology the applied per-phase 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

to-ground supply voltages and SwA , SwB , SwC are the switching functions for each phase.
Based on the phase-to-ground 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 soft-starters 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 per-phase 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 3-phase 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, non-sinusoidal, 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 3-phase load

Fig. 4.42 shows a three-phase resistive-inductive load in delta connection.

Fig. 4.42. Soft-starter-fed 3-phase R-L load in delta connection.

It can easily be demonstrated that the voltages applied to a delta-connected 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 line-to-line supply voltages (or phase voltages for the load).
T

For a delta connection the switching functions are determined based on phase-to-ground
voltages and line currents.
Considering the same load as in the precedent case the applied per-phase 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 3-phase 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, non-sinusoidal, 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

Branch-delta connected load

A soft-starter-fed 3-phase resistive-inductive load in branch-delta connection is presented in
Fig. 4.44.

Fig. 4.44. Soft-starter-fed 3-phase R-L load in branch-delta connection.

The output voltages of the soft-starter 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 branch-delta connection the switching functions are determined based on the line
voltages and line currents.
Again, a 3-phase resistive-inductive 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 3-phase branch delta connected load:
a) Mode 1 and b) Mode 3 of operation.

Using a branch-delta 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 soft-starter
The soft-starter 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 soft-starter.

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 soft-starter model in Simulink.

4.3.1.5

RMS model for soft-starter

In steady-state analysis of the soft-starter-fed 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 soft-starter 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 cut-off 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 R-L 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 R-L load as shown in Fig. 4.49

Fig. 4.49. Single-phase 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 phase-angle.
Equation (4.3.6) is valid only for a branch-delta connection. Based on the input voltage and
the phase current in complex form and a look-up table for the extinction angle the output
voltage per-phase 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 AC-controller
output voltage.
The Simulink implementation of (4.3.6) is based on a look-up 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 single-phase controller with a series R-L load.

Fig. 4.51. Simulink implementation of the RMS model for the soft-starter.

4.3.2 Voltage Source Converters

Currently, in wind turbine applications the back-to-back 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 back-to-back 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 bi-directional power flow
can be achieved.
A voltage source converter can be implemented in several ways: six-step, 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 (SV-PWM).
Independent on the implementation method the converter can be seen as a black box with
some input-output 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 input-output characteristics as a function of the switch