From 886f3f3d3a04558e7f862e9151c04ebdf08bcf6d Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Fri, 14 Jun 2019 15:18:16 -0400 Subject: [PATCH 01/46] Notes and Action Items --- Aaron_Notes.txt | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 Aaron_Notes.txt 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 From 574c54393948248f4b593861ff2c1e2bef627bea Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 18 Jun 2019 14:35:05 -0400 Subject: [PATCH 02/46] partials for mass --- motor_weight_optimization.py | 40 ++++++++++++++++++++++++++++-------- 1 file changed, 32 insertions(+), 8 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index dee46ef..5db9b06 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -218,13 +218,37 @@ def compute(self,inputs,outputs): # 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'] = + # s_d=inputs['s_d'] + + # J['sta_mass', 'rho'] = l_st * ((pi * r_m**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*s_d*1.5))) + # J['sta_mass', 'r_m'] = 2 * rho * l_st * (pi * r_m) + # 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 * r_m**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*s_d*1.5))) + # J['sta_mass', 's_d'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) + + # 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() @@ -380,4 +404,4 @@ def compute(self,inputs,outputs): } -''' \ No newline at end of file +''' From 7f4bc6424d174fb6815550f36216e4d0628d59a9 Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Tue, 18 Jun 2019 14:41:33 -0400 Subject: [PATCH 03/46] Started work on partials for motor_size --- motor_weight_optimization.py | 83 +++++++++++++++++++----------------- 1 file changed, 43 insertions(+), 40 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index dee46ef..710cfe0 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -73,50 +73,53 @@ def compute(self,inputs,outputs): # 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'] + 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) + # TODO: Get this partial working: + # Use: check_partials function to check: + def compute_partials(self, inputs, J): + gap = inputs['gap'] # Air Gap + + # 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) class torque(ExplicitComponent): @@ -238,7 +241,7 @@ def compute(self,inputs,outputs): ind.add_output('gap', val=0.001, units='m') # Stacking factor 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('r_m', 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 !! From 06ed44cf218788c12c91e8fbf633156a3d0c5731 Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Tue, 18 Jun 2019 15:08:46 -0400 Subject: [PATCH 04/46] Naming convention; slight fix --- motor_weight_optimization.py | 58 ++++++++++++++++++------------------ 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 925774a..8cce084 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -8,7 +8,7 @@ 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 @@ -48,7 +48,7 @@ def setup(self): 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,28 +69,28 @@ 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'] + # 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'] + 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) - +""" # TODO: Get this partial working: # Use: check_partials function to check: def compute_partials(self, inputs, J): gap = inputs['gap'] # Air Gap # rotor_outer_radius - r_m = inputs['r_m'] + mot_or = inputs['mot_or'] 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 + J['rot_or', 'mot_or'] = 1-w_sy-s_d-gap + J['rot_or', 's_d'] = mot_or - w_sy - 1 - gap + J['rot_or', 'w_sy'] = mot_or - 1 - s_d - gap # rotor_yoke_width rot_or = inputs['rot_or'] @@ -120,7 +120,7 @@ def compute_partials(self, inputs, J): 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) - +""" class torque(ExplicitComponent): @@ -165,14 +165,14 @@ 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') @@ -188,14 +188,14 @@ def setup(self): 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'] @@ -216,19 +216,19 @@ 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'] - # J['sta_mass', 'rho'] = l_st * ((pi * r_m**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*s_d*1.5))) - # J['sta_mass', 'r_m'] = 2 * rho * l_st * (pi * r_m) + # 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 * r_m**2)-(pi * (sta_ir+s_d)**2)+(n_s*(w_t*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'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) # rotor @@ -265,7 +265,7 @@ def compute(self,inputs,outputs): ind.add_output('gap', val=0.001, units='m') # Stacking factor 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.08, 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 !! @@ -285,12 +285,12 @@ def compute(self,inputs,outputs): 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('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','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', 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.connect('balance.rot_or', 'rot_or') @@ -318,7 +318,7 @@ def compute(self,inputs,outputs): # print('Rotor Outer Radius................', p.get_val('rot_or', 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')) @@ -337,7 +337,7 @@ def compute(self,inputs,outputs): 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')) @@ -352,7 +352,7 @@ def compute(self,inputs,outputs): rot_ir = float(p.get_val('rot_ir', units='mm')) rot_or = float(p.get_val('balance.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) From 01be6bb5dc06249653d148d1e1eab99eee1b888d Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Tue, 18 Jun 2019 15:46:27 -0400 Subject: [PATCH 05/46] rot_or (rotor outer radius) is a independent variable - no need for partial to be taken --- motor_weight_optimization.py | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 8cce084..93c4b30 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -41,7 +41,7 @@ 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') @@ -76,21 +76,12 @@ def compute(self,inputs,outputs): outputs['rot_ir'] = (rot_or - t_mag) - outputs['w_ry'] outputs['sta_ir'] = rot_or + gap 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) + outputs['j'] = 2*n*i*(2.**0.5)/(k_wb/n_s*(area-n_s*1.25*(outputs['w_t']*outputs['s_d']))*1E6) + # TODO: Better name for current density??? -""" # TODO: Get this partial working: # Use: check_partials function to check: def compute_partials(self, inputs, J): - gap = inputs['gap'] # Air Gap - - # rotor_outer_radius - mot_or = inputs['mot_or'] - s_d = inputs['s_d'] - w_sy = inputs['w_sy'] - J['rot_or', 'mot_or'] = 1-w_sy-s_d-gap - J['rot_or', 's_d'] = mot_or - w_sy - 1 - gap - J['rot_or', 'w_sy'] = mot_or - 1 - s_d - gap # rotor_yoke_width rot_or = inputs['rot_or'] @@ -120,7 +111,6 @@ def compute_partials(self, inputs, J): 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) -""" class torque(ExplicitComponent): From f1c00f652c6be0666b6c5859309baf907340b2a7 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 18 Jun 2019 15:58:56 -0400 Subject: [PATCH 06/46] small edits --- motor_weight_optimization.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 5db9b06..bcdbbf0 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -150,12 +150,12 @@ def compute(self,inputs,outputs): # 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 + # 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): From b221e626d923aab0a2e89c9bdaa4f9d570f3c49d Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 19 Jun 2019 09:55:30 -0400 Subject: [PATCH 07/46] more derivatives --- motor_weight_optimization.py | 44 ++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 127052d..0ea3e21 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -91,7 +91,7 @@ def compute_partials(self, inputs, J): 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', '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) @@ -99,7 +99,7 @@ def compute_partials(self, inputs, J): 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', '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) @@ -111,6 +111,46 @@ def compute_partials(self, inputs, J): 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'] = + J['rot_ir', 't_mag'] = + J['rot_ir', 'b_g'] = + J['rot_ir', 'n_m'] = + J['rot_ir', 'k'] = + J['rot_ir', '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'] + J['j', ''] = + J['j', ''] = + J['j', ''] = + J['j', ''] = + J['j', ''] = + J['j', ''] = + J['j', ''] = + J['j', ''] = + J['j', ''] = + J['j', ''] = class torque(ExplicitComponent): From 8de10ee1ffe04693575b6767ef025ee25c25203d Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Wed, 19 Jun 2019 10:32:26 -0400 Subject: [PATCH 08/46] Partials; Rotor Inner Radius - Motor Size --- motor_weight_optimization.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 0ea3e21..64ab6c5 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -125,12 +125,12 @@ def compute_partials(self, inputs, J): # rotor_inner_radius t_mag = inputs['t_mag'] - J['rot_ir', 'rot_or'] = - J['rot_ir', 't_mag'] = - J['rot_ir', 'b_g'] = - J['rot_ir', 'n_m'] = - J['rot_ir', 'k'] = - J['rot_ir', 'b_ry'] = + J['rot_ir', 'rot_or'] = 1 - J['w_ry', 'rot_or'] + J['rot_ir', 't_mag'] = -1 - J['w_ry', 't_mag'] + 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 From 20db8cd9a317c27929cfea4c8ce8388a008ba827 Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Wed, 19 Jun 2019 11:08:47 -0400 Subject: [PATCH 09/46] Current Density - Changed from J to j for size.J vs size.j --- motor_weight_optimization.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 64ab6c5..9c495af 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -78,7 +78,7 @@ def compute(self,inputs,outputs): 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) # TODO: Better name for current density??? - +''' # TODO: Get this partial working: # Use: check_partials function to check: def compute_partials(self, inputs, J): @@ -151,7 +151,7 @@ def compute_partials(self, inputs, J): J['j', ''] = J['j', ''] = J['j', ''] = - +''' class torque(ExplicitComponent): def setup(self): @@ -324,7 +324,7 @@ def compute(self,inputs,outputs): model.add_subsystem('torque', torque(), promotes_inputs=['rot_or','b_g','i','n_m','n','l_st'], promotes_outputs=['tq']) model.connect('balance.rot_or', 'rot_or') - model.connect('size.J', 'balance.lhs:rot_or') + model.connect('size.j', 'balance.lhs:rot_or') model.connect('balance.l_st', 'l_st') model.connect('tq', 'balance.lhs:l_st') @@ -361,7 +361,7 @@ 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')) # NOTE: Changed "J" for current density to "j" print('Stack Length......................', p.get_val('mass.l_st', units='mm')) from solid import * From fba712a581d35fcb6cbc1646d241121835c2f0eb Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Wed, 19 Jun 2019 11:09:25 -0400 Subject: [PATCH 10/46] uncommented stuff --- motor_weight_optimization.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 9c495af..bc84919 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -78,7 +78,7 @@ def compute(self,inputs,outputs): 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) # TODO: Better name for current density??? -''' + # TODO: Get this partial working: # Use: check_partials function to check: def compute_partials(self, inputs, J): @@ -151,7 +151,7 @@ def compute_partials(self, inputs, J): J['j', ''] = J['j', ''] = J['j', ''] = -''' + class torque(ExplicitComponent): def setup(self): From fcac88bd32bfce8858f180fa469fb1a8df296cef Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 19 Jun 2019 16:20:21 -0400 Subject: [PATCH 11/46] added current density derivatives --- motor_weight_optimization.py | 145 +++++++++++++++++++---------------- 1 file changed, 81 insertions(+), 64 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 0ea3e21..17e5923 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -125,12 +125,12 @@ def compute_partials(self, inputs, J): # rotor_inner_radius t_mag = inputs['t_mag'] - J['rot_ir', 'rot_or'] = - J['rot_ir', 't_mag'] = - J['rot_ir', 'b_g'] = - J['rot_ir', 'n_m'] = - J['rot_ir', 'k'] = - J['rot_ir', 'b_ry'] = + 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 @@ -141,16 +141,33 @@ def compute_partials(self, inputs, J): i = inputs['i'] k_wb = inputs['k_wb'] n_s = inputs['n_s'] - J['j', ''] = - J['j', ''] = - J['j', ''] = - J['j', ''] = - J['j', ''] = - J['j', ''] = - J['j', ''] = - J['j', ''] = - J['j', ''] = - J['j', ''] = + 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*(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): @@ -175,20 +192,20 @@ 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*.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 + 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): @@ -242,44 +259,44 @@ def compute(self,inputs,outputs): # - # 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'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) + 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'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) - # rotor - # rot_ir=inputs['rot_ir'] - # rot_or=inputs['rot_or'] - # l_st=inputs['l_st'] - # t_mag=inputs['t_mag'] + 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 + 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'] + 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 + 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 From d3d6ab5729581044919d4a60bd7074e8d609561d Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 19 Jun 2019 16:40:18 -0400 Subject: [PATCH 12/46] issue with setup --- motor_weight_optimization.py | 281 +++++++++++++++++------------------ 1 file changed, 140 insertions(+), 141 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 17e5923..a4b957b 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -81,93 +81,93 @@ def compute(self,inputs,outputs): # TODO: Get this partial working: # Use: check_partials function to check: - 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*(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'] + # 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*(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): @@ -192,20 +192,20 @@ 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*.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 + # 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): @@ -259,44 +259,44 @@ def compute(self,inputs,outputs): # - 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'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) - - 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 + # 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'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) + # + # #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 @@ -351,10 +351,9 @@ def compute_partials(self,inputs,J): model.nonlinear_solver = NewtonSolver() model.nonlinear_solver.options['maxiter'] = 100 model.nonlinear_solver.options['iprint'] = 0 - p.setup() p.final_setup() - #p.check_partials(compact_print=True) + p.check_partials(compact_print=True) p.model.list_outputs(implicit=True) p.set_solver_print(2) #view_model(p) From 7210bc9b09997567c036bf904f9bf410aed51ff1 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 19 Jun 2019 16:46:48 -0400 Subject: [PATCH 13/46] fixed j issue --- motor_weight_optimization.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index a4b957b..061e49f 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -341,7 +341,7 @@ def compute(self,inputs,outputs): model.add_subsystem('torque', torque(), promotes_inputs=['rot_or','b_g','i','n_m','n','l_st'], promotes_outputs=['tq']) model.connect('balance.rot_or', 'rot_or') - model.connect('size.J', 'balance.lhs:rot_or') + model.connect('size.j', 'balance.lhs:rot_or') model.connect('balance.l_st', 'l_st') model.connect('tq', 'balance.lhs:l_st') @@ -377,7 +377,7 @@ 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')) from solid import * From 9a305f5cec9892c320c8ae5d5464cce687a701a1 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Thu, 20 Jun 2019 09:56:13 -0400 Subject: [PATCH 14/46] working analytic derivatives --- motor_weight_optimization.py | 284 +++++++++++++++++------------------ 1 file changed, 142 insertions(+), 142 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 061e49f..41b4066 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -43,7 +43,7 @@ def setup(self): 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.declare_partials('*','*', method='fd') + self.declare_partials('*','*')#, method='fd') def compute(self,inputs,outputs): # rotor_outer_radius @@ -81,93 +81,93 @@ def compute(self,inputs,outputs): # TODO: Get this partial working: # Use: check_partials function to check: - # 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*(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'] + 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): @@ -180,7 +180,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'] @@ -192,20 +192,20 @@ 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*.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 + 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): @@ -230,7 +230,7 @@ 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 @@ -259,44 +259,44 @@ def compute(self,inputs,outputs): # - # 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'] = 2 * rho * l_st * -(pi * (sta_ir+s_d)) - # - # #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 + 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 From bf593ebf804f6dfcdce660ea84ad4b78195bc826 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Thu, 20 Jun 2019 13:19:18 -0400 Subject: [PATCH 15/46] last commit for derivs --- motor_weight_optimization.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 41b4066..60e4205 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -353,10 +353,10 @@ def compute_partials(self,inputs,J): model.nonlinear_solver.options['iprint'] = 0 p.setup() p.final_setup() - p.check_partials(compact_print=True) + #p.check_partials(compact_print=True) p.model.list_outputs(implicit=True) p.set_solver_print(2) - #view_model(p) + view_model(p) p.run_model() print('Rotor Outer Radius................', p.get_val('balance.rot_or', units='mm')) From 21848f85947142599a6a47b0532f3c244967af8f Mon Sep 17 00:00:00 2001 From: garobed1 Date: Thu, 20 Jun 2019 16:27:13 -0400 Subject: [PATCH 16/46] added efficiency class --- motor_weight_optimization.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 60e4205..20219c3 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -4,8 +4,32 @@ from openmdao.api import Problem, IndepVarComp, ExplicitComponent, ExecComp from openmdao.api import NewtonSolver, Group, DirectSolver, NonlinearRunOnce, LinearRunOnce, view_model, BalanceComp -class motor_size(ExplicitComponent): +class efficiency(ExplicitComponent): + + def setup(self): + 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', 3800, units='rpm', desc='motor speed') + + self.add_output('P_in', units='W', desc='input power') + self.add_output('P_out', units='W', desc='output power') + self.add_output('nu', desc='efficiency') + + self.declare_partials('*','*',method='fd') + def compute(self, inputs, outputs): + i = inputs['i'] + tq = inputs['tq'] + v = inputs['v'] + rm = inputs['rm'] + + outputs['P_in'] = v*i + outputs['P_out'] = tq*rm*(2*pi/60) + outputs['nu'] = outputs['P_out']/outputs['P_in'] + +class motor_size(ExplicitComponent): + def setup(self): # rotor_outer_radius self.add_input('mot_or', 0.0765, units='m', desc='motor outer radius') From 8b93f4775539ed76583ff1f64fd16a261b20db87 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Fri, 21 Jun 2019 11:37:09 -0400 Subject: [PATCH 17/46] some optimization --- motor_weight_optimization.py | 93 +++++++++++++++++++++++++++--------- 1 file changed, 70 insertions(+), 23 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 20219c3..0a9b6ac 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -2,21 +2,21 @@ 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 efficiency(ExplicitComponent): def setup(self): 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', 3800, units='rpm', desc='motor speed') + self.add_input('v', 500, units='V', desc='RMS voltage') + self.add_input('rm', 5400, units='rpm', desc='motor speed') self.add_output('P_in', units='W', desc='input power') self.add_output('P_out', units='W', desc='output power') self.add_output('nu', desc='efficiency') - self.declare_partials('*','*',method='fd') + self.declare_partials('*','*') def compute(self, inputs, outputs): i = inputs['i'] @@ -24,12 +24,34 @@ def compute(self, inputs, outputs): v = inputs['v'] rm = inputs['rm'] - outputs['P_in'] = v*i + outputs['P_in'] = (3**.5)*v*i outputs['P_out'] = tq*rm*(2*pi/60) outputs['nu'] = outputs['P_out']/outputs['P_in'] + + 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('mot_or', 0.0765, units='m', desc='motor outer radius') @@ -333,7 +355,7 @@ def compute_partials(self,inputs,J): #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('mot_or', val=0.08, units='m') # Motor outer radius @@ -351,39 +373,62 @@ def compute_partials(self,inputs,J): 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') - 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','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','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(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('efficiency', efficiency(), promotes_inputs=['tq','i'], promotes_outputs=['nu']) - model.connect('balance.rot_or', 'rot_or') - model.connect('size.j', 'balance.lhs:rot_or') + # 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.connect('balance.l_st', 'l_st') - model.connect('tq', 'balance.lhs:l_st') - - model.linear_solver = DirectSolver() + #model.linear_solver = DirectSolver() model.nonlinear_solver = NewtonSolver() model.nonlinear_solver.options['maxiter'] = 100 model.nonlinear_solver.options['iprint'] = 0 + + #optimization setup + model.driver = ScipyOptimizeDriver() + model.driver.options['optimizer'] = 'SLSQP' + model.add_objective('nu', ref=-1) + model.add_design_var('gap', lower=.001, upper=.003) + model.add_design_var('mot_or', lower=.07, upper=.10) + model.add_design_var('rot_or', lower = .05, upper=.068) + model.add_design_var('l_st', lower = .0004, upper = .0006) + model.add_constraint('tq', lower=24, scaler=1) + model.add_constraint('size.j', upper=13, scaler=1) + p.setup() + + p['gap'] = 0.001 + p['mot_or'] = 0.09 + p['rot_or'] = 0.06 + p['l_st'] = 0.01 + p.final_setup() #p.check_partials(compact_print=True) p.model.list_outputs(implicit=True) - p.set_solver_print(2) - view_model(p) - p.run_model() - - print('Rotor Outer Radius................', p.get_val('balance.rot_or', units='mm')) + p.set_solver_print() + #view_model(p) + #p.run_model() + p.run_driver() + + print('Rotor Outer Radius................', p.get_val('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')) @@ -403,6 +448,8 @@ def compute_partials(self,inputs,J): print('Mass of Magnets...................', p.get_val('mag_mass', units='kg')) print('Current Density...................', p.get_val('size.j')) print('Stack Length......................', p.get_val('mass.l_st', units='mm')) + + print('Efficiency........................', p.get_val('nu')) from solid import * from solid.utils import * @@ -420,7 +467,7 @@ def compute_partials(self,inputs,J): 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=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) From d9a5b29ea2b0b6d6723c9aceed498c4ea3cbd366 Mon Sep 17 00:00:00 2001 From: Aaron Grgurich Date: Mon, 24 Jun 2019 11:07:16 -0400 Subject: [PATCH 18/46] Model - Problem Fix - Model Driver runs now --- motor_weight_optimization.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 0a9b6ac..72a0f70 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -123,6 +123,7 @@ def compute(self,inputs,outputs): outputs['sta_ir'] = rot_or + gap 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']) # TODO: Better name for current density??? # TODO: Get this partial working: @@ -398,24 +399,25 @@ def compute_partials(self,inputs,J): #model.linear_solver = DirectSolver() - model.nonlinear_solver = NewtonSolver() - model.nonlinear_solver.options['maxiter'] = 100 - model.nonlinear_solver.options['iprint'] = 0 + #model.nonlinear_solver = NewtonSolver() + #model.nonlinear_solver.options['maxiter'] = 100 + #model.nonlinear_solver.options['iprint'] = 2 #optimization setup - model.driver = ScipyOptimizeDriver() - model.driver.options['optimizer'] = 'SLSQP' + p.driver = ScipyOptimizeDriver() + p.driver.options['optimizer'] = 'SLSQP' + p.driver.options['debug_print'] = ['desvars', 'objs'] model.add_objective('nu', ref=-1) model.add_design_var('gap', lower=.001, upper=.003) model.add_design_var('mot_or', lower=.07, upper=.10) model.add_design_var('rot_or', lower = .05, upper=.068) - model.add_design_var('l_st', lower = .0004, upper = .0006) - model.add_constraint('tq', lower=24, scaler=1) - model.add_constraint('size.j', upper=13, scaler=1) + model.add_design_var('l_st', lower = .0004, upper = .003) + model.add_constraint('tq', lower=20, scaler=1) + model.add_constraint('size.j', upper=15, scaler=1) p.setup() - p['gap'] = 0.001 + p['gap'] = 0.002 p['mot_or'] = 0.09 p['rot_or'] = 0.06 p['l_st'] = 0.01 @@ -423,7 +425,7 @@ def compute_partials(self,inputs,J): p.final_setup() #p.check_partials(compact_print=True) p.model.list_outputs(implicit=True) - p.set_solver_print() + p.set_solver_print(level=2) #view_model(p) #p.run_model() p.run_driver() From 88c51c5f411e87ec3f6c84a14a8d59ec64a55756 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 24 Jun 2019 11:13:07 -0400 Subject: [PATCH 19/46] overwrite --- motor_weight_optimization.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 0a9b6ac..8f5a094 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -396,41 +396,42 @@ def compute_partials(self,inputs,J): # model.connect('balance.l_st', 'l_st') # model.connect('tq', 'balance.lhs:l_st') - #model.linear_solver = DirectSolver() - + model.linear_solver = DirectSolver() + model.nonlinear_solver = NewtonSolver() model.nonlinear_solver.options['maxiter'] = 100 - model.nonlinear_solver.options['iprint'] = 0 + model.nonlinear_solver.options['iprint'] = 2 #optimization setup model.driver = ScipyOptimizeDriver() model.driver.options['optimizer'] = 'SLSQP' + model.driver.options['debug_print'] = ['desvars', 'objs'] model.add_objective('nu', ref=-1) model.add_design_var('gap', lower=.001, upper=.003) model.add_design_var('mot_or', lower=.07, upper=.10) model.add_design_var('rot_or', lower = .05, upper=.068) - model.add_design_var('l_st', lower = .0004, upper = .0006) - model.add_constraint('tq', lower=24, scaler=1) - model.add_constraint('size.j', upper=13, scaler=1) + model.add_design_var('l_st', lower = .0004, upper = .0030) + model.add_constraint('tq', lower=14, scaler=1) + model.add_constraint('size.j', upper=15, scaler=1) p.setup() - p['gap'] = 0.001 - p['mot_or'] = 0.09 + p['gap'] = 0.002 + p['mot_or'] = 0.08 p['rot_or'] = 0.06 - p['l_st'] = 0.01 + p['l_st'] = 0.0009 p.final_setup() #p.check_partials(compact_print=True) p.model.list_outputs(implicit=True) - p.set_solver_print() + p.set_solver_print(level=2) #view_model(p) #p.run_model() p.run_driver() - print('Rotor Outer Radius................', p.get_val('rot_or', units='mm')) + # 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('Stator Inner Radius...............', p.get_val('sta_ir', units='mm')) print('Motor Outer Radius................', p.get_val('mass.mot_or', units='mm')) From 16772c1d757a67567f978cd95d540512ee556010 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 24 Jun 2019 14:18:50 -0400 Subject: [PATCH 20/46] some opt work --- motor_weight_optimization.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index f31e96c..4f16162 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -9,7 +9,7 @@ class efficiency(ExplicitComponent): def setup(self): self.add_input('i', 30, units='A', desc='RMS current') self.add_input('tq', 25, units='N*m', desc='torque') - self.add_input('v', 500, units='V', desc='RMS voltage') + self.add_input('v', 385, units='V', desc='RMS voltage') self.add_input('rm', 5400, units='rpm', desc='motor speed') self.add_output('P_in', units='W', desc='input power') @@ -406,13 +406,13 @@ def compute_partials(self,inputs,J): #optimization setup p.driver = ScipyOptimizeDriver() p.driver.options['optimizer'] = 'SLSQP' - p.driver.options['debug_print'] = ['desvars', 'objs'] + p.driver.options['disp'] = True model.add_objective('nu', ref=-1) - model.add_design_var('gap', lower=.001, upper=.003) - model.add_design_var('mot_or', lower=.07, upper=.10) - model.add_design_var('rot_or', lower = .05, upper=.068) + 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 = .0004, upper = .003) - model.add_constraint('tq', lower=20, scaler=1) + model.add_constraint('tq', lower=12, scaler=1) model.add_constraint('size.j', upper=15, scaler=1) p.setup() @@ -421,7 +421,7 @@ def compute_partials(self,inputs,J): p['mot_or'] = 0.08 p['rot_or'] = 0.06 p['l_st'] = 0.0009 - + # p.final_setup() #p.check_partials(compact_print=True) p.model.list_outputs(implicit=True) @@ -433,7 +433,7 @@ def compute_partials(self,inputs,J): # 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('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.mot_or', units='mm')) @@ -450,6 +450,8 @@ def compute_partials(self,inputs,J): print('Mass of Magnets...................', p.get_val('mag_mass', units='kg')) print('Current Density...................', p.get_val('size.j')) print('Stack Length......................', p.get_val('mass.l_st', units='mm')) + print('Power In..........................', p.get_val('efficiency.P_in', units='W')) + print('Power Out..........................', p.get_val('efficiency.P_out', units='W')) print('Efficiency........................', p.get_val('nu')) From e424e316a260716091ce7ae0c4d4555142045b1d Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 24 Jun 2019 16:13:40 -0400 Subject: [PATCH 21/46] last working power efficiency calc --- .gitignore | 1 + motor_weight_optimization.py | 12 +++++------- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/.gitignore b/.gitignore index 6094a05..1c092bd 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ /n2.html +/old.py diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 4f16162..3c7094c 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -411,16 +411,16 @@ def compute_partials(self,inputs,J): 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 = .0004, upper = .003) - model.add_constraint('tq', lower=12, scaler=1) - model.add_constraint('size.j', upper=15, scaler=1) + 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.0009 + p['l_st'] = 0.008 # p.final_setup() #p.check_partials(compact_print=True) @@ -450,9 +450,7 @@ def compute_partials(self,inputs,J): print('Mass of Magnets...................', p.get_val('mag_mass', units='kg')) print('Current Density...................', p.get_val('size.j')) print('Stack Length......................', p.get_val('mass.l_st', units='mm')) - print('Power In..........................', p.get_val('efficiency.P_in', units='W')) - print('Power Out..........................', p.get_val('efficiency.P_out', units='W')) - + print('Motor Speed.......................', p.get_val('efficiency.rm', units='rpm')) print('Efficiency........................', p.get_val('nu')) from solid import * From 02a9a0cee759a01acf4c6ce433e3351fc4c9cf6a Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 1 Jul 2019 14:14:36 -0400 Subject: [PATCH 22/46] added temp loss models --- motor_weight_optimization.py | 186 ++++++++++++++++++++++++++++------- 1 file changed, 153 insertions(+), 33 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 3c7094c..c01ffef 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -4,49 +4,171 @@ from openmdao.api import Problem, IndepVarComp, ExplicitComponent, ExecComp from openmdao.api import NewtonSolver, Group, DirectSolver, NonlinearRunOnce, LinearRunOnce, view_model, BalanceComp, ScipyOptimizeDriver +class core_loss(ExplicitComponent): + + def setup(self): + self.add_input('rm', units='rpm', desc='motor speed') + self.add_input('alpha', 1.5, desc='core loss constant') + + self.add_output('L_core', units='W', desc='core loss') + + self.declare_partials('*','*',method=fd) + def compute(self, inputs, outputs): + rm = inputs['rm'] + alpha = inputs['alpha'] + + #need to translate motor speed to electric frequency + f = rm + 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_output('L_emag', units='W', desc='magnetic eddy loss') + + self.declare_partials('*','*',method=fd) + def compute(self, inputs, outputs): + rm = inputs['rm'] + + #need to translate motor speed to electric frequency + f = rm + 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_output('L_ewir', units='W', desc='winding eddy loss') + + self.declare_partials('*','*',method=fd) + def compute(self, inputs, outputs): + rm = inputs['rm'] + + #need to translate motor speed to electric frequency + f = rm + 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', 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', units='(m**2)/s', desc='air dynamic viscosity') + + #bearing parameters + self.add_input('muf_b', desc='bearing friction coefficient') + self.add_input('D_b', units='m', desc='bearing bore diameter') + self.add_input('F_b', 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['rm'] + F_b = inputs['F_b'] + + #air gap friction coefficient, Re number + Reg = rho_a*rm*rot_or*gap/mu_a + Cfg = .515*((gap**.3)/rot_or)/(Reg**.5) + + outputs['L_airg'] = k*Cfa*pi*rho_a*(rm**3)*(rot_or**4)*l_st + + #axial air friction coefficient, Re number + Rea = rho_a*rm*(rot_or**2)/mu_a + Cfa = .146/(Rea**2) + + outputs['L_airf'] = .5*Cfa*rho_a*(rm**3)*((rot_or**5)-(rot_ir**5)) + + outputs['L_bear'] = .5*muf_b*D_b*F_b*rm + + #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_', 5400, units='rpm', desc='motor speed') + 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('nu', desc='efficiency') + self.add_output('L_total', units='W' desc='total loss'') - self.declare_partials('*','*') + self.declare_partials('*','*',method=fd) def compute(self, inputs, outputs): i = inputs['i'] tq = inputs['tq'] v = inputs['v'] rm = inputs['rm'] - - outputs['P_in'] = (3**.5)*v*i - outputs['P_out'] = tq*rm*(2*pi/60) - outputs['nu'] = outputs['P_out']/outputs['P_in'] - - 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 + + #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 @@ -124,10 +246,8 @@ def compute(self,inputs,outputs): 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']) - # TODO: Better name for current density??? + - # TODO: Get this partial working: - # Use: check_partials function to check: def compute_partials(self, inputs, J): # rotor_yoke_width @@ -389,7 +509,7 @@ def compute_partials(self,inputs,J): # 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('efficiency', efficiency(), promotes_inputs=['tq','i'], promotes_outputs=['nu']) + model.add_subsystem('efficiency', efficiency(), promotes_inputs=['L_airg','L_airf','L_bear','L_emag','L_ewir','L_core','L_res'], promotes_outputs=['L_total']) # model.connect('balance.rot_or', 'rot_or') # model.connect('size.j', 'balance.lhs:rot_or') @@ -407,7 +527,7 @@ def compute_partials(self,inputs,J): p.driver = ScipyOptimizeDriver() p.driver.options['optimizer'] = 'SLSQP' p.driver.options['disp'] = True - model.add_objective('nu', ref=-1) + 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) @@ -451,7 +571,7 @@ def compute_partials(self,inputs,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('efficiency.rm', units='rpm')) - print('Efficiency........................', p.get_val('nu')) + print('Total Loss........................', p.get_val('L_total')) from solid import * from solid.utils import * From 91029b210718ed5f03e7c9c0c1572f4cea4cc6a1 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 2 Jul 2019 12:55:20 -0400 Subject: [PATCH 23/46] working surrogate models, optimization not useful --- motor_weight_optimization.py | 98 +++++++++++++++++++++++------------- 1 file changed, 64 insertions(+), 34 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index c01ffef..1ea47a5 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -8,47 +8,50 @@ class core_loss(ExplicitComponent): def setup(self): self.add_input('rm', units='rpm', desc='motor speed') - self.add_input('alpha', 1.5, desc='core loss constant') + 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) + self.declare_partials('*','*',method='fd') def compute(self, inputs, outputs): rm = inputs['rm'] + n_m = inputs['n_m'] alpha = inputs['alpha'] - #need to translate motor speed to electric frequency - f = rm + 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) + self.declare_partials('*','*',method='fd') def compute(self, inputs, outputs): rm = inputs['rm'] + n_m = inputs['n_m'] - #need to translate motor speed to electric frequency - f = rm + 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) + self.declare_partials('*','*',method='fd') def compute(self, inputs, outputs): rm = inputs['rm'] + n_m = inputs['n_m'] - #need to translate motor speed to electric frequency - f = rm + f = n_m*rm/120 outputs['L_ewir'] = .00040681*(f**2) class mech_loss(ExplicitComponent): @@ -56,18 +59,18 @@ class mech_loss(ExplicitComponent): def setup(self): #windage parameters self.add_input('rm', 5400, units='rpm', desc='motor speed') - self.add_input('k', desc='windage constant k') + 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', units='(m**2)/s', desc='air dynamic viscosity') + self.add_input('mu_a', 1.81e-5, units='(m**2)/s', desc='air dynamic viscosity') #bearing parameters - self.add_input('muf_b', desc='bearing friction coefficient') - self.add_input('D_b', units='m', desc='bearing bore diameter') - self.add_input('F_b', units='N', desc='bearing load') #coupled with motor mass + 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 @@ -87,22 +90,23 @@ def compute(self, inputs, outputs): mu_a = inputs['mu_a'] muf_b = inputs['muf_b'] - D_b = inputs['rm'] + D_b = inputs['D_b'] F_b = inputs['F_b'] + omega = rm*(2*pi/60) #air gap friction coefficient, Re number - Reg = rho_a*rm*rot_or*gap/mu_a + Reg = rho_a*omega*rot_or*gap/mu_a Cfg = .515*((gap**.3)/rot_or)/(Reg**.5) - outputs['L_airg'] = k*Cfa*pi*rho_a*(rm**3)*(rot_or**4)*l_st + outputs['L_airg'] = k*Cfg*pi*rho_a*(omega**3)*(rot_or**4)*l_st #axial air friction coefficient, Re number - Rea = rho_a*rm*(rot_or**2)/mu_a + Rea = rho_a*omega*(rot_or**2)/mu_a Cfa = .146/(Rea**2) - outputs['L_airf'] = .5*Cfa*rho_a*(rm**3)*((rot_or**5)-(rot_ir**5)) + 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*rm + outputs['L_bear'] = .5*muf_b*D_b*F_b*omega #def compute_partials(self, inputs, J): @@ -122,8 +126,6 @@ def setup(self): 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_', 5400, units='rpm', desc='motor speed') - 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') @@ -134,18 +136,26 @@ def setup(self): 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.add_output('P_out', units='W', desc='output power') + self.add_output('L_total', units='W', desc='total loss') - self.declare_partials('*','*',method=fd) + 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_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 @@ -496,6 +506,8 @@ def compute_partials(self,inputs,J): 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 = BalanceComp() # @@ -509,7 +521,17 @@ def compute_partials(self,inputs,J): # 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('efficiency', efficiency(), promotes_inputs=['L_airg','L_airf','L_bear','L_emag','L_ewir','L_core','L_res'], promotes_outputs=['L_total']) + + 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') @@ -528,12 +550,12 @@ def compute_partials(self,inputs,J): 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('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) + model.add_constraint('tq', lower = 24, scaler = 1) + model.add_constraint('size.j', upper = 13, scaler = 1) p.setup() @@ -570,7 +592,15 @@ def compute_partials(self,inputs,J): print('Mass of Magnets...................', p.get_val('mag_mass', units='kg')) 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('efficiency.rm', units='rpm')) + 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 ef186217a11d63b67c0b8d2887a16e2a0c1f9965 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 2 Jul 2019 14:20:49 -0400 Subject: [PATCH 24/46] added notes tex file --- .gitignore | 3 ++ notes.tex | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 notes.tex diff --git a/.gitignore b/.gitignore index 1c092bd..04d6b82 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,5 @@ /n2.html /old.py +/notes.pdf +/notes.log +/notes.aux diff --git a/notes.tex b/notes.tex new file mode 100644 index 0000000..f04a6b9 --- /dev/null +++ b/notes.tex @@ -0,0 +1,90 @@ +\documentclass[11pt]{article} +\usepackage{graphicx} +\usepackage{parskip} +\usepackage{amsmath} +\usepackage{float} +\usepackage{listings} +\usepackage{color} + +\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{black}, % 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} + +\maketitle + +\section{NEXT STEP} + +Talk with Aaron about MotorCAD manual, talk about FEM with Justin, use total mass as objective + +\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 + +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} + + + +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. + +\section{Specific Notes} + + + +\end{document} From 275a99a8eb82e78815c7ca552cd899ba67608a15 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 2 Jul 2019 14:47:42 -0400 Subject: [PATCH 25/46] fixed aarons name --- notes.tex | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/notes.tex b/notes.tex index f04a6b9..f777a6a 100644 --- a/notes.tex +++ b/notes.tex @@ -1,17 +1,18 @@ -\documentclass[11pt]{article} +\documentclass[10pt]{article} \usepackage{graphicx} \usepackage{parskip} \usepackage{amsmath} \usepackage{float} \usepackage{listings} \usepackage{color} +\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{black}, % choose the background color + 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 @@ -24,7 +25,7 @@ \begin{document} \title{Progress Notes on Motor Design Code} -\author{Garo Bedonian} +\author{Garo Bedonian, Aaron Grgurich} \maketitle @@ -79,7 +80,7 @@ \section{High-Level Notes} \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. From 5637900507c6e02040ac348f999a7bb2d2eb87ec Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 3 Jul 2019 14:10:32 -0400 Subject: [PATCH 26/46] some notes on core loss --- motor_weight_optimization.py | 4 ++-- notes.tex | 20 ++++++++++++++++++-- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/motor_weight_optimization.py b/motor_weight_optimization.py index 1ea47a5..99eee3d 100644 --- a/motor_weight_optimization.py +++ b/motor_weight_optimization.py @@ -569,8 +569,8 @@ def compute_partials(self,inputs,J): p.model.list_outputs(implicit=True) p.set_solver_print(level=2) #view_model(p) - #p.run_model() - p.run_driver() + p.run_model() + #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')) diff --git a/notes.tex b/notes.tex index f777a6a..9a396d9 100644 --- a/notes.tex +++ b/notes.tex @@ -5,6 +5,7 @@ \usepackage{float} \usepackage{listings} \usepackage{color} +\usepackage{ulem} \usepackage[margin = 1in]{geometry} \definecolor{mygreen}{rgb}{0,0.6,0} @@ -31,7 +32,9 @@ \section{NEXT STEP} -Talk with Aaron about MotorCAD manual, talk about FEM with Justin, use total mass as objective +\sout{Talk with Aaron about MotorCAD manual, talk about FEM with Justin,} use total mass as objective + +Use Tucker's resources to figure out what to do to model loss \section{High-Level Notes} @@ -51,8 +54,10 @@ \section{High-Level Notes} 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? +Build something in EMpy? Or MFEM? Haven't started on this yet @@ -86,6 +91,17 @@ \section{Notable Commits} \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 \end{document} From 1d1a76bc65d98644e7ea9426a8a2ae5f0f12f601 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Fri, 5 Jul 2019 15:26:52 -0400 Subject: [PATCH 27/46] notes on magnetic flux density equations --- notes.tex | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/notes.tex b/notes.tex index 9a396d9..9b6de78 100644 --- a/notes.tex +++ b/notes.tex @@ -34,7 +34,9 @@ \section{NEXT STEP} \sout{Talk with Aaron about MotorCAD manual, talk about FEM with Justin,} use total mass as objective -Use Tucker's resources to figure out what to do to model loss +\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 \section{High-Level Notes} @@ -103,5 +105,34 @@ \section{Specific Notes} https://www.amazon.com/Electric-Machinery-Fundamentals-Chapman/dp/007108617X +Hysteresis Losses: + +\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. \end{document} From 4b8cac2b705280af421e6ae5cb13460404d77aee Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 8 Jul 2019 16:29:24 -0400 Subject: [PATCH 28/46] added toothflux, still need to finish --- notes.tex | 4 +++- toothflux.py | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 1 deletion(-) create mode 100644 toothflux.py diff --git a/notes.tex b/notes.tex index 9b6de78..4c23626 100644 --- a/notes.tex +++ b/notes.tex @@ -38,6 +38,8 @@ \section{NEXT STEP} Hanselman book will be the best resource for this, use it to define loss and other necessary parameters. See specific notes +Current: Find $B_r$, $K_a$, and coefficients for correction factor + \section{High-Level Notes} Goal: Make strides in unifying motor design process using OpenMDAO. @@ -105,7 +107,7 @@ \section{Specific Notes} https://www.amazon.com/Electric-Machinery-Fundamentals-Chapman/dp/007108617X -Hysteresis Losses: +Hysteresis Losses (Steinmetz Equations): \begin{equation} P_h = k_hfB_{pk}^n diff --git a/toothflux.py b/toothflux.py new file mode 100644 index 0000000..a345c07 --- /dev/null +++ b/toothflux.py @@ -0,0 +1,52 @@ +''' +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 + + +def toothflux(theta, params, phi): + n_m = params['n_m'] #number of poles + n_s = params['n_s'] #number of Slots + l_st = params['l_st'] #axial motor length + sta_ir = params['sta_ir'] #stator inner radius + t_mag = params['t_mag'] #magnet thickness + rot_or = params['rot_or'] #rotor outer radius (INCLUDING MAGNETS) + gap = params['gap'] #air gaps + + nterms = 5 #flux fourier series terms + mterms = 5 #slot correction fourier series terms + + n_p = n_m/2 #number of pole pairs + + Kslm = np.zeros(mterms + 1) + Bgn = np.zeros(nterms + 1) + + phi = 0 + #flux term layer + for n in range(-nterms, nterms) + beta = n*n_p + + Ksl = 0 + #slot correction term layer + for m in range(-mterms, mterms) + + Kslm[m + mterms] = + + Ksl += Kslm[m + mterms]*sin(pi*(m + n*(n_m/(2*n_s))))/(pi*(m + n*(n_m/(2*n_s)))) + + + + Bgn[n+nterms] = -br*ka*(((sta_ir/rot_or)**(beta-1)) + ((sta_ir/rot_or)**(2*beta))*((rot_or/sta_ir)**(beta+1))) + + + phi += Bgn[n+nterms]*(2*pi*l_st*sta_ir/n_s)*Ksl*exp(1j*n*theta) + + + + + + From b1f7721b6a63ed35c1ecf08ef372249818696b3e Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 9 Jul 2019 11:44:22 -0400 Subject: [PATCH 29/46] need magnetization profile --- notes.tex | 8 +++++++- toothflux.py | 20 ++++++++++++++++---- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/notes.tex b/notes.tex index 4c23626..2d7ed18 100644 --- a/notes.tex +++ b/notes.tex @@ -38,7 +38,13 @@ \section{NEXT STEP} Hanselman book will be the best resource for this, use it to define loss and other necessary parameters. See specific notes -Current: Find $B_r$, $K_a$, and coefficients for correction factor +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. + +Very important: Need the magnetization profile of the halbach array to determine $K_{rn}$, $K_{\theta{}n}$ \section{High-Level Notes} diff --git a/toothflux.py b/toothflux.py index a345c07..884513c 100644 --- a/toothflux.py +++ b/toothflux.py @@ -12,11 +12,15 @@ def toothflux(theta, params, phi): n_m = params['n_m'] #number of poles n_s = params['n_s'] #number of Slots l_st = params['l_st'] #axial motor length - sta_ir = params['sta_ir'] #stator inner radius + r_s = params['sta_ir'] #stator inner radius t_mag = params['t_mag'] #magnet thickness - rot_or = params['rot_or'] #rotor outer radius (INCLUDING MAGNETS) + r_m = params['rot_or'] #rotor outer radius (INCLUDING MAGNETS) gap = params['gap'] #air gaps + r_r = r_m - gap + + b_r = 1.3 #magnetic remanence + nterms = 5 #flux fourier series terms mterms = 5 #slot correction fourier series terms @@ -25,6 +29,10 @@ def toothflux(theta, params, phi): Kslm = np.zeros(mterms + 1) Bgn = np.zeros(nterms + 1) + mu_r = .001 #magnet relative recoil permeability, aka relative permitivity, mu/mu_0 + #k_rn = 1 + #k_tn = 1 + phi = 0 #flux term layer for n in range(-nterms, nterms) @@ -38,12 +46,16 @@ def toothflux(theta, params, phi): Ksl += Kslm[m + mterms]*sin(pi*(m + n*(n_m/(2*n_s))))/(pi*(m + n*(n_m/(2*n_s)))) + 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))) + 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)) - Bgn[n+nterms] = -br*ka*(((sta_ir/rot_or)**(beta-1)) + ((sta_ir/rot_or)**(2*beta))*((rot_or/sta_ir)**(beta+1))) + 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))) - phi += Bgn[n+nterms]*(2*pi*l_st*sta_ir/n_s)*Ksl*exp(1j*n*theta) + phi += Bgn[n+nterms]*(2*pi*l_st*r_s/n_s)*Ksl*exp(1j*n*theta) From 9a8c942eb31142f939dd79080396b70d287ccc9c Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 10 Jul 2019 16:39:07 -0400 Subject: [PATCH 30/46] almost done with analysis --- notes.tex | 7 +++- toothflux.py | 93 +++++++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 88 insertions(+), 12 deletions(-) diff --git a/notes.tex b/notes.tex index 2d7ed18..77cbfdd 100644 --- a/notes.tex +++ b/notes.tex @@ -44,7 +44,10 @@ \section{NEXT STEP} The magnets are (for now) N48H neodymium. -Very important: Need the magnetization profile of the halbach array to determine $K_{rn}$, $K_{\theta{}n}$ +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) +Finish this up. \section{High-Level Notes} @@ -143,4 +146,6 @@ \section{Specific Notes} But I need to see if the magnetic field is sinusoidal. +Halbach array is said to have sinusoidal magnetic field profile.??? + \end{document} diff --git a/toothflux.py b/toothflux.py index 884513c..2fe1b7e 100644 --- a/toothflux.py +++ b/toothflux.py @@ -6,7 +6,7 @@ ''' import numpy as np from math import pi, sin, cos, exp - +import matplotlib.pyplot as plt def toothflux(theta, params, phi): n_m = params['n_m'] #number of poles @@ -16,8 +16,9 @@ def toothflux(theta, params, phi): t_mag = params['t_mag'] #magnet thickness r_m = params['rot_or'] #rotor outer radius (INCLUDING MAGNETS) gap = params['gap'] #air gaps + alpha = 1/30 #magnetic fraction, fraction of rotor circumference occupied by one magnet - r_r = r_m - gap + r_r = r_m - t_mag b_r = 1.3 #magnetic remanence @@ -25,27 +26,45 @@ def toothflux(theta, params, phi): mterms = 5 #slot correction fourier series terms n_p = n_m/2 #number of pole pairs + + ts = 2*pi/n_s #slot to slot angle + tt = #tooth width angle Kslm = np.zeros(mterms + 1) Bgn = np.zeros(nterms + 1) mu_r = .001 #magnet relative recoil permeability, aka relative permitivity, mu/mu_0 - #k_rn = 1 - #k_tn = 1 phi = 0 #flux term layer - for n in range(-nterms, nterms) + for n in range(-nterms, nterms): beta = n*n_p - Ksl = 0 + 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 coeffificients in Kslm + #compute using Riemann sum + #slot correction term layer - for m in range(-mterms, mterms) + for m in range(-mterms, mterms): - Kslm[m + mterms] = + #Kslm[m + mterms] = - Ksl += Kslm[m + mterms]*sin(pi*(m + n*(n_m/(2*n_s))))/(pi*(m + n*(n_m/(2*n_s)))) + slotsum += Kslm[m + mterms]*sin(pi*(m + n*(n_m/(2*n_s))))/(pi*(m + n*(n_m/(2*n_s)))) + # + # 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) + 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))) @@ -57,8 +76,60 @@ def toothflux(theta, params, phi): phi += Bgn[n+nterms]*(2*pi*l_st*r_s/n_s)*Ksl*exp(1j*n*theta) + +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): + if abs(theta) <= tt/2: + gratio = 1 + + if theta <= ts/2 and theta >= tt/2: + gratio = 1 + pi*r_s*(theta-tt/2)/(n_m*gap) + + if theta <= -tt/2 and theta >= -ts/2: + gratio = 1 - pi*r_s*(theta+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): + #func: one-dimensional function of interest + #T: period of function, will compute from center + #N: number of terms 1 and onward, this computes 2*N + 1 coefficients + + fsample = 2*N + t = np.linspace(-T/2, T/2, fsample + 2) + y = np.fft.rfft(func(t)) + + +if name == '__main__' + + gap = .003 + r_s = .07 + n_m = 20 + tt = pi/10 + ts = 1.1*pi/10 + mu_r = .001 + + period = ts - From b72a476b9a7d4048ed2a920129ce716cf2ae0c4b Mon Sep 17 00:00:00 2001 From: garobed1 Date: Thu, 11 Jul 2019 16:08:50 -0400 Subject: [PATCH 31/46] debug and test --- toothflux.py | 61 ++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 47 insertions(+), 14 deletions(-) diff --git a/toothflux.py b/toothflux.py index 2fe1b7e..240c6cb 100644 --- a/toothflux.py +++ b/toothflux.py @@ -16,6 +16,9 @@ def toothflux(theta, params, phi): t_mag = params['t_mag'] #magnet thickness r_m = params['rot_or'] #rotor outer radius (INCLUDING MAGNETS) gap = params['gap'] #air gaps + t_mag = params['t_mag'] #magnet thickness + tt = params['tt'] + ts = params['ts'] alpha = 1/30 #magnetic fraction, fraction of rotor circumference occupied by one magnet r_r = r_m - t_mag @@ -27,8 +30,7 @@ def toothflux(theta, params, phi): n_p = n_m/2 #number of pole pairs - ts = 2*pi/n_s #slot to slot angle - tt = #tooth width angle + Kslm = np.zeros(mterms + 1) Bgn = np.zeros(nterms + 1) @@ -46,11 +48,13 @@ def toothflux(theta, params, phi): #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) + #Ksl = kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, l_m) #fourier transform to obtain coeffificients 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 for m in range(-mterms, mterms): @@ -77,6 +81,8 @@ def toothflux(theta, params, phi): phi += Bgn[n+nterms]*(2*pi*l_st*r_s/n_s)*Ksl*exp(1j*n*theta) + #Bt = phi/(k_st*l_st*w_tb) + def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): gratio = gfunc(theta, gap, r_s, n_m, tt, ts) @@ -112,24 +118,51 @@ def RadialR(n, alpha): def RadialT(n, alpha): return 0 -def FourierCoeffsNumpy(func, T, N): +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 2*N + 1 coefficients + #N: number of terms 1 and onward, this computes N + 1 coefficients fsample = 2*N t = np.linspace(-T/2, T/2, fsample + 2) - y = np.fft.rfft(func(t)) + y = np.fft.rfft(func(t, gap, r_s, n_m, tt, ts, mu_r, t_mag))/t.size + c = y + for k in range(1, N): + np.insert(c, 0, np.conj(y(k))) + + return c -if name == '__main__' +if __name__ == '__main__': - gap = .003 + gap = .001 r_s = .07 n_m = 20 - tt = pi/10 - ts = 1.1*pi/10 - mu_r = .001 - - period = ts - + n_s = 24 + ts = 2*pi/n_s #slot to slot angle + tt = .6*pi/10 #tooth width angle + mu_r = 1 + t_mag = .003 + l_st = .008 + + 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 + } + + theta = ts/2 # np.linspace(-ts,ts) + + T = ts + N = 5 + + 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) + + print(x) From 07895c75a15a6de57cf516313abaa135cc981ac8 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Thu, 11 Jul 2019 16:33:45 -0400 Subject: [PATCH 32/46] small fixes --- toothflux.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/toothflux.py b/toothflux.py index 240c6cb..ba3849c 100644 --- a/toothflux.py +++ b/toothflux.py @@ -95,14 +95,16 @@ def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): def gfunc(theta, gap, r_s, n_m, tt, ts): - if abs(theta) <= tt/2: - gratio = 1 + gratio = np.zeros(len(theta)) + for x in range(0, len(theta)): + if abs(theta[x]) <= tt/2: + gratio[x] = 1 - if theta <= ts/2 and theta >= tt/2: - gratio = 1 + pi*r_s*(theta-tt/2)/(n_m*gap) + 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 <= -tt/2 and theta >= -ts/2: - gratio = 1 - pi*r_s*(theta+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 @@ -129,7 +131,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): c = y for k in range(1, N): - np.insert(c, 0, np.conj(y(k))) + np.insert(c, 0, np.conj(y[k])) return c @@ -157,7 +159,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): 'l_st':l_st } - theta = ts/2 # np.linspace(-ts,ts) + theta = np.linspace(-ts/2,ts/2) T = ts N = 5 @@ -166,3 +168,4 @@ def FourierCoeffsNumpy(func, T, N, 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) print(x) + print(y) From b6bca22bfe23cb8af1fddb6540838c374418ba45 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Fri, 12 Jul 2019 14:51:08 -0400 Subject: [PATCH 33/46] can get answer, need to debug --- toothflux.py | 92 ++++++++++++++++++++++++++++++++++------------------ 1 file changed, 61 insertions(+), 31 deletions(-) diff --git a/toothflux.py b/toothflux.py index ba3849c..ae9dbf4 100644 --- a/toothflux.py +++ b/toothflux.py @@ -8,10 +8,12 @@ from math import pi, sin, cos, exp import matplotlib.pyplot as plt -def toothflux(theta, params, phi): +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) @@ -25,21 +27,23 @@ def toothflux(theta, params, phi): b_r = 1.3 #magnetic remanence - nterms = 5 #flux fourier series terms - mterms = 5 #slot correction fourier series terms + #nterms = 10 #flux fourier series terms + mterms = 100 #slot correction fourier series terms n_p = n_m/2 #number of pole pairs - Kslm = np.zeros(mterms + 1) - Bgn = np.zeros(nterms + 1) + Kslm = np.zeros(2*mterms + 1) + Bgn = np.zeros(2*nterms + 1) + phi = np.zeros(2*nterms + 1) + Bn = np.zeros(2*nterms + 1) mu_r = .001 #magnet relative recoil permeability, aka relative permitivity, mu/mu_0 phi = 0 #flux term layer - for n in range(-nterms, nterms): + for n in range(-nterms, nterms+1): beta = n*n_p slotsum = 0 @@ -56,12 +60,14 @@ def toothflux(theta, params, phi): Kslm = FourierCoeffsNumpy(kfunc, T, mterms, gap, r_s, n_m, tt, ts, mu_r, t_mag) #slot correction term layer - for m in range(-mterms, mterms): - - #Kslm[m + mterms] = - - slotsum += Kslm[m + mterms]*sin(pi*(m + n*(n_m/(2*n_s))))/(pi*(m + n*(n_m/(2*n_s)))) - + 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 + # # add k_rn and k_tn here based on magnetic field orientation of halbach array # @@ -72,16 +78,19 @@ def toothflux(theta, params, phi): 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_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)) - + 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))) - - - phi += Bgn[n+nterms]*(2*pi*l_st*r_s/n_s)*Ksl*exp(1j*n*theta) - - - #Bt = phi/(k_st*l_st*w_tb) + + Bn[n+nterms] = Bgn[n+nterms]*(2*pi*l_st*r_s/n_s)*slotsum/(k_st*l_st*w_tb) + print(Bn[n+nterms]) + return Bn def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): @@ -126,13 +135,14 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): #N: number of terms 1 and onward, this computes N + 1 coefficients fsample = 2*N - t = np.linspace(-T/2, T/2, fsample + 2) - y = np.fft.rfft(func(t, gap, r_s, n_m, tt, ts, mu_r, t_mag))/t.size + t = np.linspace(-T/2, T/2, fsample + 1, 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): - np.insert(c, 0, np.conj(y[k])) + for k in range(1, N + 1): + c = np.insert(c, 0, np.conj(y[k])) + c = c.real return c if __name__ == '__main__': @@ -146,6 +156,9 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): mu_r = 1 t_mag = .003 l_st = .008 + k_st = 1 + w_tb = .003 + rot_or = .05 params = { 'gap':gap, @@ -156,16 +169,33 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): 'mu_r':mu_r, 't_mag':t_mag, 'sta_ir':r_s, - 'l_st':l_st + 'l_st':l_st, + 'k_st':k_st, + 'w_tb':w_tb, + 'rot_or':(r_s-gap) } - theta = np.linspace(-ts/2,ts/2) + theta = 0#np.linspace(-ts/2,ts/2) T = ts - N = 5 + N = 40 + + # 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) + + + Bn = toothflux(theta, N, params, Bn) + + #print(Bn) + + Brn = np.concatenate((Bn[40:], Bn[0:39])) + #print(Brn) + + i = np.fft.ifft(Brn) - 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) + #print(i) - print(x) - print(y) + # print(x) + # print(y) From d99dc2bc26e7888d0429a8dd1777d72b3e99aa6e Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 16 Jul 2019 16:57:26 -0400 Subject: [PATCH 34/46] little progress on conceptual issues --- notes.tex | 8 +++++--- toothflux.py | 46 ++++++++++++++++++++++++++-------------------- 2 files changed, 31 insertions(+), 23 deletions(-) diff --git a/notes.tex b/notes.tex index 77cbfdd..da5fb6b 100644 --- a/notes.tex +++ b/notes.tex @@ -44,10 +44,12 @@ \section{NEXT STEP} The magnets are (for now) N48H neodymium. -Very important: Need the magnetization profile of the halbach array to determine $K_{rn}$, $K_{\theta{}n}$. Will temporarily assume radial orientation. +\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) -Finish this up. +Very important: Need to compute fourier series of slot correction factor function (piecewise)} +Need to debug, values explode with higher truncation point + +Still an issue. Finish re-derivation, give it another shot for a day, then try an FEM approach. \section{High-Level Notes} diff --git a/toothflux.py b/toothflux.py index ae9dbf4..16c5121 100644 --- a/toothflux.py +++ b/toothflux.py @@ -17,10 +17,11 @@ def toothflux(theta, nterms, params, Bn): 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 gaps + gap = params['gap'] #air gap t_mag = params['t_mag'] #magnet thickness tt = params['tt'] ts = params['ts'] + mu_r = params['mu_r'] alpha = 1/30 #magnetic fraction, fraction of rotor circumference occupied by one magnet r_r = r_m - t_mag @@ -35,15 +36,15 @@ def toothflux(theta, nterms, params, Bn): Kslm = np.zeros(2*mterms + 1) - Bgn = np.zeros(2*nterms + 1) + Bgn = np.zeros(2*nterms + 1, dtype=np.complex_) phi = np.zeros(2*nterms + 1) - Bn = 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 + #mu_r = .001 #magnet relative recoil permeability, aka relative permitivity, mu/mu_0 phi = 0 #flux term layer - for n in range(-nterms, nterms+1): + for n in range(-nterms, nterms + 1): beta = n*n_p slotsum = 0 @@ -54,7 +55,7 @@ def toothflux(theta, nterms, params, Bn): #Ksl = kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, l_m) - #fourier transform to obtain coeffificients in Kslm + #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) @@ -73,6 +74,7 @@ def toothflux(theta, nterms, params, Bn): # ###for now, assume radial orientation k_rn = RadialR(n, alpha) + k_tn = RadialT(n, alpha) k_mc = (k_rn + 1j*beta*k_tn)/(1 - beta**2) @@ -80,17 +82,20 @@ def toothflux(theta, nterms, params, Bn): 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)) + 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) + #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(Bn[n+nterms]) - return Bn + print(Bgn[n+nterms]) + #print(Bgn) + return Bgn, Btn def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): @@ -148,12 +153,12 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): if __name__ == '__main__': gap = .001 - r_s = .07 + 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 + mu_r = 1.05 t_mag = .003 l_st = .008 k_st = 1 @@ -172,7 +177,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): 'l_st':l_st, 'k_st':k_st, 'w_tb':w_tb, - 'rot_or':(r_s-gap) + 'rot_or':(r_s - gap) } theta = 0#np.linspace(-ts/2,ts/2) @@ -184,18 +189,19 @@ def FourierCoeffsNumpy(func, T, N, 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 = toothflux(theta, N, params, Bn) + Bn, Bm = toothflux(theta, N, params, Bn) #print(Bn) - Brn = np.concatenate((Bn[40:], Bn[0:39])) + Brn = np.concatenate((Bn[N:], Bn[:N])) + #Btn = np.concatenate((Bm[10:], Bm[:10])) #print(Brn) - + #Brn = Bn[N:] i = np.fft.ifft(Brn) - #print(i) + print(i) # print(x) # print(y) From 4dc4a61acca9cd4aa870d79d9a6d55d16a6ccd40 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 17 Jul 2019 12:49:56 -0400 Subject: [PATCH 35/46] for aaron --- barntest.py | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++++ toothflux.py | 24 +++++++++++++-------- 2 files changed, 76 insertions(+), 9 deletions(-) create mode 100644 barntest.py diff --git a/barntest.py b/barntest.py new file mode 100644 index 0000000..abe351f --- /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*(((r_s/r_m)**(beta-1)) + ((r_s/r_m)**(2*beta))*((r_m/r_s)**(beta+1))) + print(Bgn) + i[n+N] = k_a + +plt.plot(i) +plt.show() diff --git a/toothflux.py b/toothflux.py index 16c5121..db383e5 100644 --- a/toothflux.py +++ b/toothflux.py @@ -22,7 +22,8 @@ def toothflux(theta, nterms, params, Bn): tt = params['tt'] ts = params['ts'] mu_r = params['mu_r'] - alpha = 1/30 #magnetic fraction, fraction of rotor circumference occupied by one magnet + alpha = params['alpha'] + #alpha = 18/30 #magnetic fraction, fraction of rotor circumference occupied by one magnet r_r = r_m - t_mag @@ -74,7 +75,7 @@ def toothflux(theta, nterms, params, Bn): # ###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) @@ -93,7 +94,7 @@ def toothflux(theta, nterms, params, Bn): 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[n+nterms]) #print(Bgn) return Bgn, Btn @@ -126,7 +127,7 @@ def gfunc(theta, gap, r_s, n_m, tt, ts): def RadialR(n, alpha): if n%2 != 0: - return alpha*(sin(n*alpha*pi/2))/(n*alpha*pi/2) + return alpha*sin(n*alpha*pi/2)/(n*alpha*pi/2) if n%2 == 0: return 0 @@ -164,6 +165,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): k_st = 1 w_tb = .003 rot_or = .05 + alpha = .85 params = { 'gap':gap, @@ -177,13 +179,14 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): 'l_st':l_st, 'k_st':k_st, 'w_tb':w_tb, - 'rot_or':(r_s - gap) + 'rot_or':(r_s - gap), + 'alpha':alpha } theta = 0#np.linspace(-ts/2,ts/2) T = ts - N = 40 + N = 200 # 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) @@ -195,13 +198,16 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): #print(Bn) - Brn = np.concatenate((Bn[N:], Bn[:N])) + #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.ifft(Brn) - + i = np.fft.irfft(Brn, N) + print(Brn) print(i) + plt.plot(i) + plt.show() # print(x) # print(y) From a5d15e91f5398584e82929df01d3c29e62fb891d Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 17 Jul 2019 16:44:34 -0400 Subject: [PATCH 36/46] bmatrix file presents a somewhat working solution --- barntest.py | 4 +-- bmatrix.py | 92 +++++++++++++++++++++++++++++++++++++++++++++++++++++ notes.tex | 3 +- 3 files changed, 96 insertions(+), 3 deletions(-) create mode 100644 bmatrix.py diff --git a/barntest.py b/barntest.py index abe351f..9973e09 100644 --- a/barntest.py +++ b/barntest.py @@ -53,9 +53,9 @@ def RadialT(n, alpha): #k_a = k_a.real #print(n) - Bgn = -b_r*k_a*(((r_s/r_m)**(beta-1)) + ((r_s/r_m)**(2*beta))*((r_m/r_s)**(beta+1))) + Bgn = -b_r*k_a*(((r_m/r_m)**(beta-1)) + ((r_s/r_m)**(2*beta))*((r_m/r_m)**(beta+1))) print(Bgn) - i[n+N] = k_a + i[n+N] = Bgn plt.plot(i) plt.show() diff --git a/bmatrix.py b/bmatrix.py new file mode 100644 index 0000000..de635fb --- /dev/null +++ b/bmatrix.py @@ -0,0 +1,92 @@ + +import numpy as np +import matplotlib.pyplot as plt + +# Da, Em, Dm +# Ea is precomputed +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 + +n_p = 10 +r_s = .051 +r_m = .05 +r_r = .045 +b_r = 1.3 +mu_r = 1.05 +mu_0 = 0.001 +alpha = .85 + +def B_fourier_terms(n): + bmat = np.zeros((3, 3)) + rhs = np.zeros(3) + + beta = n*n_p + + k_rn = RadialR(n, alpha) + k_tn = RadialT(n, alpha) + #Cm = b_r*(k_rn + 1j*beta*k_tn)/((1-beta**2)*mu_r*mu_0) + Cm = b_r*(k_rn)/((1-beta**2)*mu_r*mu_0) + + #Hatheta boundary condition + # bmat[0, 0] = (r_s**(-beta-1)) + # rhs[0] = -r_s**(beta-1) + Ea = -r_s**(2*beta) + + #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) + #print(rhs) + + bmat_inv = np.linalg.inv(bmat) + + #print(bmat_inv) + + x = bmat_inv.dot(rhs) + return x + +def B_arn(b_fourier, n, r): + Da = b_fourier[0] + #print(Da) + beta = n*n_p + Ea = -r_s**(2*beta) + b_arn = mu_0*beta*Da*(Ea*(r**(-beta-1)) - r**(beta-1)) + return b_arn + +if __name__ == '__main__': + + N = 22 + 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(B_r) + plt.show() + diff --git a/notes.tex b/notes.tex index da5fb6b..72968a9 100644 --- a/notes.tex +++ b/notes.tex @@ -49,7 +49,8 @@ \section{NEXT STEP} Very important: Need to compute fourier series of slot correction factor function (piecewise)} Need to debug, values explode with higher truncation point -Still an issue. Finish re-derivation, give it another shot for a day, then try an FEM approach. +\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. \section{High-Level Notes} From 6f73b38e4384d1680a21a4f7c8e3d789e6ff9e20 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Fri, 19 Jul 2019 16:45:36 -0400 Subject: [PATCH 37/46] need to finish slotcorrection, clean up code --- bmatrix.py | 27 +++++++++----- slotcorrection.py | 89 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 107 insertions(+), 9 deletions(-) create mode 100644 slotcorrection.py diff --git a/bmatrix.py b/bmatrix.py index de635fb..9153a73 100644 --- a/bmatrix.py +++ b/bmatrix.py @@ -4,6 +4,15 @@ # Da, Em, Dm # Ea is precomputed +def RadSinR(n, alpha): + if n == 1 or n == -1: + return .5 + else: + return 0 + +def RadSinT(n, alpha): + return 0 + def RadialR(n, alpha): if n%2 != 0: @@ -15,13 +24,13 @@ def RadialR(n, alpha): def RadialT(n, alpha): return 0 -n_p = 10 +n_p = 8 r_s = .051 r_m = .05 r_r = .045 b_r = 1.3 mu_r = 1.05 -mu_0 = 0.001 +mu_0 = 0.001 #value doesn't matter alpha = .85 def B_fourier_terms(n): @@ -38,7 +47,7 @@ def B_fourier_terms(n): #Hatheta boundary condition # bmat[0, 0] = (r_s**(-beta-1)) # rhs[0] = -r_s**(beta-1) - Ea = -r_s**(2*beta) + Ea = -(r_s)**(2*beta) #Hmtheta boundary condition bmat[0, 1] = (r_r**(-beta-1)) @@ -46,7 +55,7 @@ def B_fourier_terms(n): rhs[0] = -Cm #Htheta boundary condition - bmat[1, 0] = Ea*(r_m**(-beta-1)) + r_m**(beta - 1) + 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 @@ -68,16 +77,16 @@ def B_fourier_terms(n): return x def B_arn(b_fourier, n, r): - Da = b_fourier[0] + Da = -b_fourier[0] #print(Da) beta = n*n_p - Ea = -r_s**(2*beta) + Ea = -(r_s)**(2*beta) b_arn = mu_0*beta*Da*(Ea*(r**(-beta-1)) - r**(beta-1)) return b_arn if __name__ == '__main__': - N = 22 + N = 17 B_r = np.zeros(N) for n in range(1, N): @@ -86,7 +95,7 @@ def B_arn(b_fourier, n, r): print(B_r[n]) i = np.fft.irfft(B_r) - print(i) - plt.plot(B_r) + #print(i) + plt.plot(i) plt.show() diff --git a/slotcorrection.py b/slotcorrection.py new file mode 100644 index 0000000..07f0c99 --- /dev/null +++ b/slotcorrection.py @@ -0,0 +1,89 @@ + +import numpy as np +import matplotlib.pyplot as plt + + + +n_p = 8 +n_s = 24 +r_s = .051 +r_m = .05 +r_r = .045 +b_r = 1.3 +mu_r = 1.05 +mu_0 = 0.001 #value doesn't matter +alpha = .85 +n_m = 2*n_p +ts = 2*np.pi/n_s #slot to slot angle +tt = .6*np.pi/10 #tooth width angle +t_mag = r_m - r_r +gap = r_s - r_m + +def slot_correction(mterms): + + slotsum = 0 + T = ts + 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 + + return Kslm + +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 + 1, 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 + +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 + 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 = 20 + 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() From bb0d8cd2197609058820a3bd65c2122bc7b6e61d Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 23 Jul 2019 15:45:36 -0400 Subject: [PATCH 38/46] Working commit for computing tooth flux density --- bmatrix.py | 30 ++++++++++----- btooth.py | 94 +++++++++++++++++++++++++++++++++++++++++++++++ core_loss.py | 0 notes.tex | 6 +++ slotcorrection.py | 60 +++++++++++++++--------------- toothflux.py | 17 ++++----- 6 files changed, 160 insertions(+), 47 deletions(-) create mode 100644 btooth.py create mode 100644 core_loss.py diff --git a/bmatrix.py b/bmatrix.py index 9153a73..0071343 100644 --- a/bmatrix.py +++ b/bmatrix.py @@ -24,16 +24,25 @@ def RadialR(n, alpha): def RadialT(n, alpha): return 0 -n_p = 8 -r_s = .051 -r_m = .05 -r_r = .045 -b_r = 1.3 -mu_r = 1.05 +# n_p = 8 +# r_s = .051 +# r_m = .05 +# r_r = .045 +# b_r = 1.3 +# mu_r = 1.05 mu_0 = 0.001 #value doesn't matter -alpha = .85 +# alpha = .85 -def B_fourier_terms(n): + +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'] + b_r = params['b_r'] + mu_r = params['mu_r'] + alpha = params['alpha'] + bmat = np.zeros((3, 3)) rhs = np.zeros(3) @@ -76,7 +85,10 @@ def B_fourier_terms(n): x = bmat_inv.dot(rhs) return x -def B_arn(b_fourier, n, r): +def B_arn(b_fourier, n, params, r): + r_s = params['r_s'] + n_p = params['n_p'] + Da = -b_fourier[0] #print(Da) beta = n*n_p diff --git a/btooth.py b/btooth.py new file mode 100644 index 0000000..74f44fb --- /dev/null +++ b/btooth.py @@ -0,0 +1,94 @@ +''' +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 +import matplotlib.pyplot as plt +import bmatrix as bm +import slotcorrection as sl + +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'] + b_r = params['b_r'] + 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) + + K_slm = sl.slot_correction(M, params) + + for n in range(1, N): + + b_f = bm.B_fourier_terms(n, params) + B_gn[n] = bm.B_arn(b_f, n, params, r_s) + slotsum = 0 + #From Chapter 7 of Hanselman Book + 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 + #print(slotsum) + + B_tn[n] = B_gn[n]*(2*np.pi*r_s/(n_s*k_st*w_tb))*slotsum + print(B_gn) + 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 + + 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 + } + + B_tn = B_tooth(17, 20, params) + + i = np.fft.irfft(B_tn) + print(i) + plt.plot(i) + plt.show() + diff --git a/core_loss.py b/core_loss.py new file mode 100644 index 0000000..e69de29 diff --git a/notes.tex b/notes.tex index 72968a9..038a015 100644 --- a/notes.tex +++ b/notes.tex @@ -52,6 +52,10 @@ \section{NEXT STEP} \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. + \section{High-Level Notes} Goal: Make strides in unifying motor design process using OpenMDAO. @@ -151,4 +155,6 @@ \section{Specific Notes} 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. + \end{document} diff --git a/slotcorrection.py b/slotcorrection.py index 07f0c99..5282e55 100644 --- a/slotcorrection.py +++ b/slotcorrection.py @@ -4,47 +4,49 @@ -n_p = 8 -n_s = 24 -r_s = .051 -r_m = .05 -r_r = .045 -b_r = 1.3 -mu_r = 1.05 -mu_0 = 0.001 #value doesn't matter -alpha = .85 -n_m = 2*n_p -ts = 2*np.pi/n_s #slot to slot angle -tt = .6*np.pi/10 #tooth width angle -t_mag = r_m - r_r -gap = r_s - r_m +# n_p = 8 +# n_s = 24 +# r_s = .051 +# r_m = .05 +# r_r = .045 +# b_r = 1.3 +# mu_r = 1.05 +# mu_0 = 0.001 #value doesn't matter +# alpha = .85 +# n_m = 2*n_p +# ts = 2*np.pi/n_s #slot to slot angle +# tt = .6*np.pi/10 #tooth width angle +# t_mag = r_m - r_r +# gap = r_s - r_m -def slot_correction(mterms): - - slotsum = 0 +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'] + ts = params['ts'] + mu_r = params['mu_r'] + + gap = r_s - r_m + n_m = n_p*2 T = ts - Kslm = FourierCoeffsNumpy(kfunc, T, mterms, gap, r_s, n_m, tt, ts, mu_r, t_mag) + Kslm = FourierCoeffsNumpy(kfunc, T, M, 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 return Kslm -def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): +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 + 1, endpoint = False) + 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])) @@ -79,7 +81,7 @@ def gfunc(theta, gap, r_s, n_m, tt, ts): return gratio if __name__ == '__main__': - M = 20 + M = 30 K_slm = np.zeros(M) #for m in range(0, M): diff --git a/toothflux.py b/toothflux.py index db383e5..2040af3 100644 --- a/toothflux.py +++ b/toothflux.py @@ -30,7 +30,7 @@ def toothflux(theta, nterms, params, Bn): b_r = 1.3 #magnetic remanence #nterms = 10 #flux fourier series terms - mterms = 100 #slot correction fourier series terms + mterms = 20 #slot correction fourier series terms n_p = n_m/2 #number of pole pairs @@ -45,7 +45,7 @@ def toothflux(theta, nterms, params, Bn): phi = 0 #flux term layer - for n in range(-nterms, nterms + 1): + for n in range(1, nterms): beta = n*n_p slotsum = 0 @@ -69,7 +69,7 @@ def toothflux(theta, nterms, params, Bn): else: slotsum = 0 - + print(slotsum) # # add k_rn and k_tn here based on magnetic field orientation of halbach array # @@ -141,7 +141,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): #N: number of terms 1 and onward, this computes N + 1 coefficients fsample = 2*N - t = np.linspace(-T/2, T/2, fsample + 1, endpoint = False) + 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 @@ -155,7 +155,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): gap = .001 r_s = .051 - n_m = 20 + n_m = 16 n_s = 24 ts = 2*pi/n_s #slot to slot angle tt = .6*pi/10 #tooth width angle @@ -186,7 +186,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): theta = 0#np.linspace(-ts/2,ts/2) T = ts - N = 200 + 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) @@ -204,10 +204,9 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): #print(Brn) #Brn = Bn[N:] i = np.fft.irfft(Brn, N) - print(Brn) - print(i) + plt.plot(i) - plt.show() + #plt.show() # print(x) # print(y) From 6cc77fbc90a87122f1c4c148e93c7c36a84e8e79 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Tue, 23 Jul 2019 15:47:57 -0400 Subject: [PATCH 39/46] ignore files --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 04d6b82..d878c80 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ /notes.pdf /notes.log /notes.aux +/barntest* From 8ae0ee94b5e1104b3c6a080b34898dc2149cb803 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Wed, 24 Jul 2019 15:38:38 -0400 Subject: [PATCH 40/46] cleaned up, commented code --- barntest.py | 2 +- bmatrix.py | 57 +++++++++++++++++++++++++---------------------- btooth.py | 33 +++++++++++++++++---------- core_loss.py | 36 ++++++++++++++++++++++++++++++ notes.tex | 19 ++++++++++++++-- slotcorrection.py | 45 ++++++++++++++++--------------------- 6 files changed, 124 insertions(+), 68 deletions(-) diff --git a/barntest.py b/barntest.py index 9973e09..7d1d914 100644 --- a/barntest.py +++ b/barntest.py @@ -53,7 +53,7 @@ def RadialT(n, alpha): #k_a = k_a.real #print(n) - Bgn = -b_r*k_a*(((r_m/r_m)**(beta-1)) + ((r_s/r_m)**(2*beta))*((r_m/r_m)**(beta+1))) + 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 diff --git a/bmatrix.py b/bmatrix.py index 0071343..adf7921 100644 --- a/bmatrix.py +++ b/bmatrix.py @@ -1,9 +1,17 @@ +''' +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 -# Da, Em, Dm -# Ea is precomputed +# 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 @@ -13,6 +21,7 @@ def RadSinR(n, alpha): def RadSinT(n, alpha): return 0 +# Radial Profile def RadialR(n, alpha): if n%2 != 0: @@ -24,38 +33,35 @@ def RadialR(n, alpha): def RadialT(n, alpha): return 0 -# n_p = 8 -# r_s = .051 -# r_m = .05 -# r_r = .045 -# b_r = 1.3 -# mu_r = 1.05 -mu_0 = 0.001 #value doesn't matter -# alpha = .85 +# 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'] - b_r = params['b_r'] - mu_r = params['mu_r'] - alpha = params['alpha'] + 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) - #Cm = b_r*(k_rn + 1j*beta*k_tn)/((1-beta**2)*mu_r*mu_0) + + # 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 - # bmat[0, 0] = (r_s**(-beta-1)) - # rhs[0] = -r_s**(beta-1) + #Hatheta boundary condition (excluded from linear system) Ea = -(r_s)**(2*beta) #Hmtheta boundary condition @@ -75,25 +81,22 @@ def B_fourier_terms(n, params): bmat[2, 2] = beta*(r_m**(beta-1)) rhs[2] = (b_r*k_rn/(mu_r*mu_0)) - Cm - #print(bmat) - #print(rhs) - + #invert matrix bmat_inv = np.linalg.inv(bmat) - #print(bmat_inv) - + #solve system using inverse of system matrix x = bmat_inv.dot(rhs) 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'] - n_p = params['n_p'] + r_s = params['r_s'] #stator inner radius + n_p = params['n_p'] #number of pole pairs Da = -b_fourier[0] - #print(Da) beta = n*n_p Ea = -(r_s)**(2*beta) - b_arn = mu_0*beta*Da*(Ea*(r**(-beta-1)) - r**(beta-1)) + b_arn = mu_0*beta*Da*(Ea*(r**(-beta-1)) - r**(beta-1)) #Equation B.13 return b_arn if __name__ == '__main__': diff --git a/btooth.py b/btooth.py index 74f44fb..8a96e5f 100644 --- a/btooth.py +++ b/btooth.py @@ -1,8 +1,9 @@ ''' -This function computes magnetic flux in the stator teeth of a motor, described using a Fourier series with respect to angular position theta +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 + ''' @@ -10,6 +11,7 @@ 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 @@ -19,9 +21,9 @@ def B_tooth(N, M, params): 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'] - b_r = params['b_r'] - tt = params['tt'] + 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'] @@ -30,24 +32,25 @@ def B_tooth(N, M, params): 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 - #From Chapter 7 of Hanselman Book + # 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 - #print(slotsum) - - B_tn[n] = B_gn[n]*(2*np.pi*r_s/(n_s*k_st*w_tb))*slotsum - print(B_gn) + + # 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__': @@ -65,6 +68,7 @@ def B_tooth(N, M, params): w_tb = .005 alpha = .85 b_r = 1.3 + f = 5400/20 #rpm params = { 'gap':gap, @@ -82,13 +86,18 @@ def B_tooth(N, M, params): 'r_m':(r_s - gap), 'r_r':(r_s - gap - t_mag), 'alpha':alpha, - 'b_r':b_r + 'b_r':b_r, + 'f':f } B_tn = B_tooth(17, 20, params) 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 index e69de29..fc1e2ed 100644 --- a/core_loss.py +++ 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/notes.tex b/notes.tex index 038a015..ea6d4d0 100644 --- a/notes.tex +++ b/notes.tex @@ -44,9 +44,9 @@ \section{NEXT STEP} 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. +\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)} +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.} @@ -109,6 +109,8 @@ \section{Notable Commits} 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 @@ -157,4 +159,17 @@ \section{Specific Notes} 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. + + \end{document} diff --git a/slotcorrection.py b/slotcorrection.py index 5282e55..c1da8b8 100644 --- a/slotcorrection.py +++ b/slotcorrection.py @@ -1,24 +1,16 @@ +''' +Helps to compute the Fourier coefficients of the stator tooth flux density by computing a correction factor for the slots -import numpy as np -import matplotlib.pyplot as plt +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 +''' -# n_p = 8 -# n_s = 24 -# r_s = .051 -# r_m = .05 -# r_r = .045 -# b_r = 1.3 -# mu_r = 1.05 -# mu_0 = 0.001 #value doesn't matter -# alpha = .85 -# n_m = 2*n_p -# ts = 2*np.pi/n_s #slot to slot angle -# tt = .6*np.pi/10 #tooth width angle -# t_mag = r_m - r_r -# gap = r_s - r_m +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 @@ -26,22 +18,22 @@ def slot_correction(M, params): 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'] - ts = params['ts'] - mu_r = params['mu_r'] + 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 + - gap = r_s - r_m - n_m = n_p*2 - T = ts 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): - #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) @@ -54,6 +46,7 @@ def FourierCoeffsNumpy(func, T, N, gap, r_s, n_m, tt, ts, mu_r, t_mag): 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) @@ -63,7 +56,7 @@ def kfunc(theta, gap, r_s, n_m, tt, ts, mu_r, t_mag): return Ksl - +# Equation 7.11 def gfunc(theta, gap, r_s, n_m, tt, ts): gratio = np.zeros(len(theta)) From 952bad6f0677339b0fc6c2d72693570e68d56044 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 29 Jul 2019 16:44:25 -0400 Subject: [PATCH 41/46] progress on openmdao implementation --- notes.tex | 4 ++ radial_air_fourier_comp.py | 76 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+) create mode 100644 radial_air_fourier_comp.py diff --git a/notes.tex b/notes.tex index ea6d4d0..e8ab551 100644 --- a/notes.tex +++ b/notes.tex @@ -56,6 +56,8 @@ \section{NEXT STEP} Implement core loss equation. Consider options for determining coefficients. +Implement model in OpenMDAO + \section{High-Level Notes} Goal: Make strides in unifying motor design process using OpenMDAO. @@ -171,5 +173,7 @@ \section{Specific Notes} 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..6e72fed --- /dev/null +++ b/radial_air_fourier_comp.py @@ -0,0 +1,76 @@ +from __future__ import absolute_import +import numpy as np +from math import pi +import openmdao.api as om + +class DiffAssemble(om.ExplicitComponent()): + 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('x', desc='PDE solution coefficients') + + def compute(self, inputs, outputs): + 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) + + #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 + + #invert matrix + bmat_inv = np.linalg.inv(bmat) + + #solve system using inverse of system matrix + outputs['x'] = bmat_inv.dot(rhs) + +class MotorAirFourierComp(om.Group()): + + def setup(self): + self.add_subsystem('DiffAssemble', DiffAssemble()) + self.add_subsystem('DiffLinearSolve', DiffLinearSolve()) + self.add_subsystem('BarnComp', BarnComp()) + self.add_subsystem('CorrectSlot', CorrectSlot()) + self.add_subsystem() + + From 033f15b8827a487d923ac8c3eabf9ea34f0b78f1 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Fri, 2 Aug 2019 17:17:25 -0400 Subject: [PATCH 42/46] connection issue --- radial_air_fourier_comp.py | 363 ++++++++++++++++++++++++++++++++++--- 1 file changed, 333 insertions(+), 30 deletions(-) diff --git a/radial_air_fourier_comp.py b/radial_air_fourier_comp.py index 6e72fed..be70354 100644 --- a/radial_air_fourier_comp.py +++ b/radial_air_fourier_comp.py @@ -3,7 +3,15 @@ from math import pi import openmdao.api as om -class DiffAssemble(om.ExplicitComponent()): +N = 4 +M = 20 + +B_tna = np.zeros(N) + +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') @@ -12,65 +20,360 @@ def setup(self): 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('x', desc='PDE solution coefficients') - + + self.add_output('bmat', shape=(3,3), desc='system matrix') + self.add_output('rhs', shape=3, desc='right hand side') + 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 - + alpha = inputs['alpha'] #magnet fraction in rotor + mu_0 = 0.001 - + bmat = np.zeros((3, 3)) rhs = np.zeros(3) - beta = n*n_p - + 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) - + #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, 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 - - #invert matrix + rhs[2] = (b_r*k_rn/(mu_r*mu_0)) - Cm + outputs['bmat'] = bmat + outputs['rhs'] = rhs + +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 + +class DiffLinearSolve(om.ExplicitComponent): + + 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 compute(self, inputs, outputs): + bmat = inputs['bmat'] + rhs = inputs['rhs'] + print(bmat) bmat_inv = np.linalg.inv(bmat) - - #solve system using inverse of system matrix - outputs['x'] = bmat_inv.dot(rhs) -class MotorAirFourierComp(om.Group()): - + x = bmat_inv.dot(rhs) + print(x) + outputs['x'] = x + + # def solve_linear(self, d_outputs, d_residuals, mode): + # if mode == 'fwd': + # d_outputs['x'] = self.inv_jac * d_residuals['x'] + # elif mode == 'rev': + # d_residuals['x'] = self.inv_jac * d_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') + + 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, outputs): + # i = 1 + +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') + + 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, outputs): + # i = 1 + +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'], promotes_outputs=['bmat','rhs']) + self.add_subsystem('DiffLinearSolve', DiffLinearSolve(), promotes_inputs=['bmat','rhs'], promotes_outputs=['x']) + self.add_subsystem('BarnComp', BarnComp(n = ng), promotes_inputs=['x','sta_ir','n_m'], promotes_outputs=['B_arn']) + 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']) + + +class InvFourier(om.ExplicitComponent): + + def setup(self): + self.add_input('B_tn', shape=N) + + self.add_output('B_t', shape=(N-1)*2) + + def compute(self, inputs, outputs): + 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') + + 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'] + outputs['B_msq'] + +class CoreLoss(om.ExplicitComponent): + def setup(self): - self.add_subsystem('DiffAssemble', DiffAssemble()) - self.add_subsystem('DiffLinearSolve', DiffLinearSolve()) - self.add_subsystem('BarnComp', BarnComp()) - self.add_subsystem('CorrectSlot', CorrectSlot()) - self.add_subsystem() - - + self.add_input('B_pk') + self.add_input('B_msq') + self.add_input('f') + + self.add_output('L_core') + + 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 + + +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.051)#, units='m') + ind.add_output('sta_ir', val = 0.050)#, 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.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')) + print(p.get_val('InvFourier.B_tn')) From 795bb6e1b3f1c73d2ccba8e18a3ef403b6d3f1ec Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 5 Aug 2019 02:35:43 -0400 Subject: [PATCH 43/46] working openmdao model --- bmatrix.py | 47 +++++++++++++++++++------------------- btooth.py | 39 +++++++++++++++---------------- notes.tex | 16 +++++++------ radial_air_fourier_comp.py | 32 +++++++++++++++++--------- 4 files changed, 73 insertions(+), 61 deletions(-) diff --git a/bmatrix.py b/bmatrix.py index adf7921..fc24c08 100644 --- a/bmatrix.py +++ b/bmatrix.py @@ -1,7 +1,7 @@ -''' +''' 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 +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 ''' @@ -23,10 +23,10 @@ def RadSinT(n, alpha): # 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 @@ -35,7 +35,7 @@ def RadialT(n, alpha): # Permitivity of Free Space. Cancels out in computation, so the value doesn't matter -mu_0 = 0.001 +mu_0 = 0.001 # Compute the coefficients of the assumed solution using the boundary conditions @@ -48,69 +48,70 @@ def B_fourier_terms(n, params): 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, 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 - + 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 index 8a96e5f..40128c4 100644 --- a/btooth.py +++ b/btooth.py @@ -23,11 +23,11 @@ def B_tooth(N, M, params): 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'] + 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) @@ -40,26 +40,26 @@ def B_tooth(N, M, params): 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 + # 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 + + # 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 + 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 @@ -69,7 +69,7 @@ def B_tooth(N, M, params): alpha = .85 b_r = 1.3 f = 5400/20 #rpm - + params = { 'gap':gap, 'r_s':r_s, @@ -89,15 +89,14 @@ def B_tooth(N, M, params): 'b_r':b_r, 'f':f } - - B_tn = B_tooth(17, 20, params) - + + 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(i) print(L_core) - plt.plot(i) - plt.show() - + # plt.plot(i) + # plt.show() diff --git a/notes.tex b/notes.tex index e8ab551..8b0ade8 100644 --- a/notes.tex +++ b/notes.tex @@ -58,11 +58,13 @@ \section{NEXT STEP} 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 +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) @@ -70,7 +72,7 @@ \section{High-Level Notes} \item Menial tasks: analytic derivatives, clean up code, split across .py files \end{itemize} -1. +1. Seems to be current focus. Going to look through MotorCAD manual to replicate their models with Aaron @@ -91,13 +93,13 @@ \section{High-Level Notes} Currently using design space bounds, -Upper limit on current density +Upper limit on current density -Lower limit on torque +Lower limit on torque Objective is losses, but don't make sense to use as they are -4. +4. Do these to pass the time Especially need to clean up comments, remove deprecated lines @@ -133,7 +135,7 @@ \section{Specific Notes} P_h = k_hfB_{pk}^n \end{equation} -or +or \begin{equation} P_h = k_hfB_{pk}^{n+mB_{pk}} @@ -159,7 +161,7 @@ \section{Specific Notes} 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. +Model appears to suffer from subtractive cancellation due to exponential increase with $n$. Higher frequency terms destroy accuracy. diff --git a/radial_air_fourier_comp.py b/radial_air_fourier_comp.py index be70354..13744a7 100644 --- a/radial_air_fourier_comp.py +++ b/radial_air_fourier_comp.py @@ -8,6 +8,8 @@ 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') @@ -24,6 +26,15 @@ def setup(self): 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 @@ -50,7 +61,7 @@ def compute(self, inputs, outputs): #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)) @@ -105,9 +116,7 @@ def compute(self, inputs, outputs): rhs = inputs['rhs'] print(bmat) bmat_inv = np.linalg.inv(bmat) - x = bmat_inv.dot(rhs) - print(x) outputs['x'] = x # def solve_linear(self, d_outputs, d_residuals, mode): @@ -131,6 +140,7 @@ def setup(self): def compute(self, inputs, outputs): x = inputs['x'] + print(x) n = self.options['n'] r_s = inputs['sta_ir'] #stator inner radius n_p = inputs['n_m']/2 #number of pole pairs @@ -282,6 +292,7 @@ def setup(self): self.add_output('B_t', shape=(N-1)*2) def compute(self, inputs, outputs): + print(inputs['B_tn']) B_t = np.fft.irfft(inputs['B_tn']) outputs['B_t'] = B_t @@ -309,8 +320,8 @@ def compute(self, inputs, outputs): for t in range(1, Nl): B_msq += 2*((t*omega*B_t[t])**2) - outputs['B_pk'] - outputs['B_msq'] + outputs['B_pk'] = B_pk + outputs['B_msq'] = B_msq class CoreLoss(om.ExplicitComponent): @@ -325,7 +336,7 @@ def compute(self, inputs, outputs): B_pk = inputs['B_pk'] B_msq = inputs['B_msq'] f = inputs['f'] - + print(B_pk) # k_h, k_e, n, m are material constants and need to be fitted to real material data n = 2 m = 0.01 @@ -347,8 +358,8 @@ def compute(self, inputs, outputs): 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.051)#, units='m') - ind.add_output('sta_ir', val = 0.050)#, units='m') + 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) @@ -357,8 +368,8 @@ def compute(self, inputs, outputs): ind.add_output('w_tb', val = .005) ind.add_output('tt', val = .6*np.pi/10) - model.nonlinear_solver = om.NewtonSolver() - model.linear_solver = om.DirectSolver() + #model.nonlinear_solver = om.NewtonSolver() + #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): @@ -376,4 +387,3 @@ def compute(self, inputs, outputs): p.run_model() print(p.get_val('L_core')) - print(p.get_val('InvFourier.B_tn')) From 4bbcb2beddb1d85eda3490878873d4ccaebcf97c Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 5 Aug 2019 04:49:55 -0400 Subject: [PATCH 44/46] some partials added --- radial_air_fourier_comp.py | 167 +++++++++++++++++++++++++++++++++++-- 1 file changed, 161 insertions(+), 6 deletions(-) diff --git a/radial_air_fourier_comp.py b/radial_air_fourier_comp.py index 13744a7..cee2088 100644 --- a/radial_air_fourier_comp.py +++ b/radial_air_fourier_comp.py @@ -81,6 +81,98 @@ def compute(self, inputs, outputs): 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 + + 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 @@ -99,6 +191,14 @@ def RadialR(n, alpha): 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 @@ -138,9 +238,10 @@ def setup(self): 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'] - print(x) n = self.options['n'] r_s = inputs['sta_ir'] #stator inner radius n_p = inputs['n_m']/2 #number of pole pairs @@ -152,8 +253,23 @@ def compute(self, inputs, outputs): 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, outputs): - # i = 1 + 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): @@ -245,6 +361,11 @@ def setup(self): 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'] @@ -269,8 +390,8 @@ def compute(self, inputs, outputs): B_tna[n] = B_tn outputs['B_tn'] = B_tna - # def compute_partials(self, inputs, outputs): - # i = 1 + def compute_partials(self, inputs, J): + i = 1 class MotorAirFourierComp(om.Group): def initialize(self): @@ -291,6 +412,7 @@ def setup(self): 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']) @@ -305,6 +427,9 @@ def setup(self): 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'] @@ -323,6 +448,22 @@ def compute(self, inputs, outputs): 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): @@ -332,11 +473,13 @@ def setup(self): 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'] - print(B_pk) # k_h, k_e, n, m are material constants and need to be fitted to real material data n = 2 m = 0.01 @@ -348,6 +491,18 @@ def compute(self, inputs, outputs): 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() From 39a0d199082339f68697b5f60d9718cd33432bbe Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 5 Aug 2019 09:43:42 -0400 Subject: [PATCH 45/46] numpy solver --- radial_air_fourier_comp.py | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/radial_air_fourier_comp.py b/radial_air_fourier_comp.py index cee2088..dd5b395 100644 --- a/radial_air_fourier_comp.py +++ b/radial_air_fourier_comp.py @@ -261,7 +261,7 @@ def compute_partials(self, inputs, J): mu_0 = 0.001 Da = -x[0] - dDadx = [-1 0 0] + dDadx = [-1, 0, 0] beta = n*n_p Ea = -(r_s)**(2*beta) dEadsta = -(2*beta)*(r_s**(2*beta-1)) @@ -391,7 +391,32 @@ def compute(self, inputs, outputs): outputs['B_tn'] = B_tna def compute_partials(self, inputs, J): - i = 1 + 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) + + 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): From 08207e62105d355af8f829f4d07b3d8df54e0a72 Mon Sep 17 00:00:00 2001 From: garobed1 Date: Mon, 5 Aug 2019 11:25:43 -0400 Subject: [PATCH 46/46] working linear solve --- radial_air_fourier_comp.py | 68 +++++++++++++++++++++++++------------- 1 file changed, 45 insertions(+), 23 deletions(-) diff --git a/radial_air_fourier_comp.py b/radial_air_fourier_comp.py index dd5b395..7d24d90 100644 --- a/radial_air_fourier_comp.py +++ b/radial_air_fourier_comp.py @@ -96,6 +96,8 @@ def compute_partials(self, inputs, J): 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)) @@ -202,28 +204,29 @@ def RadialRPart(n, alpha): def RadialT(n, alpha): return 0 -class DiffLinearSolve(om.ExplicitComponent): - - 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 compute(self, inputs, outputs): - bmat = inputs['bmat'] - rhs = inputs['rhs'] - print(bmat) - bmat_inv = np.linalg.inv(bmat) - x = bmat_inv.dot(rhs) - outputs['x'] = x - - # def solve_linear(self, d_outputs, d_residuals, mode): +# 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': - # d_outputs['x'] = self.inv_jac * d_residuals['x'] + # outputs['x'] = self.inv_jac * residuals['x'] # elif mode == 'rev': - # d_residuals['x'] = self.inv_jac * d_outputs['x'] + # residuals['x'] = self.inv_jac * outputs['x'] @@ -402,6 +405,14 @@ def compute_partials(self, inputs, J): 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)))) @@ -424,11 +435,20 @@ def initialize(self): 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'], promotes_outputs=['bmat','rhs']) - self.add_subsystem('DiffLinearSolve', DiffLinearSolve(), promotes_inputs=['bmat','rhs'], promotes_outputs=['x']) - self.add_subsystem('BarnComp', BarnComp(n = ng), promotes_inputs=['x','sta_ir','n_m'], promotes_outputs=['B_arn']) + 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): @@ -549,6 +569,8 @@ def compute_partials(self, inputs, J): 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'])