Solucion de La Ecuación de Oscilación Mediante El Método de Euler y Euler Modificado

UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ FACULTAD DE INGENIERÍA ELÉCTRICA Y ELECTRÓNICA SOLUCIÓN NUMÉRICA DE LA ECUACIÓ

Views 222 Downloads 2 File size 2MB

Report DMCA / Copyright

DOWNLOAD FILE

Recommend stories

Citation preview

UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ FACULTAD DE INGENIERÍA ELÉCTRICA Y ELECTRÓNICA

SOLUCIÓN NUMÉRICA DE LA ECUACIÓN DE OSCILACIÓN MEDIANTE EL MÉTODO DE EULER Y EULER MODIFICADO CURSO:

ESTABILIDAD DE SISTEMAS DE POTENCIA

CATEDRÁTICO:

DR. ING. PEDRO NICOLÁS TORRES MAYTA

INTEGRANTES:

ARAOZ PALOMINO KLINSMAN HIDALGO HILARIO NILO KELVIN HINOSTROZA MEZA ÁNGEL CLIF MALLQUI APOLINARIO LENNIN SANDOVAL OBISPO ERIK OSORIO CHILQUILLO ELIZABETH POMA SUELDO DONATO RAMOS TORRES RYDER RAMOS YANGALI JHON SINCHE BARRA JOSELYN TORRES CLEMENTE POOL VELAZCO FLORES ANDRÉ

SEMESTRE:

X Huancayo – 2018

ÍNDICE INTRODUCCIÓN ............................................................................................................................. 2 I.

INTEGRACIÓN NUMÉRICA ................................................................................................. 3

II.

MÉTODO DE EULER Y EULER MODIFICADO ...................................................................... 7

III.

EJEMPLO DE APLICACIÓN ............................................................................................ 13

IV.

SOLUCIÓN MEDIANTE EL USO DE HERRAMIENTAS COMPUTACIONALES .................. 21

CONCLUSIONES ........................................................................................................................... 26

1

INTRODUCCIÓN En el presente trabajo se muestra las formas de solución de una ecuación de oscilación que van más allá de una solución explícita, con la solución de una formula dada y aunque es posible demostrar de forma abstracta que la mayoría de estas ecuaciones poseen solución, por lo menos de manera local, en general resulta muy difícil expresar explícitamente de que solución se trata. Principalmente se debe saber que las ecuaciones que se deben analizar no poseen solución de forma cerrada. Estos métodos son el método de Euler que, aunque presente poca precisión, permite hallar una solución aproximada, y el método de Euler Modificado cuya precisión es mayor, es decir está más próxima al valor verdadero, sin embargo, también cuenta con error. Se implementó una hoja de cálculo en MS Office Excel y un programa en Matlab para analizar la solución en diferentes puntos mostrando su aproximación al valor verdadero.

2

I.

INTEGRACIÓN NUMÉRICA

Para el análisis de la estabilidad transitoria de sistemas de potencia, se debe resolver un conjunto de ecuaciones diferenciales no lineales de 1er orden, con condiciones iniciales conocidas, de la forma: 𝑥̇ = 𝑓(𝑥, 𝑟, 𝑡) … (1) Donde: x= vector de variable de estado r= vector de variables algebraicas t= tiempo

Las variables algebraicas están relacionadas las restricciones algébricas de la forma: 0 = 𝑓(𝑥, 𝑟, 𝑡) … (2) Que deben ser satisfechas a cada instante de tiempo. En general, ecuaciones diferenciales no lineales sólo admiten solución numérica, calculada paso a paso. Diversos algoritmos de integración numérica son conocidos y la opción del “mejor” depende grandemente del problema. Siendo así, algunas características inherentes a los métodos de integración deben ser analizadas para llegar a una conclusión. Los principales métodos de integración numérica pueden ser interpretados como aproximaciones basadas en el truncamiento de la expresión de la solución del problema en términos de su serie de Taylor: 𝑥(𝑡0 + ∆𝑡) = 𝑥(𝑡0 ) +

(∆𝑡)2 𝑑 3 𝑥(𝑡) (∆𝑡)3 𝑑𝑥(𝑡) 𝑑 2 𝑥(𝑡) | ∆𝑡 + | + | + ⋯ (3) 𝑑𝑡 𝑡=𝑡 𝑑𝑡 2 𝑡=𝑡 2! 𝑑𝑡 3 𝑡=𝑡 3! 0

0

0

a) PRESICIÓN: Esta característica dada, principalmente, por dos factores: 

Error de redondeo



Error de truncamiento 3

Los errores de redondeos están relacionados a la representación de los números en un computador (aritmética finita) y pueden ser minimizados utilizando doble precisión y/u otros recursos inherentes al computador siendo utilizado para resolver el problema. Y los errores de truncamiento son debidos a la aproximación de la solución real del problema utilizada por el método de integración escogido. El error de truncamiento puede ser analizado a partir de la serie de Taylor de la función, mostrada en la ecuación (3), y será proporcional a: (∆𝑡)𝑝+1 … (4) Donde ∆t es el paso de integración escogido y p es el orden de la serie de Taylor utilizada por el método de integración como aproximación de la solución. b) ESTABILIDAD: Dos tipos de inestabilidad pueden ocurrir en la solución de ecuaciones diferenciales ordinarias: 

Inestabilidad inherente; y



Inestabilidad inducida.

La instabilidad inherente surge cuando errores numéricos son amplificados, a cada paso de la solución, hasta dominar completamente el cálculo y hacer el método divergencia de la solución real. La inestabilidad inducida está relacionada con el método de integración utilizado o, de manera más precisa, con la discretización resultante de la aplicación del método, que depende, también, del paso de integración escogido. Se puede mostrar que los métodos de integración corresponden a resolver una ecuación algebraica de la forma ∝0 𝑦𝑛+1 +∝1 𝑦𝑛 +∝2 𝑦𝑛−1 + ⋯ +∝𝑘 𝑦𝑛−𝑘+1 + 𝛽0

𝑑𝑦𝑛+1 𝑑𝑦𝑛 𝑑𝑦𝑛−𝑘+1 + 𝛽1 + ⋯ + 𝛽𝑘 𝑑𝑡 𝑑𝑡 𝑑𝑡

= 0 … (6) Donde αi y βi son constantes y se desea determinar yn+1, conociendo los valores anteriores de yi (i < n+1).

4

Esta es una ecuación discreta, lineal, y será estable si todas las raíces de su polinomio característico tuvieran módulo menor que 1. El mayor problema está en el hecho que cuanto mayor es la precisión del método, menor será su estabilidad. Esto es más grave para los métodos en que el método o paso de integración es crítico para la solución, en especial los métodos de Runge-Kutta. Estos métodos son estables sólo para algunos valores de ∆t y son dichos condicionalmente estables. Métodos con margen de estabilidad infinita (independiente de ∆t) son llamados métodos estables y los métodos de Euler Reverso y Trapezoidal están en esta categoría. c) RIGIDEZ: Un sistema de ecuaciones diferenciales es dicho rígido (“stiff”) cuando la razón entre la mayor y la menor constantes de tiempo del problema sea mucho mayor que 1. Métodos de integración tradicionales tienen su intervalo de integración definido por la menor constante de tiempo y el tiempo final de la simulación es definido por la mayor constante de tiempo. De esta forma, un gran número de puntos deberán ser determinados, aumentando sobremanera el tiempo de computación. Además de eso, si el método de integración no fuera estable, puede haber problemas de inestabilidad numérica. 1. MÉTODOS DE INTEGRACIÓN NUMÉRICA: Los métodos de integración numérica, como fue visto en la sección anterior, representan una discretización de la solución, de forma que la ecuación diferencial original (continua) es aproximada por una ecuación de diferencias (discreta) y apenas algunos valores (correspondientes a determinados instantes de tiempo) son calculados, que corresponderán a la solución aproximada obtenida. La Figura 1 ilustra este efecto, en que una función es aproximada a partir de valores discretos.

5

Fig. 1: Discretización de una Función Continua

Para ejemplificar la aplicación de los diversos métodos de integración numérica, será utilizado el siguiente sistema lineal de las ecuaciones diferenciales: 𝑥̇ = 𝑥2 { 1 𝑥̇ 2 = −𝑎1 𝑥1 − 𝑎2 𝑥2

𝑥1 (0) = 𝑥10 … (7) 𝑥2 (0) = 𝑥20

cuya solución analítica es de la forma ℷ2 𝑥10 − 𝑥20 ℷ 𝑡 𝑥20 − ℷ1 𝑥10 ℷ 𝑡 𝑒1 + 𝑒2 ℷ2 − ℷ1 ℷ2 − ℷ1 … (8) ℷ2 𝑥10 − 𝑥20 ℷ 𝑡 𝑥20 − ℷ1 𝑥10 ℷ 𝑡 1 2 𝑥 (𝑡) = ℷ1 𝑒 + ℷ2 𝑒 { 1 ℷ2 − ℷ1 ℷ2 − ℷ1 𝑥1 (𝑡) =

Donde ℷ1,2 =

−𝑎2 ± √𝑎2 2 − 4𝑎1 … (9) 2

Son las raíces del polinomio característico de la ecuación diferencial de 2do orden asociada a la ecuación (7) o, equivalentemente, son los polos de la función de transferencia de este sistema.

6

II.

MÉTODO DE EULER Y EULER MODIFICADO

1. MÉTODO DE EULER: El método de Euler es el más simple de los métodos de integración numérica y puede ser visualizado en la Figura 2. Utilizando el valor de la derivada de la función en el instante de tiempo t = t0, se puede escribir que: 𝑥(𝑡0 + ∆𝑡) = 𝑥(𝑡0 ) +

𝑑𝑥(𝑡) | … (10) 𝑑𝑡 𝑡=𝑡0

Fig.2: Interpretación Gráfica del Método de Euler

Este método corresponde a la aplicación de la serie de Taylor, mostrada en la ecuación (1), aproximada sólo por los dos primeros términos. De esta forma, la precisión de este método es del orden O(∆t)2. La discretización del sistema de ecuaciones, dada por el método de Euler, es equivalente a la ecuación (6) haciendo: {

∝0 = 1 ∝1 = −1 𝛼2 = 𝛼3 = ⋯ = 𝛼𝑘 = 0

𝛽0 = 0 𝛽1 = −∆𝑡 … … … (11) 𝛽2 = 𝛽3 = ⋯ = 𝛽𝑘 = 0

La aplicación del método de Euler para la solución del sistema dado por la ecuación (7) resulta en:

7

𝑑𝑥1 (𝑥1 )𝑛+1 = (𝑥1 )𝑛 + ( | ) ∆t 𝑑𝑡 𝑛 … (12) 𝑑𝑥2 (𝑥 ) = (𝑥2 )𝑛 + ( | ) ∆t { 2 𝑛+1 𝑑𝑡 𝑛 El método de Euler puede ser implementado de forma explícita o implícita. En su forma explícita, el valor de las derivadas mostradas en la ecuación (12) son explícitamente calculadas a cada paso del algoritmo. Ya en el caso de haber expresiones analíticas para las derivadas, estas pueden ser sustituidas, resultando un sistema de ecuaciones que no dependen explícitamente de las derivadas. Para el ejemplo dado, la ecuación (7) suministra expresiones analíticas para las derivadas, resultando en el siguiente sistema de ecuaciones: {

𝑥̇ 1 = 𝑥2 𝑥̇ 2 = −𝑎1 𝑥1 − 𝑎2 𝑥2 =[

𝑥1 (0) = 𝑥10 (𝑥 ) ⇒ [ 1 𝑛+1 ] (𝑥2 )𝑛+1 𝑥2 (0) = 𝑥20

1 −𝑎1 ∆𝑡

(𝑥 ) ∆𝑡 ] [ 1 𝑛 ] … (13) 1 − 𝑎2 ∆𝑡 (𝑥2 )𝑛

La formulación implícita permite reducir el error numérico que puede haber en la determinación de las derivadas de las funciones, que son sustituidas por expresiones analíticas equivalentes. Para el caso particular de la ecuación de oscilación se conoce que: 𝑋𝑘 = [

𝑋′𝐾

𝛿(𝑡𝑘 ) 𝑥1 𝑘 𝑘 ] = [𝜔(𝑡 )] 𝑥2 𝑘

𝜕(𝛿(𝑡𝑘 )) 𝑥 = [ 1 ′ 𝑘 ] = [ 𝜕𝑡 ] 𝜕(𝜔(𝑡𝑘 )) 𝑥2 𝜕𝑡 ′𝑘

𝜔(𝑡𝑘 ) 𝜔(𝑡𝑘 ) 𝐹(𝑋𝐾 , 𝑡𝐾 ) = [𝜕 (𝛿(𝑡𝑘 ))] = [𝜋𝑓(𝑃𝑚𝑒𝑐 + 𝑃𝑚𝑎𝑥 sin 𝛿(𝑡𝑘 ))] 𝐻 𝜕 2𝑡 2

Sustituyendo los respectivos términos:

8

𝜔(𝑡𝑘 ) 𝛿(𝑡𝑘+1 ) 𝛿(𝑡𝑘 ) [ ]=[ ] + ∆𝑡 [𝜋𝑓(𝑃𝑚𝑒𝑐 + 𝑃𝑚𝑎𝑥 sin 𝛿(𝑡𝑘 ))] 𝜔(𝑡𝑘+1 ) 𝜔(𝑡𝑘 ) 𝐻 En forma de ecuaciones se tiene: 𝛿(𝑡𝑘+1 ) = 𝛿(𝑡𝑘 ) + ∆𝑡 ∗ 𝜔(𝑡𝑘 ) 𝜔(𝑡𝑘+1 ) = 𝜔(𝑡𝑘 ) + ∆𝑡 ∗ (

𝜋𝑓(𝑃𝑚𝑒𝑐 + 𝑃𝑚𝑎𝑥 sin 𝛿(𝑡𝑘 )) ) 𝐻

Se puede suponer: 𝛿(𝑡𝑘 ) = 𝛿𝑘 𝜔(𝑡𝑘 )= 𝜔𝑘 De tal forma que se puede escribir en forma más compacta las ecuaciones: 𝛿𝑘+1 = 𝛿𝑘 + ∆𝑡 ∗ 𝜔𝑘 𝜔𝑘+1 = 𝜔𝑘 + ∆𝑡 ∗ (

𝜋𝑓(𝑃𝑚𝑒𝑐 + 𝑃𝑚𝑎𝑥 sin 𝛿𝑘 ) 𝐻

2. MÉTODO DE EULER MODIFICADO: El método de Euler modificado intenta mejorar el desempeño del método convencional a través de la mejoría de la estimativa de la derivada. Para eso, la ecuación (10) será modificada haciendo la derivada igual a la media entre la derivada en el instante t = t0 y el valor de la derivada calculada para una estimativa del estado en el instante t = t0+ ∆t. Este es los más simples de los 9

métodos de integración del tipo predictor-corrector, en que surge un paso de predicción del valor del estado en el instante t = t0+ ∆t y este valor es utilizado para la corrección y obtención de un nuevo valor del estado en ese instante. El método, por lo tanto, puede ser descrito por los siguientes pasos: Paso 1: Predictor

…(18) Paso 2: Corrector

…(19) Se debe notar que la ecuación (18) corresponde a un paso del método de Euler convencional, mostrado en la ecuación (10). Ya la ecuación (19) es muy parecida con el método trapezoidal, siendo la única diferencia el cálculo aproximado (a partir de la estimativa xp) de la derivada en el instante t = t0+ ∆t. La implementación de este método para la solución del sistema ejemplo resulta en:

…(20) Y, para el paso corrector:

…(21) Este método también puede ser implementado de forma implícita, sustituyendo las expresiones de las derivadas en las ecuaciones (20) y (21):

10

…(22)

…(23) Se puede, aún, sustituir los valores de x1p y x2p, obtenidos en las ecuaciones (22), en la ecuación (23), resultando en:

…(24) Aplicando este método a la solución a la ecuación de oscilación, se obtiene el siguiente procedimiento de solución: a) Cálculo de la potencia acelerante y las derivadas del ángulo y de la velocidad al principio del intervalo, es decir, en el instante t:

…(25) b) Cálculo del primer valor del ángulo y la velocidad al final del intervalo, es decir, en el instante 𝑡 + ∆𝑡:

…(26) c) Cálculo de la potencia acelerante y de la derivada de la velocidad al final del intervalo, es decir, en el instante 𝑡 + ∆𝑡: 11

…(27) d) Cálculo del valor corregido de la velocidad al final del intervalo, es decir, en el instante 𝑡 + ∆𝑡:

…(28) e) Cálculo de la derivada del ángulo al final del intervalo, es decir, en el instante 𝑡 + ∆𝑡:

…(29) f) Cálculo del valor corregido del ángulo al final del intervalo, es decir, en el instante 𝑡 + ∆𝑡:

…(30) El proceso iterativo se continúa hasta que se presenten nuevos cambios en el sistema o hasta que se concluya el tiempo de estudio.

12

III.

EJEMPLO DE APLICACIÓN

Un generador sincrónico a 60 Hz, posee una constante de inercia H igual a 5MJ/MVA y una reactancia transitoria de eje directo X’d = 0.3 p.u, está conectado a una barra de potencia infinita a través de un circuito puramente reactivo como se muestra en la siguiente figura:

El generador entrega una potencia real de Pelec = 0.8 p.u, y Qelec = 0.074 p.u a la barra de potencia infinita cuyo voltaje es 1.0 p.u. Una falla trifásica ocurre en la mitad de una línea, y la cual es despejada por la puesta fuera de servicio de la línea, mediante la apertura simultanea de ambos extremos de la línea. La falla es despejada en 0.3 segundos. Obtener la solución numérica de la ecuación de oscilación, empleando el método de Euler, con un paso de tiempo Δt = 0.05 segundos. Efectúe la simulación numérica hasta t = 0.5. a) En Pre falla: La impedancia de transferencia entre el voltaje generador y la barra de potencia infinita Antes de que ocurra la falla, la reactancia entre los puntos A y B puede ser encontrada por la combinación serie paralelo:

𝑋1 = 0.30 + 0.2 +

0.30 = 0.65 𝑝. 𝑢. 2 13

La corriente se calcula:

El voltaje interno de la maquina E’ es dado por:

La ecuación de potencia eléctrica antes de la perturbación queda dada por: 𝑃𝑒𝑙𝑒𝑐 =

|𝐸 ´ ||𝐸𝐵 | 1 ∗ 1.17 𝑠𝑒𝑛𝛿 = 𝑠𝑒𝑛𝛿 𝑋1 0.65 𝑃𝑒𝑙𝑒𝑐 = 1.80 ∗ 𝑠𝑒𝑛𝛿

Siendo el ángulo de operación estable inicial: 0 𝑃𝑒𝑙𝑒𝑐 = 0.8 = 1.80 ∗ 𝑠𝑒𝑛𝛿0

𝛿0 = 26.38780 = 0.46055 𝑟𝑎𝑑

b) En Falla: Durante la falla se tiene, que la línea fallada se divide en dos, resultando el diagrama de impedancias como en la figura siguiente:

Se efectúa una transformación de estrella a delta y se logra:

14

La reactancia equivalente entre el generador y la barra de potencia infinita resulta ser:

La ecuación de potencia eléctrica antes de la perturbación queda dada por:

c) En Post falla: Cuando la falla es despejada, la línea fallada es puesta fuera de servicio, resultando el diagrama de reactancias como sigue:

Donde la reactancia equivalente entre el generador y la barra de potencia infinita resulta ser: De tal modo que la potencia eléctrica después de la falla queda dada por:

Se procede a plantear las ecuaciones de oscilación, de la forma:

ANTES de la perturbación resulta:

15

DURANTE la perturbación resulta:

DESPUÉS la perturbación resulta:

Considere que la falla es despejada en 0.3 segundos. Obtener la solución numérica de la ecuación de oscilación, empleando el método de Euler, con un paso de tiempo Δt = 0.05 segundos. Efectúe la simulación numérica hasta t = 0.5. Las ecuaciones de oscilaciones a considerar resultan ser:

1. MÉTODO DE EULER:

[



𝜔(𝑡𝑘 ) 𝛿(𝑡𝑘+1 ) 𝛿(𝑡 ) ] = [ 𝑘 ] + ∆𝑡 [𝜋𝑓(𝑃𝑚𝑒𝑐 − 𝑃𝑚𝑎𝑥 sin 𝛿(𝑡𝑘 ))] 𝜔(𝑡𝑘+1 ) 𝜔(𝑡𝑘 ) 𝐻

Para un tiempo de t=0: 𝛿(0)=26.38 ; 𝜔(0)=0



Para un tiempo t=0.05: [

0 𝛿(1) 𝛿(0) ]=[ ] + 0.05 [𝜋 ∗ 60 ∗ (0.8 − 𝑃𝑚𝑎𝑥𝐼𝐼 sin 𝛿(26.38))] 𝜔(1) 𝜔(0) 5

0 𝛿(1) 𝛿(1) 26.38 26.38 [ ]=[ ] + 0.05 [𝜋∗60∗(0.8−0.65 sin 𝛿(26.38))]; [ ]=[ ] 𝜔(1) 𝜔(1) 0 0.9638 5



Para un tiempo t=0.1: [

[

𝜔(1) 𝛿(2) 𝛿(1) ]=[ ] + 0.05 [𝜋 ∗ 60 ∗ (0.8 − 𝑃𝑚𝑎𝑥𝐼𝐼 sin 𝛿(26.38))] 𝜔(2) 𝜔(1) 5

0.963 ∗ 180/𝑝𝑖 𝛿(2) 𝛿(2) 26.38 29.147 ]=[ ] + 0.05 [𝜋∗60∗(0.8−0.65 sin 𝛿(26.38))]; [ ]=[ ] 𝜔(2) 𝜔(2) 0.963 1.926 5

16



Para un tiempo t=0.15: [

[

𝜔(2) 𝛿(3) 𝛿(2) ]=[ ] + 0.05 [𝜋 ∗ 60 ∗ (0.8 − 𝑃𝑚𝑎𝑥𝐼𝐼 sin 𝛿(26.38))] 𝜔(3) 𝜔(2) 5

1.926 ∗ 180/𝑝𝑖 𝛿(3) 𝛿(3) 29.147 34.665 ]=[ ] + 0.05 [𝜋∗60∗(0.8−0.65 sin 𝛿(29.147))]; [ ]=[ ] 𝜔(3) 𝜔(3) 1.926 2.837 5



Para un tiempo t=0.2: [

[

𝜔(3) 𝛿(4) 𝛿(3) ]=[ ] + 0.05 [𝜋 ∗ 60 ∗ (0.8 − 𝑃𝑚𝑎𝑥𝐼𝐼 sin 𝛿(34.665))] 𝜔(4) 𝜔(3) 5

2.837 ∗ 180/𝑝𝑖 𝛿(4) 𝛿(4) 42.792 34.665 ]=[ ] + 0.05 [𝜋∗60∗(0.8−0.65 sin 𝛿(34.665))]; [ ]=[ ] 𝜔(4) 𝜔(4) 3.648 2.837 5



Para un tiempo t=0.25: [

[

𝜔(4) 𝛿(5) 𝛿(4) ]=[ ] + 0.05 [𝜋 ∗ 60 ∗ (0.8 − 𝑃𝑚𝑎𝑥𝐼𝐼 sin 𝛿(42.792))] 𝜔(5) 𝜔(4) 5

3.648 ∗ 180/𝑝𝑖 𝛿(5) 𝛿(5) 42.792 53.243 ]=[ ] + 0.05 [𝜋∗60∗(0.8−0.65 sin 𝛿(42.792))]; [ ]=[ ] 𝜔(5) 𝜔(5) 3.648 4.324 5



Para un tiempo t=0.3 (despeje de la falla): [

𝜔(5) 𝛿(6) 𝛿(5) ]=[ ] + 0.05 [𝜋 ∗ 60 ∗ (0.8 − 𝑃𝑚𝑎𝑥𝐼𝐼𝐼 sin 𝛿(5))] 𝜔(6) 𝜔(5) 5

4.325 ∗ 180/𝑝𝑖 𝛿(6) 𝛿(6) 53.252 65.641 [ ]=[ ] + 0.05 [𝜋∗60∗(0.8−1.4625 sin 𝛿(53.252))]; [ ]=[ ] 𝜔(6) 𝜔(6) 4.325 3.624 5

Este método se repite hasta poder observar la estabilidad o inestabilidad del sistema mediante la curva angula versus tiempo. Cuadro de resultados: t (s)

Pm

Pe

𝛿 (sexag.)

𝜔 (rad/s)

0.00

0.8

1.8000

26.38

0

0.05

0.8

0.6500

26.38

0.9638

0.10

0.8

0.6500

29.147

1.926

0.15

0.8

0.6500

34.665

2.837

0.20

0.8

0.6500

42.792

3.648

0.25

0.8

0.6500

53.243

4.324 17

0.30

0.8

1.4625

65.641

3.624

0.35

0.8

1.4625

76.022

2.62

0.40

0.8

1.4625

83.529

1.453

0.45

0.8

1.4625

87.692

0.222

0.50

0.8

1.4625

88.327

-1.024

……

……

…….

…….

………

2. MÉTODO DE EULER MODIFICADO: Considerando los cálculos previos y las condiciones iniciales de estabilidad: 𝛿 = 0.46 rad 𝜔0= 0 rad/s ∆𝑡= 0.1; tdespeje= 0.3s 

Primera iteración: t1=t0+∆𝑡= 0+0.1=0.1s Hallando parámetros previos: 𝜕(𝛿(𝑡0 )) 𝜕𝑡 𝜕(𝜔(𝑡0 )) 𝜕𝑡

=0- 377= -377 =

л x 60 5

(0.8-0.65 sen 0.46) = 19.27 rad/s²

𝛿1 (0+0.1) = 0.46 + (-377) (0.1) = -37.24 rad 𝜔1(0+0.1) = 0 + (19.27) (0.1) = 1.927 rad/s 𝜕(𝛿(0+0.1)) 𝜕𝑡 𝜕(𝜔(𝑡0 )) 𝜕𝑡

=

= 1.927 – 377 = -377.093 л x 60 5

(0.8-0.65 sen -37.24) = 19.3 rad/s²

𝛿1 (0+0.1) = 0.46 +

(−377)+(−377.073) 2

𝛿1 (0.1) = -37.14rad 𝜔1(0+0.1) = 0 +

19.27+19.3 2

x0.1

𝛿1 (0.1) = 31.83º x 0.1

𝜔1(0.1) = 1.9285 rad/s 18



Segunda iteración: t2=t1+∆𝑡= 0.1+0.1=0.2s Hallando parámetros previos: 𝜕(𝛿2(0.1)) 𝜕𝑡 𝜕(𝜔2(0.1)) 𝜕𝑡

=1.9285 - 377 = -375.07 =

л x 60 5

(0.8-0.65 sen -37.14) = 17.16 rad/s²

𝛿2 (0.1+0.1) = -37.14 + (-375.07) (0.1) = -74.65 rad 𝜔2 (0.1+0.1) = 1.9285 + (17.16) (0.1) = 3.64 rad/s 𝜕(𝛿2(0.1+0.1)) 𝜕𝑡 𝜕(𝜔2(0.1+0.1)) 𝜕𝑡

= 3.64 – 377 = -373.36 =

л x 60 5

(0.8-0.65 sen -74.65) = 13.4 rad/s²

𝛿2 (0.1+0.1) = -37.14 +

(−375.07)+(−373.36) 2

𝛿2 (0.2) = -74.5615rad

x0.1

𝛿2 (0.2) = 47.94º 17.16+13.488

𝜔2 (0.1+0.1) = 1.9285 +

2

x 0.1

𝜔2 (0.2) = 3.46 rad/s 

Tercera iteración: t3=t2+∆𝑡= 0.2+0.1=0.3s Hallando parámetros previos: 𝜕(𝛿3(0.2)) 𝜕𝑡 𝜕(𝜔3(0.2)) 𝜕𝑡

=3.46 - 377 = -373.54 =

л x 60 5

(0.8-0.65 sen -74.5615) = 11.97 rad/s²

𝛿3 (0.2+0.1) = -74.5615 + (-373.54) (0.1) = 1.18rad 𝜔3 (0.2+0.1) = 3.46 + (11.97) (0.1) = 4.657 rad/s 𝜕(𝛿2(0.2+0.1)) 𝜕𝑡 𝜕(𝜔3(0.2+0.1)) 𝜕𝑡

= 4.657 – 377 = -372.343 =

л x 60 5

(0.8-1.4625 sen 1.18) = -20.8 rad/s²

19

𝛿3 (0.2+0.1) = 0.837 +

(−373.54)+(−372.343) 2

𝛿3 (0.3) = -36.46rad 𝜔3 (0.2+0.1) = 3.46 +

x0.1

𝛿3 (0.3) = 71º 11.97+(−20.82) 2

𝛿3 (0.3) = 1.24rad

x 0.1

𝜔3 (0.3) = 3.0175 rad/s Finalmente se obtiene el cuadro de resultados: Iteración 𝛿 (sexag.) ꞷ (rad/s) 1

31.83

1.9285

2

47.94

3.46

3

71

3.0175

20

IV.

SOLUCIÓN MEDIANTE EL USO DE HERRAMIENTAS COMPUTACIONALES

1. CON EL MS OFFICE EXCEL: En la siguiente tabla se muestra los resultados transpuestos, cuando la falla fuera sostenida (sin aclaramiento), y cuando fuera aclarada en los 0.3 segundos como se suponía al principio. Iteración (i) ti (s)

ᵟi(rad)

ᵟi(sexagesimal)

wi (rad/s)

0

0

0.46055

26.3875713

0

1

0.05

0.46055

26.3875713

0.96342613

2

0.1 0.50872131

29.1475838

1.92685226

3

0.15 0.60506392

34.6676089

2.83805945

4

0.2 0.74696689

42.7980503

3.6491002

5

0.25 0.9294219

53.2519523

4.32462941

6

0.3

1.1456534

65.641103

4.8508558

7

0.35 1.3881962

79.537781

5.2426686

8

0.4

1.6503296

94.55692

5.54578141

9

0.45 1.9276187

110.44441

5.8323978

10

0.5

127.153

6.19231597

2.2192386

Ángulo vs. Tiempo 4000 3500 3000 2500 2000 1500 1000 500 0 -500 0

0,5

1

1,5 Falla despejada

2

2,5

3

3,5

Fala sostenida

21

En la figura se visualiza que el ángulo crece indefinidamente cuando la falla persiste en el sistema. La siguiente tabla muestra la evolución del ángulo a través del tiempo hasta los 1.75 segundos una vez se ha despejado la falla a los 0.3 segundos. Paso tiempo (s) Iteración (i) ti (s)

ᵟi(rad)

ᵟi(sexages.)

wi (rad/s) 0

0.05

0

0

0.46055

26.3875713

0.05

1

0.05

0.46055

26.3875713 0.96342613

0.05

2

0.1

0.05

3

0.15 0.60506392 34.6676089 2.83805945

0.05

4

0.2

0.74696689 42.7980503

0.05

5

0.25

0.9294219

53.2519523 4.32462941

0.05

6

0.3

1.14565337

65.641103

4.8508558

0.05

7

0.35 1.38819616

79.537781

3.847479

0.05

8

0.4

1.58057011

90.559997

2.64452723

0.05

9

0.45 1.71279647

98.136009

1.39587582

0.05

10

0.5

102.1349

0.17483965

0.05

11

0.55 1.79133225

102.63578 -1.01234494

0.05

12

0.6

99.735623 -2.19436058

0.05

13

0.65 1.63099697

93.449243 -3.40344242

0.05

14

0.7

1.46082485

83.699099 -4.64723162

0.05

15

0.75 1.22846327

70.385761 -5.87936182

0.05

16

0.8

0.93449518

53.54263

-6.96818202

0.05

17

0.85 0.58608608

33.580259

-7.6774674

0.05

18

0.9

11.585935

-7.6942725

0.05

19

0.95 -0.18250092 -10.45653 -6.73996616

0.05

20

0.05

21

1.05 -0.75608327 -43.32038 -1.85514151

0.05

22

1.1 -0.84884034 -48.63497

1.54416455

0.05

23

1.15 -0.77163212 -44.21126

5.12110817

1

0.50872131 29.1475838 1.92685226

1.78259026

1.740715

0.20221271

3.6491002

-0.51949923 -29.76511 -4.73168088

22

0.05

24

1.2 -0.51557671 -29.54037

8.55136932

0.05

25

1.25 -0.08800824 -5.042501

11.4185114

0.05

26

1.3

0.48291733

27.669125

13.1687793

0.05

27

1.35

1.1413563

65.394899

13.3966071

0.05

28

1.4

1.81118665

103.77335

12.3981394

0.05

29

1.45 2.43109362

139.2914

11.2286262

0.05

30

1.5

2.99252493

171.45905

10.9386065

0.05

31

1.55 3.53945526

202.79585

12.0371491

0.05

32

1.6

4.14131271

237.27974

14.6132121

0.05

33

1.65 4.87197331

279.14351

18.4404826

0.05

34

1.7

5.79399744

331.9716

22.6701658

0.05

35

1.75 6.92750573

396.91684

25.4735512

Ángulo vs. Tiempo 4000 3500 3000 2500 2000 1500 1000 500 0 -500 0

0,5

1

1,5

2

2,5

3

3,5

Falla despejada

En la figura se alcanza notar que el ángulo oscila entre valores positivos y negativos. Pero según transcurre el tiempo, el ángulo crecerá indefinidamente. Lo cual muestra un régimen de inestabilidad, a pesar del despeje de falla a los 0.3 segundos.

23

2. CON MATLAB:

El ingreso de datos arroja los siguientes resultados:

24

El sistema para un tiempo de despeje de 0.3 segundos es inestable:

25

CONCLUSIONES 1. Los métodos de Euler y Euler modificado son herramientas eficaces para la solución numérica de la ecuación de oscilación de un sistema, pudiéndose calcular el ángulo y la velocidad en función del tiempo; y de esta manera predecir la estabilidad del sistema en análisis. Para obtener resultados con mayor precisión se debe elegir un Δt pequeño, aunque incrementa el número y el tiempo de cálculo; éste inconveniente se solucionó con la implementación de una hoja de cálculo en MS Office Excel y un programa en Matlab, ambos permiten graficar la variación de la velocidad y el ángulo en función del tiempo. 2. En el ejercicio desarrollado se hace clara la inestabilidad del sistema a pesar de ser liberada la falla y también como era de esperarse la inestabilidad se hace más grave en el caso de falla sostenida.

26