EENG 577 M3 Synchronous Generator

E. Hildenbrandt D. Davis J. Brownlee

26 January 2025

PART 1

1

Fig. 1. Cross-section of generator, not to scale. The quadrature (q) axis is 90° electrically from the direct (d) axis. Because this is a 4-pole generator, there are 720° electrical for every 360° mechanical. The analysis would be the same if a, b, c, d, and q were rotated a multiple of 90° mechanical (180° electrical).

2.1: State Space, Case 1 (no damping)

Assuming the generator is balanced and rotating at f_{s}=60Hz:

L_{aa} ≅ L_{bb} ≅ L_{cc} ≅ L_{sa} + L_{sv} cos(2θ + ϕ)
where [ϕ_{a}, ϕ_{b}, ϕ_{c}] = [0, −\frac{4π}{3} , −\frac{2π}{3} ] and θ = ωt = 377 · t

Ignoring the damping windings the state space equations can be written as such:

[\begin{matrix}v_{a} \\ v_{b} \\ v_{c} \\ v_{f} \\ v_{kd} \\ v_{kq} \end{matrix}] = [\begin{matrix}r_{s} + r_{l} & 0 & 0 & 0 & 0 & 0 \\ 0 & r_{s} + r_{l} & 0 & 0 & 0 & 0 \\ 0 & 0 & r_{s} + r_{l} & 0 & 0 & 0 \\ 0 & 0 & 0 & r_{f} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}] · [\begin{matrix}i_{a} \\ i_{b} \\ i_{c} \\ i_{f} \\ i_{kd} \\ i_{kq} \end{matrix}] +

\frac{d}{dt} [\begin{matrix}L_{sa} + L_{sv} cos(2θ) + L_{l} & −L_{ma} + L_{mv} cos(2θ − 2π/3) & −L_{ma} + L_{mv} cos(2θ − 4π/3) & L_{afm} cos(2θ) & 0 & 0 \\ −L_{ma} + L_{mv} cos(2θ − 2π/3) & L_{sa} + L_{sv} cos(2θ − 4π/3) + L_{l} & −L_{ma} + L_{mv} cos(2θ) & L_{afm} cos(2θ − 2π/3) & 0 & 0 \\ −L_{ma} + L_{mv} cos(2θ − 4π/3) & −L_{ma} + L_{mv} cos(2θ) & L_{sa} + L_{sv} cos(2θ − 2π/3) + L_{l} & L_{afm} cos(2θ − 4π/3) & 0 & 0 \\ L_{afm} cos(2θ) & L_{afm} cos(2θ − 2π/3) & L_{afm} cos(2θ − 4π/3) & L_{ff} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}] · [\begin{matrix}i_{a} \\ i_{b} \\ i_{c} \\ i_{f} \\ i_{kd} \\ i_{kq} \end{matrix}]

2.2: State Space, Case 2

[\begin{matrix}v_{a} \\ v_{b} \\ v_{c} \\ v_{f} \\ v_{kd} \\ v_{kq} \end{matrix}] = [\begin{matrix}r_{s} + r_{l} & 0 & 0 & 0 & 0 & 0 \\ 0 & r_{s} + r_{l} & 0 & 0 & 0 & 0 \\ 0 & 0 & r_{s} + r_{l} & 0 & 0 & 0 \\ 0 & 0 & 0 & r_{f} & 0 & 0 \\ 0 & 0 & 0 & 0 & r_{k}d & 0 \\ 0 & 0 & 0 & 0 & 0 & r_{k}q \end{matrix}] · [\begin{matrix}i_{a} \\ i_{b} \\ i_{c} \\ i_{f} \\ i_{kd} \\ i_{kq} \end{matrix}] +

\frac{d}{dt} [\begin{matrix}L_{sa} + L_{sv} cos(2θ) + L_{l} & −L_{ma} + L_{mv} cos(2θ − 2π/3) & −L_{ma} + L_{mv} cos(2θ − 4π/3) & L_{afm} cos(2θ) & L_{adkm} cos(2θ) & −L_{akqm1} cos(2θ) \\ −L_{ma} + L_{mv} cos(2θ − 2π/3) & L_{sa} + L_{sv} cos(2θ − 4π/3) + L_{l} & −L_{ma} + L_{mv} cos(2θ) & L_{afm} cos(2θ − 2π/3) & L_{akdm} cos(2θ − 2π/3) & −L_{akqm1} cos(2θ − 2π/3) \\ −L_{ma} + L_{mv} cos(2θ − 4π/3) & −L_{ma} + L_{mv} cos(2θ) & L_{sa} + L_{sv} cos(2θ − 2π/3) + L_{l} & L_{afm} cos(2θ − 4π/3) & L_{akdm} cos(2θ − 4π/3) & −L_{akqm1} cos(2θ − 4π/3) \\ L_{afm} cos(2θ) & L_{afm} cos(2θ − 2π/3) & L_{afm} cos(2θ − 4π/3) & L_{ff} & L_{kdf} & L_{kqf} \\ L_{adkm} cos(2θ) & L_{akdm} cos(2θ − 2π/3) & L_{akdm} cos(2θ − 4π/3) & L_{kdf} & L_{kd1kd1} & 0 \\ −L_{akqm1} cos(2θ) & −L_{akqm1} cos(2θ − 2π/3) & −L_{akqm1} cos(2θ − 4π/3) & L_{kqf} & 0 & L_{kq1kq1} \end{matrix}] · [\begin{matrix}i_{a} \\ i_{b} \\ i_{c} \\ i_{f} \\ i_{kd} \\ i_{kq} \end{matrix}]

3.1

V = RI + \frac{d}{dt}LI  (1)
A = −L^{−1}R  (2)
B = L^{−1}  (3)

In compact matrix from.

X′ = AX + BU

Where without damping

A = − [\begin{matrix}1.95 + 0.05 cos (2θ) + L_{l} & −0.80 − 0.05 cos (2θ − \frac{2π}{3} ) & −0.80 − 0.05 cos (2θ − \frac{4π}{3} ) & 26.1 cos (2θ) \\ −0.80 − 0.05 cos (2θ − \frac{2π}{3} ) & 1.95 + 0.05 cos (2θ) + L_{l} & −0.80 − 0.05 cos (2θ − \frac{4π}{3} ) & 26.1 cos (2θ − \frac{2π}{3} ) \\ −0.80 − 0.05 cos (2θ − \frac{4π}{3} ) & −0.80 − 0.05 cos (2θ − \frac{2π}{3} ) & 1.95 + 0.05 cos (2θ) + L_{l} & 26.1 cos (2θ − \frac{4π}{3} ) \\ 26.1 cos (2θ) & 26.1 cos (2θ − \frac{2π}{3} ) & 26.1 cos (2θ − \frac{4π}{3} ) & 444 \end{matrix}] ^{−}1 ·

([\begin{matrix}7.4 · 10^{−4} + r_{l} & 0 & 0 & 0 \\ 0 & 7.4 · 10^{−4} + r_{l} & 0 & 0 \\ 0 & 0 & 7.4 · 10^{−4} + r_{l} & 0 \\ 0 & 0 & 0 & 0.086 \end{matrix}] + [\begin{matrix}−0.1 sin (2θ) + L_{l} & −0.1 sin (2θ − \frac{2π}{3} ) & −0.1 sin (2θ − \frac{4π}{3} ) & −52.2 sin (2θ) \\ −0.1 sin (2θ − \frac{2π}{3} ) & −0.1 sin (2θ) + L_{l} & −0.1 sin (2θ − \frac{4π}{3} ) & −52.2 sin (2θ − \frac{2π}{3} ) \\ −0.1 sin (2θ − \frac{4π}{3} ) & −0.1 sin (2θ − \frac{2π}{3} ) & −0.1 sin (2θ) & −52.2 sin (2θ − \frac{4π}{3} ) \\ −52.2 sin (2θ) & −52.2 sin (2θ − \frac{2π}{3} ) & −52.2 sin (2θ − \frac{4π}{3} ) & 0 \end{matrix}]) ·10^{−3}H^{−1}

and

B = [\begin{matrix}1.95 + 0.05 cos (2θ) + L_{l} & −0.80 − 0.05 cos (2θ − \frac{2π}{3} ) & −0.80 − 0.05 cos (2θ − \frac{4π}{3} ) & 26.1 cos (2θ) \\ −0.80 − 0.05 cos (2θ − \frac{2π}{3} ) & 1.95 + 0.05 cos (2θ) + L_{l} & −0.80 − 0.05 cos (2θ − \frac{4π}{3} ) & 26.1 cos (2θ − \frac{2π}{3} ) \\ −0.80 − 0.05 cos (2θ − \frac{4π}{3} ) & −0.80 − 0.05 cos (2θ − \frac{2π}{3} ) & 1.95 + 0.05 cos (2θ) + L_{l} & 26.1 cos (2θ − \frac{4π}{3} ) \\ 26.1 cos (2θ) & 26.1 cos (2θ − \frac{2π}{3} ) & 26.1 cos (2θ − \frac{4π}{3} ) & 444 \end{matrix}]^{−1} ·10^{−3}H^{−1}

With damping

A = − [\begin{matrix}7.4 · 10^{−4} + r_{l} & 0 & 0 & 0 & 0 & 0 \\ 0 & 7.4 · 10^{−4} + r_{l} & 0 & 0 & 0 & 0 \\ 0 & 0 & 7.4 · 10^{−4} + r_{l} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.086 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1.58 · 10^{−4} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1.227 · 10^{−4} \end{matrix}]

−0.1 sin (2θ) + L_{l} −0.1 sin (2θ − \frac{2π}{3} )  −0.1 sin (2θ − \frac{4π}{3} )  −52.2 sin (2θ)  −0.894 sin (2θ)  0.894 sin (2θ)
−0.1 sin (2θ − \frac{2π}{3} ) −0.1 sin (2θ) + L_{l}  −0.1 sin (2θ − \frac{4π}{3} ) −52.2 sin (2θ − \frac{2π}{3} ) −0.894 sin (2θ − \frac{2π}{3} ) 0.894 sin (2θ − \frac{2π}{3} )
−0.1 sin (2θ − \frac{4π}{3} ) −0.1 sin (2θ − \frac{2π}{3} )  −0.1 sin (2θ) + L_{l} −52.2 sin (2θ − \frac{4π}{3} ) −0.894 sin (2θ − \frac{4π}{3} ) 0.894 sin (2θ − \frac{4π}{3} )
−52.2 sin (2θ)  −52.2 sin (2θ − \frac{2π}{3} ) −52.2 sin (2θ − \frac{4π}{3} )  0  0  0
−0.894 sin (2θ) −0.894 sin (2θ − \frac{2π}{3} ) −0.894 sin (2θ − \frac{4π}{3} )  0  0  0
0.894 sin (2θ)  0.894 sin (2θ − \frac{2π}{3} ) 0.894 sin (2θ − \frac{4π}{3} )  0  0  0
×10^{−3}H^{−1}·

and

b = [\begin{matrix}1.95 + 0.05 cos(2θ) + L_{l} & −0.80 − 0.05 cos(2θ − 2π/3) & −0.08 − 0.05 cos(2θ − 4π/3) & 26.1 cos(2θ) & 0.447 cos(2θ) & −0.670 cos(2θ) \\ −0.80 − 0.05 cos(2θ − 2π/3) & 1.95 + 0.05 cos(2θ − 4π/3) + L_{l} & −0.80 − 0.05 cos(2θ) & 26.1 cos(2θ − 2π/3) & 0.447 cos(2θ − 2π/3) & −0.670 cos(2θ − 2π/3) \\ −0.80 − 0.05 cos(2θ − 4π/3) & 1.95 + 0.05 cos(2θ) & −0.80 − 0.05 cos(2θ − 2π/3) + L_{l} & 26.1 cos(2θ − 4π/3) & 0.447 cos(2θ − 4π/3) & −0.670 cos(2θ − 4π/3) \\ 26.1 cos(2θ) & 26.1 cos(2θ − 2π/3) & 26.1 cos(2θ − 4π/3) & 444 & 6.29 & 5.05 \\ 0.447 cos(2θ) & 0.447 cos(2θ − 2π/3) & 0.447 cos(2θ − 4π/3) & 6.29 & 0.1254 & 0 \\ −0.670 cos(2θ) & −0.670 cos(2θ − 2π/3) & −0.670 cos(2θ − 4π/3) & 5.05 & 0 & 0.3762 \end{matrix}]^{−1} ×10^{−3}H^{−1}

PART 2

1.1

The necessary excitation field current is given by i_{f} = \frac{V_{f}}{r_{f}} = \frac{338}{0.0860} = 3930[A]. The initial conditions at time t = 0 can be calculated by dividing the given phase voltages on the bus bar with the magnitude of the load resistance and reactance.

1.2

I_{a} = \frac{S}{\sqrt{3}V_{phi}} = 32950[A]  V_{a} = \frac{20 ∗ 10^{3}}{\sqrt{(}3)} \sqrt{(}2)cos(ωt) = V_{m}cos(0)  i_{a} = \frac{V_{a}}{\sqrt{R_{l}^{2} + X_{l}^{2}}}
Rl = \frac{S · pf}{I_{A}^{2}} = 0.5463[Ω]  V_{b} = V_{m}cos(0 − 2π/3)  i_{b} = \frac{V_{b}}{\sqrt{R_{l}^{2} + X_{l}^{2}}}
Xl = \frac{S · \sqrt{1 − pf^{2}}}{I_{a}^{2}} = 0.2646[Ω]  V_{c} = V_{m}cos(0 − 4π/3)  i_{c} = \frac{V_{c}}{\sqrt{R_{l}^{2} + X_{l}^{2}}}

The initial conditions are put into a column vector as follows.

I⃗ = [\begin{matrix}26904 \\ −13452 \\ −13452 \\ 3930 \\ 0 \\ 0 \end{matrix}]

Fig. 2. Current values for i_{a}, i_{b}, i_{c}, and i_{f} (left) and the torque (right)

We found that without the damping windings, the generator speeds out of control which would effectively grenade itself once the mechanical limit of the structural components is reached. This is likely because of undamped resonance in the system causing instability. But with the damping windings the generator stays at a reasonable current. The Matlab code for both models can be viewed at the end of this document.

Fig. 3. Current values for i_{a}, i_{b}, i_{c}, and i_{f} (left) i_{kd} and i_{kq} (center) and the torque (right)

1 Appendix: MATLAB Code

 1 %------------------------------------------------------------------
 2 % Code  for SS  Model  without  dampers
 3 %------------------------------------------------------------------
 4 % define  synch  gen  values
 5 fq = 60;
 6 vf = 377;
 7 nrpm = fq*120/4; % speed  of  rotor  in rev/min
 8 n = nrpm * 0.10472; % speed  of  rotor  in rad/s

9

10 S = 659*10^6;
11 vline = 20*10^3;
12 vphi = vline/sqrt(3);
13 pf = 0.9;
14 w = 377;

15

16 Ia = S/sqrt(3)/vphi;
17 P = pf*S;
18 R = P/(Ia^2);
19 Q = S*sqrt(1-pf^2);
20 Xl = Q/(Ia^2);

21 L = Xl / w;

22

23 Vm = (20*10^3/sqrt(3))*sqrt(2);
24 iao = Vm/(sqrt(R^2 + Xl^2));
25 ibo = Vm*cos(-2*pi/3)/(sqrt(R^2 + Xl^2));
26 ico = Vm*cos(-4*pi/3)/(sqrt(R^2 + Xl^2));

27

28 ifo = 338 / 0.0860;
29 % ikdo = 0;
30 % ikqo = 0;

31

32 iniCon = [iao;ibo;ico;ifo]
33 tspan = [0  300];
34 [t,x] = ode45(@myOde, tspan, iniCon);

35 36

37 plot(t,x(:,1),t,x(:,2),t,x(:,3),t,x(:,4));
38 legend(' ia ' ,' ib ' ,' ic ' ,' if ' )

39

40 plot(t,((x(:,1).*(3*Vm*pf))./(n*Xl)))
41 legend(' Torque ' )

42

43 %plot(t,((3.*x(:,1).*Vm)./(n*Xl)))

44

45 function  dx = myOde(t,x)

46

47  w = 377;
48  S = 659*10^6;
49  pf = 0.9;
50  Ia = S/ 20*10^3;
51  P = 0.9*S;
52  Rl = P/(Ia^2);
53  Q = S*sqrt(1-pf^2);
54  Xl = Q/(Ia^2);
55  L = Xl / w;

56

57  % angle  given  time t

58

59  theta = w*t;                  % (2  theta)
60  a =  -1*(2*pi/3);       % (2  theta  - 2pi/3)
61  b =  -1*(4*pi/3);       % (2  theta  - 4pi/3)
62  % define  values  for  resistance  matrix
63  rs = 7.4*10^(-4);
64  rf = 0.0860;
65  rkd1 = 1.58*10^(-4);

66  rkq1 = 1.227*10^(-4);

67

68  % R = [rs+Rl, 0,      0, 0,   0,   0;
69  %       0,      rs+Rl, 0, 0,   0,   0;
70  %       0,      0,      rs+Rl, 0,   0,   0;
71  %       0,      0,      0,      rf, 0,   0;
72  %       0,      0,      0,      0,   0,   0;
73  %       0,      0,      0,      0,   0,   0];
74  R = diag([rs+Rl rs+Rl rs+Rl rf]);

75

76  % define  values  for  inductance  matrix
77  pwr = 10^(-3);
78  Lsa = 1.95*pwr;
79  Lma = 0.80*pwr;
80  Lsv = 0.05*pwr;
81  Lmv = 0.05*pwr;

82

83  Lafm = 26.1*pwr;
84  % Lakdm1 = 0.447*pwr;
85  % Lakqm1 = 0.670*pwr;

86  %

87  % Lfkd1 = 5.05*pwr;
88  % Lkd1f = 6.29*pwr;

89

90  Lff = 444*pwr;
91  % Lkd1kd1 = 0.1254*pwr;
92  % Lkq1kq1 = 0.3762*pwr;

93

94  %------------------------------
95  Laa = Lsa + Lsv*cos(2*theta);
96  Lab = -Lma + Lmv*cos(2*theta + a);
97  Lac = -Lma + Lmv*cos(2*theta + b);
98  Laf = Lafm*cos(2*theta);
99  % Lakd = Lakdm1*cos(2*theta1);
100  % Lakq = -Lakqm1*cos(2*theta1);

101

102  % derivative  values
103  dLaa =  -2*Lsv*sin(2*theta);
104  dLab =  -2*Lmv*sin(2*theta + a);
105  dLac =  -2*Lmv*sin(2*theta + b);
106  dLaf =  -2*Lafm*sin(2*theta);
107  % dLakd =  -2*Lakdm1*sin(2*theta1);
108  % dLakq = 2*Lakqm1*sin(2*theta1);

109

110  %-----------------------------------

111  Lba = -Lma + Lmv*cos(2*theta + a);
112  Lbb = Lsa + Lsv*cos(2*theta + b);
113  Lbc = -Lma + Lmv*cos(2*theta);
114  Lbf = Lafm*cos(2*theta + a);
115  % Lbkd = Lakdm1*cos(2*theta2);
116  % Lbkq = -Lakqm1*cos(2*theta2);

117

118  % derivative
119  dLba =  -2*Lmv*sin(2*theta + a);
120  dLbb =  -2*Lsv*sin(2*theta + b);
121  dLbc =  -2*Lmv*sin(2*theta);
122  dLbf =  -2*Lafm*sin(2*theta + a);
123  % dLbkd =  -2*Lakdm1*sin(2*theta2);
124  % dLbkq = 2*Lakqm1*sin(2*theta2);

125

126  %-------------------------------
127  Lca = -Lma + Lmv*cos(2*theta + b);
128  Lcb = -Lma + Lmv*cos(2*theta);
129  Lcc = Lsa + Lsv*cos(2*theta + a);
130  Lcf = Lafm*cos(2*theta + b);
131  % Lckd = Lakdm1*cos(2*theta3);
132  % Lckq = -Lakqm1*cos(2*theta3);

133

134  % derivative
135  dLca =  -2*Lmv*sin(2*theta + b);
136  dLcb =  -2*Lmv*sin(2*theta);
137  dLcc =  -2*Lsv*sin(2*theta + a);
138  dLcf =  -2*Lafm*sin(2*theta + b);
139  % dLckd =  -2*Lakdm1*sin(2*theta3);
140  % dLckq = 2*Lakqm1*sin(2*theta3);

141

142  %---------------------------
143  Lfa = Lafm*cos(2*theta);
144  Lfb = Lafm*cos(2*theta + a);
145  Lfc = Lafm*cos(2*theta + b);

146

147  % derivative
148  dLfa =  -2*Lafm*sin(2*theta);
149  dLfb =  -2*Lafm*sin(2*theta + a);
150  dLfc =  -2*Lafm*sin(2*theta + b);

151

152  % %-------------------------------
153  % Lkda = Lakdm1*cos(2*theta1);
154  % Lkdb = Lakdm1*cos(2*theta2);
155  % Lkdc = Lakdm1*cos(2*theta3);

156  %

157  % % derivative
158  % dLkda =  -2*Lakdm1*sin(2*theta1);
159  % dLkdb =  -2*Lakdm1*sin(2*theta2);
160  % dLkdc =  -2*Lakdm1*sin(2*theta3);

161  %

162  % %-------------------------------
163  % Lkqa = -Lakqm1*cos(2*theta1);
164  % Lkqb = -Lakqm1*cos(2*theta2);
165  % Lkqc = -Lakqm1*cos(2*theta3);

166  %

167  % % derivative
168  % dLkqa = 2*Lakqm1*sin(2*theta1);
169  % dLkqb = 2*Lakqm1*sin(2*theta2);
170  % dLkqc = 2*Lakqm1*sin(2*theta3);

171

172  % phew  that  was a lot

173

174  % define  inductance  matrix  and  derivative  inductance  matrix

175

176  L = [Laa+L, Lab, Lac, Laf;
177  Lba, Lbb+L, Lbc, Lbf;
178  Lca, Lcb, Lcc+L, Lcf;
179  Lfa, Lfb, Lfc,    Lff];

180

181  dL = [dLaa, dLab, dLac, dLaf;
182  dLba, dLbb, dLbc, dLbf;
183  dLca, dLcb, dLcc, dLcf;
184  dLfa, dLfb, dLfc, 0];

185

186  % V = RI + dL*I + L*dI
187  % L*dI = V - RI - dL*I
188  % dI = -invL(R + dL)*I + invL*V

189

190  A =  -1*(R+dL)\L;
191  B = 1\L;
192  Vm = (20*10^3/sqrt(3))*sqrt(2);
193  u = [Vm*sin(w.*t+pi/2); Vm*sin(w.*t + pi/2 - 2*pi/3); Vm*sin(w.*t + pi/2 - 4*pi/3); 338];
194  %K = [Vm*cos(w.*t);Vm*cos(w.*t - (2*pi/3));Vm*cos(w.*t - (4*pi/3));338];
195  %K = [0;0;0;377];
196  %u = K.*sin(w*t + pi/2);

197

198  dx = A*x + B*u;
199 end

 1 %------------------------------------------------------------------
 2 % Code  for SS  Model  with  dampers
 3 %------------------------------------------------------------------
 4 % define  synch  gen  values
 5 fq = 60;
 6 vf = 377;
 7 nrpm = fq*120/4; % speed  of  rotor  in rev/min
 8 n = nrpm * 0.10472; % speed  of  rotor  in rad/s

9

10 S = 659*10^6;
11 vline = 20*10^3;
12 vphi = vline/sqrt(3);
13 pf = 0.9;
14 w = 377;

15

16 Ia = S/(sqrt(3)*vphi)
17 P = pf*S;
18 R = P/(Ia^2)
19 Q = S*sqrt(1-pf^2);
20 Xl = Q/(Ia^2)
21 L = Xl / w;

22

23 Vm = (20*10^(3)/sqrt(3))*sqrt(2);
24 iao = Vm/(sqrt(R^2 + Xl^2));
25 ibo = Vm*cos(-2*pi/3)/(sqrt(R^2 + Xl^2));
26 ico = Vm*cos(-4*pi/3)/(sqrt(R^2 + Xl^2));
27 ifo = 338 / 0.0860;
28 ikdo = 0;
29 ikqo = 0;

30

31 iniCon = [iao;ibo;ico;ifo;ikdo;ikqo]
32 tspan = [0  300];
33 [t,x] = ode45(@myOde, tspan, iniCon);

34

35 plot(t,x(:,1),t,x(:,2),t,x(:,3),t,x(:,4));
36 legend(' ia ' ,' ib ' ,' ic ' ,' if ' )

37

38 plot(t,x(:,5),t,x(:,6));
39 legend(' idk ' ,' idq ' )

40

41 plot(t,((x(:,1).*(3*Vm*pf))./(n*Xl)))
42 legend(' Torque ' )

43

44 function  dx = myOde(t,x)

46  w = 377;
47  S = 659*10^6;
48  vline = 20*10^3;
49  vphi = vline / sqrt(3);
50  pf = 0.9;
51  Ia = S/(sqrt(3) * vphi);
52  P = 0.9*S;
53  Rl = P/(Ia^2);
54  Q = S*sqrt(1-pf^2);
55  Xl = Q/(Ia^2);
56  L = Xl / w;

57

58  % angle  given  time t

59

60  theta = w*t;                  % (2  theta)
61  a =  -1*(2*pi/3);       % (2  theta  - 2pi/3)
62  b =  -1*(4*pi/3);       % (2  theta  - 4pi/3)
63  % define  values  for  resistance  matrix
64  rs = 7.4*10^(-4);
65  rf = 0.0860;
66  rkd1 = 1.58*10^(-4);
67  rkq1 = 1.227*10^(-4);

68

69  % R = [rs+R,0,0,0,0,0;
70  %      0,rs+R,0,0,0,0;
71  %      0,0,rs+R,0,0,0;
72  %      0,0,0,rf,0,0;
73  %      0,0,0,0,rkd1,0;
74  %      0,0,0,0,0,rkq1];
75  R = diag([rs+Rl, rs+Rl, rs+Rl, rf, rkd1, rkq1]);

76

77  % define  values  for  inductance  matrix
78  pwr = 10^(-3);
79  Lsa = 1.95*pwr;
80  Lma = 0.80*pwr;
81  Lsv = 0.05*pwr;
82  Lmv = 0.05*pwr;

83

84  Lafm = 26.1*pwr;
85  Lakdm1 = 0.447*pwr;
86  Lakqm1 = 0.670*pwr;

87

88  Lfkd1 = 5.05*pwr;
89  Lkd1f = 6.29*pwr;

 91  Lff = 444*pwr;
 92  Lkd1kd1 = 0.1254*pwr;
 93  Lkq1kq1 = 0.3762*pwr;

94

 95  %------------------------------
 96  Laa = Lsa + Lsv*cos(2*theta);
 97  Lab = -Lma + Lmv*cos(2*theta + a);
 98  Lac = -Lma + Lmv*cos(2*theta + b);
 99  Laf = Lafm*cos(2*theta);
100  Lakd = Lakdm1*cos(2*theta);
101  Lakq = -Lakqm1*cos(2*theta);

102

103  % derivative  values
104  dLaa =  -2*Lsv*sin(2*theta);
105  dLab =  -2*Lmv*sin(2*theta + a);
106  dLac =  -2*Lmv*sin(2*theta + b);
107  dLaf =  -2*Lafm*sin(2*theta);
108  dLakd =  -2*Lakdm1*sin(2*theta);
109  dLakq = 2*Lakqm1*sin(2*theta);

110

111  %-----------------------------------
112  Lba = -Lma + Lmv*cos(2*theta + a);
113  Lbb = Lsa + Lsv*cos(2*theta + b);
114  Lbc = -Lma + Lmv*cos(2*theta);
115  Lbf = Lafm*cos(2*theta + a);
116  Lbkd = Lakdm1*cos(2*theta + a);
117  Lbkq = -Lakqm1*cos(2*theta + a);

118

119  % derivative
120  dLba =  -2*Lmv*sin(2*theta + a);
121  dLbb =  -2*Lsv*sin(2*theta + b);
122  dLbc =  -2*Lmv*sin(2*theta);
123  dLbf =  -2*Lafm*sin(2*theta + a);
124  dLbkd =  -2*Lakdm1*sin(2*theta + a);
125  dLbkq = 2*Lakqm1*sin(2*theta + a);

126

127  %-------------------------------
128  Lca = -Lma + Lmv*cos(2*theta + b);
129  Lcb = -Lma + Lmv*cos(2*theta);
130  Lcc = Lsa + Lsv*cos(2*theta + a);
131  Lcf = Lafm*cos(2*theta + b);
132  Lckd = Lakdm1*cos(2*theta + b);
133  Lckq = -Lakqm1*cos(2*theta + b);

134

135  % derivative
136  dLca =  -2*Lmv*sin(2*theta + b);
137  dLcb =  -2*Lmv*sin(2*theta);
138  dLcc =  -2*Lsv*sin(2*theta + a);
139  dLcf =  -2*Lafm*sin(2*theta + b);
140  dLckd =  -2*Lakdm1*sin(2*theta + b);
141  dLckq = 2*Lakqm1*sin(2*theta + b);

142

143  %---------------------------
144  Lfa = Lafm*cos(2*theta);
145  Lfb = Lafm*cos(2*theta + a);
146  Lfc = Lafm*cos(2*theta + b);

147

148  % derivative
149  dLfa =  -2*Lafm*sin(2*theta);
150  dLfb =  -2*Lafm*sin(2*theta + a);
151  dLfc =  -2*Lafm*sin(2*theta + b);

152

153  %-------------------------------
154  Lkda = Lakdm1*cos(2*theta);
155  Lkdb = Lakdm1*cos(2*theta + a);
156  Lkdc = Lakdm1*cos(2*theta + b);

157

158  % derivative
159  dLkda =  -2*Lakdm1*sin(2*theta);
160  dLkdb =  -2*Lakdm1*sin(2*theta + a);
161  dLkdc =  -2*Lakdm1*sin(2*theta + b);

162

163  %-------------------------------
164  Lkqa = -Lakqm1*cos(2*theta);
165  Lkqb = -Lakqm1*cos(2*theta + a);
166  Lkqc = -Lakqm1*cos(2*theta + b);

167

168  % derivative
169  dLkqa = 2*Lakqm1*sin(2*theta);
170  dLkqb = 2*Lakqm1*sin(2*theta + a);
171  dLkqc = 2*Lakqm1*sin(2*theta + b);

172

173  % phew  that  was a lot

174

175  % define  inductance  matrix  and  derivative  inductance  matrix

176

177  L = [Laa+L, Lab, Lac, Laf, Lakd, Lakq;
178  Lba, Lbb+L, Lbc, Lbf, Lbkd, Lbkq;
179  Lca, Lcb, Lcc+L, Lcf, Lckd, Lckq;

180  Lfa, Lfb, Lfc,         Lff,    Lfkd1,    0;
181  Lkda, Lkdb, Lkdc,     Lkd1f, Lkd1kd1, Lkq1kq1;
182  Lkqa, Lkqb, Lkqc,     0,      0,         Lkq1kq1;];

183

184  dL = [dLaa, dLab, dLac, dLaf, dLakd, dLakq;
185  dLba, dLbb, dLbc, dLbf, dLbkd, dLbkq;
186  dLca, dLcb, dLcc, dLcf, dLckd, dLckq;
187  dLfa, dLfb, dLfc,     0,0,0;
188  dLkda, dLkdb, dLkdc, 0,0,0;
189  dLkqa, dLkqb, dLkqc, 0,0,0;];

190

191  % V = RI + dL*I + L*dI
192  % L*dI = V - RI - dL*I
193  % dI = -invL(R + dL)*I + invL*V

194

195  A =  -1*(R+dL)\L;
196  B = 1\L;
197  Vm = (20*10^(3)/sqrt(3))*sqrt(2);
198  u = [Vm*sin(w.*t+pi/2), Vm*sin(w.*t + pi/2 - 2*pi/3), Vm*sin(w.*t + pi/2 - 4*pi/3), 338, 0, 0] ';
199  %u = [Vm;Vm;Vm;377;0;0];
200  %u = K*sin(377.*t + pi/2);
201  dx = A*x + B*u;
202 end