Lecture note · Fig. 10 · control systems
Wind turbine control: catching energy from a wind you cannot choose
Numa Koudsie · for school students and first- and second-year undergraduates
01From wind to electricity
Hold a paper windmill out of a car window and it spins faster the faster you drive. A modern wind turbine is the same toy, scaled up until each blade is longer than a jumbo jet’s wing.
The turbine of this page has blades 40 m long. Their tips sweep a disc of radius \(R\) and area \(\pi R^2 = 5027\ \mathrm{m^2}\), about 70 % of a football pitch. At its leisurely-looking 15 turns per minute, an angular speed \(\omega = 1.57\ \mathrm{rad/s}\), the tips move at \(\omega R = 62.8\ \mathrm{m/s}\), over 220 km/h. The wind pushes on the blades and makes the rotor turn; a generator in the nacelle brakes the rotor and turns the work of that braking into electrical energy.
The interesting physics lies in the balance between the two torques. If the generator brakes too little, the rotor runs away. If it brakes too hard, the rotor slows until the blades stall and the wind slips past. The controller, a small computer reading the rotor speed many times a second, decides how hard to brake. It cannot choose the wind: the wind is a time-dependent drive set by someone else, exactly like the force in Fig. 9.
02The power in the wind
Count the energy that the air carries through the disc. In one second a slab of air of length \(v\) passes through the area \(A\): its mass is \(\rho A v\), with \(\rho\) the density of air, and each kilogram carries kinetic energy \(\tfrac12 v^2\). The product is the power of the wind through the disc,
\[ P_\text{wind} = \tfrac12\,\rho A v^3 . \tag{1} \]The cube is the most important fact in wind energy. Doubling the wind speed doubles the mass of air arriving each second and quadruples the energy of each kilogram: eight times the power. At \(v = 8\ \mathrm{m/s}\), with \(\rho = 1.225\ \mathrm{kg/m^3}\), the air brings \(1.58\ \mathrm{MW}\) through this rotor; at \(11\ \mathrm{m/s}\) it brings \(4.10\ \mathrm{MW}\).
A rotor takes a share \(c_p\), the power coefficient, of this flux; and since power is torque times angular speed, the wind’s torque on the rotor is that power divided by \(\omega\):
\[ P_\text{aero} = c_p\,\tfrac12\,\rho A v^3, \qquad \tau_\text{aero} = \frac{P_\text{aero}}{\omega} . \tag{2} \]03The Betz limit★
A rotor cannot take all of the wind’s energy: air that stopped completely would pile up in front of the disc and block the air behind it. Albert Betz worked out the best possible share in 1920 by replacing the rotor with an idealised disc that the air pushes against. Let the air arrive at speed \(v\), cross the disc at \(v(1-a)\) and leave in the far wake at \(v_w\). The number \(a\) measures how much the disc slows the air.
Bernoulli’s equation says that pressure plus \(\tfrac12\rho\) times the square of the speed stays the same along the flow, as long as nothing does work on the air. So it holds in front of the disc and behind it, but not across it, where the disc takes energy out. The air crosses the thin disc at one speed, \(v(1-a)\), since what flows in must flow out; only the pressure jumps there, from \(p^+\) just in front to \(p^-\) just behind. Far upstream and far downstream the pressure is that of the open air, \(p_0\). In front, \(p_0 + \tfrac12\rho v^2 = p^+ + \tfrac12\rho v^2(1-a)^2\); behind, \(p^- + \tfrac12\rho v^2(1-a)^2 = p_0 + \tfrac12\rho v_w^2\). Adding the two, \(p_0\) and the speed at the disc drop out, and what is left is the pressure jump, \(\Delta p = p^+ - p^- = \tfrac12\rho\,(v^2 - v_w^2)\).
The force on the disc, \(\Delta p\,A\), must equal the momentum the air loses per second, \(\rho A v(1-a)\,(v - v_w)\). Equating the two and dividing by \(v - v_w\) (recall \(v^2 - v_w^2 = (v - v_w)(v + v_w)\)) shows that the disc sees exactly the average of the two speeds, \(v(1-a) = \tfrac12(v + v_w)\), so \(v_w = v(1-2a)\). The power the disc extracts is force times the air’s speed through it:
\[ \begin{aligned} P &= \rho A v(1-a)\cdot 2av\cdot v(1-a) = 2\rho A v^3\,a(1-a)^2, \\ c_p &= \frac{P}{\tfrac12\rho A v^3} = 4a(1-a)^2 . \end{aligned} \tag{3} \]Setting \(\mathrm{d}c_p/\mathrm{d}a = 4(1-a)(1-3a) = 0\) gives \(a = 1/3\) (the other root, \(a = 1\), would stop the air at the disc and take nothing), and
\[ c_{p,\text{Betz}} = 4\cdot\tfrac13\left(\tfrac23\right)^2 = \frac{16}{27} \approx 0.593 . \tag{4} \]No rotor in open air can turn more than 59.3 % of the power flowing through its disc into shaft power. The best rotor slows the air to two thirds of its speed at the disc and to one third in the wake.
04Tip-speed ratio
How large \(c_p\) is depends mainly on one number, the tip-speed ratio \(\lambda = \omega R/v\), the speed of the blade tips in units of the wind speed. A rotor turning too slowly lets most of the air pass between its blades without touching them. A rotor turning too fast behaves almost like a solid wall, and the air flows around it. Behind both pictures is one mechanism: each blade meets the air at a speed and an angle set by the wind and the blade’s own motion together, and that angle of attack decides how much lift, and so how much torque, the blade takes from the flow. A piece of blade at distance \(r\) from the hub moves across the wind at \(\omega r\), so it meets the air coming from an angle whose tangent is \(v/(\omega r) = (R/r)/\lambda\): at the tip, at \(\lambda = 8.1\), the undisturbed wind comes from only 7° off the plane of the rotor, and the air the rotor has already slowed (Section 3) from about 6°. Change \(\omega\) and \(v\) together, keeping \(\lambda\), and every piece of blade meets the air at the same angle as before; that is why \(c_p\) depends on the ratio and not on \(\omega\) and \(v\) separately. Between the two extremes lies a best ratio.
This page uses an empirical fit to measured rotor curves due to Heier, with the blades at zero pitch,
\[ c_p(\lambda) = c_1\left(\frac{c_2}{\lambda_i} - c_4\right)e^{-c_5/\lambda_i} + c_6\,\lambda, \qquad \frac{1}{\lambda_i} = \frac{1}{\lambda} - 0.035, \tag{5} \]with \(c_1 = 0.5176\), \(c_2 = 116\), \(c_4 = 5\), \(c_5 = 21\) and \(c_6 = 0.0068\). Scanning \(\lambda\) in steps of 0.001, its maximum is \(c_{p,\max} = 0.4800\) at \(\lambda_\text{opt} = 8.100\), comfortably below the Betz limit (Figure 2). One widely copied reproduction of the formula writes the last term as \(c_6\lambda_i\); that version peaks at 0.502 at \(\lambda = 8.25\) and does not match the quoted optimum, so the page uses \(c_6\lambda\).
05The rotor’s equation of motion
Newton’s second law for a rotation says that the net torque on the rotor changes its angular speed, in proportion to the inverse of its moment of inertia \(J\). Two torques act: the wind’s, which drives, and the generator’s, which brakes:
\[ J\,\frac{\mathrm{d}\omega}{\mathrm{d}t} = \tau_\text{aero}(\omega, v(t)) - \tau_\text{gen} . \tag{6} \]For the rotor and drive train together, the page takes \(J = 4.0\times10^6\ \mathrm{kg\,m^2}\). Multiplying (6) by \(\omega\) turns it into an energy statement: the wind’s power \(\tau_\text{aero}\omega\) either goes into the generator, \(P_\text{gen} = \tau_\text{gen}\omega\), or changes the rotor’s kinetic energy \(\tfrac12 J\omega^2\). Integrated over the run,
\[ E_\text{aero} = E_\text{gen} + \Delta\!\left(\tfrac12 J\omega^2\right), \tag{7} \]which the page checks after every solve: the Energy balance readout.
Notice the two time dependences in (6). The wind \(v(t)\) depends on time through the weather alone. The generator torque is chosen by the controller from the measured speed \(\omega\). The first is a drive, as in Fig. 9; the second is feedback (Figure 3).
06Two controllers★
6.1Follow the wind: the optimal-torque law
To take the most energy, the rotor should always run at \(\lambda_\text{opt}\). The controller cannot measure the wind reliably at the rotor, but it can measure \(\omega\). If the rotor is at the optimum, then \(v = \omega R/\lambda_\text{opt}\), and (2) becomes a function of \(\omega\) alone:
\[ \begin{aligned} \tau_\text{aero}\big|_{\lambda_\text{opt}} &= \frac{\tfrac12\rho\pi R^2 c_{p,\max}}{\omega}\left(\frac{\omega R}{\lambda_\text{opt}}\right)^3 = K_\text{opt}\,\omega^2, \\ K_\text{opt} &= \frac{\rho\pi R^5 c_{p,\max}}{2\lambda_\text{opt}^3} . \end{aligned} \tag{8} \]The controller brakes with \(\tau_{\rm gen}=K_{\rm opt}\omega^2\), here \(K_{\rm opt}=1.78\times10^5\ \mathrm{N\,m\,s^2}\). Its intended equilibrium, the speed \(\omega_* = \lambda_\text{opt}v/R\), is locally stable at fixed wind. At a smooth maximum of \(c_p\), \(dc_p/d\lambda=0\), so a small change of \(\omega\) leaves the aerodynamic power \(P_\text{aero}\) unchanged to first order. The torque \(P_\text{aero}/\omega\) then simply falls like \(1/\omega\), with slope \(-P_\text{aero}/\omega_*^2 = -K_{\rm opt}\omega_*\), because \(P_\text{aero} = \tau_\text{aero}\omega_* = K_{\rm opt}\omega_*^3\) there. The brake’s slope is \(2K_{\rm opt}\omega_*\). The net slope is therefore \(-3K_{\rm opt}\omega_*\lt0\): a small overspeed produces a restoring braking torque, and a small underspeed the opposite.
A small speed error consequently obeys \(J\,\delta\dot\omega\approx-3K_{\rm opt}\omega_*\,\delta\omega\), with relaxation time \(J/(3K_{\rm opt}\omega_*)\approx4.62\) s at 8 m/s wind. That is the local picture. The global one follows from the net torque, which with \(v = \omega R/\lambda\) in (2) reads \(\tfrac12\rho\pi R^5\omega^2\left[c_p(\lambda)/\lambda^3 - c_{p,\max}/\lambda_\text{opt}^3\right]\). The function \(c_p(\lambda)/\lambda^3\) is not monotone (it rises between \(\lambda = 2.4\) and \(4.3\)), but it stays above its value at \(\lambda_\text{opt}\) for every \(\lambda\) below \(\lambda_\text{opt}\) and below it for every \(\lambda\) above, checked on a grid of step \(10^{-4}\). Above \(\lambda_\text{opt}\) this is automatic, because \(c_p \le c_{p,\max}\) while \(\lambda^3 > \lambda_\text{opt}^3\); below, it is a property of this fit. So in steady wind, with the torque below its limit, the rotor returns to \(\lambda_\text{opt}\) from any start at which it is already turning: too slow and the wind wins, too fast and the brake wins. (The code sets \(c_p = 0\) below \(\lambda = 0.05\), so the page gives no torque at a standstill; the fit was not made for such slow running and says nothing about starting.) What the argument does not promise is optimal capture in changing wind: the return takes a few relaxation times, and a gust shorter than that is over before the rotor has caught up.
6.2Hold speed: a PI controller
Sometimes a constant speed is wanted: for a generator tied to the grid’s frequency, for quiet operation, or to stay away from a tower resonance. The controller then compares \(\omega\) with a target \(\omega_\text{ref}\) and brakes harder the faster the rotor turns, adding a term that remembers past errors:
\[ \tau_\text{gen} = \tau_0 + K_p\,(\omega - \omega_\text{ref}) + K_i\int_0^t(\omega - \omega_\text{ref})\,\mathrm{d}t' . \tag{9} \]Here \(\tau_0 = \tau_\text{aero}(\omega_\text{ref}, v(0))\) is the torque that balances the wind at the start, so the run begins in equilibrium with nothing yet to remember. To choose the gains, linearise (6) around the target in steady wind. A small speed error \(\delta\omega = \omega - \omega_\text{ref}\) changes the aerodynamic torque by \((\partial\tau_\text{aero}/\partial\omega)\,\delta\omega\). Put this and (9) into (6): the constant torques cancel, because \(\tau_0\) balances the wind, and what is left is \(J\,\delta\dot\omega = (\partial\tau_\text{aero}/\partial\omega - K_p)\,\delta\omega - K_i\int_0^t\delta\omega\,\mathrm{d}t'\). Differentiating once in time removes the integral and gives the equation of a damped oscillator, the same as in Fig. 9:
\[ \begin{gathered} J\,\delta\ddot\omega + \left(K_p - \frac{\partial\tau_\text{aero}}{\partial\omega}\right)\delta\dot\omega + K_i\,\delta\omega = 0 \\ \Longleftrightarrow\quad \delta\ddot\omega + 2\zeta_c\omega_n\,\delta\dot\omega + \omega_n^2\,\delta\omega = 0 . \end{gathered} \tag{10} \]The proportional term reacts to the speed error now. A positive error means overspeed, so it asks for more generator braking; a negative error reduces braking. The integral term remembers an error that persists: even a small overspeed builds a larger correction if it lasts. The linearised equation (10) above fixes how quickly a small disturbance decays. Torque limits and changing wind can take the rotor outside that approximation.
Picking a natural frequency \(\omega_n = 0.6\ \mathrm{rad/s}\) and a damping ratio \(\zeta_c = 0.8\) fixes \(K_i = J\omega_n^2 = 1.44\times10^6\ \mathrm{N\,m}\) and \(K_p = 2\zeta_c\omega_n J + \partial\tau_\text{aero}/\partial\omega\). At 15 turns per minute in an 8 m/s wind, the aerodynamic slope is \(-2.48\times10^5\ \mathrm{N\,m\,s}\): the wind itself damps the rotor a little. That gives \(K_p = 3.59\times10^6\ \mathrm{N\,m\,s}\). An error then decays within a few seconds. Because \(\zeta_c\) is a little below 1 it does swing once past zero, but only by a few per cent: after the wind rises by 1.5 m/s within 2 s in the Step run, the rotor overshoots the target by 0.39 rpm and then dips below it by just 0.01 rpm. The torque is kept between 0 and \(1.2\ \mathrm{MN\,m}\), and the integral stops growing while the torque sits at a limit (anti-windup).
07Same wind, two answers
Figure 4 shows the page’s default wind, two minutes of turbulence around a mean of 8 m/s, built from sixteen sine waves between 0.005 and 0.5 Hz with amplitudes falling as \(f^{-1/3}\), and the two controllers’ responses to it.
The table collects all eight combinations of controller and wind, with the page’s default sliders (\(v_0 = 8\ \mathrm{m/s}\), size of the gust, the step or the turbulence \(\Delta v = 1.5\ \mathrm{m/s}\), target 15 rpm).
| wind | controller | energy in 120 s | mean power | rotor speed | captured share |
|---|---|---|---|---|---|
| steady | hold speed | 25.15 kWh | 754 kW | 15.0 rpm | 0.479 |
| steady | follow the wind | 25.22 kWh | 757 kW | 15.5 rpm | 0.480 |
| gust (IEC shape, 1.5 m/s) | hold speed | 25.22 kWh | 757 kW | 14.8–15.3 rpm | 0.478 |
| gust (IEC shape, 1.5 m/s) | follow the wind | 25.30 kWh | 759 kW | 15.2–16.1 rpm | 0.480 |
| step, +1.5 m/s from 30 s to 32 s | hold speed | 34.43 kWh | 1033 kW | 15.0–15.4 rpm | 0.437 |
| step, +1.5 m/s from 30 s to 32 s | follow the wind | 37.18 kWh | 1115 kW | 15.5–18.4 rpm | 0.479 |
| turbulent, σ = 1.5 m/s | hold speed | 24.73 kWh | 742 kW | 14.4–15.7 rpm | 0.443 |
| turbulent, σ = 1.5 m/s | follow the wind | 26.18 kWh | 785 kW | 12.5–19.8 rpm | 0.466 |
Two lessons stand out. In steady wind the two controllers do almost equally well, because 15 rpm is close to the optimum 15.5 rpm. When the wind changes, holding the speed costs energy: after the step to 9.5 m/s the rotor stays at 15 rpm, \(\lambda\) falls to 6.6 and \(c_p\) drops from 0.48 to 0.43, which pulls the captured share of the whole run down to 0.44. Real turbines combine the two. Below their rated wind speed they follow the wind; above it they hold speed and power constant by turning the blades out of the wind (pitch control), which this model does not include.
Choose Step, keep \(v_0 = 8\ \mathrm{m/s}\) and the step size at 1.5 m/s, and solve once with Hold speed and once with Follow the wind. Compare Energy, 120 s, and watch the operating point on the \(c_p\) curve after the step. Then compute the rotor speed that puts \(\lambda\) at 8.1 in a 9.5 m/s wind, set the target speed to the nearest slider value, solve again, and compare the energy with both earlier runs.
08Solve first, then play
As in Fig. 9, nothing is integrated while you watch. The wind is a known function of time, so pressing Solve fixes the whole run.
- Solve. RK4 advances the rotor speed \(\omega\) and the controller’s integral from \(t = 0\) to \(120\ \mathrm{s}\) in 12 000 steps of \(h = 10\ \mathrm{ms}\). The same steps carry the generator energy, the aerodynamic energy, the energy of the wind through the disc, the distance the air has travelled and the rotor angle \(\theta\). The run starts in balance: at the optimum speed when following the wind, or at the target speed with the generator torque equal to the wind’s. A step of 10 ms is small for this problem: the speed controller answers within seconds and the quickest wind wave has a period of 2 s, so each is cut into hundreds of steps, and RK4’s error shrinks like \(h^4\). Solving the defaults again with \(h = 5\ \mathrm{ms}\) changes the stored rotor speed by at most \(2\times10^{-11}\) rpm, and with \(h = 20\ \mathrm{ms}\) by 16 times as much, as a fourth-order method should. The method would lose its high order where the generator torque hits one of its limits, because the right-hand side has a corner there; none of the eight runs of section 7 reaches a limit.
- Store. Every second state goes into the dataset: 6 001 rows of wind, speed, angle, tip-speed ratio, \(c_p\), both torques and both powers. With the default settings the energy balance (7) holds to about \(10^{-12}\).
- Play. Each frame shows the stored row nearest to the playback time. The blades are drawn at the stored angle \(\theta\), so they turn exactly as fast as the solution says.
09On screen
The rotor is drawn as a disc seen at an angle, foreshortened to 50 % across, so that both its turning and the wind passing it can be seen; the burgundy blade marks the stored angle \(\theta\). The grey streaks are air: they move by the stored distance \(\int v\,\mathrm{d}t\), at 5 px per metre of air, with a constant pattern of 16 px dashes and 22 px gaps. Only the wind arrow grows with the wind speed. They pass through the rotor unchanged; the slowing of the air is not drawn. The table at the top right lists the stored values of the present row. Below, the wind, the rotor speed (with the target dashed) and the generator power are drawn over the full two minutes, dark up to the present sample and pale after it. Behind the generator power, a rose curve shows the power the wind gives the rotor, \(P_\text{aero} = \tau_\text{aero}\omega\); the gap between the two is the rotor storing or returning kinetic energy, as in (7). The second figure shows the \(c_p\) curve with the present operating point and its last three seconds as a trail: when following the wind it stays near the peak; when holding speed it slides along the curve with every gust.
10What the model leaves out
- No pitch control, so above about 11 m/s nothing limits the power except the generator’s torque limit; real turbines pitch their blades there.
- The wind is the same over the whole disc and arrives at once; real wind varies across 80 m of rotor and changes height by height.
- The power coefficient responds instantly to the tip-speed ratio. Real rotors have dynamic inflow and stall effects with delays of seconds.
- The drive train is rigid and lossless, and \(P_\text{gen}\) is mechanical power into the generator; electrical conversion loses a few per cent more.
- A natural next question for a project: hang a small pendulum, or a double pendulum, on each blade and let generators at its joints harvest its swinging. Energy conservation already settles one point: every joule the pendulums harvest comes out of the rotor, and so out of the main generator’s share. The open question is whether a small chaotic pendulum can supply useful power to sensors in the blade, and what that costs.
10.1Further reading
- T. Burton, N. Jenkins, D. Sharpe and E. Bossanyi, Wind Energy Handbook, 2nd edition, Wiley, 2011: momentum theory, the Betz limit and turbine control (chapters 3 and 8).
- J. F. Manwell, J. G. McGowan and A. L. Rogers, Wind Energy Explained, 2nd edition, Wiley, 2009.
- S. Heier, Grid Integration of Wind Energy Conversion Systems, Wiley, 1998: source of the power-coefficient fit (5); publisher’s description at the Library of Congress.
- K. Shirzadeh, H. Hangan and C. Crawford, “Experimental and numerical simulation of extreme operational conditions for horizontal axis wind turbines based on the IEC standard”, Wind Energy Science 5, 1755–1770 (2020), open access: equation (4) gives the gust shape used by the page’s Gust wind.