Loading dthydro_cfd/generator_models.py +85 −20 Original line number Diff line number Diff line Loading @@ -42,16 +42,18 @@ def parks_transform(in_vec, gamma, beta_0=BETA_1, beta_d=BETA_2, beta_q=BETA_3, return out_vec def armature_voltage(R, X_q_pp, X_d_pp, E_d_pp, E_q_dd, I_d, I_q): def armature_voltage(R, X_q_pp, X_d_pp, E_d_pp, E_q_pp, I_d, I_q): """ Armature voltage equations. eq. 11.100. Args: X_d_pp: transient reluctance d X_q_pp: transient reluctance q Returns: [V_d, V_q] """ i_arr = np.array([[I_d, I_q],]).T E_arr = np.array([[E_d_pp, E_q_dd],]).T i_arr = np.array([I_d, I_q]) E_arr = np.array([E_d_pp, E_q_pp]) tr_m = np.array([[R, X_q_pp], [-X_d_pp, R]]) v_arr = E_arr - np.matmul(tr_m, i_arr) return v_arr Loading Loading @@ -199,12 +201,23 @@ class sync_gen5(object): def E_d_pp(self): return self.x[4] def arm_v(self, R=10.0): """ armature_voltage """ X_d = self.machine_params["Xd"] X_d_p = self.machine_params["Xd1"] X_d_pp = self.machine_params["Xd2"] X_q = self.machine_params["Xq"] X_q_p = self.machine_params["Xq1"] X_q_pp = self.machine_params["Xq2"] return armature_voltage( R, X_q_pp, X_d_pp, self.E_d_pp, self.E_q_pp, self.I_d, self.I_q) def test_gen5(): """ Test of the 5th order sychronous generator model """ from scipy.integrate import solve_ivp # Driving voltage e_f = 140.0 # Driving voltage, generator main field voltage e_f = 1000.0 delta_w = 2.0 * np.pi * 1e-4 # rad / s w_0 = W_S + delta_w # Line Currents in system frame (amps) Loading @@ -212,47 +225,99 @@ def test_gen5(): #i_b = 1.0 #i_c = 1.0 i_0_0 = 0. # Must be 0 bc, (i_a + i_b + i_c) / sqrt(3) = 0 = i_0 i_d_0 = 6.0e3 i_q_0 = 6.0e3 i_a, i_b, i_c = parks_transform( np.asarray([i_0_0, i_d_0, i_q_0]), gamma=0., inv_transform=True) i_sum = np.sum((i_a, i_b, i_c)) I_0_0 = 0. # Must be 0 bc, (i_a + i_b + i_c) / sqrt(3) = 0 = i_0 I_d_0 = 6.0e3 * 0. I_q_0 = 6.0e3 * 1. I_a, I_b, I_c = parks_transform( np.asarray([I_0_0, I_d_0, I_q_0]), gamma=0., inv_transform=True) i_abc = np.array([I_a, I_b, I_c]) * np.sqrt(3.) i_sum = np.sum(i_abc) assert np.allclose((i_sum), (0.)) # apply transform to rotor frame i_0, i_d, i_q = parks_transform(np.asarray([i_a, i_b, i_c]), gamma=0.) assert np.allclose((i_0), (0.0)) assert np.allclose((i_0, i_d, i_q), (i_0_0, i_d_0, i_q_0)) gen5 = sync_gen5(w_0, W_S, I_d=i_d, I_q=i_q) I_0, I_d, I_q = parks_transform(np.asarray([I_a, I_b, I_c]), gamma=0.) i_0dq = np.array([I_0, I_d, I_q]) * np.sqrt(3.) assert np.allclose((I_0), (0.0)) assert np.allclose((I_0, I_d, I_q), (I_0_0, I_d_0, I_q_0)) gen5 = sync_gen5(w_0, W_S, I_d=I_d, I_q=I_q) x_0 = deepcopy(gen5.x) t = 0. # Mechanical power input P_m = 25.0e6 # W # Current inputs I_d = i_d I_q = i_q # Driving voltage I_0 = I_0 I_d = I_d I_q = I_q # Driving voltage from excitation system E_f = e_f # For armature_voltage calc R = 10.0 # test call to system derivative test_dx_dt = gen5(t, x_0, P_m, I_d, I_q, E_f) # test time integration of the system dt = 1e-2 dt = 1e-1 x_new = x_0 p_e_list = [] t_list = [] w_list = [] for i in range(8000): sol = solve_ivp(gen5, [t, t+dt], x_new, args=(P_m, I_d, I_q, E_f)) x_new = sol.y[:, -1] # reset rpm # gen5.x[0] = 0. t = t + dt print("Time: %0.3e, omega: %0.5e, P_e: %0.3e" % (t, gen5.w, gen5.power_e)) p_e_list.append(gen5.power_e) w_list.append(gen5.w) t_list.append(t) aV = gen5.arm_v(R=100.0) V_0, V_d, V_q = 0., aV[0], aV[1] # Eq. 11.80 v_0dq = np.array([V_0, V_d, V_q]) * np.sqrt(3.) # Eq. 11.82 i_0dq = np.array([I_0, I_d, I_q]) * np.sqrt(3.) v_abc = parks_transform(v_0dq, gamma=0, inv_transform=True) print("Voltages: ", v_abc) print("Currents: ", i_abc) p_g = np.sum(i_abc * v_abc) p_g_b = np.sum(i_0dq * v_0dq) print("Time: %0.3e, omega: %0.5e, P_e: %0.3e, P_g: %0.3e, %0.3e" % \ (t, gen5.w, gen5.power_e, p_g, p_g_b)) t_list = np.asarray(t_list) p_e_list = np.asarray(p_e_list) w_list = np.asarray(w_list) print(sol.t) print(sol.y) sol.y[:, -1][-1] sol.y[:, -1][-2] sol.y[:, -1][-3] # plots import matplotlib.pyplot as plt plt.figure() plt.plot(t_list, p_e_list, label="Gen air gap power") plt.xlabel("Time [s]") plt.ylabel("Power [W]") plt.grid(ls="--", alpha=0.5) plt.legend() plt.savefig("sync_gen5_p_e.png") plt.close() plt.figure() plt.plot(t_list, w_list / (2*np.pi) * 60., label="Rotation rate") plt.xlabel("Time [s]") plt.ylabel("RPM") plt.grid(ls="--", alpha=0.5) plt.legend() plt.savefig("sync_gen5_w.png") plt.close() if __name__ == "__main__": test_gen5() Loading
dthydro_cfd/generator_models.py +85 −20 Original line number Diff line number Diff line Loading @@ -42,16 +42,18 @@ def parks_transform(in_vec, gamma, beta_0=BETA_1, beta_d=BETA_2, beta_q=BETA_3, return out_vec def armature_voltage(R, X_q_pp, X_d_pp, E_d_pp, E_q_dd, I_d, I_q): def armature_voltage(R, X_q_pp, X_d_pp, E_d_pp, E_q_pp, I_d, I_q): """ Armature voltage equations. eq. 11.100. Args: X_d_pp: transient reluctance d X_q_pp: transient reluctance q Returns: [V_d, V_q] """ i_arr = np.array([[I_d, I_q],]).T E_arr = np.array([[E_d_pp, E_q_dd],]).T i_arr = np.array([I_d, I_q]) E_arr = np.array([E_d_pp, E_q_pp]) tr_m = np.array([[R, X_q_pp], [-X_d_pp, R]]) v_arr = E_arr - np.matmul(tr_m, i_arr) return v_arr Loading Loading @@ -199,12 +201,23 @@ class sync_gen5(object): def E_d_pp(self): return self.x[4] def arm_v(self, R=10.0): """ armature_voltage """ X_d = self.machine_params["Xd"] X_d_p = self.machine_params["Xd1"] X_d_pp = self.machine_params["Xd2"] X_q = self.machine_params["Xq"] X_q_p = self.machine_params["Xq1"] X_q_pp = self.machine_params["Xq2"] return armature_voltage( R, X_q_pp, X_d_pp, self.E_d_pp, self.E_q_pp, self.I_d, self.I_q) def test_gen5(): """ Test of the 5th order sychronous generator model """ from scipy.integrate import solve_ivp # Driving voltage e_f = 140.0 # Driving voltage, generator main field voltage e_f = 1000.0 delta_w = 2.0 * np.pi * 1e-4 # rad / s w_0 = W_S + delta_w # Line Currents in system frame (amps) Loading @@ -212,47 +225,99 @@ def test_gen5(): #i_b = 1.0 #i_c = 1.0 i_0_0 = 0. # Must be 0 bc, (i_a + i_b + i_c) / sqrt(3) = 0 = i_0 i_d_0 = 6.0e3 i_q_0 = 6.0e3 i_a, i_b, i_c = parks_transform( np.asarray([i_0_0, i_d_0, i_q_0]), gamma=0., inv_transform=True) i_sum = np.sum((i_a, i_b, i_c)) I_0_0 = 0. # Must be 0 bc, (i_a + i_b + i_c) / sqrt(3) = 0 = i_0 I_d_0 = 6.0e3 * 0. I_q_0 = 6.0e3 * 1. I_a, I_b, I_c = parks_transform( np.asarray([I_0_0, I_d_0, I_q_0]), gamma=0., inv_transform=True) i_abc = np.array([I_a, I_b, I_c]) * np.sqrt(3.) i_sum = np.sum(i_abc) assert np.allclose((i_sum), (0.)) # apply transform to rotor frame i_0, i_d, i_q = parks_transform(np.asarray([i_a, i_b, i_c]), gamma=0.) assert np.allclose((i_0), (0.0)) assert np.allclose((i_0, i_d, i_q), (i_0_0, i_d_0, i_q_0)) gen5 = sync_gen5(w_0, W_S, I_d=i_d, I_q=i_q) I_0, I_d, I_q = parks_transform(np.asarray([I_a, I_b, I_c]), gamma=0.) i_0dq = np.array([I_0, I_d, I_q]) * np.sqrt(3.) assert np.allclose((I_0), (0.0)) assert np.allclose((I_0, I_d, I_q), (I_0_0, I_d_0, I_q_0)) gen5 = sync_gen5(w_0, W_S, I_d=I_d, I_q=I_q) x_0 = deepcopy(gen5.x) t = 0. # Mechanical power input P_m = 25.0e6 # W # Current inputs I_d = i_d I_q = i_q # Driving voltage I_0 = I_0 I_d = I_d I_q = I_q # Driving voltage from excitation system E_f = e_f # For armature_voltage calc R = 10.0 # test call to system derivative test_dx_dt = gen5(t, x_0, P_m, I_d, I_q, E_f) # test time integration of the system dt = 1e-2 dt = 1e-1 x_new = x_0 p_e_list = [] t_list = [] w_list = [] for i in range(8000): sol = solve_ivp(gen5, [t, t+dt], x_new, args=(P_m, I_d, I_q, E_f)) x_new = sol.y[:, -1] # reset rpm # gen5.x[0] = 0. t = t + dt print("Time: %0.3e, omega: %0.5e, P_e: %0.3e" % (t, gen5.w, gen5.power_e)) p_e_list.append(gen5.power_e) w_list.append(gen5.w) t_list.append(t) aV = gen5.arm_v(R=100.0) V_0, V_d, V_q = 0., aV[0], aV[1] # Eq. 11.80 v_0dq = np.array([V_0, V_d, V_q]) * np.sqrt(3.) # Eq. 11.82 i_0dq = np.array([I_0, I_d, I_q]) * np.sqrt(3.) v_abc = parks_transform(v_0dq, gamma=0, inv_transform=True) print("Voltages: ", v_abc) print("Currents: ", i_abc) p_g = np.sum(i_abc * v_abc) p_g_b = np.sum(i_0dq * v_0dq) print("Time: %0.3e, omega: %0.5e, P_e: %0.3e, P_g: %0.3e, %0.3e" % \ (t, gen5.w, gen5.power_e, p_g, p_g_b)) t_list = np.asarray(t_list) p_e_list = np.asarray(p_e_list) w_list = np.asarray(w_list) print(sol.t) print(sol.y) sol.y[:, -1][-1] sol.y[:, -1][-2] sol.y[:, -1][-3] # plots import matplotlib.pyplot as plt plt.figure() plt.plot(t_list, p_e_list, label="Gen air gap power") plt.xlabel("Time [s]") plt.ylabel("Power [W]") plt.grid(ls="--", alpha=0.5) plt.legend() plt.savefig("sync_gen5_p_e.png") plt.close() plt.figure() plt.plot(t_list, w_list / (2*np.pi) * 60., label="Rotation rate") plt.xlabel("Time [s]") plt.ylabel("RPM") plt.grid(ls="--", alpha=0.5) plt.legend() plt.savefig("sync_gen5_w.png") plt.close() if __name__ == "__main__": test_gen5()
