Wind turbine design

Wind energy is the fastest and most mature resource of renewable energy. The increasing complexity of wind turbines requires advanced modeling tools for model-based design. In this tutorial, we will implement a utility-scale wind turbine model with the corresponding benchmark control system using Collimator. We will build the model developed by national renewable energy laboratory (NREL)(Jonkman et. al., 2009) following up with some modifications and simplifications provided by researchers at Aalborg university Grunnet et. al., 2010.

We begin by importing the libraries we will use:

import numpy as np
import matplotlib.pyplot as plt
import control as ctrl
import collimator as C

In the first part of this tutorial, we will build a reduced-order wind turbine model which is more suitable for model-based control system design methods. A typical wind turbine system combines a set of interacting subsystems:

Aerodynamics: A submodel that converts wind inflow into thrust force and torque.
Flexible structures: The dynamics of tower and blades oscillations.
Drivetrain: The dynamics of the shaft and gear system that transfers power from rotor to generator.
Generator: The subsystem responsible for generating electricity from mechanical power. This subsystem also lumps the corresponding power electronics.
Pitch Actuator: The submodel for the hydraulic servo that adjusts the blade pitch angle with the augmented control loop.

The wind turbine subsystem is illustrated in the following figure:

1. Aerodynamics

The aerodynamic thrust force $F_T$ and aerodynamic torque $T_r$ generated by wind inflow is defined as:

$$F_T = \frac{1}{2}\rho\pi R^2 C_T(\lambda,\beta) v_r^2$$

$$T_r = \frac{1}{2}\rho\pi R^3 \frac{C_P(\lambda,\beta)}{\lambda} v_r^3$$

where $\rho=1.2231\ kg/m^3$ is the air density, $R=63\ m$ is the rotor radius, $\lambda$ is the tip speed ratio, and $\beta$ is the collective blade pitch angle.

The relative wind speed $v_r$ is the difference between the wind speed and the nacelle speed due to tower oscillations:

$$v_r = v-\dot{q}_t $$

where $v$ is the wind inflow speed, and $q_t$ is the tower oscillatory mode position.

The thrust $C_T$ and power $C_P$ coefficients curves are highly nonlinear functions and they are obtained as lookup tables through more specialized high-fidelty modeling softwares such as WT_perf and Aerodyn.

Here, we have uploaded and read the csv files to obtain the independent variables and coefficient matrices. In the following Python code, we load the csv files and plot a contour plot for each function.

beta = np.genfromtxt('nrel5mw_beta.csv', delimiter=',')
tsr = np.genfromtxt('nrel5mw_tsr.csv', delimiter=',')
cp = np.genfromtxt('nrel5mw_cp.csv', delimiter=',')
ct = np.genfromtxt('nrel5mw_ct.csv', delimiter=',')

BETA, TSR = np.meshgrid(beta, tsr)
fig = plt.figure(figsize=(12,6))
#plt.rc('text', usetex=True)
#plt.rc('font', family='serif')
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2)
ax1.set_title('Thrust Function (ct)')
ax2.set_title('Power Function (cp)')
ax1.set_xlabel('Pitch angle (deg)')
ax2.set_xlabel('Pitch angle (deg)')
ax1.set_ylabel('Tip speed ratio')
ax2.set_ylabel('Tip speed ratio')
ax1.set_xlim(right=30)
ax1.set_xlim(left=0)
ax2.set_xlim(right=30)
ax2.set_ylim(top = 20)
c1 = ax1.contour(BETA,TSR,ct.T, levels = np.arange(0,1.5,0.005), cmap='Spectral')
c2 = ax2.contour(BETA,TSR,cp.T, levels = np.arange(0,0.5,0.001), cmap='Spectral')
fig.colorbar(c1, ax=ax1)
fig.colorbar(c2, ax=ax2)
plt.show()

The aerodynamics subsystem is modelled in Collimator as the submodel below.

The thrust and power coefficients are implemented as lookup tables. However, the elements are too many to be implemented in lookup table block by hand. Therefore, we will create a function to write the matrices and the vectors to be copied and pasted within the block.

def print_array_as_code (my_array): 
    if my_array.ndim == 1:
        code = '['
        for i in range(my_array.size-1):
            code += '{:f}, '.format(my_array[i])
        code += '{:f}]'.format(my_array[-1])   
    elif  my_array.ndim == 2:
        code = '['
        (m,n) = my_array.shape
        for i in range(m):
            code += '['
            for j in range(n):
                code += '{:f}'.format(my_array[i,j])
                if (j != n-1):
                    code += ', '
                else:
                    code += '] '
            if (i != m-1):
                code += ', '
        code += '] '
    else:
        code = 'ERROR! Unsupported ndim'
    print(code)

print_array_as_code(beta)
print_array_as_code(tsr)
print_array_as_code(ct)
print_array_as_code(cp)

The above matrices for beta and tsr are inserted as the input x array, and input y array, respectively. Similarly, we copy and paste the ct and cp matrices within the corresponding blocks.

2. Flexible Structures

The tower dynamics is modelled by a second-order linear system with a damping ration $\zeta_t = 0.08$ and eigenfrequency $f_t = 0.3210\ Hz$:

$$m_t\ddot{q}_t + 2\zeta_t\left(2\pi f_t\right)m_t\dot{q}_t+\left(2\pi f_t\right)^2 m_t q_t=F_T $$

$$F_t = -\left(2\pi f_t\right)^2 m_t $$

$$M_t = F_t h_t $$

$$M_b = \frac{2}{3} R F_T $$

where $m_t=697462\ m$ is the tower mass, $h_t=87.6\ m$, $F_t$ is the tower deflection force and $M_t$ is the tower deflection moment.

The above equations can be converted into a state-space form as:

$$\begin{bmatrix}\ddot{q}_t \\ \dot{q}_t \end{bmatrix} = \begin{bmatrix}-2\zeta_t\left(2\pi f_t\right) & -\left(2\pi f_t\right)^2 \\ 1 & 0 \end{bmatrix} \begin{bmatrix}\dot{q}_t \\ q_t \end{bmatrix} + \begin{bmatrix}\frac{1}{m_t} \\ 0 \end{bmatrix} F_T $$

$$\begin{bmatrix}\dot{q}_t \\ F_t \end{bmatrix} = \begin{bmatrix}1 & 0 \\ 0 & -\left(2\pi f_t\right)^2 m_t \end{bmatrix}\begin{bmatrix}\dot{q}_t \\ q_t \end{bmatrix} $$

The flexible structures dynamics are implemented in Collimator as shown in the following block diagram:

Flexible structure dynamics block diagram

Regarding the state-space block, we will define the following state-space matrices to be implemented in the state-space block.

m_t = 697462
damp_t = 0.08
f_t = 0.3210 
A = np.array([[-2*damp_t*2*np.pi*f_t, -(2*np.pi*f_t)**2],[1, 0]])
B = np.array([1.0/m_t, 0])
C = np.array([[1, 0],[0, -((2*np.pi*f_t)**2)*m_t]])
D = np.array([0,0])

m_t = 697462
damp_t = 0.08
f_t = 0.3210 
A = np.array([[-2*damp_t*2*np.pi*f_t, -(2*np.pi*f_t)**2],[1, 0]])
B = np.array([1.0/m_t, 0])
C = np.array([[1, 0],[0, -((2*np.pi*f_t)**2)*m_t]])
D = np.array([0,0])

3. Drivetrain

The drivetrain is assumed to have a two-mass model with torsional flexibility. The rotor-side is considered flexible whereas the generator-side is considered rigid. The rotor-side torque is defined as:

$$J_r \dot{\omega}_r = T_r - K_r \theta- B_r \dot{\theta}$$

The generator-side torque is defined as:

$$J_g \dot{\omega}_g = \frac{1}{N} \left(K_r \theta + B_r \dot{\theta}\right) - M_g $$

where $\omega_r$ is the rotor-side shaft angular velocity, $J_r=35444067\ Nm/rad$ is the rotor inertia, $\omega_g$ is the generator-side shaft angular velocity, $J_g=534.116\ kg.m^2$ is the generator inertia, $M_g$ is the generator torque, $K_r= 867637000\ Nm/rad$ is the rotor spring, $B_r= 867637000\ Nm/(rad.s)$ is the rotor damping, $N=97$ is the gear ratio, and $\dot{\theta}$ is the angular torsional deflection in the shaft which is defined as follows.

$$\dot{\theta} = \omega_r-\frac{1}{N}\omega_g $$

The above drivetrain dynamics are implemented in Collimator as the following submodel.

The torsional deflection submodel contains the diagram for the angular torsional deflection in the shaft equation: $\dot{\theta} = \omega_r-\frac{1}{N}\omega_g $

Before we move on with the next submodel, we have to initialize the angular velocity integrators for the rotor-side $\omega_r=0.5155$ and generator-side $\omega_g=50$ to ensure that the simulation converges.

Rotor and generator angular velocity integrators

4. Generator

The specification for the NREL 5MW does not contain a definition for the generator model. Therefore, we will use here a simplified 1st-order model as the electrical dynamics are much faster than the mechanical dynamics. The generator torque is modelled by the following differential equation:

$$\tau_g \dot{M}_g = T_{ref}-M_g $$

where $T_{ref}$ is the reference torque commanded by the controller, $M_g$ is the generator torque, and $\tau_g=0.1\ s$ is the generator time constant.

The generated power is defined as:

$$P_g = \eta M_g \omega_g $$

where $\eta=94.4\%$ is the generator efficiency. The Collimator model diagram is illustrated in the following figure.

5. Pitch Actuator

At this point, we will not use the original pitch servo system model described in the NREL 5MW reference. Instead, we will use the modified model defined in [Grunnet et. al., 2010]. The open-loop pitch drive system is modelled as the following second-order transfer function:

$$\frac{\beta \left(s\right)}{U\left(s\right)} = \frac{1}{s\left(\tau_{pitch}s+1\right)}$$

where $\tau_{pitch}=0.05\ s$ is the pitch servo time constant. The pitch system is controlled by proportional controller:

$$u\left(t\right) = K_p\left(\beta_{ref}-\beta\right)$$

where $K_p = 10$ is the proportional gain, and $\beta_{ref}$ is the pitch angle reference commanded by the controller.

We can have a look at the closed-loop performance of the pitch servo:

s = ctrl.tf('s')
tau_pitch = 0.05
K_p = 10
PitchServo_ol = 1/s/(tau_pitch*s+1)
PitchServo_cl = ctrl.feedback(K_p*PitchServo_ol,1)
print(PitchServo_cl)
T, yout = ctrl.step_response(PitchServo_cl,T=np.arange(0, 1, 0.01))
fig = plt.figure(figsize=(12,8))
plt.grid(which='both')
plt.ylabel('Collective Pitch Angle (deg)')
plt.xlabel('Time (s)')
plt.plot(T,yout)
plt.show()
print(ctrl.step_info(PitchServo_cl))

{'RiseTime': 0.1465281422814393, 'SettlingTime': 0.425629365674657, 'SettlingMin': 0.9022199550997985, 'SettlingMax': 1.0432137926607123, 'Overshoot': 4.321379266071235, 'Undershoot': 0, 'Peak': 1.0432137926607123, 'PeakTime': 0.3139888763173699, 'SteadyStateValue': 1.0}

The pitch system is modelled in Collimator as shown in the following diagram:

6. Composition of Submodels

The Submodel blocks are connected together as in the following aggregate wind turbine model:

Wind Turbine Controller

Variable-speed variable-pitch utility-scale wind turbine are usually controlled through two control systems:

Generator torque controller: The role of this loop is to maximize power capture below the rated operation point (partial-load region). For NREL 5MW, the partial-load region is under the rated wind speed of $11.4\ m/s$.
Blade pitch angle controller: This loop tries to regulate the generator speed at the rated operation to keep an economic operation and reduces mechanical wearing. The full-load region for the NREL 5MW is above the rated wind speed of $11.4\ m/s$.

Both of the control systems are working independently. This control configuration is a conventional control strategy for wind turbines. For more details and elaboration about this strategy, please refer to (Pao and Johnson, 2009). In this tutorial, we will implement the baseline controller developed in the NREL 5MW wind turbine definition report (Jonkman et. al., 2009).

1. Generator Torque Controller

The generator torque command is computed as a piecewise function of the generator speed based on regions of operation.

Region 1: This region is below the cut-in wind speed (For NREL 5MW, $v_{cut in} = 3\ m/s$. Operation of wind turbine in this region has no economic feasibility. Therefore, the generator torque is turned off and no power is generated from wind.
Region 2: This region is the partial-load region where the generator torque is tracking the maximum power that can be extracted. The generator torque is proportional to the generator speed squared to track the optimal tip-speed ratio.
Region 3: This is full-load region is inversely proportional to the generator speed to track the rated power.
Regions 1.5 and 2.5 are transitional regions where the generator torque is defined as linear functions of generator speed.

For the NREL 5MW wind turbine, the generator torque command function is defined as:

$$ T_{ref}= \left\{\begin{array}{} 0 & 0<\omega_{g(rpm)}<670 \\ 96.5338\ \omega_{g(rpm)}-64677.65123 & 670\leq\omega_{g(rpm)}<871 \\ 0.025576386\ \omega_{g(rpm)}^2 & 871\leq\omega_{g(rpm)}<1136.4978 \\ 412.076\ \omega_{g(rpm)}-435288.3165 & 1136.4978\leq\omega_{g(rpm)}<1161.9632 \\ \frac{50578944.12852911}{\omega_{g(rpm)}} & \omega_{g(rpm)}\geq 1161.9632\end{array} \right. $$

In addition, the generator has the following physical constraints:

$$0 \leq T_g \leq 474502.91 $$

$$-15000 \leq \dot{T}_g \leq 15000.91$$

Here, we can visualize the the torque command function:

def gen_torq_cmd(gen_speed_rpm):
    # Define Constants
    m1_5 = 96.5338         # Region 1.5 slope (Nm/rpm)
    c1_5 = -64677.65123    # Region 1.5 offset (Nm)
    K = 0.025576386        # Region 2 gain (Nm/rpm^2)
    m2_5= 412.076          # Region 2.5 slope (Nm/rpm)
    c2_5 = -435288.3165    # Region 2.5 offset (Nm)
    P0 = 5e6/0.944         # Region 3 raw po

    if (gen_speed_rpm >= 1161.9632): # Region 3  
        gen_torque_ref = 30/3.14159 * P0 / gen_speed_rpm
    elif (gen_speed_rpm >= 1136.4978): # Region 2.5
        gen_torque_ref = m2_5 * gen_speed_rpm + c2_5
    elif (gen_speed_rpm >= 871): # Region 2
        gen_torque_ref = K * gen_speed_rpm ** 2
    elif (gen_speed_rpm > 670): # Region 1.5
        gen_torque_ref = m1_5 * gen_speed_rpm + c1_5
    else: # Region 1
        gen_torque_ref = 0
    return gen_torque_ref

gen_speed_rpm = np.arange(0,1200,0.1)
vf = np.vectorize(gen_torq_cmd)
gen_torq_ref = vf(gen_speed_rpm)
fig = plt.figure(figsize=(12,8))
plt.grid(which='both')
plt.ylabel('Commanded generator torque (N.m)')
plt.xlabel('Generator speed (rpm)')
plt.plot(gen_speed_rpm,gen_torq_ref)
plt.show()

The torque controller can be modelled in Collimator directly through a Python script block:

if time == 0:
    # Define Constants
    m1_5 = 96.5338         # Region 1.5 slope (Nm/rpm)
    c1_5 = -64677.65123    # Region 1.5 offset (Nm)
    K = 0.025576386        # Region 2 gain (Nm/rpm^2)
    m2_5= 412.076          # Region 2.5 slope (Nm/rpm)
    c2_5 = -435288.3165    # Region 2.5 offset (Nm)
    P0 = 5e6/0.944         # Region 3 raw po

if (gen_speed_rpm >= 1161.9632) or (pitch_command >= 1): # Region 3  
    gen_torque_ref = 30/3.14159 * P0 / gen_speed_rpm
    region = 3 # Full-load
    
elif (gen_speed_rpm >= 1136.4978): # Region 2.5
    gen_torque_ref = m2_5 * gen_speed_rpm + c2_5
    region = 2.5 # Partial-load (Transitional)
    
elif (gen_speed_rpm >= 871): # Region 2
    gen_torque_ref = K * gen_speed_rpm ** 2
    region = 2 # Partial-load
    
elif (gen_speed_rpm > 670): # Region 1.5
    gen_torque_ref = m1_5 * gen_speed_rpm + c1_5
    region = 1.5 # Partial-load (Transitional)
    
else: # Region 1
    gen_torque_ref = 0
    region = 1

The physical constraints of the generator can be implemented by 'Saturation' and 'Rate Limiter' blocks.

2. Blade Pitch Controller

This control system regulates the rotor speed at the rated value by manipulating the blade pitch angle. The blade pitch angle command is computed by a gain-scheduled PI controller. The PI controller is defined as:

$$U(z) = K_p + K_i \frac{T_s}{z-1}$$

where $K_p = 0.196491203325120$, $K_i=0.084210515710766$, and $T_s = 0.01$ is the controller sampling time.

The control action is multiplied by a correction factor that is a function of the current blade pitch command:

$$C_f = \frac{1}{1+\frac{\beta_{ref}}{0.109997}} $$

also, the pitch servo system has the following physical constraints:

$$0 \leq \beta_{(deg)} \leq 90 $$

$$-8 \leq \dot{\beta}_{(deg)} \leq 8 $$

The PI controller is implemented in Collimator as the following diagram:

The complete gain-scheduled controller is modelled as:

Gain-Scheduled PI controller for the blade pitch angle

We can see that besides the gain scheduled PI, we have added a mechanism to switch off the pitch controller while the wind turbine is below region 3. The off switch mechanism also resets the PI integrator to have an anti wind up effect.

3. The Overall Controller

The combined control strategy that contains the generator torque controller and the pitch controller is implemented in the Model Editor diagram below:

Generator torque and pitch controller model

Wind Turbine Design Simulation Results

To test the performance of the controller, we will use three datasets for wind fields generated by specialized software's for modeling wind turbulence:

partial_load_windfield.csv: This dataset contains a windfield with mean speeds above the cut-in speed and below the rated speed. The dataset is generated by TurbSim software.
full_load_windfield.csv: This dataset contains a windfield with mean speeds above the rated speed and below the cut-off speed. The dataset is generated by TurbSim software.
windfield.csv: This dataset combines both regions and it is generated by SimWindFarm toolbox. The wind field data set is inserted into a data source block set as below: