changed basic thermal model, now realistic behavior

main.py

@@ -20,9 +20,6 @@ lon = start_lon   # deg
 
 date = Time('2020-12-13 12:00:00.000')
 
-AZ = sun_angles_astropy(45.0, 45.0, 0, date)[0]
-ELV = sun_angles_astropy(45.0, 45.0, 0, date)[1]
-
 
 # INITIALISATION
 
@@ -45,9 +42,9 @@ while t <= t_end and h >= 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