diff --git a/main.py b/main.py index 614ce41..e56564e 100644 --- a/main.py +++ b/main.py @@ -1,56 +1,274 @@ import numpy as np +import math +import xarray as xr import matplotlib.pyplot as plt +from scipy.spatial import cKDTree +from scipy import interpolate +import cartopy.crs as ccrs 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 models.simple_atmosphere import T_air, p_air, rho_air +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 u +import astropy.units as unit from models.thermal import AirMass +from netCDF4 import Dataset +import pandas as pd + +data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Mappe1.xls', sheet_name='Tabelle3') + +comp_time = pd.DataFrame(data, columns= ['Time']) +comp_height = pd.DataFrame(data, columns= ['Height']) +comp_lat = pd.DataFrame(data, columns= ['Latitude']) +comp_lon = pd.DataFrame(data, columns= ['Longitude']) + +Reynolds = [0.0, 0.002, 0.004, 0.007, 0.01] #, 0.04, 0.07, 0.1, 0.2, 0.3, 0.5, 0.7, 1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 200, +# 300, 500, 700, 1000, 2000, 3000, 5000, 7000, 10000, 20000, 30000, 50000] +drag_coeff = [20000, 12000, 6000, 3429, 2400] #, 600, 343, 240, 120, 80, 49, 36.5, 26.5, 14.4, 10.4, 6.9, 5.4, 4.1, 2.55, 2, 1.5, +# 1.27, 1.07, 0.77, 0.65, 0.55, 0.5, 0.46, 0.42, 0.4, 0.385, 0.39, 0.405, 0.45, 0.47, 0.49] + +cd_interp1d = interpolate.interp1d(Reynolds, drag_coeff) + +Re_list = [] + +def cd_func(Re): + w = np.log10(Re) + + if Re <= 0.01: + c_d = 9/2 + 24/Re + elif Re <= 20: + c_d = 24/Re * (1 + 0.1315 * Re ** (0.82 - 0.05 * w)) + elif Re <= 260: + c_d = 24/Re * (1 + 0.1935 * Re ** 0.6305) + elif Re <= 1.5 * 10e3: + c_d = 10 ** (1.6435 - 1.1242 * w + 0.1558 * w ** 2) + elif Re <= 1.2 * 10e4: + c_d = 10 ** (-2.4571 + 2.5558 * w - 0.9295 * w ** 2 + 0.1049 * w ** 3) + elif Re <= 4.4 * 10e4: + c_d = 10 ** (-1.9181 + 0.6370 * w - 0.0636 * w ** 3) + elif Re <= 3.38 * 10e5: + c_d = 10 ** (-4.3390 + 1.5809 * w - 0.1546 * w ** 2) + elif Re <= 4 * 10e5: + c_d = 29.78 - 5.3 * w + elif Re <= 10 ** 6: + c_d = 0.1 * w - 0.49 + elif Re > 10e6: + c_d = 0.19 - 8 * 10e4/Re + else: + c_d = 0.47 + + return c_d + +def cd_func_garg(Re): + if Re == 0: + c_d = 2.0 + elif Re < 2 * 10e5: + c_d = 5.4856 * 10e9 * np.tanh(4.3774 * 10e-9 / Re) + 0.0709 * np.tanh(700.6575 / Re) + 0.3894 * np.tanh(74.1539 / Re) - 0.1198 * np.tanh(7429.0843 / Re) + 1.7174 * np.tanh(9.9851 / Re + 2.3384) + 0.4744 + elif 2 * 10e5 < Re <= 10e6: + c_d = 8 * 10e-6 * ((Re/6530) ** 2 + np.tanh(Re) - 8 * np.log(Re) / np.log(10)) - 0.4119 * np.exp(-2.08 * 10 ** 43 / (Re + Re ** 2) ** 4) - 2.1344 * np.exp(-((np.log(Re ** 2 + 10.7563) / np.log(10)) ** 2 + 9.9867) / Re) + 0.1357 * np.exp(-((Re / 1620) ** 2 + 10370) / Re) - 8.5 * 10e-3 * (2 * np.log(np.tanh(np.tanh(Re))) / np.log(10) - 2825.7162) / Re + 2.4795 + elif Re > 10e6: + c_d = 0.212546 + else: + c_d = 2.0 + return 0.80 # c_d + +#cd_interp1d = interpolate.interp1d(Reynolds, drag_coeff) + +def transform(lon, lat, t): + # WGS 84 reference coordinate system parameters + A = 6378137.0 # major axis [m] + E2 = 6.69437999014e-3 # eccentricity squared + + t_s = 1000 * t + + lon_rad = np.radians(lon) + lat_rad = np.radians(lat) + # convert to cartesian coordinates + r_n = A / (np.sqrt(1 - E2 * (np.sin(lat_rad) ** 2))) + x = r_n * np.cos(lat_rad) * np.cos(lon_rad) + y = r_n * np.cos(lat_rad) * np.sin(lon_rad) + z = r_n * (1 - E2) * np.sin(lat_rad) + return x, y, z, t_s + +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'][:] # temperature in K + +vw_x = data.variables['u'][:] # v_x +vw_y = data.variables['v'][:] # v_y +vw_z = data.variables['w'][:] # v_z + +lon_era2d, lat_era2d, time_era = np.meshgrid(ERAlon, ERAlat, ERAtime) + +xs, ys, zs, ts = transform(lon_era2d.flatten(), lat_era2d.flatten(), time_era.flatten()) + +tree = cKDTree(np.column_stack((xs, ys, zs, ts))) # ! + + -# 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 +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 + + + +xt, yt, zt, tt = transform(lon, lat, t) # test coordinates + +d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8) # longitude, latitude, time in h # ! +w = 1.0 / d[0] ** 2 + +lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0] +lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1] +time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2] + +interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w) +interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w) +interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w) +interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w) +interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w) + +pressure_hPa = np.array([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]) + +pressure = 100 * pressure_hPa +# print(pressure) + + +# height_interp1d = interpolate.interp1d(interp4d_height, pressure) +temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp) +press_interp1d = interpolate.interp1d(interp4d_height, pressure) +vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x) +vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y) +vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) + + +#u = 0 +#v = 0 +#w = 0 +T_air = temp_interp1d(h) +p_air = press_interp1d(h) +rho_air = p_air/(R_air * T_air) + +u = vw_x_interp1d(h) +v = vw_y_interp1d(h) +w = -1 / g * vw_z_interp1d(h) * R_air * T_air / p_air + + +v_x = 0 +v_y = 0 v_z = 0 -p_gas = p_air(h) -T_gas = T_air(h) + +v_rel = 0 + +T_gas = T_air T_film = T_gas t_list = [] h_list = [] +v_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) + v_list.append(v_rel) + 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 + + xt, yt, zt, tt = transform(lon, lat, t_epoch) # current balloon coordinates in cartesian coordinates: x [m], y [m], z [m], time [h since 1900-01-01 00:00:00.0] + + d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8) # longitude, latitude, time in h # ! + w = 1.0 / d[0] ** 2 + + lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0] + lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1] + time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2] + + interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w) + interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w) + interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w) + interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w) + interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w) + + press_interp1d = interpolate.interp1d(interp4d_height, pressure) + temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp) + vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x) + vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y) + vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) + + + try: + p_air = press_interp1d(h) + T_air = temp_interp1d(h) + except: + if (h > interp4d_height[0]): + h = interp4d_height[0] + else: + h = interp4d_height[-1] + + p_air = press_interp1d(h) + T_air = temp_interp1d(h) + + #print("height") + #print(h) + #print("temperature") + #print(temp_interp1d(h)) + #print("pressure") + #print(press_interp1d(h)) + #print("vw_x") + #print(vw_x_interp1d(h)) + + p_gas = p_air + rho_air = p_air / (R_air * T_air) + u = vw_x_interp1d(h) + v = vw_y_interp1d(h) + w = -1/g * vw_z_interp1d(h) * R_air * T_air / p_air + + 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) 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 + 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 + 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) @@ -60,13 +278,6 @@ while t <= t_end and h >= 0: 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))) @@ -74,8 +285,8 @@ while t <= t_end and h >= 0: # 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)) + 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 @@ -95,17 +306,19 @@ while t <= t_end and h >= 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_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(h) / 273.15) ** 0.9 + 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(h) + Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air 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) + 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_zrel) * d_b * rho_air(h) / my_air + Re = np.abs(v_relz) * d_b * rho_air / my_air + + Re_list.append(Re) 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)) ** ( @@ -120,39 +333,121 @@ while t <= t_end and h >= 0: 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_ConvExt = HC_external * A_eff * (T_air - T_film) Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) - RoC = -v_z + try: + c_d = cd_func(Re) + except: + print(Re) + 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 + + # 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 - v_w = 0 - v = v_z s = h - a = F/m_virt + a_x = Drag_x/m_virt + a_y = Drag_y/m_virt + a_z = F/m_virt - s_n = s + dt * v + 0.5 * a * dt ** 2 + # 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 - v_w + dt * a + 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(s) * R_air)) + 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 = v_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 +# v_z = v t += dt # time increment - utc = Time(start_utc) + t * u.second + utc = Time(start_utc) + t * unit.second -plt.plot(t_list, h_list) + + +#print(len(lon_list), len(lat_list)) +#""" +ax = plt.axes(projection=ccrs.Mollweide()) +#ax.coastlines() +ax.stock_img() +plt.plot(start_lon, start_lat, 'rx', transform=ccrs.Geodetic()) +plt.plot(lon_list, lat_list, 'k--', transform=ccrs.Geodetic()) +plt.plot(comp_lon, comp_lat, 'r-', transform=ccrs.Geodetic()) +#plt.xticks() plt.show() +#""" + + +### WORKS! +### +###dset = xr.open_dataset("test2021.nc") +### +###print(dset['t'][1][0]) # [time] [level] +### +#fig = plt.figure() #figsize=[120,50]) +### +#ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) +### +###dset['t'][1][0].plot(ax=ax, cmap='jet', transform=ccrs.PlateCarree()) +###ax.coastlines() +###plt.show() + +#print(len(ERAlat)) +#print(len(ERAlon)) +#print(len(ERAz)) + + + + + +#""" +plt.plot(comp_time, comp_height, 'r--', label='PoGo Flight Test') +#plt.plot(t_list, Re_list, 'k-', label="Reynolds-Number") +#plt.plot(t_list, v_list, 'k--', label='Ballon v_rel') +plt.plot(t_list, h_list, 'k-', label='Balloon Altitude') +#plt.plot(t_list, p_list, 'r-', label='Air Pressure [Pa]') +#plt.plot(t_list, Temp_list, 'b-', label='Air Temperature [K]') +#plt.plot(t_list, rho_list, 'g-', label='Air Density in [kg/m$^3$]') +plt.xlabel('time in s') +plt.ylabel('Balloon Altitude in m') +plt.legend() +plt.show() +#"""