diff --git a/wetb/gtsdf/__init__.py b/wetb/gtsdf/__init__.py index 56c191e1635e05f7d8f8d8fd8ef4057447ec566d..1c26d6b949b0f7248ca38e5f1c40d03d17807d99 100644 --- a/wetb/gtsdf/__init__.py +++ b/wetb/gtsdf/__init__.py @@ -44,6 +44,7 @@ from .gtsdf import compress2statistics class Dataset(object): def __init__(self, filename): + self.filename = filename self.time, self.data, self.info = load(filename) def __call__(self, id): if isinstance(id, str): diff --git a/wetb/hawc2/at_time_file.py b/wetb/hawc2/at_time_file.py index 3819ec2b37d0bf55d3aa7ff6e65ca258b2537fa6..165ad3b81c7aebad53635b44c4efba677283d2d5 100644 --- a/wetb/hawc2/at_time_file.py +++ b/wetb/hawc2/at_time_file.py @@ -49,8 +49,9 @@ class AtTimeFile(object): self.blade_radius = bladetip_radius with open(filename, encoding='utf-8') as fid: lines = fid.readlines() - atttribute_name_line = [l.startswith("# Radius_s") for l in lines].index(True) - self.attribute_names = lines[atttribute_name_line].lower().replace("#", "").split() + atttribute_name_line = [l.strip().startswith("# Radius_s") for l in lines].index(True) + #self.attribute_names = lines[atttribute_name_line].lower().replace("#", "").split() + self.attribute_names = [n.strip() for n in lines[atttribute_name_line].lower().split("#")[1:]] data = np.array([[float(l) for l in lines[i].split() ] for i in range(atttribute_name_line+1, len(lines))]) self.data = data def func_factory(column): @@ -97,10 +98,15 @@ class AtTimeFile(object): if __name__ == "__main__": at = AtTimeFile(r"tests/test_files/at.dat", 86.3655) # load file + at = AtTimeFile(r'U:\hama\HAWC2-AVATAR\res/avatar-7ntm-scaled-rad.dat') + at.attribute_names # Attribute names at[:3,1] # first 3 twist rows - at.twist()[:3] # Twist first 3 twist rows - print (at.twist(36, curved_length=True)) # Twist at curved_length = 36 (interpolated) - print (at.twist(36)) # Twist at 36 (interpolated) + print (len(at.attribute_names)) + print ("\n".join(at.attribute_names)) + print (at.data.shape) + #at.twist()[:3] # Twist first 3 twist rows + #print (at.twist(36, curved_length=True)) # Twist at curved_length = 36 (interpolated) + #print (at.twist(36)) # Twist at 36 (interpolated) diff --git a/wetb/hawc2/tests/test_files/htcfiles/DTU_10MW_RWT.htc b/wetb/hawc2/tests/test_files/htcfiles/DTU_10MW_RWT.htc new file mode 100644 index 0000000000000000000000000000000000000000..0c044026dffb65bb5bb2433d264e83c0a7359a26 --- /dev/null +++ b/wetb/hawc2/tests/test_files/htcfiles/DTU_10MW_RWT.htc @@ -0,0 +1,702 @@ +;DTU_10MW_RWT, version 9, 25-09-2017, mhha +; +begin simulation; + time_stop 1000.0; + solvertype 1 ; (newmark) + on_no_convergence continue ; + convergence_limits 1E3 1.0 1E-7 ; + logfile ./log/DTU_10MW_RWT_ver09.log ; + begin newmark; + deltat 0.01; + end newmark; +end simulation; +; +;---------------------------------------------------------------------------------------------------------------------------------------------------------------- +begin new_htc_structure; +; beam_output_file_name ./log/DTU_10MW_RWT_beam.dat; Optional - Calculated beam properties of the bodies are written to file +; body_output_file_name ./log/DTU_10MW_RWT_body.dat; Optional - Body initial position and orientation are written to file +; body_eigenanalysis_file_name ./eig/DTU_10MW_RWT_body_eigen.dat; +; structure_eigenanalysis_file_name ./eig/DTU_10MW_RWT_strc_eigen.dat ; +;------------------------------------------------------------------------------------------------------------------------------- +;------------------------------------------------------------------------------------------------------------------------------- + begin main_body; tower 115m + name tower ; + type timoschenko ; + nbodies 1 ; + node_distribution c2_def ; + damping_posdef 0.0 0.0 0.0 4.12E-03 4.12E-03 4.5E-04 ; Mx My Mz Kx Ky Kz , M´s raises overall level, K´s raises high freguency level "tuned by Larh" + begin timoschenko_input; + filename ./data/DTU_10MW_RWT_Tower_st.dat; + set 1 1 ; + end timoschenko_input; + begin c2_def; Definition of centerline (main_body coordinates) + nsec 11; + sec 1 0 0 0.00 0 ; x,y,z,twist + sec 2 0 0 -11.50 0 ; + sec 3 0 0 -23.00 0 ; + sec 4 0 0 -34.50 0 ; + sec 5 0 0 -46.00 0 ; + sec 6 0 0 -57.50 0 ; + sec 7 0 0 -69.00 0 ; + sec 8 0 0 -80.50 0 ; + sec 9 0 0 -92.00 0 ; + sec 10 0 0 -103.50 0 ; + sec 11 0 0 -115.63 0 ; + end c2_def ; + end main_body; +; + begin main_body; + name towertop ; + type timoschenko ; + nbodies 1 ; + node_distribution c2_def ; + damping_posdef 0.0 0.0 0.0 7.00E-03 7.00E-03 7.00E-03 ; "changed by Larh" + concentrated_mass 2.0 0.0 2.6870E+00 3.0061E-01 4.4604E+05 4.1060E+06 4.1060E+05 4.1060E+06 ; Nacelle mass and inertia "corrected by Anyd 25/4/13" + begin timoschenko_input; + filename ./data/DTU_10MW_RWT_Towertop_st.dat ; + set 1 2 ; + end timoschenko_input; + begin c2_def; Definition of centerline (main_body coordinates) + nsec 2; + sec 1 0.0 0.0 0.0 0.0 ; x,y,z,twist + sec 2 0.0 0.0 -2.75 0.0 ; + end c2_def ; + end main_body; +; + begin main_body; + name shaft ; + type timoschenko ; + nbodies 1 ; + node_distribution c2_def ; + damping_posdef 0.0 0.0 0.0 4.65E-04 4.65E-04 3.983E-03 ; "tuned by Anyd 23/5/13 to 31.45 log decr. damping for free free with stiff rotor and tower" + concentrated_mass 1.0 0.0 0.0 0.0 0.0 0.0 0.0 3.751E+06 ; generator equivalent slow shaft "re_tuned by Anyd 20/2/13" + concentrated_mass 5.0 0.0 0.0 0.0 1.0552E+05 0.0 0.0 3.257E+05 ; hub mass and inertia; "re_tuned by Anyd 20/2/13" + begin timoschenko_input; + filename ./data/DTU_10MW_RWT_Shaft_st.dat ; + set 1 1 ; + end timoschenko_input; + begin c2_def; Definition of centerline (main_body coordinates) + nsec 5; + sec 1 0.0 0.0 0.0 0.0 ; Tower top x,y,z,twist + sec 2 0.0 0.0 1.5 0.0 ; + sec 3 0.0 0.0 3.0 0.0 ; + sec 4 0.0 0.0 4.4 0.0 ; Main bearing + sec 5 0.0 0.0 7.1 0.0 ; Rotor centre + end c2_def ; + end main_body; +; + begin main_body; + name hub1 ; + type timoschenko ; + nbodies 1 ; + node_distribution c2_def ; + damping_posdef 0.0 0.0 0.0 3.00E-06 3.00E-06 2.00E-05; "changed by Larh" + begin timoschenko_input; + filename ./data/DTU_10MW_RWT_Hub_st.dat ; + set 1 2 ; + end timoschenko_input; + begin c2_def; Definition of centerline (main_body coordinates) + nsec 2; + sec 1 0.0 0.0 0.0 0.0 ; x,y,z,twist + sec 2 0.0 0.0 2.8 0.0 ; + end c2_def ; + end main_body; +; + begin main_body; + name hub2 ; + copy_main_body hub1; + end main_body; +; + begin main_body; + name hub3 ; + copy_main_body hub1 ; + end main_body; +; + begin main_body; + name blade1 ; + type timoschenko ; + nbodies 10 ; + node_distribution c2_def; + damping_posdef 0.0 0.0 0.0 1.53e-3 2.55e-3 3.3e-4 ; " 3% damping tuned by tkim 23/03/13 unable to fit 3rd and higher mode" + begin timoschenko_input ; + filename ./data/DTU_10MW_RWT_Blade_st.dat; + set 1 1 ; set subset + end timoschenko_input; + begin c2_def; Definition of centerline (main_body coordinates) + nsec 27 ; + sec 1 0.00000E+00 7.00600E-05 4.44089E-16 -1.45000E+01 ; + sec 2 -2.06477E-05 -1.22119E-02 3.00000E+00 -1.45000E+01 ; + sec 3 -7.28810E-03 -2.49251E-02 6.00000E+00 -1.44851E+01 ; + sec 4 -1.89235E-02 -2.73351E-02 7.00004E+00 -1.44610E+01 ; + sec 5 -5.41282E-02 -2.82163E-02 8.70051E+00 -1.43388E+01 ; + sec 6 -1.26633E-01 -2.13210E-02 1.04020E+01 -1.40201E+01 ; + sec 7 -2.25666E-01 -1.28378E-02 1.22046E+01 -1.33904E+01 ; + sec 8 -2.88563E-01 -7.70659E-03 1.32065E+01 -1.29371E+01 ; + sec 9 -3.99194E-01 -4.88317E-03 1.50100E+01 -1.19445E+01 ; + sec 10 -5.76634E-01 -1.80296E-02 1.82151E+01 -9.98243E+00 ; + sec 11 -7.07136E-01 -5.01772E-02 2.14178E+01 -8.45147E+00 ; + sec 12 -7.91081E-01 -9.41228E-02 2.46189E+01 -7.46417E+00 ; + sec 13 -8.37195E-01 -1.48880E-01 2.78193E+01 -6.72916E+00 ; + sec 14 -8.53948E-01 -2.14514E-01 3.10194E+01 -6.08842E+00 ; + sec 15 -8.49367E-01 -2.90618E-01 3.42197E+01 -5.49322E+00 ; + sec 16 -7.93920E-01 -4.62574E-01 4.02204E+01 -4.39222E+00 ; + sec 17 -7.16284E-01 -6.88437E-01 4.66217E+01 -3.09315E+00 ; + sec 18 -6.34358E-01 -9.60017E-01 5.30232E+01 -1.75629E+00 ; + sec 19 -5.53179E-01 -1.28424E+00 5.94245E+01 -5.00650E-01 ; + sec 20 -4.75422E-01 -1.66402E+00 6.58255E+01 6.01964E-01 ; + sec 21 -4.03180E-01 -2.10743E+00 7.22261E+01 1.55560E+00 ; + sec 22 -3.30085E-01 -2.65630E+00 7.90266E+01 2.51935E+00 ; + sec 23 -3.10140E-01 -2.78882E+00 8.05267E+01 2.72950E+00 ; + sec 24 -2.86719E-01 -2.92517E+00 8.20271E+01 2.93201E+00 ; + sec 25 -2.55823E-01 -3.06577E+00 8.35274E+01 3.11874E+00 ; + sec 26 -2.07891E-01 -3.20952E+00 8.50277E+01 3.28847E+00 ; + sec 27 -8.98940E-02 -3.33685E+00 8.63655E+01 3.42796E+00 ; + end c2_def ; + end main_body; +; + begin main_body; + name blade2 ; + copy_main_body blade1; + end main_body; +; + begin main_body; + name blade3 ; + copy_main_body blade1 ; + end main_body; +;------------------------------------------------------------------------------------------------------------------------------- +; + begin orientation; + begin base; + body tower; + inipos 0.0 0.0 0.0 ; initial position of node 1 + body_eulerang 0.0 0.0 0.0; + end base; +; + begin relative; + body1 tower last; + body2 towertop 1; + body2_eulerang 0.0 0.0 0.0; + end relative; +; + begin relative; + body1 towertop last; + body2 shaft 1; + body2_eulerang 90.0 0.0 0.0; + body2_eulerang 5.0 0.0 0.0; 5 deg tilt angle + mbdy2_ini_rotvec_d1 0.0 0.0 -1.0 0.2 ; mbdy2_ini_rotvec_d1 0.0 0.0 -1.0 [init_wr]; + end relative; +; + begin relative; + body1 shaft last; + body2 hub1 1; + body2_eulerang -90.0 0.0 0.0; + body2_eulerang 0.0 180.0 0.0; + body2_eulerang 2.5 0.0 0.0; 2.5deg cone angle + end relative; +; + begin relative; + body1 shaft last; + body2 hub2 1; + body2_eulerang -90.0 0.0 0.0; + body2_eulerang 0.0 60.0 0.0; + body2_eulerang 2.5 0.0 0.0; 2.5deg cone angle + end relative; +; + begin relative; + body1 shaft last; + body2 hub3 1; + body2_eulerang -90.0 0.0 0.0; + body2_eulerang 0.0 -60.0 0.0; + body2_eulerang 2.5 0.0 0.0; 2.5deg cone angle + end relative; +; + begin relative; + body1 hub1 last; + body2 blade1 1; + body2_eulerang 0.0 0.0 0; + end relative; +; + begin relative; + body1 hub2 last; + body2 blade2 1; + body2_eulerang 0.0 0.0 0.0; + end relative; +; + begin relative; + body1 hub3 last; + body2 blade3 1; + body2_eulerang 0.0 0.0 0.0; + end relative; +; + end orientation; +;------------------------------------------------------------------------------------------------------------------------------- +begin constraint; +; + begin fix0; fixed to ground in translation and rotation of node 1 + body tower; + end fix0; +; + begin fix1; + body1 tower last ; + body2 towertop 1; + end fix1; +; + begin bearing1; free bearing + name shaft_rot; + body1 towertop last; + body2 shaft 1; + bearing_vector 2 0.0 0.0 -1.0; x=coo (0=global.1=body1.2=body2) vector in body2 coordinates where the free rotation is present + end bearing1; +; + begin fix1; + body1 shaft last ; + body2 hub1 1; + end fix1; +; + begin fix1; + body1 shaft last ; + body2 hub2 1; + end fix1; +; + begin fix1; + body1 shaft last ; + body2 hub3 1; + end fix1; +; + begin bearing2; + name pitch1; + body1 hub1 last; + body2 blade1 1; + bearing_vector 2 0.0 0.0 -1.0; + end bearing2; +; + begin bearing2; + name pitch2; + body1 hub2 last; + body2 blade2 1; + bearing_vector 2 0.0 0.0 -1.0; + end bearing2; +; + begin bearing2; + name pitch3; + body1 hub3 last; + body2 blade3 1; + bearing_vector 2 0.0 0.0 -1.0; + end bearing2; +end constraint; +; +end new_htc_structure; +;---------------------------------------------------------------------------------------------------------------------------------------------------------------- +begin wind ; + density 1.225 ; + wsp 4.0 ; + tint 0.0 ; + horizontal_input 1 ; + windfield_rotations 0.0 0.0 0.0 ; yaw, tilt, rotation + center_pos0 0.0 0.0 -119 ; hub heigth + shear_format 1 0.2 ; + turb_format 0 ; 0=none, 1=mann,2=flex + tower_shadow_method 0 ; 0=none, 1=potential flow, 2=jet + scale_time_start 0.0 ; + wind_ramp_factor 0.0 40.0 0.6 1.0 ; +; Steps ; + wind_ramp_abs 140.0 141.0 0.0 1.0 ; wsp. after the step: 5.0 + wind_ramp_abs 181.0 182.0 0.0 1.0 ; wsp. after the step: 6.0 + wind_ramp_abs 222.0 223.0 0.0 1.0 ; wsp. after the step: 7.0 + wind_ramp_abs 263.0 264.0 0.0 1.0 ; wsp. after the step: 8.0 + wind_ramp_abs 304.0 305.0 0.0 1.0 ; wsp. after the step: 9.0 + wind_ramp_abs 345.0 346.0 0.0 1.0 ; wsp. after the step: 10.0 + wind_ramp_abs 386.0 387.0 0.0 1.0 ; wsp. after the step: 11.0 + wind_ramp_abs 427.0 428.0 0.0 1.0 ; wsp. after the step: 12.0 + wind_ramp_abs 468.0 469.0 0.0 1.0 ; wsp. after the step: 13.0 + wind_ramp_abs 509.0 510.0 0.0 1.0 ; wsp. after the step: 14.0 + wind_ramp_abs 550.0 551.0 0.0 1.0 ; wsp. after the step: 15.0 + wind_ramp_abs 591.0 592.0 0.0 1.0 ; wsp. after the step: 16.0 + wind_ramp_abs 632.0 633.0 0.0 1.0 ; wsp. after the step: 17.0 + wind_ramp_abs 673.0 674.0 0.0 1.0 ; wsp. after the step: 18.0 + wind_ramp_abs 714.0 715.0 0.0 1.0 ; wsp. after the step: 19.0 + wind_ramp_abs 755.0 756.0 0.0 1.0 ; wsp. after the step: 20.0 + wind_ramp_abs 796.0 797.0 0.0 1.0 ; wsp. after the step: 21.0 + wind_ramp_abs 837.0 838.0 0.0 1.0 ; wsp. after the step: 22.0 + wind_ramp_abs 878.0 879.0 0.0 1.0 ; wsp. after the step: 23.0 + wind_ramp_abs 919.0 920.0 0.0 1.0 ; wsp. after the step: 24.0 + wind_ramp_abs 960.0 961.0 0.0 1.0 ; wsp. after the step: 25.0 +; + begin tower_shadow_potential_2; + tower_mbdy_link tower; + nsec 2; + radius 0.0 4.15 ; + radius 115.63 2.75 ; + end tower_shadow_potential_2; +end wind; +; +begin aerodrag ; + begin aerodrag_element ; + mbdy_name tower; + aerodrag_sections uniform 10 ; + nsec 2 ; + sec 0.0 0.6 8.3 ; tower bottom + sec 115.63 0.6 5.5 ; tower top + end aerodrag_element; +; + begin aerodrag_element ; Nacelle drag side + mbdy_name shaft; + aerodrag_sections uniform 2 ; + nsec 2 ; + sec 0.0 0.8 10.0 ; + sec 7.01 0.8 10.0 ; + end aerodrag_element; +end aerodrag; +; +begin aero ; + nblades 3; + hub_vec shaft -3 ; rotor rotation vector (normally shaft composant directed from pressure to sustion side) + link 1 mbdy_c2_def blade1; + link 2 mbdy_c2_def blade2; + link 3 mbdy_c2_def blade3; + ae_filename ./data/DTU_10MW_RWT_ae.dat ; + pc_filename ./data/DTU_10MW_RWT_pc.dat ; + induction_method 1 ; 0=none, 1=normal + aerocalc_method 1 ; 0=ingen aerodynamic, 1=med aerodynamic + aerosections 50 ; def. 50 + ae_sets 1 1 1; + tiploss_method 1 ; 0=none, 1=prandtl + dynstall_method 2; 0=none, 1=stig øye method,2=mhh method +; +end aero ; +;------------------------------------------------------------------------------------------------- +begin dll; +; + begin type2_dll; + name dtu_we_controller ; + filename ./control/dtu_we_controller.dll ; + dll_subroutine_init init_regulation_advanced ; + dll_subroutine_update update_regulation ; + arraysizes_init 100 1 ; + arraysizes_update 100 100 ; + begin init ; + ; Overall parameters + constant 1 10000.0 ; Rated power [kW] + constant 2 0.628 ; Minimum rotor (LSS) speed [rad/s] + constant 3 1.005 ; Rated rotor (LSS) speed [rad/s] + constant 4 15.6E+06 ; Maximum allowable generator torque [Nm] + constant 5 100.0 ; Minimum pitch angle, theta_min [deg], + ; if |theta_min|>90, then a table of <wsp,theta_min> is read ; + ; from a file named 'wptable.n', where n=int(theta_min) + constant 6 82.0 ; Maximum pitch angle [deg] + constant 7 10.0 ; Maximum pitch velocity operation [deg/s] + constant 8 0.4 ; Frequency of generator speed filter [Hz] + constant 9 0.7 ; Damping ratio of speed filter [-] + constant 10 1.80 ; Frequency of free-free DT torsion mode [Hz], if zero no notch filter used + ; Partial load control parameters + constant 11 13013100.0 ; Optimal Cp tracking K factor [Nm/(rad/s)^2], ; + ; Qg=K*Omega^2, K=eta*0.5*rho*A*Cp_opt*R^3/lambda_opt^3 + constant 12 0.683456e8 ; Proportional gain of torque controller [Nm/(rad/s)] + constant 13 0.153367e8 ; Integral gain of torque controller [Nm/rad] + constant 14 0.0 ; Differential gain of torque controller [Nm/(rad/s^2)] +; Full load control parameters + constant 15 1 ; Generator control switch [1=constant power, 2=constant torque] + constant 16 1.06713 ; Proportional gain of pitch controller [rad/(rad/s)] + constant 17 0.242445 ; Integral gain of pitch controller [rad/rad] + constant 18 0.0 ; Differential gain of pitch controller [rad/(rad/s^2)] + constant 19 0.4e-8 ; Proportional power error gain [rad/W] + constant 20 0.4e-8 ; Integral power error gain [rad/(Ws)] + constant 21 11.4 ; Coefficient of linear term in aerodynamic gain scheduling, KK1 [deg] + constant 22 402.9 ; Coefficient of quadratic term in aerodynamic gain scheduling, KK2 [deg^2] & + ; (if zero, KK1 = pitch angle at double gain) + constant 23 1.3 ; Relative speed for double nonlinear gain [-] +; Cut-in simulation parameters + constant 24 -1 ; Cut-in time [s], no cut-in is simulated if zero or negative + constant 25 1.0 ; Time delay for soft start of torque [1/1P] +; Cut-out simulation parameters + constant 26 -1 ; Shut-down time [s], no shut-down is simulated if zero or negative + constant 27 5.0 ; Time of linear torque cut-out during a generator assisted stop [s] + constant 28 1 ; Stop type [1=normal, 2=emergency] + constant 29 1.0 ; Time delay for pitch stop after shut-down signal [s] + constant 30 3 ; Maximum pitch velocity during initial period of stop [deg/s] + constant 31 3.0 ; Time period of initial pitch stop phase [s] (maintains pitch speed specified in constant 30) + constant 32 4 ; Maximum pitch velocity during final phase of stop [deg/s] +; Expert parameters (keep default values unless otherwise given) + constant 33 2.0 ; Time for the maximum torque rate = Maximum allowable generator torque/(constant 33 + 0.01s) [s] + constant 34 2.0 ; Upper angle above lowest minimum pitch angle for switch [deg], if equal then hard switch + constant 35 95.0 ; Percentage of the rated speed when the torque limits are fully opened [%] + constant 36 2.0 ; Time constant of 1st order filter on wind speed used for minimum pitch [1/1P] + constant 37 1.0 ; Time constant of 1st order filter on pitch angle used for gain scheduling [1/1P] +; Drivetrain damper + constant 38 0.0 ; Proportional gain of active DT damper [Nm/(rad/s)], requires frequency in input 10 +; Over speed + constant 39 25.0 ; Overspeed percentage before initiating turbine controller alarm (shut-down) [%] +; Additional non-linear pitch control term (not used when all zero) + constant 40 0.0 ; Rotor speed error scaling factor [rad/s] + constant 41 0.0 ; Rotor acceleration error scaling factor [rad/s^2] + constant 42 0.0 ; Pitch rate gain [rad/s] +; Storm control command + constant 43 28.0 ; Wind speed 'Vstorm' above which derating of rotor speed is used [m/s] + constant 44 28.0 ; Cut-out wind speed (only used for derating of rotor speed in storm) [m/s] +; Safety system parameters + constant 45 30.0 ; Overspeed percentage before initiating safety system alarm (shut-down) [%] + constant 46 1.5 ; Max low-pass filtered tower top acceleration level [m/s^2] +; Turbine parameter + constant 47 178.0 ; Nominal rotor diameter [m] +; Parameters for rotor inertia reduction in variable speed region + constant 48 0.0 ; Proportional gain on rotor acceleration in variable speed region [Nm/(rad/s^2)] (not used when zero) +; Parameters for alternative partial load controller with PI regulated TSR tracking + constant 49 7.8 ; Optimal tip speed ratio [-] (only used when K=constant 11 = 0 otherwise Qg=K*Omega^2 is used) +; Parameters for adding aerodynamic drivetrain damping on gain scheduling + constant 50 0.0 ; Aerodynamic DT damping coefficient at the operational point of zero pitch angle [Nm/(rad/s)] (not used when zero) + constant 51 0.0 ; Coefficient of linear term in aerodynamic DT damping scheduling, KK1 [deg] + constant 52 0.0 ; Coefficient of quadratic term in aerodynamic DT damping scheduling, KK2 [deg^2] +; Torque exclusion zone + constant 53 0.0 ; Exclusion zone: Lower speed limit [rad/s] (Default 0 used if zero) + constant 54 0.0 ; Exclusion zone: Generator torque at lower limit [Nm] (Default 0 used if zero) + constant 55 0.0 ; Exclusion zone: Upper speed limit [rad/s] (if =< 0 then exclusion zone functionality is inactive) + constant 56 0.0 ; Exclusion zone: Generator torque at upper limit [Nm] (Default 0 used if zero) + constant 57 0.0 ; Time constant of reference switching at exclusion zone [s] (Default 0 used if zero) +; DT torsion mode damper + constant 58 0.0 ; Frequency of notch filter [Hz] (Default 10 x input 10 used if zero) + constant 59 0.0 ; Damping of BP filter [-] (Default 0.02 used if zero) + constant 60 0.0 ; Damping of notch filter [-] (Default 0.01 used if zero) + constant 61 0.0 ; Phase lag of damper [s] => max 40*dt (Default 0 used if zero) +; Fore-aft Tower mode damper + constant 62 0.0 ; Frequency of BP filter [Hz] (Default 10 used if zero)\\ + constant 63 0.0 ; Frequency of notch fiter [Hz] (Default 10 used if zero)\\ + constant 64 0.0 ; Damping of BP filter [-] (Default 0.02 used if zero)\\ + constant 65 0.0 ; Damping of notch filter [-] (Default 0.01 used if zero)\\ + constant 66 0.0 ; Gain of damper [-] (Default 0 used if zero)\\ + constant 67 0.0 ; Phase lag of damper [s] => max 40*dt (Default 0 used if zero)\\ + constant 68 0.0 ; Time constant of 1st order filter on PWR used for fore-aft Tower mode damper GS [Hz] (Default 10 used if zero) + constant 69 0.0 ; Lower PWR limit used for fore-aft Tower mode damper GS [-] (Default 0 used if zero) + constant 70 0.0 ; Upper PWR limit used for fore-aft Tower mode damper GS [-] (Default 0 used if zero) +; Side-to-side Tower mode filter + constant 71 0.0 ; Frequency of Tower side-to-sede notch filter [Hz] (Default 100 used if zero) + constant 72 0.0 ; Damping of notch filter [-] (Default 0.01 used if zero) + constant 73 0.0 ; Max low-pass filtered tower top acceleration level before initiating safety system alarm (shut-down) [m/s^2] (Default 1.1 x input 46 used if zero) + constant 74 0.0 ; Time constant of 1st order filter on tower top acceleration [1/1P] (Default 1 used if zero) +; Pitch deviation monitor parameters + constant 75 1005020 ; Parameters for pitch deviation monitoring. The format is 1,nnn,mmm + ; where 'nnn' [s] is the period of the moving average and 'mmm' is threshold of the deviation [0.1 deg] (functionality is inactive if value $<$ 1,000,000) +; Gear ratio + constant 76 0.0 ; Gear ratio used for the calculation of the LSS rotational speeds and the HSS generator torque reference [-] (Default 1 if zero) + end init ; +; + begin output ; + general time ; [s] + constraint bearing1 shaft_rot 1 only 2 ; Drivetrain speed [rad/s] + constraint bearing2 pitch1 1 only 1; [rad] + constraint bearing2 pitch2 1 only 1; [rad] + constraint bearing2 pitch3 1 only 1; [rad] + wind free_wind 1 0.0 0.0 -119 ; Global coordinates at hub height + dll inpvec 2 2 ; Elec. power from generator servo .dll + dll inpvec 2 8 ; Grid state flag from generator servo .dll + mbdy state acc towertop 1 1.0 global only 1 ; Tower top x-acceleration [m/s^2] + mbdy state acc towertop 1 1.0 global only 2 ; Tower top y-acceleration [m/s^2] + end output; + end type2_dll; +; + begin type2_dll; + name generator_servo ; + filename ./control/generator_servo.dll ; + dll_subroutine_init init_generator_servo ; + dll_subroutine_update update_generator_servo ; + arraysizes_init 100 1 ; + arraysizes_update 100 100 ; + begin init ; + constant 1 20.0 ; Frequency of 2nd order servo model of generator-converter system [Hz] + constant 2 0.9 ; Damping ratio 2nd order servo model of generator-converter system [-] + constant 3 15.6E+06 ; Maximum allowable LSS torque (pull-out torque) [Nm] + constant 4 0.94 ; Generator efficiency [-] + constant 5 1.0 ; Gearratio [-] + constant 6 0.0 ; Time for half value in softstart of torque [s] + constant 7 -1 ; Time for grid loss [s] (never if lower than zero) + end init ; +; + begin output; + general time ; Time [s] + dll inpvec 1 1 ; Electrical torque reference [Nm] + constraint bearing1 shaft_rot 1 only 2; Generator LSS speed [rad/s] + mbdy momentvec shaft 1 1 shaft only 3 ; Shaft moment [kNm] (Qshaft) + end output; +; + begin actions; + mbdy moment_int shaft 1 -3 shaft towertop 2 ; Generator LSS torque [Nm] + end actions; + end type2_dll; +; + begin type2_dll; + name mech_brake ; + filename ./control/mech_brake.dll ; + dll_subroutine_init init_mech_brake ; + dll_subroutine_update update_mech_brake ; + arraysizes_init 100 1 ; + arraysizes_update 100 100 ; + begin init ; + constant 1 9.36E+06 ; Fully deployed maximum brake torque [Nm] (0.6*max torque) + constant 2 100.0 ; Parameter alpha used in Q = tanh(omega*alpha), typically 1e2/Omega_nom + constant 3 0.5 ; Delay time for before brake starts to deploy [s] + constant 4 0.6 ; Time for brake to become fully deployed [s] + end init ; +; + begin output; + general time ; Time [s] + constraint bearing1 shaft_rot 1 only 2 ; Generator LSS speed [rad/s] + dll inpvec 1 25 ; Command to deploy mechanical disc brake [0,1] + end output; +; + begin actions; + mbdy moment_int shaft 1 -3 shaft towertop 2 ; Brake LSS torque [Nm] + end actions; + end type2_dll; +; + 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.0 ; 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.0 ; Max. pitch speed [deg/s] + constant 5 15.0 ; Max. pitch acceleration [deg/s^2] + constant 6 -5.0 ; Min. pitch angle [deg] + constant 7 90.0 ; Max. pitch angle [deg] + constant 8 -1 ; Time for pitch runaway [s] + constant 9 -1 ; Time for stuck blade 1 [s] + constant 10 0.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] + dll inpvec 1 2 ; Pitch1 demand angle [rad] + 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 -- ; +;------------------------------------------- +begin type2_dll; + name towerclearance_mblade ; + filename ./control/towerclearance_mblade.dll ; + dll_subroutine_init initialize ; + dll_subroutine_update update ; + arraysizes_init 3 1 ; + arraysizes_update 15 6 ; + begin init ; Variables passed into initialization function + constant 1 4.15 ; Tower radius at tower bottom [m] + constant 2 2.75 ; Tower radius at tower top [m] + constant 3 3 ; Number of points to check [-] + end init ; + begin output; Variables passed into update function + mbdy state pos tower 1 0.0 global ; [1,2,3] global coordinates of tower base + mbdy state pos tower 10 1.0 global ; [4,5,6] global coordinates of tower top + mbdy state pos blade1 26 1.0 global ; [7,8,9] global coordinates of point 1 (blade 1 tip) + mbdy state pos blade2 26 1.0 global ; [10,11,12] global coordinates of point 2 (blade 2 tip) + mbdy state pos blade3 26 1.0 global ; [13,14,15] global coordinates of point 3 (blade 3 tip) + end output; +end type2_dll; +; +end dll; +;---------------------------------------------------------------------------------------------------------------------------------------------------------------- +; +begin output; + filename ./res/DTU_10MW_RWT_ver09 ; + data_format hawc_binary; + buffer 1 ; +; + general time; + constraint bearing1 shaft_rot 2; angle and angle velocity + constraint bearing2 pitch1 5; angle and angle velocity + constraint bearing2 pitch2 5; angle and angle velocity + constraint bearing2 pitch3 5; angle and angle velocity + aero omega ; + aero torque; + aero power; + aero thrust; + wind free_wind 1 0.0 0.0 -119; local wind at fixed position: coo (1=global,2=non-rotation rotor coo.), pos x, pos y, pos z + ; Moments: + mbdy momentvec tower 1 1 tower # tower base ; + mbdy momentvec tower 10 2 tower # tower yaw bearing ; + mbdy momentvec shaft 4 1 shaft # main bearing ; + mbdy momentvec blade1 2 2 blade1 # blade 1 root ; + mbdy momentvec blade2 2 2 blade2 # blade 2 root ; + mbdy momentvec blade3 2 2 blade3 # blade 3 root ; + mbdy momentvec blade1 13 1 local # blade 1 50% local e coo ; + mbdy momentvec blade2 13 1 local # blade 2 50% local e coo ; + mbdy momentvec blade3 13 1 local # blade 3 50% local e coo ; + ; Displacements and accellerations + mbdy state pos tower 10 1.0 global only 1 # Tower top FA displ; + mbdy state pos tower 10 1.0 global only 2 # Tower top SS displ; + mbdy state acc tower 10 1.0 global only 1 # Tower top FA acc; + mbdy state acc tower 10 1.0 global only 2 # Tower top SS acc; +; + mbdy state pos blade1 26 1.0 blade1 # blade 1 tip pos ; + mbdy state pos blade2 26 1.0 blade2 # blade 2 tip pos ; + mbdy state pos blade3 26 1.0 blade3 # blade 3 tip pos ; + mbdy state pos blade1 26 1.0 global # gl blade 1 tip pos ; +; - Monitor Aerodynamics - ; + aero windspeed 3 1 1 72.5; + aero alfa 1 72.5; + aero alfa 2 72.5; + aero alfa 3 72.5; + aero cl 1 72.5; + aero cl 2 72.5; + aero cl 3 72.5; + aero cd 1 72.5; + aero cd 2 72.5; + aero cd 3 72.5; +; DLL outputs and into HAWC2 + dll inpvec 1 1 # Generator torque reference [Nm] ; + dll inpvec 1 2 # Pitch angle reference of blade 1 [rad] ; + dll inpvec 1 3 # Pitch angle reference of blade 2 [rad] ; + dll inpvec 1 4 # Pitch angle reference of blade 3 [rad] ; + dll inpvec 1 5 # Power reference [W] ; + dll inpvec 1 6 # Filtered wind speed [m/s] ; + dll inpvec 1 7 # Filtered rotor speed [rad/s]; + dll inpvec 1 8 # Filtered rotor speed error for torque [rad/s]; + dll inpvec 1 9 # Bandpass filtered rotor speed [rad/s]; + dll inpvec 1 10 # Proportional term of torque contr. [Nm] ; + dll inpvec 1 11 # Integral term of torque controller [Nm] ; + dll inpvec 1 12 # Minimum limit of torque [Nm] ; + dll inpvec 1 13 # Maximum limit of torque [Nm] ; + dll inpvec 1 14 # Torque limit switch based on pitch [-] ; + dll inpvec 1 15 # Filtered rotor speed error for pitch [rad/s]; + dll inpvec 1 16 # Power error for pitch [W] ; + dll inpvec 1 17 # Proportional term of pitch controller [rad] ; + dll inpvec 1 18 # Integral term of pitch controller [rad] ; + dll inpvec 1 19 # Minimum limit of pitch [rad] ; + dll inpvec 1 20 # Maximum limit of pitch [rad] ; + dll inpvec 1 21 # Torque reference from DT dammper [Nm] ; + dll inpvec 1 22 # Status signal [-] ; + dll inpvec 1 23 # Total added pitch rate [rad/s] ; + dll inpvec 1 24 # Filtered Mean pitch for gain sch [rad] ; + dll inpvec 1 25 # Flag for mechnical brake [0=off/1=on] ; + dll inpvec 1 26 # Flag for emergency pitch stop [0=off/1=on] ; + dll inpvec 1 27 # LP filtered acceleration level [m/s^2] ; + dll inpvec 1 31 # Monitored average of reference pitch [rad] ; + dll inpvec 1 32 # Monitored ave. of actual pitch (blade 1) [rad] ; +; Input from generator model + dll inpvec 2 1 # Mgen LSS [Nm] ; + dll inpvec 2 2 # Pelec [W] ; + dll inpvec 2 3 # Mframe [Nm] ; + dll inpvec 2 4 # Mgen HSS [Nm] ; + dll inpvec 2 8 # Grid flag [0=run/1=stop]; +; Input from mechanical brake + dll inpvec 3 1 # Brake torque [Nm] ; +; Input from pitch servo + dll inpvec 4 1 # pitch 1 [rad]; + dll inpvec 4 2 # pitch 2 [rad]; + dll inpvec 4 3 # pitch 3 [rad]; +; Check tower clearence + dll inpvec 5 1 # Bltip tow min d [m]; +end output; +; +exit; \ No newline at end of file diff --git a/wetb/signal/filters/cy_filters.py b/wetb/signal/filters/cy_filters.py index 35fe75dba9ac05fe3ebbbdea7c6757160670e55f..1d0b2e2073c08a7698aab1499adca1a744c77c3e 100644 --- a/wetb/signal/filters/cy_filters.py +++ b/wetb/signal/filters/cy_filters.py @@ -42,6 +42,33 @@ def cy_dynamic_low_pass_filter(inp, delta_t, tau, method=1): #cpdef cy_dynamic_ output[i] = output[i - 1] * np.exp(-delta_t / tau[i]) + inp[i] * (1 - np.exp(-delta_t / tau[i])) return output +def cy_dynamic_low_pass_filter_2d(inp, delta_t, tau, method=1): #cpdef cy_dynamic_low_pass_filter_2d(np.ndarray[double,ndim=2] inp, double delta_t, np.ndarray[double,ndim=1] tau, int method=1): + #cdef np.ndarray[double,ndim=2] output, alpha + #cdef int i + + output = np.empty_like(inp, dtype=np.float) + output[0] = inp[0] + + if method == 1: + alpha = delta_t / (tau + delta_t) + for i in range(1, inp.shape[0]): + output[i] = output[i - 1] + alpha[i] * (inp[i] - output[i - 1]) # Same as output[i] = alpha*inp[i]+(1-alpha)*output[i-1] + elif method == 2: + for i in range(1, inp.shape[0]): + output[i] = (delta_t * (inp[i] + inp[i - 1] - output[i - 1]) + 2 * tau[i] * output[i - 1]) / (delta_t + 2 * tau[i]) + elif method == 3: + for i in range(1, inp.shape[0]): + output[i] = output[i - 1] * np.exp(-delta_t / tau[i]) + inp[i] * (1 - np.exp(-delta_t / tau[i])) + return output + +def cy_dynamic_low_pass_filter_test(inp): #cpdef cy_dynamic_low_pass_filter_test(np.ndarray[double,ndim=2] inp): + #cdef np.ndarray[double,ndim=2] output, alpha + #cdef int i + output = np.empty_like(inp, dtype=np.float) + output[0] = inp[0] + for i in range(1, inp.shape[0]): + output[i] = inp[i] + return output @cython.ccall @cython.locals(alpha=cython.float, i=cython.int) diff --git a/wetb/signal/filters/first_order.py b/wetb/signal/filters/first_order.py index 7d10ba70706dbf17158bcecc8ad44c06b29ea21c..f5792ec9382da81617c9b537bfcd9ccb9c9a138e 100644 --- a/wetb/signal/filters/first_order.py +++ b/wetb/signal/filters/first_order.py @@ -11,7 +11,10 @@ def low_pass(input, delta_t, tau, method=1): if isinstance(tau, (int, float)): return cy_filters.cy_low_pass_filter(input.astype(np.float64), delta_t, tau) else: - return cy_filters.cy_dynamic_low_pass_filter(input.astype(np.float64), delta_t, tau, method) + if len(input.shape)==2: + return cy_filters.cy_dynamic_low_pass_filter_2d(input.astype(np.float64), delta_t, tau, method) + else: + return cy_filters.cy_dynamic_low_pass_filter(input.astype(np.float64), delta_t, tau, method) def high_pass(input, delta_t, tau): from wetb.signal.filters import cy_filters diff --git a/wetb/signal/subset_mean.py b/wetb/signal/subset_mean.py index 7578c11ac477856a237b7ec00ea0f2e418ba6695..e753e26e13a704fc8abea28af72b6fffc62de774 100644 --- a/wetb/signal/subset_mean.py +++ b/wetb/signal/subset_mean.py @@ -137,7 +137,7 @@ def cycle_trigger(values, trigger_value=None, step=1, ascending=True, tolerance= else: return np.where((values[1:] < trigger_value - tolerance) & (values[:-1] >= trigger_value + tolerance))[0][::step] -def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1): +def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None): """Returns one index per revolution (minimum rotor position) Parameters @@ -157,19 +157,25 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1): rotor_speed = np.ones_like(rotor_position)*rotor_speed deg_per_sample = rotor_speed*360/60/sample_frq thresshold = deg_per_sample.max()*3 - drp = np.diff(rotor_position)%360 + drp = (np.diff(rotor_position)+thresshold)%360-thresshold rp = rotor_position - rp = np.array(rotor_position).copy() + rp = np.array(rotor_position).copy()%360 #filter degree increase > thresshold rp[np.r_[False, np.diff(rp)>thresshold]] = 180 upper_indexes = np.where((rp[:-1]>(360-thresshold))&(rp[1:]<(360-thresshold)))[0] lower_indexes = np.where((rp[:-1]>thresshold)&(rp[1:]<thresshold))[0] +1 + if plt: + plt.plot(rp) + plt.plot(lower_indexes, rp[lower_indexes],'.') + plt.plot(upper_indexes, rp[upper_indexes],'.') + + # Best lower is the first lower after upper best_lower = lower_indexes[np.searchsorted(lower_indexes, upper_indexes)] upper2lower = best_lower - upper_indexes diff --git a/wetb/wind/turbulence/mann_parameters.py b/wetb/wind/turbulence/mann_parameters.py index cf33153ec81b21c3280882b66447d84ad3603efb..0dd3c3f4fd351386d4cd186ad5e1ef2711dc2a62 100644 --- a/wetb/wind/turbulence/mann_parameters.py +++ b/wetb/wind/turbulence/mann_parameters.py @@ -189,13 +189,13 @@ def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2): Returns ------- - residual : float + residual : array_like + rms of each spectrum """ - _3to2list = list(np.array(logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size))) - bk1, sp_meas = _3to2list[:1] + [_3to2list[1:]] - sp_fit = np.array(logbin_spectra(k1, *get_mann_model_spectra(ae, L, G, k1)))[1:] - - return np.sqrt(((bk1 * sp_meas - bk1 * sp_fit) ** 2).mean()) + k1_sp = np.array([sp for sp in logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size) if sp is not None]) + bk1, sp_meas = k1_sp[0], k1_sp[1:] + sp_fit = np.array(get_mann_model_spectra(ae, L, G, bk1))[:sp_meas.shape[0]] + return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1)) def fit_ae2var(variance, L, G): diff --git a/wetb/wind/turbulence/spectra.py b/wetb/wind/turbulence/spectra.py index 855ea98b31154af189ea6c7df45ba5047e7e89ba..f6aa6010d78df06c5bed4da3da878b3546b1528e 100644 --- a/wetb/wind/turbulence/spectra.py +++ b/wetb/wind/turbulence/spectra.py @@ -113,7 +113,7 @@ def spectra_from_time_series(sample_frq, Uvw_lst): sample_frq : int, float or array_like Sample frequency Uvw_lst : array_like - list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...] + list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...], v and w are optional Returns ------- @@ -127,7 +127,9 @@ def spectra_from_time_series(sample_frq, Uvw_lst): For two dimensional input, the mean spectra are returned """ assert isinstance(sample_frq, (int, float)) - U,v,w = [np.array(Uvw_lst)[:,i,:].T for i in range(3)] + Uvw_arr = np.array(Uvw_lst) + Ncomp = Uvw_arr.shape[1] + U,v,w = [Uvw_arr[:,i,:].T for i in range(Ncomp)] + [None]*(3-Ncomp) k = 2 * np.pi * sample_frq / U.mean(0) if v is not None: assert np.abs(np.nanmean(v,0)).max()<1, "Max absolute mean of v is %f"%np.abs(np.nanmean(v,0)).max() @@ -138,6 +140,15 @@ def spectra_from_time_series(sample_frq, Uvw_lst): return [k1_vec] + [spectrum(x1, x2, k=k) for x1, x2 in [(u, u), (v, v), (w, w), (w, u)]] + +{'ns':( 298, 0.0455033341592, 30.7642421893, 2.32693322102, 0.00538254487239), +'nu':( 853, 0.0355868134195, 71.0190181051, 2.95576644546, 0.00268903502954), +'s' :( 367, 0.0300454276184, 24.7375807569, 2.40165270878, 0.00253338624678), +'vs':( 223, 0.0149281609201, 11.7467239752, 3.21842435078, 0.00143453350246), +'n' :( 984, 0.0494970534618, 47.3846412548, 2.83466101838, 0.00513954022345), +'u' :( 696, 0.0258456660845, 78.9511607824, 2.51128584334, 0.00250014140782), +'vu':( 468, 0.0218274041889, 73.0730383576, 2.01636095317, 0.00196337613117)} + def bin_spectrum(x, y, bin_size, min_bin_count=2): assert min_bin_count > 0 x = x / bin_size