From 728c61da3bbe3ff09847fab5438d8d92a906285e Mon Sep 17 00:00:00 2001 From: kbrown1snl <66843135+kbrown1snl@users.noreply.github.com> Date: Mon, 27 Mar 2023 14:14:27 -0600 Subject: [PATCH] Update 20_active_wake_control.py --- Examples/20_active_wake_control.py | 74 +++++++++++++++++++----------- 1 file changed, 48 insertions(+), 26 deletions(-) diff --git a/Examples/20_active_wake_control.py b/Examples/20_active_wake_control.py index bbb64244..c67555d8 100644 --- a/Examples/20_active_wake_control.py +++ b/Examples/20_active_wake_control.py @@ -3,41 +3,46 @@ 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: @@ -45,35 +50,44 @@ 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) @@ -84,7 +98,8 @@ 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] @@ -92,33 +107,40 @@ -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]. -----------------------------------------------