Modelado Euler-Lagrange del rotor de un aerogenerador tripala como sistema multicuerpo

1. Introducción

Un medio actual para el aprovechamiento de la energía eólica son las turbinas de viento (wind turbines) o aerogeneradores, siendo los de eje horizontal (HAWT, Horizontal-Axis Wind Turbine) con dos o tres palas los más comunes [1]. Uno de los principales componentes en estos sistemas es el rotor o sistema de captación de viento, el cual gira propulsado por las fuerzas que resultan de su interacción con el viento, convirtiendo la energía cinética del viento en energía mecánica útil [2]. El rotor está formado por las palas interconectadas al hub, cubo o buje, siendo deseable su alta eficiencia para una mayor extracción de energía eólica [3].

En lo que respecta al modelado y simulación de aerogeneradores un enfoque estructural explorado es la dinámica de multicuerpo (multibody dynamics), la cual consiste en considerar al aerogenerador como un conjunto de cuerpos rígidos y/o flexibles vinculados entre sí mecánicamente [4]. Sobre esta línea se han desarrollado diversos trabajos pioneros, como es el caso de [5] quien presenta el modelado de un HAWT como una cadena cinemática abierta de nueve componentes, donde la torre está rígidamente unida a la tierra y los cuerpos superiores son libres de moverse en el espacio. Más recientemente, en [6] se propone un procedimiento de diseño para el tren motriz de un aerogenerador basado en el enfoque de sistemas multicuerpo, con el propósito de estimar cargas estáticas y dinámicas empleando multiplicadores de Lagrange. Mientras en [7] se presenta el modelado de un aerogenerador marino flotante (FOWT, Floating Offshore Wind Turbine), donde las ecuaciones del cuerpo rígido acoplado y los movimientos elásticos son derivados de las ecuaciones de Lagrange. De igual manera, en [8] es modelado un FOWT como un sistema multicuerpo conformado por la plataforma flotante, la torre, la nariz, el cubo y tres palas.

En lo que respecta a la simulación, en [9] se realiza la simulación aerolástica de un HAWT empleando la dinámica de multicuerpos para representar la respuesta de movimiento a cargas aerodinámicas. Otros trabajos como el de [10] combinan un análisis multicuerpos con redes neuronales artificiales para reducir el tiempo de ejecución de la simulación de la dinámica de un aerogenerador, mientras se garantiza su precisión. Por otra parte, en [11] se aborda el modelado y simulación tridimensional (3D) de la cinemática de un aerogenerador como un sistema de cinco cuerpos (torre, eje del rotor y tres palas), empleando matrices de transformación derivadas de uso de la convención Denavit-Hartenverg (D-H) que es empleada en el modelado cinemático de robots. Software de análisis y simulación como Adams [12] y HAWC2 [13] también hacen uso del enfoque de multicuerpos.

Con base en lo anterior, en este artículo se desglosa el modelado del rotor de un aerogenerador tripala bajo el enfoque de la dinámica de multicuerpos, empleando una forma algorítmica Euler-Lagrange, y una simplificación de un modelo cinemático mediante matrices de transformación. Se describe el procedimiento para el modelado y su validación a través de una simulación por computadora en 2D y 3D.

2. Modelado y simulación

2.1. Modelado dinámico mediante el algoritmo de Euler-Lagrange

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. 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. 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

Ecuación 1

Ecuación 2

Ecuación 3

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í

Ecuación 4

Ecuación 5

Ecuación 6

Ecuación 7

Ecuación 8

Ecuación 9

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

Ecuación 10

Ecuación 11

Ecuación 12

Ecuación 13

Ecuación 14

Ecuación 15

Ecuación 16

Ecuación 17

Ecuación 18

Ecuación 19

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

Ecuación 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).

Ecuación 21

Ecuación 22

Ecuación 23

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

Ecuación 24

Ecuación 25

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

Ecuación 26

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

Ecuación 27

de forma tal que

Ecuación 28

Ecuación29

Ecuación 30

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

Ecuación 31

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

Ecuación 32

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

Ecuación 33

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:

Ecuación 34

Ecuación 35

Ecuación 36

por lo tanto

Ecuación 37

Ecuación 38

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.

Ecuación 39

Ecuación 40

Ecuación 41

Ecuación 42

Ecuación 43

Ecuación 44

Ecuación 45

Ecuación 46

Ecuación 47

Ecuación 48

Ecuación 49

Ecuación 50

Ecuación 51

Ecuación 52

Ecuación53

Ecuación 54

Ecuación 55

Ecuación 56

2.2 Simulación del modelo en Matlab®

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).

3. Discusión

La idea detrás de modelos genéricos como el presentado en este artículo es la disponibilidad de información para su implementación en diversos entornos de simulación, para lo cual solo es necesario proporcionar una parametrización del sistema real, para este caso la masa y centro de masas de las palas.

Un inconveniente de los modelos genéricos, es la posible reducción de su precisión en comparación con modelos complejos que son provistos por los fabricantes, por lo que encontrar un equilibrio entre precisión y complejidad resulta una difícil tarea. Al respecto en [4] se señala que existe una tendencia en las herramientas de análisis a disminuir el costo computacional, lo que implica el uso de modelos simplificados que si bien difícilmente capta todos los fenómenos, resultan válidos en etapas iniciales de diseño o en el análisis de pequeños aerogeneradores. En [18] se menciona que en el estudio de la simulación eléctrica un modelo aerodinámico detallado suele no ser necesario, por lo que el uso de modelos estructurales como la dinámica de multicuerpos (de menor nivel de desarrollo comparada con los de elemento finito) resulta una alternativa adecuada.

En el caso particular, la simulación del modelo genérico presentado, tiene como propósito ser una herramienta de exploración experimental y conceptual de la dinámica del rotor de un aerogenerador. La simulaciones mostradas ilustran la dinámica del sistema bajo diversas condiciones iniciales, denotando en todas ellas los efectos energéticos del movimiento de las palas sobre la posición y velocidad del rotor, precepto bajo el cual se fundamenta el modelo descrito. La simulación 3D del modelo a través de representaciones en diagramas de alambre, tal como se indica en [19] representa un medio simple para verificar su fiabilidad.

4. Conclusiones

En este artículo se expuso la obtención y simulación de un modelo del rotor de un aerogenerador tripala a partir de la descripción de su cinemática mediante matrices de transformación simplificadas. El modelo abordado desde el enfoque de la dinámica de multicuerpos es construido a partir del modelo presentado en [11] y la formulación Euler Lagrange. El modelo de características genéricas fue validado a través de una simulación 2D y 3D en el entorno de Matlab®. El método para la obtención del modelo, resulta en sí mismo una aportación de este trabajo.

El empleo del método expuesto para conformar un modelo que incluya otros cuerpos del aerogenerador como la torre, así como su aplicación en aerogeneradores FOWT son líneas de trabajo a explorar.

Agradecimientos

Al Fondo Sectorial CONACYT-Secretaría de Energía Sustentabilidad Energética y al proyecto P08: Diseño y construcción de un aerogenerador experimental con capacidad de 3 kW y desarrollo de software de simulación en realidad virtual, con fines didácticos, del Centro Mexicano de Innovación en Energía Eólica (CEMIE-Eólico), por el apoyo económico y tecnológico brindado para la realización de este trabajo.

Referencias

Información adicional

Cómo citar: F. Aguilar, AP. Matus, MA. Hernández, JJ. Arellano, S. Sánchez, D. Pacheco, “Modelado Euler-Lagrange del rotor de un aerogenerador tripala como sistema multicuerpo,” Rev. UIS Ing., vol. 19, no. 1, pp. 25-36, 2020. doi: 10.18273/revuin.v19n1-2020002

Revista UIS Ingenierías
ISSN: 1657-4583
Vol. 19
Num. 1
Año. 2020

Modelado Euler-Lagrange del rotor de un aerogenerador tripala como sistema multicuerpo

FranciscoAna PatriciaMiguel ÁngelJ. JesúsSergioDaniel Aguilar-AcevedoMatus-VicenteHernández-LópezArellano-PimentelSánchez-SánchezPacheco-Bautista
Universidad del IstmoUniversidad del IstmoUniversidad del IstmoUniversidad del IstmoUniversidad del IstmoUniversidad del Istmo,MéxicoMéxicoMéxicoMéxicoMéxicoMéxico