design freeze

air_data.py

@@ -2,7 +2,10 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import sys
+import time
+import aerosandbox.library.atmosphere as atmo
 from scipy.interpolate import interp1d
+from lmfit import Model, Parameter
 
 def get_air_data(alt):
     """
@@ -12,63 +15,171 @@ def get_air_data(alt):
     df = pd.read_csv("air_data.csv")
     h = df['h']
     a = df['a']
-    rho = df['rho/rho0'] * 1.225
-    nu = df['nu']
-    mu = nu * rho
+    rho = df['rho']
+    mu = df['mu']*1e5
     # linear for a, mu; quadratic for rho, nu
     f_rho = interp1d(h, rho, kind="quadratic")
     f_mu = interp1d(h, mu)
     f_a = interp1d(h, a)
     return np.round(f_rho(alt),4), np.round(f_mu(alt),4), np.round(f_a(alt),1)
 
-    # i = h[h == alt].index[0]
-    # return rho[i], mu[i], a[i]
+def get_atmo_data(alt):
+    """
+    Peter's fitted atmosphere model good for <40km altitude
+    ~1000x faster than interpolating
+    """
+    P = atmo.get_pressure_at_altitude(alt)
+    rho = atmo.get_density_at_altitude(alt)
+    T = atmo.get_temperature_at_altitude(alt)
+    mu = atmo.get_viscosity_from_temperature(T)*1e5
+    a = atmo.get_speed_of_sound_from_temperature(T)
+    return rho, mu, a
+
+# def fitted_air_data(alt):
+#     if alt < 11000:
+#         T = 273.15 + 15.04 - 0.00649*alt
+#         p = 101.29 * (T/288.08) ** 5.256
+#     elif alt < 25000:
+#         T = 273.15 - 56.46
+#         p = 22.65 * np.exp(1.73 - 0.000157*alt)
+#     else:
+#         T = 273.15 -131.21 + 0.00299*alt
+#         p = 2.488 * (T/216.6) ** -11.388
+#     rho = p / (.2869 * T)
+#     a = np.sqrt(1.4 * 287 * T)
+
+def fit_T(alt, b1, c1, b2, b3, c3, b4, c4):
+    T = np.zeros(alt.shape)
+    for i, h in enumerate(alt):
+        if h <= 11000:
+            T[i] = b1 + c1*h
+        elif h <= 20000:
+            T[i]= b2
+        elif h <= 32000:
+            T[i] = b3 + c3*h
+        else:
+            T[i] = b4 + c4*h
+    return T
+
+def fit_p(alt, b1, c1, exp1, b2, c2, exp2, b3, c3, exp3, b4, c4, exp4):
+    p = np.zeros(alt.shape)
+    for i, h in enumerate(alt):
+        if h <= 11000:
+            T = 288.15 - 0.0065*h
+            p[i] = b1 * (T/c1) ** exp1
+        elif h <= 20000:
+            T = 216.65
+            p[i] = b2 * np.exp(c2 + exp2*h)
+        elif h <= 32000:
+            T = 196.65 + 0.001*h
+            p[i] = b3 * (T/c3) ** exp3
+        else:
+            T = 148.236 + 0.002554*h
+            p[i] = b4 * (T/c4) ** exp4
+    return p
+
+def fit_rho(alt, crho):
+    rho = np.zeros(alt.shape)
+    for i, h in enumerate(alt):
+        if h < 11000:
+            T = 273.15 + 15.04 - 0.00649*h
+            p = 101.29 * (T/288.08) ** 5.256
+        elif h < 20000:
+            T = 273.15 - 56.46
+            p = 22.65 * np.exp(1.73 - 0.000157*h)
+        elif h < 32000:
+            T = 273.15 -131.21 + 0.00299*h
+            p = 2.488 * (T/216.6) ** -11.388
+        else:
+            pass
+        rho[i] = p / (crho * T)
+        a = np.sqrt(1.4 * 287 * T)
+    return rho
 
 def change_air_data(alt):
     """
     Rewrite qcon.def which specifies air_data for qmil and qprop to the given
     altitude
     """
-    rho, mu, a = get_air_data(alt)
+    rho, mu, a = get_atmo_data(alt)
     # print(rho, mu, a)
     data = []
     data.append(str(rho) + " ! rho  kg/m^3; alt: " + str(alt) + "m\n" )
     data.append(str(mu) + "E-5 ! mu   kg/m-s\n")
     data.append(str(a) + " ! a    m/s")
     with open('qcon.def', 'w') as file:
         file.writelines(data)
+    return rho, mu, a
 
 def plot_air_data():
     """
     Plot air density, dynamic viscosity, and speed of sound as a function of
     altitude, extrapolating data from air_data.csv
     """
-    alts = np.arange(-2000, 50000, 1000)
+    alts = np.arange(0, 50000, 1000)
     rhos = np.zeros(len(alts))
     mus = np.zeros(len(alts))
     ass = np.zeros(len(alts))
+    rhos2 = np.zeros(len(alts))
+    mus2 = np.zeros(len(alts))
+    ass2 = np.zeros(len(alts))
     for n, h in enumerate(alts):
         rhos[n], mus[n], ass[n] = get_air_data(h)
+        rhos2[n], mus2[n], ass2[n] = get_atmo_data(h)
 
     fig = plt.figure(figsize=(10,5))
     plt.title("Standard Atmosphere")
     plt.subplot(131)
     plt.plot(rhos, alts/1000)
+    plt.plot(rhos2, alts/1000)
     plt.axhline(color="0.5", linestyle=":")
     plt.axvline(color="0.5", linestyle=":")
     plt.ylabel("Altitude (km)")
     plt.xlabel(r"$\rho$ (kg/$m^3$)")
     plt.subplot(132)
     plt.plot(mus, alts/1000)
+    plt.plot(mus2, alts/1000)
     plt.xlabel(r"$\mu$")
     plt.axhline(color="0.5", linestyle=":")
     plt.subplot(133)
     plt.plot(ass, alts/1000)
+    plt.plot(ass2, alts/1000)
     plt.xlabel("a (m/s)")
     plt.axhline(color="0.5", linestyle=":")
     plt.show()
 
 if __name__ == "__main__":
-    alt = 22000 if len(sys.argv) <= 1 else sys.argv[1]
+    alt = 22000 if len(sys.argv) <= 1 else int(sys.argv[1])
     change_air_data(alt)
     # plot_air_data()
+    # df = pd.read_csv("air_data.csv")
+    # h = df['h'].to_numpy()
+    # a = df['a'].to_numpy()
+    # rho = df['rho'].to_numpy()
+    # mu = df['mu'].to_numpy()
+    # p = df['p'].to_numpy()
+    # T = df['T'].to_numpy()
+
+    # start = time.time()
+    # for i in range(10000):
+    #     get_atmo_data(20000+2*i)
+    # print(time.time()-start)
+
+    # T_model = Model(fit_T)
+    # T_params = T_model.make_params(b1=288.15, c1=-0.0065, b2=216.65, b3=196.65, c3=0.001, b4=149.4, c4=0.0025)
+    # T_result = T_model.fit(T, params, alt=h)
+    # print(T_result.fit_report())
+    # Rsquared = 1 - result.residual.var()/np.var(T)
+    # print(Rsquared)
+    # p_model = Model(fit_p)
+    # p_params = p_model.make_params(b1=101300, c1=288, exp1=5.26, b2=19900, c2=1.9, exp2=-0.00016,
+    #                                b3=2000, c3=220, exp3=-11, b4=1500, c4=180, exp4=-13)
+    # p_result = p_model.fit(p, p_params, alt=h)
+    # print(p_result.fit_report())
+    # print(1 - p_result.residual.var()/np.var(p))
+    # rho_model = Model(fit_rho)
+    # params = rho_model.make_params(crho=1)
+    # result = T_model.fit(T, params, alt=h)
+    # print(result.fit_report())
+    # Rsquared = 1 - result.residual.var()/np.var(T)
+    # print(Rsquared)

best_prop

@@ -6,41 +6,41 @@ Template prop lrp01.dat M0.3
    0.3800  6.2832    ! CL0    CL_a 
   -0.4000  1.2900    ! CLmin  CLmax
 
-   0.01140  0.03400  0.02400  0.5226    ! CD0    CD2u  CD2l  CLCD0
+   0.01140  0.03400  0.09000  0.5226    ! CD0    CD2u  CD2l  CLCD0
    100000.0  -0.500            ! REref  REexp
 
    1.0000  1.0000  1.0000   !  Rfac   Cfac   Bfac
    0.0000  0.0000  0.0000   !  Radd   Cadd   Badd
 
 #        r            c        beta
-   0.16647      0.80966E-01  65.0400
-   0.19940      0.93242E-01  60.5287
-   0.23233      0.10344      56.3969
-   0.26527      0.11167      52.6272
-   0.29820      0.11809      49.1977
-   0.33113      0.12292      46.0820
-   0.36407      0.12642      43.2516
-   0.39700      0.12879      40.6782
-   0.42993      0.13024      38.3351
-   0.46287      0.13094      36.1976
-   0.49580      0.13101      34.2436
-   0.52873      0.13055      32.4535
-   0.56167      0.12962      30.8100
-   0.59460      0.12826      29.2979
-   0.62753      0.12651      27.9034
-   0.66047      0.12440      26.6142
-   0.69340      0.12194      25.4196
-   0.72633      0.11913      24.3102
-   0.75927      0.11596      23.2777
-   0.79220      0.11240      22.3148
-   0.82513      0.10842      21.4150
-   0.85807      0.10397      20.5729
-   0.89100      0.98984E-01  19.7832
-   0.92393      0.93364E-01  19.0417
-   0.95687      0.86984E-01  18.3444
-   0.98980      0.79658E-01  17.6877
-    1.0227      0.71094E-01  17.0685
-    1.0557      0.60775E-01  16.4840
-    1.0886      0.47589E-01  15.9317
-    1.1215      0.27757E-01  15.4092
-    1.1380      0.15349E-01  15.1592
+   0.11633      0.77139E-01  72.3249
+   0.14900      0.10000      66.6475
+   0.18167      0.12314      61.3643
+   0.21433      0.14507      56.5582
+   0.24700      0.16344      52.2776
+   0.27967      0.17693      48.4979
+   0.31233      0.18518      45.1708
+   0.34500      0.18842      42.2456
+   0.37767      0.18735      39.6728
+   0.41033      0.18293      37.4053
+   0.44300      0.17618      35.3993
+   0.47567      0.16805      33.6145
+   0.50833      0.15933      32.0142
+   0.54100      0.15062      30.5648
+   0.57367      0.14231      29.2378
+   0.60633      0.13444      28.0156
+   0.63900      0.12698      26.8846
+   0.67167      0.11989      25.8330
+   0.70433      0.11312      24.8507
+   0.73700      0.10661      23.9288
+   0.76967      0.10030      23.0597
+   0.80233      0.94093E-01  22.2366
+   0.83500      0.87912E-01  21.4539
+   0.86767      0.81643E-01  20.7065
+   0.90033      0.75143E-01  19.9899
+   0.93300      0.68217E-01  19.3002
+   0.96567      0.60573E-01  18.6341
+   0.99833      0.51712E-01  17.9886
+    1.0310      0.40600E-01  17.3612
+    1.0637      0.23844E-01  16.7495
+    1.0800      0.13349E-01  16.4495

climb.py

@@ -100,8 +100,8 @@ def plot_trajectory(data_file="cycle.npz"):
     day = len(t) * 18 // 24
     res_day = np.concatenate((result[:,:night], result[:,day:]), axis=1)
     res_night = result[:,night:day]
-    t_d, h_d, v_d, p_d, thrust_d, rpm_d, Q_d, Pshaft_d, J_d, dbeta_d, eta_d, eta_tot_d = res_day
-    t_n, h_n, v_n, p_n, thrust_n, rpm_n, Q_n, Pshaft_n, J_n, dbeta_n, eta_n, eta_tot_n = res_night
+    print("Avg RPM", np.mean(rpm), "PShaft", np.mean(Pshaft), "Q", np.mean(Q))
+    print("Max RPM", max(rpm), "PShaft", max(Pshaft), "Q", max(Q))
 
     fig = plt.figure(figsize=(7,7.5))
     ax1 = plt.subplot(311)
@@ -159,19 +159,19 @@ def plot_trajectory(data_file="cycle.npz"):
     plt.show()
 
 if __name__ == "__main__":
-    num_motor = 6
+    num_motor = 4
     data =  np.load('climb_path.npz')
-    ts = data['t'][:55]/3600
-    hs = data['h'][:55]
-    vs = data['v'][:55]
-    thrusts = data['thrust'][:55]/num_motor
+    ts = data['t'][:41]/3600
+    hs = data['h'][:41]
+    vs = data['v'][:41]
+    thrusts = data['thrust'][:41]/num_motor
 
     # start = time.time()
-    eff_opt = follow_trajectory(ts, hs, vs, thrusts, npt=200, optimize=True)
-    print("Average Efficiency:", eff_opt)
-    plot_trajectory("climb.npz")
+    # eff_opt = follow_trajectory(ts, hs, vs, thrusts, npt=200, optimize=True)
+    # print("Var Pitch Average Efficiency:", eff_opt)
+    # plot_trajectory("climb.npz")
     # print(time.time() - start)
     eff_unopt = follow_trajectory(ts, hs, vs, thrusts, npt=200, optimize=False)
-    print("Average Efficiency:", eff_unopt)
+    print("Fixed Pitch Average Efficiency:", eff_unopt)
     plot_trajectory("climb_unopt.npz")
     # print(time.time() - start)

descent.py

@@ -65,16 +65,18 @@ def plot_descent(result):
 
 
 if __name__ == "__main__":
+    avg_lift = 390 * 9.8
+    num_motor = 4
     # run_qprop(32.01, -10)
     result0 = windmill(0, 20000, 0)
     result5 = windmill(0, 20000, -5)
     result10 = windmill(0, 20000, -10)
-    result15 = windmill(0, 20000, -15)
+    # result15 = windmill(0, 20000, -15)
     plt.figure()
-    plt.plot(result0[:,0]/1e3, -result0[:,7]/500, label=r"$d\beta = 0\degree$")
-    plt.plot(result5[:,0]/1e3, -result5[:,7]/500, label=r"$d\beta = -5\degree$")
-    plt.plot(result10[:,0]/1e3, -result10[:,7]/500, label=r"$d\beta = -10\degree$")
-    plt.plot(result15[:,0]/1e3, -result15[:,7]/500, label=r"$d\beta = -15\degree$")
+    plt.plot(result0[:,0]/1e3, -result0[:,7]/avg_lift*num_motor, label=r"$d\beta = 0\degree$")
+    plt.plot(result5[:,0]/1e3, -result5[:,7]/avg_lift*num_motor, label=r"$d\beta = -5\degree$")
+    plt.plot(result10[:,0]/1e3, -result10[:,7]/avg_lift*num_motor, label=r"$d\beta = -10\degree$")
+    # plt.plot(result15[:,0]/1e3, -result15[:,7]/avg_lift*num_motor, label=r"$d\beta = -15\degree$")
     plt.xlabel("Altitude [km]")
     plt.ylabel("Windmill Drag / Average Lift")
     plt.title("Total Windmill Drag for Descent")
@@ -83,18 +85,18 @@ def plot_descent(result):
     plt.plot(result0[:,0]/1e3, result0[:,2], label=r"$d\beta = 0\degree$")
     plt.plot(result5[:,0]/1e3, result5[:,2], label=r"$d\beta = -5\degree$")
     plt.plot(result10[:,0]/1e3, result10[:,2], label=r"$d\beta = -10\degree$")
-    plt.plot(result15[:,0]/1e3, result15[:,2], label=r"$d\beta = -15\degree$")
+    # plt.plot(result15[:,0]/1e3, result15[:,2], label=r"$d\beta = -15\degree$")
     plt.xlabel("Altitude [km]")
     plt.ylabel("Blade CL")
-    plt.title("Total Windmill Drag for Descent")
+    plt.title("Blade Windmill cl for Descent")
     plt.legend()
     plt.figure()
     plt.plot(result0[:,0]/1e3, result0[:,3], label=r"$d\beta = 0\degree$")
     plt.plot(result5[:,0]/1e3, result5[:,3], label=r"$d\beta = -5\degree$")
     plt.plot(result10[:,0]/1e3, result10[:,3], label=r"$d\beta = -10\degree$")
-    plt.plot(result15[:,0]/1e3, result15[:,3], label=r"$d\beta = -15\degree$")
+    # plt.plot(result15[:,0]/1e3, result15[:,3], label=r"$d\beta = -15\degree$")
     plt.xlabel("Altitude [km]")
     plt.ylabel("Blade CD")
-    plt.title("Total Windmill Drag for Descent")
+    plt.title("Blade Windmill cd for Descent")
     plt.legend()
     plt.show()

est_motor

@@ -2,6 +2,6 @@ Estimated Motor Model
 
  1     ! motor type (brushed DC)
 
- 1.1      ! Rmotor (Ohms)
- 3.6      ! Io     (Amps)
- 8.0      ! Kv     (rpm/Volt)
+ 0.3      ! Rmotor (Ohms)
+ 1.1      ! Io     (Amps)
+ 4.0      ! Kv     (rpm/Volt)