Las denominadas ecuaciones de Euler-Lagrange describen la evolución de un sistema mecánico sujeto a restricciones holonómicas (los diferentes grados están desacoplados). Con el fin de determinar las ecuaciones de Euler-Lagrange en una situación específica, se define las energías cinéticas y potenciales del sistema ^[14].

El modelo que se expone a continuación se basa en la descripción de la cinemática de un aerogenerador tripala desarrollada en ^[11] y la forma algorítmica Euler Lagrange presentada en ^[15] que usa la representación D-H. Para obtener el modelo, se realiza una primera aproximación considerando una sola pala del aerogenerador, una vez obtenido este modelo se generaliza la solución para las tres palas en base a las relaciones geométricas existentes entre las palas.

La obtención del modelo para una pala se describe a continuación.

1. 1. Se asignan a cada eslabón un sistema de referencia de acuerdo con la metodología D-H, en conformidad con el modelo expuesto en ^[11], tal como se muestra en la Figura 1.
2. 2. Obtenidos los parámetros D-H se definen las matrices de transformación respectivas, excluyendo del análisis el movimiento de la torre (θ1 = 0). De esta forma, serán el giro del eje del rotor (θ2) y el ángulo de paso de la pala (θ3), los que conformen el conjunto de variables generalizadas {q2 ,q3} que se emplean para describir la dinámica, de modo tal que,

donde: 𝑠𝑖=sen(𝑞𝑖), 𝑐𝑖=cos(𝑞𝑖), para i = 1, 2, 3; y c es la distancia a lo largo del eje x₃ que va desde la intersección del eje z₂ con el eje x₃ hasta el origen del sistema 3, se calcula como 𝑐=√𝐴2+𝐵2.

3. Enseguida se calculan las matrices 𝐔_𝑖𝑗 definidas por 𝐔_𝑖𝑗= , donde i, j = 1, 2, …, n; siendo 𝑞𝑗 la variable generalizada j. Así

4. Se obtienen las matrices 𝐔_{𝑖𝑗,𝑘} definidas por , donde i, j, k = 1, 2, …, n; siendo 𝑞𝑘 la variable generalizada k. Por tanto

5. Se definen las matrices de pseudoinercias 𝐉_𝑖 para cada elemento en base a (20).

donde: i = 1, 2, …, n. Las integrales están extendidas al elemento i considerando que (𝑥_𝑖𝑦_𝑖𝑧_𝑖) son las coordenadas cartesianas del diferencial de masa dm respecto al sistema de coordenadas del elemento. En la Figura 2 se señalan los elementos del sistema multicuerpos (torre, eje del rotor y pala).

6. Se obtienen la matriz de inercias 𝐃=[𝑑𝑖𝑗] cuyos elementos se definen por

7. Se calculan los términos ℎ𝑖𝑘𝑚 de la forma

donde i, k, m = 1, 2, …, n. Posteriormente se construye la matriz columna de fuerzas de Coriolis y centrípeta 𝐇=[ℎ_𝑖]𝑇 cuyos elementos se definen como

de forma tal que

8. Finalmente se obtiene la matriz columna de fuerzas de gravedad 𝐂=[𝑐_𝑖]^𝑇cuyos elementos están definidos por:

donde: i,j,m = 1,2,...n. Siendo g el vector de gravedad expresado en el sustema de la base {S₀} expresado por (𝑔_𝑥0, 𝑔𝑦0, 𝑔_𝑧0, 0), y 𝑖_𝑟𝑗 es el vector de coordenadas homogéneas del centro de masas del elemento j expresado en el sistema de referencias del elemento i. Por lo tanto expresado por (𝑔𝑥0,𝑔𝑦0,𝑔𝑧0,0), y 𝑟𝑗𝑖 es el vector de coordenadas homogéneas del centro de masas del elemento j expresado en el sistema de referencias del elemento i. Por lo tanto

Como se observa en (), (), (), () y (), bajo la simplificación de inicial que excluye el movimiento de la torre, la dinámica del sistema estará descrita solo por un par de ecuaciones, de modo tal que en su forma canónica 𝑟 = 𝐃𝑞̈ + 𝐇 + 𝐂 el modelo es expresado como

donde: 𝜏_𝑖 es la fuerza y torque efectivo aplicado a 𝑞_i.

El procedimiento descrito se puede generalizar para las otras dos palas (iguales entre todas ellas), obteniendo de esta forma tres conjuntos de ecuaciones con la misma estructura, las cuales pueden relacionarse a través del ángulo del rotor de la forma siguiente:

por lo tanto

Una vez que se realizan las simplificaciones adecuadas se obtiene el modelo general para el rotor de un aerogenerador descrito en las ecuaciones (39) a la (56), donde 𝑞₄ y 𝑞₅ corresponden a las variables asociadas a la rotación de las palas 2 y 3; mientras m y d, son la masa y el centro de masas (en x) comunes entre todas las palas.

Las simulaciones mostradas a continuación se realizaron bajo el entorno de Matlab® usando una fórmula explicita del método de Runge-Kutta (función ode45). En la Tabla 1 se listan los parámetros empleados. La masa y centros de masa corresponden al modelo en Solidworks® descrito en ^[16], para una pala fabricada con fibra de vidrio con un valor de punta de 33.38 m. En ^[17] se señala, que el uso de la función ode45 demanda una tolerancia de error relativa pequeña para resolver con precisión problemas en intervalos "largos", así como aquellos moderadamente inestables como es el caso del modelo que se presenta en este artículo.

En la Figura 3 se muestra la simulación del modelo de una pala, con condiciones El sistema parte del reposo y con el objeto de observar su respuesta natural el valor de los torques es establecido en cero. Para la posición inicial del eje del rotor (11/12π) la pala asumirá el papel de un péndulo, describiendo para el rotor un comportamiento sinusoidal constante en posición y velocidad, mientras la dinámica de la pala se presenta como una sinusoidal creciente, al no considerar en el modelo ningún efecto de fricción. En la Figura 4a se muestra la simulación de un diagrama de alambre que ilustra la trayectoria de movimiento en el espacio 3D del punto extremo en la pala. En la Figura 4b, se ilustra el movimiento pendular respectivo para el caso en que A→0 y B→0.

En la Figura 5 se muestra la simulación del modelo de tres palas que parte del reposo con entradas cero y condiciones . En este caso, el eje del rotor cambia lentamente (entre 0.5127 y 0.5273 rad, a una velocidad máxima de 0.0458rad/s), mientras las restantes variables muestran un comportamiento oscilatorio. Cuando el sistema no parte del reposo , el eje del rotor muestra mayor variación en posición y velocidad (ver Figura 6), describiendo un movimiento en sentido contrario respecto al caso anterior. En la Figura 7 se observan las diferencias en las trayectorias de movimiento.

En la Figura 8 se observa la respuesta natural del sistema partiendo del reposo cuando una de las palas presenta un valor angular inicial distinto a cero, esto es . En posición, las variables exhiben un comportamiento creciente muy notorio para la pala 1, y poco perceptible para el eje del rotor (menor a 0.5rad).

Para el caso en que más de una pala presente una variación angular inicial, el sistema muestra un comportamiento sinusoidal armónico (ver Figura 9) como producto del relación entre los movimientos angulares de todas las palas y sus efectos sobre el eje del rotor, que es observable en la trayectoria que describen las palas (ver Figura 10).

En la Figura 11 se muestran la simulación derivada de la conjunción de las condiciones iniciales de las pruebas anteriormente descritas, donde el eje del rotor presenta un movimiento angular sinusoidal (con una cresta de 2.73l5 rad), mientras las palas muestran un comportamiento creciente que en combinación llevan a las puntas de las palas a describir trayectorias circulares tal como se muestra en la Figura 13a. En el caso en que todas las variables presenten condiciones diferentes de cero, la dinámica propia de las palas reduce las variaciones angular (ver Figura 12), efecto que es observable en la trayectoria de movimiento (ver Figura 13b).

En la Figura 14 se muestran las trayectorias para diferentes condiciones iniciales. En la Figura 14a la pala 1 que presenta un valor inicial de 1/6π describe un movimiento creciente, mientras que el eje del rotor y las restantes palas describen trayectorias de naturaleza sinusoidal. En la Figura 14b se denotan crecimientos angulares para las tres variables, que se reflejan en remarcadas trayectorias circulares. En la Figura 14c la pala 1 muestra una mayor variación angular a razón de su valor inicial, y el impulso de la velocidad en el eje del rotor. Los efectos de velocidades iniciales en el eje y en alguna de las palas son ilustrados en la Figura 14d.

En ninguno de los casos presentados existió un cambio angular en el eje del rotor mayor a una revolución.

Finalmente en la Figura 15 se muestran los efectos de un torque constante aplicado al rotor, cuando este parte del reposo en equilibrio estático. Como se observa una mayor magnitud de torque genera menores variaciones en posición y velocidad del eje del rotor, al contrarrestar la dinámica de las palas que son libres fricción. Las variaciones en las trayectorias evidencian el impacto que tiene el incremento del torque en el rotor (ver Figura 16).