Finish final edits and documentation

README.md

@@ -1,20 +1,20 @@
 # prop-design
-Python wrapped QMil and QProp used to design propellers for HALE UAV in 16.82 Flight Vehicle Development
+Python wrapped QMil and QProp used to design propellers for HALE UAV in 16.82 Flight Vehicle Development for the Solar HALE UAV, Dawn.
 
 ## How to Use
-Important functions are sectioned off into different files.
+Important functions for prop design are sectioned off into different files.
 
 ### air_data.py
-This file is used to modify qcon.def which specifies air data for Qmil and Qprop and accepts an altitude when run in the terminal as an argument. Data is pulled from air_data.csv and interpolated.
+This file is used to modify qcon.def which specifies air data for Qmil and Qprop and accepts an altitude when run in the terminal as an argument. There are two options for data -- to pull from air_data.csv which is more accurate but slower, or to use aerosandbox curve fit atmospheric model which is faster. The latter is the default.
 
 ### qmil_design.py
-This file takes template.mil and optimizes over rpm for max efficiency to design a prop at the given conditions. The final .mil file used is saved in output.mil and the prop is saved in best_prop
+This file takes template.mil and optimizes over rpm for max efficiency to design a prop at the given conditions. The final .mil file used is saved in output.mil and the prop is saved in best_prop. Plots propeller geometry
 
 ### trajectory.py
 This file takes in a 24 hour cycle trajectory and computes various propulsive data over the cycle. The data is pulled from time_altitude_speed.npz which contains arrays for time, altitude, velocity, and thrust
 
-### qprop_sweep.py
-This file sweeps over propulsive area, altitude, airspeed, thrust, and variable pitch angle to calculate efficiency and shaft power at each condition. This data is then saved to opt_sweep.npz
+### climb.py
+This file works similar to trajectory.py. It requires a climb_path.npz input and runs QProp for the climb trajectory.
 
-### curve_fit.py
-This file takes the sweep data from qprop_sweep and attempts to fit a curve to the data.
+### descent.py
+The descend() function calculates a descent trajectory with windmilling props. See docstring and final memo for detailed information.

climb.py

@@ -96,54 +96,40 @@ def plot_trajectory(data_file="cycle.npz"):
     """
     result = np.load(data_file)["res"].T
     t, h, v, p, thrust, rpm, Q, Pshaft, J, dbeta, eta, eta_tot = result
-    night = int( len(t) * 7.5 // 24 )
-    day = len(t) * 18 // 24
-    res_day = np.concatenate((result[:,:night], result[:,day:]), axis=1)
-    res_night = result[:,night:day]
     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)
-    ax1.set_title("Climb for each Propeller")
-    # ax1.plot(t_d, h_d/304.88, ".-", label="Day", color="cornflowerblue")
-    # ax1.plot(t_n, h_n/304.88, ".-", label="Night", color="midnightblue")
+    fig = plt.figure(figsize=(12,7.5))
+    ax1 = plt.subplot(321)
+    # fig.suptitle("Climb for each Propeller")
     ax1.plot(t, h/304.88, ".-", color="cornflowerblue")
     ax1.set_ylabel('Altitude [kft]', color="midnightblue")
-    ax12 = ax1.twinx()
-    # ax12.plot(t_d, rpm_d, ".-", label="Day", color="indianred")
-    # ax12.plot(t_n, rpm_n, ".-", label="Night", color="maroon")
+    # ax12 = ax1.twinx()
+    ax12 = plt.subplot(322)
     ax12.plot(t, rpm, ".-", color="indianred")
     ax12.set_ylabel("RPM", color="maroon")
-    ax12.grid(None)
-    # ax12.set_yticks(np.linspace(ax12.get_yticks()[0], ax12.get_yticks()[-1], len(ax1.get_yticks())))
+    # ax12.grid(None)
 
-    ax2 = plt.subplot(312)
-    # ax2.plot(t_d, thrust_d, ".-", label="Day", color="cornflowerblue")
-    # ax2.plot(t_n, thrust_n, ".-", label="Night", color="midnightblue")
+    ax2 = plt.subplot(323)
     ax2.plot(t, thrust, ".-", color="cornflowerblue")
     ax2.set_ylabel('Thrust [N]', color="midnightblue")
-    ax22 = ax2.twinx()
-    # ax22.plot(t_d, eta_d, ".-", label="Day", color="indianred")
-    # ax22.plot(t_n, eta_n, ".-", label="Night", color="maroon")
+    # ax22 = ax2.twinx()
+    ax22 = plt.subplot(324)
     ax22.plot(t, eta, ".-", color="indianred")
     ax22.set_ylabel(r'$\eta$', color="maroon")
-    ax22.grid(None)
+    # ax22.grid(None)
 
-    ax3 = plt.subplot(313)
-    # plt.plot(t, dbeta)
-    # plt.ylabel(r'$d\beta$')
-    # ax3.plot(t_d, Pshaft_d, ".-", label="Day", color="cornflowerblue")
-    # ax3.plot(t_n, Pshaft_n, ".-", label="Night", color="midnightblue")
-    ax3.plot(t, Pshaft, ".-", color="cornflowerblue")
-    ax3.set_ylabel('Shaft Power [W]', color="midnightblue")
-    ax32 = ax3.twinx()
-    # ax32.plot(t_d, Q_d, ".-", label="Day", color="indianred")
-    # ax32.plot(t_n, Q_n, ".-", label="Night", color="maroon")
+    ax3 = plt.subplot(325)
+    ax3.plot(t, Pshaft, ".-", color="indianred")
+    ax3.set_ylabel('Shaft Power [W]', color="maroon")
+    # ax32 = ax3.twinx()
+    ax32 = plt.subplot(326)
     ax32.plot(t, Q, ".-", color="indianred")
     ax32.set_ylabel('Required Torque [Nm]', color="maroon")
     ax3.set_xlabel("Time after Takeoff [hr]")
-    ax32.grid(None)
+    ax32.set_xlabel("Time after Takeoff [hr]")
+    plt.tight_layout()
+    # ax32.grid(None)
 
     plt.figure()
     plt.plot(J, eta, ".-")

descent.py

@@ -41,6 +41,9 @@ def run_qprop(vel, dbeta=0, prop="best_prop", motor="est_motor"):
     return data, cl, cd
 
 def windmill(alt_min, alt_max, dbeta):
+    """
+    Calculate windmill drag at specified dbeta between alt_min and alt_max
+    """
     q = 45
     num_alts = 50
     alts = np.linspace(alt_min, alt_max, num_alts)
@@ -64,6 +67,13 @@ def plot_descent(result):
     plt.show()
 
 def descend():
+    """
+    Function to calculate descent trajectory. Descent at .85 divergence speed.
+    Recalculates CL then scales CD for new induced drag and prop windmill drag
+    is added. The descent angle is calculated and the speed is integrated
+    using Euler method at 5 minute intervals and stores data in results.
+    Plots descent trajectory
+    """
     W = 390
     g = 9.81
     V_div = 13.3 # Divergence speed at sea level
@@ -79,7 +89,7 @@ def descend():
     CD0 = CL_cruise / LD_cruise - CL_cruise**2 / (np.pi * AR * e) # 0.0152
     CD = CD0 + CL ** 2 / (np.pi * AR * e) # Descent CD with less induced drag
     drag_vehicle = CD * q * S
-    dt = 5 * 60 # 30 minutes [sec]
+    dt = 5 * 60 # 5 minutes [sec]
 
     alt, x, t = (19800, 0, 0)
     result = []

props/final_CDR_prop_0507

@@ -0,0 +1,46 @@
+
+Template prop lrp01.dat M0.3                    
+
+   2       ! Nblades
+
+   0.3800  6.2832    ! CL0    CL_a 
+  -0.4000  1.2900    ! CLmin  CLmax
+
+   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.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

qmil_design.py

@@ -73,6 +73,10 @@ def change_prop_area(area):
     design_opt_rpm()
 
 def plot_prop(propfile):
+    """
+    Plot the propeller r/c and beta/c distribution and the actual geometry
+    :params propfile: Filename or path to file containing propeller geometry
+    """
     f = open(propfile, "r")
     contents = f.readlines()
     prop_geom =  contents[-32:]
@@ -156,7 +160,7 @@ def design_opt_rpm(h=21000, plot=False, traj=False, opt=False):
     # plot_blade_cl()
     # print(-design_prop(1200, "test_prop"))
 
-
+    # Sensitivity study to design RPM
     # rpms = np.arange(600.0, 2500.0, 100.0)
     # etas = np.zeros(len(rpms))
     #

trajectory.py

@@ -105,41 +105,44 @@ def plot_trajectory(data_file="cycle.npz"):
     print("Avg RPM", rpm_n[0], "PShaft", Pshaft_n[0], "Q", Q_n[0])
     print("Max RPM", max(rpm_d), "PShaft", max(Pshaft_d), "Q", max(Q_d))
 
-    fig = plt.figure(figsize=(7,7.5))
-    ax1 = plt.subplot(311)
-    ax1.set_title("24 Hour Cycle for each Propeller")
+    fig = plt.figure(figsize=(12,7.5))
+    ax1 = plt.subplot(321)
+    fig.suptitle("24 Hour Cycle for each Propeller")
     ax1.plot(t_d, h_d/304.88, ".-", label="Day", color="cornflowerblue")
     ax1.plot(t_n, h_n/304.88, ".-", label="Night", color="midnightblue")
     ax1.set_ylabel('Altitude [kft]', color="midnightblue")
-    ax12 = ax1.twinx()
+    ax1.legend()
+    ax12 = plt.subplot(322)
     ax12.plot(t_d, rpm_d, ".-", label="Day", color="indianred")
     ax12.plot(t_n, rpm_n, ".-", label="Night", color="maroon")
     ax12.set_ylabel("RPM", color="maroon")
-    ax12.grid(None)
+    ax12.legend()
+    # ax12.grid(None)
     # ax12.set_yticks(np.linspace(ax12.get_yticks()[0], ax12.get_yticks()[-1], len(ax1.get_yticks())))
 
-    ax2 = plt.subplot(312)
+    ax2 = plt.subplot(323)
     ax2.plot(t_d, thrust_d, ".-", label="Day", color="cornflowerblue")
     ax2.plot(t_n, thrust_n, ".-", label="Night", color="midnightblue")
     ax2.set_ylabel('Thrust [N]', color="midnightblue")
-    ax22 = ax2.twinx()
+    ax22 = plt.subplot(324)
     ax22.plot(t_d, eta_d, ".-", label="Day", color="indianred")
     ax22.plot(t_n, eta_n, ".-", label="Night", color="maroon")
     ax22.set_ylabel(r'$\eta$', color="maroon")
-    ax22.grid(None)
+    # ax22.grid(None)
 
-    ax3 = plt.subplot(313)
+    ax3 = plt.subplot(325)
     # plt.plot(t, dbeta)
     # plt.ylabel(r'$d\beta$')
-    ax3.plot(t_d, Pshaft_d, ".-", label="Day", color="cornflowerblue")
-    ax3.plot(t_n, Pshaft_n, ".-", label="Night", color="midnightblue")
-    ax3.set_ylabel('Shaft Power [W]', color="midnightblue")
-    ax32 = ax3.twinx()
+    ax3.plot(t_d, Pshaft_d, ".-", label="Day", color="indianred")
+    ax3.plot(t_n, Pshaft_n, ".-", label="Night", color="maroon")
+    ax3.set_ylabel('Shaft Power [W]', color="maroon")
+    ax32 = plt.subplot(326)
     ax32.plot(t_d, Q_d, ".-", label="Day", color="indianred")
     ax32.plot(t_n, Q_n, ".-", label="Night", color="maroon")
     ax32.set_ylabel('Required Torque [Nm]', color="maroon")
     ax3.set_xlabel("Time from Solar Noon [hr]")
-    ax32.grid(None)
+    ax32.set_xlabel("Time from Solar Noon [hr]")
+    # ax32.grid(None)
 
     plt.figure()
     plt.plot(J, eta, ".-")
@@ -212,6 +215,7 @@ def rpm_trade(rpms, ts, hs, vs, thrusts, optimize):
     # plt.plot(ts, hs)
     # plt.show()
 
+    # For motor study
     # plot_motor_eff(8*np.pi/30, 0.3, 1.1)
     # Kvs = np.linspace(3, 10, 10)
     # I0s = np.linspace(1, 5, 15)
@@ -228,19 +232,6 @@ def rpm_trade(rpms, ts, hs, vs, thrusts, optimize):
     # plt.colorbar()
     # plt.show()
 
-    # rpms = np.arange(800, 1700, 50)
-    # opt_eff = rpm_trade(rpms, ts, hs, vs, thrusts, True)
-    # unopt_eff = rpm_trade(rpms, ts, hs, vs, thrusts, False)
-    # np.savez("rpm_trade.npz", opt_eff, unopt_eff)
-    # plt.figure()
-    # plt.plot(rpms, opt_eff, label="Variable Pitch")
-    # plt.plot(rpms, unopt_eff, label="Fixed Pitch")
-    # plt.title("Design RPM vs Average Efficiency")
-    # plt.xlabel("QMil Design RPM")
-    # plt.ylabel(r"Average $\eta$ over 24 Hour Cycle")
-    # plt.legend()
-    # plt.show()
-
     # start = time.time()
     # eff_opt = follow_trajectory(ts, hs, vs, thrusts, npt=200, optimize=True)
     # print("Average Efficiency Var Pitch:", eff_opt)