From d0d56c1fc966ecd78f7bc7a53310845ea7a36a43 Mon Sep 17 00:00:00 2001 From: Marcel Frommelt Date: Mon, 25 Jan 2021 20:25:41 +0900 Subject: [PATCH] changed basic thermal model, now realistic behavior --- DataRequest.py | 52 ++++++ README.md | 3 - input/Trial_2021-01-12.py | 174 ++++++++++++++++++ .../natural_constants.cpython-38.pyc | Bin 297 -> 317 bytes input/__pycache__/user_input.cpython-38.pyc | Bin 677 -> 744 bytes input/natural_constants.py | 3 +- input/user_input.py | 27 +-- main.py | 52 +++--- main_old.py | 146 +++++++++++++++ output/placeholder.txt => placeholder.txt | 0 10 files changed, 418 insertions(+), 39 deletions(-) create mode 100644 DataRequest.py delete mode 100644 README.md create mode 100644 input/Trial_2021-01-12.py create mode 100644 main_old.py rename output/placeholder.txt => placeholder.txt (100%) diff --git a/DataRequest.py b/DataRequest.py new file mode 100644 index 0000000..30372f5 --- /dev/null +++ b/DataRequest.py @@ -0,0 +1,52 @@ +import cdsapi + +c = cdsapi.Client() + +r = c.retrieve( + 'reanalysis-era5-pressure-levels', + { + 'product_type': 'reanalysis', + 'variable': [ + 'fraction_of_cloud_cover', 'geopotential', 'relative_humidity', + 'specific_humidity', 'temperature', 'u_component_of_wind', + 'v_component_of_wind', 'vertical_velocity', + ], + 'pressure_level': [ + '1', '2', '3', + '5', '7', '10', + '20', '30', '50', + '70', '100', '125', + '150', '175', '200', + '225', '250', '300', + '350', '400', '450', + '500', '550', '600', + '650', '700', '750', + '775', '800', '825', + '850', '875', '900', + '925', '950', '975', + '1000', + ], + 'year': '2016', + 'month': '07', + 'day': [ + '11', '12', '13', + '14', '15', '16', + '17', '18', + ], + 'time': [ + '00:00', '01:00', '02:00', + '03:00', '04:00', '05:00', + '06:00', '07:00', '08:00', + '09:00', '10:00', '11:00', + '12:00', '13:00', '14:00', + '15:00', '16:00', '17:00', + '18:00', '19:00', '20:00', + '21:00', '22:00', '23:00', + ], + 'area': [ + 72, -111, 67, + 22, + ], + 'format': 'netcdf', + }) +r.download('test.nc') \ No newline at end of file diff --git a/README.md b/README.md deleted file mode 100644 index 5d4db94..0000000 --- a/README.md +++ /dev/null @@ -1,3 +0,0 @@ -# BASTET - -ESBO Balloon Simulation and Trajectory Evaluation Tool \ No newline at end of file diff --git a/input/Trial_2021-01-12.py b/input/Trial_2021-01-12.py new file mode 100644 index 0000000..d5c3307 --- /dev/null +++ b/input/Trial_2021-01-12.py @@ -0,0 +1,174 @@ +import numpy as np +from astropy.time import Time + +# SIMPLE ATMOSPHERE MODEL + +def T_air(h): + if h >= 0 and h <= 11000: + res = 288.15 - 0.0065 * h + return res + elif h > 11000 and h <= 20000: + res = 216.65 + return res + elif h >= 20000: + res = 216.65 + 0.0010 * (h - 20000) + return res +def p_air(h): + if h >= 0 and h <= 11000: + res = 101325 * ((288.15 - 0.0065 * h)/288.15) ** 5.25577 + return res + elif h > 11000 and h <= 20000: + res = 22632 * np.exp(-(h - 11000)/6341.62) + return res + elif h > 20000: + res = 5474.87 * ((216.65 + 0.0010 * (h - 20000))/216.65) ** (-34.163) + return res +def rho_air(h): + res = p_air(h)/(R_air * T_air(h)) + return res + + +# USER INPUT + +start_height = 0 +start_utc = '2006-01-16 08:00:00.000' +start_lat = 0 # start latitude in [deg] +start_lon = 0 # start longitude in [deg] + + +# INITIALISATION + +t = 0 +h = start_height +lat = start_lat +lon = start_lon +utc = Time(start_utc) +v_z = 0 +p_gas = p_air(h) +T_gas = T_air(h) +T_film = T_gas +p_0 = 101325 + +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 + + + +dT_gas = (Q_g/(c_pg * m_g) - (g * M_a * m_g * T_g)/(c_pg * T_a * M_g) * dz/dt) * dt + + +V_b_new = m_gas * R_gas * T_gas_new/p_gas_new # r.h.s. with updated values! + + +v_w = 0 + +T_g1 = T_g +T_a1 = T_a +# Q_g1 = ??? +T_e1 = T_e + +# T_g2 = T_g1 + dt * (Q_g1/(c_pg * m_g1) - g * M_a * T_g1 * dh/(c_pg * T_a1 * M_g)) +# T_e2 = T_e1 + Q_e/(c_env * m_e) * dt + +T_gn = T_g + dt * (Q_g/(c_pg * m_g) - g * M_a * T_g * dh/(c_pg * T_a * M_g)) +T_en = T_e + Q_e/(c_env * m_e) * dt + +s_n = s + dt * v + 0.5 * a * dt ** 2 +v_n = v - v_w + dt * a + + + + + + + +T_g = T_gn +T_e = T_en +s = s_n +v = v_n + + +V = m_g * R_g * T_g / (p_g) + + + +dV_b = V_b_new - V_b # calculating volume increment +dT_gas = Q_ConvInt/(c_v * m_g) * dt + (gamma - 1) * T_g * (dm_g/m_g - dV_b/V_b) + + + + +V_b = V_b_new # updating current volume + + +dT_film = ((Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film)) * dt + + + + + + + + + +m_balloon = A_balloon * rho_envelope +m_gas = Vol_balloon * (M_gas * p_gas)/(R_c * T_gas) +m_gross = m_balloon + m_payload +m_tot = m_balloon + m_gas + m_payload + +GI = F_buoyant - g * (rho_a - rho_g) * Vol_balloon + +W = m_gross * g + +m_virtual = m_tot + C_virtual * rho_a * Vol_balloon + +R_z = GI + T_z - D_z - W + +dV_A = (R_z * dt)/m_virtual + +F_Drag = 0.5 * C_D * rho * np.abs(V) * V * A_cross + +Re = (rho * V * l)/my + +C_D = 24/Re + +my = my_0 * (T_0 + C)/(T + C) * (T/T_0) ** (3/2) + +C_D = 0.23175 - 0.15757 * f + 0.04744 * f ** 2 - 7.0412 * 10 ** (-3) * f ** 3 + 5.1534 * 10 ** (-4) * f ** 4 - 1.4835 * 10 ** (-5) * f ** 5 + +Vol_balloon = m_g * R_g * T_g/p_g + + +Vol_ballon_new = m_g * R_g * T_g/p_g # r.h.s. with updated values! +dVol = Vol_ballon_new - Vol_balloon # calculating volume increment +Vol_ballon = Vol_ballon_new # updating current volume + + +L_z = g * (M_atm * p_atm/(R * T_atm) - M_gas * p_gas/(R * T_gas)) * Vol_balloon + + + +dT_g = Q_ConvInt/(c_v * m_g) * dt + (gamma - 1) * T_g * (dm_g/m_g - dVol/Vol) + +dT_e = (Q_sun + Q_albedo + Q_IREarth + Q_IRsky + Q_ConvExt - Q_ConvInt + Q_IRout)/(c_e * m_e) + +Pr_g = my_g * c_pg/kappa_g +Pr_a = my_a * c_pa/kappa_a + +beta_g = 1/T_g +beta_a = 1/T_a + +Gr_g = d_b ** 3 * rho_g ** 2 * g * beta_g * (T_e - T_g)/(my_g ** 2) +Gr_a = d_b ** 3 * rho_a ** 2 * g * beta_a * (T_e - T_a)/(my_a ** 2) + + + +h = h + dt * V + 0.5 * a * dt ** 2 +V = V + dt * a + + + + + + + diff --git a/input/__pycache__/natural_constants.cpython-38.pyc b/input/__pycache__/natural_constants.cpython-38.pyc index deed592ccbd3e892d1c17c958e65690f54b2b63f..7483501f88b234370ad55a2bc4d8223baf6a4fcb 100644 GIT binary patch delta 107 zcmZ39?4ZNP9$`ezy`O?7Po4ubV c`z^-wiHp_QIDmSJI3`|*7T^K0SQrry00fW}#Q*>R diff --git a/input/__pycache__/user_input.cpython-38.pyc b/input/__pycache__/user_input.cpython-38.pyc index a6b29c8259cbbb1a60c703d20978668e3ef746e6..a9df815dd290c6da9f8030cc9ab8b125b2a7f3c4 100644 GIT binary patch delta 389 zcmZ3=`hvASl$V!_0SMTUI|pSZ=fJ`Dtf3LVl# z*112d61iue4yM;#x6c$%XJAOL28*}Yrvp6&2I*(sojmOH)!8BaSguLK_5@{zDsdwN z12bI%LtR5N1p^B!0|Ow^GcYj7JloC4;HN2hizPR{Am;vfS?AjbiSi&KC^3PTh_3S$&w3R4tQ3Ud^5 z3QH7A3TqT=3R@Ig3VReg5Obt(L~*8YMscNZMRBKaNAaZaMDeEZM)9TaMe(QbM+u|| zLjSP$obPbJk z4UH8HjjRj|fJo24zyP9RRlkFuro=6l-1vf=TWq=UX_+~>w^(!I(-Vttv4+LxCRW_y jOw1|BNQ@8j48Fx$1f(YKWc= 0: V_b = m_gas/rho_gas # calculate balloon volume from current gas mass and gas density - if V_b > V_max: - V_b = V_max - m_gas = V_max * rho_gas + if V_b > V_design: + V_b = V_design + m_gas = V_design * rho_gas else: pass @@ -55,7 +52,12 @@ while t <= t_end and h >= 0: m_tot = m_pl + m_film + m_gas m_virt = m_tot + c_virt * rho_air(h) * V_b - d_b = (6.0 * V_b / np.pi) ** (1/3) # calculate diameter of balloon from its volume + 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(h), d_b, v_z) # calculate drag force @@ -68,12 +70,6 @@ while t <= t_end and h >= 0: AZ, ELV = sun_angles_analytical(lat, lon, utc) A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) - A_surf = 4.94 * V_b ** (2/3) - - A_eff = A_surf - #A_surf1 = 4.94 * V_b - - # AirMass(x, p_0, ELV, h) # CALCULATIONS FOR THERMAL MODEL @@ -129,19 +125,31 @@ while t <= t_end and h >= 0: RoC = -v_z - dT_gas = (Q_ConvInt / (gamma * R_gas) - (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 + # 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 = 0 + v = v_z + s = h - dv_z = F/m_virt * dt # velocity increment + a = F/m_virt - v_z += dv_z * dt # velocity differential - h += v_z * dt # altitude differential + s_n = s + dt * v + 0.5 * a * dt ** 2 + dh = s_n - s + v_n = v - v_w + dt * a - p_gas = p_air(h) - T_gas = T_gas + dT_gas * dt - T_film = T_film + dT_film * dt + T_gn = T_gas + dt * (Q_ConvInt / (gamma * c_v * m_gas) - g * R_gas * T_gas * dh / (c_v * gamma * T_air(s) * 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 = v_n + + h = s + T_gas = T_g + T_film = T_e + v_z = v t += dt # time increment utc = Time(start_utc) + t * u.second diff --git a/main_old.py b/main_old.py new file mode 100644 index 0000000..79542e9 --- /dev/null +++ b/main_old.py @@ -0,0 +1,146 @@ +import numpy as np +import matplotlib.pyplot as plt +from datetime import datetime +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 models.simple_atmosphere import T_air, p_air, rho_air +from input.user_input import * +from input.natural_constants import * +from astropy.time import Time +import astropy.units as u +from models.thermal import AirMass + +# 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 +h = start_height +utc = Time(start_utc) +v_z = 0 +p_gas = p_air(h) +T_gas = T_air(h) +T_film = T_gas + +t_list = [] +h_list = [] + +while t <= t_end and h >= 0: + t_list.append(t) + h_list.append(h) + + 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(h) * 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(h), d_b, v_z) # calculate drag force + I = g * V_b * (rho_air(h) - rho_gas) # calculate gross inflation + W = g * m_gross # calculate weight (force) + F = I - W - D * np.sign(v_z) + + v_zrel = v_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(h), p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air(h), p_0, ELV, h))) + tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air(h) / p_0) + np.exp(-0.095 * p_air(h) / 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(h) ** 1.5) / (T_air(h) + 110.4) + my_gas = 1.895 * 10 ** -5 * (T_gas / 273.15) ** 0.647 + k_air = 0.0241 * (T_air(h) / 273.15) ** 0.9 + k_gas = 0.144 * (T_gas / 273.15) ** 0.7 + Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air(h) + Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas + + Gr_air = (rho_air(h) ** 2 * g * np.abs(T_film - T_air(h)) * d_b ** 3) / (T_air(h) * 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_zrel) * d_b * rho_air(h) / 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(h) - 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 + + ### + + dv_z = F/m_virt * dt # velocity increment + + v_z += dv_z * dt # velocity differential + h += v_z * dt # altitude differential + + p_gas = p_air(h) + T_gas = T_gas + dT_gas * dt + T_film = T_film + dT_film * dt + + t += dt # time increment + utc = Time(start_utc) + t * u.second + +plt.plot(t_list, h_list) +plt.show() diff --git a/output/placeholder.txt b/placeholder.txt similarity index 100% rename from output/placeholder.txt rename to placeholder.txt