1
Resumen— Los sistemas doblemente excitados son
comúnmente encontrados en los modelos de máquinas dichos
modelos pueden representarse para dos casos: con saliencia y sin
saliencia. Para las máquinas con saliencia se puede tener saliencia
en el rotor y estator liso o las condiciones invertidas, y para la
máquina sin saliencia ambos elementos son lisos. Debe recordarse
a la máquina como un transductor del que puede obtenerse
energía mecánica o energía eléctrica mediante la entrada
adecuada y la interacción de la energía interna en el modelo de la
máquina correspondiente. En este trabajo se considera un modelo
de cuarto orden, donde se encuentran dos ecuaciones de voltaje y
dos de movimiento, la solución de este dará lugar a la
determinación de diferentes parámetros que permiten el análisis
de la operación de la máquina. El modelo analizado es la máquina
con un grado de saliencia en el rotor y estator liso. El método de
integración numérica utilizado es el Runge Kutta de 5to. Orden
(Butcher).
Se modela un transductor doblemente excitado con un grado
de saliencia en el rotor propuesto en [Simmons, 1968] usando
ForTran 90 con el objetivo de estudiar los principios básicos
estudiados en la parte teórica de éste documento. Se muestran
resultados y por último, se presenta el código del programa.
I. MODELO EN FUNCIÓN DE CORRIENTES
as ecuaciones de balance de voltaje en el sistema
doblemente excitado mostrado en la figura 1 están
definidas de la siguiente forma
𝑉𝑠 = 𝑅𝑠𝑖𝑠 + 𝑝𝜓𝑠 (1)
𝑉𝑟 = 𝑅𝑟𝑖𝑟 + 𝑝𝜓𝑟 (2)
Las ecuaciones de flujo en el sistema doblemente excitado están
definidas de la siguiente forma
𝜓𝑠 = 𝐿𝑠𝑖𝑠 + 𝑀𝑠𝑟𝑖𝑟 (3)
𝜓𝑟 = 𝐿𝑟𝑖𝑟 + 𝑀𝑠𝑟𝑖𝑠 (4)
Sustituyendo (3) y (4) en (1) y (2), obtenemos las relaciones
siguientes
𝑉𝑠 = 𝑅𝑠𝑖𝑠 + 𝐿𝑠
𝑑𝑖𝑠𝑑𝑡
+ 𝑀𝑠𝑟
𝑑𝑖𝑟𝑑𝑡
+ 𝑖𝑠𝑑𝐿𝑠
𝑑𝑡+ 𝑖𝑟
𝑑𝑀𝑠𝑟
𝑑𝑡 (5)
𝑉𝑟 = 𝑅𝑟𝑖𝑟 + 𝐿𝑟
𝑑𝑖𝑟𝑑𝑡
+ 𝑀𝑠𝑟
𝑑𝑖𝑠𝑑𝑡
+ 𝑖𝑟𝑑𝐿𝑟
𝑑𝑡+ 𝑖𝑠
𝑑𝑀𝑠𝑟
𝑑𝑡 (6)
Este trabajo fue apoyado en parte por el Departamento de Ingeniería
Eléctrica de SEPI-ESIME Zacatenco, bajo el auspicio de la beca CONACYT.
Fig. 1. Transductor con campo magnético doblemente excitado [White and
Woodson, 1959].
Recordar que la matriz de coeficientes G es:
𝐺 =𝑑𝐿
𝑑𝜃
𝑖𝑟𝑝𝐿𝑟 = 𝑖𝑟𝑑𝐿𝑟
𝑑𝜃𝑟
𝑑𝜃𝑟
𝑑𝑡 𝑎𝑑𝑒𝑚𝑎𝑠 𝜔𝑟 = 𝑝𝜃𝑟 ⇒ 𝜔𝑟 =
𝑑𝜃𝑟
𝑑𝑡
𝑒𝑠 𝑙𝑜 𝑚𝑖𝑠𝑚𝑜 𝑞𝑢𝑒 𝑖𝑟𝑝𝐿𝑟 = 𝑖𝑟𝑑𝐿𝑟
𝑑𝜃𝑟
𝑑𝜃𝑟
𝑑𝑡
Rescribiendo las ecuaciones (1) y (2)
𝑉𝑠 = 𝑅𝑠𝑖𝑠 + 𝐿𝑠
𝑑𝑖𝑠𝑑𝑡
+ 𝑀𝑠𝑟
𝑑𝑖𝑟𝑑𝑡
+ 𝜔𝑟𝑖𝑠𝑑𝐿𝑠
𝑑𝜃𝑟+ 𝜔𝑟𝑖𝑟
𝑑𝑀𝑠𝑟
𝑑𝜃𝑟 (7)
𝑉𝑟 = 𝑅𝑟𝑖𝑟 + 𝐿𝑟
𝑑𝑖𝑟𝑑𝑡
+ 𝑀𝑠𝑟
𝑑𝑖𝑠𝑑𝑡
+ 𝜔𝑟𝑖𝑟𝑑𝐿𝑟
𝑑𝜃𝑟+ 𝜔𝑟𝑖𝑠
𝑑𝑀𝑠𝑟
𝑑𝜃𝑟 (8)
La forma general del sistema definida en función de las
corrientes es la siguiente
[𝑉𝑠𝑉𝑟
] = {[𝑅𝑠 00 𝑅𝑟
] + [𝐿𝑠 𝑀𝑠𝑟
𝑀𝑠𝑟 𝐿𝑟] 𝑝
+ 𝜔𝑟
𝑑
𝑑𝜃𝑟
[𝐿𝑠 𝑀𝑠𝑟
𝑀𝑠𝑟 𝐿𝑟]} [
𝑖𝑠𝑖𝑟
] (9)
[𝑉] = {[𝑅] + [𝐿]𝑝 + 𝜔𝑟[𝐺]}[𝑖] (10)
Fernando Ramos Albarrán
INSTITUTO POLITÉCNICO NACIONAL
Escuela Superior de Ingeniería Mecánica y Eléctrica, Sección de Estudios de Posgrado e Investigación
Unidad Profesional “Adolfo López Mateos”. Col. Lindavista, C. P. 07738, México D. F.
Tarea 3: Transductor Doblemente Excitado con
Saliencia
L
2
La matriz de coeficientes es
[𝐺] =𝑑
𝑑𝜃𝑟[𝐿𝑠 𝑀𝑠𝑟
𝑀𝑠𝑟 𝐿𝑟] =
𝑑[𝐿]
𝑑𝜃𝑟
La inversa de la matriz L está definida para el modelo sin
saliencia como
⇒ [𝐿]−1 = [
𝐿𝑟
𝐿𝑠𝐿𝑟 − 𝑀𝑠𝑟2−
𝑀𝑠𝑟
𝐿𝑠𝐿𝑟 − 𝑀𝑠𝑟2
−𝑀𝑠𝑟
𝐿𝑠𝐿𝑟 − 𝑀𝑠𝑟2
𝐿𝑠
𝐿𝑠𝐿𝑟 − 𝑀𝑠𝑟2
] (11)
La inversa de la matriz L está definida para el modelo con
saliencia como
⇒ [𝐿]−1
=
[
𝐿𝑟
(𝐿1 + 𝐿2𝑐𝑜𝑠2𝜃)𝐿𝑟 − 𝑀𝑚𝑎𝑥2𝑐𝑜𝑠2𝜃
−−𝑀𝑚𝑎𝑥𝑐𝑜𝑠𝜃
(𝐿1 + 𝐿2𝑐𝑜𝑠2𝜃)𝐿𝑟 − 𝑀𝑚𝑎𝑥2𝑐𝑜𝑠2𝜃
−−𝑀𝑚𝑎𝑥𝑐𝑜𝑠𝜃
(𝐿1 + 𝐿2𝑐𝑜𝑠2𝜃)𝐿𝑟 − 𝑀𝑚𝑎𝑥2𝑐𝑜𝑠2𝜃
𝐿𝑠
(𝐿1 + 𝐿2𝑐𝑜𝑠2𝜃)𝐿𝑟 − 𝑀𝑚𝑎𝑥2𝑐𝑜𝑠2𝜃 ]
(12)
A partir de la ecuación (10) se puede despejar las derivadas de
las corrientes para obtener una parte del modelo final tal y como
se muestra enseguida
[𝑉] = {[𝑅] + [𝐿]𝑝 + 𝜔𝑟[𝐺]}[𝑖]
[𝑉] = [𝑅][𝑖] + [𝐿]𝑝[𝑖] + 𝜔𝑟[𝐺][𝑖]
[𝐿]𝑝[𝑖] = −[𝑅][𝑖] − 𝜔𝑟[𝐺][𝑖] + [𝑉]
𝑝[𝑖] = −[𝐿]−1[𝑅][𝑖] − [𝐿]−1𝜔𝑟[𝐺][𝑖] + [𝐿]−1[𝑉]
𝑝[𝑖] = {−[𝐿]−1[𝑅] − [𝐿]−1𝜔𝑟[𝐺]}[𝑖] + [𝐿]−1[𝑉]
𝑝[𝑖] = {−[𝐿]−1[[𝑅] + 𝜔𝑟[𝐺]]}[𝑖] + [𝐿]−1[𝑉] (13)
Sustituyendo en la ecuación 11 los valores de inductancias para
el modelo sin saliencia se tienen las relaciones que se muestran
a continuación
⇒ 𝑝[𝑖] = {− [𝑥 𝑦𝑦 𝑧] [[
𝑅𝑠 00 𝑅𝑟
]
+ 𝜔𝑟 [0 −𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃
−𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃 0]]} [
𝑖𝑠𝑖𝑟
]
+ [𝑥 𝑦𝑦 𝑧] [
𝑉𝑠𝑉𝑟
] (12)
Sustituyendo en la ecuación 11 los valores de inductancias para
el modelo con saliencia se tienen las relaciones que se muestran
a continuación
⇒ 𝑝[𝑖] = {− [𝑥 𝑦𝑦 𝑧] [[
𝑅𝑠 00 𝑅𝑟
]
+ 𝜔𝑟 [0 −𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃
−𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃 −2𝐿2𝑠𝑒𝑛2𝜃]]} [
𝑖𝑠𝑖𝑟
]
+ [𝑥 𝑦𝑦 𝑧] [
𝑉𝑠𝑉𝑟
] (13)
Recordando la ecuación de oscilación, definida para la
operación de la máquina como motor
Para la ecuación de movimiento, el par mecánico se tomara
como constante, es necesario definir el par eléctrico, la ecuación
a definir es la siguiente
𝑇𝑒 = 𝑃. 𝑑𝑒𝑃. [1
2𝑖𝑠
2 𝑑𝐿𝑠
𝑑𝜃𝑟+
1
2𝑖𝑟
2 𝑑𝐿𝑟
𝑑𝜃𝑟+ 𝑖𝑠𝑖𝑟
𝑑𝑀𝑠𝑟
𝑑𝜃𝑟] (14)
Luego entonces la ecuación que define el par eléctrico en una
máquina sin saliencia es
𝑇𝑒 = −𝑖𝑠𝑖𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃
La ecuación que define el par eléctrico en una máquina con
saliencia es
𝑇𝑒 = −𝑖𝑠2𝐿2𝑠𝑒𝑛2𝜃 − 𝑖𝑠𝑖𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃
Recordando la ecuación de oscilación, definida para la
operación de la máquina como motor
𝐽𝑑𝜔
𝑑𝑡= 𝑇𝑒 − 𝑇𝑚 (15)
Para el modelo sin saliencia se tiene
𝑑𝜔
𝑑𝑡=
1
𝐽[−𝑖𝑠𝑖𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃] −
1
𝐽𝑇𝑚
Para el modelo con saliencia se tiene
𝑑𝜔
𝑑𝑡=
1
𝐽[−𝑖𝑠
2𝐿2𝑠𝑒𝑛2𝜃 − 𝑖𝑠𝑖𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃] −1
𝐽𝑇𝑚
El modelo final para el caso de la máquina síncrona sin
saliencia es el siguiente
𝑖�̇� = 𝑦𝑥𝑖𝑠(−𝑅𝑠) + 𝑦𝑥𝑖𝑟(𝜔𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃) + 𝑉𝑠𝑥 + 𝑉𝑟𝑦
𝑖�̇� = 𝑧𝑦𝑖𝑠(𝜔𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃) − 𝑧𝑦𝑖𝑟(𝑅𝑟) + 𝑉𝑠𝑦 + 𝑉𝑟𝑧 (16)
�̇� = −1
𝐽(−𝑖𝑠𝑖𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃) −
1
𝐽𝑇𝑚
�̇� = 𝜔
El modelo final para el caso de la máquina síncrona con
saliencia es el siguiente
𝑖�̇� = 𝑦𝑥𝑖𝑠(−𝑅𝑠 + 2𝜔𝑟𝐿2𝑠𝑒𝑛2𝜃) + 𝑦𝑥𝑖𝑟(𝜔𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃) + 𝑉𝑠𝑥+ 𝑉𝑟𝑦
𝑖�̇� = 𝑧𝑦𝑖𝑠(𝜔𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃) − 𝑧𝑦𝑖𝑟(𝑅𝑟) + 𝑉𝑠𝑦 + 𝑉𝑟𝑧 (17)
�̇� = −1
𝐽(𝑖𝑠
2𝐿2𝑠𝑒𝑛2𝜃) −1
𝐽(𝑖𝑠𝑖𝑟𝑀𝑚𝑎𝑥𝑠𝑒𝑛𝜃) −
1
𝐽𝑇𝑚
�̇� = 𝜔
El método de Runge Kutta de 5to Orden [Butcher, 1964] puede
ser escrito como:
3
(18)
Donde:
(19)
Existe una similitude entre el Método de Butcher y la regla de
Boole, éste presenta un error de truncamiento del O(h5), lo que
provee mayor exactitud. Cabe notar que se requieren 6
evaluaciones de pendientes.
II. RESULTADOS
A. Alimentando con fuentes de C.A. El estator y el rotor.
Los datos usados para los parámetros de este modelo fueron los
siguientes:
!******** DATOS DEL SISTEMA *************** Vr=15 !Magnitud de tensión en el rotor [V] Vs=100 !Magnitud de tensión en el estator [V] f=60.0 !Frecuencia (igual a cero en C.D.) [Hz] PI=3.1416 !Relación entre longitud y radio del rotor w0r= 2*PI*f !Velocidad sincrona del rotor [Rad/s] w0s= 2*PI*f !Velocidad sincrona del estator Lr= 0.09474 !Inductancia del rotor [H] Ls1= 0.09 !Inductancia en el estator [H] Ls2= 0.01 !Inductancia en el estator [H] M= 0.00136 !Inductancia mutua maxima [H] Rs= 8.8713 !Resistencia del estator [Ohms] Rr= 8.2606 !Resistencia del rotor [Ohms] J= 0.5 !Inercia [Kg m2] D=0.5 !Coeficiente de amortiguamiento [Nms/rad] Tm= 0.5 !Par mecánico [Nm] PdPolos=1.0 !Pares de Polos !*********************************************
Fig. 2. Datos de simulación en el tiempo para el caso A
Fig. 3. Gráficas de la corriente del rotor
0 1 2 3 4 5-4
-2
0
2
4
6
Tiempo (s)
Corr
iente
ir
(A)
Corriente del Rotor
0 0.5 1 1.5 2 2.5-4
-2
0
2
4
Tiempo (s)
Corr
iente
ir
(A)
Corriente del Rotor
0 2 4 6 8 10-6
-4
-2
0
2
4
6Corriente del estator
Tiempo (s)
Corr
iente
is (
A)
4
Fig. 4. Gráficas de la corriente del estator
Fig. 5. Gráficas de la velocidad angular
Fig. 6. Gráfica de la posición angular del rotor
Fig. 7. Gráficas del par eléctrico.
0 0.2 0.4 0.6 0.8 1
-4
-2
0
2
4
Corriente del estator
Tiempo (s)
Corr
iente
is (
A)
0 2 4 6 8 101036.98
1037
1037.02
1037.04
1037.06
1037.08
1037.1Velocidad angular del rotor
Tiempo (s)
Velo
cid
ad (
rad/s
)
0 0.02 0.04 0.06 0.08 0.1 0.12
1036.998
1036.999
1037
1037.001
1037.002
1037.003
Velocidad angular del rotor
Tiempo (s)
Velo
cid
ad (
rad/s
)
0 2 4 6 8 100
1000
2000
3000
4000
5000
6000
7000Posición angular del rotor
Posic
ión (
rad)
Tiempo (s)
0 2 4 6 8 10-0.2
-0.1
0
0.1
0.2
0.3Par Eléctrico
Par
elé
ctr
ico T
e (
Nm
)
Tiempo (s)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
Par Eléctrico
Par
elé
ctr
ico T
e (
Nm
)
Tiempo (s)
5
Fig. 8. Gráficas del par eléctrico - velocidad.
B. Alimentando con fuentes de C.A. El estator y el rotor.
Los datos usados para los parámetros de este modelo fueron los
siguientes, solo cambia el coeficiente de amortiguamiento D:
!******** DATOS DEL SISTEMA *************** Vr=15 !Magnitud de tensión en el rotor [V] Vs=100 !Magnitud de tensión en el estator [V] f=60.0 !Frecuencia (igual a cero en C.D.) [Hz] PI=3.1416 !Relación entre longitud y radio del rotor w0r= 2*PI*f !Velocidad sincrona del rotor [Rad/s] w0s= 2*PI*f !Velocidad sincrona del estator Lr= 0.09474 !Inductancia del rotor [H] Ls1= 0.09 !Inductancia en el estator [H] Ls2= 0.01 !Inductancia en el estator [H] M= 0.00136 !Inductancia mutua maxima [H] Rs= 8.8713 !Resistencia del estator [Ohms] Rr= 8.2606 !Resistencia del rotor [Ohms] J= 0.5 !Inercia [Kg m2] D=2.0 !Coeficiente de amortiguamiento [Nms/rad] Tm= 0.5 !Par mecánico [Nm] PdPolos=1.0 !Pares de Polos !*********************************************
Fig. 9. Datos de simulación en el tiempo para el caso B
Fig. 10. Gráficas de la corriente del rotor
1036.98 1037 1037.02 1037.04 1037.06 1037.08 1037.1-0.2
-0.1
0
0.1
0.2
0.3Par-Velocidad
Velocidad (rad/s)
Par
(Nm
)
1036.9995 1037 1037.0005 1037.001 1037.0015 1037.002-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
Par-Velocidad
Velocidad (rad/s)
Par
(Nm
)
0 5 10 15 20-4
-2
0
2
4
6Corriente del Rotor
Corr
iente
is (
A)
Tiempo (s)
0 0.2 0.4 0.6 0.8 1 1.2
-3
-2
-1
0
1
2
3
4
5
Corriente del Rotor
Corr
iente
is (
A)
Tiempo (s)
0 5 10 15 20-15
-10
-5
0
5
10
15Corriente del estator
Tiempo (s)
Corr
iente
is (
A)
6
Fig. 11. Gráficas de la corriente del estator
Fig. 12. Gráficas de la velocidad angular
Fig. 13. Gráfica de la posición angular del rotor
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
-10
-5
0
5
10
Corriente del estator
Tiempo (s)
Corr
iente
is (
A)
17.5 18 18.5 19 19.5-15
-10
-5
0
5
10
15Corriente del estator
Tiempo (s)
Corr
iente
is (
A)
0 5 10 15 203070
3080
3090
3100
3110
3120
3130
3140
Tiempo (s)
Velocidad angular del rotor
Velo
cid
ad (
rad/s
)
4.1 4.12 4.14 4.16 4.18 4.2 4.22 4.24 4.26
3091.4
3091.5
3091.6
3091.7
3091.8
3091.9
3092
Tiempo (s)
Velocidad angular del rotor
Velo
cid
ad (
rad/s
)
0 5 10 15 200
2000
4000
6000
8000
10000
Tiempo (s)
Posic
ión a
ngula
r (r
ad)
Posición angular del rotor
0 5 10 15 20-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
Tiempo (s)
Par
elé
ctr
ico (
Nm
)
Par eléctrico del rotor
7
Fig. 14. Gráficas del par eléctrico.
Fig. 15. Gráficas del par eléctrico - velocidad.
C. Alimentando con fuentes de C.D. El estator y el rotor.
Los datos usados para los parámetros de este modelo fueron
los siguientes
!******** DATOS DEL SISTEMA *************** Vr=15 !Magnitud de tensión en el rotor [V] Vs=100 !Magnitud de tensión en el estator [V] f=0.0 !Frecuencia (igual a cero en C.D.) [Hz] PI=3.1416 !Relación entre longitud y radio del rotor w0r= 2*PI*f !Velocidad sincrona del rotor [Rad/s] w0s= 2*PI*f !Velocidad sincrona del estator Lr= 0.09474 !Inductancia del rotor [H] Ls1= 0.09 !Inductancia en el estator [H] Ls2= 0.01 !Inductancia en el estator [H] M= 0.00136 !Inductancia mutua maxima [H] Rs= 8.8713 !Resistencia del estator [Ohms] Rr= 8.2606 !Resistencia del rotor [Ohms] J= 0.5 !Inercia [Kg m2] D=2.0 !Coeficiente de amortiguamiento [Nms/rad] Tm= 0.5 !Par mecánico [Nm] PdPolos=1.0 !Pares de Polos !********************************************* Además se desactivan las partes indicadas en el código del programa
Fig. 16. Datos de simulación en el tiempo para el caso C
0 0.5 1 1.5 2
-0.4
-0.2
0
0.2
0.4
0.6
Tiempo (s)
Par
elé
ctr
ico (
Nm
)Par eléctrico del rotor
17.6 17.8 18 18.2 18.4 18.6 18.8
-0.4
-0.2
0
0.2
0.4
0.6
0.8
Tiempo (s)
Par
elé
ctr
ico (
Nm
)
Par eléctrico del rotor
3070 3080 3090 3100 3110 3120 3130 3140-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
Velocidad (rad/seg)
Par
elé
ctr
ico T
e (
Nm
)
Par Eléctrico-Velocidad angular del rotor
3079 3079.2 3079.4 3079.6 3079.8 3080 3080.2 3080.4 3080.6 3080.8-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Velocidad (rad/seg)
Par
elé
ctr
ico T
e (
Nm
)
Par Eléctrico-Velocidad angular del rotor
8
Fig. 17. Gráficas de la corriente del rotor
Fig. 18. Gráficas de la corriente del estator
0 5 10 15 20 25 300
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2Corriente del rotor
Tiempo (s)
Corr
iente
ir
(A)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Corriente del rotor
Tiempo (s)
Corr
iente
ir
(A)
0 5 10 15 20 25 300
2
4
6
8
10
12
14
Tiempo (s)
Corr
iente
is (
A)
Corriente del estator
0.1 0.2 0.3 0.4 0.5 0.65
6
7
8
9
10
11
12
13
Tiempo (s)
Corr
iente
is (
A)
Corriente del estator
2 2.5 3 3.5 4
7
8
9
10
11
12
13
Tiempo (s)
Corr
iente
is (
A)
Corriente del estator
0 5 10 15 20 25 301020
1040
1060
1080
1100
1120
1140
Tiempo (s)
Veocid
ad (
rad/s
)
Velocidad angular del rotor
9
Fig. 19. Gráficas de la velocidad angular
Fig. 20. Gráfica de la posición angular del rotor
Fig. 21. Gráficas del par eléctrico.
Fig. 22. Gráficas del par eléctrico - velocidad.
6.375 6.38 6.385 6.39 6.395 6.4 6.405
1057.63
1057.64
1057.65
1057.66
1057.67
1057.68
1057.69
1057.7
1057.71
Tiempo (s)
Veocid
ad (
rad/s
)Velocidad angular del rotor
0 5 10 15 20 25 300
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000Posición angular del rotor
Posic
ión a
ngula
r (r
ad)
Tiempo (s)
0 5 10 15 20 25 30-1.5
-1
-0.5
0
0.5
1
1.5
Tiempo (s)
Par
elé
ctr
ico T
e (
Nm
)
Par Eléctrico
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5-1.5
-1
-0.5
0
0.5
1
1.5
Tiempo (s)
Par
elé
ctr
ico T
e (
Nm
)
Par Eléctrico
1020 1040 1060 1080 1100 1120 1140-1.5
-1
-0.5
0
0.5
1
1.5Par Eléctrico - Velocidad angular del rotor
Velocidad angular (rad/s)
Par
elé
ctr
ico (
Nm
)
1037.4 1037.6 1037.8 1038 1038.2 1038.4 1038.6 1038.8
-1
-0.5
0
0.5
1
1.5
Par Eléctrico - Velocidad angular del rotor
Velocidad angular (rad/s)
Par
elé
ctr
ico (
Nm
)
10
III. CONCLUSIONES Y OBSERVACIONES
Como puede observarse por ejemplo en las gráficas de las
corrientes, después de ese pequeño lapso de tiempo su
comportamiento es básicamente el mismo, es debido a eso se
hace un acercamiento en las gráficas para observar ese
comportamiento transitorio.
Además se consideró el factor de amortiguamiento [Norman
S. Nise, 2016], el cual puede observarse en comparación a
trabajos anteriores del curso de Máquinas Avanzadas I que las
gráficas tienden a ser un poco más estables y oscilar menos. En
este trabajo, en el caso B al aumentar esta constante, la
oscilación de la velocidad es básicamente nula a comparación
del caso A, al hacer el zoom en la gráfica, lo mismo sucede con
la gráfica de par-velocidad.
Se puede observar también que en el caso de alimentación en
C.D. la velocidad aumenta más rápido que en los casos de
alimentación con C.A.
IV. REFERENCIAS
[1] S. Simmons & D.O. Kelly, “IntroductiontoGeneralizedElectrical Machine
Theory”, McGrawHill1968.
[2] D. Olguín Salinas. “Apuntes de Máquinas Avanzadas I”.2012
[3] Applied Numerical Methods W/MATLAB: for Engineers & Scientists 3rd
Edition. Steven Chapra (2011)
[4] Control Systems Engineering 6th EditionNorman S. Nise (Author)
[5] Reporte de la ecuación de oscilación. Daniel Ruiz Vega.
[6] Power System Stability and Control. Kundur 1994. McGrawHill
[7]. Introduction to Programming with Fortran: With Coverage of Fortran 90,
95, 2003, 2008 and 77 3rd ed. 2015 Edition. Ian Chivers, Jane Sleightholme.
[8] Electromechanical Energy Conversion. David C. White, H.H. Woodson. 1
dic 1959.
[9] Reporte de la máquina universal. Julio Pedro V De Mata Castrejon. SEPI
ESIME ZACATENCO.
V. PROGRAMA EN FORTRAN 90
module precision ! ----------------------------------------------------- ! SP : simple precision de la norma IEEE 754 ! DP : doble precision de la norma IEEE 754 ! ----------------------------------------------------- integer, parameter :: sp = selected_real_kind(6,37) !simple precisión integer, parameter :: dp = selected_real_kind(15, 307) !doble precisión integer, parameter :: qp = selected_real_kind(33, 4931) !cuádruple precisión
end module precisión !**************************************************************************** ! ! MaquinasDoblementeExcitado.f90 ! INSTITUTO POLITÉCNICO NACIONAL ! SEPI-ESIME ZACATENCO ! PROFESOR: Daniel Olguín Salinas ! ALUMNO: Fernando Ramos Albarrán ! ! !**************************************************************************** program MaquinasDoblementeExcitado use precision implicit none real(qp), dimension(:), allocatable :: ir, is,ws,Os, t integer i, n real(qp) :: h,ti,tf !***************Tiempos y paso de integracion********************** print *, 'Tiempo inicial: ' read *, ti print *, 'Tiempo final:' read *, tf print *, 'Valor del paso de integración (h): ' read *, h n=(tf-ti)/h print *, 'n=',n allocate(ir(1:n+1),is(1:n+1),ws(1:n+1),Os(1:n+1),t(1:n+1)) t(1)=ti call procesorungekutta(n,ir,is,ws,Os,t,h) PAUSE end program MaquinasDoblementeExcitado subroutine procesorungekutta(n,ir,is,ws,Os,t,h) use precision implicit none real(qp) :: ir(1:n+1),is(1:n+1),ws(1:n+1),Os(1:n+1),t(1:n+1),Te(1:n+1) real(qp), dimension(6) :: irk,isk,wsk,Osk real(qp) :: h,ir0,is0,ws0,Os0,PRUEBA integer i, n real(qp) :: Vs, Vr,w0s,w0r,Lr,Ls1,Ls2,M,Rs,Rr,J,D,Tm,PdPolos,f,PI !************************ DATOS DEL SISTEMA ************************** Vr=15.0 !Magnitud de tensión en el rotor Vs=100.0 !Magnitud de tensión en el estator f=0.0 !Frecuencia (igual a cero en C.D.) PI=3.1416 !Relación entre longitud y radio del rotor
11
w0r= 2*PI*f ! Velocidad sincrona del rotor w0s= 2*PI*f ! Velocidad sincrona del estator Lr= 0.09474 !Inductancia del rotor Ls1= 0.09 !Inductancia en el estator Ls2= 0.01 !Inductancia en el estator M= 0.00136 !Inductancia mutua maxima Rs= 8.8713 !Resistencia del estator Rr= 8.2606 !Resistencia del rotor J= 0.5 !Inercia D=2.0 !Coeficiente de amortiguamiento Tm= 0.5 !Par mecánico PdPolos=1.0 !Pares de Polos !****************************************************************** print *, '***********************************************************' print *, ' Método de Runge Kutta 5to Orden Butcher ' print *, '***********************************************************' print *, 'valor inicial de ir(t): ' read *, ir0 print *, 'valor inicial de is(t): ' read *, is0 print *, 'valor inicial de wr(t) [3600 Valor de vel. propuesto]: ' !3600 Valor de vel. propuesto read *, ws0 print *,'wr(0)=',ws0 print *, 'valor inicial de 0r(t): ' read *, Os0 print *,'0r(0)=',Os0 is(1)=is0 print *,'ir(0)=',is(1) ir(1)=ir0 print *,'is(0)=',ir(1) ws(1)=ws0 print *,'wr(0)=',ws(1) Os(1)=Os0 print *,'0r(0)=',Os(1) print *,'t(0)=',t(1) print *,'h=',h PRUEBA= sin(0.0) print *,'SENO DE CERO', PRUEBA PAUSE 10 format (t8, 't', t22 , 'ir(t)', t36, 'is(t)',t50,'wr(t)',t64,'0r(t)',t78,'Te(t)') 20 format (t2, f14.8, t16, f14.8, t30, f14.8, t44, f14.8, t58, f14.8, t72, f14.8) do i=1,n !PRIMER VARIABLE DE ESTADO: Corriente del estator (is) isk(1)=( (-ws(i)*M*M*cos(Os(i))*sin(Os(i)) - Lr*Rs + 2*ws(i)*Lr*Ls2*sin(2*Os(i)) ) / ( (Ls1+Ls2*cos(2*Os(i)))*Lr-M*M*cos(Os(i))**2) )*is(i) &
+( ( M*cos(Os(i))*Rr + Lr*ws(i)*M*sin(Os(i)) ) / ( (Ls1+Ls2*cos(2*Os(i)))*Lr-M*M*cos(Os(i))**2) )*ir(i) & +( Lr / ( (Ls1+Ls2*cos(2*Os(i)))*Lr-M*M*cos(Os(i))**2) ) *Vs & !*sin(w0s*t(i)) & !ACTIVAR PARA CA +( -M*cos(Os(i)) /( (Ls1+Ls2*cos(2*Os(i)))*Lr -M*M*cos(Os(i))**2) )*Vr !*sin(w0r*t(i)) !ACTIVAR PARA CA isk(2)=( (-(ws(i)+(1/4)*h*isk(1))*M*M*cos(Os(i)+(1/4)*h*isk(1))*sin(Os(i)+(1/4)*h*isk(1)) - Lr*Rs + 2*(ws(i)+(1/4)*h*isk(1))*Lr*Ls2 & *sin(2*(Os(i)+(1/4)*h*isk(1))) ) / ( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*isk(1))))*Lr-M*M*cos(Os(i)+(1/4)*h*isk(1))**2) )*(is(i)+(1/4)*h*isk(1)) & +( ( M*cos(Os(i)+(1/4)*h*isk(1))*Rr + Lr*(ws(i)+(1/4)*h*isk(1))*M*sin(Os(i)+(1/4)*h*isk(1)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*isk(1))))*Lr-M*M*cos(Os(i)+(1/4)*h*isk(1))**2) )*(ir(i)+(1/4)*h*isk(1)) & +( Lr / ( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*isk(1))))*Lr-M*M*cos(Os(i)+(1/4)*h*isk(1))**2) ) *Vs & !*sin(w0s*(t(i)+(1/4)*h)) & !ACTIVAR PARA CA +( -M*cos(Os(i)+(1/4)*h*isk(1)) /( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*isk(1))))*Lr -M*M*cos(Os(i)+(1/4)*h*isk(1))**2) )*Vr !*sin(w0r*(t(i)+(1/4)*h)) !ACTIVAR PARA CA isk(3)=( (-(ws(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))*M*M*cos(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))*sin(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2)) & - Lr*Rs + 2*(ws(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))*Lr*Ls2 & *sin(2*(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))) ) / ( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))))*Lr & -M*M*cos(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))**2) )*(is(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2)) & +( ( M*cos(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))*Rr + Lr*(ws(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))*M*sin(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))))*Lr-
12
M*M*cos(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))**2) ) & *(ir(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2)) & +( Lr / ( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))))*Lr-M*M*cos(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))**2) ) & *Vs & !*sin(w0s*(t(i)+(1/4)*h)) & !ACTIVAR PARA CA +( -M*cos(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2)) /( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))))*Lr & -M*M*cos(Os(i)+(1/8)*h*isk(1)+(1/8)*h*isk(2))**2) )*Vr !*sin(w0r*(t(i)+(1/4)*h)) !ACTIVAR PARA CA isk(4)=( (-(ws(i)-(1/2)*h*isk(2)+h*isk(3))*M*M*cos(Os(i)-(1/2)*h*isk(2)+h*isk(3))*sin(Os(i)-(1/2)*h*isk(2)+h*isk(3)) & - Lr*Rs + 2*(ws(i)-(1/2)*h*isk(2)+h*isk(3))*Lr*Ls2 & *sin(2*(Os(i)-(1/2)*h*isk(2)+h*isk(3))) ) / ( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*isk(2)+h*isk(3))))*Lr & -M*M*cos(Os(i)-(1/2)*h*isk(2)+h*isk(3))**2) )*(is(i)-(1/2)*h*isk(2)+h*isk(3)) & +( ( M*cos(Os(i)-(1/2)*h*isk(2)+h*isk(3))*Rr + Lr*(ws(i)-(1/2)*h*isk(2)+h*isk(3))*M*sin(Os(i)-(1/2)*h*isk(2)+h*isk(3)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*isk(2)+h*isk(3))))*Lr-M*M*cos(Os(i)-(1/2)*h*isk(2)+h*isk(3))**2) ) & *(ir(i)-(1/2)*h*isk(2)+h*isk(3)) & +( Lr / ( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*isk(2)+h*isk(3))))*Lr-M*M*cos(Os(i)-(1/2)*h*isk(2)+h*isk(3))**2) ) & *Vs & !*sin(w0s*(t(i)+(1/2)*h)) & !ACTIVAR PARA CA +( -M*cos(Os(i)-(1/2)*h*isk(2)+h*isk(3)) /( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*isk(2)+h*isk(3))))*Lr & -M*M*cos(Os(i)-(1/2)*h*isk(2)+h*isk(3))**2) )*Vr !*sin(w0r*(t(i)+(1/2)*h)) !ACTIVAR PARA CA isk(5)=( (-(ws(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))*M*M*cos(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4)) & *sin(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4)) & - Lr*Rs + 2*(ws(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))*Lr*Ls2 &
*sin(2*(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))) ) / ( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))))*Lr & -M*M*cos(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))**2) )*(is(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4)) & +( ( M*cos(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))*Rr + Lr*(ws(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4)) & *M*sin(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))))*Lr-M*M*cos(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))**2) ) & *(ir(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4)) & +( Lr / ( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))))*Lr-M*M*cos(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))**2) ) & *Vs & !*sin(w0s*(t(i)+(3/4)*h)) & !ACTIVAR PARA CA +( -M*cos(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4)) /( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))))*Lr & -M*M*cos(Os(i)+(3/16)*h*isk(1)+(9/16)*h*isk(4))**2) )*Vr !*sin(w0r*(t(i)+(3/4)*h)) !ACTIVAR PARA CA isk(6)=( (-(ws(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) & *M*M*cos(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) & *sin(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) & - Lr*Rs + 2*(ws(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))*Lr*Ls2 & *sin(2*(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))) ) & / ( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))))*Lr & -M*M*cos(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))**2) ) &
13
*(is(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) & +( ( M*cos(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))*Rr & + Lr*(ws(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) & *M*sin(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))))*Lr & -M*M*cos(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))**2) ) & *(ir(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) & +( Lr / ( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))))*Lr & -M*M*cos(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))**2) ) & *Vs & !*sin(w0s*(t(i)+h)) & !ACTIVAR PARA CA +( -M*cos(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5)) & /( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))))*Lr & -M*M*cos(Os(i)-(3/7)*h*isk(1)+(2/7)*h*isk(2)+(12/7)*h*isk(3)-(12/7)*h*isk(4)+(8/7)*h*isk(5))**2) )*Vr !*sin(w0r*(t(i)+h)) !ACTIVAR PARA CA !SEGUNDA VARIABLE DE ESTADO: Corriente del rotor (ir) irk(1)=( (Rs*M*cos(Os(i)) - 2*ws(i)*M*cos(Os(i))*Ls2*sin(2*Os(i)) + ws(i)*(Ls1+Ls2*cos(2*Os(i)))*M*sin(Os(i)) ) & / ( (Ls1+Ls2*cos(2*Os(i)))*Lr-M*M*cos(Os(i))**2) )*is(i) & +( (-ws(i)*M*M*cos(Os(i))*sin(Os(i)) - Rr*(Ls1+Ls2*cos(2*Os(i))) ) / ( (Ls1+Ls2*cos(2*Os(i)))*Lr-M*M*cos(Os(i))**2) )*ir(i) & +( (-M*cos(Os(i)) )/( (Ls1+Ls2*cos(2*Os(i)))*Lr-M*M*cos(Os(i))**2) )*Vs & !*sin(w0s*t(i)) & !ACTIVAR PARA CA +( (Ls1+Ls2*cos(2*Os(i))) / ((Ls1+Ls2*cos(2*Os(i)))*Lr-M*M*cos(Os(i))**2) )*Vr !*sin(w0r*t(i)) !ACTIVAR PARA CA
irk(2)=( (Rs*M*cos(Os(i)+(1/4)*h*irk(1)) - 2*(ws(i)+(1/4)*h*irk(1))*M*cos(Os(i)+(1/4)*h*irk(1))*Ls2*sin(2*(Os(i)+(1/4)*h*irk(1))) & + (ws(i)+(1/4)*h*irk(1))*(Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*irk(1))))*M*sin(Os(i)+(1/4)*h*irk(1)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*irk(1))))*Lr-M*M*cos(Os(i)+(1/4)*h*irk(1))**2) )*(is(i)+(1/4)*h*irk(1)) & +( (-(ws(i)+(1/4)*h*irk(1))*M*M*cos(Os(i)+(1/4)*h*irk(1))*sin(Os(i)+(1/4)*h*irk(1)) - Rr*(Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*irk(1)))) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*irk(1))))*Lr-M*M*cos(Os(i)+(1/4)*h*irk(1))**2) )*(ir(i)+(1/4)*h*irk(1)) & +( (-M*cos(Os(i)+(1/4)*h*irk(1)) )/( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*irk(1))))*Lr-M*M*cos(Os(i)+(1/4)*h*irk(1))**2) ) & *Vs & !*sin(w0s*(t(i)+(1/4)*h)) & !ACTIVAR PARA CA +( (Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*irk(1)))) / ((Ls1+Ls2*cos(2*(Os(i)+(1/4)*h*irk(1))))*Lr-M*M*cos(Os(i)+(1/4)*h*irk(1))**2) ) & *Vr !*sin(w0r*(t(i)+(1/4)*h)) !ACTIVAR PARA CA irk(3)=( (Rs*M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)) - 2*(ws(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))*M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)) & *Ls2*sin(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))) & + (ws(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))*(Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))))*M*sin(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))))*Lr-M*M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))**2) ) & *(is(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)) & +( (-(ws(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))*M*M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))*sin(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)) & - Rr*(Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)))) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))))*Lr-
14
M*M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))**2) ) & *(ir(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)) & +( (-M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)) )/( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))))*Lr & -M*M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))**2) ) & *Vs & !*sin(w0s*(t(i)+(1/4)*h)) & !ACTIVAR PARA CA +( (Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2)))) / ((Ls1+Ls2*cos(2*(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))))*Lr & -M*M*cos(Os(i)+(1/8)*h*irk(1)+(1/8)*h*irk(2))**2) ) & *Vr !*sin(w0r*(t(i)+(1/4)*h)) !ACTIVAR PARA CA irk(4)=( (Rs*M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3)) - 2*(ws(i)-(1/2)*h*irk(2)+h*irk(3))*M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3)) & *Ls2*sin(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3))) & + (ws(i)-(1/2)*h*irk(2)+h*irk(3))*(Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3))))*M*sin(Os(i)-(1/2)*h*irk(2)+h*irk(3)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3))))*Lr-M*M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3))**2) ) & *(is(i)-(1/2)*h*irk(2)+h*irk(3)) & +( (-(ws(i)-(1/2)*h*irk(2)+h*irk(3))*M*M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3))*sin(Os(i)-(1/2)*h*irk(2)+h*irk(3)) & - Rr*(Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3)))) ) & / ( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3))))*Lr-M*M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3))**2) ) & *(ir(i)-(1/2)*h*irk(2)+h*irk(3)) & +( (-M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3)) )/( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3))))*Lr & -M*M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3))**2) ) & *Vs & !*sin(w0s*(t(i)+(1/2)*h)) & !ACTIVAR PARA CA +( (Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3)))) /
((Ls1+Ls2*cos(2*(Os(i)-(1/2)*h*irk(2)+h*irk(3))))*Lr & -M*M*cos(Os(i)-(1/2)*h*irk(2)+h*irk(3))**2) ) & *Vr !*sin(w0r*(t(i)+(1/2)*h)) !ACTIVAR PARA CA irk(5)=( (Rs*M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)) - 2*(ws(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))*M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)) & *Ls2*sin(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))) & + (ws(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))*(Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))))*M & *sin(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))))*Lr-M*M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))**2) ) & *(is(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)) & +( (-(ws(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))*M*M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))*sin(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)) & - Rr*(Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)))) ) & / ( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))))*Lr-M*M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))**2) ) & *(ir(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)) & +( (-M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)) )/( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))))*Lr & -M*M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))**2) ) & *Vs & !*sin(w0s*(t(i)+(3/4)*h)) & !ACTIVAR PARA CA +( (Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4)))) / ((Ls1+Ls2*cos(2*(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))))*Lr & -M*M*cos(Os(i)+(3/16)*h*irk(1)+(9/16)*h*irk(4))**2) ) & *Vr !*sin(w0r*(t(i)+(3/4)*h)) !ACTIVAR PARA CA
15
irk(6)=( (Rs*M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & - 2*(ws(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & *M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & *Ls2*sin(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))) & + (ws(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & *(Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))))*M & *sin(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) ) & / ( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))))*Lr & -M*M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))**2) ) & *(is(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & +( (-(ws(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))*M & *M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & *sin(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & - Rr*(Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)))) ) & / ( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))))*Lr & -M*M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))**2) ) & *(ir(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) & +( (-M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)) ) & /( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))))*Lr &
-M*M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))**2) ) & *Vs & !*sin(w0s*(t(i)+h)) & !ACTIVAR PARA CA +( (Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5)))) & / ((Ls1+Ls2*cos(2*(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))))*Lr & -M*M*cos(Os(i)-(3/7)*h*irk(1)+(2/7)*h*irk(2)+(12/7)*h*irk(3)-(12/7)*h*irk(4)+(8/7)*h*irk(5))**2) ) & *Vr !*sin(w0r*(t(i)+h)) !ACTIVAR PARA CA !TERCER VARIABLE DE ESTADO: Velocidad angular del rotor (wr) wsk(1)= (1/J)*(-Tm+PdPolos*(-M*sin(Os(i))*is(i)*ir(i)-Ls2*sin(2*Os(i))*is(i)**2)+D) wsk(2)= (1/J)*(-Tm+PdPolos*(-M*sin(Os(i)+(1/4)*h*wsk(1))*(is(i)+(1/4)*h*wsk(1))*(ir(i)+(1/4)*h*wsk(1)) & -Ls2*sin(2*(Os(i)+(1/4)*h*wsk(1)))*(is(i)+(1/4)*h*wsk(1))**2)+D) wsk(3)= (1/J)*(-Tm+PdPolos*(-M*sin(Os(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2)) & *(is(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2))*(ir(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2)) & -Ls2*sin(2*(Os(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2)))*(is(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2))**2) +D) wsk(4)= (1/J)*(-Tm+PdPolos*(-M*sin(Os(i)-(1/2)*h*wsk(2)+h*wsk(3)) & *(is(i)-(1/2)*h*wsk(2)+h*wsk(3))*(ir(i)-(1/2)*h*wsk(2)+h*wsk(3)) & -Ls2*sin(2*(Os(i)-(1/2)*h*wsk(2)+h*wsk(3)))*(is(i)-(1/2)*h*wsk(2)+h*wsk(3))**2) +D) wsk(5)= (1/J)*(-Tm+PdPolos*(-M*sin(Os(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4)) & *(is(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4))*(ir(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4)) & -Ls2*sin(2*(Os(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4)))*(is(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4))**2) +D) wsk(6)= (1/J)*(-Tm+PdPolos*(-M*sin(Os(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) &
16
*(is(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) & *(ir(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) & -Ls2*sin(2*(Os(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5))) & *(is(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5))**2) +D) !Esta sección se deja comentada para usar sin considerar valores en p.u. !wsk(1)= (1/J)*(-Tm*ws(i)+PdPolos*(-M*sin(Os(i))*is(i)*ir(i)-Ls2*sin(2*Os(i))*is(i)**2)*ws(i)+D*ws(i)) !wsk(2)= (1/J)*(-Tm*(ws(i)+(1/4)*h*wsk(1))+PdPolos*(-M*sin(Os(i)+(1/4)*h*wsk(1))*(is(i)+(1/4)*h*wsk(1))*(ir(i)+(1/4)*h*wsk(1)) & ! -Ls2*sin(2*(Os(i)+(1/4)*h*wsk(1)))*(is(i)+(1/4)*h*wsk(1))**2)*(ws(i)+(1/4)*h*wsk(1))+D*(ws(i)+(1/4)*h*wsk(1))) !wsk(3)= (1/J)*(-Tm*(ws(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2))+PdPolos*(-M*sin(Os(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2)) & ! *(is(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2))*(ir(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2)) & ! -Ls2*sin(2*(Os(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2)))*(is(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2))**2)*(ws(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2)) & ! +D*(ws(i)+(1/8)*h*wsk(1)+(1/8)*h*wsk(2))) !wsk(4)= (1/J)*(-Tm*(ws(i)-(1/2)*h*wsk(2)+h*wsk(3))+PdPolos*(-M*sin(Os(i)-(1/2)*h*wsk(2)+h*wsk(3)) & ! *(is(i)-(1/2)*h*wsk(2)+h*wsk(3))*(ir(i)-(1/2)*h*wsk(2)+h*wsk(3)) & ! -Ls2*sin(2*(Os(i)-(1/2)*h*wsk(2)+h*wsk(3)))*(is(i)-(1/2)*h*wsk(2)+h*wsk(3))**2)*(ws(i)-(1/2)*h*wsk(2)+h*wsk(3)) & ! +D*(ws(i)-(1/2)*h*wsk(2)+h*wsk(3))) !wsk(5)= (1/J)*(-Tm*(ws(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4))+PdPolos*(-M*sin(Os(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4)) & ! *(is(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4))*(ir(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4)) & ! -Ls2*sin(2*(Os(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4)
))*(is(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4))**2)*(ws(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4)) & ! +D*(ws(i)+(3/16)*h*wsk(1)+(9/16)*h*wsk(4))) !wsk(6)= (1/J)*(-Tm*(ws(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) & ! +PdPolos*(-M*sin(Os(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) & ! *(is(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) & ! *(ir(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) & ! -Ls2*sin(2*(Os(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5))) & ! *(is(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5))**2) & ! *(ws(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5)) & ! +D*(ws(i)-(3/7)*h*wsk(1)+(2/7)*h*wsk(2)+(12/7)*h*wsk(3)-(12/7)*h*wsk(4)+(8/7)*h*wsk(5))) !CUARTA VARIABLE DE ESTADO: Posición angular del rotor (0r) !Los términos comentados son para el modelo de la máquina sincrona Osk(1)= ws(i)-w0r Osk(2)= (ws(i)+(1/4)*h*Osk(1))-w0r Osk(3)= (ws(i)+(1/8)*h*Osk(1)+(1/8)*h*Osk(2))-w0r Osk(4)= (ws(i)-(1/2)*h*Osk(2)+h*Osk(3))-w0r Osk(5)= (ws(i)+(3/16)*h*Osk(1)+(9/16)*h*Osk(4))-w0r Osk(6)= (ws(i)-(7/3)*h*Osk(1)+(2/7)*h*Osk(2)+(12/7)*h*Osk(3)-(12/7)*h*Osk(4)+(8/7)*h*Osk(5))-w0r ir(i+1)=ir(i)+(7*irk(1)+32*irk(3)+12*irk(4)+32*irk(5)+7*irk(6))*(h/90) is(i+1)=is(i)+(7*isk(1)+32*isk(3)+12*isk(4)+32*isk(5)+7*isk(6))*(h/90) ws(i+1)=ws(i)+(7*wsk(1)+32*wsk(3)+12*wsk(4)+32*wsk(5)+7*wsk(6))*(h/90) Os(i+1)=Os(i)+(7*Osk(1)+32*Osk(3)+12*Osk(4)+32*Osk(5)+7*Osk(6))*(h/90) Te(i)=PdPolos*(-M*sin(Os(i))*is(i)*ir(i)-Ls2*sin(2*Os(i))*is(i)**2) t(i+1)=t(i)+h print *,'PASO',i,':' print *,'ir=',ir(i) print *,'is=',is(i) print *,'wr=',ws(i) print *,'0r=',Os(i)
17
end do open(7,file='rungekutta_log.txt') write(7,10) do i=1,n write(7,20) t(i),ir(i),is(i),ws(i),Os(i),Te(i) end do close(7) return end subroutine