|
41 | 41 | "from scipy.linalg import block_diag\n",
|
42 | 42 | "import xarray as xr\n",
|
43 | 43 | "from math import comb\n",
|
44 |
| - "from wavespectra.construct.frequency import pierson_moskowitz\n", |
45 | 44 | "\n",
|
46 | 45 | "import wecopttool as wot\n",
|
47 | 46 | "\n",
|
|
76 | 75 | "fend = 1.875\n",
|
77 | 76 | "nfreq_irreg = 150\n",
|
78 | 77 | "f1_irreg = fend / nfreq_irreg\n",
|
79 |
| - "freq_irreg = wot.frequency(f1_irreg, nfreq_irreg, False) # False -> no zero frequency\n", |
80 | 78 | "\n",
|
81 | 79 | "f1_reg = .325/2\n",
|
82 |
| - "nfreq_reg = 12\n", |
83 |
| - "freq_reg = wot.frequency(f1_reg, nfreq_reg, False) # False -> no zero frequency" |
| 80 | + "nfreq_reg = 12" |
84 | 81 | ]
|
85 | 82 | },
|
86 | 83 | {
|
|
100 | 97 | "nrealizations = 2\n",
|
101 | 98 | "\n",
|
102 | 99 | "fp = 1/Tp\n",
|
103 |
| - "efth = pierson_moskowitz(freq=freq_irreg, hs=Hs, fp=fp)\n", |
| 100 | + "spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, Hs)\n", |
| 101 | + "efth = wot.waves.omnidirectional_spectrum(f1_irreg, nfreq_irreg, spectrum, \"Pierson-Moskowitz\")\n", |
104 | 102 | "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations=nrealizations)"
|
105 | 103 | ]
|
106 | 104 | },
|
|
280 | 278 | "outputs": [],
|
281 | 279 | "source": [
|
282 | 280 | "rho = 1025. # kg/m^3\n",
|
| 281 | + "freq_reg = wot.frequency(f1_reg, nfreq_reg, False) # False -> no zero frequency\n", |
283 | 282 | "bem_data_reg = wot.run_bem(pnr_fb, freq_reg)\n",
|
284 | 283 | "omega_reg = bem_data_reg.omega.values\n",
|
285 | 284 | "\n",
|
| 285 | + "freq_irreg = wot.frequency(f1_irreg, nfreq_irreg, False) # False -> no zero frequency\n", |
| 286 | + "bem_data_irreg = wot.run_bem(pnr_fb, freq_irreg)\n", |
| 287 | + "omega_irreg = bem_data_irreg.omega.values\n", |
| 288 | + "\n", |
286 | 289 | "bem_data_irreg = wot.run_bem(pnr_fb, freq_irreg)\n",
|
287 | 290 | "omega_irreg = bem_data_irreg.omega.values\n",
|
288 | 291 | "\n",
|
|
307 | 310 | "outputs": [],
|
308 | 311 | "source": [
|
309 | 312 | "# load in the empirical data\n",
|
310 |
| - "datafile = 'data/pioneer_empirical_data.npz'\n", |
311 |
| - "\n", |
312 |
| - "with np.load(datafile) as empirical_data:\n", |
313 |
| - " # Access the data within the 'with' block\n", |
314 |
| - " omega_data = empirical_data['omega_data']\n", |
315 |
| - " exc_coeff_data = empirical_data['exc_coeff_data']\n", |
316 |
| - " Zi_data = empirical_data['Zi_data']\n", |
317 |
| - " Zi_stiffness = empirical_data['Zi_stiffness']\n", |
318 |
| - "\n", |
319 |
| - "dof_names = [\"Pitch\"]\n", |
320 |
| - "directions = [0]\n", |
| 313 | + "datafile = 'data/pioneer_empirical_data.nc'\n", |
| 314 | + "empirical_data = xr.load_dataset(datafile)\n", |
321 | 315 | "\n",
|
322 |
| - "exc_coeff_data, Zi_data = wot.utilities.create_dataarray(Zi_data, exc_coeff_data, omega_data, directions, dof_names)\n", |
| 316 | + "exc_coeff_data = empirical_data.exc_coeff_data_real + 1j*empirical_data.exc_coeff_data_imag\n", |
| 317 | + "Zi_data = empirical_data.Zi_data_real + 1j*empirical_data.Zi_data_imag\n", |
| 318 | + "Zi_stiffness = empirical_data.Zi_stiffness\n", |
323 | 319 | "\n",
|
324 | 320 | "# here extrapolation for impedance and padding with zeros for the excitation\n",
|
325 | 321 | "exc_coeff_intrp_reg = exc_coeff_data.interp(omega = omega_reg, method='linear', kwargs={\"fill_value\": 0})\n",
|
|
503 | 499 | "pto_impedance = np.array([[pto_impedance_11, pto_impedance_12],\n",
|
504 | 500 | " [pto_impedance_21, pto_impedance_22]])\n",
|
505 | 501 | "\n",
|
506 |
| - "pto = wot.pto.PTO(ndof, \n", |
507 |
| - " kinematics=[1], \n", |
508 |
| - " controller=wot.controllers.unstructured_controller(), \n", |
509 |
| - " impedance=pto_impedance, \n", |
510 |
| - " loss=None, \n", |
511 |
| - " names=None)\n" |
| 502 | + "pto = wot.pto.PTO(ndof, kinematics=[1], controller=None, impedance=pto_impedance, loss=None, names=None)\n" |
512 | 503 | ]
|
513 | 504 | },
|
514 | 505 | {
|
|
856 | 847 | "outputs": [],
|
857 | 848 | "source": [
|
858 | 849 | "fig, ax = plt.subplots(nrows=5, sharex=True, figsize=(12, 12))\n",
|
859 |
| - "t = wec_tdom.sel(realization=0).time.values\n", |
| 850 | + "t = wec_tdom[0].time.values\n", |
860 | 851 | "\n",
|
861 | 852 | "# Positions\n",
|
862 | 853 | "ax[0].plot(t, fw_pos*180/np.pi, label='Flywheel', c='k')\n",
|
|
881 | 872 | "ax1 = ax[1].twinx()\n",
|
882 | 873 | "ax1.tick_params(axis='y', labelcolor='b')\n",
|
883 | 874 | "# when using an impedance model, the excitation force cannot be split up into Froude-Krylov and diffraction\n",
|
884 |
| - "wec_tdom.sel(realization=0)['force'].sel(type='excitation').plot(ax=ax1, c='b') \n", |
| 875 | + "wec_tdom[0]['force'].sel(type='excitation').plot(ax=ax1, c='b') \n", |
885 | 876 | "ax1.set_ylabel('Excitation torque [N-m]', color='b')\n",
|
886 | 877 | "ax1.set_title('Torque')\n",
|
887 | 878 | "\n",
|
|
937 | 928 | "Zi = Zi_intrp_reg\n",
|
938 | 929 | "Rad_res = np.real(Zi.squeeze())\n",
|
939 | 930 | "\n",
|
940 |
| - "Fex = wec_fdom.sel(realization=0).force.sel(type=['excitation'])\n", |
941 |
| - "Vel = wec_fdom.sel(realization=0).vel\n", |
| 931 | + "Fex = wec_fdom[0].force.sel(type=['excitation'])\n", |
| 932 | + "Vel = wec_fdom[0].vel\n", |
942 | 933 | "\n",
|
943 | 934 | "P_max_absorbable = (np.abs(Fex)**2/(8*Rad_res) ).squeeze().sum('omega').item() # after Falnes Eq. 6.10\n",
|
944 | 935 | "P_opt_excitation = 2*P_max_absorbable # after Falnes Eq. 6.10\n",
|
|
1043 | 1034 | "pto_impedance = np.array([[pto_impedance_11, pto_impedance_12],\n",
|
1044 | 1035 | " [pto_impedance_21, pto_impedance_22]])\n",
|
1045 | 1036 | "\n",
|
1046 |
| - "pto = wot.pto.PTO(ndof, \n", |
1047 |
| - " kinematics=[1], \n", |
1048 |
| - " controller=wot.controllers.unstructured_controller(), \n", |
1049 |
| - " impedance=pto_impedance, \n", |
1050 |
| - " loss=None, \n", |
1051 |
| - " names=None)\n", |
| 1037 | + "pto = wot.pto.PTO(ndof, kinematics=[1], controller=None, impedance=pto_impedance, loss=None, names=None)\n", |
1052 | 1038 | "\n",
|
1053 | 1039 | "wec_lin = wot.WEC.from_impedance(freq_irreg, \n",
|
1054 | 1040 | " Zi_intrp_irreg,\n",
|
|
1135 | 1121 | "outputs": [],
|
1136 | 1122 | "source": [
|
1137 | 1123 | "fig, ax = plt.subplots(nrows=5, sharex=True, figsize=(12, 12))\n",
|
1138 |
| - "t = wec_tdom.sel(realization=0).time.values\n", |
| 1124 | + "t = wec_tdom[0].time.values\n", |
1139 | 1125 | "\n",
|
1140 | 1126 | "# Positions\n",
|
1141 | 1127 | "ax[0].plot(t, fw_pos*180/np.pi, label='Flywheel', c='k')\n",
|
|
1159 | 1145 | "\n",
|
1160 | 1146 | "ax1 = ax[1].twinx()\n",
|
1161 | 1147 | "ax1.tick_params(axis='y', labelcolor='b')\n",
|
1162 |
| - "wec_tdom.sel(realization=0)['force'].sel(type=['excitation']).plot(ax=ax1, c='b')\n", |
| 1148 | + "wec_tdom[0]['force'].sel(type=['excitation']).plot(ax=ax1, c='b')\n", |
1163 | 1149 | "ax1.set_ylabel('Excitation torque [N-m]', color='b')\n",
|
1164 | 1150 | "ax1.set_title('Torque')\n",
|
1165 | 1151 | "\n",
|
|
0 commit comments