import matplotlib.pyplot as plt import numpy as np H = 0.079 D = 0.090 c = 4180 rho = 1000 k_tc = 0.2 L = 3e-3 T_room = 20 V = H*np.pi*(D/2)**2 C = c*rho*V A = np.pi*D*H + 2*np.pi*(D/2)**2 G = (k_tc/L)*A dt = 0.1 t_start = 0 t_stop = 200 N_sim = int((t_stop - t_start)/dt) + 1 P_on = 1777 P_off = 0 t_array = np.zeros(N_sim) T_array = np.zeros(N_sim) T_room_array = np.zeros(N_sim) + T_room P_array = np.zeros(N_sim) T_min = 0 T_max = 100 T_k = T_init = 20.0 for k in range(0, N_sim): t_k = k*dt if (0 <= t_k <= t_start): P_k = 0 else: P_k = P_on dT_dt_k = (1/C)*(P_k + G*(T_room - T_k)) T_kp1 = T_k + dt*dT_dt_k T_kp1 = np.clip(T_kp1, T_min, T_max) t_array[k] = t_k T_array[k] = T_k P_array[k] = P_k T_k = T_kp1 T_100 = 100 # [oC] indexes_less_than_100 = np.argwhere(T_array < T_100) t_100 = t_array[np.amax(indexes_less_than_100)] print('t_100 [s] =', f'{t_100:.2f}') plt.close('all') plt.figure(num=1, figsize=(12, 9)) plt.subplot(2, 1, 1) plt.plot(t_array, T_array,'b', label='T') plt.plot(t_array, T_room_array, 'g', label='T_room') plt.legend() plt.ylabel('[deg C]') plt.grid() plt.subplot(2, 1, 2) plt.plot(t_array, P_array, 'r', label='P') plt.legend() plt.grid() plt.xlabel('t [s]') plt.ylabel('[W]') # plt.savefig('prog_sim_kettle_detect_p.pdf') plt.show()