diff --git a/.gitignore b/.gitignore index 6094a05..d878c80 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,6 @@ /n2.html +/old.py +/notes.pdf +/notes.log +/notes.aux +/barntest* diff --git a/Aaron_Notes.txt b/Aaron_Notes.txt new file mode 100644 index 0000000..a62ec3a --- /dev/null +++ b/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 diff --git a/barntest.py b/barntest.py new file mode 100644 index 0000000..7d1d914 --- /dev/null +++ b/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() diff --git a/bmatrix.py b/bmatrix.py new file mode 100644 index 0000000..fc24c08 --- /dev/null +++ b/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() diff --git a/btooth.py b/btooth.py new file mode 100644 index 0000000..40128c4 --- /dev/null +++ b/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() diff --git a/core_loss.py b/core_loss.py new file mode 100644 index 0000000..fc1e2ed --- /dev/null +++ b/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 diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index dee46ef..99eee3d 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -2,13 +2,191 @@ import numpy as np from math import pi from openmdao.api import Problem, IndepVarComp, ExplicitComponent, ExecComp -from openmdao.api import NewtonSolver, Group, DirectSolver, NonlinearRunOnce, LinearRunOnce, view_model, BalanceComp +from openmdao.api import NewtonSolver, Group, DirectSolver, NonlinearRunOnce, LinearRunOnce, view_model, BalanceComp, ScipyOptimizeDriver -class motor_size(ExplicitComponent): +class core_loss(ExplicitComponent): + + def setup(self): + self.add_input('rm', units='rpm', desc='motor speed') + self.add_input('alpha', 1.27, desc='core loss constant') + self.add_input('n_m', desc='number of poles') + + self.add_output('L_core', units='W', desc='core loss') + + self.declare_partials('*','*',method='fd') + def compute(self, inputs, outputs): + rm = inputs['rm'] + n_m = inputs['n_m'] + alpha = inputs['alpha'] + + f = n_m*rm/120 + outputs['L_core'] = .2157*(f**alpha) + +class mag_eddy_loss(ExplicitComponent): + + def setup(self): + self.add_input('rm', units='rpm', desc='motor speed') + self.add_input('n_m', desc='number of poles') + + self.add_output('L_emag', units='W', desc='magnetic eddy loss') + + self.declare_partials('*','*',method='fd') + def compute(self, inputs, outputs): + rm = inputs['rm'] + n_m = inputs['n_m'] + + f = n_m*rm/120 + outputs['L_emag'] = .0010276*(f**2) + +class wire_eddy_loss(ExplicitComponent): + + def setup(self): + self.add_input('rm', units='rpm', desc='motor speed') + self.add_input('n_m', desc='number of poles') + + self.add_output('L_ewir', units='W', desc='winding eddy loss') + + self.declare_partials('*','*',method='fd') + def compute(self, inputs, outputs): + rm = inputs['rm'] + n_m = inputs['n_m'] + + f = n_m*rm/120 + outputs['L_ewir'] = .00040681*(f**2) + +class mech_loss(ExplicitComponent): + + def setup(self): + #windage parameters + self.add_input('rm', 5400, units='rpm', desc='motor speed') + self.add_input('k', .003, desc='windage constant k') + self.add_input('rho_a', 1.225, units='kg/m**3', desc='air density') + self.add_input('rot_or', .05, units='m', desc='rotor outer radius') + self.add_input('rot_ir', .02, units='m', desc='rotor inner radius') + self.add_input('l_st', units='m', desc='length of stack') + self.add_input('gap', units='m', desc='air gap') #aka delta + self.add_input('mu_a', 1.81e-5, units='(m**2)/s', desc='air dynamic viscosity') + + #bearing parameters + self.add_input('muf_b', .3, desc='bearing friction coefficient') + self.add_input('D_b', .01, units='m', desc='bearing bore diameter') + self.add_input('F_b', 100, units='N', desc='bearing load') #coupled with motor mass + + self.add_output('L_airg', units='W', desc='air gap windage loss') #air gap + self.add_output('L_airf', units='W', desc='axial face windage loss') #axial face + self.add_output('L_bear', units='W', desc='bearing loss') + + self.declare_partials('*','*',method='fd') + + def compute(self, inputs, outputs): + + rm = inputs['rm'] + k = inputs['k'] + rho_a = inputs['rho_a'] + rot_or = inputs['rot_or'] + rot_ir = inputs['rot_ir'] + l_st = inputs['l_st'] + gap = inputs['gap'] + mu_a = inputs['mu_a'] + + muf_b = inputs['muf_b'] + D_b = inputs['D_b'] + F_b = inputs['F_b'] + + omega = rm*(2*pi/60) + #air gap friction coefficient, Re number + Reg = rho_a*omega*rot_or*gap/mu_a + Cfg = .515*((gap**.3)/rot_or)/(Reg**.5) + + outputs['L_airg'] = k*Cfg*pi*rho_a*(omega**3)*(rot_or**4)*l_st + + #axial air friction coefficient, Re number + Rea = rho_a*omega*(rot_or**2)/mu_a + Cfa = .146/(Rea**2) + + outputs['L_airf'] = .5*Cfa*rho_a*(omega**3)*((rot_or**5)-(rot_ir**5)) + + outputs['L_bear'] = .5*muf_b*D_b*F_b*omega + + #def compute_partials(self, inputs, J): + +class efficiency(ExplicitComponent): + def setup(self): + # indeps = self.add_subsystem('indeps', IndepVarComp(), promotes_inputs('*')) + # + # indeps.add_output('i', 30, units='A', desc='RMS current') + # indeps.add_output('tq', 25, units='N*m', desc='torque') + # indeps.add_output('v', 385, units='V', desc='RMS voltage') + # indeps.add_output('rm', 5400, units='rpm', desc='motor speed') + + + + self.add_input('i', 30, units='A', desc='RMS current') + self.add_input('tq', 25, units='N*m', desc='torque') + self.add_input('v', 385, units='V', desc='RMS voltage') + self.add_input('rm', 5400, units='rpm', desc='motor speed') + + self.add_input('L_airg', units='W', desc='gap windage loss') + self.add_input('L_airf', units='W', desc='face windage loss') + self.add_input('L_bear', units='W', desc='bearing loss') + self.add_input('L_emag', units='W', desc='magnetic eddy loss') + self.add_input('L_ewir', units='W', desc='wiring eddy loss') + self.add_input('L_core', units='W', desc='core loss') + self.add_input('L_res', 50, units='W', desc='resistive loss') + + self.add_output('P_in', units='W', desc='input power') + #self.add_output('P_out', units='W', desc='output power') + self.add_output('L_total', units='W', desc='total loss') + + self.declare_partials('*','*',method='fd') + + def compute(self, inputs, outputs): + i = inputs['i'] + tq = inputs['tq'] + v = inputs['v'] + rm = inputs['rm'] + L_airg = inputs['L_airg'] + L_airf = inputs['L_airf'] + L_core = inputs['L_core'] + L_bear = inputs['L_bear'] + L_emag = inputs['L_emag'] + L_ewir = inputs['L_ewir'] + L_res = inputs['L_res'] + + + outputs['P_in'] = (3**.5)*v*i + #outputs['P_out'] = tq*rm*(2*pi/60) + #outputs['nu'] = outputs['P_out']/outputs['P_in'] + outputs['L_total'] = L_airg + L_airf + L_bear + L_emag + L_ewir + L_core + L_res + # + # def compute_partials(self, inputs, J): + # i = inputs['i'] + # tq = inputs['tq'] + # v = inputs['v'] + # rm = inputs['rm'] + # + # P_in = (3**.5)*v*i + # P_out = tq*rm*(2*pi/60) + # + # J['P_in','v'] = (3**.5)*i + # J['P_in','i'] = (3**.5)*v + # + # J['P_out','tq'] = rm*(2*pi/60) + # J['P_out','rm'] = tq*(2*pi/60) + # + # J['nu','v'] = -P_out/((3**.5)*i*(v**2)) + # J['nu','i'] = -P_out/((3**.5)*v*(i**2)) + # J['nu','tq'] = rm*(2*pi/60)/P_in + # J['nu','rm'] = tq*(2*pi/60)/P_in + + + +class motor_size(ExplicitComponent): + def setup(self): # rotor_outer_radius - self.add_input('r_m', 0.0765, units='m', desc='outer radius of motor') + self.add_input('mot_or', 0.0765, units='m', desc='motor outer radius') self.add_input('gap', 0.001, units='m', desc='air gap') # rotor_yoke_width @@ -41,14 +219,14 @@ def setup(self): self.add_input('n', 16, desc='number of wire turns') self.add_input('i', 30, units='A', desc='RMS current') self.add_input('k_wb', 0.65, desc='bare wire slot fill factor') - self.add_output('J', units='A/mm**2', desc='Current density') + self.add_output('j', units='A/mm**2', desc='Current density') - self.declare_partials('*','*', method='fd') + self.declare_partials('*','*')#, method='fd') def compute(self,inputs,outputs): # rotor_outer_radius rot_or = inputs['rot_or'] - r_m = inputs['r_m'] # .0765 + mot_or = inputs['mot_or'] # .0765 gap = inputs['gap'] # rotor_yoke_width b_g = inputs['b_g'] @@ -69,55 +247,104 @@ def compute(self,inputs,outputs): b_t = inputs['b_t'] outputs['w_t'] = (2*pi*rot_or*b_g) / (n_s*k*b_t) # Exec Comps - # print(r_m ,rot_or , gap , outputs['w_sy']) - # print(r_m - rot_or - gap - outputs['w_sy']) - outputs['s_d'] = r_m - rot_or - gap - outputs['w_sy'] - #outputs['r_m'] = rot_or + gap + s_d + outputs['w_sy'] - outputs['rot_ir'] = (rot_or- t_mag) - outputs['w_ry'] + # print(mot_or ,rot_or , gap , outputs['w_sy']) + # print(mot_or - rot_or - gap - outputs['w_sy']) + outputs['s_d'] = mot_or - rot_or - gap - outputs['w_sy'] + #outputs['mot_or'] = rot_or + gap + s_d + outputs['w_sy'] + outputs['rot_ir'] = (rot_or - t_mag) - outputs['w_ry'] outputs['sta_ir'] = rot_or + gap - area = pi*(r_m-outputs['w_sy'])**2 - pi*(r_m-outputs['w_sy']-outputs['s_d'])**2 #outputs['sta_ir'] - outputs['J'] = 2*n*i*(2.**0.5)/(k_wb/n_s*(area-n_s*1.25*(outputs['w_t']*outputs['s_d']))*1E6) - - - # def compute_partials(self, inputs, J): - - # # rotor_outer_radius - # # r_m = inputs['r_m'] - # # s_d = inputs['s_d'] - # # w_sy = inputs['w_sy'] - # # J['rot_or', 'r_m'] = 1-w_sy-s_d-gap - # # J['rot_or', 's_d'] = r_m - w_sy - 1 - gap - # # J['rot_or', 'w_sy'] = r_m - 1 - s_d - gap - - # # rotor_yoke_width - # rot_or = inputs['rot_or'] - # b_g= inputs['b_g'] - # n_m= inputs['n_m'] - # k = inputs['k'] - # b_ry= inputs['b_ry'] - # J['w_ry', 'rot_or'] = (pi*b_g)/(n_m*k*b_ry) - # J['w_ry', 'b_g'] = (pi*rot_or)/(n_m*k*b_ry) - # J['w_ry', 'n_m'] = -(pi*rot_or*b_g)/(n_m**3*k*b_ry) - # J['w_ry', 'k'] = -(pi*rot_or*b_g)/(n_m*k**2*b_ry) - # J['w_ry', 'b_ry'] = -(pi*rot_or*b_g)/(n_m*k*b_ry**2) - - # # stator_yoke_width - # b_sy= inputs['b_sy'] - # J['w_sy', 'rot_or'] = (pi*b_g)/(n_m*k*b_sy) - # J['w_sy', 'b_g'] = (pi*rot_or)/(n_m*k*b_sy) - # J['w_sy', 'n_m'] = -(pi*rot_or*b_g)/(n_m**3*k*b_sy) - # J['w_sy', 'k'] = -(pi*rot_or*b_g)/(n_m*k**2*b_sy) - # J['w_sy', 'b_sy'] = -(pi*rot_or*b_g)/(n_m*k*b_sy**2) - - # # tooth_width - # n_s = inputs['n_s'] - # b_t = inputs['b_t'] - # J['w_t', 'rot_or'] = (2*pi*b_g)/(n_s*k*b_t) - # J['w_t', 'b_g'] = (2*pi*rot_or)/(n_s*k*b_t) - # J['w_t', 'n_s'] = -(2*pi*rot_or*b_g)/(n_s**2*k*b_t) - # J['w_t', 'k'] = -(2*pi*rot_or*b_g)/(n_s*k**2*b_t) - # J['w_t', 'b_t'] = -(2*pi*rot_or*b_g)/(n_s*k*b_t**2) + area = pi*(mot_or-outputs['w_sy'])**2 - pi*(mot_or-outputs['w_sy']-outputs['s_d'])**2 #outputs['sta_ir'] + outputs['j'] = 2*n*i*(2.**0.5)/(k_wb/n_s*(area-n_s*1.25*(outputs['w_t']*outputs['s_d']))*1E6) + print(outputs['j']) + + def compute_partials(self, inputs, J): + + # rotor_yoke_width + rot_or = inputs['rot_or'] + b_g= inputs['b_g'] + n_m= inputs['n_m'] + k = inputs['k'] + b_ry= inputs['b_ry'] + J['w_ry', 'rot_or'] = (pi*b_g)/(n_m*k*b_ry) + J['w_ry', 'b_g'] = (pi*rot_or)/(n_m*k*b_ry) + J['w_ry', 'n_m'] = -(pi*rot_or*b_g)/(n_m**2*k*b_ry) + J['w_ry', 'k'] = -(pi*rot_or*b_g)/(n_m*k**2*b_ry) + J['w_ry', 'b_ry'] = -(pi*rot_or*b_g)/(n_m*k*b_ry**2) + + # stator_yoke_width + b_sy= inputs['b_sy'] + J['w_sy', 'rot_or'] = (pi*b_g)/(n_m*k*b_sy) + J['w_sy', 'b_g'] = (pi*rot_or)/(n_m*k*b_sy) + J['w_sy', 'n_m'] = -(pi*rot_or*b_g)/(n_m**2*k*b_sy) + J['w_sy', 'k'] = -(pi*rot_or*b_g)/(n_m*k**2*b_sy) + J['w_sy', 'b_sy'] = -(pi*rot_or*b_g)/(n_m*k*b_sy**2) + + # tooth_width + n_s = inputs['n_s'] + b_t = inputs['b_t'] + J['w_t', 'rot_or'] = (2*pi*b_g)/(n_s*k*b_t) + J['w_t', 'b_g'] = (2*pi*rot_or)/(n_s*k*b_t) + J['w_t', 'n_s'] = -(2*pi*rot_or*b_g)/(n_s**2*k*b_t) + J['w_t', 'k'] = -(2*pi*rot_or*b_g)/(n_s*k**2*b_t) + J['w_t', 'b_t'] = -(2*pi*rot_or*b_g)/(n_s*k*b_t**2) + + # slot_depth + mot_or = inputs['mot_or'] + gap = inputs['gap'] + J['s_d', 'mot_or'] = 1 + J['s_d', 'rot_or'] = -1 - J['w_sy', 'rot_or'] + J['s_d', 'gap'] = -1 + J['s_d', 'b_g'] = -J['w_sy', 'b_g'] + J['s_d', 'n_m'] = -J['w_sy', 'n_m'] + J['s_d', 'k'] = -J['w_sy', 'k'] + J['s_d', 'b_sy'] = -J['w_sy', 'b_sy'] + + # rotor_inner_radius + t_mag = inputs['t_mag'] + J['rot_ir', 'rot_or'] = 1 - J['w_ry', 'rot_or'] + J['rot_ir', 't_mag'] = -1 + J['rot_ir', 'b_g'] = - J['w_ry', 'b_g'] + J['rot_ir', 'n_m'] = - J['w_ry', 'n_m'] + J['rot_ir', 'k'] = - J['w_ry', 'k'] + J['rot_ir', 'b_ry'] = - J['w_ry', 'b_ry'] + + # stator_inner_radius + J['sta_ir', 'rot_or'] = 1 + J['sta_ir', 'gap'] = 1 + + # current_density + n = inputs['n'] + i = inputs['i'] + k_wb = inputs['k_wb'] + n_s = inputs['n_s'] + w_sy = (pi*rot_or*b_g)/(n_m*k*b_sy) + s_d = mot_or - rot_or - gap - w_sy + w_t = (2*pi*rot_or*b_g) / (n_s*k*b_t) + area = pi*(mot_or-w_sy)**2 - pi*(mot_or-w_sy-s_d)**2 + + djdarea = -2*n*i*(2.**0.5)/(k_wb/n_s*((area-n_s*1.25*(w_t*s_d))**2)*1E6) + djds_d = djdarea*(-n_s*1.25*w_t) + djdw_t = djdarea*(-n_s*1.25*s_d) + djdn_s = 2*n*i*(2.**0.5)*area/(k_wb*((area-n_s*1.25*(w_t*s_d))**2)*1E6) + + dads_d = 2*pi*(mot_or-w_sy-s_d) + dadmot_or = 2*pi*s_d + dadw_sy = -2*pi*s_d + + + J['j', 'n'] = 2*i*(2.**0.5)/(k_wb/n_s*(area-n_s*1.25*(w_t*s_d))*1E6) + J['j', 'i'] = 2*n*(2.**0.5)/(k_wb/n_s*(area-n_s*1.25*(w_t*s_d))*1E6) + J['j', 'k_wb'] = - 2*n*i*(2.**0.5)/((k_wb**2)/n_s*(area-n_s*1.25*(w_t*s_d))*1E6) + J['j', 'mot_or'] = djdarea*(dads_d*J['s_d', 'mot_or'] + dadmot_or) + djds_d*J['s_d', 'mot_or'] + J['j', 'rot_or'] = djdarea*(dadw_sy*J['w_sy', 'rot_or'] + dads_d*J['s_d','rot_or']) + djdw_t*J['w_t','rot_or'] + djds_d*J['s_d','rot_or'] + J['j', 'gap'] = (djdarea*dads_d + djds_d)*J['s_d','gap'] + J['j', 'b_g'] = djdarea*(dadw_sy*J['w_sy','b_g']+dads_d*J['s_d','b_g']) + djdw_t*J['w_t','b_g'] + djds_d*J['s_d','b_g'] + J['j', 'n_m'] = djdarea*(dadw_sy*J['w_sy','n_m']+dads_d*J['s_d','n_m']) + djds_d*J['s_d','n_m'] + J['j', 'n_s'] = djdn_s + djdw_t*J['w_t','n_s'] + J['j', 'k'] = djdarea*(dadw_sy*J['w_sy','k'] + dads_d*J['s_d','k']) + djdw_t*J['w_t','k'] + djds_d*J['s_d','k'] + J['j', 'b_sy'] = djdarea*(dadw_sy*J['w_sy','b_sy'] + dads_d*J['s_d','b_sy']) + djds_d*J['s_d','b_sy'] + J['j', 'b_t'] = djdw_t*J['w_t','b_t'] class torque(ExplicitComponent): @@ -130,7 +357,7 @@ def setup(self): self.add_input('rot_or', .025, units='m', desc='rotor outer radius') self.add_output('tq', 25, units='N*m', desc='torque') #self.declare_partials('tq', ['n_m','n','b_g','l_st','rot_or','i']) - self.declare_partials('*','*', method='fd') + self.declare_partials('*','*')#, method='fd') def compute(self,inputs,outputs): n_m=inputs['n_m'] @@ -142,34 +369,34 @@ def compute(self,inputs,outputs): outputs['tq'] = 2*n_m*n*b_g*l_st*rot_or*i*.68 - # def compute_partials(self,inputs,J): - # n_m=inputs['n_m'] - # n= inputs['n'] - # b_g= inputs['b_g'] - # l_st= inputs['l_st'] - # rot_or = inputs['rot_or'] - # i = inputs['i'] - - # J['tq','n_m'] = 2*n*b_g*l_st*rot_or*i - # J['tq', 'n'] = 2*n_m*b_g*l_st*rot_or*i - # J['tq', 'b_g'] = 2*n_m*n*l_st*rot_or*i - # J['tq', 'l_st'] = 2*n_m*n*b_g*rot_or*i - # J['tq', 'rot_or'] = 2*n_m*n*b_g*l_st*i - # J['tq', 'i'] = 2*n_m*n*b_g*l_st*rot_or + def compute_partials(self,inputs,J): + n_m=inputs['n_m'] + n= inputs['n'] + b_g= inputs['b_g'] + l_st= inputs['l_st'] + rot_or = inputs['rot_or'] + i = inputs['i'] + + J['tq', 'n_m'] = 2*n*b_g*l_st*rot_or*i*.68 + J['tq', 'n'] = 2*n_m*b_g*l_st*rot_or*i*.68 + J['tq', 'b_g'] = 2*n_m*n*l_st*rot_or*i*.68 + J['tq', 'l_st'] = 2*n_m*n*b_g*rot_or*i*.68 + J['tq', 'rot_or'] = 2*n_m*n*b_g*l_st*i*.68 + J['tq', 'i'] = 2*n_m*n*b_g*l_st*rot_or*.68 class motor_mass(ExplicitComponent): def setup(self): # stator self.add_input('rho', 8110.2, units='kg/m**3', desc='density of hiperco-50') - self.add_input('r_m', .075, units='m', desc='motor outer radius') + self.add_input('mot_or', .075, units='m', desc='motor outer radius') self.add_input('n_s', 15, desc='number of slots') self.add_input('sta_ir', .050, units='m', desc='stator inner radius') self.add_input('w_t', units='m', desc='tooth width') self.add_input('l_st', units='m', desc='length of stack') self.add_input('s_d', units='m', desc='slot depth') self.add_output('sta_mass', 25, units='kg', desc='mass of stator') - #self.declare_partials('sta_mass', ['rho','r_m','n_s','sta_ir','w_t','l_st','s_d'], method='fd') + #self.declare_partials('sta_mass', ['rho','mot_or','n_s','sta_ir','w_t','l_st','s_d'], method='fd') # rotor self.add_input('rot_or', 0.0615, units='m', desc='rotor outer radius') self.add_input('rot_ir', 0.0515, units='m', desc='rotor inner radius') @@ -180,19 +407,19 @@ def setup(self): self.add_output('mag_mass', 0.5, units='kg', desc='mass of magnets') #self.declare_partials('rot_mass',['rho','rot_or','rot_ir'], method='fd') - self.declare_partials('*','*', method='fd') + self.declare_partials('*','*')#, method='fd') def compute(self,inputs,outputs): # stator rho=inputs['rho'] - r_m=inputs['r_m'] + mot_or=inputs['mot_or'] n_s=inputs['n_s'] sta_ir=inputs['sta_ir'] w_t=inputs['w_t'] l_st=inputs['l_st'] s_d=inputs['s_d'] - outputs['sta_mass'] = rho * l_st * ((pi * r_m**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*s_d*1.5))) + outputs['sta_mass'] = rho * l_st * ((pi * mot_or**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*s_d*1.5))) # rotor rot_ir=inputs['rot_ir'] @@ -209,22 +436,46 @@ def compute(self,inputs,outputs): # - # def compute_partials(self,inputs,J): - - # stator - # rho=inputs['rho'] - # r_m=inputs['r_m'] - # n_s=inputs['n_s'] - # sta_ir=inputs['sta_ir'] - # w_t=inputs['w_t'] - # l_st=inputs['l_st'] - - # J['sta_mass', 'rho'] = - # J['sta_mass', 'r_m'] = - # J['sta_mass', 'n_s'] = - # J['sta_mass', 'sta_ir'] = - # J['sta_mass', 'w_t'] = - # J['sta_mass', 'l_st'] = + def compute_partials(self,inputs,J): + + #stator + rho=inputs['rho'] + mot_or=inputs['mot_or'] + n_s=inputs['n_s'] + sta_ir=inputs['sta_ir'] + w_t=inputs['w_t'] + l_st=inputs['l_st'] + s_d=inputs['s_d'] + + J['sta_mass', 'rho'] = l_st * ((pi * mot_or**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*s_d*1.5))) + J['sta_mass', 'mot_or'] = 2 * rho * l_st * (pi * mot_or) + J['sta_mass', 'n_s'] = rho * l_st * (w_t*s_d*1.5) + J['sta_mass', 'sta_ir'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) + J['sta_mass', 'w_t'] = rho * l_st * (n_s*(s_d*1.5)) + J['sta_mass', 'l_st'] = rho * ((pi * mot_or**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*s_d*1.5))) + J['sta_mass', 's_d'] = rho * l_st * (-(2 * pi * (sta_ir+s_d)) + (n_s*w_t*1.5)) + + #rotor + rot_ir=inputs['rot_ir'] + rot_or=inputs['rot_or'] + l_st=inputs['l_st'] + t_mag=inputs['t_mag'] + + J['rot_mass', 'rot_ir'] = -(2 * pi * rot_ir) * rho * l_st + J['rot_mass', 'rot_or'] = (2 * pi * (rot_or - t_mag)) * rho * l_st + J['rot_mass', 'l_st'] = (pi*(rot_or - t_mag)**2 - pi*rot_ir**2) * rho + J['rot_mass', 't_mag'] = (2 * pi * (rot_or - t_mag)) * rho * l_st + J['rot_mass', 'rho'] = (pi*(rot_or - t_mag)**2 - pi*rot_ir**2) * l_st + + #magnets + rho_mag=inputs['rho_mag'] + + J['mag_mass', 'rot_or'] = ((2*pi*rot_or) - (2*pi*(rot_or-t_mag))) * rho_mag * l_st + J['mag_mass', 't_mag'] = - (2*pi*(rot_or-t_mag)) * rho_mag * l_st + J['mag_mass', 'rho_mag'] = (((pi*rot_or**2) - (pi*(rot_or-t_mag)**2))) * l_st + J['mag_mass', 'l_st'] = (((pi*rot_or**2) - (pi*(rot_or-t_mag)**2))) * rho_mag + + if __name__ == "__main__": p = Problem() @@ -235,10 +486,10 @@ def compute(self,inputs,outputs): #ind.add_output('rot_or', val=0.0615, units='m') # Outer radius of rotor, including 5mm magnet thickness ind.add_output('k', val=0.97) # Stacking factor ind.add_output('k_wb', val=0.55) # copper fill factor - ind.add_output('gap', val=0.001, units='m') # Stacking factor + ind.add_output('gap', val=0.001, units='m') # air gap ind.add_output('n', val=24) # Number of wire turns ind.add_output('i', val=33, units='A') # RMS Current - ind.add_output('r_m', val=0.0795, units='m') # Motor outer radius + ind.add_output('mot_or', val=0.08, units='m') # Motor outer radius ind.add_output('t_mag', val=.005, units='m') # Magnet thickness ind.add_output('b_g', val = 1, units='T') # Air gap flux Density !! Flux values may represent 100% slot fill !! @@ -253,45 +504,80 @@ def compute(self,inputs,outputs): ind.add_output('rho', val=8110.2, units='kg/m**3') # Density of Hiperco-50 ind.add_output('rho_mag', val=7500, units='kg/m**3') # Density of Magnets - bal = BalanceComp() + ind.add_output('rot_or', val = 0.06, units='m') + ind.add_output('l_st', val = 0.02, units='m') + + ind.add_output('rm', val = 5400, units='rpm') - bal.add_balance('rot_or', val=0.06, units='m', use_mult=False, rhs_val = 13.) - bal.add_balance('l_st', val=0.02, units='m', use_mult=False, rhs_val = 24.) + # bal = BalanceComp() + # + # bal.add_balance('rot_or', val=0.06, units='m', use_mult=False, rhs_val = 13.) + # bal.add_balance('l_st', val=0.02, units='m', use_mult=False, rhs_val = 24.) - model.add_subsystem('size', motor_size(), promotes_inputs=['n','i','k_wb','r_m','gap','rot_or','b_g','k','b_ry','n_m','b_sy','b_t','n_s'], promotes_outputs=['w_ry', 'w_sy', 'w_t','s_d','rot_ir','sta_ir']) + model.add_subsystem('size', motor_size(), promotes_inputs=['n','i','k_wb','mot_or','gap','rot_or','b_g','k','b_ry','n_m','b_sy','b_t','n_s'], promotes_outputs=['w_ry', 'w_sy', 'w_t','s_d','rot_ir','sta_ir']) # model.add_subsystem('motor_radius_prime', ExecComp('r_m_p = rot_or + .005 + .001 + s_d + w_sy',r_m_p={'units':'m'}, rot_or={'units':'m'}, s_d={'units':'m'}, w_sy={'units':'m'}), promotes_inputs=['rot_or','s_d','w_sy'], promotes_outputs=['r_m_p']) - # model.add_subsystem('mass_stator', mass_stator(), promotes_inputs=['rho','r_m','n_s','sta_ir','w_t','l_st'], promotes_outputs=['weight'] - # model.add_subsystem('stmass', ExecComp('mass = l_st * ((pi * r_m**2)-(pi * sta_ir**2)+(n_s*(w_t*1.2)))', l_st={'units':'m'},r_m={'units':'m'},sta_ir={'units':'m'},w_t={'units':'m'}), promotes_inputs=['l_st','r_m','sta_ir','n_s','w_t'], promotes_outputs=['mass'] - model.add_subsystem(name='balance', subsys=bal) - model.add_subsystem('mass', motor_mass(), promotes_inputs=['t_mag','rho_mag','rho','r_m','n_s','sta_ir','w_t','l_st','s_d','rot_or','rot_ir'], promotes_outputs=['sta_mass','rot_mass','mag_mass']) + # model.add_subsystem('mass_stator', mass_stator(), promotes_inputs=['rho','mot_or','n_s','sta_ir','w_t','l_st'], promotes_outputs=['weight'] + # model.add_subsystem('stmass', ExecComp('mass = l_st * ((pi * mot_or**2)-(pi * sta_ir**2)+(n_s*(w_t*1.2)))', l_st={'units':'m'},mot_or={'units':'m'},sta_ir={'units':'m'},w_t={'units':'m'}), promotes_inputs=['l_st','mot_or','sta_ir','n_s','w_t'], promotes_outputs=['mass'] + # model.add_subsystem(name='balance', subsys=bal) + model.add_subsystem('mass', motor_mass(), promotes_inputs=['t_mag','rho_mag','rho','mot_or','n_s','sta_ir','w_t','l_st','s_d','rot_or','rot_ir'], promotes_outputs=['sta_mass','rot_mass','mag_mass']) model.add_subsystem('torque', torque(), promotes_inputs=['rot_or','b_g','i','n_m','n','l_st'], promotes_outputs=['tq']) + + model.add_subsystem('mech_loss', mech_loss(), promotes_inputs=['rm','rot_or','rot_ir','l_st','gap'], promotes_outputs=['L_airg','L_airf','L_bear']) + model.add_subsystem('mag_eddy_loss', mag_eddy_loss(), promotes_inputs=['rm','n_m'], promotes_outputs=['L_emag']) + model.add_subsystem('wire_eddy_loss', wire_eddy_loss(), promotes_inputs=['rm','n_m'], promotes_outputs=['L_ewir']) + model.add_subsystem('core_loss', core_loss(), promotes_inputs=['rm','n_m'], promotes_outputs=['L_core']) + + #total mass + #model.add_subsystem('total_mass', ExecComp('tm = s + r + m')) + #model.connect() + + model.add_subsystem('efficiency', efficiency(), promotes_inputs=['i','rm','L_airg','L_airf','L_bear','L_emag','L_ewir','L_core'], promotes_outputs=['L_total']) - model.connect('balance.rot_or', 'rot_or') - model.connect('size.J', 'balance.lhs:rot_or') - - model.connect('balance.l_st', 'l_st') - model.connect('tq', 'balance.lhs:l_st') - - model.linear_solver = DirectSolver() + # model.connect('balance.rot_or', 'rot_or') + # model.connect('size.j', 'balance.lhs:rot_or') + # + # model.connect('balance.l_st', 'l_st') + # model.connect('tq', 'balance.lhs:l_st') - model.nonlinear_solver = NewtonSolver() - model.nonlinear_solver.options['maxiter'] = 100 - model.nonlinear_solver.options['iprint'] = 0 + #model.linear_solver = DirectSolver() + #model.nonlinear_solver = NewtonSolver() + #model.nonlinear_solver.options['maxiter'] = 100 + #model.nonlinear_solver.options['iprint'] = 2 + + #optimization setup + p.driver = ScipyOptimizeDriver() + p.driver.options['optimizer'] = 'SLSQP' + p.driver.options['disp'] = True + model.add_objective('L_total', ref=1) + model.add_design_var('gap', lower = .001, upper = .004) + model.add_design_var('mot_or', lower = .075, upper = .10) + model.add_design_var('rot_or', lower = .04, upper = .07) + model.add_design_var('l_st', lower = .004, upper = .03) + model.add_constraint('tq', lower = 24, scaler = 1) + model.add_constraint('size.j', upper = 13, scaler = 1) + p.setup() + + p['gap'] = 0.002 + p['mot_or'] = 0.08 + p['rot_or'] = 0.06 + p['l_st'] = 0.008 + # p.final_setup() #p.check_partials(compact_print=True) p.model.list_outputs(implicit=True) - p.set_solver_print(2) + p.set_solver_print(level=2) #view_model(p) p.run_model() - - print('Rotor Outer Radius................', p.get_val('balance.rot_or', units='mm')) + #p.run_driver() + + # print('Rotor Outer Radius................', p.get_val('balance.rot_or', units='mm')) print('Rotor Inner Radius................', p.get_val('rot_ir', units='mm')) - # print('Rotor Outer Radius................', p.get_val('rot_or', units='mm')) - + print('Rotor Outer Radius................', p.get_val('rot_or', units='mm')) + print('Gap...............................', p.get_val('gap', units='mm')) print('Stator Inner Radius...............', p.get_val('sta_ir', units='mm')) - print('Motor Outer Radius................', p.get_val('mass.r_m', units='mm')) + print('Motor Outer Radius................', p.get_val('mass.mot_or', units='mm')) print('Rotor Yoke Thickness..............', p.get_val('w_ry', units='mm')) print('Slot Depth........................', p.get_val('s_d', units='mm')) @@ -304,13 +590,23 @@ def compute(self,inputs,outputs): print('Mass of Stator....................', p.get_val('sta_mass', units='kg')) print('Mass of Rotor.....................', p.get_val('rot_mass', units='kg')) print('Mass of Magnets...................', p.get_val('mag_mass', units='kg')) - print('Current Density...................', p.get_val('size.J')) + print('Current Density...................', p.get_val('size.j')) print('Stack Length......................', p.get_val('mass.l_st', units='mm')) + print('Motor Speed.......................', p.get_val('rm')) + + print('Input Power.......................', p.get_val('efficiency.P_in')) + print('Gap Windage Loss..................', p.get_val('L_airg')) + print('Axial Windage Loss................', p.get_val('L_airf')) + print('Bearing Loss......................', p.get_val('L_bear')) + print('Stator Core Loss..................', p.get_val('L_core')) + print('Magnetic Eddy Loss................', p.get_val('L_emag')) + print('Wiring Eddy Loss..................', p.get_val('L_ewir')) + print('Total Loss........................', p.get_val('L_total')) from solid import * from solid.utils import * - r_m = float(p.get_val('mass.r_m', units='mm')) + mot_or = float(p.get_val('mass.mot_or', units='mm')) l_st = float(p.get_val('mass.l_st', units='mm')) sta_ir = float(p.get_val('sta_ir', units='mm')) w_sy = float(p.get_val('w_sy', units='mm')) @@ -323,9 +619,9 @@ def compute(self,inputs,outputs): fs = 0.05 rot_ir = float(p.get_val('rot_ir', units='mm')) - rot_or = float(p.get_val('balance.rot_or', units='mm')) + rot_or = float(p.get_val('rot_or', units='mm')) - stator_yolk = cylinder(r=r_m, h=l_st, center=True) - cylinder(r=r_m-w_sy, h=l_st+1, center=True) + stator_yolk = cylinder(r=mot_or, h=l_st, center=True) - cylinder(r=mot_or-w_sy, h=l_st+1, center=True) slot = cube([s_d, w_t, l_st], center=True) rotor = color("Blue");cylinder(r=rot_or, h=l_st, center=True) - cylinder(r=rot_ir, h=l_st+1, center=True) @@ -380,4 +676,4 @@ def compute(self,inputs,outputs): } -''' \ No newline at end of file +''' diff --git a/notes.tex b/notes.tex new file mode 100644 index 0000000..8b0ade8 --- /dev/null +++ b/notes.tex @@ -0,0 +1,181 @@ +\documentclass[10pt]{article} +\usepackage{graphicx} +\usepackage{parskip} +\usepackage{amsmath} +\usepackage{float} +\usepackage{listings} +\usepackage{color} +\usepackage{ulem} +\usepackage[margin = 1in]{geometry} + +\definecolor{mygreen}{rgb}{0,0.6,0} +\definecolor{mygray}{rgb}{0.5,0.5,0.5} +\definecolor{mymauve}{rgb}{0.58,0,0.82} + +\lstset{ % + backgroundcolor=\color{white}, % choose the background color + basicstyle=\footnotesize, % size of fonts used for the code + breaklines=true, % automatic line breaking only at whitespace + captionpos=b, % sets the caption-position to bottom + commentstyle=\color{mygreen}, % comment style + escapeinside={\%*}{*)}, % if you want to add LaTeX within your code + keywordstyle=\color{blue}, % keyword style + stringstyle=\color{mymauve}, % string literal style +} + +\begin{document} + +\title{Progress Notes on Motor Design Code} +\author{Garo Bedonian, Aaron Grgurich} + +\maketitle + +\section{NEXT STEP} + +\sout{Talk with Aaron about MotorCAD manual, talk about FEM with Justin,} use total mass as objective + +\sout{Use Tucker's resources to figure out what to do to model loss} + +Hanselman book will be the best resource for this, use it to define loss and other necessary parameters. See specific notes + +Current: Find \sout{$B_r$, $K_a$}, $\mu_{rn}$, $K_{rn}$, $K_{\theta{}n}$ and coefficients for correction factor + +Ask Dustin about the permanent magnets, and their magnetic remanence + +The magnets are (for now) N48H neodymium. + +\sout{Very important: Need the magnetization profile of the halbach array to determine $K_{rn}$, $K_{\theta{}n}$. Will temporarily assume radial orientation.} + +Very important: Need to compute fourier series of slot correction factor function (piecewise) +Need to debug, values explode with higher truncation point + +\sout{Still an issue. Finish re-derivation, give it another shot for a day, then try an FEM approach.} +Justin helped me solve it by solving the linear system for every $n$. Now check the answer and see what I can do with it. + +\sout{Implement slot correction factor. Clean up code.} + +Implement core loss equation. Consider options for determining coefficients. + +Implement model in OpenMDAO + +Finish computing partials for rotor outer radius, gap, alpha, w_tb, tt + +\section{High-Level Notes} + +Goal: Make strides in unifying motor design process using OpenMDAO. + +Ideas/paths for making progress +\begin{itemize} + \item Implement the lowest-fidelity analyses possible (existing surrogate models, analytic stuff) + \item Focus on improving fidelity of one, two, analyses (EM, Thermal) + \item Focus on defining optimization goals (objectives, constraints, design space) + \item Menial tasks: analytic derivatives, clean up code, split across .py files +\end{itemize} + +1. +Seems to be current focus. +Going to look through MotorCAD manual to replicate their models with Aaron + +Communicating with Tucker, will send some material on EM when he can + +Also need to consider thermal and rotordynamic analyses + +2. +Build something in EMpy? +Or MFEM? + +Haven't started on this yet + + +3. +With losses modeled as they are (dependent only on electric frequency and representing a different motor), can't use them as objective. +Want to try minimizing mass instead, see what optimizer does. + +Currently using design space bounds, + +Upper limit on current density + +Lower limit on torque + +Objective is losses, but don't make sense to use as they are + +4. +Do these to pass the time + +Especially need to clean up comments, remove deprecated lines + + + + +\section{Notable Commits} + +ef18621 2019-07-02: Added this notes.tex file + +91029b2 2019-07-02: The surrogate models for magnetic/wiring eddy and core losses work, as do the mechanical loss models. Optimization doesn't make sense as they are. + +bb0d8cd 2019-07-23: Computed magnetic flux density in the stator teeth, expandable to yoke. + +\section{Specific Notes} + +May need to look at magnetic flux density/max flux calculations (assumed constant) See pg. 156 EMF book + +Model magnetic field as a function of the properties of the coils in the motor. See pg. 161 EMF book + +MATlab code modeling three phase rotating magnetic field. See pg. 168 EMF book + +Core loss reductions, high dependence on lamination thickness and resisitivity of material (find this out for Hiperco). See pg. 219 Hanselman + +Amazon links to good books: + +https://www.amazon.com/Electric-Machinery-Fundamentals-Chapman/dp/007108617X + +Hysteresis Losses (Steinmetz Equations): + +\begin{equation} +P_h = k_hfB_{pk}^n +\end{equation} + +or + +\begin{equation} +P_h = k_hfB_{pk}^{n+mB_{pk}} +\end{equation} + +$k_h$, $n$, $m$, are material constants. $B_{pk}$ is peak flux density, + +To be useful for optimization, we need to compute $B$ in general (steady state) + +For $B_{pk}$, + +\begin{equation} +B(\theta) = \sum_{n=-\infty}^{\infty}B_ne^{jn\theta} +\end{equation} + +And $B_o={pk}$ is maximum of $B(\theta)$. Check (7.23) and (7,24) for stator tooth, (7.32) and (7.33) for yoke. (pg. 158 Hanselman) + +\begin{equation} +\phi_i(\alpha) = \sum_{n=-\infty}^{\infty}\bigg\{\frac{2L_{st}R_s}{N_m}B_{gn}\int_{-\theta_s/2}^{\theta_s/2}K_{sl}(\theta)e^{jn\theta}d\theta{}\bigg\}e^{jn\alpha} +\end{equation} + +But I need to see if the magnetic field is sinusoidal. + +Halbach array is said to have sinusoidal magnetic field profile.??? + +Model appears to suffer from subtractive cancellation due to exponential increase with $n$. Higher frequency terms destroy accuracy. + + + +Caveats to keep in mind: + +Model suffers from subtractive cancellation error for large $\beta$, limiting accuracy of solution. Should not severely impact solution. + +Model depends on knowledge of the magnetic field profile. Unsure how to characterize Halbach array, but some options are implemented. + +Need to derive material constants for Hiperco 50 so as to fit equation 9.39, then use 9.40 to generally compute loss. Material data sheets may have the data to do this, but otherwise may require experimentation. Current numbers are a guess. + +Model currently only computes core losses only within the stator teeth, but we can now quickly extend to other relevant losses. + +Datasheet for Hiperco 50, includes data for core losses +https://edfagan.com/litPDF/Hiperco%2050A%20Alloy%20Data%20Sheet.pdf + +\end{document} diff --git a/radial_air_fourier_comp.py b/radial_air_fourier_comp.py new file mode 100644 index 0000000..7d24d90 --- /dev/null +++ b/radial_air_fourier_comp.py @@ -0,0 +1,591 @@ +from __future__ import absolute_import +import numpy as np +from math import pi +import openmdao.api as om + +N = 4 +M = 20 + +B_tna = np.zeros(N) + +#need partials for rot_or, sta_ir, t_mag, alpha, w_tb, tt + +class DiffAssemble(om.ExplicitComponent): + def initialize(self): + self.options.declare('n', default=1., desc='index') + + def setup(self): + self.add_input('n_m', 16, desc='number of poles') + self.add_input('sta_ir', .51, desc='stator inner radius') + self.add_input('rot_or', .50, desc='rotor outer radius w/ magnets') + self.add_input('t_mag', .03, desc='magnet thickness') + self.add_input('b_r', 1.3, desc='magnetic remanence') + self.add_input('mu_r', 1.05, desc='relative recoil permitivity') + self.add_input('alpha', .85, desc='magnet fraction') + + self.add_output('bmat', shape=(3,3), desc='system matrix') + self.add_output('rhs', shape=3, desc='right hand side') + + self.declare_partials('bmat','sta_ir') + self.declare_partials('bmat','rot_or') + self.declare_partials('bmat','t_mag') + self.declare_partials('bmat','alpha') + self.declare_partials('rhs','sta_ir') + self.declare_partials('rhs','rot_or') + self.declare_partials('rhs','t_mag') + self.declare_partials('rhs','alpha') + + def compute(self, inputs, outputs): + n = self.options['n'] + n_p = inputs['n_m']/2 #number of pole pairs + r_s = inputs['sta_ir'] #stator inner radius + r_m = inputs['rot_or'] #rotor outer radius (INCLUDING MAGNETS) + r_r = r_m - inputs['t_mag'] #rotor outer radius + b_r = inputs['b_r'] #magnetic remanence of stator material + mu_r = inputs['mu_r'] #magnetic relative recoil permitivity + alpha = inputs['alpha'] #magnet fraction in rotor + + mu_0 = 0.001 + + 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 + outputs['bmat'] = bmat + outputs['rhs'] = rhs + + def compute_partials(self, inputs, J): + n = self.options['n'] + n_p = inputs['n_m']/2 #number of pole pairs + r_s = inputs['sta_ir'] #stator inner radius + r_m = inputs['rot_or'] #rotor outer radius (INCLUDING MAGNETS) + r_r = r_m - inputs['t_mag'] #rotor outer radius + b_r = inputs['b_r'] #magnetic remanence of stator material + mu_r = inputs['mu_r'] #magnetic relative recoil permitivity + alpha = inputs['alpha'] #magnet fraction in rotor + + drrdrot = 1 + drrdtma = -1 + + mu_0 = 0.001 + + bmat = np.zeros((3,3)) + rhs = np.zeros(3) + dbmatdsta = np.zeros((3, 3)) + dbmatdrot = np.zeros((3, 3)) + dbmatdtma = np.zeros((3, 3)) + dbmatdalp = np.zeros((3, 3)) + drhsdsta = np.zeros(3) + drhsdrot = np.zeros(3) + drhsdtma = np.zeros(3) + drhsdalp = np.zeros(3) + + beta = n*n_p + + # magnetic profile coefficients + k_rn = RadialR(n, alpha) + dkrndalp = RadialRPart(n, alpha) + + k_tn = RadialT(n, alpha) + dktndalp = 0 + + # assume output is real, ignore imaginary part (Equation B.18) + Cm = b_r*(k_rn)/((1-beta**2)*mu_r*mu_0) + dCmdkrn = b_r/((1-beta**2)*mu_r*mu_0) + + #Hatheta boundary condition (excluded from linear system) + Ea = -(r_s)**(2*beta) + dEadsta = -(2*beta)*(r_s**(2*beta-1)) + + #Hmtheta boundary condition + bmat[0, 1] = (r_r**(-beta-1)) + dbmatdrot[0, 1] = (-beta-1)*(r_r**(-beta-2))*drrdrot + dbmatdtma[0, 1] = (-beta-1)*(r_r**(-beta-2))*drrdtma + + bmat[0, 2] = (r_r**(beta-1)) + dbmatdrot[0, 2] = (beta-1)*(r_r**(beta-2))*drrdrot + dbmatdtma[0, 2] = (beta-1)*(r_r**(beta-2))*drrdtma + + rhs[0] = -Cm + drhsdalp[0] = -1*dCmdkrn*dkrndalp + + #Htheta boundary condition + bmat[1, 0] = -Ea*(r_m**(-beta-1)) - r_m**(beta - 1) + dbmatdrot[1, 0] = -Ea*(-beta-1)*(r_m**(-beta-2)) - (beta - 1)*(r_m**(beta - 2)) + dbmatdsta[1, 0] = -(r_m**(-beta-1))*dEadsta + + bmat[1, 1] = r_m**(-beta-1) + dbmatdrot[1, 1] = (-beta-1)*(r_m**(-beta-2)) + + bmat[1, 2] = r_m**(beta-1) + dbmatdrot[1, 2] = (beta-1)*(r_m**(beta-2)) + + rhs[1] = -Cm + drhsdalp[1] = -1*dCmdkrn*dkrndalp + + #Hr boundary condition + bmat[2, 0] = (beta/mu_r)*(Ea*(r_m**(-beta-1)) - r_m**(beta - 1)) + dbmatdrot[2, 0] = (beta/mu_r)*(Ea*(-beta-1)*(r_m**(-beta-2)) - (beta - 1)*(r_m**(beta - 2))) + dbmatdsta[2, 0] = (beta/mu_r)*(r_m**(-beta-1))*dEadsta + + bmat[2, 1] = -beta*(r_m**(-beta-1)) + dbmatdrot[2, 1] = -beta*(-beta-1)*(r_m**(-beta-2)) + + bmat[2, 2] = beta*(r_m**(beta-1)) + dbmatdrot[2, 2] = beta*(beta-1)*(r_m**(beta-2)) + + + rhs[2] = (b_r*k_rn/(mu_r*mu_0)) - Cm + drhsdalp[2] = (b_r/(mu_r*mu_0))*dkrndalp - 1*dCmdkrn*dkrndalp + + J['bmat','sta_ir'] = dbmatdsta + J['bmat','rot_or'] = dbmatdrot + J['bmat','t_mag'] = dbmatdtma + J['bmat','alpha'] = dbmatdalp + J['rhs','sta_ir'] = drhsdsta + J['rhs','rot_or'] = drhsdrot + J['rhs','t_mag'] = drhsdtma + J['rhs','alpha'] = drhsdalp + + +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 RadialRPart(n, alpha): + + if n%2 != 0: + return n*(np.pi/2)*np.cos(n*alpha*np.pi/2)/(n*np.pi/2) + + if n%2 == 0: + return 0 + +def RadialT(n, alpha): + return 0 + +# class DiffLinearSolve(om.ImplicitComponent): +# +# def setup(self): +# self.add_input('bmat', shape=(3,3), desc='system matrix') +# self.add_input('rhs', shape=3, desc='right hand side') +# +# self.add_output('x', shape=3) +# +# +# def apply_nonlinear(self, inputs, outputs, residuals): +# bmat = inputs['bmat'] +# rhs = inputs['rhs'] +# x = outputs['x'] +# print(bmat) +# +# residuals['x'] = np.dot(bmat, x) - rhs +# + + # def solve_linear(self, outputs, residuals, mode): + # if mode == 'fwd': + # outputs['x'] = self.inv_jac * residuals['x'] + # elif mode == 'rev': + # residuals['x'] = self.inv_jac * outputs['x'] + + + +class BarnComp(om.ExplicitComponent): + def initialize(self): + self.options.declare('n', default=1., desc='index') + + def setup(self): + self.add_input('x', shape=3, desc='assumed solution coeffs') + self.add_input('sta_ir', .51, desc='stator inner radius') + self.add_input('n_m', desc='number of poles') + + self.add_output('B_arn') + + self.declare_partials('B_arn','x') + self.declare_partials('B_arn','sta_ir') + def compute(self, inputs, outputs): + x = inputs['x'] + n = self.options['n'] + r_s = inputs['sta_ir'] #stator inner radius + n_p = inputs['n_m']/2 #number of pole pairs + mu_0 = 0.001 + + Da = -x[0] + beta = n*n_p + Ea = -(r_s)**(2*beta) + b_arn = mu_0*beta*Da*(Ea*(r_s**(-beta-1)) - r_s**(beta-1)) #Equation B.13 + outputs['B_arn'] = b_arn + + def compute_partials(self, inputs, J): + x = inputs['x'] + n = self.options['n'] + r_s = inputs['sta_ir'] #stator inner radius + n_p = inputs['n_m']/2 #number of pole pairs + mu_0 = 0.001 + + Da = -x[0] + dDadx = [-1, 0, 0] + beta = n*n_p + Ea = -(r_s)**(2*beta) + dEadsta = -(2*beta)*(r_s**(2*beta-1)) + + b_arn = mu_0*beta*Da*(Ea*(r_s**(-beta-1)) - r_s**(beta-1)) #Equation B.13 + + J['B_arn','x'] = mu_0*beta*(Ea*(r_s**(-beta-1)))*dDadx + J['B_arn','sta_ir'] = mu_0*beta*Da*((r_s**(-beta-1))*dEadsta + Ea*(-beta-1)*(r_s**(-beta-2)) - (beta - 1)*(r_s**(beta-2))) + +class SCFactor(om.ExplicitComponent): + + def setup(self): + self.add_input('n_m', 16, desc='number of poles') + self.add_input('k_st', desc='lamination stacking factor') + self.add_input('w_tb', desc='tooth body width') + self.add_input('sta_ir', .51, desc='stator inner radius') + self.add_input('rot_or', .50, desc='rotor outer radius w/ magnets') + self.add_input('t_mag', .03, desc='magnet thickness') + self.add_input('n_s', 24, desc='number of slots') + self.add_input('tt', .6*np.pi/10, desc='tooth width angle') + self.add_input('mu_r', 1.05, desc='relative recoil permitivity') + self.add_output('K_slm', shape=2*M + 2, desc='slot correction factor coeffs') + + def compute(self, inputs, outputs): + n_p = inputs['n_m']/2 #number of pole pairs + k_st = inputs['k_st'] #lamination stacking factor + w_tb = inputs['w_tb'] #tooth body width + r_s = inputs['sta_ir'] #stator inner radius + r_m = inputs['rot_or'] #rotor outer radius (INCLUDING MAGNETS) + t_mag = inputs['t_mag'] #magnet thickness + n_s = inputs['n_s'] #max angle occupied by stator tooth + tt = inputs['tt'] #tooth width angle + mu_r = inputs['mu_r'] #magnetic relative recoil permeability + + ts = 2*np.pi/n_s + gap = r_s - r_m #air gap thickness + T = ts + + outputs['K_slm'] = FourierCoeffsNumpy(kfunc, T, M, gap, r_s, inputs['n_m'], tt, ts, mu_r, t_mag) + + + # def compute_partials(self, inputs, outputs): + # i = 1 + +def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): + + fsample = 2*N + t = np.linspace(-T/2, T/2, fsample + 2, endpoint = False) + y = np.fft.rfft(func(t, gap, r_s, n_m, tt, ts, mu_r, t_mag))/fsample + + c = y + for k in range(1, N + 1): + c = np.insert(c, 0, np.conj(y[k])) + + c = c.real + return c + +# Equation 7.9 +def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): + + gratio = gfunc(theta, gap, r_s, n_m, tt, ts) + + Ksl = (1 + t_mag/(gap*mu_r))/(gratio + t_mag/(gap*mu_r)) + + return Ksl + +# Equation 7.11 +def gfunc(theta, gap, r_s, n_m, tt, ts): + + gratio = np.zeros(len(theta)) + for x in range(0, len(theta)): + if abs(theta[x]) <= tt/2: + gratio[x] = 1 + + if theta[x] <= ts/2 and theta[x] >= tt/2: + gratio[x] = 1 + np.pi*r_s*(theta[x]-tt/2)/(n_m*gap) + + if theta[x] <= -tt/2 and theta[x] >= -ts/2: + gratio[x] = 1 - np.pi*r_s*(theta[x]+tt/2)/(n_m*gap) + + + return gratio + +class CorrectSlot(om.ExplicitComponent): + + def initialize(self): + self.options.declare('n', default=1., desc='index') + + def setup(self): + self.add_input('K_slm', shape=2*M+2, desc='slot correction coeffs') + self.add_input('B_arn', desc='B_arn') + self.add_input('n_m', desc='number of poles') + self.add_input('n_s', desc='number of slots') + self.add_input('sta_ir') + self.add_input('k_st') + self.add_input('w_tb') + + self.add_output('B_tn', shape=N, desc='corrected B_arn terms') + + self.declare_partials('B_tn','K_slm') + self.declare_partials('B_tn','B_arn') + self.declare_partials('B_tn','sta_ir') + self.declare_partials('B_tn','w_tb') + + def compute(self, inputs, outputs): + n = self.options['n'] + K_slm = inputs['K_slm'] + B_arn = inputs['B_arn'] + n_m = inputs['n_m'] #number of pole pairs + n_s = inputs['n_s'] + r_s = inputs['sta_ir'] + k_st = inputs['k_st'] + w_tb = inputs['w_tb'] + + 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 = B_arn*(2*np.pi*r_s/(n_s*k_st*w_tb))*slotsum + + B_tna[n] = B_tn + outputs['B_tn'] = B_tna + + def compute_partials(self, inputs, J): + n = self.options['n'] + K_slm = inputs['K_slm'] + B_arn = inputs['B_arn'] + n_m = inputs['n_m'] #number of pole pairs + n_s = inputs['n_s'] + r_s = inputs['sta_ir'] + k_st = inputs['k_st'] + w_tb = inputs['w_tb'] + + dssdksl = np.zeros(2*M + 1) + + 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 + + if n%2 != 0: + for m in range(-M, M + 1): + dssdksl[m + M] = np.sinc(np.pi*(m + n*(n_m/(2*n_s))))#/(np.pi*(m + n*(n_m/(2*n_s)))) + else: + dssdksl = 0 + + dBtndBarn = (2*np.pi*r_s/(n_s*k_st*w_tb))*slotsum + dBtndksl = B_arn*(2*np.pi*r_s/(n_s*k_st*w_tb))*dssdksl + dBtndsta = B_arn*(2*np.pi/(n_s*k_st*w_tb))*slotsum + dBtndwtb = -B_arn*(2*np.pi*r_s/(n_s*k_st*(w_tb**2)))*slotsum + + J['B_tn','K_slm'] = dBtndBarn + J['B_tn','B_arn'] = dBtndksl + J['B_tn','sta_ir'] = dBtndsta + J['B_tn','w_tb'] = dBtndwtb + +class MotorAirFourierComp(om.Group): + def initialize(self): + self.options.declare('ng', default=1., desc='index') + + def setup(self): + ng = self.options['ng'] + self.add_subsystem('DiffAssemble', DiffAssemble(n = ng), promotes_inputs=['n_m','sta_ir','rot_or','t_mag','b_r','mu_r','alpha']) + + lingrp = self.add_subsystem('lingrp', om.Group(), promotes=['*']) + lingrp.add_subsystem('DiffLinearSolve', om.LinearSystemComp(size=3)) + self.connect('DiffAssemble.bmat', 'DiffLinearSolve.A') + self.connect('DiffAssemble.rhs', 'DiffLinearSolve.b') + + + self.add_subsystem('BarnComp', BarnComp(n = ng), promotes_inputs=['sta_ir','n_m'], promotes_outputs=['B_arn']) + self.connect('DiffLinearSolve.x','BarnComp.x') + + self.add_subsystem('CorrectSlot', CorrectSlot(n = ng), promotes_inputs=['K_slm','B_arn','n_m','n_s','sta_ir', 'k_st', 'w_tb'], promotes_outputs=['B_tn']) + + #self.linear_solver = om.DirectSolver() + +class InvFourier(om.ExplicitComponent): + + def setup(self): + self.add_input('B_tn', shape=N) + + self.add_output('B_t', shape=(N-1)*2) + + self.declare_partials('B_t','B_tn', method='fd') + def compute(self, inputs, outputs): + print(inputs['B_tn']) + B_t = np.fft.irfft(inputs['B_tn']) + outputs['B_t'] = B_t + +class MaxMeanSquared(om.ExplicitComponent): + + def setup(self): + self.add_input('B_t', shape=(N-1)*2) + self.add_input('f') + + self.add_output('B_pk') + self.add_output('B_msq') + + self.declare_partials('B_pk','B_t') + self.declare_partials('B_msq','B_t') + + def compute(self, inputs, outputs): + B_t = inputs['B_t'] + f = inputs['f'] + + omega = f*(2*np.pi/60) + + # peak magnetic flux density is the largest data point from an inverse fourier transform + B_pk = max(B_t) + + Nl = len(B_t) + #assume B_n is real, compute mean-square value of dB/dt + B_msq = 0 + for t in range(1, Nl): + B_msq += 2*((t*omega*B_t[t])**2) + + outputs['B_pk'] = B_pk + outputs['B_msq'] = B_msq + + def compute_partials(self, inputs, outputs): + B_t = inputs['B_t'] + f = inputs['f'] + + omega = f*(2*np.pi/60) + + dBpkdBt = np.zeros(len(B_t)) + dBpkdBt[B_t.index(max(B_t))] = 1 + + dBmsqdBt = np.zeros(len(B_t)) + for t in range(1, Nl): + dBmsqdBt[t] = 4*(((t*omega)**2)*B_t[t]) + + J['B_pk', 'B_t'] = dBpkdBt + J['B_msq', 'B_t'] = dBmsqdBt + +class CoreLoss(om.ExplicitComponent): + + def setup(self): + self.add_input('B_pk') + self.add_input('B_msq') + self.add_input('f') + + self.add_output('L_core') + + self.declare_partials('L_core', 'B_pk') + self.declare_partials('L_core', 'B_msq') + + def compute(self, inputs, outputs): + B_pk = inputs['B_pk'] + B_msq = inputs['B_msq'] + f = inputs['f'] + # 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)) + + outputs['L_core'] = L_core + + def compute_partials(self, inputs, J): + B_pk = inputs['B_pk'] + B_msq = inputs['B_msq'] + f = inputs['f'] + + n = 2 + m = 0.01 + k_h = 0.5 + k_e = 0.5 + + J['L_core', 'B_msq'] = k_e/(2*(np.pi**2)) + J['L_core', 'B_pk'] = k_h*f*(B_pk**(n + m*B_pk))*(m*np.ln(B_pk) + (n + m*B_pk)/B_pk) + +if __name__ == '__main__': + p = om.Problem() + model = p.model + + ind = model.add_subsystem('indeps', om.IndepVarComp(), promotes=['*']) + + ind.add_output('t_mag', val=.003)#, units='m') # Magnet thickness + ind.add_output('n_s', val=24) # Number of Slots + ind.add_output('n_m', val=16) # Number of poles + ind.add_output('rot_or', val = 0.050)#, units='m') + ind.add_output('sta_ir', val = 0.051)#, units='m') + ind.add_output('f', val = 5400/20)#, units='rpm') + ind.add_output('b_r', val = 1.3) + ind.add_output('mu_r', val = 1.05) + ind.add_output('alpha', val = .85) + ind.add_output('k_st', val = 1) + ind.add_output('w_tb', val = .005) + ind.add_output('tt', val = .6*np.pi/10) + + #model.nonlinear_solver = om.NewtonSolver() + #model.nonlinear_solver.options['solve_subsystems'] = True + + #model.linear_solver = om.DirectSolver() + model.add_subsystem('SCFactor', SCFactor(), promotes_inputs=['n_m', 'k_st', 'w_tb', 'sta_ir', 'rot_or', 't_mag', 'n_s', 'tt', 'mu_r'], promotes_outputs=['K_slm']) + + for nl in range(1, N): + model.add_subsystem('MotorAirFourierComp{}'.format(nl), MotorAirFourierComp(ng=nl), promotes_inputs=['n_m','sta_ir','rot_or','t_mag','b_r','mu_r','alpha', 'K_slm', 'n_s','w_tb','k_st']) + + model.add_subsystem('InvFourier', InvFourier(), promotes_outputs=['B_t']) + + model.connect('MotorAirFourierComp{}.B_tn'.format(nl), 'InvFourier.B_tn') + + model.add_subsystem('MaxMeanSquared', MaxMeanSquared(), promotes_inputs=['B_t', 'f'], promotes_outputs=['B_pk', 'B_msq']) + model.add_subsystem('CoreLoss', CoreLoss(), promotes_inputs=['B_pk', 'B_msq', 'f'], promotes_outputs=['L_core']) + + p.setup() + + p.run_model() + + print(p.get_val('L_core')) diff --git a/slotcorrection.py b/slotcorrection.py new file mode 100644 index 0000000..c1da8b8 --- /dev/null +++ b/slotcorrection.py @@ -0,0 +1,84 @@ +''' +Helps to compute the Fourier coefficients of the stator tooth flux density by computing a correction factor for the slots + +Based on the radial air gap solution as solved in bmatrix.py + +Based on section 7.3 of 'Brushless Permanent Magnet Motor Design' by Duane P. Hanselman +''' + + +import numpy as np +import matplotlib.pyplot as plt + +# Compute slot correction Fourier terms using model presented in Chapter 7 +def slot_correction(M, params): + 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 + r_m = params['r_m'] #rotor outer radius (INCLUDING MAGNETS) + t_mag = params['t_mag'] #magnet thickness + tt = params['tt'] #max angle occupied by stator tooth + ts = params['ts'] #angle of separation between stator teeth + mu_r = params['mu_r'] #magnetic relative recoil permeability + + gap = r_s - r_m #air gap thickness + n_m = n_p*2 #number of poles + T = ts #period of solution, equivalent to the angle of separation between teeth + + + Kslm = FourierCoeffsNumpy(kfunc, T, M, gap, r_s, n_m, tt, ts, mu_r, t_mag) + + + return Kslm + +# Compute Fouerier coefficients of a real valued function, arrange in linear order +def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): + + fsample = 2*N + t = np.linspace(-T/2, T/2, fsample + 2, endpoint = False) + y = np.fft.rfft(func(t, gap, r_s, n_m, tt, ts, mu_r, t_mag))/fsample + + c = y + for k in range(1, N + 1): + c = np.insert(c, 0, np.conj(y[k])) + + c = c.real + return c + +# Equation 7.9 +def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): + + gratio = gfunc(theta, gap, r_s, n_m, tt, ts) + + Ksl = (1 + t_mag/(gap*mu_r))/(gratio + t_mag/(gap*mu_r)) + + return Ksl + + +# Equation 7.11 +def gfunc(theta, gap, r_s, n_m, tt, ts): + + gratio = np.zeros(len(theta)) + for x in range(0, len(theta)): + if abs(theta[x]) <= tt/2: + gratio[x] = 1 + + if theta[x] <= ts/2 and theta[x] >= tt/2: + gratio[x] = 1 + np.pi*r_s*(theta[x]-tt/2)/(n_m*gap) + + if theta[x] <= -tt/2 and theta[x] >= -ts/2: + gratio[x] = 1 - np.pi*r_s*(theta[x]+tt/2)/(n_m*gap) + + + return gratio + +if __name__ == '__main__': + M = 30 + K_slm = np.zeros(M) + + #for m in range(0, M): + K_slm = slot_correction(M) + #print(K_slm) + #plt.plot(K_slm) + plt.show() diff --git a/toothflux.py b/toothflux.py new file mode 100644 index 0000000..2040af3 --- /dev/null +++ b/toothflux.py @@ -0,0 +1,212 @@ +''' +This function computes magnetic flux 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 +from math import pi, sin, cos, exp +import matplotlib.pyplot as plt + +def toothflux(theta, nterms, params, Bn): + n_m = params['n_m'] #number of poles + n_s = params['n_s'] #number of Slots + l_st = params['l_st'] #axial motor length + k_st = params['k_st'] #lamination stacking factor + w_tb = params['w_tb'] #tooth body width + r_s = params['sta_ir'] #stator inner radius + t_mag = params['t_mag'] #magnet thickness + r_m = params['rot_or'] #rotor outer radius (INCLUDING MAGNETS) + gap = params['gap'] #air gap + t_mag = params['t_mag'] #magnet thickness + tt = params['tt'] + ts = params['ts'] + mu_r = params['mu_r'] + alpha = params['alpha'] + #alpha = 18/30 #magnetic fraction, fraction of rotor circumference occupied by one magnet + + r_r = r_m - t_mag + + b_r = 1.3 #magnetic remanence + + #nterms = 10 #flux fourier series terms + mterms = 20 #slot correction fourier series terms + + n_p = n_m/2 #number of pole pairs + + + + Kslm = np.zeros(2*mterms + 1) + Bgn = np.zeros(2*nterms + 1, dtype=np.complex_) + phi = np.zeros(2*nterms + 1) + Bn = np.zeros(2*nterms + 1, dtype=np.complex_) + + #mu_r = .001 #magnet relative recoil permeability, aka relative permitivity, mu/mu_0 + + phi = 0 + #flux term layer + for n in range(1, nterms): + beta = n*n_p + + slotsum = 0 + + ###slot correction factor function, to be used to find + #first determine g(theta)/g + #gratio = gfunc(theta, gap, r_s, n_m, tt, ts) + + #Ksl = kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, l_m) + + #fourier transform to obtain coefficients in Kslm + #compute using Riemann sum + + Kslm = FourierCoeffsNumpy(kfunc, T, mterms, gap, r_s, n_m, tt, ts, mu_r, t_mag) + + #slot correction term layer + if n%2 != 0: + for m in range(-mterms, mterms + 1): + + slotsum += Kslm[m + mterms]*sin(pi*(m + n*(n_m/(2*n_s))))/(pi*(m + n*(n_m/(2*n_s)))) + + else: + slotsum = 0 + print(slotsum) + # + # add k_rn and k_tn here based on magnetic field orientation of halbach array + # + ###for now, assume radial orientation + k_rn = RadialR(n, alpha) + #print(k_rn) + 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 + + #print(k_a) + #k_a = k_a.real + #print(n) + Bgn[n+nterms] = -b_r*k_a*(((r_s/r_m)**(beta-1)) + ((r_s/r_m)**(2*beta))*((r_m/r_s)**(beta+1))) + Btn = 0 + #print(slotsum) + Bn[n+nterms] = Bgn[n+nterms]*(2*pi*l_st*r_s/n_s)*slotsum/(k_st*l_st*w_tb) + #print(Bgn[n+nterms]) + #print(Bgn) + return Bgn, Btn + +def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): + + gratio = gfunc(theta, gap, r_s, n_m, tt, ts) + + Ksl = (1 + t_mag/(gap*mu_r))/(gratio + t_mag/(gap*mu_r)) + + return Ksl + + + +def gfunc(theta, gap, r_s, n_m, tt, ts): + + gratio = np.zeros(len(theta)) + for x in range(0, len(theta)): + if abs(theta[x]) <= tt/2: + gratio[x] = 1 + + if theta[x] <= ts/2 and theta[x] >= tt/2: + gratio[x] = 1 + pi*r_s*(theta[x]-tt/2)/(n_m*gap) + + if theta[x] <= -tt/2 and theta[x] >= -ts/2: + gratio[x] = 1 - pi*r_s*(theta[x]+tt/2)/(n_m*gap) + + + return gratio + +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 + +def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): + #func: one-dimensional function of interest + #T: period of function, will compute from center + #N: number of terms 1 and onward, this computes N + 1 coefficients + + fsample = 2*N + t = np.linspace(-T/2, T/2, fsample + 2, endpoint = False) + y = np.fft.rfft(func(t, gap, r_s, n_m, tt, ts, mu_r, t_mag))/fsample + + c = y + for k in range(1, N + 1): + c = np.insert(c, 0, np.conj(y[k])) + + c = c.real + return c + +if __name__ == '__main__': + + gap = .001 + r_s = .051 + n_m = 16 + 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 + rot_or = .05 + alpha = .85 + + params = { + 'gap':gap, + 'n_m':n_m, + '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, + 'rot_or':(r_s - gap), + 'alpha':alpha + } + + theta = 0#np.linspace(-ts/2,ts/2) + + T = ts + N = 18 + + # x = kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag) + # y = FourierCoeffsNumpy(kfunc, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag) + + Bn = np.zeros(2*N + 1) + Bm = np.zeros(2*N + 1) + + Bn, Bm = toothflux(theta, N, params, Bn) + + #print(Bn) + + #Brn = np.concatenate((Bn[N:], Bn[:N])) + Brn = Bn[N:] + #Btn = np.concatenate((Bm[10:], Bm[:10])) + #print(Brn) + #Brn = Bn[N:] + i = np.fft.irfft(Brn, N) + + plt.plot(i) + #plt.show() + + # print(x) + # print(y)