Rotor interaction

# install h2lib if not present
try:
    import h2lib
except ModuleNotFoundError:
    !pip install h2lib[test]

Import modules needed for this tutorial

from h2lib import H2Lib
from h2lib_tests import tfp
import matplotlib.pyplot as plt
import numpy as np
from tqdm.notebook import tqdm
from IPython import display
model_path = tfp + 'DTU_10_MW/'
htc_path = 'htc/DTU_10MW_RWT.htc'

Diameter, position and orientation

with H2Lib() as h2:
    h2.init(htc_path, model_path)
    D = h2.get_diameter()
    rotor_pos = h2.get_rotor_position()
    rotor_ori = h2.get_rotor_orientation()
    aerosection_pos = h2.get_aerosections_position()

ax1, ax2 = plt.subplots(1,2, figsize=(12,6))[1]
x,y,z = rotor_pos
yaw, tilt, azi = rotor_ori

rotor = plt.Circle((x,-z), D/2,  fill=False)
ax1.plot(x,-z,'.k')
ax1.add_patch(rotor)

ax2.plot(y,-z,'.k')
ax2.plot([y-np.sin(tilt)*D/2, y+np.sin(tilt)*D/2], [-z-np.cos(tilt)*D/2, -z+np.cos(tilt)*D/2],'k')


for b in range(3):
    ax1.plot(aerosection_pos[b,:,0], -aerosection_pos[b,:,2],'.',ms=3, label=f'Blade {b}')
    ax2.plot(aerosection_pos[b,:,1], -aerosection_pos[b,:,2], label=f'Blade {b}')

ax1.axis('equal')
ax1.legend()
ax1.set_xlabel('x [m]')
ax2.set_ylabel('-z [m]')
ax2.axis('equal')
ax2.legend()
ax2.set_xlabel('y [m]')
ax2.set_ylabel('-z [m]')

Text(0, 0.5, '-z [m]')

DWM workflow

In a DWM (Dynamic wake meandering) workflow the wind speed including wakes but without induction is passed to HAWC2. Hawc2 calculates the induction which is passed back to the DWM which uses it to calculate the wakes.

Nxyz = (64,32,32)
class DWM():
    def __init__(self):
        self.time = 0

    def get_windfield(self):
        return np.random.random((3,) + Nxyz) # 3 components

    def set_induction(self, induction):
        pass

    def step(self):
        self.time += .2
        return self.time


dwm = DWM()

with H2Lib() as h2:
    h2.init(htc_path, model_path)
    h2.init_windfield(Nxyz=Nxyz, dxyz=(1,1,1), box_offset_yz=(-16,0), transport_speed=10)

    for _ in tqdm(range(10)):
        t = dwm.step()
        h2.run(t)
        uvw = dwm.get_windfield()
        h2.set_windfield(uvw, box_offset_x=-60)
        induc = h2.get_induction_axisymmetric()
        # induc = h2.get_induction_polargrid() # alternative if the full induction grid is needed
        dwm.set_induction(induc)

Actuator line workflow

class Flow():
    def __init__(self):
        self.time = 0

    def get_uvw(self, pos_xyz_lst):
        uvw = np.array(pos_xyz_lst) * 0
        uvw[:, :, 0] = 6
        return uvw.tolist()

    def step(self):
        self.time += .2
        return self.time

    def set_fxyz(self, pos_xyz, fxyz):
        pass


flow = Flow()
with H2Lib() as h2:
    h2.init(htc_path, model_path)
    for _ in tqdm(range(10)):
        t = flow.step()
        h2.run(t)
        pos_gl_xyz = h2.get_aerosections_position()
        uvw = flow.get_uvw(pos_gl_xyz)
        h2.set_aerosections_windspeed(uvw)
        frc_gl_xyz = h2.get_aerosections_forces()
        flow.set_fxyz(pos_gl_xyz, frc_gl_xyz)

Fully resolved blade and Miras workflow

These workflows have not been reimplemented in h2lib yet, and h2lib therefore still contains the old cpl routines.

To use these routines, the pyhawc2 package may be needed needed. This package has been updated to use the dll/so from h2lib.

Install via:

pip install pyhawc2 --index-url https://gitlab.windenergy.dtu.dk/api/v4/projects/796/packages/pypi/simple

Python controller workflow

In this example we will setup a workflow there the pitch demand value of blade 1 is updated from python with a sinusodial value.

Load htc file and find the sensor in htc.dll.type2_dll.output

from wetb.hawc2.htc_file import HTCFile
htc = HTCFile(model_path + htc_path)
servo_with_limits = htc.dll.get_subsection_by_name('servo_with_limits')
pitch1_sensor = servo_with_limits.output.sensors[1]

Modify pitch1 output from HAWC2 to servo_with_limits dll to use value of general variable with index 1 (initial value=0) instead of value from DTU controller

pitch1_sensor.type = "general"
pitch1_sensor.sensor = "variable"
pitch1_sensor.values = [1, 0]
pitch1_sensor.comments = ""
print ('Pitch sensor after:', pitch1_sensor)

Pitch sensor after: general variable    1 0;

Print full dll section, note sensor two in output

print (servo_with_limits)

begin type2_dll;
  name  servo_with_limits;
  filename      ./control/servo_with_limits.dll;
  dll_subroutine_init   init_servo_with_limits;
  dll_subroutine_update update_servo_with_limits;
  arraysizes_init       100 1;
  arraysizes_update     100 100;
  begin init;
    constant    1 3;    Number of blades [-]
    constant    2 1;    Frequency of 2nd order servo model of pitch system [Hz]
    constant    3 0.7;  Damping ratio 2nd order servo model of pitch system [-]
    constant    4 10;   Max. pitch speed [deg/s]
    constant    5 15;   Max. pitch acceleration [deg/s^2]
    constant    6 -5;   Min. pitch angle [deg]
    constant    7 90;   Max. pitch angle [deg]
    constant    8 -1;   Time for pitch runaway [s]
    constant    9 -1;   Time for stuck blade 1 [s]
    constant    10 0;   Angle of stuck blade 1 [deg] (if > 90 deg then blade is stuck at instantaneous angle)
  end init;
  begin output;
    general time;       Time                         [s]
    general variable    1 0;
    dll inpvec  1 3;    Pitch2 demand angle          [rad]
    dll inpvec  1 4;    Pitch3 demand angle          [rad]
    dll inpvec  1 26;   Flag for emergency pitch stop         [0=off/1=on]
  end output;
;
  begin actions;
    constraint  bearing2 angle pitch1;  Angle pitch1 bearing    [rad]
    constraint  bearing2 angle pitch2;  Angle pitch2 bearing    [rad]
    constraint  bearing2 angle pitch3;  Angle pitch3 bearing    [rad]
  end actions;
end type2_dll;
;
;       --- DLL for tower-blade tip distance -- ;
;-------------------------------------------

rename htc and save

htc.set_name('python_pitch_control')
htc.save()

with H2Lib() as h2:
    h2.init('htc/python_pitch_control.htc', model_path)

    # Add sensors
    h2.add_sensor('constraint bearing1 shaft_rot 1')  # rotor angle and speed
    h2.add_sensor('General variable 1 0 # $name(P1_demand) $unit(rad)')  # variable used to control the pitch servo
    h2.add_sensor('constraint bearing2 pitch1 1')  # angle and angle velocity of blade 1
    h2.add_sensor('constraint bearing2 pitch2 1')   # angle and angle velocity of blade 2

    # plot legends
    h2.step()
    c_lst = []
    sensor_idx = [1,3,4,6]
    for i in sensor_idx:
        name, unit, desc = h2.get_sensor_info(i)
        p = plt.plot([],'.', label=f'{name} [{unit}]')
        c_lst.append(p[0].get_color())
    plt.legend()


    for i in tqdm(range(25)):
        t = h2.get_time()
        # update pitch demand of blade 1
        h2.set_variable_sensor_value(1, np.sin(t)/10)

        h2.run(i/5)

        # plot sensors
        for value, c in zip(h2.get_sensor_values(sensor_idx), c_lst):
            plt.plot(t, value , '.', color=c)
        display.clear_output(wait=True)
        display.display(plt.gcf())