Skip to content

Commit

Permalink
Merge pull request #2 from kbrown1snl/patch-2
Browse files Browse the repository at this point in the history 
Update 20_active_wake_control.py
  • Loading branch information
@dzalkind
dzalkind authored Mar 27, 2023
2 parents 17ce1a5 + 728c61d commit 820aa66
Showing 1 changed file with 48 additions and 26 deletions.
74 changes: 48 additions & 26 deletions Examples/20_active_wake_control.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,77 +3,91 @@
Run openfast with ROSCO and active wake control
-----------------------------------------------
Set up and run simulation with AWC, check outputs
Active wake control (AWC) with blade pitching is implemented in this example with an adaptation into the rotating frame of the mathematical framework from the classical theory for stability of axisymmetric jets [1], which offers flexibility in specifying the forcing strategy.
Active wake control (AWC) with blade pitching is implemented in this example with two approaches as detailed below:
-----------------------------------------------
AWC_Mode = 1: Complex number method:
AWC_Mode = 1: Normal mode method:
-----------------------------------------------
The normal mode method is an adaptation into the rotating frame of the mathematical framework from the classical theory for stability of axisymmetric jets [1], which offers flexibility in specifying the forcing strategy.
The inputs to the controller are:
Name Unit Type Range Description
AWC_NumModes - Integer [0,inf] number of forcing modes
AWC_NumModes - Integer [0,inf] number of forcing modes
AWC_n - Integer [-inf,inf] azimuthal mode number(s) (i.e., the azimuthal mode number relates to the number and direction of the lobes of the wake structure according to the classical spatio-temporal Fourier decomposition of an arbitrary quantity q, sigma{sigma{q*exp(i*n*theta)*exp(i*omega*time)}}. For the case of a non-time-varying flow (i.e., where omega = 0), the azimuthal mode number specifies the number of cycles of blade pitch oscillation per one rotation around the rotor azimuth.)
AWC_clockangle deg Float [0,360] clocking angle(s) of forcing mode(s)
AWC_clockangle deg Float [0,360] clocking angle(s) of forcing mode(s)
AWC_freq Hz Float [0,inf] frequency(s) of forcing mode(s)
AWC_amp deg Float [0,inf] pitch amplitude(s) of forcing mode(s) (note that AWC_amp specifies the amplitude of each individual mode so that the total amplitude of pitching will be the sum of AWC_amp)
The latter two inputs may be specified based on the expected inflow while the former three inputs determine the type of active wake control to be used.
Readers may be familiar with several forcing strategies from literature on active wake control that can be represented as follows:
-collective dynamic induction control: AWC_NumModes = 1, AWC_n = 0, AWC_clockangle = 0
-helix clockwise [2]: AWC_NumModes = 1, AWC_n = 1, AWC_clockangle = 0
-helix counter-clockwise [2]: AWC_NumModes = 1, AWC_n = -1, AWC_clockangle = 0
-up-and-down: AWC_NumModes = 2, AWC_n = -1 1, AWC_clockangle = 0 0
-side-to-side: AWC_NumModes = 2, AWC_n = -1 1, AWC_clockangle = 90 90
-other: Higher-order modes or different combinations of the above can also be specified
Calculation methodology:
-helix clockwise [2]: AWC_NumModes = 1, AWC_n = 1, AWC_clockangle = 0
-helix counter-clockwise [2]: AWC_NumModes = 1, AWC_n = -1, AWC_clockangle = 0
-up-and-down: AWC_NumModes = 2, AWC_n = -1 1, AWC_clockangle = 0 0
-side-to-side: AWC_NumModes = 2, AWC_n = -1 1, AWC_clockangle = 90 90
-other: Higher-order modes or different combinations of the above can also be specified
These strategies are implemented using the following calculation methodology:
For each blade, we compute the total phase angle of blade pitch excursion according to:
AWC_angle(t) = 2*Pi*AWC_freq * t - AWC_n * (psi(t) + phi + AWC_clockangle*PI/180) (eq 1)
AWC_angle(t) = 2*Pi*AWC_freq * t - AWC_n * (psi(t) + phi + AWC_clockangle*PI/180) (eq 1)
where t is time
phi(t) is the angular offset of the given blade in the rotor plane relative to blade 1
psi is the angle of blade 1 in the rotor plane from top-dead center
Next, the phase angle is converted into the complex pitch amplitude:
AWC_complexangle(t) = AWC_amp*PI/180 * EXP(i * AWC_angle(t)) (eq 2)
AWC_complexangle(t) = AWC_amp*PI/180 * EXP(i * AWC_angle(t)) (eq 2)
where i is the square root of -1
Note that if AWC_NumModes>1, then eq 1 and 2 are computed for each additional mode, and AWC_complexangle becomes a summation over all modes for each blade.
Finally, the real pitch amplitude, theta(t), to be passed to the next step of the controller is calculated:
theta(t) = theta_0(t) + REAL(AWC_complexangle(t)) (eq 3)
theta(t) = theta_0(t) + REAL(AWC_complexangle(t)) (eq 3)
where theta_0(t) is the controller's nominal pitch command
Rearranging for ease of viewing:
Inserting eq 1 into eq 2, and then putting that result into eq 3 gives:
theta(t) = theta_0(t) + REAL(AWC_amp*PI/180 * EXP(i * (2*Pi*AWC_freq * t - AWC_n * (psi(t) + phi + AWC_clockangle*PI/180)))) (eq 4)
Applying Euler's formula and carrying out the REAL operator:
theta(t) = theta_0(t) + AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t - AWC_n * (psi(t) + phi + AWC_clockangle*PI/180)) (eq 5)
theta(t) = theta_0(t) + AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t - AWC_n * (psi(t) + phi + AWC_clockangle*PI/180)) (eq 5)
As an example, we can set parameters to produce the counter-clockwise helix pattern from [2] using AWC_NumModes = 1, AWC_n = -1, and AWC_clockangle = 0:
For blade 1, eq 5 becomes:
theta(t) = theta_0(t) + AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t + psi(t)) (eq 6)
theta(t) = theta_0(t) + AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t + psi(t)) (eq 6)
Note that the inverse multi-blade coordinate (MBC) transformation can also be used to obtain the same result as eq 6.
Beginning with Eq. 3 from [2], we have
/ \ / \
| theta_1(t) | | theta_0(t) |
| theta_2(t) | = T^-1(psi(t)) * | theta_tilt(t) | (eq 7)
| theta_2(t) | = T^-1(psi(t)) * | theta_tilt(t) | (eq 7)
| theta_3(t) | | theta_yaw(t) |
\ / \ /
where
/ \
| 1 cos(psi_1(t)) sin(psi_1(t)) |
T^-1(psi(t)) = | 1 cos(psi_2(t)) sin(psi_2(t)) |
| 1 cos(psi_3(t)) sin(psi_3(t)) |
\ /
Multiplying the first row of the top matrix (and dropping the subscript of blade 1) yields:
theta(t) = theta_0(t) + theta_tilt(t)*cos(psi(t)) + theta_yaw(t)*sin(psi(t)) (eq 8)
theta(t) = theta_0(t) + theta_tilt(t)*cos(psi(t)) + theta_yaw(t)*sin(psi(t)) (eq 8)
Setting theta_tilt(t) = AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t) and theta_yaw(t) = -AWC_amp*PI/180 * sin(2*Pi*AWC_freq * t) gives:
theta(t) = theta_0(t) + (AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t))*cos(psi(t)) - (AWC_amp*PI/180 * sin(2*Pi*AWC_freq * t))*sin(psi(t)) (eq 9)
theta(t) = theta_0(t) + (AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t))*cos(psi(t)) - (AWC_amp*PI/180 * sin(2*Pi*AWC_freq * t))*sin(psi(t)) (eq 9)
Applying a Ptolemy identity gives:
theta(t) = theta_0(t) + AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t + psi(t)) (eq 10)
theta(t) = theta_0(t) + AWC_amp*PI/180 * cos(2*Pi*AWC_freq * t + psi(t)) (eq 10)
which is equivlanet to eq 6 above.
-----------------------------------------------
AWC_Mode = 2: Coleman transform method:
-----------------------------------------------
A second method is the Coleman transform method.
The inputs to the controller are:
Name Unit Type Range Description
AWC_NumModes - Integer [1,2] number of modes for tilt and yaw (1: identical settings for tilt and yaw pitch angles, 2: seperate settings for tilt and yaw moments)
Expand All @@ -84,41 +98,49 @@
Using the inputs mentioned above, the user is able to specify any desired combination of sinusoidal tilt and yaw modes to be tracked by the turbine.
When a single mode is defined in the inputs, the prescribed tilt and yaw angles are assumed to be identical, except for the phase. The phase difference
between the tilt and yaw angles is taken from the input AWC_clockangle
between the tilt and yaw angles is taken from the input AWC_clockangle.
Readers may be familiar with several forcing strategies from literature on active wake control that can be represented as follows:
-collective dynamic induction control: AWC_NumModes = 1, AWC_harmonic = 0
-helix clockwise [2]: AWC_NumModes = 1, AWC_harmonic = 1, AWC_clockangle = -90 OR AWC_NumModes = 2, AWC_n = [1 1], AWC_clockangle = [0 -90]
-helix counter-clockwise [2]: AWC_NumModes = 1, AWC_harmonic = 1, AWC_clockangle = 90 OR AWC_NumModes = 2, AWC_n = [1 1], AWC_clockangle = [0 90]
-up-and-down: AWC_NumModes = 2, AWC_harmonic = [1 1], AWC_amp = [# 0] (where "#" represents the desired amplitude)
-side-to-side: AWC_NumModes = 2, AWC_harmonic = [1 1], AWC_amp = [0 #] (where "#" represents the desired amplitude)
-other: different combinations of the above can also be specified
Calculation methodology:
These strategies are implemented using the following calculation methodology:
The inputs described above enable the user to specify a desired sinusoidal signal for either the collective pitch (AWC_n = 0) or tilt and yaw pitch
angles (AWC_n = 1). These AWC pitch angles are defined as:
AWC_angle(t) = AWC_amp * sin(2*pi*AWC_freq*t + AWC_clockangle) (eq 1)
AWC_angle(t) = AWC_amp * sin(2*pi*AWC_freq*t + AWC_clockangle) (eq 1)
In case of collective pitch AWC, this signal is directly superimposed on the regular pitch control signal.
In case of IPC-based AWC, the reference tilt and yaw pitch angles theta are transformed to the rotating frame (i.e., pitch angles theta_k(t) for all
individual blades) using the inverse MBC transformation:
/ \ / \
| theta_1(t) | | theta_0(t) |
| theta_2(t) | = T^-1(psi(t)) * | theta_tilt(t) | (eq 2)
| theta_2(t) | = T^-1(psi(t)) * | theta_tilt(t) | (eq 2)
| theta_3(t) | | theta_yaw(t) |
\ / \ /
where
theta_tilt(t) = AWC_amp(1) * sin(2*pi*AWC_freq(1)*t + AWC_clockangle(1)) (eq 3)
theta_yaw(t) = AWC_amp(2) * sin(2*pi*AWC_freq(2)*t + AWC_clockangle(2)) (eq 4)
theta_tilt(t) = AWC_amp(1) * sin(2*pi*AWC_freq(1)*t + AWC_clockangle(1)) (eq 3)
theta_yaw(t) = AWC_amp(2) * sin(2*pi*AWC_freq(2)*t + AWC_clockangle(2)) (eq 4)
and
/ \
| 1 cos(psi_1(t)) sin(psi_1(t)) |
T^-1(psi(t)) = | 1 cos(psi_2(t)) sin(psi_2(t)) | (eq 5)
T^-1(psi(t)) = | 1 cos(psi_2(t)) sin(psi_2(t)) | (eq 5)
| 1 cos(psi_3(t)) sin(psi_3(t)) |
\ /
with psi_k(t) the azimuthal position of blade k at time instant t. Note that if AWC_NumModes = 1, it is assumed that:
AWC_amp(2) = AWC_amp(1)
AWC_freq(2) = AWC_freq(1)
AWC_clockangle(2) = 2*AWC_clockangle(1)
For more information on this control strategy, the user is referred to [2].
-----------------------------------------------
Expand Down

0 comments on commit 820aa66

Please sign in to comment.