Merge pull request #18 from AnttiLehikoinen/WIP

Wip

Examples/Surf-PMSM/Examples and testing/compute_inductances_2.m

@@ -0,0 +1,47 @@
+%computing dq-axis stuff
+ra = phi_bias;
+rotM = [cos(-ra) -sin(-ra);
+    sin(-ra) cos(-ra)];
+angles_p = -[0 2/3 4/3]*pi;
+Tp = 2/3*[cos(angles_p); -sin(angles_p)];
+%Phidq = rotM*Tp*simc.matrices.Ls'*dimsc.leff*simc.matrices.Cs'*A(1:simc.Np,:) * mshc.symmetrySectors;
+%Idq = rotM*Tp*Imain;
+
+%rotorAngles = wm*pars.ts;
+rotorAngles = pars.rotorAngle;
+Nsamples = numel(rotorAngles);
+
+Phidq = zeros(2, Nsamples);
+Idq = zeros(2, Nsamples);
+
+A = sim.results.Xs;
+%Imain = sim.Is;
+for k = 1:Nsamples
+    ra = phi_bias + rotorAngles(k)*dimsc.p;
+    rotM = [cos(-ra) -sin(-ra);
+        sin(-ra) cos(-ra)];
+
+    Phidq(:,k) = rotM*Tp*sim.matrices.Ls'*dimsc.leff*sim.matrices.Cs'*A(1:sim.Np,k) * mshc.symmetrySectors;
+    Idq(:,k) = rotM*Tp*Imain(:,k);
+end
+
+figure(7); clf; hold on;
+plot(rotorAngles, Phidq' );
+%plot(rotorAngles, Phi_d(:, 56), 'bx-');
+%plot(rotorAngles, Phi_q(:, 56), 'rx-');
+
+figure(8); clf;
+plot( Idq')
+
+%computing
+Id = Idq(1,:);
+Iq = Idq(2,:);
+Phi_d = Phidq(1,:);
+Phi_q = Phidq(2,:);
+c = [[ones(numel(Id),1) Id' 0*Iq']; 
+    [zeros(numel(Id),1) 0*Id' Iq']] \ [Phi_d'; Phi_q']
+%Tb = 3/2*dimsc.p * ( -Idq(1,:).*Phidq(2,:) + Idq(2,:).*Phidq(1,:) );
+%figure(6); clf; hold on; box on; grid on;
+%plot(angles/pi*180, T, 'b')
+%plot(angles/pi*180, Tb, 'r')
+

Examples/Surf-PMSM/Examples and testing/plot_Bag.m

@@ -5,8 +5,8 @@
 
 %finding nodes to plot
 angles = linspace(0, pi/4, 1000);
-%angles = linspace(-1/12, 1/12)*2*pi + pi/2;
 r = 0.55*dim.Sin + 0.45*dim.Rout;
+%r = 0.45*dim.Sin + 0.55*dim.Rout;
 x = r*cos(angles);
 y = r*sin(angles);
 [~, els] = pointSources2Elements_faster(x,y,msh_ag); angles = angles(els>0);els = els(els>0);
@@ -20,7 +20,7 @@
 Brad_all = Babs_all;
 Bcirc_all = Babs_all;
 
-for k = 1:1:Nsamples
+for k = 1:Nsamples
 
     %computing B in the airgap
     A_ag = mshc.bandData.tag_solution(sim.results.Xt(:,k));
@@ -56,5 +56,9 @@
 figure(6); clf; hold on; box on;
 %drawFluxDensity(mshc, simc.results.Xh, 'LineStyle', 'none'); colormap('jet'); colorbar;% caxis([0 2])
 drawFluxDensity(msh_ag, A_ag, 'LineStyle', 'none'); colorbar;
+msh_triplot(mshc, find(mshc.matel==2), 'k');
 
-plot(x, y, 'k.-')
+plot(x, y, 'k.-')
+
+figure(7); clf;
+surf(Babs_all)

Examples/Surf-PMSM/Examples and testing/setup.m

@@ -90,12 +90,16 @@
 
 %determining supply angles
 phi_rotor = pi/8; %angle of rotor axis
-phi_stator = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q*(1/2 + 1/2*(1 - 0.5*dimsc.c)); %angle of stator winding axis
-phi_bias = -phi_rotor + phi_stator;
+%phi_stator = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q*(1/2 + 1/2*(1 - 0.5*dimsc.c)); %angle of stator winding axis
+%phi_bias = -phi_rotor + phi_stator;
+
+angle_layer_1 = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q - 2*pi/dimsc.Qs*dimsc.c; %chorded layer
+angle_layer_2 = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q; %non-chorded layer
+phi_stator = 0.5*angle_layer_1 + 0.5*angle_layer_2;
 
 sim = MachineSimulation(mshc, dimsc);
 pars = SimulationParameters('U', 525, 'f', 200, 'slip', 0, 'N_periods', 2, 'N_stepsPerPeriod', 600, 'isDC', true, ...
-    'phi0', phi_bias-pi/2 - pi/180*15);
+    'phi0', phi_bias-pi/2 - pi/180*11);
 
 %uncomment for elegant no-load simulation
 %sim.matrices.Zew_s = eye(3)*1e9;
@@ -130,7 +134,7 @@
 
 %torque plot
 figure(5); clf; hold on; box on;
-%T_rated = sim_compute_torque(sim, pars, 'stepping');
+T_rated = sim_compute_torque(sim, pars, 'stepping');
 plot(pars.ts, T_rated);
 wm = 2*pi*pars.f/dimsc.p;
 Pave = mean(T_rated)*wm;
@@ -159,4 +163,5 @@
 plot(ts, Uplot(2, :), 'r--');
 plot(ts, Uplot(3, :), 'k--');
 
-%Idq_plot;
+%Idq_plot;
+save('data_rated.mat', 'I_rated', 'T_rated', 'E_rated', 'Phi_rated', 'Babs_all', 'Bcirc_all', 'Brad_all');

Examples/Surf-PMSM/Examples and testing/setup_currentSupply.m

@@ -93,11 +93,18 @@
 phi_stator = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q*(1/2 + 1/2*(1 - 0.5*dimsc.c)); %angle of stator winding axis
 phi_bias = -phi_rotor + phi_stator;
 
-Imain = 5*[cos(phi_bias); cos(phi_bias-2*pi/3); cos(phi_bias-4*pi/3)];
+angle_layer_1 = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q - 2*pi/dimsc.Qs*dimsc.c; %chorded layer
+angle_layer_2 = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q; %non-chorded layer
+phi_stator = 0.5*angle_layer_1 + 0.5*angle_layer_2;
+phi_bias = (phi_rotor - phi_stator)*dimsc.p;
+
+angles = linspace(0, 2*pi, 100);
+Imain = sqrt(2)*6*[cos(phi_bias+angles); cos(phi_bias-2*pi/3+angles); cos(phi_bias-4*pi/3+angles)];
+Imain = Imain + randn(size(Imain));
 
 msch.matel(mshc.rotel) = 0;
 sim = MachineSimulation(mshc, dimsc);
-pars = SimulationParameters('Is', Imain, 'slip', 1, 'rotorAngle', zeros(1, size(Imain,2)));
+pars = SimulationParameters('Is', Imain, 'slip', 1, 'rotorAngle', angles/dimsc.p);
 
 %harmonic simulation
 sim.run_static(pars);

Examples/Surf-PMSM/Examples and testing/setup_noload.m

@@ -90,8 +90,14 @@
 
 %determining supply angles
 phi_rotor = pi/8; %angle of rotor axis
-phi_stator = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q*(1/2 + 1/2*(1 - 0.5*dimsc.c)); %angle of stator winding axis
-phi_bias = -phi_rotor + phi_stator;
+%phi_stator = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q*(1/2 + 1/2*(1 - 0.5*dimsc.c)); %angle of stator winding axis
+%phi_bias = -phi_rotor + phi_stator;
+
+angle_layer_1 = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q - 2*pi/dimsc.Qs*dimsc.c; %chorded layer
+angle_layer_2 = pi/(2*dimsc.p) + pi/dimsc.Qs*dimsc.q; %non-chorded layer
+phi_stator = 0.5*angle_layer_1 + 0.5*angle_layer_2;
+
+phi_bias = (phi_rotor - phi_stator)*dimsc.p;
 
 sim = MachineSimulation(mshc, dimsc);
 pars = SimulationParameters('U', 525, 'f', 200, 'slip', 0, 'N_periods', 2, 'N_stepsPerPeriod', 600, 'isDC', true, ...
@@ -159,4 +165,5 @@
 plot(ts, Uplot(2, :), 'r--');
 plot(ts, Uplot(3, :), 'k--');
 
-%Idq_plot;
+%Idq_plot;
+%save('data_noload.mat', 'I_noload', 'T_noload', 'E_noload', 'Phi_noload', 'Babs_all', 'Bcirc_all', 'Brad_all');

SMEKlib/Classes/Elements and Shape Functions/Elements.m

@@ -39,6 +39,13 @@
             end
             bl = false;
         end
+        function bl = isPrism(type)
+            if type == Elements.prism;
+                bl = true;
+                return;
+            end
+            bl = false;
+        end
 
         function X = refPoints_nodes(type)
             switch type