First commit

air_data.csv

@@ -0,0 +1,51 @@
+h,T,p,rho/rho0,nu,k,a
+-2000,301.2,1.2778,1.2067,1.253,2.636,347.9
+-1500,297.9,1.207,1.1522,1.301,2.611,346
+-1000,294.7,1.1393,1.0996,1.352,2.585,344.1
+-500,291.4,1.0748,1.0489,1.405,2.56,342.2
+0,288.15,1.01325,1,1.461,2.534,340.3
+500,284.9,0.9546,0.9529,1.52,2.509,338.4
+1000,281.7,0.8988,0.9075,1.581,2.483,336.4
+1500,278.4,0.8456,0.8638,1.646,2.457,334.5
+2000,275.2,0.795,0.8217,1.715,2.431,332.5
+2500,271.9,0.7469,0.7812,1.787,2.405,330.6
+3000,268.7,0.7012,0.7423,1.863,2.379,328.6
+3500,265.4,0.6578,0.7048,1.943,2.353,326.6
+4000,262.2,0.6166,0.6689,2.028,2.327,324.6
+4500,258.9,0.5775,0.6343,2.117,2.301,322.6
+5000,255.7,0.5405,0.6012,2.211,2.275,320.5
+5500,252.4,0.5054,0.5694,2.311,2.248,318.5
+6000,249.2,0.4722,0.5389,2.416,2.222,316.5
+6500,245.9,0.4408,0.5096,2.528,2.195,314.4
+7000,242.7,0.4111,0.4817,2.646,2.169,312.3
+7500,239.5,0.383,0.4549,2.771,2.142,310.2
+8000,236.2,0.3565,0.4292,2.904,2.115,308.1
+8500,233,0.3315,0.4047,3.046,2.088,306
+9000,229.7,0.308,0.3813,3.196,2.061,303.8
+9500,226.5,0.2858,0.3589,3.355,2.034,301.7
+10000,223.3,0.265,0.3376,3.525,2.007,299.8
+10500,220,0.2454,0.3172,3.706,1.98,297.4
+11000,216.8,0.227,0.2978,3.899,1.953,295.2
+11500,216.7,0.2098,0.2755,4.213,1.952,295.1
+12000,216.7,0.194,0.2546,4.557,1.952,295.1
+12500,216.7,0.1793,0.2354,4.93,1.952,295.1
+13000,216.7,0.1658,0.2176,5.333,1.952,295.1
+13500,216.7,0.1533,0.2012,5.768,1.952,295.1
+14000,216.7,0.1417,0.186,6.239,1.952,295.1
+14500,216.7,0.131,0.172,6.749,1.952,295.1
+15000,216.7,0.1211,0.159,7.3,1.952,295.1
+15500,216.7,0.112,0.147,7.895,1.952,295.1
+16000,216.7,0.1035,0.1359,8.54,1.952,295.1
+16500,216.7,0.09572,0.1256,9.237,1.952,295.1
+17000,216.7,0.0885,0.1162,9.99,1.952,295.1
+17500,216.7,0.08182,0.1074,10.805,1.952,295.1
+18000,216.7,0.07565,0.0993,11.686,1.952,295.1
+18500,216.7,0.06995,0.09182,12.639,1.952,295.1
+19000,216.7,0.06467,0.08489,13.67,1.952,295.1
+19500,216.7,0.0598,0.0785,14.784,1.952,295.1
+20000,216.7,0.05529,0.07258,15.989,1.952,295.1
+22000,218.6,0.04047,0.05266,22.201,1.968,296.4
+24000,220.6,0.02972,0.03832,30.743,1.985,297.7
+26000,222.5,0.02188,0.02797,42.439,2.001,299.1
+28000,224.5,0.01616,0.02047,58.405,2.018,300.4
+30000,226.5,0.01197,0.01503,80.134,2.034,301.7

air_data.py

@@ -0,0 +1,35 @@
+import pandas as pd
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+import numpy as np
+import sys
+
+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
+    # 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 change_air_data(alt):
+    rho, mu, a = get_air_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)
+
+if __name__ == "__main__":
+    alt = 25000 if len(sys.argv) <= 1 else sys.argv[1]
+    change_air_data(alt)

best_prop

@@ -0,0 +1,46 @@
+
+Template prop                                   
+
+   2       ! Nblades
+
+   0.4000  6.2832    ! CL0    CL_a 
+  -0.4000  1.3000    ! CLmin  CLmax
+
+   0.02000  0.00800  0.00600  0.3000    ! 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.16917      0.26251      73.7697
+   0.20750      0.31523      70.0227
+   0.24583      0.36321      66.4530
+   0.28417      0.40570      63.0727
+   0.32250      0.44223      59.8872
+   0.36083      0.47267      56.8969
+   0.39917      0.49711      54.0980
+   0.43750      0.51588      51.4837
+   0.47583      0.52941      49.0450
+   0.51417      0.53820      46.7718
+   0.55250      0.54277      44.6533
+   0.59083      0.54362      42.6787
+   0.62917      0.54121      40.8371
+   0.66750      0.53594      39.1182
+   0.70583      0.52816      37.5123
+   0.74417      0.51815      36.0102
+   0.78250      0.50614      34.6033
+   0.82083      0.49229      33.2839
+   0.85917      0.47671      32.0448
+   0.89750      0.45943      30.8794
+   0.93583      0.44045      29.7817
+   0.97417      0.41969      28.7464
+    1.0125      0.39699      27.7683
+    1.0508      0.37206      26.8431
+    1.0892      0.34451      25.9666
+    1.1275      0.31366      25.1352
+    1.1658      0.27845      24.3454
+    1.2042      0.23690      23.5941
+    1.2425      0.18476      22.8786
+    1.2808      0.10742      22.1963
+    1.3000      0.59302E-01  21.8676

est_motor

@@ -0,0 +1,7 @@
+Estimated Motor Model
+
+ 1     ! motor type (brushed DC)
+
+ 0.2      ! Rmotor (Ohms)
+ 0.5      ! Io     (Amps)
+ 20.0     ! Kv     (rpm/Volt)

output.mil

@@ -0,0 +1,25 @@
+
+Template prop
+
+   2       ! Nblades
+
+   0.4000  6.2832    ! CL0    CL_a
+  -0.4000  1.3000    ! CLmin  CLmax
+
+   0.02000   0.008  0.006  0.30  ! CD0    CD2u   CD2l  CLCD0
+   100000.0 -0.500               ! REref  REexp
+
+
+   0.0  0.5  1.0   ! XIdes   (r/R locations where design cl is specified)
+   0.6  0.5  0.4   ! CLdes   (specified cl)
+
+  0.15    !  hub radius(m)
+  1.30    !  tip radius(m)
+  35.0    !  speed(m/s)
+  713.1864009684666 ! rpm
+
+  0.0      !  Thrust(N)   ( 0 if power  specified )
+    1333.0 !  Power(W)    ( 0 if thrust specified )
+
+ 0    0.0   ! Ldes   [ KQdes ]
+   30       !  Nout     number of output stations (optional)

plot_mil.py

@@ -0,0 +1,39 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy import radians as r
+
+CL0 = 0.4
+CL_a = 6.2832
+CLmin = -0.4
+CLmax = 1.3
+CD0 = 0.02
+CD2u = 0.008
+CD2l = 0.006
+CLCD0 = 0.30
+REref = 100000
+REexp = -0.5
+alphas = np.linspace(r(-10), r(10), 200)
+v = 20
+a = 298.4
+Re = 150000
+beta = (1-(v/a)**2)**(1/2)
+CL = np.zeros(200)
+CD = np.zeros(200)
+for i, alp in enumerate(alphas):
+    CL[i] = (alp*CL_a + CL0)/beta
+    stall = False
+    if CL[i] > CLmax:
+        stall = True
+        CL[i] = CLmax
+    elif CL[i] < CLmin:
+        stall = True
+        CL[i] = CLmin
+    CD2 = CD2u if CL[i] > CLCD0 else CD2l
+    CD[i] = (CD0 + CD2*(CL[i]-CLCD0)**2) * (Re/REref)**REexp
+    if stall:
+        CD[i] += 2*(np.sin(alp - (CLCD0-CL0)/CL_a))**2
+
+plt.plot(alphas, CL)
+plt.show()
+plt.plot(alphas, CD)
+plt.show()

qcon.def

@@ -0,0 +1,3 @@
+0.0401 ! rho  kg/m^3; alt: 25000m
+1.4486E-5 ! mu   kg/m-s
+298.4 ! a    m/s

qmil_design.py

@@ -0,0 +1,110 @@
+import subprocess
+import numpy as np
+import re
+import matplotlib.pyplot as plt
+from scipy.optimize import minimize_scalar
+
+def get_num(line):
+    """
+    Regex function to find efficiency number in qmil output
+    """
+    nums = re.findall("\d+\.\d+", line)
+    try:
+        num = float(nums[0])
+    except:
+        # catch error of efficiency not found
+        print("no number found here")
+        return 0
+    if len(nums) > 1 or num > 1 or num < 0:
+        # check if too many values found or efficiency out of range
+        print("something wrong")
+        return 0
+    return num
+
+def get_eff(rpm, out_file=""):
+    """
+    Rewrites template qmil file and replaces desired rpm and returns
+    efficiency of new designed prop
+    If out_file is specified, that prop geometry is saved to the given filename
+    """
+    with open('template.mil', 'r') as file:
+        data = file.readlines()
+    data[18] = "  " + str(rpm) + " ! rpm\n"
+    with open('output.mil', 'w') as file:
+        file.writelines(data)
+
+    filename = 'output.mil'
+
+    cmd = ['qmil', filename, out_file]
+    process = subprocess.Popen(cmd, stdout=subprocess.PIPE)
+    output = process.communicate()[0]
+    for line in output.splitlines():
+        l = str(line)
+        if "eff = " in l:
+            eff = get_num(l)
+    process.wait()
+    return -eff
+
+def plot_prop(propfile):
+    f = open(propfile, "r")
+    contents = f.readlines()
+    prop_geom =  contents[-32:]
+    P = len(prop_geom)
+    r = []
+    c = []
+    beta = []
+    R = 1.3
+    for i in range(1,P):
+        line = prop_geom[i].split()
+        # print(line)
+        r.append(float(line[0])/R)
+        c.append(float(line[1])/R)
+        beta.append(float(line[2]))
+    f = plt.figure()
+    plt.subplot(211)
+    plt.plot(r,c,'b')
+    plt.ylabel("c/R")
+    plt.title("Propeller, R = 1.3m")
+    plt.subplot(212)
+    plt.plot(r, beta, 'b', label = "QMIL")
+    plt.ylabel(r"$\beta$")
+    plt.xlabel('r/R')
+
+    f1 = plt.figure()
+    c = np.array(c)
+    plt.plot(r, c/2)
+    plt.plot(r, -c/2)
+    plt.xlim(0, 1)
+    plt.ylim(-0.5, 0.5)
+    plt.title("Propeller, R = 1.3m")
+    plt.ylabel("x/R")
+    plt.xlabel('y/R')
+    plt.gca().set_aspect('equal', adjustable='box')
+    plt.show()
+
+def opt_rpm():
+    """
+    Optimize for efficiency on rpm given an airspeed
+    """
+    res = minimize_scalar(get_eff, bounds=(200, 1500), method='bounded')
+    if res.success:
+        get_eff(res.x, "best_prop")
+        print(-res.fun, res.x)
+        plot_prop("best_prop")
+        return res.x
+    else:
+        print("Unsuccessful optimization", res.x, res.fun)
+
+opt_rpm()
+# plot_prop("best_prop")
+# rpms = np.arange(200.0, 1500.0, 50.0)
+# etas = np.zeros(len(rpms))
+#
+# for i, rpm in enumerate(rpms):
+#     eff = get_eff(rpm)
+#     etas[i] = eff
+#
+# print(etas)
+# best_rpm_i = np.argmax(etas)
+# print(best_rpm_i, rpms[best_rpm_i], etas[best_rpm_i])
+# get_eff(rpms[best_rpm_i], "best_prop")