diff --git a/explore_wave.py b/explore_wave.py new file mode 100644 index 0000000..9d12d9b --- /dev/null +++ b/explore_wave.py @@ -0,0 +1,68 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 17 17:42:06 2024 + +@author: Yihan Liu +""" + +def pierson_moskowitz_spectrum(U19_5, zeta, eta, t, random_phases): + """ + This function generates the Pierson-Moskowitz spectrum for a given wind speed U10 and frequency f. + + parameters + ---------- + U19_5 : float + the average wind speed at 19.5m above the sea surface + zeta : float + the x component to evaluate + eta : float + the y component to evaluate. (Note: the coordinate system here is different + from the Betti model. The downward is negative + in this case) + t: float + the time to evaluate. + random_phase : Numpy Array + the random phase to generate wave. Should be in [0, 2*pi) + + Returns + ------- + wave_eta : float + The wave elevation + [v_x, v_y, a_x, a_y]: list + The wave velocity and acceleration in x and y direction + """ + g = 9.81 # gravitational constant + alpha = 0.0081 # Phillips' constant + + f_pm = 0.14*(g/U19_5) # peak frequency + + N = 400 + + cutof_f = 3*f_pm # Cutoff frequency + + f = np.linspace(0.1, cutof_f, N) # Array + omega = 2*np.pi*f # Array + delta_f = f[1] - f[0] # Array + + S_pm = (alpha*g**2/((2*np.pi)**4*f**5))*np.exp(-(5/4)*(f_pm/f)**4) # Array + + a = np.sqrt(2*S_pm*delta_f) + k = omega**2/g + + # Generate random phases all at once + + + # Perform the calculations in a vectorized manner + sin_component = np.sin(omega*t - k*zeta + random_phases) + cos_component = np.cos(omega*t - k*zeta + random_phases) + exp_component = np.exp(k*eta) + + wave_eta = np.sum(a * sin_component) + + v_x = np.sum(omega * a * exp_component * sin_component) + v_y = np.sum(omega * a * exp_component * cos_component) + + a_x = np.sum((omega**2) * a * exp_component * cos_component) + a_y = -np.sum((omega**2) * a * exp_component * sin_component) + + return wave_eta, [v_x, v_y, a_x, a_y] \ No newline at end of file diff --git a/gen_torque_controller.py b/gen_torque_controller.py index 6793928..e27577e 100644 --- a/gen_torque_controller.py +++ b/gen_torque_controller.py @@ -263,8 +263,8 @@ def pierson_moskowitz_spectrum(U19_5, zeta, eta, t, random_phases): a_x = np.sum((omega**2) * a * exp_component * cos_component) a_y = -np.sum((omega**2) * a * exp_component * sin_component) - #return wave_eta, [v_x, v_y, a_x, a_y] - return 0, [0,0,0,0] + return wave_eta, [v_x, v_y, a_x, a_y] + #return 0, [0,0,0,0] def structure(x_1, beta, omega_R, t, Cp_type, performance, v_w, v_aveg, random_phases): @@ -572,6 +572,9 @@ def WindTurbine(omega_R, v_in, beta, T_E, t, Cp): rho = 1.225 # (kg/m^3) Density of air A = 12469 # (m^2) Rotor area eta_G = 97 # (-) Speed ratio between high and low speed shafts + P_0 = 5296610 # rated power + + T_E = P_0 / (eta_G*omega_R) tildeJ_R = eta_G**2*J_G + J_R tildeT_E = eta_G*T_E @@ -583,6 +586,8 @@ def WindTurbine(omega_R, v_in, beta, T_E, t, Cp): T_A = P_A/omega_R domega_R = (1/tildeJ_R)*(T_A - tildeT_E) + #print(T_E*eta_G*omega_R) + return domega_R @@ -666,6 +671,7 @@ def rk4(Betti, x0, t0, tf, dt, beta_0, T_E, Cp_type, performance, v_w, v_wind, s """ d_BS = 37.550 # (m) The position of center of weight of BS (platform and tower) + eta_G = 97 # (-) Speed ratio between high and low speed shafts n = int((tf - t0) / dt) + 1 t = np.linspace(t0, tf, n) @@ -686,11 +692,12 @@ def rk4(Betti, x0, t0, tf, dt, beta_0, T_E, Cp_type, performance, v_w, v_wind, s def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i): - - eta_G = 97 # (-) Speed ratio between high and low speed shafts J_G = 534.116 # (kg*m^2) Total inertia of electric generator and high speed shaft J_R = 35444067 # (kg*m^2) Total inertia of blades, hub and low speed shaft tildeJ_R = eta_G**2*J_G + J_R + + + rated_omega_R = 1.26711 # The rated rotor speed is 12.1 rpm #rated_omega_R = 1.571 @@ -704,6 +711,9 @@ def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i): K_p = 0.0765*(2*tildeJ_R*rated_omega_R*zeta_phi*omega_phin*GK)/(eta_G*(-dpdbeta_0)) K_i = 0.013*(tildeJ_R*rated_omega_R*omega_phin**2*GK)/(eta_G*(-dpdbeta_0)) K_d = 0.187437 + #K_p = (2*tildeJ_R*rated_omega_R*zeta_phi*omega_phin*GK)/(eta_G*(-dpdbeta_0)) + #K_i = (tildeJ_R*rated_omega_R*omega_phin**2*GK)/(eta_G*(-dpdbeta_0)) + #K_d = 0 error_omega_R = omega_R - rated_omega_R error[i] = error_omega_R @@ -738,6 +748,7 @@ def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i): h_waves = [] for i in range(n - 1): betas.append(beta) + k1, h_wave, Q_t = Betti(x[i], t[i], beta, T_E, Cp_type, performance, v_wind[i], v_w, random_phases) k2 = Betti(x[i] + 0.5 * dt * k1, t[i] + 0.5 * dt, beta, T_E, Cp_type, performance, v_wind[i], v_w, random_phases)[0] k3 = Betti(x[i] + 0.5 * dt * k2, t[i] + 0.5 * dt, beta, T_E, Cp_type, performance, v_wind[i], v_w, random_phases)[0] @@ -1169,7 +1180,7 @@ def plot_fft(wave_eta, t): -#reproduce_save_driver(seeds) +reproduce_save_driver(seeds) diff --git a/reproduce_results.py b/reproduce_results.py index a7025f2..b7b6cc0 100644 --- a/reproduce_results.py +++ b/reproduce_results.py @@ -12,6 +12,7 @@ import matplotlib.pyplot as plt from scipy.signal import spectrogram from matplotlib.lines import Line2D +from gen_wind_Van_Der_Hoven import generate_wind def process_rotor_performance(input_file = "Cp_Ct.NREL5MW.txt"): """ @@ -100,6 +101,64 @@ def CpCtCq(TSR, beta, performance): # Get the C_p value at the index return C_p[TSR_index][pitch_index], C_t[TSR_index][pitch_index] +def gen_turbulence(v_bar, L, k_sigma_v, T_s, N_t, white_noise, + delta_omega = 0.002, M = 5000, N = 100): + """ + Generate turbulencec component for wind speed + + Parameters + ---------- + v_bar : int + Average wind speed + L : int + Turbulence length + k_sigma_v : float + Slope parameter + T_s : int + Time step + N_t : int + Number of step + white_noise : np.array + the white noise with mean = 0 and std = 1, has length + + Returns + ------- + Array of wind speed with turbulence + + """ + + # Step 1: Update the current values of the parameters in + # the turbulence component model + T_F = L / v_bar + sigma_v = k_sigma_v * v_bar + K_F = np.sqrt((2 * np.pi * T_F)/(4.20654 * T_s)) + + # Step 2: Calculate the discrete impulse response of the filter + delta_omega = 0.002 # Frequency step size + M = 5000 # Number of frequency points + N = 100 # Numerical parameters for convolution integration, divide + # finite integral from 0 to t to N regions + + # Discrete frequency domain P(omega) + P = np.zeros(M + 1) + for r in range(M + 1): + P[r] = np.real(K_F / (1 + 1j * r * delta_omega * T_F)**(5/6)) + + # Discrete impulse response h(k) === h(T_s*k), k range from 0 to N + h = np.zeros(N + 1) + for k in range(N + 1): + h[k] = T_s * delta_omega * (2/np.pi) * np.sum(P * np.cos(k * np.arange(M + 1) * T_s * delta_omega)) + + # Step 3: Generate the turbulence component in the interval using convolution + v_t = np.zeros(N_t + 1) + + # Zero-pad the white noise + white_noise_padded = np.pad(white_noise, (0, N), 'constant') + + for m in range(N_t + 1): + v_t[m] = T_s * np.sum(h * white_noise_padded[m : m + N + 1]) + + return v_bar + sigma_v * v_t def genWind(v_w, end_time, time_step, seed): """ @@ -263,8 +322,8 @@ def pierson_moskowitz_spectrum(U19_5, zeta, eta, t, random_phases): a_x = np.sum((omega**2) * a * exp_component * cos_component) a_y = -np.sum((omega**2) * a * exp_component * sin_component) - #return wave_eta, [v_x, v_y, a_x, a_y] - return 0, [0,0,0,0] + return wave_eta, [v_x, v_y, a_x, a_y] + #return 0, [0,0,0,0] def structure(x_1, beta, omega_R, t, Cp_type, performance, v_w, v_aveg, random_phases): @@ -772,18 +831,18 @@ def PI_blade_pitch_controller(omega_R, dt, beta, integral, error, i): # dicard data for first 500s discard_steps = int(500 / 0.5) - np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}_{seeds[2]}_wind_wave_time.npz', + np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}_wind_wave_time.npz', t=t, v_wind=v_wind, wave_eta=wave_eta) - t_sub = t[::steps][discard_steps:] - x_sub = x[::steps][discard_steps:] - v_wind_sub = v_wind[:len(t)][::steps][discard_steps:] - wave_eta_sub = np.array(wave_eta)[::steps][discard_steps:] - h_wave_sub = np.array(h_waves)[::steps][discard_steps:] - betas_sub = betas[::steps][discard_steps:] - Qt_list_sub = Qt_list[::steps][discard_steps:] + t_sub = t[::steps][:-discard_steps] + x_sub = x[::steps][:-discard_steps] + v_wind_sub = v_wind[:len(t)][::steps][:-discard_steps] + wave_eta_sub = np.array(wave_eta)[::steps][:-discard_steps] + h_wave_sub = np.array(h_waves)[::steps][:-discard_steps] + betas_sub = betas[::steps][:-discard_steps] + Qt_list_sub = Qt_list[::steps][:-discard_steps] return t_sub-t_sub[0], x_sub, v_wind_sub, wave_eta_sub, h_wave_sub, betas_sub, Qt_list_sub @@ -818,9 +877,25 @@ def main(end_time, v_w, x0, seeds, seed_wave, time_step = 0.05, Cp_type = 0): # modify this to change initial condition #[zeta, v_zeta, eta, v_eta, alpha, omega, omega_R] #v_wind = genWind(v_w, end_time, time_step, seeds) - v_wind = np.load(f'reproduced_results/turbsim_output/{seeds[0]}_{seeds[1]}.npy') + #v_wind = np.load(f'reproduced_results/turbsim_output/{seeds[0]}_{seeds[1]}.npy') #v_wind = genWind_seeds(seeds) #v_wind = np.full(45000, 20) + def insert_elements(arr, num_elements): + + result = [] + for i in range(len(arr) - 1): + start, end = arr[i], arr[i + 1] + result.extend(np.linspace(start, end, num_elements + 2)[:-1]) + result.append(arr[-1]) + return np.array(result) + + # generate a random seed + state_before = np.random.get_state() + np.random.seed(seeds) + white_noise = np.random.randn(end_time + 1) + np.random.set_state(state_before) + + v_wind = insert_elements(gen_turbulence(v_w, 180, 0.13, 1, end_time, white_noise), int(1/time_step)) # modify this to change run time and step size #[Betti, x0 (initial condition), start time, end time, time step, beta, T_E] @@ -882,8 +957,8 @@ def reproduce_save_driver(seeds): v_w = 20 end_time = 2000 #end_time < 3000 - seeds_wind = [seeds[0], seeds[1]] - seed_wave = seeds[2] + seeds_wind = seeds[0] + seed_wave = seeds[1] x0 = np.array([-2.61426271, @@ -896,7 +971,7 @@ def reproduce_save_driver(seeds): t, x, v_wind, wave_eta, h_wave, betas, Q_t = main(end_time, v_w, x0, seeds_wind, seed_wave) end_time -= 500 - np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}_{seeds[2]}.npz', + np.savez(f'reproduced_results/data/{seeds[0]}_{seeds[1]}.npz', t=t, x=x, v_wind=v_wind, @@ -916,25 +991,25 @@ def load_data(seeds): load the percentile and extreme value ''' - output_file_name = f'{seeds[0]}_{seeds[1]}_{seeds[2]}.npz' + output_file_name = f'{seeds[0]}_{seeds[1]}.npz' data = np.load(f'reproduced_results/data/{output_file_name}', allow_pickle=True) # Extracting the data t = data['t'][:-1] state = data['x'][:-1] - #beta = np.rad2deg(data['betas']) + beta = np.rad2deg(data['betas']) x = data['x'][:-1] wind_speed = data['v_wind'][:-1] wave_eta = data['wave_eta'][:-1] data.close() - + ''' pitch_rate = x[:, 5] pitch_acceleration = np.diff(pitch_rate) last_acceleration = pitch_acceleration[-1][None] pitch_acceleration = np.concatenate((pitch_acceleration, last_acceleration), axis=0)[:, None] state = np.concatenate((x[:, :6], pitch_acceleration), axis=1) - + ''' # Extracting percentile data @@ -985,8 +1060,8 @@ def plot_helper(ax): ax[1].set_xlim(0, t[-1]) # plot 7 states - for j in range(7): - #for j in range(6): + #for j in range(7): + for j in range(6): ax[j+2].plot(t, max_state[:,j], alpha=0.6, color='green', linewidth=0.5) ax[j+2].plot(t, min_state[:,j], alpha=0.6, color='orange', linewidth=0.5) @@ -1004,7 +1079,7 @@ def plot_helper(ax): ax[j+2].tick_params(axis='both', labelsize=16) - ''' + ax[8].plot(t, state[:, -1], color='black', linewidth=0.5) ax[8].set_xlabel('Time (s)', fontsize=12) #ax[j+2].set_ylabel(f'{state_names[j]}') @@ -1037,7 +1112,7 @@ def plot_helper(ax): Line2D([0], [0], color='orange', lw=1, alpha=0.6, label='The Minimum at Each Time Step')] ax[9].legend(handles=legend_elements, loc='center', fontsize=17.5) - + ''' # for 8 states including pitch acceleration: @@ -1112,7 +1187,7 @@ def plot_fft(wave_eta, t): #seeds = [966870, 8017384, 8986784] #surge events #seeds = [-4966241, 2503014, 9071735] -seeds = [ -402337, -6699134, 7762480] +#seeds = [ -402337, -6699134, 7762480] #surge velocity #2eeds = [-3919707, 4205247, 211083] #wave short corr @@ -1154,7 +1229,7 @@ def plot_fft(wave_eta, t): #seeds = [ 4200242, -3440907, 5773012] #std: 0.31226751594476654 #std = load_data(seeds) - +seeds =[7619098, 6844915] seed = [[-8615404, 1149694, 9191470], [3973823, 4556159, 3377501], @@ -1169,7 +1244,7 @@ def plot_fft(wave_eta, t): -#reproduce_save_driver(seeds) +reproduce_save_driver(seeds) diff --git a/reproduced_results/data/-402337_-6699134_7762480.npz b/reproduced_results/data/-402337_-6699134_7762480.npz index 4ecf139..2133f98 100644 Binary files a/reproduced_results/data/-402337_-6699134_7762480.npz and b/reproduced_results/data/-402337_-6699134_7762480.npz differ diff --git a/reproduced_results/data/-402337_-6699134_7762480_wind_wave_time.npz b/reproduced_results/data/-402337_-6699134_7762480_wind_wave_time.npz index fb748c3..32bfbf7 100644 Binary files a/reproduced_results/data/-402337_-6699134_7762480_wind_wave_time.npz and b/reproduced_results/data/-402337_-6699134_7762480_wind_wave_time.npz differ diff --git a/reproduced_results/data/2186117_5245701.npz b/reproduced_results/data/2186117_5245701.npz new file mode 100644 index 0000000..52ce438 Binary files /dev/null and b/reproduced_results/data/2186117_5245701.npz differ diff --git a/reproduced_results/data/2186117_5245701_wind_wave_time.npz b/reproduced_results/data/2186117_5245701_wind_wave_time.npz new file mode 100644 index 0000000..f4e46c1 Binary files /dev/null and b/reproduced_results/data/2186117_5245701_wind_wave_time.npz differ diff --git a/reproduced_results/data/609551_3823752.npz b/reproduced_results/data/609551_3823752.npz new file mode 100644 index 0000000..76bf3e9 Binary files /dev/null and b/reproduced_results/data/609551_3823752.npz differ diff --git a/reproduced_results/data/609551_3823752_wind_wave_time.npz b/reproduced_results/data/609551_3823752_wind_wave_time.npz new file mode 100644 index 0000000..5281348 Binary files /dev/null and b/reproduced_results/data/609551_3823752_wind_wave_time.npz differ diff --git a/reproduced_results/data/6699134_7762480.npz b/reproduced_results/data/6699134_7762480.npz new file mode 100644 index 0000000..75eae16 Binary files /dev/null and b/reproduced_results/data/6699134_7762480.npz differ diff --git a/reproduced_results/data/6699134_7762480_wind_wave_time.npz b/reproduced_results/data/6699134_7762480_wind_wave_time.npz new file mode 100644 index 0000000..5a6eec9 Binary files /dev/null and b/reproduced_results/data/6699134_7762480_wind_wave_time.npz differ diff --git a/reproduced_results/data/7619098_6844915.npz b/reproduced_results/data/7619098_6844915.npz new file mode 100644 index 0000000..41ad137 Binary files /dev/null and b/reproduced_results/data/7619098_6844915.npz differ diff --git a/reproduced_results/data/7619098_6844915_wind_wave_time.npz b/reproduced_results/data/7619098_6844915_wind_wave_time.npz new file mode 100644 index 0000000..9f83918 Binary files /dev/null and b/reproduced_results/data/7619098_6844915_wind_wave_time.npz differ diff --git a/reproduced_results/data/8236358_3518863.npz b/reproduced_results/data/8236358_3518863.npz new file mode 100644 index 0000000..9ca887f Binary files /dev/null and b/reproduced_results/data/8236358_3518863.npz differ diff --git a/reproduced_results/data/8236358_3518863_wind_wave_time.npz b/reproduced_results/data/8236358_3518863_wind_wave_time.npz new file mode 100644 index 0000000..6d975e9 Binary files /dev/null and b/reproduced_results/data/8236358_3518863_wind_wave_time.npz differ diff --git a/reproduced_results/data/8973511_9490074.npz b/reproduced_results/data/8973511_9490074.npz new file mode 100644 index 0000000..cf5f592 Binary files /dev/null and b/reproduced_results/data/8973511_9490074.npz differ diff --git a/reproduced_results/data/8973511_9490074_wind_wave_time.npz b/reproduced_results/data/8973511_9490074_wind_wave_time.npz new file mode 100644 index 0000000..44bf5a1 Binary files /dev/null and b/reproduced_results/data/8973511_9490074_wind_wave_time.npz differ diff --git a/reproduced_results/figure/large_surge_1_rotor.png b/reproduced_results/figure/large_surge_1_rotor.png index aa4bdb3..502f65f 100644 Binary files a/reproduced_results/figure/large_surge_1_rotor.png and b/reproduced_results/figure/large_surge_1_rotor.png differ