Gbopt

.gitignore

@@ -1 +1,6 @@
 /n2.html
+/old.py
+/notes.pdf
+/notes.log
+/notes.aux
+/barntest*

Aaron_Notes.txt

@@ -0,0 +1,29 @@
+Motor Sizing Program
+--------------------
+Using OpenMDAO to optimize motor geometry based on given parameters and constraints.
+
+Main Goal:
+Create more accurate motor sizing code and then cross-check the results with an optimized motor design.
+
+Create different branches in git for different features.
+i.e.
+- Torque Calculation
+- Losses Calculation
+-etc
+
+--------------------
+
+Known Issues:
+- Flux Density calculation is incorrect
+	- Research and fix
+
+- Torque calculation is incorrect
+	- Research and fix
+
+
+Wanted Features:
+- Torque vs Speed curve plot
+
+- Efficiency Table (Using Interpolation)
+	- Create Table / .csv of motor efficiencies (from Motor-CAD or Python code???)
+	- Feedback into motor-sizing program

barntest.py

@@ -0,0 +1,61 @@
+import numpy as np
+from math import pi, sin, cos, exp
+import matplotlib.pyplot as plt
+
+
+def RadialR(n, alpha):
+
+    if n%2 != 0:
+        return alpha*(sin(n*alpha*pi/2))/(n*alpha*pi/2)
+
+    if n%2 == 0:
+        return 0
+
+def RadialT(n, alpha):
+    return 0
+
+
+gap = .001
+r_s = .051
+n_m = 20
+n_s = 24
+ts = 2*pi/n_s #slot to slot angle 
+tt = .6*pi/10 #tooth width angle
+mu_r = 1.05
+t_mag = .003
+l_st = .008
+k_st = 1
+w_tb = .003
+r_m = r_s - gap
+b_r = 1.3
+alpha = 1/30 #magnetic fraction, fraction of rotor circumference occupied by one magnet
+n_p = n_m/2 #number of pole pairs
+r_r = r_m - t_mag
+
+N = 21
+i = np.zeros(2*N+1, dtype = complex)
+for n in range(-N, N+1):
+    beta = n*n_p
+    k_rn = RadialR(n, alpha)
+
+    k_tn = RadialT(n, alpha)
+
+    k_mc = (k_rn + 1j*beta*k_tn)/(1 - beta**2)
+
+    del_r = (mu_r + 1)*((r_r/r_m)**(2*beta) - (r_s/r_m)**(2*beta)) + (mu_r - 1)*(1 - ((r_r/r_m)**(2*beta))*((r_s/r_m)**(2*beta)))
+
+    #print(del_r)
+    if n != 0:
+        k_a = (k_mc*k_rn/del_r)*((1-(r_r/r_m)**(2*beta))*((beta + 1)*((r_r/r_m)**(2*beta)) - (2*beta)*((r_r/r_m)**(beta+1)) + beta - 1))
+    else:
+        k_a = 0
+
+
+    #k_a = k_a.real
+    #print(n)
+    Bgn = -b_r*k_a*(((.0505/r_m)**(beta-1)) + ((r_s/r_m)**(2*beta))*((r_m/.0505)**(beta+1)))
+    print(Bgn)
+    i[n+N] = Bgn
+
+plt.plot(i)
+plt.show()

bmatrix.py

@@ -0,0 +1,117 @@
+'''
+Compute the coefficients of the solution to the air gap and magnet region flux density
+
+Also compute the Fourier coefficients specifically for the radial component of the air gap
+
+Based on Appendix B of 'Brushless Permanent Magnet Motor Design' by Duane P. Hanselman
+'''
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Definitions of magnetic field profiles, which determine k_rn and k_tn. Use as appropriate
+
+# Radial Sinusoidal Profile
+def RadSinR(n, alpha):
+    if n == 1 or n == -1:
+        return .5
+    else:
+        return 0
+
+def RadSinT(n, alpha):
+    return 0
+
+# Radial Profile
+def RadialR(n, alpha):
+
+    if n%2 != 0:
+        return alpha*np.sin(n*alpha*np.pi/2)/(n*alpha*np.pi/2)
+
+    if n%2 == 0:
+        return 0
+
+def RadialT(n, alpha):
+    return 0
+
+
+# Permitivity of Free Space. Cancels out in computation, so the value doesn't matter
+mu_0 = 0.001
+
+
+# Compute the coefficients of the assumed solution using the boundary conditions
+# Solves a size 3 linear system for every nth Fourier series term
+def B_fourier_terms(n, params):
+    n_p = params['n_p'] #number of pole pairs
+    r_s = params['r_s'] #stator inner radius
+    r_m = params['r_m'] #rotor outer radius (INCLUDING MAGNETS)
+    r_r = params['r_r'] #rotor outer radius
+    b_r = params['b_r'] #magnetic remanence of stator material
+    mu_r = params['mu_r'] #magnetic relative recoil permitivity
+    alpha = params['alpha'] #magnet fraction in rotor
+
+    bmat = np.zeros((3, 3))
+    rhs = np.zeros(3)
+
+    beta = n*n_p
+
+    # magnetic profile coefficients
+    k_rn = RadialR(n, alpha)
+    k_tn = RadialT(n, alpha)
+
+    # assume output is real, ignore imaginary part (Equation B.18)
+    Cm = b_r*(k_rn)/((1-beta**2)*mu_r*mu_0)
+
+    #Hatheta boundary condition (excluded from linear system)
+    Ea = -(r_s)**(2*beta)
+    print(r_s)
+    #Hmtheta boundary condition
+    bmat[0, 1] = (r_r**(-beta-1))
+    bmat[0, 2] = (r_r**(beta-1))
+    rhs[0] = -Cm
+
+    #Htheta boundary condition
+    bmat[1, 0] = -Ea*(r_m**(-beta-1)) - r_m**(beta - 1)
+    bmat[1, 1] = r_m**(-beta-1)
+    bmat[1, 2] = r_m**(beta-1)
+    rhs[1] = -Cm
+
+    #Hr boundary condition
+    bmat[2, 0] = (beta/mu_r)*(Ea*(r_m**(-beta-1)) - r_m**(beta - 1))
+    bmat[2, 1] = -beta*(r_m**(-beta-1))
+    bmat[2, 2] = beta*(r_m**(beta-1))
+    rhs[2] = (b_r*k_rn/(mu_r*mu_0)) - Cm
+
+    print(bmat)
+    #invert matrix
+    bmat_inv = np.linalg.inv(bmat)
+
+    #solve system using inverse of system matrix
+    x = bmat_inv.dot(rhs)
+    print(x)
+    return x
+
+# Compute Fourier Coefficients of radial component air gap magnetic flux Density
+def B_arn(b_fourier, n, params, r):
+    r_s = params['r_s'] #stator inner radius
+    n_p = params['n_p'] #number of pole pairs
+
+    Da = -b_fourier[0]
+    beta = n*n_p
+    Ea = -(r_s)**(2*beta)
+    b_arn = mu_0*beta*Da*(Ea*(r**(-beta-1)) - r**(beta-1)) #Equation B.13
+    return b_arn
+
+if __name__ == '__main__':
+
+    N = 17
+    B_r = np.zeros(N)
+
+    for n in range(1, N):
+        b_f = B_fourier_terms(n)
+        B_r[n] = B_arn(b_f, n, r_s)
+        print(B_r[n])
+
+    i = np.fft.irfft(B_r)
+    #print(i)
+    plt.plot(i)
+    plt.show()

btooth.py

@@ -0,0 +1,102 @@
+'''
+This function computes magnetic flux density in the stator teeth of a motor, described using a Fourier series with respect to angular position theta
+
+Taken from Brushless Permanent Magnet Motor Design, Second Edition by Dr. Duane Hanselman, Chapter 7
+Magnetic Flux Density Coefficient formula taken from Appendix B, pg 349
+
+'''
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import bmatrix as bm
+import slotcorrection as sl
+import core_loss as cl
+
+def B_tooth(N, M, params):
+    n_s = params['n_s'] #number of Slots
+    n_p = params['n_p'] #number of pole pairs
+    k_st = params['k_st'] #lamination stacking factor
+    w_tb = params['w_tb'] #tooth body width
+    r_s = params['r_s'] #stator inner radius
+    t_mag = params['t_mag'] #magnet thickness
+    r_m = params['r_m'] #rotor outer radius (INCLUDING MAGNETS)
+    r_r = params['r_r'] #rotor outer radius
+    b_r = params['b_r'] #magnetic remanence
+    tt = params['tt']
+    ts = params['ts']
+    mu_r = params['mu_r']
+    alpha = params['alpha']
+
+    n_m = n_p*2
+    B_gn = np.zeros(N)
+    B_tn = np.zeros(N)
+
+    # Compute slot correction factor terms
+    K_slm = sl.slot_correction(M, params)
+
+    for n in range(1, N):
+        # Compute solution coefficients for a particular n
+        b_f = bm.B_fourier_terms(n, params)
+        # Compute radial air gap magnetic flux density nth Fourier Coefficient
+        B_gn[n] = bm.B_arn(b_f, n, params, r_s)
+
+        slotsum = 0
+        # Compute characteristic summation of slot correction factor for the nth Fourier Coefficient
+        if n%2 != 0:
+            for m in range(-M, M + 1):
+                slotsum += K_slm[m + M]*np.sinc(np.pi*(m + n*(n_m/(2*n_s))))#/(np.pi*(m + n*(n_m/(2*n_s))))
+        else:
+            slotsum = 0
+
+        # Modify air gap solution with slot correction
+        B_tn[n] = B_gn[n]*(2*np.pi*r_s/(n_s*k_st*w_tb))*slotsum
+    return B_tn
+
+if __name__ == '__main__':
+
+    gap = .001
+    r_s = .051
+    n_p = 8
+    n_s = 24
+    ts = 2*np.pi/n_s #slot to slot angle
+    tt = .6*np.pi/10 #tooth width angle
+    mu_r = 1.05
+    t_mag = .003
+    l_st = .008
+    k_st = 1
+    w_tb = .005
+    alpha = .85
+    b_r = 1.3
+    f = 5400/20 #rpm
+
+    params = {
+        'gap':gap,
+        'r_s':r_s,
+        'n_p':n_p,
+        'n_s':n_s,
+        'tt':tt,
+        'ts':ts,
+        'mu_r':mu_r,
+        't_mag':t_mag,
+        'sta_ir':r_s,
+        'l_st':l_st,
+        'k_st':k_st,
+        'w_tb':w_tb,
+        'r_m':(r_s - gap),
+        'r_r':(r_s - gap - t_mag),
+        'alpha':alpha,
+        'b_r':b_r,
+        'f':f
+    }
+
+    B_tn = B_tooth(4, 20, params)
+    print(B_tn)
+    i = np.fft.irfft(B_tn)
+
+    L_core = cl.core_loss(i, params)
+
+    # print(i)
+    print(L_core)
+    # plt.plot(i)
+    # plt.show()

core_loss.py

@@ -0,0 +1,36 @@
+'''
+Compute core loss for either stator teeth or stator yoke
+
+Depends on 4 constants characteristic of the lamination material, in this case Hiperco 50
+Will need to determine them, either through data sheets or experimentation
+
+From section 9.4 of Hanselman book
+'''
+
+import numpy as np
+
+def core_loss(B_n, params):
+    f = params['f'] # motor speed in rpm
+
+    # convert motor speed to rad/s
+    omega = f*(2*np.pi/60)
+
+    # peak magnetic flux density is the largest data point from an inverse fourier transform
+    B_pk = max(B_n)
+
+    N = len(B_n)
+    #assume B_n is real, compute mean-square value of dB/dt
+    B_msq = 0
+    for t in range(1, N):
+        B_msq += 2*((t*omega*B_n[t])**2)
+
+    # k_h, k_e, n, m are material constants and need to be fitted to real material data
+    n = 2
+    m = 0.01
+    k_h = 0.5
+    k_e = 0.5
+
+    # compute core loss
+    L_core = k_h*f*(B_pk**(n + m*B_pk)) + B_msq*k_e/(2*(np.pi**2))
+
+    return L_core