import numpy as np import matplotlib.pyplot as plt from datetime import datetime start = datetime.now() print(start) begin_time = datetime.now() from models.sun import sun_angles_analytical from models.sun import sun_angles_astropy from models.drag import drag, c_d from input.user_input import * from input.natural_constants import * from astropy.time import Time import astropy.units as unit from models.thermal import AirMass from netCDF4 import Dataset data = Dataset("test2021.nc", "r", format="NETCDF4") ERAtime = data.variables['time'][:] # time ERAlat = data.variables['latitude'][:] # latitude ERAlon = data.variables['longitude'][:] # longitude ERAz = data.variables['z'][:]/g # geopotential to geopotential height ERApress = data.variables['level'][:] # pressure level ERAtemp = data.variables['t'][:] vw_x = data.variables['u'][:] # v_x vw_y = data.variables['v'][:] # v_y vw_z = data.variables['w'][:] # v_z # t: simulation time (0 .. t_max) in [s] # utc = Time(start_utc) + t * u.second # current time (UTC) lat = start_lat # deg lon = start_lon # deg # date = Time('2020-12-13 12:00:00.000') # INITIALISATION t = 0 # simulation time in seconds h = start_height utc = Time(start_utc) epoch_diff = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 t_epoch = epoch_diff sq_diff_lat = (ERAlat - lat) ** 2 sq_diff_lon = (ERAlon - lon) ** 2 sq_diff_time = (ERAtime - t_epoch) ** 2 min_index_lat = sq_diff_lat.argmin() min_index_lon = sq_diff_lat.argmin() min_index_time = sq_diff_time.argmin() sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 min_index_h = sq_diff_h.argmin() h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] p_grid = ERApress[min_index_h] * 100 T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid v = v_grid u = u_grid w = w_grid T_air = T_grid p_air = p_grid rho_air = p_air/(R_air * T_air) v_x = 0 v_y = 0 v_z = 0 T_gas = T_air T_film = T_gas t_list = [] h_list = [] lat_list = [] lon_list = [] p_list = [] Temp_list = [] rho_list = [] while t <= t_end and h >= 0: t_list.append(t) h_list.append(h) lat_list.append(lat) lon_list.append(lon) p_list.append(p_air) Temp_list.append(T_air) rho_list.append(rho_air) t_epoch = epoch_diff + t/3600 sq_diff_lat = (ERAlat - lat) ** 2 sq_diff_lon = (ERAlon - lon) ** 2 sq_diff_time = (ERAtime - t_epoch) ** 2 min_index_lat = sq_diff_lat.argmin() min_index_lon = sq_diff_lat.argmin() min_index_time = sq_diff_time.argmin() sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 min_index_h = sq_diff_h.argmin() h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] p_grid = ERApress[min_index_h] * 100 T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid p_air = p_grid p_gas = p_air T_air = T_grid rho_air = p_air / (R_air * T_air) u = u_grid v = v_grid w = w_grid v_relx = u - v_x v_rely = v - v_y v_relz = w - v_z v_rel = (v_relx ** 2 + v_rely ** 2 + v_relz ** 2) ** (1/2) p_gas = p_air rho_gas = p_gas/(R_gas * T_gas) # calculate gas density through ideal(!) gas equation V_b = m_gas/rho_gas # calculate balloon volume from current gas mass and gas density if V_b > V_design: V_b = V_design m_gas = V_design * rho_gas else: pass m_gross = m_pl + m_film m_tot = m_pl + m_film + m_gas m_virt = m_tot + c_virt * rho_air * V_b d_b = 1.383 * V_b ** (1/3) # calculate diameter of balloon from its volume L_goreB = 1.914 * V_b ** (1/3) h_b = 0.748 * d_b A_surf = 4.94 * V_b ** (2/3) A_surf1 = 4.94 * V_design ** (2/3) * (1 - np.cos(np.pi * L_goreB/L_goreDesign)) A_eff = 0.65 * A_surf + 0.35 * A_surf1 A_top = np.pi/4 * d_b ** 2 D = drag(c_d, rho_air, d_b, v_rel) # calculate drag force if v_rel == 0: Drag_x = 0 Drag_y = 0 Drag_z = 0 else: Drag_x = D * v_relx/v_rel Drag_y = D * v_rely/v_rel Drag_z = D * v_relz/v_rel I = g * V_b * (rho_air - rho_gas) # calculate gross inflation W = g * m_gross # calculate weight (force) F = I - W + Drag_z AZ, ELV = sun_angles_analytical(lat, lon, utc) A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) # CALCULATIONS FOR THERMAL MODEL if ELV >= -(180 / np.pi * np.arccos(R_E / (R_E + h))): tau_atm = 0.5 * (np.exp(-0.65 * AirMass(p_air, p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air, p_0, ELV, h))) tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air / p_0) + np.exp(-0.095 * p_air / p_0)) else: tau_atm = 0 tau_atmIR = 0 doy = int(utc.doy) MA = (357.52911 + 0.98560028 * (utc.jd - 2451545)) % 360 # in degree, reference: see folder "literature" TA = MA + 2 * e * np.sin(np.deg2rad(MA)) + 5 / 4 * e ** 2 * np.sin(np.deg2rad(2 * MA)) I_Sun = 1367.5 * ((1 + e * np.cos(np.deg2rad(TA))) / (1 - e ** 2)) ** 2 I_SunZ = I_Sun * tau_atm q_sun = I_SunZ q_IRground = epsilon_ground * sigma * T_ground ** 4 q_IREarth = q_IRground * tau_atmIR if ELV <= 0: q_Albedo = 0 else: q_Albedo = Albedo * I_Sun * np.sin(np.deg2rad(ELV)) my_air = (1.458 * 10 ** -6 * T_air ** 1.5) / (T_air + 110.4) my_gas = 1.895 * 10 ** -5 * (T_gas / 273.15) ** 0.647 k_air = 0.0241 * (T_air / 273.15) ** 0.9 k_gas = 0.144 * (T_gas / 273.15) ** 0.7 Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas Gr_air = (rho_air ** 2 * g * np.abs(T_film - T_air) * d_b ** 3) / (T_air * my_air ** 2) Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25 HC_free = Nu_air * k_air / d_b Re = np.abs(v_relz) * d_b * rho_air / my_air HC_forced = k_air / d_b * (2 + 0.41 * Re ** 0.55) HC_internal = 0.13 * k_gas * ((rho_gas ** 2 * g * np.abs(T_film - T_gas) * Pr_gas) / (T_gas * my_air ** 2)) ** ( 1 / 3) HC_external = np.maximum(HC_free, HC_forced) HalfConeAngle = np.arcsin(R_E / (R_E + h)) ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 Q_Sun = alpha_VIS * A_proj * q_sun * (1 + tau_VIS / (1 - r_VIS)) Q_Albedo = alpha_VIS * A_surf * q_Albedo * ViewFactor * (1 + tau_VIS / (1 - r_VIS)) Q_IREarth = alpha_IR * A_surf * q_IREarth * ViewFactor * (1 + tau_IR / (1 - r_IR)) Q_IRfilm = sigma * epsilon * alpha_IR * A_surf * T_film ** 4 * 1 / (1 - r_IR) Q_IRout = sigma * epsilon * A_surf * T_film ** 4 * (1 * tau_IR / (1 - r_IR)) Q_ConvExt = HC_external * A_eff * (T_air - T_film) Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) # RoC = -v_z # dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * g / (rho_gas * R_gas) * RoC) * dt # dT_film = ((Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film)) * dt # v_w = w # v = v_z s = h a_x = Drag_x/m_virt a_y = Drag_y/m_virt a_z = F/m_virt # DIFFERENTIAL EQUATIONS dx = dt * v_x + 0.5 * a_x * dt ** 2 # lon dy = dt * v_y + 0.5 * a_y * dt ** 2 # lat lon = lon + np.rad2deg(np.arctan(dx/((6371229.0 + h) * np.cos(np.deg2rad(lat))))) lat = lat + np.rad2deg(np.arctan(dy/(6371229.0 + h))) v_xn = v_x + dt * a_x v_yn = v_y + dt * a_y s_n = s + dt * v_z + 0.5 * a_z * dt ** 2 dh = s_n - s v_n = v_z + dt * a_z T_gn = T_gas + dt * (Q_ConvInt / (gamma * c_v * m_gas) - g * R_gas * T_gas * dh / (c_v * gamma * T_air * R_air)) T_en = T_film + (Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film) * dt T_g = T_gn T_e = T_en s = s_n v_x = v_xn v_y = v_yn v_z = v_n h = s T_gas = T_g T_film = T_e # v_z = v t += dt # time increment utc = Time(start_utc) + t * unit.second end = datetime.now() print(end) plt.plot(t_list, h_list, 'k-') plt.plot(t_list, p_list, 'r-') plt.plot(t_list, Temp_list, 'b-') plt.plot(t_list, rho_list, 'g-') plt.show()